Dynamical Phase Transition from Nonequilibrium Dynamics of Dark Solitons

By holographic duality, we identify a novel dynamical phase transition which results from the temperature dependence of non-equilibrium dynamics of dark solitons in a superfluid.For a non-equilibrium superfluid system with an initial density of dark solitons, there exists a critical temperature $T_d$,above which the system relaxes to equilibrium by producing sound waves, while below which it goes through an intermediate phase with a finite density of vortex-antivortex pairs. In particular, as $T_d$ is approached from below, the density of vortex pairs scales as $(T_d - T)^\gamma$ with the critical exponent $\gamma = 1/2$.


Introduction
In this paper we discuss a novel dynamical phase transition for superfluids. Non-equilibrium superfluid states obtained from external driving or phase/density imprinting often contain a large number of dark solitons [1][2][3][4][5][6][7], which are regions of low densities in a superfluid. The subsequent evolution of such a system and its relaxation to equilibrium are of great interest. We argue that the dynamical relaxation process undergoes a "phase transition" at some critical temperature T d . For T > T d , the solitons directly decay into sound waves which subsequently relax to equilibrium, while for T < T d an intermediate phase with a finite density of vortices emerges. See FIG. 1. The intermediate phase is unstable but could in principle be long-lived enough for interesting dynamics including turbulence to happen. The density n v of vortices, which may be considered as the "order parameter" of the intermediate phase, is proportional to the density of original dark solitons and increases as the temperature is decreased. Near T d , it has critical behavior (1) While we are dealing with a non-equilibrium process, nevertheless a sensible temperature can be defined which measures the fractions of normal and superfluid components, and characterizes the dissipation of the system. This is also usually how temperature is determined in experimental situations (see e.g. [8,9]). The phase transition is deduced by examining in detail the temperature dependence of the decay of dark solitons in two classes of theories, which encompasses most systems currently under study in the literature. One class is the dissipative Gross-Pitaevskii equation (GPE) [10][11][12], which is a phenomenological theory with a dissipative parameter κ added "by hand" to the GPE. This theory is rather crude and requires significant modeling. The other class is holographic superfluid. Holographic duality equates certain strongly correlated systems of quantum arXiv:1810.11424v4 [hep-th] 22 Jan 2020 matter without gravity to classical gravitational systems in a curved spacetime with one additional spatial dimension. The gravity description provides a first-principle description of finite temperature dissipative effects; as we will see below once the "microscopic" theory is fixed, all aspects of the superfluid phase are determined. There is no phenomenological modeling involved. Remarkably, despite that these two classes of theories are governed by completely different sets of equations, and have very different physical origins, we find qualitatively similar results.
A dark soliton can decay by self-accelerating uniformly and then turning into sound waves, or by fragmenting into vortex pairs or filaments, a process called snake instability. We find that for T > T d it decays dominantly by self-acceleration but by snake instabilities for T < T d . We show that this pattern persists with a finite density of dark solitons, which results in the phase transition. We also present a general analytic argument for the "critical exponent" 1 2 in (1). Instabilities of dark solitons in various superfluids were studied extensively in the literature at theoretical level (see e.g. [13][14][15][16][17]), mainly based on the GPE or Bogoliubov-de Gennes (BdG) equations [18], which are suitable for zero or very low temperatures. A main technical advance which enables us to find the transition is a systematic study of temperature dependence of the dynamics and decay mechanisms of dark solitons (see also [19][20][21][22][23][24][25][26][27] for other finite temperature studies).

