Acta Polytechnica


Acta Polytechnica 53(3):259–267, 2013 © Czech Technical University in Prague, 2013
available online at http://ctn.cvut.cz/ap/

STATIONARY AND DYNAMICAL SOLUTIONS OF THE
GROSS-PITAEVSKII EQUATION FOR A BOSE-EINSTEIN

CONDENSATE IN A PT SYMMETRIC DOUBLE WELL

Holger Cartarius∗, Dennis Dast, Daniel Haag, Günter Wunner,
Rüdiger Eichler, Jorg Main

Institut für Theoretische Physik 1, Universität Stuttgart, Pfaffenwaldring 57, 70550 Stuttgart, Germany
∗ corresponding author: Holger.Cartarius@itp1.uni-stuttgart.de

Abstract. We investigate the Gross-Pitaevskii equation for a Bose-Einstein condensate in a PT
symmetric double-well potential by means of the time-dependent variational principle and numerically
exact solutions. A one-dimensional and a fully three-dimensional setup are used. Stationary states
are determined and the propagation of wave function is investigated using the time-dependent Gross-
Pitaevskii equation. Due to the nonlinearity of the Gross-Pitaevskii equation the potential depends
on the wave function and its solutions decide whether or not the Hamiltonian itself is PT symmetric.
Stationary solutions with real energy eigenvalues fulfilling exact PT symmetry are found as well as
PT broken eigenstates with complex energies. The latter describe decaying or growing probability
amplitudes and are not true stationary solutions of the time-dependent Gross-Pitaevskii equation.
However, they still provide qualitative information about the time evolution of the wave functions.

Keywords: Bose-Einstein condensates, PT symmetry, Gross-Pitaevskii equation, stationary states,
dynamics.

1. Introduction
Bose-Einstein condensates can, at extremely low tem-
peratures, be described by the Gross-Pitaevskii equa-
tion [1, 2], which reads in a particle number scaled
form and in appropriate units

iψ̇(x, t) =
(
−∆ + V (x) −g|ψ(x, t)|2

)
ψ(x, t). (1)

This is the Hartree approximation of the correspond-
ing many-particle equation for the dilute atomic gas,
where the assumption that all atoms are in the quan-
tum mechanical ground state is used. An external
potential, in which the atom cloud is trapped, is de-
scribed by V (x). Additionally, the atoms interact via
the short-range van der Waals force. In the dilute gas
it can be described sufficiently exact by an s-wave scat-
tering process, which leads in the mean-field approx-
imation to the nonlinear contribution −g|ψ(x, t)|2.
The strength g is determined by the s-wave scattering
length. Using magnetic fields acting on the hyperfine
levels of the atoms and exploiting Feshbach resonances,
the scattering length can be tuned and g is a true
parameter of the system, which can be varied.
Bose-Einstein condensates are a promising candi-

date for a first experimental realization of a special
class of non-Hermitian Hamiltonians in quantum me-
chanics, viz. systems which possess a PT symmetry
in spite of their non-Hermitian nature. As was shown
first by Bender and Boettcher, these Hamiltonians
exhibit remarkable properties [3] such as stationary
states with real eigenvalues in the presence of loss and
gain terms in the potential.

With the operators of spatial reflection P: x →−x,
p → −p, and of time reversal T : x → x, p → −p,
i →−i the PT symmetry of the Hamiltonian can be
expressed in terms of the commutator relation

[PT ,H] = 0. (2)

Since the kinetic energy term in the Hamiltonian is
always PT symmetric one obtains in a simple calcu-
lation from (2) the necessary condition

V ∗(−x) = V (x) (3)

for the potential.
Whereas the experimental realization of PT sym-

metry has successfully been achieved [4, 5] for optical
wave guides [6–12], a verification in a genuine quan-
tum system is still lacking. However, a proposal for
a Bose-Einstein condensate with a PT symmetric
external potential was given by Klaiman et al. [11],
who suggested a double-well setup, where atoms are
coherently incoupled in one well and outcoupled from
the other.
In this article, we will solve the Gross-Pitaevskii

equation of such a system and show that it indeed
supports PT symmetric solutions. Since the model
describes a true quantum system, our investigations
provide a good starting point for an experimental re-
alization of this quantum system with PT symmetry.
However, in the case of the nonlinear Gross-Pitaevskii
equation (1) one has to address one additional critical
question. The wave function has an influence on the
nonlinear Hamiltonian’s symmetry. Along with the

259

http://ctn.cvut.cz/ap/


Holger Cartarius et al. Acta Polytechnica

external potential the interaction term −g|ψ(x, t)|2
has to fulfill the condition (3). The Hamiltonian is
only PT symmetric if the square modulus of its so-
lution |ψ(x, t)|2 is a symmetric function of x. This
is always fulfilled for PT symmetric wave functions.
Thus, one may state that the nonlinear Hamiltonian
of the Gross-Pitaevskii equation preserves its own
symmetry in the case of exact PT symmetry. But
does this happen also with nonlinearity? Wave func-
tions describing PT broken eigenstates usually do not
possess a symmetric square modulus and destroy also
the PT symmetry of the Hamiltonian.
Previous studies of PT symmetric systems with

