Behavior of eigenvalues in a region of broken-PT symmetry

TL;DR

This paper investigates the eigenvalue behavior of PT-symmetric quantum Hamiltonians in the broken symmetry region, revealing new phenomena such as an infinite-order exceptional point and spectral transitions.

Contribution

It provides the first detailed numerical and analytical analysis of eigenvalues in the broken PT symmetry region for <, uncovering novel spectral transitions and critical points.

Findings

Discovery of an infinite-order exceptional point at 

Transition from discrete to partially continuous spectrum at 

Eigenvalue behavior near the conformal limit 

Abstract

PT-symmetric quantum mechanics began with a study of the Hamiltonian $H=p^2+x^2(ix)^\varepsilon$. When $\varepsilon\geq0$, the eigenvalues of this non-Hermitian Hamiltonian are discrete, real, and positive. This portion of parameter space is known as the region of unbroken PT symmetry. In the region of broken PT symmetry $\varepsilon<0$ only a finite number of eigenvalues are real and the remaining eigenvalues appear as complex-conjugate pairs. The region of unbroken PT symmetry has been studied but the region of broken PT symmetry has thus far been unexplored. This paper presents a detailed numerical and analytical examination of the behavior of the eigenvalues for $-4<\varepsilon<0$. In particular, it reports the discovery of an infinite-order exceptional point at $\varepsilon=-1$, a transition from a discrete spectrum to a partially continuous spectrum at $\varepsilon=-2$, a transition at the Coulomb value $\varepsilon=-3$, and the behavior of the eigenvalues as $\varepsilon$ approaches the conformal limit $\varepsilon=-4$.

Figures

Tables

Table 1. Table 1: Comparison of the real parts of the eigenvalues of the differential equation ( 15 ) at δ = 0.01 𝛿 0.01 \delta=0.01 with the asymptotic approximation in ( 29 ). The rate at which the accuracy increases with increasing k 𝑘 k is similar to the increase in accuracy of the standard WKB approximation to the eigenvalues of the quartic anharmonic oscillator R21 .

$k$	Numerical value of	$O\left({\delta }^2\right)$ calculation	relative
	$\mathrm{Re}{E}_k$ at $\delta =0.01$	of $\mathrm{Re}{E}_k$ in (29)	error
0	$-1.0352$	$-1.0414$	$8.70\%$
2	$-1.0426$	$-1.0461$	$0.33\%$
4	$-1.0469$	$-1.0493$	$0.30\%$
6	$-1.0499$	$-1.0518$	$0.18\%$
8	$-1.0523$	$-1.0538$	$0.15\%$
10	$-1.0542$	$-1.0555$	$0.12\%$
12	$-1.0559$	$-1.0569$	$0.10\%$

Table 2. Table 2: Comparison of the imaginary parts of the eigenvalues of the differential equation ( 15 ) at δ = 0.01 𝛿 0.01 \delta=0.01 with the asymptotic approximation in ( 29 ).

$k$	Numerical value of	$O\left({\delta }^2\right)$ calculation	relative
	$\mathrm{Im}{E}_k$ at $\delta =0.01$	of $\mathrm{Im}{E}_k$ in (29)	error
0	$0.03397$	$0.03210$	$5.3\%$
2	$0.03352$	$0.03220$	$3.8\%$
4	$0.03339$	$0.03224$	$3.4\%$
6	$0.03334$	$0.03226$	$3.2\%$
8	$0.03332$	$0.03228$	$3.1\%$
10	$0.03332$	$0.03229$	$3.0\%$
12	$0.03333$	$0.03231$	$3.0\%$

Table 3. Table 3: First three real discrete eigenvalues as a function of δ 𝛿 \delta , where ε = − 4 + δ 𝜀 4 𝛿 \varepsilon=-4+\delta . All the eigenvalues approach 0 as δ → 0 → 𝛿 0 \delta\to 0 . In fact, Fig.  19 indicates that they approach zero in a linear fashion.

$\delta$	First real	Second real	Third real
	eigenvalue	eigenvalue	eigenvalue
0.15	$0.173$	$0.440$	$0.807$
0.12	$0.114$	$0.321$	$0.628$
0.08	$0.080$	$0.230$	$0.454$
0.06	$0.060$	$0.177$	$0.351$
0.04	$0.035$	$0.116$	$0.236$
0.02	$0.012$	$0.049$	$0.106$

Equations

$$H=p^2+x^2(ix{)}^{\epsilon }$$

$$-y^{\prime \prime }(x)+x^2(ix{)}^{\epsilon }y(x)=Ey(x).$$

$${\theta }_{\mathrm{right}\mathrm{wedge},\mathrm{center}}$$

$${\theta }_{\mathrm{left}\mathrm{wedge},\mathrm{center}}$$

$${\theta }_{\mathrm{right}\mathrm{wedge},\mathrm{upper}\mathrm{edge}}$$

$${\theta }_{\mathrm{left}\mathrm{wedge},\mathrm{upper}\mathrm{edge}}$$

$${\theta }_{\mathrm{right}\mathrm{wedge},\mathrm{lower}\mathrm{edge}}$$

$${\theta }_{\mathrm{left}\mathrm{wedge},\mathrm{lower}\mathrm{edge}}$$

$$-ix\big{[}1+\delta\ln(ix)+{\rm O}\big{(}\delta^{2}\big{)}\big{]}.$$

$$\psi(x)=y_{0}(x)+\delta y_{1}(x)+{\rm O}\big{(}\delta^{2}\big{)}.$$

