M E Brachet

M Abid

C Huepe

S Metens

C Nore

C T Pham

L S Tuckerman

Gross-Pitaevskii Dynamics of Bose-Einstein Condensates and Superfluid Turbulence

The Gross-Pitaevskii equation, also called the nonlinear Schrödinger equation (NLSE), describes the dynamics of low-temperature superflows and Bose-Einstein Condensates (BEC). We review some of our recent NLSE-based numerical studies of superfluid turbulence and BEC stability. The relations with experiments are discussed.

Introduction

The present paper is a review of results, obtained by our group during the last 10 years, by numerically studying the nonlinear Schrödinger equation (NLSE). Direct numerical simulations (DNS) and branch-following methods were extensively used to investigate the dynamics and stability of NLSE solutions in 2 and 3 space dimensions.

Much work has been devoted to the determination of the critical velocity at which superfluidity breaks into a turbulent regime [START_REF] Donnelly | Quantized Vortices in Helium II[END_REF]. A mathematical model of superfluid 4 He, valid at temperatures low enough for the normal fluid to be negligible, is the nonlinear Schrödinger equation (NLSE), also called the Gross-Pitaevskii equation [START_REF] Gross | Structure of a quantized vortex in boson systems[END_REF][START_REF] Pitaevskii | Vortex lines in an imperfect Bose gas[END_REF][START_REF] Landau | Fluid Mechanics[END_REF]. In a related context, dilute Bose-Einstein condensates (BEC) have been recently produced experimentally [START_REF] Anderson | Observation of Bose-Einstein condensation in a dilute atomic vapor[END_REF][START_REF] Davis | Bose-Einstein condensation in a gas of sodium atoms[END_REF][START_REF] Bradley | Bose-Einstein condensation of lithium: Observation of limited condensate number[END_REF]. The dynamics of these compressible nonlinear quantum fluids is accurately described by the NLSE allowing direct quantitative comparison between theory and experiment [START_REF] Dalfovo | Theory of Bose-Einstein condensation in trapped gases[END_REF].

Excitations of superfluid 4 He are described by the famous Landau spectrum which includes phonons in the low wave number range, and maxons and rotons in the high (atomic-scale) wave number range. In contrast, the standard NLSE (the equation used in the present paper) used only phonon excitations. It therefore incompletely represents the atomic-scale excitations in superfluid 4 He. However, note that there exist generalizations of the NLSE [START_REF] Pomeau | Model of superflow with rotons[END_REF][START_REF] Roberts | The nonlinear Schrödinger equation as a model of superfluidity[END_REF] that do reproduce the correct excitation spectrum, at the cost of introducing a spatially non-local interaction potential.

Several problems pertaining to superfluidity and BEC can thus be studied in the framework of the NLSE. In this review, we concentrate on two such problems: (i) low-temperature superfluid turbulence [START_REF] Nore | Decaying kolmogorov turbulence in a model of superflow[END_REF][START_REF] Nore | Kolmogorov turbulence in lowtemperature superflows[END_REF][START_REF] Abid | Experimental and numerical investigations of low-temperature superfluid turbulence[END_REF] and (ii) stability of BEC in the presence of a moving obstacle [START_REF] Huepe | Solutions de nucléation tourbillonnaires dans un modèle d'écoulement superfluide[END_REF][START_REF] Huepe | Scaling laws for vortical nucleation solutions in a model of superflow[END_REF][START_REF] Nore | Subcritical dissipation in threedimensional superflows[END_REF] or an attractive interaction [START_REF] Huepe | Decay rates in attractive Bose-Einstein condensates[END_REF][START_REF] Huepe | Stability and decay rates of non-isotropic attractive Bose-Einstein condensates[END_REF].

The authors recognize that this paper relates for the most part to their own work and does not include important and relevant contributions by other authors. However, we now give a short (and partial) list of references to provide a starting point to the reader motivated to undertake a deeper exploration of the field. A recent review of superfluid turbulence can be found in [START_REF] Vinen | Quantum turbulence[END_REF] and a proceeding devoted to the same subject is [START_REF] Barenghi | Quantized vortex dynamics and superfluid turbulence[END_REF]. Alternate simulations (by Biot-Savart vortex methods) of low-temperature superfluid turbulence can be found in [START_REF] Araki | Energy spectrum of superfluid turbulence with no normal-fluid component[END_REF]. Reconnection and acoustic emission are studied in [START_REF] Ogawa | Study of reconnection and acoustic emission of quantized vortices in superfluid by the numerical analysis of the Gross-Pitaevskii equation[END_REF]. A standard reference on the subject of vortex reconnection is [START_REF] Koplik | Vortex reconnection in superfluid Helium[END_REF]. Cascade processes and Kelvin waves are investigated in [START_REF] Araki | Cascade process of vortex tangle dynamics in superfluid 4 He without mutual friction[END_REF][START_REF] Leadbeater | Decay of superfluid turbulence via Kelvin-wave radiation[END_REF][START_REF] Kivotides | Kelvin waves cascade in superfluid turbulence[END_REF][START_REF] Svistunov | Superfluid turbulence in the low-temperature limit[END_REF] and [START_REF] Vinen | Decay of turbulence at a very low temperature: the radiation of sound from a Kelvin wave on a quantized vortex[END_REF]. Recent (very) low-temperature experiments in Helium are described in [START_REF] Davis | Decay of quantized vorticity in superfluid 4 He at mK temperature[END_REF]. The experimental field of Bose-Einstein condensation is in rapid evolution. Recent results include the observation of an isolated quantum vortex [START_REF] Matthews | Vortices in a Bose-Einstein condensate[END_REF][START_REF] Inouye | Observation of vortex phase singularities in Bose-Einstein condensates[END_REF] and the nucleation of several vortices [START_REF] Madison | Vortex formation in a stirred Bose-Einstein condensate[END_REF]. Details of vortex dynamics [START_REF] Rosenbush | Dynamics of a single vortex line in a Bose-Einstein condensate[END_REF] and even Kelvin waves [START_REF] Bretin | Quadrupole oscillation of a single-vortex condensate: evidence for Kelvin modes[END_REF] are now being observed.

The present paper is organized as follows: in section 2 the basic definitions and properties of the model of superflow are given. A short presentation of the hydrodynamic form, through Madelung's transformation, of NLSE with an arbitrary nonlinearity is derived. Simple solutions are discussed. Section 3 is devoted to superfluid turbulence. The basic tools that are needed to numerically study 3D turbulence using NLSE are developed and validated in Section 3.1. The NLSE numerical results are given in section 3.2. Experimental results are given in section 3.3.

The stability of BEC is studied in section 4. Exact 1D results are given in section 4.1 and a general formulation of stability is given in section 4.2. Numerical branch-following methods are explained in section 4.3. The stability of a superflow around a cylinder is studied in section 4.4. The stability of an attractive Bose-Einstein condensate is studied in section 4.5. Finally, section 5 is our conclusion.

Hydrodynamics using the NLSE

The hydrodynamical form of NLSE with an arbitrary nonlinearity, corresponding to a barotropic fluid with an arbitrary equation of state, is introduced in this section. Basic hydrodynamical features such as acoustic propagation and stationary vortex solutions are also discussed.

Madelung's transformation

The connection between the NLSE and fluid dynamics can be obtained directly using the following action [START_REF] Spiegel | Fluid dynamical form of the linear and nonlinear Schrödinger equations[END_REF]:

A = 2α dt d 3 x i 2 ψ ∂ψ ∂t -ψ ∂ψ ∂t -F (1)
with

F = d 3 x α|∇ψ| 2 + f (|ψ| 2 ) , (2)
where ψ( x, t) is a complex wave field and ψ its complex conjugate, α is a positive real constant and f is a polynomial in |ψ| 2 ≡ ψψ with real coefficients:

f (|ψ| 2 ) = -Ω|ψ| 2 + β 2 |ψ| 4 + f 3 |ψ| 6 + . . . + f n |ψ| 2n . ( 3
)
The NLSE is the Euler-Lagrange equation of motion for ψ corresponding to [START_REF] Donnelly | Quantized Vortices in Helium II[END_REF]. It reads

∂ψ ∂t = -i δF δψ , or ∂ψ ∂t = i(α∇ 2 ψ -ψf (|ψ| 2 )). (4)
Madelung's transformation [START_REF] Donnelly | Quantized Vortices in Helium II[END_REF][START_REF] Spiegel | Fluid dynamical form of the linear and nonlinear Schrödinger equations[END_REF]

ψ = √ ρ exp i ϕ 2α (5)
maps the nonlinear wave dynamics of ψ into equations of motion for a fluid of density ρ and velocity v = ∇ϕ. Indeed, using equation [START_REF] Anderson | Observation of Bose-Einstein condensation in a dilute atomic vapor[END_REF], equation ( 1) can be written as:

A = -dtd 3 x ρ ∂ϕ ∂t + 1 2 ρ(∇ϕ) 2 + 2αf (ρ) + 1 2 (2α∇( √ ρ)) 2 (6)
and the corresponding Euler-Lagrange equations of motion become:

∂ρ ∂t + ∇ • (ρ v) = 0, ( 7
)
∂ϕ ∂t + 1 2 (∇ϕ) 2 + 2αf (ρ) -2α 2 ∆ √ ρ √ ρ = 0. ( 8
)
These equations are the continuity and Bernoulli equation [START_REF] Landau | Fluid Mechanics[END_REF] for an isentropic, compressible and irrotational fluid if one drops the last term of ( 8). This term is called the "quantum pressure".

Using this identification, one can define the following "thermodynamical functions"1 . First, by inspecting the Bernoulli equation, the fluid's enthalpy per unit mass is given by:

h = 2αf (ρ). (9)
Second, noting that 1 2 ρ(∇ϕ) 2 corresponds to kinetic energy in equation ( 6), the fluid's internal energy per unit mass reads:

e = 2αf (ρ) ρ . ( 10
)
The general thermodynamic relation

h = e + p/ρ, (11)
gives the expression

p = 2α(ρf (ρ) -f (ρ)) (12)
for the fluid's pressure.

The units of the variables used in (2) and (3) can be recovered as follows: Madelung's transformation [START_REF] Anderson | Observation of Bose-Einstein condensation in a dilute atomic vapor[END_REF]

