数理解析研究所講究録 1226 巻 2001 年 176-218

176
Superfluid turbulence and stability of BEC:
using the Gross-Pitaevskii equation
Marc Brachet aa Laboratoire dePhysique Statistique, CNRS UMR 8550, \’Ecole Nor maleSup\’erieure, 24 rue Lhomond, 75231 Paris Cedex 05, France–
Abstract
The Gross-Pitaevskii equation, also called the nonlinear Schr\"odinger equation (NLSE),
describes the dynamics low-temperature superflows and Bose-Einstein Condensates(BEC). We review some recent NLSE based numerical studies of superfluid turbulence and BEC stability. The relations with experiments are discussed.
Key words:–
Contents

1Introduction

2

2Hydrodynamics using the NLSE

3

2.1

3

Madelung transformation

2.2 Sound waves

5

2.3 Vortices in 2and 3D

6

3Superfluid turbulence

7

3.1

8

Tools for vortex dynamics

3.2 Numerical results

13

3.3

Experimental results

13

4

Stability of stationary solutions

16

4.1

Exact solution in 1D
ID

116
fi

4.2

General formulation
ulation119$\mathrm{O}$177
4.3 Branch following methods
superflow around acylinder

21

24

4.4

Stability of

4.5

Stability of attractive Bose-Einstein condensates

33

5Conclusion

39

6Acknowledgments

39

References

39

References40$2\mathrm{D}$1Introduction
The present paper is areview of results, obtained in the last 10 years or so,
by numerically studying the nonlinear Schrodinger equation (NLSE). Direct
numerical simulations (DSN) and branch-following methods were extensively
used to investigate the dynamics and stability of NLSE solutions in 2and 3
space dimensions.

Much work has been devoted to the determination of the critical velocity at
which superfluidity breaks into aturbulent regime [1]. Amathematical model, valid at temperatures low enough for the normal fluid toof superfluid
be negligible, is the nonlinear Schrodinger equation (NLSE), also called theGross-Pitaevskii equation [2-4]. In arelated context, dilute Bose-Einstein condensates (BEC) have been recently produced experimentally [5-7]. The dynamics of these compressible nonlinear quantum fluids is accurately describedby the NLSE allowing direct quantitative comparison between theory and experiment [8].
$4\mathrm{H}\mathrm{e}$Several problems pertaining to superfluidity and BEC can thus be studied
in the framework of the NLSE. In this review, we concentrate on two suchproblems: (i) low-temperature superfluid turbulence [9-11] and (ii) stabilityof BEC in the presence of amoving obstacle [12-14] or an attractive interaction[15]. The paper is organized as follows: in section 2the basic definitions andproperties of the model of superflow are given. Ashort presentation of thehydrodynamic form, through Madelung’s transformation, of NLSE with anarbitrary nonlinearity is derived. Simple solutions are discussed.
Section 3is devoted to superfluid turbulence. The basic tools that are needed
turbulence using NLSE are developed and validated in
to numerically study
Section 3.1. The NLSE numerical results are given in section 3.2. Experimenta$3\mathrm{D}$178
results are given in section 3.3.
The stability of BEC is studied in section 4. Exact ID results are given in
section 4.1 and ageneral formulation of stability is given in section 4.2. Numerical branch following methods are explained in section 4.3. The stabilityof asuperflow around acylinder is studied in section 4.4. The stability of an
attractive Bose-Einstein condensate is studied in section 4.5. Finally, section
5is our conclusion.

2Hydrodynamics using the NLSE
The hydrodynamical form of NLSE with an arbitrary nonlinearity, corresponding to abarotropic fluid with an arbitrary equation of state is introduced inthis section. Basic hydrodynamical features such that acoustic propagationand vortex time independent solutions are also discussed.

2.1

Madelung transformation

The most direct way to understand the generality of the connexion between
the NLSE and fluid dynamics is to consider the following action [16] :$A=2 \alpha\int dt\{d^{3}x\{$(1)$\frac{i}{2}(\overline{\psi}\frac{\partial\psi}{\partial t}-\psi\frac{\partial\overline{\psi}}{\partial t}))-\mathcal{F}\}$

with
(2)

$\mathcal{F}=\int d^{3}x(\alpha|\nabla\psi|^{2}+f(|\psi|^{2}))$

1

where
is acomplex wave field andits complex conjugate,
is apositive real constant and is apolynomial inwith real coefficients:$\psi(\vec{x}, t)$$\overline{\psi}$$\alpha$$f$$|\psi|^{2}\equiv\overline{\psi}\psi$(3)$f(| \psi|^{2})=-\Omega|\psi|^{2}+\frac{\beta}{2}|\psi|^{4}+f_{3}|\psi|^{6}+\ldots+f_{n}|\psi|^{2n}$The NLSE is the Euler-Lagrange equation of motion for(1), it reads$\frac{\partial\psi}{\partial t}=-i\frac{\delta \mathcal{F}}{\delta\overline{\psi}}$,$\psi$

corresponding to

179

uq)

(4)

$\mathrm{i}(\mathrm{o}\mathrm{V}^{2}\mathrm{Q}-\mathrm{e}f’]\mathrm{e}|^{2}))$$\ovalbox{\tt\small REJECT}\ovalbox{\tt\small REJECT}$

Madelung’s transformation $[1,16]$(5)$\psi=\sqrt{\rho}\exp(i\frac{\varphi}{2\alpha})$maps the nonlinear wave dynamics of into equations of motion for afluid ofdensity and velocity. Indeed with the help of (5), (1) can be written$\psi$$\vec{v}=\nabla\varphi$$\rho$$A=- \int dtd^{3}x(\rho\frac{\partial\varphi}{\partial t}+\frac{1}{2}\rho(\nabla\varphi)^{2}+2\alpha f(\rho)+\frac{1}{2}(2\alpha\nabla(\sqrt{\rho}))^{2})$(6)

and the corresponding Euler-Lagrange equations of motion read.

$(\rho\vec{v})=0$(7)$\frac{\partial\varphi}{\partial t}+\frac{1}{2}(\nabla\varphi)^{2}+2\alpha f’(\rho)-2\alpha^{2}\frac{\triangle\sqrt{\rho}}{\sqrt{\rho}}=0$(8)$\frac{\partial\rho}{\partial t}+\nabla$

Without the last term of (8) (the s0-called “quantum pressure” term), theseequations are the continuity and Bernoulli equations [4] for an isentropic,
compressible, irrotational fluid.
It is possible to use this identification to define the corresponding “therm0dynamical functions”. Being isentropic $(S=0)$ , the fluid is barotropic, andthere is only one independent thermodynamical variable. First, the Bernoulli
equation readily gives the fluid’s enthalpy per unit mass ash$=2\alpha f’(\rho)$.

(9)

Second, theterm of (6) corresponds to kinetic energy. Thus the fluid’s
internal energy per unit mass is given by$\frac{1}{2}\rho(\nabla\varphi)^{2}$e$= \frac{2\alpha f(\rho)}{\rho}$

.

(10)

The general thermodynamical identity

h $=e+p/\rho$ ,(11)gives the expression
p $=2\alpha(\rho f’(\rho)-f(\rho))$

(1)

180
for the fluid’s pressure.
The physical dimensions of the variables used in (2) and (3) are fixed by thefollowing considerations. Madelung’s transformation (5) imposes that$[\rho]=ML^{-3}$ and $[\alpha]=L^{2}T^{-1}$ . Using (10), one gets $[f(\rho)/\rho]=T^{-1}$ and thus,from (3), $[\Omega]=T^{-1}$ , $[\beta]=T^{-1}\rho^{-1}$ and $[f_{i}]=T^{-1}\rho^{1-i}$ . Note that, in the caseof aBose condensate of particles of mass $m$ , has the value $\hslash/2m[17]$ .
$[|\psi|^{2}]=$$\alpha$2.2 Sound waves
Dispersion relation

2.2.1

The nature of the extra quantum pressure term in (8) can be understood
through the dispersion relation corresponding to acoustic (density) waves(withpropagating around aconstant density level . Settingin (7) and in the gradient of (8), one gets (keeping only,
the linear terms) :$\rho_{0}$$\mathrm{f}(\mathrm{p}\mathrm{O})=0)$$\rho=\rho_{0}+\delta\rho$$\nabla\varphi=\delta u$$\partial_{t}\delta\rho+\rho_{0}\nabla\delta u=0$$\partial_{t}\delta u+2\alpha f’(\rho_{0})\nabla\delta\rho-2\alpha^{2}\Delta\frac{\nabla\delta\rho}{2\rho_{0}}=0$or.

$\partial_{t^{2}}\delta\rho=2\alpha\rho_{0}f’(\rho_{0})\triangle\delta\rho-\alpha^{2}\Delta^{2}\delta\rho$

