On the eigenvalues of spectral gaps of matrix-valued Schrödinger operators

This paper presents a method for calculating eigenvalues lying in the gaps of the essential spectrum of matrix-valued Schrödinger operators. The technique of dissipative perturbation allows eigenvalues of interest to move up the real axis in order to achieve approximations free from spectral pollution. Some results of the behaviour of the corresponding eigenvalues are obtained. The effectiveness of this procedure is illustrated by several numerical examples.


Introduction
This paper is a sequel to two articles of the second author, written over a period of more than 25 years. The first, [12], presents numerical methods for self-adjoint Sturm-Liouville-type equations with matrix coefficients, while the second, [13], analyses the application of a dissipative barrier scheme to a Schrödinger equation on a half-line. In the years since [12] appeared, there has been a lot of activity on Schrödinger-type equations with matrix-valued coefficients-see, e.g., Clark and Gesztesy [4] and Clark, Gesztesy, Holden and Levitan [5], together with the substantial bibliographies therein. Some results from the scalar case carry across to the matrix case in a straightforward way; some require new proofs; and some are simply no longer true. As a simple example, the usual spectral data only determines the Titchmarsh-Weyl coefficient, and hence the matrix-valued potential, up to unitary equivalence. Our concern in this article is to examine which of the results in [13] are still true in the case of a matrix-valued potential, and which not. At the end of the article, we indicate some directions for future work on quantum waveguide problems.
To fix notation, we consider the dissipative matrix Schrödinger equation on the half-line [0, ∞), with a regular self-adjoint-type boundary condition at the origin. The precise form of this condition is not important for our results, so we shall use where α ∈ [0, π), even though this is not the most general form. Here u is a vectorvalued function in a subspace of L 2 ([0, ∞)) n , the parameter γ is a non-zero real, and the coefficients Q, S satisfy the following hypotheses: (A1): Q(x) is a Hermitian-valued, integrable over compact subsets of [0, ∞), and is eventually periodic with period a > 0 i.e, there exists R 0 ≥ 0 such that (A2): S is a cutoff function with support in [0, R] for some R ≥ R 0 : there exists 0 < c < 1 such that When x ∈ (cR, R), we assume that elements of S are measurable and their values lie in [0, 1]. The hypothesis (A2) in particular is stronger than is really needed for most of our results, but sufficient to analyse most dissipative barrier schemes.
We define an operator L 0 by: with domain: Our aim is to present a substantial analysis of the interaction between the dissipative barrier approach to the problem of numerical approximation of the spectrum of L 0 , and interval truncation methods. Our methods will be based on the Floquet theory and Weyl-Titchmarsh functions.

Summary of results
We investigate the following results for an eigenvalue λ γ of our problem with the dissipative term iγ S(·) which develops from the eigenvalue λ when γ = 0 : 1. For our non-truncated problem, if λ is an isolated eigenvalue of multiplicity ν, where 1 ≤ ν ≤ n, Q is a compactly supported perturbation of a Hermitian periodic function, and for a sufficiently small γ > 0, the approximation λ γ,j for each j = 1, . . . , ν satisfies the bound: where C 1 and C 2 are positive constants. 2. If our problem is truncated to some interval [0, X], X > R, then by imposing a boundary condition at x = X, any eigenvalue λ γ,X,good of the truncated problem which converges to λ γ as X → ∞ satisfies: where C 3 and C 4 are positive constants and depend on λ γ . Moreover, the total algebraic multiplicities of all λ γ,X,good converging to λ γ is equal to the algebraic multiplicity of λ γ . 3. If our problem is truncated to some interval [0, X], X > R, then by imposing a boundary condition at x = X, then any eigenvalue λ γ,X,bad of the truncated problem which converges as X → ∞ to a point which is not in the spectrum of where C 5 and C 6 are positive constants. 4. One crucial difference between the operators considered here and the scalarcoefficient operators in [13] concerns the behaviour of eigenvalues of L 0 + iγ S as γ 0. In [13], it is shown that if an eigenvalue λ γ of L 0 + iγ S evolves continuously to become an interior point of a spectral band of L 0 + iγ S when γ = 0, then a threshold effect occurs: there exists γ crit > 0 such that as γ γ crit , λ γ converges to the interior point of the spectral band. We give an example to show that this is not generally true for the operators considered in the present article.
The paper is organised as follows. Section 3 is devoted to the representation of the Floquet theory and Glazman decomposition for the main equation. Section 4 and Section 5 contain the analysis of the method for the truncated and non-truncated problem respectively. Finally, Section 6 represents some numerical experiments to illustrate our results.