nonlinearity indicate that the nonlinearity does not
destroy the relevant effects known from linear PT
symmetric Hamiltonians. An investigation of a Bose-
Einstein condensate in a two-mode approximation us-
ing a non-Hermitian Bose-Hubbard model was made
by Graefe et al. [13–15]. Furthermore, the combina-
tion of PT symmetry and nonlinearity has been stud-
ied in quantum mechanical model potentials [10, 16],
in optics [17] or for a Bose-Einstein condensate with
an idealized double-δ trap [18], a system of which the
linear variant already attracted much interest [19–23].
To investigate Bose-Einstein condensates in a PT

symmetric double well we first introduce our numerical
approach in Section 2. Then we present and discuss
the numerical results for the energy eigenvalues, the
stationary wave functions and the time evolution of
non-stationary initial wave packets in Section 3. Con-
clusions are drawn in Section 4.

2. Numerical approach to
Bose-Einstein condensates in a
PT symmetric double well

2.1. Gross-Pitaevskii equation
We consider a Bose-Einstein condensate of particles
with mass m in a double-well setup described by the
external potential

V (x) =
m

2
ω2xx

2 +
m

2
ω2y,z(y

2 + z2) + v0e−σx
2

+ iΓxe−ρx
2
, (4)

where a three-dimensional harmonic trap is superim-
posed by a Gaussian barrier in the x direction (cf.
Figure 1). The trapping frequencies are ωx for the
direction of the double-well structure and ωy,z for
the two remaining directions. The barrier has its
maximum at x = 0, the height v0, and its width
is determined by σ. An outcoupling (incoupling) of
atoms is reflected by a negative (positive) imaginary
potential contribution in the left (right) well. Its
strength is determined by the gain-loss parameter Γ.
Since it affects the probability amplitude of the whole
condensate, the physical interpretation is a coherent
out-/incoupling. With this ansatz we are not con-
sidering individual atoms but the macroscopic wave
function of the condensed phase. The potential (4)

1

2

3

4

5

6

7

8

9

−6 −4 −2 0 2 4 6
−4

−3

−2

−1

0

1

2

3

4

R
e
V
(x
)

Im
V
(x
)

x

Re V (x)

Im V (x)

Figure 1. PT symmetric external potential. The
real part (solid line) defines the confinement of the
condensed atom cloud and the imaginary part (dashed
line) describes the gain-loss contributions due to the
coherent in- and outcoupling of atoms.

has been chosen such that it fulfills the condition (3),
i.e. its real part is a symmetric function of x, while
its imaginary part is antisymmetric. Thus, the linear
external potential (4) is PT symmetric.
Introducing the length scale a0 =

√
~/mωy,z de-

fined by the trap frequency in the direction perpendic-
ular to the double-well shape and the unit of energy
E0 = ~2/2ma20 the dimensionless potential reads

V (x) = ω2xx
2 + y2 + z2 + v0e−σx

2
+ iΓxe−ρx

2
, (5)

and the evolution of the condensate is described by
the Gross-Pitaevskii equation (1). With the chemical
potential µ and the usual separation ansatz ψ(x, t) =
φ(x)e−iµt in the units defined above one obtains the
time-independent Gross-Pitaevskii equation(

−∆ + V (x) −g|φ(x)|2
)
φ(x) = µφ(x). (6)

In our calculations we use the potential parameters
ωx = 0.5, v0 = 4, and σ = 0.5. The width parameter

ρ =
σ

2 ln(v0σ/ω2x)
(7)

of the gain-loss potential is chosen such that the ex-
trema of the real and imaginary potential parts coin-
cide, cf. Figure 1. If only the x direction is considered
and all y and z terms are removed from the potential
and the wave function, we can reduce the model to one
dimension that contains the relevant PT symmetric
information.

2.2. Two methods: Variational Gaussian
and numerically exact

We use two methods to solve the time-dependent
and time-independent Gross-Pitaevskii equations (1)
and (6). The Gaussian variational method has been
shown to provide highly precise solutions with low
numerical effort [24–26]. Since in our work it is applied

260



vol. 53 no. 3/2013 Stationary and Dynamical Solutions of the Gross-Pitaevskii Equation

for the first time to Bose-Einstein condensates in a
PT symmetric complex potential, we compare it to
numerically exact solutions of the Gross-Pitaevskii
equation in the one-dimensional case.
The idea of the Gaussian variational method con-

sists of the restriction to a Gaussian shaped wave
function

ψ(z, x) =
2∑
k=1

e−[A
k
x(x−q

k
x )

2+Aky,z (y
2+z2)]

× eip
k
x(x−q

k
x )−ϕ

