Exponential asymptotics for the eigenvalues in the broken PT-symmetric region

TL;DR

This paper applies exponential asymptotics to accurately predict eigenvalues in PT-symmetric quantum models, especially in the broken region where traditional WKB methods fail, extending understanding of non-Hermitian Hamiltonians.

Contribution

It introduces exponential asymptotics as a new approach to predict eigenvalues in broken PT-symmetric regions, surpassing the limitations of traditional WKB methods.

Findings

Excellent agreement with numerical eigenvalues across the parameter range

Method extends to a wider class of PT-symmetric problems

Provides insight into eigenvalue behavior in broken symmetry region

Abstract

Stemming from the seminal work of Bender & Boettcher in 1998 (Phys. Rev. Lett. vol. 80 pp. 5243-5246), there has been great interest in the study of PT-symmetric models of quantum mechanics, where the primary focus is with the study of non-Hermitian Hamiltonians that nevertheless produce countably infinite sets of real-valued eigenvalues. One of the fundamental models of such a system is governed by the Hamiltonian $H = \hat{p}^2 + x^2(ix)^{\varepsilon}$. In their work, Bender & Boettcher proposed a WKB methodology for the prediction of the discrete eigenvalues in the so-called unbroken region of $\varepsilon > 0$. However, the authors noted that this methodology fails to predict those 'broken' eigenvalues for $\varepsilon < 0$. Here, we shall explain why the traditional WKB methodology fails, and we shall demonstrate how eigenvalues for all relevant values of $\varepsilon$ can be predicted using techniques in exponential asymptotics. These predictions provide excellent agreement to exact numerical results over nearly the entire range of values. Moreover, such techniques can be extended to a much wider range PT-symmetric problems.

Equations

$$-\frac{d{}^2\psi }{d{x}^2}+x^2(\mathrm{i}x{)}^{\epsilon }\psi =E\psi ,$$

$$p=2+\epsilon ,$$

$$E=E_n={\left[\frac{\sqrt{\pi }\left(n+\frac{1}{2}\right)\Gamma \left(\frac{3}{2}+\frac{1}{p}\right)}{\Gamma \left(1+\frac{1}{p}\right)\sin \left(\frac{\pi }{p}\right)}\right]}^{\frac{2p}{p+2}}\text{for n→∞}.$$

$$ϵ=E^{-\frac{p+2}{2p}},$$

$$2\mathrm{i}\exp\left[\frac{2R(p)\cos\bigl{(}\frac{\pi}{p}\bigr{)}}{\epsilon}\right]\cos\left[\frac{2R(p)\sin\bigl{(}\frac{\pi}{p}\bigr{)}}{\epsilon}\right]=0,$$

$$R(p)=\frac{\sqrt{\pi }\Gamma \left(1+\frac{1}{p}\right)}{2\Gamma \left(\frac{3}{2}+\frac{1}{p}\right)}.$$

$$2\mathrm{i}\exp\left[\frac{2R(p)\cos\bigl{(}\frac{\pi}{p}\bigr{)}}{\epsilon}\right]\cos\left[\frac{2R(p)\sin\bigl{(}\frac{\pi}{p}\bigr{)}}{\epsilon}\right]-\frac{2\pi\mathrm{i}\epsilon^{p}}{2^{p+2}\Gamma(-p)}=0.$$

$$x=E^{1/p}z\text{and}ϵ=E^{-\frac{p+2}{2p}},$$

$$-{ϵ}^2{f}^{\prime \prime }(z)-(\mathrm{i}z{)}^{p}f(z)=f(z),$$

$$\text{with f→0 as z→∞eiπ[−1/2∓2/(p+2)]},$$

$${\theta }_{\text{left/right}}=\pi \left(-\frac{1}{2}\mp \frac{2}{p+2}\right),$$

$$f(z)\sim e^{\mathrm{i}\varphi (z)/ϵ}\sum _{n=0}^{\infty }{ϵ}^n{A}_n(z).$$

$$({\varphi }^{\prime }{)}^2-(\mathrm{i}z{)}^p$$

$$2\mathrm{i}{\varphi }^{\prime }{A}_0^{\prime }+\mathrm{i}{\varphi }^{\prime \prime }{A}_0$$

$$\varphi (z)$$

$$A_0(z)$$

$$1+(\mathrm{i}z{)}^p=0.$$

$$z_A=-\mathrm{i}{\mathrm{e}}^{-\mathrm{i}\pi /p}\text{and}{z}_B=-\mathrm{i}{\mathrm{e}}^{\mathrm{i}\pi /p}.$$

$$\begin{array}{c}{\Phi }_{\text{I}}(z)={\int }_{z_A}^z[1+(\mathrm{i}t{)}^p{]}^{1/2}dt,{\Phi }_{\text{III}}(z)={\int }_{z_B}^z[1+(\mathrm{i}t{)}^p{]}^{1/2}dt,\\ \Phi (z)={\int }_a^z[1+(\mathrm{i}t{)}^p{]}^{1/2}dt,\end{array}$$

$$f(z)\sim\begin{cases}\displaystyle\biggl{[}\frac{a_{1}}{(\phi^{\prime})^{1/2}}\biggr{]}\mathrm{e}^{\mathrm{i}\Phi_{\textrm{I}}/\epsilon}+\biggl{[}\frac{b_{1}}{(\phi^{\prime})^{1/2}}\biggr{]}\mathrm{e}^{-\mathrm{i}\Phi_{\textrm{I}}/\epsilon}&\text{for $z$ in Region I},\\[11.38092pt] \displaystyle\biggl{[}\frac{a_{2}}{(\phi^{\prime})^{1/2}}\biggr{]}\mathrm{e}^{-\mathrm{i}\Phi/\epsilon}+\biggl{[}\frac{b_{2}}{(\phi^{\prime})^{1/2}}\biggr{]}\mathrm{e}^{\mathrm{i}\Phi/\epsilon}&\text{for $z$ in Region II},\\[11.38092pt] \displaystyle\biggl{[}\frac{a_{3}}{(\phi^{\prime})^{1/2}}\biggr{]}\mathrm{e}^{-\mathrm{i}\Phi_{\textrm{III}}/\epsilon}+\biggl{[}\frac{b_{3}}{(\phi^{\prime})^{1/2}}\biggr{]}\mathrm{e}^{\mathrm{i}\Phi_{\textrm{III}}/\epsilon}&\text{for $z$ in Region III}.\end{cases}$$

