Acta Polytechnica Vol. 51 No. 4/2011 Exceptional Points for Nonlinear Schrödinger Equations Describing Bose-Einstein Condensates of Ultracold Atomic Gases G. Wunner, H. Cartarius, P. Köberle, J. Main, S. Rau Abstract The coalescence of two eigenfunctions with the same energy eigenvalue is not possible in Hermitian Hamiltonians. It is, however, a phenomenon well known from non-hermitian quantummechanics. It can appear, e.g., for resonances in open systems, with complex energy eigenvalues. If two eigenvalues of a quantum mechanical system which depends on two or more parameters pass through such a branch point singularity at a critical set of parameters, the point in the parameter space is called an exceptional point. We will demonstrate that exceptional points occur not only for non-hermitean Hamiltonians but also in the nonlinear Schrödinger equations which describe Bose-Einstein condensates, i.e., the Gross- Pitaevskii equation for condensates with a short-range contact interaction, and with additional long-range interactions. Typically, in these condensates the exceptional points are also found to be bifurcation points in parameter space. For condensates with a gravity-like interaction between the atoms, these findings can be confirmed in an analytical way. 1 Introduction In 1924, Satyendra Nath Bose and Albert Einstein predicted that when the thermal de Broglie wave- length becomes of the order of the interparticle dis- tance, bosons begin to “condense” into their ground state. All bosons have the same energy and quan- tum characteristics, similar to the way all photons in a laser share the same quantum state. Superfluidity of 4He atoms below a temperature of approximately T = 2.17 K was the first experimental realisation in 1937 of such a macroscopic quantum state [1]. It took, however, until 1995 for Bose-Einstein condensa- tion to be observed in dilute atomic vapours of alkali atoms, namely, by Wieman et al. [2] in 87Rb and by Ketterle et al. [3] in 23Na (in both cases the number of nucleons and electrons add up to an even number, to form a boson). In these experiments, advanced optical techniques (capture in a magneto-optical trap and application of laser cooling and evaporative cool- ing) had to be employed to cool the atoms to tem- peratures very near to absolute zero (170 nK in the rubidium experiment and 2 μK in the sodium ex- periment). Meanwhile Bose-Einstein condensates of many more alkaline and alkaline earth atoms have been produced (for a review see, e.g., ref. [4]). Al- kalis, however, are not ideal Bose gases: at small distances, the atoms still interact via the short-range van der Waals force, whose action can be simulated by a short-range isotropic s-wave scattering contact potential Vs(r,r ′) = 4πah̄2 m δ(r − r′), where a is the s-wave scattering length. The latter can be tuned, using Feshbach resonances of hyperfine levels and changing the applied magnetic field, from positive values (repulsive interaction) to negative values (at- tractive interaction), until the interaction becomes so attractive that the condensates collapse. In this paper we will be concerned with Bose- Einstein condensates where an additional long-range atom-atom interaction is active. An example is Bose- Einstein condensation of 52Cr atoms (Griesmeier et al. [5], Koch et al. [6]), with a large magnetic mo- ment of μ = 6μB, so that the atoms interact via the dipole-dipole interaction over large distances. In these condensates the possibility of tuning the rela- tive strengths of the short-range interaction and the long-range interaction has opened the way to a vari- ety of new phenomena (see, e.g., the review article by Lahaye et al. [7]). Furthermore, the stability of these condensates depends strongly on the trap geometry. The atoms are aligned by an external magnetic field in the z direction. If the confinement in the z direc- tion is stronger than perpendicular to it (if the aspect ratio λ of the trapping frequencies ωz/ω� is greater than one), the condensate assumes an oblate shape and the dipole-dipole interaction acts predominantly repulsive, while for λ < 1 the shape is prolate, the atoms are arranged in the favoured head-to-tail con- figuration, and the interaction acts attractive. As an alternative system, O’Dell et al. [8] have proposed Bose-Einstein condensation of neutral atoms with an electromagnetically induced attrac- tive long-range 1/r interaction. In their proposal, 6 “triads” of intense off-resonant laser beams aver- age out 1/r3 interactions in the near-zone limit of the retarded dipole-dipole interaction of the atoms produced by the irradiation, while the weaker 1/r 95 Acta Polytechnica Vol. 51 No. 4/2011 interaction is retained. The novel feature of these condensates is that the atoms can be self-trapped, without an external trap. On the theoretical side, the advantage is that for self-trapping analytical vari- ational calculations are possible, which can serve as a guide for investigations of more complex situations and condensates with other interatomic interactions. A theoretical description of condensates of weakly interacting bosons at ultracold temperatures can be performed within a mean-field picture by the Gross- Pitaevskii equation. This equation can be consid- ered as the Hartree equation for the single particle orbital ψ(r) which all atoms occupy. It is a nonlin- ear Schrödinger equation, and therefore phenomena not present in linear quantum mechnics appear. For example, the superposition principle is no longer ap- plicable to the equation. For Bose-Einstein condensates with gravity-like long-range interaction Vu(r,r ′) = −u/|r−r ′| (“monopolar condensates”, u describes the strength of the interaction) the Gross-Pitaevskii equation reads{ −Δ + γ2r2 + N 8π a au |ψ(r)|2 − (1) 2N ∫ |ψ(r′)|2 |r−r′| d3r′ } ψ(r) = ih̄ ∂ψ(r) ∂t . Here γ = h̄ω0/Eu is the trapping frequency measured in the time scale given by the “Rydberg” energy Eu associated with the strength of the 1/r interaction, N is the particle number, and a/au the scattering length in units of the “Bohr” radius au. It can be shown [9] that the equation effectively depends only on two physical parameters, viz. the particle number scaled quantities N2a/au and γ/N 2. For condensates with dipole-dipole interactions (“dipolar condensates”) in which the atoms are spin- polarized in the z direction, the Gross-Pitaevskii equation reads { − Δ + γ2ρ ρ 2 + γ2z z 2 + N 8π a ad |ψ(r)|2 + N ∫ |ψ(r ′)|2 (1 − 3 cos2 ϑ′) |r−r ′|3 d3r ′ } ψ(r) = (2) ih̄ ∂ψ(r) ∂t , where we have used as units of length, energy, and frequency the quantities ad = μ0μ 2m/(2πh̄2) (μ is the magnetic moment), Ed = h̄ 2/(2ma2d), and ωd = Ed/h̄, respectively. Again it can be shown [10] that the equation only depends on the three scaled param- eters N2γ̄, λ, a/ad (or alternatively N 2γ̄ = γ2/3ρ γ 1/3 z , λ = γz /γρ). Equations (1) and (2) are the starting point for the investigations below. 2 Monopolar condensates 2.1 Variational results To find stationary solutions of the Gross-Pitaevskii equation (1) we can enter into the corresponding mean-field energy functional E[ψ] = N ∫ d3r ψ∗(r) ( −Δr + γ2r2 + 4πN a au |ψ(r)|2 − N ∫ d3r′ |ψ(r′)|2 |r−r′| ) ψ(r) with the Gaussian variational ansatz ψ(r) = A exp ( −k2r2/2 ) , with A = ( k/ √ π )3/2 . All integrals, including the gravity-like interaction in- tegral, can be evaluated analytically. Requiring (3) to become stationary with respect to variation of the width leads to a quintic equation for k, which, how- ever, reduces to a simple quadratic equation in the special case γ = 0, i.e., for a self-trapping condensate. The two solutions read [11, 13] k± = 1 2 √ π 2 1 N a au ( ± √ 1 + 8 3π N2 a au − 1 ) , (3) where the plus sign corresponds to a local minimum (stable ground state of the condensate) and the minus sign to a local maximum (collectively excited state of the condensate) of the energy. From (3) it can be seen that these solutions only exist for N2a/au ≥ −3π/8. Accurate numerical stationary solutions of (1) can be determined by direct outward integration start- ing with suitable initial conditions at r = 0 [9]. Fi- gure 1 compares the variational and numerically ex- act results for the mean-field energy and the chemi- cal potential. It can be seen that the simple Gaus- sian ansatz qualitatively well describes the overall be- haviour: both in the variational calculations and in the numerical calculations two solutions are born in a tangent bifurcation at a critical negative value of the scattering strength, with the value being only slightly underestimated by the Gaussian ansatz (N2a/au = −3π/8 = −1.178 0 variationally, and −1.025 1 nu- merically). We note that similar tangent bifurcation behaviour is also present for non-vanishing trapping potentials [9]. At the bifurcation point, the energies and the wave functions ψk+ and ψk− are identical, which is characteristic of exceptional points known to appear in non-hermitian quantum mechanics. To confirm that the bifurcation point is an exceptional point, we have to perform a circle around the bifurcation point. If the two eigenvalues permute after travers- ing the full circle, we have an exceptional point. This requires a two-dimensional parameter space and thus 96 Acta Polytechnica Vol. 51 No. 4/2011 Fig. 1: Comparison of variational and numerically accurate results for themean-field energy (left) and chemical potential (right) of a self-trapped monopolar condensate as a function of the particle number scaled scattering length Fig. 2: Paths of the real and imaginary parts of the mean-field energy and the chemical potential in the variational calculation with a Gaussian-type orbital for self-trapped monopolar condensates if a circle of radius 10−8 is traversed around the bifurcation point N2a/au = −3π/8 an extension of the scattering length to complex val- ues [11, 13]. We demonstrate this first for the analytical so- lution obtained with the Gaussian ansatz for self- trapping condensates. Figure 2 shows the resulting paths of the values of the mean-field energy and the chemical potential in the complex plane if a small circle reiϕ, ϕ = 0 . . . 2π, with radius r = 10−8 is per- formed around the critical scattering length N2a/au = −3π/8. A surprising result is that while for the chemical potential a clear permutation of the two solutions is visible, after a full circle the two mean-field energy values do not permute. This can be explained by looking at the analytical fractional power series ex- pansion of the mean-field energy around the critical value [11] Ẽ±(ϕ) = − 4 9π + 0 · √ reiϕ/2 + 32 27π2 · √ r 2 eiϕ ±( 4 9π − 32 9π2 ) · √ r 3 e(3/2)iϕ + O (√ r 4 ) . Evidently the first-order term with the phase factor eiϕ/2 vanishes, the lowest non-vanishing order has the phase factor eiϕ, which does not lead to a permuta- tion of the energies, and it is only the third-order term which can yield the permutation. By contrast, in the fractional power series expansion of the chem- ical potential μ̃±(ϕ) = − 20 9π ± 8 3π · √ reiϕ/2 −( 4 3π + 128 27π2 ) · √ r 2 eiϕ ±( 8 9π − 64 9π2 ) · √ r 3 e(3/2)iϕ + O (√ r 4 ) the first-order term does not vanish and leads to a permutation of the chemical potential values. Thus the permutation of the eigenvalues appears for a small radius, in the same way that it occurs for ex- ceptional points in simple non-hermitian linear model systems. 97 Acta Polytechnica Vol. 51 No. 4/2011 Fig. 3: Same as Figure 2, but for a larger radius of r =10−3 Fig. 4: Same as Figure 2, but for a still larger radius of r =10−1 For increasing radius, the higher-order terms in the expansions become more and more important, and the permutation of the mean-field energy values after a full circle becomes recognizable, as can be seen in Figure 3. The permutation behaviour of the chem- ical potential solutions is not altered, as is expected from the expansion. For the case of an ever increasing radius of r = 10−1, shown in Figure 4, the higher-order terms lead to a deformation of the circle of the mean-field energy (as in linear model systems). The shape of the chem- ical potential circle does not change, but the spacing between the points is no longer uniform. For the Gaussian ansatz, analytical expansions can also be performed for circles around the bifur- cation point with large radii, r � 1. The expression for the mean-field energy reads [11] Ẽ±(ϕ) = ± ( π 16 − 1 2 ) · e−iϕ/2 √ r + 1 2 · e−iϕ √ r 2 ±( 3π2 64 − 3π 8 ) · e−(3/2)iϕ √ r 3 + O ( 1 √ r 4 ) while for the chemical potential we obtain μ̃±(ϕ) = ± (π 8 − 1 ) · e−iϕ/2 √ r + ( 1 − 3π 16 ) · e−iϕ √ r 2 ±( 3π2 32 − 3π 8 ) · e−(3/2)iϕ √ r 3 + O ( 1 √ r 4 ) . From the expressions it is evident that a clear semi- circle structure should appear for large radii, as can be seen in Figure 5 for r = 108. 2.2 Numerical analysis of the bifurcation point for self-trapped condensates For accurate numerical solutions of the Gross- Pitaevskii equation a similar investigation of the ex- ceptional point behaviour can be carried out. The task again is to perform a circle in the complex scat- tering length plane around the numerically accurate bifurcation point at N2a/au = −1.025 1 and to nu- merically solve the equation for each complex scatter- ing length on the circle. This yields not only complex values for the chemical potential and the mean-field energy but also complex wave functions. 98 Acta Polytechnica Vol. 51 No. 4/2011 Fig. 5: Same as Figure 2, but for the very large radius r =108 Fig. 6: Comparison of variational andnumerically accurate results proving thebranchpoint singularity at the bifurcation point of the Gross-Pitaevskii equation for self-trapped Bose-Einstein condensates with an attractive 1/r long-range interaction. Panels (a) and (b) show the circles of radius r = 10−3 in the complex scattering length plane around the critical values, (c) and (d) the reactions of the chemical potential values as the circles are traversed, (e) and (f) those of the mean-field energies (variational results: left column, accurate numerical results: right column). (Adapted from [11]) 99 Acta Polytechnica Vol. 51 No. 4/2011 The problem is that the Gross-Pitaevskii equation contains the square modulus of the wave function ψ and is therefore a non-analytic function of ψ. We use the following procedure for complex wave functions. Any complex wave function can be written as ψ(r) = eα(r)+iβ(r) , (4) where the real functions α(r) and β(r) determine the amplitude and phase of the wave function, respec- tively. The complex conjugate and the square modulus of ψ(r) read ψ∗(r) = eα(r)−iβ(r) , |ψ(r)|2 = e2α(r) . (5) The Gross-Pitaevskii system can then be transformed into two coupled nonlinear differential equations for the real functions α(r) and β(r), however, without any complex conjugate or square modulus. These equations can now be continued analytically by al- lowing for complex valued functions α(r) and β(r). This implies that ψ∗(r) = eα(r)−iβ(r) , |ψ(r)|2 = e2α(r) . without complex conjugation of α(r) and β(r) is for- mally used for the calculation of ψ∗ and |ψ(r)|2. As a consequence, the square modulus of ψ can become complex, leading to a complex absorbing potential in the Gross-Pitaevskii equation. Figure 6 shows a comparison between the results of the variational and the numerically accurate cal- culation. The first row shows the circles of radius r = 10−3 around the respective critical values of the scattering length where the bifurcation appears. The second row compares the results for the chemical po- tential, and the third row for the mean-field energy as the circles in the complex scattering length plane are traversed. Obviously the permutation behaviour is identical in both approaches, proving that the bifur- cation point is indeed an exceptional point also in the full numerical solution of the nonlinear Schrödinger equation which describes the Bose-Einstein conden- sation of ultracold atomic gases with a long-range attractive 1/r interaction. 3 Dipolar condensates To find stationary solutions of the Gross-Pitaevskii equation (2) of dipolar gases, we could of course again extremize the mean-field energy functional for a given variational ansatz. A more sophisticated way is to exploit the time-dependent variational principle [12] ||iφ(t) − Hψ(t)||2 != min with respect to φ (6) and afterwards set ψ̇ ≡ φ. If a complex parametriza- tion of the trial wave function ψ(t) = χ(λ(t)) is as- sumed, the variation leads to equations of motion for these parameters λ(t) [13]. 〈 ∂ψ ∂λ ∣∣∣iψ̇ − Hψ〉 = 0 ↔ (7) Kλ̇ = −ih with K = 〈 ∂ψ ∂λ ∣∣∣∂ψ ∂λ 〉 ,h = 〈 ∂ψ ∂λ ∣∣∣H∣∣∣ψ〉 . 3.1 Variational results We assume an axisymmetric trapping potential, and take the time-dependent Gaussian trial wave function ψ(ρ, z, t) = ei(Aρρ 2+Az z 2+γ), (8) Aρ = Aρ(t), Az = Az(t), γ = γ(t) , where all three parameters are complex functions of time. The equations of motion that follow for the real and imaginary parts of the variational param- eters from the time dependent variational principle read Ȧrρ = −4((A r ρ) 2 − (Aiρ) 2) + fρ(A i ρ, A i z , γ i) Ȧiρ = −8A r ρA i ρ Ȧrz = −4((A r z) 2 − (Aiz ) 2) + fz(A i ρ, A i z , γ i) Ȧiz = −8A r zA i z γ̇r = −4Aiρ − 2A i z + fγ (A i ρ, A i z , γ i) γ̇i = 4Arρ + 2A r z and are solved with the initial values Arρ = 0, A i ρ > 0, Arz = 0, A i z > 0 and γi = 1 2 ln π3/2 2 √ 2Aiρ √ Aiz . (9) This results in four remaining coupled ordinary differ- ential equations for Ȧrρ, Ȧ i ρ, Ȧ r z , Ȧ i z. Variational sta- tionary solutions of the Gross-Pitevskii equation (2) are then obtained by searching for the fixed points of these differential equations. As in the monopo- lar case two stationary points are found, one stable, representing the ground state, and one unstable, rep- resenting a collectively excited state. Figure 7 shows the results for the chemical po- tential and the mean-field energy as functions of the scattering length for different aspect ratios. It is evi- dent that again the two stationary solutions are born at a critical scattering length in a tangent bifurcation, and below the critical scattering length no stationary solutions exist. Are the bifurcation points found in the varia- tional calculations with an axisymmetric Gaussian wave function also exceptional points? To check 100 Acta Polytechnica Vol. 51 No. 4/2011 Fig. 7: Chemical potential (left) and mean-field energy (right) as functions of the scattering length for an axisymmetric variational Gaussian ansatz for dipolar condensates for N2γ̄ =3.4 × 104, and different aspect ratios Fig. 8: Chemical potentials (middle) and mean-field energies (right) for the two stationary solutions that emerge in a tangent bifurcation when a circle with radius r = 10−3 around the critical scattering length acrit is traversed in the complex extended parameter plane (left) (for N2γ̄ =3.4×104, and two different aspect ratios). The permutation of the quantities after one cycle proves that the bifurcation points are exceptional points this, an extension of the scattering length to com- plex values is again necessary. In contrast to the self-trapped monopolar case analytical investigations are no longer possible, and the analysis has to be carried out numerically. Figure 8 shows the results for the behaviour of the chemical potential and the mean-field energy when the scattering lengths cor- responding to the bifurcation points in Figure 7 for the aspect ratios λ = 1 and λ = 6 are circled in the complex plane with a radius r = 10−3. It is evident that both the values of the chemical potentials and the mean-field energies of the two solutions permute as one full cycle is traversed, an unambiguous sign that the bifurcation points found in the solution of the nonlinear Schrödinger equation which describes the Bose-Einstein condensation of ultracold atomic gases with a long-range dipole-dipole interaction are also exceptional points. 3.2 Numerical results More sophisticated variational calculations for dipo- lar condensates can be performed using the method of coupled Gaussian wave packets. For the Gross- Pitaevskii equation with dipolar interaction the method has been shown [14, 15] to be a full-fledged alternative to numerical computations on a grid with the method of imaginary time evolution. In contrast to the latter, the method allows us to find both sta- ble and unstable solutions. The idea is to take as trial functions superpositions of N different Gaus- sians centered at the origin Ψ(r) = N∑ k=1 ei(a k x(t) x 2+aky(t) y 2+akz(t) z 2+γk(t)) ≡ N∑ k=1 gk(yk (t),r) , 101 Acta Polytechnica Vol. 51 No. 4/2011 Fig. 9: Mean-field energy of a dipolar condensate for (particle number scaled) trap frequencies N2γz = 25200 and N 2 γρ = 3600 as a function of the scattering length in units of ad. In the variational calculation with one Gaussian a stable ground state (gN=1) and an unstable excited state (eN=1) emerge in a tangent bifurcation. In the calculation with coupled Gaussians two unstable states emerge (labeled ucoupled), of which the lower one turns into a stable ground state (gcoupled) in the region of the small rectangle in a pitchfork bifurcation. (Adapted from [14]) where the ak(t) are complex width parameter func- tions (3N ) and the γk(t) fix the weights and phases of the individual Gaussians (N ). Equations of mo- tion of the variational parameters again follow from the time-dependent variational principle. Figure 9 shows the results for a dipolar conden- sate as a function of the scattering length for trap frequencies N2γz = 25 200 and N 2γρ = 3 600. In the calculation with five Gaussians the results from a grid calculation are well reproduced, and the bifur- cation of the mean-field energy is confirmed. Again the variational calculation with a single Gaussian underestimates the critical scattering length where the bifurcation appears. However, a more compli- cated bifurcation scenario evolves: As the scattering length is decreased the stable ground state solution already turns unstable before reaching the tangent bifurcation point, indicating a pitchfork bifurcation in the region of the small rectangle shown in the fi- gure. Neither the tangent nor the pitchfork bifurca- tion point have been conclusively analyzed so far as to their properties as exceptional points, but the analy- sis should reveal new features, in particular when the tangent bifurcation point is encircled with a radius that also encompasses the pitchfork bifurcation. 4 Summary We have demonstrated that “nonlinear versions” of exceptional points appear in bifurcating solutions of the (extended) Gross-Pitaevskii equations describing the Bose-Einstein condensation of ultracold atomic gases with attractive 1/r interaction and with dipole- dipole interaction. The identification as exceptional points required a complex extension of the scattering length. In summary, Bose-Einstein condensates near the collapse point can be regarded as realizations of physical systems close to exceptional points. References [1] Kapitza, P., Allen, J. F., Misener, A. D.: Nature 141, (1938), 74, 75. [2] Anderson, M. H., Ensher, J. R., Matthews, M. R., Wieman, C. E., Cornell, E. A.: Science 269 (1995), 198. [3] Davis, K. B., Mewes, M. O., Andrews, M. R., van Druten, N. J., Durfee, D. S., Kurn, D. M., Ketterle, W.: Phys. Rev. Lett. 75 (1995), 3 969. [4] Ueda, M.: Fundamentals and New Frontiers of Bose-Einstein Condensation, World Scientific Publishing, 2010. [5] Griesmaier, A., Werner, J., Hensler, S., Stuh- ler, J., Pfau, T.: Phys. Rev. Lett. 94 (2005), 160 401. [6] Koch, T., Lahaye, T., Metz, J., Fröhlich, B., Griesmaier, A., Pfau, T.: Nature Physics 4 (2008), 218. [7] Lahaye, T., Menotti, C., Santos, L., Lewen- stein, M., Pfau, T.: Rep. Prog. Phys. 72 (2009), 126 401. 102 Acta Polytechnica Vol. 51 No. 4/2011 [8] O’ Dell, D., Giovanazzi, S., Kurizki, G., Akulin, V. M.: Phys. Rev. Lett. 84 (2000), 5 687. [9] Papadopoulos, I., Wagner, P., Wunner, G., Main, J.: Phys. Rev. A 76 (2007), 053 604. [10] Köberle, P., Cartarius, H., Fabčič, T., Main, J., Wunner, G.: New J. Phys. 11 (2009), 023 017. [11] Cartarius, H., Main, J., Wunner, G.: Phys. Rev. A 77 (2008), 013 618. [12] McLachlan, A. D.: Mol. Phys. 8 (1964), 39. [13] Cartarius, H., Fabčič, T., Main, J., Wunner, G.: Phys. Rev. A 78 (2008), 013 615. [14] Rau, S., Main, J., Köberle, P., Wunner, G.: Phys. Rev A 81 (2010), 031 605(R). [15] Rau, S., Main, J., Köberle, P., Wunner, G.: Phys. Rev A 82 (2010), 046 201. Guenter Wunner E-mail: wunner@itp1.uni-stuttgart.de H. Cartarius P. Köberle J. Main S. Rau Institute for Theoretical Physics 1 University of Stuttgart 70550 Stuttgart, Germany 103