Holographic boiling and generalized thermodynamic description beyond local equilibrium

Tuning a very simple two-component holographic superfluid model, we can have a first order phase transition between two superfluid phases in probe limit. Inspired by the potential landscape discussion, an intuitive physical picture for systems with first order phase transitions is provided. We stress that holography perfectly offers a generalized thermodynamic description of certain strongly coupled systems even out of local equilibrium, which enables us to carefully study domain wall structures of the system under first order phase transitions, either static or in real time dynamics. We numerically construct the 1D domain wall configuration and compute the surface tension of the domain wall from its generalized grand potential. We also numerically simulate the real time dynamics of a 2D bubble nucleation process (holographic boiling). The surface tension of the 1D domain wall nicely matches the final state of the 2D bubble nucleation process when the bubble radius is large enough. This study can be regarded as a fist step towards a holographic description of the Helium-3 dynamics, which is a strongly coupled system and has two superfluid phases with a first order phase transition.

First order phase transitions are not only quite common in our daily life, but are also important phenomena in superfluids such as Helium-3. For example, the phase transition between the two different kinds of superfluid phases, namely the A-phase and B-phase of Helium-3, is a first order one [17]. An important universal feature accompanied by the first order phase transition is the creation and evolution of bubbles and domain walls.
Such structures also play an important role in the evolution of the early universe and in gravitational wave physics [18][19][20].
However, it is quite difficult to accurately describe the first order phase transitions, because dynamical process of first order phase transitions are usually far from (local) equilibrium. Such a dynamical process in quantum many body systems is extremely complicated in general, so certain theoretical description of a general dynamical process is very important and helpful. Such kind of efforts for holographic systems are made recently in [21], where some generalized thermodynamic (hydrodynamic) description beyond local equilibrium turns out to be available and a generalized grand potential (free energy) can be well defined accordingly. We then would like to see some concrete application of this description in, e.g. superfluid systems with first order phase transitions mentioned above, where it can be verified and hopefully even developed further.
In a recent paper [22], the authors realized a dynamic process of domain wall formation holographically. There are also some other progresses in this topic [23,24]. However, in these studies the holographic systems are effectively 1D, so more realistic structures like bubbles (at least 2D) cannot be considered. In order to study the more realistic and more complicated structures in first order phase transition, it is wise to consider a holographic model with first order phase transition in probe limit, where back-reaction of the matter field to the bulk space-time is neglected and complex numerical relativity problems can be avoided. Such a system with first order phase transition in the probe limit can be easily built from a holographic model with two competing orders, as in Ref. [25].
In this paper, we study the nontrivial configurations and dynamical processes in the first order phase transition between superfluid phases holographically. In fact, our model here is probably the simplest for holographic study of first order phase transitions, and it turns out that essential physics can already be observed quite well in our system. We first realize the first order phase transition between two s-wave superfluid phases holographically, and give a concrete picture to explain how the first order phase transition occurs from a "grand potential landscape" point of view.
This grand potential landscape point of view is supposed to be true not only in the holographic context, but also in general thermodynamics, offering an intuitive picture of the stability of phases and the phase structures. Actually, in thermodynamic systems with first order phase transitions (like the van der Waals gas-liquid system), as well as in many holographic models including ours here, there are typically three equilibrium states involved, the stable one, the meta-stable one and the unstable one. It can be argued that the unstable state is a saddle point in the landscape and acts as the minimal potential barrier between the stable and meta-stable states, which will also be verified in our discussions of inhomogeneous configurations (mixture states of different phases).
For the inhomogeneous configurations with local structures like domain wall, which go beyond local equilibrium, we can use the grand potential of the unstable state to estimate the maximum of the grand potential density of the domain wall (bubble) structures, from the grand potential landscape point of view. A more quantitative discussion of such configurations out of local equilibrium needs to be based on the generalized thermodynamic description introduced in [21], where it turns out that the potential (including the generalized free energy and grand potential) landscape point of view is also crucial. To show the practicability of this point of view and the generalized thermodynamic description more concretely, we numerically build some typical inhomogeneous configurations. In the numerical work, we first construct the domain wall between the two superfluid phases in a simple 1D setup of our model. We can calculate the surface tension of the domain wall from its generalized grand potential. We then study the dynamical process of bubble nucleation after a local quench to a homogeneous meta-stable state in 2D and find that the surface tension with "generalized balance conditions" for the bubble configuration matches well with the above calculation for the 1D domain wall configuration. This paper is organized as follows. We firstly realize a first order phase transition between two s-wave superfluid phases holographically, and give the concrete picture of first order phase transitions from a "grand potential landscape" point of view in Sec. II. We then give a domain wall configuration in Sec. III, and study the dynamical process of bubble nucleation in Sec. IV. We finally take some conclusions and discussions in Sec. V.

