Dynamical evolution of spinodal decomposition in holographic superfluids

,

First (1st) order phase transitions, including for example liquid-gas (evaporation) and solid-liquid (melting) phase transitions, are ubiquitous in nature.Commonly, a 1st order phase transition point is defined as the location at which two phases share the same value of the thermodynamic potential (e.g., free energy).1st order phase transitions usually occur either from meta-stable states by jumping barriers, or from unstable states following a Cahn-Hilliard linear instability at finite wave-vector.The dynamics related to jumping potential barriers from meta-stable states are known as "bubble nucleation".They involve a thermodynamic barrier to phase separation, and they are usually described by Classical Nucleation Theory (CNT) [1,2].On the other hand, the dynamics arising around unstable states, and the consequent time-dependent processes, do not involve any thermodynamic barrier nor nucleation step and they give rise to the so called "spinodal decomposition" [3].In both processes, inhomogeneous configurations appear and phase separation is often observed with the production and evolution of bubbles, localized regions of the favoured phase inside the thermodynamically disfavoured phase (see Figure 1 for a schematic representation).
The main subject of this work is spinodal decomposition, defined as the phase separation triggered by a linear instability in the spinodal region of a first-order phase transition.On general grounds, spinodal decomposition is well described by Cahn-Hilliard theory [4].In order to explain the essence of this framework, let us consider the evolution of a system with a conserved charge density c within the mean-field Ginzburg-Landau approximation.At leading order, the free energy can be written as a function of the charge density c (or component concentration) and gradient, as follows [5]: where V is the volume of the system and f a scalar thermodynamic function.Here, κ controls the cost in the free energy to produce spatial variations in the concentration c, i.e., the energetic cost of inhomogeinities.With further constraint from the conservation of the total charge, Q = c dV , spinodal decomposition can be described by a generalized diffusion equation: where M is the mobility of the system and the chemical potential µ is defined as µ = ∂F/∂c.The variation of this free energy functional yields to the Cahn-Hilliard equation for the charge density perturbation which can be expressed in the following form [5]: By going to Fourier space e −i(ωt−kx) , the linearized Cahn-Hilliard equation can be easily solved in terms of the imaginary part of the frequency: Let us note that the susceptibility χ is the curvature of the thermodynamic potential f .We can see from this formula that, when the susceptibility χ is negative (e.g., a maximum of the potential f ), there is a critical wave-vector: FIG. 1.The routes towards phase separation in a first-order phase transition.On the left, the dynamics of spinodal decomposition from a local maximum of the thermodynamic potential F and the corresponding time-dependent dynamics, well described by the Cahn-Hilliard equation.On the right, the dynamics from a metastable local minimum of the potential which evolve via the jump of a thermodynamic barrier indicated with the vertical dashed line.The picture shows an example of bubble nucleation in which bubbles of the thermodynamic favoured phase build inside the thermodynamic disfavoured one, producing phase separation.
below which the solution in Eq.( 4) results in positive values of Im [ω(k)], indicating a linear instability at large length scales (or small k).If the system is large enough, the solution in the region with χ < 0 will spontaneously undergo the growth of these inhomogeneous perturbations -spinodal decomposition -which ultimately results in phase separation.The Cahn-Hilliard theory has been widely used to model spinodal decomposition phenomena from unstable initial states [6][7][8] at which (see Fig. 1) the thermodynamic susceptibility is negative.
In recent years, holography has been found to be a convenient and effective tool for studying strongly coupled field theories at finite temperatures [9][10][11], especially out of equilibrium [12].Through holography, numerical experiments were conducted to study the non-equilibrium physics and phase transition processes in several setups (for a non complete list see Refs.[13][14][15][16][17][18][19][20][21][22][23][24]).Among the various setups, the holographic superfluid model [25] has emerged as one of the most attractive platforms to study non-equilibrium physics.To mention some of the main results within this model, non equilibrium processes were realized in Ref. [13], which also confirmed dynamically the onset of the phase transition indicated by the linear stability analysis.The dynamics of vortexes and turbulence in the holographic superfluid are studied in Ref. [14,26], which found that the superfluid kinetic energy spectrum obeys the Kolmogorov -5/3 scaling law.The Kibble-Zurek mechanism describing the formation of vortexes is also well tested in the holographic superfluid model [17,27], and the breaking of its universality in the fast quench processes has interestingly been reported in a recent study [28].
Back to our main topic, the spinodal decomposition phenomenon has been recently discussed in holographic models [29][30][31][32][33][34][35][36][37].In Refs.[29,30], the quasinormal modes (QNMs) of a holographic system with 1st order phase transitions are calculated, providing solid evidence for the linear instability related to spinodal decomposition.The full dynamical evolution in the same model has been then considered in Ref. [31].The spinodal decomposition triggered by the Gregory-Laflamme instability has been studied in Ref. [32] and the relation with hydrodynamics was discussed.In Ref. [33], it has been proposed that the time evolution of spinodal decomposition consists of four generic stages: a first, linear stage in which the instability grows exponentially; a second, non-linear stage in which peaks and/or phase domains are formed; a third stage in which these structures merge; and a fourth stage in which the system finally relaxes to a static, phase-separated configuration.Further studies [34][35][36][37] also realized spinodal decomposition in various holographic models and studied problems such as the formation and collision of domain walls as well as effects of different values of surface tension.Importantly, the above studies only considered an effective one-dimensional problem, whereas the realization of spinodal decomposition with two effective spacial dimensions was studied only in Ref. [38] with a special focus on the emission of gravitational waves.In all these studies, the holographic setup for the 1st order phase transition was an AdS-Einstein-scalar model, which is usually applied to study the 1st order QCD phase transition [39].On the contrary, the bubble nucleation process, involving a finite potential barrier, has been realized in both one and two spacial dimensions in the holographic superfluid model with two competing s-wave orders in Ref. [15], in which the surface tension was calculated and verified.However, the full dynamical processes which accompany spinodal decomposition have not been studied in a holographic superfluid model yet.This will be our main motivation and the main novelty of this work.
Very recently, the linear stability of a simple holographic superfluid model with various kinds of phase transitions including 2nd order, 1st order and 0th order, was studied in detail within the probe limit [40].There, linear instabilities at large length scale, with 0 < k < k c , were discovered in the spinodal region of the 1st order superfluid phase transition.However, differently from the diffusive dynamics in the Cahn-Hilliard theory presented above, the instabilities in this superfluid model turned out to have a different origin.Indeed, it is well known that a superfluid presents dynamics which are richer than the conservation of a single charge, due to the emergence of Goldstone degrees of freedom arising because of the spontaneous symmetry breaking of the U(1) symmetry.In particular, its hydrodynamic description cannot be correctly formulated ignoring this fact.In a relativistic superfluid, the U(1) charge does not diffuse anymore, as envisaged in Eq.( 2), but it rather couples to the Goldstone mode (i.e., the fluctuations of the phase of the superfluid order parameter) creating a propagating wave-like excitation known as 2nd sound [41].In the broken phase, the dispersion of this mode at low wave-vector is given by and has been studied in detail in the holographic superfluid model [42][43][44][45][46]. Hence, in a superfluid, the instability which leads to spinodal decomposition arises because the second sound speed square c 2 s becomes negative.When that happens, the sound dispersion takes the following form which implies the presence of a critical wave-vector: below which: A negative c 2 s is a consequence of a negative thermodynamic susceptibility χ < 0. Indeed, from superfluid hydrodynamics it is well known that c 2 s ∝ χ −1 [42,44].Although this instability arises from the growth of the superfluid sound mode, rather than a diffusive mode as considered in Eq.( 2), the similarities with Cahn-Hilliard theory are striking.First, both cases concern an inhomogeneous instability which develops only below a critical wave-vector.Moreover, in both cases the instability is in a sense thermodynamic since it corresponds to a thermodynamic susceptibility becoming negative.
In summary, in the first-order holographic superfluid model, we do expect an instability with Im(ω) > 0 for 0 < k < k c .This behavior is reminiscent of the well-known Gregory-Laflamme type instability of black strings [47], and it is qualitatively similar to the sound mode instabilities observed in the previous holographic studies which do not involve a superfluid system [30][31][32][33][34][35][36].Based on these progresses, our aim is to extend the previous analyses by considering the real-time dynamics of spinodal decomposition in the holographic superfluid model of [40].For simplicity, we will work in the probe limit, where there are less numerical challenges and the formation and evolution of 2-dimensional round bubbles can be observed and tracked in an easier way.
The structure of the paper is as follows: in Section II, we introduce the holographic model and briefly review the static solutions and the QNMs of the model.In Section III, we perform the "numerical experiments".We test the linear stability of homogeneous initial states in the spinodal region, and we confirm the theoretical expression for the critical length scale.Moreover, we perform quenching experiments to show the full dynamical processes of spinodal decomposition in the dual (2+1)-dimensional field theory.Finally, we provide some comments on the question of whether the final state is uniform or not.In Section IV, we summarize and present an outlook of our results.

