Complex eigenvalue splitting for the Dirac operator

We analyze the eigenvalue problem for the semiclassical Dirac (or Zakharov-Shabat) operator on the real line with general analytic potential. We provide Bohr-Sommerfeld quantization conditions near energy levels where the potential exhibits the characteristics of a single or double bump function. From these conditions we infer that near energy levels where the potential (or rather its square) looks like a single bump function, all eigenvalues are purely imaginary. For even or odd potentials we infer that near energy levels where the square of the potential looks like a double bump function, eigenvalues split in pairs exponentially close to reference points on the imaginary axis. For even potentials this splitting is vertical and for odd potentials it is horizontal, meaning that all such eigenvalues are purely imaginary when the potential is even, and no such eigenvalue is purely imaginary when the potential is odd.


Introduction
Consider the eigenvalue problem (1.1) P (h)u = λu on the real line for the Dirac (or Zakharov-Shabat) operator given by the 2 × 2 non-selfadjoint system where u is a column vector, h a small positive parameter, λ a spectral parameter, and V a real-valued analytic function on R. Solving (1.1) constitutes an essential step in the treatment of many important nonlinear evolution equations by means of the inverse scattering transform, including the focusing nonlinear Schrödinger (NLS) equation, the sine-Gordon equation and the modified Korteweg-de Vries equation [7]. Among the numerous applications of these equations are nonlinear wave propagation in plasma physics, nonlinear fiber optics, hydrodynamics and astrophysics. The operator P (h) is the massless Dirac operator on the real line with antiselfadjoint potential. In the selfadjoint case, resonances have been studied in various settings by many authors, see for example [17] for a historical account of the massless case, and [18] for the massive case. Certain types of massless Dirac operators have also been shown to be effective models for twistronics, such as in Twisted Bilayer Graphene (TBG), see e.g. [30]. In fact, the twist-angles producing superconductive properties in TBG can be characterized in terms of the spectrum of the related Dirac operator [2,3], which in this case is highly non-normal.
Here, we shall focus on the connection between P (h) and the NLS equation, which is one of the most fundamental nonlinear evolution equations in physics. In the focusing semiclassical case one is interested in the asymptotic behavior of ψ = ψ(t, x; h) in the semiclassical limit h → 0, where ψ is the solution to the initial value problem (1.2) ih ∂ψ ∂t + h 2 2 and V is a real-valued function independent of h. In the inverse scattering method the initial data is substituted by the soliton ensembles data, defined by replacing the scattering data for ψ(0, x) = V (x) with their formal WKB approximation. The focusing NLS equation (1.2) is then solved with this new set of h-dependent initial data, and the asymptotic behavior of the obtained approximate solution is analyzed in the limit h → 0. However, it is a priori not clear how such an h-dependent approximation of initial data affects the behavior of ψ as h → 0, or if it is even justified at all. For this a rigorous semiclassical description of the spectrum of the corresponding Dirac operator P (h) is required, which has so far only been provided in a few cases such as for periodic potentials by Fujiié and Wittsten [12], and for bell-shaped, even potentials by Fujiié and Kamvissis [9]. Both of the mentioned articles employ the exact WKB method which we describe in Section 2 below. For an in-depth discussion on the necessity (as well as effects) of a precise description of the semiclassical spectral data of P (h) in the context of inverse scattering and the focusing NLS equation we refer to the second paper mentioned above. The interest in the spectrum of the operator P (h) and its relatives dates back to Zakharov and Shabat [27]. Since P (h) is not selfadjoint the eigenvalues are not expected to be real in general. These complex eigenvalues directly determine the energy and speed of the soliton (solitary wave) solutions of (1.2); the energy, or amplitude, given by the imaginary part and the speed by the real part of the eigenvalue. Early on it was realized that there are examples of potentials V (x) for which all the complex eigenvalues are in fact purely imaginary, thus giving rise to soliton pulses with zero velocity in the considered frame of reference. (In the defocusing case, obtained from (1.2) by changing sign of the nonlinear term, no such solutions exist in general. In fact, the corresponding Dirac operator is then selfadjoint, and the first author has shown that it has real spectrum even under small non-selfadjoint perturbations [16].) In 1974, Satsuma and Yajima [26] studied P (h) with V (x) = V 0 sech(x), V 0 > 0, and solved (1.1) by reducing it to the hypergeometric equation. They found that if h = h N = V 0 /N , there are exactly N purely imaginary eigenvalues λ k given by For many years thereafter, the literature was filled with erroneous statements about eigenvalues being confined to the imaginary axis whenever the potential V is realvalued and symmetric. In the nonsemiclassical regime (h = 1) the question was given rigorous consideration in a series of papers by Klaus and Shaw [20,21,22] who established that (a) if V is of Klaus-Shaw type, that is, a "single-lobe" potential defined by a non-negative, piecewise smooth, bounded L 1 function on the real line which is nondecreasing for x < 0 and nonincreasing for x > 0, then all eigenvalues are purely imaginary (symmetry not being a factor); (b) there are examples of real-valued, even, piecewise constant or piecewise quadratic potentials with two or more "lobes" giving rise to eigenvalues that are not purely imaginary; (c) if V ∈ L 1 is an odd function, there are no purely imaginary eigenvalues at all. We shall consider these questions in the semiclassical setting and analytic category, and show that a counterpart of (a) holds for eigenvalues near λ 0 = iµ 0 ∈ iR even if one only assumes that V locally has the shape of a single-lobe 1 potential near the "energy level" µ 0 . We will also derive precise conditions for eigenvalues when V locally has the shape of a double-lobe potential near the energy level µ 0 , and show that when V is symmetric, this leads to an exponentially small splitting of the eigenvalues akin to the well-known splitting phenomenon observed for eigenvalues of the selfadjoint Schrödinger operator with a double-well potential. We prove that when V is even and h > 0 is sufficiently small, this splitting is vertical from reference points on the imaginary axis; in particular, all eigenvalues are purely imaginary then. (This is in contrast to the examples in (b) which of course do not satisfy the analyticity assumption, and we believe this might help explain the confusion witnessed in the literature prior to the mentioned papers by Klaus and Shaw.) We also show that when V is odd and h > 0 is sufficiently small, the splitting is horizontal from reference points on the imaginary axis; in particular, in accordance with (c) there can be no purely imaginary eigenvalues in this case. Here we note that for fixed h, (1.1) can be formally interpreted as a non-semiclassical Zakharov-Shabat eigenvalue problem with potential q(x) = h −1 V (x) and spectral parameter ζ = h −1 λ, so it makes sense to compare results between the two settings. In particular, the eigenvalue formation threshold ∞ −∞ |q(x)| dx > π/2 established by Klaus and Shaw [22] is always reached as h → 0. We also wish to mention that some of the examples in (b) together with the corresponding focusing NLS equation have been further analyzed by Desaix, Andersson, Helczynski, and Lisak [5], and Jenkins and McLaughlin [19], among others.

