A study of two PT -symmetric quantum mechanical systems

(1)

Theoretical Physics

A study of two PT -symmetric quantum mechanical

systems

Kristoffer Aronsen

karonsen@kth.se

Simon Sandell

simsan@kth.se

SA114X Degree Project in Engineering Physics, First Level

Department of Theoretical Physics

KTH Royal Institute of Technology

Supervisor: Tommy Ohlsson

(2)

Abstract

In this report we examine the concept of PT -symmetric quantum mechanics and analyze two such systems. PT -symmetric quantum mechanics is a new branch of quantum me-chanics, first proposed by Bender and Boettcher in 1998 and has been the subject of much research since then. Various important aspects of PT -symmetric quantum mechanics is explained, such as asymptotic behaviour of solutions of the Schrödinger equation, as well as the concept of the C operator that is unique to PT -symmetric quantum mechanics. We briefly discuss both the concept of broken and unbroken PT -symmetry, and the relation between PT -symmetric and Hermitian quantum mechanics through similarity transforms. We specifically direct our attention to two theoretical quantum mechanical systems described by Hamiltonians with the property of being PT -symmetric. For these systems we determine the energy levels, one by analytical means and the other by use of numerical methods. To determine the energy levels numerically we have used the built-in ODE-solver ode45 built-in MATLAB built-in conjunction with the equation solver fzero. Our results are compared to those of previous studies and our methods and results are discussed.

(3)

Sammanfattning

I denna uppsats undersöker vi begreppet PT -symmetrisk kvantmekanik och analyserar två sådana system. PT -symmetrisk kvantmekanik är en ny gren av kvantmekaniken, först framförd av Bender och Boettcher 1998, och har sedan dess varit ämnet i många rappor-ter och studier. Olika viktiga aspekrappor-ter av den PT -symmetriska kvantmekaniken förklaras, såsom asymptotiskt beteende för lösningar av Schrödingerekvationen, liksom begreppet C-operator som är unik för PT -symmetrisk kvantmekanik. Vi diskuterar kortfattat be-greppet brutnen och obruten PT -symmetri, och relationen mellan den PT -symmetriska och Hermitska kvantmekanik genom likhets-transformer. Vi riktar särskilt vår uppmärk-samhet mot två teoretiska kvantmekaniska system som beskrivs av Hamiltonianer med PT -symmetri. För dessa system bestämmer vi energinivåer, en genom analytiska metoder och den andra genom användning av numeriska metoder. För att bestämma energinivåer numeriskt vi har använt den inbyggda ODE-lösaren ode45 i MATLAB i samband med ekvationslösaren fzero. Våra resultat jämförs med tidigare studier och våra metoder och resultat diskuteras.

(4)

Chapter 1 Introduction

1.1 Briefly about classical quantum mechanics

Quantum mechanics is an area of physics that applies to the microscopic scale. The fun-damental idea of quantum mechanics is that energy is not actually continuous but instead need to be transferred in packages of certain sizes, called quanta. The notion that the energy states of a physical system could be discrete was first suggested by Ludwig Boltz-mann in 1877. Max Planck proposed in an attempt to solve the problem with the black body radiation that energy was emitted and absorbed in discrete quanta. This hypothesis was consistent with observations. In 1905 Albert Einstein used a quantum-based model to describe the photoelectric effect. Numerous papers were published during the first half of the twentieth century on the subject of quantization in physics, and the founda-tions of the theory of quantum mechanics was established. The mathematically rigorous formulation of quantum mechanics was proposed by Dirac, von Neumann, Hilbert and Weyl [1, 2, 3].

In quantum mechanics a physical system is described by an equation called the “Schrödinger equation”. The Schrödinger equation is the quantum mechanical analog of Newton’s second law of motion and describes the time evolution of a quantum me-chanical system. In the most general form the Schrödinger equation can be written as

i~∂

∂tΨ (t, r) = HΨ, (1.1) where Ψ is called the “wave function” of the system and H is a mathematical operator called the Hamilton operator or the Hamiltonian. If one solves the eigenvalue problem

Hψ = Eψ (1.2)

one can determine the energy levels E of the system, which are given by the eigenvalues of the operator H, i.e. its spectrum.

A key feature of the formulations made by Dirac et al. was that for the system to be physically viable the Hamiltonian needs to be Hermitian. This means that the Hamiltonian H is its own adjoint, or in other words, its own complex conjugate and transpose. The Hermiticity of H is taken as an axiom but the reasons for this is that it guarantees that the eigenvalues of H are real and that the time evolution operator, usually denoted U , is unitary, which means that its transpose is equal to its inverse, UT _{= U}−1_{. The unitarity of U is in turn necessary for the probability of the system to}

(6)

1.2 Introduction of PT -symmetric quantum mechanics

In 1998 Bender and Boettcher [4] published an article investigating the prospect that the demand for Hermiticity of a Hamiltonian in a quantum mechanical system might be unnecessarily strict. In earlier investigations made by Wu and Brower et al. [5, 6], non-Hermitian Hamiltonians were used to describe theoretical systems, and real energy eigenvalues were obtained despite the Hamiltonians being non-Hermitian. These publi-cations were at their time criticized for this fact. Bender and Boettcher realized that the reason that the spectra of the Hamiltonians used in these articles were real was because the Hamiltonians were invariant under the inversion of parity followed by the inversion of time. In other words they were PT -symmetric.

If a Hamiltonian is PT -symmetric we write H = HPT, where PT is defined as the successive application of parity P and time reversal T :

P : ˆp → −ˆp, ˆx → −ˆx

T : ˆp → −ˆp, i → −i (1.3) where ˆp and ˆx are the momentum and position operators and i denotes the imaginary unit.