k

(8)

described by a small set of variational parameters, viz.

z(t) =
{
Akx(t),A

k
y,z(t),q

k
x(t),p

k
x(t),ϕ

k
}
. (9)

In the case of the PT symmetric double-well setup it
is reasonable to start with two Gaussian wave func-
tions, where each of them is located in one of the wells.
The widths of the Gaussians are determined by the
complex parameters A1x, A2x, A1y,z and A2y,z. Since
the relevant dynamics only affects the x coordinate
and the trap was assumed to be symmetric in y and
z directions, we include in our ansatz the same sym-
metry for the Gaussian wave functions. The positions
and the corresponding momenta of both Gaussians
are determined by the real coordinates q1x, q2x, p1x and
p2x. The amplitudes and phases are introduced via the
complex variables ϕ1 and ϕ2. This leads in total to
16 real parameters completely defining the condensate
wave function. Reducing the model to one dimension
we end up with twelve real variables.

Certainly, the ansatz (8) cannot solve the time-
dependent Gross-Pitaevskii equation (1) exactly. One
way to find the “best” approximative solution with
a wave function restricted to the Gaussian form (8)
is the application of the McLachlan time-dependent
variational principle [27],

δI = δ
∥∥iχ(z(t), x) −Hψ(z(t), x)∥∥2 != 0. (10)

In this procedure the variation with respect to the
parameters in the wave function χ is performed such
that the functional I is minimized. In a second step
one sets ψ̇ ≡ χ and obtains the equations of motion,
which in our case are of the form

Ȧkx = −4i
(
(Akx)

2 + (Aky,z)
2) + iV k2;x, (11a)

Ȧky,z = −4i
(
(Akx)

2 + (Aky,z)
2) + iV k2;y,z, (11b)

q̇kx = 2p
k
x + s

k
x, (11c)

ṗkx = −Re v
k
1;x − 2 Im A

k
xs
k
x − 2 Re V

k
2;xq

k
x, (11d)

ϕ̇k = ivk0 + 2i(A
k
x + A

k
y,z) − i(p

k
x)

2

− ipkxs
k
x + iq

k
xv
k
1;x + iq

k
xV

k
2;xq

k
x, (11e)

with

skx =
1
2

(Re Akx)
−1(Im vk1;x + 2 Im V

k
2;xq

k
x). (11f)

The equations of motion (11a)-(11f) contain effec-
tive potential terms v = (v10, . . . ,v11;x, . . . ,V 12;x, . . . ),
which are obtained from a system of linear equations,
Kv = r, where the matrix K contains weighted over-
lap integrals of the Gaussians and the vector r consists
of weighted Gaussian averages of all potential terms
including the nonlinearity. The detailed form can
be found in [28]. The dynamics of the condensate
wave function is found by solving the ordinary dif-
ferential equations (11a)-(11f) with a Runge-Kutta
algorithm [29].
Stationary states or solutions of the time-indepen-

dent Gross-Pitaevskii equation are found by the re-
quirements Ȧkx = Ȧky,z = q̇kx = ṗkx = 0 (12 conditions
for real numbers), ϕ̇1 = ϕ̇2 (2 conditions). These
states have to be normalized, ‖ψ‖ = 1, which is im-
portant in the nonlinear Gross-Pitaevskii equation
and adds an additional constraint. Due to the arbi-
trary global phase one of the 16 Gaussian parameters
introduced above is free and 15 parameters must be
varied to fulfill the 15 conditions. This is done with
a 15-dimensional root search by applying a Powell
hybrid method [29]. If we consider a one-dimensional
condensate the root search reduces to 11 conditions,
which have to be fulfilled by 11 appropriately cho-
sen parameters. This small difference exemplifies the
high scalability of the variational Gaussian method.
An increase of the dimension or the flexibility of the
wave function leads only to a moderate increase of the
numerical effort.
The numerically exact integrations of the Gross-

Pitaevskii equations (1) and (6) presented in this
paper are carried out for the one-dimensional setup,
where the computational costs are reasonable. The
stationary wave functions are integrated outward from
x = 0 in positive and negative direction using a Runge-
Kutta algorithm with initial values Re ψ(0),ψ′(0) ∈ C,
and µ ∈ C. They have to be chosen such that the five
conditions ψ(∞) → 0, ψ(−∞) → 0, and ‖ψ‖ = 1
are satisfied. The conditions define square-integrable
and normalized wave functions, i.e. the stationary
states we are interested in. Note that the arbitrary
global phase has been exploited by setting Im ψ(0) = 0
for the initial values of the integration. For dynamical
calculations the split operator method is used, as
explained in [28].