Statement of results.
To be more precise, we shall view P (h) as a densely defined operator on L 2 and study the eigenvalue problem (1.1) for spectral parameters λ = iµ close to λ 0 = iµ 0 ∈ i(0, V 0 ), where V 0 = max x∈R |V (x)|, for which the potential is either a single or double lobe in a sense to be specified below. We assume that the potential satisfies the following assumptions: is real-valued on R and analytic in a complex domain D ⊂ C containing an open neighborhood of the real line, and (ii) lim sup x→±∞ |V (x)| < µ 0 . Examples of D are tubular neighborhoods of R, or more generally, domains {x ∈ C : |Im x| < δ(Re x)} where δ : R → R + is a positive continuous function which is allowed to decay as |Re x| → ∞. Note that the spectrum of P (h) is symmetric with respect to reflection in R (as well as with respect to reflections in the imaginary 1 Here lobe is terminology adopted from Klaus and Shaw referring to a projecting or hanging part of something, like in earlobe, or the lobe of a leaf. axis), so it is not necessary to treat λ 0 ∈ i(−V 0 , 0) separately. We will also not consider spectral parameters close to the real line. In fact, if (ii) is strengthened to a decay condition of the form then it is known that the continuous spectrum of P (h) consists of the entire real axis, and that away from the origin there are no real eigenvalues. For potentials of Klaus-Shaw type satisfying (ii) , a precise description of the reflection coefficients as well as the eigenvalues close to zero has recently been obtained by Fujiié and Kamvissis [9]. Finally, it is not necessary to consider eigenvalues away from R i[−V 0 , V 0 ] since the spectrum of P (h) accumulates on this set in the limit as h → 0. In fact, if has no spectrum in Ω as long as h is sufficiently small, see Dencker [4,Section 2] or Fujiié and Wittsten [12, Proposition 2.1]. After obtaining the necessary properties in Section 2 of the exact WKB solutions needed for our analysis, we shall therefore in Section 3 study the spectrum of P (h) near λ 0 = iµ 0 when the potential locally, near the energy level µ 0 , corresponds to a single lobe in the following sense. Definition 1.1. Let 0 < µ 0 < V 0 and assume that the equation V (x) 2 − µ 2 0 = 0 has exactly two real solutions α l (µ 0 ) and α r (µ 0 ) with α l < α r and V (α • ) = 0, • = l, r. We then say that V is a single-lobe potential near µ 0 .
We may without loss of generality assume that the roots of the equation V (x) 2 − µ 2 = 0 (called turning points) are roots to V (x) = µ 0 (so that V (α l ) > 0 and V (α r ) < 0) since the case when they are roots to V (x) = −µ 0 can be studied by replacing the potential V with −V and reducing the resulting eigenvalue problem to the original one. 2 (It is not possible that one turning point is a root to V (x) = µ 0 and the other to V (x) = −µ 0 since this would give four solutions to V (x) 2 − µ 2 0 = 0 when V (α l ), V (α r ) = 0.) Figure 1 illustrates two stereotypical examples of this situation. Of course, any potential V of Klaus-Shaw type is a single-lobe potential near µ 0 ∈ (0, V 0 ). Note also that the turning points depend continuously (even analytically) on µ as long as the multiplicity is constant. In particular, Definition 1.1 cannot hold at µ 0 = 0 or µ 0 = V 0 (or at any local extreme values of V ) because the situation degenerates then, which explains why these values are excluded. It also explains why it makes sense to say that V is a single-lobe potential near µ 0 , since there is then an ε-neighborhood B ε (µ 0 ) ⊂ C around µ 0 such that if µ ∈ B ε (µ 0 ) then the equation V (x) 2 − µ 2 = 0 has exactly two solutions α l (µ) and α r (µ) with Re α l < Re α r , Re V (α l ) > 0 and Re V (α r ) < 0. For such µ we define the action This can of course also be realized by noting that if ψ solves (1. where the determination of the square root is chosen so that I(µ) is real and positive for real µ. In this case, we prove in Section 3 that there are constants ε, h 0 > 0 such that if µ ∈ B ε (µ 0 ) and 0 < h ≤ h 0 then λ = iµ is an eigenvalue if and only if the Bohr-Sommerfeld quantization condition is satisfied for some integer k, see Theorem 3.4. Here r(µ, h) is a function defined on B ε (µ 0 ) × (0, h 0 ] with r = O(1) as h → 0. In particular, if µ sl k (h) is the unique root of I(µ) = (k + 1 2 )πh near µ 0 (where the superscript sl refers to single lobe), and λ sl k = iµ sl k , then there is a unique eigenvalue λ k = iµ k such that |λ k (h) − λ sl k (h)| = O(h 2 ), see Remark 3.5. We also obtain the following refinement of [9,Theorem 2.2] showing that for single-lobe potentials, the semiclassical eigenvalues are confined to the imaginary axis: If V is a single-lobe potential near µ 0 = −iλ 0 , then there exist positive constants h 0 and ε such that the point spectrum of P (h) satisfies Section 4 studies the eigenvalue problem for potentials assumed to locally have the features of a double lobe. Definition 1.3. Let 0 < µ 0 < V 0 and assume that the equation V (x) 2 −µ 2 0 = 0 has exactly four real solutions α l (µ 0 ), β l (µ 0 ), β r (µ 0 ) and α r (µ 0 ) with α l < β l < β r < α r and V (α • ), V (β • ) = 0, • = l, r. We then say that V is a double-lobe potential near µ 0 . Figure 2 shows two stereotypical examples of double-lobe potentials. In the first example, V (β l ) = V (β r ) > 0, whereas V (β l ) = −V (β r ) > 0 in the second. As indicated, it suffices to consider these two situations (i.e., peak-peak and peak-valley) since the other two cases can be obtained, as for single-lobe potentials, by replacing the potential V by −V and reducing the corresponding eigenvalue problem to the original one. By continuity there is an ε > 0 such that for µ ∈ B ε (µ 0 ), the equation V (x) 2 − µ 2 = 0 still has exactly four solutions α l (µ), β l (µ), β r (µ) and α r (µ) with Re α l < Re β l < Re β r < Re α r and the signs of For such µ we introduce the action integrals where the determinations of the square roots are chosen in such a way that each action integral is real-valued and positive for real µ. We show that there are positive constants ε, h 0 such that if µ ∈ B ε (µ 0 ) and 0 < h ≤ h 0 then λ = iµ is an eigenvalue in the case when V (β l ) = ±V (β r ) if and only if as h → 0, and * denotes the operation γ * l (µ) = γ l (μ), see §2.5. From the quantization condition (1.7) we see that modulo an exponentially small error the eigenvalues λ = iµ for µ ∈ B ε (µ 0 ) are given in terms of the roots to the equation This is equivalent to the two Bohr-Sommerfeld quantization conditions corresponding to each potential lobe, i.e., These may be rewritten in the form are both bounded when h tends to 0. Thus we conclude that the set of eigenvalues produced by a double-lobe potential is exponentially close to the union of the sets of eigenvalues produced by each potential lobe (cf. (1.4)). This is a well-known fact for the Schrödinger equation, see [15] and [13]. Remark. For readers familiar with the time-independent Schrödinger equation we wish to mention that "inside" the lobe(s) (the projection of the shaded regions in Figures 1-2 onto the real axis), solutions to (1.1) are oscillating, while they are exponential in character "outside" the lobe(s). In this sense, the lobes can thus be said to correspond to potential wells (rather than to barriers) in the terminology of quantum mechanics.
Section 5 considers the special case of double-lobe potentials V such that V (x) is either an even or an odd function of x ∈ R. If this assumption holds, the quantization condition (1.7) can be rewritten in the case when V (x) = ±V (−x) as for some integer k = k(h). If µ dl k (h) is the unique root of (1.9) near µ 0 (where the superscript dl stands for double lobe), it turns out that iµ dl k is purely imaginary. Now, eigenvalues λ = iµ of the Dirac operator (where µ satisfies the quantization condition (1.8)) split in pairs symmetrically about the reference points iµ dl k . Theorem 1.4. Suppose that V is a double-lobe potential near µ 0 such that V (x) is either an even or an odd function of x ∈ R, and let µ dl k (h) be the unique root of (1.9) near µ 0 . Then iµ dl k ∈ iR and the two eigenvalues have the following asymptotic behavior as h → 0: (2) If V (x) is an odd function, then Moreover, the eigenvalues split precisely vertically in the even case, whereas they split precisely horizontally in the odd case. Thus, for 0 < h ≤ h 0 , all eigenvalues are purely imaginary when V is even, and no eigenvalue is purely imaginary when V is odd.
The proof relies on the explicit exponential error term in (1.8) which we obtain by using a novel method, inspired by recent work due to Mecherout, Boussekkine, Ramond and Sjöstrand [24], to refine the WKB analysis for the Dirac operator by introducing carefully chosen WKB solutions defined "between" the lobes. As already mentioned, the results are reminiscent of the well-known splitting of eigenvalues for the linear Schrödinger operator with a symmetric double-well potential, going back to the work of Landau and Lifshitz [23] and studied mathematically by, among others, Simon [28], Helffer and Sjöstrand [15] and Gérard and Grigis [13]. This type of tunneling effect has recently also been observed for a system of semiclassical Schrödinger operators by Assal and Fujiié [1]. For more on this topic we refer to the mentioned works and the references therein.
In the literature a common focus of study is the appearance and location of purely imaginary eigenvalues as the L 1 norm of the potential increases, for example by taking q(x) = h −1 V (x) and letting h decrease. Potentials of the form consisting of two separated sech-shaped pulses have been numerically investigated by Desaix, Anderson and Lisak [6] for different separations x 0 . They found that at the first critical amplitude h −1 = 1/4, a purely imaginary eigenvalue ζ 1 appears, and for h −1 < 1/4 there are no eigenvalues (consistent with the threshold of Klaus and Shaw [22]). For small separations, q behaves almost like a single-lobe potential, and the second critical amplitude h −1 = 3/4 also gives rise to a purely imaginary eigenvalue. However, for larger separations such as x 0 = 5, two complex eigenvalues ζ 2,3 = ±ξ + iη with nonzero real parts are created already in the vicinity of h −1 = 4/10. As the amplitude h −1 increases, the real parts decrease while η increases until the two eigenvalues meet and then separate along the imaginary axis (both now purely imaginary, ζ 2 with increasing and ζ 3 with decreasing imaginary part). As h −1 reaches the second critical amplitude 3/4, ζ 3 is destroyed and only ζ 1 and ζ 2 remain. Since ζ = h −1 λ, we should be able to see a similar type of behavior for semiclassical eigenvalues of P (h) as h decreases, which is something we hope to investigate in a future paper. Of course, as h becomes sufficiently small, our results show that for a potential consisting of two separated sech-pulses, all eigenvalues λ = iµ are purely imaginary as long as µ = 0 is not close to a local extreme value of V , see Figure 3. It would also be interesting to see if this property persists as µ tends to local extreme values of V ; the exact WKB method does not work then so other methods are needed.

Exact WKB analysis
2.1. Exact WKB solutions. Here we recall the construction of a solution of the Dirac system in a complex domain as a convergent series, known as an exact WKB solution. Such solutions were first introduced by Ecalle [8] and later used by Gérard and Grigis [13] to study the Schrödinger operator. We shall follow the construction for systems due to Fujiié, Lasser and Nédélec [10].  Figure 3. An even potential V (x) with a local minimum at x = 0. Away from the shaded region V is either a single-lobe or a doublelobe potential. For sufficiently small h, any eigenvalue λ of P (h) with imaginary part away from the shaded region must therefore be purely imaginary by Theorems 1.2 and 1.4.
The system (1.1) can be written in the form Recall (see [10]) that the exact WKB solutions of systems of type (2.1) are of the form where the function z(x) is the complex change of coordinates for some choice of phase base point x 0 in the strip D where V is assumed to be analytic, while Q is the matrix valued function Here z(x) and Q(z) are defined on the Riemann surfaces of (V 2 +λ 2 ) 1/2 and H(z(.)) over D, respectively. These Riemann surfaces are defined by introducing branch cuts emanating from the zeros of x → det(M (x, λ)), i.e., of iV ± λ (the turning points of the system (2.1)), see §2.4. The amplitude vectors w ± in (2.2) are defined as the (formal) series , where w ± 0 (z) ≡ 1, while w ± j (z) for j ≥ 1 are the unique solutions to the scalar transport equations with prescribed initial conditions w ± n (z) = 0 for some choice of amplitude base point z = z(x) wherex is not a turning point. When we want to signify the dependence on the base pointz = z(x) we write Recall that if Ω is a simply connected open subset of D which is free from turning points then z = z(x) is conformal from Ω onto z(Ω). For fixed h > 0, the formal series (2.4) converges uniformly in a neighborhood of the amplitude base pointx, and w ± even (x, h) and w ± odd (x, h) are analytic functions in Ω, see [10, Lemma 3.2]. As a consequence, the functions u ± given by (2.2) are exact solutions of (2.1) and when we wish to indicate the particular choice of amplitude base pointx ∈ Ω and phase base point x 0 ∈ D we will write u ± (x; x 0 ,x). We remark that these solutions are defined for example everywhere on R, although some of the expressions involved are only defined on Riemann surfaces of (V 2 + λ 2 ) 1/2 or H(z(.)).
For fixedx ∈ Ω, let Ω ± be the set of points x for which there is a path fromx to x along which t → ± Re z(t) is strictly increasing. In other words, x ∈ Ω ± if there is a path which intersects the the level curves of t → Re z(t) transversally in the appropriate direction. The level curves of t → Re z(t) are called Stokes lines.
uniformly on compact subsets of Ω ± as h → 0, see [10, Proposition 3.3]. In particular, 2.2. The Wronskian formula. For vector-valued solutions u and v of (2.1), let W(u, v) be the Wronskian defined by Since the trace of the matrix M (x, λ) is zero, it follows that W(u, v) is in fact independent of x. If x 0 is a phase base point in D andx,ỹ are different amplitude base points in Ω, a straightforward calculation shows that , where the solutions u ± are given by (2.2). Recalling the initial conditions of the transport equations (2.5)-(2.6) and evaluating at x =ỹ we get . We may of course also choose x =x, which gives In particular, we see that if there is a path fromx toỹ along which the function showing that such a pair of solutions is linearly independent if h is sufficiently small. We also recall the Wronskian formula for pairs of solutions of the same type: 2.3. Stokes geometry. We now describe the configuration of Stokes lines for single-lobe and double-lobe potentials.
2.3.1. Single-lobe potentials. Suppose that V is a single-lobe potential near µ 0 and let µ ∈ B ε (µ 0 ). Fix determinations of H(z(x)) given by (2.3) and of Note that this is in accordance with (1.3). The Stokes lines (level curves of t → Re z(t; α • )) are then found by taking the union of for x 0 = α l , α r . When µ is real it is known that there are three Stokes lines emanating from α l ∈ R having arguments 0, 2π/3, 4π/3, while the Stokes lines emanating from α r ∈ R have arguments π/3, π, 5π/3, see Gérard and Grigis [13]. We define the Riemann surfaces of z(x) and H(z(x)) by introducing branch cuts along the Stokes line with argument 2π/3 at α l and the Stokes line with argument 5π/3 at α r . For real µ ∈ B ε (µ 0 ) there is a bounded Stokes line lying on R starting at α l and ending at α r . Hence, the Stokes lines separate the complex domain D into four sectors (called Stokes regions). In the top and bottom sectors the function z(x) takes the form (2.10). By continuing the chosen determination of z(x) through rotation clockwise around the turning points (thus avoiding the branch cuts) it is easy to see that for x belonging to the left and right sector when • = l and • = r, respectively. For general µ ∈ B ε (µ 0 ) the picture is slightly perturbed; as iµ is rotated off the imaginary axis α l and α r start migrating in opposite directions along paths in the upper and lower half plane, and the bounded Stokes line connecting α l and α r is broken into two unbounded curves, see Figure 4. (We refer to [11] for a detailed explanation of this phenomenon.) However, for small ε the arguments of the Stokes lines at the turning points are almost unchanged so for µ ∈ B ε (µ 0 ) we may still place branch cuts as described above. Note that there are now three Stokes regions around the left turning point and three around the right, and (2.12) is still valid if interpreted in this sense. However, we will avoid introducing notation for the  different Stokes regions, and simply say (informally) that x is near the lobe if x is not in the Stokes region to the left of α l or to the right of α r . We also remark that if x 0 (µ) is a turning point satisfying V (x 0 (µ)) = µ, then x 0 (−µ) is also a solution to V (x) 2 − µ 2 = 0; hence the original Stokes configuration is reached again already when iµ has traversed half a circuit around the origin, see the left panel of Figure  6 below.
In Figure 4 we have also indicated that Re z(x) increases as one travels from top to bottom and left to right, while not passing through a branch cut. This is realized in the following way: For x in the regions between turning points we have by (2.10) and Taylor's formula that where g 1 is analytic and g 1 (x 0 ) = 0. Since V (x 0 ) > µ 0 if α l (µ 0 ) < x 0 < α r (µ 0 ) we see by picking x 0 real that the square root is approximately real when µ ∈ B ε (µ 0 ), so Re z(x) increases as Im x decreases. On the other hand, for x in the Stokes region left of α l or right of α r we have by (2.12) and Taylor's formula that where g 2 is analytic and g 2 (x 0 ) = 0. By picking increases as Re x increases. This also shows that Re z(x) is constant along lines which are essentially vertical near R when |Re x| is large.

2.3.2.
Double-lobe potentials. Suppose now that V is a double-lobe potential near µ 0 and let µ ∈ B ε (µ 0 ). Again, fix determinations of H(z(x)) and z(x) in accordance with (1.5)-(1.6); the obtained configuration of Stokes lines will essentially be two side-by-side copies of the configuration for single-lobe potentials with an appropriate gluing in the region between the two middle turning points β l and β r .
Indeed, the Stokes lines are given by the union of (2.11) for x 0 = α l , β l , β r , α r . When µ is real there are three Stokes lines emanating from α l and three from β r having arguments 0, 2π/3, 4π/3, while the Stokes lines emanating from β l , α r ∈ R have arguments π/3, π, 5π/3, see Gérard and Grigis [13]. As iµ is rotated off the imaginary axis the turning points start migrating in alternating, opposing directions along paths in the upper and lower half plane, so that α l moves in the direction opposite from β l but similar to β r . We place branch cuts along the Stokes lines which for real µ have arguments 2π/3 at α l , β r and the Stokes lines with arguments 5π/3 at β l , α r . Performing the same analysis as above shows that in the sectors to the left of α l and to the right of α r , and in the intersection of the sectors to the right of β l and to the left of β r (i.e., between β l and β r ), z(x; α • ) takes the form (2.12). When x is in the other sectors (between α l and β l or between β r and α r ), z(x; α • ) is given by (2.10), and as for single-lobe potentials we shall informally say that x is near the lobes in this case. Using Taylor's formula as in (2.13)-(2.14) then shows that Re z(x) increases as one travels from top to bottom and left to right, while not passing through a branch cut, see Figure 5. The right panel of Figure 6 shows an example of how the turning points of a double-lobe potential migrate as iµ is rotated off the imaginary axis.
2.4. The Riemann surface. Let R(x 0 , θ) denote the operator acting through rotation around x 0 by θ radians, so that, e.g., R(0, θ)x = e iθ x. Since V − µ is analytic and V (α l ) − µ = 0 it follows that Figure 6. The migration paths of turning points (solutions to V (x) 2 − µ 2 = 0) of the potentials in Figure 4 (left) and Figure 5 (right) as iµ is rotated π radians in the positive direction from the starting value µ = 0.2 until µ = −0.2 when the original Stokes geometry is recovered. Black dots and circles mark starting and finishing locations, respectively. Rotation in the opposite direction reverses the direction of migration. Note that since sech(x + iπ) = − sech(x) the turning points appear periodically in C with complex period iπ for both potentials. Examples of the domain D (gray) are shown to indicate that only small rotations of iµ are of interest for the problem under consideration here.
i.e., when t is rotated 2πk radians anticlockwise around α l then V (t) − µ is rotated 2πk radians anticlockwise around the origin. (Negative k results in clockwise rotation by 2π|k| radians.) We of course have similar behavior near the other turning points of the same type, as well as for V + µ in the case when e.g. V (β r ) + µ = 0. Definition 2.2. Suppose that V is a single-lobe (double-lobe) potential near µ 0 and let y be a point in the upper half plane with Re α l < Re y < Re α r (Re α l < Re y < Re β l ). The point over y that is obtained when rotating y anticlockwise once around α l will be denoted byŷ, i.e., y = R(α l , 2π)y.
More generally, the sheet reached (from the usual sheet) by entering the cut starting at α l from the right will be referred to as thex-sheet. The point over y that is obtained when rotating y clockwise once around α l will be denoted byy, i.e., y = R(α l , −2π)y.
The sheet reached (from the usual sheet) by entering the cut starting at α l from the left will be referred to as thex-sheet.
Note that this definition is in accordance with [12, Definition 5.2]. When winding this way around a turning point we always assume that the path is appropriately deformed so that it is not obstructed by other branch cuts. Informally, we think ofx as lying in the sheet "above" the usual sheet, andx as lying in the sheet "below" the usual sheet. It is straightforward to check that thex-sheet is also reached (from the usual sheet) whenever we rotate anticlockwise once around the other zeros of V − µ (i.e., around β l , β r and α r if V (β l ) = V (β r ), and around β l if V (β l ) = −V (β r )). Similarly, thex-sheet is reached (from the usual sheet) by rotating clockwise once around zeros of V − µ. The directions are reversed when rotating around zeros of V +µ, i.e., when rotating around β r and α r if V (β l ) = −V (β r ). For a proof of these facts we refer to [12,Lemma 5.3]. We also record the following identities describing how WKB solutions are transformed when switching sheets. Lemma 2.3. [12,Lemma 5.4] Letx andx be defined as above and in accordance with Definition 2.2. Let x 0 be any of the turning points α l , β l , β r , α r , and let y be an amplitude base point. Then

2.5.
Symmetry. For constants c = c(λ) depending on the spectral parameter λ we shall simply write c(µ) with the convention that µ is always defined via λ = iµ. We then write c(μ) to represent the value of c at the reflection of λ in the imaginary axis, i.e., at iμ = −λ. We let c * (µ) = c(μ).
Recall that we fixed a determination of H(z(x)) so that if µ ∈ R then atx = (α l + β l )/2 ∈ R we have It is straightforward to check that for µ ∈ R, this determination implies that H(z(x)) ∈ R + when α l < x < β l , (2.15) H(z(x)) ∈ e iπ/4 R + when x < α l , β l < x < β r , α r < x, When V (β l ) = V (β r ) > 0 this is in accordance with the fact that for some constant c. Using the determination above we find that for µ ∈ R and x < α l , e iπ/4 R + H(−x) = c/H(x) ∈ ce −iπ/4 R + which implies that c = i. The same conclusion can also be drawn from the observation that if µ ∈ R and α l < x < β l then iR + H(−x) = c/H(x) ∈ cR + by (2.17), again showing that c = i, that is, These observations will be used to prove two symmetry properties: one with respect to reflection of the spectral parameter in the imaginary axis, and one with respect to parity.
Proposition 2.4. Let µ ∈ B ε (µ 0 ) and let x 0 (µ) ∈ C be a solution to V (x) 2 −µ 2 = 0. Then x 0 (µ) = x 0 (μ). Let y be an amplitude base point independent of µ. If V is a single-lobe, or a double-lobe with V (β l ) = V (β r ) then near the lobe(s) we have If V is a double-lobe with V (β l ) = −V (β r ) then (2.20) holds near the left lobe while near the right lobe. In the Stokes region to the left of α l or to the right of α r , Proof. Since V is real-analytic we have V (x) = V (x), which implies that V (α l (µ))− µ = 0. Since α l (μ) also satisfies this equation it follows that α l (µ) = α l (μ), for α l (µ 0 ) ∈ R and the turning points depend continuously on µ ∈ B ε (µ 0 ). Hence α * l (µ) = α l (µ). The same arguments show that x * 0 (µ) = x 0 (µ) when x 0 is any of the other three turning points.
Next, if V is a single-lobe and x lies in the domain between the turning points, or if V is a double-lobe and x lies in either the domain between the left pair or in the domain between the right pair of turning points, then z(x, µ) = i (V 2 − µ 2 ) 1/2 dt with real integrand when x, µ ∈ R. It is then easy to check that z(x,μ) = −z(x, µ). (In particular, when x and µ are real, z(x, µ) is purely imaginary, as expected.) If V is a single-lobe or a double-lobe with V (β l ) = ±V (β r ) then, using (2.15) or (2.17), one checks that H(z(x,μ)) = H(z(x, µ)) near the left lobe and H(z(x,μ)) = ±H(z(x, µ)) near the right lobe, which implies that with c = 1 near the left lobe, and c = ±1 near the right lobe, with sign determined according to V (β l ) = ±V (β r ). Since z (x,μ) = −z (x, µ), inspection of the governing equations for the amplitude function w ± (x, h; y, µ) shows that which gives (2.20)-(2.21). Finally, if x lies in the domain left of α l or right of α r then z(x, µ) = (µ 2 − V 2 ) 1/2 dt with real integrand when x, µ ∈ R, so z(x,μ) = z(x, µ). Using (2.16) one checks that H(z(x,μ)) = −iH(z(x, µ)) and Q(z(x,μ)) = iQ(z(x, µ)). Inspection of the governing equations for the amplitude function w ± (x, h; y, µ) shows that w ± (x, h; y,μ) = w ± (x, h;ȳ, µ). This proves the last statement of the proposition and the proof is complete.
Proof. Since we are only concerned with symmetry with respect to x → −x we will omit µ from the notation.
and H(z(−x)) = H(z(x)) by (2.18). The governing equations for the amplitude function w ± (x, h; y) imply that Noting that Q(z(−x)) = Q(z(x)) and − 0 1 1 0 and that squaring the right-most matrix gives the identity, we obtain the first formula.
If V (x) = −V (−x) for x ∈ R then z satisfies the same relations as above while H(z(−x)) = i/H(z(x)) by (2.19). The governing equations for w ± (x, h; y, µ) now give ) . Since the second formula therefore follows by checking that with Q(z(−x)) described above. This straightforward verification is left to the reader. Proof. We adapt the arguments in the proof of [16,Lemma IV.2]. Since We end this section with a result that will be used to determine the location of the reference points µ sl k and µ dl k mentioned in the introduction. In the statement, we let for brevity I(µ) denote either the action integral (1.3), or one of the action integrals I l , I r given by (1.5). It will be convenient to allow an error term which can be made exponentially small for any fixed h.
Proof. Note that since α l and β l depend analytically on µ and are roots to V (x) 2 − µ 2 = 0. At µ 0 ∈ R, this is a real integral with positive integrand. Hence, I l (µ 0 ) < 0, where prime denotes differentiation with respect to µ, and we can ensure that I l (µ) = 0 for µ ∈ B ε (µ 0 ) by choosing ε sufficiently small. The same arguments show that h), a contradiction. By assumption we have I * = I, so I(μ k , h) = I(µ k , h) = y k = y k = I(µ k , h) since y k is real. Since I is injective, we conclude that µ k =μ k .