leads to [|ψ| 2 ] = [ρ] = M L -3 and [α] = L 2 T -1 .Using (10), one gets [f (ρ)/ρ] = T -1 and thus, from (3), [Ω] = T -1 , [β] = T -1 ρ -1 and [f i ] = T -1 ρ 1-i .
Note that, in the case of a Bose condensate of particles of mass m, α has the value /2m [START_REF] Nozières | The Theory of Quantum Liquids[END_REF].

Sound waves

Dispersion relation

The effect of the quantum pressure term in (8) can be found, at least to linear order, by means of the dispersion relation of acoustic (density) waves propagating in a constant density background ρ 0 . Writing ρ = ρ 0 + δρ (with f (ρ 0 ) = 0), ∇ϕ = δu in [START_REF] Bradley | Bose-Einstein condensation of lithium: Observation of limited condensate number[END_REF] and in the gradient of ( 8), one obtains (keeping only the linear terms):

∂ t δρ + ρ 0 ∇δu = 0 ∂ t δu + 2αf (ρ 0 )∇δρ -2α 2 ∆ ∇δρ 2ρ 0 = 0 or ∂ t 2 δρ = 2αρ 0 f (ρ 0 )∆δρ -α 2 ∆ 2 δρ.
The dispersion relation for an acoustic wave δρ = (exp(i(ωt

-k • x)) + c.c.) (with 1) is thus ω = 2αρ 0 f (ρ 0 ) k 2 + α 2 k 4 . ( 13
)
It is clear from this relation that the quantum pressure has a noticeable dispersive effect for large wave numbers. For small wave numbers, the usual propagation, with a constant speed of sound given by:

c = ∂p ∂ρ 1 2 = 2αρ 0 f (ρ 0 ),
is recovered. The length scale ξ = α/(ρ 0 f (ρ 0 )) at which dispersion becomes noticeable is known as the coherence length.

Nonlinear acoustics

The description given by linear acoustics can be somewhat improved by including the dominant nonlinear effects. Such an equation was derived in [START_REF] Nore | Numerical study of hydrodynamics using the nonlinear Schrödinger equation[END_REF].

Numerical simulations of NLSE in one space dimension using a standard Fourier pseudo-spectral method [START_REF] Gottlieb | Numerical Analysis of Spectral Methods[END_REF] can be used to study the acoustic regime triggered by an initial disturbance of the form :

ψ(x) = 1 + ae -x 2 l 2 .
Such simulations were performed in ref. [START_REF] Nore | Numerical study of hydrodynamics using the nonlinear Schrödinger equation[END_REF] where it was found that the shocks which would have appeared under compressible Euler dynamics (i.e. following equation ( 8) without the last term in l.h.s.) are regularized by the dispersion.

There was no evidence of finite-time singularity in our numerics: the spectrum of the solution was well resolved, with a conspicuous exponential tail.

Vortices in 2 and 3D

Stationary solutions of the equations of motion can give more insight into the connection between the NLSE and fluid dynamics. Indeed, stationary solutions of NLSE ( 4) are also solutions of the Real Ginzburg-Landau Equation (RGLE)

∂ψ ∂t = - δF δψ = (α∇ 2 ψ -ψf (|ψ| 2 )). (14)
They are thus extrema of the free energy F.

The simplest solution of this type corresponds to a constant density fluid at rest. Therefore, ψ is constant in space and (14) reads:

f (|ψ| 2 ) = -Ω + β|ψ| 2 + 3f 3 |ψ| 4 + . . . + nf n |ψ| 2n-2 = 0. ( 15
)
This equation, for given values of the coefficients β and f i , i = 3, . . . , n, fixes the Ω term of f by the fluid's density |ψ| 2 . However, Ω could be removed from the Bernoulli equation ( 8) by the change of variable ϕ → ϕ + 2αΩt that corresponds to a change of phase ψ → ψe iΩt in NLSE (4). Thus, Ω does not play an important role in the NLSE dynamics. It is however a matter of convention not to perform these changes of variable in order that stationary solutions of ( 14) coincide with those of (4).

Vortex solutions are another important kind of stationary solutions of NLSE. They are topological defects, or singularities, of Madelung's transformation when ρ = 0 (i.e. when both (ψ) = 0 and (ψ) = 0). These two conditions localize singularities into points in two dimensions and lines in three dimensions. The circulation of v around such a generic singularity is ±4πα. Therefore, they are known in the framework of superfluidity as quantum vortices [START_REF] Donnelly | Quantized Vortices in Helium II[END_REF]. Solutions of ( 14) with cylindrical symmetry are obtained numerically in [START_REF] Kawatra | Quantized vortices in imperfect Bose gas[END_REF]. It is found that the density profile of a vortex admits a horizontal tangent near the core while the velocity diverges as the inverse of the core distance. The momentum density ρ v is thus a regular quantity. It is important to realize that such vortex solutions are regular solutions of the NLSE (4), the singularity stemming only from Madelung's transformation [START_REF] Anderson | Observation of Bose-Einstein condensation in a dilute atomic vapor[END_REF].

Superfluid turbulence

The mathematical description of superfluid flows (i.e. laboratory 4 He flows) is based on Landau's two-fluid model [START_REF] Landau | Fluid Mechanics[END_REF]. The interaction of normal fluid and superfluid vortices is called mutual friction and must be taken into account as pioneered by Schwarz [START_REF] Schwarz | Three-dimensional vortex dynamics in superfluid 4 He: line-line and line-boundary interactions[END_REF]. At sufficiently low temperatures, one can neglect the normal fluid (below T = 1 • K for helium at normal pressure) and another mathematical description is given by the Nonlinear Schrődinger Equation NLSE (or Gross-Pitaevskii equation [START_REF] Gross | Structure of a quantized vortex in boson systems[END_REF][START_REF] Pitaevskii | Vortex lines in an imperfect Bose gas[END_REF]). Note that it is difficult to estimate the precise temperature below which the normal component can be neglected. There remains an urgent need for more experiments at much lower temperatures such as those reported in [START_REF] Davis | Decay of quantized vorticity in superfluid 4 He at mK temperature[END_REF].

In this section, we will use the simplest form for f , corresponding to a cubic nonlinearity in the NLSE (4). The NLSE, with convenient normalization, then reads:

∂ t ψ = (ic/ √ 2ξ)(ψ -|ψ| 2 ψ + ξ 2 ∇ 2 ψ). ( 16
)
Madelung's transformation [START_REF] Anderson | Observation of Bose-Einstein condensation in a dilute atomic vapor[END_REF] takes the form

ρ = |ψ| 2 (17)
ρv j = (icξ/ √ 2)(ψ∂ j ψ -ψ∂ j ψ) ( 18
)
where ξ and c are the coherence length defined above and speed of sound (when ρ 0 = 1 [START_REF] Donnelly | Quantized Vortices in Helium II[END_REF]) respectively. The superflow is irrotational everywhere but near the lines ψ = 0 (topological defects). There, the flow evolves under Eulerian dynamics [START_REF] Neu | Vortices in complex scalar fields[END_REF][START_REF] Lund | Defect dynamics for the nonlinear Schrödinger equation derived from a variational principle[END_REF]. The topological defect lines are the superfluid vortices, whose velocity circulation is automatically correct in this model [START_REF] Nozières | The Theory of Quantum Liquids[END_REF].

The present section is devoted to the analogy between turbulence in lowtemperature superfluids and classical turbulence in incompressible viscous fluids. This is done numerically by conducting numerical simulations of NLSE for the Taylor-Green (TG) vortex [START_REF] Taylor | Mechanism of the production of small eddies from large ones[END_REF] and comparing the results with prior Navier-Stokes simulations for the same vortex. The well-documented TG vortex [START_REF] Brachet | Small-scale structure of the Taylor-Green vortex[END_REF][START_REF] Brachet | Géométrie des structures à petite échelle dans le vortex de Taylor-Green[END_REF][START_REF] Domaradzki | An analysis of sugrid-scale interactions in numerically simulated isotropic turbulence[END_REF] is the solution of the Navier-Stokes equations with initial velocity field:

v T G = sin(x) cos(y) cos(z), -cos(x) sin(y) cos(z), 0 .

