Supersonic optical tunnels for Bose-Einstein condensates

We propose a method for the stabilisation of a stack of parallel vortex rings in a Bose-Einstein condensate. The method makes use of a hollow laser beam containing an optical vortex. Using realistic experimental parameters we demonstrate numerically that our method can stabilise up to 9 vortex rings. Furthermore we point out that the condensate flow through the tunnel formed by the core of the optical vortex can be made supersonic by inserting a laser-generated hump potential. We show that long-living immobile condensate solitons generated in the tunnel exhibit sonic horizons. Finally, we discuss prospects of using these solitons for analogue gravity experiments.


Introduction
Ultracold quantum gases -Bose-Einstein condensates (BECs) -provide a uniquely accessible testbed for an astonishing range of effects occuring in diverse physical systems ranging from antiferromagnetic materials [1] and nonlinear photonic structures [2] to astrophysical black holes [3,4]. The recent observation of a quantum three-body bound state in an ultracold cesium gas [5] is only one example of a fundamental effect that was predicted decades ago, but had not lent itself to experimental observation in any other physical system. The detailed exploration of these distinct physical phenomena in condensates is enabled by well-developed experimental techniques that provide unprecedented control over the interatomic interactions and the potentials confining the atoms. Advanced trapping schemes have recently been demonstrated by constructing box-like traps [6] and circular waveguides [7].
In this article we study a Bose-Einstein condensate with repulsive interatomic interactions confined in a harmonic trap and additionally exposed to a laser beam containing an optical vortex [8]. By tightly focussing the laser beam blue-detuned from an atomic resonance, it is possible to create a finite length optical tunnel for the BEC, defined by the "dark" region of the beam in the vortex core. This dark region is located between the central vertical white lines in Fig. 1(a) and corresponds to the potential valley between the potential peaks in Fig. 2(a).
The main purpose of this paper is to show that the proposed trapping configuration can enable the controlled study of a complex topological structure consisting of multiple stationary vortex rings in a BEC. Unsupported ring vortices are energetically unstable against contraction [9] and have so far been only observed in non-stationary situations, for example as the decay product of dark solitons or in soliton collisions [10,11]. In contrast, our method could enable the stabilisation of a stack of up to nine parallel vortex rings under realistic experimental conditions. The creation of stationary vortex rings would also aid studies of their interactions [12] and represent a stepping-stone towards the creation of even more complex matter-wave structures such as skyrmions [13,14,15,16].
A skyrmion is a topologial soliton in a two component BEC, in which the contraction of the ring vortex in one hyperfine component is prevented by filling it with a second hyperfine state carrying a line vortex. This also represents a tool to stabilise multiple ring vortices [26]. The optical tunnel system presented here results from a skyrmion if one replaces the hyperfine component carrying a line vortex by an optical field. Our system can thus be thought of as "atom-light-skyrmion".
Furthermore, in the proposed trapping configuration the BEC is flowing through the optical tunnel with velocities close to the speed of sound. Therefore we discuss the prospects of using our system for studies of analogue gravity effects. It has been discovered by W. Unruh [3] that within the hydrodynamic approximation [17], the equations of motion for quantised sound waves in a flowing BEC are analogous to those of a scalar quantum field in curved space time. Within this analogy, a sonic horizon, defined as a surface on which the normal component of the fluid flow becomes supersonic, corresponds to the event horizon of a black hole. Unruh proposed that a sonic horizon created in a BEC flow could enable the experimental observation of Hawking radiation [18,19] of sound waves in a BEC.
Here, we show that the BEC flow in the optical tunnel can be driven into the supersonic regime by inserting a laser generated "hump" potential into the tunnel. The perturbation can be engineered to result in the development of a single grey soliton in the condensate flow, which is immobile in the lab-frame. The BEC flow past the soliton is supersonic.
The rest of this paper is organised as follows: In section 2.1 we introduce our model and describe the optical vortex potential. In section 2.2 we explain how the total potential has to be adjusted to allow stabilised ring flows, while we present the initial state of our imaginary time method in section 2.3. In section 3.1 we show some exemplary stable ring flows found numerically, followed by a brief description of a scheme for their experimental realisation in section 3.2. Section 3.3 contains a discussion of the dynamics that arise, when the flow is driven into the supersonic regime. Finally, we devote section 3.4 to an analysis of supersonic flow past a grey soliton in the context of analogue gravity.

