Phase separation and exotic vortex phases in a two-species holographic superfluid

At finite temperature, the stable equilibrium states of coupled two-component Bose Einstein condensations (BECs) with the same conformal mass in both non-rotating and rotating condition can be obtained by studying its real dynamics via holography. The equilibrium state is the state where the free energy reaches the minimum and does not change any more. In the case of no rotation, the spatial phase separated states of the two components become more stable than the miscible condensates state when the direct repulsive inter-component coupling constant η>ηc>0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta>\eta _c>0$$\end{document}. Under rotation, the quantum fluid reveals four equilibrium structures of vortex states by varying the η\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta $$\end{document} from the miscible region to the phase separated region. Among the four structures, the vortex sheet solution is the most exotic one that appears in the phase separated region.


Introduction
The gauge-gravity duality [1][2][3] which relates strongly interacting quantum field theories to theories of classic gravity in higher dimensions has provided a new scheme, which can not only study strongly interacting condensed matter systems in equilibrium [4,5], but also study the real time dynamics when the system is far away from equilibrium [6][7][8][9][10]. Also, Schwinger-Keldysh propagators from AdS/CFT correspondence [11][12][13][14] has been studied in details in order to compute the real-time n-point functions from the partition function in the bulk. The first proposed model of a single component holographic superfluid/BEC was given in [15][16][17]. The array of vortices is known to happen in a rotating superfluid [18], which has been obtained in holography by studying the Wei-Can Yang and Chuan-Yin Xia contributed equally to this work. a e-mail: hbzeng@yzu.edu.cn (corresponding author) full dynamics of the single component holographic superfluid/BEC in a rotating disk as the final time-independent solutions [19,20], while single vortex has been obtained in the same system as a static solution in the presence of rotation without the dynamic process [21][22][23][24][25][26]. Besides the superfluid/BEC with only one order parameter, a two-component superfluids/BECs with two coupled order parameters have become one of the most concerned topics in condensed matter physics, since it demonstrates novel quantum states which can not be observed in a single component system as confirmed in the frame work of two-component G-P equations [27]. The two components Ψ 1 and Ψ 2 are coupled through the direct coupling term η|Ψ 1 | 2 |Ψ 2 | 2 . In the case of no rotation, the two-component BECs will enter a spatial separation state due to the repulsive interaction η > 0 between the two components [28][29][30][31]. In the presence of rotation, by solving the time dependent two-component Gross-Pitaevskii (G-P) equations where the two scalars have the same mass and the same intracomponent interaction, the equilibrium vortex structures can be obtained after a long time revolution, which reveal rich structures by varying the direct coupling between the two components [32]. As η increase the interlocked vortex lattices undergo phase transitions from triangular to square, to double-core lattices, and eventually develop vortex sheets [33]. Vortex sheet is a state that the vortices are aligned to make up winding chains of single quantized vortices, and the chains of two components are interwoven alternately, or form as the cylindrical vortex sheets. Such a solution was proposed for the first time by Landau and Lifshitz [34] and has been observed experimentally in 3 He A [35]. Vortex sheet has proved to be an important physical object with nontrivial topology, which has many physical applications [36].
Though the rotating single component superfluid has been studied thoroughly in holography, the holographic study on a rotating two component superfluid is still lacking. The pur-pose of this article is to study the equilibrium properties of strongly coupled two component BECs with and without rotation by using the holographic method. Indeed, continuous phase transition can not happen in 1D and 2D since the divergent fluctuation will destroy the long range correlation. However, the holographic model we used is in the large N limit while the fluctuations are suppressed, also the size of the system is much larger than the correlation length, we are not thinking of the holographic finite size effect [37] in the present work. The equilibrium states can be obtained by evolving the system from an initial homogeneous superfluid state under the disturbance of rotation or small fluctuation. Our treatment is very similar to the calculation of the time-dependent two-component G-P equation, where the static solution is also obtained as the final state that no longer changes with time anymore [32,33]. When there is no rotation, the phase separation states was studied in detail for different temperatures and different values of coupling constants. We found that the phase separation states tend to occur in the case of low temperature and relatively high inter-components direct repulsive interaction, while Josephson coupling always prevents the system from phase separation. Under rotation, by turning off Josephson coupling, the holographic superfluid reveals four vortex structures corresponding to different direct coupling values between two components. The four structures include triangle lattice, square lattice, vortex stripe and vortex sheet. Vortex sheet solution is the result of phase separation. The vortex diagram in the intercomponent direct coupling η versus rotation-frequency Ω is also be obtained, which is similar to the simulation results of the two-component G-P equation [33]. Because of the nonlinear nature of the bulk theory, we will have to rely on numerical methods by solving a highly nonlinear partial differential bulk equations of motion. This paper is organized as follows. Section 2 introduces the holographic superfluid model with two condensations, the equations of motion (EoMs) of the model and also the numerical method we adapted to solve the EoMs. Section 3 discusses the properties of the two-component system without rotation and the phase separation state is obtained with a large repulsive inter-components coupling. Section 4 presents the appearance of exotic vortex phases of the two-component system under rotation, and the vortex phase diagram is also obtained. The dynamic process of vortex evolution and corresponding free energy properties are also studied. Section 5 is devoted to the some discussions.