Eigenvalues for a single-lobe potential
Here we suppose that V is a single-lobe potential near µ 0 , and let B ε (µ 0 ) be a small neighborhood as described in connection with Definition 1.1. We will consider eigenvalues λ = iµ with µ ∈ B ε (µ 0 ) with the purpose of deriving the quantization condition (1.4) and proving Theorem 1.2. We ask the reader to recall the relevant Stokes geometry described in §2.3 and Figure 4.
It is known that there are solutions u 0 = u 0 (µ) ∈ L 2 (R + ) and v 0 = v 0 (µ) ∈ L 2 (R − ) of (1.1) (unique modulo constant factors) such that λ = iµ is an eigenvalue of P (h) if and only if u 0 = cv 0 for some c = c(µ, h), thus As in the Schrödinger case, this can be shown by following the program of Olver [25] -see [14] for a detailed presentation in this direction.
easy to see that λ = iµ is an eigenvalue of P (h) if and only if u 0 = cc * v 0 . Since u * 0 = i u 0 , v * 0 = i v 0 and (cc * ) * = cc * , the claim follows.
To calculate the Wronskian (3.1), we shall follow [13, §2] and first show that modulo an exponentially small error, u 0 and v 0 are each multiples of exact WKB solutions. Pick real numbers x l ,x l andx r , x r such that x l <x l < Re α l and Re α r <x r < x r , and pick y in the upper half plane such that Re α l < Re y < Re α r , see Figure 7. These may be chosen independent of µ ∈ B ε (µ 0 ) if ε is small enough. Define two pairs of linearly independent exact WKB solutions u l , u l and u r , u r by setting with u ± given by (2.2). By (a slight modification of) the arguments of Gérard and Grigis [13, §2.2] we obtain the following representation formulas, where we use similar notation to make comparison easier. Note that by Proposition 2.4, u * r (x) = iu r (x) when x belongs to the Stokes domain containing x r , and the same is true for u r . By Remark 3.1 it is also true for u 0 . In the domain containing x l , the same relation holds for u l , u l , v 0 . A simple calculation then shows that (3.6) l * ± (µ, h) = l ± (µ, h), m * ± (µ, h) = m ± (µ, h). Next, introduce two pairs u + l , u − l and u + r , u − r of intermediate exact WKB solutions given as where I(µ) is the action integral (1.3). Represent v 0 and u 0 as the linear combinations v 0 = c 11 u + l + c 12 u − l , u 0 = c 21 u + r + c 22 u − r , where the coefficients c jk depend on the parameters µ and h. For x near the lobe we get for some symbols c l and c r .
Lemma 3.3. Let µ ∈ B ε (µ 0 ). For any A > 0 we may choose x r x r and x l x l so that the symbols c l and c r in (3.10)-(3.11) are given by Proof. Using (3.4)-(3.5) we see that .
For W(u l , u − l ), W(u + l , u − l ), W(u + r , u r ) and W(u + r , u − r ), we can directly apply the Wronskian formula (2.7), and obtain W(u l , u − l ) = 4iw + even (ȳ, h; x l ), W(u + r , u r ) = 4iw + even (x r , h; y), W(u + l , u − l ) = 4iw + even (ȳ, h; y), W(u + r , u − r ) = 4iw + even (ȳ, h; y). In particular, we can easily find curves such that each amplitude function appearing in these expressions has an asymptotic expansion described by Remark 2.1. Indeed, this just requires being able to connect the relevant points (e.g., x l andȳ in w + even (ȳ, h; x l )) through curves along which Re z(x) is increasing, which is clearly possible in view of the discussion connected to Figure 4 (see the figure for comparison). Hence, = w + even (x r , h; y) w + even (ȳ, h; y) have the stated asymptotic properties as h → 0 by Remark 2.1.
Since Re z(x) → ±∞ as x → ±∞, the asymptotic properties of R • follow from Lemma 3.2 and Remark 2.1.
These intermediate WKB solutions will allow us to prove the following quantization condition.
To see that r * = r, we use (3.16) and a logarithmic identity and get which completes the proof.
Proof of Theorem 1.2. Let r be given by Theorem 3.4, and letr be defined by the logarithm on the right of (3. 16), so that r =r + O(e −A/h ) for any A > 0. Note that the amplitude functions w + even are so-called analytic symbols with respect to the spectral parameter λ = iµ and h > 0. This means that ∂w + even (µ)/∂µ = O(h) uniformly for µ ∈ B ε (µ 0 ), see [13] or [29]. Using the definition ofr together with (3.14) it is then easy to see that h∂r(µ, h)/∂µ = O(h).
Then I * = I so we may apply Lemma 2.7 (with a in the lemma given by a(µ, h) = −hr(µ, h)) to conclude that if h is sufficiently small then there is precisely one µ k which solves I(µ, h) = (k + 1 2 )πh. Moreover, µ k ∈ R. By Theorem 3.4 this means that eigenvalues λ = iµ of P (h) are purely imaginary for µ near µ 0 .
Remark 3.5. By Lemma 2.7 (with a(µ, h) ≡ 0) there is precisely one solution µ sl k to I(µ) = (k + 1 2 )πh, and µ sl k ∈ R. From the previous proof we then infer that |µ k − µ sl k | = O(h 2 ) by the aid of Taylor's formula, where λ k = iµ k is the eigenvalue of P (h) satisfying (3.15). Moreover, similar arguments also show that where C is an upper bound of ∂I(µ)/∂µ for µ ∈ B ε (µ 0 ). Hence, if λ j = iµ j is an eigenvalue such that µ j solves (3.15) with k replaced by j = k, then showing that there is a unique eigenvalue O(h 2 )-close to µ sl k .
Remark. As shown by Theorem 1.2, eigenvalues of P (h) are purely imaginary for single-lobe potentials. In particular, the Stokes geometry depicted in the right panel of Figure 4 is never realized in the occasion of an eigenvalue. Heuristically this can be explained by the fact that there would otherwise be a curve transversal to the Stokes lines which connects the Stokes sector to the left of α l with the sector to the right of α r . Hence, the exact WKB solution u l above, which can be written as u l (x) = e z(x)ũ for someũ, could be continued into this right sector along a curve where Re z(x) is increasing. Letting x → ∞ along R would yield a contradiction to the fact that u l is collinear (modulo an exponentially small error) with the function u 0 ∈ L 2 (R + ).