Bender has since the first publication written multiple articles on the subject of PT -symmetric quantum mechanics and has also extended his theories into the fields of classi-cal mechanics as well as quantum field theory. He has shown that the systems investigated by Wu and Brower et al. while not being Hermitian, are in fact PT -symmetrical, and that this is the cause for their real spectra.

His theory has inspired extensive research in multiple areas of physics, optics and condensed matter physics amongst others. He argues that PT -symmetric quantum me-chanics has all the qualities required to replace conventional Hermitian quantum mechan-ics as a generalization. While there have been reports of observations of PT -symmetric systems in optics and condensed matter physics, according to Bender and others [7, 8] there is as of yet no experimentally observed quantum mechanical system that can be described by a PT -symmetric Hamiltonian.

1.3 Objective

In this report we try to illuminate the concept of PT -symmetric quantum mechanics and what theoretically differs it from classical Hermitian quantum mechanics. We also present two different cases of PT -symmetric quantum mechanical systems and determine the energy eigenvalues. The systems have been thoroughly analyzed in previous works, and we will in large be using the same methods as previous authors.

The first of the systems we will be analyzing was initially proposed by Bender and Boettcher [4]. This system is defined as

H = p2− (ix)N_, _(1.4)

where N is a real parameter that can be varied. Different values of N will give different energy eigenvalues and for some values of N the system will transition from having a fully real spectrum to a spectrum consisting of real and complex eigenvalues. When N is

(7)

equal to 2, this system becomes the quantum harmonic oscillator which is a well-known system in quantum mechanics.

The second system was proposed by Cannata et al. [9] and is analytically solvable. It is defined as

H = p2+ exp (2ix) (1.5)

and it can be constructed by superposition of the Bender-Boettcher potentials as

∞

X

N =0

(ix)N

(8)

Chapter 2 Background material

In this section we will introduce some techniques that will be used in the following inves-tigations of the systems, such as the concept of natural units and asymptotic relations of solutions to differential equations. We will also to a lesser extent explain some concepts of the theory of PT -symmetric quantum mechanics, namely the C operator, and bro-ken and unbrobro-ken PT -symmetry. We explain how PT -symmetric Hamiltonians belong to a broader class of so-called pseudo-Hermitian Hamiltonians and can be brought to Hermitian form by similarity transforms.

2.1 Natural units and mass

In this report we will express the Hamiltonians using natural units. Natural units are commonly used in various areas of physics to increase the readability of mathematical formulas and expressions. By assigning certain natural constants like Planck’s reduced constant, ~, or the speed of light, c, to the value of 1, many expressions take a simpler form without losing any information, since we can linearly rescale back to any relevant units. We will let ~

2m = 1 as it simplifies the Hamiltonians that we will investigate below.

2.2 Parity and time inversions in quantum mechanics

Parity transformation, P, is a reflection in position. In quantum mechanics this means that ˆx → −ˆx. Since ˆp can be written as ˆ_{p = −i~}_∂x∂ we also find that ˆp → −ˆp under parity transformation. Time inversion, T , takes t → −t as expected. This has the effect that ˆ

p → −ˆp since ˆp can be defined as ˆp = m∂ ˆ_∂tx. It also has a less obvious effect of switching the sign of the imaginary unit i → −i, which can be loosely shown by studying the commutation relation between ˆx and ˆp, which is also called the canonical commutation relation. The commutator in quantum mechanics is defined as

[ Â, ˆB]ϕ = Â ˆB − ˆB Âϕ,

where ˆA and ˆB are two operators and ϕ is some test function which we will omit below. The canonical commutation relation [ˆx, ˆ_{p] = i~ must be conserved under T :}

(9)

which leads to the equation

−i~ = (i~)T . (2.2) The last equation shows that T needs to change the sign of the right hand side of the equality. An operator that has the effect that i → −i is said to be anti-unitary and it was rigorously shown by Eugene Wigner that the quantum mechanical time inversion operator must be one such operator. It is also a known fact that an anti-unitary operator changes the sign of the energy to be negative, and in order to ensure positive energy, the operator must also exchange t → −t.

2.3 The C operator

For a quantum system to be physical, the state vectors of the Hilbert space needs to be positive under an inner product. An initial problem with the theory of PT -symmetric quantum mechanics was the fact that the modulus of the wave functions of PT -symmetric systems were not always positive. Bender showed that using the intuitive choice of inner product

hf |gi = Z

C

(f (x))PT g(x)dx, (2.3)

the eigenfunctions of eq. 1.4 in fact followed the rule

hψn|ψmi = (−1)nδnm (2.4)

under this inner product. The path C is generally a complex path chosen so that the state vectors or wave functions fulfill some boundary condition. To solve this problem, a new operator which commutes with H and PT and is a root of unity distinct from P must be introduced so that the inner product becomes

hf |gi = Z

C

(f (x))CPT g(g)dx. (2.5)

Bender showed that in position space C must have the form

C(x, y) =

∞

X

n=0

ψn(x)ψn(y), (2.6)

which implies that we need to know the wave functions of the Hamiltonian before we can determine the inner product. In classical Hermitian quantum mechanics the inner product is defined as

hf |gi = Z

f (x)∗g(x)dx, (2.7)

for any state vectors of the physical system. In this inner product the path is constrained to integrate along the real line and we are not free to analytically continue into the complex domain, in contrast to the inner product of PT -symmetric quantum mechanics, which allows analytical continuation as long as we stay within the Stokes wedges of the

(10)

eigenvalue problem. The definition of Stokes wedges will be explained in the following section.

To calculate eigenvalues of a quantum mechanical PT -symmetric system, knowing the expression for the C operator is not necessary. Numerical and analytical methods can be applied directly to the time-independent Schrödinger equation as usual. This is fortunate, since explicitly writing down the C operator for a system is non-trivial. For many systems the eigenfunctions are not known in analytical form and can only be numerically approximated.

2.4 Asymptotic relations