3. Numerical solutions
3.1. Energy eigenvalues
Figure 2 shows a typical example for the eigenvalues
of the time-independent Gross-Pitaevskii equation (6).
First, one observes that the linear system (g = 0)
reveals the known features of complex Hamiltonians
with PT symmetry. Two real eigenvalue solutions
are found below a value ΓEP, where they merge in
an exceptional point. For larger values of Γ two com-
plex and complex conjugate solutions are found. The
agreement between the Gaussian approximation (solid

261



Holger Cartarius et al. Acta Polytechnica

2.34

2.36

2.38

2.40

2.42

2.44

2.46

2.48

2.50

2.52

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08

R
e
µ

Γ

(a)

−0.08
−0.06
−0.04
−0.02
0.00

0.02

0.04

0.06

0.08

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08

Im
µ

Γ

(b)
g = 0

g = 0.1

g = 0.2

g = 0.3

Figure 2. Eigenvalues of the time-independent Gross-
Pitaevskii equation (6) for different values of the non-
linearity g in dependence on the gain-loss parameter
Γ. With increasing g the real part of the energies
decreases. The Gaussian approximation (solid lines)
and the numerically exact solutions (dashed lines) are
shown. Vanishing imaginary parts are not plotted.
Two solutions with real eigenvalues are obtained up
to a value ΓEP ≈ 0.04, where they merge in an excep-
tional point. Two additional solutions with complex
eigenvalues are obtained, starting at a critical value
Γc, where Γc < ΓEP for g 6= 0.

lines) and the numerically exact solution (dashed lines)
is excellent.

An important result is the persistence of real eigen-
value solutions for nonvanishing values of the gain-loss
parameter Γ in the case g 6= 0. This demonstrates
that the Gross-Pitaevskii equation with a nonlinearity
in the potential still supports real eigenvalue solutions,
i.e. non-decaying states. There are, however, a few
crucial differences between the linear and the nonlin-
ear system. The complex eigenvalue solutions are now
born at a value Γc < ΓEP, whereas at the exceptional
point ΓEP only the real eigenvalue states vanish and
new complex solutions do not appear. The bifurcation
scenario of the linear system seems to be split between
the emergence of complex eigenvalue solutions at Γc
and the disappearance of the real eigenvalue states
at ΓEP 6= Γc. It is possible to explain this unusual
characteristics by the non-analyticity of the Gross-
Pitaevskii equation [28]. For sufficiently high values
of g the complex eigenvalues exist down to Γ → 0 and
only lose their imaginary contribution for Γ = 0.
The one-dimensional model already contains the

whole important potential shape for the analysis of
a PT symmetric condensate. Thus, one can expect
that it includes all important features in the eigen-
value structure. This is confirmed by the comparison

4.34

4.36

4.38

4.40

4.42

4.44

4.46

4.48

4.50

4.52

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08

R
e
µ

Γ

(a)

−0.08
−0.06
−0.04
−0.02
0.00

0.02

0.04

0.06

0.08

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08
Im

µ
Γ

(b)
g3D = 0
g3D = 0.2π
g3D = 0.4π
g3D = 0.6π

Figure 3. Real (a) and nonvanishing imaginary (b)
parts of the energy eigenvalues of the time-independent
Gross-Pitaevskii equation in the one- (dashed lines)
and three-dimensional (solid lines) models. The com-
parison shows that, by applying an appropriate rescal-
ing of the nonlinearity parameter g and an energy
shift for the energy contributions from the two addi-
tional directions, the one-dimensional model already
provides a very good quantitative description of the
fully three-dimensional setup. It is almost impossible
to identify the differences in the graph.

of the eigenvalues obtained with the one- (dashed)
and three-dimensional (solid) potentials in Figure 3.
The three-dimensional setup contains two additional
directions in which an external harmonic oscillator
potential is present. This leads to the assumption
that, in the ground state, the energy is shifted by a
value of ∆µ = 2 in the units introduced above. Fur-
thermore, the spatial extension in y and z directions
leads to additional energy contributions from the con-
tact interaction. The difference can be estimated from
the expectation value of the contact energy. We de-
mand that the expectation values of both models are
identical,∫

R3
dx dy dz g3D|ψ3D(x)|4

!=
∫
R

dxg1D|ψ1D(x)|4

(12)
and obtain the condition

g3D = 2πg1D (13)

for a rescaled g1D of the one-dimensional model
which leads to the same contact energy as a three-
dimensional wave function with g3D. In this cal-
culation we used a product ansatz ψ3D(x) ≈
ψ1D(x)ψ0(y)ψ0(z) of the wave function, where ψ0 de-
scribes the harmonic oscillator ground state [28]. This
assumption is only correct in the linear case without

262



vol. 53 no. 3/2013 Stationary and Dynamical Solutions of the Gross-Pitaevskii Equation

−0.3
−0.2
−0.1

0

0.1

0.2

0.3

0.4

0.5

−6 −4 −2 0 2 4 6

φ

x

(a)

−0.5
−0.4
−0.3
−0.2
−0.1

0
0.1
0.2
0.3
0.4
0.5

−6 −4 −2 0 2 4 6

φ

x

(b)