Setup and dark soliton solutions
The dissipative GPE for a superfluid can be written as where we have done some rescalings. The parameter µ specifies the equilibrium value of the condensate ψ 0 = v ≡ √ µ. For definiteness we will consider two spatial dimensions with spatial coordinates x = (x, y). κ > 0 is a dissipative parameter. To see its effect, note that where F is the free energy functional of the system Both κ and µ should be considered as temperature dependent. We expect κ to increase monotonically with temperature and µ to decrease with temperature. In particular µ should approach 0 as the critical temperature for a superfluid is approached from below. Below we will take κ as a proxy for temperature, i.e. when we say increasing the temperature for the dissipative GPE, we mean increasing κ. At the moment there is no first principle to decide the precise κ-dependence of µ, which may be considered as a major defect of this model.
The static dark soliton solution to (2), which we will denote as ψ S , can be readily obtained analytically. For definiteness we will take it to be translation invariant along y direction and centered at x = 0, i.e., ψ S (x) has the following properties: Let us now turn to holographic superfluids in 2 + 1 dimensions, which can be described by an Abelian-Higgs model in a (3+1)-dimensional anti-de Sitter (AdS 4 ) black hole spacetime [28,29]. The Lagrangian can be written as with the background spacetime metric In (9), L is the curvature radius of AdS and f (z) = 1 − ( z z h ) 3 with z = z h the location of an event horizon and z = 0 the AdS boundary. The black hole spacetime (9) has a Hawking temperature T = 3 4πz h , which is identified with the temperature of the dual boundary system. In (8) the U (1) gauge field A a is dual to a conserved current J a for a U (1) global symmetry in the boundary system, and the complex scalar field Ψ is dual to a boundary order parameter ψ charged under the U (1) symmetry. The system is in a superfluid phase below some critical temperature T c , when Ψ develops a normalizable profile in the bulk spacetime which corresponds to ψ developing a nonzero expectation value [28,29]. We will ignore the backreaction of A a and Ψ to the background black hole geometry, an approximation which works well when the temperature is not too low [30]. We are mostly interested in the regime T /T c > ∼ 0.3 where the approximation is sufficient.
For definiteness we will take the mass square of Ψ to be m 2 = − 2 L 2 . In this case, there are two possible boundary conditions for Ψ, leading to two different types of superfluids which will be denoted respectively as ψ ± . Dark soliton solutions ψ S in the ψ ± -superfluid phases can be found by solving equations of motion following from (8) with appropriate boundary conditions [31][32][33][34]. They also satisfy (6)- (7). Furthermore, it was observed there that a dark soliton in the ψ + -superfluid (ψ − -superfluid ) has features resembling those of a soliton in the BCS regime (BEC regime) of ultracold atomic gases. See FIG. 2 and its caption. While holographic superfluids cannot be directly mapped to either BEC or BCS regimes, the parallels are nevertheless striking. For ease of comparison, we will below refer to ψ ± -superfluids as BCS and BEClike superfluids, although one should keep the caveats in mind. 48. The behavior of dark solitons and vortices has been known to exhibit many differences in the BEC and BCS regimes of ultracold atomic gases [35][36][37]. Interestingly, these different features also appear in ψ± holographic superfluids [31][32][33][34].
We emphasize that once the parameters and boundary conditions of the gravity description (8) are fixed, the system is fully specified at all scales, including the full set of properties of the superfluid phase such as the critical temperature, temperature dependence of the order parameter, dynamics of dark solitons, and so on.