II. HOLOGRAPHIC SETUP AND LINEAR STABILITY A. Holographic setup and the 1st order phase transition
We consider the holographic s-wave superfluid model with nonlinear self-interaction terms in (3+1) dimensional asymptotic AdS spacetime [40], which is a simple extension of the HHH model [25].The bulk action of this model is given by Here, F µν = ∇ µ A ν − ∇ ν A µ is the field strength of the U(1) gauge field, and D µ Ψ = ∇ µ Ψ − iqA µ Ψ the standard covariant derivative.Λ = −3/L 2 is the negative cosmological constant, where L is the AdS radius.Ψ is a complex scalar field charged under the U(1) gauge symmetry.The additional two terms appearing in the matter Lagrangian, Eq.( 12), λ(Ψ * Ψ) 2 and τ (Ψ * Ψ) 3 , introduce nonlinear self-interactions for the bulk scalar field Ψ.
The influence of the nonlinear term λ(Ψ * Ψ) 2 has been investigated analytically in Ref. [48] and numerically in a holographic s+p model in asymptotic AdS 4 spacetime in Ref. [49].In Ref. [40], both the two nonlinear terms λ(Ψ * Ψ) 2 and τ (Ψ * Ψ) 3 were considered and the corresponding QNMs spectrum analyzed.Thanks to these additional terms, the phase diagram of the model becomes very rich and it includes 2nd, 1st as well as 0th order phase transitions.
In the rest of this paper, for simplicity, we fix m 2 = −2, q = 1, z h = 1, L = 1.The values of the two nonlinear parameters are set to λ = −2, τ = 0.8 and λ = −0.8,τ = 0.28, in order to focus on two typical 1st order phase transitions.In the rest of this section, we briefly review previous results concerning the thermodynamics and the linear stability of this model for the choice of parameters λ = −2, τ = 0.8, which were already presented in Ref. [40].
In the probe limit, the background geometry is given by a (3+1) dimensional black brane In these notations, z is the radial coordinate in the bulk, with z = 0 the location of the AdS boundary and z = z h the horizon.The Hawking temperature of this black brane solution is given by For homogeneous static solutions, we consider the following ansatz for the matter fields: which results in the following equations of motion where the prime denotes the derivative with respect to z.
To solve these coupled equations, we need to specify the boundary conditions.The expansions for the bulk fields near the horizon z → z h are given by On the contrary, near the AdS boundary z → 0 the bulk fields are of the form The free energy difference between the normal solution and the superfluid solution.The black triangle denotes the "critical point" at which the broken phase connects to the normal one.The green cross denotes the turning point of the condensate curve in the canonical ensemble, while the red cross denotes the turning point of the condensate curve in the grand canonical ensemble.The black square denotes the location where the amplitude mode collides with another pure imaginary mode (see [40] for details).The solid red line is the unstable branch of solutions, while the solid blue is either stable or meta-stable.Finally, the dotted blue part is stable under homogeneous perturbations but unstable under inhomogeneous ones.
where µ and ρ are the chemical potential and charge density of the boundary system, respectively.We follow the standard quantization scheme, in which ψ (1) is regarded as the source term for the boundary operator while ψ (2)  is regarded as the vacuum expectation value.We set the source free condition ψ (1) = 0 to obtain the solutions corresponding to the spontaneous symmetry breaking of the boundary global U (1) symmetry.
Using the AdS/CFT dictionary, the free energy coincides with the Euclidean on-shell action with proper boundary terms, This thermodynamic potential involves only the matter contribution because of the probe limit and it is valid only for homogeneous solutions.Next, we show the results for the choice of parameters λ = −2, τ = 0.8, whose linear dynamical stability has already been presented in Ref. [40].The superfluid condensate, the imaginary part Im(ω) of the amplitude mode as well as the free energy with respect to the charge density ρ are plotted in the left, central, and right panels of Figure 2, respectively.The other case, λ = −0.8,τ = 0.28, will be consider later on.
The condensate at the critical point first grows leftwards up to the turning point (green cross) and then changes direction, rather than directly growing rightwards at the critical point as in second order phase transitions.Notice that the imaginary part of the frequency of the amplitude mode is positive along the red branch, implying the linear instability of those solutions.The free energy curve forms a "swallow" tail shape and has discontinuous first derivatives at the phase transition point, which is defined as the intersection point of the black and blue curves in the free energy.All these features clearly demonstrate the presence of a first order phase transition in this system.

