Stability of Bose-Einstein condensates in a $\mathcal{PT}$ symmetric double-$\delta$ potential close to branch points

A Bose-Einstein condensate trapped in a double-well potential, where atoms are incoupled to one side and extracted from the other, can in the mean-field limit be described by the nonlinear Gross-Pitaevskii equation (GPE) with a $\mathcal{PT}$ symmetric external potential. If the strength of the in- and outcoupling is increased two $\mathcal{PT}$ broken states bifurcate from the $\mathcal{PT}$ symmetric ground state. At this bifurcation point a stability change of the ground state is expected. However, it is observed that this stability change does not occur exactly at the bifurcation but at a slightly different strength of the in-/outcoupling effect. We investigate a Bose-Einstein condensate in a $\mathcal{PT}$ symmetric double-$\delta$ potential and calculate the stationary states. The ground state's stability is analysed by means of the Bogoliubov-de Gennes equations and it is shown that the difference in the strength of the in-/outcoupling between the bifurcation and the stability change can be completely explained by the norm-dependency of the nonlinear term in the Gross-Pitaevskii equation.


Introduction
In quantum mechanics non-Hermitian Hamiltonians with imaginary potentials have become an important tool to describe systems with loss or gain effects [1]. Non-Hermitian PT symmetric Hamiltonians, i.e. Hamiltonians commuting with the combined action of the parity (P: x → −x, p → −p) and time reversal (T : x → x, p → −p, i → −i) operators, possess the interesting property that, in spite of the gain and loss, they can exhibit stationary states with real eigenvalues [2]. When the strength of the gain and loss contributions is increased typically pairs of these real eigenvalues pass through an exceptional point, i.e. a branch point at which both the eigenvalues and the wave functions are identical, and become complex and complex conjugate.
A promising candidate for the realisation of a PT symmetric quantum system are Bose-Einstein condensates. At sufficiently low temperatures and densities they can in the mean-field limit be described by the nonlinear Gross-Pitaevskii equation [3,4]. If a condensate is trapped in a double-well potential it is possible to add atoms to one side of the double well and remove atoms from the other. This leads to a gain or loss to the condensate's probability amplitude. If the strength of both contributions is equal, the process can effectively be described by a complex external potential V (−x) = V * (x) rendering the Hamiltonian PT symmetric [5]. The experimental realisation of an open quantum system with PT symmetry will be an important step since the experimental verification has only been achieved in optics so far [6,7].
From the mathematical point of view a new type of bifurcation appearing in the nonlinear PT symmetric Gross-Pitaevskii equation is of special interest [8][9][10][11][12][13][14][15][16][17]. As in linear quantum systems two PT symmetric stationary states merge in an exceptional point if the strength of a parameter describing the gain and loss processes is increased. However, in contrast to its linear counterpart the Gross-Pitaevskii equation possesses no PT broken states emerging at this exceptional point. They already appear for lower strengths of the gain/loss parameter and bifurcate from one of the PT symmetric states in a pitchfork bifurcation. The new bifurcation point has been identified to be a third-order exceptional point [10]. For attractive nonlinearities one finds that the PT broken solutions bifurcate from the ground state. In this scenario the PT symmetric ground state is the only state which arXiv:1401.2354v1 [quant-ph] 10 Jan 2014 exists on both sides of the bifurcation and always possesses a real energy eigenvalue. The pitchfork bifurcation is expected to entail a change of its stability. However, it is observed that this stability change does not occur exactly at the bifurcation but at a slightly different value of the gain/loss parameter.
The discrepancy between the points of bifurcation and stability change seems to be surprising and does not appear in all similar systems. The mean-field limit of a two mode approximation with a Bose-Hubbard Hamiltonian [14][15][16][17] does not show this effect. The model has, however, two crucial differences to the treatment of Bose-Einstein condensates with the full Gross-Pitaevskii equation in Ref. [11]. The latter system contains a harmonic trap in which infinitely many stationary states can be found, whereas the nonlinear two-mode system exhibits only four states, viz. the two PT symmetric and the two PT broken states mentioned above. Furthermore the nonlinearity derived in [14,15] is slightly different from that of the Gross-Pitaevskii equation. It has the form ∝ |ψ| 2 /||ψ|| 2 , and hence does not depend on the norm of the wave function. Thus, there might be two reasons for the appearance of the discrepancy. It could have its origin in the existence of higher modes influencing the ground state's stability or in the norm-dependency of the Gross-Pitaevskii nonlinearity.
It is the purpose of this article to clarify this question. To do so we study a Bose-Einstein condensate in an idealised double-δ trap [8][9][10], a system of which already its linear counterpart helped to understand basic properties of PT symmetric structures [29][30][31][32][33]. This system is described by the Gross-Pitaevskii equation, i.e. the contact interaction has the norm-dependent form ∝ |ψ(x, t)| 2 . However, it exhibits only four stationary states of which two are PT symmetric and two are PT broken as in the two-mode model [14][15][16][17]. Additionally, in a numerical study the structure of the nonlinearity can easily be changed such that the system's mathematical properties can be brought in agreement with the mean-field limit of the Bose-Hubbard dimer.
The article is organised as follows. We will introduce and solve the Gross-Pitaevskii equation of a Bose-Einstein condensate in a double-δ trap for an attractive atom-atom interaction in Section 2. Some properties of the stationary solutions which are important for the following discussions are recapitulated. Then we will investigate the ground state's stability in the vicinity of the bifurcation in Section 3. The Bogoliubov-de Gennes equations are solved for both types of the nonlinearity, and the origin of the discrepancy between the bifurcation and the stability change is discussed. Conclusions are drawn in Section 4.