Gross-Pitaevskii equation
We describe the dynamics of the BEC wavefunction ψ with the Gross-Pitaevskii equation (GPE) [20]: Using cylindrical coordinates r, z, ϕ, we model a BEC cloud in a cigar-shaped trap with frequencies ω r radially and ω z longitudinally. Throughout the paper we assume that ψ does not depend on ϕ, hence ∆ = ∂ 2 ∂r 2 + 1 r ∂ ∂r + ∂ 2 ∂z 2 . The interaction coefficient U is related to the s-wave scattering length a s : U = 4π 2 a s /m, where m is the atomic mass. Both parameters are taken for 87 Rb: m = 1.44 × 10 −25 kg and a s = 5.67nm. If the wavefunction is written in the polar form ψ = √ ρ exp (iϑ), we can extract the atom number density ρ and velocity v = ∇ϑ/m. For later reference we also define the flow speed v = |v|, the local speed of sound c = Uρ/m, the Mach number M = v/c and the healing length ξ = /( √ 2mc). Finally, V = V t + V v + V h denotes the external potential whose individual components are: The repulsive potential V v can be generated by a blue-detuned optical vortex laser beam [21]. The charge of the optical vortex l is taken to be equal 2. In the dynamical simulations we employ a "hump" potential V h , generated by one or more additional laser beams without any vorticity. They will be used to create small scale modifications of the overall potential experienced by the atoms within the core of the optical vortex. The propagation direction of the laser beams coincides with the long trap axis z. Further, w 0,1,2 are parameters describing the waist of the laser beams at their focus, k 0,1,2 are the wave numbers of the light and the z  Figure 2 shows a slice in the rz plane of the combined potential V of the harmonic trap and optical vortex V = V t + V v . An equal energy contour is indicated by the white lines. In three dimensions the tight focus of the laser beam generates a circular local maximum, whose intersect with the rz plane yields the peaks situated symmetrically around the origin at |r| ∼ 5µm, see figure 2. For a suitable choice of parameters the combined potential exhibits saddle points, marked S in figure 2(a). If the value of the chemical potential is above the potential value at the saddle points, the BEC can classically "wrap" around the intense region of the optical vortex.

Potential design and trial states
We find parameter sets for which the saddle points exist using the following method: First we determine the position of the maximum of the vortex potential V v Eq. (3) with respect to r, keeping z fixed. We thus solve for r. We ignore the harmonic trap in this step. The solution of Eq. (6) is: Now taking V t into account, the height of the potential "ridge" visible in figure 2(a) is approximately given by V (r m (z), z), where V t is now taken into account. Figure 2 (d) shows V (r m (z), z). If saddle points exist, they must be minima of this function.
For V (r m (z), z) to have extrema besides the obvious maximum at z = 0, the quartic expression must have real roots. After selecting the optical vortex parameters and ω r , equation (8) yields a useful approximation for the strongest axial confinement ω z allowed.

Energetically stable flows
In order to find energetically stable ring flows we solve Eq. (1) for V = V t + V v in imaginary time [22] and obtain stationary states of the system. The initial trial function for the imaginary time algorithm must have the same topological structure as the target ring flow. To this end, we multiply an amplitude ρ(r, z) obtained from the Thomas-Fermi approximation [17] with a phase factor exp (iη(r, z)): ψ(r, z, φ) = ρ(r, z) exp (iη(r, z)).
The phase is given by The quantity measures the distance to the single circular ring singularity implicit in Eq. (10). On any closed loop threaded through this singularity, the phase η varies by 2πq, where q is the ring vortex charge. The singularity is located at the point r = D, z = 0 in our coordinates. An energetically stable excited state with high-charge ring flow corresponds to a local minimum of energy. If the trial condensate function is chosen to have the same topology as the energy minimum, we obtain the stabilised flow after a sufficiently long imaginary time [22].