$$a_2{\mathrm{e}}^{\mathrm{i}\Phi (B)/ϵ}-\mathrm{i}{b}_2{\mathrm{e}}^{-\mathrm{i}\Phi (B)/ϵ}$$

$$a_2{\mathrm{e}}^{\mathrm{i}\Phi (A)/ϵ}+\mathrm{i}{b}_2{\mathrm{e}}^{-\mathrm{i}\Phi (A)/ϵ}$$

$$\left|\begin{array}[]{cc}\mathrm{e}^{\mathrm{i}\Phi(B)/\epsilon}&-\mathrm{i}\mathrm{e}^{-\mathrm{i}\Phi(B)/\epsilon}\\ \mathrm{e}^{\mathrm{i}\Phi(A)/\epsilon}&\mathrm{i}\mathrm{e}^{-\mathrm{i}\Phi(A)/\epsilon}\end{array}\right|=\mathrm{i}\mathrm{e}^{\mathrm{i}[\Phi(B)-\Phi(A)]/\epsilon}+\mathrm{i}\mathrm{e}^{-\mathrm{i}[\Phi(B)-\Phi(A)]/\epsilon}=0.$$

$$\mathrm{i}{\mathrm{e}}^{2\mathrm{i}\Phi (A)/ϵ}+\mathrm{i}{\mathrm{e}}^{2\mathrm{i}\Phi (B)/ϵ}=0⟹{\mathrm{e}}^{2\mathrm{i}[\Phi (B)-\Phi (A)]/ϵ+\mathrm{i}\pi }=1,$$

$$\frac{2}{\epsilon}\Bigl{[}\phi(B)-\phi(A)\Bigr{]}=\frac{2}{\epsilon}\int_{z_{A}}^{z_{B}}(1+(\mathrm{i}t)^{p})^{1/2}\,\operatorname{d\!}{}{t}=(2n+1)\pi\qquad\text{for $n\in\mathbb{Z}$}.$$

$$\varphi (z)$$

$$A_0(z)$$

$$f(z)\sim \frac{1}{[1+(\mathrm{i}z{)}^p{]}^{1/4}}\exp \left[\frac{1}{ϵ}{\int }_a^z[1+(\mathrm{i}t{)}^p{]}^{1/2}dt\right].$$

$$ϵA^{\prime \prime }+2\mathrm{i}{\varphi }^{\prime }{A}^{\prime }+\mathrm{i}{\varphi }^{\prime \prime }A=0.$$

$$A(z)\sim \sum _{n=0}^{\infty }{ϵ}^n{A}_n(z).$$

$$A_{n-1}^{\prime \prime }+2\mathrm{i}{\varphi }^{\prime }{A}_n^{\prime }+\mathrm{i}{A}_n{\varphi }^{\prime \prime }=0.$$

$$A_n(z)\sim \frac{B(z)\Gamma (n+\gamma )}{[\chi (z){]}^{n+\gamma }}\text{as n→∞}$$

$$-2\mathrm{i}{\chi }^{\prime }{\varphi }^{\prime }+({\chi }^{\prime }{)}^2=0,$$

Full text

1 Introduction

In classical or quantum mechanics, the equations governing the time-evolution of a system can be derived from a Hamiltonian, $H$. It is a standard axiom in quantum mechanics that $H$ must be Hermitian, and thus its eigenvalues real. This is in connection with the assumption that measurements of the system must correspond to eigenvalues of $H$, and hence the Hermitian property guarantees real-valued outcomes from the model.

The question, however, is whether the Hermitian property of $H$ is a necessity or whether there exist more fundamental restrictions on the Hamiltonian that would nevertheless result in real-valued eigenvalues. Largely beginning with the work of Bender & Boettcher [4] in 1998 (though similar ideas proliferated the literature going back to Dyson [11]), there has been a great deal of interest in studying so-called $\mathcal{PT}$-symmetric theories, which generalise the Hermitian property of Hamiltonians. Such theories posit that as a fundamental assumption, $H$ should be invariant under a parity transformation, $\mathcal{P}$, and time reversal, $\mathcal{T}$.

The canonical model of a non-Hermitian $\mathcal{PT}$-symmetric Hamiltonian is given by $H=\hat{p}^{2}+x^{2}(\mathrm{i}x)^{\varepsilon}$ and corresponds to the time-independent Schrödinger eigenvalue problem for $\psi=\psi(x)$, given by

[TABLE]

where $E$ is the eigenvalue. The boundary conditions for the above problem require that $\psi\to 0$ exponentially rapidly as $|x|\to\infty$ in a pair of adjacent wedges in the complex $x$-plane [see later (9)]. The $\mathcal{PT}$-symmetric nature of (1) is equivalent to the fact that the equation is invariant under the transformation $x\mapsto-x$ and $\mathrm{i}\mapsto-\mathrm{i}$.

As explained in e.g. [2], the remarkable significance of (1) is that for $\varepsilon>0$, its eigenvalues, $E$, are real, despite the fact that the Hamiltonian is non-Hermitian. In particular, the following is now known about the spectrum of $H$, shown in Fig. 1.

For $\varepsilon>0$, the eigenvalues, $E=E_{n}(\varepsilon)$, are real and form a discrete countably infinite set. These eigenvalues match with the case of the harmonic oscillator, $E_{n}=2n+1$ at $\varepsilon=0$. An asymptotic approximation of $E_{n}$ in the limit $n\to\infty$ was developed by Bender & Boettcher [4], and the reality of the spectrum was proved by Dorey et al. [10]. 2. 2.