(

It admits symmetries that are used to speed up computations: rotation by π about the axis (x = z = π/2), (y = z = π/2) and (x = y = π/2) and reflection symmetry with respect to the planes x = 0, π, y = 0, π, z = 0, π. The velocity is parallel to these planes which form the sides of the impermeable box which confines the flow.

Tools for vortex dynamics

It is well known that compressible fluid dynamics, with an arbitrary chosen initial condition, leads to a flow dominated by acoustic radiation. We must thus generate an initial data with as small acoustic emission as possible if we are to use NLSE to study vortex dynamics.

Preparation method

We now present a method for generating a vortex array whose NLSE dynamics is similar to the classical vortex dynamics of the large scale flow v T G . Our method has two steps. In the first, we exhibit a global Clebsch representation of v T G . The acoustic wave emission is minimized in the second step [START_REF] Nore | Simulation numérique d'écoulements cisaillés tridimensionnels à l'aide de l'équation de Schrödinger non linéaire[END_REF].

The Clebsch potentials:

λ(x, y, z) = cos(x) 2 | cos(z)| (20) µ(x, y, z) = cos(y) 2 | cos(z)| sgn(cos(z)) (21)
(where sgn gives the sign of its argument) correspond to the TG flow in the sense that ∇ × v T G = ∇λ × ∇µ and λ and µ are periodic functions of (x, y, z). These Clebsch potentials map the physical space (x, y, z) into the (λ, µ) plane.

The complex field ψ c , corresponding to the large scale TG flow circulation, is given by ψ c (x, y, z) = (ψ 4 (λ, µ)) [γ d /4] with γ d = 2 √ 2/(πcξ) ([ ] denotes the integer part of a real) and

ψ 4 (λ, µ) = ψ e (λ -1/ √ 2, µ)ψ e (λ, µ -1/ √ 2) × ψ e (λ + 1/ √ 2, µ)ψ e (λ, µ + 1/ √ 2) ( 22
)
where The second step of our procedure consists of integrating, to convergence, the Advective Real Ginzburg-Landau Equation (ARGLE):

ψ e (λ, µ) = (λ + iµ)tanh( √ λ 2 + µ 2 / √ 2ξ)/ √ λ 2 + µ 2 .
∂ t ψ = c/( √ 2ξ)(ψ -|ψ| 2 ψ + ξ 2 ∇ 2 ψ) -iv T G • ∇ψ - (v T G ) 2 /(2 √ 2cξ)ψ (23)
with initial data ψ = ψ c .

Using the TG symmetries we expand ψ(x, y, z, t), solution of the ARGLE and NLSE equations, as [START_REF] Nore | Decaying kolmogorov turbulence in a model of superflow[END_REF]:

ψ(x, y, z, t) = N/2 m=0 N/2 n=0 N/2 p=0 ψ(m, n, p, t) cos mx cos ny cos pz ( 24
)
where N is the resolution and ψ(m, n, p, t) = 0, unless m, n, p are either all even or all odd integers. Furthermore ψ(m, n, p, t) satisfies the additional conditions ψ(m, n, p, t) = (-1) r+1 ψ(n, m, p, t) where r = 1 when m, n, p are all even and r = 2 when m, n, p are all odd. Implementing this expansion in a pseudospectral code yields a saving of a factor 64 in computational time and memory size when compared to general Fourier expansions.

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

Energy spectra

The total energy of the vortex array, conserved by NLSE dynamics, can be decomposed into three parts E tot = (1/2π) 3 d 3 x(E kin + E int + E q ), with kinetic energy E kin = ρv j v j /2, internal energy E int = (c 2 /2)(ρ -1) 2 and quantum energy E q = c 2 ξ 2 (∂ j √ ρ) 2 . Each of these parts can be defined as the integral of the square of a field, for example, E kin = ( √ ρv j ) 2 /2. In order to separate the kinetic energy corresponding to compressibility effects, E kin can be further decomposed into a compressible and incompressible parts using √ ρv j = ( √ ρv j ) c + ( √ ρv j ) i with ∇.( √ ρv j ) i = 0. Using Parseval's theorem, the angle-averaged kinetic energy spectrum is written as:

E kin (k) = 1 2 k 2 sin θdθdφ 1 (2π) 3 d 3 re ir j k j √ ρv j 2 which verifies E kin = (1/2π) 3 d 3 xE kin = ∞ 0 dkE kin (k).
The angle-averaged spectrum is obtained by summation over shells in Fourier space. A mode (m, n, p) belongs to the shell numbered

k = [ √ m 2 + n 2 + p 2 + 1/2].
Note that the radius of curvature of the vortex lines in Fig. 1 is large compared to their radii. Thus these 3D lines can be considered as straight, and then compared to the 2D axisymmetric vortices which are exact solutions to the 2D NLSE. A 2D vortex at the origin can be written as ψ vort (r) = ρ(r) exp(imϕ), m = ±1, where (r, ϕ) are polar coordinates. The vortex profile admits different limits ρ(r) ∼ r as r → 0 and ρ(r) = 1 + O(r -2 ) for r → ∞. It can be computed numerically using mapped Chebychev polynomials and an appropriate functional [START_REF] Nore | Decaying kolmogorov turbulence in a model of superflow[END_REF]. The corresponding velocity field is azimuthal and is given by v(r) = √ 2cξ/r. Using the mapped Chebychev polynomials expansion for ρ(r), the angle averaged spectrum of √ ρv j can then be computed with

the formula E vort kin (k) = (c 2 ξ 2 /2πk) ∞ 0 drJ 0 (kr)∂ r √ ρ 2 [11]
, where J 0 is the zeroth order Bessel function.

The incompressible kinetic energy spectrum E i kin of the ARGLE converged vortex array of Fig. 1 is displayed on Fig. 2. For large wavenumbers, the spectrum is well represented by extending a collection of 2D vortices into 3D vortex lines via E line kin (k) ≡ const. × E vort kin (k). (We will see that the constant of proportionality is related to the length l of vortex lines by const. = l/(2π) = 175 at time t = 0.) In contrast, the small wavenumber region cannot be represented by E line kin . This stems from the average separation distance between the vortex lines in Fig. 1. Denoting this distance d bump = k -1 bump = 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. Due to constructive interference, the energy spectrum at k = 2 has a value close to its corresponding value in TG viscous flow (namely 0.125), which greatly exceeds the value of E line kin (k = 2). In contrast, for 2 < k ≤ k bump , destructive interference decreases E i kin below E line kin .

Numerical results

The evolution in time via NLSE ( 16) of the incompressible kinetic energy is shown in Fig. 3. The main quantitative result is the excellent agreement of the energy dissipation rate, -dE i kin /dt, with the corresponding data in the incompressible viscous TG flow (see reference [START_REF] Brachet | Small-scale structure of the Taylor-Green vortex[END_REF], and reference [START_REF] Frisch | Turbulence, the legacy of A. N. Kolmogorov[END_REF], figure 5.12). Both the moment t max ∼ 5 -10 of maximum energy dissipation (the inflection point of Fig. 3) and its value (t max ) ∼ 10 -2 at that moment are in quantitative agreement. Furthermore, both t max and (t max ) depend weakly on ξ.

Another important quantity studied in viscous decaying turbulence is the scaling of the kinetic energy spectrum during time evolution and, especially, at the moment of maximum energy dissipation, where a k -5/3 range can be observed (see reference [START_REF] Brachet | Small-scale structure of the Taylor-Green vortex[END_REF]). Fig. 2 (b) shows the energy spectrum at t = 5.5. A least-square fit over the interval 2 ≤ k ≤ 16 with a power law E i kin (k) = Ak -n gives n = 1.70 (solid line). For 5 < t < 8, a similar fit gives n = 1.6 ± 0.2 (data not shown). Note that recent numerical simulations [START_REF] Araki | Energy spectrum of superfluid turbulence with no normal-fluid component[END_REF] using incompressible Eulerian dynamics of vortex filaments also show evidence of a k -5/3 energy spectrum range.

The full Kolmogorov law is given by E(k) = Cε 2/3 k -5/3 with E(k) the energy spectrum (per unit mass), C the Kolmogorov constant and ε the energy dissipation rate (per unit mass). The viscous energy spectrum E(k) is comparable (see ref. [START_REF] Brachet | Small-scale structure of the Taylor-Green vortex[END_REF]), in the inertial range, to the superfluid incompressible kinetic energy E i kin (k). As stated above, the superfluid dissipation rate -dE i kin /dt is also comparable to the viscous dissipation rate ε. Therefore, the Kolmogorov constant C is of the same order of magnitude in both viscous and superfluid flows. Note that the viscous Taylor-Green vortex, because of the non-homogenous character of the flow, has a Kolmogorov constant that is larger (by a factor 1.5) [START_REF] Brachet | Small-scale structure of the Taylor-Green vortex[END_REF] than that of homogeneous turbulence.

Fitting E i kin (k) in the interval 30 ≤ k ≤ 170 with l/(2π) × E vort kin (k) leads to l/2π = 452, roughly three times the t = 0 length of the vortex lines. The time evolution of l/2π obtained by this procedure is displayed in Fig. 3, showing that the length continues to increase beyond t max . The computations were performed with c = 2 corresponding to a root-mean-square Mach number

M rms ≡ |v T G rms |/c = 0.25.
As it is very costly to decrease M rms , we checked [START_REF] Nore | Decaying kolmogorov turbulence in a model of superflow[END_REF] that compressible effects were non-dominant at this value of M rms .

The vortex lines are visualized in physical space in Figs. 4 and5 at time t = 4 and t = 8. At t = 4, no reconnection has yet taken place while a complex vortex tangle is present at t = 8. Detailed visualizations (data not shown) demonstrate that reconnections occur for t > 5. Note that the viscous TG vortex also undergoes a qualitative (and quantitative) change in vortex dynamics around t ∼ 5.

Experimental results

The TG flow is related to an experimentally studied swirling flow [START_REF] Douady | Direct observation of the intermittency of intense vorticity filaments in turbulence[END_REF][START_REF] Zocchi | Measurement of the scaling of the dissipation at high Reynolds numbers[END_REF][START_REF] Fauve | Pressure fluctuations in swirling turbulent flows[END_REF]. The relation between the experimental flow and the TG vortex is a similarity in overall geometry [START_REF] Douady | Direct observation of the intermittency of intense vorticity filaments in turbulence[END_REF]: a shear layer between two counter-rotating eddies. The TG vortex, however, is periodic with free-slip boundaries while the experimental flow is contained inside a tank between two counter-rotating disks.

The spectral behavior of NLSE can be compared to standard (viscous) turbulence only for k ≤ k bump . It is thus of interest to estimate the scaling of k bump in terms of the characteristic parameters of the large scale flow and of the fluid. As seen above, k bump ∼ d -1 bump , where d bump is the average distance between neighboring vortices. Consider a flow with characteristic integral scale l 0 and large scale velocity u 0 (in the case of the TG flow, l 0 ∼ 1 and u 0 ∼ 1). The fluid characteristics are the velocity of sound c and the coherence length ξ (with corresponding wavenumber k ξ ∼ ξ -1 ). The number n d of vortex lines crossing a typical large-scale l 2 0 area is given by the ratio of the large-scale flow circulation l 0 u 0 to the quantum of circulation Γ = 4πcξ/ √ 2, i.e n d ∼ l 0 u 0 /cξ.

On the other hand, the assumption that the vortices are uniformly spread over the large scale area gives n d ∼ l 2 0 /d 2 bump . Equating these two evaluations of n d yields the relation d bump ∼ l 0 (cξ)/(l 0 u 0 ). Note that this argument assumes that the large-scale vorticity is coherent. It therefore yields a maximum possible value of d bump , and thus a minimum for k bump .

In the case of helium, the viscosity at the critical point ( T = 5.174 • K, P = 2.2 10 5 P a) is ν cp = 0.27 × 10 -7 m 2 s -1 while the quantum of circulation, Γ = h/m He has the value 0.99 × 10 -7 m 2 s -1 . Thus, ν cp ∼ 0.25Γ. The order of magnitude for d bump is thus d bump l 0 / R cp ∼ l λ where R cp is the integral scale Reynolds number at the critical point and l λ the Taylor micro-scale. In other words, the value of d bump in a superfluid helium experiment at T = 1 • K is of the same order as the Taylor micro-scale in the same experimental set-up run with viscous helium at the critical point.

The experimental set-up is similar to that described in [START_REF] Zocchi | Measurement of the scaling of the dissipation at high Reynolds numbers[END_REF]. To work with T ∼ 1.2K some modifications are, however, necessary. A cylinder, 8 cm in diameter and 12 cm high, limited axially by two counter-rotating disks limits the flow. One disk is flat and 8 radial blades, forming an angle of 45 o between each other, are fixed on the other one. To stabilize the turbulent shear region a stator is mounted halfway the height of the container. The two disks are driven by two DC motors rotating from 1 to 30 Hz. The whole system is enclosed in a liquid Helium bath used as the experimental fluid and this is the main difference with the set-up described in [START_REF] Zocchi | Measurement of the scaling of the dissipation at high Reynolds numbers[END_REF]. The pressure above the liquid bath is adjusted by a pumping system and this fixes the temperature of the fluid.

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 a quartz pressure transducer WHM 112 A22 from PCB. Details are given in [START_REF] Abid | Experimental and numerical investigations of low-temperature superfluid turbulence[END_REF].

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:

P meas (t) = P (t) + ρU 2 (t)/2 (25)
In the flow region where the probe is immersed, a well established axial mean flow U exists so that, after removing the mean parts of equation ( 25), one gets:

p meas (t) = p(t) + ρU u(t) (26)
where p meas , p and u are the fluctuations of the measured pressure, the actual pressure, and the local velocity respectively. It is currently admitted that, in ordinary turbulent situations, and at low fluctuation rates, equation ( 26) is dominated by the dynamic term, so that, by measuring the pressure fluctuations at the total head tube, one has a direct access to the velocity fluctuations u(t).

The situation is less clear when the probe is immersed in the superfluid. It is, however, possible to write an equation similar to [START_REF] Kivotides | Kelvin waves cascade in superfluid turbulence[END_REF]. Details can be found in [START_REF] Abid | Experimental and numerical investigations of low-temperature superfluid turbulence[END_REF].

The analysis of the pressure fluctuations obtained with the total head tube placed at 2 cm above the mid plane and 2 cm from the cylinder axis, yields interesting conclusions. Figures 6 and7 show the spectra of the pressure fluctuations above and below T λ (i.e respectively at 2.3K and 1.4K). Fig. 6 clearly shows, as expected, that such fluctuations follow a Kolmogorov regime between the injection scale (signaled by the peak at 25 Hz) and the largest resolved frequency, i.e 900 Hz. The spectrum obtained at 1.4K is similar to that obtained at T = 2.3K (see Fig. 7). A clear Kolmogorov like regime exists for the same range of frequencies. The corresponding Kolmogorov constant turns out also to be indistinguishable from the classical value. We have further analyzed the deviations from Kolmogorov in the superfluid regime. The striking result is that they have the same magnitude as in classical turbulence. More details are given in [START_REF] Maurer | Local investigation of superfluid turbulence[END_REF]. These observations -both on global and local quantities -agree well with the theoretical approach developed in the previous section. In particular, it seems clear that the Kolmogorov cascade survives in the superfluid regime.

Stability of stationary solutions

This section is devoted to the stability of BEC. Exact 1D bifurcation results are given in section 4.1. A general formulation of stability is presented in section 4.2. Numerical branch following methods are explained in section 4.3.

The stability of a superflow around a cylinder is investigated in section 4.4 and attractive Bose-Einstein condensates are studied in section 4.5.

Exact solution in 1D

Definition of the system

We consider a point impurity moving within a 1D superflow. In the frame of the moving impurity, the system can be described by the following action functional

A[ψ, ψ] = dt i 2 dx (ψ∂ t ψ -ψ∂ t ψ) -K . ( 27
)
In this expression, ψ is a complex field, ψ its conjugate and the energy functional K reads

K = E -vP + v R 2 (+∞)φ(+∞) -R 2 (-∞)φ(-∞) , (28)
with

E = dx |∂ x ψ| 2 + 1 2 (|ψ| 2 -1) 2 + gδ(x)(|ψ| 2 -1) , ( 29
)
P = dx 1 2i ψ(∂ x ψ) -ψ(∂ x ψ) , ( 30
)
ψ = R exp(iφ) . ( 31
)
The Dirac (pseudo) potential gδ(x) in ( 29) represents the impurity and the last term in [START_REF] Vinen | Decay of turbulence at a very low temperature: the radiation of sound from a Kelvin wave on a quantized vortex[END_REF] imposes the appropriate boundary conditions for the phase φ [START_REF] Hakim | Nonlinear Schrödinger flow past an obstacle in one dimension[END_REF]. R obeys the boundary conditions R 2 (±∞) = 1.

The Euler-Lagrange equation associated to [START_REF] Svistunov | Superfluid turbulence in the low-temperature limit[END_REF], δA/δψ = 0, is the nonlinear Schrödinger equation (NLSE)

i∂ t ψ = -∂ xx ψ + iv ∂ x ψ -ψ + |ψ| 2 ψ + gδ(x)ψ , ( 32
)
where the jump condition

∂ x ψ(0 + , t) -∂ x ψ(0 -, t) = gψ(0, t) (33)
is imposed in order to balance the gδ(x)ψ singularity with the -∂ xx ψ term for all times t.

Stationary solutions

Time-independent solutions of the NLSE [START_REF] Madison | Vortex formation in a stirred Bose-Einstein condensate[END_REF] are best studied by performing the change of variables defined above in [START_REF] Inouye | Observation of vortex phase singularities in Bose-Einstein condensates[END_REF]: ψ = R exp(iφ). Using these variables, the NLSE reads

∂ t R = v∂ x R -R∂ xx φ -2∂ x R∂ x φ , ( 34
)
∂ t φ = v∂ x φ -(∂ x φ) 2 + 1 -R 2 -gδ(x) + ∂ xx R R , ( 35
)
and the jump condition [START_REF] Rosenbush | Dynamics of a single vortex line in a Bose-Einstein condensate[END_REF] reads

∂ x R(0 + , t) -∂ x R(0 -, t) = gR(0, t) , ( 36
) ∂ x φ(0 + , t) -∂ x φ(0 -, t) = 0 . ( 37
)
Note that equations ( 34) and ( 35) can be respectively interpreted as the continuity and Bernoulli equations for a fluid of density ρ = R 2 (x) and velocity u = 2∂ x φ (as done in section 2).

Explicit time-independent solutions of equations ( 34) and [START_REF] Spiegel | Fluid dynamical form of the linear and nonlinear Schrödinger equations[END_REF] were found by Hakim [START_REF] Hakim | Nonlinear Schrödinger flow past an obstacle in one dimension[END_REF], using what are called gray solitons in the nonlinear optics terminology. Gray solitons [START_REF] Tsuzuki | Nonlinear waves in the Pitaevskii-Gross equation[END_REF][START_REF] Zakharov | Interaction between solitons in a stable medium[END_REF] are stationary solutions of equations ( 34) and [START_REF] Spiegel | Fluid dynamical form of the linear and nonlinear Schrödinger equations[END_REF], without the potential term gδ(x). They are localized density depletion of the form

R 2 GS (x) = v 2 /2 + (1 -v 2 /2) tanh 2 [ 1/2 -v 2 /4 x] , ( 38
)
φ GS (x) = arctan v √ 2 -v 2 exp[ √ 2 -v 2 x] + v 2 -1 . ( 39
)
Patching together pieces of gray solitons, Hakim found the following ξ-indexed stationary solutions of equations ( 34) and [START_REF] Spiegel | Fluid dynamical form of the linear and nonlinear Schrödinger equations[END_REF], including the potential term gδ(x)

R ξ (x) = R GS (x ± ξ) , x ≷ 0 (40) φ ξ (x) = φ GS (x ± ξ) -φ GS (±ξ) , x ≷ 0 ( 41
)
where the jump conditions ( 36) and (37) impose the relation

g(ξ) = √ 2(1 -v 2 /2) 3/2 tanh[ 1/2 -v 2 /4 ξ] v 2 /2 + sinh 2 [ 1/2 -v 2 /4 ξ] . ( 42
)
The function g(ξ) reaches a maximum [START_REF] Pham | Dynamical scaling laws in two types of extended hamiltonian systems at dissipation onset[END_REF] The two stationary solutions of (32) corresponding to ξ + (g) > ξ c and ξ -(g) < ξ c obtained by inverting (42) for g < g c thus disappear, merging in a saddlenode bifurcation at a critical strength g c . Note that the bifurcation can also be obtained by varying v and keeping g constant. In the following, the strength g of the delta function is used as the control parameter of our system, keeping v constant.

g c = g(ξ c ) at ξ c = argcosh 1+ √ 1+4v 2 2 √ 2-v 2 with g c = 4(1 -v 2 /2) [ √ 1 + 4v 2 -(1 + v 2 )] 1/2 2v 2 -1 + √ 1 + 4v 2 . ( 43
)
Fig. 8a shows the energetically unstable and stable solutions (K(ξ

-(g)) > K(ξ + (g)))
. The bifurcation diagram corresponding to the energy K (see equation [START_REF] Vinen | Decay of turbulence at a very low temperature: the radiation of sound from a Kelvin wave on a quantized vortex[END_REF]) is also displayed on the figure as an insert. In fig. 8b, note that the phase φ ξ (x), as defined in equation [START_REF] Neu | Vortices in complex scalar fields[END_REF], differs from that considered in [START_REF] Hakim | Nonlinear Schrödinger flow past an obstacle in one dimension[END_REF] by an (x-independent) constant. The phase in [START_REF] Hakim | Nonlinear Schrödinger flow past an obstacle in one dimension[END_REF] is set to 0 at x = +∞, whereas [START_REF] Neu | Vortices in complex scalar fields[END_REF] is antisymmetric in x. This difference is unimportant because equations ( 34) and ( 35) are invariant under the constant phase shift

φ(x) → φ(x) + ϕ . ( 44
)
20

General formulation

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

Consider the following action functional associated with the NLSE

A = d t dx i 2 ψ ∂ψ ∂ t -ψ ∂ψ ∂ t -F , ( 45
)
where ψ is a complex field, ψ its conjugate and F is the energy of the system.

Here, x and t correspond to nondimensionalized space and time variables respectively.

The Euler-Lagrange equation corresponding to [START_REF] Brachet | Géométrie des structures à petite échelle dans le vortex de Taylor-Green[END_REF] leads to the NLSE in terms of the functional

F ∂ψ ∂ t = -i δF δψ . ( 46
)
This equation obviously has stationary solution ψ S if δF/δψ| ψ=ψ S = 0. Thus, stationary solutions of ( 46) are extrema of F. In general, we are looking for an extremum of an energy functional E under some constraint Q[ψ] = cst. The usual Lagrange multiplier trick consists in introducing a control parameter ν and, rather than solving for extrema of E[ψ], searching for extrema of the new functional

F[ψ] = E[ψ] -νQ[ψ].
We thus solve for

δF δψ ν=cst. = 0. ( 47
)
We now turn to the precise definitions, corresponding to the two systems considered in this section : superflows and Bose-Einstein condensates.

Superflows

In the problem of a superflow past an obstacle, E is the hydrodynamic energy and ν ≡ U is the flow velocity with respect to the obstacle [START_REF] Huepe | Solutions de nucléation tourbillonnaires dans un modèle d'écoulement superfluide[END_REF][START_REF] Nore | Subcritical dissipation in threedimensional superflows[END_REF]. This implies that Q ≡ P is the flow momentum. Functionals F, E and P are given by the expressions

F = E -P • U (48) E = c 2 d 3 x [-1 + V ( x)]|ψ| 2 + 1 2 |ψ| 4 + ξ 2 |∇ψ| 2 (49)
P = √ 2cξ d 3 x i 2 ψ∇ψ -ψ∇ψ . ( 50
)
Here, c and ξ are the physical parameters characterizing the superfluid. They correspond to the speed of sound (c) for a fluid with mean density ρ 0 = 1, and to the coherence length (ξ). The potential

V ( x) = (V o /2)(tanh[4(r - D/2)/∆] -1
) is used to represent a cylindrical obstacle of diameter D. The computations are performed with V o = 10 and ∆ = ξ for which the density |ψ| ∼ 0 in the disk. The NLSE reads

∂ψ ∂t = - i √ 2cξ δF δψ = i c √ 2ξ [1 -V ( x)]ψ -|ψ| 2 ψ + ξ 2 ∇ 2 ψ + U • ∇ψ. ( 51
)
We will be interested in the solutions of δF/δψ = 0, for a given value of U . According to equation ( 47), these solutions are extrema of E at constant momentum P.

Bose-Einstein condensates

We consider a condensate of N particles of mass m and effective scattering length a in a radial confining harmonic potential V (r) = mω 2 r 2 /2 [START_REF] Huepe | Decay rates in attractive Bose-Einstein condensates[END_REF]. Quantities are rescaled by the natural quantum harmonic oscillator units of time τ 0 = 1/ω and length L 0 = /mω, thus obtaining the nondimensionalized variables t = t/τ 0 , x = x/L 0 and ã = 4πa/L 0 . The control parameter ν becomes in this context the chemical potential µ. The total number of particles in the condensate is therefore given by Q ≡ N . Functionals F, E and N are given, in terms of rescaled variables, by

F = E -µN (52) E = d 3 x 1 2 |∇ xψ| 2 + V (x)|ψ| 2 + ã 2 |ψ| 4 (53)
N = d 3 x|ψ| 2 . ( 54
)
Two different situations are possible, depending on the sign of the (rescaled) effective scattering length ã. When ã is positive the particles interact repulsively. A negative ã corresponds to an attractive interaction. The dynamical equation is

∂ψ ∂ t = -i δF δψ = i 1 2 ∇ 2 xψ - 1 2 |x| 2 ψ -ã|ψ| 2 -µ ψ . ( 55
)
We will be interested in the solutions of δF/δψ = 0, for a given value of µ. According to equation ( 47), these solutions are extrema of E at constant particle number N .

Branch following methods

When the extremum of F is a local minimum, the stationary solution ψ S of ( 51) can be reached by a relaxation method. If the extremum is not a minimum, Newton's iterative method is used to solve for ψ S .

Relaxation 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 the superflow past an obstacle:

∂ψ ∂t = -i δF δψ = i α∇ 2 ψ + [Ω -V ( x)] ψ -β|ψ| 2 ψ + U • ∇ψ. ( 56
)
When the extremum of F is a local minimum, the stationary solution ψ S of ( 56) can be reached by integrating to relaxation the associated real Ginzburg-Landau equation (RGLE)

∂ψ ∂t = - δF δψ = α∇ 2 ψ + [Ω -V ( x)] ψ -β|ψ| 2 ψ -i U • ∇ψ. ( 57
)
Indeed, ( 56) and ( 57) have the same stationary solutions.

In our numerical computations, equation ( 57) is integrated to convergence by using the Forward-Euler/Backwards-Euler time stepping scheme

ψ(t + σ) = Θ -1 1 -i σ U • ∇ + σ [Ω -V ( x)] -β|ψ(t)| 2 ψ(t) (58) with Θ = 1 -σ α∇ 2 . ( 59
)
The advantage of this method is that it converges to the stationary solution of (56) independently of the time step σ.

Newton's method

We use Newton's method [START_REF] Seydel | From Equilibrium to Chaos: Practical Bifurcation and Stability Analysis[END_REF] to find unstable stationary solutions of the RGLE.

In order to work with a well-conditioned system [START_REF] Mamun | Asymmetry and Hopf bifurcation in spherical Couette flow[END_REF], we search for the fixed points of [START_REF] Seydel | From Equilibrium to Chaos: Practical Bifurcation and Stability Analysis[END_REF]. These can be found as the roots of 60) for the problem of a superflow past a cylinder with ξ/D = 1/10 and a field ψ (j) discretized into n = 128 × 64 = 8190 collocation points. The error measure is given by n j=1 f 2 (j) (ψ)/n. The convergence is faster than exponential, as expected for a Newton method. where Θ -1 was already introduced in equation ( 58).