B. Quasinormal modes (QNMs) and linear stability
Taking advantage of the gauge/gravity correspondence, the linear stability of the dual quantum system under infinitesimal perturbations can be studied by analyzing the QNMs of the classical gravity dual.In this subsection, we review the results for the 1st order phase transition with λ = −2 and τ = 0.8 obtained in Ref. [40].In the QNM analysis, if the imaginary part of the frequency ω of one mode is positive, this mode will grow exponentially in time ∼ exp(|Im(ω)|t) destabilizing the original background solution.Therefore, only when the imaginary part of ω of all the modes is negative, the background solution is linearly stable.In other words, all the possible linear perturbations will not grow up in time but eventually relax back to the background solution.In general, the complex frequency ω also depends on the real wave vector k.Homogeneous perturbations corresponds to k = 0 and constitute only a sub-set of all possible kind of linear perturbations.In order to test completely the linear stability of the system around a certain background solution, the full spectrum of ω(k) has to be considered.
The results for the imaginary frequency of the amplitude mode with k = 0 are shown in the central panel of Figure 2.They show that under homogeneous perturbations, the superfluid solution is unstable below the turning point of the condensate curve (red branch in the left panel of Figure 2), and stable above (blue branch in the left panel of Figure 2).This fits perfectly with the analysis from the free energy landscape, right panel of Figure 2.
The results for the QNMs at finite k are more interesting.We select six different solutions marked by different colored points in the left panel of Figure 3.In the corresponding central panel, we show the imaginary part of the lowest mode as a function of the wave-vector k for each of these six solutions.We observe that the solutions marked with black and blue points in the left panel are stable for any value of the wave-vector k, since the imaginary part of this mode (and indeed of all other modes) is negative.On the contrary, the solutions marked by cyan, magenta, red, and violet colors are unstable.More precisely, for the red and purple points, the corresponding solution is stable under homogeneous k = 0 perturbations but unstable under perturbations with finite wave-vector.Finally, the solutions marked with magenta and cyan colors in the left panel of Figure 3 are unstable under both k = 0 and k ̸ = 0 perturbations.As already anticipated in [40], this is closely related to the difference between the canonical and grand canonical ensembles.To better illustrate this point, we plot the condensate curves in the two ensembles with solid and dashed curves respectively in the left panel of Figure 3, where we use dashed green lines to denote the two turning points for the two different ensembles.Importantly, these horizontal dashed lines correspond the onset of linear instability at k = 0 and k ̸ = 0 respectively.The canonical ensemble agrees well with the linear stability analysis at k = 0 but not with the stability analysis under inhomogeneous perturbations, which is on the contrary captured by the grand canonical ensemble.We can see this explicitly by looking at the purple solution between the two different curves.In the canonical ensemble the derivative of the condensate is positive, indicating thermodynamic stability.On the contrary, in the grand-canonical ensemble the derivative is negative hinting towards a thermodynamic instability.At the same time, the QNM analysis reveals that such a solution is stable at k = 0, as one could derive from the canonical ensemble, but unstable at k ̸ = 0, as hinted by the grand-canonical ensemble.For more details on this point see [40].
It is worth noticing that the imaginary part of the unstable mode crosses the horizontal axes at a critical value k = k c and the mode becomes stable again at k > k c .Using the numerical data, we derive a critical length scale l c = 2π/k c .This scale sets the onset of the instability driven by the inhomogeneous perturbations.Importantly, these unstable modes are sound modes described by Eq.( 6), which share the same qualitative feature with the Gregory-Laflamme instability of black rings.The study of the QNMs also tell us that, at large enough length scales, the solutions between the two turning points (between the two green dashed lines in Figure 3) will be destabilized by the growth of unstable inhomogeneous modes -spinodal decomposition.However, the QNM analysis does not tell us anything about the non-linear time dependent dynamics nor the endpoint of the instability.This will be the topic of next section.

III. NUMERICAL EXPERIMENTS AND NONLINEAR DYNAMICAL EVOLUTION
In this section, we simulate the full dynamical processes.Let us start with a brief overview of the setup used and the numerical methods.