$$0=\psi_{\rm R}(x)\psi_{\rm L}^{\prime}(x)-\psi_{\rm L}(x)\psi_{\rm R}^{\prime}(x)\bigg{|}_{x=0}.$$

$$y_0^{\prime \prime }(x)+ix{y}_0(x)+E{y}_0(x)=0.$$

$$y_{0,(\mathrm{L},\mathrm{R})}^{\prime \prime }(r)-\left(r-E{e}^{\mp i\pi /3}\right)y_{0,(\mathrm{L},\mathrm{R})}(r)=0,$$

$$y_{0,(\mathrm{L},\mathrm{R})}(x)$$

$$\displaystyle\psi_{0,{\rm R}}(x)\psi_{0,{\rm L}}^{\prime}(x)-\psi_{0,{\rm L}}(x)\psi_{0,{\rm R}}^{\prime}(x)\big{|}_{x=0}$$

$$=-C_{\mathrm{L}}{C}_{\mathrm{R}}[e^{-i\pi /6}\mathrm{Ai}\left(E{e}^{-2i\pi /3}\right)A{i}^{\prime }\left(E{e}^{2i\pi /3}\right)$$

$$+e^{i\pi /6}\mathrm{Ai}\left(E{e}^{2i\pi /3}\right)A{i}^{\prime }\left(E{e}^{-2i\pi /3}\right)]$$

$$=-i{C}_{\mathrm{L}}{C}_{\mathrm{R}}\mathrm{W}\left[\mathrm{Ai}\left(E{e}^{2i\pi /3}\right),\mathrm{Ai}\left(E{e}^{-2i\pi /3}\right)\right]$$

$$=\frac{1}{2\pi }{C}_{\mathrm{L}}{C}_{\mathrm{R}}\ne 0.$$

$$y_1^{\prime \prime }(x)+ix{y}_1(x)+ix\ln (ix)y_0(x)+E{y}_1(x)=0.$$

$${\left[y_0^2(x)Q^{\prime }(x)\right]}^{\prime }=-ix\ln (ix)y_0^2(x).$$

$$Q_{\mathrm{L},\mathrm{R}}^{\prime }(x)$$

$$0$$

$$\mathrm{Ai}(x)\sim \frac{1}{2\sqrt{\pi }}{x}^{-1/4}\exp \left(-\frac{2}{3}{x}^{3/2}\right)(|x|\to \infty ,|\arg x|<\pi ).$$

$$\textstyle{\frac{2}{\delta}}\sim\textstyle{\frac{1}{\sqrt{E}}}\exp\left(\mbox{$\textstyle{\frac{4}{3}}$}E^{3/2}\right)\left[Q_{\rm R}^{\prime}(0)-Q_{\rm L}^{\prime}(0)\right].$$

$$\frac{2}{\delta }$$

$$\arg \frac{2}{\delta }=\arctan \left(-\frac{\beta }{2\alpha }\right)-2{\alpha }^{1/2}\beta =2m\pi ,$$

$$\mathrm{Re}E\sim {\left(-\frac{3}{4}\ln |\delta |\right)}^{2/3},\mathrm{Im}E\sim \frac{n\pi }{2\sqrt{\mathrm{Re}E}},$$

$$\epsilon =-2+\delta$$

$$-y^{\prime \prime }(x)-(ix{)}^{\delta }y(x)=Ey(x).$$

$$\displaystyle-y^{\prime\prime}(x)-\delta\ln(ix)y(x)-\mbox{$\textstyle{\frac{1}{2}}$}\delta^{2}[\ln(ix)]^{2}y(x)$$

$$=(E+1)y(x).$$

$$\displaystyle-y^{\prime\prime}(t)+\delta\ln(-t)y(t)+\mbox{$\textstyle{\frac{1}{2}}$}\delta^{2}[\ln(-t)]^{2}y(t)$$

$$=-(E+1)y(t).$$

$$\displaystyle-y^{\prime\prime}(t)+\delta[\ln(t)\pm i\pi]y(t)+\mbox{$\textstyle{\frac{1}{2}}$}\delta^{2}[\ln(t)\pm i\pi]^{2}y(t)$$

Full text

I Introduction

$\mathcal{PT}$-symmetric quantum theory has its roots in a series of papers that proposed a new perturbative approach to scalar quantum field theory: Instead of a conventional expansion in powers of a coupling constant, it was proposed that a perturbation parameter $\delta$ be introduced that measures the nonlinearity of the theory. Thus, to solve a $g\phi^{4}$ field theory one studies a $g\phi^{2}(\phi^{2})^{\delta}$ theory and treats $\delta$ as a small parameter. After developing a perturbation expansion in powers of $\delta$, the parameter $\delta$ is set to one to obtain the results for the $g\phi^{4}$ theory. This perturbative calculation is impressively accurate and does not require the coupling constant $g$ to be small T1 ; T2 . A crucial technical feature of this idea is that $\phi^{2}$ and not $\phi$ be raised to the power $\delta$ in order to avoid raising a negative number to a noninteger power and thereby generating complex numbers as an artifact of the procedure.

