# Wave packet pseudomodes of variable coefficient differential operators

Lloyd N Trefethen

## Abstract

The pseudospectra of non-selfadjoint linear ordinary differential operators with variable coefficients are considered. It is shown that when a certain winding number or twist condition is satisfied, closely related to Hörmander's commutator condition for partial differential equations, ϵ-pseudoeigenfunctions of such operators for exponentially small values of ϵ exist in the form of localized wave packets. In contrast to related results of Davies and of Dencker, Sjöstrand & Zworski, the symbol need not be smooth. Applications in fluid mechanics, non-hermitian quantum mechanics and other areas are presented with the aid of high-accuracy numerical computations.

## 1. Introduction

Certain differential and pseudodifferential operators with variable coefficients have exponentially large resolvent norms, with respect to a parameter, in regions of the complex plane far from the spectrum. These large norms are explained by the existence of pseudoeigenfunctions in the form of localized wave packets, which, although they may not satisfy the eigenvalue equation or the boundary conditions exactly, satisfy them approximately with an exponentially small error. Another way to describe this phenomenon is to say that these operators have ϵ-pseudospectra that extend over regions of the complex plane far from the spectrum, even for exponentially small values of ϵ. Applications in which such effects may be important include ‘ghost solutions’ of ordinary differential equations (Domokos & Holmes 2003), non-Hermitian quantum mechanics (Bender & Boettcher 1998; Bender et al. 1999; Davies 1999a,b; Benilov et al. 2003), the theory of dichotomy for ordinary differential equations and their numerical approximations (Massera & Schäfer 1966; Daleckii & Krein 1974; Coppel 1978; Ascher et al. 1995; Chicone & Latushkin 1999), fluid dynamics (Reddy et al. 1993; Trefethen et al. 1993; Cossu & Chomaz 1997; Chapman 2002; Benilov et al. 2003) and the Lewy/Hörmander phenomenon of non-existence of solutions to certain linear partial differential and pseudodifferential equations (Lewy 1957; Hörmander 1960; Zworski 2003).

The fact that variable coefficient non-selfadjoint differential operators may have extended pseudospectra with wave packet pseudomodes was pointed out by Davies (1999a,b). Shortly thereafter, Zworski observed that Davies's discoveries could be related to long-established solvability results in the theory of PDE by Hörmander and others (Hörmander 1960; Zworski 2001). Later publications stimulated by Davies's work include Aslanyan & Davies (2000), Boulton (2002), Pravda-Starov (2003), Zworski (2003), Dencker et al. (2004) and Davies & Kuijlaars (2004). The most convenient framework for all of these results is a semi-classical formulation, in which differentiated terms contain a small parameter h or its powers (Dimassi & Sjöstrand 1999; Martinez 2002). Depending on the problem and the method of proof, the meaning of ‘exponentially small’ for such an operator may involve terms of order O(M−1/h) for some M>1 or of order O(hN) for all N>0 as h→0.

The standard method of deriving results about wave packet pseudomodes is WKBJ or microlocal analysis, in which a wave packet is constructed that is localized with respect to both the space variables and their dual, the wavenumber vector. This approach applies to both differential and pseudodifferential operators in any number of space dimensions, provided that the coefficients are smooth; the typical assumption is C. The most important presentation of this approach to date is the recent article of Dencker et al. (2004).

For the special case of ordinary differential operators, however, an alternative method of proof is available, based on consideration of intersections of certain subspaces, which requires no assumption of smooth dependence on the space variable. The purpose of this article is to apply this method to establish a fundamental theorem about the existence of exponentially good pseudomodes for such problems. Instead of smoothness, this theorem depends upon a condition involving winding numbers of the symbol. One might think that this replacement of one condition by another would simply be a technical matter, of little significance for applications, since the conclusions reached about exponentially good pseudomodes are in the end the same. However, we shall show that the difference is a substantial one (§7). For at least some problems with smooth coefficients that violate the winding number condition, the pseudospectral effects are structurally unstable and vanish when the coefficients are perturbed in a non-smooth manner. By contrast, the effects associated with operators that satisfy the winding number condition are robust.

This article also has another purpose—to relate these phenomena to analogous effects that have been investigated in the highly developed field of Toeplitz matrices and operators and their generalizations with variable or ‘twisted’ coefficients. The consideration of winding numbers is standard in the Toeplitz literature, as is the focus of attention on the case of one space dimension (although there are generalizations to block Toeplitz matrices and operators and to Berezin–Toeplitz operators on manifolds; Borthwick & Uribe 2003). The analogy can be portrayed at a high level by considering four fundamental classes of non-normal operators:

View this table:

In all four cases, the resolvent norm grows exponentially for certain problems, as an appropriate parameter is increased or decreased, throughout a region of the complex plane determined by the symbol curves of the operator and their winding numbers. For each, we have listed the two most important published references. The present contribution belongs to the upper-right cell of this 2×2 array. Mathematically, it is an analogue for differential operators of the earlier work for twisted Toeplitz matrices by Trefethen & Chapman (2004). Both our formulation of the main theorem and its proof stay as close as possible to those in Trefethen & Chapman (2004), even to the point of repeating wordings where possible. However, the reader who compares the two papers will find that the two developments differ in many ways.

