text:
# Bose-Einstein condensates with a bent vortex in rotating traps

## I Introduction

Several experimental groups have now produced vortices in Bose condensates of atomic gases by a rotation of the trapping potential ens ; MIT ; Foot . It is therefore important to characterize the steady state vortex configurations for parameters relevant to the experiments. In particular the regime of the ENS and MIT still deserves some investigation. It corresponds to a condensate trapped in a cigar-shaped harmonic potential, that is with a harmonic oscillation frequency along the rotation axis $`z`$ smaller by one order of magnitude than the (quasi-degenerate) oscillation frequencies in the $`xy`$ plane. Furthermore, a configuration with a single vortex can be produced in the experiment in a reproducible way ens . The single vortex configuration in a cigar-shaped trap is the main object of the present work.

Why is the single vortex cigar-shaped regime so interesting? First the weak harmonic confinement along $`z`$ makes the system very different from previously studied configurations. If the harmonic confinement along $`z`$ was stronger than in the $`xy`$ plane the vortex configurations would be similar to what happens in a 2D rotating Bose gas, a well studied limiting case Fetter\_review ; Rokhsar ; castindum ; du . If the harmonic confinement along $`z`$ was absent, e.g. if the condensate was a cylinder, one would face a situation typical of superfluid helium II, a well studied subject Donnelly .

But now the gas keeps its 3D character while being spatially inhomogeneous along $`z`$. This has the amusing consequence that a single vortex line has a tendency to bend. First the steady state condensate can have a bent vortex line, as shown numerically by a minimization of the Gross-Pitaevskii energy functional in garcia  and as obtained also numerically with an approximate vortex line energy functional in AftaRiv . Second, even a condensate with a straight vortex line has low energy Bogoliubov modes corresponding to large fluctuations of the end points of the vortex line Fetter\_dynamics ; feder  so that the vortex line can easily bend in presence of thermal fluctuations.

After a presentation of the model considered in this article, see section II, we give a systematic numerical study of the steady state bending as function of the trap aspect ratio, see section III, an analytic understanding of the bending of the steady state vortex line in a cigar-shaped trap, see section IV, and the description of the finite temperature fluctuations of the vortex line, see section V, including a discussion of the effect of vortex line bending on the experimental absorption images.

## II Model and basic assumptions

Let us consider an almost pure Bose-Einstein condensate of $`N`$ atoms confined in a trapping potential rotating along axis $`z`$ at angular velocity $`\mathrm{\Omega }`$. The thermodynamically metastable equilibrium configurations of the system correspond to local minima of the Gross-Pitaevskii energy functional in the rotating frame castindum

$$E[\varphi ,\varphi ^{}]=d^3\stackrel{}{r}\left[\varphi ^{}\left(H_0\mathrm{\Omega }L_z\right)\varphi +\frac{Ng}{2}|\varphi |^4\right],$$
(1)
where the condensate wave function $`\varphi `$ obeys the normalization condition

$$\varphi ^2\varphi |\varphi =d^3\stackrel{}{r}|\varphi |^2=1.$$
(2)
The Hamiltonian operator $`H_0`$ in Eq. (1) contains the kinetic and trapping potential terms

$$H_0=\frac{\mathrm{}^2^2}{2m}+U(\stackrel{}{r})$$
(3)
and $`L_z=i\mathrm{}(x_yy_x)`$ is the angular momentum operator along $`z`$. Here we will consider harmonic trapping potentials $`U(\stackrel{}{r})`$ with an adjustable slight anisotropic deformation in the $`xy`$ plane:

$$U(\stackrel{}{r})=\frac{1}{2}m\omega _{}^2\left[(1ϵ)x^2+(1+ϵ)y^2\right]+\frac{1}{2}m\omega _z^2z^2.$$
(4)
For the choice of parameters we refer to the typical values of the recent ENS experiments with 87Rb atoms (scattering length $`a=100a_05.29\times 10^9m`$) ens : total number of atoms $`N=1.4\times 10^5`$, axial frequency $`\omega _z=2\pi \times 11.7\mathrm{Hz}`$ and anisotropy parameter $`ϵ1`$. In this paper the radial frequency $`\omega _{}`$ is varied in order to explore a wide range of trap anisotropies: $`\lambda ^1\omega _{}/\omega _z[1,15]`$. In what follows energies and lengths are given in units of $`\mathrm{}\omega _{}`$ and $`a_{}=\sqrt{\mathrm{}/m\omega _{}}`$ respectively, $`m`$ being the atomic mass of 87Rb.

## III Local minima of energy: numerical results

In this section we discuss the method and the algorithm employed to minimize numerically the energy functional (1). Then we present the results obtained for the stationary configurations with and without vortices in a long cigar trap, and finally we discuss the existence domain for a single vortex configuration for several trap geometries.

### III.1 Method

First of all, we reformulate the energy functional in order to account more easily for the normalization constraint (2): we define $`\varphi =\psi /\psi `$ so that the value of $`E`$ can be obtained for condensate wave functions $`\psi `$ with a norm different from unity. This corresponds to dividing the terms of $`E[\psi ,\psi ^{}]`$ quadratic in $`\psi `$ by $`\psi ^2`$, and the interaction term quartic in $`\psi `$ by $`\psi ^4`$. The modified energy functional reads

$$E[\psi ,\psi ^{}]=d^3\stackrel{}{r}\frac{\psi ^{}\left[H_0\mathrm{\Omega }L_z\right]\psi }{\psi ^2}+\frac{Ng}{2}\frac{|\psi |^4}{\psi ^4}.$$
(5)
Then we discretize $`\psi `$ on a three-dimensional grid with periodic boundary conditions in position and in momentum space. The number of grid points that we use ranges from 64 to 256 for each lattice side, depending on the trap geometry (a typical choice is $`64\times 64\times 192`$ for long cigar configurations and $`96\times 96\times 64`$ for spherical geometries with vortices). The minimization is performed by using the conjugate gradient algorithm described in numrec . This algorithm is in general much more efficient that the usual steepest descent method (see Ref. numrec  for a comprehensive discussion). One ingredient of the conjugate gradient method is a line minimization of the energy functional, that is the minimization of $`E`$ along the line $`\psi =\psi _0+\lambda \delta \psi `$, $`\lambda `$ being a real parameter, where $`\psi _0`$ is the current trial wave function and $`\delta \psi `$ is the proposed direction (or gradient) along which to move the trial wave function. An important issue here is to find the first minimum encountered when moving downhill in energy along the line: if the algorithm can jump to another minimum on the line, corresponding for example to an energy valley with a different number of vortices, one gets wrong predictions, in the sense that one does not fully explore the stability domain of a branch of solution with a given number of vortices. This issue is usually not considered as important in the textbook implementation of the conjugate gradient method: e.g., in numrec  the line minimization does not necessarily find the first minimum. What we have done is to use the fact that $`E`$ is a rational function of $`\lambda `$ so that it is completely characterized by the coefficients of the monomials in $`\lambda `$ in the numerator and the denominator. We then easily find the roots of $`dE/d\lambda `$ and the first local minimum of $`E`$ encountered when one moves along the line downhill in $`E`$ starting from $`\lambda =0`$. Afterwards we proceed with another line minimization along a conjugate direction numrec , and so forth until we find a local minimum of the energy functional (5) in the full configuration space spanned by the wave function $`\psi `$. As convergence criterion we use

$$H_{\mathrm{gp}}\psi \mu \psi <\gamma \mu \psi $$
(6)
where $`H_{gp}`$ is the Gross-Pitaevskii Hamiltonian

$$H_{\mathrm{gp}}H_0+Ng|\psi |^2/\psi ^2$$
(7)
and $`\mu `$ is given by

$$\mu =d^3\stackrel{}{r}\frac{\psi ^{}H_{\mathrm{gp}}\psi }{\psi ^2}$$
(8)
which eventually converges to the chemical potential. The parameter $`\gamma `$ is a small parameter ($`10^8÷10^{10}`$) controlling the quality of the convergence.

### III.2 Results for a fixed trap geometry: case of a long cigar

We start by considering a cigar shaped trap geometry typical of the ENS experiments in which vortex configurations are nucleated ens ; feder . In particular we set $`ϵ=0.03`$ and $`\omega _{}=2\pi \times 174`$Hz, giving $`\lambda ^114.9`$. To explore the configuration space we start from a trial wave function for a fixed value of the rotation frequency $`\mathrm{\Omega }`$ until an equilibrium configuration (i.e. a local minimum of the energy functional (5)) is found. Then we increase or decrease slightly $`\mathrm{\Omega }`$, and then find the new local minimum. In this way we can follow continuously branches of configurations with or without vortices, depending on the rotation frequency and on the path followed. In Figure 1 we show the energy $`E`$ and the angular momentum per particle $`L_z`$ along the rotation axis of each configuration, as function of the rotation frequency $`\mathrm{\Omega }`$, up to 4-vortex configurations.

We notice that, for this trap geometry, when a single vortex appears (i.e. when it becomes a stable configuration) it has immediately a lower energy than the 0-vortex configuration, differently to what happens in 2D castindum . Moreover this vortex configuration is characterized by a bending of the core line, a phenomenon already obtained with a different numerical technique in garcia , and this accounts for the fact that the angular momentum per particle $`L_z`$ is lower than $`\mathrm{}`$ ($`L_z`$ is equal to $`\mathrm{}`$ for a straight centered vortex).

To better investigate this aspect we have studied the deformation of the vortex line as function of the rotating frequency. If we start with a single vortex configuration and we decrease the rotation frequency $`\mathrm{\Omega }`$ the bending of the core line becomes more pronounced, until we reach a certain critical frequency $`\mathrm{\Omega }_1`$ and the system jumps to the 0-vortex configuration. In the opposite direction, when $`\mathrm{\Omega }`$ increases, we find that the vortex line tends to become more straight as expected (see Figure 2), but then, whereas a small bending is still present at the extremities of the condensate, it “decays” into a three vortex configuration. The value of $`\mathrm{\Omega }`$ at which this happens defines a second critical rotation frequency $`\mathrm{\Omega }_2`$.

Therefore for the long cigar trap and the parameters considered here, one never obtains a single, straight vortex: when $`\mathrm{\Omega }`$ is too large, other vortices come in. As we will show in the following, this is due to destabilization of high angular momentum surface modes dalfovo , and corresponds to the fact that the single vortex configuration ceases to be a local minimum of energy when the energy of these modes becomes negative (i.e. there is at least one direction along which the energy is going down in the functional space). We discuss this effect in details in the next subsection.

### III.3 Destabilization of surface modes

We consider in this subsection the case of a cylindrically symmetric trapping potential (with a vanishing asymmetry parameter $`ϵ=0`$) so that eigenstates of the $`N`$-body Hamiltonian have a well defined angular momentum. The expression for the energy $`_l`$ of an excitation mode of the condensate with angular momentum $`\mathrm{}l`$ in the rotating frame is

$$_l=\mathrm{}\omega (l)l\mathrm{}\mathrm{\Omega },$$
(9)
where $`\omega (l)`$ is the excitation energy of the same mode in the laboratory frame castindum . This determines the rotation frequency $`\stackrel{~}{\mathrm{\Omega }}_l\omega (l)/l`$ at which the energy of such mode becomes negative in the rotating frame. In particular the minimum value over $`l`$ of $`\stackrel{~}{\mathrm{\Omega }}_l`$ defines the Landau critical rotation frequency $`\stackrel{~}{\mathrm{\Omega }}`$ at which the vortex becomes thermodynamically unstable due to destabilization of these surface modes

$$\stackrel{~}{\mathrm{\Omega }}=\underset{\{l\}}{\mathrm{min}}\left(\frac{\omega (l)}{l}\right).$$
(10)
The value $`\stackrel{~}{\mathrm{\Omega }}`$ can be estimated analytically by using the sum rules approach which provides an estimate for $`\omega (l)`$, as discussed in Refs. dalfovo ; zambelli . In order to describe the surface modes, we consider the excitation operators for $`l0`$

$$F_\pm =\underset{k=1}{\overset{N}{}}(x_k\pm iy_k)^l.$$
(11)
Then one defines the moments $`m_p^\pm `$

$$m_p^\pm _0^+\mathrm{}𝑑E\left[S_+(E)\pm S_{}(E)\right]E^p$$
(12)
where the strength distribution functions $`S_\pm `$ are given by

$$S_\pm (E)\underset{n}{}|n|F_\pm |0|^2\delta (E\mathrm{}\omega _n).$$
(13)
The states $`|n`$ form a complete set of eigenstates of the Hamiltonian operator $`H`$ for our system of $`N`$ interacting particles confined by the trapping potential:

$$H=\underset{k=1}{\overset{N}{}}\left[\frac{p_k^2}{2m}+U(\stackrel{}{r}_k)+g\underset{j=1}{\overset{k1}{}}\delta (\stackrel{}{r}_j\stackrel{}{r}_k)\right].$$
(14)
Here we assume the two mode approximation, as discussed in Ref. zambelli :

$$S_\pm (E)=\sigma \delta (E\mathrm{}\omega _\pm ),$$
(15)
where the two modes are equally weighted (they have the same strength $`\sigma `$) due to the vanishing of the $`m_0^{}`$ momentum zambelli . This is justified by the fact that in the large $`N`$ limit the strength distributions are dominated by the contribution of two modes with energy $`\mathrm{}\omega _\pm `$ and angular momentum $`\pm \mathrm{}l`$ excited respectively by the operators $`F_\pm `$. With this ansatz it is straightforward to obtain the relation between the frequencies $`\omega _\pm `$ of such modes and the first three moments in Eq. (12); in particular we have

$`\delta `$ $``$ $`\omega _+\omega _{}={\displaystyle \frac{m_2^{}}{m_1^+}}`$ (16)
$`\overline{\omega }^2`$ $``$ $`\left({\displaystyle \frac{\omega _++\omega _{}}{2}}\right)^2={\displaystyle \frac{m_3^+}{m_1^+}}{\displaystyle \frac{3}{4}}\delta ^2.`$ (17)
Notice that only the mode with frequency $`\omega _+=\overline{\omega }+\delta `$ is relevant to our discussion, since the energy of the mode with angular momentum $`\mathrm{}l`$ increases for increasing rotation frequency $`\mathrm{\Omega }`$, and can never become negative.

The next step is to calculate the moments $`m_p^\pm `$ from Eqs. (12) and (13) by using closure relations and then explicitly evaluating the commutators

