Generation of vortices and stabilization of vortex lattices in holographic superfluids

Within the simplest holographic superfluid model and without any ingredient put by hand, it is shown that vortices can be generated when the angular velocity of rotating superfluids exceeds certain critical values, which can be precisely determined by linear perturbation analyses (quasi-normal modes of the bulk AdS black brane). These vortices appear at the edge of the superfluid system first, and then automatically move into the bulk of the system, where they are eventually stabilized into certain vortex lattices. For the case of 18 vortices generated, we find (at least) five different patterns of the final lattices formed due to different initial perturbations, which can be compared to the known result for such lattices in weakly coupled Bose-Einstein condensates from free energy analyses.


Introduction
Vortices play a central role in the non-trivial dyamics of superfluids. The study of vortex generation and vortex lattice formation in rapidly rotating superfluids, both experimentally and theoretically, has a long history, from that in the old liquid helium-4 or helium-3 to that in the modern Bose-Einstein condensates (BEC) or superfluid ultracold fermions. Quantized vortices and vortex lattices are firstly discovered in helium experiments where the superfluid is driven by a rotating bucket [1,2] (see [3] for a review). The interest for vortices in BEC comes from the fact that the theories are more tractable and size of vortex core in BEC is much larger than that in liquid helium [4]. As a result, a large quantity of experiments on vortices and vortex lattices are performed in BEC (see [5] for a review), among which the critical angular velocity is determined in [4] and formation of lattices with large numbers of vortices are performed in [6,7]. Most of the related theoretical works are based on the Gross-Pitaevskii (GP) equation [8][9][10] as well as its generalizations, which is a mean field description of weakly coupled BEC at zero temperature, for example, formations of vortex lattices are studied in [11][12][13].
The GP equation based method is relatively simple, and can reveal some features of vortex generation and vortex lattice formation observed in experiments. In [12], for example, after entering the BEC, vortices move towards the central area of the rotating JHEP02(2020)104 condensate and finally form a regular lattice, which can reveal the basic features of the lattice formation process. However, it is not difficult to find several disadvantages in these GP equation based methods. First of all, there is no dissipation in the original GP equation, so in order to yield the dissipation effect at finite temperature, a classical thermal cloud has to be introduced [12]. In the numerical simulations by this method, vortices seem to nucleate in the thermal cloud and then enter the BEC to form the lattice. But an explanation how the classical vortices transform to the quantum objects is absent. Second, the GP equation can describe the weakly coupled BEC well, but is not expected to be a proper desciption of liquid helium systems, which have strong interparticle interactions, as well as general ultracold fermion systems. 1 It is thus valuable to look for alternative methods to cope with the problem of vortex generation and vortex lattice formation. Applied AdS/CFT duality, or holography, just provides one of such alternative methods, which has been proved as a very useful framework to study superfluid physics at finite temperature (see e.g. [14,15]). Holography equates certain quantum systems without gravity to gravitational systems in a curved spacetime with one additional spatial dimension [16][17][18]. It describes strongly coupled systems by design and naturally incorporates dissipations by including a black hole on the gravity side. Most interestingly, even the simplest holographic superfluid model [20] shows some features of superfluid ultracold fermions [21,22]. All the above facts make holography a very promising tool to describe superfluids other than the zero temperature weakly coupled BEC. In fact, a number of holographic studies of the vortex are performed, where single vortex solutions are considered in [23][24][25][26], and vortex lattice solutions obtained from static equations can be found in [27][28][29].
In this paper, we use the simplest holographic superfluid model to study the vortex generation and vortex lattice formation in rapidly rotating superfluids. Compared with a requisite thermal cloud introduced in [12,13] or other dissipation terms in [11], no additional ingredient is required in our model. Actually, we first investigate the critical angular velocity from the linear instability indicated by critical quasi-normal modes (QNM) of the bulk AdS black brane, which is a good comparison with experimental work [4] and theoretical work [30]. Next, in the nonlinear time evolution, it is shown that the vortices appear at the edge of the superfluid system at first, and then automatically move into the bulk of the system, where they are eventually stabilized into certain vortex lattices. In particular, we carefully examine the final configurations when there are 18 vortices generated in total. In this case, we find (at least) five different patterns of the final lattices formed due to different initial perturbations, which can be compared to the known result for such lattices in weakly coupled Bose-Einstein condensates from free energy analyses. Moreover, the influences of radius, temperature and angular velocity on both the linear analysis and nonlinear evolutions are discussed. Experimentally, the setup in [7] is most relevant to our theoretical study, where superfluids already have angular velocities at the initial stage when the condensate forms.