Holographic model and the time evolution equations of motion
The holographic model we use is a bottom-up construction contains two coupled charged scalar fields living in an Ad S black hole background, which is an extension of the single component holographic superfluid model proposed in [16]. The action include two charged scalar fields and one U (1) gauge field [38] the inter-component coupling potential between the two charged scalar fields takes the form where η is called the direct coupling. Note that the interaction potential takes the same form as the G-P equations for coupled two-component superfluid [32,39,40]. In the holography, there is another more complex model dual to a twocomponent superfluid with two gauge fields corresponding to two charged scalar fields respectively [41,42]. Also in the action, The metric of the Ad S 4 black hole background adopts Eddington-Finkelstein coordinates, in which is the AdS radius, z is the AdS radial coordinate of the bulk and f (z) = 1 − (z/z h ) 3 . Thus, z = 0 is the AdS boundary while z = z h is the horizon; the Hawking temperature is T = 3/(4π z h ). r and θ are respectively the radial and angular coordinates of the dual 2 + 1 dimensional boundary, which is a disk that suitable to study the properties of the two-component superfluid under rotation. The axial gauge A z = 0 is adopted as in [7,17]. Near the boundary z = 0, by choosing that m 2 1 = m 2 2 = −2, the general solutions take the asymptotic form as, j = 1, 2, the coefficients a r,θ can be regarded as the superfluid velocity along r, θ directions while b r,θ as the conjugate currents [21]. Coefficients a t and b t are respectively interpreted as chemical potential μ and charge density ρ in the boundary field theory. Moreover, Ψ 1 j is a source term which is always set to be zero, Ψ 2 j is the vacuum expectation value O j of the dual scalar operator in the boundary in the spontaneous symmetry broken phase. As the simplest case we set m 2 1 = m 2 2 , there is a discrete symmetry upon interchange of the two charged scalars. Such a discrete symmetry was also considered in the two-component G-P equations simulations [32,33].
We take the superfluid interpretation of the holographic model, where the rotation is introduced by imposing the Dirichlet boundary condition for gauge fields on the boundary as [17,22] where Ω is the constant angular velocity of the disk and a θ is the superfluid velocity along the θ direction. Since the system is only given rotational speed, and that the superfluid is confined to the container, so there is no radial velocity (a r = 0). In another word, the disk in the simulation should be static in the inertial frame, but the superfluid rotates relative to the disk with an angular velocity. Then the picture is that the disk is static while the superfluid is rotating, or equivalently we treat the superfluid as static while the disk is rotating. The radius of the boundary disk is set as r = R. The Neumann boundary conditions ∂ r h i = 0 are adopted both at r = R and r = 0, where h i represents all the fields except a θ . Without loss of generality we rescale = z h = q = 1. Therefore, by scaling Ψ j = ψ j z and using the axial gauge that A z = 0, the equations of motion (EoMs) can be written as The EoMs are solved numerically by the Chebyshev spectral method in the z, r direction, while Fourier decomposition is adopted in the θ direction. The time evolution is simulated by the fourth order Runge-Kutta method. The GPU computing is used to speed up the calculation. The initial state at t = 0 is always chosen to be a homogeneous solution for fixed η when there is no rotation, which can be obtained by solving the time independent EoMs by fixing a chemical potential μ with the Newton-Raphson method. The solutions are considered to be the equilibrium states when the norm of changes of all fields become smaller than 10 −5 for sufficient long time interval. Practically, a solutions at later time t is used to minus the solutions at t − δt, when the maximum of the change of the solution become smaller than 10 −5 for a time interval (δt = 100). In order to verify whether its stability, random perturbation of the order ψ * 10 −2 is added to the order parameter. Then the system starts to evolve and will return to the same stable state in a short time, which is the state before the perturbation was added. The solution can remain unchanged at t = 20,000 and beyond, so we think we attained a stable solution.
There is critical value μ c (η) above which the two scalar fields will condense with the same value as a result of the discrete symmetry between the two components. From numerics we found that μ c (0) ∼ 4.07, then the dimensionless critical temperature is T 0 c = 3 4πμ c (0) = 0.0587. Turning on the direct coupling η from − 0.5 to 0.5 will not alter the critical temperature, though a positive/negative η will reduce/increase the value of order parameter a bit. In all the dynamic simulation we set T = 0.82T 0 c , at which for every combination of − 0.5 ≤ η ≤ 0.5, the system always have a stable homogeneous solution with the two components have the same value of order parameter.