A. Holographic setup for the dynamical tests
To simplify the numerical calculation, it is convenient to use in-going Eddington coordinates In this choice of coordinates, the component A z of the U(1) gauge field is nonzero.Nevertheless, we can set A z = 0 using the gauge symmetry, at the cost of having a complex valued scalar field Ψ.For general time dependent solutions, it is consistent to use the following ansatz where ⃗ x includes space directions x and y, and ψ(z, ⃗ x, t) is a complex valued function.
From the lagrangian, Eq.( 12), we get the time dependent equations of motion In order to solve these equations, we consider a finite volume boundary system in a square box with periodic boundary conditions in the x and y directions.The dependence on these two directions is expanded in Fourier spectrum with ∆x = 0.25.The z dependence of the functions is expanded using a Chebyshev spectrum with 21 grid points.Finally, we apply the fourth-order Runge-Kutta (RK4) method to promote the time evolution with a time step ∆t = 0.05.
Besides the periodic boundary conditions in the x and y directions, we also need to specify the boundary conditions along the z direction.Eq.( 27) is a constraint equation; we take its value on the conformal boundary as the boundary condition for the bulk field ϕ along the z direction.This condition is equivalent to the conservation of U(1) charge on the boundary system.Finally, we impose Dirichlet boundary conditions for ψ and A µ at z = 0.After specifying the equations of motion and boundary conditions, we are ready to run the numerical experiment and observe the full time dependent evolution of this holographic superfluid system along the spinodal decomposition.

B. Testing linear stability with dynamical processes
In the first stage of the numerical experiments, we test the linear stability results obtained from QNMs, by setting an initial state and see whether it stays still or go away to another configuration in a very long time.From this setup, the critical scale defined in the previous section from the wave-vector at which the mode becomes unstable is also confirmed.2, and other points can be found in the left panel of Figure 3.

Homogeneous perturbations
We start testing the stability of initial homogeneous states under homogeneous perturbations, which should verify the results of the QNMs with k = 0.In this group of numerical experiments, we set the time dependent fields as well as perturbations to only depend on t and z.We set the initial condition to be a static solution and let it evolve following the non-linear dynamics.
We choose six typical static homogeneous solutions marked by the different colored points in the left panel of Figure 3 as the initial states of the dynamic evolution.We show the value of condensate as a function of time in the right panel of Figure 3.For convenience, we also show the detailed information for the six typical solutions used as initial states in the "numerical experiments", as well as the turning point in the canonical ensemble and the grand canonical ensemble marked by the green cross and red cross and the critical point marked by the black triangle in Table I.
We can see from the right panel of Figure 3 that the red, violet, blue, and black curves are horizontal lines with no time dependence.This implies that the four related solutions marked by the red, violet, blue, and black points in the left panel of Figure 3 are therefore stable upon linear perturbations.For the solutions marked by the magenta and cyan points, the dynamical evolution shows that they are unstable under very small homogeneous perturbations.Indeed, after a finite amount of time, the two solutions run away from the initial state.It is interesting to notice that there are two different paths for each initial state, which also end in different final states.This is because the initial unstable states are saddle points in the potential landscape, where on one side is the way "down the hill" towards the (meta-)stable normal phase solution, while on the other side is the way "down the hill" towards the (meta-)stable superfluid solution on the upper branch of the condensate curve (see inset in the right panel of Figure 3 for a cartoon).Indeed, at late time, one of the two branches reaches a final state with zero condensate, while the other a final state with a constant and finite condensate, different from that of the initial state.We have further confirmed that the difference between the two paths, and the corresponding final states, is related to the sign of the perturbations of the scalar field around the initial states.In addition, we observe that when the final state is a superfluid state with finite condensate, the cyan and magenta curves show an oscillating behavior while approaching the equilibrium final state.This phenomenon is consistent with the results of QNMs of the final state.In particular, it is simply the sign that the late time relaxational dynamics are dominated by the lowest QNM -the superfluid sound -(see [50] for similar findings).These oscillations are not present when the time evolution ends in the normal phase, since the latter does not contain any sound mode in the spectrum (in the probe limit).
The above numerical experiments show that the solutions marked by the red, violet, blue, and black points are stable under homogeneous perturbations, while the the solutions marked by the magenta and cyan points are unstable, which is consistent with the results of the k = 0 quasinormal modes, as well as the potential landscape analysis in the canonical ensemble.

