Towards quantum turbulence in finite temperature Bose-Einstein condensates

Motivated by the various indications that holographic superfluid is BCS like at the standard quantization but BEC like at the alternative quantization, we have implemented the alternative quantization in the dynamical holographic superfluid for the first time. With this accomplishment, we further initiate the detailed investigation of quantum turbulence in finite temperature BEC by a long time stable numerical simulation of bulk dynamics, which includes the two body decay of vortex number caused by vortex pair annihilation, the onset of superfluid turbulence signaled by Kolmogorov scaling law, and a direct energy cascade demonstrated by injecting energy to the turbulent superfluid. All of these results share the same patterns as the holographic superfluid at the standard quantization, thus suggest that these should be universal features for quantum turbulence at temperatures order of the critical temperature.


Introduction and motivation
Turbulence is defined loosely as sort of temporally and spatially complex state of fluid motion. There are two types of turbulence. One is classical turbulence (CT) in normal fluid, where the eddies can have arbitrary circulation. The other is quantum turbulence (QT) in superfluid, where the involved vortices are instead quantized by quantum mechanics. Due to the presence of quantized vortices in QT, the hydrodynamical description of turbulence breaks down in superfluid. Thus some new approaches are developed in condensed matter physics to address QT, such as Gross-Pitaevskii equation (GPE) for zero temperature superfluid [1][2][3], and Zaremba-Nikuni-Griffin (ZNG) formalism to take into account the finite temperature effect [4]. However, both of them are dealing with the weakly interacting systems. Here comes holography as a complementary tool to attack QT at strongly coupled regime. Holographic duality provides us with a powerful framework, in which a complete description of a strongly coupled quantum many-body system, valid at all scales, can be encoded in a classical gravitational system with one extra dimension [5][6][7]. In particular, the finite temperature effect can be beautifully implemented by putting a black hole in the dual bulk.
In order to address QT, one is required to first construct a holographic model for superfluid phase. Actually this has been successfully achieved by putting a U(1) gauge field minimally coupled to a complex scalar field with m 2 L 2 = −2 in AdS Schwarzschild, where L is AdS curvature radius [8,9]. As a result, when the temperature decreases to a certain critical value, the black hole will become hairy with the complex scalar field condensed, which corresponds exactly to the phase transition to a superfluid from a normal fluid on the boundary. Based on this simplest holographic superfluid model, the two dimensional QT at temperatures of order the critical temperature is investigated by holography in JHEP07(2016)092 the seminal work [10]. It is found that although QT shares the same Kolmogorov energy spectrum as CT, QT exhibits a direct energy cascade from the large scales to the small scales, which is distinct from the inverse energy cascade in CT. Further work shows that the finite temperature leads to the decrease of vortex number, the behavior of which is well captured by vortex pair annihilation to sound waves [11,12].
It is noteworthy that the above work on holographic superfluid turbulence focuses exclusively on the standard quantization [13]. But as investigated from various aspects in [14,15], this standard quantization might give rise to a BCS-like superfluid while the alternative quantization might correspond to a BEC-like superfluid. For example, the density depletion in the core of both a single dark solition and vortex seems to be saturating at 40% for the standard quantization and at 100% for the alternative quantization when the temperature is lowered. This is reminiscent of a loosely bound BCS-like fermion pair condensed system and a more tightly bound BEC-like fermion pair condensed system, respectively. With this in mind, one is tempted to ask whether the aforementioned patterns for QT are universal, or dependent on whether the superfluid is BCS-like or BEC-like. As such, it is significant for one to explore what happens to QT for the alternative quantization. At first sight, the alternative quantization appears to be technically a routine generalization of the standard one by simply changing Dirichlet boundary condition to Neumann one at the AdS boundary. This is actually the case for the situation without the temporal evolution, where the Schwarzschild coordinates are preferred. But it turns out to be a non-trivial numerical challenge to implement the alternative quantization for the situation involving the temporal evolution like QT, where the infalling Eddington coordinates are favorably adopted. The main purpose of this paper is to conquer this numerical challenge and explore holographic superfluid turbulence at the alternative quantization, where the focus is put on the possible universal features such as the vortex number decay behavior, the scaling law, and energy cascade direction. Various aspects, especially the time evolution of vortex number and the scaling behavior in the kinetic energy spectrum, are investigated more carefully than in previous works.
The paper is organized as follows. In the next section, we shall develop the holographic superfluid model at the alternative quantization in the infalling Eddington coordinates and elaborate on the relevant numerics for our numerical simulation of the bulk dynamics. Then we present our numerical results on the decrease of vortex number by vortex pair annihilation, the scaling behaviors in the kinetic energy spectrum, and energy cascade direction separately in section 3, section 4 and section 5. In the end, we summarize what we have founded with some discussions.

