A Bose-Einstein Condensate with $\mathcal{PT}$-Symmetric Double-Delta Function Loss and Gain in a Harmonic Trap: A Test of Rigorous Estimates

We consider the linear and nonlinear Schr\"odinger equation for a Bose-Einstein condensate in a harmonic trap with $\cal {PT}$-symmetric double-delta function loss and gain terms. We verify that the conditions for the applicability of a recent proposition by Mityagin and Siegl on singular perturbations of harmonic oscillator type self-adjoint operators are fulfilled. In both the linear and nonlinear case we calculate numerically the shifts of the unperturbed levels with quantum numbers $n$ of up to 89 in dependence on the strength of the non-Hermiticity and compare with rigorous estimates derived by those authors. We confirm that the predicted $1/n^{1/2}$ estimate provides a valid upper bound on the the shrink rate of the numerical eigenvalues. Moreover, we find that a more recent estimate of $\log(n)/n^{3/2}$ is in excellent agreement with the numerical results. With nonlinearity the shrink rates are found to be smaller than without nonlinearity, and the rigorous estimates, derived only for the linear case, are no longer applicable.


Introduction
Bose-Einstein condensates with PT -symmetric loss and gain have been proposed [1] as a first experimental realisation of PT symmetry in a real quantum system. The idea is to confine the condensate in a double well potential, and to create a PT -symmetric situation by coherently injecting atoms into one well and removing them from the other. A particular difficulty in the theoretical treatment is that, on account of the s-wave scattering of the atoms, the Gross-Pitaevskii equation, which describes the condensates, contains a term g|ψ| 2 , and thus is nonlinear in the wanted wave function. In a series of papers both for realistic set-ups as well as for delta-function models of the double wells we have shown [2][3][4][5][6][7] that the nonlinearity introduces new features in the evolution of the eigenvalue spectrum as the non-Hermiticity is increased but yet PT symmetry of the wave function is preserved if both nonlinearity and non-Hermiticity are not too strong.
Bose-Einstein condensates are usually trapped in harmonic potentials produced by counterpropagating laser beams. Therefore condensates with an additional PT -symmetric double well potential can be regarded as a "perturbation" of the harmonic oscillator. Aducci and Mityagin [8] and Siegl and Mityagin [9] have recently analysed perturbations of harmoniclike operators from a mathematical point of view. The second paper in particular also allows for singular perturbations, such as delta-functions.
These authors have proved that the eigenvalues of the perturbed operator become eventually simple and the root system forms a Riesz basis. Their results are valid if the following criterion is fulfilled: whereB is the perturbation operator, M is a constant that depends on the type of the perturbation, and the ψ m are the harmonic oscillator eigenfunctions. A further prediction is that the shrink rate of the disks into which the eigenvalues can be shifted is proportional to 1/n 2α , with α the exponent appearing in (1).
It is the purpose of this paper to test the estimates numerically for the example of a Bose-Einstein condensate in a double well potential confined by a harmonic trap. To this end we consider the model of two PT -symmetric delta-function wells, since in this case simple analytical estimates can be obtained, whereas in the case of the realistic double well discussed in Refs. [4,5] complicated estimates in terms of hypergeometric functions result.
It must be mentioned that of course because of their simplicity delta functions have widely been used in the literature in the context of PT symmetry. Spectral properties of scattering and bound states in PTsymmetric double-and multiple-delta function potentials have been investigated e.g. in Refs [2,3,[11][12][13][14][15][16][17][18]. In all these papers no external potential was present, in additon to the delta potentials.
A paper in which PT -symmetric point interactions were studied embedded in an external potential is that by Jakubský and Znojil [19] who positioned the delta functions in an infinitely high square well and analysed the spectrum in dependence on the position of the delta functions within the well. Their work was extended by Krejčiřík and Siegl [20,21] who replaced the delta functions by PT -symmetric Robin boundary conditions at the edges of the square well. We note that the spectrum of a harmonic oscillator perturbed by two identical real-valued point interactions has been analysed by Fassari and Rinaldi [22]. But to our knowledge the situation of two PT -symmetric delta functions in an external harmonic potential has not yet been investigated. While in the square well the delta functions can be placed only within the well, the harmonic oscillator potential has the advantage that the delta functions in principle can be shifted to any position on the real axis.

Bose-Einstein condensate in a
PT -symmetric harmonic trap At low temperatures and densities Bose-Einstein condensates are well described by the Gross-Pitaevskii equation [23,24] − Here ψ denotes the condensate wave function, the eigenvalue µ is the chemical potential, and V (x) is the trapping potential to confine the condensate. The nonlinear term in (2) arises from the s-wave scattering interaction of the atoms; g is a measure for the strength of this interaction. We consider a harmonic trapping potential and model a PT -symmetric double well with equilibrated loss and gain by imaginary delta functions. Thus the Hamiltonian we consider is given bŷ (3) In (3) ±b denotes the position of the imaginary deltas and γ the strength of the non-Hermiticity We will later consider the effects of the nonlinearity on the spectrum, but for the time being assume that the nonlinearity is negligible in order to be in a position to compare with the predictions of Mityagin and Siegl [9].
The eigenvalues of the unperturbed spectrum are given by µ n = 2n + 1 (n = 0, 1, 2, . . . ). Figure 1 shows the unperturbed spectrum together with the wave functions of the lowest five states. The dashed vertical lines designate different positions at which the delta functions are placed. From the figure it is already obvious that only such states will be significantly affected by the perturbation which are within the classically allowed region at the positions the delta functions. By contrast, states for which the delta functions lie in the classically forbidden (exponentially decaying) regime will not be affected. This means that as the delta functions are shifted further and further out an increasing number of low-lying eigenvalues will not be changed by the perturbation. It is easy to show that the Mityagin-Siegl criterion (1) is fulfilled for the imaginary delta function perturbation in (3) since (4) The oscillator eigenfunctions have either even or odd parity, therefore | ψ m |B|ψ n | = 0 if m and n are both even or odd, while for m even and n odd, and vice versa, we have In the last line we have exploited an inequality given by Mityagin and Siegl [9] which is valid for 2 The prediction then is that eigenvalues can only move to a distance M , where M is a constant, uniform for all eigenvalues, and in particular, that there exists an n 0 such that for n ≥ n 0 all eigenvalues stay in disjoint disks with shrinking radii, with the shrink rate being bounded from above by Cn −1/2 . In fact, one of the authors of Ref. [9] has pointed out [10] that for the case of two delta potentials the shrink rate can be estimated even more precisely to behave as log(n)/n 3/2 . Note that no such statements can be made from their theorems for the case that the nonlinearity is also included as a perturbation.

Eigenvalue spectra
The (real and complex) eigenvalues µ of the Hamiltonian (3) and its eigenstates ψ are obtained, both in the linear and the nonlinear case, by integrating the wave functions outward from x = 0, and varying the initial values of Re ψ(0), ψ (0) ∈ C, and µ ∈ C to find square-integrable normalised solutions (the arbitrary global phase is exploited by choosing Im ψ(0) = 0). Figure 2 shows for the lowest 30 levels the real and imaginary parts of the eigenvalues of (3) (for g = 0) as functions of the strength γ of the non-Hermiticity, for a position of the delta functions at b = 0.2, close to the centre of the oscillator. One recognises that successive pairs of eigenvalues coalesce at branch points, from whereon they turn into complex conjugate pairs. The branch points are shifted to larger values of γ as one goes up in the spectrum. What is surprising is that both the real and imaginary part of the complex eigenvalues emerging from the branch point of the eighth and ninth excited state experience huge shifts, but eventually saturate, like the other levels. For the latter, the real parts remain approximately constant beyond the branch points, and the imaginary parts quickly tend to zero. In Figure 3 we move the delta functions further away from the centre of the harmonic oscillator to the classical turning points b = 1 and √ 7 of the unperturbed levels with n = 0 and n = 3, respectively. At the turning points the unperturbed wave functions enter into the classically forbidden, i.e. exponentially decreasing region. It is therefore no surprise that for b = 1 the ground state no longer "feels" the delta functions, and no longer unites with the first excited state at a branch point. Rather its eigenvalue remains real for any strength of the non-Hermiticity. For b = √ 7 it is the states with n = 0, 1, 2, 3 which exhibit this property.
We note that similar behaviour was found by Jakubský and Znojil [19] in their PT -symmetric square well model. In their terminology, energy levels which coalesce at a branch point and turn complex are called "fragile", while energy levels whose eigenvalues remain real for any strength of the perturbation are called "robust". The latter also include states where the delta functions happen to be at or close to a node of the wave function, and therefore remain unaffected by the perturbation. Examples of this can be seen in Figure  3.
In Figure 2 the ground state coalesces with the first excited state, while in Figure 3, at b = 1 (top panel), it has become a single real level for any value of γ, and the first excited state coalesces with the second excited one. The question arises how the transition between the different coalescence behaviour occurs. This is illustrated in Figure 4 where the real and imaginary parts of the eigenvalues emerging from the ground state and the first two excited levels are shown as functions of γ, for three positions of the delta functions around b ≈ 0.9. It is evident that at b = 0.897 the ground and first excited state still coalesce, giving rise to a pair of complex conjugate eigenvalues which "collides" with the real eigenvalue of the second excited level. The latter then turns into another pair of complex conjugate eigenvalues. The behaviour is similar for b = 0.915 but here the pair of complex eigenvalues resulting from the merger of the ground and first excited state disappears by splitting into two real eigenvalues the lower of which remains real for any γ, while the higher after a small interval of γ coalesces with the real eigenvalue of the second excited  state and a new pair of complex eigenvalues is born. Finally, at b = 0.925 the transition has occured, the ground state has become a single real level, and the two first excited states come together. We note again that similar behaviour was found by Krejčiřík and Siegl [20] in their studies of eigenvalues in a square well with Robin boundary conditions (cf. Figure 6 in [20]). We now proceed to results for nonvanishing nonlinearity. Figure 5 shows the spectrum for a value of g = 5 with the delta functions placed at b = 1. The overall behaviour is similar (many states coalesce at branch points), but there are significant differences. Like in our previous studies of PT -symmetric Hamiltonians with nonlinearity [2][3][4][5][6]25], pairs of complex conjugate eigenvalues (imaginary parts not shown) appear before the branch points are reached. Again there exist "robust" levels whose eigenvalues remain real for any value of γ.

