Acta Polytechnica Vol. 51 No. 4/2011

Perturbation Theory for P T-Symmetric Sinusoidal Optical Lattices at
the Symmetry-Breaking Threshold

H. F. Jones

Abstract

The P T symmetric potential V0[cos(2πx/a)+ iλsin(2πx/a)] has a completely real spectrum for λ ≤ 1, and begins to
develop complex eigenvalues for λ > 1. At the symmetry-breaking threshold λ = 1 some of the eigenvectors become
degenerate, giving rise to a Jordan-block structure for each degenerate eigenvector. In general this is expected to give
rise to a secular growth in the amplitude of the wave. However, it has been shown in a recent paper by Longhi, by
numerical simulation and by the use of perturbation theory, that for an initial wave packet this growth is suppressed,
giving instead a constant maximum amplitude. We revisit this problem by developing the perturbation theory further.
We verify that the results found by Longhi persist to second order, and with different input wave packets we are able to
see the seeds in perturbation theory of the phenomenon of birefringence first discovered by El-Ganainy et al.

Keywords: pseudo-Hermitian quantum mechanics, optical lattices, perturbation theory.

1 Introduction
The study of quantum mechanical Hamiltonians that
are P T -symmetric but not Hermitian [1–6] has re-
cently found an unexpected application in classical
optics [7–15], due to the fact that in the paraxial ap-
proximation the equation of propagation of an elec-
tromagnetic wave in a medium is formally identical
to the Schrödinger equation, but with different in-
terpretations for the symbols appearing therein. It
turns out that propagation through such a medium
exhibits many new and interesting properties, such as
power oscillations and birefringence. The equation of
propagation takes the form

i
∂ψ

∂z
= −

(
∂2

∂x2
+ V (x)

)
ψ, (1)

where ψ(x, z) represents the envelope function of the
amplitude of the electric field, z is a scaled propa-
gation distance, and V (x) is the optical potential,
proportional to the variation in the refractive index
of the material through which the wave is passing.
A complex V corresponds to a complex refractive
index, whose imaginary part represents either loss
or gain. In principle the loss and gain regions can
be carefully configured so that V is P T symmetric,
that is V ∗(x) = V (−x). There is also a non-linear
version of this equation, arising from sufficiently in-
tense beams, where there is an additional term pro-
portional to |ψ|2ψ. However, for the purposes of this
paper we shall limit ourselves to the linear case.

A model system exemplifying some of the novel
features of beam propagation in P T -symmetric opti-

cal lattices uses the sinusoidal potential

V = V0 [cos(2πx/a) + iλ sin(2πx/a)]

This model has been studied numerically and the-
oretically in Refs. [9, 12, 13]. The propagation in z
of the amplitude ψ(x, z) is governed by the analogue
Schrödinger equation (1), which for an eigenstate of
H, with eigenvalue β and z-dependence ψ ∝ e−iβz
reduces to the eigenvalue equation

− ψ′′ − V0 [cos(2πx/a) + iλ sin(2πx/a)] ψ = βψ .(2)

It turns out that these eigenvalues are real for λ ≤ 1,
which corresponds to unbroken P T symmetry, where
the eigenfunctions respect the (anti-linear) symmetry
of the Hamiltonian. Above λ = 1 complex eigenval-
ues begin to appear, and indeed above λ ≈ 1.776 87
all the eigenvalues are complex [16]. Clearly one
would expect oscillatory behaviour of the amplitude
below the threshold at λ = 1 and exponential be-
haviour above the threshold, but the precise form
of the evolution at λ = 1 is less obvious. At first
sight one would expect linear growth because of the
appearance of Jordan blocks associated with the de-
generate eigenvalues that merge at that value of λ,
but, as Longhi [12] has emphasized, this behaviour
can be significantly modified depending on the na-
ture of the initial wave packet.

In a previous paper [17] we approached this
problem by explicitly constructing the Bloch wave-
functions and the associated Jordan functions cor-
responding to the degenerate eigenvalues and then
using the method of stationary states to construct
the z-dependence. We found that the explicit lin-
ear dependence arising from the Jordan associated

21



Acta Polytechnica Vol. 51 No. 4/2011

functions is indeed cancelled by the combined contri-
butions from the non-degenerate wave-functions and
were able to understand how this cancellation came
about. In the present paper we approach the problem
from a different point of view by revisiting the com-
plementary perturbative calculation of Longhi [12].
In Section 2 we briefly recapitulate how the spec-
trum and eigenfunctions are calculated. Then in Sec-
tion 3, which forms the main body of the paper, we
give an explicit expression for the first-order contri-
bution and carry out the second-order calculation in
detail. This enables us to investigate the saturation
phenomenon for a variety of different inputs. Finally,
in Section 4 we give a brief discussion of our results.