Eigenvalues for a double-lobe potential
In this section we suppose that V is a double-lobe potential near µ 0 and consider eigenvalues λ = iµ with µ ∈ B ε (µ 0 ), where B ε (µ 0 ) is a small neighborhood as described in connection with Definition 1.3. The goal is to derive a quantization condition for such eigenvalues, which will then be used to prove the eigenvalue splitting occurring for symmetric potentials described in the introduction. For this reason, we will repeatedly include additional statements resulting from imposing the assumption that The Stokes geometry has been described in §2.3 and Figure 5. As in Section 3 we choose real numbers x l ,x l andx r , x r such that x l <x l < Re α l and Re α r < x r < x r . We also choose points y l and y r in the upper half-plane such that Re α l < Re y l < Re β l , Re β r < Re y r < Re α r , see Figure 8. All points are chosen independent of µ ∈ B ε (µ 0 ). In the case when V (x) = ±V (−x) for x ∈ R we choose y • and x • ,x • so that y l = −ȳ r , and x l = −x r ,x l = −x r .
Let u 0 ∈ L 2 (R + ) and v 0 ∈ L 2 (R − ) be the functions described in Section 3 such that λ = iµ is an eigenvalue of P (h) if and only if u 0 = cv 0 for some constant c = c(µ, h). We choose u 0 and v 0 in accordance with Remark 3.1 so that v * 0 = iv 0 , u * 0 = iu 0 and c = c * . When V is even we can choose v 0 as above and define Then u * 0 = iu 0 and, using the fact that v 0 solves (1.1), it is easy to check that u 0 also solves (1.1). When V is odd we instead define u 0 as Note that these are defined as in Section 3 except that y has now been replaced by y • , compare with (3.7). The reason for this is of course that we now have two lobes instead of one. Represent v 0 and u 0 as the linear combinations v 0 = c 11 u + l + c 12 u − l , u 0 = c 21 u + r + c 22 u − r . As for single-lobes one checks that c 12 = −ic * 11 using Proposition 2.4 near the left lobe, see (3.9). When V (β l ) = ±V (β r ) we use Proposition 2.4 near the right lobe and obtain for some symbols c l and c r .
by Proposition 2.5. In view of the definition (4.1) of u 0 for even V we get Comparing the right-hand side with (4.4) we see that c l = c r since V (β l ) = V (β r ) when V is even. If V is odd then by Proposition 2.5. In view of the definition (4.2) of u 0 for odd V we get We again see that c l = c r in view of (4.4) since V (β l ) = −V (β r ) when V is odd.
Let u l , u l and u r , u r be given by (3.2)-(3.3). Then Lemma 3.2 holds also for double-lobe potentials, and by replacing y with y • in the proof of Lemma 3.3 we find that the symbols c l , c r in (4.3)-(4.4) can also be written as in (3.12)-(3.13), that is, When V is symmetric we have u l (−x) = u r (x) and u l (−x) = u r (x), so using the representations (3.4)-(3.5) one checks that m − = l − (while m + = ±l + when V (x) = ±V (−x)) as in Remark 4.1. Since c l = c r it follows that τ l = τ r .

