Excitation of knotted vortex lines in matter waves

We study the creation of knotted ultracold matter waves in Bose-Einstein condensates via coherent two-photon Raman transitions with a $\Lambda$ level configuration. The Raman transition allows an indirect transfer of atoms from the internal state $\left| a \right\rangle$ to the target state $\left| b \right\rangle$ via an excited state $\left| e \right\rangle$, that would be otherwise dipole-forbidden. This setup enables us to imprint three-dimensional knotted vortex lines embedded in the the probe field to the density in the target state. We elaborate on experimental feasibility as well as on subsequent dynamics of the matter wave.


Introduction
Waves can exhibit wavefront dislocations or vortices, i.e. phase singularities of the complex wavefuntion, where the modulus of the wave vanishes, and around which the phase of the wave changes by a multiple of 2π [1]. In [2], it has been shown that Maxwell's equations in fact admit solutions where these lines of dislocation or vortex lines are tied into knots embedded in light fields. Here, the term "knot" refers to an embedding of a circle S 1 into a three-sphere S 3 [3]. Both theoretical [4] as well as experimental advances [5][6][7][8] in three-dimensional light shaping allowed to the experimental realization of knotted vortex lines embedded in optical fields. Apart from optics, the study of knotted topological defect lines and their dynamics has fascinated scientists from diverse settings, including classical fluid dynamics [9,10], excitable media [11,12], chiral nematic colloids [13] to semiconductors [14,15].
Lord Kelvin speculated in 1867 [16] that atoms can be described by vortex tubes in the aether. In fact, almost two decades ago it was suggested that (nontrivial) knots might exist as stable solitons or Hopfions in three-dimensional field theories [17]. A Skyrme model served as an example for a field theory which admits stable knot solitons [18].
Contrasting these studies, numerical investigations on the robustness and centre of mass motion of knotted vortex lines embedded in a single component BEC have been undertaken [26,27]. Such formations are typically unstable, since there is no topological stabilization mechanism and reconnections of vortex lines are allowed due to the occurrence of the quantum stress tensor in the hydrodynamic formulation of the Gross-Pitaevskii Equation. Nucleation [28] and reconnection of vortex lines [29][30][31] and vortex line bundles [32] and subsequent emission of sound waves [33] have been theoretically studied extensively.
In this paper, we follow up on the idea formulated in [34], and discuss a general experimental scheme to create a two-component BEC which contains a knotted vortex line in one of its components. We suggest using a light field containing a knotted vortex line as probe field of a Raman-pulse that drives a coherent two-photon Raman transition of three-level atoms with Λ-level configuration [cf. Fig. 1(a)]. Previously, similar methods have been used to create dark solitons and vortices [34][35][36][37][38] and have been experimentally realized using microwave [39,40] and more recently optical coupling [41]. The pump beam will be the mentioned knotted (stationary) light beam, whereas the control beam is a co-propagating plane wave to remove fast oscillations in the z-direction [34]. The large controllability of the pulse parameters (i.e. strength and duration of the Rabi-pulse) allows for a large controllability of the excitation. The dynamics of two-component BECs can be monitored in real time via in situ measurements without ballistic expansion [42]. Using numerical methods, we study excitation and subsequent dynamics of specific examples of knotted matter waves for experimentally feasible parameters. This paper thus paves the way for experimentally accessing many of the phenomena discussed only theoretically in, e.g. [26,27,43], and if such system were realized experimentally, it would give controlled experimental access to reconnection of vortex lines, subsequent emission of sound waves and more generally quantum turbulence.