## 2. An elementary example

We can see the essence of the matter in a simple example. Consider the operator(2.1)where h is a small parameter, acting on a dense subspace of L2[−1,1] of differentiable functions satisfying the boundary conditions u(−1)=u(1)=0. The equation Ahu=λu cannot be solved for any λC; the spectrum of Ah is empty. However, the function(2.2)satisfies Ahu=0, and it satisfies the boundary conditions up to an error e−1/2h. Alternatively, we could note that satisfies the boundary conditions and satisfies Ahu=0 up to an error no greater than e−1/2h.

Thus, we can say that the Gaussian (2.2) of width O(h1/2) (or more precisely, the same function shifted by a constant) is an ϵ-pseudoeigenfunction corresponding to the ϵ-pseudoeigenvalue 0 for a value of ϵ of size O(M−1/h) for some M>1 as h→0.1 Moreover, the same is true for any number λ with −1<Re λ<1, as is shown by the pseudoeigenfunction(2.3)The situation is summarized in figure 1 for h=0.02. This and the other figures in this article were computed numerically by spectral discretization (Trefethen 2000) and the EigTool pseudospectra plotting system (Wright 2002); we shall not give details. We see that the pseudospectra of Ah approximate the strip −1<Re λ<1, and for the particular value λ=1/2+i, the optimal pseudoeigenfunction comes very close to the predicted form: a wave packet centred at x=1/2 with wavenumber 1/h=50 (i.e. wavelength 2π/50≈0.13).

Figure 1

Above, ϵ-pseudospectra of the operator (2.1) with h=0.02 (the spectrum is empty). From outside in, the contours correspond to ϵ=10−1,…,10−9. The figure confirms that the resolvent norm ‖(λAh)−1‖ grows exponentially as h→0 for any λ in the strip −1<Re λ<1, marked by the dashed lines. Below, an optimal pseudoeigenfunction for λ=1/2+i (marked by a cross in the top image) with x*=1/2 and k*=−1. (Both real part and envelope are shown.)

Once we have proved our main theorem and expressed it in terms of winding numbers, this example will be interpretable as follows. We associate the operator Ah of (2.1) with the x-dependent symbol(2.4)which maps the real k-axis onto the negatively oriented vertical line Re f=x in the complex plane. We define the winding number of this symbol curve about any point λC by completing it by a large semicircle traversed counter-clockwise in the right half-plane. Thus for each x, the winding number is 0 if Re λ<x, 1 if Re λ>x, and undefined if Re λ=x. As x increases from −1 to 1, each λ with −1<Re λ<1 experiences a decrease in winding number when x passes through the value Re λ. Theorem 3.1 ensures that each such value of λ is an ϵ-pseudoeigenvalue of Ah for a value of ϵ of size O(M−1/h) for some M>1 as h→0.

## 3. Theorems

Let an interval [a,b] be given, a<b, and for a small parameter h>0, let Dh be the scaled derivative operatorFor some integer n≥0, let continuous coefficient functionsbe defined on (a,b), which may or may not be smooth. We consider the family of linear operators {Ah}, h>0, defined bytogether with arbitrary homogeneous boundary conditions at x=a and x=b (the details of the boundary conditions will not matter), acting in a suitable dense domain in L2[a,b].

Given kC and h, consider the function in L2[a,b]. We can write Ahv explicitly asIn other words, we havewhere f is the symbol of {Ah}, defined for x∈(a,b) and kC by(3.1)Our aim is to build wave packets localized near a particular x*∈(a,b). Let x* and k* be given, and define λ=f(x*,k*). If f were independent of x, then the function(3.2)would satisfy the eigenfunction equation for Ah with eigenvalue λ,(3.3)although, in general, it would not be an eigenfunction because of the boundary conditions. If f varies with x, however, then we shall look for solutions to (3.3) near x* with the form of wave packets (exact solutions, not approximate), and if these decay exponentially away from x*, then they can be extended smoothly to zero so as to make exponentially good pseudoeigenfunctions, regardless of the boundary conditions.

The symbol f=f(x,k) satisfies the twist condition at x=x*∈(a,b), k=k*R if at this point it is differentiable with respect to x withf/∂k≠0 and(3.4)It satisfies the antitwist condition if it has the same properties with (3.4) replaced by(3.5)Here is our main theorem.

Let {Ah} be a family of variable coefficient differential operators on [a,b] with symbol f(x,k) as described above. Given x*∈(a,b) and k*R, define λ=f(x*,k*), and suppose that the twist condition is satisfied. Moreover, suppose that f(x*,k)≠λ for all real kk* and that cn(x*)≠0. Then, there exist constants C1, C2>0 and M>1 such that for all sufficiently small h, there exists a non-zero pseudoeigenfunction v(h) that is exponentially good,(3.6)and localized,(3.7)