In this report we will solve differential equations with the boundary condition that the solution ψ satisfies

|ψ| → 0 as |x| → ∞. (2.8) However, some of the equations we consider will not have solutions that can satisfy this condition if we restrict the variable x to be real only, so in order to find any solutions at all we will have to let x be complex. In section 3.1 we consider a class of potentials which generally do not have analytical solutions, but we use numerical methods to find the eigenvalues and the approximate shapes of the solution functions. In order to solve the equation numerically we still need to know beforehand where we can expect to find solutions that satisfy eq. 2.8 as the method is iterative and integrates the equation along a single contour in the complex plane. Fortunately it is possible to determine these regions, which are defined as Stokes wedges, by studying the asymptotic relations of the equation. The specific procedure that will be described now only works when the problem is a homogeneous ordinary linear differential equation with an irregular singularity.

The differential equations we study in section 3.1 have the form

y00(x) + p(x)y(x) = 0 (2.9)

and if p(x) → ∞ as x → x0 faster than (x − x0) −1

the equation is said to have an irregular singularity at x = x0. It should be mentioned that x0 = ∞ is allowed.

When working with differential equations with irregular singularities, one can make the substitution

y(x) = exp [S(x)] (2.10) to study the behaviour of the solution near the singularity. With this substitution the original problem will transform to a differential equation in S(x), as the original dif-ferential equation is homogeneous, and thus let us divide the entire equation by eS(x)_.

To simplify the problem further, we assume that S00(x) (S0(x))2 near the irregular singularity. Depending on the given differential equation, additional approximations can sometimes be made, such as dropping constant terms or terms that are dominated by some other. After these alterations have been made to the original equation we are left with a much easier equation that can be solved analytically. This new equation will naturally not yield a solution to the first problem, but in most cases it can be used to approximate the behaviour of the real solution in a neighborhood of the singularity. The expression exp [S(x)] is called the asymptotic relation for y(x) near the irregular singular-ity. This procedure and its validity are discussed in further detail in the book “Advanced Mathematical Methods for Scientists and Engineers” by Bender and Orszag [10].

(11)

Since the singularity of the systems we study in section 3.1 are located at |x| = ∞ we can use the asymptotic relation to find the regions that satisfy the boundary condition 2.8. These regions will be bounded by what are defined as Stokes lines. Such lines are the complex paths along which the asymptotic relation of y(x) is oscillatory. The regions between these Stokes lines can either be regions where the asymptotic behaviour of y(x) is exponentially growing or declining. Regions with exponentially declining behaviour will contain some complex path where the solution declines most rapidly. These lines are referred to as the anti-Stokes lines. For the first system we discuss in this report, these regions will be wedge-shaped and can therefore be referred to as Stokes wedges.

2.5 Broken and unbroken PT -symmetry

While the Hermiticity property of a Hamiltonian is always a guarantee for a real spectrum, this is not the case for the PT -symmetry of a Hamiltonian. In some cases the spectrum of a Hamiltonian with PT -symmetry will not be real, but it can be proven that if the eigenstate functions are the same for H as for PT the spectrum of H will be entirely real. For a linear operator it holds that if the operator commutes the Hamiltonian of the system the linear operator will have the same eigenfunctions as the Hamiltonian. If PT would have been linear this would have been a useful relation, but the T operator is anti-linear and it follows that PT is as well. In PT -symmetric quantum mechanics it therefore becomes necessary to show that the eigenstate functions of H and PT are the same, and this is not trivial. When the eigenstate functions of H and PT are equal it is said that the Hamiltonian H has an unbroken PT -symmetry and elsewise H has broken PT -symmetry.

2.6 PT -symmetry and pseudo-hermiticity

In a series of articles by Mostafazadeh [11, 12, 13] it is shown that PT -symmetric Hamiltonians are categorized in a broader class of Hamiltonians called pseudo-Hermitian. He also show that for every pseudo-Hermitian Hamiltonian there exists a similarity trans-formation that transforms the Hamiltonian to Hermitian form. The received Hermi-tian Hamiltonian will have the same energy eigenvalues that the previous non-HermiHermi-tian Hamiltonian had. This led him to make the argument that in fact PT -symmetric Hamil-tonians are physically equivalent to Hermitian HamilHamil-tonians and states that the only advantage to study PT -symmetric systems over their corresponding Hermitian systems is that the Hermitian system might be more difficult to work with [14]. Regarding this Bender argues that while every PT -symmetric quantum mechanical system is just another representation of some Hermitian quantum mechanical system, the similarity transformation is always difficult to determine and that so far only one such transforma-tion has been found in practice. Bender also state that while the energy eigenvalues of the two systems may be equal, they are not physically equal [7].

(12)

Chapter 3 Investigation

In this section we will focus on two examples of PT -symmetric quantum mechanics. We will apply analytical and numerical methods to determine energy eigenvalues as well as eigenstate functions for the two systems. The two systems have been studied previously in articles by Bender and Cannata et al. but the calculations presented in these papers are often neglecting several steps, making comprehension harder for those new to the subject of PT -symmetry and complex differential equations. We will try to repeat the calculations and present them in a form more adapted to newer readers.

3.1 The Bender-Boettcher potentials

The Hamiltonian operator is H = p2 − (ix)N_{, and thus, the potential for the system is}

V (x) = − (ix)N where N is a positive number. In particular when N = 2 this system reduces to the quantum harmonic oscillator. This system was investigated thoroughly in [4] and subsequent papers by Bender. Here we aim to reproduce the quantum mechanical calculations and present them in more detail. As a formality, we start by checking that the potential is PT -symmetric:

V (x)P = −(ix)NP = − (i(−x))N, (3.1)

V (x)PT =(−ix))N

T

= − (ix)N = V (x). (3.2)

3.1.1 Analytical calculations