Eigenvalue shifts
For vanishing nonlinearity, the upper part of Figure 6 shows, for γ = 1 and different positions of the delta functions, the shifts of the eigenvalues as a function of the quantum number n in comparison with the predicted shrink rates of n −1/2 of [9] and the improved estimate of log(n)/n 3/2 [10]. It can be recognized that the n −1/2 dependence is indeed an upper bound on the shrink rate, and, moreover, that the improved  estimate is in excellent agreement with the numerical data, irrespective of the position of the deltas! The picture changes when the nonlinearity is switched on. The bottom part of Figure 6 shows eigenvalue shifts for g = 2, and different positions of the imaginary delta potentials. From the comparison of the scales of the vertical axes in Figure 6 it can be seen that with nonlinearity the shifts even for the highest states are still bigger than 0.1, while they have dropped already well below 10 −2 in the linear case. Furthermore, the shrink rate is found to be slower (approximately proportional to n −0.37 ) than predicted by both rigorous mathematical estimates for the linear case. We therefore find significant differences in the shrink rates with and without nonlinearity.
For small γ, all eigenvalues are still real. The question therefore suggests itself how the situation changes when eigenvalues are involved that are shifted into the complex plane, which happens for increasing γ. Figure 7 shows for the eigenvalue spectra with b = 0.2, vanishing nonlinearity and different values of γ the distances |∆µ| of the real or complex eigenvalues from their original values (cf. also Figure 2). The outlier at n = 8 and 9 is caused by the giant change of the real and imaginary parts of the complex eigenvalues beyond the branch point of the corresponding states observed already in Figure 2. A similar outlier occurs around n = 75, and a monotonous decrease of the eigenvalue shifts in our calculations only sets in for n ≈ 80.
Finally in Figure 8 we investigate the effect of the nonlinearity on the eigenvalue shifts for γ = 1.5 and b = 0.5 and b = 1. We find that the eigenvalue shifts oscillate around straight lines with slopes -0.37, irrespective of the strength of the nonlinearity. The amplitude of the oscillations, however, and their numbers depend on the position of the delta functions. Again the slope of the estimate proportional to n −1/2 , also shown in the Figure, is steeper than the actual slopes found in the numerical results.