$`m_1^+`$ $`=`$ $`[F_{},[H,F_+]]={\displaystyle \frac{2N\mathrm{}^2}{m}}l^2r_{}^{2l2}`$ (18)
$`m_2^{}`$ $`=`$ $`[[F_{},H],[H,F_+]]={\displaystyle \frac{4N\mathrm{}^3}{m^2}}l^2(l1)r_{}^{2l4}L_z`$ (19)
$`m_3^+`$ $`=`$ $`[F_{},H],[H,[H,F_+]]={\displaystyle \frac{2N\mathrm{}^4\omega _{}^2}{m}}l^3`$ (20)
$`\left[r_{}^{2l2}+{\displaystyle \frac{l1}{m^2\omega _{}^2}}{\displaystyle \frac{p_+r_{}^{2l2}p_{}+2\mathrm{}(l2)r_{}^{2l6}L_z+2(l2)r_{}^{2l6}L_z^2/l}{r_{}^{2l2}}}\right]`$
where $`\mathrm{}`$ stands for the expectation value in the state $`|0`$. These expressions generalize those in dalfovo  (case of a vortex free condensate) and in zambelli  (case of excitations of $`l=2`$ of a condensate with vortices). We note that the use of $`p_\pm p_x\pm ip_y`$ and $`r_\pm x\pm iy`$ rather than the Cartesian coordinates $`p_x`$, $`p_y`$, $`x`$, $`y`$, e.g. by writing the kinetic energy as $`p_+p_{}`$, considerably simplifies the calculations of the commutators.

For the case of a straight vortex with charge $`q`$ the key quantities to be inserted in Eqs. (16) and (17) are

$`{\displaystyle \frac{m_2^{}}{m_1^+}}`$ $`=`$ $`2q(l1)\omega _{}\left[{\displaystyle \frac{r_{}^{2l4}}{r_{}^{2l2}}}\right]`$ (21)
$`{\displaystyle \frac{m_3^+}{m_1^+}}`$ $`=`$ $`\omega _{}^2l\left[1+(l1){\displaystyle \frac{p_+r_{}^{2l4}p_{}}{r_{}^{2l2}}}+{\displaystyle \frac{2q(l+q)(l2)(l1)}{l}}{\displaystyle \frac{r_{}^{2l6}}{r_{}^{2l2}}}\right],`$ (22)
where quantities within square brackets are expressed in dimensionless units, as defined in Sect. II.

In Figure 3 we plot the curve $`\omega (l)/l`$ for a cylindrically symmetric trap ($`ϵ=0`$), for the case $`q=0`$ (ground state without vortices) and $`q=+1`$ (straight vortex). In the same figure we also indicate the critical frequencies $`\mathrm{\Omega }_2^0`$ and $`\mathrm{\Omega }_2`$ which define the upper bound of the existence domain respectively for states without vortices and with a bent vortex, as found from the numerical minimization of the energy functional in presence of a slight anisotropy as considered in the previous subsection ($`ϵ=0.03`$). We notice that the sum rule prediction gives a value of $`\mathrm{min}\left(\omega (l)/l\right)`$ which is quite close to our numerical result. This supports the statement that the surface modes are indeed responsible for the destabilization of the single-vortex state and the vortex-free state.

### III.4 Effect of a ballistic expansion

In general, in the experiments, the condensate density is imaged after a time-of-flight: the trapping potential is switched off and the cloud expands freely for some time before being imaged by the absorption of a laser beam. This has the advantage of increasing the diameter of the vortex core, making it larger than the optical wavelength. Some theoretical work is however needed to relate the images obtained after the expansion to the density profile of the initial trapped condensate.

To calculate the expansion of the condensate after the release from the trap, we propagate the initial wavefunction according to the full 3D Gross-Pitaevskii equation in real time after having abruptly switched off the trapping potential, and after having performed the scaling and gauge transform of Castin-Dum ; Kagan . In terms of the rescaled spatial coordinates $`\stackrel{~}{r}_ir_i/\lambda _i(t)`$, the modified wave function $`\stackrel{~}{\psi }(\stackrel{}{\stackrel{~}{r}},t)`$ satisfies the equation

$`i_t\stackrel{~}{\psi }(\stackrel{}{\stackrel{~}{r}},t)`$ $`=`$ $`[{\displaystyle \frac{\mathrm{}^2}{2m}}{\displaystyle \underset{j}{}}{\displaystyle \frac{1}{\lambda _j^2}}_{\stackrel{~}{r}_j}^2`$
$`+{\displaystyle \frac{1}{_j\lambda _j}}(U(\stackrel{}{\stackrel{~}{r}})+g|\stackrel{~}{\psi }|^2\mu )]\stackrel{~}{\psi }(\stackrel{}{\stackrel{~}{r}},t)`$
and the scaling parameters $`\lambda _j(t)`$ are solutions of Castin-Dum

$`\ddot{\lambda }_j`$ $`=`$ $`{\displaystyle \frac{\omega _j^2(0)}{\lambda _j\lambda _1\lambda _2\lambda _3}}`$ (24)
with initial values equal to unity and with vanishing initial derivatives. The use of these scaling equations allows us to calculate numerically the expansion of any vortex configuration by using the same grid as the one used for the stationary trapped state as most of the evolution of $`\psi `$ due to the ballistic expansion is absorbed by the scaling and gauge transform. In a 2D model with an axially symmetric trapping potential the net effect of the ballistic expansion on the density corresponds exactly to a scaling transformation and a finite angle rotation of the condensate eigenaxes castindum . Here we expect this property to remain approximately true as the condensate is very elongated along $`z`$ and expands much more rapidly radially than axially.

A first application of this calculation is the analysis of the column density of the condensate, that is the density integrated along $`z`$, after ballistic expansion. There is indeed an important issue concerning the mechanisms which are responsible for the reduction of contrast for the vortex configurations which are observed in the experiments ens . As pointed out already in garcia , in the case of a bent vortex before ballistic expansion, an important contribution can arise due to the bending itself, as shown in Figure 4. In this picture we plot a cut of the column density along the $`y`$ axis where the vortex bends, for three values of the rotation frequency $`\mathrm{\Omega }`$. In Figure 5 we compare the column density for a bent vortex at $`\mathrm{\Omega }=0.59\omega _{}`$, before and after the release from the trap (time of flight $`30`$ms). The $`y`$ coordinate is given in rescaled units $`\stackrel{~}{y}y/\lambda _2(t)`$. From this figure we can see that the expansion has no relevant effects on the contrast; the same result holds also for the other configurations shown in Figure 2.

A second application is the analysis of the shape of the vortex line after ballistic expansion. This is also an important issue, as the shape of the vortex line can be measured using transverse absorption imaging rather than axial absorption imaging vu\_a\_ENS . Focusing on the configuration $`\mathrm{\Omega }=0.37\omega _{}`$ exhibiting the strongest bending, we have calculated the shape of the vortex line before and after expansion by looking in each horizontal plane for the local minimum of density closest to the vertical axis $`z`$. Initially the vortex line is contained in the half plane $`(y)z`$ which is at an angle $`\theta (0)=\pi /2`$ with respect to $`x`$ axis. After a time $`t_e30`$ms of ballistic expansion, we find that the vortex line remains almost planar, but in a vertical half plane at an angle $`\theta (t_e)1.14`$rad with respect to $`x`$ axis, and not passing through the origin. We show in Figure 6 the vortex line before and after expansion, projected in the vertical planes in which it is contained at $`t=0`$ and $`t=t_e`$. One then sees that the vortex line essentially preserves its initial shape after time of flight, in the frame of the rescaled spatial coordinates.

### III.5 Existence domain of a single vortex configuration

In the last part of this section we discuss the thermodynamic existence domain of a single vortex configuration, that is the domain where it is a local minimum of energy. By varying the radial frequency $`\omega _{}`$ we have explored a wide region of trap geometries, ranging from the spherical case to long cigar traps ($`\lambda ^1[1,15]`$).

In Figure 7 we show this domain in the $`\omega _{}/\omega _z\mathrm{\Omega }/\omega _{}`$ plane. Here we consider an almost axially symmetric trap, with a very small asymmetry ($`ϵ=10^4`$), which fixes the direction along which the vortex bends. By minimizing numerically the energy functional (1) with the conjugate gradient method we find that thermodynamically stable configurations with one vortex lie in the interval between the lines with empty and full circles, representing respectively $`\mathrm{\Omega }_1`$ and $`\mathrm{\Omega }_2`$. The black squares represent the predictions of the sum rules, which give a good estimate for $`\mathrm{\Omega }_2`$, as discussed in §III.3. In the same picture we also show the results from the Bogoliubov approach for a straight vortex, as described in §V. Up-triangles correspond to the critical rotation frequency which stabilize a straight vortex ($`\mathrm{\Omega }_1`$); this line separates the existence domains for straight and bent vortices. Down-triangles correspond to the frequency at which the straight vortex is no longer thermodynamically stable due to the destabilization of the surface modes ($`\mathrm{\Omega }_2`$). Notice that these values are in very good agreement with the result of the conjugate gradient in the whole range of existence of the straight vortex.

Finally, the solid line represent a 2D ansatz based on the analytical model in Ref. castindum , as described in the next section. If we imagine our 3D condensate as a collection of 2D slices, we can suppose that the 3D vortex is a stable configuration (i.e. the vortex line remains close to $`z`$ axis and does not get away) if the rotation frequency $`\mathrm{\Omega }`$ at least exceeds the 2D critical rotation frequency $`\mathrm{\Omega }_c^{2D}`$ at which a 2D vortex in the central slice ($`z=0`$) becomes energetically favorable with respect to the solution without vortices, that is: $`\mathrm{\Omega }>\mathrm{\Omega }_c^{2D}(z=0)`$. The expression for $`\mathrm{\Omega }_c^{2D}`$, which will be derived in the next section, is

$$\stackrel{~}{\mathrm{\Omega }}_c^{2D}=\frac{1}{\stackrel{~}{\mu }}\mathrm{ln}\left(C^{}\left(\stackrel{~}{\mu }+\eta \right)\right)$$
(25)
where the tilda indicates rescaling by $`\omega _{}`$, e.g. $`\stackrel{~}{\mu }\mu /\mathrm{}\omega _{}`$, $`C^{}1.8011`$, and $`\eta =1.9378`$. Then, by expressing $`\stackrel{~}{\mu }`$ in terms of $`\lambda =\omega _z/\omega _{}`$ in the Thomas-Fermi approximation:

$$\stackrel{~}{\mu }=\frac{1}{2}\left(15N\frac{a}{a_z}\right)^{2/5}\lambda ^{1/5}13.1\lambda ^{1/5}$$
(26)
we have

$$\stackrel{~}{\mathrm{\Omega }}_c^{2D}=\frac{b}{\lambda ^{1/5}}\mathrm{ln}\left(c\lambda ^{1/5}+d\right)$$
(27)
with $`b0.0763`$, $`c23.6`$ and $`d3.49`$. As we can see from Figure 7 this simple ansatz gives a good estimate of the critical frequency $`\mathrm{\Omega }_1`$. This ansatz will be justified in the next section.

## IV Understanding the bending of the vortex line analytically

The previous numerical results show that the vortex line in a steady state configuration is not necessarily straight when the condensate is cigar shaped along the rotation axis $`z`$, in accordance with previous numerical results based on a different algorithm garcia . This however does not explain physically why the vortex bends. To get the required physical understanding we derive an approximate energy functional for the vortex line in the Thomas-Fermi limit in the spirit of AftaRiv , and we minimize this energy functional with a simple variational ansatz. We reach a very simple prediction for the minimal rotation frequency required to stabilize a bent vortex, which is in good agreement with the full numerics when the condensate is cigar-shaped.

### IV.1 Deriving a simple energy functional

We restrict in what follows to the interesting regime of a cigar shaped condensate, where $`\omega _z\omega _{}`$. The first step is to transform the Gross-Pitaevskii energy functional (1) into a functional of the shape of the vortex line only. This assumes that both the modulus and the phase of the condensate wavefunction $`\varphi `$ can be expressed in terms of the shape of the vortex line only. This is in general a formidable task, as the condensate density is not uniform in a harmonic potential pas\_uni , but it is greatly simplified if we restrict to the Thomas-Fermi limit $`\mu \mathrm{}\omega _{}`$. We present a rather detailed and pedestrian derivation in the appendix A, we give here only the main ideas.

In the Thomas-Fermi regime there is a clear separation of spatial scales between the vortex core radius, of the order of the healing length $`\xi =(\mathrm{}^2/m\mu )^{1/2}`$, and the transverse Thomas-Fermi radius of the condensate $`R_{}=(2\mu /m\omega _{}^2)^{1/2}`$, where $`\mu `$ is the chemical potential of the gas. The total density can then be written as the product of a slowly varying envelope function and of a narrow ‘hole’ function defining the vortex core Fetter\_review ; castindum . We further assume that the rotation frequency $`\mathrm{\Omega }`$ is of the order of $`\mathrm{}\omega _{}^2/\mu `$. As a consequence the rotational velocity term $`\stackrel{}{\mathrm{\Omega }}\stackrel{}{r}`$, at most of the order of $`\mathrm{\Omega }R_{}`$, is much smaller than the typical velocity field in the lab frame at a distance $`\xi `$ from the vortex core, $`v\mathrm{}/m\xi `$, in the Thomas-Fermi limit, so that the structure of the vortex line is not distorted by rotation. Another consequence is that the envelope function is also not destabilized by the rotation personal\_estimate  and is close to the usual Thomas-Fermi expression.

Expressing the phase of the condensate as function of the vortex line is made difficult by the spatial inhomogeneity of the density profile of the condensate pas\_uni . In principle this phase has to be determined everywhere in the condensate if one wants to calculate the kinetic energy term of (1). Fortunately, using the fact that the condensate is in a steady state, one can replace the volume integral giving the kinetic energy stored in the condensate phase by a line integral along the vortex line, using the same type of techniques as in AftaRiv . It is then possible to rely on approximations for the condensate phase valid close to the vortex line. We use in the appendix A the simplifying hypothesis that the vortex line is weakly curved, with a radius of curvature of the order of $`R_{}`$ or larger, which allows to approximate the condensate phase close to the vortex line by the one of a straight vortex.

We finally obtain the following energy functional of the vortex line, taking the vortex free configuration as the zero of energy:

$`E_\mathrm{v}{\displaystyle 𝑑s}`$ $`{\displaystyle \frac{g_{2\mathrm{D}}(z_0(s))}{g}}[E_{2\mathrm{D}}^{\mathrm{\Omega }=0}(r_0(s);z_0(s))`$ (28)
$`+\mathrm{cos}(\alpha (s))E_{2\mathrm{D}}^{\mathrm{rot}}(r_0(s);z_0(s))].`$
In this expression the vortex line is parametrized by the curvilinear abscissa $`s`$. At the point of abscissa $`s`$ the vortex line is at the elevation $`z_0(s)`$ and at a distance $`r_0(s)`$ from the rotation axis, and makes an angle $`\alpha (s)`$ with respect to $`z`$. A remarkable feature of (28) is that it is expressed in terms of the energy functionals of a vortex core in a 2D condensate, $`E_{2\mathrm{D}}^{\mathrm{\Omega }=0}`$ for the energy in the absence of rotation and $`E_{2\mathrm{D}}^{\mathrm{rot}}`$ for the energy due to the $`\mathrm{\Omega }L_z`$ term. This is physically plausible considering the cigar shaped nature of the condensate, and this allows to view the 3D condensate as a collection of 2D horizontal slices. The slice of elevation $`z`$ constitutes a 2D Bose condensate with a chemical potential $`\mu m\omega _z^2z^2/2`$, where $`\mu `$ is the chemical potential of the 3D condensate, and has a Thomas-Fermi radius coinciding with the local 3D one. The coupling constant $`g_{2\mathrm{D}}(z)`$ of the 2D gas can be expressed in terms of the 3D coupling constant, see (115). We arrive at the simple formula

$$\frac{g_{2\mathrm{D}}(z)}{g}=\frac{15}{16R_z}\left(1\frac{z^2}{R_z^2}\right)^2$$
(29)
where $`R_z`$ is the Thomas-Fermi radius of the condensate along $`z`$. An interesting remark is that the rotational energy term in (28) is multiplied by $`\mathrm{cos}\alpha (s)`$. As $`E_{2\mathrm{D}}^{\mathrm{rot}}`$ is proportional to the rotation frequency $`\mathrm{\Omega }`$, this means that $`\mathrm{cos}\alpha (s)E_{2\mathrm{D}}^{\mathrm{rot}}`$ is the rotational energy of a vortex core in a 2D condensate rotating at the effective frequency

$$\mathrm{\Omega }_{2\mathrm{D}}(s)=\mathrm{\Omega }\mathrm{cos}\alpha (s).$$
(30)
### IV.2 Minimizing numerically the simple energy functional

In a first stage we have to check that the energy functional derived in the appendix A correctly reproduces the results of the minimization of the full Gross-Pitaevskii energy functional. We perform this check numerically: we discretize the vortex line in little segments having all the same length $`dl`$ much smaller than the transverse Thomas-Fermi radius $`R_{}`$ of the condensate. As the vortex line is symmetric with respect to $`z`$ reflexion and lies in the $`xz`$ plane, the left extremity of the first segment in the calculation moves along $`x`$ axis only, with an abscissa $`x_0`$. The $`z>0`$ part of the vortex line is discretized in $`k`$ segments and its shape is parametrized by the angles $`\alpha _i`$, $`i=1,\mathrm{},k`$ at which the $`k`$ segments are with respect to the axis $`z`$. The energy functional (28) in its discretized version is then a function of $`k+1`$ coordinates, that is of $`x_0`$ and of the $`k`$ angles $`\alpha _i`$. Starting with a straight vortex line at some small angle with respect to $`z`$ axis we move $`x_0`$ and the $`\alpha _i`$’s according to the simple gradient method or imaginary time evolution method, that is we move the parameters by a small step in the direction opposite to the local gradient of the energy functional.

It is known that this simple gradient method is not efficient when the desired minima are at the bottom of a very elongated valley in the coordinate space numrec . This potential problem is minimized by a rescaling of the coordinates by their natural units $`x_0^{\mathrm{typ}}=R_{}`$ and $`\alpha ^{\mathrm{typ}}=2\pi `$, so that our specific algorithm is to iterate the following small coordinate changes

$`dx_0`$ $`=`$ $`d\tau \left(x_0^{\mathrm{typ}}\right)^2_{x_0}E_\mathrm{v}/E_v^{\mathrm{typ}}`$ (31)
$`d\alpha _i`$ $`=`$ $`d\tau \left(\alpha ^{\mathrm{typ}}\right)^2_{\alpha _i}E_\mathrm{v}/E_v^{\mathrm{typ}}`$ (32)
where the typical energy scale is $`E_v^{\mathrm{typ}}=\mathrm{}^2\omega _{}^2/\mu `$ and the dimensionless ‘time’ step $`d\tau `$ is small as compared to unity (typically $`0.1`$). The iteration stops when $`|dx_0|/(d\tau x_0^{\mathrm{typ}})`$ and $`|d\alpha _i|/(d\tau \alpha ^{\mathrm{typ}})`$ are below some threshold, taken here to be $`10^6`$.

The resulting prediction for the minimum rotation frequency required to stabilize the vortex line is shown in Figure 8, deliberately restricted to the domain of cigar shaped condensates. The agreement with the minimization of the full Gross-Pitaevskii energy functional is remarkable, considering the fact that the points of Figure 8 are moderately in the Thomas-Fermi regime ($`\mu 7.6\mathrm{}\omega _{}`$ for $`\omega _{}/\omega _z=15`$).

### IV.3 Why the vortex line bends in a cigar shaped condensate

The actual goal of this section is to understand physically why the vortex line bends in a cigar shaped condensate. This can be achieved intuitively thanks to the very suggestive form of the energy functional (28). One just needs to have in mind the following characteristics of the 2D vortex problem:

* if the effective rotation frequency $`\mathrm{\Omega }_{2\mathrm{D}}`$ is too small the 2D energy functional has a maximum for the vortex core at the center of the trap and is a purely repulsive potential, see the dotted line in figure 9: the vortex core cannot be stabilized inside the condensate and its equilibrium position is at infinity.
* if $`\mathrm{\Omega }_{2\mathrm{D}}`$ is above the stabilization frequency $`\mathrm{\Omega }_{\mathrm{stab}}^{2\mathrm{D}}`$ and below the critical rotation frequency $`\mathrm{\Omega }_\mathrm{c}^{2\mathrm{D}}`$, the 2D energy functional has a local but not global minimum for the vortex core at the center of the trap (dashed line in Figure 9). In this situation, the vortex core is stabilized at the trap center with an energy larger than the vortex free condensate.
* for $`\mathrm{\Omega }_{2\mathrm{D}}>\mathrm{\Omega }_\mathrm{c}^{2\mathrm{D}}`$, the energy minimum at the trap center is now below the energy of the vortex free configuration (solid line in Figure 9).

The important feature of the 2D case is that the equilibrium positions of the vortex core are either the trap center or the infinity. Another point, crucial for the 3D case, is that both $`\mathrm{\Omega }_{\mathrm{stab}}^{2\mathrm{D}}`$ and $`\mathrm{\Omega }_\mathrm{c}^{2\mathrm{D}}`$ are decreasing functions of the chemical potential $`\mu _{2\mathrm{D}}`$.

Let us now follow the vortex line travelling through the 3D cigar shaped condensate, in the case where the 2D critical rotation frequency in the central slice $`z=0`$ is smaller than the actual rotation frequency $`\mathrm{\Omega }`$. Let us call $`z_\mathrm{c}`$ the elevation of the 2D slice with a local critical frequency $`\mathrm{\Omega }_\mathrm{c}^{2\mathrm{D}}`$ equal to $`\mathrm{\Omega }`$, see Figure 10a.

It is clear that the vortex line will be straight along the rotation axis for $`|z|<z_\mathrm{c}`$: the vortex line is there in a valley corresponding to the global minimum of energy of each local 2D slice. When the vortex line reaches the domain of elevation $`z>z_\mathrm{c}`$, having a vortex core on the rotation axis costs more energy than having the vortex core at infinity in each local 2D slice. The tempting strategy then offered to the vortex line is to bend and leave the condensate. Assume that the vortex line leaves the condensate radially, as shown in Figure 10a. The corresponding horizontal vortex line has an energy

$$E_\mathrm{v}^{\mathrm{horiz}}=\frac{g_{2\mathrm{D}}(z_c)}{g}_0^+\mathrm{}𝑑xE_{2\mathrm{D}}^{\mathrm{\Omega }=0}(x;z_c).$$
(33)
The corresponding integral can be calculated exactly, but it is sufficient to give here an order of magnitude: $`g_{2\mathrm{D}}/g`$ scales as $`1/R_z`$, $`E_{2\mathrm{D}}^{\mathrm{\Omega }=0}`$ is of the order of $`\mathrm{}^2\omega _{}^2/\mu _{2\mathrm{D}}`$ and the integral over $`x`$ converges over a distance given by the Thomas-Fermi radius $`R_{}`$, so that

$$E_\mathrm{v}^{\mathrm{horiz}}\frac{R_{}}{R_z}\frac{(\mathrm{}\omega _{})^2}{\mu }\left(1\frac{z_c}{R_z}\right)^{3/2}.$$
(34)
We have included some approximate $`z_c`$ dependence relevant for the extreme case of $`z_c`$ close to $`R_z`$. What would be the energy cost for the vortex line to remain on the rotation axis from $`z_c`$ to $`R_z`$? The energy of the corresponding vertical segment is

$`E_\mathrm{v}^{\mathrm{vert}}`$ $`=`$ $`{\displaystyle _{z_c}^{R_z}}𝑑z{\displaystyle \frac{g_{2\mathrm{D}}(z)}{g}}\left[E_{2\mathrm{D}}^{\mathrm{\Omega }=0}(0;z)+E_{2\mathrm{D}}^{\mathrm{rot}}(0;z)\right]`$
$`=`$ $`{\displaystyle \frac{15}{16R_z}}{\displaystyle _{z_c}^{R_z}}𝑑z\left(1{\displaystyle \frac{z^2}{R_z^2}}\right)^2\left[\mathrm{}\mathrm{\Omega }_\mathrm{c}^{2D}(z)\mathrm{}\mathrm{\Omega }_\mathrm{c}^{2D}(z_c)\right]`$
where we have used the fact that the 2D rotational energy of a vortex core in the center of a trap rotating at frequency $`\mathrm{\Omega }`$ is $`\mathrm{}\mathrm{\Omega }`$, also equal to $`\mathrm{}\mathrm{\Omega }_\mathrm{c}^{2D}(z_c)`$ by definition of $`z_c`$. The integrand has the same order of magnitude as in $`E_\mathrm{v}^{\mathrm{horiz}}`$ but the integration length is now of the order of $`R_z`$ so that

$$E_\mathrm{v}^{\mathrm{vert}}\frac{\mathrm{}^2\omega _{}^2}{\mu }\left(1\frac{z_c}{R_z}\right)^2$$
(36)
is typically $`R_z/R_{}`$ times larger than $`E_\mathrm{v}^{\mathrm{horiz}}`$. This proves that in the limit of a cigar-shaped condensate $`R_zR_{}`$, the strategy of bending is more favorable than the strategy of following the rotation axis, except for a $`z_c`$ very close to the end point of the condensate, $`1z_c/R_zR_{}^2/R_z`$ that is $`\mathrm{\Omega }(\mathrm{}^2\omega _{}^2/\mu )(R_z/R_{})^2`$.

The bending at the value $`z_c`$ in the above reasoning can be justified variationally as follows. We perform a simple minimization of the energy functional (28) with the following linear piecewise variational ansatz for the vortex line in the plane $`y=0`$: an horizontal line linking $`(x=+\mathrm{},z=Z)`$ to $`(x=0,z=Z)`$, then a vertical segment of length $`2Z`$ along the rotation axis, then an horizontal line from $`(x=0,z=Z)`$ to $`(x=+\mathrm{},z=Z)`$. The energy of the ansatz depends on the single variational parameter $`Z`$:

$$E_\mathrm{v}^{\mathrm{ans}}(Z)=E_\mathrm{v}^{\mathrm{segm}}(Z)+E_\mathrm{v}^{\mathrm{line}}(Z)$$
(37)
that is the sum of the energies of the vertical segment and of the horizontal lines. One has to extremize this function over $`Z`$. In the case where the extremum is in the interior of the interval $`(0,R_z)`$, one has to solve

$$\frac{d}{dZ}E_\mathrm{v}^{\mathrm{ans}}=0.$$
(38)
This non-trivial task becomes simple in the limit of a very elongated condensate along $`z`$ proviso . As shown in appendix A, the vortex 2D energy function $`E_{2\mathrm{D}}^{\mathrm{\Omega }=0}(r_{};z)`$ depends on $`r_{}`$ only via the ratio $`r_{}/R_{}(z)`$ where $`R_{}(z)`$ is the local Thomas-Fermi radius. As a consequence $`E_\mathrm{v}^{\mathrm{line}}(Z)`$ depends on $`Z`$ through $`Z/R_z`$ only, and its derivative is $`R_{}/R_z`$ times smaller than the derivative of $`E_\mathrm{v}^{\mathrm{segm}}(Z)`$. In the limit $`R_z/R_{}`$ tending to infinity, equation (38) reduces to

$$\frac{d}{dZ}E_\mathrm{v}^{\mathrm{segm}}(Z)=0.$$
(39)
Using the explicit expression of $`E_\mathrm{v}^{\mathrm{segm}}(Z)`$ similar to (IV.3) we obtain the condition

$$\mathrm{\Omega }_\mathrm{c}^{2\mathrm{D}}(Z)=\mathrm{\Omega }$$
(40)
that is $`Z=z_\mathrm{c}`$, the vortex starts bending at the elevation where the 2D critical rotation frequency equals the trap rotation frequency $`\mathrm{\Omega }`$. As a consequence the minimal rotation frequency to stabilize the bent vortex line is given in the $`\omega _z/\omega _{}0`$ limit by the 2D critical frequency in the central slice $`z=0`$, see (125):

$$\mathrm{\Omega }_1=\frac{\mathrm{}\omega _{}^2}{\mu }\mathrm{log}\left[e^{C+1/2}\left(\frac{\mu }{\mathrm{}\omega _{}}+\eta \right)\right]$$
(41)
with $`\eta 1.938`$ and $`C0.0884`$. This asymptotic prediction is plotted as a solid line in Figure 8 and is in good agreement both with the numerical minimization of the vortex energy functional and with the numerical minimization of the full Gross-Pitaevskii energy functional general .

To conclude this subsection we point out two striking properties of the result (41). Firstly, it explains why in a cigar-shaped condensate, the bent vortex is first stabilized with an almost vanishing energy gap with respect to the 0-vortex configuration, at least much smaller than $`\mathrm{}^2\omega _{}^2/\mu `$, see Figure 1: at $`\mathrm{\Omega }=\mathrm{\Omega }_1`$, the variational ansatz predicts an energy scaling as $`(\mathrm{}^2\omega _{}^2/\mu )R_{}/R_z`$, see the expression for $`E_v^{\mathrm{horiz}}`$ with $`z_c/R_z0`$. This is very different from the 2D case, where the vortex, when first stabilized, has a large and positive energy $`\mathrm{}^2\omega _{}^2/\mu `$ with respect to the vortex free configuration. Secondly, it is remarkable that the 2D stabilization frequency plays no role in the 3D case for the cigar-shaped condensate. This means that being a local minimum of energy in 2D slices does not imply that the vortex line is a local minimum of energy in 3D. As an example, it is possible in 3D to shorten the vertical segment of the vortex line, that is to introduce the bending at a slightly lower elevation, whereas this infinitesimal transformation has no equivalent in a purely 2D case.