Subsequently, the $\delta$ expansion was used to solve an array of nonlinear classical differential equations taken from various areas of physics: The Thomas-Fermi equation (nuclear charge density) $y^{\prime\prime}(x)=[y(x)]^{3/2}/\sqrt{x}$ is modified to $y^{\prime\prime}(x)=y(x)[y(x)/x]^{\delta}$; the Lane-Emdon equation (stellar structure) $y^{\prime\prime}(x)+2y^{\prime}(x)/x+[y(x)]^{n}=0$ is modified to $y^{\prime\prime}(x)+2y^{\prime}(x)+[y(x)]^{1+\delta}$; the Blasius equation (fluid dynamics) $y^{\prime\prime\prime}(x)+y^{\prime\prime}(x)y(x)=0$ is modified to $y^{\prime\prime\prime}(x)+y^{\prime\prime}(x)[y(x)]^{\delta}=0$; the Korteweg-de Vries equation (nonlinear waves) $u_{t}+uu_{x}+u_{xxx}=0$ is modified to $u_{t}+u^{\delta}u_{x}+u_{xxx}=0$. In each of these cases the quantity raised to the power delta is positive and when $\delta=0$ the equation becomes linear. Just a few terms in the $\delta$ expansion gives an accurate numerical result T3 .

The breakthrough of $\mathcal{PT}$-symmetric quantum theory was the surprising discovery that to avoid the appearance of spurious complex numbers it is actually not necessary to raise a positive quantity to the power $\delta$ so long as the quantity is symmetric under combined space and time reflection. This fact is highly nontrivial and was totally unexpected. For example, a quantum-mechanical potential of form $x^{2}(ix)^{\varepsilon}$ does not necessarily lead to complex eigenvalues because the quantity $ix$ is $\mathcal{PT}$ invariant. Indeed, the non-Hermitian $\mathcal{PT}$-symmetric Hamiltonian

[TABLE]

has the property that its eigenvalues are entirely real, positive, and discrete when $\varepsilon\geq 0$ (see Fig. 1). The reality of the spectrum was noted in Refs. R1 ; R2 and was attributed to the $\mathcal{PT}$ symmetry of $H$. Dorey, Dunning, and Tateo proved that the spectrum is real when $\varepsilon>0$ R3 ; R4 . Following the observation that the eigenvalues of non-Hermitian $\mathcal{PT}$-symmetric Hamiltonians could be real, many papers were published in which various $\mathcal{PT}$-symmetric model Hamiltonians were studied R5 .

A particularly interesting feature of $\mathcal{PT}$-symmetric Hamiltonians is that they often exhibit a transition from a parametric region of unbroken $\mathcal{PT}$ symmetry in which all of the eigenvalues are real to a region of broken $\mathcal{PT}$ symmetry in which some of the eigenvalues are real and the rest of the eigenvalues occur in complex-conjugate pairs. The $\mathcal{PT}$ transition occurs in both the classical and the quantized versions of a $\mathcal{PT}$-symmetric Hamiltonian R2 and this transition has been observed in numerous laboratory experiments R6 ; R7 ; R8 ; R9 ; R10 ; R11 ; R12 ; R13 ; R14 ; R15 ; R16 ; R17 .

There have been many studies of the real spectrum of $H$ in (1) but essentially nothing has been published regarding the analytic behavior of the complex eigenvalues as functions of $\varepsilon$ in the region of broken $\mathcal{PT}$ symmetry. However, it is known that there is a sequence of negative-real values of $\varepsilon$ lying between $-1$ and [math] at which pairs of real eigenvalues become degenerate and split into pairs of complex-conjugate eigenvalues. These special values of $\varepsilon$ are often called exceptional points R18 . In general, eigenvalues usually have square-root branch-point singularities at exceptional points.

Exceptional points in the complex plane, sometimes called Bender-Wu singularities, explain the divergence of perturbation expansions R19 ; R20 . The appearance of exceptional points is a generic phenomenon. In these early studies of coupling-constant analyticity it was shown that the energy levels of a Hamiltonian, such as the Hamiltonian for the quantum anharmonic oscillator $H=p^{2}+x^{2}+gx^{4}$, are analytic continuations of one another as functions of the complex coupling constant $g$ due to the phenomenon of level crossing at the exceptional points. Thus, the energy levels of a quantum system, which are discrete when $g$ is real and positive, are actually smooth analytic continuations of one another in the complex-$g$ plane R21 . A simple topological picture of quantization emerges: The discrete energy levels of a Hamiltonian for $g>0$ are all branches of a multivalued energy function $E(g)$ and the distinct eigenvalues of this Hamiltonian correspond with the sheets of the Riemann surface on which $E(g)$ is defined. Interestingly, it is possible to vary the parameters of a Hamiltonian in laboratory experiments and thus to observe experimentally the effect of encircling exceptional points R13 ; R22 ; R23 .

The purpose of this paper is to study the analytic continuation of the real eigenvalues shown in Fig. 1 as $\varepsilon$ moves down the negative-$\varepsilon$ axis. In Sec. II we show that there is an infinite-order exceptional point at $\varepsilon=-1$ where there is an elaborate logarithmic spiral (a double helix) of eigenvalues. The real part of each complex-conjugate pair of eigenvalues that is formed at exceptional points between $\varepsilon=-1$ and $\varepsilon=0$ approaches $+\infty$ like $|\ln(\varepsilon+1)|^{2/3}$ as $\varepsilon$ approaches $-1$. In contrast, the imaginary parts of each pair of eigenvalues vanish logarithmically at $\varepsilon=-1$. As $\varepsilon$ goes below $-1$, the real parts of the eigenvalues once again become finite and the imaginary parts of the eigenvalues rise up from 0. As $\varepsilon$ goes from just above to just below $-1$, the imaginary parts of the eigenvalues appear to undergo discrete jumps but in fact they vary continuously as functions of $\varepsilon$.