JHEP02(2020)104
This paper is organized as follows. In the next section, we briefly introduce the holographic superfluid model that will be used. In section 3, we construct our initial configurations and consider the linear instability from the QNM on this initial background. Then, in section 4, the dynamical process of vortex generation and vortex lattice formation is studied numerically by the nonlinear time evolution, and the final vortex lattices with 18 vortices are investigated. Finally, the conclusion with some discussion provided in the last section.

Action of the model
The simplest holographic model to describe superfluids is given in [20], which consists of a complex scalar field Ψ coupled to a U(1) gauge field A M in the (3 + 1)D gravity with a cosmological constant related to the AdS radius as Λ = −3/L 2 . The action is: where D M Ψ = (∇ M − iA M ) Ψ, and κ 2 = 8πG with G Newton's gravitational constant.
A simple way to solve the model is considering the probe limit, where backreaction of the matter fields is neglected. To do so, we take the limit e → ∞. As a result, the matter fields live in a bulk spacetime that is fixed to be the Schwarzschild-AdS balck brane for the finite temperature system, where f (z) = 1 − ( z z h ) 3 is the blackening factor. The corresponding Hawking temperature is which is identified with the temperature for the boundary system. The equations of motion are Without losing generality, we choose m 2 L 2 = −2, so the asymptotic behaviors of the matter fields are In addition, the axial gauge A z = 0 is adopted in what follows. According to the dictionary of AdS/CFT duality, µ and ρ are the chemical potential and particle number density of the JHEP02(2020)104 boundary system, respectively. To describe superfluidity as a spontaneous U(1) symmetry breaking, we will switch off the source ψ − , and then ψ + corresponds to the condensate (expectation value of the dual operator).
The boundary system is parametrized by T /µ. When T is fixed, a phase transition happens at a critical chemical potential µ c , above which ψ + does not vanish and the system is in the superfluid phase [20]. We will work in the superfluid phase, so in the following we will choose µ > µ c .
3 Initial configurations and quasi-normal modes

Initial configurations
We will work in polar coordinates, so the metric of the Schwarzschild-AdS balck brane (2.2) becomes For numerical simplicity, we set z h = 1 and L = 1 from now on. Hence, µ c = 4.06 in our numerical computations [20]. We will prepare a rotating BEC as our initial state, which is implemented by taking all fields real and independent of θ: With ansatz given above, we can choose A r to be 0, which actually requires no radial current in the bulk [24]. Then the simplified equations of motion are where we have made the replacement Ψ = zψ. At the AdS conformal boundary z = 0, the boundary conditions are as usual with µ the chemical potential. For the superfluid system, a rigid rotation can be introduced by letting where Ω is the angular velocity. The system should have a finite size in the r direction, i.e. the radius, and in what follows we use R to denote this radius. The boundary conditions at the edge r = R of the system are taken as For detailed discussions of vortices in holographic systems, see [24,25]. In particular, the differences between superfluid vortices and superconductor vortices in holographic models JHEP02(2020)104  are discussed in [25]. At the center r = 0 of the polar coordinates, we impose the boundary conditions from the rotation invariance of the initial configuration. From the above equations of motion and boundary conditions, we can obtain the specific initial configuration numerically, once the chemical potential µ and the angular velocity Ω are given. For illustration, we plot the distribution of order parameter as a function of space coordinates (x, y) when µ = 5.5 and Ω = 0.168 for our initial states in figure 1, where we utilise different colors to denote values of order parameter. To make a good comparison, this color scheme will be preserved from now on for all figures of distributions of order parameters in figure 2, figure 4, figure 5 and figure 11. From now on, we choose R = 12 in all our numerical studies except for figure 3d, figure 3e, figure 6 and figure 10.
For simplicity, we will omit space coordinates axes when plotting distributions of order parameters in all following plots for order parameters. We replot the distribution of order parameter in figure 2a when Ω = 0.168 with coordinates ommited, and figure 2b shows distribution of order parameter when Ω = 0.22. Here, the distribution in figure 2a will be used as initial state for generation of lattice with 7 vortices, and the distribution in JHEP02(2020)104 figure 2b is initial state for generation of lattices with 18 and 19 vortices. Figure 2 shows that with the increasement of Ω, boundary values of order parameters are depressed, and this is the expected result when the rotating superfluid consists of no vortex.