II. HOLOGRAPHIC SUPERFLUID MODEL WITH A FIRST ORDER PHASE TRANSITION
We use the holographic s+s model to realize a first order phase transition in the probe limit, under which the gravitational background is fixed, and the numerical work can be greatly simplified. The gravitational background can be taken as the 4 dimensional asymptotic AdS black brane with The boundary superfluid system is in a thermal bath with the same temperature as the Hawking temperature of the bulk black brane.
The Lagrangian for the matter fields is The probe limit can be attained consistently by taking the large charge limit, and a scaling symmetry implies that only the ratio e 1 /e 2 is important. Therefore we set e 2 = 1 in the rest of this paper.
With the standard procedure, we can get the solutions dual to superfluids on the boundary field theory. Because we have two orders dual to Ψ 1 and Ψ 2 respectively, we can get three Solution S1+S2 respectively. The dashed lines denote the condensate value in unstable or meta-stable sections. and set Ψ 1− = 0 for Ψ 1 and Ψ 2+ = 0 for Ψ 2 . In this way, we have chosen different quantizations for the two scalar fields, the standard quantization for Ψ 1 and the alternative quantization for Ψ 2 . We define the boundary operator dual to Ψ 1 as O 1 and the boundary operator dual to Ψ 2 as O 2 . In this way, we get a dimension 2 operator O 1 and a dimension For the two scalar operators with different conformal dimension, the condensed solution will get "parallel" thermodynamic potential curves if the two orders also have the same charge. We can therefore tune the ratio e 1 /e 2 to make the two thermodynamic potential curves have an intersection point. Near this intersection point, an s+s coexistent solution appears. For convenience, we can choose this ratio to be e 1 /e 2 = e = 4.5.
Without considering the interaction term (or set λ 12 = 0 equivalently), the s+s solution will emerge near the intersection point of the two thermodynamic potential curves, and connect the two solutions with single condensate. It is the one with lowest thermodynamic potential in the region where it exists. In this case, the system undergoes two second order phase transitions from the stable region of Solution S1 to the stable region of Solution S2, showing a typical "x-type" figure for condensates as in the left panel of Fig. 1.
In order to get a first order phase transition, we can tune the value of λ 12 . In Ref. [25], the influence of a such kind of interacting parameter in a holographic system with multiple orders is already shown in the s+d system. Here in the holographic s+s model, the qualitative S1 phase S2 phase S1 S2 phase  Figure 2: λ 12 − µ phase diagram. The right plot is the enlarged version of the dashed rectangular region in the left one.
influence of λ 12 on the phase diagram is the same. It is easy to see that this term will not change the condensate as well as the thermodynamic potential of the superfluid solutions with single condensate. It will only change Solution S1+S2 and has a quite obvious influence.
Generally, if we decrease the value of λ 12 , Solution S1+S2 will exist in a larger region, and become stabler, while if we increase the value of λ 12 , Solution S1+S2 will firstly shrink to be in a smaller region, and become less stable. When we keep increasing λ 12 , at a certain point, Solution S1+S2 will become totally unstable. After that, the region of Solution S1+S2 will then increase, but its thermodynamic potential is still increasing, making this coexisting solution more and more unstable. With some numerical work, we have confirmed that when λ 12 > λ c ≈ 0.0962 , the Solution S1+S2 becomes totally unstable, and there is a first order phase transition between S1 and S2.
To make the influence of λ 12 more concrete, we draw a λ 12 − µ phase diagram in Figure 2.
We choose µ rather than ρ as the horizontal axis because it is more convenient to work in the grand canonical ensemble. In this phase diagram, the red and blue lines for second order phase transitions are made up of critical points of second order phase transitions. It is a bit more complex to get the green and black lines for first order phase transitions, because we need to compare the thermodynamic potential of two solutions to get the first order phase transition point [26].
In the phase diagram, we can see that there is a triple critical point among the S1 phase, S2 phase and the S1+S2 phase. Above the triple critical point, there is a first order phase transition between the S1 phase and S2 phase. Below the triple point, an S1+S2 phase exist between S1 phase and S2 phase. Further numerical calculation tells that in this system, the phase transition between S1+S2 phase and S2 phase is always second order, while the phase transition between S1+S2 phase and S1 phase is second order at a lower value of λ 12 (related to the solid red line) and is first order at a higher value of λ 12 (related to the solid green line).
To study the detail of a first order phase transition, we can also choose λ 12 to be in the region related to the green line. However, we choose λ 12 > λ c in order to make the transition clearer. The stable phase on one side is Solution S1 with only order O 1 nonzero, and the stable phase on the other side is Solution S2 with only order O 2 nonzero. In the next subsection, we give the condensate behavior and grand potential of the three homogeneous superfluid solutions, which is helpful to understand the picture of first order phase transition.