In Sec. III we discuss the Stokes wedges that characterize the eigenvalue problem as $\varepsilon$ goes below $-1$. We give plots of the eigenvalues in the region $-2<\varepsilon<-1$ and perform an asymptotic analysis of the eigenvalues near $\varepsilon=-2$. As $\varepsilon$ approaches $-2$, the entire spectrum becomes degenerate; the real parts of all the eigenvalues approach $-1$ and the imaginary parts coalesce to [math].

Section IV presents a numerical study of the eigenvalues in the region $-4<\varepsilon<-2$. We show that a transition occurs at $\epsilon=-2$ in which the eigenspectrum goes from being discrete to becoming partially discrete and partially continuous. The continuous part of the spectrum lies on complex-conjugate pairs of curves in the complex-$\varepsilon$ plane. Another transition occurs at $\varepsilon=-3$ (the $\mathcal{PT}$-symmetric Coulomb potential); below $\varepsilon=-3$ some of the discrete eigenvalues become real. As $\varepsilon$ approaches the conformal point $\varepsilon=-4$, the eigenvalues collapse to the single value 0. Section V gives brief concluding remarks.

II Eigenvalue behavior as $\varepsilon\to-1$

II.1 Stokes wedges

The time-independent Schrödinger eigenvalue problem for the Hamiltonian $H$ in (1) is characterized by the differential equation

[TABLE]

The boundary conditions imposed on the eigenfunctions require that $y(x)\to 0$ exponentially rapidly as $|x|\to\infty$ in a pair of Stokes wedges in the complex-$x$ plane. This subsection explains the locations of these Stokes wedges.

As has been previously discussed at length, the potential $x^{2}(ix)^{\varepsilon}$ has a logarithmic singularity in the complex-$x$ plane when $\varepsilon$ is not an integer. Thus, it is necessary to introduce a branch cut. This branch cut is chosen to run from [math] to $\infty$ in the complex-$x$ plane along the positive-imaginary axis because this choice respects the $\mathcal{PT}$ symmetry of the Hamiltonian. This is because $\mathcal{PT}$ symmetry translates into left-right symmetry in the complex-$x$ plane (that is, mirror symmetry with respect to the imaginary-$x$ axis) R1 ; R2 . The argument of $x$ on the principal sheet (sheet 0 of the Riemann surface) runs from $-3\pi/2$ to $\pi/2$. On sheet 1, $\pi/2<{\rm arg}\,x<5\pi/2$, on sheet $-1$, $-7\pi/2<{\rm arg}\,x<-3\pi/2$, and so on.

As explained in Refs. R1 ; R2 , the Stokes wedges in which the boundary conditions on $y(x)$ are imposed are located in the complex-$x$ plane in a $\mathcal{PT}$-symmetric fashion. If $\varepsilon=0$, the Stokes wedges have angular opening $\pi/2$ and are centered about the positive-$x$ and negative $x$-axes on the principal sheet of the Riemann surface. As $\varepsilon$ increases from [math], the wedges get narrower and rotate downwards; as $\varepsilon$ decreases from [math], the Stokes wedges get wider and rotate upwards. WKB analysis provides precise formulas for the location of the center line of the Stokes wedges,

[TABLE]

the upper edges of the Stokes wedges,

[TABLE]

and the lower edges of the Stokes wedges,

[TABLE]

The locations of the Stokes wedges for eight values of $\varepsilon$ are shown in Fig. 2. As $\varepsilon$ decreases to $-1$, the opening angles of the wedges increase to $120^{\circ}$ and the upper edges of the wedges touch. At the special value $\varepsilon=-1$ the logarithmic Riemann surface collapses to a single sheet; the wedges fuse and are no longer separated. As a result there are no eigenvalues at all (the spectrum is null) R25 . When $\varepsilon$ goes below $-1$, the wedges are again distinct and no longer touch; the left wedge rotates in the negative direction and enters sheet $-1$ while the right wedge rotates in the positive direction and enters sheet $1$.

II.2 Numerical behavior of the eigenvalues as $\varepsilon$ decreases below [math]

Previous numerical studies of the (real) eigenvalues for $\varepsilon\geq 0$ were done by using the shooting method. However, when the eigenvalues become complex, the shooting method is not effective and we have used the finite-element method and several variational methods. We have checked that the eigenvalues produced by these different methods all agree to at least five decimal places.

Figure 1 may seem to suggest that the real eigenvalues disappear pairwise at special isolated values of $\varepsilon$. However, the eigenvalues do not actually disappear; rather, as each pair of real eigenvalues fuse, these eigenvalues convert into a complex-conjugate pair of eigenvalues. At this transformation point both the real and the imaginary parts of each pair of eigenvalues vary continuously; the real parts remain nonzero and the imaginary parts move away from zero as $\varepsilon$ goes below the transition point. A more complete plot of the eigenvalues in Fig. 3 shows that the real parts of each pair of eigenvalues decay slightly as $\varepsilon$ decreases towards $-1$, while the imaginary parts grow slowly in magnitude. However, just as $\varepsilon$ reaches $-1$ the real parts of the eigenvalues suddenly diverge logarithmically to $+\infty$ and the imaginary parts of the eigenvalues suddenly vanish logarithmically. Below $\varepsilon=-1$ the real parts of the eigenvalues rapidly descend from $+\infty$ and the imaginary parts of the eigenvalues rise up from [math]. This behavior is depicted in Fig. 3 and a detailed description of the region $-1.05<\varepsilon<-0.95$ is shown in Fig. 4.

II.3 Asymptotic study of the eigenvalues near $\varepsilon=-1$

