Dispersive shocks and spectral analysis for linearized Quantum Hydrodynamics

TL;DR

This paper analyzes the spectral properties and stability of dispersive shocks in a 1-D quantum hydrodynamics model, focusing on quantum effects via the Bohm potential and using spectral analysis techniques.

Contribution

It provides a spectral analysis of linearized quantum hydrodynamics with dispersive effects, including stability assessment using the Evans function method.

Findings

Spectral stability of dispersive shocks established.

Linearized operator computed explicitly.

Quantum effects influence shock stability characteristics.

Abstract

In this paper we perform the analysis of spectral properties of the linearized system around constant states and dispersive shock for a 1-D compressible Euler system with dissipation--dispersion terms. The dispersive term is originated by the quantum effects described through the Bohm potential, as customary in Quantum Hydrodynamic models. The analysis performed in this paper includes the computation of the linearized operator and the spectral stability through the Evans function method.

Equations

$$\left\{\begin{array}[]{ll}\rho_{t}+m_{x}=0,\\ m_{t}+\Big{(}\frac{m^{2}}{\rho}+\rho^{\gamma}\Big{)}_{x}=\epsilon\mu m_{xx}+\epsilon^{2}k^{2}\rho\Big{(}\frac{(\sqrt{\rho})_{xx}}{\sqrt{\rho}}\Big{)}_{x},\end{array}\right.$$

$$\rho=P\Big{(}\frac{x-st}{\epsilon}\Big{)}\mbox{, }m=J\Big{(}\frac{x-st}{\epsilon}\Big{)}\mbox{. }$$

$$\underset{y\to \pm \infty }{\lim }P(y)=P^{\pm }and\underset{y\to \pm \infty }{\lim }J(y)=J^{\pm }$$

$$J^{+}-J^{-}$$

$$\displaystyle\Big{(}\frac{J^{2}}{P}+P^{\gamma}\Big{)}^{+}-\Big{(}\frac{J^{2}}{P}+P^{\gamma}\Big{)}^{-}$$

$$\rho=P\Big{(}\frac{x-st}{\epsilon}\Big{)}\mbox{, }m=J\Big{(}\frac{x-st}{\epsilon}\Big{)}$$

$$\rho\Big{(}\frac{(\sqrt{\rho})_{xx}}{\sqrt{\rho}}\Big{)}_{x}=\frac{1}{2}\Big{(}\rho(\ln\rho)_{xx}\Big{)}_{x}\mbox{. }$$

$$-s{P}^{\prime }+J^{\prime }$$

$$\displaystyle-sJ^{\prime}+\Big{(}\frac{J^{2}}{P}+P^{\gamma}\Big{)}^{\prime}$$

$$-s{\int }_{-\infty }^y{P}^{\prime }(x)dx+{\int }_{-\infty }^y{J}^{\prime }(x)dx=0,$$

$$J(y)-sP(y)=J^{-}-s{P}^{-}.$$

$$J(y)-sP(y)=J^{+}-s{P}^{+},$$

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

$$A=s{P}^{\pm }-J^{\pm },$$

$$P^{\prime \prime }=\frac{2}{k^2}f(P)-\frac{2s\mu }{k^2}{P}^{\prime }+\frac{P^{\prime 2}}{P}.$$

$$f(P)$$

$$=P^{\gamma }-(As+B)+\frac{A^2}{P}$$

$$B=-sJ^{\pm}+\Big{(}\frac{J^{2}}{P}+P^{\gamma}\Big{)}^{\pm},$$

$$P^{\prime }=Q$$

$$Q^{\prime }=\frac{2f(P)}{k^2}-\frac{2\mu s}{k^2}Q+\frac{Q^2}{P}$$

$$f(P)$$

$$-\frac{(P^{+}{)}^{\gamma +1}-(P^{-}{)}^{\gamma +1}}{P^{+}-P^{-}}.$$

$$F(P)=\frac{P^{\gamma +1}}{\gamma +1}-(As+B)P+A^2\ln P,$$

$$\displaystyle F_{1}(P)=F(P)-\Big{(}\frac{s\mu P}{k}\Big{)}^{2}-F(P^{-}),$$

$$\displaystyle\tilde{F}_{1}(\xi)=F(\xi)-\Big{(}\frac{s\mu\xi}{k}\Big{)}^{2}-F(P^{+})$$

$$\frac{s\mu }{k}<\sqrt{-2f^{\prime }(P^{+})},$$

$$-\frac{s\mu }{k}<\sqrt{-2f^{\prime }(P^{-})},$$

$$J=\left[\begin{array}{cc}0& a\\ f^{\prime }(P^{-})& -b\end{array}\right].$$

$${\mu }_{1,2}(P^{-})=(-b\pm \sqrt{b^2+4a{f}^{\prime }(P^{\pm })})/2$$

$$v_2=-\left[\begin{array}{c}\frac{(b+\sqrt{b^2+4a{f}^{\prime }(P^{-})})}{2f^{\prime }(P^{-})}\\ 1\end{array}\right]$$

$$L\left[\begin{array}{c}\tilde{\rho }\\ \tilde{J}\end{array}\right]=\left[\begin{array}{c}s{\tilde{\rho }}_y-{\tilde{J}}_y\\ s{\tilde{J}}_y+(\frac{J^2}{P^2}\tilde{\rho }{)}_y-(\frac{2J}{P}\tilde{J}{)}_y-\gamma (P^{\gamma -1}\tilde{\rho }{)}_y+\mu {\tilde{J}}_{yy}+L_V\tilde{\rho }\end{array}\right],$$

$$L_{V}\tilde{\rho}=\frac{k^{2}}{2}\tilde{\rho}_{yyy}-2k^{2}\Big{(}(\sqrt{P})_{y}\Big{(}\frac{\tilde{\rho}}{\sqrt{P}}\Big{)}_{y}\Big{)}_{y},$$

$$\lambda \left[\begin{array}{c}\tilde{\rho }\\ \tilde{J}\end{array}\right]=L\left[\begin{array}{c}\tilde{\rho }\\ \tilde{J}\end{array}\right].$$

Full text

1. Introduction

The aim of this paper is to study the following Euler system with dissipation–dispersion terms:

[TABLE]

where $t\geq 0$, $x\in\mathbb{R}$, $\rho=\rho(t,x)>0$, $m=m(t,x)$, $\gamma\geq 1$, $0<\epsilon\ll 1$, $\mu>0$, $k>0$. The positive coefficients $\epsilon\mu$ and $\epsilon^{2}k^{2}$ stand for the viscosity and dispersive coefficients, respectively, and $\rho^{\gamma}$ is the pressure. The particular shape of the dispersion terms is due to the Bohm potential, and the resulting system is referred to as the quantum hydrodynamics system, being used for instance in superfluidity or to model semiconductor devices. We are in particular interested in the analysis of the linearized version of (1) around special solutions, as constant states and dispersive shocks, namely solutions of (1) written as travelling waves

[TABLE]

As customary, the speed $s\in\mathbb{R}$ of the travelling wave and its limiting end states

[TABLE]

are assumed to verify the Rankine–Hugoniot conditions:

[TABLE]

The existence of such solution is studied in full details in the companion paper [19], where in particular the interplay between the diffusive and dispersive effects is analyzed in connections with the existence, monotonicity, and stability of such solutions. The notion and study of the effect of dispersive terms have been first considered by [13, 21], see also [12, 14, 20, 23], while a fairly complete analysis, via the Whitham modulation theory, has been investigated in [15], which also includes a wide bibliography on these topics. The first attempt to analyze the spectral theory of the linearized operator around dispersive shocks has been discussed in [16] regarding the case of $p$-system with real viscosity and linear capillarity, but only in the case of monotone shocks, while the mathematical theory of the quantum hydrodynamic systems has been developed in [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11].

In the present paper, we shall set up the dynamical systems solved by these profiles, obtained as heteroclinic connections, and give some numerical computations of them in the next two sections. Moreover Section 4 is devoted to the study of the linearization of (1), both around a profile and around a constant state, the latter case being relevant for the subsequent analysis, which strongly relies on the linearized, constant coefficient, asymptotic operators, obtained as the parameter of the profile tends to $\pm\infty$. In Section 5 we treat the analysis of the spectrum of linearized operators, in particular giving a resolvent estimate for the constant coefficient case, namely, for the linearized operator about a constant state. To investigate the point spectrum of the linearized operator around a profile, one needs to supplement the resolvent estimate (available for large eigenvalues) with the study of the Evans function (see, for instance [22]); the corresponding numerical results about it are collected in Section 5.2.

2. The equations for the profile

In order to analyze travelling wave profiles, we plug

[TABLE]

in (1), and we rewrite the Bohm potential in conservative form as follows:

[TABLE]

After substituting the profiles $P$ and $J$ in the system (1) and multiplying by $\epsilon$ we obtain

[TABLE]

where ′ denotes $d/dy$ and $P=P(y)$, $J=J(y)$. Integrating equation (4), we get

[TABLE]

that is

[TABLE]

We can also integrate (4) from $y$ to $+\infty$ to get

[TABLE]

that is

[TABLE]

with

[TABLE]

thanks to the Rankine-Hugoniot condition (2).

Similarly, we can integrate (5) up to $\pm\infty$, after substituting the expression for $J(y)$ into it, to obtain the planar ODE

[TABLE]

In (7),

[TABLE]

and the constant $B$ is given by

[TABLE]

as follows from (3).

Introducing the new variable $P^{\prime}=Q$, we rewrite (7) as a first-order system of ODEs

[TABLE]

Finally, it is worth to observe that the constants $A,B$ in $f(P)$ can be expressed solely in terms of $P^{\pm}$ as follows:

[TABLE]

3. Numerical computation of the profiles

The existence of profiles under different assumptions on the asymptotic states is proved in [19]; here we collect some numerical computations about such profiles in the parameters’ range where the existence is proved. In what follows, we briefly recall the case of large shocks, which in particular includes the possibility of non monotone profiles.

To this aim, let us denote

[TABLE]

and

[TABLE]

Then the next theorem holds.

Theorem 1**.**

Let $s>0$ and $P^{+}<P^{-}$. If $F_{1}(P)\leq 0$ for $0<P<P^{+}$, then there is a traveling wave profile, connecting $[P^{-},0]$ to $[P^{+},0]$. If in addition

[TABLE]

*then the traveling wave profile is non-monotone.

Let $s<0$ and $P^{-}<P^{+}$. If $\tilde{F}_{1}(\xi)\leq 0$ for $0<\xi<P^{-}$, then there is a traveling wave profile, connecting $[P^{-},0]$ to $[P^{+},0]$. If in addition*

[TABLE]

then the traveling wave profile is non-monotone.

The condition $s>0$, $P^{+}<P^{-}$ is implied by the condition for a 2-shock with a supersonic right state, and $u^{+}>c_{s}(P^{+})$, where $u^{+}=J^{+}/P^{+}$, and $c_{s}(P^{+})$ is the speed of sound.

We describe the result of Theorem 1 in the particular case of $\gamma=1$, in Figure 1 (a), and $\gamma=\frac{3}{2}$, in Figure 1 (b). The condition of Theorem 1 can be easily checked numerically for the profiles presented in Figure 1. To give a numerical computation for such profile, we analyze the jacobian of (9)-(10) at the steady-state $[P^{-},0]$. With the notation $a=\frac{2}{k^{2}}>0$, $b=\frac{2s\mu}{k^{2}}>0$, the jacobian is given by

[TABLE]

a direct computation shows that his eigenvalues are

[TABLE]

and, since $f^{\prime}(P^{-})>0$, we have $\sqrt{b^{2}+4af^{\prime}(P^{-})}>b$, and we conclude $\mu_{1}(P^{-})<0<\mu_{2}(P^{-})$ and therefore the steady-state $[P^{-},0]$ is a saddle. The eigenvector corresponding to the positive eigenvalue $\mu_{2}(P^{-})$ with negative second component is

[TABLE]

and, as customary, it is tangent to the unstable subspace of the equilibrium $[P^{-},0]$. Then, to numerically compute the traveling wave profile, we numerically integrate (9)-(10) with the initial condition $[P^{-},0]+hv_{2}/|v_{2}|$, with $h>0$, sufficiently small, for instance $h=10^{-5}$; see Figure 1. The global existence Lemma from [19] applies to the parameters of these profiles. The jacobian at $[P^{+},0]$ for the parameters in Figure 1 (b) has imaginary eigenvalues, which shows that the profile is non-monotone.

Using $P(y)$, we obtain $J(y)$ from equation (6).

4. Linearization

Using the change of variables $\tau=t/\epsilon$, $y=(x-st)/\epsilon$, we get the full linearized operator around the profile for (1):

[TABLE]

where

[TABLE]

with associated eigenvalue problem given by

[TABLE]

For the analysis of the eigenvalue problem (12) we will need the Evans function and to this end it is also important to re–express the above linearized systems 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 [16]. To this end, consider

[TABLE]

Integrating the equation (12) it follows that for $\lambda\neq 0$ the integrated variables $\tilde{\rho}$ and $\tilde{J}$ decay exponentially as $|x|\rightarrow+\infty$. Expressing $\tilde{\rho}$ and $\tilde{J}$ in terms of $\hat{\rho}$ and $\hat{J}$, and integrating (12) from $-\infty$ to $x$ we get the system in integrated variables:

[TABLE]

with

[TABLE]

We can rewrite (13)-(14) as $V^{\prime}=\hat{M}(x,\lambda)V$, where $V=[\hat{\rho},\hat{J},u_{1},u_{2}]^{T}$ and

[TABLE]

The limit of $\hat{M}(x,\lambda)$ as $x\rightarrow\pm\infty$ is given by

[TABLE]

which is the same asymptotics one shall obtain from the original linearization (11) after rewriting the latter as a first order $4\times 4$ system. This will lead to a linear, constant coefficient operator which is crucial in the spectral analysis of the linerazed operator around the profile.

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

[TABLE]

Then the operator, corresponding to the linearization around the constant steady–state $(\bar{P},\bar{J})$ is

[TABLE]

where ${}^{\prime}=d/dy$. The asymptotic operators at $\pm\infty$ for (11) are given by

[TABLE]

where

[TABLE]

We may rewrite the equation

[TABLE]

as a first order system $V^{\prime}=M^{\pm}V$, with $V=[\tilde{\rho},\tilde{J},u_{1},u_{2}]^{T}$, and the matrices $M^{\pm}$ are exactly (16).

5. Spectral analysis

In this section we shall locate the spectrum of our linearized operators, starting from the constant steady–case. Indeed, there are no eigenfunctions which decay at $\pm\infty$ for the constant coefficient linear operator $L_{c}$, thus the spectrum of the latter reduces to the essential one. In addition, we shall perform a resolvent estimate in that case valid in the whole unstable half plane $\Re(\lambda)>0$, to present in particular a simple calculation which is useful to understand the behavior of the spectrum for the linearized operator about the profile $L$; more details in this direction are presented in [19].

5.1. Constant steady state

To obtain a resolvent estimate in $H^{1}$, we need to invert in that space the relation

[TABLE]

and prove $(\lambda I-L_{c})^{-1}:H^{1}\to H^{1}$ is bounded for $\Re(\lambda)>0$. With the notation (17), the system at hand rewrites

[TABLE]

and denote $f=[f_{1},f_{2}]^{T}$.

The speed of sound is $c_{s}(\bar{P})=\sqrt{\gamma\bar{P}^{\gamma-1}}$ and the flow velocity is $\bar{u}=\bar{J}/\bar{P}$. The condition for a subsonic steady-state is $|\bar{u}|<c_{s}(\bar{P})$, which becomes $\alpha<0$ after squaring. As an example, we present here below a resolvent estimate in the above framework in $H^{1}$.

Proposition 2**.**

If the steady-state is subsonic, then for $\Re(\lambda)>0$ we have the resolvent estimate

[TABLE]

where $C=\max\{\frac{k^{2}}{2},-\alpha,1\}$ and $h(\Re(\lambda)):=\min\{-\alpha\Re(\lambda),\Re(\lambda),\frac{k^{2}}{2}\Re(\lambda),\mu\}$.

Proof.

We multiply equation (18) by $-\alpha\overline{\tilde{\rho}}$ and equation (19) by $\overline{\tilde{J}}$. After taking the real part we get

[TABLE]

Let $\tilde{\rho}=\tilde{\rho_{r}}+i\tilde{\rho_{i}}$ and $\tilde{J}=\tilde{J_{r}}+i\tilde{J_{i}}$. We have

[TABLE]

Therefore

[TABLE]

After integration, the derivative term on the right hand side of (20) disappears and we get

[TABLE]

By integration by parts we obtain

[TABLE]

We have

[TABLE]

hence by integration by parts

[TABLE]

Let $f_{1}=f_{1,r}+if_{1,i}$. Substituting $\tilde{J_{r}}^{\prime}$ and $\tilde{J_{i}}^{\prime}$ from equation (18) yields

[TABLE]

where

[TABLE]

Using integration by parts we get

[TABLE]

and again by integration by parts

[TABLE]

which yields

[TABLE]

Substituting (22) and (23) into (21) we get

[TABLE]

Using integration by parts we get

[TABLE]

Equation (25) yields

[TABLE]

With the notation $h(\Re(\lambda)):=\min\{-\alpha\Re(\lambda),\Re(\lambda),\frac{k^{2}}{2}\Re(\lambda),\mu\}$, we estimate the left-hand side of (24) by

[TABLE]

Let $C=\max\{\frac{k^{2}}{2},-\alpha,1\}$. Taking into account (26) we estimate the right-hand side of (24) by

[TABLE]

From inequalities (27) and (28) we get

[TABLE]

This yields the proof of the proposition. ∎

For completeness, we underline that the stability of the essential spectrum can be obtained by Lemma 3 below, which is a restatement for the present case involving $L_{c}$ of the corresponding result in [19] for the asymptotic operators $L_{\pm\infty}$. The statement refers to the asymptotic operators at $x\to\pm\infty$, as it is used to locate the essential spectrum also for the linearized operator around the profile [17], but clearly can be applied to any linearized operator about a constant state. To this end, we report the dispersion relation referred to $L_{c}$ here below:

[TABLE]

Lemma 3**.**

If $\alpha\leq 0$, then the essential spectrum is in the closed left half-plane. Also, if $\xi\neq 0$, then $\Re(\lambda_{1,2})<0$.

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, and in particular the resolvent estimate obtained above is proved in the companion paper [19] only for $\Re(\lambda)$ 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 $|\lambda|$, which gives a quantitative version of the asymptotic results of [22], excluding the presence of eigenvalues for $|\lambda|>C$ with an explicit bound for the constant. This result is completed with the numerical study of the winding number of the Evans function in a sufficiently big contour; all details can be found in [19]. In the next section we shall briefly recall the main ingredients/results for the latter numerical analysis, thus yielding to a numerical evidence of point spectrum stability.

5.2. The Evans function and its numerical evaluation

To define the Evans function, let us consider the equation $Y^{\prime}=\hat{M}(y,\lambda)Y$, where $\hat{M}(y,\lambda)$ is defined in (15). As it is manifest, its limits at $\pm\infty$ are given by the matrices $M^{\pm}$, defined by (16), corresponding to limit states $P^{\pm}$, 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 $\nu^{-}_{1},\dots,\nu^{-}_{k}$ (i.e. $\Re(\nu^{-}_{i})>0$), and $M^{+}$ has $n-k$ stable eigenvalues $\nu^{+}_{1},...,\nu^{+}_{n-k}$ (i.e. $\Re(\nu^{+}_{i})<0$), and denote the corresponding (normalized) eigenvectors by $v^{\pm}_{i}$. In our case $n=4$ and $k=2$. Let $Y^{-}_{i}$ be a solution of $Y^{\prime}=M(y,\lambda)Y$, satisfying $exp(\nu^{-}y)Y^{-}(y)$ tends to $v^{-}_{i}$ as $y\rightarrow-\infty$ and $exp(\nu^{+}y)Y^{+}(y)$ tends to $v^{+}_{i}$ as $y\rightarrow+\infty$. Then, the Evans function can be defined by

[TABLE]

As a consequence, a point $\lambda\in\mathbb{C}$ is in the point spectrum of $L$ if and only if $E(\lambda)=0$. In contrast, to compute the Evans function numerically, we use the compound matrix method; for instance, see [16] . This method is used in order to get a stable numerical procedure, because the system $Y^{\prime}=M(y,\lambda)Y$ is numerically stiff. Specifically, the compound matrix $B(y,\lambda)$ is given by:

[TABLE]

We integrate the equation $\phi^{\prime}=(B(y,\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. Given a profile $[P(y),J(y)]$ at discrete points $(y_{k})_{k=1}^{N}$, we obtain the matrix $B(y,\lambda)$, by linear interpolation. Similarly we integrate the equation $\phi^{\prime}=(B(y,\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:

[TABLE]

Specifically, for our calculations we choose $L_{1}=40$.

Under the assumption that $E(\lambda)$ is analytic in the region surrounded by a closed contour $\Gamma$, and it does not vanish on the contour, we can use the winding number

[TABLE]

to count the number of zeros inside the contour. The remaining part of this paper is devoted to provide numerical evidence that the integral (29) is indeed zero in a sufficiently big contur $\Gamma$ lying in the unstable half plane. We present the calculations on Section 5.3 in integrated variables.

5.3. Numerical evaluation of $E(\lambda)$ for integrated variables

In this subsection we compute the Evans function for parameters $A=1$, $B=1.1$, $s=1$, $\gamma=3/2$, $\mu=1$, $k=\sqrt{2}$. These parameters in particular correspond to a non-monotone shock, and Theorem 1 applies to them. Here we use integrated variables to avoid the smallness of $E(\lambda)$ near zero, namely, we solve the ODEs $\phi^{\prime}=(\hat{B}(y,\lambda)-\mu^{\pm})\phi$, where the compound matrix $\hat{B}$ is constructed from $\hat{M}$ defined in (15).

In the present situation, the essential spectrum touches the origin $\lambda=0$ and therefore, to rigorously use the Evans function tool in order to locate eigenvalues with small positive real part, one should first extend it across the essential spectrum. This rigorous analysis goes beyond the aim of the present discussion; however we numerically check the Evans function is well defined and different from zero at $\lambda=0$ by showing that it is almost constant (and non zero) on a small semi-circular contour without a vertical segment in the unstable half-plane with radius $10^{-6}$ and with center at $\lambda=0$; see Figure 2.

To compute the initial conditions we integrate the reduced Kato ODE

[TABLE]

with the algorithm from [24], that is $|r^{1}|=1$ eigenvector as before (referring to maximal/minimal decay/growt rate of $B^{\pm}$) and for $k>0$

[TABLE]

Using a careful analysis of the Evans function at infinity one can show numerically that there are no eigenvalues for $|\lambda|\geq 1.9\cdot 10^{4}$. Here we complement this information with a numerical evidence of absence of eigenvalues inside that circle.

We have to initialize the computation on the real axis and for stability reasons for a not very large values of $\lambda$. For this we use two contours. One semi-circular contour with radius 10, center at $\lambda=0$ and vertical segment on the imaginary axis. We don’t evaluate the Evans function at [math], but evaluate it till $i10^{-6}$. Along this contour we integrate the Kato ODE using $10^{4}$ points. We also use another countour which surrounds a semi-annular region in the right half-plane with two semi-circles with radii $5$ and $1.9\cdot 10^{4}$, center at $\lambda=0$ and vertical segment on the imaginary axis. Along this contour we use $10^{7}$ points to integrate the Kato ODE. We use higher density near the origin.

Then using these initial conditions we compute $E(\lambda)$ with the stiff solver ode15s in matlab, with relative tolerance $10^{-4}$ and we set $L_{1}=40$. Finally, we apply the symmetry of $E(\bar{\lambda})=\overline{E(\lambda)}$. The union of the areas, surrounded by the two contours covers the whole region of the unstable half-plane inside $|\lambda|\leq 1.9\cdot 10^{4}$. The Evans function $E(\lambda)$ is plotted in Figure 3 and its winding number is (approximately) 0, giving a numerical evidence of point spectrum stability.

Finally, we note that if $\lambda\in\mathbb{R}$, then $E(\lambda)\in\mathbb{R}$. Our numerics agrees with this simple observation, because we get $\Im{E(\lambda)}\approx 0$ for $\lambda=1.9\cdot 10^{4}$. Moreover, to corroborate this accuracy, we also use the Cauchy integral formula

[TABLE]

for $a$ inside the contour $\Gamma$ surrounding the semi-annular region, and, for $a=1.9\cdot 10^{4}-20$, we get a relative error less than $5\cdot 10^{-4}$.

Moreover, we present a computation of the Evans function along a contour, surrounding a semi-annular region with radii $10^{-6}$ and $1.9\cdot 10^{4}$, center at $\lambda=0$ and vertical segment on the imaginary axis. Along the contour we integrate the Kato ODE with $10^{7}$ points. The winding number of the Evans function is zero. See Figure 4. Using the Cauchy integral formula with $a=1.9\cdot 10^{4}-20$, we observe a relative error less than $5\cdot 10^{-4}$.