A similar result applies if the antitwist condition is satisfied. In this case the pseudoeigenfunctions are localized not at x* but at one or more other points in the interior of [a,b] or at the boundaries.

Suppose the same assumptions hold as in theorem 3.1 except that the antitwist condition (3.5) is satisfied instead of the twist condition. Assume also that λ is not in the spectrum of Ah. Then, there exists M>1 such that for all sufficiently small h, there exists a non-zero pseudoeigenfunction v(h) satisfying (3.6).

The remainder of this article is devoted to proving and illustrating theorem 3.1. We shall not spell out the proof of theorem 3.2 in full, which is analogous to that of theorem 7.3 of Trefethen & Chapman (2004), but in outline, it follows from theorem 3.1 by considering the operator Bh defined as the complex conjugate of the adjoint of Ah; that is,together with boundary conditions adjoint to those of Ah. If f(x,k) is the symbol of {Ah}, then the symbol of {Bh} has the form f(x,−k)+O(h), an O(h) perturbation of a function that satisfies the twist condition. It can be shown that the arguments used to prove theorem 3.1 remain valid in the presence of this perturbation, implying that for all sufficiently small h. If ‖(λAh)−1‖ is finite, then this implies (3.6) as required. If ‖(λAh)−1‖ is infinite—that is, λ is in the spectrum of Ah—then we know that λAh lacks a bounded inverse, but it does not necessarily have eigenfunctions or good pseudoeigenfunctions. For example, suppose we negate the sign in (2.1) and consider the operator(3.8)acting on the space of absolutely continuous functions in L2[−1,1] with u(−1)=u(1)=0. Then, Bh is the operator (2.1) but with no boundary conditions. Its spectrum is all of C, so the same is true of Ah, but whereas Bh has an eigenfunction (and hence good pseudoeigenfunctions) for each λC, Ah has ‘too many boundary conditions’ and no eigenfunctions.

## 4. Winding number interpretation

The condition f(x*,k)≠λ for kk* of theorems 3.1 and 3.2 might easily be overlooked, but it is important. It is the price we pay for not requiring f to depend smoothly on x. This condition has a geometric interpretation. Given a symbol f=f(x,k) and a particular value x=x*, define the symbol curve of f at x* to be the curve f(x*,R), a subset of the complex plane. Given a number , we wish to define the winding number of f with respect to λ. Because f(x*,R) is not a closed curve, this requires some care. A suitable approach is to make use of contours ΓR consisting of the real interval [−R,R] closed by a semicircle of radius R in the upper half complex plane. We consider the winding number of f about λ associated with this closed contour traversed in the usual anticlockwise direction, and we define the winding number to be the limiting winding number of f(x*,ΓR) about λ obtained for all sufficiently large R. If λf(x*,R), is undefined.

It is readily verified that the conditions of theorem 3.1 and 3.2 have the following geometric interpretations:

1. the symbol curve at x* passes just once through λ; and

2. as x increases through x*, the winding number about λ decreases by 1 (theorem 3.1) or increases by 1 (theorem 3.2).

One can use these conditions to assess quickly where the pseudospectra of a variable coefficient differential operator will lie (figure 2). One imagines the symbol curve moving about the complex plane as x increases from a to b. The pseudospectra fill the region swept out by this moving curve.

Figure 2

Winding number interpretation of theorem 3.1. If the curve f(x,R) crosses λ as x increases through x* in such a way that the winding number about λ is decreased, then there is an exponentially good wave packet pseudomode localized at x* with pseudoeigenvalue λ.

In fact, we can make an explicit connection between the winding number at a point x where it is defined and the eigenvalues of the associated transfer matrix C(x) introduced below in (5.1) and (5.2). By the principle of the argument of complex analysis, the symbol polynomial f at x has exactly I zeros in the upper half-plane. In lemmas 5.1 and 5.2, this will correspond to a space of amplified functions of dimension I and a space of attenuated functions of dimension nI.

## 5. Proof of theorem 3.1

In this section, we prove theorem 3.1. The essential idea, involving a non-empty intersection of subspaces of functions decaying to the left and the right that satisfy the eigenvalue equation exactly in an interval, is summarized at the end. The assumption of differentiability in the theorem is used in the estimates of lemma 5.2 associated with the constant γ.

Given an n-times differentiable function v(x), defineandThen, (3.3) can be written as(5.1)or, more compactly, as(5.2)where for 1≤jn and . Thus, we have a first-order vector-valued linear ordinary differential equation (ODE) with continuous coefficients defined throughout the neighbourhood of x* where cn(x)≠0, with n linearly independent solutions throughout this region (e.g. sect. 3.6 of Coddington & Levinson (1955)). We may associate equation (5.1) frozen at the point x* with the characteristic polynomial(5.3)For κC, the frozen-coefficient ODE has a solution with(5.4)if and only if κ satisfies p(κ)=0. The polynomial p has exactly n roots {κ} counted with multiplicity for x sufficiently close to x*. Each root with Imκ>0 corresponds to exponential growth as x increases and exponential decay as x decreases, and for Imκ<0, the pattern is reversed. As mentioned at the end of the §4, if the winding number of f about λ at x* is well-defined, then there will be I roots of the former type and nI roots of the latter type.