Re φ

Im φ

|φ|2

Figure 4. Wave functions of the ground (a) and ex-
cited (b) eigenstates with real eigenvalues in the case
g = 0.2 and Γ = 0.03. Only the Gaussian solutions are
shown since they are almost identical with the numer-
ically exact values. The square moduli of both wave
functions are symmetric functions of x, preserving the
PT symmetry of the nonlinear Hamiltonian.

contact interaction. The remarkable agreement of
both calculations in Figure 3 shows, however, that it
still approximates the nonlinear case very well even
for considerable nonlinearities g1D ≈ 0.3. Due to the
excellent quantitative agreement between the one- and
three-dimensional calculations we will only consider
the one-dimensional variant in the following sections.

3.2. Wave functions
We have already seen that real eigenvalue solutions
do appear in the nonlinear Gross-Pitaevskii equation.
But, we still do not know whether we have found
truly PT symmetric solutions. As explained in the
introduction, the Hamiltonian is only PT symmetric
if the square modulus of the eigenstate is a symmetric
function of x. This question will be answered in this
section.

Figure 4 shows the wave functions belonging to the
real eigenvalues for g = 0.2 and Γ = 0.03. Obviously,
the square modulus of both wave functions is a sym-
metric function of x. This preserves the Hamiltonian’s
PT symmetry, it is not destroyed by the nonlinearity.
Or in other words, the nonlinear Hamiltonian picks
as eigenstates wave functions which render itself PT
symmetric.
The situation is, however, completely different for

the wave functions of the states with complex eigen-
values. In linear PT symmetric models the complex
eigenvalue solutions belong to the PT broken case,
in which the wave functions do not reflect the PT

−0.2
−0.1

0

0.1

0.2

0.3

0.4

0.5

0.6

−6 −4 −2 0 2 4 6

φ

x

(a)

−0.2
−0.1

0

0.1

0.2

0.3

0.4

0.5

0.6

−6 −4 −2 0 2 4 6

φ
x

(b)
Re φ

Im φ

|φ|2

Figure 5. Wave functions of the complex eigenvalue
solutions with negative (a) and positive (b) imaginary
part in the case g = 0.2 and Γ = 0.03. Again only
the Gaussian wave functions are shown. The square
moduli are not symmetric functions of x and destroy
the Hamiltonian’s PT symmetry.

symmetry of the Hamiltonian. This is also the case for
the solutions of the Gross-Pitaevskii equation as can
be seen in Figure 5. The wave functions of the solu-
tions with complex eigenvalues are not PT symmetric.
Their square moduli are not symmetric functions of x,
and thus even the PT symmetry of the Hamiltonian
is destroyed, a circumstance which is not possible in a
linear quantum system. The wave function belonging
to the eigenvalue with negative imaginary part has a
higher amplitude in the left well with loss, whereas
the probability amplitude of the state with positive
imaginary part of the energy is shifted more to the
right well with gain.

A PT symmetric Bose-Einstein condensate will only
be observable if the eigenstates are stable with respect
to quantum fluctuations. We performed a linear sta-
bility analysis to check the behavior of the eigenstates
with real eigenvalues. A linearization of the equations
of motion (11a)–(11f) of the Gaussian parameters
around the stationary states and the solution of the
Bogoliubov-de Gennes equations for the numerically
exact wave functions has been done [28]. We found
that the excited state is stable for all combinations
of g and Γ as long as it exists. The ground state,
however, becomes unstable as soon as the PT bro-
ken branches with complex eigenvalues emerge in the
energy diagram of Figure 2.

263



Holger Cartarius et al. Acta Polytechnica

0 60 120 180 240 300

t

−4

−2

0

2

4

x

0

0.1

0.2

0.3

0.4

0.5

0.6
(a)

0 60 120 180 240 300

t

−4

−2

0

2

4

x

0

0.1

0.2

0.3

0.4
(b)

Figure 6. Evolution of the wave packet (14) in the
PT symmetric double well for nonvanishing nonlinear-
ity g = 0.2, ϕ = π/2, and the gain-loss contributions
Γ = 0 (a) and Γ = 0.02 (b).

3.3. Temporal evolution and the
significance of “stationary”
solutions with complex eigenvalues

We first want to investigate the temporal evolution
of condensate wave functions close to the stationary
real eigenvalue solutions for values of Γ below the
appearance of the PT broken states. Two examples
of numerically exact propagations using the split oper-
ator method are given in Figure 6 for the initial wave
packet

ψ(x,t = 0) =
1
√

2
(
φGS(x) + eiϕφES(x)

)
(14)

with the ground state φGS(x) and the excited state
φES(x). In this case one expects for linear systems
the behavior that an oscillation of the probability am-
plitude between the two wells sets in. The frequency
decreases with increasing Γ until the oscillation stops
at Γ = ΓEP [11]. This behavior is reproduced in Figure
6, i.e. the qualitative pattern of the motion is con-
served in the nonlinear Gross-Pitaevskii equation (1).
A quantitative analysis shows that the nonlinearity
leads to slightly higher oscillation frequencies.