B. Homogeneous solutions
When λ 12 > λ c , the s+s solution becomes totally unstable, and we can get first order phase transition between S1 phase and S2 phase. In order to get the numerical solution with a thin bubble wall, we wish to get a larger value of potential barrier, this is related to a larger value of λ 12 . Thus in the next section, we set e 1 = 4.5, e 2 = 1, to get the inhomogeneous numerical solution dual to a mixed state on the boundary field theory.
In the rest of this section, we show the condensate behavior and grand potential of the static homogeneous solutions with the above choice of λ 12 . The condensate behavior is plotted in the right panel of Figure 1, and the grand potential curves are plotted in Figure 3.
In the AdS/CFT dictionary, the grand potential of the three solutions can be calculated from the Euclidean on-shell action on the bulk side. Because we work in the probe limit, the contribution from the gravitational part is the same for different solutions, and we only need to calculate the contribution from the matter part (4) to the grand potential [27]: where S ME denotes the Euclidean on-shell action of matter fields in the black brane back-  Figure 3: Grand potential of Solution S2 (blue) and Solution S1+S2 (green) with respect to that of Solution S1 (red). The orange line is for comparison, which denotes the grand potential of Solution S1+S2 with λ 12 = 0.
ground. With the ansatz and equations of motion, we can finally get with V the total volume of the space. Since we only consider the probe limit in this paper, we will suppress the subscript m of the grand potential (and other thermodynamic quantities) in the rest of this paper.

C. A picture for first order phase transitions
In this section, we describe a concrete picture for first order phase transitions, which is consistent with the swallow tail shape of thermodynamic potential curves as well as our numerical results for the inhomogeneous domain wall and bubble configurations in the following sections.
The general picture of first order phase transitions tell us that when a first order phase transition occurs, the system changes from a meta-stable state (which is a local minimum in the configuration space) to a stable state (which is a global minimum in the configuration space). We must stress that in order to better understand the first order phase transition, we need also involve the third state, which is an unstable state and is a saddle point in the configuration space.
Another feature of the first order phase transition is that the thermodynamic potential usually form a "swallow" tail shape as in Figure 3, where we can see that the "swallow tail" is formed by three part: the red curve, the blue curve and the green curve. The red and blue curves are related to the stable and meta-stable states, with the lower one more stable, while the green curve is related to the unstable state.
In a homogeneous solution of a thermodynamic system, we have the (static) equation In order to clarify the physical picture, we can consider the full configuration space, which includes not only the states that satisfies the equation of state (equations of motion in the bulk), but also the non-equilibrium or more general states that do not satisfy the equation of state. In this configuration space the stable solution, the meta-stable solution and the unstable solution can all be connected continuously through some general states. If we draw the thermodynamic potential in this space, we can get a landscape of the potential (see Sec. II E for details), where the stable and meta-stable solutions are the local minimum points, and the unstable solution is a col (saddle point).
We draw the thermodynamic potential curve with swallow tail shape in Figure 4, where we mark five typical points and for each point we show a special curve in the grand potential landscape as in Figure 5. For Point 1 , the S1+S2 solution just emerge from the S2 solution.
We can see that the blue and green points coincide at a stationary point. At Point 2 , the S1+S2 solution, which is a col point, is separated away from the S2 solution, which is a local minimum point. Point 3 denote the transition point for the first order phase transition, where the thermodynamic potential for the two minimum are equal. Point 4 is similar to Point 2 , where the S1+S2 solution is also a col point. The difference is that in this situation, the S2 solution has a lower value of grand potential. Point 5 is similar to Point 1 , which  are both at the tip of the swallow tail, where the S1+S2 merge with the S1 solution.
With the above picture, we can see that to make a first order phase transition occur, the meta-stable state need to skip over a barrier and go to the stable state. The unstable state is right at the point through which the thermodynamic potential barrier is lowest. The reason is the following. Imagine that we draw all lines connecting the stable state and the unstable that there is only one saddle point in the landscape, 1 we see that this saddle point is just the unstable state, which is by design the lowest barrier between the stable and unstable states. Therefore the unstable state is very useful to estimate the (nonlinear) stability of meta-stable state qualitatively.