Equations of motion for light-matter wave coupling
Consider a three-level atom in Λ-configuration with ground states |a and |b , which are off-resonantly coupled to an excited state |e with detuning ∆ and spatially dependent Rabi-frequencies Ω a (r) and Ω b (z). Assuming ∆ ≫ Ω a , Ω b , the excited state can be adiabatically eliminated, and the dynamics of the condensate wave function confined in a trap V (r) can be described [34,36] with g kj = 4π 2 a kj /m, a kj denoting the scattering length between the species and Fig. 1(a) depicts the level scheme of three-level atoms in Λ-type configuration. The Rabi-frequencies Ω k , k = a, b relate to their corresponding electric fields via where the transition dipole element d has been introduced. As electric fields, we consider a monochromatic, quasi-linearly polarized light field with wavevector k 0 . Then, its slowly varying envelope E is defined by where c.c denotes the complex conjugate. Within the paraxial approximation, the slowly varying envelope E of the beam is governed by where r ⊥ = (x, y), and ∇ 2 ⊥ denotes the transverse Laplacian. Such a system can be realized by considering the two hyperfine ground states |a = |S 1/2 , F = 2, M F = −1 , |b = |S 1/2 , F = 2, M F = 1 of 87 Rb, that are off-resonantly coupled to an excited state manifold. Then, the scattering lengths between the species a kj are given by a ba = 5.5nm, and the ratio a aa :a ba :a bb is given by 1.03:1:0.97 [44]. The motivation for assuming cylindrical symmetry of the trapping potential V (r) is associated with the specific form of the paraxial wave equation Eq. (6). Whereas we (color online) (a) Schematic sketch of a Raman-type transition with Λ-configuration. The state |a is off-resonantly coupled to state |b with two-photon detuning δ. The knotted light field E proportional to Ω a (r) off-resonantly couples state |a to |e with detuning ∆. The final state |b then reflects the involved structure of the light field E associated with Ω a (r). (b) depicts an isosurface of the intensity (purple) of the light-field E at a low isointensity value and a slice of its phase in the (x, y)-plane. Additional to the vortex lines in the center of the beam, the usual diffraction cones become visible. (c) Illustration of the situation at t = 0 for the parameters discussed in 3.2. The isodensity of the cigar-shaped wave function ψ a (yellow) is subject to the knotted light field (purple) from Fig. (b). use different beams for trapping and Raman transition, the aspect ratio ω z /ω for our otherwise independent cigar shaped trap cannot be chosen arbitrarily. Instead, in order to accommodate the optical knot in the BEC (see Fig. 1(c)), we have to make sure that extent of the BEC in the z-direction is large enough. The aspect ratio of the optical beam (and the optical knot) can be expressed using the ratio between the Rayleigh length z r = k 0 σ 2 and the width σ of the light field in the (x, y)-plane: The basic idea towards imprinting knotted vortex lines into BECs is depicted in Fig. 1(a).
Consider the case where all atoms are initially in state |a [ψ b (t = 0) = 0]. Then, assuming Ω i ≪ ∆ and δ ≈ 0 the Rabi pulse coherently transfers population from |a to |b as illustrated in according to Eqs. (1)- (2). The light field is depicted in Fig. 1(b) and together with the density of ψ a in Fig. 1(c). The laser associated with Ω b is chosen to be a simple plane wave, which coherently co-propagates with E, so that the fast variations ∼ e ik 0 z in Eq. (5) are canceled whenever products Ω a Ω * b , as in Eqs. (1)-(2), occur. The explanation as to why such a setup allows to imprint the phase of the structured light field Ω a onto the condensate is the following. For short time and two-photon detuning δ = 0, the small change δψ b to ψ b is given by Thus, we may conclude that at least for small times, it should be possible to populate state |b with a given phase profile using our light field. Once the phase profile θ has been imprinted (E = |E|e iθ ) onto the condensate, the atoms will display motion according to v = ∇θ/m [34]. To realize such a scenario, the intensity of the involved laser beams must be sufficiently strong, such that the time-scales of the imprinting are small compared to time-scales of the dynamics of the condensate. We will use numerical methods and realistic experimental parameters to extend this simple idea beyond the perturbative limit and study its subsequent dynamics.

