## Realizing PT -symmetric non-Hermiticity with ultra-cold atoms and Hermitian multi-well potentials

Manuel Kreibich, Jörg Main, Holger Cartarius, and Günter Wunner Institut für Theoretische Physik 1, Universität Stuttgart, 70550 Stuttgart, Germany

We discuss the possibility of realizing a non-Hermitian , i. e. an open two-well system of ultra-cold atoms by enclosing it with additional time-dependent wells that serve as particle reservoirs. With the appropriate design of the additional wells PT -symmetric currents can be induced to and from the inner wells, which support stable solutions. We show that interaction in the mean-field limit does not destroy this property. As a first method we use a simplified variational ansatz leading to a discrete nonlinear Schrödinger equation. A more accurate and more general variational ansatz is then used to confirm the results.

## I. INTRODUCTION

Since the first realization of PT -symmetric gain and loss in optical wave guides [1] much effort has been made to realize analogous systems in various fields of physics, e. g. lasers [2], electronics [3], microwave cavities [4], and to make use of new effects arising from nonlinearity, viz. the Kerr nonlinearity in optical wave guides [5]. PT symmetry originates from quantum mechanics and stands for a combined action of parity and time-reversal. Despite the non-Hermiticity of PT -symmetric Hamiltonians, in a certain range of parameters entirely real eigenvalue spectra exist [6]. Due to a formal analogy between quantum mechanics and electromagnetism, the formulation of PT symmetry has spread to those systems, leading to the above mentioned realizations. However, the experimental verification in a genuine quantum system has not been achieved so far.

In a PT -symmetric system there is a balanced gain and loss of probability density (or electromagnetic field in analogue systems). According to a proposal in [7] such a quantum mechanical system could be realized with a Bose-Einstein condensate (BEC) in a double-well potential where particles are injected into one well and removed from the other. Indeed, it could be shown that the system is ideally suited for a first-experimental observation of PT symmetry in a quantum system [8]. There is progress in coupling two BECs and at the same time eject particles [9], but an actual realization of PT symmetry using a BEC is still missing.

On the other hand much theoretical work has been done on PT -symmetric BECs, including the proof of existence of stable states of interacting systems [10], a microscopic treatment based on the Bose-Hubbard model [11], and a thorough investigation on dynamics of stable and unstable regimes [12]. In contrast to injecting and removing particles from a double-well potential, one could think of a double-well included in a tilted optical lattice, with an incoming and outgoing transport of particles, so called Wannier-Stark systems [13, 14]. These incoming and outgoing particle currents could in principle serve as the necessary currents for realizing a PT -symmetric system.

In a recent Rapid Communication [15] we followed the

Figure 1. (a) Two-mode model with an imaginary potential V , where - in this case - particles are removed from the left well and injected into the right. (b) Such a system could be realized by coupling two additional wells to the system, where the tunneling currents to and from the outer wells are used as an implementation of the imaginary potential. The outer wells act as (limited) particle reservoirs and have to be chosen time-dependent.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0158v1-figure-6bcf9c88-18b7-4eb4-bdfc-4b3ac150c904.jpg)

simplest approach to this idea, namely to couple two additional wells to the double well potential (see Fig. 1), and use the arising currents as a theoretical realization of PT symmetry in a quantum mechanical system. We showed that for a specific time-dependent choice of the well depths and coupling strengths of the additional wells, the inner wells behave exactly as the wells of a PT -symmetric two-well system, thus serving as a possible experimental realization.

It is the purpose of this paper investigate more deeply our approach, especially for the case of interacting atoms, and to formulate a more accurate method to confirm the results. In Sec. II we derive the few-mode model or discrete nonlinear Schrödinger equation (DNLSE) from the Gross-Pitaevskii equation (GPE) using a simple variational ansatz, where each well is populated by a single time-independent Gaussian function. The results are applied in Sec. III to show that a Hermitian fourmode model can be chosen such that the middle wells behave exactly as the wells of the PT -symmetric twomode model. This is exemplified for some parameters. In Sec. IV we use an extended variational ansatz for the realization of PT symmetry. Still one Gaussian function per well is used, but each Gaussian function is fully dynamical in its parameters, thus giving a more accurate description of the system. The method is compared with

the simply Gaussian ansatz. In Sec. V we conclude our work.

## II. FROM THE GROSS-PITAEVSKII EQUATION TO A FEW-MODE MODEL

In this section we apply the steps that lead from the GPE, which describes the dynamics of a BEC at absolute zero temperature, to a DNLSE or few-mode-model. At first we presents our ansatz and calculate the necessary integrals. To derive the equations of motion (EOM) in a Schrödinger form we use the method of symmetric orthogonalization. In Sec. III the results are used and applied to different scenarios.

## A. Simplified variational ansatz

We start with the well-known GPE,

$$i \hbar { \partial } _ { t } \psi ( { \boldsymbol r }, t ) = \left [ - \frac { \hbar { ^ { 2 } } } { 2 m } \Delta + V ( { \boldsymbol r } ) + g \left | \psi ( { \boldsymbol r }, t ) \right | ^ { 2 } \right ] \psi ( { \boldsymbol r }, t ). \quad \text{Our } \colon _ { \substack { ( 1 ) \ \text{Schröd} \\ \text{modif} } }$$

The interaction strength is given by g = 4 π glyph[planckover2pi1] 2 Na/m with the s-wave scattering length a and particle number N . For the external multi-well potential V ( r ) we use a Gaussian profile for each well,

$$V ( r ) = \sum _ { k = 1 } ^ { N _ { w } } V ^ { k } \exp \left [ - \frac { 2 x ^ { 2 } } { w _ { x } ^ { 2 } } - \frac { 2 y ^ { 2 } } { w _ { y } ^ { 2 } } - \frac { 2 ( z - s _ { z } ^ { k } ) ^ { 2 } } { w _ { z } ^ { 2 } } \right ], \ \ ( 2 )$$

with a total number of N w wells. Each well k has the potential depth V k &lt; 0 and is displaced along the z -direction by s k z . Such a potential could, e. g., be created experimentally by the method presented in [16].

There is various work on how to transform the GPE (1) into a DNLSE [17, 18] by integrating out the spatial degrees of freedom. In general, an ansatz of the form

$$\psi ( r, t ) = \sum _ { k } d ^ { k } ( t ) g ^ { k } ( r ),$$

is used, where d k ( ) t is the (complex) amplitude and g k ( r ) the localized wave function in well k . Usually, a set of Wannier functions is applied [19]. In this work [Sec. III], we use a single Gaussian function for each well. With this ansatz we can derive analytical formulae for the matrix elements of the DNLSE or few-mode model, and thus open an easy way to actually compute the localized wave functions by a simple energy minimization process. In Sec. IV we extend the variational ansatz to obtain more accurate results.

The simplified variational ansatz is given by the superposition of N G Gaussian functions, one for each well, N w = N G , with

$$g ^ { k } ( r ) = e ^ { - A _ { x } ^ { k } x ^ { 2 } - A _ { y } ^ { k } y ^ { 2 } - A _ { z } ^ { k } ( z - q _ { z } ^ { k } ) ^ { 2 } }. \quad \ ( 4 )$$

The parameters A ,A ,A k x k y k z ∈ ❘ describe the width of each Gaussian function in each direction, and q k z ∈ ❘ the displacement along the z -direction. We treat these parameters in this section as time-independent. The timedependence of the wave function comes from the amplitude d k ( ) t in Eq. (3).

We insert ansatz (3) into the GPE i glyph[planckover2pi1] ∂ ψ t = ˆ Hψ , multiply from the left with ( g l ) ∗ , and integrate over ❘ 3 . Then, we are left with

$$\mathrm i \hbar { K } \dot { d } = \mathrm H d.$$

Here, the amplitudes d k are written as a vector d = ( d , . . . , d 1 N G ) T ∈ ❈ N G . The matrix K describes the overlap, and H the Hamiltonian matrix elements. They are given by the integrals (see Sec. II B for their explicit calculation)

$$K _ { l k } = \left \langle g ^ { l } | g ^ { k } \right \rangle,$$

$$H _ { l k } = \left \langle g ^ { l } \right | \hat { H } \left | g ^ { k } \right \rangle. \text{ \quad \ \ } ( 6 b )$$

Our aim is to derive EOM for the amplitudes in a Schrödinger form, i. e. i glyph[planckover2pi1] ˙ d eff = H eff d eff , with possibly modified amplitudes and matrix elements. Equation (5) would be of that form if K was proportional to the identity matrix. Since the Gaussian functions g k are nonorthogonal this is not the case. For that reason, we perform a symmetric orthogonalization in Sec. II C.

## B. Evaluation of the integrals

The occurring integrals in Eqs. (6) are solely based on the Gaussian integral

$$\intertext { a d } \int _ { \substack { \infty \\ \mathcal { d } = - \infty } } ^ { \infty } \int _ { \substack { \mathcal { d } x \, x ^ { n } \, e ^ { - A x ^ { 2 } + p x } } } & = \frac { \partial ^ { n } } { \partial p ^ { n } } \int _ { \substack { \mathcal { d } x \, e ^ { - A x ^ { 2 } + p x } } } \\ & = \frac { \partial ^ { n } } { \partial p ^ { n } } \sqrt { \frac { \pi } { A } } \, e ^ { p ^ { 2 } / 4 A }, \quad \text{Re } A > 0. \\ \intertext { a d } \int _ { \substack { \mathcal { d } x \, - \infty \\ \mathcal { d } = - \infty } } \cdots \cdots \cdots \cdots$$

For convenience we define the abbreviations

$$A _ { \alpha } ^ { k l } = A _ { \alpha } ^ { k } + ( A _ { \alpha } ^ { l } ) ^ { * }, \quad \alpha = x, y, z, \quad ( 8 \text{a} )$$

$$\kappa _ { \alpha } ^ { k l } = A _ { \alpha } ^ { k } ( A _ { \alpha } ^ { l } ) ^ { * } / A _ { \alpha } ^ { k l },$$

$$\beta _ { \alpha } ^ { k l } = \sqrt { \frac { A _ { \alpha } ^ { k l } w _ { \alpha } ^ { 2 } } { A _ { \alpha } ^ { k l } w _ { \alpha } ^ { 2 } + 2 } }, \quad \quad \quad ( 8 c )$$