Inhomogeneous case
To continue, we now focus on the numerical experiments regarding the stability of the initial homogeneous states under inhomogeneous perturbations, which will be compared with the results from the QNMs at finite k.In this set of numerical experiments, for simplicity, we still turn off the y dependence, which means the dynamical fields as well as perturbations depends only on t, z, x.
As before, we take the six typical static homogeneous solutions marked by different colored points in the left panel of Figure 3 as the initial states of the dynamic evolution (see Table I for details), and show the average value of condensate as a function of time in Figure 4.The left and right panels correspond respectively to l x = 4 and l x = 10, where l x is the size of the system along the x direction.
We can see that the left panel of Figure 4 is almost the same as the right panel of Figure 3.This is because the length scale of the system l x = 4 is very small and indeed below the onset of instability shown in the central panel of Fig. 3.In order to make this clearer, we use the relation between the wave vector and the wave length, k = 2π/l, to calculate the two values of k x related to l x = 4 and l x = 10.Here, k x is the minimum wave-vector allowed in a box of size l x .These two values are marked with dashed red vertical lines in the central panel of Figure 3.The dashed right line corresponds to l x = 4 and the left one to l x = 10.We clearly see that the value of k x corresponding to l x = 4 is larger than all the critical wave-vectors.In other words, since only modes with k > k x are allowed in this system, no unstable modes can appear.As a consequence, the solutions marked by the red and violet points are stable.However, the finite size does not rule out the homogeneous mode with k = 0 since all allowed wave-vectors are given by 2nπ/l x , (n = 0, 1, 2...).Hence, the solutions marked by the cyan and magenta points are still unstable.Now, let us look at the right panel of Figure 4, where the length scale l x is larger and part of the finite k unstable modes for the solutions marked by the cyan, magenta, and red points are allowed.The purple solution is stable because the system is too small to contain the corresponding unstable modes (see central panel of Figure 3).As shown by the red curve, the solution marked by the red point is not stable and it decays into the normal phase.The magenta curve is the same as the left panel of Figure 4, while one of the cyan curves shows a new process approaching a final state with a smaller average value of condensate.We expect that this new final state is inhomogeneous.In order to confirm this, we plot the local value of condensate as a function of the x coordinate in the inset of the right panel in Figure 4, where we see a clear x dependence in this final state.The inset of the left panel in Figure 4 shows the case of the final state of the upper cyan curve in the left plot of Figure 4 as a comparison.
The 1-dimensional non uniform final state shown in the inset in the right panel of Figure 4 is not a typical 1dimensional bubble state.More precisely, the region inside of the bubble still presents a finite condensate and therefore is not in the normal phase.In order to better explore the formation and evolution of bubbles, in the next subsection, we try different initial states and we study the dynamical quenching processes in 2-dimensional boundary space.Before that, let us discuss in more detail the critical length and perform more dynamical tests to confirm the value of this critical length l c calculated from the critical wave vector k c of the QNMs.
In order to do so, we select 11 solutions in the range of ⟨O⟩ ρ t < ⟨O⟩ < ⟨O⟩ µ t as initial states and repeat this second set of numerical experiments.For each initial solution, we set the size of the box along the x direction to different values, and see whether the infinitesimally deformed initial solution stays stable or run away to another final state in a long time t max = 150000.The results are shown in Figure 5.
In the left panel of Figure 5, we show the dependence of the imaginary part of the unstable mode on the wave vector k for the 11 solutions.The intersection of the curves with the horizontal axes indicates the critical value k c , which is related to the critical length scale l c = 2π/k c of the system, beyond which the solution becomes unstable.Our results show that the solution which is closer to the turning point in the grand canonical ensemble has a larger value of the critical wave vector k c , and therefore a smaller value of the critical length scale l c .At exactly the turning point in the canonical ensemble, the critical wave-vector diverges and the finite size system becomes stable.
In the right panel of Figure 5, we show the results of the dynamical test.The magenta triangles indicate the length scale calculated from the critical wave vector of the unstable mode from QNMs shown in the left panel of Figure 5, and the red curve is an interpolation of the numerical data.The blue asterisks show the length scale at and above which the initial solution runs away before t max , while the blue cross show the length scale at and below which the FIG. 6.The critical scale lc near the turning point ρ µ t in the grand canonical ensemble.The magenta triangle represents the critical scale calculated from QNMs, while the red curve is the fitting result performed using Eq.(32).initial solution does not evolve in time in the range 0 < t < t max .For practical purposes, we set the time scale t max = 150000.Near the critical length, the growth of the unstable modes takes an extremely long time.We can see that both the asterisk and cross are very close to the magenta triangle (inset in the right panel of Figure 5), showing a perfect match between the results from quasinormal modes (QNMs) and this set of dynamical tests.
When the solution approaches the turning point in the grand canonical ensemble, ρ → ρ µ t , the critical value of the wave vector becomes zero while the critical length scale diverges.It is interesting to fit this divergent behavior with a power-law function In order to get better results, we add additional 29 data points by calculating the QNMs in the region close to the turning point ρ µ t , and present these data points (the magenta triangles) as well as the fitting curve (the red line) in Figure 6.
Fitting the data points with the three free parameters R, ρ µ t and η, as in Eq.( 32), we obtain Interestingly, η is approximately 1/2, which is the typical value of the critical exponent in mean field theory.

Inhomogeneous final states with phase separation
In the previous numerical experiments, we have confirmed the non uniform instability as well as the existence of a critical length scale which can be predicted by the quasinormal mode dispersion.Another interesting question left to be investigated is whether the final state of the time dependent evolution triggered by the non uniform instability is also non uniform.
In principle, in order to answer this question, we would have to follow the full dynamical process up to the final state and then check whether the latter is homogeneous or not.However, because of charge conservation, it is possible to use the charge density of the uniform initial state to locate the position of the final state.With the full information of the QNM spectrum, we then know whether the possible uniform final state is stable.In case all the uniform solutions with the value of charge density equal to the initial one are unstable, then the real final state must necessarily be inhomogeneous.We will show such an example of this sort in the next subsection.
However, if we find that the final states are linearly stable, we are not able to conclude that the final state is homogeneous.Indeed, the final state might be a meta stable non-uniform phase which is linearly stable but nonlinearly unstable.In this case we have to run the full time dependent evolution to reach the final state and confirm this explicitly.Fortunately, within the holographic setup, we are able to run the full non-equilibrium evolution and analyze the nature of the final state.This is what the next subsection will be about.We anticipate that, only in some cases, the initial solutions with inhomogeneous instability end up into a uniform final state.In fact, we also find examples where the system reaches a non uniform final state at the end of the dynamical evolution.
We start by considering as initial state the solution marked by the purple point in Fig. 3 and set the box size to l x = l y = 100.We then add an infinitesimal random perturbation to the initial state and let it evolve.For completeness, we have considered both 1D and 2D systems.The behavior of the condensate as a function of time during the full nonlinear time evolution is shown in Figure 7 for the two cases.There, we clearly observed that during the time evolution bubbles are created.They then decrease in number and grow in size until the final state is reached.The final state contains a large symmetric bubble whose interior region is in a superfluid state with a very large value of the condensate.This implies that the final state is not homogeneous but displays evident phase separation between two regions (indicated in blue and yellow colors): one, inside the bubble, in the superfluid phase; the other, outside the bubble, in the normal phase with vanishing condensate.
From the real time evolution of the 1D and 2D systems, we see that these phase separation process, known as spinodal decomposition, is naturally triggered by the inhomogeneous linear instability appearing in the dispersion relation of the superfluid second sound mode.Interestingly, the whole dynamical evolution can be divided into four stages.In the first stage (or the early stage), the inhomogeneous unstable perturbations grow up exponentially as described by the linear analysis.In the second stage, the now large perturbations deform the homogeneous initial state and induce the formation of initial bubbles and domain wall structures.We see that the number of initial bubbles are determined by the scale of the wavelength of the most unstable mode.The allowed discrete values for the wave-vector k in this finite system with l x = 100 are The maximum peak of the purple curve in the central panel of Figure 3 is located at k m = 0.3119 and indicates the most unstable wave-vector.We see that k m ≈ k 5 = 5k 1 and therefore the most unstable mode is the one with wavevector equal to k 5 .The characteristic scale of the initial inhomogeneous state right after the onset of the instability is then given by 2π/k 5 , which is 1/5 of the length scale of the whole system.Therefore, in the second stage we see approximately 5 initial bubbles in the 1D system and 5 × 5 = 25 in the 2D one.This agrees very well with what observed in Figure 7.In the third stage (or middle stage), the initial bubbles undergo a very diversified evolution.Some small bubbles shrink and disappear, while some big ones expand.Bubbles can also merge with each other creating a larger one.However, in most cases, when the chemical potential become balanced at each point, there is only one big bubble left in the final state.When only one big bubble is left, the system enters into the fourth stage (the final stage), in which the evolution of the big bubble becomes very slow.We notice that, occasionally, there might be two or more bubbles with the same radius left in the final state and when the radius of all bubbles become equal, the system goes into the fourth stage.The four stages just described are consistent with what suggested previously in Ref. [33].