Knotted Light Field
Laguerre-Gaussian (LG) functions LG l,p form a basis set for solutions for the paraxial wave equation Eq. (6), and are given by the expression with r 2 ⊥ = x 2 + y 2 and L |l| p being the associated Laguerre polynomials. We seek a linear superposition of LG-modes, l,p a lp LG l,p , that describe lightbeams containing knotted vortex lines. To this end, we will review the method proposed in [7], which uses Milnor polynomials [45] as an ansatz for complex light fields in the shape of torus knots to determine appropriate amplitudes a lp , and rescale for application to our setup.
The basic idea [7,46] is to parametrize an N-strand braid as the roots of the polynomial Here, h denotes the height of the periodic braid and N denotes the number of strands or roots of p h (u) in u. Let us choose N = 2 and [46] The projection of the braid onto the (h = 0)-plane leads to a circle. The parameter n represents the number of braid crossings. We will consider the cases n = 2, 3 in the following. A small computation allows us to find an explicit expression for the Milnor polynomial p h in the variables u and exp(ih) =: v, We can now imagine a cylinder containing the braid, and "glue" top and bottom surfaces together to obtain a knot. In fact, one can show [47] that any knot can be represented as a closure of a braid. An easy way to deform our cylinder into a torus and to make it explicitly dependent on r ′ = (x ′ , y ′ , z ′ ) is to write u and v as an inverse stereographic projection from three-dimensional space to a unit three-sphere, Here, r ′ = x ′2 + y ′2 + z ′2 and the units r ′ are non-dimensional. One can easily show that For n = 2, 3, Eq. (15) describes a Hopf link and a trefoil knot, respectively, which represent the simplest non-trivial examples of a link and a knot, respectively. The latter two will be used as exemplary fields in the following. In order to use Eq. (15) to find appropriate coefficients a l,p , let us rescale the light field Eq. (9) by R s to nondimensional units (x ′ , y ′ , z ′ ), such that to equate the light field to our knot in the same dimensionless units. Effectively, R s scales the abstract unit sphere to have a transverse extent of approximately R s with respect to the laser beam. Hence, there are two length scales of our system, R s and σ, that are associated with the nodal lines of the knot and the transverse extent of the beam. As we will see, the ratio w = σ/R s will play a crucial role as an important degree of freedom, that allows us to change the width of the beam relative to the positions of the nodal lines of the knot. To find appropriate superpositions of LG-modes, it is sufficient to restrict considerations to the (x, y)-plane only. Let us equate Eq. (15) with a superposition of the above-mentioned rescaled version of Eq. (9): Comparing different powers in r ′ ⊥ allows us to determine the finite number of coefficients a l,p uniquely, which depend only on the real number w. Whereas this ansatz can be used for some knots, there are counterexamples, and there is no rigorous general proof as to why and in which cases this ansatz leads to success. A large value of w ensures that the vortex lines are actually embedded in the beam, and not chopped off. Furthermore, for finite w, additional vortex lines in the shape of hairpins appear for larger z-values, which leads to the fact that w must be chosen large enough. On the other hand, if w is too large, the polynomial increase in the polynomials describing the knots will not be attenuated quickly enough by the Gaussian, and thus intensity variations become huge, which is undesirable for our setup. It is possible [7] to further optimize these coefficients to separate vortices with regions of larger intensity. Once appropriate coefficients a l,p have been found, the light field can be written down as superposition of LG modes by rescaling back into physical units. For the sake of clarity, let us introduce an auxiliary function f defined as f (r ⊥ /σ, z/z r )/σ := LG σ,zr l,p (r ⊥ , z). Then, we find Here, the amplitude A of the light field has been introduced, which gives the intensity of the beam the right value in appropriate units. Note, that Eq. (21) is no longer dependent on the choice of R s .