Floquet theory and Glazman decomposition for matrix Schrödinger operators
The essential spectrum of L 0 can be described using the Floquet theory, studying the solutions of the differential (1) over just one period. We shall review the elements of the Floquet theory below, primarily to introduce the notation which we require for an analysis of the domain truncation technique.
For the point spectrum of L 0 + iγ S, we shall apply the Glazman decomposition technique [1, Appendix 2].
Recall the parameter R > 0 from hypothesis (A2). For a fixed λ ∈ C and any non-zero constant vector h, consider the following two boundary value problems: If these problems can be solved uniquely for every h 1 , then the maps v(R) −→ v (R) ∈ C n and w(R) −→ −w (R) ∈ C n are linear operators (Weyl-Titchmarsh operators or Dirichlet to Neumann maps) which admit representations by n × n matrices: . The matrices M left (λ) and M right (λ) are analytic functions of λ on suitable subsets of C. In particular M left is meromorphic with poles in C + when γ > 0, and on the real axis when γ = 0; the function M right is analytic outside Spec ess (L 0 ) except at a set of poles (see [4,9]).
If λ is an eigenvalue of L 0 + iγ S, then there exists an eigenfunction u; we can Assuming that u(R) is not zero, this leads to the condition that ker(M left (λ) + M right (λ)) is not trivial. In fact, if u(R) were zero, then both M left (·) and M right (·) would be undefined at λ, so the condition that u(R) be non-zero is satisfied automatically if M left (λ) and M right (λ) are well defined.
Conversely, suppose there exists μ ∈ C such that Take h ∈ ker(M left (μ) + M right (μ)) and define a non-trivial vector u by: where v and w are the solutions of P left and P right respectively for the case λ = μ. Then u is a solution for the differential equation −u +(Q+iγ S)u = μu on both the intervals (0, R) and (R, ∞), which is continuous. Moreover, it has a continuous first derivative at x = R which follows from (12) since h ∈ ker(M left (μ) + M right (μ)) and using the definitions of M left (μ) and M right (μ). This means that u is an eigenfunction of L 0 + iγ S with eigenvalue μ. We have proved the following result.
where U is the n × n matrix-valued solution of the initial value problem In order to find M right , we bear in mind that the Q(x) is periodic for x ≥ R 0 and hence for x ≥ R ≥ R 0 ; also the dissipative perturbation S(x) = 0 for x ≥ R. Therefore, we can apply the Floquet theory [6] for the system of differential equations. We rewrite (1) as a first-order differential system: Let (x, λ) be the fundamental matrix of this equation, i.e.
Suppose that A(λ) has a canonical Jordan form, i.e.
Proof 1. This statement holds because det(A(λ)) = 1, which follows from the fact that det( (x, λ)) is a non-zero constant (the coefficient matrix on the right-hand side of (16) having zero trace). 2. The (1) is in the limit-point case at infinity, and hence for λ outside the essential spectrum it has precisely an n-dimensional space of solutions in L 2 (0, ∞). None of the Floquet multipliers 1 (λ), . . . , 2n (λ) has absolute value 1, for otherwise it is possible to construct a Weyl singular sequence of oscillatory solutions from the corresponding Floquet solution; this is impossible as λ lies outside the essential spectrum. Thus, precisely n of the Floquet multipliers must have absolute value strictly less than 1, precisely n have absolute value strictly greater then 1, and we can order them so that It follows from Lemma 2 that if λ lies outside the essential spectrum, then the Jordan matrix J (λ) in (18) decomposes into n × n blocks as where J 1 (λ) corresponds to 1 (λ), . . . , n (λ) and J 2 (λ) corresponds to n+1 (λ), . . ., 2n (λ). In this case, the matrix has columns which span the n-dimensional space of square-integrable solutions; moreover, for N ∈ N : We can partition (x, λ) as where (x, λ) is an n × n solution of whose columns (as we mentioned above) span the space of all square integrable solutions of (1). Hence, by direct verification, if (R, λ) is invertible, then the function We immediately have the following corollary to Lemma 1.