D. Domain walls and bubbles
In first order phase transitions, a characteristic phenomenon is phase separation, where there is an interface between different phases, the domain wall. In higher (larger than one) dimensions, the domain wall can have various topologies and shapes, typically a bubble. The two homogeneous regions separated by the wall are in different stable (or meta-stable) states, and the wall region is where the system change from one state to another. Therefore, the wall region is supposed to have higher thermodynamic potential density than the homogeneous region, and one should take into account this excess amount of thermodynamic potential, which is well known as the surface tension and will be further discussed in the following sections, to understand the mechanics of bubbles.
With the above picture of first order phase transition, we can see that through the wall region, the local state must go up the thermodynamic potential barrier. Because the lowest barrier is right on the unstable state, the point with maximum thermodynamic potential in the wall would just be very similar to the unstable state on the swallow tail, which means that they are very close in the configuration space mentioned in the last subsection. The route in configuration space related from the points on the bubble wall may not pass the col point because of the complex shape of the landscape. Another difference is that the points on the wall involve additional spacial gradient.
We shall obtain configurations with bubble wall in two different ways. First we build a domain wall solution with only one inhomogeneous spatial dimension. And later, we obtain two dimensional bubble configurations using a quench method.
To make our previous arguments for the first order phase transition clearer and quantitative, we shall discuss here the landscape picture of the free energy (grand potential) of holographic systems, taking the simplest holographic superfluid model as an example.
Although numerical evolution of holographic dynamics is usually (and most simply) done under the Eddington-Finkelstein coordinates (68), in order to make our discussion as general as possible, we will consider a class of coordinates (including the Eddington-Finkelstein one as a special case), under which the time evolution of holographic systems can be done in principle. In this class of coordinates, ξ = ∂ t is a time-like Killing vector field. The futuredirected orthonormal co-vector to a constant t surface is and so the total flux of the energy current across a constant t surface Σ (null for the Eddington-Finkelstein case or space-like for more general cases) is Note that in holographic systems at finite temperature, the surface Σ is bounded by the horizon and the asymptotic AdS conformal boundary.
Let us first illustrate the picture by a scalar field, the simplest case. The Lagrangian density is so we have which is obviously a sum of the kinetic energy density and the potential density. Note that here we achieve this nice expression under a coordinate system that is generally not time orthogonal. Actually, the potential density and we can define the static free energy (grand potential) as On the other hand, the generalized free energy (grand potential) is just the total flux (16): We see that due to the positivity of the kinetic term, where the equality holds when the configuration is static or g tt = 0 (i.e. the Eddington-Finkelstein case). Locally, in holographic systems we can define as the static and generalized free energy (grand potential) densities, respectively.
The full equation of motion of the scalar field can be written as while the static equation of motion, i.e. the equation of motion for static configurations, is just the extreme of the variation The time evolution of the generalized free energy (grand potential) can be computed as Under the ingoing Eddington-Finkelstein coordinates or more general ingoing coordinates the only non-vanishing g zµ | z h is g zt | z h < 0, so the generalized free energy (grand potential) decreases monotonically if there is no source on the AdS conformal boundary, which can be proved in more general holographic systems with the null energy condition [21]. Physically, this decrease (the energy flux across the horizon) is interpreted as the energy dissipation in this dynamical process [13,28,29].
For a Maxwell field, the situation is subtler. The Lagrangian density is so we have On the other hand, the flux of the Noether current [21] (or called the canonical current) is given by the Hamiltonian density. The canonical momentum density is so the Hamiltonian density is Then it is easy to check that which is consistent with the general fact that these two currents are equal on shell (the last step) up to a total divergence [21]. Holographically, this total divergence gives a boundary contribution µρ to the flux across Σ, which in the equilibrium case is just the difference between the free energy and the grand potential. Specifically, the flux of the energy current is the generalized free energy F and that of the Noether current is the generalized grand potential Ω: where ∼ = means equality on shell.
In the above case, the potential density is i.e. the part of L independent of any time derivative. Then the kinetic energy density is which is positive definite. The static grand potential Ω s is still defined by (20), as well as the corresponding density which again satisfy and (22), respectively, due to the positivity of the kinetic energy density (36). A caveat here is that both the potential density (35) and the kinetic energy density (36) are gauge dependent, so certain gauge fixing will be helpful. In particular, under the Eddington-Finkelstein coordinates, since g tt = 0 and g ti = 0 unless i = z, if we take the radial gauge A z = 0 as usual, it is obvious that the kinetic energy density (36) completely vanishes even for non-static configurations, similar to the scalar field case.
For the coupled Maxwell-scalar case similarly we have The potential density is and the Hamiltonian density Then the kinetic energy density is again positive definite. Under the Eddington-Finkelstein coordinates and the radial gauge, this kinetic energy density again vanishes even for non-static configurations. For our double scalar field model, the generalization of the above discussion is straightforward.
However, in holography, the above quantities generally diverge due to the asymptotic behavior approaching the AdS boundary, which need to be renormalized. In the standard procedure of holographic renormalization (see, e.g. [30]), there will be a boundary counter added to the original action in order for the on-shell action and other quantities to be finite when the boundary tends to the AdS conformal boundary. This counter term leads to a "boundary energy momentum tensor" which gives a boundary contribution to the grand potential density ω. The static part of (45) gives a boundary contribution to the static grand potential density ω s , too.
Anyway, the landscape picture for holographic systems is just like that in the scalar field case. However, in the case involving a gauge field, it turns out that only the holographic boundary condition of fixed total particle number is possible for the dynamic evolution [21], i.e the holographic systems can be sourceless for the free energy F instead of the grand potential Ω. Accordingly, only the free energy has the property of monotonic decreasing during a dynamical process without driving. The full equations of motion have a form similar to (24), where all the time derivative terms are on the left hand side and the right hand side is a functional gradient of Ω s . But the potential landscape should be considered as the landscape of the static free energy 2 with the static equations of motion Actually, the form of the above equation is the same as δΩ s = 0, while the only difference is the boundary condition (fixing the total particle number (63) versus fixing µ). So the dynamics is just a point particle rolling on the landscape with friction (dissipation). Eventually, it will tend to stop at some local minimum of the landscape, i.e. a static state satisfying (50) and δ 2 F s ≥ 0. In particular, under the Eddington-Finkelstein coordinates, since the kinetic energy always vanishes, the dynamics is extremely clear as the point particle keeping rolling down the landscape with all the decrease of its potential energy directly dissipated.
In the above discussion, we have considered the most general inhomogeneous configurations. To make things simpler, we can also consider the landscape of free energy (grand potential) for the homogeneous configurations, where all the spatial derivatives (other than z) are dropped from the quantities in the above discussion, and then we do not need to distinguish between F s (or Ω s ) and their densities f s = ω s + µρ (or ω s ). In particular, sometimes we would like to explore the spacial variation of the local state in an inhomogeneous static configuration while ignoring the spacial gradient, e.g. as in the intuitive discussion of the domain wall configuration in Sec. II D. In this case, since the chemical potential µ 2 Note that F s is not F (∂ t = 0), though they are equal for static, on shell configurations. Generally, from the relation (34) and the definition (49) of F s we have with T the total kinetic energy, which is the integral of (43) for the coupled Maxwell-scalar case.
is the same everywhere in this static system (see the discussion around (62) below), it is especially convenient to consider this simpler landscape of grand potential density ω s . So our discussion here offers a concrete description of the physical picture in the previous two subsections.
In either case, we can explore the landscape by the following two ways: • The global way: finding extremal points on the landscape by solving the static equations of motion For simplicity, we use the scalar field case to illustrate our argument, keeping in mind that it also works for more general cases. We start from an extremal point on the free energy landscape, which means δF s = 0 there, and consider the linear perturbation around that point. Suppose the dominant mode (the quasi-normal mode with largest imaginary part) is ω, so the perturbation δφ of the real scalar field takes the form From the previous discussion we know and where we have defined the second order static free energy as an inner product on the functional space of δφ. 4 (In the case involving gauge fields, the relations similar to (52) and (53) only hold on shell, as mentioned in Footnote 2, but that is not a problem because the linear perturbations here do be on shell.) So we have where we have used the fact that the inner product is symmetric. Thus, δ 2 F cannot be Conversely, if (, ) s is not positive definite (but is of full rank), we can choose an initial δφ such that (δφ, δφ) s < 0.
In this case, at least one of the quasi-normal frequencies ω should have a positive imaginary part, since otherwise the expansion clearly shows that δφ will tend to zero eventually, again contradicting the fact (47). In sum, there is or is not an unstable mode if (, ) s is indefinite (but of full rank) or positive definite, respectively. In between, (, ) s will be degenerate (actually positive semi-definite), where it has a zero mode and ω = 0 becomes a quasi-normal frequency. Therefore, in this mechanism, the system loses its stability always by a quasi-normal mode crossing the real axis at exactly zero.
Since the holographic superfluid model provides a description, valid at all scales, of certain strongly coupled superfluid system, it serves as a very powerful tool to study dynamic (far from equilibrium) processes of inhomogeneous configurations, even out of local equilibrium [13][14][15]. In general cases without local equilibrium, it is expected that most thermodynamic quantities cannot be well defined for the system, as well as most (even local) thermodynamic relations cannot hold [31]. But in the context of holographic investigations, there have been hints that some generalization of non-equilibrium thermodynamics (hydrodynamics) can be applied to such systems [21]. Systems with first order phase transitions have inhomogeneous (mixture) configurations of different phases, which should be the ideal arena to carefully investigate to what extent a generalized thermodynamic description works for such systems. Therefore, we shall try to make it clear in this section. First, we propose a general formalism as follows.
In principle, a complete formalism should involve backreaction of the matter fields onto the bulk geometry, but that will render the problem extremely complicated [21]. Here, for simplicity, we stay in the probe limit, tailoring the discussion to the holographic superfluid model in this paper. Roughly speaking, a holographic model in the probe limit corresponds to the limit that certain quantities, e.g. the specific heat capacity, of the boundary system tends to infinity, so that its temperature and some other parameters cannot be changed.
For a static configuration at finite temperature, we have with I bulk the (renormalized) bulk on-shell Euclidean action, β the inverse temperature and Ω s the static grand potential defined in (20), which is the grand potential of the boundary system by the holographic duality. Besides the inverse temperature β that is fixed to be homogeneous in the probe limit, the chemical potential µ and the grand potential are originally defined in global equilibrium, which can be generalized to the case of local equilibrium if the spatial variation of the configuration is gentle enough. In order to cope with local structures like domain walls, which obviously break local equilibrium due to the small scales characterizing them, now Ω s is taken as the (generalized) grand potential of the boundary system even for arbitrary static configurations, i.e. beyond local equilibrium, where the (generalized) chemical potential µ( x) = A t ( x)| z=0 can vary acutely in space. For the variation of the bulk on-shell action, we have with ρ( x) the local particle number density, so the grand potential Ω s satisfies as a generalization of the usual thermodynamic relation. Then the (generalized) free energy F s is given by the following functional Legendre transform (as the time-independent version of (49)): with the transformed on-shell action which yields the variation Under the boundary condition that there is no particle exchange between the system and the environment or there is no boundary, e.g. periodic boundary conditions, which means that the total particle number is a constant, a static configuration should satisfy (50). Combined with (62) and the condition (50) gives the result that µ is a constant (independent of x). We call it the chemical balance condition for static configurations without local equilibrium, which is a generalization of the corresponding condition in the ordinary thermodynamics and will be confirmed in our numerical evolution of configurations containing domain walls.
In the absence of local equilibrium, e.g. at the domain wall, the standard local thermodynamic relation (the Gibbs-Duhem relation) with ω the local grand potential density (i.e. the grand potential per unit volume) and p the pressure does not hold. However, away from the wall, the relation (64) should hold due to the restoration of local equilibrium. Now consider the surface tension of the domain wall. The surface tension coefficient σ is defined as the external work W done (under certain conditions) to enlarge the domain wall by a unit area. Besides the isothermal condition, the certain conditions for our static inhomogeneous configurations here also include the chemical balance condition described above. Under those conditions, the external work is just the increase of the grand potential due to the appearance of the domain wall: with Ω 0 the grand potential of the corresponding homogeneous system (without the domain wall). In the 1D domain wall case, the pressures on different sides of the wall should be equal to the same value p, which can be called the mechanical balance condition, otherwise a variation of the position of the wall would result in a change of the total grand potential and so the system cannot be in a static state. Suppose ω 1 and ω 2 are the grand potential densities of the two phases separated by (and far away from) the domain wall, respectively.
Then we have Actually, we see Ω 0 = −pV . At the position of the wall, ω will be significantly different from −p, which is the real contribution to (65). So σ is reduced to the following integral: with x 1 and x 2 the left and right edges of the domain wall, respectively. The above arguments will be verified in the following numerical evolution.