Bose-Einstein condensates in the PT symmetric double-δ trap
We assume the condensate to be trapped in an idealised trap of two delta functions, i.e. the potential has the shape [8] where the units are chosen such that the real part due to the action of the δ functions has the value -1. It describes two symmetric infinitely thin potential wells at positions ±b. In the left well we describe an outflux of atoms by a negative imaginary contribution −iγ, and on the right side an influx of particles is described by a positive imaginary part +iγ of the same strength. This leads to the time-independent Gross-Pitaevskii equation in dimensionless form [8], with the energy eigenvalue or chemical potential µ = −κ 2 . The parameter g is determined by the s-wave scattering length, which effectively describes the van der Waals interaction for low temperatures and densities. Physically it can be tuned via Feshbach resonances. Throughout this article we assume g to be negative, i.e. the atom-atom interaction is attractive. Solutions to the Gross-Pitaevskii equation are found with a numerical exact integration. The wave function is integrated outward from x = 0 in positive and negative direction. To do so the initial values ψ(0), ψ (0), and κ have to be chosen. The arbitrary global phase is exploited such that ψ(0) is chosen to be real. Together with ψ (0) ∈ C and κ ∈ C we have five parameters which have to be set such that physically relevant wave functions are obtained. These have to be square-integrable and normalised in the nonlinear system. The five conditions ψ(∞) → 0, ψ(−∞) → 0, and ||ψ|| = 1 ensure that these conditions are fulfilled, and, together with the five initial parameters, define a five-dimensional root search problem, which is solved numerically. Figure 1 shows a typical example for all stationary states found in the case g = −1 (red solid lines). Two states with purely real eigenvalues vanish for increasing γ in an exceptional point. Since κ is plotted and µ = −κ 2 is the energy eigenvalue, the upper line corresponds to the ground state. Two complex and complex conjugate eigenvalues bifurcate from this ground state in a pitchfork bifurcation at a critical value γ κ ≈ 0.3071. The same plots for g = −0.5 (blue dotted lines) and g = 0 (green dashed lines) demonstrate how the pitchfork bifurcation is introduced by the nonlinearity.   (2) for the case g = −1 (red solid lines), which is compared with its linear counterpart g = 0 (green dashed lines) and the weaker interaction g = −0.5 (blue dotted lines). In the nonlinear case we observe two complex conjugate states bifurcating from the ground state in a pitchfork bifurcation. For g = −1 this occurs at γ κ ≈ 0.3071.
This eigenvalue structure is very generic for PT symmetric systems with a quadratic term in the Hamiltonian. It is found for Bose-Einstein condensates in a true spatially extended double well in one and three dimensions [11,12] which contain the nonlinearity ∝ |ψ| 2 as well as for the mean-field limit of the Bose-Hubbard dimer [14][15][16][17] with a term ∝ |ψ| 2 /||ψ| 2 . The difference in the norm-dependency of the nonlinearity between both systems does not lead to different eigenvalues, which has already been mentioned [17]. The difference appears, however, for the dynamics, where the norm plays a crucial role in the presence of gain and loss [13].