Approximation of spectral-gap eigenvalues using truncated problems
In this section, we truncate the problem over [0, ∞) to a problem on [0, X] for some X > R. At x = X we impose, for some β ∈ R, a self-adjoint artificial boundary condition The operator L 0 is replaced by L 0,X defined by: with domain: In this case, the spectra of L 0,X and L 0,X + iγ S are purely discrete. To characterise the eigenvalues of L 0,X + iγ S, we replace P right in (11) by: Let (x, λ) be a 2n × n matrix of non-L 2 (0, ∞) solutions corresponding to the eigenvalues n+1 , . . . , 2n of (17), thus, where J 2 (λ) is the corresponding Jordan matrix introduced in (19). Moreover, for N ∈ N : We may partition (λ, x) as construct the matrix: choosing C X (λ) so that cos(β) X (X)−sin(β) X (X) = 0: (32) Hence, the solution w of the boundary value problem (28) exists if the corresponding Thus, the eigenvalues of L 0,X + iγ S may be characterised by an analogue of Lemma 1 if we replace M right by Corollary 1 also has the following analogue.
We analyse the effect of interval truncation through a sequence of intermediate results and technical lemmas. Furthermore, in particular, where C is a positive constant, Proof Equation (35) follows directly from (20, 29) upon using the definitions (21, 30). Substituting (35) into (32) yields (36). In order to estimate the norm of C X (λ), we need to find the norm of J 1 (λ) N . The expression for J 1 (λ) N may have terms, of worst case scenario, like: A similar analogue would be followed to estimate the norm of J 2 (λ) −N . Therefore, Finally, (37) follows from the expressions of the norms of J 1 (λ) N and J 2 (λ) −N . (22) and suppose (λ) > 0. Then (R, λ) is invertible.

Lemma 4 Let (x, λ) be defined as in
Proof Assume that (R, ·) is not invertible: then we can find a non-zero vector c ∈ C n such that (R, λ)c = 0. Define a function u(x, λ) = (x, λ)c, which satisfies the differential equation for all x. Also, u(R, λ) = 0. However, the fact that with Dirichlet condition at R. This problem is self-adjoint, so (λ) = 0, which is a contradiction.
is an eigenfunction of the problem on [R, X] with Dirichlet condition at R and the boundary condition (25) at X with β ∈ R. This is also a self-adjoint problem, so again we have the contradiction (λ) = 0.
Finally, we have a quantitative lemma on continuity of determinants, which will be needed in the proof of Theorem 1. We shall use this lemma with X = R + Na and large N.
where b is a positive constant and 0 < τ < 1. Then whereb is a positive constant.
Then there exist approximations λ γ,X,good to λ γ , whose total algebraic multiplicity is equal to the algebraic multiplicity of λ γ , obtained as eigenvalues of the operator L 0,X + iγ S defined in (26, 27), which satisfies Here C 3 and C 4 are positive constants which depend on λ γ .
Proof Without loss of generality, it is sufficient to check the cases X = R + Na where N ∈ N is sufficiently large. The other cases follow by exploiting the freedom in the choice of the constants c and R in (26) and (27). For example, if X = R + Na + b, with 0 < b < a, then we can replace R by R + b and use a smaller constant c in (26). First, we observe that for γ > 0, λ γ has strictly positive imaginary part. If u γ is the corresponding normalised eigenfunction, then a standard integration by parts yields: is the Hermitian conjugate of u γ (x). Next, we observe consequently from Lemma 4 and Lemma 5, (R, ·) and X (R, ·) are invertible for all λ in a neighbourhood of λ γ . It follows that M right (·) and M right,X (·) are well defined in a neighbourhood of λ γ .

Remark 1
When M left (λ) has a pole at λ, it is still possible for λ to be an eigenvalue of L 0 + iγ S. However, M right (λ) and M right,X (λ) cannot have poles off the real axis, as Lemma 4 and Lemma 5 show that (R, λ) and X (R, λ) are invertible for (λ) = 0.