f (ψ) = Θ -1 1 -i σ U • ∇ + σ [Ω -V ( x)] -β|ψ(t)| 2 ψ(t) -ψ(t), (60)
We denote by ψ (j) the value of the field ψ at the j-th collocation point. The roots of f (ψ) are found by solving the linear problem

k df (j) dψ (k) δψ (k) = -f (j) (ψ). ( 61
)
for δψ and then incrementing ψ by ψ (j) = ψ (j) + δψ (j) [START_REF] Klein | Self-stretching of perturbed vortex filaments[END_REF] and iterating the Newton process ( 61)-( 62) to convergence.

The solution to ( 61) is obtained by an iterative bi-conjugate gradient method either BCGM [START_REF] Press | Numerical Recipees in Fortran[END_REF] or BiCGSTAB [START_REF] Van Der | Bi-cgstab: A fast and smoothly converging variant of bi-cg for the solution of nonsymmetric linear systems[END_REF]. These methods require only the ability to act repeatedly with [df (j) /dψ (k) ] on an arbitrary field ϕ to obtain an approximative solution of [START_REF] Van Der | Bi-cgstab: A fast and smoothly converging variant of bi-cg for the solution of nonsymmetric linear systems[END_REF]. Note that since the convergence of the time step [START_REF] Seydel | From Equilibrium to Chaos: Practical Bifurcation and Stability Analysis[END_REF] does not depend on σ, the roots found through this Newton iteration are also independent of σ. Therefore, σ becomes a free parameter that can be used to adjust the pre-conditioning of the system in order to optimize the convergence of the BCGM [START_REF] Mamun | Asymmetry and Hopf bifurcation in spherical Couette flow[END_REF].