### IV.4 Saddle points of the simple energy functional

Using a simple energy functional rather than the full Gross-Pitaevskii energy functional has allowed to understand the bending of the vortex line physically with the help of a simple piecewise variational ansatz. But much more can be done: as we now show, one can investigate not only the energy minimum but also possible saddle points of the energy functional. This will explain an intriguing feature of the numerical results of Figure 1: the angular momentum of the condensate has a discontinuous jump from the one-vortex branch to the no-vortex branch, which opens a gap in the allowed values of $`L_z`$.

The idea is to consider now a more general ansatz for the vortex line than the one of the previous subsection, in order to allow the vortex line to move towards the border of the condensate. We introduce the two parameter linear piecewise ansatz in the $`y=0`$ plane shown in Figure 10b: a horizontal line linking $`(x=+\mathrm{},z=Z)`$ to $`(x=X,z=Z)`$, then a vertical segment of length $`2Z`$ parallel to $`z`$, and finally a horizontal line from $`(x=X,z=Z)`$ to $`(x=+\mathrm{},z=Z)`$. In other words, the vortex line has a C shape of adjustable distance $`X`$ to the rotation axis and adjustable height $`2Z`$.

We consider the parameters of the ENS experiment, with $`\lambda =\omega _z/\omega _{}=1/15`$. Eq.(26) leads to a Thomas-Fermi chemical potential of $`\mu =7.62\mathrm{}\omega _{}`$. We choose a rotation frequency $`\mathrm{\Omega }=0.448\omega _{}`$, which is a factor 1.2 above the threshold value $`\mathrm{\Omega }_1`$ of (41). Figure 11 shows a contour plot of the simple energy functional as function of the coordinates $`(X,Z)`$ of the bending point of the vortex line. One recognizes first the global minimum of energy, marked with a M, with a position of the C shaped vortex line similar to the one of Figure 10a: the vertical part of the vortex line is very close to the rotation axis, with $`X`$ less than 0.05 times the transverse Thomas-Fermi radius $`R_{}`$. Second, one finds a saddle point in the energy, marked with a S, corresponding to a quite different position of the C shape: the vortex line is now far from the rotation axis, with $`X0.45R_{}`$, and the half elevation of the $`C`$ shape is of the order of $`R_{}`$.

What is the physical meaning of this saddle point ? As it corresponds to an extremum of the energy functional (vanishing first order derivatives of the functional), it represents a stationary shape for the vortex line in the rotating frame. In addition to this dynamic property, it has the following interesting energetic aspect: it corresponds to an energy minimum for a fixed value of the angular momentum per particle $`L_z`$, as one can check explicitly for the two-parameter ansatz.

These properties are not specific to the two-parameter ansatz. They apply for arbitrary shapes of the vortex line, but also for the exact Gross-Pitaevskii energy functional. Consider indeed the subspace of condensate wavefunctions with a fixed angular momentum per particle $`l_z`$. The Gross-Pitaevskii energy functional in the absence of rotational term $`\mathrm{\Omega }L_z`$ is bounded from below (for $`g>0`$) in this subspace so it has a minimum for some wavefunction $`\psi _{l_z}`$, with an energy $`E^{\mathrm{\Omega }=0}(l_z)`$. One can then show that the wavefunction $`\psi _{l_z}`$ is a stationary point for the full Gross-Pitaevskii energy functional, that is including the rotational energy term $`\mathrm{\Omega }L_z`$, without imposing a fixed the angular momentum, provided that one takes the rotation frequency equal to $`\mathrm{\Omega }=dE^{\mathrm{\Omega }=0}/dl_z`$. One can show that it is however not necessarily a minimum of the full energy functional: it is a minimum if

$$\frac{d\mathrm{\Omega }}{dl_z}\frac{d^2E^{\mathrm{\Omega }=0}}{dl_z^2}>0$$
(42)
otherwise it is a saddle point.

We can exemplify these properties with our simple two-parameter ansatz. We have plotted in Figure 12 the mean angular momentum per particle of the energy minimum and of the saddle point as function of the rotation frequency $`\mathrm{\Omega }`$. One finds that the branch of minimum has an angular momentum ranging from 0.42$`\mathrm{}`$ to $`\mathrm{}`$, whereas the branch of saddle point fills the gap of angular momentum from 0 to 0.42$`\mathrm{}`$. One also finds that the angular momentum is an increasing function of the rotation frequency on the branch of minimum, whereas it is a decreasing function of $`\mathrm{\Omega }`$ on the saddle point branch, in accordance with the general condition (42).

In an experiment where the angular momentum of the condensate is gradually decreased, starting from a value close to $`\mathrm{}`$, the vortex first follows the minimum energy branch: the distance of the vortex line to the rotation axis $`X`$ remains much smaller than $`R_{}`$, whereas the vortex half-height $`Z`$ and the rotation frequency $`\mathrm{\Omega }`$ of the line in the lab frame decrease with time. Then the angular momentum reaches a critical value: the minimum in the $`(X,Z)`$ parameter space merges with the saddle point. When the angular momentum is further decreased, the vortex line follows the saddle point branch: the distance $`X`$ increases and becomes of the order of $`R_{}`$, as the vortex moves towards the boundary of the condensate; $`Z`$ takes values as small as $`R_{}`$, in which case the validity of our simple analytical approach becomes marginal, and the rotation frequency of the vortex line increases. This scenario was very recently observed at ENS vu\_a\_ENS .

To conclude this subsection, we can emphasize again the analogy of the 3D bent vortex problem and the 2D single vortex problem. In 2D, the minimum of energy in the rotating frame with no constrainst on the angular momentum corresponds to the vortex at the trap center, with an angular momentum $`\mathrm{}`$, or to the vortex at infinity with a vanishing angular momentum. If one minimizes the energy for a fixed angular momentum $`l_z`$, one finds an off-center vortex core, shifting from the trap center to infinity as $`l_z`$ is varied from $`\mathrm{}`$ downto zero; these off-center solutions do not satisfy the criterion (42) Rokhsar  and are very analogous to our saddle point solutions.

## V Finite temperature fluctuations of an otherwise straight vortex line

### V.1 Method: Bogoliubov approach around the straight vortex steady state

In order to evaluate the contribution of thermal fluctuations on the contrast, we study the limiting case of a straight vortex. For the ENS trap parameters, this means that $`1/\lambda <7`$ and in all the following, we illustrate our discussion using the case $`(\lambda =1/5,\mathrm{\Omega }=0.5\omega _{})`$. In this regime, the system has a rotational symmetry around the $`z`$ axis so that the numerical problem is effectively 2D in cylindrical coordinates. We use the $`U(1)`$-symmetry preserving Bogoliubov approach described in U1 . In this way, the problem of the spurious mode of the condensate is avoided. The Bose field is expanded as

$$\widehat{\psi }(\stackrel{}{r})=\varphi (\stackrel{}{r})\widehat{a}_\varphi +\widehat{A}_\varphi \underset{k}{}\widehat{b}_ku_k(\stackrel{}{r})+\widehat{b}_k^{}v_k^{}(\stackrel{}{r}).$$
(43)
$`\varphi `$ is the condensate wave function normalized to unity. The operator $`\widehat{a}_\varphi `$ is the annihilation operator in the condensate mode and the almost unitary operator $`\widehat{A}_\varphi =(\widehat{a}_\varphi ^{}\widehat{a}_\varphi +1)^{1/2}\widehat{a}_\varphi `$ gives the phase factor of the field in the condensate mode Girardeau ; Wigner . The modal functions $`u_k,v_k`$ are normalized like $`d^3\stackrel{}{r}|u_k|^2|v_k|^2=1`$. They are obtained from the usual modal functions of the Bogoliubov-de Gennes equations after orthogonalization of $`u_k`$ and of $`v_k`$ with respect to $`\varphi `$. The Bose operators $`\widehat{b}_k`$ annihilate one quasi-particle in the mode $`k`$ but conserve the total number of particles. The index $`k=\{n,l,s\}`$ denotes the quantum numbers of the mode linked to the symmetry of the system: $`s=1`$ (or $`1`$) for symmetric (or antisymmetric) modes with respect to the plane $`z=0`$, $`l\mathrm{}`$ is the angular momentum with respect to the condensate and the integer $`n`$ is the radial quantum number.

Concerning the transverse direction $`xy`$, we expand $`\varphi ,u_k,v_k`$ on the harmonic oscillator basis $`\left\{\mathrm{\Phi }_{\mathrm{ho}}^{m^{},n^{}}\right\}`$ of pulsation $`\omega _{}`$ ($`m^{}\mathrm{}`$ is the angular momentum and $`n^{}`$, the radial quantum number). For example, in this basis we have

$$\varphi (\stackrel{}{r})=\underset{n^{}=0}{\overset{\mathrm{}}{}}c_n^{}(z)\mathrm{\Phi }_{\mathrm{ho}}^{m=1,n^{}}(x,y).$$
(44)
Numerically, the basis is truncated: the spatial grid along $`z`$ is surrounded by infinite walls, and also the number of wave functions $`\mathrm{\Phi }_{\mathrm{ho}}^{m^{},n^{}}`$ is limited (for the value $`\lambda =1/5`$, the harmonic oscillator basis in our computation contains all the wave functions of energy less than $`42`$ $`\mathrm{}\omega _{}`$). The choice of the grid and the discretization of the Laplacian along $`z`$ has been made so that the first 160 energy levels $`ϵ_n`$ of the pure 1D harmonic oscillator are recovered with an error $`|ϵ_nϵ_n^{\mathrm{exact}}|<10^2`$. In this part we have computed the condensate wave function using an imaginary time method and the convergence criterion (6) (typical values of $`\gamma `$ are of the order of $`10^{13}`$).

### V.2 Expectation values of some observables

The mean density of atoms out of the condensate is obtained straightforwardly with the usual expression

$$\rho _{\mathrm{exc}}(\stackrel{}{r})=\underset{k}{}|v_k|^2(\stackrel{}{r})+\underset{k}{}n_k\left(|u_k|^2(\stackrel{}{r})+|v_k|^2(\stackrel{}{r})\right),$$
(45)
where

$$n_k=1/(\mathrm{exp}(ϵ_k/k_BT)1)$$
(46)
is the Bose occupation factor and $`ϵ_k`$ is the energy of the Bogoliubov mode $`k`$ in the rotating frame. Because of numerical constraints we have limited the sum over states of energy less than $`20`$ $`\mathrm{}\omega _{}`$ in the rotating frame, and we have studied configurations with temperatures less than $`11.4`$ $`\mathrm{}\omega _{}/k_B`$ that is less than $`0.4T_c`$, where $`k_BT_c=\mathrm{}\left(\omega _z\omega _{}^2N/\zeta (3)\right)^{1/3}`$ is the ideal Bose gas critical temperature. For comparison with experimental results, we extract the mean column density $`n`$ in the $`xy`$ plane:

$$n(x,y)=𝑑z\left[(NN_{\mathrm{exc}})|\varphi |^2(\stackrel{}{r})+\rho _{\mathrm{exc}}(\stackrel{}{r})\right],$$
(47)
where $`N_{\mathrm{exc}}=d^3r\rho _{\mathrm{exc}}(\stackrel{}{r})`$ is the mean number of atoms out of the condensate. We plot this column density as function of the distance to the $`z`$ axis for different values of the temperature in Figure 13.

At zero temperature the effect of quantum depletion at the center of the vortex line is not observable on the scale of the figure. Hence quantum fluctuations of the vortex line are clearly not responsible for the low contrast measured in ENS experiments.

To be more quantitative, we denote $`n^{\mathrm{min}}`$ the value of the mean column density at the center of the trap and $`n^{\mathrm{max}}`$ the maximal value of the mean column density. Then, the mean contrast is defined as

$$\overline{𝒞}=1\frac{n^{\mathrm{min}}}{n^{\mathrm{max}}}.$$
(48)
The modal functions are of the form

$`u_k(\stackrel{}{r})=U_k(r,z)\mathrm{exp}\left[i(l+1)\theta \right]`$ (49)
$`v_k(\stackrel{}{r})=V_k(r,z)\mathrm{exp}\left[i(l1)\theta \right],`$ (50)
If $`l+10`$ the function $`U_k`$ vanishes in $`r=0`$ because of the centrifugal barrier. For the same reason $`V_k(r=0,z)`$ vanishes for $`l10`$. As a consequence only modes with $`l=\pm 1`$ contribute to $`n^{\mathrm{min}}`$. Furthermore, in the rotating frame the energy of the quasi-particles is given by Eq.(9) and the most populated modes are the lowest energy kelvons, characterized by $`l=1`$, which have a negative energy in the absence of rotation Fetter\_dynamics ; feder . In the example presented in this section, the lowest kelvon mode localized at the extremities of the condensate has an energy in the rotating frame given by $`E0.018\mathrm{}\omega _{}`$. Figure 14 shows that the dependence of $`\overline{𝒞}`$ as a function of temperature is almost linear. This behavior can be understood from the fact that the lowest energy modes have a semi-classical character ($`ϵ_kk_BT`$), with an occupation number depending linearly on temperature $`n_kk_BT/ϵ_k`$.

Figure 14 shows that temperature contributes significantly to the observed contrasts. In particular, for a temperature of the order of the chemical potential, which is a typical situation in experiments, we find that the contrast is 65%, not far from the contrast value of $``$ 50% observed at ENS. Note that some zero temperature bending of the vortex line is expected to occur for the parameters of ENS, see Figure 4, which adds up to the thermal fluctuations and further reduce the contrast in the experiment.

### V.3 Mimicking a single experimental run: a Glauber $`P`$ method

Let us insist now on the fact that expectations values calculated in the previous section are not sufficient for a quantitative understanding of experiments. They just give an order of magnitude of what is observed. Indeed, for example, for the particular value of $`\lambda `$ that we have chosen, the straight vortex is at the edge of the existence domain and one expects large fluctuations of the phase and of the density due to the soft core mode and also to the emergence of low energy surface modes. As a consequence, observables obtained from a single measurement of the system may differ notably from their mean value. Hence, the aim of this section is to evaluate the fluctuations induced by the soft modes. Calculations of the previous subsection have shown that quantum fluctuations do not contribute to the contrast at finite temperature, so that semi-classical field approaches are good tools to answer this problem. We use for that purpose the Glauber $`P`$ method of quantum optics gardiner .