It is also known, primarily through numerical solutions of the eigenvalue problem (1) that as $\varepsilon$ decreases below zero the eigenvalues move into the complex plane, forming complex-conjugate pairs. The ‘fingers’ in the bifurcation diagram begin to close off. 3. 3.

There are further interesting behaviours in regards to the complex-valued eigenvalues and eigenfunctions for $\varepsilon<0$. For example, there is an infinite-order exceptional point at $\varepsilon=-1$ where $|\operatorname{Re}E|\to\infty$ and $|\operatorname{Im}E|\to 0$ in the form of a logarithmic spiral.

Our interest relates to the so-called broken region $\varepsilon<0$, and we highlight two open questions of significance. First, is there a simple explanation of why the spectrum must change at $\varepsilon=0$ and why, for $\varepsilon<0$, the ‘fingers’ of the bifurcation diagram terminate? Second, is there an asymptotic approximation of the eigenvalues that remains valid in the broken region?

1.1 Key idea of why the eigenvalues terminate

We find it convenient to set

[TABLE]

so that the exceptional point of $\varepsilon=0$ corresponds to $p=2$.

In Bender & Boettcher [4] it was shown that in the large eigenvalue limit, $|E|\to\infty$, the eigenvalues are given by [cf. their eqn (5)]

[TABLE]

Bender & Boettcher’s derivation of (3) relies upon the development of distinct WKB approximations to (1), each of which is valid as $|x|\to\infty$ in separate sectors of the complex plane. Matching the WKB solutions from one sector to the next involves a solvability condition that produces the eigenvalue conditions. If we introduce the small parameter

[TABLE]

so that $\epsilon\to 0$ as $E\to\infty$, then in our notation, the solvability condition is

[TABLE]

where we have defined

[TABLE]

Setting the argument of the cosine term in (5) to $(n+1/2)\pi$ for $n\in\mathbb{Z}$ results in the Bender & Boettcher asymptotic result of (3). The WKB analysis that produces (5) mirrors the well-known semi-classical approach [12], with the distinction of requiring integration of the differential equation in the complex plane. However, Bender & Boettcher [4] note that for $p<2$, the path of continuation along which the WKB solutions are typically matched must proceed through a branch cut in the complex plane; hence apparently condition (5) is no longer valid. We clarify this remark in §3.

In this work, we shall demonstrate that for $p<2$ (or $\varepsilon<0)$, the eigenvalue condition (5) should instead include an additional term,

[TABLE]

When $p>2$, then $R(p)\cos(\pi/p)>0$ and the contribution from the second term of (7) is exponentially small in comparison to the first; thus (3) remain valid. However, for $1<p<2$ the dominance exchanges with $R(p)\cos(\pi/p)<0$, and the second term now exponentially dominates. As $\epsilon\to 0$, there is no way to zero this second term, and hence in the large eigenvalue limit, there are no real eigenvalues that satisfy (7).

Consequently the new solvability condition (7) provides a formal argument of why the discrete family of eigenvalues must terminate for $p<2$, with the ‘fingers’ in the bifurcation diagram of Figure 1 closing firstly for the largest eigenvalues. Moreover, note that for finite values of $\epsilon$, it is still possible to satisfy (7). The solution of this transcendental equation then allows for an asymptotic prediction of the eigenvalues. The agreement between asymptotically calculated eigenvalues and exact numerical eigenvalues is excellent, even at moderate values of the eigenvalue (see the later Figure 5).

The rest of this paper will be devoted to presenting the methodology that leads to the above results; in particular, we wish to emphasise that this methodology, which uses exponential asymptotics, is considerably general, and can likely be applied to the range of linear eigenvalue problems of interest to studies in $\mathcal{PT}$-symmetry. In fact, as we discuss in Sec. 6, similar methodologies are applied to much more difficult problems in nonlinear differential equations where standard WKB approaches are not valid.

2 Mathematical formulation

Recall our choice of $p=\varepsilon+2$, and in anticipation of studying the large eigenvalue limit of (1), with $|E|\to\infty$, we re-scale the independent variable in (1) by setting

[TABLE]

The $\mathcal{PT}$-symmetric eigenvalue problem for $\psi(x)=f(z)$ where $z\in\mathbb{C}$ is now given by

[TABLE]

where we are primarily interested in values of real values of $p\geq 1$. Above and in the remainder of this paper, primes (′) denote differentiation with respect to $z$.

Let us explain the boundary conditions. The path on which we solve the above equation is shown in Fig. 2 and corresponds to the boundary conditions (9b). Note that for large $|z|$, the dominant balance of (9a) involves $\epsilon^{2}f^{\prime\prime}\sim-(\mathrm{i}z)^{p}f$. Consequently a WKB solution can be developed in the limit $|z|\to\infty$ and contains the exponential factor $\exp[\pm(\mathrm{i}z)^{\alpha}/(\alpha\epsilon)]$ where $\alpha=p/2+1$. When $\alpha$ is non-integral, consider the positive branch of $(\mathrm{i}z)^{\alpha}$ with the branch cut taken vertically upwards. With $z=|z|\mathrm{e}^{\mathrm{i}\theta}$ the exponential corresponding to this branch decays within two wedges (known as the Stokes wedges) centred at angles

[TABLE]

as shown in Fig. 2. The size of each wedge is $2\pi/(p+2)$. The exponential argument that appears in the boundary condition (9b) corresponds to the centre of the wedge, where the WKB solution is decaying most rapidly at infinity.