2 Band structure at threshold
At the threshold λ = 1, the potential V in Eq. (2) be-
comes the complex exponential V = V0 exp(2iπx/a),
for which the Schrödinger equation is

− ψ′′ − V0 exp(2iπx/a)ψ = βψ. (3)

This is a form of the Bessel equation, as is seen
by the substitution y = y0 exp(iπx/a), where y0 =
(a/π)

√
V0, giving

y2
d2ψ

dy2
+ y

dψ

dy
− (y2 + q2)ψ = 0, (4)

�1 1
k.

a

Π

2

4

9

16

Β.
Π2

a2

Fig. 1: Band structure for λ = 1 in the reduced zone
scheme. The Bloch momentum k is plotted in units of
π/a and the eigengvalue β in units of (a/π)2

where q2 = β(a/π)2. Thus the spectrum is that of
a free massive particle, shown in the reduced zone
scheme in Fig. 1, and for q ≡ ka/π not an integer the

solutions ψk(x) = Iq (y) and ψ−k(x) = I−q (y) are lin-
early independent, and have exactly the correct peri-
odicity, ψk(x + a) = e

ikaψk(x), to be the Bloch wave-
functions. It is important to note, however, that be-
cause the original potential is P T -symmetric rather
than Hermitian, these functions are not orthogonal
in the usual sense, but rather with respect to the P T
inner product, namely∫

dxψ−k (x)ψk′ (x) = δkk′
∫

dxψ−k(x)ψk (x) (5)

However, for q = n, a non-zero integer, In(y)
and I−n(y) are no longer independent. In that case
the Bloch eigenfunctions do not form a complete set,
and we must search for other functions, still with
the same periodicity, to supplement them. These are
the Jordan associated functions, which we denote by
ϕk(x) ≡ χn(y). They may be defined as derivatives
of the eigenfunctions with respect to β, and satisfy
the generalized eigenvalue equation[

y2
d2

dy2
+ y

d
dy

− (y2 + n2)
]

χn(y) = In(y) (6)

The crucial feature of the Jordan functions is that be-
cause of this latter equation they naturally give rise
to linear growth in z, provided that they are excited:

e−iHz ϕr = e
−iβr ze−i(H−βr)zϕr

= e−iβr z(ϕr − izψr). (7)

However, as was found numerically in Ref. [12],
and explored further in Ref. [17], this natural lin-
ear growth may become saturated due to the contri-
butions of neighbouring Bloch functions, which are
closely correlated with those of the Jordan functions.

3 Perturbation Theory

The analysis of Ref. [17] approached the problem
from one point of view, in which the interplay be-
tween the contributions of the Bloch eigenfunctions
and the Jordan associated functions was made ex-
plicit. A complementary way of looking at things,
which does not separate these two contributions, is
to use the perturbative expansion, which instead em-
phasizes the contributions of the free propagation and
the corrections brought about by the potential. The
general framework for an expansion of ψ(x, z) in pow-

ers of V0, namely ψ(x, z) =
∞∑

r=0

V r0 ψr(x, z), has been

given in Ref. [12], along with an approximate form
of the first-order term ψ1(x, z) for the case q0 = −1
and w large. In this section we generalize this calcu-
lation by obtaining analytic expressions for both the
first- and second-order terms for general q0 and w.
Of course, this can only be used as a guide because

22



Acta Polytechnica Vol. 51 No. 4/2011

�400 �200 200 400
y

0.5

1

1.5

2
�a�

�400 �200 200 400
y

0.1

0.2

0.3

0.4

0.5

�b�

Fig. 2: Characteristic behaviour of error functions whose arguments are (a) real or (b) have a large imaginary part. In

(a) we plot erf(4+ y/w)+erf(4 − y/w) with w =80. In (b) we plot w e−w
2
|erf(4+ y/w − iw)+erf(4− y/w + iw)|

there is no guarantee that the expansion converges,
nor that the large-z behaviour of the complete am-
plitude can be extracted from the behaviour of the
truncated series. We will take as our input a Gaus-
sian profile of the form

ψ(x, 0) = f (x) ≡ e−(x/w)
2+ik0x, (8)

with offset k0 and width w. The zeroth-order term,
ψ0(x, z), is just the freely-propagating wave-packet

ψ0(x, z) =
w

√
(w2 + 4iz)

·

eik0(x−k0z)e−(x−2k0z)
2/(w2+4iz),

(9)

while the first-order term, ψ1(x, z), is given by

ψ1(x, z) = −i
∫