To obtain the physical quantized energy eigenvalues for different values of N , only relying on solutions to the Schrödinger equation defined along the real line is not enough. The analysis must be continued into the complex plane in order to obtain quantized energy levels when dealing with most PT -symmetric systems. In our case, when N > 4 the eigenfunctions do not fulfill the boundary condition |φ| → 0 when |x| → ∞, and even when 1 ≤ N < 4 we can gain precision in our numerical methods by integrating along a complex contour instead of the real line as in regular quantum mechanics. The ideal path in the complex plane to follow is defined by the anti-Stokes lines, which are the lines where the wave function declines most quickly. Derivation of the Stokes lines of the system 1.4 is therefore the main subject of this section.

(13)

We have the time-independent Schrödinger equation, −∂ 2_ψ(x) ∂x2 − (ix) N ψ(x) = Eψ(x), (3.3)

and want to find the asymptotic form of the solutions near the irregular singularity at |x| → ∞. We follow the procedure presented in section 2.4. Substituting ψ(x) = exp(S(x)) into eq. 3.3 gives a new differential equation

−S00

(x) − (S0(x))2 − (ix)N_{exp [S(x)] = E exp [S(x)],} _(3.4)

− S00(x) − (S0(x))2 = E + (ix)N. (3.5) This is a nonlinear second-order ordinary differential equation and cannot be solved easily by analytical means. Following the procedure we make a second assumption that S00(x) (S0(x))2_{, which brings the differential equation down to first order and makes}

it very easy to solve analytically

S0(x) = ±p− (E + (ix)N) ≈ ±p−(ix)N ₌_i2 _{= −1 = ±x}N₂_iN +2₂ _. _(3.6)

Since |x| is very large, we can set E ≈ 0 to simplify the expression. Now we can integrate to obtain the two solutions

S(x) = ±x N 2+1i N +2 2 N 2 + 1 . (3.7)

We throw away the arbitrary integration constant since |x| 0. We want our solutions to be the analytical continuation of the harmonic oscillator and the wave functions of the harmonic oscillator have asymptotic form ψn,HO ∼ exp

−x2 2 . So when we fix N = 2 we can demand S(x)|N =2= −x 2

2 . Thus from eq. 3.7 we have

S(x)|N =2 = ± x22+1i 2+2 2 2 2 + 1 = ±(−1)x 2 2 . (3.8)

The correct solution is the one with positive sign,

S(x) = +x N 2+1i N +2 2 N 2 + 1 . (3.9)

This is the expression for the asymptotic form of the solutions ψ(x).

The Stokes lines, where the solutions converge most rapidly, are found by letting Im S(x) = 0, or equivalently, arg S(x) = π ± 2nπ, where n is an integer. Note that the angles Arg S(x) = 0 ± 2nπ are not allowed, as those are the angles ensuring the fastest divergence of the asymptotic form exp [S(x)].

Making the substitution x = rexp(iθ) and using n = 1 we have the equation for the angle θ(N 2 + 1) + π 2( N 2 + 1) = π, (3.10) θ(N 2 + 1) = π 2(2 − 1 − N 2), (3.11) θ = π 2 1 − N₂ 1 + N₂ = − π 2 N − 2 N + 2 = θr. (3.12)

(14)

This angle is located in the lower right quadrant for N > 2 and is zero for N = 2. As it turns out this system have two distinct Stokes lines. The other can be found by letting n = −1. Then we have θ(N 2 + 1) + π 2( N 2 + 1) = −π, (3.13) which leads to θ = −_Nπ 2 + 1 − π 2 = −π 2 + N₂ + 1 N + 2 ! = −π − π 1 − N 2 N + 2 ! = −π + π 2 N − 2 N + 2 = θl. (3.14)

This angle is in the lower right quadrant for N > 2 and just as θr is the right half of

the real line, this angle is equal to −π for N = 2. The expression we stopped at makes it easy to see that the two angles are indeed each other’s reflection through the complex axis for arbitrary N . We now know the contours to integrate along for our numerical analysis, we start at the complex infinity at the angle θl and move towards the origin, and

then integrate outwards toward the infinity at the angle θr. This will ensure maximum

precision in the numerical calculations, and the speed of convergence determines how far we have to integrate before we can approximate the solution to be zero.A longer path of integration will lead to a larger accumulative error. The accumulative error is not significant when small, if it is not dominating the wave function, the solution will still converge. However, it will reduce the accuracy of the calculated eigenvalues, as the calculated solution will not tend to zero exactly, but will have the magnitude of the error. If the error grows to large, it will have a leading role in the solution, and therefore the method will fail completely.

One can also find the angles on which the solutions become oscillatory by letting Re S(x) = 0, or equivalently, arg S(x) = π₂ ± nπ, where again n is an integer. These angles are the Stokes lines, but they will not be of much use in our further analysis. However, if calculated, they will show us that the eigenvalue problem cannot be solved on the real number line for values of N ≥ 4 as mentioned above.

3.1.2 Numerical analysis

When the Stokes wedges have been found it is possible to determine the energy eigenvalues of the system through numerical methods. The eigenvalues are calculated by following a contour through the complex plane and the choice of contour is arbitrary as long as it is contained within the Stokes wedges when they are far away from the origin.

In previous papers by Bender et al. [4] the energy eigenvalues were calculated nu-merically by separating the real part and the imaginary part of the differential equation into a coupled system of two real second-order differential equations. The system was then solved using a version of the Runge-Kutta method. Explicit details of how this was performed has been left out from Bender’s publications, forcing us to construct our own method.

The second-order differential equation was rewritten as a system of first order equa-tions. The eigenfunctions could then be numerically calculated using MATLAB’s built-in function ode45 for solving systems of ordinary differential equations. The method works

(15)

by taking starting values of the wave function and its derivative, and transforming the system of ordinary differential equations into a system of difference equations, which are then integrated along a given contour.