Theorem 2 Suppose that assumptions (A1) and (A2) hold. For γ > 0, let λ γ,X,bad
be an eigenvalue of the non-self-adjoint Schrödinger operator L 0,X + iγ S defined in (5,6) which converges, as X → +∞, to a point which is not in the spectrum of L 0 + iγ S. Then for some positive constants C 5 and C 6 : Proof Similarly to the proof of Theorem 1, it is sufficient to consider the case X = R + Na where N ∈ N. Since λ γ,X,bad has the property: ker M left (λ γ,X,bad ) + M right,X (λ γ,X,bad ) = {0}, in particular M left (λ γ,X,bad ) + M right,X (λ γ,X,bad ) has an eigenvalue 0. It is known (see, e.g., [3]) that spectral pollution must lie on the real axis, since L 0 + iγ S is a relatively compact perturbation of the semi-bounded self-adjoint operator L 0 : thus, λ γ,X,bad → μ ∈ R ∩ {Spectral Gap} as X → ∞. Additionally, since μ is not in the spectrum of L 0 + iγ S, the matrix M left (μ) + M right (μ) is invertible if it is well defined, i.e. if (R, μ) and U(R, μ) are invertible (see (14,23)). We assume without loss of generality that this is indeed so, for if not then one may simply use a greater value of R for the Glazman decomposition. Hence, there is a compact neighbourhood of μ, say B(μ, r); r > 0 such that M left (λ) + M right (λ) is invertible ∀λ ∈ B(μ, r). Following the reasoning which led to (42), we deduce that provided C R (μ) is well defined, then M right,X will converge locally uniformly to M right in a neighbourhood of μ. Since such a uniform convergence excludes spectral pollution, it follows that C R (μ) must not be well defined, i.e. cos(β) (R, μ) − sin(β) (R, μ)is not invertible.