Implementation

We use standard Fourier pseudo spectral methods [START_REF] Gottlieb | Numerical Analysis of Spectral Methods[END_REF]. Typical convergences of the Newton and bi-conjugate gradient iterations are shown in figures 9 and 10.

In the case of the radially symmetric Bose condensate, ψ(r, t) is expanded as ψ(r, t) = N R /2 n=0 ψ2n ( t)T 2n (r/R), where T n is the n-th order Chebychev polynomial and ψN R is fixed to satisfy the boundary condition ψ(R, t) = 0.

The NLSE is integrated in time by a fractional step (operator-splitting) method [START_REF] Klein | Self-stretching of perturbed vortex filaments[END_REF].

Stability of a superflow around a cylinder

In this section, following references [START_REF] Huepe | Solutions de nucléation tourbillonnaires dans un modèle d'écoulement superfluide[END_REF][START_REF] Huepe | Scaling laws for vortical nucleation solutions in a model of superflow[END_REF][START_REF] Nore | Subcritical dissipation in threedimensional superflows[END_REF], we investigate the stationary stable and unstable (nucleation) solutions of the NLSE describing the superflow around a cylinder, using the numerical methods developed in section 4.3. We study a disc of diameter D, moving at speed U in a two-dimensional (2D) superfluid at rest. The NLSE (51) can be mapped into two hydrodynamical equations by applying Madelung's transformation [START_REF] Donnelly | Quantized Vortices in Helium II[END_REF][START_REF] Spiegel | Fluid dynamical form of the linear and nonlinear Schrödinger equations[END_REF]:

ψ = √ ρ exp iφ √ 2cξ . ( 63
)
The real and imaginary parts of the NLSE produce for a fluid of density ρ and velocity

v = ∇φ -U , ( 64
)
the following equations of motion

∂ρ ∂t + ∇ (ρ v) = 0 (65) ∂φ ∂t -U • ∇φ + 1 2 (∇φ) 2 + c 2 [ρ -(1 -V ( x))] -c 2 ξ 2 ∇ 2 √ ρ √ ρ = 0. ( 66
)
In the coordinate system x that follows the obstacle, these equations correspond to the continuity equation and to the Bernoulli equation [START_REF] Landau | Fluid Mechanics[END_REF] (with a supplementary quantum pressure term c 2 ξ 2 ∇ 2 √ ρ/ √ ρ) for an isentropic, compressible and irrotational flow. Note that, in the limit where ξ/D → 0, the quantum pressure term vanishes and we recover the system of equations describing an Eulerian flow.