For values of E that are not eigenvalues, the integration will diverge, and since the wave functions are complex-valued, generally to some complex infinity. By running the solver for different values of E and looking at sign-changes of the divergence of the wave functions, one can get a rough estimate of intervals that contain eigenvalues. When the rough locations are known, more accurate methods can be used to determine the eigenvalues with more precision. This can be performed either by hand, or by utilizing other built-in methods such as fzero.

This procedure is common practice in quantum mechanics and is described in [15]. A notable difference is that the eigenvalue problem is not confined to the real line and this caused some difficulties when constructing the numerical method.

3.2 Periodic PT -symmetric potential

Another PT -symmetric system has been examined by F. Cannata et al. [9]. In contrast to eq. 1.4 this potential is analyticaly solvable and no numerical analysis will be necessary in determining the levels of this system. Cannata et al. presented analysis of this system along with two other analytically solvable systems in [9].

The Hamiltonian of the system is given by

H = p2+ exp (2ix) (3.15)

and again if we operate with PT on this Hamiltonian

HPT =(−p)2+ exp (−2ix)T = p2+ exp (2ix), (3.16) we find that H = HPT. Note that this potential function V (x) = eix can be viewed as a superposition of the group of potentials V (x) = −(ix)N _{we study in section 3.1, by}

noting the relation

∞

X

N =0

(ix)N

N ! = exp (ix). (3.17)

It is also periodic with period π,

exp (2i(x + π)) = exp (2ix) exp (2πi) = exp (2ix). (3.18)

3.2.1 Analytical calculations

The time-independent Schrödinger equations for this system is

−∂

2_ψ(x)

∂x2 + exp (2ix)ψ(x) = Eψ(x). (3.19)

The first step in analyzing this system is to realize that this equation can be brought to the form of Bessel’s differential equation [16] by making the substitution

(16)

In the equation, we have a second order partial derivative with regard to x, which we need to transform. By using the product rule for derivatives

∂2 ∂x2 = ∂z ∂x ∂ ∂z ∂z ∂x ∂ ∂z = ∂z ∂x ∂ ∂z ∂z ∂x ∂ ∂z + ∂z ∂x ∂2 ∂z2 =By definition: ∂z ∂x = ∂ ∂x[exp (ix)] = iz = iz ∂ ∂z (iz) ∂ ∂z + iz ∂2 ∂z2 = −z ∂ ∂z − z 2 ∂2 ∂z2. (3.21)

Substituting equations 3.20 and 3.21 into equation 3.19 gives

z2∂ 2_ψ(z) ∂z2 + z ∂ψ(z) ∂z + z 2 _{− α}2_{ψ(x) = 0,} (3.22)

where α2 = E. This is Bessel’s differential equation with the general solutions

ψ1(x) = J√E(e ix_), _(3.23) ψ2(x) = Y√E(e ix_), _(3.24) ψ3(x) = H√(1)_E(eix), (3.25) ψ4(x) = H (2) √ E(e ix_). _(3.26) The functions Jν, Yν, H (1)

ν , Hν(2) are Bessel functions of first and second kind and the

Hankel functions of first and second kind respectively. Any two of these four solutions can be chosen as the linearly independent basis of the solution space. We want our wave functions to be normalizable, so we need to know the asymptotic behaviour of the solutions for |x| → ∞. These formulas have already been calculated, but we will go through the procedure of arriving at the asymptotic forms, since the mathematics is non-trivial.

To find the asymptotic behaviour of the solutions we will use the method of steepest descent described in [17]. There are four conditions a function f (t) must fulfill for this method to be of use:

• f (t) can be represented by an integral on the generic form R

CF (z, t)dz.

• The integration path can be smoothly varied so that for large t the dominant contribution to the integral arises from a small range of z in the neighbourhood of the point z0 where |F (z0, t)| is a maximum on the path.

• The integration path passes through the point z0 in the orientation that causes the

most rapid decrease in |F | on departure from z0 in either direction along the path.

• In the limit of large t the contribution to the integral from the neighborhood of z0

asymptotically approaches the exact value of f (t).

(17)

in the complex plane: H_ν(1)(t) = 1 πi Z C1 exp(t/2)(z − z−1 ) dz zν+1, (3.27) H_ν(2)(t) = 1 πi Z C2 exp(t/2)(z − z−1 ) dz zν+1, (3.28) (3.29)

where the contours C1 and C2 are any contours that go from z0 = 0+ to zend = ∞ and

arg zend = +π for C1, whereas arg zend = −π for C2. Thus the first and second conditions

are fulfilled. The third condition can be fulfilled by studying the integrand in formula 3.27 and making appropriate constraints on the path. The fourth condition is realized in much the same way by studying the limit t → ∞ at the saddle point.

The method of steepest descent tells us that the asymptotic behaviour of a function can be obtained by Z C g(z, t)expw(z, t)]dz ≈ g(z0, t)exp " w(z0, t) + iθ s 2π |w00_(z 0, t)| # , (3.30)

where z = z0 is a saddle point of w(z, t) and

θ = −arg(w 00_(z 0, t)) 2 + π 2 or 3π 2 . (3.31)

The equation for θ ensures that the third condition is fulfilled if we know a saddle-point of the integrand. It is easy to realize that the Hankel functions can be brought on the form of the integral above. By comparing equations 3.27 and 3.30 we realize that if we make the following definitions

g(z, t) = 1 πiz

−(ν+1)

, (3.32)

w(z, t) = (t/2)(z − z−1), (3.33)

the integral formulas are precisely on the form required to use the method of steepest descent. The saddle points are found by

w0(z0) = (t/2) 1 + 1 z2 0 = 0, (3.34) z0 = ±i. (3.35)

For H(1) we can deform the contour C1 to pass through z = i and for H(2) the point is