dkf̃ (k)ei(k+kB)x ·[
−i
∫ z
0

dy e−ik
2yei(k+kB)

2(y−z)
]

,
(10)

where

f̃ (k) =
w

2
√

π
e−(k−k0)

2w2/4 (11)

is the Fourier transform of f (x) of Eq. (8) and kB =
2π/a is the width of the first Brillouin zone. We can
reverse the order of integration in the expression for
ψ1, performing the (Gaussian) k integration first, to
obtain

ψ1(x, z) =

−i
w

√
(w2 + 4iz)

e
1
2 ikB x−

1
4 ik

2
B z−

1
4(k0+

1
2 kB)

2w2 × (12)∫ z
0

dy e−(2kB y−(x+kBz)+
1
2 iw

2(k0+
1
2 kB))

2
/(w2+4iz).

The y integration is then also a Gaussian integration
over a finite range, giving the result

ψ1(x, z) = −i
w
√

π

4kB
·

e
1
2
ikB x−14 ik

2
B z−

1
4
(k0+

1
2

kB)
2w2 (erf(η1) + erf(η2)) ,

(13)

where

η1 =
kB z + x − 12iw

2(k0 +
1
2kB )√

(w2 + 4iz)
(14)

η2 =
kB z − x + 12iw

2(k0 +
1
2kB )√

(w2 + 4iz)
. (15)

For the purposes of considering large w, it is conve-
nient to rewrite these in the form

η1 =
x − 2k0z√
(w2 + 4iz)

−

1
2
i(k0 +

1
2
kB )

√
(w2 + 4iz) (16)

η2 =
2z(kB + k0) − x√

(w2 + 4iz)
+

1
2
i(k0 +

1
2
kB )

√
(w2 + 4iz) . (17)

The case k0 = −
1
2
kB , i.e. q0 = −1, is clearly very

special, since in this case the second terms in the
contributions to η1 and η2 vanish, so that we get
the simple expressions η1 = (x + kB z)/

√
(w2 + 4iz)

and η2 = −(x − kBz)/
√

(w2 + 4iz). In that case, as
long as w2 � 4z the arguments may be treated as
effectively real, and each error function behaves like
a sign function of its argument (see Figure 2(a)), so
that the sum of the two behaves like the step func-
tion θ(kB z − |x|). This is the function Φ(x/(kB z)) of
Ref. [12]. In this case the qualitative features of the
perturbative calculation are in complete agreement
with the spreading of the wave-function in Figure 3
of that paper, and the saturation1 of ψmax. How-
ever, in this treatment there is of course no mention
of whether or not any Jordan functions are excited.

1ψ1(x, z) grows initially like z for small z.

23



Acta Polytechnica Vol. 51 No. 4/2011

-400 -200 200 400

0.2

0.4

0.6

0.8

1

-400 -200 200 400

0.002

0.004

0.006

0.008

0.01

0.012

0.014

Fig. 3: |ψ2(x, z)| versus x for z =50. (a) q0 = −1, (b) q0 =0. The parameters are: a =1, V0 =2 and w =6π

The same expressions in Eqs. (13) and (16) can
also be used for the cases q0 = 0 and q0 = 1 in the
limit of large w. In each case the arguments η1 and
η2 now have a large imaginary part. In that situation
the modulus of the erf has a narrow peak where the
real part vanishes (see Figure 2(b)). Thus the result
consists of two narrow rays, which are centered on
x = 2k0z and x = 2(kB + k0)z. Here we have the
seeds of the birefringence first observed in Ref. [7].
For the case k0 = q0 = 0, the two rays are centered
on x = 0 and x = 2kBz, while for the case q0 = 1, or

k0 =
1
2
kB , the two rays are centered on x = kBz and

x = 3kBz.
We now go on to second-order perturbation the-

ory to investigate the behaviour of ψ2(x, z), which is
given by [12]

ψ2(x, z) = −
∫

dkf̃ (k)ei(k+2kB)
∫ z
0

dη · (18)∫ z−η
0

dξ e−ik
2z−4ikb(k+kB)η+ikB(kB+2k)ξ

Again the k integration is a Gaussian, which leaves
finite-range Gaussian integrations over ξ and η. It
is convenient to change the integration variable ξ to
y ≡ 2η − ξ. The integration over y then yields the
expression

ψ2(x, z) = −
w

2kB

∫ z
0

dη

√
π

2
(erf(a) − erf(b)) · (19)

e−2ik
2
B ηe

3
2 ikB x−

1
4 ik

2
B t−

1
4 w

2Δ2,

where we have written k0 = −
(

1
2
kB + Δ

)
, and a

and b are given by

b =
4kBη − x − kBz − 12iw

2Δ
√