Quasi-normal modes
Before performing the nonlinear real time evolution, we first make a linear analysis of the system, which corresponds to quasi-normal modes of the bulk system. From now on, we will switch to the ingoing Eddington-Finkelstein coordinates, so that the ingoing boundary conditions at the horizon are easily incorporated by regularity. Under the Eddington-Finkelstein coordinates, the metric (3.1) becomes Note that the original axial gauge A z = 0 in the Schwarzschild coordinates (3.1) is violated after this coordinate transformation. In order to preserve A z = 0 in the new coordinates, a U(1) gauge transformation should be performed after the coordinate transformation: As a result, A r = ∂ r λ (z, r) and so it no longer vanishes, while ψ should not be taken real anymore.
Since the background configurations constructed in the last subsection are time translation invariant and rotation invariant, the linear perturbations can be expanded as [31] where the coefficients p, q and a µ are all complex functions of (z, r). In particular, δψ is complex, so q (z, r) is independent of p (z, r). Equations for quasi-normal modes are shown in appendix A. At the AdS conformal boundary z = 0, Dirichlet boundary conditions should be imposed [31]: At r = 0, boundary conditions can be obtained by asymptotic analyses of the perturbation equations: Boundary conditions at r = R should be consistent with static case: On top of a background configuration constructed in the last subsection and for a given m, the quasi-normal frequencies ω are determined by the condition that the perturbation JHEP02(2020)104 equations together with the boundary conditions (3.14) and (3.15) have nontrivial solutions. Actually, for every m, we can obain a series of ω, among which the one with the maximal imaginary part determines the instability (if there is any ω with a positive imaginary part) of the m mode: this mode is more unstable if the imaginary part of this ω is larger. Figure 3a shows our results for the background configuration in figure 2a with Ω = 0.168, where we can see that the system is unstable against the perturbation of the 4 ≤ m ≤ 29 modes and the most unstable mode has m = 19. 2 A critical angular velocity Ω m can be determined for every m, which is just the lowest value of Ω that makes the m mode unstable. The critical angular velocities for different m are shown in figure 3b. We can see that the critical angular velocity for the system is Ω c ≈ 0.16, which is the lowest value among all Ω m s. If the angular velocity Ω < Ω c , the rotating superfluid system is stable and there will be no vortex generation. In addition, relations between Ω c and chemical potential µ as well as R are studied. In figure 3c, the critical angular velocity Ω c is calculated as a function of µ. Our result shows that Ω c grows with µ, which indicates that themal fluctuations contribute to instability of the system. In figure 3d, the critical angular velocity Ω c is calculated with R varying from 5 to 30, which shows that the critical angular velocity Ω c decreases with the rise of R. In order to understand this relation better, we calculate the linear velocity v at r = R in figure 3e. We find that compared with the scale of the system, the change of linear velocity ∆v is rather small. This rather small ∆v, as well as the result from nonlinear evolutions in the next section that the most obvious change of the system first appears at the edge, indicates that the instability of the system should be mainly determined by the linear velocity at the edge r = R.

Formation of vortex lattices
When the angular velocity Ω is above the critical value Ω c , the system is unstable. In real systems, there are random perturbations from environmental noises, thermal and quantum fluctuations, imperfect preparation of the initial configurations and so on. If the system is unstable, the tiny perturbations of the unstable modes will exponentially grow up at the early stage and then go beyond the linear regime. Will the instability when Ω > Ω c finally result in vortex generation and vortex lattice formation? To answer this question, which is very important but outside the scope of the linear analyses, we will slightly perturb the system from the initial configurations constructed in section 3, and see what will happen by numerically simulating the nonlinear dynamics of this holographic system. To be specific, we will set µ = 5.5 from now on.
Basically, we use the same strategy as in [15,32,33] for the numerical evolution here. The only difference and complication is the polar coordinates for this system, which has never been coped with in nonlinear time evolution of holographic superfluids. The full equations of motion in the polar coordinates are shown in appendix B.
We also impose the boundary conditions at z = 0 as At the edge r = R of the system, we set We do spectral expansion in the (z, r, θ) directions, but skip the coordinate singularity at r = 0, so no boundary conditions are needed there [38]. As well, under the ingoing Eddington-Finkelstein coordinates, no boundary conditions are needed at the horizon. To evolve the system, initial perturbations are necessary. In our time evolution, perturbations at t = 0 are imposed merely on ψ with the following form:

JHEP02(2020)104
Here, A m is the amplitude which should be small enough, and z preserves Dirichlet boundary condition for ψ. sin m πr 2R behaves as r m when r → 0, which is required by the analytic property of ψ, and preserves the Neumann boundary condition of ψ at r = R.
With the setup above, we can do the time evolution by the standard fourth order Runge-Kuta method. Here, we will mainly show our numerical results for a typical process of the vortex lattice formation when Ω = 0.168. The results are as follows: when t < 60, some distortions appear at the edge of the superfluid. These distortions come from instabilities we studied in section 3.2 and will lead to the generation of vortices (see for example at t = 30 in figure 4a). When t ≈ 60, vortices appear at the edge (figure 4b).
These vortices then move toward the center of the superfluid with the ongoing generation of new vortices at the edge. When t ≈ 120, the generation process ends, while the movement of vortices continues (figure 4c). Vortices' movement basically ends when t ≈ 360, at which the lattice forms (figure 4d). We continue the time evolution and find that the structure of the lattice does not change anymore (see figure 4e), which means that the lattice gets stabilized.
To make a good comparison, we show another typical evolution that consists of 19 vortices in figure 5. 3 The evolution process is similar, except for the final background order parameters. In both cases, vortices are arranged triangularly. A more complex lattice will be obtained if we study the case with more vortices in a larger system.
It is important to study the relation between the number of vortices and other conditions or parameters, such as the angular velocity and chemical potential, of the system. In order to mimic best the experimental situations, we choose our initial perturbations as eq.   vortex number is witnessed when Ω grows. This relation is consistent with the Feynman rules [3]. Moreover, the relation between the vortex number and µ is plotted in figure 6b, with fixed Ω = 0.11 and R = 18. We see the tendency that the vortex number decreases with the increase of µ. Overall, a larger Ω and smaller µ can generally increase the number of vortices.

Properties of dissipation
To make a vortex lattice eventually stabilized, energy dissipation is necessary. In holography, the energy dissipated in this process corresponds to the energy falling into the black hole [14,34,35]. To be precise, the energy dissipation in holographic systems is just the energy flux across the bulk black hole horizon: where T M N is the energy momentum tensor. Figure 7a shows the energy dissipation for the whole evolution. We can see that nearly all nonzero values of Q lies in the period 15 < t < 250, which includes our main process from the initial unstable configurations to the formation of vortex lattices. When t > 250, the energy dissipation decreases rapidly    and tends to zero. Figure 7b shows energy dissipation when the lattice consists of 19 vortices, which shares a similar behavior. As shown in figure 7, the maximum value of dissipation differs in different evolutions, so it is interesting to study its dependence on other physical quantities, such as Ω and µ. With other quantities, except for Ω and µ, fixed, we perform several different evolutions and plot the results in figure 8. In figure 8a, we fix µ = 5.5 and calculate the maximum dissipation when Ω changes from 0.17 to 0.22. In this plot, the maximum dissipation decreases with the increase of Ω. In figure 8b, we fix Ω = 0.168, and calculate the maximum dissipation when µ increases from 5 to 5.5. It is obvious that the maximum dissipation increases with the increase of µ. These results seem to be interesting, but the physical mechanism behind them is yet to be understood further.
In our time evolutions, the time when the first vortex appears can be viewed as the start of evolution, 4 while the time when dissipation reaches its maximum can be regarded   to mark the starting point that the system begins to stabilize. As a result, the period between these two moments, which we name as ∆t from now on, can reflect efficiency of the evolution. Similar to studies of maximum dissipation, we calculate ∆t when either Ω or µ is not fixed, and the results are shown in figure 9. In figure 9a, the tendency that ∆t grows with the increase of Ω is witnessed. In figure 9b, ∆t drops with the increase of µ, implying that a system with lower temperature involves more quickly.
As discussed in the first paragraph of this subsection, when the lattices stabilize, the dissipations generally reach rather small values, so it is reasonable to estimate the equilibration time by the period between the moment when the first vortex appears and that when the dissipation is small enough. In order to get a better result, we calculate the equlibration time by averaging over 12 evolutions with random initial perturbations for every fixed set of parameters [33], setting R = 6. More specifically, we take the time when the dissipation becomes 0.06 times its maximum as the equilibration moment to best match our numerical precision, and make the initial perturbations random by introducing random phases into different modes in the initial perturbation (27) (with all unstable modes of the same tiny enough amplitude) [15]. Our results are plotted in figure 10. In figure 10a, the equilibration time grows with the increase of Ω when Ω < 0.48, and then remains roughly constant when Ω > 0.48. In figure 10b, the equilibration time drops with the increase of µ, which is similar to the result in figure 9b.