z = −i. The angle that ensures that the third condition is fulfilled is given by eq. 3.31 and is for the first Hankel function

θ(1) = −arg(w 00_{(i, t))} 2 + π 2 or 3π 2 = π 4 + π 2 or 3π 2 = 3π 4 or 7π 4 . (3.36)

We must choose the first value θ(1) ₌ 3π

4 since if we choose to cross the saddle point in the

(18)

to cross the imaginary axis once more at a location further away from the origin where the absolute value of the integrand is larger than at z0. For the second Hankel function

we have θ(2) = −π 4 + π 2 or 3π 2 = π 4 or −3π 4 , (3.37)

and must choose θ(2) ₌ −3π 4 .

We are now ready to use the method of steepest descent. We insert the obtained angles and saddle points into eq. 3.30 and obtain

H_ν(1)(t) ≈ 1 πiexp iπ 2 (−ν − 1) exp 3iπ 4 exp(it) r 2π t = r 2 πtexp h it − νπ 2 − π 4 i (3.38) H_ν(2)(t) ≈ r 2 πtexp h −it − νπ 2 − π 4 i . (3.39)

When deriving these formulas, we assumed real values of t, but fortunately these formulas can be analytically continued to the complex domain. We have from [17]

H_ν(1)(z) ≈ r 2 πzexp h iz − νπ 2 − π 4 i , − π < arg z < 2π, (3.40) H_ν(2)(z) ≈ r 2 πzexp h −iz − νπ 2 − π 4 i , − 2π < arg z < π. (3.41) Finally we make the substitution 3.20, z = exp(ix), and set ν = √E to obtain the final expressions for the asymptotic behaviours of the wave functions to the studied problem, namely ψ3(x) ∼ r 2 π exp −ix 2 exphi exp(ix) − iπ 4 2√E + 1i, (3.42) ψ3(x) ∼ r 2 π exp −ix 2

exph−i exp(ix) + iπ 4

2√E + 1i. (3.43) Since Jν and Yν are the real and imaginary part of Hν(1), respectively, albeit for real

z, but again possible to analytically continue into the complex domain, the asymptotic behaviour of the final two wave functions are

ψ1(x) ∼ r 2 πexp −ix 2 coshi exp(ix) − iπ 4 2√E + 1i, (3.44) ψ2(x) ∼ r 2 πexp −ix 2 sinhi exp(ix) − iπ 4 2√E + 1i, (3.45) both valid for −π < arg z < π.

Knowing the asymptotic behaviour of the solutions, we start to study which solutions are physical. Rewriting x = Re x + i Im x makes the analysis easier. We restrict our analysis to positive energies, E > 0. In the complex upper half-plane, Im x > 0, we have

(19)

which tends to zero for large |x|. Bessel’s function of second kind is known to be divergent for arguments equal to 0, so the solutions two to four will not be normalizable along any contour in the complex upper half-plane and we are left with the solution ψ1, which

is normalizable, but with no condition on E, and therefore giving rise to a continuous eigenvalue spectra without a lower bound. Thus it is not a physical solution.

In the lower complex half-plane,

z = exp (iRe x) exp (|Im x|), (3.47) and we observe that we cannot make the same argument as for the upper half-plane, as the variable of the Bessel functions will go to some complex infinity. The asymptotic forms will now come in handy to tell us which solutions are physical. Looking at eqs. 3.42 to 3.43 we see that the first exponential will tend to zero when Im x− > −∞, but it will be dominated by the tetrated exponential term. The second term in eq. 3.42 gives

ψ3(x) ∼ exp

i exp (iRe x + |Im x|) − iπ 4 2√E + 1 = exp

exp (|Im x|) (i cos (Re x) − sin (Re x)) − iπ 4 2 √ E + 1 , (3.48)

and thus ψ3(x) is normalizable when sin (Re x) > 0. It is convenient to define

Sl = {x ∈ C|Im x < 0, lπ ≤ Re x < (l + 1)π} . (3.49)

Now we note that of ψ3(x) is normalizable when x ∈ S2l, where l is some integer. Since the

asymptotic form of ψ4(x) only differs from that ψ3(x) by a sign in the second exponential

term, we realize that ψ4(x) will be normalizable in the regions S2l+1. At the boundaries

Re x = nπ, the second exponential becomes oscillatory for both ψ3(x) and ψ4(x), and

thus, the first exponential will ensure normalizability. Remembering the recollecting that

sin (t) = exp (it) − exp (−it)

2i , (3.50)

cos (t) = exp (it) + exp (−it)

2 , (3.51)

we can quickly convince ourselves that ψ1(x) and ψ2(x) will not constitute normalizable

solutions in the lower complex half-plane.

Now that we know the normalizable solutions, we seek constraints on the energy. We are free to choose any contour in the lower complex half-plane. If we choose a contour confined to for example S0, we know that ψ3 is the solution, but we have no constraint

on the energy, leading again to a continuous and unbounded spectra. We can study the case of a contour starting at Im x = −∞ in S0, but ending in another sector of the

lower complex half-plane, Sm, m > 0, without loss of generality. From eq. 3.40 we have

the constraint that −π < arg z < 2π, equivalent to −π < Re x < 2π, so in the general case when the end sector is not within this small region of the complex plane, we must analytically continue the solutions. From [16] the analytical continuations of the Hankel functions are given by

H_ν(1)[z exp (imπ)] = −sin [(m − 1)πν] sin (πν) H (1) ν (z) − exp (−iπν) sin (πνm) sin (πν) H (2) ν (z), (3.52)

H_ν(2)[z exp (imπ)] = sin [(m + 1)πν] sin (πν) H (2) ν (z) + exp (iπν) sin (πνm) sin (πν) H (1) ν (z). (3.53)

(20)

Now, if m = 2l, we must demand that the Hν(2)(z) term disappears, as it is exponentially