Bifurcation diagram and scaling in 2D

In this section, varying the ratio of the coherence length ξ to the cylinder diameter D, we obtain scaling laws in the ξ/D → 0 limit.

Bifurcation diagram

We present results for ξ/D = 1/10 which are representative of all ratios we computed. The functional E and energy F of the stationary solutions are shown in Fig. 11 We can relate branches in Fig. 11 to the presence of vortices in the solution.

When M n ≤ M ≤ M c , solutions are irrotational (M n ∼ 0.405 as indicated in Fig. 11). For M ≤ M n the stable branch (a) remains irrotational (Fig. 12A) while the unstable branch (b) corresponds to a one vortex solution (Fig. 12B) and the unstable branch (c), to a two vortex solution (Fig. 12C). The distance between the vortices and the obstacle in branches (b) and (c) increases when M is decreased. Branch (c) is precisely the situation described in [START_REF] Frisch | Transition to dissipation in a model of superflow[END_REF]. Furthermore, the value M c ≈ 0.4286 is close to the predicted value 2/11. Figure 12D shows the result of integrating the NLSE forward in time with, as initial condition, a slightly perturbed unstable symmetric stationary state (Fig. 12C). The perturbation drives the system over the nucleation barrier The point where vortices cross the surface of the disc (see Fig. 12) is labelled by M n . The total fluid momentum is given by -dF/dU (see text).

and cycles it, after the emission of two vortices, back to a stationary stable solution. This shows that the branch (c) corresponds to hyperbolic fixed points of NLSE.

Figures 12E,F show the phase of the field at the surface of the disc (r = D/2 and θ ∈ [0, 2π]) for four different flow speeds. In both unstable branches, 2πdiscontinuities, a diagnostic of vortex crossing, appear between M = 0.40 and M = 0.41.

Scaling laws

We now characterize the dependence on ξ/D of the main features of the bifurcation diagram. When ξ/D is decreased, M c and M pf become indistinguishable. In the limit where ξ/D = 0, the critical Mach number M c will be that of an Eulerian flow M c Euler .

Figure 13 shows the convergence of M c to the Eulerian critical velocity. This convergence can be characterized by fitting the polynomial law

M c = K 1 (ξ/D) K 2 + M c
Euler to M c (ξ/D). This fit is shown on Fig. 13 as a dotted line, yielding K 1 = 0.322, K 2 = 0.615 and M c Euler = 0.35.

Dynamical solutions

The stationary solutions obtained in the above subsection provide adequate initial data for the study of dynamical solutions. Indeed, after a small perturbation, their integration in time will generate a dynamical evolution with very small acoustic emission. Therefore, this procedure is an efficient way to initialize vortical dynamics in a controlled manner.

0 π/2 ¡ π ¢ -π/2 ¢ -π ¢ θ £ 0 π/2 π -π/2 -π phase ¤ 12 F 0 π/2 ¡ π ¢ -π/2 ¢ -π ¢ θ £ 0 π/2 π -π/2 -π phase
Starting from a two-vortex unstable stationary solution at a supercritical Mach number M c = 0.9, the evolution of the NLSE time integration shows a clearly periodic emission of vortex pairs (see Fig. 12). This emission conserves total circulation. sp law for the frequency is equivalent to that expected for a dissipative system.

We have studied the behavior of the frequency of vortex emission close to the bifurcation for these symmetric wakes with different supercritical velocities (characterized by δ sp = (M -M c )/M c = -δ > 0). Our results, plotted on Fig. 14, are consistent with a δ 1/2 sp scaling.

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 method

We used the 2D laminar stationary solution ψ 0V (x, y) (corresponding to branch (a) of the preceding section) and the one-vortex unstable stationary solution ψ 1V (x, y) (branch (b)) to construct the 3D initial condition

ψ 3D (x, y, z) = f I (z)ψ 1V (x, y) + [1 -f I (z)]ψ 0V (x, y). ( 67
)
The function f I (z), defined by

f I (z) = (tanh[(z -z 1 )/∆ z ] -tanh[(z -z 2 )/∆ z ])/2,
takes the value 1 for z 1 ≤ z ≤ z 2 and 0 elsewhere, with ∆ z an adaptation length.

Figure 15 represents a 3D initial data prepared with this method for ξ/D = 0.025,

| U |/c = 0.26 and ∆ z = 2 √ 2ξ in the [L x × L y × L z ] periodicity box (L x /D = 2.4 √ 2π, L y /D = 1.2 √ 2π and L z /D = 0.4 √ 2π
). The surface |ψ 3D | = 0.5 defines the cylindrical surface and the initial condition vortex line, with both ends pinned to the right side of the cylinder.

Short time dynamics

Starting from the initial condition (67), the evolution of the NLSE time integration shows different short-time and long-time dynamics.

During the short-time dynamics, the initial pinned vortex line rapidly contracts, evolving through a decreasing number of half-ring-like loops, down to a single quasi-stationary half-ring (see Figs. 16a,16b,16c). This evolution takes place mainly on the plane perpendicular to the flow, provided that the initial vortex is long enough to contract to a quasi-stationary half-ring as shown on Fig. 16c. Otherwise, the vortex line collapses against the cylinder while moving upstream.

Note that this quasi-stationary half-ring has been used by Varoquaux [START_REF] Avenel | The nucleation of vortices in superfluid 4 He: Answers and questions[END_REF][START_REF] Ihas | Quantum nucleation of vortices in the flow of superfluid He-4 through an orifice[END_REF] to estimate the nucleation barrier in a 3D experiment.

The dynamics of the half-ring situation (Fig. 16c) is very slow and can be shown to be close to a stationary field. Indeed, the local flow velocity v in an Eulerian flow around a cylindrical obstacle is known to vary from v = | U | at infinity to v = 2| U | at both sides of its surface. Moreover, the diameter d of a stationary vortex ring in an infinite Eulerian flow with no obstacle is given by [START_REF] Donnelly | Quantized Vortices in Helium II[END_REF]:

| U |/c = ( √ 2ξ/d) [ln (4d/ξ) -K] , (68)
where | U | is the flow velocity at infinity and the vortex core model constant K ∼ 1 is obtained by fitting the numerical results in [START_REF] Jones | Motions in a bose condensate: Iv. axisymmetric solitary waves[END_REF]. Therefore, for the values used on Figs. [START_REF] Nore | Subcritical dissipation in threedimensional superflows[END_REF]

Vortex stretching as a subcritical drag mechanism

A small 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 stretched while the pinning points move towards the back of the cylinder. log-log plot. The squares stand for our numerical stretching cases while the crosses correspond to non-stretching cases. There is a frontier between the 3D numerical dissipative and non-dissipative cases [START_REF] Nore | Subcritical dissipation in threedimensional superflows[END_REF]. For 1/30 < ξ/D < 1/20, the frontier corresponds to the expression R s = 5.5 with

R s ≡ | U |D/cξ = M D/ξ. ( 69
)
This superfluid 'Reynolds' number is defined in the same way as the standard (viscous) Reynolds number Re ≡ | U |D/ν (with ν the kinematic viscosity). It has been shown, in the superfluid turbulent (R s 1) regime, that R s is equivalent to the standard (viscous) Reynolds number Re [START_REF] Nore | Kolmogorov turbulence in lowtemperature superflows[END_REF][START_REF] Nore | Decaying kolmogorov turbulence in a model of superflow[END_REF][START_REF] Abid | Experimental and numerical investigations of low-temperature superfluid turbulence[END_REF]. Note that, for a Bose condensate of particles of mass m, the quantum of velocity circulation around a vortex, Γ = 2π √ 2cξ, has the Onsager-Feynman value Γ = h/m (h is Planck's constant) and the same physical dimensions L 2 T -1 as ν.

The value of R s divides the space of parameters into a laminar flow zone and [START_REF] Wilks | The properties of liquid and solid helium[END_REF]. Squares stand for our numerical stretching cases while crosses correspond to non-stretching cases [START_REF] Nore | Subcritical dissipation in threedimensional superflows[END_REF]. Note the wide difference in values of parameter ξ/D 10 -1 for our numerics (corresponding to current experiments in BEC) and ξ/D down to 10 -8 older experiments in liquid Helium). a recirculating flow zone, very much like in the problem of a circular disc in a viscous fluid in which this frontier is also found to be around Re ∼ 5. Thus, there seems to exist some degree of universality between viscous normal fluids and superfluids modeled by NLSE as discussed in [START_REF] Nore | Kolmogorov turbulence in lowtemperature superflows[END_REF][START_REF] Nore | Decaying kolmogorov turbulence in a model of superflow[END_REF][START_REF] Abid | Experimental and numerical investigations of low-temperature superfluid turbulence[END_REF]. In the context of superfluid 4 He flow, the experimental critical velocity is known to depend strongly on the system's characteristic size D. It is often found to be well below the Landau value (based on the velocity of roton excitation) except for experiments in which ions are dragged in liquid helium. Feynman's alternative critical velocity criterion R s ∼ log(D/ξ) is based on the energy needed to form vortex lines. It produces better estimates for various experimental settings, but does not describe the vortex nucleation mechanism [START_REF] Donnelly | Quantized Vortices in Helium II[END_REF].

In a recent experiment, Raman et al. have studied dissipation in a Bose-Einstein condensed gas by moving a blue detuned laser beam through the condensate at different velocities [START_REF] Raman | Evidence for a critical velocity in a Bose-Einstein condensed gas[END_REF]. In their inhomogeneous condensate, they observed a critical Mach number for the onset of dissipation M c 2D /1.6.

Our computations were performed for values of ξ/D comparable to those in Bose-Einstein condensed gas experiments. They demonstrate the possibility of a subcritical drag mechanism, based on 3D vortex stretching. It would be very interesting to determine experimentally the dependence of the critical Mach number on the parameter ξ/D and the nature (2D or 3D) of the excitations [START_REF] Nore | Subcritical dissipation in threedimensional superflows[END_REF].

Stability of attractive Bose-Einstein condensates