B. Numerical evolution of the mixture configuration and domain wall in 1D
In order to evolve the system nonlinearly, we switch to the Eddington-Finkelstein coordinates, under which the metric of the bulk black hole becomes We have adopted the radial gauge A z = 0 when working in the Schwarzschild coordinates, which is not satisfied after the coordinate transformation. To preserve the radial gauge, a U (1) gauge transformation is required: As a result, Ψ 1 and Ψ 2 should not be taken real anymore. Moreover, we should recover the spatial direction x, so the equations of motion are Here, we have made replacements Ψ i = zψ i (i = 1, 2). We will adopt the same nonlinear time evolution scheme as in (for example) [32], where (76) is taken to be a constraint equation.
We impose Dirichlet boundary conditions for ψ 1 and A µ at z = 0, while a Neumann-like boundary condition is adopted for ψ 2 [14]: This boundary condition can be easily derived from the coordinate transformation of D z ψ 2 | z=0 = 0 in the Schwarzschild coordinates to the Eddington-Finkelstein coordinates.
The two different boundary conditions for the two scalar fields correspond to the two different quantizations, respectively, taken in our model in Sec. II A. For the x direction, we impose periodic boundary conditions for all fields.
Formation of domain wall structures can be observed from nonlinear time evolution after perturbing the coexisting state by varying ρ. The form of the initial perturbation is our nonlinear time evolution, the coexisting state breaks down even when A is very small (≈ 10 −12 ), which confirms us with the fact that this coexisting state is unstable. After a long time evolution, the 1D domain wall structure forms finally. Fig. 6 shows a typical result when λ 12 = 0.4 and µ = 1.47. In this plot, distributions of order parameters O 1 , O 2 , particle number density ρ and grand potential density ω are given as functions of x, where the explicit form of grand potential density can be found in Appendix C. In all these figures, it is obvious to find sharp structures at the center of domain wall. By employing the numerical integral, σ defined in equation (67) can be calculated: Similarly, we can evolve the system nonlinearly with different values of λ 12 . Our results are shown in Fig. 7, where µ = 1.47, the same as the last paragraph. In Fig. 7a, the tension σ is calculated, in which an approximately linear relation can be witnessed between σ and λ 12 . In order to estimate the thickness of a domain wall, we calculate the relative height of the grand potential density ω, where the (λ 12 -independent) background value of ω far from the domain wall has been subtracted, as shown in Fig. 7b. To calculate the thickness of the domain wall, we define the region where relative heights are greater than 30% of their maximum to be wall areas. The results are shown in Fig. 7c, from which we find that the thickness of the wall decreases with λ 12 .