Now write
Recall that λ = iµ is an eigenvalue precisely when det(v 0 u 0 ) = 0. Since v l and v r are linearly independent by Lemma 4.2, a straightforward computation using (4.9)-(4.10) shows that det(v 0 u 0 ) = 0 is equivalent to i.e., We rewrite this as 0 = e iI l /h d 12 + e −iI l /h d 11 e iIr/h d 22 + e −iIr/h d 21 (4.14) − e −2J/h sin(I l /h) sin(I r /h). Lemma 4.3. Let µ ∈ B ε (µ 0 ). In the case when V (β l ) = ±V (β r ) we have Proof. We first note that since v * l = iv l , u * 1 = iu 2 we have d * 11 = d 12 in view of (4.12). Similarly, d * 21 = d 22 by (4.13). Using the arguments in Remark 4.1 it is easy to check that in the symmetric case V (x) = ±V (−x) we have d 12 = ±d 22 .
Proof. We start with the proof of (i) and note that the phase base points of u + (x; β l , y l ) and u − (x; β r ,ȳ r ) differ. We therefore rewrite u + (x; β l , y l ) as (A.1) u + (x; β l , y l ) = e J/h u + (x; β r , y l ), where we have used the fact that This identity is straightforward to check; in fact it can be established using the proof of [12, Lemma 5.5] with obvious modifications. Since we can find a curve from y l toȳ r along which Re z(x) is strictly increasing we can evaluate the Wronskian at y r (see (2.7)) which gives (i), with w + even (ȳ r ; y l ) = 1 + O(h) by Remark 2.1. We now prove (ii). By Lemma 2.3 we have u + (x; β r , y r ) = iu − (x; β r ,y r ) foř x neary r in case 1 • and u + (x; β r , y r ) = −iu − (x; β r ,ŷ r ) forx nearŷ r in case 2 • . In each case, take the function on the right and continue it through the branch cut starting at β r into the domain in the usual sheet containing y l . Note that at y l these functions take the values iu − (y l ; β r ,y r ) and −iu − (y l ; β r ,ŷ r ), respectively. Using (A.1) and evaluating the Wronskian atŷ r (see (2.7)) gives W(u + (x; β l , y l ), u + (x; β r , y r )) = 4ie J/h iw + even (y r ; y l ) in case 1 • , −iw + even (ŷ r ; y l ) in case 2 • .
In view of (i) this proves (iv) since J * = J. Since W(u 3 , u 2 ) = (W(u 1 , u 4 )) * anď y =ŷ, it is easy to check that (iii) follows from (ii) in a similar manner.