Holographic setup and relevant numerics
A simple holographic model for the two dimensional superfluid is a gravitational system in asymptotically AdS 4 spacetime, where dual to the conserved current j a and the condensate operator O in the superfluid, a dynamical U(1) gauge field A a and a complex scalar field Ψ with mass m are introduced respectively and minimally coupled by the charge q. After

JHEP07(2016)092
rescaling the matter fields, the corresponding bulk action can be written as [8,9] Here G is the Newton's constant, and the matter Lagarangian reads where D = ∇ − iA with ∇ the covariant derivative compatible to the metric. We are required to work in the regime with L 2 16πG 1 such that classical gravity is reliable, which corresponds to the large N limit of the dual boundary system. A black hole solution corresponds to placing the boundary system at the Hawking temperature. The bulk gauge field A t = µ at the AdS boundary corresponds to the chemical potential of the dual system. If the charge q and conformal dimension of the scalar operator O ± lie in certain range, taking the chemical potential sufficiently large will drive the bulk scalar field Ψ to condense via Higgs mechanism. As a result, the black hole carries a scalar hair outside the horizon, which corresponds to a superfluid phase on the boundary. The variation of action above will give rise to 16 equations of motion, including 10 from Einstein equation, 4 from Maxwell equation and 2 from Klein-Gordon equation. For simplicity, in this paper we shall ignore the backreaction of matter fields onto the metric, which can be implemented by taking the large q limit. Thus the background geometry is simply required to be a solution to vacuum Einstein equation with a negative cosmological constant. Since we are considering the superfluid at finite temperature, we take the Schwarzschild black brane solution as our background geometry, which can be expressed in the infalling Eddington coordinates as where the blackening factor f (z) = 1 − ( z z h ) 3 with z = z h the horizon and z = 0 the AdS boundary. As alluded to above, the dual boundary system is placed at Hawking temperature, which is given by On top of this background geometry, the equations of motion for the matter fields can then be written as In what follows, we will take the units in which L = 1, 16πGq 2 = 1, and z h = 1, and work with the case of m 2 = −2, in which the conformal dimension of the dual condensate operators O ± can be ∆ = 2, 1, corresponding to the standard and alternative quantization JHEP07(2016)092 respectively [16]. In the axial gauge A z = 0, the asymptotic solution of A and Ψ near the AdS boundary can be expanded as Then according to the holographic dictionary, the expectation value of j µ and O ± can be obtained explicitly as the variation of renormalized bulk on-shell action with respect to the sources. For the standard quantization, we have [17] where the dot denotes the time derivative,and the renormalized action is given by with the counter term added to make the original action finite. For the alternative quantization, one may naively think of ψ as the source, with φ the expectation value of the dual condensate operator. But this is not the case because we are working with the infalling Eddington coordinates rather than the Schwarzschild ones. Instead, it isψ = O + that plays the role of the external source for the dual condensate operator. With this in mind, we thus have [18] Although the resulting expectation value of condensate operator is still given by φ up to an unimportant overall minus sign, it is compulsory to keep in mind that for the alternative quantization the spontaneously broken superfluid phase is accomplished via Higgs mechanism not with the boundary condition ψ = 0, but with the source free boundary conditionψ at the AdS boundary. It is this boundary condition that we shall employ to evolve the bulk scalar at the AdS boundary. Vortices in superfluid are topological excitations with the quantized circulation. To see this, let us recall the definition of superfluid velocity