IV. BUBBLE NUCLEATION AND STABILIZATION IN 2D
In the previous section, we have investigated inhomogeneous holographic configurations from the generalized thermodynamic point of view, and then concretely calculated the surface tension of the domain wall from the final mixture state of a simple 1D evolution of our holographic superfluid model.
In this section, we will test in the 2D evolution the validity of our generalized thermodynamic picture and, in particular, the concrete value of the surface tension calculated in the 1D setup. In the 2D case, the domain wall is an effectively 1D object separating different phases, which can have certain shapes (and, in higher dimensions, also topologies) as bubbles. Then the surface tension just manifests itself as making the domain wall of least with r the radius of the bubble. The total volume V and particle number N are both fixed, so we have The equations of state read Here ρ 1 (µ) and ρ 2 (µ) are given functions characterizing the properties of the two phases, respectively, which can be obtained numerically for our model, as shown in Appendix B.
The chemical potential µ is the same for the phases in and outside the bubble, due to the chemical balance condition proved in the last section. As usual in textbooks, we also have the mechanical balance condition for the bubble: Here p 1 (µ) and p 2 (µ) are also equations of state in our model, which are actually plotted in unknown variables (r, V 1 , V 2 , µ, ρ 1 , ρ 2 ), which can be solved straightforwardly.