Linear instability analysis
Now consider small perturbations around a dark soliton solution. The calculations are very different for the dissipative GPE and holographic superfluids, but the conclusions are qualitatively similar. We will describe the main results below, leaving calculation details to Appendix B.
Since the system is translationally invariant along the y-direction, at the level of linear analysis, we can decompose perturbations in terms of Fourier modes e iqy in the y-direction. More explicitly, we write with a small parameter. Due to the parity symmetry y → −y, q and −q have identical behavior, so we will restrict to q ≥ 0. The equations for f ω,q (x) can be reduced to an eigenvalue equation with the eigenvalue proportional to ω. Since ψ S is odd in x, one finds that the even and odd parts of f ω,q (x) decouple. The eigenvalues ω e (q), ω o (q) (respectively for the even and odd parts with appropriate boundary condition on f ωq ) have a discrete spectrum and are all complex. Note that an ω-eigenvalue with a positive imaginary part leads to exponential time growth in (10) and thus corresponds to an unstable mode. At a given temperature T , one finds that there exists a q c (T ) such that there is exactly one unstable mode for each q ∈ [0, q c (T )) and the mode is even in x. The unstable mode turns out to be pure imaginary ω (0) e (q) = iλ q and λ q > 0, which leads to e −iω (0) e (q)t = e λqt . We will denote the eigenfunction f ω,q (x) for the unstable mode as f The upper value q c (T ) for the unstable mode decreases as one increases the temperature, and for holographic superfluids q c (T ) → 0 as T → T c . More explicitly, we find that as See FIG. 3. For the dissipative GPE, due to lack of understanding how µ and κ depend on T , a precise plot is not possible. Nevertheless, assuming that We will now show that: (i) the unstable mode at q = 0 corresponds to self-acceleration; (ii) the unstable mode at q = 0 corresponds to snake instability; (iii) when a soliton decays via snake instability, there is a vortex-anti vortex pair to be produced for each wavelength. For example, consider the system in a finite periodic box with length R y along the y-direction. The allowed qs are thus of the form q = 2πN Ry with N an integer. A snake instability with q = 2πN Ry will then create N pairs of vortices. To see (i), let us note that the center x c of the dark soliton can be identified as being located at the minimum of the condensate, i.e. ∂|ψ| 2 ∂x | xc = 0. Now consider (10) with f ω,q given by f q=0 (x)e λ0t , which leads to It can readily seen from the above equation that We thus see that the dark soliton accelerates exponentially. Now let us consider (10) for a general unstable q, with equation (12) becoming We then have (15) and θ = argf (15), x c depends on y sinusoidally, leading to a snake instability. In particular, for each wave length 2π q there are two points of zero velocity with opposite circulations around them, corresponding to the locations of a vortex and an anti-vortex to be created. This thus demonstrates (ii) and (iii) above.

Switch of dominant decay channel
Among all the unstable modes q < q c (T ), the one with the maximal λ q grows fastest and is thus expected to be the dominant decay channel of a soliton. In FIG. 4, we plot λ q as a function q for both holographic superfluids at two different temperatures. We notice that the plots for the higher temperature has maximal λ q at q = 0, so the dominant decay channel is the self-acceleration, while for the lower temperature plots, the maximal value of λ q occurs at some nonzero q, so the dominant decay channel is the snake instability. In FIG. 5 we plot the value of q max (T ) for which maximal value of λ q occurs as a function of temperature. We see that q max (T ) decreases with temperature monotonically until a value T d < T c after which q max becomes identically zero.
To study this phenomenon for the dissipative GPE, we now have to make assumptions on the κ-dependence of µ. We find that the qualitative behavior does not depend sensitively on the choice of µ(κ). For example, by choosing µ to be independent of κ or to be a linear decreasing function of κ, the very similar behavior is obtained. We present results for the latter in FIG. 6 and FIG. 7. We thus find that the sharp transition also occurs for the dissipative GPE.
As T → T d from below we find the following "critical behavior" with exponent γ in all cases (holographic superfluids and dissipative GPE) numerically around 1 2 . While our current numerical method is not sufficient to make a very accurate determination of γ, there is a simple Ginzburg-Landau type argument which shows that γ should be given by 1 2 generically. Consider expanding λ q (T ) in small q Since for T ≈ T d , q max is close to zero, and the above expansion should suffice for determining the behavior of q max near T d . For T > ∼ T d , q max is zero, we thus should have a(T ) > 0, while for T < ∼ T d , we should generically have a(T ) < 0, b(T d ) > 0 in order to have q max = 0. We thus conclude that near T = T d we can expand a(T ) as It then follows that From our earlier discussion of self-acceleration and snake instability, we thus conclude that there is a sharp transition at T d in how dark solitons decay: for T c > T > T d , dark solitons decay by self-acceleration, while for T < T d , they decay by snake instability. In particular, in the snake instability regime, the number of vortex and anti-vortex pairs created by the decay of a dark soliton, which is proportional to q max , increases as we decrease the temperature.
The value λ qmax characterizes the growth rate of instabilities and thus its inverse may be considered as giving the time scale for the "lifetime" of a dark soliton. In  FIG. 8 we plot the temperature dependence of λ qmax for the three types of superfluids. It is interesting to note the higher temperature, the instabilities grow slower and thus the longer lifetime of a soliton.