However, with the emergence of the additional PT
broken states we observe a qualitative change between
the linear and nonlinear systems. As was mentioned
above, the ground state becomes unstable in this
regime. The relevance of this unstable character can
be seen in Figure 7a, where the initial wave packet
(14) is evolved for g = 0.2, Γ = 0.03, and the phase
ϕ = π. Already during the first oscillation a complete
deformation of the typical pattern is observed. Then
the oscillation stops and changes into an unrestricted
growth of the probability amplitude in the right well,

0 60 120

t

−4

−2

0

2

4

x

0

0.5

(a)

0 60 120 180 240 300

t

−4

−2

0

2

4

x

0

0.1

0.2

0.3
(b)

Figure 7. Evolution of the wave packet (14) in the
PT symmetric double well for nonvanishing nonlinear-
ity g = 0.2 and the gain-loss contributions Γ = 0.03
with phase ϕ = π (a) and Γ = 0.04 with ϕ = π/2
(b).

the wave function “explodes”. This behavior is not
surprising. One of the superimposed states in (14)
is unstable and the two PT broken solutions with
complex eigenvalues exist. One can expect that, in
the nonlinear system, there is a considerable overlap
of the time evolved wave function (14) with the PT
broken eigenstates, cf. also Figures 4 and 5. Then,
the eigenstate with positive imaginary part of the
energy will always dominate for long times since it
increases, and determines the further evolution. In a
realistic situation an infinite growth of the probability
amplitude will not appear. It has its origin in the PT
symmetric potential (4) and requires an infinite reser-
voir of atoms. At some point in time this description
will break down.

The instability does not necessarily lead immedi-
ately to a destruction of the original oscillating be-
havior, as is shown for the example with g = 0.2 and
Γ = 0.04 in Figure 7b, i.e. very close to ΓEP. Even
for a gain-loss influence larger than that chosen for
Figure 7a it was possible to find a stable oscillation
for the phase ϕ = π/2. The behavior known from
linear systems re-emerges. The probability amplitude
almost pulsates in both wells with a low frequency.
From the dynamical point of view the solutions of

the time-independent Gross-Pitaevskii equation (6)
with complex eigenvalues have to be considered with
some care. Strictly speaking, they lose their physical
relevance. Due to the decay or growth of the proba-
bility amplitude mentioned above the nonlinear term
−g|ψ|2 in the Hamiltonian becomes explicitly time
dependent. Thus, the states are not true stationary
solutions of the time-dependent Gross-Pitaevskii equa-
tion (1). They are, however, still useful to indicate the

264



vol. 53 no. 3/2013 Stationary and Dynamical Solutions of the Gross-Pitaevskii Equation

0
1
2
3
4
5
6
7
8
9
10

0 10 20 30 40 50

N
2

t

actual time developement

exp(−2 Im µt)

Figure 8. Temporal evolution of the norm of an
initial complex eigenvalue state with Im µ > 0. The
parameters are g = 0.2 and Γ = 0.03.

temporal evolution of the condensate. For example,
it was possible above to explain the explosion of the
condensate wave function in Figure 7a by the presence
of the growing PT broken state.
Furthermore, it is possible to show that the com-

plex eigenvalue solutions can still indicate the actual
temporal evolution of an initial wave packet. In par-
ticular, for small times the prediction from the “sta-
tionary” complex eigenvalue solutions for the norm
N2 =

∫
|ψ|2 dx describes very well the true onset of

the growth, as can be seen for an initial state with
Im µ > 0 in Figure 8.
To investigate the correspondence of the complex

eigenvalue solutions with the exact time integration
for longer times we introduce the difference

D =
1
N

(√∫
right well |ψ|

2 dx

−
√∫

left well |ψ|
2 dx

)
, (15)

of the wave function’s norm in the right and left wells.
It tells us how the probability density is distributed in
both wells. A positive value signalizes a higher proba-
bility density in the right well with gain, whereas a
negative value indicates a concentration of the con-
densate’s probability density in the left well with loss.

In Figure 9 we compare the norm difference D of the
correct temporal evolution with that of the complex
eigenvalue solution with positive imaginary part, i.e.
the growing state for the same parameters. Since the
norm of the time evolved wave changes due to the
gain-loss part of the potential we have to adapt the
effective g in the stationary Gross-Pitaevskii equation
(6),

g → gN2, (16)

such that the two norm differences are comparable.
In Figure 9a we plot D for the situation depicted in
Figure 7. With increasing time the state first decays
and the norm drops. At low values of the norm the
effective g (16) assumes values at which only the true
stationary states with real eigenvalues and D = 0
exist. The latter property can be seen in the figure.
However, at t ≈ 12 the overlap with the growing

−0.8
−0.6
−0.4
−0.2
0.0