(w2 + 4iz)

a =
3kB (2η − z) − x − 12iw

2Δ
√

(w2 + 4iz)
. (20)

The final η integrations are then of the form

Iη ≡
∫

dη erf(c1η + c2)e
c3η

=
1
c3

[
ec3ηerf(c1η + c2) − (21)

ec3(c3−4c1c2)/(4c
2
1)erf (c1η + c2 − c3/(2c1))

]
.

Thus in principle ψ2 is expressible in terms of eight
error functions. However, it turns out in practice that
only six are involved.

In the case k0 = −kB/2, or q0 = −1, when
Δ = 0, the arguments of the error functions are such
as to give a plateau in |ψ2| between x = −3kBz
and x = −kBz, i.e. a widening of the beam to
the left of ψ0, and a much smaller peak, centered
around x = 3kBz, that is, a second weak beam to
the right. More importantly, the second-order con-
tribution again shows no sign of the linear growth2

in z näıvely expected from the excitation of Jordan
associated functions.

In the other case, q0 = 0, corresponding to
Δ = −kB/2, there are three peaks, centered around
x = 0, x = −2kBz and x = 4kBz, representing a
further splitting of the initial beam. Both cases are
illustrated in Figure 3.

4 Discussion
Of course one needs to treat the results of pertur-
bation theory with caution. In general terms we
have no proof that the perturbation series converges,
and in particular the asymptotic behaviour in z of
a few terms of the series does not necessarily give
the correct asymptotic behaviour of the entire sum.
Nonetheless this series does appear to give reliable re-
sults. For the parameters used by Longhi in Ref. [12],
with V0 = 0.2, first-order perturbation theory already
reproduces the numerical results very well, and the
second-order term gives only an extremely small cor-
rection.

2In fact ψ2(x, z) is proportional to z
2 for small z.

24



Acta Polytechnica Vol. 51 No. 4/2011

In all of our calculations we have not allowed z
to become too large, restricting it by the condition
z � w2/4, in which case saturation is an inbuilt fea-
ture of perturbation theory. In fact it was shown in
Ref. [17] using the method of stationary states that
if one goes to much larger values of z the amplitude
so calculated begins to grow again but this ultimate
resumption of linear growth is not physical, because
it corresponds to the situation where the beam has
widened beyond lateral limits of the optical lattice.

It is an elegant feature of the perturbative expan-
sion that the different types of possible behaviour of
the beam — spreading or splitting into two or more
beams — arise from very simple properties of the er-
ror function depending crucially on the offset k0. In
the first case the arguments are essentially real, and
the error functions behave like sign functions, while
in the second case there is a large imaginary part,
and the moduli of the error functions behave instead
like narrow peaks.

References

[1] Bender, C. M., Boettcher, S.: Phys. Rev. Lett.
80, 5243 (1998).

[2] Bender, C. M.: Contemp. Phys. 46, 277 (2005);
Rep. Prog. Phys. 70, 947 (2007).

[3] Mostafazadeh, A.: arXiv:0810.5643.

[4] Bender, C. M., Brody, D. C., Jones, H. F.:
Phys. Rev. Lett. 89, 270401 (2002); 92,
119902(E) (2004).

[5] Bender, C. M., Brody, D. C., Jones, H. F.:
Phys. Rev. D 70, 025001 (2004); 71, 049901(E)
(2005).

[6] Mostafazadeh, A.: J. Math. Phys. 43, 205
(2002); J. Phys. A 36, 7081 (2003).

[7] El-Ganainy, R. et al.: Optics Letters 32, 2632
(2007).

[8] Musslimani, Z. et al.: Phy. Rev. Lett. 100,
030402 (2008).

[9] Makris, K. et al.: Phy. Rev. Lett. 100, 103904
(2008).

[10] Klaiman, S., Günther, U., Moiseyev, N.:
Phys. Rev. Lett. 101 080402 (2008).

[11] Guo, A. et al.: Phy. Rev. Lett. 103, 093902
(2009).

[12] Longhi, S.: Phys. Rev. A 81, 022102 (2010).

[13] Makris, K. et al.: Phy. Rev. A 81, 063807
(2010).

[14] Rüter, C., et al.: Nature Physics 6, 192 (2010).

[15] Ramezani, H. et al.: arXiv:1005.5189 (2010).

[16] Midya, N., Roy, B., Choudhury, R.: Phys.
Lett. A 374, 2605 (2010).

[17] Graefe, E. M., Jones, H. F.: arXiv:1104.2838
(2011).

Hugh F. Jones
E-mail: h.f.jones@imperial.ac.uk
Physics Department
Imperial College
London SW7 2BZ, UK

25