Rescaling
Let us rescale Eqs. (1)-(2) by introducing spatial r ′ and time t ′ coordinates rescaled by oscillator length a 0 = /mω and trapping frequency ω, respectively, Then, after multiplying Eqs. (1)-(2) by 1/mω 2 √ a 0 , Eqs. (1)-(2) become where we left away the prime for our non-dimensional time and space variables and introduced γ z = ω z /ω, and κ kj = 4πa kj N/a 0 . We can express as κ kj = a 2 0 /2a 2 h using the healing length a h = 8πa kj N/a 3 0 . Since the aspect ratio γ z of the trapping frequencies is typically small, it is more convenient for numerical studies to rescale the elongated z-axis as follows: This rescaling finally yields our equations of motion: where we have dropped the primes for convenience. Furthermore, applying the same rescaling to the paraxial wave equation [Eq. (6)], we find the following expression in our rescaled units: Here, k 0 has been redefined and stands for k 0 = 2πa 0 γ z /λ. Then, the functions defined in Eq. (9) remain solutions to Eq. (31) in the new coordinates, and using our modified k 0 , we can still write the Rayleigh length as z r = k 0 σ 2 , where z r is measured in units of a 0 /γ z and σ in units of a 0 . We have m = 1.44×10 −25 kg for 87 Rb, so that a trapping frequency of ω = 2π×10Hz leads to an order of magnitude of a 0 = 3.4µm for the typical transverse extent of the BEC. An order of magnitude estimate for the off-diagonal coupling terms is given by for Ω i ≈ 1-10MHz and Ω i /∆ = 0.01-0.1. The required tight focusing of the beam leads to the fact, that we only need moderate beam powers to achieve adequate Rabifrequencies.
In the following, we set the two-photon detuning δ = 0 to achieve optimal transfer. Finally, we need to find an estimate for the Rabi pulse duration t d . To this end, consider the simple case, where Ω a and Ω b describe driving by plane waves. We seek to find the optimal time for maximal transfer of atoms from |a to |b . Given that the "off-diagonal" terms in the coupled equations Eqs. (29)-(30) are much larger than the diagonals, we can approximate the population N b in state |b as This expression allows us to find the optimal (non-dimensional) time duration of the Rabi pulse t d that leads to a complete transfer of population: On the other hand, this consideration does not apply to our case due to the spatial dependence of the Rabi frequency Ω a (r). Due to spatial dependence of the intensity variations, we do not really have an optimal pulse duration t d , and our choice of t d is a compromise between on one hand having sufficiently large transfer of atoms to state |b and on the other hand avoiding population being transferred back from state |b to state |a in positions where the intensity is large and thus the transfer dynamics is faster. Hence, we have to choose shorter pulse durations than what Eq. (34) suggests. Furthermore, due to the nodal lines in the pump beam, it is never possible to use this arrangement for a total conversion of atomic population from |a to |b . With these aspects kept in mind, we may still use Eq. (34) as a rough estimate for the order of magnitude for our choice of t d .

Numerical Investigation
In the previous section, we elaborated on how to inscribe complex knotted vortex lines into matter waves. Whereas the simple argument of Eq. (8) allows us to conclude that our setup should work at least in the perturbative regime for short time-scales and small amounts of atoms in state |b , we cannot infer what happens beyond this perturbative regime. In this section, we aim to numerically study excitation and decay of our highly excited states into more elementary unknots or vanishing of the latter by collision and annihilation of vortex lines. We will do so by considering specific examples. However, these examples are by no means an exhaustive treatment of the the huge potential manifold of realizations possible. To thoroughly understand the dynamics remains an elusive goal, which goes beyond the scope of this paper.