$$c ^ { k l } = \, e ^ { - \kappa _ { z } ^ { k l } ( q _ { z } ^ { k } - q _ { z } ^ { l } ) ^ { 2 } },$$

The integrals of K can now easily be calculated with Eq. (7) for the case n = 0 , the result is

$$K _ { l k } = \sqrt { \frac { \pi } { A _ { x } ^ { k l } } } \sqrt { \frac { \pi } { A _ { y } ^ { k l } } } \sqrt { \frac { \pi } { A _ { z } ^ { k l } } } c ^ { k l }.$$

For the calculation of the Hamiltonian matrix H , we separate it into a kinetic term T , an external potential term V , and an interaction term W .

For the kinetic term we need integrals of the form shown in Eq. (7) with n = 2 and we obtain

$$T _ { l k } = K _ { l k } \left [ \kappa _ { x } ^ { k l } + \kappa _ { y } ^ { k l } + \kappa _ { z } ^ { k l } - 2 ( \kappa _ { z } ^ { k l } ) ^ { 2 } ( q _ { z } ^ { k } - q _ { z } ^ { l } ) ^ { 2 } \right ]. \ ( 1 0 ) \quad \text{with } i \,.$$

To calculate the integrals necessary for the external potential, the parameters A and p in Eq. (7) are shifted by the parameters of the external potential. The result is

$$& V _ { l k } = K _ { l k } \beta _ { x } ^ { k l } \beta _ { y } ^ { k l } \beta _ { z } ^ { k l } \sum _ { m } V ^ { m } \\ & \times \exp \left \{ - \frac { 2 \left [ A _ { z } ^ { k } ( s _ { z } ^ { m } - q _ { z } ^ { k } ) + ( A _ { z } ^ { l } ) ^ { * } ( s _ { z } ^ { m } - q _ { z } ^ { l } ) \right ] ^ { 2 } } { A _ { z } ^ { k l } ( A _ { z } ^ { k l } w _ { z } ^ { 2 } + 2 ) } \right \}. \quad ( 1 1 ) \quad \text{D.} \quad \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{                                                                                                                                                                                                        
d _ { \quad } } \quad \quad \quad$$

Since the operator of the interaction potential W is nonlinear, its matrix elements depend on the amplitudes d k . With the standard Gaussian integral we obtain

$$W _ { l k } = \sum _ { i, j } \tilde { W } _ { l k j i } ( d ^ { j } ) ^ { * } d ^ { i }, \text{ \quad \ \ } ( 1 2 ) \text{ \ \ } \text{ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
 \text{$\text{$\text{no}$} } }$$

where we defined the four-rank tensor

$$& \text{where we defined the four-rank tensor} \\ & \bar { W } _ { l k j i } = \frac { 4 \pi \hslash ^ { 2 } N a } { m } \sqrt { \frac { \pi } { A _ { z } ^ { i j k l } } } \sqrt { \frac { \pi } { A _ { z } ^ { i j k l } } } \sqrt { \frac { \pi } { A _ { z } ^ { i j k l } } } \\ & \quad \times \exp \left [ - \frac { A _ { z } ^ { i } ( A _ { z } ^ { j } ) ^ { * } ( q _ { z } ^ { i } - q _ { z } ^ { j } ) ^ { 2 } + A _ { z } ^ { i } ( A _ { z } ^ { l } ) ^ { * } ( q _ { z } ^ { i } - q _ { z } ^ { l } ) ^ { 2 } } { A _ { z } ^ { i j k l } } \right ] \\ & \quad \times \exp \left [ - \frac { A _ { z } ^ { k } ( A _ { z } ^ { j } ) ^ { * } ( q _ { z } ^ { k } - q _ { z } ^ { j } ) ^ { 2 } + A _ { z } ^ { k } ( A _ { z } ^ { l } ) ^ { * } ( q _ { z } ^ { k } - q _ { z } ^ { l } ) ^ { 2 } } { A _ { z } ^ { i j k l } } \right ] \\ & \quad \times \exp \left [ - \frac { A _ { z } ^ { i } A _ { z } ^ { k } ( q _ { z } ^ { i } - q _ { z } ^ { k } ) ^ { 2 } + ( A _ { z } ^ { j } ) ^ { * } ( A _ { z } ^ { l } ) ^ { * } ( q _ { z } ^ { j } - q _ { z } ^ { l } ) ^ { 2 } } { A _ { z } ^ { i j k l } } \right ] \\ & \quad \text{and} \\ & \text{and the new abbreviations}$$

and the new abbreviations

$$A _ { \alpha } ^ { i j k l } = A _ { \alpha } ^ { i j } + A _ { \alpha } ^ { k l }, \quad \alpha = x, y, z. \quad \ ( 1 4 ) \quad \text{whe}$$

## C. Symmetric orthogonalization

We calculated all necessary integrals for setting up the EOMs (5), which are not of the form of a Schrödinger equation (cf. Sec. II A). To achieve this, we now transform the EOM into a Schrödinger form with the method of symmetric orthogonalization [20]: Since K is Hermitian there exists a unitary transformation U such that D = UKU † is diagonal with real entries. Then we can construct the Hermitian matrix

$$\mathbf x = \mathbf U ^ { \dagger } \mathbf D ^ { - 1 / 2 } \mathbf U. \quad \quad \ \ ( 1 5 )$$

This matrix has the properties that XKX = ✶ , and furthermore, if H is Hermitian, so is XHX . Now, we can use this matrix to transform Eq. (5) to

$$i \hbar { ( } X K X ) ( X ^ { - 1 } \dot { d } ) = ( X H X ) ( X ^ { - 1 } d ). \quad ( 1 6 )$$

With the definitions d eff = X -1 d and H eff = XHX , and the above property of X we arrive at the EOM for d eff in Schrödinger form

$$i \hbar { d } _ { \text{eff} } = \mathbf H _ { \text{eff} } d _ { \text{eff} }$$

with a Hermitian Hamiltonian matrix H eff (if H is Hermitian). We note, that the normalization condition for the wave function ψ transforms as

$$\text{by} \quad \ 1 & \stackrel {! } { = } \int \text{d} ^ { 3 } r \, \left | \psi \right | ^ { 2 } = d ^ { \dagger } \text{K} d = ( X ^ { - 1 } d ) ^ { \dagger } ( \text{XXKX} ) ( X ^ { - 1 } d ) \\ & \quad = d ^ { \dagger } _ { \text{eff} } d _ { \text{eff} }.$$

## D. Application of symmetric orthogonalization and nearest-neighbor approximation

By means of symmetric orthogonalization we can transform the EOM (5) into a Schrödinger form since we know K analytically. Without any further approximations the analytical expressions can become quite complicated. At this point, we introduce as usual the nearest neighbor approximation to obtain analytical and simple approximate expressions for X and H eff .

The structure of the matrix elements K lk [see Eq. (9)] allows for a natural approximation since the function c kl drops exponentially with increasing distance of wells | k - | l . We say a term is of order n when | k - | l = n . As an approximation we only consider terms up to order n = 1 . Hence, we have K = K (0) + K (1) with

$$K _ { l k } ^ { ( 0 ) } = \left ( \frac { \pi } { 2 } \right ) ^ { 3 / 2 } \frac { \delta _ { k l } } { \sqrt { A _ { x, R } ^ { k } A _ { y, R } ^ { k } A _ { z, R } ^ { k } } } } \quad \ \ ( 1 9 )$$

and

$$K _ { l k } ^ { ( 1 ) } = \sqrt { \frac { \pi } { A _ { x } ^ { k l } } } \sqrt { \frac { \pi } { A _ { y } ^ { k l } } } \sqrt { \frac { \pi } { A _ { z } ^ { k l } } } c ^ { k l } \left ( \delta _ { k, l + 1 } + \delta _ { k + 1, l } \right ), \ \ ( 2 0 )$$

where A k α,R = Re A k α , and - for later use A k α,I = Im A k α . Using this expansion, we calculate the zeroth order of X by the requirement X (0) K (0) X (0) = ✶ . We obtain

$$X _ { l k } ^ { ( 0 ) } = \left ( \frac { 2 } { \pi } \right ) ^ { 3 / 4 } \, \delta \sqrt { A _ { x, R } ^ { k } A _ { y, R } ^ { k } A _ { z, R } ^ { k } } \delta _ { k l }. \quad ( 2 1 )$$

To calculate the first order we start with