JHEP07(2016)092
With this, the winding number ω of a vortex is further defined as where c denotes a counterclockwise oriented path surrounding a single vortex. Note that close to the core of a single vortex with winding number ω, the condensate φ ∝ (z − z 0 ) ω for ω > 0 and φ ∝ (z − z 0 ) −ω for ω < 0 with z the complex coordinate and z 0 the location of the core. Whence the winding number must be an integer due to the single valued condensate. In addition, besides the vanishing magnitude of condensate at the core of a vortex, the corresponding phase shift around the vortex is given precisely by 2πω. This characteristic property will be exploited as an efficient way to identify quantized vortices in our later vortex counting.
To address QT at the alternative quantization by holography, we would first like to impose the following boundary conditions onto the bulk fields, i.e., on top of (2.14). In particular, we set the chemical potential µ = 2 > µ c with µ c ≈ 1.12 the critical chemical potential for the onset of a homogeneous and isotropic superfluid phase. Next we are required to prepare initial bulk configurations for the gauge field and complex scalar field at the Eddington time t = 0. We shall investigate superfluid turbulence in an 100 × 100 periodic box, where the initial 32 vortex-antivortex pairs with the winding number ω = ±6 are alternately placed at the square periodic lattice with lattice spacing b = 100/8. As detailed in appendix and illustrated in figure 1, it is not hard to construct an initial bulk configuration for the complex field Φ = Ψ z , dual to the aforementioned initial boundary state. Regarding the spatial components of gauge field A, for simplicity but without loss of generality, we shall simply take A = 0 as the bulk initial configuration. Furthermore, the initial configuration for A t can be determined by the constraint equation once we prescribe the second boundary condition ∂ z A t | z=0 = −ρ with ρ the boundary charge density. As such, we take the initial charge density ρ to be the value for the homogeneous and isotropic superfluid phase at the chemical potential µ.
With the above initial data and boundary conditions, the later time behavior of bulk fields can be obtained by the following evolution equations But in practice, to make our numerical simulation more reliable after a long time evolution, we adopt an alternative numerical scheme. Although we still evolve the scalar field Φ and JHEP07(2016)092 the spatial components of gauge field A in time by (2.19) and (2.20), we no longer use (2.21) to achieve the time evolution of A t except at z = 0, where (2.21) reduces to the conservation law for charge density as follows Then we can keep using the constraint equation (2.18) to solve A t from Φ and A at every step of the later time evolution. We numerically solve these non-trivial equations by employing the pseudo-spectral method with 28 Chebyshev modes in the z direction and 361 Fourier modes in the x direction, as well as the fourth order Runge-Kutta method in time direction with the time step ∆t = 0.05. We evolve the system for a total period of time t = 450, because it is expected that the intrinsic dynamics keeps unchanged without new interesting phenomenon popping up at later times. Finally, the vortex dynamics can be extracted from the near boundary behavior of bulk fields (2.7) by the holographic dictionary (2.11) and (2.12). Furthermore, we locate vortices by calculating the winding number (2.16) at each Fourier collocation point, which can be implemented by computing the total phase difference of φ around each plaquette formed by the nearest neighboring collocation points.

Vortex pair annihilation and decrease of vortex number
We have obtained 18 groups of data for the temporal evolution of superfluid turbulence at the chemical potential µ = 2. For the purpose of demonstration, we plot the configuration of turbulent superfluid at t = 400 in figure 2. The typical behavior can be described as follows. Initially, the lattice of vortices with winding number ±6 decays directly to the vortices with winding number ±1 in the absence of vortices with intermediate winding numbers. Meanwhile, due to the repulsive force between the vortices with like-winding number, not only do the resulting clusters of six ω = ±1 vortices rotate around the cluster centers but also expand away from them. After a very short period, such an expansion JHEP07(2016)092 makes the clusters of ω = 1 vortices meet the clusters of ω = −1 vortices, which triggers the vortex-antivortex pair annihilation. The coalescence of vortex and antivortex cores leads universally to the formation of a crescent-shaped depletion of condensate, which then converts into a shock wave, dissipated gradually in the superfluid. Thus the vortex number gets decreased. This cluster collision induced vortex pair annihilation proceeds till the neighboring clusters are well mixed. Then it takes some time for these well mixed clusters to get dissolved before the superfluid is full of randomly distributed vortices. During this cross-over process, the vortex pair annihilation rate increases gradually. When the vortices get randomly distributed, the vortex pair annihilation is driven essentially by the thermal effect, thus will be absent at zero temperature. At this moment, the vortex pair annihilation rate should also arrive at its maximal value and will keep constant till the end of our evolution.
To be quantitative, we would like to plot the temporal evolution of averaged vortex number over 18 groups of data in figure 3. As we see, the vortex number keeps unchanged when t < 25, and then gets decreased with time. Since the decrease of vortex number is induced completely by vortex pair annihilation, we expect to have where the vortex pair annihilation rate Γ(t), as alluded to above, is supposed to be time dependent. To see this explicitly, we plot the time variation of inverse of averaged vortex number in figure 4. Note that (3.1) leads to Thus figure 4 tells us that the vortex pair annihilation rate keeps invariant in the 25 < t < 110 region, which corresponds to the aforementioned cluster collision induced vortex pair annihilation. In the cross-over region, namely 110 < t < 240, the vortex pair annihilation rate increases gradually, which is followed by the purely thermal fluctuation driven vortex pair annihilation with a constant rate all the way to the end of our evolution.

