Dissipative dynamics of superfluid vortices at non-zero temperatures
Dissipative dynamics of superfluid vortices at non-zero temperatures
Natalia G. Berloff and Anthony J. Youd
Department of Applied Mathematics and Theoretical Physics,
University of Cambridge, Cambridge, CB3 0WA
We consider the evolution and dissipation of vortex rings in a condensate at non-zero temperature,
in the context of the classical field approximation, based on the defocusing nonlinear Schrödinger
equation. The temperature in such a system is fully determined by the total number density and
the number density of the condensate. A vortex ring is introduced into a condensate in a state
of thermal equilibrium, and interacts with non-condensed particles. These interactions lead to a
gradual decrease in the vortex line density, until the vortex ring completely disappears. We show
that the square of the vortex line length changes linearly with time, and obtain the corresponding
universal decay law. We relate this to mutual friction coefficients in the fundamental equation of
vortex motion in superfluids.
PACS numbers: 03.65.Sq, 03.75.Kk, 05.65.+b, 67.40.Vs, 67.57.De
The processes of self-organization, formation of large-
scale coherent localized structures and interactions of
these structures with small-scale fluctuations are at the
heart of nonlinear sciences, ranging from classical tur-
bulence, superfluids, ultracold gases and Bose–Einstein
Condensates (BECs), to the formation of the early Uni-
verse. Turbulence is characterised by the co-existence of
motions with many length and time scales described by
many degrees of freedom. The key to our understand-
ing of turbulence is to elucidate physics of interactions
between large scales (eg. large eddies) and small scales
(eg. turbulent fluctuations), and to develop mathemat-
ical models that account for the effects of small scales
without actually solving for them. According to the Lan-
dau description, superfluid 4He consists of the ground
state and the excitations — quasiparticles drifting on top
of the ground state. In the language of relativistic quan-
tum fields this corresponds to vacuum and matter (eg.
gravity waves interacting with a vacuum). Indeed, there
is a close relationship between superfluid hydrodynamics
and quantum gravity, so that at some level of hierarchy
of parameters the interactions of the quantum vacuum
and matter can be described by the defocusing nonlinear
Schrödinger (NLS) equation [1]. The dynamics of Bose
condensates depends on the energy exchange between the
condensed and non-condensed parts of the gas. Again,
the NLS equation (reformulated as the Gross–Pitaevskii
(GP) equation [2]) describes equilibrium and dynamical
properties of BEC as well as the formation of BEC from
a strongly degenerate gas of weakly interacting bosons
[3, 4]. The formation of the large-scale coherent localized
ground state (condensate) from a non-equilibrium initial
state has been studied in a number of papers addressing
different stages of the formation. Weak turbulence the-
ory has been used to predict the self-similar evolution of
the field in the regime of random phases of Fourier ampli-
tudes [5, 6], the transition from the regime of weak tur-
bulence to superfluid turbulence via states of strong tur-
bulence in the long-wavelength region of energy space [7],
and the final stage, resulting in the formation of a genuine
condensate [8, 9]. The related question about the effect
of finite temperature on the BEC dynamics has also been
addressed recently [10]. For instance, it was shown that
the presence of the thermal cloud in a trapped conden-
sate creates an effective dissipation that forces a single
vortex to move away from the centre and disappear [11].
The problem of the vortex tangle interacting with the
normal fluid (thermal cloud) is the key question in su-
perfluid turbulence. The Landau two-fluid theory of su-
perfluidity pre-dated the discovery of quantised vortex
lines and therefore omitted significant dynamical effects.
This was remedied— in the limit in which the mean spac-
ing between the vortex lines is small compared with any
other length scale of interest — by HVBK theory [12, 13].
In this limit, the superfluid vorticity is treated as a con-
tinuum, but the discrete nature of the vorticity gives rise
to an extra force on the superfluid component, arising
from the tension in the vortex lines. This term is absent
from the classical Euler equation of motion for an inviscid
fluid. The vortex lines also create a force of mutual fric-
tion between superfluid and normal fluid in addition to
the mutual friction included by Landau in his equations,
and represents the effects of collisions of the quasiparti-
cles with the vortex cores. Such forces were introduced
into the Landau model in an ad hoc way. This Letter
is the first attempt to study the effect of these collisions
quantitatively: we shall find the vortex line decay law
at non-zero temperature in the context of the defocus-
ing NLS equation. The NLS equation is a good starting
point, as the non-dissipative Landau two-fluid model can
be obtained from the equations of conservation of mass
and momentum for a one-component barotropic fluid us-
ing a general expression for the internal energy functional
of the density [14]. Through the Madelung transforma-
tion the NLS equation can be written in that form. Anal-
ogously, the transport coefficients in the Landau model
http://arxiv.org/abs/0704.1971v2
have been obtained directly from the NLS equation by
following the Chapman–Engskog expansion [15]. Note
that the separation of scales needed to carry out the
derivation of the Landau two-fluid model from the NLS
equation does not allow the inclusion of vortices as part
of the ground state. It is natural, therefore, to attempt
to derive the corresponding effects of the interactions of
vortices with the quasiparticles directly from the NLS
equation.
We consider the normalised defocusing NLS equation
for the complex function ψ [2]:
i∂tψ = −∇2ψ + |ψ|2ψ. (1)
The dynamics conserves the total number of parti-
cles N =
|ψ|2dx, and the total energy E =∫ (
|∇ψ|2 + 1
dx. We consider the uniform discrete
system of volume V = N 3, which is a periodic box on a
computational grid with 1283 discrete points.
Our goal is to determine the universal decay law for
the vortex line density in the entire range of tempera-
tures from 0 to the critical temperature of condensation,
Tλ. Our approach consists of three essential steps. We
aimed to: (1) achieve the thermal equilibrium state for
the given number of particles and given energy, starting
from a non-equilibrium stochastic initial condition for the
wavefunction ψ; (2) introduce a vortex ring into this state
and follow its decay via interactions with non-condensed
quasiparticles; (3) relate the decay rate to the tempera-
ture at equilibrium, where we derive the expression for
the relative temperature, T/Tλ, as a function of the total
number density, ρ = N/V , and the number density of
the condensate, ρ0.
We performed large scale numerical simulations of
Eq. (1) starting from a strongly non-equilibrium initial
condition[7], where the phases of the complex Fourier am-
plitudes ak(t) =
ψ(x, t)e−ip·x dx are distributed ran-
domly at t = 0. Here the momentum p takes quantised
values p = (2π/N )n with n = (0, 0, 0), (±1, 0, 0), · · ·.
This initial state describes a weakly interacting Bose
gas that is so rapidly cooled below the critical BEC tem-
perature that the particles remain in a strongly non-
equilibrium state. The kinetics of the initial weak tur-
bulent state has been analysed in [5, 6, 16] discovering a
quasiparticle cascade from high energies to low energies
in the wave number space. The ordering of the system
and the violation of the assumptions of weak turbulence
occurs very rapidly in a low-energy part of the spectrum,
with the formation of a quasi-condensate consisting of a
tangle of quantised vortices. The vortex tangle decays as
the system reaches a state of thermal equilibrium with
some portion (ρ0 ≡ |a0|2/V ) of particles occupying the
zero momentum state (genuine condensate) and the rest
of the non-condensed particles being distributed accord-
ing to the Rayleigh–Jeans equilibrium distribution [17],
modified by the presence of nonlinear interactions with
the condensate [9]:
p 6=0
|2 = T
ωB(p)
, (2)
where T is the temperature and ωB(p) is the Bogoli-
ubov dispersion relation (see below). An ultraviolet cut-
off for this distribution appears naturally through the
spatial discretization of the NLS equation. The numer-
ical scheme consists of fourth-order finite difference dis-
cretization in space and fourth-order Runge–Kutta in
time, so it is globally fourth-order accurate. This scheme
corresponds to the Hamiltonian system in the discrete
variables ψjkn, such that
iψ̇jkn =
, j, k, n = 1, ...,N , (3)
where
ψ∗jkn[
Ψ2 − 43Ψ1 +
ψjkn] +
|ψjkn|4
with Ψ2 = ψj+2,k,n + ψj−2,k,n + ψj,k+2,n + ψj,k−2,n +
ψj,k,n+2 + ψj,k,n−2 and Ψ1 = ψj+1,k,n + ψj−1,k,n +
ψj,k+1,n + ψj,k−1,n + ψj,k,n+1 + ψj,k,n−1.
The thermodynamic description of the condensation
process has been obtained in [9] by adapting the Bogoli-
ubov theory of a weakly interacting Bose gas [18] to the
classical system (1). We follow the same basic idea to de-
rive expressions for the energy and non-condensed part of
the discretised energy (4) written in terms of the Fourier
amplitudes ap as
K2(p)a
pap +
p1,p2,p3,p4
a∗p1a
ap3ap4δp1+p2−p3−p4 ,
where δp is the Kronecker delta symbol and
K2(p) =
sin2(pi/2)(7− cos(pi)). (6)
The Bogoliubov transformation bp = upap − vpa∗−p,
such that up = 1/
1−Q2p and vp = Qp/
1−Q2p with
Qp = [−K2−2ρ0+ωB(p)]/ρ0 diagonalises the term in (5),
which is quadratic in a0, to
ωB(p) b
pbp, where
excludes the p = 0 mode. Here ωB(p) =
K22 + 2ρ0K2
is the Bogoliubov-type dispersion relation.
Using the equilibrium distribution of the uncondensed
particles (2) the non-condensed number density can then
be expressed in terms of the basis used in this diagonal-
isation as
ρ− ρ0 =
′K2(p) + ρ0
. (7)
The discretised energy density H/V in the new basis
takes the form
ρ2 + (ρ− ρ0)2
1. (8)
The Eqs. (7)–(8) are analogous to Eqs. (8)–(9) of [9] but
modified for the discrete Hamiltonian discretization (4).
Given the energy density, H/V , and the total number
density, ρ, one can determine the temperature, T , at
equilibrium and the number density of the condensed
particles, ρ0, from Eqs. (7) and (8). The condensate
fraction ρ0/ρ as a function of the energy density H/V
is shown in FIG.1. This figure can be compared with
FIG.2 of [9] for the spectral representation of the total
energy. The analytical formulae (7)–(8) predict the sub-
critical behaviour of condensation, whereas the numerics
does not support this conclusion, as shown in the insert
of FIG.1. We use a linear approximation for small ρ0 to
determine the critical maximum energy for condensation
as shown in the insert. This energy is then used to deter-
mine the critical temperature for condensation Tλ (= T
for minH/V for which ρ0 = 0) from (7)–(8). We found
a phenomenological formula that determines T/Tλ as a
function of ρ0 and ρ as
1− α√ρ
− α√ρ
, (9)
where α is the only fitting parameter that we found as
α = 0.227538. The insert in FIG.1 shows the graph of
T/Tλ as a function of ρ0/ρ for ρ = 1/2. Eq. (9) gives an
excellent fit to the values computed from (7)–(8) across
all the values of ρ0 and ρ.
In order to analyze the decay of the vortex line length
at non-zero temperatures, we insert a vortex ring into a
state of thermal equilibrium and follow its decay due to
the interactions with the non-condensed particles. The
condensate healing length, which determines the size
of the vortex core, is calculated based on the density
of the condensate, and in our non-dimensional units is
ξ = 1/
ρ0. In healing lengths, the radius of the ring
is set to R0 = 10. The new initial state is ψv(t =
0) = ψeq ∗ ψvortex, where ψeq is the equilibrium state
and ψvortex is a wavefunction of the vortex ring [19]. The
vortex line length, L, is calculated as a function of time
with high frequencies being filtered out from the field ψ,
according to ãp = ap ∗ max(
1− p2/p2c, 0), where the
cut-off wavenumber is chosen as pc = 10(2π/N ) [20].
The first important conclusion of our numerical simula-
tions is that at all temperatures, the square of the vortex
line length decays linearly with time,
= −γ(ρ, T/Tλ), (10)
where γ does not depend on t. FIG. 2 shows this depen-
dence for various temperatures. The actual isosurfaces
of the decaying vortex line are shown in the inserts.
FIG. 1: (colour online) Condensate fraction, ρ0/ρ, as a func-
tion of the energy density as obtained from the numerical
simulations (points) and from the analytical expressions (7)–
(8) (solid line). The inserts show (a) the plot of T/Tλ as
a function of ρ0/ρ, obtained using Eqs. (7)–(8) or Eq. (9)
and (b) subcritical condensation predicted by Eqs. (7)–(8)
(black line), the linear approximation used to obtain the crit-
ical temperature of condensation (gray (red) line), and nu-
merical calculations (blue dots). The total number density is
ρ = 1/2.
0.5 1 1.5 2 H/V
2.55 2.56 2.57 2.58 2.59 H/V
0 0.2 0.4 0.6 0.8
FIG. 2: (colour online) The decay of the square of the vortex
line length as a function of time for various T indicated next
to the graphs. The fit to the linear function is shown by
the gray (red) lines. The inserts show isosurface plots of the
vortex line (for filtered fields ψ; see text) for T = 0.52Tλ at
time=130 (left) and time=1300 (right); between these two
times the vortex line length is reduced by a factor of 2. The
perturbations to the vortex line due to collisions with non-
condensed particles are clearly seen on the left insert. These
collisions generate Kelvin waves that also radiate energy to
sound. The total number density is ρ = 1/2.
100 200 300 400 500 600
0.44Tλ
0.27Tλ
0.63Tλ
This result agrees with predictions of the HVBK theory
for superfluid helium [12] according to which the funda-
mental equation of the motion of a vortex line, vL, is
given by (see also page 90, Eq. (3.17) of [21])
vL = vsl+αs
′×(vn − vsl)−α′s′×[s′×(vn − vsl)], (11)
where vsl is the local superfluid velocity that consists of
the ambient superfluid flow velocity and the self-induced
vortex velocity ui, vn is the normal fluid velocity, s is
a position vector of a point on the vortex and s′ is the
unit tangent at that point. Mutual friction parameters α
and α′ are ad hoc coefficients in the HVBK theory that
are functions of ρn, ρ, and T only. Eq. (11) is a gen-
eral and universal equation used to follow the evolution
of three-dimensional vortex motion in an arbitrary flow.
When formulated for a single vortex ring Eq. 11 reads
dR/dt = −αui, where ui = κ[log(8R/ξ) − δ + 1]/(4πR)
and δ is the vortex core parameter. For the GP vortices
δ ≈ 0.38 [2]. In dimensionless units used in our paper
ui = [log(8R) − δ + 1]/R. After integration of the equa-
tion for Ṙ we get αt = (R20 − R2)/[2(log(8R̂) + δ − 1)],
where R̂ is the mean radius of the ring. When this is com-
pared with (10) we get the following relationship between
γ and α: γ = 8π2(log(8R̂) + δ − 1)α. From our numer-
ics we obtained a general result valid across all ranges of
temperatures and total densities: γ ≈ Kρ(T/Tλ)2, where
K ≈ 68. Note that for a GP condensate T/Tλ ≈ ρn/ρ
to the first order (see insert (a) of FIG.1), so alterna-
tively, we can write γ ≈ K1ρn(T/Tλ) FIG. 3 shows the
comparison of the numerically calculated γ/ρ and the
quadratic fit K(T/Tλ)
2. Thus, we found that the mutual
friction coefficient in condensate superfluids is given by
α ≈ K2ρn(T/Tλ).
The existence of the transverse force on superfluid vor-
tices which is parametrised by the parameter α′ has been
a subject of much debate in mid-1990s, when calcula-
tions of the classical Magnus force applied to superfluid
vortices have been offered and argued about [22]. The
criticism is based on the observation that the classical
hydrodynamic equations are inapplicable in the vortex
core. Whether or not the details of the non-classical vor-
tex dynamics are crucial to the existence of the trans-
verse force is still an open question. The estimate of α′
can be obtained from our numerical procedure as follow-
ing. Eq. (11) written for a distance travelled by a single
vortex ring takes form (see Eq. (3.53) on page 107 of [21])
dz/dt = (1−α′)ui. We compared the distances travelled
by a vortex ring at various temperatures obtained nu-
merically with the distances travelled by a vortex ring in
the absence of the transverse force according to the ana-
lytical formula dz/dt = ui, where ui = ui(R(t)) and R(t)
varies with time according to (10). The insert of FIG.3
shows these distances for T/Tλ = 0.27. Our calculations
fail to detect any significant presence of the transverse
force for any temperature considered: the deviation from
the analytical curve is insignificant within the accuracy
of (10). We plan to perform a more thorough analyti-
cal and numerical study of transverse force from a single
FIG. 3: (Color online) Values of γ/ρ as a function of tempera-
ture T/Tλ for various values of the total number density ρ de-
picted in various shades of gray (in various colours): ρ = 1/2
(dark (red)), ρ = 1/4 (light (green)) and ρ = 3/4 (medium
(blue)). The plot of the quadratic fit γ/ρ = 68(T/Tλ)
given by the dashed line. The relative temperature is calcu-
lated using Eq. (9). The result is not sensitive to whether
we use the values of ρ and ρ0 that correspond to the state of
thermodynamical equilibrium before the introduction of the
vortex ring or after the vortex ring disappears and the sys-
tem equilibrates. The insert shows the distance travelled by
a vortex ring as a function of time for T/Tλ = 0.27 (red dots
– distances calculated using dz/dt = ui, black line using nu-
merics). Curves depart when the vortex ring becomes small
in radius and the analytical formula is no longer accurate ap-
proximation of the vortex velocity.
0 0.2 0.4 0.6 0.8
0 200 400 600 800
phonon acting on a single vortex in context of the GP
model in future.
In summary, we considered the effect of temperature
on the decay of vortex line density via interactions with
non-condensed particles in the context of the defocus-
ing NLS equation. We obtained a simple expression for
the temperature at equilibrium as a function of the to-
tal number density and the number density of the con-
densed particles. Depending on these two parameters,
a vortex ring introduced into the condensate shows dif-
ferent decay rates with time. We identified this decay
law as linear for the square of the vortex line length and
showed the universal dependence of the decay rate on
temperature and total density. It has been suggested
that the emission of sound by vortex reconnections and
vortex motion is the only active dissipation mechanism
responsible for the decay of superfluid turbulence. The
decay of superfluid turbulence via Kelvin wave radiation
and vortex reconnections was studied in the framework
of the GP equation [23] at near zero temperature, via
collision of two vortex rings, and confirmed that in the
Kelvin wave cascade, where energy is transferred to much
shorter wavelengths with a cut-off below a critical wave-
length, the vortex line density can be described by the
famous Vinen equation [24] d(L/V )/dt = −χ(L/V )2. It
has also been shown [25] that the presence of localized
finite amplitude sound waves greatly enhances the dissi-
pation of the vortex tangle, essentially changing the de-
cay law to exponential decay. This Letter complements
the existing Kelvin wave cascade scenario by consider-
ing an opposite limit when there are no reconnections,
and the decay mechanism depends only on the energy
exchange with non-condensed particles. This mechanism
exceeds the energy transfer via the Kelvin wave cascade.
Finally, we related our results about a single vortex ring
to the mutual friction coefficients in the general equation
of vortex motion.
NGB acknowledges the support from EPSRC-UK. She
is also very grateful to Professor Joe Vinen for several
illuminating discussions about superfluid turbulence and
his suggestion to look into the decay of vortex rings in
the context of the NLS equation.
[1] G. Volovik, gr-qc/0612134.
[2] V. L. Ginzburg and L. P. Pitaevskii, Sov. Phys. JETP
7, 858 (1958); L. P. Pitaevskii, Sov. Phys. JETP 13, 451
(1961); E. P. Gross J. Math. Phys. 4, 195 (1963).
[3] E. Levich and V. Yakhot, J. Phys. A: Math. Gen. 11,
2237 (1978).
[4] Yu. Kagan and B.V. Svistunov, Phys. Rev. Lett. 79, 3331
(1997).
[5] V.E.Zakharov, S.L. Musher, and A.M.Rubenchik, Phys.
Rep. 129, 285 (1985) and S. Dyachenko, A.C. Newell, A.
Pushkarev and V.E.Zakharov, Physica D 57, 96 (1992).
[6] B.V. Svistunov, J. Moscow Phys. Soc. 1, 373 (1991); Yu.
Kagan and B.V. Svistunov, Zh. Eksp. Theor. Fiz. 105,
353 (1994) [Sov. Phys. JETP 78, 187 (1994)].
[7] N. G. Berloff and B. V. Svistunov, Phys. Rev. A 66,
013603 (2002)
[8] M.J.Davis, S.A. Morgan, and K. Burnett, Phys. Rev. Lett
87, 160402 (2001) and Phys. Rev. A 66, 053618 (2002)
[9] C. Connaughton et al, Phys. Rev. Lett. 95, 26901 (2005).
[10] M.Brewczyk et al J. Phys. B: At. Mol. Opt. Phys. 40,
R1-R37 (2007) and reference within.
[11] H. Schmidt et al J.Opt.B: Quantum Simiclass. Opt. 5
S96 (2003).
[12] H.E. Hall and W. F. Vinen, Proc. R. Soc. Lond., A238,
215 (1956); I.L. Bekharevich and I.M. Khalatnikov Soviet
Phys., JETP, 13, 643 (1961).
[13] R.N. Hills, and P.H. Roberts, Archiv. Rat. Mech. &
Anal., 66, 43 (1977a); Int. J. Eng. Sci., 15, 305 (1977b);
J. Low Temp. Phys., 30, 709 (1978a); J. Phys. C11, 4485
(1978b).
[14] S.J. Putterman and P.H.Roberts, Physica, 117A, 369
(1983).
[15] T.R.Kirkpatrick and J. R. Dorfman J. Low Temp. Phys.
58, 301 (1985); 399 (1985).
[16] Yu. Kagan, B.V. Svistunov, and G.V. Shlyapnikov, Zh.
Eksp. Teor. Fiz. 101, 528 (1992) [Sov. Phys. JETP 75,
387 (1992)]; Yu. Kagan and B.V. Svistunov, Zh. Eksp.
Theor. Fiz. 105, 353 (1994) [Sov. Phys. JETP 78, 187
(1994)].
[17] V. E. Zakharov, V. S. L’vov and G. Falkovich, Kol-
mogorov Spectra of Turbulence I (Springer, Berlin, 1992);
A. C. Newell, S. Nazarenko and L. Biven, Physica D 152,
520 (2001).
[18] N. N. Bogoliubov, Journal of Physics 11, 23 (1947).
[19] N.G. Berloff J. Phys. A: Math. and Gen., 37(5), 1617
(2004).
[20] Various choices of pc taken from the interval (5, 15) ×
(2π/N ) give different L, but similar vortex line decay
rates. A more accurate, but more computationally inten-
sive, way to calculate the vortex line length in a quasi-
condensate (a condensate with vortices) is through a time
averaging of the field that removes the high frequencies,
see the discussion in [7].
[21] R.J.Donnelly “Quantized Vortices in Helium II”, Cam-
bridge University Press, Cambridge 1991.
[22] S.V.Iordanskii, Sov. Phys. JETP 22, 160 (1966); E. B.
Sonin, Sov. Phys. JETP 42 469 (1975); Phys. Rev. B 55,
485 (1997), Ao and Thouless, Phys. Rev. Lett. 70, 2158
(1993).
[23] M. Leadbeater, D.C. Samuels, C.F. Barenghi and C.S.
Adams, Phys. Rev. A 67 015601 (2003).
[24] W.F. Vinen, Proc. R. Soc. London. Ser. A 242, 493
(1957)
[25] N.G. Berloff Phys. Rev A, 69 053601 (2004)
http://arxiv.org/abs/gr-qc/0612134