In this section, following reference [START_REF] Huepe | Decay rates in attractive Bose-Einstein condensates[END_REF], 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 a critical value N c [START_REF] Bradley | Bose-Einstein condensation of lithium: Observation of limited condensate number[END_REF]. Various physical processes compete to determine the lifetime of attractive condensates. Among them one can distinguish macroscopic quantum tunneling (MQT) [START_REF] Ueda | Macroscopic quantum tunneling of a Bose-Einstein condensate with attractive interaction[END_REF][START_REF] Stoof | Macroscopic quantum tunneling of a Bose condensate[END_REF], inelastic two and three body collisions (ICO) [START_REF] Shi | Bose-Einstein condensation in an atomic gas with attractive interactions[END_REF][START_REF] Dodd | Role of attractive interactions on Bose-Einstein condensation[END_REF] and thermally induced collapse (TIC) [START_REF] Stoof | Macroscopic quantum tunneling of a Bose condensate[END_REF][START_REF] Sackett | Growth and collapse of a Bose-Einstein condensate with attractive interactions[END_REF]. We compute the life-times, using both a variational Gaussian approximation and the exact numerical solution for the condensate wave-function.

Computations of stationary states

Gaussian approximation

A Gaussian approximation for the condensate density can be obtained analytically through the following procedure.

Inserting

ψ(r, t) = A( t) exp -r 2 /2r 2 G ( t) + ib( t)r 2 (70)
into the action [START_REF] Brachet | Géométrie des structures à petite échelle dans le vortex de Taylor-Green[END_REF], where F is given by ( 52), yields a set of Euler-Lagrange equations for r G ( t), b( t) and the (complex) amplitude A( t). The stationary solutions of the Euler-Lagrange equations produce the following values [START_REF] Huepe | Bifurcations et instabilités dans les condensats de Bose-Einstein et les écoulements superfluides[END_REF]:

N (µ) = 4 √ 2π 3 -8µ + 3 √ 7 + 4µ 2 7|ã| -2µ + √ 7 + 4µ 2 3/2 , ( 71
) E = N (µ) -µ + 3 7 + 4µ 2 /7. ( 72
)
N is found to be maximal at

N G c = 8 √ 2π 3 /|5 5/4 ã|. The corresponding value of the chemical potential is µ = µ G c = 1/2 √ 5.
Linearizing the Euler-Lagrange equations around the stationary solutions, yields the following expression for the eigenvalues [START_REF] Huepe | Bifurcations et instabilités dans les condensats de Bose-Einstein et les écoulements superfluides[END_REF]:

λ 2 (µ) = 8µ 2 -4µ 7 + 4µ 2 + 2 (73)
This qualitative behavior is the generic signature of a Hamiltonian saddle-node (HSN) bifurcation defined, at lowest order, by the normal form [START_REF] Guckenheimer | Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields[END_REF]

m ef f Q = δ -βQ 2 , ( 74
)
where δ = (1 -N /N c ) is the bifurcation parameter. The critical amplitude Q is related to the radius of the condensate [START_REF] Huepe | Bifurcations et instabilités dans les condensats de Bose-Einstein et les écoulements superfluides[END_REF]. We can relate the parameters β and m ef f to critical scaling laws, by defining the appropriate energy

E = E 0 + m ef f Q2 /2 -δ Q + βQ 3 /3 -γδ. (75)
From [START_REF] Huepe | Bifurcations et instabilités dans les condensats de Bose-Einstein et les écoulements superfluides[END_REF] it is straightforward to derive, close to the critical point δ = 0, the universal scaling laws

E ± = E c -E l δ ± E ∆ δ 3/2 , ( 76
) λ 2 ± = ±λ 2 ∆ δ 1/2 , ( 77
)
where

E c = E 0 , E l = γ, E ∆ = 2/3 √ β and λ 2 ∆ = 2 √ β/m ef 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 × 10 -3 , that corresponds to experiments with 7 Li atoms in a radial trap [START_REF] Bradley | Evidence of Bose-Einstein condensation in an atomic gas with attractive interactions[END_REF][START_REF] Bradley | Bose-Einstein condensation of lithium: Observation of limited condensate number[END_REF].

As is apparent on Fig. 21, the exact critical N E c = 1258.5 is smaller than the Gaussian value N G c = 1467.7 [START_REF] Ruprecht | Time-dependent solution of the nonlinear Schrödinger equation for Bose-condensed trapped neutral atoms[END_REF][START_REF] Ueda | Macroscopic quantum tunneling of a Bose-Einstein condensate with attractive interaction[END_REF]. The critical amplitudes corresponding to the Gaussian approximation can be computed from ( 71) and [START_REF] Dodd | Role of attractive interactions on Bose-Einstein condensation[END_REF]. One finds E c = 4 √ 2π 3 /|5 3/4 ã|, E ∆ = 64 √ π 3 /|5 9/4 ã| and λ 2 ∆ = 4 √ 10. For the exact solutions, we obtain the critical amplitudes by performing fits on the data.

One finds E ∆ = 1340 and λ 2 ∆ = 14.68. Thus, the Gaussian approximation captures the bifurcation qualitatively, but with quantitative 17% error on N c [START_REF] Ruprecht | Time-dependent solution of the nonlinear Schrödinger equation for Bose-condensed trapped neutral atoms[END_REF], 24% error on E ∆ and 14% error on λ 2 ∆ . Fig. 22 shows the physical origin of the quantitative errors in the Gaussian approximation. By inspection it is apparent that the exact solution is well approximated by a Gaussian only for small N on the stable (elliptic) branch.

Estimation of life-times

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 Γ T is estimated using the formula [START_REF] Gardiner | Handbook of Stochastic Methods[END_REF] Γ

T ω = |λ + | 2π exp -ω (E + -E -) k B T (78)
where ω(E + -E -) is the (dimensionalized) height of the nucleation energy barrier, T is the temperature of the condensate and k B is the Boltzmann constant. Note that the prefactor characterizes the typical decay time which is controlled by the slowest part of the nucleation dynamics: the top-of-the-barrier saddle point eigenvalue λ + . The behavior of Γ T can be obtained directly from the universal saddle-node scaling laws [START_REF] Bradley | Evidence of Bose-Einstein condensation in an atomic gas with attractive interactions[END_REF] and [START_REF] Ruprecht | Time-dependent solution of the nonlinear Schrödinger equation for Bose-condensed trapped neutral atoms[END_REF]. Thus the exponential factor and the prefactor vanish respectively as δ 3/2 and δ 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 quantum action path integral [START_REF] Stoof | Macroscopic quantum tunneling of a Bose condensate[END_REF][START_REF] Ueda | Macroscopic quantum tunneling of a Bose-Einstein condensate with attractive interaction[END_REF]. This trajectory is approximated as the solution of d 2 q(t) dt 2 = -dΦ(q) dq .

(79) Φ(q) is a polynomial such that -Φ(q) reconstructs the Hamiltonian dynamics, and is determined by the relations q m q f q m q b q f Fig. [START_REF] Koplik | Vortex reconnection in superfluid Helium[END_REF]. The bounce trajectory is shown as dashes, above the potential Φ(q).

Φ(q

m ) = -E + (80) Φ(q f ) = -E - (81) ∂ 2 q Φ(q m ) = |λ + (N )| (82) ∂ 2 q Φ(q f ) = -|λ -(N )|. ( 83
)
The bounce trajectory is displayed on Fig. 23 (dashed line) above the potential Φ(q). The MQT rate is estimated as

Γ Q ω = |λ -|v 2 0 4π exp -4 √ 2 q b q f Φ(q) -Φ(q f )dq , ( 84
)
where v 0 is defined by the asymptotic form of the bounce trajectory q(t) [START_REF] Stoof | Macroscopic quantum tunneling of a Bose condensate[END_REF]:

q(τ ) ∼ q f + (v 0 /|λ -|) exp[-|λ -τ |].
Universal scaling laws can be derived close to criticality from ( 74), [START_REF] Bradley | Evidence of Bose-Einstein condensation in an atomic gas with attractive interactions[END_REF] and [START_REF] Ruprecht | Time-dependent solution of the nonlinear Schrödinger equation for Bose-condensed trapped neutral atoms[END_REF]. The exponential factor in (84) follows the same scaling than |E + -E -|dq. It therefore vanishes as δ 5/4 . ¿From the asymptotic form of q(t), dq follows the same law as v 0 /|λ -|. Thus v 0 ∼ δ 3/4 and the prefactor vanishes as δ 7/8 .

Inelastic collision

The inelastic collision rate (ICO) is estimated using the relation

dN dt = f C (N ) (85)
with

f C (N ) = K |ψ| 4 d 3 x + L |ψ| 6 d 3 x, (86)
where K = 3.8 × 10 -4 s -1 and L = 2.6 × 10 -7 s -1 . The ICO rate can be evaluated from the stable branch alone. In order to compare the particle decay rate f C (N ) to the condensate collective decay rates obtained for TIC and MQT, we compute the condensate ICO half-life as:

τ 1/2 (N ) = N N /2 dn/f C (n) ( 87
)
and plot τ -1 1/2 on Fig. 24.

Discussion

It is apparent by inspection of Fig. 24 that for a given value of N the exact and Gaussian approximate rates are dramatically different. We now compare the relative importance of the different exact decay rates. At T ≤ 1 nK, the MQT effect becomes important compared to the ICO decay in a region very close to N E c (δ ≤ 8 × 10 -3 ) as was shown in [START_REF] Ueda | Macroscopic quantum tunneling of a Bose-Einstein condensate with attractive interaction[END_REF] using Gaussian computations but evaluating them with the exact maximal number of condensed particles N E c . Considering thermal fluctuations for temperatures as low as 2 nK, it is apparent on Fig. 24 (see insert) that the MQT will be the dominant decay mechanism only in a region extremely close to N c (δ < 5 × 10 -3 ) where the condensates will live less than 10 -1 s. Thus, in the experimental case of 7 Li atoms, the relevant effects are ICO and TIC, with cross-over determined in Fig. 24.