Figure 3 shows that the eigenvalues are singular at $\varepsilon=-1$ and suggests that this singularity is more complicated than the square-root branch-point singularities that occur at standard exceptional points R26 . To identify the singularity we perform a local asymptotic analysis about the point $\varepsilon=-1$. We begin by letting $\varepsilon=-1+\delta$ and we treat $\delta$ as small ($\delta<<1$). This allows us to approximate the potential $x^{2}(ix)^{\varepsilon}$ in (1) as

[TABLE]

We also expand the eigenfunctions in powers of $\delta$:

[TABLE]

Because we are treating $\delta$ as small, the Stokes wedges have an angular opening close to $2\pi/3$ and are approximately centered about the angles $\theta_{\rm L}=-7\pi/6$ and $\theta_{\rm R}=\pi/6$. We construct solutions $\psi_{\rm L}(x)$ and $\psi_{\rm R}(x)$ in the left and right Stokes wedges. We then patch together these eigenfunctions and their first derivatives at the origin $x=0$. The patching condition is

[TABLE]

To zeroth order in powers of $\delta$ the Schrödinger eigenvalue equation $H\psi(x)=E\psi(x)$ reads

[TABLE]

Substituting $x=re^{i\theta_{\rm L,R}}$ reduces this equation to an Airy equation R24 for the zeroth-order eigenfunctions $y_{0,{\rm(L,R)}}(x)$ in the left and right wedges:

[TABLE]

where the derivatives are now taken with respect to $r$.

The boundary conditions on the eigenfunctions in each wedge require that $y_{0,{\rm(L,R)}}(r)\to 0$ as $r\to\infty$, so the solutions to (7) are Airy functions R24 :

[TABLE]

where $C_{\rm L,R}$ are multiplicative constants.

The right side of the patching condition (6) for the zeroth-order solutions is calculated from the Wronskian identity for Airy functions R24 :

[TABLE]

When $\delta$ is exactly [math], the potential is linear in $x$ and $y_{0,{\rm(L,R)}}(x)$ are the exact solutions to the Schrödinger equation. The above calculation shows that these solutions cannot be patched, and thus there are no eigenvalues at all when $\varepsilon=-1$ ($\delta=0$). This conclusion is consistent with Fig. 3, which shows that the real parts of all of the eigenvalues become infinite as $\varepsilon$ approaches $-1$. The fact that the spectrum is empty at $\varepsilon=-1$ is not a new result; the absence of eigenvalues of a linear potential was established in Ref. R25 .

Next, we perform a first-order {\rm O}\big{(}\delta^{1}\big{)} analysis. We set $y_{1}(x)=Q(x)y_{0}(x)$. (This substitution is motivated and explained in detail in Ref. R21 .) The first-order Schrödinger equation now reads

[TABLE]

We multiply this equation by the integrating factor $y_{0}(x)$ and insert the leading-order approximation to the eigenfunctions and obtain

[TABLE]

We then integrate this equation along the center ray of each Stokes wedge:

[TABLE]

Thus, to first order in $\delta$ with $\psi(x)=y_{0}(x)[1+\delta Q(x)]$ the patching condition (6) becomes

[TABLE]

where we have used the zeroth-order patching condition (6) and the leading-order eigenfunction (8). Note that because the Schrödinger equation is linear we are free to choose $Q_{\rm L}(0)+Q_{\rm R}(0)=0$.

For large $E$, we use the asymptotic expansion of the Airy function R24

[TABLE]

Thus, the patching condition for $|E|\to\infty$ becomes

[TABLE]

Note that because we are treating $\delta$ as small, the difference $Q_{\rm R}^{\prime}(0)-Q_{\rm L}^{\prime}(0)$ is approximately a positive real number. For real $E$ this difference is exactly real because $Q_{\rm R}^{\prime}(0)$ and $-Q_{\rm L}^{\prime}(0)$ are complex conjugates.

We expand the right side of (12) to first order in $\beta/\alpha$, where $\alpha={\rm Re}\,E>0$ and $\beta={\rm Im}\,E$. This expansion is justified because, as we can see in Fig. 3, the imaginary parts are small compared with the real parts near $\varepsilon=-1$. The patching condition (12) then becomes

[TABLE]

Hence, when $\delta$ is positive, we obtain the condition

[TABLE]

where $m$ is an integer. This result simplifies because the arctangent term is small; to leading-order we obtain $2\alpha^{1/2}\beta=2m\pi$. Similarly when $\delta<0$, we find that $2\alpha^{1/2}\beta=(2m+1)\pi$.

We conclude that for either sign of $\delta$ we obtain a simple formula for the real part of the eigenvalues. Specifically, if we combine the above three equations, we obtain \textstyle{\frac{2}{|\delta|}}\sim\alpha^{-1/2}\exp\left(\mbox{\textstyle{\frac{4}{3}}}\alpha^{3/2}\right). Hence, in the neighborhood of $\varepsilon=-1$ (that is, when $\delta$ is near [math]), the real parts of the eigenvalues are logarithmically divergent while the imaginary parts of the eigenvalues remain finite:

[TABLE]

where $n$ is an even integer for $\delta>0$ and $n$ is an odd integer for $\delta<0$. Evidently, the imaginary parts of the eigenvalues vary rapidly as $\varepsilon$ passes through $-1$ because there is a logarithmic singularity at $\varepsilon=-1$. A blow-up of the region $-1.05<\varepsilon<1.05$ is given in Fig. 4.