Phase separation and the effect of η
Phase separation is an inhomogeneous solution that the two condensates do not overlap spatially. The immiscible state can be more stable than the miscible state with a lower free energy, which is a result of the repulsive interaction between the two condensations when the positive η is large enough and can be seen from the potential Eq. (2). Such an inhomogeneous stable solution can be obtained by solving the full time dependent equation. At initial time the system is in a homogeneous and miscible state with Ψ 1 = Ψ 2 at fixed η, such a homogeneous state might be a metastable state which can evolve to a final stable inhomogeneous state, which will not change in time anymore, the real ground state. Without rotation, it is more convenient to adopt the Cartesian coordinate rather than the polar coordinate in the boundary theory, the metric is then To study the phase separation just in x direction, we take the assumption that all fields are functions of t, z, x. At the boundary, we use periodic boundary conditions. To quantitatively describe the spatial overlap of the two condensations, we define an integral where O h is the corresponding homogeneous order parameter at t = 0 and the length of the x-axis is L = 50. The system size L has no effect on the properties. If χ = 1, system shows no phase separation, and it stays in the initial homogeneous state. Otherwise, if χ 1, it would be fair to say the system shows phase separation. In the intermediate case, the system is partially phase-separated and partially phase-mixed. In Fig. 1a, we show that at three different temperatures, the overlapping integral χ as a function of η. The increase of χ along temperature when η = 0.5 indicates that the system is harder to enter the phase separation phase, agrees with the increased correlation length of both condensations when T is approaching T c . Two samples of O 1,2 of a separated phase are illustrated in Fig. 1b, c. We tried to interchange the solutions of Ψ 1 and Ψ 2, and then substituted them into the equations to continue the simulation. We found that the solutions are still stable, so the solutions we obtained satisfy the commutative symmetry in the action. Such an exotic immiscible state has been experimentally observed in a two- [43] and a three- [44] component quantum fluid. In addition, the phase separated state can also be obtained from a holographic first order phase transition in an inhomogeneous black holes [45,46].