Of course, the problem at hand involves variable coefficients. Our choice of λ will ensure that one root has Imκ=0 and passes from above to below the real axis as x increases through x*, and thus corresponds to exponential decay in both directions. We begin in lemma 5.1, however, with a result based on the assumption that the roots are separated from the axis. We work with a fixed range of values of x around x*, chosen narrow enough so that the coefficients {c(x)} vary sufficiently little and so that cn(x)≠0. We consider the variable coefficient ODE(5.5)with C=C(x*), where E(x) will be taken to be small and to depend continuously on x, and let Sh(x) denote its solution matrix or matrizant: if v(x*) is an n-vector, then is the corresponding unique solution of (5.5) near x*. (If E(x)≡0, .) We first show that under suitable hypotheses, the solutions to this equation attenuate and amplify certain vectors exponentially. This is a basic result with the flavour of the stable manifold theorem in dynamical systems. Afterwards, we shall refine it to the more specialized lemma 5.2, analogous to the centre manifold theorem, which is actually the one needed to prove theorem 3.1 (Shub 1987; Wiggins 1990).

Let C be an (η+ν)×(η+ν) matrix (η, ν≥0,η+ν≥1) that has η eigenvalues κ with Imκ<ρ<0 and ν eigenvalues with Imκ>R>0, let E(x) be a continuous matrix function withE(x)‖≤ϵ, let Sh(x) be the solution matrix for (5.5) for x near x*=0, and let Δx>0 be fixed. There exist ϵ>0, independent of h, such that the solution matrices {Shx)} separate Cη+ν into exponentially amplified and attenuated subspaces in the following sense: for each sufficiently small h>0, there is an η-dimensional subspace such that(5.6)and a ν-dimensional subspace such that(5.7)

We shall see that ν can be chosen independently of h, since almost all vectors lead to exponential growth, but η must depend on h, since exponential decay is more delicate.

The first step is to reduce C to block-diagonal form. There exists a non-singular matrix X such that(5.8)where Gν has dimension ν and amplifies all ν-vectors,(5.9)and Gη has dimension η and shrinks all η-vectors,(5.10)we may choose any T>R which is smaller than the imaginary parts of all the eigenvalues of C in the upper half-plane, and any τ<ρ which is larger than the imaginary parts of all the eigenvalues of C in the lower half-plane. If (5.5) is left-multiplied by X, then we obtain a new differential equation for the variable Xv involving the matrix (5.8); the norms of X and X−1 introduce fixed constants and thus have no effect on the small-h assertion. Because E(x) becomes XE(x)X−1 in this process, the norm bound ϵ has to be adjusted by the condition number . All this is straightforward, and rather than encumber the rest of the argument with X and associated details, let us assume from now on, without loss of generality, that C itself has block diagonal form to begin with,(5.11)with Gν and Gη satisfying (5.9) and (5.10).

We identify first a growing space Sν as in (5.7). Define δ, ϵ>0 byand let denote the cone of all vectors satisfying

Consider what happens when Sh(x) acts on a vector :assuming (without loss of generality) that ‖v‖=1. At x=0, we have, by (5.5) and (5.11),from which we calculate, using (5.9),and, using (5.10),Thus, locally Sh(x) increases the v component of a vector at a rate greater than R and decreases the w component, if the latter has norm greater than δ/2. These observations imply that Sh(x) maps into itself and that (5.7) holds for sufficiently small h for all u. Finally, we note that contains the ν-dimensional invariant subspace associated with Gν; that is, the subspace spanned by the ν vectors . Thus, this subspace is a suitable choice for Sν.

To find a decaying space Sη for (5.6), we consider the inverse solution operator Shx)−1, withBy the sine argument as before, Shx)−1 has a growing space of dimension η, and we take . ▪

The next lemma is a variant of lemma 5.1 designed for systems that have a root passing through real axis, as occurs at the centre of a wave packet pseudomode.

Let {C(x)} be a family of (η+ν+1)×(η+ν+1) matrices (η,ν≥0) depending differentiably on x at x=x*. Assume that C*=C(x*) has η eigenvalues with Imκ<ρ<0, ν eigenvalues with Imκ>R>0, and one eigenvalue μ* with Imμ*=0. Assume further that dImμ/dx<0 at x=x*, where μ=μ(x) is the eigenvalue of C(x) that converges to μ* as xx*. For x near x* and h>0, consider the solution matrix Sh(x) as before. There exist Δx>0 and M>0, independent of h, such that these operators separate Cη+ν+1 into exponentially amplified and attenuated subspaces in the following sense: for each sufficiently small h, there is an (η+1)-dimensional subspace such that(5.12)and a ν-dimensional subspace such that(5.13)

It is only (5.12) we shall need, not (5.13), so this is the claim we shall prove. As in the final paragraph of the proof of lemma 5.1, to show that Shx) has a decaying subspace of dimension η+1, we shall show that Shx)−1 has a growing subspace of this dimension. Our first step, as in the proof of lemma 5.1, is to assume without loss of generality that −C(x*) has the block-diagonal formwhere Gη has dimension η and amplifies all η-vectors as in (5.9), Gν has dimension ν and shrinks all ν-vectors as in (5.10) and Imμ*=0, with τ<0<T.