0.2

0.4

0.6

0.8

1.0

0 20 40 60 80 100 120 140 160

D

t

(a)

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

0 10 20 30 40 50 60 70 80

D

t

(b)

actual time developement

adapted broken eigenstate

Figure 9. Comparison of the norm difference D of
the correct temporal evolution with that of complex
eigenvalue solutions. In (a) the same situation as in
Figure 7 is depicted, whereas (b) shows the evolution of
the complex eigenvalue state with positive imaginary
part. In both cases the norm difference follows the
growing complex eigenvalue state for long times. Γ
was chosen to be 0.03 and the initial g was 0.2.

eigenstate wins and the norm starts to grow. After
some oscillations of D, one observes for times t > 120
that the norm differences of the two calculations agree
more and more. This indicates that the wave function
initially prepared in a superposition of two eigenstates
evolves into the shape of the growing eigenstate for
long times. That is, the complex eigenvalue solution
with positive imaginary part still has a meaning for
long times although it is not a true stationary solu-
tion of the time-dependent Gross-Pitaevskii equation
(1). A similar behavior can be observed in Figure
9b, where the initial state was the growing complex
eigenvalue solution of (6). Since we start with the
increasing norm solution the agreement between both
calculations sets in earlier.

4. Conclusion
The work presented in this article shows that PT
symmetric eigenfunctions exist in nonlinear quan-
tum systems, in particular, in Bose-Einstein con-
densates described by the Gross-Pitaevskii equa-
tion. The solutions render the Hamiltonian it-
self PT symmetric. A comparison of numerically
exact calculations and a Gaussian approximation
demonstrated that the Gaussian ansatz provides
quantitatively well converged numerical results. A
fully three-dimensional calculation of the conden-
sate is possible, but is not necessary to describe
the effects appearing due to the non-Hermitian

265



Holger Cartarius et al. Acta Polytechnica

gain-loss potential. It is even possible to extract
quantitative values from a one-dimensional descrip-
tion of the condensate containing only the relevant
direction, in which the PT symmetric potential
acts.

Solutions of the time-independent Gross-Pitaevskii
equation with complex energy eigenvalues are found
as well, and belong to eigenstates with broken
PT symmetry, destroying the Hamiltonian’s sym-
metry. They have no direct physical meaning,
since they are not true stationary states of the
time-dependent Gross-Pitaevskii equation due to
their complex energy eigenvalues indicating a growth
or decay of the probability amplitude. However,
one can observe that they influence the ground
state. For nonvanishing nonlinearity g the com-
plex eigenvalue solutions bifurcate from the ground
state when a critical value of the gain-loss parame-
ter Γ is exceeded. At this point the ground state
becomes unstable, whereas the excited real eigen-
value solution is not affected by the appearance of
the new states and stays stable as long as it ex-
ists.
The time evolution of the condensate showed that

the eigensolutions of the time-independent Gross-
Pitaevskii equation with complex eigenvalues may
help to estimate the true temporal behavior of a con-
densate wave function. For small times, the imaginary
part of the energy correctly describes the onset of
the wave function’s growth or decay. For large times,
all initial wave functions tend to the state with pos-
itive imaginary part of the energy, which is located
predominantly in the well with gain.

Certainly, there is a number of questions which still
have to be answered. A better understanding of the
nonlinearity’s influence on the solutions is important.
Analytically solvable matrix models might help to get
a better insight. Furthermore, the relation of the three-
and one-dimensional results should be investigated for
higher values of the nonlinearity g.

A full analytical continuation of the non-analytical
Gross-Pitaevskii equation should help to understand
the change in the number of solutions observed
in the eigenvalue diagrams. In a proper exten-
sion, the critical values of Γ at which solutions ap-
pear or vanish should turn out to be bifurcation
points.

It would be desirable to understand how a coherent
in- or outcoupling of atoms can be understood on
a microscopic level. This will be important for a
realistic description of experimental situations. It is
also possible to extend the model such that it can be
understood as an embedding in a chain of potential
wells. Here the gain-loss contributions can result from
an effective description of a transport effect when
only two wells somewhere in the middle of the chain
are taken into account. Furthermore, it would be
interesting to see how a gain-loss potential interacts
with long-range inter-atomic interactions such as the

dipole-dipole interaction leading often to qualitatively
new effects [30].

References
[1] Gross, E. P. Structure of a Quantized Vortex in Boson
Systems. Nuovo Cimento 20, 454 (1961).

[2] Pitaevskii, L. P. Vortex Lines in an Imperfect Bose
Gas. Sov. Phys. JETP 13, 451 (1961).

[3] Bender, C. M. and Boettcher, S. Real spectra in
non-Hermitian Hamiltonians having PT symmetry.
Phys. Rev. Lett. 80, 5243 (1998).

[4] Guo, A., et al. Observation of PT -Symmetry
Breaking in Complex Optical Potentials. Phys. Rev.
Lett. 103, 093902 (2009).