This approach is made simple by the fact that the $`N`$-body density operator $`\widehat{D}`$ is the exponential of $`H_{\mathrm{Bog}}/(k_BT)`$, where the Bogoliubov Hamiltonian is a sum of decoupled harmonic oscillator terms, $`H_{\mathrm{Bog}}=E_0+_kϵ_k\widehat{b}_k^{}\widehat{b}_k`$. We introduce the coherent state of quasi-particles $`|\{\beta _k\}=|\beta _{k_1},\beta _{k_2},\beta _{k_3}\mathrm{}`$, where $`|\beta _k`$ is the normalized eigenstate of the annihilation operators $`\widehat{b}_k`$ with the eigenvalue $`\beta _k`$. Then one defines the Glauber $`P`$ representation of $`\widehat{D}`$:

$$\widehat{D}=\underset{k}{}\frac{d\text{Re}\beta _kd\text{Im}\beta _k}{\pi }P(\{\beta _k^{}\})|\{\beta _k^{}\}\{\beta _k^{}\}|.$$
(51)
The explicit calculation of $`P`$ is possible, as $`D`$ is Gaussian in the $`b_k`$’s gardiner , this gives the following Gaussian distribution:

$$P(\{\beta _k^{}\})=\underset{k}{}\frac{1}{n_k}\mathrm{exp}\left(\frac{\beta _k^{}\beta _k}{n_k}\right),$$
(52)
where $`n_k`$ is the mean number of quasi-particles in the mode $`k`$ given in Eq. (46). This Gaussian distribution is easily sampled. As $`P`$ is positive one can imagine that a given experimental realization of the gas is in the state $`|\{\beta _k\}`$ where the complex numbers $`\beta _k`$ vary randomly from one realization to the other. In the following, we want to determine the density and the velocity field of the atomic gas for a given Monte Carlo realization of the $`\{\beta _k\}`$.

The $`N`$-body distribution function corresponding to a single term of the statistical mixture (51)

$`\rho (\stackrel{}{r}_1,\mathrm{},\stackrel{}{r}_N)`$ $`=`$ $`\{\beta _k^{}\}|\widehat{\psi }^{}(\stackrel{}{r}_1)\mathrm{}\widehat{\psi }^{}(\stackrel{}{r}_N)`$ (53)
$`\times \widehat{\psi }(\stackrel{}{r}_N)\mathrm{}\widehat{\psi }(\stackrel{}{r}_1)|\{\beta _k^{}\}`$
is not easy to calculate for the interacting Bose gas as $`\widehat{\psi }`$ is a superposition of $`\widehat{b}_k`$ and $`\widehat{b}_k^{}`$ so that the product of field operators in (53) is not normally ordered in terms of the $`\widehat{b}_k`$. We perform the following approximation

$$\widehat{b}_k^{}|\{\beta _k^{}\}\beta _k^{}|\{\beta _k^{}\}.$$
(54)
The error has a root mean square norm equal to unity, hence the approximation is good for modes with a large occupation number, bad for empty modes. In this way we describe correctly the fluctuations due to finite temperature, but not the quantum fluctuations existing even at zero temperature. As we have shown previously, this is fine in the dilute limit $`(\rho a^3)^{1/2}1`$, where quantum depletion is small. The approximation amounts to taking

$$\widehat{\psi }(\stackrel{}{r})|\{\beta _k^{}\}\psi (\stackrel{}{r})|\{\beta _k^{}\}$$
(55)
with a classical field

$$\psi (\stackrel{}{r})=\sqrt{N_0}\varphi (\stackrel{}{r})+\underset{k}{}\beta _ku_k(\stackrel{}{r})+\beta _k^{}v_k^{}(\stackrel{}{r}),$$
(56)
where $`N_0`$ is such that the norm squared of $`\psi `$ is equal to the total number of particles $`N`$ phase\_globale . In this approximation, a single realization is now in the coherent state $`|\{\psi \}`$ for the field $`\widehat{\psi }`$. In this case the $`N`$-body distribution function of the atoms for a single experimental realization is factorized:

$$\rho (\stackrel{}{r}_1,\mathrm{},\stackrel{}{r}_N)=\underset{k}{}|\psi (\stackrel{}{r}_k)|^2.$$
(57)
In this paragraph, we describe temperature effects in a single stochastic realization of $`\psi `$ at $`k_BT=0.2k_BT_C\frac{1}{2}\mu `$. First, we have extracted the shape of the vortex line: Figure 15 shows the $`x`$ and $`y`$ coordinates of the vortex line as function of $`z`$.

The soft core modes localized at the two extremities Fetter\_dynamics ; feder  of the condensate are thermally excited and as a result push the core away from the rotation axis. This is a temperature induced bending. Second, for the same stochastic realization, Figure 16 represents the location of the vortices in the velocity field of the gas in the plane $`z=0`$. It shows that satellite vortices, which are not easily observable in the density profile, appear at the border of the condensate. The presence of these vortices is due to the excitation of surface modes of low energy and relatively high angular momentum (recall that for $`\lambda =1/5,\mathrm{\Omega }=0.5\omega _{}`$, the system is at the border of the thermodynamic stability domain).

We address now the problem of the integrated density contrast for a single stochastic realization. In the previous single stochastic realization at a temperature $`T=0.2T_c`$, the value of the contrast is around $`90\%`$ which does not coincide with the mean contrast $`\overline{𝒞}85\%`$ derived in the previous subsection. This disagreement suggests strong fluctuations of density along the $`z`$axis. To confirm this idea we have computed 10000 stochastic realizations of the column density at the center of the trap. We have reported in Figure 17 the corresponding histogram.

This figure shows effectively that the probability law is far from being Gaussian: the long tail on the right side is due to 1D character of the excitation modes of the vortex line, the so-called kelvons, similar indeed to the fluctuations of the number of condensate particles in a 1D Bose gas cartago . This is another indication of the importance of fluctuations in the properties of this system.

To be complete we now consider the effect on the contrast of the ballistic expansion of the cloud performed in the experiment. As in subsection III.4 we integrate numerically the rescaled Gross-Pitaevskii equation, with a classical field $`\psi `$ including thermal fluctuations. As we see in Figure 18 the rescaling now absorbs to a lower extent the effect of the ballistic expansion, but the density contrast is only weakly changed.

## VI Conclusion

We have understood why the vortex line can bend in a steady state cigar-shaped condensate rotated at frequency $`\mathrm{\Omega }`$ around its long axis $`z`$. In the Thomas-Fermi regime the cigar-shaped condensate can be viewed as a collection of slices parallel to the $`xy`$ plane, each slice corresponding formally to a 2D rotating condensate. For each 2D condensate one defines as usual the critical rotation frequency $`\mathrm{\Omega }_\mathrm{c}^{2\mathrm{D}}`$ above which it is energetically more favorable to have the vortex core at the trap center rather than at infinity. As the 3D condensate density is inhomogeneous along $`z`$ axis due to the harmonic confinement the local $`\mathrm{\Omega }_\mathrm{c}^{2\mathrm{D}}`$ is minimal in the plane $`z=0`$ and maximal at the end points of the cigar.

The vortex line then uses the following strategy to minimize its energy: it follows the rotation axis $`z`$ in the elevation interval where the trap rotation frequency is larger than the local $`\mathrm{\Omega }_\mathrm{c}^{2\mathrm{D}}`$’s, and moves away from the rotation axis to infinity where $`\mathrm{\Omega }`$ becomes smaller than the local $`\mathrm{\Omega }_\mathrm{c}^{2\mathrm{D}}`$.

This leads to the analytical prediction that the minimal rotation frequency $`\mathrm{\Omega }_1`$ required to stabilize the bent vortex in a cigar-shaped condensate is equal to $`\mathrm{\Omega }_\mathrm{c}^{2\mathrm{D}}(z=0)`$, that is the 2D critical rotation frequency corresponding to the central slice $`z=0`$.

Another analytical prediction is that, for $`\mathrm{\Omega }>\mathrm{\Omega }_1`$, the energy presents a saddle point, in addition to the previously discussed minimum. This saddle point is a stationary configuration of the vortex line in the rotating frame, with an angular momentum smaller than the one of the energy minimum. The vortex line in the saddle point has also a different position with respect to $`z`$ axis than in the minimum energy configuration: the vortex line, rather than being on axis, is shifted away from the rotation axis by an amount of the order of the transverse Thomas-Fermi radius $`R_{}`$, and it has a much smaller length along $`z`$, also on the order of $`R_{}`$. This description matches very recent experimental observations of bent vortices at ENS vu\_a\_ENS .

To test the analytical predictions for the energy minimum we have performed a full numerical minimization of the Gross-Pitaevskii energy functional with the efficient conjugate gradient method for typical parameters of the ENS experiment. We have constructed in this way a phase diagram giving the existence domain of a single vortex configuration as a local minimum of energy, as function of the trapping potential aspect ratio and of the rotation frequency. The numerics confirm the analytical prediction for $`\mathrm{\Omega }_1`$, and give the rotation frequency $`\mathrm{\Omega }_2`$ above which surface modes of the condensate are destabilized and several vortices enter the condensate.

We have also studied the effect of thermal fluctuations of an otherwise straight vortex line in a cigar-shaped condensate. We find values of the mean density contrast of the vortex hole close to the experimental results of ens . Using the Glauber P representation of the density operator of the gas for the field of Bogoliubov quasi-particles, we can predict what happens in a single realization of the experiment: the vortex line experiences some bending, due to the thermal population of low energy kelvon modes localized close to the end points of the condensate. The 1D nature of the kelvon modes leads to a remarkable non Gaussian character of the vortex line fluctuations.

###### Acknowledgements.

We acknowledge useful discussions with Bernard Bernu, Alice Sinatra, Gora Shlyapnikov, Subhasis Sinha, Amandine Aftalion, Tristan Rivière, Jean Dalibard and the members of his group. This work was supported by the BEC2000+ Programme of ESF. M.M. would like to thank the École normale supérieure in Paris for hospitality. Laboratoire de Physique Théorique des Liquides is the Unité Mixte de Recherche 7600 of Centre National de la Recherche Scientifique. Laboratoire Kastler Brossel is a research unit of École normale supérieure and of Université Pierre et Marie Curie, associated to Centre National de la Recherche Scientifique.

## Appendix A Simple derivation of a vortex line energy functional

We present here a derivation of an approximate energy functional for a single vortex line in a Thomas-Fermi condensate subject to a harmonic potential rotating around one of its eigenaxes. Such a derivation is available in the literature AftaRiv  but it is rather complex as it directly includes the effect of curvature of the vortex line in the energy functional, and it applies for an arbitrary aspect ratio of the condensate. The derivation that we propose has the advantage of simplicity as (i) it is specialized to the case of a very elongated condensate along the rotation axis $`z`$ in a cylindrically symmetric trap, that is the atomic oscillation frequency $`\omega _z`$ along $`z`$ is much smaller than the oscillation frequency $`\omega _{}`$ in the transverse $`xy`$ plane, and (ii) it immediately supposes that the explicit dependence of the energy functional on the curvature of the steady state bent vortex line can be neglected. Point (ii) allows to assume locally that the vortex line is straight, which greatly simplifies the derivation of the energy functional. By comparison with the more general derivation of AftaRiv  we have checked the validity of point (ii) verif\_curv .

Furthermore we assume in our derivation that the rotation frequency $`\mathrm{\Omega }`$ is of the order of

$$\mathrm{\Omega }\frac{\mathrm{}\omega _{}^2}{\mu }$$
(58)
where $`\omega _{}`$ is the oscillation frequency of the atoms in the $`xy`$ plane and $`\mu `$ is the chemical potential. It is indeed in this range of rotation frequencies that the single vortex configuration is first stabilized, and this also greatly simplifies the derivation of the energy functional. One then finds an energy difference between the vortex free configuration (of energy $`E_0`$) and a configuration with a vortex of the order of

$$EE_0\frac{\mathrm{}^2\omega _{}^2}{\mu }.$$
(59)
The first step in the derivation of the approximate energy functional is to rewrite the full Gross-Pitaevskii energy functional (1) in terms of the modulus and the phase of the condensate wavefunction $`\varphi (\stackrel{}{r})`$:

$$\varphi (\stackrel{}{r})=n^{1/2}(\stackrel{}{r})e^{iS(\stackrel{}{r})}$$
(60)
where $`n`$ is the probability density normalized to unity. At this stage it is convenient to introduce the so-called velocity field of the condensate:

$$\stackrel{}{v}(\stackrel{}{r})=\frac{\mathrm{}}{m}\stackrel{}{\text{grad}}S(\stackrel{}{r}).$$
(61)
Note that $`\stackrel{}{v}`$ corresponds to the velocity field in the lab frame, the velocity field in the rotating frame being $`\stackrel{}{v}\stackrel{}{\mathrm{\Omega }}\stackrel{}{r}`$. Inserting (60) into the Gross-Pitaevskii energy functional (1) leads to an expression that we split for convenience in two contributions, one linked to the modulus and the other one linked to the phase:

$`E`$ $`=`$ $`E_{\mathrm{mod}}+E_{\mathrm{phase}}`$ (62)
$`E_{\mathrm{mod}}`$ $`=`$ $`{\displaystyle d^3\stackrel{}{r}\left[\frac{\mathrm{}^2}{2m}(\stackrel{}{\text{grad}}\sqrt{n})^2+Un+\frac{1}{2}Ngn^2\right]}`$ (63)
$`E_{\mathrm{phase}}`$ $`=`$ $`{\displaystyle d^3\stackrel{}{r}n\left[\frac{1}{2}m\stackrel{}{v}^2\stackrel{}{\mathrm{\Omega }}(\stackrel{}{r}m\stackrel{}{v})\right]}.`$ (64)
As already mentioned, the trapping potential $`U`$ is axi-symmetric with respect to $`z`$:

$$U=\frac{1}{2}m\left[\omega _{}^2(x^2+y^2)+\omega _z^2z^2\right].$$
(65)
We also give the conditions on $`n`$ and $`\stackrel{}{v}`$ ensuring that $`\varphi `$ is a local extremum of the Gross-Pitaevskii energy functional:

$`0`$ $`=`$ $`\text{div}\left[n(\stackrel{}{v}\stackrel{}{\mathrm{\Omega }}\stackrel{}{r})\right]`$ (66)
$`\mu `$ $`=`$ $`{\displaystyle \frac{1}{2}}m\stackrel{}{v}^2+U+Ngn{\displaystyle \frac{\mathrm{}^2}{2m}}{\displaystyle \frac{\mathrm{\Delta }\sqrt{n}}{\sqrt{n}}}m\stackrel{}{\mathrm{\Omega }}\stackrel{}{r}\stackrel{}{v}`$
where $`\mu `$ is the chemical potential. These conditions are simply the time independent Gross-Pitaevskii equation written in the modulus-phase representation.