Patterns of vortex lattices with 18 vortices
Vortex lattice has a perfect triangular shape when the number of vortices is infinite, while there are distortions if we consider a finite vortex number case. There are two possibilities. The first case is that the number is so special that interior vortices can form a triangular structure. In this case, there is only one kind of vortex lattice.   It has been calculated that when the lattice consists of 18 vortices, there are 7 different patterns from free energy analyses based on the GP equation [36]. Still, we obtain these lattices from time evolution with perturbations added as equation (4.4). It needs to be emphasized that when Ω is specified, different groups of A m s can result in different evolution processes. In our time evolutions, we obtain 5 patterns with fixed angular velocity (chosen to be 0.22) but different groups of A m s (see figure 11). In this plot, we name these different lattices "18-N", where N is the order of the corresponding pattern in [36].
We have not found the 3rd and 5th patterns from our time evolution. There exist two possible reasons. First, we differentiate lattices by comparing shapes of our lattices with results in [36]. In our model, these 2 absent lattices may not be so distinct that we can merely recognize them by calculating their free energies, which cannot be accomplished because of our current numerical accuracy. The second reason is that these absent lattices are not easy to be realized from time evolution with random initial perturbations, i.e. one may need very special initial perturbations to produce them.

JHEP02(2020)104 5 Conclusion and discussion
By holography, we have shown that vortices can be generated when the angular velocity of rotating superfluids exceeds certain critical values. In numerical simulations of the nonlinear dynamics, these vortices appear at the edge of the superfluid system first, and then automatically move into the bulk of the system, where they are eventually stabilized into certain vortex lattices. For the case of 18 vortices generated, we have found (at least) five different patterns of the final lattices formed due to different initial perturbations. Actually, these patterns can be recognized as five of seven such patterns obtained from free energy analyses based on the GP equation.
Compared with previous works or traditional methods, our study in this paper has the following advantages. First, in holography, the linear instability that can precisely determine the critical angular velocity is naturally captured by QNM of the bulk black hole. But, as far as we know, such efficient linear analysis in traditional methods is still lacking, though it can be done in principle. More importantly, our study yields all the essential physics of rotating superfluid systems just using the simplest holographic superfluid model without any ingredient put by hand, while the traditional methods cannot incorporate dissipation or finite temperature effects easily and naturally.
For simplicity, we do not consider the backreaction of the matter fields onto the bulk geometry in our holographic superfluid model. But for a full holographic duality, such backreaction should be taken into account [37]. It will be also interesting to investigate the vortex generation and vortex lattice formation in holographic superconductors. These topics are left for future exploration.

Acknowledgments
YT would like to thank Hua-Bi Zeng, Hai-Qing Zhang, Hong Liu and Shan-Quan Lan for useful discussions. He would also like to thank the Center for Theoretical Physics, Massachusetts Institute of Technology for the hospitality. This work is partially supported by NSFC with Grant No.11675015 and No.11975235. YT is partially supported by the grants (No. 14DZ2260700) from Shanghai Key Laboratory of High Temperature Superconductors. He is also supported by the "Strategic Priority Research Program of the Chinese Academy of Sciences" with Grant No.XDB23030000. HZ is also supported by the Vrije Universiteit Brussel through the Strategic Research Program "High-Energy Physics", and he is also an individual FWO Fellow supported by 12G3515N.

A Equations for quasi-normal modes
We solve problems of quasi-normal modes in Eddington-Finkelstein coordinates, so metric (3.10) is adopted. Linealize equations (B.1)-(B.5), substitute (3.12) and (3.13) into the JHEP02(2020)104 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.