Stationary states
The low density region of a vortex ring in a BEC has a toroidal shape, which matches the shape of the high intensity region of a focussed optical vortex laser beam. In the presence of the optical vortex potential it becomes energetically favourable for the BEC ring singularity to be located at the maximum of laser light intensity. This prevents the contraction of the vortex ring, a mechanism similar to the stabilisation of line vortices within a harmonically trapped BEC cloud [25]. Examples of stationary states resulting from the imaginary time evolution are shown in figure 1. In a harmonic trap with ω r = 7.8 × 2πHz and ω z = 0.5 × 2πHz the BEC contains 2.2 × 10 6 atoms. The subplots are for two different sets of optical vortex parameters and ring flow charges: case (i) and case (ii), as given in the caption. Case (ii) represents the more intense optical vortex. It can be seen in figure 1(a) that the single charge 7 ring singularity that was present in the initial state has broken up into a regular stack of seven unit charge ring singularities. For case (i) they are responsible for the multi-peak structure of the Mach number shown in figure 1(c). If the parameters are altered towards those of case (ii), the influence of the singularities on the flow in the tunnel is reduced. This corresponds to the transition from the solid to the dashed lines in figure 1(b-c).
The integral along the z axis of the flow velocity is directly connected to the ring vortex charge q. Due to the quantisation of circulation in a BEC, its velocity must obey C v · dl = h m q, where C denotes any closed contour threading through all the ring singularities. As the BEC has a much smaller cross-sectional area available when it flows through the tunnel than when it returns on the outer shell of the cloud, the velocity in the tunnel is much higher. The above equation thus approximately becomes: where L denotes the length of the optical tunnel in the z direction.
The presented stability analysis of stationary states is not completely general, as we impose cylindrical symmetry throughout the imaginary time evolution. An important class of ring vortex excitations, the Kelvin modes (see Refs. [23,24] and references therein), could in principle break this cylindrical symmetry. However, a full three-dimensional (3D) stability analysis of skyrmions [15] indicates that the 2D model should give a very good indication of the stability properties: In Ref. [15], we have not found an indication of energetically unstable Kelvin modes of skyrmion's ring singularities. Nonetheless we point out that 3D studies of the optical tunnel would be an interesting, numerically challenging extension of our present work.

Experimental creation scheme
The method of stabilising the ring singularities by pinning, proposed in the present article, should be experimentally realisable with present technology. The procedure to form the ring flow consists of initially creating a BEC at rest, already in the presence of the optical vortex. Then a ring flow structure with high vorticity is seeded by means of phase imprinting. In Ref. [13] J. Ruostekoski and J. R. Anglin show in detail how the phase structure of the vortex ring, similar to the one given by Eq. (10), can be created. Their method uses a superposition of coherent light fields. In the presence of dissipation the BEC will, after the seeding, evolve towards stationary states like those in figure 1.
Compared to the previously known method to stabilise the ring vortex as a part of a skyrmion, our proposal is more flexible. A skyrmion is a topological soliton in a two component BEC, composed of ring vortex in one hyperfine component, which is filled with a second hyperfine component carrying a line vortex. This filling prevents the contraction of the ring vortex in the skyrmion. [14,15]. However, even the stabilisation of a charge two ring flow in a skyrmion requires large atom numbers of the order 1 ×10 7 , and more for higher charges [26]. In comparison, case (ii) represents a charge 9 ring flow in a condensate of only 2.2 × 10 6 atoms.