The dispersion relation for an acoustic wave
) is thus(with$\delta\rho=\epsilon(\exp(i(\omega t -\vec{k}\cdot\vec{x}))+c.c.)$$\epsilon\ll 1$(13)$=\sqrt{2\alpha\rho \mathrm{o}f’(\rho_{0})\vec{k}^{2}+\alpha^{2}\vec{k}^{4}}$$\omega$

This relation shows that the quantum pressure has adispersive effect that becomes important for large wave numbers. For small wavenumbers, one recovers
the usual propagation, with asound velocity given byc$=( \frac{\partial p}{\partial\rho})^{\frac{1}{2}}=\sqrt{2\alpha\rho_{0}f’(\rho_{0})}$.

The length scale
is known as the “coherence length”$\xi=\sqrt{\alpha/(\rho_{0}f’(\rho_{0}))}$

at which dispersion becomes noticeable181
2.2.2 Nonlinear acoustics
The description given by linear acoustic can be somewhat improved by including the dominant nonlinear effects. Such an equation was derived in [18].

Numerical simulations of NLSE in one space dimension using astandard
Fourier pseud0-spectral method [19] can be used to study the acoustic regime
triggered by an initial disturbance of the form :$\psi(x)=1+ae^{-\frac{x^{2}}{\iota^{2}}}$.

Such simulations were performed in ref. [18] where it was found that the shocks
which would have appeared under compressible Euler dynamics (i.e. following
Eq. (8) without the last term in r.h.s.) are regularized by the dispersion. There
was no evidence of finite-time singularity in our numerics: the spectrum of thesolution was well resolved, with aconspicuous exponential tail.

2.3

Vortices in 2and 3D

Further insight on the connexion between the NLSE and fluid dynamics canbe obtained by considering stationary solutions of the equations of motion.
Indeed, by inspection of (1), time independent solutions of NLSE (4), are also
solutions of the Real Ginzburg-Landau Equation (RGLE)$\frac{\partial\psi}{\partial t}=-\frac{\delta \mathcal{F}}{\delta\overline{\psi}}=$$(\alpha\nabla^{2}\psi-\psi f’(|\psi|^{2}))$

(14)

.

They are thus extrema of the free energy$\mathcal{F}$.

The simplest solution of this type corresponds to aconstant density fluid atrest. In this simple case, is constant in space and (14) reads$\psi$.

(15)$f’(|\psi|^{2})=-\Omega+\beta|\psi|^{2}+3f_{3}|\psi|^{4}+\ldots+nf_{n}|\psi|^{2n-2}=0$, relates
This equation, for given values of the coefficients and , $i=3$ ,Q.
Note that the term of does notto the value ofthe fluid density
play acrucial role in the NLSE dynamics. Indeed, it could be removed fromthat
the Bernoulli equation (8) by the change of variable
better,byin NLSE (4). It is however
amounts to achange of phaseconvention, not to perform these changes of variable. With this convention,
stationary solutions of (14) coincide with stationary solutions of (4). The
term of is thus be fixed by the fluid’s density through (15).
$\beta$$f_{i}$$\Omega$$|\psi|^{2}$$\ldots$$n$$f$$\varphiarrow\varphi+2\alpha\Omega t$$\psi$$arrow\psi e^{i\Omega t}$$\Omega$$f$

Another important type of time-independent solutions of NLSE are the vortex
solutions. Madelung’s transformation is singular when $\rho=0$ (i.e. when both

182
and $\propto s(\psi)=0$ . As two conditions are required, the singularities
generically happen on points in two dimensions and lines in three dimensions.
. These topologThe circulation of around such ageneric singularity is
ical defects are known in the context of superfluidity as “quantum vortices”[1]. Solutions of (14) with cylindrical symmetry can be obtained numerically[20]. The density profile of avortex admits ahorizontal asymptote near thecore while the velocity diverges as the inverse of the core distance. The m0is thus aregular quantity. It is important to realize that
mentum density
such vortex solutions are regular solutions of the NLSE (4), the singularity
stemming only from Madelung’s transformation (5).
$\Re(\psi)=0$$\pm 4\pi\alpha$$\vec{v}$$\rho\vec{v}$3Superfluid turbulenceSuperfluid flows ( $i.e$ . laboratory $4He$ flows) are described mathematically interms of Landau’s tw0- luid model [4]. When both normal fluid and superfluid
vortices are present, their interaction, called “mutual friction”, must be taken
into account as pioneered by Schwarz [21]. At temperatures low enough for thefor helium at normal
normal fluid to be negligible (in practice below
pressure), an alternative mathematical description is given by the Nonlinear
Schrodinger Equation (NLSE), also called the Gross-Pitaevskii equation $[2,3]$ .
In this section, we will use the simplest form for , corresponding to acubic
nonlinearity in the NLSE (4). The NLSE, with convenient normalization, reads$T=1^{\mathrm{o}}K$$f$$\partial_{t}\psi=ic/(\sqrt{2}\xi)(\psi-|\psi|^{2}\psi+\xi^{2}\nabla^{2}\psi)$

.

(16)

Madelung’s transformation (5) takes the form
(17)$\rho=|\psi|^{2}$(18)$\rho v_{j}=ic\xi/\sqrt{2}(\psi\partial_{j}\overline{\psi}-\overline{\psi}\partial_{j}\psi)$

where 4is the s0-called “coherence length” and is the velocity of sound (whenthe mean density $\rho_{0}=1[1])$ . The superflow is irrotational, except near thenodal lines of which are known to follow Eulerian dynamics $[22,23]$ . These
topological defects correspond to the superfluid vortices that appear naturally,
with the correct velocity circulation, in this model [17].
$c$$\psi$The basic goal of the present section is to qualify the degree of analogy between
turbulence in low-temperature superfluids and incompressible viscous fluids.
We will do this by comparing numerical simulations of NLSE with existing
numerical simulations of the Navier-Stokes equations, in particular the TaylorGreen (TG) vortex [24]. The TG vortex is the solution of the Navier-Stokes

183
equations with initial velocity field$\mathrm{v}^{TG}=($$\sin(x)\cos(y)\cos(z),$ $-\cos(x)\sin(y)\cos(z)$ ,

0)

(19)

. This flow is well documented in the literature [25-27]. It admits symmetries
that are used to speed up computations: rotation by about the axis $(x=$$z=\pi/2)$ , $(y=z=\pi/2)$ and $(x=y=\pi/2)$ and reflection symmetry withrespect to the planes $x=0$ , , $y=0$ , , $z=0$ , . The velocity is parallel tothese planes which form the sides of the impermeable box which confines the
flow.
$\pi$

$\pi$

$\pi$

$\pi$

Tools for vortex dynamics

3.1

Under compressible fluid dynamics, an arbitrary chosen initial condition will
generally lead to aregime dominated by acoustic radiation. In order to study
vortex dynamics using NLSE, we thus need to prepare the initial data in such
away that the acoustic emission is as small as possible.

3.1.1

Preparation method

We now show how to construct avortex array whose NLSE dynamics mimics. The first step of our method isthe vortex dynamics of the large scale flow
based on aglobal Clebsch representation ofand the second step minimizes
the emission of acoustic waves [28].
$\mathrm{v}^{TG}$$\mathrm{v}^{TG}$The Clebsch potentialsy,$z)=\cos(x)\sqrt{2|\cos(z)|}$

(20)

y,

(21)

$\lambda(x,$$\mu(x,$$z)=\cos(y)\sqrt{2|\cos(z)|}\mathrm{s}\mathrm{g}\mathrm{n}(\cos(z))$(where sgn gives the sign of its argument) correspond to the TG flow in thesense thatand Aand are periodic functions of $(x, y, z)$ .
These Clebsch potentials are used to map the physical space $(x, y, z)$ into the, corresponding to the large scale TG flow
plane. The complex field
circulation, is given bywith $\gamma_{d}=2\sqrt{2}/(\pi c\xi)([]$
denotes the integer part of areal) and$\nabla\cross \mathrm{v}^{TG}=\nabla\lambda$$(\lambda, \mu)$$\cross\nabla\mu$$\mu$$\psi_{c}$$\psi_{c}(x, y, z)=(\psi_{4}(\lambda, \mu))^{[\gamma d/4]}$$\psi_{4}(\lambda, \mu)=\psi_{e}(\lambda-1/\sqrt{2}, \mu)\psi_{e}(\lambda, \mu-1/\sqrt{2})\cross$(22)$\psi_{e}(\lambda+1/\sqrt{2}, \mu)\psi_{e}(\lambda, \mu+1/\sqrt{2})$

where

.

$\psi_{e}(\lambda, \mu)=(\lambda+i\mu)\tanh(\sqrt{\lambda^{2}+\mu^{2}}/\sqrt{2}\xi)/\sqrt{\lambda^{2}+\mu^{2}}$184

for the Tayvisualization of the vector field
Fig. 1.
$\xi=0.1/(8\sqrt{2})$, sound velocity
lor-Green flow at time $t=0$ with coherence length$c=2$ and $N=512$ in the impermeable box $[0, \pi]$.
$\nabla\cross(\rho\vec{v})$$\mathrm{T}\mathrm{h}\mathrm{r}\mathrm{e}\triangleright \mathrm{d}\mathrm{i}\mathrm{m}\mathrm{e}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{a}\mathrm{l}$$\cross[0, \pi]\cross[0, \pi]$

The second step of our procedure consists of integrating to convergence theAdvective Real Ginzburg-Landau Equation (ARGLE):

Ae$=c/(\sqrt{2}\xi)(\psi-|\psi|^{2}\psi+\xi^{2}\nabla^{2}\psi)-i\mathrm{v}^{TG}\cdot\nabla\psi-$$(\mathrm{v}^{TG})^{2}/(2\sqrt{2}ae)\psi$(23). It is shown in [9] that the TG symmetries can bewith initial data
used to expand $\psi(x,$y, z, t), solution of the ARGLE and NLSE equations as:$\psi=\psi_{c}$$\psi(x,$y, z,n,p, t) $\cos mx\cos ny\cos pz$$t)= \sum_{m=0}^{N/2}\sum_{n=0}^{N/2}\sum_{p=0}^{N/2}\hat{\psi}(m,$(24)

where $N$ is the resolution and $\hat{\psi}(m, n,p, t)=0$ , unless $m$ , $n,p$ are either all evensatisfies the additional conditions
or all odd integers. Furthermore
$r=1$
where
when $m$ , $n,p$ are all even and$r=2$ when $m$ , $n,p$ are all odd. Implementing this expansion in apseud0spectral code yields asaving of afactor 64 in computational time and memory
size when compared to general Fourier expansions.
$\hat{\psi}(m, n,p, t)$$\hat{\psi}(m, n,p, t)=(-1)^{\mathrm{r}+1}\hat{\psi}(n, \mathrm{n},\mathrm{p}, t)$

The ARGLE converged periodic vortex array obtained in this manner is displayed on Fig. 1. with coherence length $=0.1/(8\sqrt{2})$ , sound velocity $c=2$
and resolution $N=512$
$\xi$

185

. The bottom
Fig. 2. Plot of the incompressible kinetic energy spectrum,
curve (a) (circles) corresponds to time $t=0$ . The spectrum of asingle axisymmetricvortex multiplied by $(l/2\pi)=175$ is shown as the bottom solid line. The top
curve (b) (pluses) corresponds to time $t=5.5$ . Aleast-square fit over the interval$2\leq k\leq 16$ with apower law $E_{kin}^{i}(k)=Ak^{-}n$ gives $n=1.70$ (top solid line).
$E_{kin}^{i}(k)$$2\mathrm{D}$

3.1.2 Energy spectra
The total energy of the vortex array, conserved by NLSE dynamics, can be, with kinetic
decomposed into three parts$\mathcal{E}_{int}=(c^{2}/2)(\rho-1)^{2}$and quantum, internal energy
energyas the intebe
definedcan. Each of these parts
energy. In order togral of the square of afield, for example,
canseparate the kinetic energy corresponding to compressibility effects,
using
parts
incompressible
be further decomposed into acompressible and. Using Parseval’s theorem, thewithangle-averaged kinetic energy spectrum is defined as:$E_{tot}=1/(2 \pi)^{3}\int d^{3}x(\mathcal{E}_{kin}+\mathcal{E}_{int}+\mathcal{E}_{q})$$\mathcal{E}_{kin}=1/2\rho v_{j}v_{j}$$\mathcal{E}_{q}=c^{2}\xi^{2}(\partial_{j}\sqrt{\rho})^{2}$$\mathcal{E}_{kin}=1/2(\sqrt{\rho}v_{j})^{2}$$\mathcal{E}_{kin}$$\nabla.(\sqrt{\rho}v_{j})^{i}=0$$\sqrt{\rho}v_{j}=(\sqrt{\rho}v_{j})^{c}+(\sqrt{\rho}v_{j})^{i}$$E[in(k)= \frac{1}{2}\int k^{2}\sin\theta d\theta d\phi$$| \frac{1}{(2\pi)^{3}}\int d^{3}re^{ir_{j}k_{j}}\sqrt{\rho}v_{j1^{2}}$. The incompressible
which satisfies
spectrum computing
single-averaged, is thekinetic energy spectrum,
over shells in Fourier space. Amode $(m, n,p)$ belongs to the shell numbered
as $k=[\sqrt{m^{2}+n^{2}+p^{2}}+1/2]$ .
$E_{kin}=1/(2 \pi)^{3}\int d^{3}x\mathcal{E}_{kin}=\int_{0}^{\infty}dkE_{kin}(k)$$E_{kin}^{i}(k)$The radius of curvature of the vortex lines in Fig. 1is large compared tolines can be considered as straight, and then
their radius. Thus theseaxisymmetric vortices which are exact solutions to thecompared to the$3\mathrm{D}$$2\mathrm{D}$

10

186

”

$\mathrm{a}\leq$, plotted versus time forFig. 3. Total incompressible kinetic energy,( $=0.1/(2\sqrt{2})$ , $N=128$ (long-dash line); $\xi=0.1/(4\sqrt{2})$ , $N=256$ (dash);$\xi=0.1/(6.25\sqrt{2})$ , $N=400$ (dot) and $\xi=0.1/(8\sqrt{2})$ , $N=512$ (solid line). Allruns are realized with $c=2$ . The evolution of the total vortex filament length di(crosses) for the $N=512$ run is also shown (scale given on the rightvided by$E_{kin}^{i}$$2\pi$

y-axis).

Fig. 4. Same visualization as in Fig. 1but at time t $=4$ .

11

187

Fig. 5. Same visualization as in Fig. 1but at time t $=8$ .

$2\mathrm{D}$NLSE. A$2\mathrm{D}$vortex at the origin is given by $\psi^{vort}(r)=\sqrt{\rho(r)}\exp(im\varphi)$ ,where
are polar coordinates. The vortex profile
asand $\sqrt{\rho(r)}=1+O(r^{-2})$ for $rarrow\infty$ . It can be computed numerically using mapped Chebychev polynomials and an appropriate functional [9].
The corresponding velocity field is azimuthal and is given by $v(r)=\sqrt{2}c\xi/r$ .
Using the mapped Chebychev polynomials expansion for, the angle averaged spectrum of
can then be computed with the formula, whereis the zeroth order Bessel func$m=\pm 1$ ,$\sqrt{\rho(r)}\sim r$$(r, \varphi)$$rarrow \mathrm{O}$$\sqrt{\rho(r)}$$E_{kin}^{vot\iota}(k)=$$\sqrt{\rho}v_{j}$$c^{2} \xi^{2}/(2\pi k)(\int_{0}^{\infty}drJ_{0}(kr)\partial_{r}\sqrt{\rho})^{2}[9]$$J_{0}$

tion.

The incompressible kinetic energy spectrumof the ARGLE converged
vortex array of Fig. 1is displayed on Fig. 2. For large wavenumbers, the spectrum is well represented by extending acollection ofvortices into
vortex$E_{kin}^{line}(k)\equiv
const$
lines via.
. (We will see that the constant of proportionality is related to the length of vortex lines by const. $=l/(2\pi)=175$ attime $t=0.$ ) In contrast, the small wavenumber region cannot be representedby. This stems from the average separation distance between the vortex
lines in Fig. 1. Calling this distance $d_{bump}=k_{bump}^{-1}=1/16$ , the wavenumber range between the large-scale wavenumber $k=2$ and the characteristic
separation wavenumber $k_{bump}$ can be explained by interference effects. Dueto constructive interference, the energy spectrum at $k=2$ has avalue close
to its corresponding value in TG viscous flow (namely 0.125), which is much
above the value of $E_{kin}^{line}(k=2)$ . In contrast, for $2<k\leq k_{bump}$ , destructive
interference decreases.
below$E_{kin}^{i}$$2\mathrm{D}$$\cross E_{kin}^{vort}(k)$$l$$E_{kin}^{line}$$E_{kin}^{i}$$E_{kin}^{line}$

12$3\mathrm{D}$188
3.2 Numerical results
The evolution in time via NLSE (16) of the incompressible kinetic energy isshown in Fig. 3. The main quantitative result is the remarkable agreement, with the corresponding data in theof the energy dissipation rate,incompressible viscous TG flow (see reference [25], and reference [29], figure5.12). Both the moment $t_{\max}\sim 5-10$ of maximum energy dissipation (theat that moment are in
inflection point of Fig. 3) and its value
depend weakly
andquantitative agreement. Furthermore, bothon .
$-dE_{k\dot{\iota}n}^{\dot{l}}/dt$$\epsilon(t_{\max})\sim 10^{-2}$$t_{\max}$$\epsilon(t_{\max})$$\xi$

Another important quantity studied in viscous decaying turbulence is the scaling of the kinetic energy spectrum during time evolution and, especially, at
range can be obthe moment of maximum energy dissipation, where aspectrumat $t=5.5$ . A
energyserved (see reference [25]). Fig. 2(b) shows theleast-square fit over the interval $2\leq k\leq 16$ with apower law $E_{kin}^{i}(k)=Ak^{-n}$gives $n=1.70$ (solid line). For $5<t<8$ , asimilar fit gives $n=1.6\pm 0.2$ (datain the interval $30\leq k\leq 170$ with $l/(2\pi)$ times
not shown). Fitting
leads to $l/2\pi=452$ , roughly three times the $t=0$ length of the vortex lines. The time evolution of $l/2\pi$ obtained by this procedure is displayed. The
in Fig. 3, showing that the length continues to increase beyond
computations were performed with $c=2$ corresponding to aroot-mean-square. As it is very costly to decrease $M_{rms}$ ,Mach number
we checked [9] that compressible effects were non-dominant at this value of$M_{rms}$ .
$k^{-5/3}$$E_{kin}^{l}(k)$$E_{kin}^{vort}(k)$$t_{\max}$$M_{tms}\equiv|\mathrm{v}_{tms}^{TG}|/c=0.25$

The vortex lines are visualized in physical space in Figs. 4and 5at time
$t=4$ and $t=8$ . At $t=4$ , no reconnection has yet taken place while acomplex vortex tangle is present at $t=8$ . Detailed visualizations (data notshown) demonstrate that reconnection occurs for $t>5$ . Note that the viscous
TG vortex also undergoes aqualitative (and quantitative) change in vortex
dynamics around $t\sim 5$ .

3. 3Experimental results
The TG flow is related to an experimentally studied swirling flow [30-32]. The
relation between the experimental flow and the TG vortex is asimilarity inoverall geometry [30]: ashear layer between two counter-rotating eddies. The
TG vortex, however, is periodic with free-slip boundaries while the experimental flow is contained inside atank between two counter-rotating disks.
The spectral behavior of NLSE can be compared to standard (viscous) turbu. It is thus of interest to estimate the scaling oflence only for k$k_{bump}$$\leq k_{bump}$

13

189
in terms of the characteristic parameters of the large scale flow and of thefluid. As seen above,, whereis the average distance between neighboring vortices. Consider aflow with characteristic integral scaleand large scale velocity(in the case of the TG flow,).
andThe fluid characteristics are the velocity of sound and the coherence length4(with corresponding wavenumber). The numberof vortex lines
crossing atypical large-scale area is given by the ratio of the large-scale flow
circulation
to the quantum of circulation,.
On the other hand, the assumption that the vortices are uniformly spread overthe large scale area gives. Equating these two evaluations ofyields the relation.
$k_{bump}\sim d_{bump}^{-1}$$l_{0}$$d_{bump}$$l_{0}\sim 1$$u_{0}$$u_{0}\sim 1$$c$$k_{\xi}\sim\xi^{-1}$$n_{d}$$l_{0}^{2}$$\Gamma=4\pi c\xi/\sqrt{2}$$l_{0}u_{0}$$\mathrm{i}.\mathrm{e}n_{d}\sim l_{0}u_{0}/c\xi$$n_{d}\sim l_{0}^{2}/d_{bump}^{2}$$n_{d}$$d_{bump}\sim l_{0}\sqrt{(c\xi)/(l_{0}u_{0})}$In the case of helium, the viscosity at the critical point $(T=5.174^{\mathrm{o}}K$ , $P=$
$2.210^{5}Pa)$ is $\nu_{cp}=0.27\cross 10^{-7}m^{2}s^{-1}$ while the quantum of circulation,$\Gamma=h/m_{He}$ has the value $0.99\cross 10^{-7}m^{2}s^{-1}$ . Thus,. The order of$d_{bump}$
magnitude foris thuswhere
is the integral
scale Reynolds number at the critical point and the Taylor micr0-scale. In
other words, the value ofin asuperfluid helium experiment at
is of the same order as the Taylor micr0-scale in the same experimental set-up
run with viscous helium at the critical point.
$0.25\nu_{cp}\sim\Gamma$$d_{bump}=l_{0}/\sqrt{R_{cp}}\sim l_{\lambda}$$R_{cp}$$l_{\lambda}$$T=1^{\mathrm{o}}K$$d_{bump}$

The experimental set-up is similar to the one described in [31]. Some modifications have been made to work down to. The flow is produced in acylinder, 8cm in diameter and 12 cm high, limited axially by two counterrotating disks. One disk is flat and, on the other one, are fixed 8radial blades,
forming an angle ofbetween each other. Astator is mounted at half thetotal height of the cell in order to stabilize the turbulent shear region. The
two disks are driven by two DC motors rotating from 1to 30 Hz. The wholesystem is enclosed directly in aliquid Helium bath which is used as the experimental fluid, the main difference with the set-up described in [31]. The
temperature of the fluid is fixed by the pressure above the liquid bath, which
is itself controlled by the pumping system.
$1.2\mathrm{A}\mathrm{T}$$45^{o}$Local pressure fluctuations are measured by using small total-head pressure
tubes, immersed in the flow. The pressure sensors are hollow metallic tubes,
connected to aquartz pressure transducer WHM 112 A22 from PCB. Details
are given in [11].
In normal fluids, the pressure measured at the tip of the total-head tube can be
related to the upstream flow $U(t)$ and the local pressure $P(t)$ using Bernoulli
theorem:(25)$P_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}}(t)=P(t)+\rho U^{2}(t)/2$

14

190

1
1

.
$\vee\wedge--\cdot$.
.
$\cdot 3$$1|^{1}$1 .1$\check{\mathrm{r}}\mathrm{c}$11$\cdot\cdot 1$

111$\cdot$.

frequency (Hz)1$\cdots$

Fig. 6. Experimental pressure fluctuation spectrum (in non-dimensional units) measured with atotal head pressure tube immersed in the flow at $T=2.3K$ .

In the flow region where the probe is immersed, awell established axial mean
flow exists so that, after removing the mean parts of Equation (1), one gets:$U$

Pmeas $(t)=p(t)+\rho Uu(t)$(26)

where pmeas, and are the fluctuations of the measured pressure, the actualpressure, and the local velocity respectively. It is currently admitted that, inordinary turbulent situations, and at low fluctuation rates, Equation (2) isdominated by the dynamic term, so that, by measuring the pressure fluctuations at the total head tube, one has adirect access to the velocity fluctuations$u(t)$ .
$p$$u$

The situation is less clear when the probe is immersed in the superfluid. It is
however possible to write an equation similar to (2). Details can be found in[11].

The analysis of the pressure fluctuations obtained with the total head tube
placed at 2cm above the mid plane and 2cm from the cylinder axis, yields
interesting informations. Figures 6and 7shows the spectra of the pressure(i.e respectively at $2.3K$ and $1.4K$ ). Fig. 6
fluctuations above and below
clearly shows, as expected, that such fluctuations follow aKolmogorov regime
between the injection scale (signaled by the peak at 25 Hz) and the largest
Hz. The spectrum obtained at 1AK is similar toresolved frequency,(see Fig. 7). Aclear Kolmogorov like regime existsthat obtained at
for the same range of frequencies. The corresponding Kolmogorov constantturns out also to be indistinguishable from the classical value. We have further analyzed the deviations from Kolmogorov in the superfluid regime. The$T_{\lambda}$$\mathrm{i}.\mathrm{e}900$$T=2.\mathrm{S}\mathrm{K}$

15

191

.

$\vee-\wedge-\cdot$$\check{\mathrm{E}}\circ$1$0^{*}$$1$$0^{3}$

11$0^{l}$$0^{1}$

10
1$0^{\cdot}1$

10100

freq uen cy$1\cdot 00$$(\mathrm{H}\mathrm{z})$Fig. 7. Experimental pressure fluctuation spectrum (in non-dimensional units) measured with atotal head pressure tube immersed in the flow at T $=1.4$ K.

striking result is that they have the same magnitude as in classical turbulence.
More details are given in [33].

These observations -both on global and local quantities -agree pretty well
the theoretical approach developed in the previous section. In particular, it
seems rather clear that Kolmogorov cascade survives in the superfluid regime.

4Stability of stationary solutions
This section is devoted to the stability of BEC. Exact ID bifurcation results are
given in section 4.1. Ageneral formulation of stability is given in section 4.2.
Numerical branch following methods are explained in section 4.3. The stabilityof asuperflow around acylinder is studied in section 4.4 and attractive BoseEinstein condensates are studied in section 4.5.

4.1

Exact solution in ID

4.1.1 Definition of the system
We consider apunctual impurity moving within aID superflow. In the frameof the moving impurity, the system can be described by the following action

16

192
functional(27)$A[ \psi, \overline{\psi}]=\int \mathrm{d}t[\frac{i}{2}\int \mathrm{d}x(\overline{\psi}\partial_{t}\psi-\psi\partial_{t}\overline{\psi})-\mathcal{K}]$In this expression,tional reads$\psi$is acomplex field,$\overline{\psi}$

its conjugate and the energy func-$\mathcal{K}$$\mathcal{K}=\mathcal{E}-v\mathcal{P}+v[R^{2}(+\infty)\phi(+\infty)-R^{2}(-\infty)\phi(-\infty)]$

(28)

,

with
$\mathcal{E}=\int \mathrm{d}x[|\partial_{x}\psi|^{2}+\frac{1}{2}(|\psi|^{2}-1)^{2}+g\delta(x)(|\psi|^{2}-1)]$

(29)
(30)

,

$\mathcal{P}=\int \mathrm{d}x\frac{1}{2i}[\overline{\psi}(\partial_{x}\psi)-\psi(\partial_{x}\overline{\psi})]$$\psi=R\exp(i\phi)$

,

.

(31)

in (29) represents the impurity and the
The Dirac (pseudo) potential
boundary conditions for the phaseappropriate
imposes(28)
thelast term in$[34]$ .
obeys the boundary conditions $R^{2}(\pm\infty)=1$ .
$g\delta(x)$$R$$\phi$The Euler-Lagrange equation associated to (27),Schr\"odinger equation (NLSE)$\delta A/\delta\overline{\psi}=0$$i\partial_{t}\psi=-\partial_{xx}\psi+iv\partial_{x}\psi-\psi+|\psi|^{2}\psi+g\delta(x)\psi$, is the nonlinear

(32)

,

where the discontinuity condition
(33)

t)

$\partial_{x}\psi(0^{+}, t)-\partial_{x}\psi(0^{-}, t)=g\psi(0,$is imposed in order to balance thefor all times t.

4.1.2$g\delta(x)\psi$

singularity with the$-\partial_{xx}\psi$

term

Stationary solutions

Time-independent solutions of the NLSE (32) are best studied by performing
the change of variables defined above in (31): $\psi=R\exp(i\phi)$ . Using thesevariables, the NLSE reads$\partial_{t}R=v\partial_{x}R-R\partial_{xx}\phi-2\partial_{x}R\partial_{x}\phi$,$\partial_{t}\phi=v\partial_{x}\phi-(\partial_{x}\phi)^{2}+1-R^{2}-g\delta(x)+\frac{\partial_{xx}R}{R}$

(34)
,

(35)

and the jump condition (33) readst),$\partial_{x}R(0^{+}, t)-\partial_{x}R(0^{-}, t)=gR(0,$$\partial_{x}\phi(0^{+}, t)-\partial_{x}\phi(0^{-}, t)=0$

.

17

(36)
(37)

193
Note that Eqs. (34) and (35) can be respectively interpreted as the continuity
and velocity u
and Bernoulli equations for afluid of density p(as done in section 2).
$\ovalbox{\tt\small REJECT}$$\#^{2}(\ovalbox{\tt\small REJECT} \mathrm{z})$$\ovalbox{\tt\small REJECT}$$2\mathrm{a}_{x}$Explicit time-independent solutions of Eqs. (34) and (35) were found by Hakim[lusing the s0-called gray solitons (a nonlinear optics terminology). Gray solitons$[35,36]$ are stationary solutions of Eqs. (34) and (35), without the potential. They are localized density depletion of the form
term
$g\delta(x)$$R_{\mathrm{G}\mathrm{S}}^{2}(x)=v^{2}/2+(1-v^{2}/2)\tanh^{2}[\sqrt{1/2-v^{2}/4}x]$$\phi_{\mathrm{G}\mathrm{S}}(x)=\arctan(\frac{v\sqrt{2-v^{2}}}{\exp[\sqrt{2-v^{2}}x]+v^{2}-1})$

(38)

,

.

(39)

Patching together pieces of gray solitons, Hakim found the following 4-indexed
stationary solutions of Eqs. (34) and (35), including the potential term
$g\delta(x)$$R_{\xi}(x)=R_{\mathrm{G}\mathrm{S}}(x\pm\xi)$

x $<0>$,$\phi_{\xi}(x)=\phi_{\mathrm{G}\mathrm{S}}(x\pm\xi)-\phi_{\mathrm{G}\mathrm{S}}(\pm\xi)$

,

(40)
(41)

$x><0$

where the jump conditions (36) and (37) impose the relation$g( \xi)=\sqrt{2}(1-v^{2}/2)^{3/2}\frac{\tanh[\sqrt{1/2-v^{2}/4}\xi]}{v^{2}/2+\sinh^{2}[\sqrt{1/2-v^{2}/4}\xi]}$

The function $g(\xi)$ reaches amaximum $g_{c}=g(\xi_{c})$ at$g_{c}=4(1-v^{2}/2) \frac{[\sqrt{1+4v^{2}}-(1+v^{2})]^{1/2}}{2v^{2}-1+\sqrt{1+4v^{2}}}$

.

(42)

$\xi_{c}=\frac{\arg\cosh(\frac{1+\sqrt{1+4v^{2}}}{2}}{\sqrt{2-v^{2}}}$

.

with
(43)

The two stationary solutions of (32) corresponding to $\xi_{+}(g)>\xi_{c}$ and $\xi_{-}(g)<$obtained by inverting (42) for $g<g_{c}$ thus disappear, merging in asaddlenode bifurcation at acritical strength . Note that the bifurcation can also be
obtained by varying and keeping constant. In the following, the strengthof the delta function is used as the control parameter of our system, keeping
constant.

$\xi_{c}$$g_{c}$

$v$

$g$

$g$

$v$

(see Eq. (28)) is
The bifurcation diagram corresponding to the energy
shown on fig.8, The energetically unstable and stable solutions, as deare also displayed on the figure. Note that the phasefined in eq. (41), differs from that considered in [34] by an (x-independent)
constant. The phase in [34] is set to 0at $x=+\infty$ , whereas (41) is antisymmetric in . This difference is unimportant because Eqs. (34) and (35) are$\mathcal{K}$$(\mathcal{K}(\xi_{-}(g))>$$\phi_{\xi}(x)$$\mathcal{K}(\xi_{+}(g)))$

$x$

18

194
(b)

(a)

$R$

$\phi$of the stable (–) and unstable (—) stationary solutions ofof theeq.(32) (see eq.(40)) for $g=1.250$ and $v=0.5$ ;insert, energy functional
stationary solutions versus for $v=0.5$ (see eq.(28)); lower branch: energeticallystable branch, upper branch: energetically unstable branch. The bifurcation occurs
Phase of the stable (–) and unstable (—) stationary solutions
at $g=1.5514$(see eq.(41)), same conditions as in (a).
Fig. 8. (a) Modulus$R$$\mathcal{K}$$g$$\phi$$(\mathrm{b})$invariant under the constant phase shift$\phi(x)\mapsto*\phi(x)+\varphi$

4.2

.

(44)

General formulation

In this section we define and test the numerical tools needed to obtain thestationary solutions of the NLSE.

Consider the following action functional associated to the NLSE
$A= \int d\tilde{t}\{\int d\tilde{\mathrm{x}}\frac{i}{2}(\overline{\psi}\frac{\partial\psi}{\partial\tilde{t}}-\psi\frac{\partial\overline{\psi}}{\partial\tilde{t}})-\mathcal{F}\}$

(45)

,

is the energy of the system.
and correspond to adequately nondimensionalized space and time
variables respectively.
where
Here,

$\psi$is acomplex field,$\overline{\psi}$

its conjugate and$\mathcal{F}$$\tilde{t}$$\mathrm{x}$The Euler-Lagrange equation corresponding to (45) leads to the NLSE interms of the functional$\mathcal{F}$$\frac{\partial\psi}{\partial\tilde{t}}=-\dot{\iota}\frac{\delta \mathcal{F}}{\delta\overline{\psi}}$

(46)

.

.
as astationary solutions if
This equation obviously admits
lookingIn
weare
general,(46)
extremaofare
Thus, stationary solutions ofT.
.
for an extremum of an energy functional under some constraint
The usual Lagrange multiplier trick consists in introducing acontrol parameter$\psi_{S}$$\delta \mathcal{F}/\delta\psi|_{\psi=\psi_{S}}=0$$\mathcal{E}$19$Q[\psi]=\mathrm{c}st$195$7/\mathrm{a}\mathrm{n}\mathrm{d}$, rather than solving for extrema of

new functional$\ovalbox{\tt\small REJECT} \mathrm{F}[\mathrm{e}]\ovalbox{\tt\small REJECT} \mathrm{f}[\mathrm{e}]$$\frac{\delta \mathcal{F}}{\delta\psi}|_{\nu=\mathrm{c}\mathrm{s}\mathrm{t}}.=0$$\ovalbox{\tt\small REJECT} 5[\mathrm{e}]$, searching for extrema of the–&(2[e]. We thus solve for

.

(47)

We now turn to the precise definitions, corresponding to the two systems
considered in this section :(a) Bose-Einstein Condensates and (b) Superflows.

4.2. 1Superflows
In the problem of asuperflow past an obstacle, is the hydrodynamic energy
andis the flow velocity with respect to the obstacle $[12,14]$ . This implies
thatis the flow momentum. Functional, and are given by theexpressions$\mathcal{E}$$\nu\equiv\vec{U}$$Q$$\equiv\vec{P}$$\mathcal{F}$$\mathcal{E}$$\vec{P}$$\mathrm{s}$(48)$\mathcal{F}=\mathcal{E}-P^{\prec}\cdot U^{\prec}$(49)$\mathcal{E}=c^{2}\int d^{3}x([-1+V(\vec{x})]|\psi|^{2}+\frac{1}{2}|\psi|^{4}+\xi^{2}|\nabla\psi|^{2})$$\vec{P}=\sqrt{2}c\xi\int d^{3}x\frac{i}{2}(\psi\nabla\overline{\psi}-\overline{\psi}\nabla\psi)$

.

(50)

Here, and 4are the physical parameters characterizing the superfluid. They
correspond to the speed of sound (c) for afluid with mean density $\rho_{0}=1$ ,and to the coherence length (4). The potential
is used to represent a$D$cylindrical obstacle of diameter . The NLSE reads$c$$V(x)\prec$.

$\frac{\partial\psi}{\partial t}=-\frac{i}{\sqrt{2}c\xi}\frac{\delta \mathcal{F}}{\delta\overline{\psi}}=i\frac{c}{\sqrt{2}\xi}([1-V(\vec{x})]\psi-|\psi|^{2}\psi+\xi^{2}\nabla^{2}\psi)+\vec{U}\cdot\nabla\psi$(51)

We will be interested in the solutions of, for agiven value of. According to equation (47), these solutions are extrema of at constantmomentum .
$\delta \mathcal{F}/\delta\overline{\psi}=0$$\vec{U}$$\mathcal{E}$$\vec{P}$4.2.2 Bose-Einstein Condensates
We consider acondensate of $N$ particles of mass $m$ and effective scattering
length in aradial confining harmonic potential $V(r)=m\omega^{2}r^{2}/2[15]$ . Quantities are rescaled by the natural quantum harmonic oscillator units of timeand length, thus obtaining the nondimensionalized
variablesand. The control parameter becomes in this context the chemical potential . The total number of particles$a$$\tau_{0}=1/\omega$$L_{0}=\sqrt{\hslash/m\omega}$$\tilde{t}=t/\tau_{0},\tilde{\mathrm{x}}=\mathrm{x}/L_{0}$$\tilde{a}=4\pi a/L_{0}$

$\mu$

20

$\nu$

196
in the condensate is therefore given by q
given, in terms of rescaled variables, by. Functionals ye and$\ovalbox{\tt\small REJECT}/\mathrm{V}^{\ovalbox{\tt\small REJECT}}$$\ovalbox{\tt\small REJECT}$$\ovalbox{\tt\small REJECT} \mathrm{V}$

are

(52)

$\mathcal{F}=\mathcal{E}-\mu \mathcal{N}$(53)$\mathcal{E}=\int d^{3}\tilde{\mathrm{x}}(\frac{1}{2}|\nabla_{\overline{\mathrm{x}}}\psi|^{2}+V(\tilde{\mathrm{x}})|\psi|^{2}+\frac{\tilde{a}}{2}|\psi|^{4})$$N= \int d^{3}\tilde{\mathrm{x}}|\psi|^{2}$

.

(54)

Two different situations are possible, depending on the sign of the (rescaled)effective scattering length . When is positive the particles interact repulsively. Anegative corresponds to an attractive interaction. The dynamical
equation is$\tilde{a}$$\tilde{a}$$\tilde{a}$$\frac{\partial\psi}{\partial\tilde{t}}=-i\frac{\delta \mathcal{F}}{\delta\overline{\psi}}=i[\frac{1}{2}\nabla\frac{2}{\mathrm{x}}\psi-\frac{1}{2}|\tilde{\mathrm{x}}|^{2}\psi-(\tilde{a}|\psi|^{2}-\mu)\psi]$

.

(55)

, for agiven value of
We will be interested in the solutions of. According to equation (47), these solutions are extrema of at constantparticle number $N$ .
$\delta \mathcal{F}/\delta\overline{\psi}=0$$\mathcal{E}$

$\mu$

4.3

Branch following methods

is alocal minimum, the stationary solution
When the extremum ofof (51) can be reached by arelaxation method. If the extremum is not a.
minimum, Newton’s iterative method is used to solve for$\psi_{S}$$\mathcal{F}$$\psi s$4.3.1

Relasation method

In what remains of this section, we will write the NLSE under the following
generic form, which is valid for both the Bose-Einstein condensates and thesuperflow past an obstacle:$\frac{\partial\psi}{\partial t}=-i\frac{\delta \mathcal{F}}{\delta\overline{\psi}}=i(\alpha\nabla^{2}\psi+[\Omega-V(\vec{x})]\psi-\beta|\psi|^{2}\psi)+\vec{U}\cdot\nabla\psi$

.

(56)

of
When the extremum of $T$ is alocal minimum, the stationary solution(56) can be reached by integrating to relaxation the associated real GinzburgLandau equation (RGLE)$\psi_{S}$$\frac{\partial\psi}{\partial t}=-\frac{\delta \mathcal{F}}{\delta\overline{\psi}}=\alpha\nabla^{2}\psi+[\Omega-V(\vec{x})]\psi$$-\beta|\psi|^{2}\psi-iU^{\prec}\cdot\nabla\psi$

Indeed, (56) and (57) have the same stationary solutions.

21

.

(57)

197

$10^{0}$$\overline{\mathrm{u}^{=}\mathrm{o}_{1}}10^{-10}$$10^{-20}$$10^{-30}$

Fig. 9. Two typical examples of the Newton method convergence towards the solution of equation (60) for the problem of asuperflow past acylinder with $\xi/D=1/10$
discretized into $n=128\cross 64=8190$ collocation points. The error
and afield. The convergence is faster than exponential, as
measure is given by
expected for a Newton method.
$\psi_{(j)}$$\sum_{j=1}^{n}f_{(j)}^{2}(\psi)/n$

In our numerical computations, equation (57) is integrated to convergence byusing the Forward-Euler/Backwards-Euler time stepping scheme$\psi(t+\sigma)=^{-1}[(1-i\sigma\vec{U}\cdot\nabla)+\sigma([\Omega-V(_{X}^{\prec})]-\beta|\psi(t)|^{2})]\psi(t)$

(58)

with
$=[1-\sigma\alpha\nabla^{2}]$

(59)

.

The advantage of this method is that it converges to the stationary solutionof (56) independently of the time step .
$\sigma$

4.3.2

Newton method

We use Newton’s method [37] to find unstable stationary solutions of the
RGLE.
In order to work with awell-conditioned system [38], we search for the fixed
points of (58). These can be found as the roots of,(60)$f(\psi)=^{-1}[(1-i\sigma\vec{U}\cdot\nabla)+\sigma([\Omega-V(\vec{x})]-\beta|\psi(t)|^{2})]\psi(t)-\psi(t)$the value of
where $-1$ was already introduced in equation (58). Callingthe field over the -th collocation point, finding the roots of $f(\psi)$ is equivalent$\psi_{(j)}$$\psi$$\mathrm{j}$

22

198

Fig. 10. Two typical examples of a -conjugate gradient method convergence corresponding to the case shown on figure 9. The convergence of the relative error
achieved for the solution ofis given by, where,.
and$\mathrm{b}\mathrm{i}$$\mathrm{A}\mathrm{x}=\mathrm{b}$$\mathrm{x}$$\mathrm{b}=-f_{(j)}(\psi)$$\mathrm{A}=[df_{(j)}/d\psi_{(k)}]$$|\mathrm{A}\mathrm{x}-\mathrm{b}|/|\mathrm{b}|$$\mathrm{x}=\delta\psi_{(k)}$iterating the Newton step
(61)$\psi_{(j)}=\psi_{(j)}+\delta\psi_{(j)}$up to convergence. Every Newton step (61) requires the solution for$\sum_{k}[\frac{df_{(j)}}{d\psi_{(k)}}]\delta\psi_{(k)}=-f_{(j)}(\psi)$.

$\delta\psi_{(k)}$

of

(62)

This solution is obtained by an iterative -conjugate gradient method (BCGM)[39]. The BCGM uses the direct application of $[df_{(j)}/d\psi(k)]$ over an arbitrary
field to obtain an approximative solution of (62). Note that since the convergence of the time step (58) does not depend on , the roots found throughthis Newton iteration are also independent of . Therefore, becomes afree
parameter that can be used to adjust the pre-conditioning of the system inorder to optimize the convergence of the BCGM [38].
$\mathrm{b}\mathrm{i}$$\varphi$$\sigma$$\sigma$

4.3.3$\sigma$

Implementation

We use standard Fourier pseudo spectral methods [19]. Typical convergencesof the Newton and -conjugate gradient iterations are shown in figures 9and
10.
$\mathrm{b}\mathrm{i}$In the case of the radially symmetric Bose Condensate,$\psi(r,\tilde{t})$is expanded$\mathrm{p}\mathrm{o}11\mathrm{y}\mathrm{n}\mathrm{o}\mathrm{m}\mathrm{i}\mathrm{a}1\hat{\psi}_{N_{R}}\mathrm{i}\mathrm{s}\mathrm{f}\mathrm{i}\mathrm{x}\mathrm{e}\mathrm{d}\mathrm{t}\mathrm{o}\mathrm{s}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{s}\mathrm{f}\mathrm{y}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{b}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{a}\mathrm{a}\mathrm{r}\mathrm{y}\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{d}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\psi(R,\tilde{t})=0\mathrm{a}\mathrm{s}\psi(r,\tilde{t})=\sum_{\mathrm{a}\mathrm{n}\mathrm{d}}n=0\hat{\psi}_{2n}(\tilde{t})T_{2n}(r/R)\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}T_{n}\mathrm{i}\mathrm{s}\mathrm{t}\mathrm{h}\mathrm{e}n- \mathrm{t}\mathrm{h}\mathrm{o}\mathrm{r}\mathrm{d}\mathrm{e}\mathrm{r}\mathrm{C}\mathrm{h}\mathrm{e}\mathrm{b}\mathrm{y}\mathrm{c}\mathrm{h}.\mathrm{e}\mathrm{v}N_{R}/2$

23

199
The time integration of the NLSE is done by using afractional step (OperatorSplitting) method [40].

4.4

Stability of 2D superflow around a cylinder

In this section, following references [12-14], we investigate the stationary stable and unstable (nucleation) solutions of the NLSE describing the superflow
around acylinder, using the numerical methods developed in section 4.3. We
study adisc of diameter $D$ , moving at speed
in atw0-dimensional
superfluid at rest. The NLSE (51) can be mapped into two hydrodynamical
equations by applying Madelung’s transformation $[1,16]$ :$\vec{U}$$\psi$$= \sqrt{\rho}\exp(\frac{i\phi}{\sqrt{2}c\xi})$$(2\mathrm{D})$

.

(63)

The real and imaginary parts of the NLSE produce for afluid of density
velocity$\vec{v}=\nabla\phi-\vec{U}$$\rho$

and
(64)

,

the following equations of motion(65)$\frac{\partial\rho}{\partial t}+\nabla(\rho v)\prec=0$$[ \frac{\partial\phi}{\partial t}-\vec{U}\cdot\nabla\phi]+\frac{1}{2}(\nabla\phi)^{2}+c^{2}[\rho-\Omega(\vec{x})]-c^{2}\xi^{2}\frac{\nabla^{2}\sqrt{\rho}}{\sqrt{\rho}}=0$

.

(66)

In the coordinate system that follows the obstacle, these equations correspond to the continuity equation and to the Bernoulli equation [4] (with asupplementary quantum pressure term
) for an isentropic, compressible and irrotational flow. Note that, in the limit where, thequantum pressure term vanishes and we recover the system of equations describing an Eulerian flow.
$\vec{x}$$c^{2}\xi^{2}\nabla^{2}\sqrt{\rho}/\sqrt{\rho}$$\xi/Darrow \mathrm{O}$4.4.1 Bifurcation diagram and scaling in 2D
In this section, varying the ratio of the coherence length 4to the cylinder
diameter D, we obtain scaling laws in thelimit.
$\xi/Darrow \mathrm{O}$

Bifurcation diagram
We present results for $\xi/D=1/10$ which are representative of all ratios wecomputed. The functionaland energyof the stationary solutions are$\mathcal{F}$$\mathcal{E}$

24

200
shown in Fig. 11 as afunction of the Mach number $(M=|\overline{U}|/c)$ . The stablebranch (a) disappears with the unstable solution (c) at asaddle-node bifurcation when $M=M^{c}\approx 0.4286$ . The energy $T$ has acusp at the bifurcation
point, which is the generic behavior for asaddle-node. There are no stationary solutions beyond this point. When $M\approx 0.4282$ , the unstable symmetricbranch (c) bifurcates at apitchfork to apair of asymmetric branches (6). Their
which is roughly half of thenucleation energy barrier is given by.
barrier for the symmetric branch$(\mathcal{F}_{\mathcal{U}}-\mathcal{F}_{a’})$$(\mathcal{F}_{d}-\mathcal{F}_{a’})$We can relate branches in Fig. 11 to the presence vortices in the solution.
as indicated, solutions are irrotational (When)
the stable branch (a) remains irrotational (Fig.
in Fig. 11). For)
while the unstable branch (b) corresponds to aone vortex solution (Fig.
). The disand the unstable branch (c), to atwo vortex solution (Fig.
tance between the vortices and the obstacle in branches (b) and (c) increases
when $M$ is decreased. Branch (c) is precisely the situation described in [41].
.
Furthermore, the value $M^{c}\approx 0.4286$ is close to the predicted value
with,
in
time
forwardthe
NLSE
shows the result of integrating
Figureas initial condition, aslightly perturbed unstable symmetric stationary state). The perturbation drives the system over the nucleation barrier(Fig.
and cycles it, after the emission of two vortices, back to astationary stable s0lution. This shows that the branch (c) corresponds to hyperbolic fixed points
of NLSE.
$M^{\mathrm{p}\mathrm{f}}\sim 0.405$$M^{\mathrm{p}\mathrm{f}}\leq M\leq M^{c}$$12\mathrm{A}$$M\leq M^{\mathrm{p}\mathrm{f}}$$12\mathrm{B}$$12\mathrm{C}$$\sqrt{2/11}$$12\mathrm{D}$$12\mathrm{C}$show the phase of the field at the surface of the disc $(r=D/2$
for four different flow speeds. In both unstable branches,
discontinuities, adiagnostic of vortex crossing, appear between $M=0.40$ and$M=0.41$ .

Figuresand$12\mathrm{E},\mathrm{F}$$2\pi-$$\mathit{0}\in[0,2\pi])$

Scaling laws
We now characterize the dependence on $\xi/D$ of the main features of the bifurbecome indistinguishandcation diagram. When $\xi/D$ is decreased,
will be thatMachnumber,thecriticalable. In the limit where $\xi/D=0$.
of an Eulerian flow$M^{c}$$M^{\mathrm{p}\mathrm{f}}$$M^{c}$$M_{\mathrm{E}\mathrm{u}1\mathrm{e}\mathrm{r}}^{c}$to the Eulerian critical velocity.
Figure 13 shows the convergence ofThis convergence can be characterized by fitting the polynomial law $M^{c}=$to $M^{c}(\xi/D)$ . This fit is shown on Fig. 13 as adotted line,.
yielding $K_{1}=0.322$ , $K_{2}=0.615$ and$M^{c}$$K_{1}(\xi/D)^{K_{2}}+M_{\mathrm{E}\mathrm{u}1\mathrm{e}\mathrm{r}}^{c}$$M_{\mathrm{E}\mathrm{u}1\mathrm{e}\mathrm{r}}^{c}=0.35$

Dynamical solutions The stationary solutions obtained in the above subsection provide adequate initial data for the study of dynamical solutions.
Indeed, after asmall perturbation, their integration in time will generate adynamical evolution with very small acoustic emission. Therefore, this proce-

25

201

$\mathrm{u}^{\mathrm{c}_{1}}\frac{}{\Phi}\frac{\Phi}{\Phi}$

versus Mach number $(M=|\vec{U}|/c)$ ,, and functional
Fig. 11. Plot of the energy
with $D=10\xi$ . Stable state (a). Nucleation solutions: asymmetric branch (b) andsymmetric branch (c). The diagram shows asaddle-node and apitchfork bifurcation.
The point where vortices cross the surface of the disc (see fig. 12) is gabelled by(see text).
. The total fluid momentum is given by$(\mathcal{E})$$(\mathcal{F})$$-d\mathcal{F}/dU$$M_{\mathrm{n}}$dure corresponds to an efficient way to start vortical dynamics in acontrolledmanner.
Starting from atw0-vortex unstable stationary solution at asupercritical Mach
number $M^{c}=0.9$ , the evolution of the NLSE time integration shows aclearly
periodical emission of vortex pairs (see Fig. 12). This emission conserves total
circulation.

We have studied the behavior of the frequency of vortex emission close to thebifurcation for such symmetrical wakes with different supercritical velocities(characterized by $\delta_{\mathrm{s}\mathrm{p}}=(M-M^{c})/M^{c}=-\delta>0$). Our results, plotted on
scaling.
Fig. 14, are consistent with a$\delta_{\mathrm{s}\mathrm{p}}^{1/2}$

4.4.2

Subcriticality and vortex-stretching in 3D

In this section, using a 3D version of our code to integrate the NLSE, we study
3D instabilities of the basic 2D superflow.
Preparation methody) (corresponding to branch
We used the 2D laminar stationary solution,(a) of precedent section) and the one-vortex unstable stationary solution$\psi_{0V}(x,$$\psi_{1V}(x$

26

202

.

.

$2\mathrm{E}_{*^{\mu}4}\mathrm{r}2_{|,|}^{k}|,--$$\mathrm{r}\alpha:$$0|$$2\mathrm{F}_{*}\dot{L}0,\lceil_{1}^{--_{n2\prime}}\mathrm{a}- \mathrm{r}^{\underline{\triangleleft}}\mathrm{r}2--’- n20^{1\mathrm{I}}|$$\}|$$- \mathrm{r}\iota_{\mathrm{I}}- \mathrm{z}_{\mathrm{I}}^{1}||$$-|*$$-0?|$

0$*|$Fig. 12. Stationary states :stable (A), one vortex unstable (B), two vortices unstable (C). The surface indicates the fluid density around the cylinder$(M=0.24, \xi/d=0.1)$ . (D) Shows the result of the NLSE integration, startingfrom aslightly perturbed stationary (C) state. Figures (E) and (F) display thephase of the complex field at the surface of the cylinder versus the polar angle
0. Asymmetric branch (A), symmetric branch (B). $M=0.4286$, $M=0.41(\square )$ ,$M=0.40(+)$ , $M=0.30(\mathrm{x})$ . The crossing out of the vortex produces aphase
discontinuity at.
$\psi$$(\circ)$$M^{\mathrm{p}\mathrm{f}}\sim 0.405$(branch (b)) to construct the 3D initial condition$\psi_{3\mathrm{D}}(x,$y,$z)=f_{1}(z)\psi_{1V}(x, y)+[1-f_{1}(z)]\psi_{0V}(x,$ y).

The function$f_{1}(z)$(67), defined by,$f_{i}(z)=(\tanh[(z-z_{1})/\Delta_{z}]-\tanh[(z-\mathrm{z}2)/\mathrm{A}\mathrm{z}])/2$

takes the value 1for $z_{1}\leq z\leq z_{2}$ and 0elsewhere, withlength.

27$\Delta_{z}$an adaptation

203
05 ——————-$+\cross$

0.45

$\Xi$

$\dot{\mathrm{M}_{1}}------------------*_{----------------------}$$0.4\dot{\mathrm{q}}-----------------------------------------**$$\mathrm{x}$

0.350

0.06

0.15

Q.i

0.2

$\xi/_{\mathrm{d}}$Fig. 13. Saddle-node bifurcation Mach number $M^{c}(+)$ and pitchfork bifurcation
Mach number, as afunction of $\xi/D$ . The dotted curve corresponds to a
fit to the polynomial law
with $K_{1}=0.322$ , $K_{2}=0.615$
and. The dashed lines
andcorrespondrespectively to first and second order compressible corrections to the $M^{c}=0.5$
critical velocity computed using alocal sonic criterion for an incompressible flow$M^{\mathrm{p}\mathrm{f}}(\cross)$$M^{c}=K_{1}(\xi/D)^{K_{2}}+M_{\mathrm{E}\mathrm{u}1\mathrm{e}\mathrm{r}}^{c}$$M_{1}^{*}\approx 0.4264$$M_{\mathrm{E}\mathrm{u}1\mathrm{e}\mathrm{r}}^{c}=\mathrm{Q}.35$$M_{2}^{*}\approx 0.3903$(see text).
0.024\prime\prime$\prime\prime\prime$$,\acute{\mathrm{o}}$

’

0.020

$’\phi’$”$’\sigma’$

0.016

”
’

$’\varphi’$$>$\prime\prime\prime

”

0.012

”
$\prime\prime\prime\phi$”$’\varphi’\prime\prime\prime\prime\prime s$

0.008
0.02

0.04

6

0.06

0.08

$\mathrm{s}\mathrm{p}$Fig. 14. Vortex emission frequency as afunction of(At $-M^{c}$ ) $/M^{c}\ll 1$ (with$M^{c}=0.3817)$ , for asymmetric wake and $\xi/D=1/20$ . The dashed line shows afit ofapolynomialwith $K_{1}=0.081$ . The obtainedlaw for the frequency
is equivalent to the one expected for adissipative system.
$\delta_{\mathrm{s}\mathrm{p}}=$$\nu=K_{1}\delta_{\mathrm{s}\mathrm{p}}^{1/2}$$\delta_{\mathrm{s}\mathrm{p}}^{1/2}$

Figure 15 represents a
initial data prepared with this method for $\xi/D=$
0.025, $|U^{\prec}|/c=0.26$ andin theperiodicity box$L_{x}/D=2.4\sqrt{2}\pi$$L_{z}/D=0.4\sqrt{2}\pi$(, $L_{y}/D=1.2\sqrt{2}\pi$ and). The surface$0.5$ draws the cylinder surface and the initial condition vortex line,with bothends pinned to the right side of the cylinder.
$3\mathrm{D}$$\triangle_{z}=2\sqrt{2}\xi$$[L_{x}\cross L_{y}\cross L_{z}]$$|\psi_{3\mathrm{D}}|=$

28

204

Fig. 15. Initial condition of avortex pinned to the cylinder generated by eq.(67).
is shown for $\xi/D=0.025$ , $|\vec{U}|/c=0.26$ and
The surface
periodicity box ( $L_{x}/D=2.4\sqrt{2}\pi$ , $L_{y}/D=1.2\sqrt{2}\pi$ andin the$\Delta_{z}=2\sqrt{2}\xi$$|\psi_{3\mathrm{D}}|=0.5$$[L_{x}\cross L_{y}\cross L_{z}]$$L_{z}/D=0.4\sqrt{2}\pi)$.

Short time dynamics
Starting from the initial condition (67), the evolution of the NLSE time integration shows ashort-time and along-time dynamics.
During the short-time dynamics, the initial pinned vortex line rapidly contracts, evolving through adecreasing number of half-ring-like loops, down to,). This evolution,
asingle quasi-stationary half-ring (see Figs.
flow,
providedthat the inito
thehappens mainly on the plane perpendiculartial vortex is long enough to contract to aquasi-stationary half-ring as shown. Otherwise, the vortex line collapses against the cylinder while
on Fig.
moving upstream.
$16\mathrm{a}$$16\mathrm{b}$$16\mathrm{c}$$16\mathrm{c}$Note that this quasi-stationary half-ring has been used by Varoquaux [42,43]
to estimate the nucleation barrier in a 3D experiment.
) is very slow and can be
The dynamics of the half-ring situation (Fig.
Indeed,local flow velocity in an
thefield.
astationary
to
to
be
close
shown
Eulerian flow around acylindrical obstacle is known to vary from $v=|\vec{U}|$ at
infinity to $v=2|\vec{U}|$ at both sides of its surface. Moreover, the diameter ofastationary vortex ring in an infinite Eulerian flow with no obstacle is givenby [1]:$16\mathrm{c}$

$v$

$d$

$|\vec{U}|/c=(\sqrt{2}\xi/d)[\ln(4d/\xi)-K]$

(68)

,

is the flow velocity at infinity and the vortex core model constant$K\sim 1$ is obtained by fitting the numerical results in [44]. Therefore, for thevalues used on Figs. 16, we expect that local velocities range from $v=0.25$ to$v=2\cross 0.25$ . Equation (68) thus implies that the diameter of an hypotheticalstationary half-ring should be bounded by $d(v=|\vec{U}|=0.25)=18.8\xi$ and$d(v=2|\vec{U}|=0.5)=6.3\xi$ . The diametermeasured on the half-ringis consistent with its quasi-stationary behavior. Similarly
observed on Fig.
is also found to be
the diameter of the half-ring shown on Fig. 18
between the corresponding bounds $d(0.35)=11.4\xi$ and $d(2\cross 0.35)=3\xi$ .
where$|\vec{U}|$$d\approx 9\xi$$16\mathrm{c}$$d\approx 7.6\mathrm{f}$

29

205

Fig. 16. Short-time dynamics for $\xi/D=1/40$ and $|\vec{U}|/c=0.25$ starting from Fig. 15:
A $(t=5\xi/c)$ , $\mathrm{B}(t=10\xi/c)$ and $\mathrm{C}(t=15\xi/c)$ . The contraction of the initial vortex
line occurs in the plane perpendicular to the flow. The half-rings have adiametercompatible with that of aquasi-stationary half-ring (see text).

Vortex stretching as asubcritical drag mechanism Asmall perturbation
over the half-ring solution can drive the system into two opposite situations
where the half-ring either starts moving upstream or downstream.

When driven upstream, the half-ring eventually collapses against the cylinder,
dissipating its energy as sound waves. Otherwise, the vortex loop is stretchedwhile the pinning points move towards the back of the cylinder. Figures 17
show the long-time dynamics for astretching case with $\xi/D=1/40$ and$|U|\prec/c=0.25$ starting from Fig.
. Figure 19 shows alater situation for$\xi/D=1/20$ and $|\vec{U}|/c=0.35$ starting from Fig. 18. As the vortex loop grows,
its backmost part remains oblique to the flow. The described vortex stretchingmechanism consumes energy, thus generating drag. It can be responsible forthe appearance of drag in experimental superflows if fluctuations are strong
enough to nucleate the initial vortex loop (which is imposed extrinsically insubcritical velocities.
our numerical system). Note that it takes place for$16\mathrm{c}$$2\mathrm{D}$

30

206

Fig. 17. Long-time dynamics for $\xi/D=1/40$ and $|\vec{U}|/c=0.25$ starting from
Fig. 16c. The half-ring moves downstream while growing.

Fig. 18. Quasi-nucleation solution for $|\vec{U}|/c=0.35$ and $\xi/D=1/20$ at time
t $=15\xi/c$ .

31

207

Fig. 19. Vortex stretching at $t=150\xi/c$ with $|U^{\prec}|/c=0.35$ and $\xi/D=1/20$ . The
vortex line is oblique to the flow.
$10^{0}$$\cross_{\mathrm{R}_{\mathrm{I}\mathrm{Q}}}$$10^{-2}$$\mathrm{o}$$\circ$$\mathrm{o}$$> 0^{-4}0\tilde{\mathrm{o}}$$\mathrm{o}_{\mathrm{o}\mathrm{q}}\mathrm{o}$8$\{0^{-6}$$0_{8}$$10^{-8}$$0^{0}$$10^{2}$$10^{6}$$10^{8}$$\{0_{\mathrm{d}/\xi}^{4}$

Fig. 20. Critical Mach number $Vc/C$ versus scale ratio of numerical and experimental
data $D/\xi$ . Circles correspond to several experiments from [45]. Squares stand forour numerical stretching cases while crosses correspond to non-stretching cases [14].

Figure 20 displays several numerical and experimental [45] critical Mach numbers $(V_{c}/C)$ with respect to $D/\xi$ , which seem to follow a(-1) slope in a
log-log plot. The squares stand for our numerical stretching cases while thecrosses correspond to non-stretching cases. There is afrontier between the
numerical dissipative and non-dissipative cases [14]. For $1/30<\xi/D<1/20$ ,$3\mathrm{D}$

32

208
the frontier corresponds to the expression $R_{s}=5.5$ with
$R_{s}\equiv|\tilde{U}|D/c\xi=MD/\xi$

(69)

.

This superfluid ‘Reynolds’ number is defined in the same way as the standard(viscous) Reynolds number $Re\equiv|\tilde{U}|D/\nu$ (with the kinematic viscosity).
It has been shown, in the superfluid turbulent $(R_{s}\gg 1)$ regime, thatis equivalent to the standard (viscous) Reynolds number $Re[10,9,11]$ . Note
that, for aBose condensate of particles of mass $m$ , the quantum of velocity
circulation around avortex, $\Gamma=2\pi\sqrt{2}c\xi$ , has the Onsager-Feynman valueas$\Gamma=h/m$ (is Planck’s constant) and the same physical dimensions.
$\nu$$R_{s}$$L^{2}T^{-1}$

$h$

$\nu$

divides the space of parameters into alaminar flow zone
The value ofand arecirculating flow zone, very much like in the problem of acircular
disc in aviscous fluid in which this frontier is also found to be around $Re\sim$. It therefore seems to exist some degree of universality between viscous
normal fluids and superfluids modeled by NLSE as discussed in [10,9,11]. In
flow, the experimental critical velocity is known
the context of superfluid
to depend strongly on the system’s characteristic size $D$ . It is often found to bewell below the Landau value (based on the velocity of roton excitation) except
for experiments where ions are dragged in liquid helium. Feynman’s alternative
critical velocity criterion $Rs\sim\log(D/\xi)$ is based on the energy needed to form
vortex lines. It produces better estimates for various experimental settings, butdoes not describe the vortex nucleation mechanism [1].
$R_{s}$$5$$4\mathrm{H}\mathrm{e}$In arecent experiment, Raman et al. have studied dissipation in aBoseEinstein condensed gas by moving ablue detuned laser beam through thecondensate at different velocities [46]. In their inhomogeneous condensate,.
they observed acritical Mach number for the onset of dissipation$M_{2D}^{c}/1.6$

Our computations were performed for values of $\xi/D$ comparable to those in
Bose-Einstein condensed gas experiments. They demonstrate the possibility ofvortex stretching. It would be very
asubcritical drag mechanism, based on$3\mathrm{D}$interesting to determine experimentally the dependence of the critical Mach
) of the excitationsornumber on the parameter $\xi/D$ and the nature ([14].
$2\mathrm{D}$4.5$3\mathrm{D}$Stability of attractive Bose-Einstein condensates

In this section, following reference [15], we study condensates with attractive
interactions which are known to be metastable in spatially localized systems,
provided that the number of condensed particles is below acritical value[7]. Various physical processes compete to determine the lifetime of attractive$N_{c}$

33

209
condensates. Among them one can distinguish macroscopic quantum tunneling (MQT) [47,48], inelastic two and three body collisions (ICO) [49,50] andthermally induced collapse (TIC) [48,51]. We compute the life-times, using
both avariational Gaussian approximation and the exact numerical solution
for the condensate wave-function

4.5.1

Computations of stationary states

Gaussian approximation

AGaussian approximation for the condensate density can be obtained analytically through the following procedure.
Inserting(70)$\psi(r,\tilde{t})=A(\tilde{t})\exp(-r^{2}/2r_{\mathrm{G}}^{2}(\tilde{t})+ib(\tilde{t})r^{2})$into the action (45), where is given by (52), yields aset of Euler-Lagrange. The stationary,
equations forand the (complex) amplitude
solutions of the Euler-Lagrange equations produce the following values [52]:$\mathcal{F}$$r_{\mathrm{G}}(\tilde{t})$$A(\tilde{t})$$b(\tilde{t})$$N( \mu)=\frac{4\sqrt{2\pi^{3}}(-8\mu+3\sqrt{7+4\mu^{2}})}{7|\tilde{a}|(-2\mu+\sqrt{7+4\mu^{2}})^{3/2}},$$\mathcal{E}=N(\mu)(-\mu+3\sqrt{7+4\mu^{2}})/7$

N is found to be maximal atof the chemical potential is

(71)

,

(72)

.

$N_{c}^{G}=8\sqrt{2\pi^{3}}/|5^{5/4}\tilde{a}|$$\mu=\mu_{c}^{G}=1/2\sqrt{5}$. The corresponding value.

Linearizing the Euler-Lagrange equations around the stationary solutions,yields the following expression for the eigenvalues [52]:(73)$\lambda^{2}(\mu)=8\mu^{2}-4\mu\sqrt{7+4\mu^{2}},+2$This qualitative behavior is the generic signature of aHamiltonian Saddle
Node (HSN) bifurcation defined, at lowest order, by the normal form [53]$m_{\mathrm{e}ff}\dot{Q}=\delta-\beta Q^{2}$

(74)

,

where $\delta=(1-N/N_{c})$ is the bifurcation parameter. The critical amplitude
is related to the radius of the condensate [52]. We can relate the parametersand meff to critical scaling laws, by defining the appropriate energy$Q$$\beta$.

$\mathcal{E}=\mathcal{E}_{0}+m_{\mathrm{e}fJ}\dot{Q}^{2}/2-\delta Q+\beta Q^{3}/3-\gamma\delta$

34

(75)

210

Fig. 21. Stationary solutions of the NLSE equation versus particle number $N$ . Left:value of the energy functionalon the stable (elliptic) branch andon theunstable (hyperbolic) branch. Right: square of the bifurcating eigenvalue
is the energy of small excitations around the stable branch. Solid lines: exact solutionof the NLSE equation. Dashed lines: Gaussian approximation.
$\mathcal{E}_{-}$$\mathcal{E}_{+}$$(\lambda_{\pm}^{2});|\lambda_{-}|$

From (74) it is straightforward to derive, close to the critical point
universal scaling laws$\mathcal{E}_{\pm}=\mathcal{E}_{c}-\mathcal{E}_{l}\delta\pm \mathcal{E}_{\Delta}\delta^{3/2}$$\lambda_{\pm}^{2}=\pm\lambda_{\Delta}^{2}\delta^{1/2}$

where$\mathcal{E}_{c}=\mathcal{E}_{0}$,$\mathcal{E}_{l}=\gamma$,,$\delta$$=0$ ,

the

(76)

,

$\mathcal{E}_{\Delta}=2/3\sqrt{\beta}$

(77)

and

$\lambda_{\Delta}^{2}=2\sqrt{\beta}/m_{\mathrm{e}[f}$.

Numerical branch following

Using the branch-following method described in section 4.3, we have computed
the exact stationary solutions of the NLSE. We use the following value$-5.74\cross 10^{-3}$ ,that corresponds to experiments withatoms in aradial trap$[54,7]$ .
$\tilde{a}=$$7\mathrm{L}\mathrm{i}$As apparent on Fig. 3, the exact critical $N_{c}^{E}=1258.5$ is smaller than theGaussian one $N_{c}^{G}=1467.7$ $[55,47]$ . The critical amplitudes corresponding
to the Gaussian approximation can be computed from (71) and (72). One
,
finds
and
. For the exact
solutions, we obtain the critical amplitudes by performing fits on the data.
One finds
and. Thus, the Gaussian approximation
captures the bifurcation qualitatively, but with quantitative 17% error on[55], 24% error on
and 14% error on. Fig. 4shows the physical originof the quantitative errors in the Gaussian approximation. By inspection it isapparent that the exact solution is well approximated by aGaussian only forsmall Aon the stable (elliptic) branch.
$\mathcal{E}_{c}=4\sqrt{2\pi^{3}}/|5^{3/4}\tilde{a}|$$\mathcal{E}_{\Delta}=1340$$\mathcal{E}_{\Delta}=64\sqrt{\pi^{3}}/|5^{9/4}\tilde{a}|$$\lambda_{\Delta}^{2}=4\sqrt{10}$$\lambda_{\Delta}^{2}=14.68$$N_{c}$$\mathcal{E}_{\Delta}$$\lambda_{\Delta}^{2}$

35

211

20
$|\psi|^{2}$

10

0
$r$

versus radius , in reduced units (see text). Solid
Fig. 22. Condensate density
lines: exact solution of the NLSE equation. Dashed lines: Gaussian approximation.
$N$
Stable (elliptic) solutions are shown for particle number $N$ $=252(\mathrm{a})$ and $=1132$$N$$=1132$ (see insert).
(b). (c) is the unstable (hyperbolic) solution for$r$$|\psi|^{2}$

Estimation of life-time

4.5.2

$s$

In this section, we estimate the decay rates due to thermally induced collapse,
macroscopic quantum tunneling and inelastic collisions.

Thermally induced collapse

The thermally induced collapse (TIC) rate
[56]$\Gamma_{T}$is estimated using the formula(78)$\frac{\Gamma_{T}}{\omega}=\frac{|\lambda_{+}|}{2\pi}\exp[\frac{-\hslash\omega(\mathcal{E}_{+}-\mathcal{E}_{-})}{k_{B}T}]$

is the (dimensionalized) height of the nucleation energy barwhereis the Boltzmann constant.
rier, is the temperature of the condensate andtime which is controlled
decay
Note that the prefactor characterizes the typicalby the slowest part of the nucleation dynamics: the top-0f-the-barrier saddle
can be obtained directly from the. The behavior ofpoint eigenvalue
universal saddle-node scaling laws (76) and (77). Thus the exponential factor.
andand the prefactor vanish respectively as$\hslash\omega(\mathcal{E}_{+}-\mathcal{E}_{-})$$k_{B}$$T$$\lambda_{+}$$\Gamma_{T}$$\delta^{3/2}$$\delta^{1/4}$

Macroscopic quantum tunneling

We estimate the MQT decay rate using an instanton technique that takes into
account the semi classical trajectory giving the dominant contribution to the

36

212

$\mathrm{V}(\mathrm{q}_{\mathrm{f}})$$\mathrm{V}(\mathrm{q}_{\mathrm{n}}\}$

Fig. 23. The bounce trajectory is shown as dashes, above the potential $V(q)$ .

quantum action path integral [48,47]. This trajectory is approximated as thesolution of$\frac{d^{2}q(t)}{dt^{2}}=-\frac{dV(q)}{dq}$

,

(79)

is a polynomial such that $-V(q)$ reconstructs the Hamiltonian dynamics.
is determined by the relations$V(q)$$V(q)$

(80)
(81)
(82)

$V(q_{m})=-\mathcal{E}_{+}$$V(q_{f})=-\mathcal{E}_{-}$$\partial_{q}^{2}V(q_{m})=|\lambda_{+}(N)|$$\partial_{q}^{2}V(q_{f})=-|\lambda_{-}(N)|$

.

(83)

The bounce trajectory is displayed on Fig. 23 (dashed line) above the potential$V(q)$ . The MQT rate is estimated as$\frac{\Gamma_{Q}}{\omega}=\sqrt{\frac{|\lambda_{-}|v_{0}^{2}}{4\pi}}\exp[\frac{-4}{\sqrt{2}}\int_{q_{f}}^{qb}\sqrt{V(q)-V(q_{f})}dq]$

where

,

(84)

is defined by the asymptotic form of the bounce trajectory $q(t)[48]$ :. Universal scaling laws can be derived close
to criticality from (74), (76) and (77). The exponential factor in (84) follows
the same scaling than
. It therefore vanishes as. From theasymptotic form of $q(t)$ , $dq$ follows the same law as. Thusand the prefactor vanishes as.
$v_{0}$$q(\tau)\sim qf+(v_{0}/|\lambda_{-}|)\exp[-|\lambda_{-}\tau|]$$\sqrt{|\mathcal{E}_{+}-\mathcal{E}_{-}|}dq$$\delta^{5/4}$$v_{0}/|\lambda_{-}|$$\delta^{7/8}$

Inelastic collision37$v_{0}\sim\delta^{3/4}$

213–$\underline{\omega|}$$\frac{\mathrm{q})}{oe\varpi}$$[mathring]_{\mathrm{O}}0\iota \mathfrak{o}^{\backslash }>$,$\frac{\Phi}{\varpi,\{\emptyset\subset}$$\mathfrak{D}q)$$\mathrm{O}\mathrm{o}\subset$

Fig. 24. Condensate decay rates versus particle number. ICO: inelastic collisions.
MQT: macroscopic quantum tunneling TIC: thermally induced collapse at temperatures $1nK(1),2nK(2)$ , $50nK(3)$ , 100nK (4), $200nK(5)$ , $300nK(6)$ and $400nK$(7). The insert shows the details of the cross-0ver region between quantum tunnelingand thermal decay rate. Solid lines: exact solution of the NLSE equation, dashedlines: Gaussian approximation.

The inelastic collision rate (ICO) is estimated using the relation(85)$\frac{dN}{dt}=f_{C}(N)$with$f_{C}(N)=K \int|\psi|^{4}d^{3}\tilde{\mathrm{x}}+L\int|\psi|^{6}d^{3}\tilde{\mathrm{x}}$

(86)

,

where $K=3.8\cross 10^{-4}\mathrm{s}^{-1}$ and $L=2.6\cross 10^{-7}\mathrm{s}^{-1}$ . The ICO rate can beevaluated from the stable branch alone. In order to compare the particle decayrate $f_{C}(N)$ to the condensate collective decay rates obtained for TIC and
MQT, we compute the condensate ICO half-life as:(87)$\tau_{1/2}(N)=\int_{N/2}^{N}dn/f_{C}(n)$and plot$\tau_{1/2}^{-1}$

Discussion
It is apparent by inspection of Fig. 6that for agiven value of $N$ the exact andGaussian approximate rates are dramatically different. We now compare thethe MQT
relative importance of the different exact decay rates. At$T\leq 1\mathrm{n}\mathrm{K}$

38

214
effect becomes important compared to the ICO decay in aregion very closetoas it was shown in [47] using Gaussian computations
but evaluating them with the exact maximal number of condensed particles. Considering thermal fluctuations for temperatures as low as 2
, it isapparent on Fig. 6(see insert) that the MQT will be the dominant decay
mechanism only in aregion extremely close to( $\delta<5$ x.1O -3) where thecondensates will live less than 10 . Thus, in the experimental case of
atoms, the relevant effects are ICO and TIC, with cross-0ver determined inFig. 6.
$N_{c}^{E}(\delta\leq 8\cross 10^{-3})$$N_{c}^{E}$$\mathrm{n}\mathrm{K}$$N_{c}$$7\mathrm{L}\mathrm{i}$$\mathrm{s}$5Conclusion
The main result of the NLSE simulations presented in section 3.2 is thattwo diagnostics of Kolmogorov’s regime in decaying turbulence are satisfied.
These diagnostics are, at the time of the maximum of energy dissipation: (i) aparameter-independent kinetic energy dissipation rate and (ii) a
spectralscaling in the inertial range. Thus, the NLSE simulations were shown to bevery similar as far as energetics is concerned, with the viscous simulations. The
experimental results shown in section 3.3 show that the Kolmogorov cascadesurvives in the superfluid regime.
$k^{-5/3}$

We have seen that the numerical tools developed in section 4.3 can be used inpractice to obtain the stationary solutions of the NLSE. These methods haveallowed us to find the full bifurcation diagrams of Bose-Einstein condensateswith attractive interactions and superflows past acylinder. Furthermore, thestationary solutions have given us efficient way to start vortical dynamics (in) in acontrolled manner.
and$2\mathrm{D}$$3\mathrm{D}$6Acknowledgments

The work reviewed in this paper was performed in collaboration with M.
Abid, G. Dewel, C. Huepe, J. Maurer, S. Metens, C. Nore, -T. Pham and P.
Tabeling. Computations were performed at the Institut du Developpement etdes Ressources en Informatique Scientifique.
$\mathrm{C}$

39

215
References[1] R. J. Donnelly. Quantized Vortices in Helium $II$. Cambridge Univ. Press,
Cambridge, 1991.

[2]. Gross. Structure of aquantized vortex in boson systems. Nuovo Cimento,$\mathrm{E}.\mathrm{P}$20(3), 1961.

[3]. Pitaevskii. Vortex lines in an imperfect Bose gas. $Sov$ . Phys.-JETP, 13(2),1961.

$\mathrm{L}.\mathrm{P}$[4] L. Landau and E. Lifchitz. Fluid Mechanics. Pergamon Press, Oxford, 1980.

[5] M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A.
Cornell. Science, 269:198, 1995.
[6] K. B. Davis, M. O. Mewes, M. R. Adrews, N. J. van Druten, D. S. Durfee,
D. M. Kurn, and W. Ketterle. Phys. Rev. Lett., 75:3969, 1995.

[7] C. C. Bradley, C. A. Sackett, and R. G. Hulet. Bose-einstein cond.e nsation oflithium: Observation limited condensate number. Phys. Rev. Lett., 78(6):985,

1997.

[8] Franco Dalfovo, Stefano Giorgini, Lev P. Pitaevskii, and Sandro Stringari.
Theory of bose-einstein condensation in trapped gases. Reviews of Modern
Physics, 71(3), 1999.
[9] C. Nore, M. Abid, and M. Brachet. Decaying kolmogorov turbulence in amodelof superflow. Phys. Fluids, 9(9):2644, 1997.
[10] C. Nore, M. Abid, and M. E. Brachet. Kolmogorov turbulence in lowtemperature superflows. Phys. Rev. Lett., 78(20):3896-3899,1997.

[11] M. Abid, M. Brachet, J. Maurer, C. Nore, and P. Tabeling. Experimental
and numerical investigations of low-temperature superfluid turbulence. $Eur$ . .
Mech.
Fluids, 17(4):665-675, 1998.
$J$

$B$

[12] C. Huepe and M.-E. Brachet. Solutions de nucleation tourbillonnaires dans
un modele d’ecoulement superfluide. C. R. Acad. Sci. Paris, 325(11):195-202,

1997.

[13] C. Huepe and M. E. Brachet. Scaling laws for vortical nucleation solutions inamodel of superflow. Physica $D$ , 140:126-140, 2000.
[14] C. Nore, C. Huepe, and M. E. Brachet. Subcritical dissipation in threedimensional superflows. Phys. Rev. Lett., 84(10):2191, 2000.
[15] C. Huepe, S. Metens, G. Dewel, P. Borckmans, and M.-E. Brachet. Decay rates
in attractive bose-einstein condensates. Phys. Rev. Lett., 82(2):1999.
[16] E. A. Spiegel. Fluid dynamical form of the linear and nonlinear schrodinger
equations. Physica $D$ , 1:236, 1980.

40

216
[17] P. Nozieres and D. Pines. The Theory of Quantum Liquids. Addison Wesley,
New York, 1990.
[18] C. Nore, M. Brachet, and S. Fauve. Numerical study of hydrodynamics using
the nonlinear schr\"odinger equation. Physica D, 65:154-162, 1993.

[19] D. Gottlieb and S. A. Orszag. Numerical Analysis of Spectral Methods. SIAM,
Philadelphia, 1977.
[20] M. P. Kawatra and R. K. Pathria. Quantized vortices in imperfect bose gas.
Phys. Rev., 151:1, 1966.

[21] K. W. Schwarz. Three-dimensional vortex dynamics in superfluid
and line-boundary interactions. Phys. Rev. , 31:5782, 1985.

$4\mathrm{h}\mathrm{e}$: line-line$B$[22] J. C. Neu. Vortices in complex scalar fields. Physica$D$, 43:385, 1990.

[23] F. Lund. Defect dynamics for the nonlinear schrodinger equation derived fromavariational principle. Phys. Rev.Lett., A159:245, 1991.
[24] G. I. Taylor and A. E. Green. Mechanism of the production of small eddies
from large ones. Proc. Soc. Lond. A, 158:499, 1937.
[25] M. E. Brachet, D. I. Meiron, S. A. Orszag, B. G. Nickel, R. H. Morf, and
U. Frisch. Small-scale structure of the taylor-green vortex. J. Fluid Mech.,$130:411\triangleleft 52$,1983.

[26] M. Brachet. Geometrie des structures \‘a petite \’echelle dans le vortex de taylorgreen. C.R.A.S II, 311:775, 1990.
[27] J. Domaradzki, W. Liu, and M. Brachet. An analysis of sugrid-scale interactions
in numerically simulated isotropic turbulence. Phys. Fluids A, 5:1747, 1993.
[28] C. Nore, M. Abid, and M. Brachet. Simulation num\’erique d’\’ecoulements
cisailles tridimensionnels \‘a Paide de l’\’equation de schr\"odinger non lineaire.
C.R.A. , 31911(7):733, 1994.
$S$[29] U. Frisch. Turbulence, the legacy of A. N. Kolmogorov. Cambridge Univ. Press,
Cambridge, 1995.

Direct observation of the[30] S. Douady, Y. Couder, and M. E. Brachet.
Phys. Rev. Lett,
turbulence.
inintermittency of intense vorticity filaments

67:983, 1991.
[31] G. Zocchi, P. Tabeling, J. Maurer, and H. Willaime. Measurement of the scalingof the dissipation at high reynolds numbers. Phys. Rev. E, 50:3693, 1994.
[32] S. Fauve, C. Laroche, and B. Castaing.
turbulent flows. J. Phys. If,3:271, 1993.

Pressure fluctuations in swirling[33] F. Belin, P. Tabeling, and H. Willaime. Exponents of the structure functions
in alow temperature helium experiment. Physica D, 93:52-63, 1996.

[34] V. Hakim. Nonlinear Schr\"odinger flow past an obstacle in one dimension. Phys.
Rev. E, 55(3):2835-2845, 1997.

41

217
[35] T. Tsuzuki. Nonlinear waves in the Pitaevskii-Gross equation. J. Low Temp.
Phys., $4(4)$ , 1971.
[36]. Shabat. Interaction between soliton in astable. Zakharov andmedium. $Sov$ . Phys.-JETP, 37(5):823-828, 1973.
$\mathrm{A}.\mathrm{B}$$\mathrm{V}.\mathrm{E}$[37] R. Seydel. From Equilibrium to Chaos: Practical Bifurcation and Stability
Analysis. Elsevier, New York, 1988.
[38] C. Mamun and L. Tuckerman. Asymmetry and hopf bifurcation in sphericalcouette flow. Phys. Fluids, 7(1):80, 1995.

[39] W. Press, S. Teukolsky, W. Vetterling, and B.Flannery. Numerical Recipres in
Fortran. Cambridge Univ. Press, Cambridge, 1994.

[40] R. Klein and A. J. Majda. Self-stretching of perturbed vortex filaments. Physica, 53:267, 1991.
$D$[41] T. Prisch, Y. Pomeau, and S. Rica. Transition to dissipation in amodel ofsuperflow. Phys. Rev. Lett., 69:1644, 1992.

[42] O. Avenel,$4\mathrm{h}\mathrm{e}$[43]$\mathrm{G}.\mathrm{G}$. Ihas, and E. Varoquaux. The nucleation of vortices in superfluid:Answers and questions. J. Lout Temp.Phys., 93:1031-1057,1993.

$\mathrm{G}.\mathrm{G}$. Ihas, O. Avenel, R. Aarts,$\mathrm{R}$Salmelin, and E. Varoquaux. Phys. Rev.

Lett, 69(2):327, 1992.

[44] C. A. Jones and P. H. Roberts. J. Phys.

$A$, 15:2599, 1982.

[45] J. Wilks. The properties of liquid and solid helium. Clarendon Press, Oxford,

1967.

. Kuklewicz,. Durfee and[46] C. Raman, M. Kohl, R. Onofrio,in abosevelocity
Z. Hadzibabic, and W. Ketterle. Evidence for acriticaleinstein con densed gas. Phys. Rev. Lett, 83(13):2502, 1999.
$\mathrm{D}.\mathrm{S}$$\mathrm{C}.\mathrm{E}$[47] M. Ueda and A. J. Leggett. Macroscopic quantum tunneling of abose-einsteincondensate with attractive interaction. Phys. Rev. Lett, 80(8):1576, 1998.

[48] H. T. C. Stoof. Macroscopic quantum tunneling of abose condensate. J. Stat.
Phys., 87:1353, 1997.

[49] H. Shi and W.-M. Zheng. Bose-einstein condensation in an atomic gas withattractive interactions. Phys. Rev. , 55(4) $:2930$ ,1997.
$A$[50] R. J. Dodd, Mark Edwards, C. J. Williams, C. W. Clark, M. J. Holland, P. A.
Ruprecht, and K. Burnett. Role of attractive interactions on bose-einsteincondensation. Phys. Rev. , 54(1) $:661$ ,1996.
$A$[51] C. A. Sackett, H. T. C. Stoof, and R. G. Hulet. Growth and collapse of aboseeinstein condensate with attractive interactions. Phys. Rev. Lett, 80(10):2031,
1998.
[52] C. Huepe. Bifurcations et instabilites dans les condensats de Bose-Einstein etthesis, Ecole Normale Supe\’erieure, 1999.
les icoulements superfluides.
$\mathrm{P}\mathrm{h}\mathrm{D}$

42

218
[53] J. Guckenheimer and P. Holmes. Nonlinear Oscillations, Dynamical Systems
and Bifurcations of Vector Fields. Springer-Verlag, Berlin, 1983.
[54] C. C. Bradley, C. A. Sackett, J. J. Tollett, and R. G. Hulet. Evidence of boseeinstein condensation in an atomic gas with attractive interactions. Phys. Rev.
Lett, 75(9):1687, 1995.

[55] P. A. Ruprecht, M. J. Holland, K. Burnett, and Mark. Edwards. Timedependent solution of the nonlinear schrodinger equation for bose-condensed
trapped neutral atoms. Phys. Rev. , 51(6) $:4704$ ,1995.
$A$[56] C.W. Gardiner. Handbook of Stochastic Methods. Springer-Verlag, Berlin, 1985.

43