To visualize the behavior of the eigenvalues near $\varepsilon=-1$ more clearly, we have plotted the imaginary and real parts of the eigenvalues in the *complex *-$\varepsilon$ plane in the left and right panels of Fig. 5. Observe that the imaginary parts of the eigenvalues lie on a helix and that the real parts of the eigenvalues lie on a it double helix as $\varepsilon$ winds around the logarithmic singularity at $\varepsilon=-1$. This logarithmic singularity is an infinite-order exceptional point, which one discovers only very rarely in studies of the analytic structure of eigenvalue problems.

III Eigenvalue behavior as $\varepsilon\to-2$

In Fig. 6 we plot the first three complex-conjugate pairs of eigenvalues in the range $-2.0\leq\varepsilon\leq-1.1$. Note that the eigenvalues $E_{k}$ coalesce to the value $-1$ as $\varepsilon$ approaches $-2$. As $\varepsilon$ decreases towards $-2$ the real part of $E_{k}$ becomes more negative as $k$ increases, and the spectrum becomes inverted; that is, the higher-lying real parts of the eigenvalues when $\varepsilon$ is near $-1.7$ (for example) decrease as $\varepsilon$ decreases and they cross when $\varepsilon$ is near $-1.3$. This crossing region is shown in detail in Fig. 7.

The objective of this section is to explain the behavior of the eigenvalues as $\varepsilon$ approaches $-2$ by performing a local analysis near $\varepsilon=-2$. To do so we let

[TABLE]

and treat $\delta$ as small ($\delta<<1$) and positive. With this change of parameter (2) becomes

[TABLE]

The boundary conditions on $y(x)$, which we can deduce from Fig. 2, are that the eigenfunctions $y(x)$ must vanish asymptotically at the ends of a path that originates at $e^{-3\pi i/2}\infty$ in the complex-$x$ plane, goes down to the origin along the imaginary axis, encircles the origin in the positive direction, goes back up the imaginary axis, and terminates at $e^{\pi i/2}\infty$. The eigenfunctions are required to vanish at the endpoints $e^{-3\pi i/2}\infty$ and $e^{\pi i/2}\infty$.

We now make the crucial assumption that it is valid to expand the potential term in (15) as a series in powers of $\delta$. To second order in $\delta$ we then have

[TABLE]

In this form one can see that to every order in powers of $\delta$ the potential terms in the Schrödinger equation are singular at $x=0$. As a consequence, the solution $y(x)$ vanishes at $x=0$. (One can verify that $y(0)=0$ by examining the WKB approximation to $y(x)$; the prefactor $[V(x)-E]^{-1/4}$ vanishes logarithmically.)

We then make the change of independent variable $t=-ix$. In terms of $t$ (16) becomes

[TABLE]

This eigenvalue equation is posed on a contour on the real-$t$ axis that originates at $t=+\infty$, goes down the positive-real $t$ axis, encircles the origin in the positive direction, and goes back up to $e^{2\pi i}\infty$, and $y(t)$ is required to vanish at the endpoints of this contour. We then replace $\ln(-t)$ with $\ln(t)\pm i\pi$:

[TABLE]

Next, we make the scale change

[TABLE]

This converts (18) into the Schrödinger equation

[TABLE]

where the energy term $F$ is given by

[TABLE]

and the order $\delta$ term in the potential is given by

[TABLE]

Our procedure will be as follows. First, we neglect the $U(s)$ term in (19) because $\delta$ is small and we use WKB theory to solve the simpler Schrödinger equation

[TABLE]

Second, we find the energy shift $\Delta F$ due to the $U(s)$ term in (19) by using first-order Rayleigh-Schrödinger theory R21 ; to wit, we calculate the expectation value of $U(s)$ in the WKB approximation to $y_{0}(s)$ in (22). Having found $F=F_{0}+\Delta F$, we obtain the energy $E$ from (20):

[TABLE]

This approach gives a very good numerical approximation to the energies shown in Fig. 6.

The standard WKB quantization formula for the eigenvalues $F_{0}$ in a single-well potential $V(s)$ (the two-turning-point problem) is

[TABLE]

For (22) the potential $V(s)$ is $\ln(s)$ and the boundary conditions on $y_{0}(s)$ are given on the positive half line: $y_{0}(s)$ vanishes at $s=0$ and at $s=+\infty$. In order to apply (24) we extend the differential equation to the whole line $-\infty<s<+\infty$ by replacing $\ln(s)$ with $\ln(|s|)$ and consider only the odd-parity solutions. Thus, we must replace the integer $n$ in (24) with $2k+1$, where $k=0,\,1,\,2,\,\ldots$. The turning points are given by $s_{1}=-e^{F_{0}}$ and $s_{2}=e^{F_{0}}$. Hence, the WKB formula (24) becomes

[TABLE]

The substitution $s=ue^{F_{0}}$ simplifies this equation to

[TABLE]

and the further substitution $v=-\ln(u)$ reduces the integral to a Gamma function:

[TABLE]

Thus, the WKB approximation to the eigenvalues $F_{0}$ is

[TABLE]

which is valid for large $k$.

Next, we calculate the order-$\delta$ correction $\Delta F$ to (25) due to the potential $U(s)$ in (19). To do so we calculate the expectation value of $U(s)$ in the WKB eigenfunction $y_{0}(s)$ of (22):

[TABLE]

where $U(s)$ is given in (21).

Integrals of this type are discussed in detail in Chap. 9 of Ref. R21 . To summarize the procedure, in the classically-forbidden region beyond the turning point, $y_{0}(s)$ is exponentially small, and the contribution to the integral from this region is insignificant. In the classically-allowed region the square of the eigenfunction has the general WKB form