B. Numerical evolution of bubble nucleation and stabilization
We will quench a meta-stable state locally. From numerical evolution, we find that proper quenches can trigger a (local) first order phase transition, eventually resulting in a stable bubble configuration under certain conditions as described in the previous subsection.
We will consider 2-dimensional configurations, so both the x and y directions must be included into the equations of motion. Boundary conditions for the z direction are the same as the domain wall formation in the 1D case. (A y has same boundary condition as A x .) For both the x and y directions, we impose the periodic boundary conditions. Without loss of generality, the meta-stable state that we are going to quench is the state with a single ψ 2 condensation. The amplitude of this local quench should be large enough to overcome the potential barrier, as described by the physical picture in Sec. II.
The time evolution can be divided into two stages. The first stage t < 1000 corresponds to the bubble formation process. To characterize this formation process more precisely, we introduce a locally defined energy dissipation [13]: where T M N is energy momentum tensor. In (holographic) superfluids at finite temperature, dissipations occur around moving local structures [13], so we expect that this dissipation reaches its maximum at the domain wall position. As a result, the moment when this dissipation distribution becomes round can be defined as the bubble formation time, since we begin with a non-round quench here. Moreover, we can see that all the fields including the order parameters and particle density have round distributions. The second stage t > 1000 corresponds to expansion of the bubble. In order to describe the speed of expansion, distance between the position of maximum energy dissipation and the origin can be adopted to characterize the size of the bubble. Fig. 8a shows the evolution of radius r calculated from the dissipation. In this plot, the bubble expands all the time when t < 6000 with its speed slowing down gradually. Then when t > 6000, the radius r tends to a constant. As the system approaches a stable state, the dissipation becomes rather small (< 3 × 10 −9 at the final stage of our evolution), which in practice makes it not a proper method to determine the bubble size by the maximum of dissipation. However, as shown in Fig. 8a, it is possible to obtain a rough size of the bubble from the dissipation even when t < 1000 (in contrast  to the grand potential density discussed below), so this method enables us to know more about the bubble radius at the initial stage.
Aside from the dissipation, the evolution of the radius r can be determined by studying the dynamical grand potential density (see Appendix C for the concrete form of the grand potential). In this way, the radius r is calculated by the distance between positions where the grand potential density reaches its maximum and the center of the system. Fig. 8b shows the result, where we can find a similar shape to that calculated by dissipation in Fig. 8a when t > 1000. The curve in Fig. 8b extends smoothly to the final state when the system stabilizes, so it is a better choice to describe the evolution of radius when t is large. And we are able to read the radius of the bubble from Fig. 8b at the final stage: Moreover, configurations of the final stable bubble are plotted in Fig. 9, where, similar to  potential can be defined locally as the gauge invariant expression As has been discussed in Sec. III A, the chemical potential must satisfy the chemical balance condition, if the system has stabilized. Numerically speaking, the static state is achieved if µ approaches a constant independent of t and x. In our time evolution, this condition is specified as that the standard deviation of µ is small enough (< 10 −6 ). Moreover, we can read off the chemical potential for the final static bubble configuration: µ b = 1.4292.
Substituting µ b and σ calculated from (67) into (81-86), we can obtain another radius of the bubble: which is close to the result in (88). There are mainly two factors for the deviation between r b and r. The first one is that we have taken the surface tension σ at the critical chemical potential, which is not the true tension of the bubble here. The second one is the curvature effect of the domain wall, which can be estimated to be of the order which is consistent with the order of our deviation here.