Theorem 3 Suppose that assumption (A2) holds (see
holds for some positive constants C and C 2 , then there exists C 1 > 0, independent of R, such that for all R > 0, where c ∈ (0, 1) is the constant appearing in assumption (A2).
Proof In fact, we prove more than stated: only the estimate (48) depends on (A2); the rest of the theorem requires only that S to be bounded independently of R as a multiplication operator. The existence of λ γ,j with |λ γ,j − λ| → 0 as γ → 0 is a consequence of results in [11] on analytic families. Since γ is sufficiently small, let be a contour which encloses the spectral point λ of L 0 and suppose that λ γ,j , j = 1, . . . , ν are the only spectral points of L 0 + iγ S inside . Clearly, S = 1 independently of R and since L 0 is a self-adjoint operator then |λ − λ γ,j | ≤ γ independently of R; thus, the contour can be chosen independently of R. Suppose u γ,j , j = 1, . . . , ν are eigenvectors of L 0 +iγ S, linearly independent with u γ,j = 1. Following Kato [11,VII,§3], let P (γ ) be the projection onto the eigenspace of L 0 + iγ S spanned by the u γ,j , j = 1, . . . , ν, and P (0) be the projection onto the eigenspace of L 0 associated with λ; the projection P (γ ) is analytic as a function of γ, so that Since taking the norm of (49) and using (48), we conclude that P (0)u γ,j is bounded away from zero uniformly with respect to R and γ for sufficiently small γ . Now, since (L 0 + iγ S)u γ,j = λ γ,j u γ,j and (L 0 − iγ I )u k = (λ − iγ )u k , using the inner product for the first equation with u k , we obtain: Similarly, using the inner product for the second equation with u γ,j , we obtain: Because L 0 and S are self-adjoint and u k and u γ,j are in the domain of L 0 and it is contained in the domain of S then from (51), we have: From (50) and (52), we obtain: |λ + iγ − λ γ,j | u γ,j , u k = iγ (I − S)u γ,j , u k = iγ u γ,j , (I − S)u k .
Since P (0)u γ,j is bounded away from zero, we may choose k (possibly depending on γ ) such that u γ,j , u k is bounded away from zero for small γ , uniformly with respect to R; furthermore, from the assumption (46) and (A2), we deduce for some positive constants C and C 2 . The result is proved.

Remark 2
We may also ask what happens when an eigenvalue λ γ merges into the essential spectrum with decreasing γ . For the scalar problem, Theorem 5 in [13] states that in this situation, there is a strictly positive value γ crit > 0 of the the coupling constant at which λ γ merges into the essential spectrum. In our current setup, this is false. Consider the following system: and S(x) = I 2 in [0, 50], 0 in (50, ∞). The system can be seen as two scalar problems with q 1 (x) = 0 and q 2 (x) = −40 1+x 2 χ [0,40] (x) + sin(x). The problem with potential q 2 and γ = 0 has a spectral gap which is approximately (−0.340363, 0.595942), with infinitely many eigenvalues accumulating from below at the top end of the gap [15]. However, these eigenvalues all lie in the essential spectrum [0, ∞) for the problem with potential q 1 and hence also lie in the essential spectrum of the matrix Schrödinger system. Nevertheless, they emerge from the real axis with positive speed, into the upper half-plane, as soon as γ is increased from zero, following the estimate in Theorem 3, (47).

Numerical examples
In this section, we present some numerical examples to demonstrate our theoretical results. We generally apply the 3-point finite difference scheme to our system; the number of steps-per-period is fixed when the number of periods increases. The main advantage of using 'low tech' finite differences rather than a spectral method or Galerkin method is that there exists a Floquet theory for periodic Jacobi matrices, so the theorems of the previous section may be expected to have analogues for the (infinite) discretised system. Example 1 Consider the following matrix Schrödinger systems: problem P satisfies hypothesis (A1) while problem P' satisfies (A1'). The fact that P' satisfies (A1') can be proved by an ODE version of the Agmon-type results in Janas, J. et al [10]. In both these equations, we use: Both problems have exactly the same essential spectrum, though their point spectra are slightly different due to the compactly supported potential well in P. The first two spectral bands are approximately [13]: For problem P', we expect an eigenvalue close to λ ≈ −0.1076 in a gap, and an eigenvalue embedded in a spectral band near λ ≈ 0.6336. For both problems, we use the perturbation: and perceive the resulting eigenvalues with γ = 1/4. Figure 1 shows eigenvalue computations for both P and P'. They show that our estimate (47) holds with C 1 ≈ 1293.6 and C 2 ≈ 0.5386 for (P'), and C 1 ≈ 1291.9 and C 2 ≈ 0.5384 for (P). These figures were calculated using X = 100 and a step-size h = 0.1, both of which were chosen to ensure that the effects of discretisation error would not be seen in the estimated constants. For the embedded eigenvalues, Fig. 2 shows the effects of interval truncation. In particular we observe that for P the prediction of Theorem 1 holds with C 3 ≈ 0.0039 and C 4 ≈ 0.3842. Theorem 1 has not been proved for P', which does not have an eventually periodic Q. However, the experiments indicate that the result still holds, with C 3 ≈ 0.0047 and C 4 ≈ 0.3842. For these experiments, the step-size was h = 0.1.
Finally, we shall show that the prediction of Theorem 2 for a spurious eigenvalue also holds; in fact, Fig. 3 indicates that it holds not only for P, but also for P', for which it is not proved. In these figures, we consider the value λ bad ≈ −0.1847, which lies in a spectral gap but is not an eigenvalue of either problem. We fixed X = 55 and varied R from 19 to 43 in steps of 4. The value of h was again chosen small enough to suppress effects of discretisation error. The constants C 5 and C 6 predicted by Theorem 2 are C 5 ≈ 0.0017 and C 6 ≈ 0.5345; it seems also that C 5 ≈ 0.0017 and C 6 ≈ 0.5345 for P'.
with the following perturbed periodic potential: This rather unwieldy formula was obtained by an expression in whichQ(x) = −6k(x)I 5 + diag(n 2 /4, n = 1, . . . , 5) and T is the matrix of orthonormal eigenvectors of a randomly chosen real, symmetric 5×5 matrix. According to [2], one eigenvalue for the scalar problem with q(x) = −6k(x) is −1; from this, we know that one of the eigenvalues of our Hamiltonian system is 1.25 which is an embedded eigenvalue in the essential spectrum [1/4, ∞) of this multi-channel problem. We expect that if the dissipative barrier method is working well, then we should see an eigenvalue close to 1.25+0.25i when γ = 1/4. Figure 4 shows the plot of the corresponding eigenfunction, for which the calculated eigenvalue was approximately 1.24998 + 0.24993i. This illustrates the usefulness of the dissipative barrier method for lifting embedded eigenvalues out of the essential spectrum and making them easier to calculate, even bearing in mind Remark 2.

Example 3 We consider an equation of the form
Here W is a continuous, periodic matrix, strictly positive definite, so that the corresponding weighted L 2 space is equivalent to L 2 (0, ∞). The matrix S(x) was chosen to be compactly supported and bounded. The matrix D is a diagonal matrix with D(j, j ) = j 2 , j = 1, . . . , n; qualitatively this is a second-order differentiation matrix. The main differences compared to the cases considered in our theorems are, firstly, the weight W (x) is no longer the identity; and secondly, the compactly supported dissipative barrier S(x) is now also multiplied by the spectral parameter. The form of this problem is inspired by work in optics with complex index of refraction (see, e.g., [7]).
It is not difficult to show that the Floquet theory still holds for the problem on [R, ∞), and so M right,X (λ) still converges exponentially fast to M right (λ) as X → ∞. For the problem on [0, R], there is now an additional λ-dependence of the barrier term, but M left (λ) is still meromorphic. We therefore expect to see fast convergence of the eigenvalues lying well away from the real axis. Figure 5 shows the results of computations in the purely diagonal case with all matrices being of dimension 5 × 5. These results were computed using the Numerov discretisation [14], with a uniform mesh of 80 intervals per period (mesh size 2π/80). Because the values of λ in Fig. 5 are not large, this mesh size is sufficient to ensure that the points plotted in Fig. 5 will not move in the 'eyeball norm' if the mesh size is halved. In Fig. 5, we see that the asterisks (shorter interval approximations) and circles (longer interval approximations) are essentially coincident for the eigenvalues well away from the real axis. We expected this, due to the exponentially small error which interval truncation causes. The more interesting parts are the 'loops' in the upper half plane, one of which starts at approximately λ = 1.75 and returns to the real axis around λ = 2.6. Here the asterisk loop (approximating interval [0, 20π ]) is approximately twice as far from the real axis as the circle loop (approximating interval [0, 40π ]). These loops are approximations to a spectral band. Though we have not proved this, they appear to converge at a rate 1/X, as the interval [0, X] goes to infinity. Note that there are also other approximations to (parts of) spectral bands, due to the fact that the spectral multiplicity of the higher bands can be greater than 1. It seems that, in this picture, only the bands near 0.5, and from 0.75 to just below 1, are simple. In Fig. 6, we repeated the experiments using non-diagonal W (x) given by The same phenomena are noted as in the diagonal coefficient case, though the different scale on the vertical axis makes the slow convergence to the essential spectrum more stark.

Conclusion and future work
We have introduced a method to calculate eigenvalues in gaps of matrix-valued Schrödinger operators. Theoretically, we have shown that the relatively compact dissipative perturbation technique together with domain truncation obtain approximations of isolated eigenvalues close to the ones of the original problem. Moreover, spurious eigenvalues can be predicted using this method and are characterised by exponentially small imaginary parts. These approximations have been computed using finite difference schemes for some numerical examples and have shown excellent agreement with the theoretical part. An additional remark on this procedure is that the approximating results of the implementations of both fast decaying periodic potentials and compactly supported periodic potentials (see, e.g., Example 1) are mostly the same.
We have also observed the effectiveness of the presented method when the weighted matrix is different from the identity and the dissipative barrier is multiplied by the spectral parameter as in Example 3. The only caveat for these cases is that the approximations to the spectral bands do not converge fast, so if very high accuracy is required then one should use λ-dependent non-reflecting boundary conditions [7].
One of the main sources of Hamiltonian eigenvalue problems is the use of semidiscretisation for PDE problems on waveguides. In future work, we will consider the PDE: − u + Q(x, y)u = λu, on a semi-infinite waveguide. These give rise to a Hamiltonian system in which Q(x) in our problem is an operator in l 2 (N); they can also be studied directly in PDE form.