Stability analysis of the ground state
The linear stability is analysed with the Bogoliubov-de Gennes equations. They are derived under the assumption that a stationary state ψ 0 (x, t) is perturbed by a small fluctuation θ(x, t), i.e. where With this ansatz and a linearisation in the small quantities u and v one obtains from the Gross-Pitaevskii equation the coupled system of the Bogoliubov-de Gennes differential equations, In Equation (3b) it can be seen that ω decides on the temporal evolution of the fluctuation. Real values of ω describe stable oscillations, whereas imaginary parts lead to a growth or decay of the fluctuation's amplitude. Thus, ω measures the stability of the stationary solution ψ 0 against small fluctuations.
Numerically the Bogoliubov-de Gennes equations are solved with the same method as the stationary states, i.e. the wave functions u and v are integrated outward from x = 0. It can easily be seen that the Bogoliubov-de Gennes equations (4) are invariant under the transformation u(x) → u(x)e iχ , v(x) → v(x)e iχ with a real phase χ. Similarly to the procedure for the integration of the stationary states this symmetry can be exploited. The remaining initial values with which the integration has to be started are In a nine-dimensional root search they have to be chosen such that the conditions u(±∞) → 0, v(±∞) → 0, and ∞ −∞ |u(x) + v * (x)| 2 dx = 1 are fulfilled [11].
Two further symmetries can be exploited to reduce the number of independent stability eigenvalues which have to be calculated. The replacement (u, v, ω) → (v * , u * , −ω * ) leaves the ansatz (3b) invariant. Thus, if ω is a stability eigenvalue also −ω * is a valid solution. Furthermore for every eigenvalue ω there is also one solution with the eigenvalue ω * if the stationary state ψ 0 is PT symmetric. This can be verified by applying the PT operator to the Bogoliubov-de Gennes equations (4). Due to these two symmetries it is sufficient to search for stability eigenvalues with Re(ω) > 0 and Im(ω) > 0.

Stability in the vicinity of the bifurcation
The relevant question which has to be answered by our calculation is whether or not the discrepancy between the γ values of the pitchfork bifurcation and the stability change appears for the double-δ potential. Thus we calculated the stability eigenvalue with Re(ω) > 0, Im(ω) > 0 for a range of γ around the bifurcation, which is shown in Figure 2 for g = −1. For increasing γ the eigenvalue ω switches from real to imaginary at γ ω ≈ 0.3138 marking the stability change. The pitchfork bifurcation is visible in the real parts of κ of all stationary states of the system. It is marked by the black dashed-dotted line. Obviously there is a discrepancy between γ κ and γ ω . The difference is ∆γ ≈ −0.0067.
The system does not possess any further stationary states besides those shown in Figure 1. Three of these four states are participating in the pitchfork bifurcation. Only the excited PT symmetric solution would be able to influence the dynamics of the ground state at a value γ = γ ω . However, it stays real for all parameters γ shown in Figure 2 and cannot cause any qualitative different behaviour of the ground state's dynamics. Thus, an influence of further states can be ruled out to be the reason for the discrepancy in the double-well system of Refs. [11,12].
The remaining difference between the Gross-Pitaevskii equation (2) and the two mode system of Refs. [14][15][16][17] is the norm-dependency of the nonlinearity. Indeed, an influence of the norm is already present if the study done in Figure 2 is repeated for different values of the nonlinearity parameter g. Figure  3 shows ∆γ as a function of g. A strong dependency is visible. Even the sign changes. For g 0.4 the ground state becomes unstable at γ ω < γ κ . For g → 0 the discrepancy vanishes as expected.