[TABLE]

where $C$ is a multiplicative constant and $\phi$ is a constant phase shift.

Making the replacement \sin^{2}\theta=\mbox{\textstyle{\frac{1}{2}}}-\mbox{\textstyle{\frac{1}{2}}}\cos(2\theta), we observe that because of the Riemann-Lebesgue lemma, the cosine term oscillates to zero for large quantum number $k$, and we may replace $[y_{0}(s)]^{2}$ in the integrals in (26) by the simple function \mbox{\textstyle{\frac{1}{2}}}[F_{0}-V(s)]^{-1/2}. Thus, the shift in the eigenvalues is given by

[TABLE]

After making the previous changes of variable $s=e^{F_{0}}u$ followed by $v=-\ln(u)$, we obtain

[TABLE]

which evaluates to

[TABLE]

Finally, we substitute $F=F_{0}+\Delta F$ in (23) to obtain the eigenvalues $E_{k}$:

[TABLE]

where $F_{0}$ is given in (25).

To verify these results, in Tables 1 and 2 we compare our numerical calculation of ${\rm Re}\,E_{k}$ and ${\rm Im}\,E_{k}$ with the asymptotic prediction in (29).

IV Eigenvalue behavior for $-4<\varepsilon<-2$

This section reports our numerical calculations of the eigenvalues for $\varepsilon$ between $-2$ and $-4$. We rotate $x$ in (2) by $90^{\circ}$ by making the transformation $s=ix$. In the $s$ variable the eigenvalue equation (2) becomes

[TABLE]

In the $x$ variable the center-of-wedge angles (3) are $-\pi+\varepsilon\pi/(8+2\varepsilon)$ and $-\varepsilon\pi/(8+2\varepsilon)$ but in the $s$ variable these angles are simply $\mp 2\pi/(4+\varepsilon)$. Thus, the integration contour makes $2/(4+\varepsilon)$ loops around the logarithmic branch-point at the origin in the complex-$s$ plane.

For example, if $\varepsilon=-3$ (this is the complex $\mathcal{PT}$-symmetric version of the Coulomb potential for which $H=p^{2}+i/x$ R27 ), the contour loops around the origin exactly twice; it goes from an angle $-2\pi$ to the angle $2\pi$. Looping contours for other complex eigenvalue problems have been studied in the past and have been called “toboggan contours” R28 . In the $\mathcal{PT}$-symmetric Coulomb case the contour is shown in Fig. 8. Figure 9 shows the contours for the cases $\varepsilon=-2.5$ and $\varepsilon=-3.5$.

To solve these eigenvalue problems with looping contours we introduce the change of variable

[TABLE]

which parametrizes the looping path in the complex-$s$ plane in terms of the real variable $t$. As $t$ ranges from $-1$ to $+1$, the path in the complex-$s$ plane comes in from infinity in the center of the left Stokes wedge, loops around the logarithmic branch-point singularity at the origin, and goes back out to infinity in the center of the right Stokes wedge. In terms of the $t$ variable the eigenvalue equation (30) has the form

[TABLE]

where $\psi(t)$ satisfies $\psi(-1)=\psi(1)=0$.

To solve this eigenvalue problem we use the Arnoldi algorithm, which has recently come available on Mathematica R29 . This algorithm finds low-lying eigenvalues, whether or not they are real. We apply the Arnoldi al-

gorithm to (32) subject to the homogeneous Dirichlet boundary conditions $\psi(-1+\eta)=\psi(1-\eta)=0$ and let $\eta\to 0^{+}$. There are two possible outcomes: (i) In this limit, some eigenvalues rapidly approach limiting values; these eigenvalues belong to the discrete part of the spectrum. (ii) Other eigenvalues become dense on curves in the complex plane as $\eta\to 0^{+}$; these eigenvalues belong to the continuous part of the spectrum.

IV.1 $\varepsilon$ slightly below $-2$

As soon as $\varepsilon$ goes below $-2$, the eigenvalues explode away from the value $-1$ (shown at the left side of Fig. 6). In Fig. 10 we plot about 100 eigenvalues for $\varepsilon=-2.0001$ and $-2.001$. In each plot we see both discrete and con-

tinuous eigenvalues. The continuous eigenvalues lie on a complex-conjugate pair of curves in the left-half plane; the discrete eigenvalues also lie in the left-half plane but closer to the real axis.

IV.2 Discrete and continuous eigenvalues

While the purpose of Fig. 10 is to show that the eigenvalues explode away from $-1$ as $\varepsilon$ goes below $-2$, it is also important to show how to distinguish between discrete and continuous eigenvalues. To illustrate this we apply the Arnoldi algorithm at $\varepsilon=-2.6$. Our results are given in Fig. 11 for $\eta=0.01$. The spectrum in this case is qualitatively different from the spectrum near $\varepsilon=-2$; there are now two pairs of curves of continuous eigenvalues, and these curves are now in the right-half complex plane. The discrete eigenvalues are still in the left-half complex plane but further from the negative real axis. There is an elaborate spectral structure near the origin and this is shown in Fig. 12. (We do not investigate this structure in this paper and reserve it for future research.)

We emphasize that when the Arnoldi algorithm is used to study a spectrum, it can only return discrete values. Thus, one must determine whether an Arnoldi eigenvalue belongs to a discrete or a continuous part of the spectrum. To distinguish between these two possibilities we study the associated eigenfunctions and observe how they obey the boundary conditions. Plots of discrete and continuous eigenfunctions associated with eigenvalues shown in Fig. 11 are given in Figs. 13 and 14.