Without loss of generality we assume x*=0, as in lemma 5.1.

The proof rests on an observation of linear algebra: if a diagonal matrix is perturbed by O(ϵ) in the off-diagonal positions, its eigenvalues of multiplicity 1 change by only O(ϵ2). (One can prove this by considering the characteristic polynomial.) We exploit this phenomenon as follows. Our task is to consider the solution operator associated with a set of matriceswith ‖E(x)‖=O(x), where E(x) is constructed to have a zero entry in the middle position (η+1,η+1) and c(x) is defined accordingly. Because , it follows from the fact of linear algebra just mentioned that . Choose γ to be a constant in the range 0<γ<−dImμ/dx, and collect the upper-left (η+1)×(η+1) block of C(x) into a single matrix Gη+l(x). Then, for some σ>0 and xmax>0, we have(5.14)and(5.15)(5.16)for some and all xxmax. (The reasoning that establishes (5.15) is as follows: if G is a matrix and α is the minimal imaginary part of the field of values of G, i.e. the minimal eigenvalue of i(GG*)/2, then for any v; here, by the O(ϵ2) observation above, α increases faster than γx as x increases from 0.) We now follow estimates as in the proof of lemma 5.1. DefineLet be the cone of all vectors satisfyingConsider how a vector evolves locally near a point x∈[0,Δx], assuming, without loss of generality, that ‖v‖=1. We havefrom which we calculateandThese estimates show that Sh(x)−1 maps into itself and increases the norm of the upper part of the vector at a rate at least γx/2h. It follows that if u, then and for all sufficiently small h, for some constant MγΔx/4. Finally, contains the (η+1)-dimensional invariant subspace associated with Gη+1, which is accordingly a suitable choice of the subspace needed. ▪

To establish theorem 3.1, we construct a wave packet in an interval [x*−Δx,x*x] that satisfies the eigenvalue equation (3.3) exactly in that interval and is also localized in the sense of (3.7) there. We do this by applying lemma 5.2 in both directions, moving both right and left from x=x*, taking advantage of the equivalence between (3.3) and (5.1).

First, let us move from x* to the right. The hypothesis of theorem 3.1 that f(x*,k)≠λ for all real kk* ensures that the matrix C(x*) has just a single eigenvalue μ* on the real axis, as required in lemma 5.2. The twist condition ensures that dImμ/dx<0 at x=x*, and the hypothesis cn(x*)≠0 ensures that C(x) is well defined for all x sufficiently close to x*. The matrices Sh(x) of lemma 5.2 are then the solution operators that transfer an n-vector of data to the right a distance x. According to lemma 5.2, for some η≥0 there is an (η+1)-dimensional subspace of such data vectors that generate solutions in [x*,x*x] that are exponentially small in the sense of (5.12) at the end of that interval. By the same argument, although this is not spelled out in the statement of lemma 5.2, these solutions have the Gaussian decay behaviour described by (3.7).

Similarly, by an obvious symmetry, we can apply lemma 5.2 moving from x* to the left. We conclude that with ν=m+n−1−η there is a (ν+1)-dimensional subspace of data vectors that generate solutions in [x*−Δx,x*] with the appropriate Gaussian decay.

Now, we take the intersection of these two subspaces. Their dimensions are η+1 and ν+1=m+nη, and they are subspaces of Cm+n. It follows that is a subspace of dimension at least 1 (in fact, it will be exactly 1) of vectors that generate solutions of the eigenvalue equation in [x*−Δx,x*x] that decay appropriately on both sides of x*.

Thus there exists a vector that exactly satisfies the eigenvalue equation in [x*−Δx,x*x] and has the necessary decay there. Near x*−Δx and x*x, such a vector will be exponentially small. By extending it smoothly to 0 outside this interval, we obtain a pseudoeigenfunction v(h) satisfying both (3.6) and (3.7). ▪

## 6. Examples

We shall present six numerical examples, the first three fitting the framework of theorem 3.1 and the others having a double-crossing symbol curve, so that only the theorems of Dencker et al. (2004) can be applied. The next section will comment on the significance of this distinction.

Our first example is a variable coefficient advection-diffusion operator from a paper of Cossu & Chomaz (1997). The operator (after some simplification) is(6.1)for xR, and the eigenvalues and eigenfunctions are known explicitly:(6.2)for n=1,2,3,…, where Hn is the nth Hermite polynomial (Chomaz et al. 1987). Pseudospectra for the case h=0.02 are shown in figure 3, and we see that the resolvent norm is exponentially large in the region of the complex plane bounded by the parabola Re λ=1/4−(Imλ)2. To explain this behaviour, we note that the symbol is(6.3)which implies that the symbol curves are parabolas adjusted by the variable offset 1/4−x2. The winding number about a value λ decreases from 1 to 0 as this curve crosses λ from left to right, which occurs for a negative value of x, and this explains why the pseudoeigenfunction in the figure sits in the left half of the domain.