growing in these sectors. Thus we have the conditions

sin (πν2l) = 0, (3.54) sin (πν) 6= 0. (3.55)

However ν = √E, so this is an eigenvalue problem. In fact, we have infinitely many eigenvalue problems as we can choose m arbitrarily large. The fist condition gives

2πl√E = nπ, (3.56) √

E = n

2l, (3.57)

whereas the second condition gives √

Eπ 6= kπ, (3.58)

√

E /_{∈ N.} (3.59)

These two restrictions completely describe the energy levels for a given l. As an example, when l = 42, we have the spectrum

En =

n2

7056, n ∈ {N \ {84, 168, 336, . . . }} . (3.60) Letting instead m = 2l + 1, we must demand that the Hν(1) term in eq. 3.52 goes to zero.

This yields the exact same spectrum. There exists however a special case, when l = 0 the first term in eq. 3.52 vanishes and the second term can be reduced so that the equation becomes

H_ν(1)[z exp (iπ)] = − exp (−iπν)sin (πν) sin (πν)H

(2)

ν (z) = − exp (−iπν)H (2)

ν (z). (3.61)

From this it is impossible to derive a constraint on E, the spectrum is again continuous, and this case does not constitute a physical system.

3.3 Results

In this section we present the results of our analytical and numerical analyses of the two PT -symmetric quantum systems subject of this report. In the figures we present the numerically calculated eigenvalues for a range of values of N , as well as the eigenfunction for a specific case.

3.3.1 Bender-Boettcher potentials

Several different numerical methods were applied, to varying degree of success. Most successful in recalculation of the energy levels was a method using ode45. In figure 3.1 the resulting values are plotted and in table 3.1 numerical values for some values of N are listed.

In figure 3.2 one of the numerically calculated wave function for N = 5 and E = 8.59 is presented and in figure 3.3 the wave function for the same energy and N , but along a path much closer to one of the Stokes lines is presented.

(21)

1 2 3 4 5 6 0 5 10 15 20 25 N E n e rg y

Figure 3.1: The eigenvalues of the system in section 3.1. The vertical perforated lines mark N = 1 and N = 2. When N = 1 there is no solution to the eigenvalue problem of the system. In the region between N = 1 and N = 2 the PT -symmetry is broken and it can be seen how neighboring eigenvalues float closer together with decreasing N .

N E Eref 2 1 -3 -5 -3 1.1563 1.1562 4.1093 4.1092 7.5632 7.5621 4 1.4771 1.4771 6.0039 6.0033 11.807 11.8023 5 1.9083 -8.5877 -17.711

-Table 3.1: In this table some choice values of E are presented along with their respective values of N . For comparison some values from ref. [4] are included in the column labeled Eref.

(22)

Re x

-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5

R

e

ψ

×104 -2 -1 0 1 2

Re x

-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5

Im

ψ

×104 -2 -1 0 1 2

Figure 3.2: The numerically determined wave function for N = 5 and E = 8.587795. The equation has been integrated along the anti-Stokes line for the corresponding value of N . Note that the span of the y-axes are from −2 · 104 _{to 2 · 10}4 _{as the wave function is not}

normalized. This has no effect on the eigenvalue as the eigenfunction that corresponds to it is determined up to a scaling constant.

(23)

Re x

-6

-4

-2

0

2

4

6 R

e

ψ

-20

-10

0

10

20 Re x

-6

-4

-2

0

2

4

6 Im

ψ

-20

-10

0

10

20

Figure 3.3: Wave function for N = 5 and E = 8.587795 along a contour contained within the Stokes wedges, but with different angles than the anti-Stokes lines. The function displays a more oscillatory behaviour than it does along the real line and it also converges more slowly than for the case in figure 3.2.

(24)

Re

x

-6

-4

-2

0

2

4

6 Im

x

-6

-4

-2

0

2

4

6 θ

Stokes

= π/7

θ = 0.14π

θ = 0

Figure 3.4: The blue contour is the contour of the wave function of figure 3.2 and the red contour is for figure 3.3. The angle θ is the angle between anti-Stokes lines and the contour for which it represents. Additionally θStokes is the angle between the anti-Stokes

lines and the Stokes lines. The blue contour follows the anti-Stokes lines completely, while the red contour lies very close to the Stokes lines. The perforated beams represent the Stokes and anti-Stokes lines.

(25)

3.3.2 Periodic potential

Our analysis of the PT -symmetrical quantum system 1.5 resulted in an infinite amount of eigenvalue-problems, characterized by choice of ending region of the complex path. The energy-levels were, in the general case of staring in the region S0 and ending in

either S2l or S2l+1, calculated to be

En =

n2

4l2, n ∈ {N \ {2l, 4l, 6l, . . . }} . (3.62)

This formula is identical to that presented by Cannata et al. [9], disregarding the different definition of natural units used.

3.4 Discussion

In this section we discuss both the theoretical aspects of PT -symmetry and our own analytical and numerical analysis presented in this report.

3.4.1 Difference between PT -symmetric and Hermitian quantum

mechanics

One notable difference between Hermitian and PT -symmetric quantum mechanics is the fact that we have allowed the positional parameter x to be complex. The physical interpretation of this is not yet established nor is it discussed in any literature we have found. From a theoretical standpoint PT -symmetric quantum mechanics introduces more layers of difficulty in treating systems since it is not possible to know beforehand if a system has a real spectrum or not.

3.4.2 The difficulties of the numerical analysis