Dynamical phase transition
Now consider a superfluid system in a non-equilibrium state with a finite initial density of dark solitons. Then there is a non-equilibrium dynamical phase transition at T d in how the system relaxes back to equilibrium. For T > T d , the solitons directly decay to sound waves which subsequently equilibrate. For T < T d there exists an intermediate phase with an initial density n v of vortexantivortex pairs, which is proportional to q max . The lower temperature, the larger q max and n v . Furthermore, from (16) as T d is approached from below n v has critical behavior To confirm our expectation that the decay of a dark soliton is controlled by the mode q max with the largest λ q , and the dynamical phase transition, we now consider the full nonlinear evolution of initial states with dark solitons.
Let us first consider the case of a single initial dark soliton under the following two types of perturbations: (i) a single-wave length perturbation δψ ∝ e iqy with q ∈ [0, q c (T )), and (ii) a general perturbation of the form δψ ∝ q e iαq e iqy with random phases α q . For (i), with q = 0, one finds that different parts of the dark soliton accelerate uniformly. During its acceleration, the dark soliton also broadens, and eventually dissolves into sound waves. For q = 0, the acceleration pattern for different parts of the soliton shows sinusoidal behavior with wave length 2π q , consistent with the prediction of (15). In particular, the vortex-antivortex formation as well as their numbers are precisely as predicted below (15). See FIG. 9  (a) (b). For a general perturbation (ii), the resulting evolutions are shown in FIG. 9 (c) (d). We indeed find that at a temperature above T d , the soliton accelerates without any vortex formation (FIG. 9 (c)). At a temperature below T d , the evolution is dominated by snake instability, evidenced by formation of vortex pairs. Furthermore, the number of vortex pairs is consistent with linear analysis. For instance, for ψ + -superfluid at T Tc = 0.78, the linear analysis tells us that the mode q = 2πN Ry with N = 3 is the most unstable mode. We find the non-linear evolution indeed produces 3 vortex pairs (see FIG. 9 (d)).
Finally we consider the full nonlinear evolution of a system with a finite initial density of dark solitons. On physical ground, we expect our results regarding the stability of a single soliton should apply to a finite density of solitons as far as the density is not too high, i.e. as far as the average distances between solitons are larger than the healing length of a soliton (which is often treated as a microscopic scale). This expectation is indeed confirmed by full nonlinear simulations, see FIG. 10 for contrast of two temperatures below and above T d . Thus nonlinear evolutions fully confirm the existence of a dynamical phase transition.