Figure 3

Above, eigenvalues and ϵ-pseudospectra of the operator (6.4) of Cossu & Chomaz (1997) with h=0.02, ϵ=10−2,…,10−15. By theorem 3.1, the resolvent norm ‖(λA)−1‖ grows exponentially as h→0 for any λ lying in the region to the left of the dashed parabola. Below, an optimal pseudoeigenfunction for λ=0.05i (marked by a cross in the top image) with central position and wavenumber x*=−5−1/2≈−0.45 and k*=0.

This example illustrates some behaviour of widespread physical importance. The operator (6.1) can be interpreted as a standard advection–diffusion process coupled with a factor that introduces exponential amplification for |x|<1/2 and exponential attenuation for |x|>1/2. It is clear that the associated time-dependent process ut=Ahu must be susceptible to transient growth of order O(C1/h) on a time-scale O(h−1) for some C>1, for a pulse will grow exponentially during the time of order O(h−1) that it spends passing through the amplification region. Cossu and Chomaz relate this transient effect to the local convective instability of fluid flows in unbounded domains, including wakes, jets, and boundary layers. Alternatively, it is implied by the appearance of exponentially large resolvent norms in the right half-plane (see theorem 5 of Trefethen 1997).

Our second example is a more exotic advection–diffusion operator investigated by Benilov et al. (2003). These authors consider the instability of a thin viscous liquid film on the inner surface of a rotating cylinder in an approximation in which gravitational effects are included but inertial and capillary effects are ignored. They reduce their problem (again after some simplification) to the operator(6.4)with periodic boundary conditions on [−π,π], with symbol . An unusual feature here is that for two values of x, the coefficient of the quadratic term passes through zero. For each x, the symbol curve is the parabola described in the direction of decreasing imaginary part. Completing this curve in the usual manner by a semicircle at infinity, we see that its winding number is 1 about points to its right and 0 about points to its left. For any λ in the domain bounded by the two parabolas Re z=±(Imz)2, the winding number accordingly decreases by 1 at some value of x in the interval (−π,−π/2) or (π/2,π), leading to an exponentially good wave packet pseudoeigenfunction. Figure 4 confirms this geometry for the case h=0.1.

Figure 4

Above, eigenvalues and ϵ-pseudospectra of the operator (6.4) of Benilov et al. (2003) with h=0.1 and ϵ=10−1,…,10−7. By theorem 3.1, ‖(λAh)−1‖ grows exponentially as h→0 for all λ in the region between the two dashed parabolas. Below, an optimal pseudoeigenfunction for λ=1+1.5i (marked by a cross) with central position and wavenumber and k*=−1.5.

From a dynamical point of view, a distinctive feature of (6.4) is that although the pseudospectra fill unbounded expanses of the right half-plane, all the eigenvalues lie on the neutrally stable imaginary axis. Benilov et al. (2003) argue that the underlying fluid mechanical problem is susceptible to a phenomenon of ‘explosive instability’.

For a higher-order example, consider the fourth-order differential operator,(6.5)with periodic boundary conditions and symbol . For x=−π/2, the symbol curve is the quartic k4−ik, enclosing each point λ inside with winding number 3 (once by the quartic itself, twice more by the fourth power of a large semicircle at ∞). As x increases, for any λ in this region, the winding number decreases to 2 when the curve crosses once and then to 1 as it crosses a second time. Thus, for each such λ, we expect exponentially good pseudomodes consisting of a pair of wave packets (figure 5). In the special case Imλ=0, both crossings occur at the same value of x*. (Theorem 3.1 as written does not apply in this case, but that is an accident of wording, for in fact its proof is valid in such cases of multiple crossings so long as there is a net decrease in winding number.) In this special case, there will be pseudomodes consisting of two wave packets superimposed at the same x* and with opposite values of k*. This explains the lack of a smooth envelope in the figure.

Figure 5

Above, eigenvalues and ϵ-pseudospectra of the fourth-order operator (6.5) with h=0.04, ϵ=10−2,…,10−10. By theorem 3.1, ‖(λAh)−1‖ grows exponentially as h→0 for any λ lying in the quartic region marked by the dashed line. Below, an optimal pseudomode for λ=1 (cross), with x*=0 and k* taking both values 1 and −1.

We now move to a different type of example, in which the symbol curve at a particular value x* passes through the value λ of interest for two different wavenumbers k*. In such cases, theorem 3.1 does not apply, but since our coefficients depend analytically on x, we can appeal instead to theorem 1.2 of Dencker et al. (2004), which ensures as before the existence of pseudoeigenfunctions, now localized in k as well as x, that are exponentially good in the sense of (3.6).