C. Quench experiments
We have already successfully realized phase separation from spinodal decomposition and analyzed its evolution in real time.However, in the previous "experiments", we have always started from unstable initial states.The condensate curve in the canonical ensemble (see left panel in Fig. 3) reveals that it is possible to quench the system from a stable initial state towards an unstable state in the spinodal region.This might reflect a more realistic scenario which can be possibly compared with experimental setups.
Therefore, in this section, we perform additional numerical simulations to show the time evolution in these quench "experiments".In the first case with λ = −2, τ = 0.8, we quench the charge density of the system uniformly from a stable uniform initial state with ρ i = 3.5005 to the charge density value of the solution marked by the violet point in Figure 3.The results of the 1-dimensional numerical simulation are shown in Figure 8, while those of the 2-dimensional numerical simulation are presented in Figure 9.
In Figure 8, we plot the time dependent value of the average charge density (solid black line) and average condensate (solid purple line) in the left panel, while we plot the time dependent value of local condensate in the right panel.We can see that in this quench experiment, the quench speed is fast, and the system first goes through the boundary of stability uniformly.Quickly after the end of the quench, the system first stays nearly uniform for a while, and then the spinodal decomposition appears.The domain walls emerge in the system due to the uneven disturbance.They then either disappear or grow until reaching a final equilibrium state.
In Figure 9, we plot the 2D snapshot of the local condensate at four typical instants, where the creation and evolution of bubbles are evident.The quench speed is the same the 1D and 2D cases and the quenches are fast.Therefore, the qualitative behavior in this case is almost the same as the dynamical process from the unstable initial state in Figure 7. Similar to the results presented in the previous section, the evolution can be also divided into four stages.In summary, the spinodal decomposition induced by a fast quench is qualitatively the same as that induced by considering an initial unstable state.
From the QNM spectrum of the static uniform states we can conclude that the final uniform normal solution is still stable under homogeneous and inhomogeneous perturbations.This solution also corresponds to the lowest value of free energy at this value of charge density.Nevertheless, we still cannot exclude the existence of a non-uniform meta-stable final state with lower value of the free energy.Therefore, in order to confirm the nature of final state, the full time evolution is necessary.
To compare with the previous case, we now consider a different system with λ = −0.8 and τ = 0.28, in which we are able to confirm the non-uniform final state without running the full dynamical evolution.In this case, the uniform solutions with the same value of charge density are all unstable under linear perturbations.In the rest of this section, we present the condensate curves for this case and introduce the logic that leads to the confirmation a non-uniform final state.The previous analysis for the critical length is also valid here.Nevertheless, since the results are analogous, to avoid clutter, we only give results for the 1-dimensional and 2-dimensional quench "experiments" and focus on the non-uniform final state.
We plot the condensate curves in the canonical and grand canonical ensembles in the left panel of Figure 10, where we can see that the condensate curve in the canonical ensemble show a typical 2nd order phase transition, while the condensate curve in the grand canonical ensemble show a typical 1st order phase transition [51].From the analysis of the QNMs performed in previous studies [40], we know that the solutions below the turning point in the grand canonical ensemble are stable under homogeneous perturbations, but unstable under inhomogeneous perturbations.Therefore, if we choose the final value of the averaged charge density equal to that of a solution in this region, and quench the system from a stable normal phase, we will certainly get a non-uniform final state exhibiting phase separation.
The main reason for this is that in the canonical ensemble, at this value of average charge density ρ, the normal phase is unstable, and the superfluid solution has lower value of the free energy.However, this superfluid solution is unstable when inhomogeneous perturbations are included and the total charge kept fixed.As a result, this homogeneous state with lowest value of free energy at this value of average charge density ρ is also unstable.In summary, the final stable state must be some new non-uniform configuration.
To confirm this simple argument, we run a 1-dimensional quench "experiment" and plot the average value of condensate versus time t in the central panel of Figure 10.At the same time, we provide a 3D figure to show the evolution of the local value of condensate in the right panel of Figure 10.As expected, we observe that the final state is not homogeneous as it presents a large one-dimensional bubble.
Finally, we run a 2-dimensional version of the same quench process and show the evolution of local condensate in four snapshots at different time in Figure 11.It is clear that the system runs into a non-uniform final state with a big bubble as well.Notice that the final bubble is a bubble of normal state inside a superfluid region, similar to a superfluid vortex, and different from the previous cases, e.g., the last panel in Figure 9.
The spatial distribution of the grand potential provides valuable insights, enabling us to conduct a more concrete analysis of this dynamical behavior.First, we show the local value of the grand potential density for the 1D and 2D quench experiments in Figure 12.The detailed formula for the grand potential density in in-going Eddington coordinate is presented in the Appendix, Eq.(A6).In Figure 12, the first and second panels depict the grand potential density configurations of the 1-dimensional and 2-dimensional final states for the case λ = −2 and τ = 0.8, respectively.The third and fourth panels show the relative results for the case λ = −0.8 and τ = 0.28.These configurations clearly show the grand potential density barrier along the bubble walls.The height of the wall is larger in the first case and much lower in the second case, which is important in explaining the difference in the phase separation dynamics.In the first case, small bubbles emerge quickly in the second stage of time evolution.While in the second case, small bubbles do not show up clearly, and some connected narrow region emerge instead.The reason behind this difference is that in the first case the energy barrier of the bubble wall is large, and the decrease of the length of bubble wall is efficient in minimizing the total grand potential.However, in the second case, the height of the energy barrier on the bubble wall is much lower, and the width of the bubble wall is larger.In this case the phase separation process is not likely to produce many small bubbles.The detailed relation between the height of the bubble wall and the creation and evolution of small bubble is worth to be studied more carefully in the future.
Another important issue is that during the 1-dimensional evolution, the chemical potential inside and outside the bubble has to be balanced.In the 2-dimensional case, because of the different dependence of the surface and volume terms on the radius of the bubble, the grand potential density inside a stable bubble should be lower than the grand potential density outside [15].This feature is also evident in the four plots in Figure 12.
For convenience, in the numerical simulations, we fix the gauge by setting A t | z=0 to a constant independent of t and ⃗ x, instead of using the usual gauge choice A t | z=1 = 0 as in the static case.For any final equilibrium state, including the non-uniform states, the chemical potential, defined as should be balanced everywhere -the so-called "chemical balance condition" -as proved in [15].However, during the non-equilibrium evolution, the chemical potential varies in space.Therefore, the variance of the chemical potential defined as is a good measure of the degree of inhomogeneity during the time evolution -dynamical heterogeneity.We plot the spatial variance of chemical potential δµ for the 1-dimensional and 2-dimensional quench experiments in Figure 13.The left two panels are for the 1-dimensional (top panel) and 2-dimensional (bottom panel) quench experiments for the first case with (λ = −2, τ = 0.8).The two panels on the right side are for the 1-dimensional (top panel) and 2-dimensional (bottom panel) quench experiments for the second case with (λ = −0.8,τ = 0.28).In the two panels for the 2-dimensional processes, we use four dashed vertical red lines to indicate the times corresponding to the four snapshots of local condensate presented above in Figure 9 and Figure 11.We see that the variance of chemical potential δµ reaches an initial maximum, which roughly corresponds to the exponential growth of the unstable modes.
After the initial peak, there are some following smaller peaks, which corresponds to the regime in which small bubbles merge and disappear.The sub peaks are sharper in the 1-dimensional cases than in the 2-dimensional cases, because in the 1-dimensional case, the domain wall region expands into a larger portion of the whole volume.Finally, the system gradually approaches the stable final state forming a big bubble, and the variance of chemical potential δµ decreases to 0 at that point.