Conclusions
We have carried out a numerical analysis of a PTsymmetric double delta perturbation of the harmonic oscillator. We have also considered the case were in addition a Gross-Pitaevskii nonlinearity proportional to |ψ| 2 is present. With the latter, the system can be considered as a model of a Bose-Einstein condensate in a double well with loss and gain of atoms.
We have checked that the Mityagin-Siegl criterion for the perturbed eigenfunctions to form a Riesz basis is fulfilled, and compared rigorous mathematical estimates for the shrink rate of the eigenvalue shifts in dependence on the harmonic oscillator quantum number n for various strengths of the non-Hermiticity and the nonlinearity. We have verified that in the linear case the mathematical prediction for the shrink rates proportional to 1/n 1/2 is a valid estimate, and that the improved estimate proportional to log(n)/n 3/2 is in excellent agreement with the behaviour of the shrink rates found in the numerical results. By contrast, with nonlinearity we find slopes of approximately −0.37, less steep than both mathematical estimates. Evidently, to derive estimates also for the nonlinear case remains a mathematical challenge. A peculiarity found is the occurrence of outliers in the eigenvalue spectra which appear beyond branch points with unusually large real and imaginary parts of their complex conjugate eigenvalues. They are the reason why for growing strength of the non-Hermiticity the asymptotic shrink rate behaviour is attained only for high values of n. Certainly the nature of these outliers and their mathematical importance should be clarified in future studies.