Our fourth example is Davies's non-selfadjoint harmonic oscillator (1999a,b). Written with a small parameter h, we have(6.6)with symbol . For any fixed x*R, the symbol curve is the half-line in the complex plane traversed from ∞ to and back again to ∞. For each λ along this half-line and corresponding choice of x*, there are two values of k*, one of which satisfies the twist condition (the one whose sign is the same as that of x*). We can see this either by calculating the twist ratio , or by thinking of sections of the symbol curve. We conclude that every λC with Re λ>0, Imλ>0 is an ϵ-pseudoeigenvalue of Ah for an exponentially small value of ϵ, as shown in figure 6. There are two values of x* for each λ, which explains why the optimal pseudomode in the figure consists of two wave packets rather than one.

Figure 6

Above, eigenvalues and ϵ-pseudospectra of the Davies example (6.6) with h=1/10, ϵ=10−1,…,10−13. By the results of Dencker et al. (2004), ‖(λAh)−1‖ grows exponentially as h→0 for any λ lying in the first quadrant of the complex plane. Below, an optimal pseudoeigenfunction for λ=3+2i (cross), with and .

It is worth commenting further on the significance of a case like figure 6 in which two wave packets appear in a computed pseudomode. The arguments in this paper and in (Dencker et al. 2004) construct exponentially good pseudomodes in the form of single wave packets, not double ones. The present case is special because there are two values of x* that are equally good for this construction. In some sense, the multiplicity of the pseudoeigenvalue is two rather than the usual value of one. Thus, it would be equally valid to show a pseudomode with just one wave packet on the left, or on the right; except that the optimal pseudomode, less than 0.003% better than these, is the odd function with two bumps. (Orthogonal to this would be another function with two bumps, but even instead of odd, and little different to the eye.)

Davies & Kuijlaars (2004) have analysed the operator (6.6) in detail (with i replaced by an arbitrary complex constant), basing their arguments on the theory of polynomials orthogonal with respect to a complex weight function. Among other results, their theorem 3 implies that the condition numbers κ(λn) of the eigenvalues of (6.6) grow exponentially at the ratewith eigenvalues indexed with increasing distance from the origin. Such a precise estimate goes beyond the results presented here or in Dencker et al. (2004), where the constants are not specified.

Our next example is closely related to Davies's, the only difference being that the coefficient ix2 is replaced by ix3. This ‘complex cubic oscillator’ is a representative of a class of operators that have been discussed by Bender and others, starting from unpublished work of D. Bessis in 1995, for applications in non-Hermitian quantum mechanics (Bender & Boettcher 1998; Bender et al. 1999; Delabaere & Trinh 2000; Mezincescu 2000; Handy 2001). The equation is(6.7)with symbol . Mathematically, this is much the same as the Davies example, but the pseudospectra fill the right half-plane instead of the first quadrant since x3 ranges over all of R rather than just [0,∞) (figure 7). Most of this literature is concerned with establishing properties of the eigenvalues of (6.7) and related operators and does not question their physical significance. Again, the results of Dencker et al. apply to this operator, but theorem 3.1 does not.

Figure 7

Above, eigenvalues and ϵ-pseudospectra of the complex cubic oscillator (6.7) of Bender & Boettcher (1998) and Bender et al. (1999) with h=0.1, ϵ=10−1,…,10−13. By the results of Dencker et al., ‖(λAh)−1‖ grows exponentially as h→0 for any λ lying in the right half-plane. Below, an optimal pseudoeigenfunction for λ=2+i (cross), with x*=1 and .

Another operator with the same flavour as those of Davies and Bender is(6.8)where α and γ are parameters. This has no particular physical significance, but the pseudospectra, for α=3+3i and γ=1/16, are explored in Trefethen (1999). As with the Cossu–Chomaz operator, there are significant transient effects, with all eigenvalues in the left half-plane but large resolvent norms in the right half-plane.

None of the examples we have presented were the first examples of variable coefficient differential operators whose pseudospectra were computed numerically. That distinction belongs to the Orr–Sommerfeld and Airy operators, whose pseudospectra were computed together by Reddy et al. (1993). The Orr–Sommerfeld operator is in generalized eigenvalue form and does not quite fit the framework of this paper. The Airy operator is(6.9)with boundary conditions u(−1)=u(1)=0 and symbol . Figure 8 presents results corresponding to h=0.02. This operator has also been investigated by Stoller et al. (1991), Shkalikov (1997) and Redparth (2001).

Figure 8

Above, eigenvalues and ϵ-pseudospectra of the Airy operator (6.9) with h=0.02, ϵ=10−1,…,10−13. By the results of Dencker et al., ‖(λAh)−1‖ grows exponentially as h→0 for any λ lying in the dashed half-strip. Below, an optimal pseudoeigenfunction for λ=−1+0.5i (cross), with x*=0.5 and k*=−1.

## 7. Robustness and structural stability

For all of our six examples, we have a theorem to explain why exponentially good wave packet pseudoeigenfunctions had to appear: theorem 1.2 of Dencker et al. (2004) for all six, and theorem 3.1 of this paper for just the first three, because the others had double-crossing symbol curves. Now, the theorem of Dencker et al. requires smoothness of the symbol, whereas theorem 3.1 does not. Therefore, one might expect that the first three examples should be robust in the sense that the exponentially good pseudomodes persist if the operator coefficients undergo a non-smooth perturbation, whereas the other three may be fragile.