Discussion
To summarize, by studying the decay of dark solitons, we identified a novel dynamical phase transition which results from temperature dependence of non-equilibrium dynamics of solitons. The phase transition involves the appearance of an intermediate (unstable) phase with the order parameter given by the density of vortex pairs and a critical exponent γ = 1 2 . This vortex gas/liquid phase could have important implications for understanding non-equilibrium dynamics of a superfluid, for example, it could be long-lived enough to exhibit vortex and wave turbulent states. It would be extremely interesting to search for such a dynamical phase transition in experiments, and to measure the critical temperature T d as well as the critical exponent γ. The fact that we observe the phase transition in such vastly different theories as the dissipative GPE and holographic systems gives much confidence in that it is a universal phenomenon. A natural place to look for such a phase transition is ultracold atomic gases which we expect should exhibit this phenomenon both in the BEC and BCS regimes. The transition should be already detectable in the usual experimental setup of a conventional harmonic potential. Recent experimental advances in the realization of box-like optical traps [38] can make the comparisons with our theoretical results even more direct.
At a technical level, our approach to decay of dark solitons has a few important new features compared with earlier studies [13][14][15][16][17][19][20][21][22][23][24][25][26][27], Our linear instability analyses not only provide a unified treatment for self-acceleration 1 and snake instabilities, but also give a wealth of information on temperature dependence of dominant decay channels, the number of vortex pairs produced, as well as time scale for the lifetime of a soliton. Our full nonlinear evolution of the decay of a dark soliton at finite temperatures provide strong support for the conclusions of linear analyses. Nonlinear evolution of a dark soliton was studied previously in [17] using the timedependent BdG equations of the BCS-BEC crossover at zero temperature, where snake instability was observed, but not self-acceleration. This is consistent with our general conclusion that at low temperatures snake instabilities should dominate. For simplicity but without loss of generality, we will focus only on the case of m 2 L 2 = −2 in the axial gauge A z = 0, in which the equations of motion on top of the background (9) can be written in an explicit way as [40] [41]. The asymptotic solution of A and Φ can be expanded near the AdS boundary as Then the expectation value of J and ψ ± (the condensate of ψ ± -superfluids) can be explicitly obtained by holography as the variation of the renormalized bulk on-shell action with respect to the sources, i.e., where the dot denotes the time derivative, and S ± is the renormalized action, obtained by adding the corresponding counter term to the original action to make it finite and well posed for the variational principle as [42] In [31,32], the dark soliton was found in the Schwarzschild coordinates by the finite difference scheme with the relaxation method. Here we will instead reconstruct such a dark soliton solution in Eddington coordinates (9) by the pseudo-spectral scheme with the Newton-Raphson method. This turns out to be much more efficient, and enables us to study lower temperatures. To achieve this, we make the following ansatz for the non-vanishing bulk fields, Then the equations of motion reduce to Associated with the superfluid condensate ψ ± , the dark soliton solution can thus be obtained by imposing the boundary conditions at the AdS boundary, as well as the Neumann boundary conditions at x = ± Rx 2 , where R x is set to be much larger than the characteristic length of the dark soliton in consideration, so that the finite size effect can be removed.
On the other hand, for our purpose, the fully non-linear simulation starts with the initial data where the subscript S denotes the corresponding profile for the static dark soliton solution. For simplicity but without loss of generality, we take δΦ = ze iqy or a random superposition of modes with different qs, and δA = 0 [43]. Then A t can be solved by the constraint equation (A3) subject to the boundary conditions at the AdS boundary. Next, we can evolve Φ and A through Eq.(A1) and Eq.(A2) subject to the source free boundary condition at the AdS boundary, and the boundary conditions at x = ± Rx 2 . Furthermore, we can evolve the second boundary condition in Eq.(A17) by evaluating Eq.(A4) at the AdS boundary, i.e., which is essentially the conservation law of the boundary particle current. The later time behavior for A t can be obtained in the same way as described before. The above evolution scheme is implemented numerically by the fourth order Runge-Kutta method along the time direction, together with the pseudo-spectral method along the space directions, where Chebyshev Polynomials are used in the z and x directions, while Fourier modes are used in the y direction. This gives a computational advantage over the previous holographic investigations, which deal exclusively with periodic boundary conditions along both the x and y directions [40,[44][45][46][47] and therefore cannot properly accomodate the dynamics of dark solitons. free boundary condition at the AdS boundary. In addition, we impose Dirichlet boundary condition for δA and Neumann boundary condition for δΦ at x = ± Rx 2 . Then these boundary conditions, together with the above other perturbation equations, can be cast into the form [A(q) + B(q)ω]δ = 0 with δ the perturbation fields evaluated at the grid points by the pseudo-spectral method. Then the quasi-normal modes can be obtained numerically by solving this generalized eigenvalue problem.
Note that the dark soliton background solution has even parity for A t and odd parity for (Φ, A x ) with respect to the x axis. By inspection of the above linear perturbation equations on top of the dark soliton background, one can see that the quasi-normal modes with odd parity for (δA t , δA y ) and even parity for (δΦ, δA x ) decouple from the modes with opposite parity. So to reduce our computational resource for quasi-normal modes, we can restrict our computational domain onto [0, Rx 2 ] in the x direction by imposing Dirichlet or Neumann boundary condition at x = 0, depending on whether the parity of the perturbed fields is odd or even. As a demonstration, we plot the spectrum of low lying quasi-normal modes in FIG. 11. As we can see, only the even parity branch of δΦ gives rise to the unstable mode, which is purely imaginary.
The strategy for the linear perturbation analysis of dissipative GPE is exactly the same, but much simpler.
Here we present only the corresponding linear perturbation equation where we have also decomposed the condensate wave function to its real and imaginary parts as ψ = a + ib and assumed that the perturbation behaves like δ(x)e −iωt+iqy .