### A.1 Energy from the modulus

In the present Thomas-Fermi regime the contribution depending only on the modulus can be evaluated along the lines of our previous work castindum . One splits the density in a slowly varying envelope $`n_{\mathrm{slow}}`$ and a function $`f^2`$ representing the density hole due to the vortex core:

$$n(\stackrel{}{r})=n_{\mathrm{slow}}(\stackrel{}{r})f^2(\stackrel{}{r}).$$
(68)
The envelope $`n_{\mathrm{slow}}`$ varies at the scale of the transverse Thomas-Fermi radius of the condensate $`R_{}=(2\mu /m\omega _{}^2)^{1/2}`$ whereas $`f`$ varies at the scale of the diameter of a vortex core, that is the healing length $`\xi =(\mathrm{}^2/m\mu )^{1/2}`$ convenience .

We now determine the function $`f`$ from the requirement that it deviates significantly from unity only at a distance at most a few $`\xi `$’s from the vortex core. As such a length scale we neglect the spatial variation of $`U`$ and of $`n_{\mathrm{slow}}`$; we further check that in the range (58) of rotation frequencies, the rotational velocity term $`\stackrel{}{\mathrm{\Omega }}\stackrel{}{r}`$ is negligible as compared to the vortex velocity field $`\stackrel{}{v}`$, which allows to neglect the $`\stackrel{}{\mathrm{\Omega }}\stackrel{}{r}`$ terms close to the vortex core both in (66) and (LABEL:eq:sta2):

$$\frac{|\stackrel{}{\mathrm{\Omega }}\stackrel{}{r}|}{v}\frac{\mathrm{\Omega }R_{}}{\mathrm{}/m\xi }\frac{\mathrm{\Omega }}{\omega _{}}\frac{\mathrm{}\omega _{}}{\mu }1.$$
(69)
Furthermore the minimal radius of curvature of the vortex line is found in the subsequent calculations to be of the order of $`R_{}`$ much larger than $`\xi `$. After all these simplifications the function $`f`$ is found locally to solve the well-known Gross-Pitaevskii equation for a straight vortex line in an infinite spatially homogeneous condensate Donnelly , this fictitious homogeneous condensate having a particle density given by $`Nn_{\mathrm{slow}}`$ evaluated on the vortex line. We therefore take for $`f`$:

$$f(\stackrel{}{r})=F(d/\xi _{\mathrm{loc}})$$
(70)
where $`d`$ is the distance of $`\stackrel{}{r}`$ to the vortex line and

$$\xi _{\mathrm{loc}}=\frac{\mathrm{}}{\sqrt{mgNn_{\mathrm{slow}}(\stackrel{}{r}_0)}}$$
(71)
is the local healing length at the position of the vortex line convenience . The function $`F(u)`$ does not depend on any physical parameter. It can be obtained from a numerical solution of the reduced Gross-Pitaevskii equation for a vortex in a homogeneous condensate (see equation (2.84) of Donnelly ). In the large $`u`$ limit, its deviation from unity tends to zero as $`O(1/u^2)`$. For moderate values of $`u`$ it is well approximated by castindum

$$F(u)\mathrm{tanh}(0.7687u).$$
(72)
The slowly varying envelope $`n_{\mathrm{slow}}`$ is obtained from (LABEL:eq:sta2) by removing the short range $`mv^2/2`$ term already included in $`f`$ and by neglecting the quantum pressure term $`\mathrm{\Delta }\sqrt{n}/\sqrt{n}`$ in the spirit of the Thomas-Fermi approximation. We write $`n_{\mathrm{slow}}`$ as $`n_{\mathrm{TF}}+\delta n`$ where

$$n_{TF}(\stackrel{}{r})=\frac{\mu _0U(\stackrel{}{r})}{Ng}$$
(73)
is the usual Thomas-Fermi approximation for the probability density in the absence of vortex and

$$Ng\delta n=\delta \mu +\stackrel{}{\mathrm{\Omega }}\stackrel{}{r}m\stackrel{}{v}.$$
(74)
The density correction $`\delta n`$ includes the rotational term and the deviation $`\delta \mu `$ of $`\mu `$ from the vortex free Thomas-Fermi chemical potential $`\mu _0`$. Both contributions are of the same small order. Using the estimate $`v\mathrm{}/mR_{}`$ we find that the rotational term leads to

$$\frac{\delta n^{\mathrm{rot}}}{n_{\mathrm{TF}}}\frac{\mathrm{}\mathrm{\Omega }}{\mu _0}\left(\frac{\mathrm{}\omega _{}}{\mu _0}\right)^2.$$
(75)
We can also estimate $`\delta \mu `$ by multiplying (LABEL:eq:sta2) by $`n`$ and integrating over the whole space: we find a contribution to $`\delta \mu `$ involving the rotational term,

$$\delta \mu ^{\mathrm{rot}}=\mathrm{\Omega }L_z$$
(76)
of the order of $`\mathrm{}\mathrm{\Omega }`$. A second contribution comes from the fact that the vortex line digs an empty tube of volume $`R_z\xi ^2`$ in the condensate of volume $`R_{}^2R_z`$, where $`R_z`$ is the Thomas-Fermi radius along $`z`$, so that $`\mu `$ has to differ from $`\mu _0`$ by a relative amount

$$\frac{\delta \mu ^{\mathrm{norm}}}{\mu _0}\frac{R_z\xi ^2}{R_{}^2R_z}\left(\frac{\mathrm{}\omega _{}}{\mu _0}\right)^2$$
(77)
to ensure that $`n`$ is normalized to unity.

We now proceed with the calculation of $`E_{\mathrm{mod}}`$, inserting the ansatz (68) in (63). We calculate first the harmonic plus interaction potential energy part of $`E_{\mathrm{mod}}`$, then the kinetic energy part of $`E_{\mathrm{mod}}`$.

In the harmonic plus interaction potential energy terms, we use the identity $`f^2=(f^21)+1`$ and we collect the terms in powers of $`f^21`$:

$`E_{\mathrm{pot}}`$ $`=`$ $`E_{\mathrm{pot}}^{(0)}+E_{\mathrm{pot}}^{(1)}+E_{\mathrm{pot}}^{(2)}`$ (78)
$`E_{\mathrm{pot}}^{(0)}`$ $`=`$ $`{\displaystyle d^3\stackrel{}{r}\left(\frac{Ng}{2}n_{\mathrm{slow}}^2+Un_{\mathrm{slow}}\right)}`$ (79)
$`E_{\mathrm{pot}}^{(1)}`$ $`=`$ $`{\displaystyle d^3\stackrel{}{r}(Ngn_{\mathrm{slow}}+U)n_{\mathrm{slow}}(f^21)}`$ (80)
$`E_{\mathrm{pot}}^{(2)}`$ $`=`$ $`{\displaystyle d^3\stackrel{}{r}\frac{Ng}{2}n_{\mathrm{slow}}^2(f^21)^2}.`$ (81)
The zeroth degree term is of the order of magnitude of $`\mu `$ so that one has to include the deviation of $`n_{\mathrm{slow}}`$ from $`n_{\mathrm{TF}}`$ to first order in order to get the leading term (59) in the vortex energy. Using the fact that $`Ngn_{\mathrm{TF}}+U=\mu _0`$ we obtain

$$E_{\mathrm{pot}}^{(0)}=E_{\mathrm{TF}}+\mu _0d^3\stackrel{}{r}\delta n+o(\mathrm{}^2\omega _{}^2/\mu ),$$
(82)
where we have introduced the Thomas-Fermi approximation to the energy $`E_0`$ of the vortex free configuration:

$$E_0E_{\mathrm{TF}}=d^3\stackrel{}{r}\left(\frac{Ng}{2}n_{\mathrm{TF}}^2+Un_{\mathrm{TF}}\right).$$
(83)
The first degree term $`E_{\mathrm{pot}}^{(1)}`$ and the second degree term $`E_{\mathrm{pot}}^{(2)}`$ are of the order of $`\mathrm{}^2\omega _{}^2/\mu _0`$ as $`|f^21|`$ is close to unity in a cylinder of volume $`R_z\xi ^2`$ and is negligible outside. We can therefore approximate $`n_{\mathrm{slow}}`$ by $`n_{\mathrm{TF}}`$ in $`E_{\mathrm{pot}}^{(1)}`$ and $`E_{\mathrm{pot}}^{(2)}`$. We then get for $`E_{\mathrm{pot}}^{(1)}`$:

$`E_{\mathrm{pot}}^{(1)}`$ $``$ $`\mu _0{\displaystyle d^3\stackrel{}{r}n_{\mathrm{TF}}(f^21)}`$ (84)
$`=`$ $`\mu _0{\displaystyle d^3\stackrel{}{r}\delta nf^2}\mu _0{\displaystyle d^3\stackrel{}{r}\delta n}`$
where we have used the normalization of $`n_{\mathrm{TF}}`$ and of $`n_{\mathrm{slow}}f^2=(n_{\mathrm{TF}}+\delta n)f^2`$ to unity. We thus see that at the present order of the calculation $`E_{\mathrm{pot}}^{(1)}`$ compensates the term linear in $`\delta n`$ in $`E_{\mathrm{pot}}^{(0)}`$. We finally get:

$$E_{\mathrm{pot}}E_{\mathrm{TF}}d^3\stackrel{}{r}\frac{Ng}{2}n_{\mathrm{TF}}^2(f^21)^2.$$
(85)
In this integral we introduce a local system of cylindrical coordinates $`(\rho ,\theta ,Z)`$ with a vertical axis $`Z`$ tangent to the vortex line, see Figure 19, the first coordinate $`\rho `$ being the distance to the vortex line. We then approximate $`n_{\mathrm{TF}}`$ by its local value on the vortex line. The angular integral over $`\theta `$ gives a factor $`2\pi `$. The function $`f^21`$ depends on the distance $`\rho `$ to the vortex line, see equation (70), and the corresponding integral over $`\rho `$ can be extended to infinity as $`(f^21)^2`$ tends rapidly to zero far from the vortex line. We parametrize the vortex line by its curvilinear abscissa $`s`$. We then realize that in the integral over $`Z`$ we can make the reinterpretation $`dZ=ds`$ so that

$`E_{\mathrm{pot}}E_{\mathrm{TF}}`$ (86)
$`{\displaystyle 𝑑s\frac{Ng}{2}n_{\mathrm{TF}}^2(\stackrel{}{r}_0(s))_0^+\mathrm{}2\pi \rho 𝑑\rho \left[F^2(\rho /\xi _{\mathrm{loc}}(s))1\right]^2}`$
where $`\stackrel{}{r}_0(s)`$ is the position of the vortex line at abscissa $`s`$. Finally rescaling $`\rho `$ by $`\xi _{\mathrm{loc}}`$ in the integral over $`\rho `$ and replacing $`n_{\mathrm{slow}}`$ by $`n_{\mathrm{TF}}`$ in (71) leads to

$$E_{\mathrm{pot}}E_{\mathrm{TF}}A_0𝑑s\frac{\pi \mathrm{}^2}{m}n_{\mathrm{TF}}(\stackrel{}{r}_0(s))$$
(87)
where the constant factor $`A_0`$ is

$$A_0=_0^+\mathrm{}u𝑑u\left[F(u)^21\right]^2.$$
(88)
In the kinetic energy term of $`E_{\mathrm{mod}}`$, we neglect the spatial derivative of the slowly varying envelope $`n_{\mathrm{slow}}`$ and of the local healing length $`\xi `$, as they both vary on a length scale $`R_{}`$. We calculate the gradient of $`f`$ in the local system of cylindrical coordinates $`(\rho ,\theta ,Z)`$ of Figure 19:

$$(\stackrel{}{\text{grad}}f)^2=\frac{1}{\xi _{\mathrm{loc}}^2}F^{}_{}{}^{}2(\rho /\xi _{\mathrm{loc}}).$$
(89)
This function vanishes at $`\rho \xi _{\mathrm{loc}}`$ from the vortex line so that we can replace $`n_{\mathrm{slow}}`$ by its value on the vortex line and extend the integral over $`\rho `$ to infinity. The integral over the angle $`\theta `$ gives a factor $`2\pi `$. As in the previous paragraph we use $`dZ=ds`$ where $`s`$ is the curvilinear abscissa on the vortex line, we rescale $`\rho `$ by $`\xi _{\mathrm{loc}}`$ in the integral over $`\rho `$ and we replace $`n_{\mathrm{slow}}`$ by the vortex free Thomas-Fermi expression $`n_{\mathrm{TF}}`$. We then obtain for the kinetic energy part of $`E_{\mathrm{mod}}`$:

$$E_{\mathrm{mod}}^{\mathrm{kin}}A_1𝑑s\frac{\pi \mathrm{}^2}{m}n_{\mathrm{TF}}(\stackrel{}{r}_0(s))$$
(90)
where the constant factor $`A_1`$ is given by

$$A_1=_0^+\mathrm{}u𝑑uF^{}_{}{}^{}2(u).$$
(91)
To summarize the vortex energy stored in the modulus of the condensate wavefunction is approximated by

$$E_{\mathrm{mod}}E_{\mathrm{TF}}(A_0+A_1)𝑑s\frac{\pi \mathrm{}^2}{m}n_{\mathrm{TF}}(\stackrel{}{r}_0(s)).$$
(92)
We estimate $`E_{\mathrm{mod}}E_{\mathrm{TF}}`$ to be indeed of the order of magnitude of $`\mathrm{}^2\omega _{}^2/\mu _0`$, as expected from (59), by taking a length $`R_z`$ for the part of the vortex line inside the condensate and a typical value for the Thomas-Fermi envelope $`n_{\mathrm{TF}}1/R_zR_{}^2`$.

### A.2 Energy involving the velocity field

The energy term $`E_{\mathrm{phase}}`$ is quite difficult to evaluate in its present form (64): since the velocity field tends to zero slowly away from the vortex line (as $`\mathrm{}/md`$ where $`d`$ is the distance to the vortex line), all the parts of the Thomas-Fermi volume of the condensate give approximately the same contribution. A direct 3D integration cannot be performed as the velocity field far from the vortex line is not known explicitly en\_principe , with the notable exception of a straight vortex line on axis $`z`$.

The trick that we borrow from AftaRiv  is to transform the volume integral (64) in a line integral, by ‘integrating by part’ and using the fact that $`\stackrel{}{\text{curl}}\stackrel{}{v}`$ is known exactly:

$$\stackrel{}{\text{curl}}\stackrel{}{v}(\stackrel{}{r})=\frac{2\pi \mathrm{}}{m}𝑑s\stackrel{}{t}(s)\delta (\stackrel{}{r}\stackrel{}{r}_0(s))$$
(93)
where $`s`$ is the curvilinear abscissa along the vortex line and $`\stackrel{}{t}`$ is the unit vector tangent to the vortex line and oriented in the direction of the local vorticity ($`\stackrel{}{t}`$ has a positive component along $`z`$ for a positive rotation frequency $`\mathrm{\Omega }`$). Physically one can use an electromagnetic analogy: $`\stackrel{}{\text{curl}}\stackrel{}{v}`$ corresponds to a linear current, that is a perfectly filiform charge current, of intensity $`2\pi \mathrm{}/m`$ circulating in the vortex line.

The ‘integration by part’ to be performed relies on the following vectorial identity

$$\text{div}(\stackrel{}{A}\stackrel{}{B})=\stackrel{}{B}\stackrel{}{\text{curl}}\stackrel{}{A}\stackrel{}{A}\stackrel{}{\text{curl}}\stackrel{}{B}.$$
(94)
Integrating this identity over whole space and using Ostrogradski’s formula we arrive at the desired identity

$$d^3\stackrel{}{r}\stackrel{}{A}\stackrel{}{\text{curl}}\stackrel{}{B}=d^3\stackrel{}{r}\stackrel{}{B}\stackrel{}{\text{curl}}\stackrel{}{A}$$
(95)
where we have assumed that the total flux of the vector $`\stackrel{}{A}\stackrel{}{B}`$ vanishes through a surface at infinity.

We apply (95) first to transform the contribution $`E_{\mathrm{phase}}^{\mathrm{kin}}`$ of $`\stackrel{}{v}^2`$ to (64):

$$E_{\mathrm{phase}}^{\mathrm{kin}}=d^3\stackrel{}{r}\frac{1}{2}mn\stackrel{}{v}^2.$$
(96)
We have to choose $`\stackrel{}{A}=\stackrel{}{v}`$ so that $`\stackrel{}{\text{curl}}\stackrel{}{v}`$ appears in the left hand side of (95). We then set different :

$$\stackrel{}{\text{curl}}\stackrel{}{B}n\stackrel{}{v}.$$
(97)
This is an approximation as $`n\stackrel{}{v}`$ does not have strictly speaking a vanishing divergence. The vector which has an exactly vanishing divergence is the probability current in the rotating frame, $`n(\stackrel{}{v}\stackrel{}{\mathrm{\Omega }}\stackrel{}{r})`$, see (66). Fortunately we can repeat the reasoning of subsection (A.1). At a distance of up to a few healing lengths $`\xi `$ from the vortex core, the rotation term $`\stackrel{}{\mathrm{\Omega }}\stackrel{}{r}`$ is negligible as compared to the velocity field $`\stackrel{}{v}`$, see (69) so that (66) reduces to $`\text{div}(n\stackrel{}{v})0`$. At a distance much larger than $`\xi `$ from the vortex core the density $`n`$ can be approximated by the vortex free Thomas-Fermi density $`n_{\mathrm{TF}}`$ so that (66) reduces to

$$\text{div}\left[n_{\mathrm{TF}}(\stackrel{}{v}\stackrel{}{\mathrm{\Omega }}\stackrel{}{r})\right]0.$$
(98)
As the Thomas-Fermi density $`n_{\mathrm{TF}}`$ is rotationally symmetric with respect to the rotation axis $`z`$, we have pourquoi

$$\text{div}\left[n_{\mathrm{TF}}\stackrel{}{\mathrm{\Omega }}\stackrel{}{r}\right]=0$$
(99)
so that (98) reduces to $`\text{div}(n\stackrel{}{v})0`$. This finally leads to the desired line integral reformulation for $`E_{\mathrm{phase}}^{\mathrm{kin}}`$:

$$E_{\mathrm{phase}}^{\mathrm{kin}}\pi \mathrm{}𝑑s\stackrel{}{B}(\stackrel{}{r}_0(s))\stackrel{}{t}(s).$$
(100)
The calculation of $`\stackrel{}{B}`$ remains a challenge. Formally $`\stackrel{}{B}`$ can be considered as the static magnetic field created by a current proportional to $`n\stackrel{}{v}`$ so that we have the Biot and Savart formula Jackson :

$$\stackrel{}{B}(\stackrel{}{r})=\frac{1}{4\pi }d^3\stackrel{}{r}^{}\frac{n\stackrel{}{v}(\stackrel{}{r}^{})(\stackrel{}{r}\stackrel{}{r}^{})}{|\stackrel{}{r}\stackrel{}{r}^{}|^3}$$
(101)
but this requires in principle the knowledge of the velocity field $`\stackrel{}{v}`$ everywhere terme\_de\_bord . In practice the problem is simplified by the fact that one needs to know $`\stackrel{}{B}`$ on the vortex line only, $`\stackrel{}{r}=\stackrel{}{r}_0(s)`$ for any fixed $`s`$, and by the fact that the integrand (101) tends rapidly to zero for increasing $`|\stackrel{}{r}_0(s)\stackrel{}{r}^{}|`$ so that it is sufficient to know the velocity field $`\stackrel{}{v}`$ close to the vortex line. As mentioned in the beginning of this appendix we assume that the vortex line is locally straight around the point $`\stackrel{}{r}_0(s)`$. We then approximate $`\stackrel{}{v}`$ by the velocity field of a straight vortex in a homogeneous medium. In the local system of cylindrical coordinates $`(\rho ,\theta ,Z)`$ defined in Figure 19 we thus write

$$\stackrel{}{v}(\stackrel{}{r}^{})\frac{\mathrm{}}{m}\frac{\stackrel{}{e}_\theta }{\rho }.$$
(102)
In this local frame the vector $`\stackrel{}{r}^{}\stackrel{}{r}_0(s)`$ is equal to $`\rho \stackrel{}{e}_\rho +Z\stackrel{}{e}_Z`$ and the unit vector $`\stackrel{}{e}_Z`$ tangent to the vortex line actually coincides with $`\stackrel{}{t}(s)`$ so that one has

$$\stackrel{}{t}(s)\left[\stackrel{}{v}(\stackrel{}{r}^{})(\stackrel{}{r}_0(s)\stackrel{}{r}^{})\right]\frac{\mathrm{}}{m}.$$
(103)
This leads to the rather explicit expression

$$E_{\mathrm{phase}}^{\mathrm{kin}}\frac{\mathrm{}^2}{4m}𝑑sd^3\stackrel{}{r}^{}\frac{n(\stackrel{}{r}^{})}{|\stackrel{}{r}_0(s)\stackrel{}{r}^{}|^3}.$$
(104)
To calculate (104) we write $`n=n_{\mathrm{slow}}f^2n_{\mathrm{TF}}(f^21)+n_{\mathrm{TF}}`$ as in subsection A.1. This leads to a splitting of $`E_{\mathrm{phase}}^{\mathrm{kin}}`$ in two pieces.

The piece involving $`f^21`$ is re-expressed in terms of the local cylindrical coordinates $`(\rho ,\theta ,Z)`$ of the local $`XYZ`$ frame of Figure 19:

$`E_{\mathrm{phase}}^{\mathrm{kin}}(\text{I})`$ $``$ $`{\displaystyle \frac{\mathrm{}^2}{4m}}{\displaystyle 𝑑sd^3\stackrel{}{r}^{}\frac{n_{\mathrm{TF}}(\stackrel{}{r}^{})[f^2(\stackrel{}{r}^{})1]}{\left(|\stackrel{}{r}_0(s)\stackrel{}{r}^{}|^2+ϵ^2\right)^{3/2}}}`$
$`=`$ $`{\displaystyle \frac{\mathrm{}^2}{4m}}{\displaystyle 𝑑sd^3\stackrel{}{R}\frac{n_{\mathrm{TF}}(\stackrel{}{r}^{})[F^2(\rho /\xi _{\mathrm{loc}}(s))1]}{\left(\rho ^2+Z^2+ϵ^2\right)^{3/2}}}`$
where we introduced an arbitrarily small $`ϵ`$ in the denominator to prevent a divergence of the integral. The integrand as a function of $`Z`$ is tending to zero as soon as $`|Z|`$ exceeds a few times $`\rho `$. As $`\rho `$ is limited to a few times $`\xi _{\mathrm{loc}}(s)`$ by the function $`F^21`$, the integrand becomes negligible as soon as $`|Z|`$ exceeds a few times $`\xi _{\mathrm{loc}}(s)`$, which allows to approximate the slowly varying Thomas-Fermi envelope $`n_{\mathrm{TF}}`$ by its value on the vortex line and to extend the integration over $`Z`$ to infinity. Using

$$_{\mathrm{}}^+\mathrm{}𝑑Z\frac{1}{(Z^2+\rho ^2+ϵ^2)^{3/2}}=\frac{2}{\rho ^2+ϵ^2}$$
(105)
and extending in the resulting integral the integration over $`\rho `$ to $`\mathrm{}`$ we obtain the result

$$E_{\mathrm{phase}}^{\mathrm{kin}}(\text{I})\frac{\pi \mathrm{}^2}{m}𝑑sn_{\mathrm{TF}}(\stackrel{}{r}_0(s))_0^+\mathrm{}\rho 𝑑\rho \frac{F^2(\rho /\xi _{\mathrm{loc}}(s))1}{\rho ^2+ϵ^2}.$$
(106)
A more explicit form will be given in the next subsection.

The remaining piece of $`E_{\mathrm{phase}}^{\mathrm{kin}}`$ involves the Thomas-Fermi envelope only:

$$E_{\mathrm{phase}}^{\mathrm{kin}}(\text{II})=\frac{\mathrm{}^2}{4m}𝑑sd^3\stackrel{}{r}^{}\frac{n_{\mathrm{TF}}(\stackrel{}{r}^{})}{\left(|\stackrel{}{r}_0(s)\stackrel{}{r}^{}|^2+ϵ^2\right)^{3/2}}.$$
(107)
We simplify this expression by taking advantage of the cigar shaped nature of the condensate. We first integrate (104) over $`z^{}`$. As in the previous paragraph we use the fact that the integral

$$_{z_0(s)\zeta }^{z_0(s)+\zeta }𝑑z^{}\frac{1}{|\stackrel{}{r}_0(s)\stackrel{}{r}^{}|^3}$$
(108)
converges to its $`\zeta =+\mathrm{}`$ value as soon as $`\zeta `$ exceeds a few times $`|\stackrel{}{r}_0(s)\stackrel{}{r}_{}^{}|`$, where $`\stackrel{}{r}_{}`$ is the projection of the vector $`\stackrel{}{r}`$ in the $`xy`$ plane. The range of $`|\stackrel{}{r}_{}|`$ is limited to the transverse Thomas-Fermi radius of the condensate by the presence of $`n_{\mathrm{TF}}(\stackrel{}{r}^{})`$ in the integrand. The range of $`|\stackrel{}{r}_0(s)|`$ exceeds $`R_{}`$ for a bent vortex line as the vortex line gets out of the Thomas-Fermi profile of the condensate; however the contribution to the energy of the vortex line segments at a distance exceeding a few $`R_{}`$’s becomes much smaller than (58) and is hence negligible, see the appendix B. We can therefore assume that $`|\stackrel{}{r}_0(s)\stackrel{}{r}_{}^{}|`$ is at most a few times $`R_{}`$ in (107) so that the integral over $`z^{}`$ converges over a distance of $`R_{}`$. At such a length scale along the rotation axis $`z`$, $`n_{\mathrm{TF}}(x^{},y^{},z^{})`$ is almost constant for a cigar shaped condensate and can be approximated by $`n_{\mathrm{TF}}(x^{},y^{},z_0(s))`$. We then extend the integration over $`z^{}`$ to infinity and we use (105) to get

$$E_{\mathrm{phase}}^{\mathrm{kin}}(\text{II})\frac{\mathrm{}^2}{2m}𝑑s𝑑x^{}𝑑y^{}\frac{n_{\mathrm{TF}}(x^{},y^{},z_0(s))}{|\stackrel{}{r}_0(s)\stackrel{}{r}^{}|^2+ϵ^2}.$$
(109)
The resulting integral over $`x^{},y^{}`$ can be calculated exactly and was already encountered in the 2D calculation of castindum . We will give the result in the next subsection.

Finally we apply the ‘integration by part’ technique (95) to the last term of $`E_{\mathrm{phase}}`$, the rotational energy term:

$$E_{\mathrm{phase}}^{\mathrm{rot}}d^3\stackrel{}{r}n\left[\stackrel{}{\mathrm{\Omega }}\stackrel{}{r}m\stackrel{}{v}\right]=md^3\stackrel{}{r}\stackrel{}{v}(n\stackrel{}{\mathrm{\Omega }}\stackrel{}{r}).$$
(110)
Note that this term is simply $`\mathrm{\Omega }L_z`$, where $`L_z`$ is the angular momentum per particle along $`z`$. Calculation of this term is considerably simpler than $`E_{\mathrm{phase}}^{\mathrm{kin}}`$. First one can approximate $`n`$ by the Thomas-Fermi envelope $`n_{\mathrm{TF}}`$, neglecting in particular the density hole due to the vortex line plus\_facile . Then one realizes that $`n_{\mathrm{TF}}\stackrel{}{\mathrm{\Omega }}\stackrel{}{r}`$, having a vanishing divergence, see (99), can be written as the curl of some vectorial field $`\stackrel{}{B}`$. One finds inside the Thomas-Fermi condensate AftaRiv

$$\stackrel{}{B}=\frac{Ng}{2m\omega _{}^2}n_{\mathrm{TF}}^2\stackrel{}{\mathrm{\Omega }}$$
(111)
and one takes $`\stackrel{}{B}=\stackrel{}{0}`$ out of the Thomas-Fermi condensate. One then uses (95) and (93) to obtain

$$E_{\mathrm{phase}}^{\mathrm{rot}}\frac{\pi \mathrm{}Ng\mathrm{\Omega }}{m\omega _{}^2}𝑑s\stackrel{}{t}(s)\stackrel{}{e}_zn_{\mathrm{TF}}^2(\stackrel{}{r}_0(s))$$
(112)
where $`\stackrel{}{e}_z`$ is the unit vector along $`z`$. Note that $`\stackrel{}{t}(s)\stackrel{}{e}_z=\mathrm{cos}\alpha (s)`$ where $`\alpha `$ is the angle between the vortex line and the rotation axis $`z`$, so that the integration element $`ds\mathrm{cos}\alpha (s)`$ is simply $`dz`$, the length of the projection of the vortex line along $`z`$. From the estimate $`n_{\mathrm{TF}}1/R_{}^2R_z\mu _0/Ng`$ and for a vortex length $`R_z`$ inside the Thomas-Fermi envelope one gets $`E_{\mathrm{phase}}^{\mathrm{rot}}\mathrm{}\mathrm{\Omega }`$, as expected from $`L_z\mathrm{}`$ for a single vortex configuration. This is in the energy range (59) for a rotation frequency of the order of (58).