In Fig. 13 we plot the absolute values of the eigenfunctions corresponding to the complex-conjugate pair of eigenvalues $E=-1.79\pm 4.31i$ for $\varepsilon=-2.6$. Observe that as $t$ approaches the boundaries $-1$ and $1$, the eigenfunctions decay to [math] exponentially. We conclude from this that the eigenvalues are discrete. This result can then be verified by taking finer cell sizes in the Arnoldi algorithm. As the cell size decreases, the numerical values of $E$ are stable. In contrast, in Fig. 14 in which the absolute values of the eigenfunctions corresponding to the pair of eigenvalues $E=-0.01\pm 0.18i$ are plotted, we see that the eigenfunctions vanish exponentially at one endpoint but vanish sharply at the other endpoint. We therefore identify these eigenvalues as belonging to the continuous spectrum. Decreasing the Arnoldi cell size results in a denser set of eigenvalues along the same curve.

IV.3 Complex Coulomb potential $\varepsilon=-3$

For the Coulomb potential $\varepsilon=-3$, (30) becomes

[TABLE]

which is a special case of the Whittaker equation

[TABLE]

with \mu^{2}=\mbox{\textstyle{\frac{1}{4}}} R24 . The boundary conditions are unusual (they differ from those in conventional atomic physics) in that $\psi(s)\to 0$ as $|s|\to\infty$ with ${\rm arg}(s)=\pm 2\pi$. Rather than performing an analytic solution to the eigenvalue problem, we simply present the numerical results, which are obtained by solving (32) with $\varepsilon=-3$. Figure 15 displays about 100 eigenvalues, which lie on two pairs of complex-conjugate curves in the left-half plane. These eigenvalues are part of the continuous spectrum. A blow-up of the region around the origin is shown in Fig. 16.

The Coulomb case $\varepsilon=-3$ is a transition point between the regions $\varepsilon>-3$ and $\varepsilon<-3$. In the first region the discrete eigenvalues occur in complex-conjugate pairs and there are no real discrete eigenvalues (as we see in Fig. 11). In the region $\varepsilon<-3$ the discrete spectrum includes both real and complex-conjugate pairs of eigenvalues in addition to the continuous spectrum. Figure 17 illustrates the typical distribution of eigenvalues in the latter region for the choice $\varepsilon=-3.8$. In Fig. 18 we display the eigenfunction for the real discrete eigenvalue $E=0.0804$. Unlike the eigenfunctions in Figs. 13 and 14 this eigenfunction is symmetric in $t$.

IV.4 Conformal limit $\varepsilon\to-4$

The limit $\varepsilon\to-4$ is the conformal limit of the theory and thus the behavior of the eigenvalues in this limit is interesting to determine. It is difficult to study this limit because the eigenvalue equation in the complex-$s$ plane follows a contour that loops around the origin many times when $\varepsilon$ is near $-4$. Indeed, the number of loops approaches $\infty$ as $\varepsilon\to-4$ and, as a consequence, we are less confident about the dependability of the Arnoldi algorithm that we are using to obtain our numerical results. Nevertheless, we have studied the spectrum for values of $\varepsilon$ that are slightly greater than $-4$ and examine the trend as $\varepsilon$ moves closer to $-4$. We find that in this limit the entire spectrum collapses to the origin. It is not easy to demonstrate this by studying the continuous part of the spectrum; these points merely become denser in the vicinity of the origin. However, the discrete eigenvalues move toward the origin as $\varepsilon\to-4$. In Table 3 we show the behavior of the first three real eigenvalues as $\delta\to 0$, where $\varepsilon=-4+\delta$. These data are plotted in Fig. 19. This figure suggests that the eigenvalues vanish linearly with $\delta$.

V Conclusions

In this paper we have studied the eigenvalues of $H$ in (1) for $-4<\varepsilon<0$ and we have shown that there is a rich analytic structure as a function of the parameter $\varepsilon$. We have identified transition points at the integer values $\varepsilon=0,\,-1,\,-2,\,-3$. Just above $\varepsilon=0$ the eigenvalues are all real and positive but below $\varepsilon=0$ the eigenvalues split sequentially into complex-conjugate pairs and all of the eigenvalues but one are complex below about $\varepsilon=-0.58$. At $\varepsilon=-1$ the real parts of the eigenvalues approach $\infty$ but the imaginary parts of the eigenvalues all vanish.

Below $\varepsilon=-1$ the eigenvalues are once again finite, but as $\varepsilon$ approaches $-2$ the entire spectrum coalesces to the value $-1$. Below $\varepsilon=-2$ the eigenvalues explode away from the value $-1$ and a new feature of the spectrum arises: The spectrum is partly continuous and partly discrete. The continuous part of the spectrum lies along complex-conjugate pairs of lines in the complex plane that begin near the origin and run off to $\infty$. By contrast, the eigenvalues belonging to the discrete part of the spectrum have negative real parts.

At the Coulomb value $\varepsilon=-3$ the continuous parts of the spectrum swing around to the negative complex plane and the discrete eigenvalues disappear. Below the Coulomb transition the discrete eigenvalues reappear and some of the discrete eigenvalues are now real. As $\varepsilon$ approaches the conformal point $-4$, the spectrum appears to implode to the origin.

Acknowledgements.

We thank S. Sarkar for initial discussions of this investigation. CMB thanks the Heidelberg Graduate School of Fundamental Physics for its hospitality and CS and ZW thank the Physics Department at Washington University for its hospitality.