Note that a non-trivial solution to (9a) can be obtained by enforcing the condition that $f$ decays in any two non-adjacent Stokes sectors. The two particular sectors chosen in (9) are such that the solution of the system is an analytic continuation of the case of the harmonic oscillator. For $p=2$ there are two wedges of size $\pi/2$ and the solution is required to tend to zero as $|z|\to\infty$ along the real axis. For values of $p$ larger-than $p=2$, the two wedges move into the lower-half plane, and coalesce along the negative imaginary axis when $p\to\infty$. For values of $p$ smaller than $p=2$, the wedges move into the upper half-plane.

3 Failure of the traditional WKB approach

We first give an explanation of the traditional WKB approach, as applied in e.g. Bender & Boettcher [4], and explain why this approach apparently fails for the case of the broken eigenvalue region $p<2$.

In the limit $\epsilon\to 0$, we approximate the solution using the WKB ansatz,

[TABLE]

Substitution into (9) gives an eikonal equation for $\phi$ at leading order and an amplitude equation for $A_{0}$ at next order,

[TABLE]

These are solved to give

[TABLE]

From (13), the WKB approximation fails at those turning points where

[TABLE]

Of the turning points in (14) only two lie in the appropriate sectors of the complex plane and are relevant for the analysis:

[TABLE]

The fact that the other turning points are not involved at leading-order will become clear in the exponential asymptotics methodology presented later.

We may write down a composite WKB approximation. From (13a), we choose the positive sign of $\phi$ and define

[TABLE]

where $a$ can be taken to be any point where the integral is defined (e.g. $a=0$). The principal branches are chosen throughout. Then we form the following composite solution:

[TABLE]

The three regions are shown in Figure 3 for the case of $p=3$. From each of the two turning points, marked A and B, there are three solid lines, indicating curves where $\operatorname{Im}(\mathrm{i}\Phi_{\textrm{I}})=0$ and $\operatorname{Im}(\mathrm{i}\Phi_{\textrm{III}})=0$, respectively. Along PA, $(\mathrm{i}\Phi_{\textrm{I}})$ is purely real and negative and similarly along BQ, $(\mathrm{i}\Phi_{\textrm{III}})$ is real and negative. Thus, examining the composite solution (17), in order for $f$ to decay as $|z|\to\infty$ along PA and BQ, we require $b_{1}=a_{3}=0$.

This leaves four unknown constants. A local analysis near the turning points, $z=z_{A}$, lets us relate solutions in Region I with II and hence $\{a_{1},\,b_{2}\}$ to $\{a_{2},\,b_{2}\}$. Similarly solutions in Region II and III are matched from a local analysis near $z=z_{B}$, and this relates $\{a_{2},\,b_{2}\}$ to $\{a_{3},b_{3}\}$. The details of this matching procedure are given in Appendix A. This local asymptotic analysis re-scales the independent variable, $z$, such that locally near $z=z_{B}$ say, the solution, $f$, behaves as an Airy function to leading order. The asymptotic behaviour of the (real-valued) Airy function at positive and negative infinity then allows the solutions in Regions II and III to be matched. Crucially, this matching occurs along the segment AB show in Figure 3 where $\operatorname{Re}(\mathrm{i}\Phi_{\textrm{I}})=\operatorname{Re}(\mathrm{i}\Phi_{\textrm{III}})$.

The result is a pair of homogeneous linear equations

[TABLE]

for the two unknowns $a_{2}$ and $b_{2}$, and where we have used the shorthand $\Phi(A)=\Phi(z_{A})$ and $\Phi(B)=\Phi(z_{B})$. Thus, in order for a nonzero solution to exist we need

[TABLE]

This gives the eigenvalue condition

[TABLE]

or in terms of $\phi$ in (13a),

[TABLE]

Evaluating the integral explicitly, and remembering that $\epsilon=E^{1/2+1/p}$ from (8), returns the Bender & Boettcher [4] result of (3).

As we had noted earlier in the introduction, the eigenvalues given by this formula are illustrated in Figure 1, along with a numerical calculation of the exact eigenvalues. We see that, where the eigenvalues exist, the approximation (3) is good. However, it fails to find that there are only a finite number of eigenvalues for some values of $p$.