Supersonic flow
Supersonic flow is expected to be energetically unstable [27]. This manifests itself in the imaginary time evolution as a decay of the total ring flow charge from q 1 to q 2 < q 1 , whenever the parameters are such that a final state with charge q 1 would breach the speed of sound. Despite the energetic instability, it is not generally true that the supersonic flow is also dynamically unstable [28]. Previously, it has been shown that dynamically stable flows with a supersonic region exist in a 1D ring system [29,30], and that they can be obtained from the subsonic ones by an adiabatic change of trapping potential parameters.
Here we use the following sequence for the creation of supersonic flows: Initially the BEC is prepared in the stable state with subsonic ring flow corresponding to case (ii). Then an additional potential as in Eq. (4) is added dynamically. We found that for all potential widths, locations and ramp-up times employed, the generic response of the BEC was dynamical instability and the emission of a single or multiple grey solitons, travelling in the narrow confines of the optical tunnel. Refs. [31,32] report on similar results in 1D simulations. For our 2D simulations (that pertain to a 3D situation) we have employed an adaptive Runge-Kutta-Fehlberg method within the high-level programming language XMDS [33].
According to hydrodynamic theory, a hump potential with a single peak does not allow a quasi-stationary flow that is subsonic in a finite interval and subsonic on either end of this interval [34,35]. Such a flow requires the realisation of a double de-Laval nozzle, using two hump potentials. The first of these renders the flow supersonic, while the second decelerates it back to subsonic velocities. We have numerically examined a variety of dynamical sequences, inserting double de-Laval potentials into the existing ring flows, and found that the result was always the same: The speed of sound is breached initially at the hump farthest in the direction of BEC flow, thus not realising a double de-Laval structure. Nonetheless, we can not exclude that with more careful adjustment of the potential the stationary ring flows presented here could be used as a staging point for double de-Laval nozzles. These have been suggested for the creation of sonic horizons in the context of analogue gravity [35,36].
We now present our results achieved with the dynamical addition of a particularly shaped single peak potential. For the sudden addition of a hump potential, shaped as in figure 3(d), we find the following: Once the flow exceeds the speed of sound a single grey soliton is emitted, and for time spans as long as 0.58s it remains located near z = 0, see figure 3. Only after that time does the grey soliton begin to accelerate. The potential was obtained using Eq. (4) with parameters V h1 = −8.5 ω r , V h2 = 12.5 ω r , w 1 = 1.2µm and w 2 = 0.7µm.
We observe that at all times the flow speed in the tunnel is a function of z only, i.e. it remains quasi one dimensional. It is thus justified to apply 1D soliton theory. We find excellent agreement between the functional form for 1D grey solitons [37,38] and the structures generated dynamically in our simulations. Further the dynamical behaviour of the soliton after its creation can be well understood with the aid of a 1D Lagrangian approach [37,38]. To this end, we apply the Ansatz [37,38] ψ(z, t) = ρ bg (z)e −iµt/ u(z, t), Here µ is the chemical potential, u(z, t) the soliton wavefunction and ρ bg (z) denotes the background condensate density. We work in a frame of reference that moves with the speed of the background flow (v bg ) at the soliton centre (z 0 ). This speed is allowed to vary slowly in time. The parameters in the wavefunction A(t) and B(t) are related to the solitons width and depth and ξ(t) denotes the codensates healing length at z 0 (t). Further we have A(t) + B(t) = 1. The leading order result for the soliton equation of motion following [37,38] is: We find that this equation describes the trajectories of solitons in our simulations well, within uncertainties arising from the numerical extraction of ρ bg and v bg . In particular, Eq. (15) explains the immobility of the soliton shown in figure 3(a-b). In the vicinity of the soliton's initial position the effective potential for the solitons, given by the background density, is approximately flat. Equation (15) suggests that the final travel direction of the soliton can be manipulated between positive and negative z by using potentials which are slightly asymmetric along that axis, which we have confirmed numerically. The potential shape in such a case is still similar to figure 3 (d), but the wells adjacent to the local maximum are of different depths. The BEC density increases on the side with the deeper well, making it the favoured escape direction for the grey soliton.