[5] Rüter, C. E., et al. Observation of parity-time
symmetry in optics. Nat. Phys. 6, 192 (2010).

[6] Makris, K. G., et al. Beam Dynamics in PT Symmetric
Optical Lattices. Phys. Rev. Lett. 100, 103904 (2008).

[7] Makris, K. G., et al. PT -symmetric optical lattices.
Phys. Rev. A 81, 063807 (2010).

[8] Ruschhaupt, A., et al. Physical realization of PT
-symmetric potential scattering in a planar slab
waveguide. J. Phys. A 38, L171 (2005).

[9] El-Ganainy, R., et al. Theory of coupled optical
PT-symmetric structures. Opt. Lett. 32, 2632 (2007).

[10] Musslimani, Z., et al. Optical solitons in PT periodic
potentials. Phys. Rev. Lett. 100, 30402 (2008).

[11] Klaiman, S., et al. Visualization of Branch Points in
PT-Symmetric Waveguides. Phys. Rev. Lett. 101,
080402 (2008).

[12] Driben, R. and Malomed, B. A. Stability of solitons
in parity-time-symmetric couplers. Opt. Lett. 36, 4323
(2011).

[13] Graefe, E. M., et al. A non-Hermitian PT symmetric
Bose-Hubbard model: eigenvalue rings from unfolding
higher-order exceptional points. J. Phys. A 41, 255206
(2008).

[14] Graefe, E. M., et al. Mean-field dynamics of a
non-Hermitian Bose-Hubbard dimer. Phys. Rev. Lett.
101, 150408 (2008).

[15] Graefe, E.-M., et al. Quantum-classical
correspondence for a non-Hermitian Bose-Hubbard
dimer. Phys. Rev. A 82, 013629 (2010).

[16] Musslimani, Z. H., et al. Analytical solutions to a
class of nonlinear Schrödinger equations with PT -like
potentials. J. Phys. A 41, 244019 (2008).

[17] Ramezani, H., et al. Unidirectional nonlinear
PT -symmetric optical structures. Phys. Rev. A 82,
043803 (2010).

[18] Cartarius, H. and Wunner, G. Model of a
PT -symmetric Bose-Einstein condensate in a δ-function
double-well potential. Phys. Rev. A 86, 013612 (2012).

[19] Jakubský, V. and Znojil, M. An explicitly solvable
model of the spontaneous PT-symmetry breaking.
Czech. J. Phys. 55, 1113 (2005).

[20] Mostafazadeh, A. Delta-function potential with a
complex coupling. J. Phys. A 39, 13495 (2006).

266



vol. 53 no. 3/2013 Stationary and Dynamical Solutions of the Gross-Pitaevskii Equation

[21] Mostafazadeh, A. and Mehri-Dehnavi, H. Spectral
singularities, biorthonormal systems and a
two-parameter family of complex point interactions. J.
Phys. A 42, 125303 (2009).

[22] Mehri-Dehnavi, H., et al. Application of
pseudo-Hermitian quantum mechanics to a complex
scattering potential with point interactions. J. Phys. A
43, 145301 (2010).

[23] Jones, H. F. Interface between Hermitian and
non-Hermitian Hamiltonians in a model calculation.
Phys. Rev. D 78, 065032 (2008).

[24] Rau, S., et al. Pitchfork bifurcations in
blood-cell-shaped dipolar Bose-Einstein condensates.
Phys. Rev. A 81, 031605(R) (2010).

[25] Rau, S., et al. Variational methods with coupled
Gaussian functions for Bose-Einstein condensates with

long-range interactions. I. General concept. Phys. Rev.
A 82, 023610 (2010).

[26] Rau, S., et al. Variational methods with coupled
Gaussian functions for Bose-Einstein condensates with
long-range interactions. II. Applications. Phys. Rev. A
82, 023611 (2010).

[27] McLachlan, A. D. A variational solution of the time-
dependent Schrödinger equation. Mol. Phys. 8, 39 (1964).

[28] Dast, D., et al. A Bose-Einstein condensate in a PT
symmetric double well. Fortschr. Phys. 61, 124 (2013).

[29] Press, W. H., et al. Numerical Recipes in Fortran 77.
Cambridge University Press, New York, second edition
(1992).

[30] Lahaye, T., et al. The physics of dipolar bosonic
quantum gases. Rep. Prog. Phys. 72, 126401 (2009).

267


	Acta Polytechnica 53(3):259–267, 2013
	1 Introduction
	2 Numerical approach to Bose-Einstein condensates in a PT symmetric double well
	2.1 Gross-Pitaevskii equation
	2.2 Two methods: Variational Gaussian and numerically exact

	3 Numerical solutions
	3.1 Energy eigenvalues
	3.2 Wave functions
	3.3 Temporal evolution and the significance of ``stationary'' solutions with complex eigenvalues

	4 Conclusion
	References