Kolmogorov scaling law and onset of quantum turbulence
Kolmogorov's −5/3 scaling law in the kinetic energy spectrum is derived by dimensional analysis under the assumption that there is an inertial range during which the energy is passed from one side to another without substantial loss. Thus not only is it supposed to be obeyed by CT but also by QT. Actually this scaling law can be used as a characteristic signal for the onset of turbulence. To this end, we define the kinetic energy density as where ν = −φu. By a spatial Fourier transform, the total kinetic energy can be written as an integral over momentum, i.e., with (k, θ) the polar coordinates for the momentum k space.  Figure 5. The kinetic energy spectrum kin (t, k) at time t = 50, 100, 120, 150, 300, 450 for one of the 18 groups data at chemical potential µ = 2. The horizontal axis denotes log 10 k and the vertical axis denotes log 10 kin (t, k). The blue and orange lines are fitting results of the spectrum.
We plot the energy spectrum kin (t, k) at the different times for one of 18 groups of data in figure 5, where the horizontal axis is Log 10 k and the vertical axis is Log 10 kin (t, k). 1 As shown in figure 5, Kolmogorov scaling k −5/3 shows up in the cross-over region and persists till the end of our simulation. This implies that only in the cross-over regime does the system get driven to a turbulent state, which is consistent with our intuition that the early alternately distributed vortices should not be regarded as a turbulent state. We also list the inertial range for Kolmogorov scaling law and averaged vortex spacing at different times in table 1. As we see, the averaged vortex spacing falls well within the corresponding inertial region, which implies that Kolmogorov scaling law should be accounted for by the collective dynamics of turbulent vortices.
On the other hand, besides Kolmogorov scaling law, we also see k −3 scaling behavior, which appears after the initial vortices of winding number ±6 finish decaying into the JHEP07(2016)092 vortices of winding number ±1, and lasts all the way to the end of our evolution. This scaling has nothing to do with QT, but arises from the global rotation of vortices around the cluster centers at early times as well as the local rotation of superfluid around the vortex cores at all times. Note that the IR end of this scaling region corresponds to the scale for the above regular rotation. As time passes, this scale moves from the cluster size k ≈ 1, asymptotically to the vortex size k ≈ 3, signaling the transition from the dynamics dominated by the initial clusters of vortices at early times to the late state featured by randomly moving vortices.

Energy injection and direct energy cascade
To see the energy cascade direction for our superfluid, we shall inject energy into the system at the IR and UV scales respectively and see how the energy flows in the kinetic energy spectrum. In particular, we would like to inject energy by switching on the electric field on the boundary. Note that the injected power is given by [10,19,20] Thus we first choose the external source to be where η is a small constant and t 0 = 180 and g(t) is a gaussian function with width 8. We turn on the sources at t = 150 and turn them off at t = 210. Although the driven sources are not functions of x and y, the injected energy is not restricted at k = 0, instead distributed mainly at low ks, and even a little bit at high ks. 2 This is because there is a −∂ z A i term in the current of eq. (5.1), which depends on x and y in general. Also due to this term, the injected power is not positive everywhere. As a demonstration, we plot the space distribution of injected power at t = 170 and t = 190 in figure 6. As we see, not only is the injected power inhomogeneous, but also positive at some places and negative JHEP07(2016)092 at other places. Actually before t 0 the total injected power is positive while after t 0 the total injected power is negative. For example, the total power at t = 170 is 1928 while it is −1583 at t = 190. As a result, there is a net positive energy injected into the system by the electric field. Now we can see how the energy flows by comparing the evolution of the kinetic energy spectra between the driven and undriven systems with the same initial conditions, which is depicted in figure 7. As we expect, there is a positive energy injected to the turbulent superfluid at large scales, which then propagates deep into small scales and eventually results in an overall downwards shift of the energy spectrum at large scales and an overall upwards shift of the energy spectrum at small scales. This is a smoking gun for a direct energy cascade from the IR to the UV.
But nevertheless, to be reassuring, we also inject a net positive energy into the turbulent superfluid at high ks by the following external sources a x (t) = ηg(t − t 0 ) cos 250π 100 x + cos 320π 100 x , a y (t) = −ηg(t − t 0 ) cos 250π 100 y + cos 320π 100 y . (5.3) As depicted in figure 8, the injected energy at small scales dissipates away rapidly without affecting the kinetic energy spectrum at the IR. This is conceivable because the dominant dissipation mechanisms are expected to be vortex drag and vortex annihilation, both of which occur at small scales [10].