Prospects of analogue gravity
The purpose of the present section is to consider whether and how grey solitons could play a role in studies of analogue gravity. Our goal here is not to give a rigorous and complete account, but to present basic ideas and stimulate further interest.
The grey solitons obtained in section 3.3 possess a sonic horizon where the BEC flow past them becomes supersonic, termed black hole horizon (BH). They also carry a white hole horizon (WH), the location where the flow returns to subsonic velocities once it has passed the soliton. The separation d between the black-hole and while-hole horizons depends on the soliton size, which is determined by the healing length. We numerically find that d ∼ 2ξ. The conventional derivation of analogue Hawking radiation [4] requires the length scale of variations in the BEC flow, and hence the separation between BH and WH, to be much greater than the condensate healing length. The reason for this is that only modes with phononic wave numbers k such that kξ ≫ 1 propagate in direct analogy to scalar fields in curved space-time [3,4,29,40]. ‡.
There exist interesting phenomena in analogue gravity, which crucially involve white holes in addition to black holes (e.g. the black hole laser effect [43]). However, BEC flow configurations with a widely separated black-and white-hole pair are usually dynamically unstable [29,44], resulting in grey soliton emission from the white-hole [30]. This is also the origin of the solitons discussed in section 3.3. As a grey soliton already represents the end-product of dynamical instability, the black and white holes contained within it are stable, for as long as the soliton exists. Given the instabilities associated with well separated BH-WH pairs, compared to the ease with which one can obtain grey solitons, it would be interesting to enquire, which features of the Hawking or black-hole laser effects persist in a rapidly varying flow. Further motivation for research on quantum emission properties of rapidly varying condensate flows arises from the strict limits on the achievable analogue Hawking temperatures as long as the flow varies slow enough to be considered in the hydrodynamic regime [45]. It was shown that these arise due to three-body losses.
It has been stated that the core ingredients of Hawking radiation are exclusively the existence of an apparent horizon and the validity of the eikonal approximation there [46]. The Mach number around the long time immobile soliton becomes as high as 600 §. The group velocity v g of Bogoliubov excitations with wave number k is given by [40]: Inserting the speed of sound near the soliton, we see that only wave packets comprised of wavelengths shorter than gξ with g = 0.015 have a group velocity in excess of 600c. There is thus a window of wave numbers, π/ξ ≪ k ≪ 2π/(gξ), for which a point of no return exists and the eikonal approximation should be valid in its vicinity. The location of this "horizon" varies depending on the major wavelength of the wave packet. Due to the steep Mach number profile, the "horizons" for different wavelengths ‡ This issue is however a subtle one: All excitations experience a formally infinite blueshift i.e. increase in wave number, when traced back in time towards the horizon. They thus necessarily enter the nonphononic part of the Bogoliubov spectrum. However, it was argued in Refs. [41,42] that the analogue Hawking effect persists despite this fact. § The Mach number is truncated at M = 2 in Fig. 3(a). are spatially close. Quantum emission of quasi-particles thus appears conceivable for wave numbers within the window, however further research is required to determine the existence of an effect.
Another interesting analogue gravity system, for which a BH-WH pair associated with the propagation of a soliton was previously discussed is fermionic superfluid 3 He-A [47]. There the soliton of choice is a domain wall between regions of parallel and anti-parallel orientation of the spin and orbital angular momentum of Cooper-pairs. In such a system the separation between BH and WH can significantly exceed the coherence length [48]. Owing to the structure of the superfluid order parameter near the domain wall, the scenario discussed in Ref. [47] realises a rotating analogue BH, with an ergoregion boundary distinct from the apparent horizon [47], in contrast to the simpler setup under consideration here.
Returning to our condensate scenario, the velocity of the bulk BEC beyond the tunnel is effectively zero, see Fig. 1(c). This provides a region of flat space for the detection of quasi-particles emitted by the horizon. The steering of the soliton's escape direction could be an advantage for spectroscopic detection of horizon radiation [49]. The horizon-soliton can be made to travel towards z > 0, while the generated quantum excitations might be spectroscopically detected in the bulk BEC at z < 0 before the soliton reaches the cloud edge.
One of the advantages of the analogue gravity program over the astrophysical original, is the ability to push the studied system to the validity limits of existing derivations of horizon radiation. This allows the segregation of essential premises from merely analytically convenient assumptions. In this sense our system might one day complement more conventional analogue gravity setups along the lines proposed in Refs. [29,50,35,36,30].

Conclusions
We have introduced a novel arrangement for the studies of persistent ring flows in a BEC. Our method utilises optical vortices, laser beams with a phase singularity at the centre and hence a zero intensity hollow core. If tightly focussed, the highest intensity region of the optical vortex can fit into the low density region of the matter-wave ring vortex in the BEC, stabilising the ring singularity against contraction. We numerically showed that energetically stable BEC states with ring flows of charges up to q = 9 exist in our system.
We have further analysed the response of the ring flows to dynamically varied potentials which drive the flow into the supersonic regime, realising a supersonic optical tunnel for the BEC. In this case, our simulations show the creation of quasi one dimensional grey solitons in the tunnel. The solitons appear due to the dynamical instability once the speed of sound is reached by the velocity of the BEC flow. A single grey soliton that is immobile in the lab-frame for a long time was found for a particular choice of potential.
Finally, we have discussed the prospects of using this single, stabilised soliton as a tool for analogue gravity. We point out that whilst the Mach number around the realised horizon varies too rapidly to expect the usual form of analogue Hawking radiation, the soliton's sonic horizons might be suited to probe the physical boundaries for quantum emission from a horizon.