Norm-independent variant of the Gross-Pitaevskii equation
An even clearer identification of ∆γ with the normdependency of the nonlinearity in the Gross-Pitaevskii equation (2) can be given with a small modification. The replacement makes the Gross-Pitaevskii equation (2) normindependent. Note that this is exactly the form of the mean-field limit of Refs. [14][15][16][17]. Since the stationary states are normalised to 1 they are not influenced by the replacement. However, it influences the dynamics and also the linear stability. The Bogoliubov-de Gennes equations have to be adapted. Assuming again a small perturbation of the form (3a) and (3b) a linearisation in u and v leads us to For a numerical solution of the modified Bogoliubovde Gennes equations the value of S is included in the root search. i.e. for the integration of the Equations (7a) and (7c) a value for S is guessed and subsequently compared with the result of Equation (7c). Since S is in general a complex value this increases the dimension of the root search to 11.
An example for a typical result is shown in Figure  4. The discrepancy between the γ values at which the pitchfork bifurcation and the stability change occur vanishes. Both values agree as it is expected for a pitchfork bifurcation. This is true for all values of g. Thus, we conclude that the discrepancy appearing in

Discussion
The reason why the stability change does not occur exactly at the bifurcation can also be understood intuitively. Figure 1 indicates that the value γ κ of the pitchfork bifurcation is not equal for all values of g. For g = 0 there is no pitchfork bifurcation. The two real eigenvalues κ vanish in a tangent bifurcation and two complex eigenvalues emerge. Only for nonvanishing g the pitchfork bifurcation exists and moves to smaller values of γ for increasing |g|.
If now the norm of a wave function changes, this can also be understood as a variation of g, i.e. a wave function with a norm N = ||ψ|| 2 has the same effect as a wave function with the norm 1 and a modified valuẽ g = N g. The values γ κ of the pitchfork bifurcation is obviously different for g andg. With this relation in mind it is not surprising that a fluctuation changing the norm of the wave function may cause a qualitative change of the condensate's stability properties in the vicinity of γ κ .

Conclusions
We investigated the origin of the discrepancy between the value of the gain/loss parameters, at which the ground state of a Bose-Einstein condensate in a doublewell trap passes through a pitchfork bifurcation and at which its stability changes. In a naive expectation these γ values should be identical. However, it is found that this is not exactly fulfilled. Since this discrepancy does not occur for a similar system, viz. the mean-field limit of a Bose-Hubbard dimer [14][15][16][17], we investigated the differences between the equations describing both systems. It was found that the normdependency of the nonlinearity in the Hamiltonian ∝ |ψ| 2 is unambiguously the origin of the discrepancy. It can be completely removed with the replacement |ψ| 2 → |ψ| 2 /||ψ|| 2 . An intuitive explanation can be given. Fluctuations which change the norm of a stationary state are able to shift the position γ κ of the bifurcation.
Since the dynamical properties of the wave functions are crucial for an experimental observability of a PT symmetric Bose-Einstein condensate it is important to know about all processes introducing possible instabilities. As has been shown in this article the stability relations are nontrivial close to branch points in condensate setups with gain and loss. The Gross-Pitaevskii has a norm-dependent nonlinearity, and therefore it should be clarified in future work, which type of fluctuation influences the wave function's norm such that additional instabilities appear. In particular, it would be interesting to see, how the amplitude of a fluctuation is related to the size of the difference ∆γ. Also a deeper understanding, how this effect influences realistic setups generating the PT symmetric external potential [34], would be of high value.