The WKB approximation fails due to the fact that when $p$ is not an integer there is a branch point at the origin, with a branch cut up the positive imaginary axis, say. For $p>2$ the path between $z_{A}$ and $z_{B}$ on which $\phi$ is real misses this branch cut (it passes below the origin, as shown in Figure 3, and the WKB analysis gives a good approximation to the eigenvalues. However, when $p<2$ the path between $z_{A}$ and $z_{B}$ on which $\phi$ is real passes through the branch cut, as shown in Figure 4. Thus the turning points which have been connected via the WKB analysis (and thus the points at infinity at which the boundary conditions are imposed) lie on different Riemann sheets of the solution.

In principle, this traditional WKB approach can be modified by taking the first WKB approximation, analytically continuing it around the branch point at $z=0$, and then matching with the second WKB approximation. However, to do so would already require a knowledge of Stokes switchings in the WKB approximation, so that it is more straightforward to simply keep track of the Stokes lines on the original Riemann sheet.

4 An alternative approach using Stokes lines

In the traditional WKB approach, solutions are developed on different subregions of the plane, and then matched together. In this case, the eigenvalues $\epsilon=\epsilon_{n}$ emerge as a result of solvability conditions on the constants of integration. However, we saw that this approach does not allow a prediction of the eigenvalues in the region $p<2$. In this section, we present an approach that uses exponential asymptotics. In this approach, the discrete spectrum (3) arises in connection with the Stokes phenomenon, and we shall find an additional contribution that causes the bifurcation curves to close.

We return to the WKB approximation in (11) and this time, we choose a single exponential argument from (13a), with

[TABLE]

where $a$ can be freely chosen. For simplicity, the constant in (13b) is taken to be unity. Thus we consider the leading-order approximation

[TABLE]

Above in (22a) $\phi$ has been defined with the positive square-root branch in mind (say with both branch cuts from $z=z_{A},z_{B}$ taken upwards). We may verify that as $|z|\to\infty$ along the centre of the Stokes wedges, the WKB ansatz (23) decays with $f\to 0$. As we know, the WKB solution is singular at the two turning points, and this necessitated the matched procedure presented in the last section.

However, it is still possible to travel from one wedge to the other without passing through the two singularities. If $z$ remains below the curve PABQ in Figure 3, then there is no apparent indication that (23) should fail to approximate the true solution, $f$, over the desired path and indeed, (23) seems to be valid whatever the value of $\epsilon$. One arrives at the erroneous conclusion that there is no restriction on $\epsilon$ and a continuous spectrum results.

The apparent paradox is resolved by taking into account the extra terms that are switched-on due to the Stokes Phenomenon. In what follows, we shall see that as (23) is analytically continued from $-\infty$ to $+\infty$, two subdominant exponentials are switched-on—one due to a Stokes line from $z=-1$ and another due to a Stokes line from $z=1$. The subsequent restriction of these exponentials in order to satisfy the boundary conditions is what produces the discrete set of eigenvalues.

4.1 Exponential asymptotics

We require exponential asymptotics to determine those exponentially small terms and the conditions for their switching, and this can be done using a procedure of optimal truncation and Stokes-line smoothing [7]. For convenience, let us write $f(z)=\mathrm{e}^{\mathrm{i}\phi/\epsilon}A(z)$ so that the differential equation (9a) is now

[TABLE]

We then proceed further by expanding fully

[TABLE]

At $\mathcal{O}(1)$, $A_{0}$ is given by (22b) as we know. At $\mathcal{O}(\epsilon^{n})$ in (24), the equation for $A_{n}$ is given by

[TABLE]

Notice that $A_{0}$ in (22b) has singularities at the turning points, $1+(\mathrm{i}z)^{p}=0$, as well as at $z=0$ if $p$ is not an integer. However, solving for $A_{n}$ involves differentiation of the previous order, $A_{n-1}$. Thus generically, the power of each singularity must grow at each subsequent order and consequently $A_{n}$ diverges in the form of a factorial over a power as $n\to\infty$. Such divergence is typical of singular perturbation problems [7, 9]. We assume the late terms diverge in the form of

[TABLE]

with $\chi=0$ at the singularities of the early terms. Note that a separate factorial-over-power ansatz of the form (27) is required for each singularity in the late terms, but only those singularities that are associated with active Stokes lines need to be accounted for. Substituting (27) into (26) gives, for the first two orders in the limit $n\to\infty$,

[TABLE]

The first equation in (28a) yields $\chi(z)$, a quantity known as the singulant [9]:

[TABLE]

and by assumption, $\chi=0$ at those singularities $z=z_{*}$ that cause $A_{n}$ to diverge. By linearity of the asymptotic procedure no new singularities are introduced apart from those in the early terms. Hence $z_{*}$ must correspond either to the turning points given by (14) or, for non-integral powers $p$, the branch point at $z=0$.

At next order (28b) is solved, giving

[TABLE]

Hence the late terms in the WKB expansion of $f$ are

[TABLE]

As noted by e.g. Dingle [9], this divergence of the expansion is associated with subdominant exponentials that are switched-on via the Stokes Phenomenon.

4.2 Inner-matching procedure

As it will be shown from analysis of the Stokes lines in §4.3, only three of the singularities are important for specification of the eigenvalues: the two former points introduced as $z_{A}$ and $z_{B}$ in (15) and the branch point $z=0$. There are thus three factorial/power ansatzes of the form (31) to consider, and values of $\gamma$ and $\Lambda$ must be determined for each case.

First, the value of $\gamma$ can be found by ensuring that the order of the singularity in $A_{n}$ (27) is consistent with that of $A_{0}$ in (22b). Next, the value of $\Lambda$ can be found by developing an inner solution valid near the singularities and matching with the WKB solution. The matching procedure detailed in Appendix B shows how these are derived. In summary, they are given by

[TABLE]

We have now completely determined the behaviour of the late terms, of which the relevant divergence is driven by a sum of three factorial/power ansatzes of the form (27) or alternatively (31). The values of $\chi(z)$, $B(z)$, $\gamma$, and $\Lambda$ are given by (29), (30), and (32).

4.3 Exponential switchings

The connection between the late terms (31) and the exponentials switched-on from the Stokes Phenomenon can be understood based on a procedure of optimal truncation and Stokes-line smoothing [6, 7, 9]. Briefly, the idea is as follows. First, the solution to (24) is expressed as a truncated expansion plus a remainder term,

[TABLE]

so that the equation for the remainder satisfies

[TABLE]

In the limit $\epsilon\to 0$, the optimal truncation point, $N\to\infty$. Thus the behaviour of the remainder is predicated by the behaviour of the late terms, thus establishing a connection between the exponentials switched on and the divergence of the series. If the series is optimally truncated, then $R_{N}$ is exponentially small. Moreover, the exponentially small remainder is switched-on when the solution is analytically continued across critical curves in the complex plane known as Stokes lines in a process known as the Stokes Phenomenon. Analysis of (34) indicates two crucial facts that have been demonstrated by e.g. [6, 7, 9].

First, there are Stokes lines where successive late terms of (25) have the same phase, i.e. those points $z\in\mathbb{C}$ where $u$ is positive and real. Thus,

[TABLE]

where $u$ is given by (29).

Second, the exponentials that are switched-on are given by

[TABLE]

and a derivation of this result is given in Appendix C.

Thus, our goal is to consider the paths of continuation shown in Figure 2 and determine which Stokes lines from which singularities are crossed. As the base asymptotic series (25) is analytically continued across each Stokes line, it must switch-on the exponentially small term (36) where $\chi$, $B$, and $\gamma$ correspond to the singularity from which the Stokes line originates.

As shown in Figure 3, for general non-integral values of $p$, there are three sets of relevant Stokes lines. Moving through the three Regions I, II, and III thus switches on the three contributions,

[TABLE]

where we have substituted values of $\gamma$ and $\Lambda$ from (32). Notice that the exponential arguments are proportional to $\mathrm{e}^{-2\mathrm{i}\phi}$. Following the discussion of (17), we note that $\operatorname{Re}(-2\mathrm{i}\phi)$ is real and negative as $|z|$ tends to infinity along the (Anti-Stokes) line marked BQ in Figures 3 and 4. Thus, although the terms in (37) exponentially small at the point they are switched on, they become exponentially large as $|z|\to\infty$ in Region III. For the solution to decay it is necessary for (37) to be zero. Using $\chi$ from (32), we have

[TABLE]

Note that by direct integration,

[TABLE]

where $R$ is defined in (6), and moreover $\phi(z_{B})=-\overline{\phi(z_{A})}$. Combination of the two first terms in (38) then yields the same equation as first stated in (7), or

[TABLE]

The eigenvalue condition (7) is the main result.

Notice that for $p>2$, $R(p)\cos(\pi/p)>0$ and the first term exponentially dominants the second. Setting the cosine term to zero hence yields the Bender & Boettcher [4] result of (3). In other words, this is equivalent to the condition that the exponential switched-on across the Stokes line from $z_{A}$ is switched off due to the Stokes line from $z_{B}$. The classical quantisation condition has been re-derived.

However, for for $p<2$, $R(p)\cos(\pi/p)<0$ and the exponential dominance switches, with the second term in (40) now dominant. In the limit $\epsilon\to 0$, there are no solutions; thus there are only a finite number of real eigenvalues for $p<2$. The roots of equation (38) are shown in Figure 5, along with the numerical calculation of the eigenvalues. The agreement is very good.

Extension to complex eigenvalues

As the branches of the bifurcation diagram coalesce (as in Fig. 5), the real-valued eigenvalues merge in pairs; thereafter, for smaller values of $p$, there are two associated complex-conjugate eigenvalues. For the case of the two merging braches going through $E=1$ and $E=3$ and $p=2$, this produces the curves shown in Fig. 6.

In a recent paper, Bender et al. [5] have presented numerical and asymptotic results for the eigenvalues of (1) (with $p=2+\varepsilon$) near the special values of $p=\{1,-1,-2\}$. In fact, the exponential asymptotics we have presented §4.1 can be used to reproduce many of their results in a more unified way.

By way of example, let us investigate the limit of $p\to 1^{+}$, where $|\operatorname{Im}E|\to 0$ and $|\operatorname{Re}E|\to\infty$. In this limit, the two turning points, labeled A and B, in Fig. 4 coalesce at $z=\mathrm{i}$ and consequently, the width of Region II, bounded by the two Stokes lines, must shrink. Although turning points A and B coalesce, the previous Stokes-line contributions we have derived for points A, B, and C remains valid and it may be verified that the inner-region asymptotics of Appendix B remain unaltered in this limit. Thus we may let $p=1+\delta$ for $\delta\to 0^{+}$ directly in the eigenvalue condition (40). Then $R(p)\sim 2/3$ and $1/\Gamma(-p)\sim\delta$, and

[TABLE]

Since $\epsilon\sim E^{-3/2}$ by (8), we have in terms of $E$,

[TABLE]

Fig. 7 provides verification that the exponential scaling predicted in (42) agrees with the numerical calculation of the eigenvalues in the limit $p\to 1^{+}$ and $\delta\to 0$. Similarly Fig. 8 provides verification that the lower-order algebraic dependence of $E^{3/2}$ is as expected. Note that in providing the exact pre-factor and algebraic scaling, our (42) provides a more accurate prediction of the eigenvalues than asymptotic formula in eqn (12) of [5].

5 Application to other $\mathcal{PT}$-Ssymmetric problems

The exponential asymptotic techniques we have presented provide a more powerful and general framework than the traditional WKB analysis of (3). Thus, this idea of locating singularities in the asymptotic expansions, and examining the switching-on of exponentials as Stokes lines are crossed can be used to study the broken and unbroken eigenvalues of wider range of $\mathcal{PT}$-symmetric problems than the case of (9).

As another illustrative example, let us consider the Hamiltonian $H=\hat{p}^{2}-x^{4}-\mathrm{i}\mathcal{A}x$ that corresponds to the Schrödinger equation

[TABLE]

where $\mathcal{A}$ is a real parameter. The eigenvalues are shown in Fig.9 and they illustrate the familiar unbroken ($\mathcal{A}$ small) and broken ($\mathcal{A}$ large) regions we have observed in our previous example. Interestingly, the prediction of the eigenvalues in both regions were performed by Bender et al. [3] using traditional WKB techniques. Here, we shall show how the same analysis can be done through a straightforward application of the theory we have already developed.

We shall seek to describe the bahaviour in the limit $|E|\to\infty$. Under the re-scaling $x=E^{1/4}z$, we have

[TABLE]

where $a=\mathcal{A}E^{-3/4}$ and $\epsilon=E^{-3/4}$. We set $f(z)=\mathrm{e}^{\mathrm{i}\phi/\epsilon}A(z)$ and find that (taking the positive branch)

[TABLE]

where $b$ is an arbitrary point of integration.

The equation for $A(z)$ is identical to that of (24). Thus, written in terms of $\phi$, the development of the late terms (27) and Stokes switching relation in (36) are identical. This time, the only singularities of the leading-order problem $A_{0}\sim\text{const.}/(\phi^{\prime})^{1/2}$ are those corresponding to turning points, where

[TABLE]

The four turning points are shown for the typical value of $a=1$ in Figure 10. For all real values of $a$, there are two roots along the imaginary axis, say $z_{C}$ and $z_{D}$, and two roots located at constant $\operatorname{Im}z$, say $z_{A}$ and $z_{C}$. The Stokes lines, where $\operatorname{Im}\chi=0$ and $\operatorname{Re}\chi\geq 0$, computed from (29) are plotted solid in the figure. The path of continuation considered by [3] is shown dashed. Thus the WKB solution is analytically continued across the three Stokes lines from $z_{A}$, $z_{C}$, and $z_{B}$, the switching is given by an expression analogous to (38), with

[TABLE]

It is convenient to choose the point of integration in (45) as $b=z_{C}$, and consequently $\phi(z_{B})=-\overline{\phi(z_{A})}$. Dividing by the third exponential, we now have

[TABLE]

where we have defined

[TABLE]

so that the notation is consistent to that of eqn (13) in [3]. From numerical integration, it can be verified that the function $V(a)$ in (49) begins at $V\approx 0.87402$ and decreases monotonically as $a$ increases, passing through $V=0$ at $a=a^{*}\approx 1.18384$.

Thus for $a<a^{*}$, the first term of (48) is exponentially large and the eigenvalues are predicted by those values of $\epsilon$ where $2U(a)/\epsilon\sim(2\mathbb{Z}+1)/2$. This produces the countably infinite set seen in Fig. 9. However, when $a>a^{*}$, the second term in (48) is now dominant and the bifurcation curves are expected to close-off. The fit between asymptotic and numerical results are shown in Fig. 9 and we see that even through the entire range of $E$-values, the agreement is nearly visually indistinguishable. Thus our results in this section have duplicated those of Bender et al. [3] in their study of the quartic oscillator (44), but here we have used the exponential asymptotic framework of Sec. 4 instead of the traditional WKB methodology in Sec. 3.

There is an important distinction between the main $\mathcal{PT}$-symmetric problem of this paper, given in (9) corresponding to the Hamiltonian $H_{1}=\hat{p}^{2}-(\mathrm{i}x)^{p}$, and the sub-problem of this section, given in (44), with Hamiltonian $H_{2}=\hat{p}^{2}+x^{4}+\mathrm{i}\mathcal{A}x$. The case of $H_{1}$ is unexpectedly challenging on account of the fact that the closing-off of eigenvalues is due to a singularity of a different (non-turning-point) type than those responsible for generating the main discrete set of eigenvalues. This necessitated the more powerful exponential asymptotic techniques we have proposed, and was perhaps the reason why the prediction of the eigenvalues of $H_{1}$ remained unresolved, in contrast to those of $H_{2}$.

6 Conclusion

In this paper, we have shown how the eigenvalues of the Bender & Boettcher [4] problem can be predicted in the broken region of $\varepsilon<0$ (or $p<2)$. Previous asymptotic analyses have relied on a traditional WKB framework of matching between turning points, and we have shown that this approach is inadequate. Instead, we have proposed a methodology that use exponential asymptotics to derive the eigenvalues in both broken and unbroken regions. Thus, as the relevant parameter decreases through the critical value of $\varepsilon=0$, a previously subdominant contribution switches on and dominance exchanges. The match between the numerical and asymptotic values is excellent over the entire range of parameters.

More generally, we have shown that the large eigenvalue asymptotics of linear eigenvalue problems may be understood in terms of Stokes lines in the complex plane, and that this interpretation provides a means of calculating the eigenvalues. A subdominant exponential is turned on across Stokes lines. This exponential grows at infinity, so that its coefficient must eventually be zero, which gives the eigenvalue condition. This may be thought of in terms of the contributions from two Stokes lines being exactly out of phase, and so destructively interfering. In the classical case (e.g. corresponding to the harmonic oscillator) the eigenvalues form a series for which the corresponding eigenfunctions have an integer number of oscillations between the two Stokes lines, mirroring the picture on the real axis between the turning points, but in terms exponentially small “beyond all orders”. However, for the case of general Hamiltonians, their associated singularity and Stokes-line structures may yield additional contributions that change the previous selection mechanism of eigenvalues. The Bender & Boettcher [4] is one such case where such an additional contribution causes the eigenvalues to transition from unbroken (countably infinite) to broken (finitely many). It is expected that many other problems in $\mathcal{PT}$-symmetric quantum mechanics can be studied using this framework.

Appendix A Determination of constants

We first explain how (18a) is derived by matching about the turning point near $z=z_{B}$, although much of the presentation applies to an arbitrary turning point at $z=z_{*}$. We-scale in the inner region with

[TABLE]

where $\mathrm{i}p(\mathrm{i}z_{*})^{p-1}K^{3}=-1$. The branch of $K$ is chosen so that $\phi$ in the later expression (65) corresponds to a decaying exponential, $\mathrm{e}^{\mathrm{i}\phi/\epsilon}$, as $y\to\infty$. Under this scaling,

[TABLE]

where $P$ is a constant that depends on $K$, but whose exact value we do not require.

Under the substitution (50), the main equation (9a) now yields the standard Airy equation

[TABLE]

The solution of (53) that matches with the requisite decay condition at infinity is $f=C\operatorname{Ai}(y)$ for a constant $C$. From [1], we may write the full expansion of the Airy function as

[TABLE]

valid as $y\to\infty$. This solution is then matched with the inner limit of the WKB solution $f=f_{\textrm{III}}$ in (17). In this limit,

[TABLE]

Matching (54) and (55) gives $C/(2\sqrt{\pi})=b_{3}/(P\epsilon^{1/6})$.

Next, we turn to the matching of the inner solution with the solution of Region II in (17). Firstly, from [1], we have as $y\to-\infty$,

[TABLE]

This should be matched with the inner limit of $f=f_{\textrm{II}}$ in (17) as $z\to z_{B}$ and in fact, this is the only step in the above presentation that requires the use of $z_{B}$ instead of a generic $z_{*}$ for the turning point (since the definition of the solution in region III shifts the exponential argument to be zero at $z_{B}$). We have

[TABLE]

and combining (56) with (57) gives the pair of equations

[TABLE]

for which $C$ can be eliminated to give

[TABLE]

which is the desired (18a). Similarly, through a local analysis at $z=z_{A}$ we find

[TABLE]

which is the desired (18b). In fact, the analysis near $z=z_{A}$ results in replacing the argument of $\Phi$ in (59) with the turning point $z_{A}$ and morever swapping $a_{2}\mapsto b_{2}$ and $b_{2}\mapsto a_{2}$, which is a result of our definition of $f_{\textrm{I}}$ in (17).

Appendix B Inner matching procedure for general $p$

The exponential asymptotics procedure requires deriving the form of the factorial/power ansatz of the late terms. From (27), (29), (30) we have

[TABLE]

where $\chi(z_{*})=0$ at the respective singularities given by $z=z_{*}$. The two numerical constants $\gamma$ and $\Lambda$ are determined by taking the limit of $z\to z_{*}$ and matching with the inner solution.

B.1 Derivation of $\gamma$ and $\Lambda$ for the turning points

It is somewhat easier to see the necessary re-scalings by examining the inner-region first. We let $z_{*}$ be a turning point and re-scale in the inner region with

[TABLE]

where $\mathrm{i}p(\mathrm{i}z_{*})^{p-1}K^{3}=-1$. The branch of $K$ is chosen so that $\phi$ in the later expression (65) corresponds to a decaying exponential, $\mathrm{e}^{\mathrm{i}\phi/\epsilon}$, as $y\to\infty$. Under the substitution (62), the main equation (9a) now yields the standard Airy equation

[TABLE]

The solution of (63) that matches with the requisite decay condition at infinity is $f=C\operatorname{Ai}(y)$ for a constant $C$. From [1], we may write the full expansion of the Airy function as

[TABLE]

valid as $y\to\infty$.

Now with the outer limit of the inner solution, we return to the outer expansion and take the inner limit. First, using the scalings (62) in (22a) and (29), we have

[TABLE]

and where $P$ is an $\mathcal{O}(1)$ constant that we do not need to specify. Next, from (22), the leading-order outer solution, $A_{0}$, written in the inner coordinates, satisfies

[TABLE]

while substituting of the scalings (65) into the general late-orders expression (61) gives

[TABLE]

For the power of the singularity in (66) to match (67) at $n=0$, we require $\gamma=0$.

To determine $\Lambda$ we match with the outer solution with the inner solution (64) in the vicinity of the turning point. First, notice that the behaviour in (64) must match in the limit $y\to\infty$ with the leading-order outer WKB solution given by $A_{0}\mathrm{e}^{\mathrm{i}\phi/\epsilon}$. With $\phi$ from (65) and $A_{0}$ from (66), we have that

[TABLE]

We now match the $n$th term of the outer-to-inner limit with the $n$th term of the inner-to-outer limit. This gives

[TABLE]

Substitution of $C$ from (68) into the above yields the much more compact expression,

[TABLE]

Notice that had we chosen the constant of integration in $A_{0}$ in (22b) to be a general value of $a_{0}$ (instead of $a_{0}=1$), then $\Lambda=a_{0}/(2\pi)$ above.

B.2 Derivation of $\gamma$ and $\Lambda$ for the singularity $z=0$

In order to determine the values of $\gamma$ and $\Lambda$ that correspond to the divergence due to $z=0$, we apply a similar procedure to that of Sec. B.2. Here, the main difference is that the full solution in the inner region does not exist in terms of special functions.

In the limit $x\to 0$, we have that $u\sim 2\mathrm{i}x$, and thus from (13b), the leading-order prefactor behaves as

[TABLE]

[TABLE]

The order of the singularity is correct when $n=0$ only if $\gamma=-p$. To determine $\Lambda$ we need to match with an inner region in the vicinity of $x=0$. The inner scaling is $A=1+\epsilon^{p}g$, $x=-i\epsilon y$, giving

[TABLE]

at leading order, with the condition that $g\sim-y^{p}/4$ as $y\to\infty$. We set $g=-y^{p}h/4$ to give

[TABLE]

with $h\to 1$ as $y\to\infty$. To match with the outer WKB solution expand as $y\to\infty$ as

[TABLE]

giving after simplification,

[TABLE]

with solution $h_{n}=\Gamma(n-p)/[2^{n}\Gamma(-p)]$. Hence as $y\to\infty$, the outer limit of the leading-order inner solution is

[TABLE]

The above outer limit of the inner solution must be matched with the inner limit of the outer solution, $A\sim\sum\epsilon^{n}A_{n}$, whose $n$th term is

[TABLE]

Thus by the Van Dyke matching rule [13], matching (75) with (76) gives

[TABLE]

Appendix C Derivation of the Stokes switching

In order to derive (36), we substitute

[TABLE]

into the equation for the remainder (34), where $B$ is given in (30) and $\chi$ is given in (29). The pre-factor $\mathcal{S}$ is the switching function, which is expected to rapidly vary across Stokes lines in the limit $\epsilon\to 0$. Using the late terms (31), this gives the expression for the switching function as

[TABLE]

Writing $\chi^{\prime}=2\mathrm{i}\phi^{\prime}$ and converting $\mathcal{S}^{\prime}=\chi^{\prime}\frac{\operatorname{d\!}{}\mathcal{S}}{\operatorname{d\!}{}\chi}$, we have

[TABLE]

By expanding $\Gamma(N+\gamma+1)$ in the limit $N\to\infty$, it can be shown that the right hand-side of (80) is algebraically small except if $N$ is chosen optimally, i.e. at the point where adjacent terms of the expansion are equal,

[TABLE]

or $N\sim|\chi|/\epsilon$. At optimal truncation, it can be shown that there exists a boundary layer near Stokes lines specified by (35) where $\mathcal{S}$ incurs a jump of magnitude $2\pi\mathrm{i}/\epsilon^{\gamma}$. Returning to (80), we note this is identical to eqn (4.4) of Chapman & Vanden-Broeck [8] with a negated right-hand side and where their $N$ is our $N-1$. The derivation of the jump conditions is identical otherwise and yields (36). More details of the optimal truncation and Stokes switching procedure can be found in [7].