During the numerical analysis of the system discussed in section 3.1 we encountered difficulties finding a proper numerical algorithm to use. Many algorithms for solving differential equations numerically seem implicitly to rely on the fact that the variable is a real number. In the systems we have considered the differential equations were defined in the complex plane. The built-in ODE solver in MATLAB could only handle values of N < 6. When we attempted to solve the differential equation for N ≥ 6 we received messages saying that the interval needed to be either strictly increasing or decreasing for the algorithm to work. This is a strange message as for complex numbers one cannot really speak of increasing or decreasing since they have no inherent internal ordering. This most likely is due to that the solver is not explicitly designed for complex differential equations. It would be interesting to investigate this further and to find another method or modify the existing code to work for N ≥ 6, for consistency. There is no theoretical reason for why the eigenvalue problem could not be solved for N ≥ 6, except for the limit where N → ∞.

Another problem with the numerical analysis was that the convergence or divergence of the wave functions were dependent on the choice of starting distance from the origin. If the distance was chosen too large, the accumulative error inherent to most numerical methods for solving ordinary differential equations starts to dominate the behaviour of

(26)

the solution, and thus the energies became less precise or even impossible to determine. As an example of that, when N = 2 the system is equal to the quantum harmonic oscillator and the energies to this are known (in our choice of natural units) to be

En= 1 + 2n. (3.63)

We can solve the Schrödinger equation numerically as described above starting somewhere on the negative real axis and using one of the energies according to the formula 3.63. We would then expect the wave function to converge as we know that the energies are correct. However if the starting point is too far off from the origin the correct energy has to be replaced by an incorrect energy for the wave function to converge. If the starting point is chosen even further away convergence even becomes completely unattainable. We used the distance that gave best result for the harmonic oscillator to calculate the energies for the other values of N . As a rule, larger value of N required a larger distance of integration.

3.4.3 Observations of PT -symmetric systems

This report has been dedicated to PT -symmetric quantum mechanics. We have not been able to find any literature on experimental observations of strictly quantum mechanical systems with this property, but there are multiple articles concerning observations of PT -symmetric optical systems. However, we are uncertain as to if the optical systems described in these articles can be considered to be quantum mechanical or not. We have not spent much time on these articles.

(27)

Chapter 4 Summary and conclusions

In this report we have presented and studied PT -symmetric quantum mechanics, in par-ticular we have analytically and numerically analyzed two PT -symmetric Hamiltonians to obtain energy eigenvalues and wave functions.

PT -symmetric quantum mechanics is a new branch of quantum mechanics where the classical demand for Hermiticity of the Hamiltonian is replaced by the condition of invariance under parity and time reversal. The replacement of Hermiticity with the less restrictive PT -symmetry admits physical systems that were previously forbidden. This relaxation does however come at the cost of some additional complications in the form of a need for a new operator C, and the proof of unbrokenness of the PT -symmetry.

In our study of the two systems we have made an effort of presenting all of the steps in the calculations in a clear way so that the reasoning can be easily followed for readers with less knowledge about complex differential equations. In previous papers these steps were not often easy to follow nor very clear.

The first system studied required both analytical and numerical analyses to obtain energy-levels and approximate eigenfunctions. We used well defined analytical methods for calculating the asymptotic forms of the solutions to the Schrödinger equation, and utilized these results in our numerical methods to obtain the spectra for the family of PT -symmetrical potentials studied. We presented calculated eigenvalues for a range of N and a plot of an approximate eigenfunction for a specific N and energy.

The second system was completely analytically solvable. By a variable transformation the Schrödinger equation could be rewritten in a form of Bessel’s equation, and then solved in terms of Bessel functions. Regional studies then revealed that the system harbours infinitely many eigenvalue problems defined by the choice of complex contour. A general formula for the systems energy-levels was presented.

(28)

References

[1] Dirac PAM. The Principles of Quantum Mechanics. Oxford University Press; 1930.

[2] von Neumann J. Mathematische Grundlagen der Quantenmechanik. Springer; 1932.

[3] Weyl H. Quantenmechanik und Gruppentheorie. Z. Phys. 1927;46.

[4] Bender CM, Boettcher S. Real Spectra in Non-Hermitian Hamiltonians Having PT -symmetry. Phys. Rev. Lett. 1998;80.

[5] Wu TT. Ground State of a Bose System of Hard Spheres. Phys. Rev. 1959 Sep;115.

[6] Brower RC, Furman MA, Moshe M. Critical exponents for the Reggeon quantum spin model. Phys. Lett. B. 1978;76.

[7] Bender CM. Making sense of non-Hermitian Hamiltonians. Rep. Prog. Phys. 2007;70.

[8] Brody DC. Consistency of PT -symmetric quantum mechanics. J. Phys. A. 2016;49.

[9] Cannata F, Junker G, Trost J. Schrödinger operators with complex potential but real spectrum. Phys. Lett. A. 1998;246.

[10] Bender CM, Orszag SA. Advanced Mathematical Methods for Scientists and Engi-neers I. Springer New York; 1999.

[11] Mostafazadeh A. Pseudo-Hermiticity versus PT symmetry: The necessary condition for the reality of the spectrum of a non-Hermitian Hamiltonian. J. Math. Phys. 2002;43.

[12] Mostafazadeh A. Pseudo-Hermiticity versus PT -symmetry II: A complete charac-terization of non-Hermitian Hamiltonians with a real spectrum. J. Math. Phys. 2002;43.

[13] Mostafazadeh A. Pseudo-Hermiticity versus PT -symmetry III: Equivalence of pseudo-Hermiticity and the presence of antilinear symmetries. J. Math. Phys. 2002;43.

[14] Mostafazadeh A. Exact PT -symmetry is equivalent to Hermiticity. J. Phys. A. 2003;36.

[15] Griffiths DJ. Introduction to Quantum Mechanics. Pearson Education Limited; 2013.

(29)

[16] Magnus W, Oberhettinger F, Soni RP. Formulas and theorems for the special func-tions of mathematical physics. Berlin; 1966.

[17] Harris GB, Arfken HJ, Weber FE. Mathematical Methods for Physicists. Academic Press; 2013.

A study of two PT -symmetric quantum mechanical systems