V. CONCLUSION AND DISCUSSION
In this paper, we have studied the first order phase transition between two superfluid In order to perform the analysis of quasi-normal modes numerically, we start from the equations of motion in Eddington-Finkelstein coordinates, under which equations (73-77) are obtained as we referred before. In contrast to Sec. III B, when studying the stabilities of the homogeneous solutions, we do not need to consider the x direction, so equations (73-77) get simplified: Similar to the case when studying the time evolution, all fields above, including ψ 1 , ψ 2 and A t , are related to the fields in equations (6-8) through coordinates transformation and the U (1) gauge transformation (69-72). Linearize equations (A1-A4), we will be left with condition for a t : − iω∂ z a t + (µψ * 2 + i∂ z ψ * 2 ) p 2 + (µψ 2 − i∂ z ψ 2 ) q 2 | z=0 = 0. (A16) Here, as the reductions of equations (A8-A11) at z = 1 contribute another four boundary conditions (natural boundary conditions), we have gotten enough boundary conditions to solve equations (A8-A12).
Transform the three homogeneous solutions, Solution S1, Solution S2 and Solution S1+S2, into their proper forms as we mentioned in the last paragraph, and then submit them into equations (A8-A12), ω can be solved numerically. And, from perturbations (A5-A7), it is obvious that if there exists an ω whose imaginary part Im(ω) > 0, the corresponding solution is unstable, otherwise it is stable. Our results are shown in Fig. 10, where we choose ρ = 3.1 and λ 12 = 0.4. In both Fig. 10a and Fig. Fig. 10b, imaginary parts of all solutions of ω are equal or less than 0, so it confirms us that these two solutions are stable or metastable.
In Fig. 10c, the Solution S1+S2 is studied, which seems to be stable as the former two solutions. However, when looking at the enlarged drawing of its central area (see Fig. 10d), it is obvious that an ω with imaginary part greater than 0 (≈ 0.0026) exists, so the Solution S1+S2 is unstable.

Appendix B: Equations of state
Once the chemical potential µ is given, we are able to solve static equations (6-8) numerically, from which density ρ (= − ∂ z A t | z=0 ) is obtained. In this way, we express ρ in terms of µ, and thus obtain numerical solutions to equations of state (84) and (85). Our results are shown in Fig. 11, where the blue, red and green curves correspond to equations of state for S1, S2 and S1+S2.
where L ξ is the Lie derivative with respect to ξ = (∂/∂t) and L m the Lagrangian in (5).
The grand potential is defined as where the gauge A t | z=1 = 0 should be adopted, the last two terms are counter terms, and γ µν is the induced metric. Expanding all the indices in (C2) and apply boundary equations for all the fields, we can reach the final form of the grand potential: where we have adopted ψ i = Ψ i /z.