### A.3 In terms of a 2D energy functional

We now reexpress the 3D vortex line energy functional that we have obtained in terms of the energy functional of a 2D condensate with a vortex. Physically we associate a 2D fictitious condensate to each horizontal slice of the 3D condensate. Each 2D condensate is stored in a isotropic harmonic potential $`m\omega _{}^2(x^2+y^2)/2`$ with a rotation frequency $`\mathrm{\Omega }`$, and has the same number $`N`$ of particles as the 3D condensate. The fictitious 2D condensate at elevation $`z`$ has a vortex at a position given by the intersection of the 3D vortex line with the horizontal plane of elevation $`z`$. The 2D condensate has a Thomas-Fermi chemical potential

$$\mu _{2\mathrm{D}}(z)=\mu _0\frac{1}{2}m\omega _z^2z^2$$
(113)
where $`\mu _0`$ is the 3D vortex free Thomas-Fermi chemical potential. The Thomas-Fermi radius of the 2D condensate therefore coincides with the one $`R_{}(z)`$ of the 3D condensate at elevation $`z`$:

$$R_{}(z)=\left(\frac{2\mu _{2\mathrm{D}}(z)}{m\omega _{}^2}\right)^{1/2}.$$
(114)
Using the Thomas-Fermi value of the chemical potential for a 2D condensate we get the effective 2D coupling constant of the 2D condensates:

$$g_{2\mathrm{D}}(z)=\frac{\pi \mu _{2\mathrm{D}}^2(z)}{Nm\omega _{}^2}.$$
(115)
We can then relate the 2D Thomas-Fermi density of the 2D condensate to the Thomas-Fermi 3D density of the true condensate:

$$\mu _0U(\stackrel{}{r})=gn_{\mathrm{TF}}(\stackrel{}{r})=g_{2\mathrm{D}}(z)n_{2\mathrm{D}}(\stackrel{}{r}_{};z).$$
(116)
Collecting all the energy terms of this section and after lengthy calculations we produce the result

$$EE_{\mathrm{TF}}𝑑s\frac{g_{2\mathrm{D}}(z_0(s))}{g}\left[E_{2\mathrm{D}}^{\mathrm{\Omega }=0}(r_0(s);z_0(s))+\mathrm{cos}(\alpha (s))E_{2\mathrm{D}}^{\mathrm{rot}}(r_0(s);z_0(s))\right]$$
(117)
where $`\alpha (s)`$ is the angle between the vortex line and the axis $`z`$. The first 2D energy functional term contains the kinetic energy and the harmonic plus interaction potential energy. For a vortex core inside the Thomas-Fermi radius of the condensate, that is for $`\stackrel{~}{r}_{}r_{}/R_{}(z)<1`$, it is given by

$$E_{2\mathrm{D}}^{\mathrm{\Omega }=0}(r_{};z)=\frac{(\mathrm{}\omega _{})^2}{\mu _{2\mathrm{D}}(z)}\left\{\frac{1}{2}+\left(1\stackrel{~}{r}_{}^2\right)\left[C+\mathrm{log}\left[\frac{\mu _{2\mathrm{D}}(z)}{\mathrm{}\omega _{}}\left(1\stackrel{~}{r}_{}^2\right)\right]\right]\right\}$$
(118)
where the value of the constant $`C`$ is

$$C=A_0+A_11+\frac{1}{2}\mathrm{log}2+_0^1𝑑u\frac{F^2(u)}{u}+_1^+\mathrm{}\frac{F^2(u)1}{u}.$$
(119)
For $`\stackrel{~}{r}_{}>1`$ the vortex line at elevation $`z`$ is out of the Thomas-Fermi condensate and only $`E_{\mathrm{phase}}^{\mathrm{kin}}`$ has a non-vanishing contribution to $`E_{2\mathrm{D}}^{\mathrm{\Omega }=0}`$:

$$E_{2\mathrm{D}}^{\mathrm{\Omega }=0}(r_{};z)=\frac{(\mathrm{}\omega _{})^2}{2\mu _{2\mathrm{D}}(z)}\left[1+\left(1\stackrel{~}{r}_{}^2\right)\mathrm{log}\frac{\stackrel{~}{r}_{}^2}{\stackrel{~}{r}_{}^21}\right].$$
(120)
The second energy functional term in (117) is the rotational energy; it is given for $`\stackrel{~}{r}_{}<1`$ by

$$E_{2\mathrm{D}}^{\mathrm{rot}}(r_{};z)\mathrm{}\mathrm{\Omega }(1\stackrel{~}{r}_{}^2)^2$$
(121)
and it vanishes out of the condensate, that is for $`\stackrel{~}{r}_{}>1`$. If one uses the hyperbolic tangent estimate of castindum  for $`F`$, see (72), one gets $`C0.0884`$. The 2D energy functional $`E_{2\mathrm{D}}^{\mathrm{\Omega }=0}+E_{2\mathrm{D}}^{\mathrm{rot}}`$ then coincides exactly with the one of equation (64) in castindum .

### A.4 Improving the 2D energy functional

By a more detailed analysis than in castindum  of the properties of the 2D energy functional $`E_{2\mathrm{D}}^{\mathrm{\Omega }=0}`$ we have identified a pathology that lead to problems in the full 3D energy functional minimization: the derivative of $`E_{2\mathrm{D}}^{\mathrm{\Omega }=0}(r_{};z)`$ with respect to $`\stackrel{}{r}_{}`$ presents a logarithmic singularity on the border of the Thomas-Fermi condensate $`r_{}=R_{}(z)`$, that is it diverges logarithmically to $`+\mathrm{}`$ in $`\stackrel{~}{r}_{}=1^{}`$ and to $`\mathrm{}`$ in $`\stackrel{~}{r}_{}=1^+`$. As a consequence $`E_{2\mathrm{D}}^{\mathrm{\Omega }=0}(r_{};z)`$ has a local minimum inside the Thomas-Fermi condensate, at a distance of the order of the healing length $`\xi _{2\mathrm{D}}=\mathrm{}/(m\mu _{2\mathrm{D}}(z))^{1/2}`$ from the border, see Figure 20. This local minimum of energy is an artifact of the approximations used in the appendix: it is found by a numerical minimization of the full Gross-Pitaevskii energy functional with imaginary time propagation, that a vortex is going out of the condensate to infinity in the absence of trap rotation.

This artifact was not relevant in the 2D case, one just had to restrict to the Thomas-Fermi limit $`\mu _{2\mathrm{D}}\mathrm{}\omega _{}`$ with vortex cores inside the Thomas-Fermi condensate (excluding a thin layer of thickness $`\xi _{2\mathrm{D}}`$ around the boundary). In the 3D case this artifact cannot easily be avoided: the bent vortex necessarily tries to cross the boundary of the condensate, where it encounters the logarithmic singularity and can remain trapped. This problem is also worse in 3D because the 3D vortex line can explore the extremities of the cigar shaped condensate, where the local 2D chemical potential $`\mu _{2\mathrm{D}}(z)`$ can be of the order of $`\mathrm{}\omega _{}`$; the local minimum artifact in $`E_{2\mathrm{D}}^{\mathrm{\Omega }=0}`$ then becomes very pronounced, see Figure 20b, which seriously prevents the vortex line from bending and leaving the condensate.

Heuristically we perform a modification to the energy functional $`E_{2\mathrm{D}}^{\mathrm{\Omega }=0}`$ eliminating the spurious trapping of the vortex line in the absence of rotation. In the region out of the Thomas-Fermi condensate, $`\stackrel{~}{r}_{}>1`$, we realized that $`E_{2\mathrm{D}}^{\mathrm{\Omega }=0}`$ is well approximated by a small $`1/\stackrel{~}{r}_{}`$ expansion of the logarithmic term in (120) so that we take

$$\overline{E}_{2\mathrm{D}}^{\mathrm{\Omega }=0}(r_{};z)=\frac{(\mathrm{}\omega _{})^2}{2\mu _{2\mathrm{D}}(z)}\left(\frac{1}{2\stackrel{~}{r}_{}^2}+\frac{1}{2\stackrel{~}{r}_{}^4}\right).$$
(122)
In the inner Thomas-Fermi region, $`\stackrel{~}{r}_{}<1`$, we eliminate the logarithmic singularity of the derivative at the inner border of the condensate by adding a constant term $`\eta `$ to the argument of the logarithm of (118):

$$\overline{E}_{2\mathrm{D}}^{\mathrm{\Omega }=0}(r_{};z)=\frac{(\mathrm{}\omega _{})^2}{\mu _{2\mathrm{D}}(z)}\left\{\frac{1}{2}+\left(1\stackrel{~}{r}_{}^2\right)\left[C+\mathrm{log}\left[\eta +\frac{\mu _{2\mathrm{D}}(z)}{\mathrm{}\omega _{}}\left(1\stackrel{~}{r}_{}^2\right)\right]\right]\right\}.$$
(123)
These definitions lead to a continuous function $`\overline{E}_{2\mathrm{D}}^{\mathrm{\Omega }=0}`$. The value of $`\eta `$ is adjusted to further ensure continuity of the derivatives of $`\overline{E}_{2\mathrm{D}}^{\mathrm{\Omega }=0}`$ in $`\stackrel{~}{r}_{}=1`$:

$$\eta =e^{3/4C}1.938.$$
(124)
As seen on Figure 20 the regularized function $`\overline{E}_{2\mathrm{D}}^{\mathrm{\Omega }=0}`$ has no local minimum, even for a chemical potential smaller than $`\mathrm{}\omega _{}`$.

In 2D the physical predictions derived from $`\overline{E}_{2\mathrm{D}}^{\mathrm{\Omega }=0}`$ slightly differ in the Thomas-Fermi regime from the ones from $`E_{2\mathrm{D}}^{\mathrm{\Omega }=0}`$. For example the critical rotation frequency $`\mathrm{\Omega }_\mathrm{c}^{2\mathrm{D}}`$ such that the single vortex configuration has the same energy as the vortex free configuration is, from $`\overline{E}_{2\mathrm{D}}^{\mathrm{\Omega }=0}`$:

$$\overline{\mathrm{\Omega }}_\mathrm{c}^{2\mathrm{D}}=\frac{\mathrm{}\omega _{}^2}{\mu _{2\mathrm{D}}}\mathrm{log}\left[e^{C+1/2}\left(\frac{\mu _{2\mathrm{D}}}{\mathrm{}\omega _{}}+\eta \right)\right]$$
(125)
whereas the original one in castindum  is

$$\mathrm{\Omega }_\mathrm{c}^{2\mathrm{D}}=\frac{\mathrm{}\omega _{}^2}{\mu _{2\mathrm{D}}}\mathrm{log}\left[e^{C+1/2}\left(\frac{\mu _{2\mathrm{D}}}{\mathrm{}\omega _{}}\right)\right].$$
(126)
In the Thomas-Fermi limit $`\overline{\mathrm{\Omega }}_\mathrm{c}^{2\mathrm{D}}`$ and $`\mathrm{\Omega }_\mathrm{c}^{2\mathrm{D}}`$ differ by a term scaling as $`\eta \mathrm{}^3\omega _{}^3/\mu _{2\mathrm{D}}^2`$ which is beyond accuracy of the energy functional derivation of the present appendix. The stabilization rotation frequency $`\mathrm{\Omega }_{\mathrm{stab}}^{2\mathrm{D}}`$ such that the single vortex core is a local minimum of energy is also modified:

$$\overline{\mathrm{\Omega }}_{\mathrm{stab}}^{2\mathrm{D}}=\frac{\mathrm{}\omega _{}^2}{2\mu _{2\mathrm{D}}}\mathrm{log}\left[e^{C+1/(1+\eta \mathrm{}\omega _{}/\mu _{2\mathrm{D}})}\left(\frac{\mu _{2\mathrm{D}}}{\mathrm{}\omega _{}}+\eta \right)\right]$$
(127)
but this varies only to second order in $`\eta \mathrm{}\omega _{}/\mu _{2\mathrm{D}}`$.

## Appendix B Energy of the vortex line segments far from the condensate

We give here the behaviour of the field $`\stackrel{}{B}`$ defined by (101) far from the Thomas-Fermi condensate, that is at a distance $`r`$ much larger than the Thomas-Fermi radii of the condensate. This problem is formally equivalent to the calculation in magnetostatics of the magnetic field $`\stackrel{}{B}`$ very far from a localized distribution of current $`\stackrel{}{j}=n\stackrel{}{v}/\mu _0`$ where $`\mu _0`$ is the magnetic permeability of vacuum. This calculation is performed for example in Jackson . From the property $`\text{div}\stackrel{}{j}=0`$ one finds that the leading term in the $`1/r`$ expansion corresponds to a dipolar field created by the magnetic dipole moment $`\stackrel{}{M}`$ of the current distribution:

$$\stackrel{}{B}(\stackrel{}{r})\frac{\mu _0}{4\pi r^3}\left[3\stackrel{}{u}(\stackrel{}{u}\stackrel{}{M})\stackrel{}{M}\right]$$
(128)
where $`\stackrel{}{u}=\stackrel{}{r}/r`$ and the moment $`\stackrel{}{M}`$ is here proportional to the mean angular momentum per particle in the condensate:

$$\stackrel{}{M}\frac{1}{2}d^3\stackrel{}{r}^{}\stackrel{}{r}^{}\stackrel{}{j}(\stackrel{}{r}^{})=\frac{1}{2m\mu _0}\stackrel{}{L}.$$
(129)
As a consequence $`\stackrel{}{B}`$ tends to zero as $`1/r^3`$ at infinity, so that the total flux of $`\stackrel{}{v}\stackrel{}{B}`$ vanishes at infinity, as it was assumed in (95). When the condensate wavefunction is symmetric under reflection with respect to the $`xy`$ plane, as expected for a planar bent vortex contained in the $`xz`$ plane, the $`x`$ and $`y`$ component of $`\stackrel{}{L}`$ vanish and we get

$$\stackrel{}{B}(\stackrel{}{r})\frac{L_z}{8\pi m}\frac{3\stackrel{}{u}(\stackrel{}{u}\stackrel{}{e}_z)\stackrel{}{e}_z}{r^3}.$$
(130)