Calculations similar to those described here have also been generalized to an attractive Bose-Einstein condensate confined by a cylindrically symmetric potential V (r) = m(ω 2 r r 2 + ω 2 z z 2 )/2, which lead to cigar-shaped or pancakeshaped condensates. Newton's method was used to calculated stationary states, as described above. The inverse Arnoldi method was used to calculate the leading eigenvalues λ 2 ± , yielding decay rates that are similar to those computed for the spherically symmetric case. This was found to be due to spontaneous isotropization of the condensates. See [START_REF] Huepe | Stability and decay rates of non-isotropic attractive Bose-Einstein condensates[END_REF] for more details.

Conclusion

The main result of the NLSE simulations presented in section 3.2 is that two diagnostics of Kolmogorov's regime in decaying turbulence are satisfied. These diagnostics are, at the time of the maximum of energy dissipation: (i) a parameter-independent kinetic energy dissipation rate and (ii) a k -5/3 spectral scaling in the inertial range. Thus, the NLSE simulations were shown to be very similar, as far as energetics is concerned, with the viscous simulations. The experimental results presented in section 3.3 prove that the Kolmogorov cascade survives in the superfluid regime.

We have seen that the numerical tools developed in section 4.3 can be used in practice to obtain the stationary solutions of the NLSE. These methods have allowed us to find the full bifurcation diagrams of Bose-Einstein condensates with attractive interactions and superflows past a cylinder. Furthermore, the stationary solutions have given us efficient way to start vortical dynamics (in 2D and 3D) in a controlled manner.

Fig. 1 .

 1 Fig. 1. Three-dimensional visualization of the vector field ∇ × (ρ v) for the Taylor-Green flow at time t = 0 with coherence length ξ = 0.1/(8 √ 2), sound velocity c = 2 and N = 512 in the impermeable box [0, π] × [0, π] × [0, π].

Fig. 2 .

 2 Fig. 2. Plot of the incompressible kinetic energy spectrum, E i kin (k). The bottom curve (a) (circles) corresponds to time t = 0. The spectrum of a single axisymmetric 2D vortex multiplied by (l/2π) = 175 is shown as the bottom solid line. The top curve (b) (pluses) corresponds to time t = 5.5. A least-square fit over the interval 2 ≤ k ≤ 16 with a power law E i kin (k) = Ak -n gives n = 1.70 (top solid line).

Fig. 3 .

 3 Fig. 3. Total incompressible kinetic energy, E i kin , plotted versus time for ξ = 0.1/(2 √ 2), N = 128 (long-dash line); ξ = 0.1/(4 √ 2), N = 256 (dash); ξ = 0.1/(6.25 √ 2), N = 400 (dot) and ξ = 0.1/(8 √ 2), N = 512 (solid line). All runs are realized with c = 2. The evolution of the total vortex filament length divided by 2π (crosses) for the N = 512 run is also shown (scale given on the right y-axis).

Fig. 4 .

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

Fig. 5 .

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

Fig. 6 .

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

Fig. 7 .

 7 Fig. 7. Experimental pressure fluctuation spectrum (in non-dimensional units) measured with a total head pressure tube immersed in the flow at T = 1.4 K.

Fig. 8 .

 8 Fig. 8. (a) Modulus R of the stable (-) and unstable (---) stationary solutions of equation (32) (see equation (40)) for g = 1.250 and v = 0.5; insert, energy functional K of the stationary solutions versus g for v = 0.5 (see equation (28)); lower branch: energetically stable branch, upper branch: energetically unstable branch. The bifurcation occurs at g = 1.5514 (b) Phase φ of the stable (-) and unstable (---) stationary solutions (see equation (41)), same conditions as in (a).

Fig. 9 .

 9 Fig. 9. Two typical examples of the Newton method convergence towards the solution of equation (60) for the problem of a superflow past a cylinder with ξ/D = 1/10 and a field ψ (j) discretized into n = 128 × 64 = 8190 collocation points. The error measure is given by n j=1 f 2 (j) (ψ)/n. The convergence is faster than exponential, as expected for a Newton method.

Fig. 10 .

 10 Fig. 10. Two typical examples of a bi-conjugate gradient method convergence corresponding to the case shown on figure 9. The convergence of the relative error achieved for the x solution of Ax=b is given by |Ax-b|/|b|, where A= [df (j) /dψ (k) ], b= -f (j) (ψ) and x= δψ (k) .

  as a function of the Mach number (M = | U |/c). The stable branch (a) disappears with the unstable solution (c) at a saddle-node bifurcation when M = M c ≈ 0.4286. The energy F has a cusp at the bifurcation point, which is the generic behavior for a Hamiltonian saddle-node bifurcation, as described in section 4.5.1. There are no stationary solutions beyond this point. When M pf ≈ 0.4282, the unstable symmetric branch (c) bifurcates at a pitchfork to a pair of asymmetric branches (b). Their nucleation energy barrier is given by (F b -F a ) which is roughly half of the barrier for the symmetric branch (F c -F a ).

Fig. 11 .

 11 Fig. 11. Plot of the energy (F), and functional (E) versus Mach number (M = | U |/c), with D = 10ξ. Stable state (a). Nucleation solutions: asymmetric branch (b) and symmetric branch (c). The diagram shows a saddle-node and a pitchfork bifurcation.The point where vortices cross the surface of the disc (see Fig.12) is labelled by M n . The total fluid momentum is given by -dF/dU (see text).

¤Fig. 12 .

 12 Fig. 12. Stationary states : stable state (A), one vortex unstable state (B), two vortex unstable state (C). The surface indicates the fluid density around the cylinder (M = 0.24, ξ/D = 0.1). Figure (D) shows the result of the NLSE integration, starting from a slightly perturbed stationary (C) state. Figures (E) and (F ) display the phase of the complex field ψ at the surface of the cylinder versus the polar angle θ. Symmetric branch (E), asymmetric branch (F ). M = 0.4286 (•), M = 0.41 ( ), M = 0.40 (+), M = 0.30 (×). The crossing out of the vortex produces a phase discontinuity at M n ∼ 0.405.

Fig. 13 .Fig. 14 .

 1314 Fig. 13. Saddle-node bifurcation Mach number M c (+) and pitchfork bifurcation Mach number M pf (×), as a function of ξ/D. The dotted curve corresponds to a fit to the polynomial lawM c = K 1 (ξ/D) K 2 + M cEuler with K 1 = 0.322, K 2 = 0.615 and M c Euler = 0.35.

Fig. 15 .

 15 Fig. 15. Initial condition of a vortex pinned to the cylinder generated by equation (67). The surface |ψ 3D | = 0.5 is shown for ξ/D = 0.025, | U |/c = 0.26 and ∆ z = 2 √ 2ξ in the [L x × L y × L z ] periodicity box (L x /D = 2.4 √ 2π, L y /D = 1.2 √ 2π and L z /D = 0.4 √ 2π).

  , we expect that local velocities range from v = 0.25 to v = 2 × 0.25. Equation (68) thus implies that the diameter of an hypothetical stationary half-ring should be bounded by d(v = | U | = 0.25) = 18.8ξ and d(v = 2| U | = 0.5) = 6.3ξ. The diameter d ≈ 9ξ measured on the half-ring observed on Fig. 16c is consistent with its quasi-stationary behavior. Similarly the diameter of the half-ring shown on Fig. 18 d ≈ 7.6ξ is also found to be between the corresponding bounds d(0.35) = 11.4ξ and d(2 × 0.35) = 3ξ.

  Figures 17 show the long-time dynamics for a stretching case with ξ/D = 1/40 and | U |/c = 0.25 starting from Fig. 16c. Figure 19 shows a later situation for ξ/D = 1/20 and | U |/c = 0.35 starting from Fig.18. As the vortex loop grows, its backmost part remains oblique to the flow. This vortex stretching mechanism consumes energy, thus generating drag. It could be responsible for the appearance of drag in experimental superflows if fluctuations are strong enough to nucleate the initial vortex loop (which is imposed extrinsically in our numerical system). Note that it takes place for 2D subcritical velocities.

Figure 20 UFig. 16 .

 2016 Figure 20 displays several numerical and experimental[START_REF] Wilks | The properties of liquid and solid helium[END_REF] critical Mach numbers (V c /C) with respect to D/ξ, which seem to follow a (-1) slope in a

Fig. 17 .Fig. 18 .Fig. 19 .

 171819 Fig. 17. Long-time dynamics for ξ/D = 1/40 and | U |/c = 0.25 starting from Fig. 16c. The half-ring moves downstream while growing.

Fig. 20 .

 20 Fig.[START_REF] Barenghi | Quantized vortex dynamics and superfluid turbulence[END_REF]. Critical Mach number V c /C versus scale ratio of numerical and experimental data D/ξ. Circles correspond to several experiments from[START_REF] Wilks | The properties of liquid and solid helium[END_REF]. Squares stand for our numerical stretching cases while crosses correspond to non-stretching cases[START_REF] Nore | Subcritical dissipation in threedimensional superflows[END_REF]. Note the wide difference in values of parameter ξ/D 10 -1 for our numerics (corresponding to current experiments in BEC) and ξ/D down to 10 -8 older experiments in liquid Helium).

Fig. 21 .

 21 Fig. 21. Stationary solutions of the NLSE equation versus particle number N . Left: value of the energy functional E + on the stable (elliptic) branch and E -on the unstable (hyperbolic) branch. Right: square of the bifurcating eigenvalue (λ 2 ± ); |λ -| is the energy of small excitations around the stable branch. Solid lines: exact solution of the NLSE equation. Dashed lines: Gaussian approximation.

Fig. 22 .

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

Fig. 24 .

 24 Fig.[START_REF] Araki | Cascade process of vortex tangle dynamics in superfluid 4 He without mutual friction[END_REF]. 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[START_REF] Bradley | Bose-Einstein condensation of lithium: Observation of limited condensate number[END_REF]. The insert shows the details of the cross-over region between quantum tunneling and thermal decay rate. Solid lines: exact solution of the NLSE equation, dashed lines: Gaussian approximation.

Being isentropic (S = 0), the fluid is barotropic, and only one independent thermodynamic variable is needed.

Acknowledgments

The work reviewed in this paper was performed in collaboration with G. Dewel, J. Maurer and P. Tabeling. It was supported by ECOS-CONICYT program no. C01E08 and by the NSF DMR 0094569. Computations were performed at the Institut du Développement et des Ressources en Informatique Scientifique.