We found confirming this prediction numerically to be challenging, for the computations underlying figures 1–5 are based on spectral methods (Trefethen 2000), a technology that relies on smooth functions for its power, whereas if one reverts to simpler finite differences or finite elements, based on weaker smoothness assumptions, the accuracy is usually too low to resolve ϵ-pseudospectra for small values of ϵ.2 The compromise we eventually reached was to continue to use spectral methods but to choose perturbations that have five continuous derivatives—enough smoothness so that the numerical method is still quite accurate, but enough roughness to excite the effects of interest. In a typical case, we perturbed one coefficient of the ODE by multiplying it by the variable coefficientand the other coefficient by multiplying it byThe results appear in figure 9. As predicted, three of the cases shown are robust and three are fragile, with pseudospectra distorted almost beyond recognition by this C5 perturbation. Theorem 3.1 ensures that the robust cases would in fact stand up to far rougher perturbations; it is just difficult to verify such cases numerically.

Figure 9

A repetition of our six examples with each operator coefficient modified by a C5 multiplicative perturbation. The first three have simply crossing symbol curves, and the perturbation has little effect on the pseudospectra (theorem 3.1). For the next three, it changes them completely.

The same distinction between robust and fragile pseudospectra arises for twisted Toeplitz matrices and is discussed in Trefethen & Chapman (2004; see figs 7.4 and 8.1 of that paper). We must emphasize, however, that the experiments of this section only go far enough to suggest that this distinction is a real one for variable coefficient differential operators. Because of the challenges in these numerical simulations as indicated, further study would be needed to elucidate more of the details of these effects and to establish whether they are indeed generic.

## 8. Discussion

The results we have presented are restricted to ordinary differential operators in one space dimension. If the symbol is smooth, the WKBJ/microlocal approach of Dencker et al. (2004) generalizes to partial differential or pseudodifferential operators, and their results are presented in this framework. On the other hand, it is not clear whether our results for non-smooth symbols can be generalized in this way.

Another kind of generalization would concern discontinuities of coefficients or boundaries of the domain. For twisted Toeplitz matrices, such generalizations are considered at length in Trefethen & Chapman (2004), where it is pointed out that earlier results of Reichel & Trefethen (1992) on pseudomodes pinned at boundaries can be derived as a special case of the wave packet theory for discontinuous coefficients. For differential operators, one could presumably work out analogous results, deriving the results of Reddy (1993) and Davies (2000) for boundary pseudomodes from a wave packet theory. Such an approach would make it clear that the existence of boundary pseudomodes does not depend upon constant coefficients, a matter first taken up in Davies (in press).

Related to boundaries and discontinuities is the matter of the antitwist condition, which we have mentioned in theorem 3.2 but not illustrated in our examples. Suppose one has a family of operators such that as x increases, the symbol curves sweep across λC to decrease the winding number but do not sweep back again to increase the winding number. Then λ will be an exponentially good pseudoeigenvalue with wave packet pseudoeigenfunctions, whereas for the complex conjugates of the adjoint operators, the pseudoeigenfunctions will be localized elsewhere, often at the boundary.

We believe that the existence of exponentially good pseudoeigenmodes will prove to have applications to many non-selfadjoint dynamical problems in the mathematical sciences involving variable coefficients, not just those listed in the first paragraph of this article. We hope that the theory and examples presented here will help to bring these phenomena to wider attention.

## Acknowledgments

Much of this article was completed during a visit to the National University of Singapore in December 2003; I thank Toh Kim-Chuan and Lee Seng-Luan and the other stimulating mathematicians at the NUS for supporting this visit. Concerning the mathematics of wave packet pseudomodes, I have benefited from interactions with many people, of whom three have been especially important: Jon Chapman of Oxford University, Brian Davies of King's College London and Maciej Zworski of the University of California, Berkeley.

## Footnotes

• An ϵ-pseudoeigenvalue and corresponding ϵ-pseudoeigenfunction of an operator A are a scalar λ and a non-zero function u satisfying ‖(Aλ)u‖<ϵu‖, and the ϵ-pseudospectrum of A is the union of its spectrum and the set of all of its ϵ-pseudoeigenvalues; that is, the set of all λC with resolvent norm satisfying ‖(λA−1)‖>ϵ−1, if we adopt the convention of writing ‖(λA)−1‖=∞ if λ is in the spectrum (Trefethen 1997; Trefethen & Embree 2005).

• Indeed, a further difficulty arises here. Suppose one discretizes a smooth differential equation by a finite difference approximation on a uniform grid. The result is a twisted Toeplitz matrix as in Trefethen & Chapman (2004). However, if the symbol curve for the differential equation has no crossings f(x*,k)=f(x*,k*) for kk*, so that theorem 3.1 is applicable, then this does not imply the same property for the matrix approximation. Thus, even when a differential equation has pseudospectra that are robust with respect to perturbations, those of its finite difference approximations will often be fragile.

• Received August 7, 2004.
• Accepted April 20, 2005.

View Abstract