Exotic vortex phases under rotation and phase diagram
In this section we are going to discuss the dynamics of the two-component BECs under rotation and also the phase diagram of final equilibrium vortex. We think of a physical process that begin with a stable homogeneous configuration without rotation corresponding to a fixed temperature (which can be obtained by solving the time independent EOMs by the Newton-Raphson method). Then, by applying the boundary condition Eq. (6), we increase the superfluid velocity along the θ direction to make the system suddenly rotate. The added velocity will drive the system to evolve and finally enter a stable state. To illustrate the stability of the long time static static states, one can also compute the free energy F of the system. We are working in the probe limit, which means that the fixed back ground black hole plays the role of a heat reservoir.
The temperature of the superfluid is fixed, then the free energy of the probe superfluid can evolve to the lowest value. At a fixed temperature, F can be computed from the renormalized on-shell action S ren , which consists of two parts is the term to remove the divergence near the boundary z = 0, in which γ is the determinant of the reduced metric. Therefore, the form of the renormalized onshell action is In the presence of rotation, focusing on the case when the holographic two-component superfluid is at a temperature T = 0.82T 0 c , development of vortex lattices from t = 0 to the final equilibrium state at t = 2000 for R = 20 , Ω = 0.1 and η = −0.2 are shown in the above two rows of Fig. 2a. The Fig. 2b shows the dynamics for the average value of order parameters O 1 and O 2 . After t 1500, the two averages are equal and do not change any more. The bottom row of Fig.2 shows the time evolution of the corresponding free energy (F − F n )/T , in which F n is the free energy in the normal state, i.e., Ψ j = 0. In the regime 30 t 100 the system is in a high free energy metastable state. At time t 100 the free energy begins to decrease, and at the same time the vortices begin to form at edge of the disk. After t 200, when vortices have fully formed and gradually move into the inner of the disk and align to the lattice structure, we can see the system is in a equilibrium state with constant lower free energy. Until the final square lattice is formed, the free energy always stays at this minimum, which confirms the stability of the final structure. Also in [47], the time dependent "free energy" was also plotted in the vortex dynamics in a type II superconductor under a periodic magnetical field.
Through the dynamic evolution of vortex under different coupling constants, we found several kinds of vortex phases. By increasing the inter-component interaction η from − 0.6, the system will approach the phase separated regime. Then the interlaced lattices are expected, where vortices in either component are filled by the other component. The sum of O 1 and O 2 is found to be an approximate constant with small fluctuation, this is the holographic realization of Thomas-Fermi distribution. In Fig. 3a-c, by increasing η from − 0.6 we see triangular lattices (− 0.6 ≤ η ≤ − 0.45), square lattices (− 0.45 ≤ η ≤ − 0.1) and vortex stripes (− 0.1 < η < 0.05). The square lattice is stable, presumably due to the fact that each vortex in one component can have all its nearest-neighbor vortices to be in the other component [48].These various structures are similar to that obtained by two component G-P equations [33].
In a classical turbulence, the vortex sheet is a thin interface across which the tangential component of the flow velocity is discontinuous. In quantum fluid, Landau and Lifshitz firstly proposed the vortex sheets scenario in rotating superfluid [34], almost at the same time when Feynman published his paper on quantized vortices in superfluid. A quantum vortex sheet solution is that the vorticity concentrated in line with the irrotational circulating flow stay between them. A typical picture of vortex sheet can be found in Fig. 1 of a review paper [36], where the vortices concentrated in circles with a uniform distance between the circles. As a novel quantum state, the vortex sheet had never been observed in superfluid 4 He. However, vortex sheet has been observed in chiral superfluid 3 He-A since it can be stable due to the confinement of the vorticity within the topologically stable solitons [35], which may has many physical applications. The other candidate to demonstrate vortex sheet is a quantum fluid with repulsive two order parameters, because the phase separation naturally provides a region of vanishing order parameters of one com-  Fig. 1a), we find the expected exotic vortex sheet solutions. The vortices of the O 1 are located at the region of the domains of O 2 component. This can be understood from the fact that the condensate of one component works as a pinning potential for the vortices in the other component due to the phase separation nature [33]. By forming vortex sheets, the condensate achieves remarkable phase separation compared to a lattice. Furthermore, the vortex sheets nearly uniformly fill the disk, and the distance d between the layers are equal.
According to the calculation by Landau and Lifshitz [34], the distance d between sheets is determined by the surface tension σ of the soliton and the kinetic energy of the counterflow (v n − v s ) outside the sheet, where v n is the normal fluid velocity and v s is the vortex-free superfluid velocity. In unit volume, the counterflow energy is 1 2 6 , and the surface energy is σ d , where ρ s is the superfluid density, which can be obtained from the conjugate current where A θ = 1 2 ρ 2 Ω and the denominator (∇arg[Ψ i ]) θ − A φ is the gauge-invariant superfluid velocity along the angular direction.
By the minimization of energy, one obtains We confirm that the formula also hold in a two-component superfluid, a sample result when η = 0.2 is plotted in Fig. 5. We show in Fig. 6 the profile of ρ s in the x-direction for the vortex sheet case η = 0.5. ρ s of one component has a constant value where the other one vanishes, we take the constant value of superfluid density in the Landau's equation, then σ can be written as σ = 1 3 ( ρ s 0.95 ) 3 . As another important properties of superfluid, the Feynman linear relation [18] between the excited vortex numbers and angular velocity in one component superfluid may can naturally generalizes to a two-component superfluid as where N represents the number density of vortex and M is the atomic mass of BECs. Since the discrete symmetry upon interchange of the two charged scalars, the vortices numbers of the two components are expected to be the same. In the holographic model we also investigated the validity of the Feynman relation in the two-component system and found no deviations from it, also the numbers for the two components are the same, N 1 ≈ N 2 . A sample result is given in Fig. 7 for η = 0.2, R = 20. By fitting we can get the value of the mass M 1 = M 2 = 0.6820. We also computed the cases where R is 15, 20, 22, 25, the linear Feynman relations between rotating velocity and vortices still holds and we got that M is equal to 0.67, 0.68, 0.7, 0.69 respectively. This means that the radius has no effect on the vortices number density as confirmed in the single component superfluid under rotation [19]. At last, we investigate the phase diagram of the vortex structures in the intercomponent direct coupling η versus rotation-frequency Ω as plotted in Fig. 8. The upper limit of the rotation frequency is set by Ω = 0.1 while the bottom limit is set by Ω = 0.05. The diagram show a transition from triangular lattices to square lattices and then stripe and sheet for all angular velocities. The triangle lattices locate at the region when the coupling constant η is negative, where the two components prefer to attract each other. Keep increasing η to zero, square lattices appear and the vortices arrange in line then the vortex stripes appear at η = 0. When η > 0.05, where is the region of phase separation, it always presents a sheet solution, that shows that the sheet solution under rotation is a result of phase separation.