Summary and discussion
We have explored QT in finite temperature Bose-Einstein condensates by working with the alternative quantization in the simplest holographic model for the superfluid. By successfully accomplishing the long time stable numerical simulation of bulk dynamics in the infalling Eddington coordinates, we first investigate the temporal evolution of vortex number after shaking the lattice of vortices with winding number ±6. As a result, the vortex number keeps unchanged during the decay into the vortices of winding number ±1,  Figure 7. The comparison of the kinetic energy spectra kin (t, k) at time t = 165, 215, 240, 265 between the driven and undriven systems. The external sources are turned on at t = 150 and turned off at t = 210. The blue curve is the kinetic energy spectrum for the undriven superfluid while the red one is that for the driven superfluid. The injected energy at the IR is propagating to UV. and then decreases due to the cluster collision induced vortex pair annihilation, which is further followed by a crossover to thermal fluctuation driven vortex pair annihilation. By fitting the data with the N 2 decay formula, we find that the vortex pair annihilation rate keeps constant in the cluster collision induced vortex pair annihilation, increases in the crossover process, and keeps a constant again in the thermal fluctuation driven vortex pair annihilation. This process is also exhibited somehow by the right movement of the IR end of the region for the scaling k −3 in the kinetic energy spectrum. We further see Kolmogorov scaling law appearing in the crossover process, which can be regarded as the benchmark for the occurrence of the superfluid turbulence. At last we show that our superfluid has a direct energy cascade by injecting a positive energy into the system at both large and small scales.
Remarkably these results share the same features with those obtained in [10][11][12] for the standard quantization case, although it is suggested in [14,15] that the former resembles BEC while the latter resembles BCS. Consequently, it is natural to conjecture that the N 2 decay of vortex number, Kolmogorov scaling law and direct energy cascade should be universal patterns for the two dimensional superfluid turbulence at least at temperatures order of the critical temperature. On the other hand, it follows from the effective description of superfluid turbulence that QT at very low temperatures should have an inverse energy cascade with the N 5/3 decay behavior for the vortex number [21]. So it is intriguing to see how such a transition occurs when the temperature is lowered. But in order to achieve this The blue curve is the kinetic energy spectrum for the undriven superfluid while the red one is that for the driven superfluid. The injected energy at the UV is dissipated away rapidly, having no effect onto the kinetic energy spectrum at the IR.
by holography, one is required to take into account the backreaction of matter fields onto the metric. This is highly non-trivial and goes beyond the scope of this paper. We hope to address this issue in the near future and report our result somewhere else. namely we put the vortices onto the homogeneous and isotropic equilibrium configuration Φ eq (z) background, which can be easily obtained by numerically solving the ordinary differential equations for scalar field Φ eq (z) and time component of gauge field A t (z). The resultant Φ eq (z) is plotted in figure 9 for µ = 2.
For the single vortex with winding number ω, one can take the bulk configuration as Φ ω = g(r)e iωθ Φ eq (z), (A.1) where (r, θ) are polar coordinates. We require that the initial profile function g(r) → 1 for r → ∞, and g(r) → r ω for r → 0 to guarantee the smoothness of Φ ω at the core of vortex. The concrete form of the initial profile function g(r) will not affect the intrinsic vortex dynamics in superfluid, because the unwanted modes will be dissipated rapidly into the black hole. Then the bulk configuration for a variety of vortices can be readily constructed by multiplying Φ eq (z) simultaneously with g(r − r i )e iω i θ i . Regarding the lattice of vortices, one is generically required to further multiply the above bulk configuration with a random phase e iχ(x,y) to make vortices get moved. For the lattice of vortices in our periodic box, we take χ(x, y) = Re γ nΛ kx=−nΛ nΛ ky=−nΛ α(k)e ik·x , where γ is a small constant, n is a small integer with Λ = 2π 100 , and α(k) is a set of O(1) random complex coefficients.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.