IV. OUTLOOK
In this work, we have performed several numerical experiments on the full dynamical evolution triggered by the spinodal instability in first order phase transitions in a holographic superfluid model.
In the first part of this work, we focused on confirming the linear dynamical stability analysis (based on the QNM spectrum) in the spinodal region.This critical value of wave vector k c identifies a critical value for the size of the system, l c = 2π/k c , below which the system is stable.This prediction has also been confirmed explicitly by our numerical "experiments" involving the fully nonlinear time evolution of the system.
In the second part, we presented the results of the full time dependent evolution in order to analyze in more detail the out-of-equilibrium dynamics, as well as the nature of the final state.We have considered both 1-dimensional and 2-dimensional quench "experiments" for two different configurations of the potential in the model.In the real time evolution from both the unstable initial state and the quench processes, the inhomogeneous instabilities trigger the spinodal decomposition.The following time evolution can be divided into four typical stages as suggested in Ref. [33].
The local value of grand potential for the final stable state is consistent with the general picture of bubbles.We also find that the variance of chemical potential δµ is a good indicator to describe the dynamical heterogeneity of the system during the real time evolution.There, a clear peak in the variance is observed at early times and is related to the sudden growth of the inhomogeneous instabilities.Following sub peaks are accompanied by the disappearance and merge of small bubbles.
Our study provides a concrete example of the dynamical evolution of a system under spinodal decomposition, which is triggered by the inhomogeneous linear instability in a region with negative susceptibility ∂ρ/∂µ < 0, when the size of the system is larger than some critical length.The initial and final regimes of the time evolution perfectly agree with the results of the linear dynamical stability based on quasi-normal modes [40].Our study also shows that the two thermodynamic ensembles (canonical and grand-canonical) are both necessary in order to understand the dynamical processes behind the spinodal decomposition.Finally, we also find that quenching the extensive variable ρ of the system is a perhaps more realistic way to reach the unstable spinodal region and study the full dynamics of the phase separation triggered by the spinodal decomposition in such a setup. .The top two panels denote the evolution of the variance of µ in one-dimensional quenching processes while the bottom two panels denote the relation in two-dimensional quenching processes.The two left panels are for the case with (λ = −2, τ = 0.8) and the right two panels are for the case with (λ = −0.8,τ = 0.28).We find that the peaks of the curves correspond to the generation or disappearance of domain walls and bubbles.In the bottom plots, we use four dashed vertical red lines to denote the times at which we plotted the snapshots of local condensate in Figure 9 and Figure 11.
There are many open questions left for the future.In the quenching experiment from a normal state to the final superfluid phase, the coexistence of two different topological defects, the bubbles and the superfluid vortexes, is expected in a large system.The detailed relation and effective interaction between the two is not known and can be studied using holographic methods in the future.Also, a more detailed analysis of the evolution of the energy power spectrum during phase separation might disclose important information about gravitational wave production from cosmological first order phase transitions.Finally, it would also be interesting to study the phase separation with a conserved total magnetic flux, and investigate its influence on the bubbles and vortexes.We are planning to pursue some of these directions in the near future.
setup and linear stability A. Holographic setup and the 1st order phase transition B. Quasinormal modes (QNMs) and linear stability III.Numerical experiments and nonlinear dynamical evolution A. Holographic setup for the dynamical tests B. Testing linear stability with dynamical processes 1. Homogeneous perturbations I. INTRODUCTION