Hopf-link
One of the simplest torus knots is the so-called Hopf-link, which corresponds to setting n = 2 in Eq. (15) and is depicted e.g. by the green isodensity surface in Fig. 2(b). Using the method outlined in 2.2 and Eq. (17), we find the following superposition coefficients of Laguerre-Gaussian polynomials for a Hopf-link: To this end, let us consider the dynamics for κ ab = 4πa ab N/a 0 = 1887, γ z = 0.125. We set the two-photon detuning to δ = 0. This can be achieved by using typical parameters of 87 Rb with ω = 2π × 10Hz, ω z = 2.5πHz, N = 93000, a ab = 5.5nm, λ = 780nm, which amounts to typical units of length scales of a 0 = 3.4µm in the (x, y)plane and a 0 /γ z = 27.2µm in the z-direction. Thomas-Fermi like dynamics, where the nonlinearity is large compared to the broadening due to the Laplacian, is preferable to ensure robustness of the vortex cores and avoid the trivial broadening of the latter. For that reason, we have to choose a sufficiently large value for κ ij . For the other a ij , we assume the ratio a aa :a ab :a bb to be given by 1.03:1:0.97 [44]. With respect to our Rabi pulse, we choose A = 450, w = 1.5, σ = 0.675 corresponding to 2.3µm, and a pulse duration of t d = 0.002 corresponding to 31.8µs. We use the common Fourier split-step method [48] and adaptive Runge-Kutta algorithm to fourth order for the time-step for numerical computation. Parts of the code were written using [49]. To find the appropriate state ψ a (t = 0), we used imaginary time-evolution (e.g. [50]). It is crucial to use a variable time-step, since the dynamics during the time duration of the Rabi-pulse t d has to be temporally fully resolved, and thus the time-step has to be much smaller than t d within t d . After t > t d , the time-step can be chosen much larger, since it only needs to resolve the characteristic time-scale of the BEC dynamics. We used time-steps varying between ∆t = 10 −8 -5 × 10 −4 and a spatial resolution of ∆x = 0.05 − 0.06 with 220 3 points.  Fig. (a), the green surface represents a low-value isodensity surface of |ψ b | 2 around the vortex line imprinted on the BEC. (c) Projection of Fig (b) into the (x, y)-plane and associated indices according to Eq. (36). Expected dynamics should thus be, that the central region moves faster in the z-direction relative to its neighbouring regions. (d) shows again the isosurface from Fig. (b), however, the coloring of the isosurface illustrates the value of the phase at each position of the isosurface. Fig. 2(a-d) depict an isodensity surface of the condensate wave function ψ b with small value, which serves as an illustration of the knotted vortex core, and its phase profile in different ways, shortly after the imprinting occurred. Here, the initial wave function ψ a (t = 0) was computed as the ground state of the trap and ψ b (t = 0) = 0.
A basic qualitative anticipation of the dynamics can be found by projecting the knot orthogonally to the beam-or z-axis and considering the momentum associated with the (unit-)areas, as elaborated in [51]. Consider the projection in the (x, y)-plane, as shown in Fig. 2(c). The index I j assigned to each of the areas R j enclosed by the vortex lines γ can be found by evaluating the sum [51] which runs over all intersections of the chosen vector ρ with the projected vortex line γ for a given region R j . Here, ρ denotes an arbitrary vector pointing from the inside to the outside of the enclosed area and sign is the usual sign-function taking the values sign( * ) = ±1. Furthermore, the vector t is a tangent vector to the vortex line curve (denoted as γ) at the point where ρ crosses γ. The direction of t is given by the orientation of the arcs, which is determined by the phase [see Fig. 2(c)]. The values assigned to the regions in Fig. 2(c) correspond to the values computed by Eq. (36). The associated momenta of a region R j can be found by [51] ( where ω denotes the vorticity. What one can deduce from these arguments is a movement in the positive zdirection with the central region traveling fastest [ Fig. 2(c)]. Let us now look at the numerically computed dynamics. Snapshots of the latter are shown in Fig. 3 (see also movies hopf 1.mp4 and hopf 2.mp4 ). Clearly, there is an overall center-of-mass mass motion towards the positive z-direction, as expected. However, reconnection of 1 (corresponding to ∆t = 1.6ms) starting from t = 0.05 (corresponding to t = 0.8ms). Both rows use the same isosurface for |ψ b | 2 , however, in the lower row we cropped the box at a smaller value, so that the "boundary" of the condensate does not obfuscate the dynamics of the vortex lines. Additional to that, we used the phase to color the isosurface in the lower row. We can see an overall drift the positive zdirection, which can be understood from Fig. 2. See also the movies hopf 1.mp4 and hopf 2.mp4 for the full propagation dynamics.
vortex lines as well as finite size effects of the condensate leading to sharp gradients in the density lead to a very involved dynamics, which we were not able to predict or understand in simple terms. Surprisingly, the formation decays into two unknots that propagate into the positive and negative y-direction, respectively.

Trefoil
In this section we will study the dynamics of a trefoil imprinted on a BEC. The dynamics is generic, quite different initial conditions lead to similar results, the basics of which can be understood in the simple terms described before.
In this case, let w = 1.2. Instead of using the prefactors a l,p (w) from Eq. (17), we use the optimized prefactors given in [7,46] of the Laguerre-Gaussian modes for the light field.
To this end, let us use σ = 0.78 corresponding to 2.65µm and A = 700. Furthermore, let us choose κ = 3753, γ z = 0.05. This can be realized using again ω = 2π × 10Hz, ω z = 1πHz, N = 1.85 × 10 5 , a ab = 5.5nm and λ = 780nm. With respect to our Rabi pulse, we use a pulse duration of t d = 0.004 corresponding to 63.7µs. Fig. 4(a) illustrates the phase of ψ b and (b) its vortex core shortly after imprinting. Projecting onto the (x, y)-plane and assigning values according to Eq. (36) leads to Fig. 4(c). Again, we can expect that the overall knot will propagate in the positive z-direction, and the central part carries the largest momentum. Fig. 4(d) illustrates the expected reconnection dynamics according to the rules found in e.g. [29]. The arrows  Fig (b) into the (x, y)-plane and associated indices. Since the indices reflect the momenta of the regions, we expect that the central region moves faster in the z-direction relative to its neighboring regions. (d) shows the expected breakup of the green vortex lines into the solid black lines. The black arrows illustrate the orientation, which is determined by the phase. The grey dashed inset illustrates the expected product of collision and reconnection of the vortex lines according to, e.g. [29].
indicate the direction of t, which is again determined by the phase. The inset (grey dashed box) illustrates how the vortex lines should reconnect. If we apply this simple rule to all three crossings, we can expect that the green vortex line evolves into what is depicted by the solid black line, i.e. a decay into two unkots.
Let us now consider the dynamics, which is shown in Fig. 5 (see also the movie trefoil 1.mp4 and trefoil 2.mp4 for full dynamics). We see that the vortex lines of this specific trefoil knot first reconnect [see Fig. 5], and the central regions travels fastest into the positive z-direction according to our expectation Fig. 4(b-c). After that, the central region expels an unknot or vortex ring which decouples from the rest of the knot, as shown in Fig. 5. The single unknot propagates to the positive z-direction faster than the remaining part of the knot. Interestingly, our dynamics differs from the usual dynamics of a clear breakup into two unknots as expected from [ Fig. 4(c), solid black line] and what has been found in [26,43]. Differences to previous observations are due to finite size effects of the cloud, that in our case the knot has a large aspect ratio, which with our rescaling basically means that the dynamics in the z-direction is much slower, smallness of the knot, and finally that we actually consider two coupled fields (dynamics of ψ a not shown).

Conclusions
Recently there has been a growing interest in the dynamics of knotted vortex lines in BECs. However, thus far there has been no actual proposal on how to excite such matter waves in a controlled fashion. In this work, we have presented a setup with physically realistic parameters to create such knotted vortex lines in ultracold matter waves. We combined recent theoretical and experimental results from complex light shaping which allowed us to create knotted vortex lines embedded in light fields. We use the latter (color online) Dynamics of the condensate wavefunction ψ b , whose vortex lines form a trefoil knot. Both upper and lower row show the same dynamics in timesteps of ∆t = 0.1 (corresponding to ∆t = 0.8ms) starting from t = 0.05 (corresponding to t = 0.8ms). Both rows use the same isosurface for |ψ b | 2 , however, in the lower row we cropped the box at a smaller value, so that the "boundary" of the condensate does not obfuscate the dynamics of the vortex lines. Additional to that, we used the phase to color the isosurface in the lower row. Upon evolution, vortex cores reconnect, which leads to a decay into a vortex ring being expelled from the central region, which can be understood from Fig. 4). See also the movie trefoil 1.mp4 and trefoil 2.mp4 for the full propagation dynamics.
as a probe field for three-level atoms in a Λ-type setup to inscribe the nodal lines to BECs. The setup is quite generic, and allows investigation of a large variety of potential dynamics. The finite time span to populate state ψ b reduces production of soundwaves. The finiteness of the condensate, the smallness of the knot as well as the reconnections of the vortex lines give rise to a very involved dynamics. The velocity in the z-direction of the vortex knot should not be large. This velocity can be controlled by the ratio between the trapping frequencies γ z = ω z /ω relative to the condensate. The probe field which determines Ω a should not be too spatially broad, otherwise the value of γ z required to fully embed the knotted vortex line becomes impractical. We have assumed that the light field can be treated within the paraxial approximation, although we note that, in principle, this condition can be relaxed to non-paraxial knotted light fields [52]. Due to the tight focusing of the optical beam, only moderate powers of the light fields that determine the Rabi-frequencies Ω i are required. Thus, the setup should work generically for a large range of parameters. We have shown two illustrative examples of a trefoil knot and a Hopf-link, and discussed the dynamics using already well-established techniques. The data presented in this paper is available online at [53].