Discussions
In order to study the properties of rotating multi-components BECs, a two-components Einstein-Maxwell-Higgs model defined in a fixed black hole background is used while the charged scalars and gauge fields are just a probe to the gravity background. The probe limit corresponds to neglecting the effect of the bulk stress-energy due to the scalar field and gauge fields on the space-time geometry. In the dual fluid, this corresponds to neglecting the effect of the superfluid condensate on the normal component of the fluid [7,8]. We don't consider the effect of the superfluid onto the normal fluid. So our hypothesis is similar to the G-P equation, in which vortex lattice and phase separation do not occur in normal fluid components. Furthermore, it is extremely important to extend present work to the back-reaction case, as long as we are going to study the how the dynamics change as we incorporate the back reaction of the superfluid on the normal component and also the change of condition for "phase separation".
The vortex phase diagram obtained from AdS/CFT correspondence in the simulation are very similar to that shown in the two-component GP equations [32,33], the similarities include also four vortex states and also the phase transitions from triangular lattices to square lattices, to double-core lattices, and eventually develop vortex sheets. The results can also related to some experiments, for example, in ref [48], the interlaced square lattice similar to Fig. 3b was observed. Also in the same experiment, the vortex core size of the interlocked are bigger than the one of single component, which we can also observe in Fig. 3. The sheet solution is expected in the highly separated region, which may can be observed experimentally in a two-component BECs with tunable intercomponent interactions, which can be deeply in a phase separate region [49,50]. The two species are of the same mass, then the vortex core size are the same as shown in Figa. 3 and 4. Using of different masses for the two components then the discrete symmetry upon interchange of the two charged scalars will not be reserved, one will realize a coexistence system of vortices with different vortex-core sizes, then the lattice structure shown in this work may be changed, which deserves to be studied in future.