Realizing PT-symmetric non-Hermiticity with ultra-cold atoms and Hermitian multi-well potentials
RRealizing PT -symmetric non-Hermiticity with ultra-cold atoms and Hermitianmulti-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-coldatoms by enclosing it with additional time-dependent wells that serve as particle reservoirs. Withthe appropriate design of the additional wells PT -symmetric currents can be induced to and fromthe inner wells, which support stable solutions. We show that interaction in the mean-field limitdoes not destroy this property. As a first method we use a simplified variational ansatz leading toa discrete nonlinear Schrödinger equation. A more accurate and more general variational ansatz isthen used to confirm the results. I. INTRODUCTION
Since the first realization of PT -symmetric gain andloss in optical wave guides [1] much effort has been madeto realize analogous systems in various fields of physics,e. g. lasers [2], electronics [3], microwave cavities [4], andto make use of new effects arising from nonlinearity, viz.the Kerr nonlinearity in optical wave guides [5]. PT sym-metry originates from quantum mechanics and stands fora combined action of parity and time-reversal. Despitethe non-Hermiticity of PT -symmetric Hamiltonians, in acertain range of parameters entirely real eigenvalue spec-tra exist [6]. Due to a formal analogy between quantummechanics and electromagnetism, the formulation of PT symmetry has spread to those systems, leading to theabove mentioned realizations. However, the experimen-tal verification in a genuine quantum system has not beenachieved so far.In a PT -symmetric system there is a balanced gainand loss of probability density (or electromagnetic fieldin analogue systems). According to a proposal in [7] sucha quantum mechanical system could be realized with aBose-Einstein condensate (BEC) in a double-well poten-tial where particles are injected into one well and removedfrom the other. Indeed, it could be shown that the systemis ideally suited for a first-experimental observation of PT symmetry in a quantum system [8]. There is progressin coupling two BECs and at the same time eject parti-cles [9], but an actual realization of PT symmetry usinga BEC is still missing.On the other hand much theoretical work has beendone on PT -symmetric BECs, including the proof of ex-istence of stable states of interacting systems [10], a mi-croscopic treatment based on the Bose-Hubbard model[11], and a thorough investigation on dynamics of stableand unstable regimes [12]. In contrast to injecting and re-moving particles from a double-well potential, one couldthink of a double-well included in a tilted optical lattice,with an incoming and outgoing transport of particles, socalled Wannier-Stark systems [13, 14]. These incomingand outgoing particle currents could in principle serveas the necessary currents for realizing a PT -symmetricsystem.In a recent Rapid Communication [15] we followed the V ( x ) ψ ( x )Im V (a) V ( x ) ψ ( x ) (b) tunneling Figure 1. (a) Two-mode model with an imaginary potential V , where – in this case – particles are removed from the leftwell and injected into the right. (b) Such a system could berealized by coupling two additional wells to the system, wherethe tunneling currents to and from the outer wells are usedas an implementation of the imaginary potential. The outerwells act as (limited) particle reservoirs and have to be chosentime-dependent. simplest approach to this idea, namely to couple two ad-ditional wells to the double well potential (see Fig. 1),and use the arising currents as a theoretical realiza-tion of PT symmetry in a quantum mechanical system.We showed that for a specific time-dependent choice ofthe well depths and coupling strengths of the additionalwells, the inner wells behave exactly as the wells of a PT -symmetric two-well system, thus serving as a possi-ble experimental realization.It is the purpose of this paper investigate more deeplyour approach, especially for the case of interacting atoms,and to formulate a more accurate method to confirm theresults. In Sec. II we derive the few-mode model ordiscrete nonlinear Schrödinger equation (DNLSE) fromthe Gross-Pitaevskii equation (GPE) using a simple vari-ational ansatz, where each well is populated by a sin-gle time-independent Gaussian function. The resultsare applied in Sec. III to show that a Hermitian four-mode model can be chosen such that the middle wellsbehave exactly as the wells of the PT -symmetric two-mode model. This is exemplified for some parameters.In Sec. IV we use an extended variational ansatz for therealization of PT symmetry. Still one Gaussian functionper well is used, but each Gaussian function is fully dy-namical in its parameters, thus giving a more accuratedescription of the system. The method is compared with a r X i v : . [ qu a n t - ph ] A ug the simply Gaussian ansatz. In Sec. V we conclude ourwork. II. FROM THE GROSS-PITAEVSKIIEQUATION TO A FEW-MODE MODEL
In this section we apply the steps that lead from theGPE, which describes the dynamics of a BEC at abso-lute zero temperature, to a DNLSE or few-mode-model.At first we presents our ansatz and calculate the neces-sary integrals. To derive the equations of motion (EOM)in a Schrödinger form we use the method of symmetricorthogonalization. In Sec. III the results are used andapplied to different scenarios.
A. Simplified variational ansatz
We start with the well-known GPE, i (cid:126) ∂ t ψ ( r , t ) = (cid:20) − (cid:126) m ∆ + V ( r ) + g | ψ ( r , t ) | (cid:21) ψ ( r , t ) . (1)The interaction strength is given by g = 4 π (cid:126) N a/m withthe s-wave scattering length a and particle number N .For the external multi-well potential V ( r ) we use a Gaus-sian profile for each well, V ( r ) = N w (cid:88) k =1 V k exp (cid:20) − x w x − y w y − z − s kz ) w z (cid:21) , (2)with a total number of N w wells. Each well k has thepotential depth V k < and is displaced along the z -direction by s kz . Such a potential could, e. g., be createdexperimentally 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 de-grees of freedom. In general, an ansatz of the form ψ ( r , t ) = (cid:88) k d k ( t ) g k ( r ) , (3)is used, where d k ( t ) is the (complex) amplitude and g k ( r ) the localized wave function in well k . Usually, a set ofWannier functions is applied [19]. In this work [Sec. III],we use a single Gaussian function for each well. With thisansatz we can derive analytical formulae for the matrixelements of the DNLSE or few-mode model, and thusopen an easy way to actually compute the localized wavefunctions by a simple energy minimization process. InSec. IV we extend the variational ansatz to obtain moreaccurate results.The simplified variational ansatz is given by the su-perposition of N G Gaussian functions, one for each well, N w = N G , with g k ( r ) = e − A kx x − A ky y − A kz ( z − q kz ) . (4) The parameters A kx , A ky , A kz ∈ R describe the width ofeach Gaussian function in each direction, and q kz ∈ R thedisplacement along the z -direction. We treat these pa-rameters in this section as time-independent. The time-dependence of the wave function comes from the ampli-tude d k ( t ) in Eq. (3).We insert ansatz (3) into the GPE i (cid:126) ∂ t ψ = ˆ Hψ , multi-ply from the left with ( g l ) ∗ , and integrate over R . Then,we are left with i (cid:126) K ˙ d = H d . (5)Here, the amplitudes d k are written as a vector d =( d , . . . , d N G ) T ∈ C N G . The matrix K describes the over-lap, and H the Hamiltonian matrix elements. They aregiven by the integrals (see Sec. II B for their explicit cal-culation) K lk = (cid:10) g l (cid:12)(cid:12) g k (cid:11) , (6a) H lk = (cid:68) g l (cid:12)(cid:12)(cid:12) ˆ H (cid:12)(cid:12)(cid:12) g k (cid:69) . (6b)Our aim is to derive EOM for the amplitudes in aSchrödinger form, i. e. i (cid:126) ˙ d eff = H eff d eff , with possiblymodified amplitudes and matrix elements. Equation (5)would be of that form if K was proportional to the iden-tity matrix. Since the Gaussian functions g k are non-orthogonal this is not the case. For that reason, we per-form a symmetric orthogonalization in Sec. II C. B. Evaluation of the integrals
The occurring integrals in Eqs. (6) are solely based onthe Gaussian integral ∞ (cid:90) −∞ d x x n e − Ax + px = ∂ n ∂p n ∞ (cid:90) −∞ d x e − Ax + px = ∂ n ∂p n (cid:114) πA e p / A , Re A > . (7)For convenience we define the abbreviations A klα = A kα + ( A lα ) ∗ , α = x, y, z, (8a) κ klα = A kα ( A lα ) ∗ /A klα , (8b) β klα = (cid:115) A klα w α A klα w α + 2 , (8c) c kl = e − κ klz ( q kz − q lz ) . (8d)The integrals of K can now easily be calculated withEq. (7) for the case n = 0 , the result is K lk = (cid:114) πA klx (cid:114) πA kly (cid:114) πA klz c kl . (9)For the calculation of the Hamiltonian matrix H , we sep-arate 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 formshown in Eq. (7) with n = 2 and we obtain T lk = K lk (cid:2) κ klx + κ kly + κ klz − κ klz ) ( q kz − q lz ) (cid:3) . (10)To calculate the integrals necessary for the external po-tential, the parameters A and p in Eq. (7) are shifted bythe parameters of the external potential. The result is V lk = K lk β klx β kly β klz (cid:88) m V m × exp (cid:40) − (cid:2) A kz ( s mz − q kz ) + ( A lz ) ∗ ( s mz − q lz ) (cid:3) A klz ( A klz w z + 2) (cid:41) . (11)Since the operator of the interaction potential W is non-linear, its matrix elements depend on the amplitudes d k .With the standard Gaussian integral we obtain W lk = (cid:88) i,j ˜ W lkji ( d j ) ∗ d i , (12)where we defined the four-rank tensor ˜ W lkji = 4 π (cid:126) N am (cid:114) πA ijklx (cid:115) πA ijkly (cid:114) πA ijklz × exp (cid:20) − A iz ( A jz ) ∗ ( q iz − q jz ) + A iz ( A lz ) ∗ ( q iz − q lz ) A ijklz (cid:21) × exp (cid:20) − A kz ( A jz ) ∗ ( q kz − q jz ) + A kz ( A lz ) ∗ ( q kz − q lz ) A ijklz (cid:21) × exp (cid:20) − A iz A kz ( q iz − q kz ) + ( A jz ) ∗ ( A lz ) ∗ ( q jz − q lz ) A ijklz (cid:21) (13)and the new abbreviations A ijklα = A ijα + A klα , α = x, y, z. (14) C. Symmetric orthogonalization
We calculated all necessary integrals for setting up theEOMs (5), which are not of the form of a Schrödingerequation (cf. Sec. II A). To achieve this, we now trans-form the EOM into a Schrödinger form with the methodof symmetric orthogonalization [20]: Since K is Hermi-tian there exists a unitary transformation U such that D = UKU † is diagonal with real entries. Then we canconstruct the Hermitian matrix X = U † D − / U . (15)This matrix has the properties that XKX = , and fur-thermore, if H is Hermitian, so is XHX . Now, we canuse this matrix to transform Eq. (5) to i (cid:126) ( XKX )( X − ˙ d ) = ( XHX )( X − d ) . (16) With the definitions d eff = X − d and H eff = XHX , andthe above property of X we arrive at the EOM for d eff in Schrödinger form i (cid:126) ˙ d eff = H eff d eff (17)with a Hermitian Hamiltonian matrix H eff (if H is Her-mitian). We note, that the normalization condition forthe wave function ψ transforms as ! = (cid:90) d r | ψ | = d † K d = ( X − d ) † ( XKX )( X − d )= d † eff d eff . (18) D. Application of symmetric orthogonalization andnearest-neighbor approximation
By means of symmetric orthogonalization we cantransform the EOM (5) into a Schrödinger form since weknow K analytically. Without any further approxima-tions the analytical expressions can become quite com-plicated. At this point, we introduce as usual the nearestneighbor approximation to obtain analytical and simpleapproximate 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 (0) lk = (cid:16) π (cid:17) / δ kl (cid:113) A kx,R A ky,R A kz,R (19)and K (1) lk = (cid:114) πA klx (cid:114) πA kly (cid:114) πA klz c kl ( δ k,l +1 + δ k +1 ,l ) , (20)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 (0) lk = (cid:18) π (cid:19) / (cid:113) A kx,R A ky,R A kz,R δ kl . (21)To calculate the first order we start with (cid:16) X (0) + X (1) (cid:17) (cid:16) K (0) + K (1) (cid:17) (cid:16) X (0) + X (1) (cid:17) = , (22)ignore second order terms, insert the known quantities,and solve for X (1) , which yields X (1) lk = − (cid:18) π (cid:19) / c kl (cid:112) A klx (cid:113) A kly (cid:112) A klz ( δ k,l +1 + δ k +1 ,l ) × (cid:113) A kx,R A ky,R A kz,R A lx,R A ly,R A lz,R (cid:113) A kx,R A ky,R A kz,R + (cid:113) A lx,R A ly,R A lz,R . (23)Now we are able to calculate analytical approximationsfor the transformed Hamiltonian H eff . We observe thatfor the linear parts, i. e. the kinetic and external potentialpart, the matrix elements of H are proportional to K , sowe can write H lin lk = K lk h 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 Hamiltonianwe obtain H lin , (0) eff ,lk = (cid:88) m,n X (0) lm K (0) mn X (0) nk h lin mn = h lin kk δ kl . (24)The zeroth order contributes to the linear diagonal ele-ments of H lineff , so they can be identified with the onsiteenergy E k = h lin kk . With Eqs. (10) and (11) we can write E k = (cid:126) m (cid:0) A kx,R + A ky,R + A kz,R (cid:1) + V k β kkx β kky β kkz exp (cid:2) − β kkz ) ( s kz − q kz ) /w z (cid:3) , (25)where we considered only the term with m = k in thesum in Eq. (11) to be consistent with zeroth order.The first order terms of H lineff are given by the first orderterms of the expression ( X (0) + X (1) )( H lin , (0) + H lin , (1) )( X (0) + X (1) ) , (26)which gives, after inserting all known quantities, H lin , (1) eff ,lk = − √ (cid:113) A kx,R A ky,R A kz,R A lx,R A ly,R A lz,R (cid:113) A kx,R A ky,R A kz,R + (cid:113) A lx,R A ly,R A lz,R × E k − h lk (cid:113) A kx,R A ky,R A kz,R + E l − h lk (cid:113) A lx,R A ly,R A lz,R × c kl (cid:112) A klx (cid:113) A kly (cid:112) A klz ( δ k,l +1 + δ k +1 ,l ) . (27)These elements can be identified as the tunneling ele-ments of the few-mode model, i. e. J lk = − H lin , (1) eff ,lk , sincethey are the super- and sub-diagonal entries. The quan-tity h lin lk is given by Eqs. (10) and (11), where in v lk weonly consider nearest-neighbor contributions, v lk = β klx β kly β klz × (cid:32) V k exp (cid:40) − (cid:2) A kz ( s kz − q kz ) + ( A lz ) ∗ ( s kz − q lz ) (cid:3) A klz ( A klz w z + 2) (cid:41) + V l exp (cid:40) − (cid:2) A kz ( s lz − q kz ) + ( A lz ) ∗ ( s lz − q lz ) (cid:3) A klz ( A klz w z + 2) (cid:41)(cid:33) . (28)So far we have calculated the matrix elements of thelinear transformed Hamiltonian H lineff . We now turn ourattention 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, weobtain W eff ,lk = (cid:88) i,j,m,n,p,q X lm ˜ W mnji X nk ( X jp d p eff ) ∗ ( X iq d q eff )= (cid:88) p,q (cid:88) i,j,m,n X lm X pj ˜ W mnji X nk X iq (cid:124) (cid:123)(cid:122) (cid:125) ≡ ˜ W eff ,lkpq ( d p eff ) ∗ d q eff . (29)For the interaction term, we only consider onsite-elements, i. e. elements with i = j = k = l in ˜ W lkji , sincemainly these terms contribute to the interaction part.For the transformed four-rank tensor, we then obtain ˜ W eff ,kkkk = 4 (cid:126) N a √ πm (cid:113) A kx,R A ky,R A kz,R . (30)These elements are the nonlinear coupling constants, c k = ˜ W eff ,kkkk . The whole action of the transformedHamiltonian then reads W eff ,lk = c k | d k | δ kl . (31)We have all necessary matrix-elements in nearest-neighbor approximation. Using this knowledge we cancalculate these elements from the Gaussian variationalapproach, or – vice verse – determine the parametersof the realistic potential (2) by knowing the matrix el-ements. As an application this is done in Sec. III. Theparameters A kα and q kz can be computed by minimizingthe mean-field energy of the system for a given initialcondition with E mf = (cid:88) k,l ( d l ) ∗ ( T lk + V lk ) d k + 12 (cid:88) i,j,k,l ( d l ) ∗ d k W lkji ( d j ) ∗ d i . (32)For this minimization, all parameters A kα , q kz and d k haveto be varied with the norm constraint (18), which is com-putationally much cheaper than a grid calculation. Theresults of this section could also be used to calculate theparameters of the Bose-Hubbard model.Now that we have all matrix elements, we can begin ouractual investigation of a realization of a non-Hermitiansystem, first by means of few-mode models. In the fol-lowing we will adopt the usual notation and write for thediscrete wave function ψ k instead of d k . III. FEW-MODE MODELSA. Equivalence of the two- and four-mode models
We now return to the idea mentioned in the intro-duction and sketched in Fig. 1, i. e. the embedding of anon-Hermitian double-well with a closed four-well struc-ture. Before analyzing the possibility of realizing a non-Hermitian system with a Hermitian model, we discuss thebasic properties of the non-Hermitian two-mode model,which is given by H = (cid:18) E + c | ψ | − J − J E + c | ψ | (cid:19) (33)with complex onsite energies E k ∈ C , which make theHamiltonian non-Hermitian. This system has intensivelybeen investigated in [21]. The time derivatives of the ob-servables, that is particle number n k = ψ ∗ k ψ k and particlecurrent j kl = i J kl ( ψ k ψ ∗ l − ψ ∗ k ψ l ) , can be easily calculated,which gives ∂ t n = − j + 2 n Im E , (34a) ∂ t n = j + 2 n Im E , (34b) ∂ t j = 2 J ( n − n ) + (Im E + Im E ) j + J (Re E − Re E + c n − c n ) C , (34c) ∂ t C = (Im E + Im E ) C − (Re E − Re E + c n − c n )˜ j , (34d)where we defined C kl = ψ k ψ l + ψ ∗ k ψ ∗ l and ˜ j kl = i( ψ k ψ ∗ l − ψ ∗ k ψ l ) = j /J . In Eqs. (34a) and (34b) we see thatthe imaginary parts of the onsite energies act as addi-tional sources or sinks to the particle flow, the quantities j e = 2 n Im E and j e = − n Im E can be identi-fied as currents to and from the environment, thus, thisHamiltonian describes an open quantum system.The general solutions for the nonlinear model (33) arecalculated in [21]. For convinience we only discuss theresults of the linear, i. e. non-interacting case ( c = c =0 ). Then, the eigenvalues are given by E ± = E + E ± (cid:114) J + ( E − E ) . (35)Requirement for PT symmetry leads to the condition E = E ∗ , which is for instance fulfilled by E = iΓ and E = − iΓ , a situation with balanced gain and loss. Then,the eigenvalues are E ± = (cid:112) J − Γ . For Γ < J theyare real, whereas for Γ > J the PT symmetry is brokenand the eigenvalues are purely imaginary.It is now our main purpose to investigate whether thebehavior of the non-Hermitian two-mode model (33) canbe described by a Hermitian four-mode model, which isgiven by H ( t ) = E ( t ) − J ( t ) 0 0 − J ( t ) E − J − J E − J ( t )0 0 − J ( t ) E ( t ) + diag (cid:16) c | ψ | , c | ψ | , c | ψ | , c | ψ | (cid:17) . (36)Two additional wells are coupled to the system, they actas a particle reservoir for the inner wells. The onsite energies and tunneling elements of the outer wells maybe time-dependent. Since this model shall be Hermitian,the onsite energies are real, E k ∈ R .To derive a relationship between the two- and four-mode model we calculate the time derivatives of the sameobservables as in the two-mode model, which yields ∂ t n = j − j , (37a) ∂ t n = j − j , (37b) ∂ t j = 2 J ( n − n ) − J ( J C − J C )+ J ( E − E + c n − c n ) C , (37c) ∂ t C = J ˜ j − J ˜ j − ( E − E + c n − c n ) ˜ j . (37d)Comparing with Eqs. (34) we find that for the two-modemodel the following condition must hold, Im E = − Im E ⇒ E = E ∗ . (38)This means that only the PT -symmetric two-modemodel, as discussed above, may be simulated by thefour-mode model. From now on we write E = iΓ and E = − iΓ with Γ ∈ R . Furthermore, the real parts ofthe onsite energies E and E and the coupling element J have to agree between the two models. We are thenleft with the following conditions, j = 2Γ n , (39a) j = 2Γ n , (39b) J C − J C , (39c) J ˜ j − J ˜ j . (39d)To summarize at this point, if these conditions are ful-filled at every time, the two middle wells of the Hermi-tian four-mode model behave exactly as the wells of the PT -symmetric two-mode model. These four conditionsseem to be independent, but we shall see in App. A thatEq. (39d) follows if Eqs. (39a)–(39c) are fulfilled.We now have to make clear that these conditions canindeed be fulfilled by giving suitable time-dependenciesto the elements of the four-mode Hamiltonian. Equa-tion (39c) can simply be fulfilled by choosing J ( t ) = dC ( t ) , J ( t ) = dC ( t ) , (40)where d (cid:54) = 0 is a real, time-independent quantity. Thetunneling elements are then time-dependent. With J and J determined there are no free parameters left tofulfill Eqs. (39a) and (39b). Instead, we can set the time-derivatives of these equations with the help of the remain-ing onsite energies E and E . We calculate the time-derivatives of j and j using the Schrödinger equation,which yields ∂ t j = ( ∂ t J ) ˜ j + 2 J ( n − n ) + J J C + J ( E − E + c n − c n ) C = ∂ t j tar , (41a) ∂ t j = ( ∂ t J ) ˜ j + 2 J ( n − n ) − J J C + J ( E − E + c n − c n ) C = ∂ t j tar , (41b)with the target currents determined by Eqs. (39a)and (39b). The time derivatives of the tunneling ele-ments J and J are given by ∂ t J = d∂ t C = d (cid:0) J ˜ j + J ˜ j − J ˜ j (cid:1) + d ( E − E + c n − c n ) ˜ j , (42a) ∂ t J = d∂ t C = d (cid:0) J ˜ j − J ˜ j − J ˜ j (cid:1) + d ( E − E + c n − c n ) ˜ j . (42b)Inserting this into Eqs. (41), the onsite energies E and E are determined by a linear set of equations, (cid:18) C C ˜ j ˜ j − ˜ j ˜ j − C C (cid:19) (cid:18) E E (cid:19) = (cid:18) v v (cid:19) , (43)with the entries v = ∂ t j tar /d − J ˜ j − d (cid:0) C ˜ j − C ˜ j (cid:1) ˜ j − dC ( n − n ) − J C C − ( c n − c n ) ˜ j ˜ j − ( c n − c n ) C C , (44a) v = ∂ t j tar /d + J ˜ j − d (cid:0) C ˜ j − C ˜ j (cid:1) ˜ j − dC ( n − n ) + J C C − ( c n − c n ) ˜ j ˜ j − ( c n − c n ) C C , (44b)where we set E = E = 0 , since – as discussed above– the real parts of these elements have to agree betweenthe two- and four-mode model.We have to note that by means of the onsite energies E and E we do not fix the currents, but their time-derivatives. For this reason the initial probabilities andcurrents have to be chosen such that they fulfill the condi-tions (39). In the following section we apply these resultsto obtain solutions for different starting conditions. B. Results
Before presenting our results we give the necessary ini-tial conditions. We insert the definition of the particlecurrent and the solutions for the tunneling elements (40) n n (b)0.40.50.6 n , n (a)0.40.50.6 0 10 20 30 40 t [¯ h/J ] j , j , j (c) -8-6-4-202 0 10 20 30 40 t [¯ h/J ] J , J E , E (d) Figure 2. Realization of a stationary PT -symmetric solutionwithin the four-mode model. (a) Particles numbers in well 1(solid) and well 2 (dashed), both equal to / . (b) Particlenumbers in well 0 (solid) and well 3 (dashed). (c) Particle cur-rents j (solid), j (dashed) and j (dash-dotted) in unitsof J / (cid:126) . (d) Tunneling elements J (solid), J (dashed),and onsite energies E (solid), E (dashed), in units of J . into the two Eqs. (39a) and (39b) and express everythingin terms of the wave function ψ k . We are allowed tochoose the global phase of the initial wave function, sowe choose ψ I = 0 . We assume the real parts ψ R and ψ R to be arbitrary but fixed. We are left with (after simplerearrangement) ψ I = Γ2 dψ R , (45a) ψ I = ψ R ψ I ψ R − Γ n dψ R (cid:0) ψ R ψ R + ψ I ψ I (cid:1) . (45b)Now, the wave function ψ k fulfills the necessary condi-tions and with the time-dependencies of the matrix ele-ments from Sec. III A we can simulate the PT -symmetrictwo-mode model.First we consider a stationary solution in the non-interacting case [cf. Eq. (35)] for Γ /J = 0 . . Fig. 2shows the results. As expected the particle numbers inthe two middle wells are both equal and constant in time(with normalization n + n = 1 ), as well as the particlecurrents. For this reason, the number of particles de-creases (left well) or increases linearly (right well), withthe slope ˙ n = − Γ / (cid:126) and ˙ n = Γ / (cid:126) , respectively. Thematrix elements of the four-mode model vary slightly intime. Thus, we have shown that the conditions (39) canbe fulfilled for a finite time.As a second example we prepare a non-stationary so-lution at t = 0 , with ψ (0) = √ . and ψ (0) = √ . (see Fig. 3). In this case the particle numbers in thetwo middle wells oscillate in time with a phase differ-ence ∆ φ < π leading to a non-constant added numberof particles, which is a typical feature of PT -symmetric n n (b)00.511.52 n , n n + n (a)00.511.5 0 5 10 15 20 25 t [¯ h/J ] j , j , j (c) -15-10-50 0 5 10 15 20 25 t [¯ h/J ] J , J E , E (d) Figure 3. Same as Fig. 2, but for oscillatory dynamics. Ad-ditionally, the added number of particles n + n is plottedas crosses in (a). As the particle number in the reservoirwell 0 decreases and gets close to zero, the conditions cannotbe fulfilled anymore and the simulation breaks down. n n (b)024681012 n n n + n (a)-40481216 0 1 2 3 4 5 t [¯ h/J ] (c) j , j , j J , J -15015 0 1 2 3 4 5 t [¯ h/J ] E , E (d) Figure 4. Same as Fig. 2, but for the interacting case ( c = − ), which leads to a collapsing dynamics, indicating a changeof stability. The reservoir is emptying in a short time domainsuch that the time available for PT symmetry is quite short. systems. Since well 0 acts as a particle reservoir the num-ber of particles decreases and gets close to zero up to thepoint where the linear system of equations (43) cannotbe fulfilled anymore (determinant of coefficient matrixequals zero). This shows that the time available for asimulation of PT symmetry with a reservoir is limited.The matrix elements in this case show a quasi-oscillatorybehavior.Next, we consider a system with attractive interaction c = − . In such a system it is known that the groundstate changes its stability in an additional bifurcation (forrepulsive interaction, the excited state changes its stabil-ity) [12]. This phenomenon is highly correlated to the n n (b)450452454456458460 (a) n n t [ t ] j , j , j (c) -65-60-55 V V (d)-0.10.00.1 0 100 200 300 t [ t ] δ δ Figure 5. Same as Fig. 2, but for an adiabatic current ramp( Γ f /J = 0 . , t f /t = 60 ) for realistic parameters of the exter-nal potential. Instead of the matrix elements, in (d) we showthe parameters of the potential, δ k marks the displacementfrom equidistant potential wells. There is a total number ofparticles of N = 10 . See text for units. occurrence of self-trapping states. In our example × system [Eq. (33)] the ground state is already unstable forthe given interaction strength, so that the ground statewith a small perturbation is expected to collapse. Thecalculation is shown in Fig. 4. The particle number n in-creases exponentially, which marks the beginning of col-lapse. To support the exponential increase, the particlenumber in the reservoir n is decreasing exponentially,so that the simulation breaks down after a quite shorttime interval. The matrix elements now show a rathercomplex time-dependency.As the last example we calculate an experimentallymore realistic scenario. Up to now – as mentioned above– we had to start with appropriate wave functions inorder to obtain PT symmetry. This is nearly impossiblein 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-modemodel (with Γ = 0 ) and increase Γ to its target value Γ f .If this change is slow enough, we arrive approximatelyat the PT -symmetric ground state. The specific timedependence we use is Γ( t ) = Γ f [1 − cos( πt/t f )] / for t ∈ [0 , t f ] .In the following we simulate a condensate of N = 10 atoms of Rb. We use units based on the potential widthin z -direction, w z , which we set w z = 1 µ m . The basicunity of energy is E = (cid:126) /mw z , which yields E /h =116 Hz . The basic unit of time is t = mw z / (cid:126) = 1 .
37 ms .The wells have widths of w x = w y = 4 w z = 4 µ m , andare initially positioned equidistantly with a distance of . w z = 1 . µ m . The potential depths are given by V = V = − E and V = V = − E . We usea scattering length of a = 2 . a B with a B being the Bohrradius.We calculate the adiabatic ramp for Γ f /J = 0 . and t f /t = 60 . The results are shown in Fig. 5. The to-tal particle number and currents are scaled to a particlenumber of N = 10 , the currents are given in units of /t , the potential depths in units of E . Instead of the positions of the potential wells we plot the difference ofthe position from its initial value, δ k = s kz − s kz ( t = 0) .By means of the relations derived in Sec. II we can trans-form from the calculated matrix elements to parametersof the external potential.After time t f an approximate PT -symmetric state isobtained. This can be seen by the constant number ofparticles in the middle wells, despite an increasing anddecreasing number in wells 3 and 0, respectively, whichis a clear sign of PT symmetry, also in a possible exper-iment. The potential depths of the outer wells have tobe adjusted in a simple way, the positions vary within afew percent of the distance of the wells. This shows thatthe simple four-mode model can be used to describe arealistic scenario and to obtain the necessary parametersof the external potential.We have applied our results from Sec. III A to differentinitial conditions and system parameters and have shownthat it is indeed possible to realize a PT -symmetric sys-tem with a Hermitian four-mode model. In the followingSection we investigate this scenario in the framework ofthe GPE and a variational ansatz with Gaussian func-tions. In particular we verify that the results from thissection agree. IV. VARIATIONAL ANSATZ WITH GAUSSIANFUNCTIONSA. Variational ansatz and equations of motion
After the simple approach in the previous section wenow investigate the possible realization in the frameworkof the GPE with an extended variational ansatz. We usea superposition of Gaussian functions, where – in con-trast to the simpler approach of Sec. II A – all variationalparameters are considered as time-dependent, ψ = (cid:88) k e − ( r − q k ) T A k ( r − q k )+i p k ( r − q k ) − γ k . (46)The matrices A k = diag( A kx , A ky , A kz ) ∈ C × describethe width of each Gaussian function in each direction,the vectors q k = (0 , , q kz ) T ∈ R and p k = (0 , , p kz ) T ∈ R the position and momentum, respectively, and thescalar quantities γ k ∈ C the amplitude and phase. Thisvariational ansatz can describe a much richer dynamicsthan ansatz (4) and thus we expect more accurate results.The equations of motion follow from the time-dependent variational principle in the formulation ofMcLachlan [22], which states that the quantity I = || i (cid:126) φ − ˆ Hψ || (47) V ext n n n n j j j x z Figure 6. External potential V ext ; the darker, the deeper thepotential. The wells are separated into regions by walls atpositions right in the middle between two minima. The in-tegral over | ψ | in a region k is the particle number n k ; thecurrent through a wall is denoted by j k,k +1 . See text for ex-act definitions. With this separation we obtain a formulationanalogous to the four-mode model (Sec. III). shall be minimized with respect to φ and φ = ˙ ψ is setafterwards. Details of the necessary calculation can befound in Ref. [23], where the ansatz (46) has been usedto study the collision of anisotropic solitons in a BEC.The EOM are given by i (cid:126) (cid:88) k (cid:28) ∂ψ∂z l (cid:12)(cid:12)(cid:12)(cid:12) ∂ψ∂z k (cid:29) ˙ z k = (cid:28) ∂ψ∂z l (cid:12)(cid:12)(cid:12)(cid:12) ˆ H (cid:12)(cid:12)(cid:12)(cid:12) ψ (cid:29) , (48)where z = ( A x , . . . , γ N G ) T is the vector of all varia-tional parameters. Stationary solutions are given by fixedpoints of (48), whereas the dynamics can be calculatedby integrating Eq. (48) with a numerical integrator.In Sec. III we started our considerations by calculatingthe time-derivatives of the observables. In the case ofthe continuous GPE with a non-Hermitian potential thisleads to the modified equation of continuity ˙ ρ + div j = 2 ρ Im V, (49)with ρ = | ψ | and j = (cid:126) ( ψ ∗ ∇ ψ − ψ ∇ ψ ∗ ) / m i . The imag-inary part of the potential acts as an additional sourceor sink to the probability flow. Instead of adjusting thematrix elements, we need to adjust the parameters of theexternal potential. Since Eq. (49) is given on the wholespace R , one would need to change the external poten-tial 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. IIIas a roadmap for this continuous system. For that, weseparate the individual wells by walls between two min-ima of the potential at positions z k = ( s kz + s k +1 z ) / for k = 1 , , , z → −∞ and z → ∞ (see Fig. 6). Toobtain time derivatives of discrete observables as in thefew-mode model, we integrate the continuity equation forthe Hermitian four-well potential over a volume V k = { ( x, y, z ) T | − ∞ < x, y < ∞ , z k − < z < z k } . (50)This yields ∂ t (cid:90) (cid:90) (cid:90) V k d r ρ = − (cid:90) (cid:90) (cid:90) V k d r div j = − (cid:90) (cid:90) ∂V k d A · j . (51)The integral on the left-hand side can be interpreted asthe number of particles in well k . On the right-handside we used Gauss’s theorem. The surface integrals goesover a box, where the areas in the x - z - and y - z -plane goto infinity. We assume the current j to also vanish atinfinity, so only the integral over the two surfaces in the x - y -plane contribute. Thus, we can write ∂ t n k = − ∞ (cid:90) −∞ d x ∞ (cid:90) −∞ d y [ j z ( x, y, z k ) − j z ( x, y, z k − )] . (52)The integral over the z -component of j evaluated at z k represents the current from well k to k +1 , which we writeas j k,k +1 , for z k − this represents the current from well k − to k , j k − ,k . We finally obtain from the continuityequation ∂ t n k = j k − ,k − j k,k +1 , (53)where j − , = j , = 0 .Therefore, with Eq. (53) we have an analogous equa-tion as for the few-mode-model (37). We can use thesame method to realize a PT -symmetric system. At ev-ery time step of the numerical integration of the equa-tions of motion (48) we vary the parameters of the ex-ternal potential such that the observables show the de-sired behavior. This is done via a nonlinear root search,in contrast to the analytical solutions of Sec. III A. InSec. III B we showed that the positions of the outer wellsbarely needed to be adjusted. As a simplification we onlyvary the depths of the outer wells V , V and demandthe outer currents j , j to reach the desired values. Weare left with a two-dimensional root search. In the follow-ing section we determine the realization of PT symmetrywithin the GPE and compare the results with those ofthe few-mode-model. B. Results and comparison with few-mode models
We use the variational ansatz to simulate the adia-batic 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 . Fig. 7 showsthe results. It is not a priori clear whether the quan-tity n k = ( d k eff ) ∗ d k eff – resulting from the transformed am-plitudes of the simplified variational approach – can becompared with the integrated probability density n k of n n n n -62-61-60-59-58 (c) V V t [ t ] (d) j Figure 7. Adiabatic current ramp calculated with variationalansatz (lines) and four-mode model (crosses) for the samesystem. Shown are (a) particle numbers in the middle wells,(b) particle numbers in outer wells, (c) potential depths ofthe outer wells in units of E , and (d) particle current j .There are small deviations, but overall the four-mode modelis a valid approximation compared to the variational ansatz. the extended variational approach. We show in App. Bthat this comparison is indeed possible and reasonable.Overall, there is the same qualitative behavior betweenthe two approaches. The particle number in the middlewells is slightly under estimated in the four-mode modeldue to its origin as an approximation. Because of thefixed positions of the outer wells for the variational ap-proach, n and n are not exactly equal, but the differ-ence is small compared to the absolute number. We canconclude that we can now also create an (approximate) PT -symmetric state within the variational approach andverify that the results from the four-mode model are agood approximation.The variational approach gives us access to the wavefunction at each time step (see Fig. 8). As indicated bythe particle numbers in Fig. 7, the amplitudes in the mid-dle wells stay almost constant, whereas the amplitudes inthe outer wells increase or decrease in time. This couldalso be measured in an experiment and could serve as aproof of realized PT symmetry.0 t = 0 z [ w z ] t = 80 t t = 120 t -2.7 -0.9 0.9 2.7 z [ w z ] t = 170 t Figure 8. Modulus squared of the wave function | ψ | at x, y =0 for different times t . The amplitudes of the outer partsincrease or decrease, whereas the amplitudes in the middlewells stay constant, a clear sign of PT symmetry, also in anexperiment. V. CONCLUSION
By means of a simple variational ansatz – with lo-calized time-independent Gaussian functions – and themethod of symmetric orthogonalization we could trans-form the GPE for a multi-well potential to a simple few-mode model. The parameters of the variational ansatzare obtained by an energy minimization process in a com-putationally cheap way. With the results one can, on theone hand, calculate the elements of the few-mode modelfor a specific system, and on the other determine the pa-rameters of the external potential by knowing the time-dependent matrix elements.These results are applied to the PT -symmetric two-mode model and to the Hermitian four-mode model.With this model we could show that the four-mode modelcan be designed in such a way such that the middle wellsbehave exactly as the two-wells of the PT -symmetrictwo-mode model. Thus, this Hermitian (closed) systemcan be used as a possible realization of a quantum me-chanical PT -symmetric system. With an extended Gaus-sian variational ansatz – with full time-dependent Gaus-sian functions, which can describe a much richer dynam-ics – we confirmed the qualitative results of the few-modemodel.In this paper we applied the idea to embed a two-wellsystem into a closed system for the simplest case. Itis possible for future work to use a greater embeddingwhich could result in a simpler time-dependence of theadditional parameters. The ideas presented in this papercould also serve as a possible road map for an investi-gation of the Bose-Hubbard model and its many-particleeffects. ACKNOWLEDGMENTS
This work was supported by DFG. M. K. is gratefulfor support from the Landesgraduiertenförderung of theLand Baden-Württemberg.
Appendix A: Analytical considerations
In Sec. III A we derived four conditions [Eqs. (39)] suchthat the Hermitian four-mode model can simulate the PT -symmetric two-mode model. However, only threeof them are independent, which will be shown in thisappendix. For that, we need the relations C jj ˜ j ik = C ij ˜ j jk + ˜ j ij C jk , (A1a) C jj C ik = C ij C jk − ˜ j ij ˜ j jk , (A1b)which can be proved by simply inserting the definitionsof these quantities.We start with Eqs. (39a) and (39b), insert the solu-tions (40), take the square and use Eqs. (A1). We get ˜ j ˜ j + 2Γ d ˜ j ˜ j ˜ j − C n n ˜ j + 4Γ d n n = 0 , (A2a) ˜ j ˜ j + 2Γ d ˜ j ˜ j ˜ j − C n n ˜ j + 4Γ d n n = 0 . (A2b)Taking the difference yields ˜ j = n n n n ˜ j ⇔ ˜ j = s (cid:114) n n n n ˜ j , (A3)where s = ± gives the sign of ˜ j compared to ˜ j .Inserting this result into Eq. (A2a) gives a fourth-orderpolynomial equation for ˜ j . The result is ˜ j = s (cid:114) n n (cid:16) − α + s (cid:112) (1 − α ) − β (cid:17) , (A4)with the definitions α = γ (cid:16) β + γ (cid:17) , (A5a) β = s Γ d √ n n , (A5b) γ = ˜ j √ n n . (A5c)Using Eqs. (A1) we can find all other quantities ˜ j kl and C kl , especially those needed for the investigation of the1fourth condition J ˜ j − J ˜ j = dC ˜ j − dC ˜ j , C = s sign d (cid:114) n n (cid:16) − α − s (cid:112) (1 − α ) − β (cid:17) , (A6a) C = s sign d (cid:114) n n (cid:16) − α − s (cid:112) (1 − α ) − β (cid:17) , (A6b) ˜ j = s (cid:114) n n (cid:16) α + s (cid:112) (1 − α ) − β (cid:17) , (A6c) ˜ j = s (cid:114) n n (cid:16) α + s (cid:112) (1 − α ) − β (cid:17) . (A6d)All signs s i are determined by the phases of the initialwave function. We can now calculate the fourth conditionand find J ˜ j − J ˜ j = 0 . (A7)Thus, the first three conditions of Eqs. (39) imply thevalidity of the fourth one. With the results of the ma-trix elements of Sec. III A, namely J , J , E and E ,we find an exact equivalence of the two- and four-modemodels also for interacting atoms. Furthermore, with theresults of this appendix, we can calculate these matrix el-ements, once the quantities n , n and ˜ j are known. Appendix B: Comparison of probabilities forfew-mode-model and Gaussian functions
As discussed in Sec. IV B it is not a priori clearwhether we can compare the particle numbers obtainedfrom the four-mode model and the extended variationalapproach. In this appendix we take a closer look at bothquantities.The transformation of the amplitudes from the four-mode model, resulting from symmetric orthogonaliza-tion, is given by the inverse of the matrix X (cf. Sec. II C), d l eff = (cid:88) k ( X − ) lk d k = (cid:88) k (cid:104) ( X − ) (0) lk + ( X − ) (1) lk (cid:105) d k . (B1)The individual orders of the inverse matrix can easily becalculated, we get ( X − ) (0) lk = (cid:114) π (cid:113) A kx,R A ky,R A kz,R δ lk , ( X − ) (1) lk = √ π (cid:113) A kx,R A ky,R A kz,R A lx,R A ly,R A lz,R (cid:113) A kx,R A ky,R A kz,R + (cid:113) A lx,R A ly,R A lz,R × c kl (cid:112) A klx (cid:113) A kly (cid:112) A klz ( δ k,l − + δ k,l +1 ) . (B2a) The main contribution to the particle number n k =( d k eff ) ∗ d k eff is then given by n k = (cid:114) π (cid:113) A kx,R A ky,R A kz,R (cid:12)(cid:12) d k (cid:12)(cid:12) , (B3)plus terms involving neighboring amplitudes.The extended variational ansatz is given by ψ = (cid:88) k d k e − A kx x − A ky y − A kz ( z − q kz ) +i( z − q kz ) p kz , (B4)with the amplitude d k = exp( − γ k ) . The probability den-sity can then be written as ρ = (cid:88) k,l d k ( d l ) ∗ e − A klx x − A kly y − A klz z +i˜ p klz z − ˜ γ kl , (B5)with the definitions ˜ p klz = 2[ A kz q kz + ( A lz ) ∗ q lz ] + i( p kz − p lz ) , (B6a) ˜ γ kl = A kz ( q kz ) + ( A lz ) ∗ ( q lz ) + i( q kz p kz − q lz p lz ) . (B6b)We integrate ρ over the region discussed in Sec. IV A.The integral over x and y yields ρ z = (cid:88) k,l d k ( d l ) ∗ π (cid:112) A klx (cid:113) A kly e − A klz z +˜ p klz z − ˜ γ kl . (B7)The integral over z goes over a finite interval, say a ≤ z ≤ b , we can express the result in terms of the errorfunction, n ( a, b ) = 12 (cid:88) k,l d k ( d l ) ∗ π (cid:112) A klx (cid:113) A kly (cid:112) A klz e (˜ p klz ) / A klz − ˜ γ kl × (cid:34) erf (cid:32) A klz b − ˜ p klz (cid:112) A klz (cid:33) − erf (cid:32) A klz a − ˜ p klz (cid:112) A klz (cid:33)(cid:35) . (B8)To obtain the number of particles in well i , we evaluatethis expression at points a and b between the wells, wecan write a = q iz − (cid:96)/ and b = q iz + (cid:96)/ with (cid:96) thedistance between two wells. As in the case for the four-mode model, the main contribution of the sum results forthe summand with i = k = l . We get approximately n i = (cid:114) π (cid:16)(cid:113) A iz,R (cid:96) / (cid:17)(cid:113) A ix,R A iy,R A iz,R (cid:12)(cid:12) d i (cid:12)(cid:12) . (B9)For typical values of A iz,R (cid:96) the error function is close tounity. Thus, it is reasonable to compare the particle num-bers obtained from the four-mode model [Eq. (B3)] withthose from the extended variational ansatz [Eq. (B9)].Since the particle currents are given by the derivatives ofthe particle numbers, this is also the case for the particlecurrents.2 [1] C. E. Rüter, K. G. Makris, R. El-Ganainy, D. N.Christodoulides, M. Segev, and D. Kip, Nat. Phys. ,192 (2010).[2] Y. D. Chong, L. Ge, and A. D. Stone, Phys. Rev. Lett. , 093902 (2011).[3] J. Schindler, A. Li, M. C. Zheng, F. M. Ellis, and T. Kot-tos, Phys. Rev. A , 040101 (2011).[4] S. Bittner, B. Dietz, U. Günther, H. L. Harney, M. Miski-Oglu, A. Richter, and F. Schäfer, Phys. Rev. Lett. ,024101 (2012).[5] H. Ramezani, T. Kottos, R. El-Ganainy, and D. N.Christodoulides, Phys. Rev. A , 043803 (2010).[6] C. M. Bender and S. Boettcher, Phys. Rev. Lett. , 5243(1998).[7] S. Klaiman, U. Günther, and N. Moiseyev, Phys. Rev.Lett. , 080402 (2008).[8] D. Dast, D. Haag, H. Cartarius, G. Wunner, R. Eichler,and J. Main, Fortschritte der Physik , 124 (2013).[9] Y. Shin, G. B. Jo, M. Saba, T. A. Pasquini, W. Ketterle,and D. E. Pritchard, Phys. Rev. Lett. , 170402 (2005).[10] H. Cartarius and G. Wunner, Phys. Rev. A , 013612(2012).[11] E. M. Graefe, H. J. Korsch, and A. E. Niederle, Phys. Rev. Lett. , 150408 (2008).[12] D. Haag, D. Dast, A. Löhle, H. Cartarius, J. Main, andG. Wunner, Phys. Rev. A , 023601 (2014).[13] A. Zenesini, C. Sias, H. Lignier, Y. Singh, D. Ciampini,O. Morsch, R. Mannella, E. Arimondo, A. Tomadin, andS. Wimberger, New J. Phys. , 053038 (2008).[14] C. Elsen, K. Rapedius, D. Witthaut, and H. J. Korsch,J. Phys. B , 225301 (2011).[15] M. Kreibich, J. Main, H. Cartarius, and G. Wunner,Phys. Rev. A , 051601(R) (2013).[16] K. Henderson, C. Ryu, C. MacCormick, and M. G.Boshier, New J. Phys. , 043030 (2009).[17] A. Trombettoni and A. Smerzi, Phys. Rev. Lett. , 2353(2001).[18] A. Smerzi and A. Trombettoni, Phys. Rev. A , 023613(2003).[19] I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. , 885 (2008).[20] P. Löwdin, J. Chem. Phys. , 365 (1950).[21] E.-M. Graefe, J. Phys. A , 444015 (2012).[22] A. D. McLachlan, Mol. Phys. , 39 (1964).[23] R. Eichler, D. Zajec, P. Köberle, J. Main, and G. Wun-ner, Phys. Rev. A86