$$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{	\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{ $\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{#\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{s\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{
$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text
$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\test{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\test$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text	$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text`$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text}$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text'$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text<$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$
$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$$$

ignore second order terms, insert the known quantities, and solve for X (1) , which yields

$$\text{---- ---- ---- ---- ---- ---- } \,, \, \cdot \, \dots \, \cdot \, \dots \, \cdot \, \cdot \, \dots \\ \text{fural} \,, & X ^ { ( 1 ) } _ { l k } = & - \left ( \frac { 8 } { \pi } \right ) ^ { 3 / 4 } \frac { c ^ { k l } } { \sqrt { A ^ { k l } _ { x } } \sqrt { A ^ { k l } _ { y } } \sqrt { A ^ { k l } _ { z } } } \left ( \delta _ { k, l + 1 } + \delta _ { k + 1, l } \right ) \\ \text{can} \, & \times \frac { \sqrt { A ^ { k } _ { x, R } A ^ { k } _ { y, R } A ^ { k } _ { z, R } } ^ { k } A ^ { l } _ { x, R } A ^ { l } _ { y, R } A ^ { l } _ { z, R } } { \sqrt { A ^ { k } _ { x, R } A ^ { k } _ { y, R } A ^ { k } _ { z, R } } + \sqrt { A ^ { l } _ { x, R } A ^ { l } _ { y, R } A ^ { l } _ { z, R } } }. \quad ( 2 3 )$$

Now we are able to calculate analytical approximations for the transformed Hamiltonian H eff . We observe that for the linear parts, i. e. the kinetic and external potential part, the matrix elements of H are proportional to K , so we can write H lin lk = K h lk lin lk [cf. Eqs. (10) and (11)]. Thus, the orders of H lin are given by the orders of K . Therefore, for the zeroth order transformed Hamiltonian we obtain

$$H _ { \text{eff}, l k } ^ { \text{lin}, ( 0 ) } = \sum _ { m, n } X _ { l m } ^ { ( 0 ) } K _ { m n } ^ { ( 0 ) } X _ { n k } ^ { ( 0 ) } h _ { m n } ^ { \text{lin} } = h _ { k k } ^ { \text{lin} } \delta _ { k l }. \quad ( 2 4 )$$

The zeroth order contributes to the linear diagonal elements of H lin eff , so they can be identified with the onsite energy E k = h lin kk . With Eqs. (10) and (11) we can write

$$& E _ { k } = \frac { \hslash ^ { 2 } } { 2 m } \left ( A _ { x, R } ^ { k } + A _ { y, R } ^ { k } + A _ { z, R } ^ { k } \right ) \\ & \quad + V ^ { k } \beta _ { x } ^ { k k } \beta _ { y } ^ { k k } \beta _ { z } ^ { k k } \exp \left [ - 2 ( \beta _ { z } ^ { k k } ) ^ { 2 } ( s _ { z } ^ { k } - q _ { z } ^ { k } ) ^ { 2 } / w _ { z } ^ { 2 } \right ], \ \ ( 2 5 ) \quad \text{These}$$

where we considered only the term with m = k in the sum in Eq. (11) to be consistent with zeroth order.

The first order terms of H lin eff are given by the first order terms of the expression

$$( X ^ { ( 0 ) } + X ^ { ( 1 ) } ) ( \mathbf H ^ { \mathrm i \mathrm n, ( 0 ) } + \mathbf H ^ { \mathrm i \mathrm n, ( 1 ) } ) ( X ^ { ( 0 ) } + X ^ { ( 1 ) } ), \ \ ( 2 6 ) \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

which gives, after inserting all known quantities,

$$H _ { \text{eff}, l k } ^ { \text{lin}, ( 1 ) } = & - 2 \sqrt { 2 } \frac { \sqrt { A _ { x } ^ { k }, R } ^ { k } A _ { y, R } ^ { k } A _ { z, R } ^ { l } A _ { x, R } ^ { l } A _ { y, R } ^ { l } A _ { z, R } ^ { l } } { \sqrt { A _ { x, R } ^ { k } } A _ { y, R } ^ { k } A _ { z, R } ^ { k } } + \sqrt { A _ { x, R } ^ { l } A _ { y, R } ^ { l } A _ { z, R } ^ { l } } } { \sqrt { A _ { x, R } ^ { k } } A _ { y, R } ^ { k } A _ { z, R } ^ { k } } + \frac { \text{$_{l}$} } { \sqrt { A _ { x, R } ^ { k } } A _ { y, R } ^ { k } A _ { z, R } ^ { l } } \right ) } \quad \text{condit} \\ & \times \left ( \frac { E _ { k } - h _ { l k } } { \sqrt { A _ { x, R } ^ { k } } A _ { y, R } ^ { k } A _ { z, R } ^ { k } } + \frac { E _ { l } - h _ { l k } } { \sqrt { A _ { x, R } ^ { l } } A _ { y, R } ^ { l } A _ { z, R } ^ { l } } \right ) } \right ) } { \sqrt { A _ { x } ^ { k l } } \sqrt { A _ { y } ^ { k l } } \sqrt { A _ { z } ^ { k l } } } \left ( \delta _ { k, l + 1 } + \delta _ { k + 1, l } \right ). \quad ( 2 7 ) \\ \text{These elements can be identified as the tunneling ele-}$$

These elements can be identified as the tunneling elements of the few-mode model, i. e. J lk = -H lin , (1) eff ,lk , since they are the super- and sub-diagonal entries. The quantity h lin lk is given by Eqs. (10) and (11), where in v lk we only consider nearest-neighbor contributions,

$$v _ { l k } & = \beta _ { x } ^ { k l } \beta _ { y } ^ { k l } \beta _ { z } ^ { k l } \\ & \times \left ( V ^ { k } \exp \left \{ - \frac { 2 \left [ A _ { z } ^ { k } ( s _ { z } ^ { k } - q _ { z } ^ { k } ) + ( A _ { z } ^ { l } ) ^ { * } ( s _ { z } ^ { k } - q _ { z } ^ { l } ) \right ] ^ { 2 } } { A _ { z } ^ { k l } ( A _ { z } ^ { k l } w _ { z } ^ { 2 } + 2 ) } \right \} \quad \text{sister} \quad \text{lowin} \quad \text{discre} \\ & + V ^ { l } \exp \left \{ - \frac { 2 \left [ A _ { z } ^ { k } ( s _ { z } ^ { l } - q _ { z } ^ { k } ) + ( A _ { z } ^ { l } ) ^ { * } ( s _ { z } ^ { l } - q _ { z } ^ { l } ) \right ] ^ { 2 } } { A _ { z } ^ { k l } ( A _ { z } ^ { k l } w _ { z } ^ { 2 } + 2 ) } \right \} \right ). \\ & \quad \infty \, \text{for two calculated the matrix alternate of the} \quad \text{A}.$$

So far we have calculated the matrix elements of the linear transformed Hamiltonian H lin eff . We now turn our attention to the interaction part. Since the operator W is nonlinear and depends on the amplitudes d k , they must be transformed, too. For the transformed operator, we obtain

$$\intertext { a }, \sigma _ { 1 } ^ { \infty } \intertext { tial } \ W _ { \text{eff}, l k } & = \sum _ { i, j, m, n, p, q } X _ { l m } \tilde { W } _ { m n j i } X _ { n k } ( X _ { j p } d _ { \text{eff} } ^ { p } ) ^ { * } ( X _ { i q } d _ { \text{eff} } ^ { q } ) \\ \intertext { K } _ { \text{ian} } & = \sum _ { p, q } \underbrace { \sum _ { i, j, m, n } X _ { l m } X _ { p j } \tilde { W } _ { m n j i } X _ { n k } X _ { i q } } _ { \equiv \tilde { W } _ { \text{eff}, l k p q } } ( d _ { \text{eff} } ^ { p } ) ^ { * } d _ { \text{eff} } ^ { q }. \\ & \intertext { r } _ { \text{ r } } \cdots \cdots \cdots \cdots \cdots \cdots.$$

For the interaction term, we only consider onsiteelements, i. e. elements with i = j = k = l in ˜ W lkji , since mainly these terms contribute to the interaction part. For the transformed four-rank tensor, we then obtain

$$\tilde { W } _ { \text{eff}, k k k k } = \frac { 4 \hbar { ^ { 2 } } N a } { \sqrt { \pi } m } \sqrt { A _ { x, R } ^ { k } A _ { y, R } ^ { k } A _ { z, R } ^ { k } }. \quad ( 3 0 )$$

These elements are the nonlinear coupling constants, c k = ˜ W eff ,kkkk . The whole action of the transformed Hamiltonian then reads

$$W _ { \text{eff}, l k } = c _ { k } | d ^ { k } | ^ { 2 } \delta _ { k l }. \text{ \quad \ \ } ( 3 1 )$$

We have all necessary matrix-elements in nearestneighbor approximation. Using this knowledge we can calculate these elements from the Gaussian variational approach, or - vice verse - determine the parameters of the realistic potential (2) by knowing the matrix elements. As an application this is done in Sec. III. The parameters A k α and q k z can be computed by minimizing the mean-field energy of the system for a given initial condition with

$$\int _ { \substack { E _ { \text{mf} } = \sum _ { k, l } ( d ^ { l } ) ^ { * } \left ( T _ { l k } + V _ { l k } \right ) d ^ { k } \\ + \frac { 1 } { 2 } \sum _ { i, j, k, l } ( d ^ { l } ) ^ { * } d ^ { k } W _ { l k j i } ( d ^ { j } ) ^ { * } d ^ { i }. \quad ( 3 2 ) \\ \text{el} \text{-} }$$

For this minimization, all parameters A k α , q k z and d k have to be varied with the norm constraint (18), which is computationally much cheaper than a grid calculation. The results of this section could also be used to calculate the parameters of the Bose-Hubbard model.

Now that we have all matrix elements, we can begin our actual investigation of a realization of a non-Hermitian system, first by means of few-mode models. In the following we will adopt the usual notation and write for the discrete wave function ψ k instead of d k .

## III. FEW-MODE MODELS

## A. Equivalence of the two- and four-mode models

We now return to the idea mentioned in the introduction and sketched in Fig. 1, i. e. the embedding of a

non-Hermitian double-well with a closed four-well structure. Before analyzing the possibility of realizing a nonHermitian system with a Hermitian model, we discuss the basic properties of the non-Hermitian two-mode model, which is given by

$$\mathbf H _ { 2 } = \begin{pmatrix} E _ { 1 } + c _ { 1 } \left | \psi _ { 1 } \right | ^ { 2 } & - J _ { 1 2 } \\ - J _ { 1 2 } & E _ { 2 } + c _ { 2 } \left | \psi _ { 2 } \right | ^ { 2 } \end{pmatrix} \quad ( 3 3 ) \quad \ \ \ \ a _ { n } ^ { \prime }$$

with complex onsite energies E k ∈ ❈ , which make the Hamiltonian non-Hermitian. This system has intensively been investigated in [21]. The time derivatives of the observables, that is particle number n k = ψ ψ ∗ k k and particle current j kl = i J kl ( ψ ψ k ∗ l -ψ ψ ∗ k l ) , can be easily calculated, which gives

$$\partial _ { t } n _ { 1 } = - j _ { 1 2 } + 2 n _ { 1 } \, \text{Im} \, E _ { 1 }, & & ( 3 4 a ) & &.$$

$$\partial _ { t } n _ { 2 } = j _ { 1 2 } + 2 n _ { 2 } \, \text{Im} \, E _ { 2 }, & & ( 3 4 b ) \quad \text{Com} _ { \dots \, 1 \,. }$$

$$\partial _ { t } j _ { 1 2 } & = 2 J _ { 1 2 } ^ { 2 } ( n _ { 1 } - n _ { 2 } ) + ( \text{Im} \, E _ { 1 } + \text{Im} \, E _ { 2 } ) j _ { 1 2 } \\ & \quad + J _ { 1 2 } ( \text{Re} \, E _ { 1 } - \text{Re} \, E _ { 2 } + c _ { 1 } n _ { 1 } - c _ { 2 } n _ { 2 } ) C _ { 1 2 }, \ ( 3 4 c ) \\ \cdot \cdot \cdot \cdot \cdot \cdot$$

$$( 3 4 c )$$

$$\partial _ { t } C _ { 1 2 } & = ( \text{Im} \, E _ { 1 } + \text{Im} \, E _ { 2 } ) C _ { 1 2 } \\ & \quad - ( \text{Re} \, E _ { 1 } - \text{Re} \, E _ { 2 } + c _ { 1 } n _ { 1 } - c _ { 2 } n _ { 2 } ) \tilde { j } _ { 1 2 }, \quad ( 3 4 d ) \quad \text{This}$$

where we defined C kl = ψ ψ k l + ψ ψ ∗ k ∗ l and ˜ j kl = i( ψ ψ k ∗ l -ψ ψ ∗ k l ) = j 12 /J 12 . In Eqs. (34a) and (34b) we see that the imaginary parts of the onsite energies act as additional sources or sinks to the particle flow, the quantities j e 1 = 2 n 1 Im E 1 and j e 2 = -2 n 2 Im E 2 can be identified as currents to and from the environment, thus, this Hamiltonian describes an open quantum system.

The general solutions for the nonlinear model (33) are calculated in [21]. For convinience we only discuss the results of the linear, i. e. non-interacting case ( c 1 = c 2 = 0 ). Then, the eigenvalues are given by

$$E _ { \pm } = \frac { E _ { 1 } + E _ { 2 } } { 2 } \pm \sqrt { J _ { 1 2 } ^ { 2 } + \frac { ( E _ { 1 } - E _ { 2 } ) ^ { 2 } } { 2 } }. \quad ( 3 5 ) \quad \text{To s} \quad.$$

Requirement for PT symmetry leads to the condition E 1 = E ∗ 2 , which is for instance fulfilled by E 1 = iΓ and E 2 = -iΓ , a situation with balanced gain and loss. Then, the eigenvalues are E ± = √ J 2 12 -Γ 2 . For Γ &lt; J 12 they are real, whereas for Γ &gt; J 12 the PT symmetry is broken and the eigenvalues are purely imaginary.

It is now our main purpose to investigate whether the behavior of the non-Hermitian two-mode model (33) can be described by a Hermitian four-mode model, which is given by

$$\mathbf H _ { 4 } ( t ) = \begin{pmatrix} E _ { 0 } ( t ) & - J _ { 0 1 } ( t ) & 0 & 0 \\ - J _ { 0 1 } ( t ) & E _ { 1 } & - J _ { 1 2 } & 0 \\ 0 & - J _ { 1 2 } & E _ { 2 } & - J _ { 2 3 } ( t ) \\ 0 & 0 & - J _ { 2 3 } ( t ) & E _ { 3 } ( t ) \end{pmatrix} \quad \text{where} \\ + \text{diag} \left ( c _ { 0 } \left | \psi _ { 0 } \right | ^ { 2 }, c _ { 1 } \left | \psi _ { 1 } \right | ^ { 2 }, c _ { 2 } \left | \psi _ { 2 } \right | ^ { 2 }, c _ { 3 } \left | \psi _ { 3 } \right | ^ { 2 } \right ). \quad ( 3 6 ) \text{ fullfill}$$

Two additional wells are coupled to the system, they act as a particle reservoir for the inner wells. The onsite energies and tunneling elements of the outer wells may be time-dependent. Since this model shall be Hermitian, the onsite energies are real, E k ∈ ❘ .

To derive a relationship between the two- and fourmode model we calculate the time derivatives of the same observables as in the two-mode model, which yields

$$\partial _ { t } n _ { 1 } = j _ { 0 1 } - j _ { 1 2 },$$

$$\text{the} \quad \partial _ { t } n _ { 2 } = j _ { 1 2 } - j _ { 2 3 },$$

$$\text{ly} _ { \mathbb { b } _ { \mathbb { c } } } \quad \partial _ { t } j _ { 1 2 } & = 2 J _ { 1 2 } ^ { 2 } \left ( n _ { 1 } - n _ { 2 } \right ) - J _ { 1 2 } \left ( J _ { 0 1 } C _ { 0 2 } - J _ { 2 3 } C _ { 1 3 } \right ) \\ \text{cle} _ { \mathbb { c } } & \quad + J _ { 1 2 } \left ( E _ { 1 } - E _ { 2 } + c _ { 1 } n _ { 1 } - c _ { 2 } n _ { 2 } \right ) C _ { 1 2 }, \quad \left ( 3 7 c \right )$$

$$( 3 7 c )$$

$$\text{ated}, \quad \partial _ { t } C _ { 1 2 } = J _ { 0 1 } \tilde { j } _ { 0 2 } - J _ { 2 3 } \tilde { j } _ { 1 3 } - ( E _ { 1 } - E _ { 2 } + c _ { 1 } n _ { 1 } - c _ { 2 } n _ { 2 } ) \tilde { j } _ { 1 2 }.$$

Comparing with Eqs. (34) we find that for the two-mode model the following condition must hold,

$$\text{Im} \, E _ { 1 } = - \text{Im} \, E _ { 2 } \Rightarrow E _ { 1 } = E _ { 2 } ^ { * }. \text{ \quad \ \ } ( 3 8 )$$

This means that only the PT -symmetric two-mode model, as discussed above, may be simulated by the four-mode model. From now on we write E 1 = iΓ and E 2 = -iΓ with Γ ∈ ❘ . Furthermore, the real parts of the onsite energies E 1 and E 2 and the coupling element J 12 have to agree between the two models. We are then left with the following conditions,

$$j _ { 0 1 } = 2 \Gamma n _ { 1 },$$

$$j _ { 2 3 } = 2 \Gamma n _ { 2 },$$

$$0 = J _ { 0 1 } C _ { 0 2 } - J _ { 2 3 } C _ { 1 3 }, \text{ \quad \ \ } ( 3 9 c )$$

$$0 = J _ { 0 1 } \tilde { j } _ { 0 2 } - J _ { 2 3 } \tilde { j } _ { 1 3 }. \text{ \quad \ \ } ( 3 9 d )$$

To summarize at this point, if these conditions are fulfilled at every time, the two middle wells of the Hermitian four-mode model behave exactly as the wells of the PT -symmetric two-mode model. These four conditions seem to be independent, but we shall see in App. A that Eq. (39d) follows if Eqs. (39a)-(39c) are fulfilled.

We now have to make clear that these conditions can indeed be fulfilled by giving suitable time-dependencies to the elements of the four-mode Hamiltonian. Equation (39c) can simply be fulfilled by choosing

$$J _ { 0 1 } ( t ) = d C _ { 1 3 } ( t ), \quad J _ { 2 3 } ( t ) = d C _ { 0 2 } ( t ), \quad ( 4 0 )$$

glyph[negationslash]

where d = 0 is a real, time-independent quantity. The tunneling elements are then time-dependent. With J 01 and J 23 determined there are no free parameters left to fulfill Eqs. (39a) and (39b). Instead, we can set the timederivatives of these equations with the help of the remaining onsite energies E 0 and E 3 . We calculate the timederivatives of j 01 and j 23 using the Schrödinger equation,

which yields

$$\partial _ { t } j _ { 0 1 } & = ( \partial _ { t } J _ { 0 1 } ) \tilde { j } _ { 0 1 } + 2 J _ { 0 1 } ^ { 2 } \left ( n _ { 0 } - n _ { 1 } \right ) + J _ { 0 1 } J _ { 1 2 } C _ { 0 2 } \\ & \quad + J _ { 0 1 } \left ( E _ { 0 } - E _ { 1 } + c _ { 0 } n _ { 0 } - c _ { 1 } n _ { 1 } \right ) C _ { 0 1 } \\ & \stackrel {! } { = } \partial _ { t } j _ { 0 1 } ^ { \text{tar} }, \\ \wedge \cdot \quad \tau \cdot \tilde { \lambda } \cdot \tilde { \lambda } \cdot \tilde { \lambda } \cdot \tilde { \lambda } \cdot \tilde { \lambda } \cdot \tilde { \lambda } \cdot \tilde { \lambda } \cdot \tilde { \lambda } \cdot \tilde { \lambda }$$

$$\begin{array} { c } \bar { \ } _ { \L v _ { 0 1 } } \,, & \vee \cdots \,, \\ \partial _ { t } j _ { 2 3 } = & ( \partial _ { t } J _ { 2 3 } ) \, \tilde { j } _ { 2 3 } + 2 J _ { 2 3 } ^ { 2 } \, ( n _ { 2 } - n _ { 3 } ) - J _ { 2 3 } J _ { 1 2 } C _ { 1 3 } \\ & + J _ { 2 3 } \, ( E _ { 2 } - E _ { 3 } + c _ { 2 } n _ { 2 } - c _ { 3 } n _ { 3 } ) \, C _ { 2 3 } \\ \overset {! } { = } & \partial _ { t } j _ { 2 3 } ^ { \text{tar} }, & ( 4 1 b ) \end{array}$$

$$& \stackrel {. } { = } & \partial _ { t } j _ { 0 1 } ^ { \quad. } \quad ( 4 1 \tt a ) \quad ( 4 1 \tt a ) \quad \\ _ { 3 } & = ( \partial _ { t } J _ { 2 3 } ) \, \tilde { j } _ { 2 3 } + 2 J _ { 2 3 } ^ { 2 } \, ( n _ { 2 } - n _ { 3 } ) - J _ { 2 3 } J _ { 1 2 } C _ { 1 3 } \quad & 0. \\ & \quad + J _ { 2 3 } \, ( E _ { 2 } - E _ { 3 } + c _ { 2 } n _ { 2 } - c _ { 3 } n _ { 3 } ) \, C _ { 2 3 } \quad & 0. \\ & \stackrel {! } { = } & \partial _ { t } j _ { 2 3 } ^ { \quad \tt a r } \quad ( 4 1 b ) \quad & 0.$$

with the target currents determined by Eqs. (39a) and (39b). The time derivatives of the tunneling elements J 01 and J 23 are given by

$$\partial _ { t } J _ { 0 1 } & = d \partial _ { t } C _ { 1 3 } \\ & = d \left ( J _ { 0 1 } \tilde { j } _ { 0 3 } + J _ { 1 2 } \tilde { j } _ { 2 3 } - J _ { 2 3 } \tilde { j } _ { 1 2 } \right ) \\ & \quad + d \left ( E _ { 3 } - E _ { 1 } + c _ { 3 } n _ { 3 } - c _ { 1 } n _ { 1 } \right ) \tilde { j } _ { 1 3 }, \quad \left ( 4$$

$$\begin{array} { c } \kappa \, \langle \nu _ { 1 1 } \nu _ { 0 0 } \cdot \, \kappa _ { 1 2 } \nu _ { 2 0 } \cdot \, \kappa _ { 2 0 1 1 } \rangle \\ & + d \left ( E _ { 3 } - E _ { 1 } + c _ { 3 } n _ { 3 } - c _ { 1 } n _ { 1 } \right ) \tilde { j } _ { 1 3 }, \quad & ( 4 2 \tt a ) \quad & \text{Figur} \\ & + d \left ( E _ { 3 } - E _ { 0 } + c _ { 2 } \\ & = d \left ( J _ { 0 1 } \tilde { j } _ { 1 2 } - J _ { 1 2 } \tilde { j } _ { 0 1 } - J _ { 2 3 } \tilde { j } _ { 0 3 } \right ) \\ & + d \left ( E _ { 2 } - E _ { 0 } + c _ { 2 } n _ { 2 } - c _ { 0 } n _ { 0 } \right ) \tilde { j } _ { 0 2 }. \quad & ( 4 2 b ) \quad & \text{of } J _ { 1 1 }. \\ & \text{and } o \end{array}$$

Inserting this into Eqs. (41), the onsite energies E 0 and E 3 are determined by a linear set of equations,

$$\begin{pmatrix} C _ { 0 1 } C _ { 1 3 } & \tilde { j } _ { 0 1 } \tilde { j } _ { 1 3 } \\ - \tilde { j } _ { 0 2 } \tilde { j } _ { 2 3 } & - C _ { 0 2 } C _ { 2 3 } \end{pmatrix} \begin{pmatrix} E _ { 0 } \\ E _ { 3 } \end{pmatrix} = \begin{pmatrix} v _ { 0 } \\ v _ { 3 } \end{pmatrix}, \quad ( 4 3 ) \quad \text{cho} \\ \quad \text{ we c}$$

with the entries

$$v _ { 0 } & = \partial _ { t } j _ { 0 1 } ^ { \text{tar} } / d - J _ { 1 2 } \tilde { j } _ { 2 3 } - d \left ( C _ { 1 3 } \tilde { j } _ { 0 3 } - C _ { 0 2 } \tilde { j } _ { 1 2 } \right ) \tilde { j } _ { 0 1 } \\ & \quad - 2 d C _ { 1 3 } ^ { 2 } \left ( n _ { 0 } - n _ { 1 } \right ) - J _ { 1 2 } C _ { 0 2 } C _ { 1 3 } \\ & \quad - \left ( c _ { 3 } n _ { 3 } - c _ { 1 } n _ { 1 } \right ) \tilde { j } _ { 0 1 } \tilde { j } _ { 1 3 } - \left ( c _ { 0 } n _ { 0 } - c _ { 1 } n _ { 1 } \right ) C _ { 0 1 } C _ { 1 3 }, \\ & \quad \sim \text{tar} \, \cdots \text{ } \tilde { \cdot } \text{ } \tilde { \cdot } \text{ } \tilde { \cdot } \text{ } \tilde { \cdot } \text{ } \tilde { \cdot } \text{ } \tilde { \cdot } \text{ } \tilde { \cdot } \text{ } \tilde { \cdot } \text{ } \tilde { \cdot } \text{ } \tilde { \cdot } \text{ } \tilde { \cdot } \text{ } \tilde { \cdot } \text{ } \tilde { \cdot } \text{ } \tilde { \cdot }$$

$$v _ { 3 } & = \partial _ { t } j _ { 2 3 } ^ { \text{tar} } / d + J _ { 1 2 } \tilde { j } _ { 0 1 } - d \left ( C _ { 1 3 } \tilde { j } _ { 1 2 } - C _ { 0 2 } \tilde { j } _ { 0 3 } \right ) \tilde { j } _ { 2 3 } \\ & \quad - 2 d C _ { 0 2 } ^ { 2 } \left ( n _ { 2 } - n _ { 3 } \right ) + J _ { 1 2 } C _ { 0 2 } C _ { 1 3 } & \quad \text{Now,} \\ & \quad - \left ( c _ { 2 } n _ { 2 } - c _ { 0 } n _ { 0 } \right ) \tilde { j } _ { 0 2 } \tilde { j } _ { 2 3 } - \left ( c _ { 2 } n _ { 2 } - c _ { 3 } n _ { 3 } \right ) C _ { 0 2 } C _ { 2 3 }, & \quad \text{ments} \\ & \quad - \text{(44b)} & \quad \text{two-n}$$

where we set E 1 = E 2 = 0 , since - as discussed above - the real parts of these elements have to agree between the two- and four-mode model.

We have to note that by means of the onsite energies E 0 and E 3 we do not fix the currents, but their timederivatives. For this reason the initial probabilities and currents have to be chosen such that they fulfill the conditions (39). In the following section we apply these results to obtain solutions for different starting conditions.

## B. Results

Before presenting our results we give the necessary initial conditions. We insert the definition of the particle current and the solutions for the tunneling elements (40)

0.6

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0158v1-figure-77c1af6e-74f3-4c67-8736-adb7ac6ffa67.jpg)

40

Figure 2. Realization of a stationary PT -symmetric solution within the four-mode model. (a) Particles numbers in well 1 (solid) and well 2 (dashed), both equal to 1 2 / . (b) Particle numbers in well 0 (solid) and well 3 (dashed). (c) Particle currents j 01 (solid), j 12 (dashed) and j 23 (dash-dotted) in units of J 12 / glyph[planckover2pi1] . (d) Tunneling elements J 01 (solid), J 23 (dashed), and onsite energies E 0 (solid), E 3 (dashed), in units of J 12 .

into the two Eqs. (39a) and (39b) and express everything in terms of the wave function ψ k . We are allowed to choose the global phase of the initial wave function, so we choose ψ I 2 = 0 . We assume the real parts ψ R 0 and ψ R 3 to be arbitrary but fixed. We are left with (after simple rearrangement)

$$\psi _ { 3 } ^ { I } & = \frac { \Gamma } { 2 d \psi _ { 0 } ^ { R } }, \\ & \quad \,.$$

$$\psi _ { 0 } ^ { I } = \frac { \psi _ { 0 } ^ { R } \psi _ { 1 } ^ { I } } { \psi _ { 1 } ^ { R } } - \frac { \Gamma n _ { 1 } } { 2 d \psi _ { 1 } ^ { R } \left ( \psi _ { 1 } ^ { R } \psi _ { 3 } ^ { R } + \psi _ { 1 } ^ { I } \psi _ { 3 } ^ { I } \right ) }. \quad ( 4 5 b )$$

Now, the wave function ψ k fulfills the necessary conditions and with the time-dependencies of the matrix elements from Sec. III A we can simulate the PT -symmetric two-mode model.

First we consider a stationary solution in the noninteracting case [cf. Eq. (35)] for Γ /J 12 = 0 5 . . Fig. 2 shows the results. As expected the particle numbers in the two middle wells are both equal and constant in time (with normalization n 1 + n 2 = 1 ), as well as the particle currents. For this reason, the number of particles decreases (left well) or increases linearly (right well), with the slope ˙ n 0 = -Γ / glyph[planckover2pi1] and ˙ n 3 = Γ / glyph[planckover2pi1] , respectively. The matrix elements of the four-mode model vary slightly in time. Thus, we have shown that the conditions (39) can be fulfilled for a finite time.

As a second example we prepare a non-stationary solution at t = 0 , with ψ 1 (0) = √ 0 6 . and ψ 2 (0) = √ 0 4 . (see Fig. 3). In this case the particle numbers in the two middle wells oscillate in time with a phase difference ∆ φ &lt; π leading to a non-constant added number of particles, which is a typical feature of PT -symmetric

35

Figure 3. Same as Fig. 2, but for oscillatory dynamics. Additionally, the added number of particles n 1 + n 2 is plotted as crosses in (a). As the particle number in the reservoir well 0 decreases and gets close to zero, the conditions cannot be fulfilled anymore and the simulation breaks down.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0158v1-figure-48939cde-e691-4ecf-8b41-0f6a272880c6.jpg)

12

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0158v1-figure-17a5c7bd-716c-4b16-be84-510b344e4c50.jpg)

20

Figure 4. Same as Fig. 2, but for the interacting case ( c = -1 ), which leads to a collapsing dynamics, indicating a change of stability. The reservoir is emptying in a short time domain such that the time available for PT symmetry is quite short.

systems. Since well 0 acts as a particle reservoir the number of particles decreases and gets close to zero up to the point where the linear system of equations (43) cannot be fulfilled anymore (determinant of coefficient matrix equals zero). This shows that the time available for a simulation of PT symmetry with a reservoir is limited. The matrix elements in this case show a quasi-oscillatory behavior.

Next, we consider a system with attractive interaction c = -1 . In such a system it is known that the ground state changes its stability in an additional bifurcation (for repulsive interaction, the excited state changes its stability) [12]. This phenomenon is highly correlated to the

Figure 5. Same as Fig. 2, but for an adiabatic current ramp ( Γ f /J 12 = 0 5 . , t /t f 0 = 60 ) for realistic parameters of the external potential. Instead of the matrix elements, in (d) we show the parameters of the potential, δ k marks the displacement from equidistant potential wells. There is a total number of particles of N = 10 4 . See text for units.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0158v1-figure-726743f5-bf86-4534-9fa4-f8e95606b3a1.jpg)

occurrence of self-trapping states. In our example 2 × 2 system [Eq. (33)] the ground state is already unstable for the given interaction strength, so that the ground state with a small perturbation is expected to collapse. The calculation is shown in Fig. 4. The particle number n 1 increases exponentially, which marks the beginning of collapse. To support the exponential increase, the particle number in the reservoir n 0 is decreasing exponentially, so that the simulation breaks down after a quite short time interval. The matrix elements now show a rather complex time-dependency.

As the last example we calculate an experimentally more realistic scenario. Up to now - as mentioned above - we had to start with appropriate wave functions in order to obtain PT symmetry. This is nearly impossible in an experiment. To overcome this, we use an adiabatic change of the PT parameter from zero to its target value. We start at the ground state of the Hermitian four-mode model (with Γ = 0 ) and increase Γ to its target value Γ f . If this change is slow enough, we arrive approximately at the PT -symmetric ground state. The specific time dependence we use is Γ( ) = Γ [1 t f -cos( πt/t f )] / 2 for t ∈ [0 , t f ] .

In the following we simulate a condensate of N = 10 5 atoms of 87 Rb. We use units based on the potential width in z -direction, w z , which we set w z = 1 µ m . The basic unity of energy is E 0 = glyph[planckover2pi1] 2 /mw 2 z , which yields E /h 0 = 116 Hz . The basic unit of time is t 0 = mw / 2 z glyph[planckover2pi1] = 1 37 ms . . The wells have widths of w x = w y = 4 w z = 4 µ m , and are initially positioned equidistantly with a distance of 1 8 . w z = 1 8 . µ m . The potential depths are given by V 0 = V 3 = -60 E 0 and V 1 = V 2 = -45 E 0 . We use a scattering length of a = 2 83 . a B with a B being the Bohr radius.

We calculate the adiabatic ramp for Γ f /J 12 = 0 5 . and t /t f 0 = 60 . The results are shown in Fig. 5. The total particle number and currents are scaled to a particle number of N = 10 5 , the currents are given in units of 1 /t 0 , the potential depths in units of E 0 . Instead of the positions of the potential wells we plot the difference of the position from its initial value, δ k = s k z -s k z ( t = 0) . By means of the relations derived in Sec. II we can transform from the calculated matrix elements to parameters of the external potential.

After time t f an approximate PT -symmetric state is obtained. This can be seen by the constant number of particles in the middle wells, despite an increasing and decreasing number in wells 3 and 0, respectively, which is a clear sign of PT symmetry, also in a possible experiment. The potential depths of the outer wells have to be adjusted in a simple way, the positions vary within a few percent of the distance of the wells. This shows that the simple four-mode model can be used to describe a realistic scenario and to obtain the necessary parameters of the external potential.

We have applied our results from Sec. III A to different initial conditions and system parameters and have shown that it is indeed possible to realize a PT -symmetric system with a Hermitian four-mode model. In the following Section we investigate this scenario in the framework of the GPE and a variational ansatz with Gaussian functions. In particular we verify that the results from this section agree.

## IV. VARIATIONAL ANSATZ WITH GAUSSIAN FUNCTIONS

## A. Variational ansatz and equations of motion

After the simple approach in the previous section we now investigate the possible realization in the framework of the GPE with an extended variational ansatz. We use a superposition of Gaussian functions, where - in contrast to the simpler approach of Sec. II A - all variational parameters are considered as time-dependent,

$$\psi = \sum _ { k } \, \mathrm e ^ { - ( r - q ^ { k } ) ^ { \mathrm T } \mathrm A ^ { k } ( r - q ^ { k } ) + \mathrm i p ^ { k } ( r - q ^ { k } ) - \gamma ^ { k } }. \quad ( 4 6 ) \quad \text{exter} \quad \text{space} \quad \text{$space} \quad \text{$total_s$}$$

The matrices A k = diag( A ,A ,A k x k y k z ) ∈ ❈ 3 × 3 describe the width of each Gaussian function in each direction, the vectors q k = (0 0 , , q k z ) T ∈ ❘ 3 and p k = (0 0 , , p k z ) T ∈ ❘ 3 the position and momentum, respectively, and the scalar quantities γ k ∈ ❈ the amplitude and phase. This variational ansatz can describe a much richer dynamics than ansatz (4) and thus we expect more accurate results.

The equations of motion follow from the timedependent variational principle in the formulation of McLachlan [22], which states that the quantity

$$I = | | \text{i} \hbar { \phi } - \hat { H } \psi | | ^ { 2 }$$

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0158v1-figure-d007c20a-787f-4c02-b350-3c7ed4f0b164.jpg)

z

Figure 6. External potential V ext ; the darker, the deeper the potential. The wells are separated into regions by walls at positions right in the middle between two minima. The integral over | ψ | 2 in a region k is the particle number n k ; the current through a wall is denoted by j k,k +1 . See text for exact definitions. With this separation we obtain a formulation analogous to the four-mode model (Sec. III).

shall be minimized with respect to φ and φ = ˙ ψ is set afterwards. Details of the necessary calculation can be found in Ref. [23], where the ansatz (46) has been used to study the collision of anisotropic solitons in a BEC. The EOM are given by

$$\mathrm i \hbar { \sum } _ { k } \left \langle \frac { \partial \psi } { \partial z _ { l } } \left | \frac { \partial \psi } { \partial z _ { k } } \right \rangle \dot { z } _ { k } = \left \langle \frac { \partial \psi } { \partial z _ { l } } \right | \hat { H } \left | \psi \right \rangle, \quad ( 4 8 )$$

where z = ( A ,. . . , γ 1 x N G ) T is the vector of all variational parameters. Stationary solutions are given by fixed points of (48), whereas the dynamics can be calculated by integrating Eq. (48) with a numerical integrator.

In Sec. III we started our considerations by calculating the time-derivatives of the observables. In the case of the continuous GPE with a non-Hermitian potential this leads to the modified equation of continuity

$$\dot { \rho } + \text{div} \, j = 2 \rho \text{ Im } V,$$

with ρ = | ψ | 2 and j = glyph[planckover2pi1] ( ψ ∗ ∇ ψ -ψ ∇ ψ ∗ ) / 2 m i . The imaginary part of the potential acts as an additional source or sink to the probability flow. Instead of adjusting the matrix elements, we need to adjust the parameters of the external potential. Since Eq. (49) is given on the whole space ❘ 3 , one would need to change the external potential at every point in space, which is not manageable, neither in theory nor in an experiment.

To overcome this problem we use the results of Sec. III as a roadmap for this continuous system. For that, we separate the individual wells by walls between two minima of the potential at positions z k = ( s k z + s k +1 z ) / 2 for k = 1 2 3 , , , z 0 → -∞ and z 4 → ∞ (see Fig. 6). To obtain time derivatives of discrete observables as in the few-mode model, we integrate the continuity equation for the Hermitian four-well potential over a volume

$$V ^ { k } = \{ ( x, y, z ) ^ { \mathrm T } | - \infty < x, y < \infty, z _ { k - 1 } < z < z _ { k } \}.$$

This yields

$$\mu _ { \nu } \dots \quad & \partial _ { t } \iiint _ { V ^ { k } } d ^ { 3 } r \, \rho & = - \iiint _ { V ^ { k } } d ^ { 3 } r \, \text{div} \, j & \quad \begin{array} { c c c c } 5 0 \\ 4 8 \\ 4 4 \\ 4 2 \\ 4 0 \\ 4 0 \\ \cdot \, & - \iiint _ { \partial V ^ { k } } d A \cdot j. & \quad ( 5 1 ) & \quad \begin{array} { c c c c } 5 0 \\ 4 0 \\ 5 5 0 \\ 5 0 0 \\ 5 0 0 \\ \cdot \, & \cdot \, & \cdot \, & \cdot \, & \cdot \, & \cdot \, & \cdot \, & \cdot \, & \cdot \, & \cdot \, & \cdot \, & \cdot \, & \cdot \, & \cdot \, & \cdot \, & \cdot \, & \cdot \, & \cdot \, & \cdot \, & \cdot \, & \cdot \, & \cdot \, & \cdot \, & \cdot \, & \cdot \, & \cdot \, & \cdot \, & \cdot \, & \cdot \, & \cdot \, & \cdot \, & \cdot \, & \cdot \, & \cdot \, & \end{array}$$

The integral on the left-hand side can be interpreted as the number of particles in well k . On the right-hand side we used Gauss's theorem. The surface integrals goes over a box, where the areas in the x z -- and y z --plane go to infinity. We assume the current j to also vanish at infinity, so only the integral over the two surfaces in the x y --plane contribute. Thus, we can write

$$\partial _ { t } n _ { k } = - \underset { - \infty } { \overset { \infty } { \int } } \, \mathrm d x \underset { - \infty } { \overset { \infty } { \int } } \, \mathrm d y \left [ j _ { z } ( x, y, z _ { k } ) - j _ { z } ( x, y, z _ { k - 1 } ) \right ]. \, \underset { - \infty } { \overset { \infty } { \int } } \, \mathrm d x \underset { 3 } { \overset { \infty } { \int } } \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \cdots } { \overset { \infty } { \int } } \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \,  \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm d x \, \mathrm \cdots }$$

The integral over the z -component of j evaluated at z k represents the current from well k to k +1 , which we write as j k,k +1 , for z k -1 this represents the current from well k -1 to k , j k -1 ,k . We finally obtain from the continuity equation

$$\partial _ { t } n _ { k } = j _ { k - 1, k } - j _ { k, k + 1 }, \quad \ \ ( 5 3 ) \quad ^ { \text{sy$} } _ { \infty \, \vee }$$

where j -1 0 , = j 3 4 , = 0 .

Therefore, with Eq. (53) we have an analogous equation as for the few-mode-model (37). We can use the same method to realize a PT -symmetric system. At every time step of the numerical integration of the equations of motion (48) we vary the parameters of the external potential such that the observables show the desired behavior. This is done via a nonlinear root search, in contrast to the analytical solutions of Sec. III A. In Sec. III B we showed that the positions of the outer wells barely needed to be adjusted. As a simplification we only vary the depths of the outer wells V 0 , V 3 and demand the outer currents j 01 , j 23 to reach the desired values. We are left with a two-dimensional root search. In the following section we determine the realization of PT symmetry within the GPE and compare the results with those of the few-mode-model.

## B. Results and comparison with few-mode models

We use the variational ansatz to simulate the adiabatic current ramp of Sec. III B (with same parameters) for a BEC in a four-well potential. For simplification, we do not give a specific value of the PT parameter Γ but a target particle current of j tar = 5 /t 0 . Fig. 7 shows the results. It is not a priori clear whether the quantity n k = ( d k eff ) ∗ d k eff - resulting from the transformed amplitudes of the simplified variational approach - can be compared with the integrated probability density n k of

Figure 7. Adiabatic current ramp calculated with variational ansatz (lines) and four-mode model (crosses) for the same system. Shown are (a) particle numbers in the middle wells, (b) particle numbers in outer wells, (c) potential depths of the outer wells in units of E 0 , and (d) particle current j 12 . There are small deviations, but overall the four-mode model is a valid approximation compared to the variational ansatz.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0158v1-figure-b4fb14e1-55c1-41c1-9c88-825e5cf035cc.jpg)

the extended variational approach. We show in App. B that this comparison is indeed possible and reasonable.

Overall, there is the same qualitative behavior between the two approaches. The particle number in the middle wells is slightly under estimated in the four-mode model due to its origin as an approximation. Because of the fixed positions of the outer wells for the variational approach, n 1 and n 2 are not exactly equal, but the difference is small compared to the absolute number. We can conclude that we can now also create an (approximate) PT -symmetric state within the variational approach and verify that the results from the four-mode model are a good approximation.

The variational approach gives us access to the wave function at each time step (see Fig. 8). As indicated by the particle numbers in Fig. 7, the amplitudes in the middle wells stay almost constant, whereas the amplitudes in the outer wells increase or decrease in time. This could also be measured in an experiment and could serve as a proof of realized PT symmetry.

0.20

Figure 8. Modulus squared of the wave function | ψ | 2 at x, y = 0 for different times t . The amplitudes of the outer parts increase or decrease, whereas the amplitudes in the middle wells stay constant, a clear sign of PT symmetry, also in an experiment.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0158v1-figure-2d315197-e3dc-412c-93f5-15a38480de7e.jpg)

## V. CONCLUSION

By means of a simple variational ansatz -with localized time-independent Gaussian functions - and the method of symmetric orthogonalization we could transform the GPE for a multi-well potential to a simple fewmode model. The parameters of the variational ansatz are obtained by an energy minimization process in a computationally cheap way. With the results one can, on the one hand, calculate the elements of the few-mode model for a specific system, and on the other determine the parameters of the external potential by knowing the timedependent matrix elements.

These results are applied to the PT -symmetric twomode model and to the Hermitian four-mode model. With this model we could show that the four-mode model can be designed in such a way such that the middle wells behave exactly as the two-wells of the PT -symmetric two-mode model. Thus, this Hermitian (closed) system can be used as a possible realization of a quantum mechanical PT -symmetric system. With an extended Gaussian variational ansatz - with full time-dependent Gaussian functions, which can describe a much richer dynamics - we confirmed the qualitative results of the few-mode model.

In this paper we applied the idea to embed a two-well system into a closed system for the simplest case. It is possible for future work to use a greater embedding which could result in a simpler time-dependence of the additional parameters. The ideas presented in this paper could also serve as a possible road map for an investigation of the Bose-Hubbard model and its many-particle effects.

## ACKNOWLEDGMENTS

This work was supported by DFG. M.K. is grateful for support from the Landesgraduiertenförderung of the Land Baden-Württemberg.

## Appendix A: Analytical considerations

In Sec. III A we derived four conditions [Eqs. (39)] such that the Hermitian four-mode model can simulate the PT -symmetric two-mode model. However, only three of them are independent, which will be shown in this appendix. For that, we need the relations

$$C _ { j j } \tilde { j } _ { i k } = C _ { i j } \tilde { j } _ { j k } + \tilde { j } _ { i j } C _ { j k }, \text{ \quad \ \ } ( \text{A1a} )$$

$$C _ { j j } C _ { i k } = C _ { i j } C _ { j k } - \tilde { j } _ { i j } \tilde { j } _ { j k }, \text{ \quad \ \ } ( \text{A1b} )$$

which can be proved by simply inserting the definitions of these quantities.

We start with Eqs. (39a) and (39b), insert the solutions (40), take the square and use Eqs. (A1). We get

$$\tilde { j } _ { 0 1 } ^ { 2 } \tilde { j } _ { 2 3 } ^ { 2 } + \frac { 2 \Gamma } { d } \tilde { j } _ { 1 2 } \tilde { j } _ { 0 1 } \tilde { j } _ { 2 3 } - C _ { 1 2 } ^ { 2 } \frac { n _ { 3 } } { n _ { 1 } } \tilde { j } _ { 0 1 } ^ { 2 } + \frac { 4 \Gamma ^ { 2 } } { d ^ { 2 } } n _ { 1 } n _ { 2 } = 0, \\ \intertext { e } _ { \omega }$$

$$\ e s ^ { \quad \tilde { j } _ { 0 1 } ^ { 2 } \tilde { j } _ { 2 3 } ^ { 2 } } + \frac { 2 \Gamma } { d } \tilde { j } _ { 1 2 } \tilde { j } _ { 0 1 } \tilde { j } _ { 2 3 } - C _ { 1 2 } ^ { 2 } \frac { n _ { 0 } } { n _ { 2 } } \tilde { j } _ { 2 3 } ^ { 2 } + \frac { 4 \Gamma ^ { 2 } } { d ^ { 2 } } n _ { 1 } n _ { 2 } = 0. \\ \ t z ^ { ( \text{A2b} ) }$$

Taking the difference yields

$$\tilde { j } _ { 2 3 } ^ { 2 } & = \frac { n _ { 2 } n _ { 3 } } { n _ { 0 } n _ { 1 } } \tilde { j } _ { 0 1 } ^ { 2 } \\ \Leftrightarrow \tilde { j } _ { 2 3 } & = s _ { 3 } \sqrt { \frac { n _ { 2 } n _ { 3 } } { n _ { 0 } n _ { 1 } } } \tilde { j } _ { 0 1 }, \quad \quad ( \Lambda 3 )$$

where s 3 = ± 1 gives the sign of ˜ j 23 compared to ˜ j 01 . Inserting this result into Eq. (A2a) gives a fourth-order polynomial equation for ˜ j 01 . The result is

$$\underset { \mathfrak { m } _ { 1 } } { \mathfrak { e } } _ { \mathfrak { m } _ { 2 } } \quad \tilde { j } _ { 0 1 } = s _ { 2 } \sqrt { 2 n _ { 0 } n _ { 1 } \left ( 1 - \alpha + s _ { 1 } \sqrt { ( 1 - \alpha ) ^ { 2 } - \beta ^ { 2 } } \right ) }, \ \ ( A 4 )$$

with the definitions

$$\alpha = \frac { \gamma } { 2 } \underset { \varepsilon \quad. } { \left ( \beta + \frac { \gamma } { 2 } \right ) } \,, \quad \quad \quad ( A 5 \tt a )$$

$$\beta = s _ { 3 } \frac { \Gamma } { d } \frac { 1 } { \sqrt { n _ { 0 } n _ { 3 } } }, \quad \quad ( A 5 b )$$

$$\gamma = \frac { \tilde { j } _ { 1 2 } ^ { \cdot } } { \sqrt { n _ { 1 } n _ { 2 } } }. \quad \quad \quad ( A 5 c )$$

Using Eqs. (A1) we can find all other quantities ˜ j kl and C kl , especially those needed for the investigation of the

fourth condition J 01 ˜ j 02 -J 23 ˜ j 13 = dC 13 ˜ j 02 -dC 02 ˜ j 23 ,

$$C _ { 0 2 } = s _ { 2 } \sin g \, d \sqrt { 2 n _ { 0 } n _ { 2 } \left ( 1 - \alpha - s _ { 1 } \sqrt { ( 1 - \alpha ) ^ { 2 } - \beta ^ { 2 } } \right ) }, \\ ( A 6 a )$$

$$C _ { 1 3 } = s _ { 2 } \, \text{sign} \, d \sqrt { 2 n _ { 1 } n _ { 3 } \left ( 1 - \alpha - s _ { 1 } \sqrt { ( 1 - \alpha ) ^ { 2 } - \beta ^ { 2 } } \right ) }, \quad \text{plus t} \\ \frac { ( A 6 b ) } { \sqrt { \ } \text{Th} }$$

$$\tilde { j } _ { 0 2 } = s _ { 6 } \sqrt { 2 n _ { 0 } n _ { 2 } \left ( 1 + \alpha + s _ { 1 } \sqrt { ( 1 - \alpha ) ^ { 2 } - \beta ^ { 2 } } \right ) }, \ ( A 6 c ) \quad \psi$$

$$\tilde { j } _ { 1 3 } = s _ { 6 } \sqrt { 2 n _ { 1 } n _ { 3 } \left ( 1 + \alpha + s _ { 1 } \sqrt { ( 1 - \alpha ) ^ { 2 } - \beta ^ { 2 } } \right ) }. \ ( A 6 d ) \quad \text{with $1$}$$

All signs s i are determined by the phases of the initial wave function. We can now calculate the fourth condition and find

$$J _ { 0 1 } \tilde { j } _ { 0 2 } - J _ { 2 3 } \tilde { j } _ { 1 3 } = 0. \quad \quad ( A 7 ) \quad \text{wii}$$

Thus, the first three conditions of Eqs. (39) imply the validity of the fourth one. With the results of the matrix elements of Sec. III A, namely J 01 , J 23 , E 0 and E 3 , we find an exact equivalence of the two- and four-mode models also for interacting atoms. Furthermore, with the results of this appendix, we can calculate these matrix elements, once the quantities n 1 , n 2 and ˜ j 12 are known.

## Appendix B: Comparison of probabilities for few-mode-model and Gaussian functions

As discussed in Sec. IVB it is not a priori clear whether we can compare the particle numbers obtained from the four-mode model and the extended variational approach. In this appendix we take a closer look at both quantities.

The transformation of the amplitudes from the fourmode model, resulting from symmetric orthogonalization, is given by the inverse of the matrix X (cf. Sec. II C),

$$d _ { \text{eff} } ^ { l } & = \sum _ { k } ( \mathbf X ^ { - 1 } ) _ { l k } d ^ { k } \\ & = \sum _ { k } \left [ ( \mathbf X ^ { - 1 } ) _ { l k } ^ { ( 0 ) } + ( \mathbf X ^ { - 1 } ) _ { l k } ^ { ( 1 ) } \right ] d ^ { k }. \quad ( B 1 ) \begin{pmatrix} \text{$\text{can}$} \\ \text{dist} \\ \text{mod} \\ \end{pmatrix}.$$

The individual orders of the inverse matrix can easily be calculated, we get

$$\text{calculated, we get} \\ ( X ^ { - 1 } ) _ { l k } ^ { ( 0 ) } & = \sqrt { \frac { \pi ^ { 3 } } { 8 } } \frac { 1 } { \sqrt { A _ { x, R } ^ { k } A _ { y, R } ^ { k } A _ { z, R } ^ { k } } } \delta _ { l k }, \\ ( X ^ { - 1 } ) _ { l k } ^ { ( 1 ) } & = \sqrt { 8 \pi ^ { 3 } } \frac { \sqrt { A _ { x, R } ^ { k } A _ { y, R } ^ { k } A _ { z, R } ^ { k } A _ { x, R } ^ { l } A _ { y, R } ^ { l } A _ { z, R } ^ { l } } } { \sqrt { A _ { x, R } ^ { k } A _ { y, R } ^ { k } A _ { z, R } ^ { k } } } + \sqrt { A _ { x, R } ^ { l } A _ { y, R } ^ { l } A _ { z, R } ^ { l } } } + \sqrt { \text{unify.} } \text{bersc} \\ & \quad \times \frac { c ^ { k l } } { \sqrt { A _ { x } ^ { k l } } \sqrt { A _ { y } ^ { k l } } \sqrt { A _ { z } ^ { k l } } } ( \delta _ { k, l - 1 } + \delta _ { k, l + 1 } ). \ ( B 2 a ) \ \text{Since} \ \text{the} \ p \varepsilon$$

The main contribution to the particle number n k = ( d k eff ) ∗ d k eff is then given by

$$n _ { k } = \sqrt { \frac { \pi ^ { 3 } } { 8 } } \frac { 1 } { \sqrt { A _ { x, R } ^ { k } A _ { y, R } ^ { k } A _ { z, R } ^ { k } } } \left | d ^ { k } \right | ^ { 2 }, \quad \text{(B3)}$$

plus terms involving neighboring amplitudes.

The extended variational ansatz is given by

$$\psi = \sum _ { k } d ^ { k } e ^ { - A _ { x } ^ { k } x ^ { 2 } - A _ { y } ^ { k } y ^ { 2 } - A _ { z } ^ { k } ( z - q _ { z } ^ { k } ) ^ { 2 } + i ( z - q _ { z } ^ { k } ) p _ { z } ^ { k } }, \quad ( B 4 )$$

with the amplitude d k = exp( -γ k ) . The probability density can then be written as

$$\underset { \rho } { \pi } = \sum _ { k, l } d ^ { k } ( d ^ { l } ) ^ { * } \, e ^ { - A _ { x } ^ { k l } x ^ { 2 } - A _ { y } ^ { k l } y ^ { 2 } - A _ { z } ^ { k l } z ^ { 2 } + i \tilde { p } _ { z } ^ { k l } z - \tilde { \gamma } ^ { k l } }, \quad ( B 5 )$$

with the definitions

$$\text{he} \quad \tilde { p } _ { z } ^ { k l } = 2 [ A _ { z } ^ { k } q _ { z } ^ { k } + ( A _ { z } ^ { l } ) ^ { * } q _ { z } ^ { l } ] + \text{i} ( p _ { z } ^ { k } - p _ { z } ^ { l } ), \quad \text{(B6a)}$$

$$\frac { \iota \sigma } { 2 }, \quad \tilde { \gamma } ^ { k l } = A _ { z } ^ { k } ( q _ { z } ^ { k } ) ^ { 2 } + ( A _ { z } ^ { l } ) ^ { * } ( q _ { z } ^ { l } ) ^ { 2 } + i ( q _ { z } ^ { k } p _ { z } ^ { k } - q _ { z } ^ { l } p _ { z } ^ { l } ). \ \ ( B 6 b )$$

We integrate ρ over the region discussed in Sec. IVA. The integral over x and y yields

$$\rho _ { z } = \sum _ { k, l } d ^ { k } ( d ^ { l } ) ^ { * } \frac { \pi } { \sqrt { A _ { x } ^ { k l } } \sqrt { A _ { y } ^ { k l } } } \, e ^ { - A _ { z } ^ { k l } z ^ { 2 } + \tilde { p } _ { z } ^ { k l } z - \tilde { \gamma } ^ { k l } }. \quad ( B 7 )$$

The integral over z goes over a finite interval, say a ≤ z ≤ b , we can express the result in terms of the error function,

$$\text{clear} \quad \text{uncount}, \\ \text{ional} \quad \text{$n(a,b)=\frac{1}{2}\sum_{k,l}d^{k}(d^{l})^{*}\frac{\pi^{3}}{\sqrt{A_{x}^{k}}\sqrt{A_{y}^{k}}\sqrt{A_{z}^{k}}^{l}-\tilde{\gamma}^{k}l}$} \\ \text{four} \quad \text{$\times\left[\text{erf}\left(\frac{2A_{z}^{k}b-\tilde{p}_{z}^{k}l}{4\sqrt{A_{z}^{k}} \right) -\text{erf} \left(\frac{2A_{z}^{k}a-\tilde{p}_{z}^{k}l}{4\sqrt{A_{z}^{k}} \right) \right].} \\ \text{$\times\left[\text{erf}\left(\frac{2A_{z}^{k}b-\tilde{p}_{z}^{k}l}{4\sqrt{A_{z}^{k}} \right) \right].} \\ \text{$\times\left[\text{erf}\left(\frac{2A_{z}^{k}b-\tilde{p}_{z}^{k}l}{4\sqrt{A_{z}^{k}} \right) \right].}$$

To obtain the number of particles in well i , we evaluate this expression at points a and b between the wells, we can write a = q i z -glyph[lscript]/ 2 and b = q i z + glyph[lscript]/ 2 with glyph[lscript] the distance between two wells. As in the case for the fourmode model, the main contribution of the sum results for the summand with i = k = l . We get approximately

$$n _ { i } = \sqrt { \frac { \pi ^ { 3 } } { 8 } } \frac { \text{erf} \left ( \sqrt { A _ { z, R } ^ { i } \ell ^ { 2 } / 2 } \right ) } { \sqrt { A _ { x, R } ^ { i } A _ { y, R } ^ { i } A _ { z, R } ^ { i } } } \left | d ^ { i } \right | ^ { 2 }.$$

For typical values of A i z,R glyph[lscript] 2 the error function is close to unity. Thus, it is reasonable to compare the particle numbers obtained from the four-mode model [Eq. (B3)] with those from the extended variational ansatz [Eq. (B9)]. Since the particle currents are given by the derivatives of the particle numbers, this is also the case for the particle currents.

- [1] C. E. Rüter, K. G. Makris, R. El-Ganainy, D. N. Christodoulides, M. Segev, and D. Kip, Nat. Phys. 6 , 192 (2010).
- [3] J. Schindler, A. Li, M. C. Zheng, F. M. Ellis, and T. Kottos, Phys. Rev. A 84 , 040101 (2011).
- [2] Y. D. Chong, L. Ge, and A. D. Stone, Phys. Rev. Lett. 106 , 093902 (2011).
- [4] S. Bittner, B. Dietz, U. Günther, H. L. Harney, M. MiskiOglu, A. Richter, and F. Schäfer, Phys. Rev. Lett. 108 , 024101 (2012).
- [6] C. M. Bender and S. Boettcher, Phys. Rev. Lett. 80 , 5243 (1998).
- [5] H. Ramezani, T. Kottos, R. El-Ganainy, and D. N. Christodoulides, Phys. Rev. A 82 , 043803 (2010).
- [7] S. Klaiman, U. Günther, and N. Moiseyev, Phys. Rev. Lett. 101 , 080402 (2008).
- [9] Y. Shin, G. B. Jo, M. Saba, T. A. Pasquini, W. Ketterle, and D. E. Pritchard, Phys. Rev. Lett. 95 , 170402 (2005).
- [8] D. Dast, D. Haag, H. Cartarius, G. Wunner, R. Eichler, and J. Main, Fortschritte der Physik 61 , 124 (2013).
- [10] H. Cartarius and G. Wunner, Phys. Rev. A 86 , 013612 (2012).
- [11] E. M. Graefe, H. J. Korsch, and A. E. Niederle, Phys.
- Rev. Lett. 101 , 150408 (2008).
- [13] A. Zenesini, C. Sias, H. Lignier, Y. Singh, D. Ciampini, O. Morsch, R. Mannella, E. Arimondo, A. Tomadin, and S. Wimberger, New J. Phys. 10 , 053038 (2008).
- [12] D. Haag, D. Dast, A. Löhle, H. Cartarius, J. Main, and G. Wunner, Phys. Rev. A 89 , 023601 (2014).
- [14] C. Elsen, K. Rapedius, D. Witthaut, and H. J. Korsch, J. Phys. B 44 , 225301 (2011).
- [16] K. Henderson, C. Ryu, C. MacCormick, and M. G. Boshier, New J. Phys. 11 , 043030 (2009).
- [15] M. Kreibich, J. Main, H. Cartarius, and G. Wunner, Phys. Rev. A 87 , 051601(R) (2013).
- [17] A. Trombettoni and A. Smerzi, Phys. Rev. Lett. 86 , 2353 (2001).
- [19] I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80 , 885 (2008).
- [18] A. Smerzi and A. Trombettoni, Phys. Rev. A 68 , 023613 (2003).
- [20] P. Löwdin, J. Chem. Phys. 18 , 365 (1950).
- [22] A. D. McLachlan, Mol. Phys. 8 , 39 (1964).
- [21] E.-M. Graefe, J. Phys. A 45 , 444015 (2012).
- [23] R. Eichler, D. Zajec, P. Köberle, J. Main, and G. Wunner, Phys. Rev. A 86 , 053611 (2012).
Manuel Kreibich Jörg Main Holger Cartarius Günter Wunner
Realizing PT-symmetric non-Hermiticity with ultra-cold atoms and Hermitian multi-well potentials
We discuss the possibility of realizing a non-Hermitian, i.e. an open
two-well system of ultra-cold atoms by enclosing it with additional
time-dependent wells that serve as particle reservoirs. With the appropriate
design of the additional wells PT-symmetric currents can be induced to and from
the inner wells, which support stable solutions. We show that interaction in
the mean-field limit does not destroy this property. As a first method we use a
simplified variational ansatz leading to a discrete nonlinear Schr\"odinger
equation. A more accurate and more general variational ansatz is then used to
confirm the results.