4 FIG. 3 .
FIG. 3. Left: The condensate as a function of the charge density ρ (dashed line) or chemical potential µ (solid line).The two lines correspond respectively to the canonical ensemble and the grand canonical ensemble.The black triangle represents the critical point.The two dashed green horizontal lines denote the turning points in the two different ensembles.Center: The dependence of Im(ω) on k for the most unstable mode around different solutions.The dashed black line is the horizontal axes with Im(ω) = 0.The two dashed red vertical lines indicate the two spatial values of the wave-vector which are related to the two spacial scales lx=10 and 4, respectively, used in the dynamical tests.The two dashed blue vertical lines indicate k1 = 2π/100 and k5 = 5 * 2π/100 for the system with lx = 100, respectively.Right: The time dependent condensate of the uniform initial solution in the real-time dynamical homogeneous evolution.The uniform initial solution is marked with the same color in the left panel.The inset is a cartoon to illustrate the two different evolutions for the magenta and cyan curves.

FIG. 4 .
FIG. 4. The time dependent value of the average condensate from different initial solutions when the inhomogeneous states are involved in the time evolution.The size of the system along the x direction is set to lx = 4 for the left panel and lx = 10 for the right panel.The lines with different colors represent the time evolution from different initial solutions.The colors are the same as in the left plot of Figure 3.The condensate values are the average in the x direction.The insets show the configuration of the local condensate in the final state for the cyan solution.

FIG. 7 .
FIG. 7. Local value of condensate at different time slices during the real-time evolution from a homogeneous initial state in one (top panel) and two (bottom panel) spatial dimensions with system size lx = ly = 100.The initial state corresponds to the violet point in Figure 3 with charge density ρ = 3.00.

5 FIG. 8 .
FIG.8.Time dependent quantities for a real time process with the average charge density quenched from ρi = 3.5005 in the initial stable state to ρ f = 3.00, which is a value for the superfluid solution located in the spinodal region with ⟨O⟩ ρ t < ⟨O⟩ < ⟨O⟩ µ t .The spatial size of the boundary system is lx = 100.Left: The solid violet curve is the average value of condensate as a function of time t while the solid black curve is the average charge density ρ as a function of t, with dashed black vertical line denoting the end time of quenching the charge density.Right: Local condensate as a function of time t and space x.

FIG. 9 .
FIG.9.Local condensate at four different time slices of a real time process with the average charge density quenched from ρi = 3.5005 of the initial stable state to ρ f = 3.00 which is a value for the superfluid solution located in the spinodal region with ⟨O⟩ ρ t < ⟨O⟩ < ⟨O⟩ µ t .In this quench process the two parameters for nonlinear terms are set to λ = −2, τ = 0.8.

8 FIG. 10 .
FIG. 10.The condensate curves for uniform static superfluid solutions for λ = −0.8 and τ = 0.28 in two different ensembles and time dependent condensate in a one dimensional quench process from the initial stable normal state with ρi = 3.50 to ρ f = 4.71 which is a value for the superfluid solution in the spinodal region.The size of the system is set to lx=100.Left: The condensate as a function of the charge density ρ in the canonical ensemble denoted by the solid blue line and the condensate as a function of chemical potential µ in the grand canonical ensemble denoted by the dashed red line.The brown point is the end point of the quench process with ρ = 4.71.Center: The solid brown curve is the averaged value of the condensate as a function of time t.The solid black curve is the average charge density ρ as a function of t and the dashed vertical black line indicate the end time of quenching ρ.Right: Local condensate as a function of time t and space x.

FIG. 13 .
FIG.13.Dynamical heterogeneity during the out-of-equilibrium evolution.The averaged value of the variance of the local value of chemical potential (µ − μ)2 .The top two panels denote the evolution of the variance of µ in one-dimensional quenching processes while the bottom two panels denote the relation in two-dimensional quenching processes.The two left panels are for the case with (λ = −2, τ = 0.8) and the right two panels are for the case with (λ = −0.8,τ = 0.28).We find that the peaks of the curves correspond to the generation or disappearance of domain walls and bubbles.In the bottom plots, we use four dashed vertical red lines to denote the times at which we plotted the snapshots of local condensate in Figure9and Figure11.

TABLE I .
The different initial states considered in the "numerical experiments" and the turning points.Specific values of ρ, µ and ⟨O⟩2/ρ for the different colored points in Figures2 and Figures 3. The black triangle is used to mark the "critical point" in both the two figures.The red and green crosses can be found in Figure The critical instability scale.Left: Im(ω) as a function of k for the most unstable mode.The initial states of these colored lines are superfluid solutions in the range ⟨O⟩ ρ t < ⟨O⟩ < ⟨O⟩ µ t .Right: Numerical test of the critical scales.Magenta triangles represent the critical scale extracted from the dispersion relations of the unstable modes presented in the left panel.The blue asterisk shows the length scale at and above which the initial solution runs away before tmax, while the blue cross show the length scale at and below which the initial solution does not evolve in time up to tmax.