Dynamics of Phase Separation from Holography

TL;DR

This paper uses holography to study the real-time evolution of phase separation in a strongly-coupled gauge theory, identifying four stages of the process and comparing hydrodynamic models' effectiveness.

Contribution

It provides a detailed holographic analysis of phase separation dynamics, including stages, stability, and hydrodynamic descriptions, with new insights into the limitations of existing models.

Findings

Four stages of phase separation identified via holography.

Second-order hydrodynamics accurately describes all stages.

Israel-Stewart hydrodynamics fails near zero sound speed points.

Abstract

We use holography to develop a physical picture of the real-time evolution of the spinodal instability of a four-dimensional, strongly-coupled gauge theory with a first-order, thermal phase transition. We numerically solve Einstein's equations to follow the evolution, in which we identify 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. On the gravity side the latter is described by a static, stable, inhomogeneous horizon. We conjecture and provide evidence that all static, non-phase separated configurations in large enough boxes are dynamically unstable. We show that all four stages are well described by the constitutive relations of second-order hydrodynamics that include all second-order gradients that are purely spatial in the local rest frame. In contrast, a M\"uller-Israel-Stewart-type formulation of hydrodynamics fails to provide a good description for two reasons. First, it misses some large, purely-spatial gradient corrections. Second, several second-order transport coefficients in this formulation, including the relaxation times $\tau_\pi$ and $\tau_\Pi$, diverge at the points where the speed of sound vanishes.

Figures

Equations

$$S=\frac{2}{{\kappa }_5^2}\int d^{5}x\sqrt{-g}\left[\frac{1}{4}\mathcal{R}-\frac{1}{2}{\left(\nabla \varphi \right)}^2-V(\varphi )\right],$$

$${\ell }^{2}V(\varphi )=-3-\frac{3{\varphi }^2}{2}-\frac{{\varphi }^4}{3}-\frac{{\varphi }^6}{3{\varphi }_M^2}+\frac{{\varphi }^6}{2{\varphi }_M^4}-\frac{{\varphi }^8}{12{\varphi }_M^4},$$

$$\ell W(\varphi )=-\frac{3}{2}-\frac{{\varphi }^2}{2}-\frac{{\varphi }^4}{4{\varphi }_M^2},$$

$$V(\varphi )=-\frac{4}{3}W(\varphi {)}^2+\frac{1}{2}{W}^{\prime }(\varphi {)}^2.$$

$${\mathcal{E}}_{\text{high}}\simeq 5.9\times \frac{{\Lambda }^4}{10^2},{\mathcal{E}}_{\text{low}}\simeq 9.4\times \frac{{\Lambda }^4}{10^5}.$$

$$(\mathcal{E},P_L,P_T)=\frac{{\kappa }_5^2}{2{\ell }^3}(-T_t^t,T_z^z,T_{x_{\perp }}^{x_{\perp }}),$$

$${\mathcal{E}}_{\text{high}}/{\mathcal{E}}_{\text{low}}\simeq 0.6\times 10^3.$$

$$c_s^2=\frac{s}{c_v}.$$

$$\mathcal{A}\sim e^{-i\omega (k)t},$$

$${\omega }_{\pm }(k)=\pm \sqrt{c_s^2{k}^2-k^4{f}_L-k^4\frac{1}{4}{\Gamma }^2}-\frac{i}{2}\Gamma k^2.$$

$$\Gamma =\frac{1}{T}\left(\frac{4}{3}\frac{\eta }{s}+\frac{\zeta }{s}\right)$$

$$f_L=-\frac{c_s^2}{\mathcal{E}+P}\left(\zeta {\tau }_{\Pi }+\frac{4}{3}\eta {\tau }_{\pi }\right).$$

$$\gamma (k)\equiv -i\omega (k).$$

$$\gamma (k)=\left|c_s\right|k-\frac{1}{2}\Gamma k^2+k^3\frac{4f_L+{\Gamma }^2}{8\left|c_s\right|}+O(k^5).$$

$$k_{\ast }\simeq \frac{2\left|c_s\right|}{\Gamma }.$$

$$k_n=\frac{2\pi n}{L}.$$

$$\{{\mathcal{E}}_1,{\mathcal{E}}_2,{\mathcal{E}}_3,{\mathcal{E}}_4,{\mathcal{E}}_5\}\simeq \{2.6,1.9,1.6,1.2,1.0\}\times \frac{{\Lambda }^4}{10^2}.$$

$$f_L({\mathcal{E}}_2)\simeq f_L({\mathcal{E}}_3)\simeq -0.136/{\Lambda }^2.$$

$$\gamma \simeq 2.5\times 10^{-4}\Lambda .$$

$$10^7\times \frac{{\mathcal{E}}_3}{{\Lambda }^4}\simeq a_{-}{e}^{{\gamma }_{-}t}+a_{+}{e}^{{\gamma }_{+}t}$$

$$a_{-}=-4.17,a_{+}=3.26,{\gamma }_{-}=-3.5\times 10^{-4}\Lambda ,{\gamma }_{+}=2.0\times 10^{-4}\Lambda .$$

$$\mathcal{E}-{\mathcal{E}}_{\text{high}}\propto \exp [-i({\omega }_R+i{\omega }_I)t].$$

$$\omega (k)\simeq \left|c_s\right|k-\frac{i}{2}\Gamma k^2,$$

$$\mathcal{E}(t,z)\simeq\tilde{\mathcal{E}}_{\textrm{\tiny final}}(z,\Delta z(t))+A_{0}\,e^{{\rm{Im}}\,\omega(t-t_{0})}\,\cos\big{[}{\rm{Re}}\,\omega(t-t_{0})+\psi_{0}\big{]}\,\cos\big{[}k(z-z_{0})\big{]}\,.$$

$$\tilde{\mathcal{E}}_{\textrm{\tiny final}}(z,\Delta z(t))\equiv\left\{\begin{array}[]{lcc}\mathcal{E}_{\textrm{\tiny final}}(z+\Delta z(t))&\mbox{if}&z\leq z_{0}\,,\\ \mathcal{E}_{\textrm{\tiny final}}(z-\Delta z(t))&\mbox{if}&z\geq z_{0}\,,\end{array}\right.$$

$$2\mathcal{E}_{\text{high}}\Delta z(t)+\int_{z_{0}-\lambda/4}^{z_{0}+\lambda/4}dzA_{0}\,e^{{\rm{Im}}\,\omega(t-t_{0})}\,\cos\big{[}{\rm{Re}}\,\omega(t-t_{0})+\psi_{0}\big{]}\,\cos\big{[}k(z-z_{0})\big{]}=0\,,$$

$$\Delta z(t)=-\frac{1}{k\mathcal{E}_{\text{high}}}A_{0}\,e^{{\rm{Im}}\,\omega(t-t_{0})}\,\cos\big{[}{\rm{Re}}\,\omega(t-t_{0})+\psi_{0}\big{]}\,.$$

$$A_0/{\Lambda }^4=0.0006775,{\psi }_0=-0.031,z_0\Lambda =45.86.$$

$$\mathcal{E}(t,z)$$

$$+$$

$$A_0^{\prime }/{\Lambda }^4=1.125\times 10^{-6},{\psi }_0^{\prime }=2.85.$$

Full text

1. 1 Introduction

2. 2 Model and instability

3. 3 Time evolution

4. 3.1 Initial state

5. 3.2 The role of the box

6. 3.3 Linear regime

7. 3.4 End of the linear regime

8. 3.5 Reshaping

9. 3.6 Mergers

10. 3.7 Unstable static configurations

11. 3.8 Domain relaxation

12. 4 Phase separation

13. 5 Hydrodynamics

14. 6 Discussion

15. A Dispersion relation of sound modes

1 Introduction

Holography has provided numerous insights into the out-of-equilibrium properties of hot, strongly-coupled, non-Abelian plasmas. Examples include [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 18, 19, 17, 20, 21, 22, 23, 24] (see e.g. [25] for a review). In most of these cases the final state of the plasma at asymptotically late times is a homogeneous state. The purpose of this paper is to analyse a case in which the final configuration is expected to be an inhomogeneous state exhibiting phase separation. Specifically, we will study a four-dimensional gauge theory with a first-order, thermal phase transition. We will place the theory in an initial homogeneous state with an energy density in the unstable, spinodal region (see Fig. 1). If this state is perturbed, the system will evolve to a final state that will necessarily be inhomogeneous. Following this real-time evolution with conventional quantum-field theoretical methods for an interacting gauge theory is challenging. Therefore we will use holography, in which case following the evolution can be done by solving the time-dependent Einstein’s equations on the gravity side. A first study of this system was presented in [26]. Subsequently, [27] analysed a case in which the gauge theory is three-dimensional. Here we will extend the analysis of [26] in several directions and we will develop a detailed physical picture of the entire evolution of the system.

The real-time dynamics of phase separation in a strongly-coupled, non-Abelian gauge theory might be relevant to understand the physics of future heavy ion collision (HIC) experiments such as the beam energy scan at RHIC, the compressed baryonic matter experiment at FAIR and other experiments at NICA. These experiments will open an unprecedented window into the properties of the phase diagram of Quantum Chromodynamics (QCD) at large baryon chemical potential, which is expected to contain a line of first order phase transitions ending at a critical point [28, 29, 30]. If this scenario is realised then the real-time dynamics of the spinodal instability may play an important role, which provides one motivation for our work.

Hydrodynamics has been extremely successful at describing the quark-gluon plasma formed in HICs so far. If the future HIC experiments do explore a region with phase transitions, the applicability of the formulation of hydrodynamics used in most numerical codes might need to be reconsidered. In [26, 24] we initiated a study in this direction. We found that near a critical point, due to the slowing down of the dynamics, the system is accurately described by the constitutive relations of a formulation of hydrodynamics that includes all second-order gradients that are purely spatial in the local rest frame. However, this formulation is an acausal theory for which the initial-value problem is not well posed. A cure that is vastly used in hydrodynamic codes consists of exchanging the terms with second-order purely spatial derivatives in the local rest frame for terms with one time and one spatial derivative (see [31] for a review). We found in [26, 24] that the resulting theory, which we will call a Müller-Israel-Stewart-type (MIS) formulation, failed to provide a good description in the region near the phase transition. In this paper we find similar conclusions and we elaborate on the reasons for the failure of the MIS-type formulation. In particular, we show that some second-order transport coefficients, including the relaxation times $\tau_{\pi}$ and $\tau_{\Pi}$, diverge at the points where the speed of sound vanishes.

2 Model and instability

Our gravity model is described by the Einstein-scalar action

[TABLE]

with potential

[TABLE]

where $\ell$ is the asymptotic curvature radius of the corresponding AdS geometry and $\phi_{M}$ is a parameter that we will set to $\phi_{M}=2.3$. This potential can be derived from the superpotential

[TABLE]

via the usual relation

[TABLE]

The potential (2.2) is the same as in [26, 24]. The dual gauge theory is a Conformal Field Theory (CFT) deformed with a dimension-three scalar operator with source $\Lambda$. On the gravity side this scale appears as a boundary condition for the scalar $\phi$.

Our motivation to choose this model is simplicity. The presence of the scalar breaks conformal invariance. The first two terms in the superpotential are fixed by the asymptotic AdS radius and by the dimension of the dual scalar operator. The quartic term in the superpotential is the simplest addition that results in a thermal first-order phase transition in the gauge theory (for $\phi_{M}\leq 2.521$), which we extract from the homogeneous black brane solutions on the gravity side. In particular, the gauge theory possesses a first-order phase transition at a critical temperature $T_{c}=0.247\Lambda$, as illustrated by the multivalued plot of the free energy density as a function of the temperature in Fig. 1(left). Fig. 1(right) shows the corresponding energy density, where the high- and low-energy phases at $T_{c}$ have energy densities

[TABLE]

Notice that we work with the rescaled quantities

[TABLE]

where ${{P}}_{L}$ and ${{P}}_{T}$ are the longitudinal and transverse pressures with respect to the $z$-spatial direction along which the dynamics will take place. For an $SU(N_{\textrm{\tiny c}})$ gauge theory the prefactor on the right-hand side typically scales as $N_{\textrm{\tiny c}}^{2}$. Note the large hierarchy between the energy densities:

[TABLE]

States on the dotted red curves of Fig. 1 are locally thermodynamically unstable since the specific heat is negative, $c_{v}=d{\mathcal{E}}/dT<0$. As we now explain, these states are also dynamically unstable. This is in agreement with the Gubser-Mitra conjecture [32, 33], which states that a black brane with a non-compact translational symmetry is classically stable if, and only if, it is locally thermodynamically stable. We will comment again on this conjecture at the end of Sec. 3.7.

The connection with the dynamic instability was pointed out in an analogous context in [36] and it arises as follows. The speed of sound is related to the specific heat $c_{v}$ and the entropy density $s$ through

[TABLE]

Therefore $c_{s}^{2}$ is negative on the dashed red curves of Fig. 1, as shown in Fig. 2(left), and consequently $c_{s}$ is purely imaginary. The amplitude of long-wave length, small-amplitude sound modes behaves as

[TABLE]

with a dispersion relation given by (see Appendix A for the derivation)

[TABLE]

The plus sign corresponds to an unstable mode, while the minus sign leads to a stable mode. In this expression

[TABLE]

is the sound attenuation constant, $\eta$ and $\zeta$ are the shear and bulk viscosities, and $f_{L}$ is a second-order transport coefficient related to the coefficients $\tau_{\pi}$ and $\tau_{\Pi}$ in [48] through (see Sec. 5)

[TABLE]

In our model $\eta/s=1/4\pi$ [34], we compute $\zeta$ numerically following [35], and we obtained $f_{L}$ in [26, 24].

An imaginary value of $c_{s}$ leads to a purely real value of the growth rate

[TABLE]

For small momenta (2.10) yields for the unstable mode

[TABLE]

The first two terms alone give the familiar parabolas corresponding to the curves in Fig. 3. Note that these curves depend on the energy density $\mathcal{E}$ of the state under consideration because both $c_{s}$ and $\Gamma$ depend on $\mathcal{E}$.

We see that the growth rate is positive for momenta in the range $0\leq k<k_{*}$ with

[TABLE]

The corrections to these parabolas coming from the inclusion of the $k^{3}$ term in (2.14) or from evaluation of the full square root in (2.10) are small if we use the values of $f_{L}$ obtained in [26, 24]. Nevertheless, because $4f_{L}$ is very close to $\Gamma^{2}$ in magnitude but has the opposite sign, fitting these corrections provide an accurate method to extract this coefficient. In subsection 3.3 we will use this method to obtain a value of $f_{L}$ in good agreement with [26, 24].

In conclusion, states on the red dashed curves of Fig. 1 are afflicted by a dynamical instability, known as spinodal instability, whereby long-wave length, small-amplitude perturbations in the sound channel grow exponentially in time. The corresponding statement on the gravity side is that the black branes dual to the states on the dashed red curve are afflicted by a long-wave length instability. Although this is similar [36, 37, 38] to the Gregory-Laflamme (GL) instability of black strings [39], there is an important difference: In the GL case all strings below a certain mass density are unstable, whereas in our case only states on the red dashed curves of Fig. 1 are unstable.

3 Time evolution

3.1 Initial state

To investigate the fate of the spinodal instability we compactify the $z$-direction on a circle of length $L$ in the range $L\Lambda\in(107,213)$. For comparison, Ref. [26] considered $L\Lambda\simeq 57$. This infrared cut-off reduces the number of unstable sound modes to a finite number, since modes along the $z$-direction must have quantized momenta

[TABLE]

For simplicity, we impose homogeneity along the other two gauge theory directions $x_{\perp}$. We then consider a set of homogeneous, unstable initial states with energy densities $\{\mathcal{E}_{1},\mathcal{E}_{2},\mathcal{E}_{3},\mathcal{E}_{4},\mathcal{E}_{5}\}$ in the spinodal region given by

[TABLE]

For comparison, we have also listed the energy density $\mathcal{E}_{5}$ of the state considered in [26]. These states are indicated by the dashed, horizontal grey lines in Fig. 1(right). To trigger the instability, Ref. [26] introduced a small $z$-dependent perturbation in the energy density corresponding to a specific Fourier mode on the circle. However, this is not indispensable since numerical noise alone is enough to trigger the instability. Therefore in this paper we will consider both cases. Note that for the instability to play a role the size of the box must be large enough to fit at least one unstable mode. In other words we must have $L>2\pi/k_{*}$.

We follow the instability by numerically evolving the Einstein-plus-scalar equations as in [14, 18]. From the dynamical metric we extract the boundary stress tensor. We have performed 15 runs in which we observe phase separation at late times. We have published the corresponding boundary data extracted from the evolutions as open data [40] with a script to visualize each of them. The results for the energy density for two representative runs are shown in Fig. 4.

3.2 The role of the box

In order to avoid confusion in the discussion below, it is important to realize from the beginning in which ways the physics of the system can depend on the size of the box or, more precisely, on the dimensionless quantity $L\Lambda$. As we said above, the finite value of $L$ implements and IR cut-off on the allowed modes. Thus a first potential effect is that, if the box has size $L<2\pi/k_{*}$, then homogeneous states that would be dynamically unstable in bigger boxes become dynamically stable because no unstable modes fit in the box.

The dynamics of inhomogeneous states can also be crucially affected. Suppose for example that we have a configuration with a single domain (to be defined more precisely in Sec. 3.5) like that at late times in Figs. 4(left) or (right). Here we must clarify an issue of terminology. Because of our periodic boundary conditions, the number of domains is only unambiguously defined in relation to the size of the box. In other words, given a static configuration with one domain in a box of size $L$ we can take $n$ copies of it and generate a new configuration with $n$ domains in a box of size $nL$. The crucial point is that, although these two configurations are equivalent as static configurations, their dynamics once slightly perturbed may be radically different. Physically, the reason is that the configurations in Figs. 4(left) or (right) are phase-separated and hence stable, whereas $n$ copies of them in a box of size $nL$ are unstable towards merging of the different domains into a single one, as we will see in Sec. 3.7. Technically, the reason is again that some modes that are unstable in the box of size $nL$ do not fit in a box of size $L$. In fact, even the final stable domain in the $nL$-sized box will not fit in the smaller box if $n$ is large enough. This illustrates that, ultimately, the different dynamics is due to the fact that, while any configuration in a box of size $L$ can be viewed as a configuration in a box of size $nL$, the reverse is not true. The space of possible configuration and the dynamics in a bigger box are richer.

In summary, in each simulation we will specify and keep fixed the size of the box. We will also refer to the number of domains in the system, or speak of multi-peak configurations, or count the number of maxima of the energy profile, etc. with the understanding that these are meaningful, unambiguous concepts because we have a fixed, specific box size in mind.

3.3 Linear regime

Since the initial perturbation is small, the first stage of the evolution is well described by a linear analysis around the initial homogeneous state. Linear theory predicts a behavior which is the sum of two exponentials, precisely the two solutions of the sound mode (2.10). In the spinodal region, one of these modes decays with time while the other one grows. After some time the latter dominates. Figs. 5 and 6 show the time evolution of the amplitudes of several Fourier modes corresponding to the runs in Fig. 4. The straight lines in the logarithmic plots of Figs. 5 and 6 correspond to the regime of exponential growth.

In Fig. 7 we compare the growth rates predicted by linear analysis up to order $k^{2}$ with those extracted from a fit to the slopes of the straight lines in Figs. 5 and 6. We obtain good agreement, except for some particular cases. These correspond to resonant behavior, i.e. to the coupling between two modes which contributes to the growth of a third mode. For example, in Fig. 7(left) the 11th mode corresponds to a resonant behavior of the 1st and 10th modes. Other modes are also affected by resonant behavior, for example the 14th mode in Fig. 5(top) changes its slope at $t\Lambda\sim 650$ from the growth rate given by the sound mode to a new growth rate given by a resonance. As these may change in time, we have obtained the dots of Fig. 7 from early times, when the resonant behaviors are minimal.

The continuous black curves in Fig. 7 show the prediction of linear analysis to order $k^{2}$. As explained above, we can consider the full non-linear expression (2.10) to extract a value for $f_{L}$ by performing a fit. The dotted red curves in Fig. 7 show the results of these fits, from which we obtain the values

[TABLE]

These agree well with the values obtained in [26, 24].

Eqs. (2.9)-(2.10) determine the time evolution of each Fourier mode once two initial conditions are specified, for example its amplitude and its derivative. We have illustrated this with the dotted black curve shown in Fig. 5(top), which we obtained by fitting the initial conditions at $t=0$ for the $n=1$ mode, and which falls on top of the exact $n=1$ red curve. In principle, this could be done for every Fourier mode and obtain the full description of the system along the linear evolution. In practice, it is not possible to specify the precise initial conditions for the modes that are excited from the noise, but an estimate can be given by recalling that it is white noise. For example, in Fig. 6 a reasonable estimate would be obtained by assuming that the initial amplitudes are equal for all modes.

3.4 End of the linear regime

After the initial period of linear evolution, eventually the system enters into the non-linear regime. Of course, this time is not sharply defined. We choose to define it as the time at which the slope of the leading Fourier mode in the log plots of Figs. 5 and 6 deviates from the straight line predicted by the linear analysis by more than 10%. The resulting times are indicated by the grey vertical lines in Figs. 5 and 6. The corresponding energy profiles at those times are shown in Fig. 8 in solid blue. As can be seen in Figs. 5 and 6, the subleading modes can deviate earlier from the exponential growth predicted by the linear analysis due to resonant behavior.

A priori one may think that the linear regime ends simply when the amplitude of the inhomogeneous modes becomes a significant fraction of, but is still smaller than, the amplitude of the initial, homogeneous zero mode. However, our analysis indicates that this is too simplistic. On the one hand, in some cases the linear regime can persist until the amplitude in certain modes is so large that it is actually larger than that of the homogeneous mode. For example, this is the case in Fig. 6(top). Note that this does not mean that the energy density becomes negative in some regions, since the large negative contribution of the leading mode in these regions is compensated by the positive contributions of the other modes, as is clear from Fig. 8(right).

On the other hand, in some other cases we expect the linear regime to end when the amplitude of the inhomogeneous modes is still arbitrarily small. The reason for this is that, in generic circumstances, we expect the final state of the evolution to be a phase-separated configuration that should maximise the total entropy given the total energy available in the box. This implies that, at the latest, the exponential growth of the inhomogeneous modes should cease around the time when the energy density profile reaches $\mathcal{E}_{\text{high}}$ or $\mathcal{E}_{\text{low}}$, since in the final state no region should have energy higher than $\mathcal{E}_{\text{high}}$ or lower than $\mathcal{E}_{\text{low}}$. Note that this is a global condition in the sense that it does not apply to an individual mode but to the full energy density, which is the sum of all the modes. In some cases this condition should cut off the growth of the inhomogeneous modes at arbitrarily small amplitudes. For example, in the case of a gauge theory with a first order phase transition with arbitrarily small latent heat the growth stops because the energy profile quickly reaches both $\mathcal{E}_{\text{high}}$ and $\mathcal{E}_{\text{low}}$. Similarly, in a generic theory in a homogeneous initial state with energy very close to the upper or to the lower endpoints of the unstable branch, the growth stops because the profile of the energy density quickly reaches $\mathcal{E}_{\text{high}}$ or $\mathcal{E}_{\text{low}}$, respectively. We leave a more detailed investigation of the precise mechanism that cuts off the exponential growth for future work.

At the end of the linear regime, the exact number of maxima and minima of the energy profile is given by the leading Fourier modes at that time. These depend on the initial amplitude of the modes and on the growth rates, and can be determined from the initial conditions. Consider first the case in which the initial modes correspond to numerical noise. Since this is assumed to be white noise all modes start with similar initial amplitudes. Therefore in this case the modes with the largest growth rates will dominate. For example, we can see this in Fig. 6(top), where the $n=4$ mode clearly dominates, which is why we have 4 maxima and minima in Fig. 8(right). Now consider instead the case in which some initial mode is excited by hand with a large amplitude. If this amplitude is large enough then this mode will be the dominant one at the end of the linear regime. However, in some cases other modes with larger growth rates may still overtake this mode and become dominant at the end of the linear regime. This is illustrated in Fig. 5, in which the $n=1$ cosine mode is overtaken by the faster $n=6,8$ cosine modes and also by the $n=7$ sine mode. The latter is actually the dominant mode at the end of the linear regime, which is the reason why there are 7 peaks in Fig. 8(left).

3.5 Reshaping

By “reshaping” we mean a stage of non-linear evolution immediately after the end of the linear regime in which energy keeps being redistributed in the system in such a way that the structures formed during the linear regime keep adjusting their shape. This adjustment may or may not include a change in the number of maxima and minima. As we will see in Sec. 3.6, the reshaping period results in the formation of some structures that are either static or move with respect to one another with slowly varying velocity and almost constant shape.

There are two qualitative possibilities for the type of structures that can be formed at the end of the reshaping period: peaks and domains. By peaks we mean Gaussian-looking profiles as those around the maxima of the energy density in Fig. 8. By domains we mean plateaus in which the energy density is approximately constant and equal to either $\mathcal{E}_{\text{high}}$ or $\mathcal{E}_{\text{low}}$. If we need to distinguish we will refer to these as “high-energy domains” and “low-energy domains”, respectively. If we simply use the term “domain” we mean a high-energy domain. Some low-energy domains are present in Fig. 8, and examples of both types are shown in Fig. 9.

The distinction between a peak and a domain is of course not a sharp one, since the size of the domain can be reduced continuously until it turns into a peak.

We define the end of the reshaping period as the time beyond which the energy contained in each peak and domain, defined as the integral between the inflection points in their profiles, does no longer change by more than 5%. This is illustrated in Fig. 10, where we plot these integrated energies

In our model the reshaping period takes a few hundred times $\Lambda^{-1}$, as can be seen in Figs. 8 and 9. After this time, the peaks and/or domains move rigidly with slowly varying velocity with almost no distortion of their shape. This can be seen in Fig. 11, which is a zoom on early times of Fig. 4. The 7 initial peaks in Fig. 11(left) correspond to the 7 peaks shown in dashed green in Fig. 8(left). Similarly, the 4 peaks formed at early times in Fig. 11(right) correspond to the 4 peaks in dashed green of Fig. 8(right). To confirm that the profile of the peaks moves almost undeformed, in Fig. 12(left) we plot the profile of one of these peaks at different times. In Fig. 12(right) we plot the same profiles shifted by a constant amount to check that they coincide. The first snapshot in Fig. 12 corresponds exactly to the peak indicated with an arrow in Fig. 8(left) at the end of the reshaping period. We see that, on the scale of this plot, the shape is indeed “frozen” after this moment, and that it moves rigidly. However, a finer analysis reveals that neither the shape nor the velocity of the peaks are strictly constant in time, hence our use of the terms “almost constant shape” and “slowly varying velocity” above. Specifically, we have verified that the maximum local energy density of a peak stays constant to within 1% over very long periods of time of order $\Lambda\Delta t\sim 15000$. Over a similar amount of time the velocity of a peak can decrease by a factor of two. Presumably, the almost constancy of the shape and the relatively slow variation of the velocity rely on two features. One is the fact that the initial velocity of the peaks and/or domains is not too high, and the other is the fact that the hierarchy between the stable equilibrium energy densities is large, $\mathcal{E}_{\text{low}}/\mathcal{E}_{\text{high}}\ll 1$, which implies that a peak or domain can move through the cold phase with little “friction”. Presumably the distortion in the shape of the moving peaks and/or domains and their slowing down will become more pronounced as the ratio $\mathcal{E}_{\text{low}}/\mathcal{E}_{\text{high}}$ increases. We leave further investigation of this point for future work.

3.6 Mergers

In generic situations the different structures formed at the end of the reshaping period move with different velocities, so they eventually collide with one another. In Fig. 11 we can see several of these collisions. In all the cases that we have considered in this paper, a collision of any two structures leads to a merger, i.e. the result of the collision is a single, larger structure.111However, preliminary investigations indicate that for sufficiently high velocities the result of the collision consists of two objects moving away from each other. The merger between two peaks may result in a larger peak or a domain. The merger between a peak and a domain or between two domains results in a larger domain. A particularly clear illustration is shown in Fig. 13, which has the largest box that we have considered, $L\Lambda\simeq 213$, and an initial $n=1$ mode. The large box allows several peaks to coalesce to form two separate domains that move with non-vanishing velocity, and that finally collide forming a larger domain.

The final structure formed after a collision relaxes to equilibrium through damped oscillations, as can seen in e.g. Fig. 11. Although these are oscillations around an inhomogeneous state, we will show in the next subsection that, if the final structure is a phase domain, then at late times these oscillations correspond to linearised sound mode oscillations around the high energy phase.

In the case of a collision between a peak and a domain the dynamics can be qualitatively understood even from early times due to the separation of scales provided by the different sizes of the colliding structures. Indeed, in this case the peak creates a perturbation on top of the domain. This perturbation suffers attenuation and widening as it travels from one side of the domain to the other, as illustrated in Fig. 14. It would be interesting to verify if the attenuation and the widening are the same as for a perturbation traveling on an infinite plasma with the same energy density, as one would expect for large domains. When the perturbation reaches the end of the domain it bounces off entirely in the sense that no energy escapes the domain. In other words, the interface between the domain and the cold phase acts as a rigid wall with negligible transmission coefficient. This is illustrated in Fig. 15, where we plot the same curves as in Fig. 14 shifted by a constant amount in order to show that the shape of the interface on the left side is hardly modified by the bouncing of the perturbation. In fact, the total energy comprised between the midpoints of the two interfaces is practically constant in time (within $0.1\%$) once the peak has merged with the domain.

After the perturbation has bounced back and forth a few times, the system is well described by the linearised sound mode, as we will verify in Sec. 3.8.

3.7 Unstable static configurations

In non-generic situations the result at the end of the reshaping period may be an almost-static configuration. Typically this happens when a single mode (and multiples thereof) completely dominates the configuration. The reason for this is that the positions of the maxima and minima of a single Fourier mode are time-independent. An example of this kind of situation with a dominant $n=3$ mode is illustrated by Fig. 2 of [26], and two further examples with $n=5$ and $n=2$ are shown in Fig. 16. Although these configurations seem static at late times on the scale shown in these plots, they are actually not, as we will now see.

Consider first the case of the simulation of Fig. 16(right), whose Fourier modes at very late times are shown in Fig. 17.

We see that the cosine modes approach constant values at late times whereas some sine modes, although very small in amplitude, are growing exponentially with growth rate

[TABLE]

Note that this is two orders of magnitude smaller than the typical growth rates of the unstable modes around the initial homogeneous state (see Fig. 7), which is the reason why the late-time part of Fig. 16(right) looks approximately static. This situation should be contrasted with that of phase-separated configurations. For example, the Fourier modes at very late times of the simulation of Fig. 4(left) are displayed in Fig. 18. In this case all modes approach constant values, showing that this configuration becomes truly static at late times.222We have checked on the gravity side that at asymptotically late times the vector field $\partial_{t}$ becomes both Killing and hypersurface-orthogonal. In other words, the configuration at late times is truly static and not just stationary.

Since this configuration is also periodic, by taking multiple copies in a bigger box we conclude that truly static configurations with multiple, equally spaced domains also exist (recall the discussion in Sec. 3.2), and the same is true for multi-peak configurations. Our analysis then suggests that all multi-domain configurations are actually in unstable equilibrium towards merging into a single domain. Presumably on the gravity side this is due to the gravitational attraction between the high-energy regions of the horizon at different points on the circle.

The physics of the simulation in Fig. 16(right) can therefore be understood as follows. The system starts in a homogeneous state perturbed by an unstable $n=2$ mode. The initial amplitude of this mode is very small compared to the homogeneous, $n=0$ mode, but very large compared to any other mode, including of course numerical noise. Therefore time evolution initially takes the system very close to a static configuration with two domains on antipodal points of the $z$-circle. Note that in the truly static configuration all the sine modes would be exactly zero by symmetry. However, although they are very small, they are non-zero in our case because they are generated by numerical noise in the initial, homogeneous state. These modes can thus be viewed as perturbations of the exactly static configuration with two antipodal domains. Since some of these modes are unstable, they drive the system away from the antipodal configuration.

The evolution proceeds by moving the two domains towards each other, by simultaneously compressing them while keeping the shape of the interfaces fixed, and by increasing the energy density in the low-energy regions. To illustrate this in Fig. 19(left) we plot the energy profiles of the central domain at two late and widely separated times.

The small but appreciable relative displacement between the two curves shows that the central domain is moving towards the right. A similar plot shows that the other domain is moving towards the left, as indicated by the arrows in Fig. 20(top).

The compression of the domains and the rigidity of the interfaces are illustrated by the middle and bottom rows of Fig. 20, which are produced as follows. In the time interval between the two times shown in Figs. 19(left) and 20(top) the central domain moves to the right an amount $\Lambda\Delta z=1.380$ (defined as the average motion of the two interfaces) and the size of the domain decreases by an amount $\Lambda\Delta\ell=0.0264$ (defined as the relative motion of the two interfaces). Thus to produce Fig. 20(middle) we shift the red, dashed curve to the left by an amount $\Delta z-\Delta\ell/2$, so that the inflection points of the interfaces on the right-hand side of the central domain are on top of one another. Then we plot the difference between the shifted curve and the continuous, black curve. We see that the result vanishes in the centre of the domain and also on the location of the right interface. This shows that the value of the energy density in the centre of the domain has remained constant and that the shape of the right interface has not changed. The fact that the result is negative at the location of the left interface shows that the domain has decreased in size. To produce Fig. 20(bottom) we repeat the procedure except that we shift the red, dashed curve to the left by an amount $\Delta z+\Delta\ell/2$, so that now it is the left inflection points of the central domain that fall on top of one another. The result shows that the left interface has also moved with constant shape. We obtain analogous results for the other domain. Since both domains decrease in size the energy that they carry also decreases. Figs. 19(left) and 20(top) show that this excess energy is transferred to the low-energy region between the domains, whose energy density clearly increases as the domains approach each other. Due to the different velocities of the interfaces the cold region shape differ between being squeezed or opened up, where minimas appear. The direction of motion of the domains is consistent with the fact that the longitudinal pressure is lower in the low-energy region towards which the domains are moving than in the region that they leave behind, as shown in Fig. 19(right). The minimum of the longitudinal pressure deviates by more than $10\%$ from the critical pressure, $P_{c}$. Mechanically speaking, the high pressure in one region pushes the domains towards the low-pressure region.

The picture above is also consistent with several features of the unstable modes in Fig. 17. First, all cosine modes are stable, because adding a cosine perturbation to the antipodal configuration does not displace the domains towards each other but instead changes the distribution of energy between the two domains. While a large perturbation of this type could potentially take the system towards a single-domain configuration, this does not seem to lead to an instability at the linear level. Second, sine modes with even mode numbers are also stable. In this case this is due to the fact that a perturbation of this type shifts the position of the antipodal points simultaneously in the same direction and hence does not change the relative distance between them. Third, all the sine modes with odd mode number grow exponentially with the same growth rate (3.4). The reason why these modes are unstable is that they do change the relative distance between the domains. The reason why they all grow with the same growth rate is a consequence of the rigid motion of the two domains.

We now turn to Fig. 16(left), whose Fourier modes are shown in Fig. 21.

Since the modes are a combination of one growing and one decaying exponential and the times are not long enough, the unstable modes do not yet appear in the figure as straight lines. For example, Fig. 22 shows a fit to the unstable, $n=3$ cosine mode of the form

[TABLE]

with

[TABLE]

Note that the growth rates are of the same order of magnitude as (3.4).

Fig. 2 of ref. [26] shows an $n=3$ analog of the $n=5$ and $n=2$ runs of Fig. 16, namely a seemingly static, triple-peak configuration. The exactly static configuration has been constructed by directly solving a manifestly static problem in [41]. Fig. 23 shows that, contrary to what was stated in [26], this static configuration is unstable, as expected from the discussion in this section. Indeed, on the scale shown in the figure, the triple-peak configuration appears static only in the range of times $300\lesssim t\Lambda\lesssim 4500$. At $t\Lambda\sim 5500$ the first and second peaks merge, while the third peak is slowly moving left. At $t\Lambda\sim 12500$ all three peaks have merged into a single peak. The instability of the triple-peak configuration was missed in [26] because that reference only explored times of order $t\Lambda\lesssim 500$.

In summary, we conclude that the simulations shown in Fig. 16 are slowly evolving but not static because the corresponding static configurations are unstable. Therefore for sufficiently long times the different structures will presumably merge. The configurations in Fig. 16 are dominated by an $n=5$ and $n=2$ mode, respectively. We have performed an analogous analysis for configurations with $n=3$, which is the case of Fig. 2 of [26], and with $n=4,6,7$. In all cases the conclusion is the same. This leads us to conjecture that, given a fixed-size box (recall the discussion in Sec. 3.2), the only stable configurations are those with a single structure. In particular, all static, non-phase separated configurations in large enough boxes should be dynamically unstable. For large enough boxes the only stable states should be phase-separated configurations, which we will study in detail in Sec. 4, whereas for smaller boxes they would correspond to configurations with a single peak.

We close this section by pointing out that the conclusion about the instability of multi-peak configurations is outside the scope of the Gubser-Mitra conjecture [32, 33]. Recall that this states that a black brane with a non-compact translational symmetry is classically stable if, and only if, it is locally thermodynamically stable. The non-compactness assumption is violated since the $z$-direction is periodically identified, and the translational symmetry assumption is violated because the multi-peak configurations are inhomogeneous. The system is still translationally invariant along the transverse directions, but these are simply spectator directions and one could periodically identify them, thus making them also compact.

3.8 Domain relaxation

After two structures merge and form a phase domain, or after the latter forms directly at the end of the reshaping period, the domain oscillates and relaxes to equilibrium. Although the merger is a non-linear process, we will show that the subsequent relaxation can be very well described by linear theory, in particular by the linear sound mode perturbations around the high-energy phase. This may seem surprising given that the full configuration is inhomogeneous. However, as we have already seen in Fig. 14 and as we will further confirm below, the interfaces at the end of the domain behave as rigid walls, thus effectively confining the oscillations to the interior of the domain.

Soon after its formation, the relaxation of the domain is controlled by the largest sound mode that fits within it. To see this, in Fig. 24 we examine the two mergers that take place in Fig. 4(left) between the central domain and the two peaks that hit it from the right (i.e. from larger values of $z$).

Specifically, in Fig. 24(left) we show the time evolution of the energy density at a constant position $z\Lambda=46$, which corresponds to the center of the final domain. We see two regions of exponentially damped oscillations corresponding to the relaxation of the domain after each of the two hits. In each of them we extract the imaginary part of the frequency from the dampening coefficient, i.e. from the slope of the straight lines in the figure, and the real part of the frequency from the period of the oscillations. In other words, in each region we perform a fit to an amplitude of the form

[TABLE]

We then compare the result to that predicted by the dispersion relation (2.10) expanded to quadratic order:

[TABLE]

with $c_{s}$ and $\Gamma$ evaluated at $\mathcal{E}=\mathcal{E}_{\text{high}}$. Once (3.8) is assumed, the comparison of the real and imaginary parts of $\omega$ with the values extracted from the fit determines two independent values of $k$. These agree excellently with one another (within 0.7%), as illustrated in Fig. 24(right), so we will not distinguish between them. Each value of $k$ has an associated wavelength given by $\lambda=2\pi/k$. Remarkably, in each case the size of the domain measured between midpoints of the interface agrees almost exactly with 1/2 of the corresponding $\lambda$. For the second merger this is illustrated in Fig. 25(left), where we see that the vertical lines, which we have drawn at a distance $\lambda_{2}/2=50/\Lambda$ from one another, intersect the interfaces at their midpoints. In Fig. 25(right) we illustrate once more the rigidity of the interfaces by shifting by a constant amount the curves of Fig. 25(left). This rigidity implies that, effectively, the oscillations obey Dirichlet boundary conditions at the ends of the domain, which is the reason why the size of the domain equals $\lambda/2$ as opposed to $\lambda$.

In the previous analysis we have shown that the oscillations are associated with the longest wave-length sound mode of the high energy phase that fits within the phase domain. This suggests the possibility of describing not just the time dependence of the energy density at the center of the domain but the full spacetime evolution of the oscillations by using the analytical expression for the sound mode. Mathematically, this means that the form of the energy density in the domain should be of the following form:

[TABLE]

The second term on the right-hand side describes the oscillation in time and space of the sound mode, where $t_{0}$ is an arbitrarily chosen time, $\psi_{0}$ is the phase at $t_{0}$, $z_{0}$ is the center of the domain and $A_{0}$ is the amplitude at $t_{0}$. We now explain the meaning of the first term on the right-hand side.

We have observed that the interfaces move rigidly left and right as the phase domain oscillates, see Fig. 25(right). Physically, this motion is a consequence of the fact that the total energy inside the phase domain remains constant during the oscillations. Thus, when the cosine of the sound mode oscillates downwards, the two interfaces must move outwards to keep the energy constant, and vice versa, as shown in Fig. 25(left). If we call $\Delta z(t)$ the displacement of each interface, then mathematically this means that the oscillations happen on top of a domain with size $\lambda/2+2\Delta z(t)$ whose energy density can be written approximately as

[TABLE]

with $\mathcal{E}_{\textrm{\tiny final}}(z)$ the static domain profile at asymptotically late times. This formula is just a simple way of “stretching” (for positive $\Delta z$) or “compressing” (for negative $\Delta z$) the domain profile by “gluing in” or “cutting out” a small piece at the centre of the domain, taking advantage of the fact that the energy density is almost exactly constant there. In order to determine $\Delta z$ at each time, we simply impose conservation of energy, namely that the energy change associated to the rigid shift of the interfaces is exactly compensated by the energy change associated to the oscillations of the sound mode:

[TABLE]

where we recall that $\lambda=2\pi/k$. The integral is trivial, and solving for $\Delta z$ we find

[TABLE]

Note that, strictly speaking, the change in time of the size of the domain implies that the value of $k$ in (3.9), and through the dispersion relation also the value of $\omega$, depend on time. However, correcting these values would result in second-order effects since the second term in (3.9) is itself small to begin with.

For concreteness, let us consider applying (3.9), with $\Delta z(t)$ given by (3.12), to the case of the final phase domain of Fig. 4(left). The parameters $\rm{Re}\,\omega_{2}$, $\rm{Im}\,\omega_{2}$ and $k_{2}$ were already computed. We fit the different parameters at late times $t_{0}\Lambda=6000$ and obtain

[TABLE]

With these values we find a very good description of the full profile at earlier times, specifically within $1\%$ from $t\Lambda\simeq 4860$ (or $0.1\%$ from $t\Lambda\simeq 5430$) to the end of the evolution. We illustrate these results in Fig. 26 for two specific times.

Above we have only included the mode with the longest wavelength that fits in the domain. This mode is also the slowest mode to decay, so it is the dominant one at late times. At earlier times the description can be improved by including higher modes. Let us illustrate this by including the second mode. In this case (3.9) is replaced by

[TABLE]

The second mode has wavelength $\lambda_{2}/2$ and, by virtue of (3.8), double $\rm{Re}\,\omega_{2}$ and quadruple $\rm{Im}\,\omega_{2}$. It also oscillates in $z$ as a sine as opposed to a cosine, as expected on general grounds. We fit the new parameters with the result

[TABLE]

With this we obtain a good description of the system, specifically within $1\%$ from $t\Lambda=4690$ (or within $0.1\%$ from $t\Lambda=5210$). Note that these times precede those below (3.13), meaning that adding the second mode improves the description at early times. This improvement is also visible in Fig. 27, where we can see that at early times the second mode captures the odd (under $z\to-z$) component of the oscillation. Notice that the second mode does not contribute to (3.11) because the integral of the sine over $z$ vanishes.

In summary, we conclude that soon after the merger the system is well described by linearized hydrodynamics. In Sec. 5 we will verify that, in fact, the full evolution, from the initial spinodal instability to the final state, is well described by non-linear hydrodynamics.

4 Phase separation

Provided that the box is large enough and that the initial conditions are generic, the end state of the spinodal instability is a fully phase-separated configuration consisting of one high-energy domain, one low-energy domain and the interfaces between them, as in Figs. 4 and 9(left). This configuration is expected to maximise the entropy given the available total energy and the box size. The entire system is at rest since the net momentum in the initial configuration was zero.

In Fig. 28(left) we plot the energy profiles at late times of several simulations, together with the high- and low-energy phases obtained from the thermodynamics of homogeneous configurations. The good agreement between the latter and the energy densities of the domains confirms that this is a phase-separated configuration. Moreover, from the surface gravity of the horizon on the gravity side we obtain a temperature that is constant and equal to $T_{c}$ (within $0.01\%$) across the entire configuration, as expected from phase coexistence. Also as expected, we find that the interface that separates one phase from the other is universal, meaning that it is a property of the theory, independent of the initial conditions and of the size of the box. This is clearly demonstrated in Fig. 28(right), where we show that shifting each curve by a constant amount all the interfaces agree with one another.

As shown in Fig. 29(top), the shape of this universal interface is very well approximated by the function

[TABLE]

where $\Delta{\mathcal{E}}=\mathcal{E}_{\text{high}}-\mathcal{E}_{\text{low}}$, $z_{0}$ is the point at which the energy density is exactly half way between $\mathcal{E}_{\text{high}}$ and $\mathcal{E}_{\text{low}}$, and $b\Lambda\simeq 2.75$ can be taken as a definition of the size of the interface.

The universality of the interface implies that, in a phase-separated configuration, the size of each domain is fixed by the size of the box and by the total energy in it.

The surface tension of the interface is defined as the excess free energy in the system, per unit area in the transverse directions $x_{\perp}$, due to the presence of the interface. In a homogeneous system the free energy density per unit volume is constant and equal to minus the pressure, $F=-P$. Our system is only homogeneous along $x_{\perp}$, so it is the transverse pressure that appears in this relation (see e.g. [42]), i.e. we have $F(z)=-P_{T}(z)$, and moreover both densities are $z$-dependent. The transverse pressure in the final phase-separated configuration is shown in Fig. 29(bottom). By definition, at $T=T_{c}$ the homogeneous, stable, high-energy and low-energy phases have the same free energy density $F_{c}$, and hence the same transverse pressure, $P_{c}/\Lambda^{4}\simeq 7.5\times 10^{-6}$. This is the value to which the transverse pressure $P_{T}(z)$ asymptotes in Fig. 29(bottom) away from the interface. In the absence of the interface the free energy density per unit transverse area in the box would be $LF_{c}=-LP_{c}$. The excess free energy per unit transverse area due to the presence of the interface is therefore

[TABLE]

where the factor of 1/2 is due to the fact that there are two interfaces in the box. Note that this surface tension is positive because $P_{T}(z)<P_{c}$ or, equivalently, because $F(z)>F_{c}$. This is consistent with the fact that the presence of gradients associated to the interface increases the free energy of the system, as expected on general grounds.

If the box is not large enough then the final state will certainly not be a phase-separated configuration. To be precise, the box must be able to fit at least two interfaces plus two domains of sizes at least as large as the interfaces so that they can be distinguished from the interfaces themselves. Fine-tuned initial conditions may also prevent the final state from being a phase-separated configuration, as in the examples shown in Fig. 16.

5 Hydrodynamics

We saw in Sec. 3.8 that the relaxation of a domain can be described by linearised hydrodynamics. We will now show that non-linear, second-order hydrodynamics describes the entire spacetime evolution of the system from the beginning of the spinodal instability to the final, phase-separated configuration. For concreteness we will illustrate this for the evolution shown in Fig. 4(left). Our discussion parallels that in Sec. 4 of [24] with additional details included.

In modern language we define hydrodynamics as a gradient expansion around local equilibrium that, at any given order, includes all possible gradients of the hydrodynamic variables that are purely spatial in the local rest frame. Let us refer to this as the purely spatial formulation. To second order the hydrodynamic stress tensor takes the form

[TABLE]

with

[TABLE]

and

[TABLE]

In these expressions $P_{\textrm{\tiny eq}}$ is the equilibrium pressure, $u^{\mu}$ is the fluid four-velocity, $\Delta^{\mu\nu}=\eta^{\mu\nu}+u^{\mu}u^{\nu}$ is the projector onto spatial directions in the local rest frame, and $\Pi_{(1)}^{\mu\nu}$ contains the first-order corrections, with $\eta$ and $\zeta$ the shear and bulk viscosities, respectively. The shear tensor is $\sigma^{\mu\nu}=\nabla^{<\mu}u^{\nu>}$, where $\nabla^{\mu}\equiv\Delta^{\mu\nu}\partial_{\nu}$, and $A^{<\mu\nu>}$ denotes the symmetric, transverse and traceless part of any rank-two tensor. As in other holographic models as e.g. [45], the bulk viscosity remains finite at the points where the speed of sound vanishes as a consequence of the large-$N_{c}$ approximation implicit in the holographic set-up [46].

All the second-order terms are contained in $\Pi_{(2)}^{\mu\nu}$. For the case of interest here of fluid motion in flat space in 1+1 dimensions its tensor and scalar parts may be expanded as

[TABLE]

In order to make contact with [26] we chose the basis of operators to be

[TABLE]

Part of the notation above is chosen to make contact with [47, 48] below. The coefficients $\tilde{c}_{1},\tilde{c}_{2},\tilde{b}_{2},\tilde{b}_{3}$ are known because they are related to the coefficients $c_{L},c_{T},f_{L},f_{T}$ of [26] through

[TABLE]

These four coefficients are shown in Fig. 30 as a function of the energy density. Note that they are finite and smooth at the points where the speed of sound vanishes. We will come back to this point below.

We have not computed the coefficients $\tilde{c}_{7},\tilde{b}_{4}$ but, as we will see, they are not needed in order to obtain a good hydrodynamic description of our system. As in [26, 24], the reason is that the operators $\tilde{\mathcal{O}}_{7}^{\mu\nu},\tilde{\mathcal{S}}_{4}$ are suppressed in the dynamical situation under consideration because they are quadratic in the fluid velocity. In the case of the evolution shown in Fig. 4 the absolute value of this velocity is everywhere smaller than 0.2 and typically no larger than 0.1, as can be seen in Fig. 31.

In Figs. 32, 33 and 34 we compare the longitudinal and transverse pressures $P_{L}$ and $P_{T}$ that we read off from the simulation on the gravity side with the second-order hydrodynamic pressures $P_{L}^{\textrm{\tiny hyd}},P_{T}^{\textrm{\tiny hyd}}$.

To obtain the latter we read off the energy density and the fluid velocity from gravity and we apply the constitutive relations (5.1) with $\Pi_{(2)}^{\mu\nu}$ given by (5.5) omitting the contributions of $\tilde{\mathcal{O}}_{7}^{\mu\nu}$ and $\tilde{\mathcal{S}}_{4}$. In Fig. 32 we see an excellent agreement with the exact pressures in the final, phase-separated configuration. Two aspects of this agreement are particular remarkable. First, the interface between the high- and the low-energy domains is well reproduced. Second, the tiny hydrodynamic longitudinal pressure, which is of order $10^{-6}\Lambda^{4}$, results from a huge cancellation between the equilibrium terms and the second-order gradient corrections, both of which are several orders of magnitude larger in the high-energy region. Presumably, the small differences between $P_{L}^{\textrm{\tiny hyd}}$ and the exact longitudinal pressure may be reduced by increasing the precision in the determination of the second-order transport coefficients.

Note that the configuration in Fig. 32 cannot possibly be described by first-order hydrodynamics. The reason is that all first-order terms are linear in the velocity, see (5.4a). Since this vanishes in the final configuration, in this case first-order hydrodynamics reduces to ideal hydrodynamics. In turn, this implies that the profile of the longitudinal pressure follows that of the energy density through point-wise application of the equation of state, as is clear in Fig. 32, in contradiction with the conservation of the stress-energy tensor, which for this case implies that $P_{L}$ must be constant. In contrast, in the static limit the constitutive relations for the pressures at second order become

[TABLE]

These expressions motivated the definition of the $c_{L},c_{T},f_{L},f_{T}$ coefficients in [26]. We conclude that it is the contribution of the second-order terms with purely spatial gradients that brings about the agreement between the exact pressures and the hydrodynamic pressures in the phase-separated configuration.

In fact, second-order hydrodynamics describes well not just the final spatial profile but the entire spacetime evolution of the system. This is illustrated in Figs. 33 and 34 which show, respectively, the spatial dependence at an intermediate time of the evolution, $t\Lambda=2000$, and the time dependence of the entire evolution at a particular point $z\Lambda=21$.

In Figs. 32, 33 and 34 we have also plotted the ideal (equilibrium) pressure, as well as the hydrodynamic pressures obtained by including only the first-order viscous corrections. These agree rather well with one another almost everywhere but fail to describe the exact pressures. This shows that the first-order terms are suppressed not just in the final, phase-separated configuration but all along the evolution of the system, and also that the second-order terms with purely spatial gradients are as large as the ideal terms.

The purely spatial formulation of hydrodynamics is an acausal theory for which the initial-value problem is not well posed. For second-order hydrodynamics, a cure that is vastly used in hydrodynamic codes consists of using the first-order equations of motion to exchange the terms with second-order purely spatial derivatives in the local rest frame for terms with one time and one spatial derivative (see [31] for a review). This results in what we will call a Müller-Israel-Stewart-type (MIS) formulation. We emphasize that, strictly speaking, what is known as the MIS formulation is the phenomenological approach introduced in [49, 50, 51], which is not second-order accurate. Building on it, different second-order-accurate formulations have been constructed [47, 48, 52], to which we will collectively refer as MIS-type formulations. The key point is that, while they differ from MIS and they may also differ from one another in certain details, all these formulations share the common property that, as a first step to make the initial-value problem well posed, a second-order spatial derivative is replaced by one time and one spatial derivative. Since these two sets of second-order terms differ by higher-order terms, the purely spatial formulation and the MIS-type formulations are equivalent if all gradients are small. Different deviations from the MIS-type formulations were first reported near the critical point of $\mathcal{N}=2^{*}$ [53] and for fluids with small viscosities [54]. Since second-order gradients are large in our situation, one may expect that the two formulations will differ, as we will now verify. We follow [48], which is completely general for a non-conformal neutral fluid (see [47] for the conformal case).

In 3+1 dimensions the tensor and the scalar parts of $\Pi_{(2)}^{\mu\nu}$ can be expanded in a basis of eight tensor operators $\mathcal{O}_{i}^{\mu\nu}$ and seven scalar operators $\mathcal{S}_{j}$, respectively [48]. For the case of fluid motion in flat space in 1+1 dimensions only the following operators of the basis chosen in [48] do not vanish identically:333Our definition of $\mathcal{O}_{1}$ differs from that in [48] but agrees with the one in [47].

[TABLE]

Note that the ${\mathcal{O}}_{7}^{\mu\nu}$ and $\mathcal{S}_{4}$ operators are the same in both basis. Moreover, in 1+1 dimensions we have

[TABLE]

showing that the number of independent operators is the same as in (5.5). The first-order equations of motion imply the following identities

[TABLE]

where $D=u^{\mu}\nabla_{\mu}$ is the time derivative in the local rest frame and the equal signs here mean equality up to third- and higher-order terms. These identities may be used to replace $\mathcal{O}_{1}^{\mu\nu}$ and $\mathcal{S}_{3}$ in the expansions of $\pi_{(2)}^{\mu\nu}$ and $\Pi_{(2)}$ in favor of the left-hand sides of (5.11) [48], thus replacing terms with two spatial derivatives in the local rest frame with terms with one time and one spatial derivative. Upon these replacement the expansions read

[TABLE]

where we have made use of (5.10) and we have labelled the second-order coefficients as in [48]. The coefficients in the expansion (5.12) can be related to those in (5.5) by changing from one basis of operators to the other, with the result

[TABLE]

where

[TABLE]

is the enthalpy. Upon this change the fact that in our dynamical situation

[TABLE]

translates into

[TABLE]

Using (5.16) in (5.12) we finally arrive at the MIS-type constitutive relations

[TABLE]

As shown in Figs. 32, 33 and 34 the second-order hydrodynamic pressures determined from these constitutive relations, $P_{L}^{\textrm{\tiny hydMIS}}$ and $P_{T}^{\textrm{\tiny hydMIS}}$, fail to describe the exact pressures. It is interesting that there are two different reasons for this. The first one is the fact that the coefficients $\eta\tau_{\pi},\lambda_{4},\zeta\tau_{\Pi},\xi_{4}$ entering the MIS-type formulation (5.17) diverge at the points where the speed of sound vanishes, as can be seen in Fig. 35.

This can be traced back to their relation (5.13) with the $\tilde{c}_{1},\tilde{c}_{2},\tilde{b}_{2},\tilde{b}_{3}$ coefficients that enter the purely spatial formulation (5.5). The fact that the latter are smooth and finite at the points where $c_{s}^{2}=0$ (see Fig. 30) shows that $\eta\tau_{\pi},\lambda_{4},\zeta\tau_{\Pi},\xi_{4}$ diverge as inverse powers of $c_{s}^{2}$ at those points. In turn, this results in the divergences of the MIS hydrodynamic pressures at the spacetime points in the evolution at which the energy density goes through a value such that $c_{s}^{2}({\mathcal{E}})=0$, as can be seen in Figs. 32, 33 and 34. In contrast, the hydrodynamic pressures predicted by the purely-spatial formulation are smooth and finite everywhere.

The second reason why the MIS-type formulation fails to reproduce the evolution of the system is that it does not include all the independent, spatial, second-order gradient corrections. Indeed, the operators in the purely spatial formulation (5.5)-(5.6) contain $\nabla_{z}^{2}{\mathcal{E}}$ and $(\nabla_{z}{\mathcal{E}})^{2}$ terms. Both types of terms are necessary in order to describe the evolution correctly. Instead, in the MIS-type formulation the $\nabla_{z}^{2}{\mathcal{E}}$ terms are absent because the only operators that contain them, $\tilde{\mathcal{O}}_{1}^{\mu\nu}$ and $\tilde{\mathcal{S}}_{3}$, have been eliminated in favor of the left-hand sides of (5.11), which contain crossed $\nabla_{t}\nabla_{z}$ derivatives. The effect of this replacement is most clearly illustrated by the late-time, static, phase-separated configurations. In these states the fluid velocity and all time derivatives vanish, so all first-order terms are zero and all second-order gradient corrections reduce to a combination of terms of the form $\partial_{z}^{2}{\mathcal{E}}$ and $(\partial_{z}{\mathcal{E}})^{2}$. Both of these are correctly captured by the purely-spatial formulation, as shown in (5.8). In contrast, in the MIS-type formulation the only non-vanishing second-order operators are $\mathcal{O}_{8}^{\mu\nu}$ and $\mathcal{S}_{6}$ in (5.17), which only include $(\partial_{z}{\mathcal{E}})^{2}$ terms. Incidentally, this also shows that the divergences in the MIS-type pressures in the phase-separated configuration at the points where $c_{s}^{2}=0$ are not due to the divergences of the relaxation times but to those of the $\lambda_{4},\xi_{4}$ coefficients.

6 Discussion

We have used holography to develop a detailed physical picture of the real-time evolution of the spinodal instability of a four-dimensional, strongly-coupled, non-Abelian gauge theory with a first-order, thermal phase transition.

We have identified several characteristic stages in the dynamics of the system. In the first, linear stage the instability grows exponentially. In the second stage the evolution is non-linear and leads to the formation of peaks and/or domains. In the third stage these structures move towards each other with approximately constant shapes and slowly varying velocities until they merge, forming larger structures. In the fourth stage the system relaxes to equilibrium through damped oscillations which can be described in terms of linearised sound modes. For large enough boxes the final state after all mergers have taken place is a phase-separated configuration with one high- and one low-energy domain. As expected on general grounds, the interface separating them is universal in the sense that it is a property of the theory and does not depend on the initial conditions that led to the phase-separated configuration. We computed the surface tension of the interface in terms of the microscopic scale of the theory with the result (4.2). We noted that this interface moved with little deformation in the final relaxation stage towards the phase separated configuration. It would be interesting to understand in detail the precise conditions behind this “rigidity” property.

If the perturbation of the initial, homogenous state is dominated by a single mode with mode number $n\geq 2$ then the system evolves through an intermediate, almost static state with $n$ peaks or domains. Exactly static solutions with these number of structures do exist and, in principle, upon time evolution the system comes arbitrarily close to them provided that numerical noise is sufficiently suppressed. However, since these multi-structure static configurations are unstable, the evolution eventually deviates away. We have shown that this instability precisely pushes the different peaks or domains towards each other so that the final configuration at asymptotically late times is a phase-separated configuration if the box is large enough, or a configuration with a single peak otherwise.

Remarkably, along the entire spacetime evolution of the system the pressures are well described by the constitutive relations of a formulation of second-order hydrodynamics in which all the gradient terms that are purely spatial in the local rest frame are included. In particular, the interface in the final phase-separated configuration is well described by this formulation. It is therefore interesting to place our system in the context of the dynamics of fluids with boundaries. A good discussion of this topic in modern language can be found in [43, 44]. The general idea is that one can formulate hydrodynamics in the presence of an interface or phase boundary. This has an associated stress tensor that can be derivatively-expanded just like the stress tensor for the bulk of the fluid. At the zeroth, non-derivative order the time-time component of this stress tensor is the surface energy, namely the energy per unit area associated to the interface. Similarly, the diagonal space-space components give the pressure of the interface. General considerations imply that this equals minus the surface tension that we computed in (4.2). At higher orders in the derivative expansion the stress tensor of the interface is characterised by a set of transport coefficients called surface transport coefficients. In general, these coefficient are completely independent from those in the bulk of the fluid. Our case is an exception because the fact that the entire phase-separated configuration is well described by the hydrodynamics of the bulk fluid implies that the surface transport coefficients could be computed in terms of the bulk transport coefficients.

Purely spatial hydrodynamics is known to be acausal. This was not an issue for us since we did not evolve in time the hydrodynamic equations but simply verified the constitutive relations, but it is an issue in situations in which hydrodynamics is the only available description. For this reason we also investigated the validity of an MIS-type formulation, in which acausality is remedied by replacing terms with second-order spatial derivatives in the local rest frame by terms with one time and one space derivative. In the limit of small gradients this produces an equivalent formulation at long wavelengths. However, in our system the spatial gradients are large and the result is not equivalent. This is one reason why the MIS-type formulation fails to reproduce the correct pressures. The other reason is that, unlike in the purely spatial formulation, several second-order coefficients in the MIS-type formulation, in particular some relaxation times, diverge at the points where the speed of sound vanishes, leading to a divergent prediction for the hydrodynamic pressures at those points.

Since the conclusions in the paragraph above are important, let us state them in slightly different terms. Specifically, some readers may wonder what is the merit of using two different basis (5.6) and (5.9) given that they are supposed to be equivalent and hence give the same result by construction. The point is that the last statement is only true provided two conditions hold: (1) that all the gradients are small, because the equivalence is only accurate up to third- and higher-order terms, and (2) that the transformation between the two basis is non-singular. In our case both conditions can be violated. The fact that the gradients are not small follows from the fact that second-order terms are as large as ideal terms, as we saw in Sec. 5. The failure of this condition means that the replacement (5.11) is not justified, since it neglects third- and higher-order terms that are large. The discrepancy between the left- and the right-hand sides of (5.11) becomes most extreme in the static, phase-separated configurations that we have considered, since in these the left-hand side vanishes identically whereas the right-hand side is generically non-zero. Condition (2) is violated at the points where $c_{s}^{2}=0$. At these points the two basis are not equivalent even if the gradients are small, because the transformation between the two basis becomes singular. This is most clearly exhibited by the fact the relations (5.13) between the corresponding coefficients in each base become singular due to the inverse factors of $c_{s}^{2}$. This feature is the one that is ultimately responsible for the divergences seen in the hydrodynamic pressures in the MIS formulation.

Although our model is a specific bottom-up model our analysis suggests that the qualitative physics that we have just described may be quite universal in situations in which the gradients are large and/or the speed of sound is small. The latter property is guaranteed to hold near a critical point, so our results may have important implications for experimental searches of the QCD critical point.

It would be interesting to allow for dynamics in the $x_{\perp}$ directions. While some details will change, we expect that some of the qualitative lessons that we have extracted will remain true. For example, unstable modes in the initial homogeneous state will grow exponentially with the same growth rates as in our analysis. The details of the reshaping period will of course be more complicated as they will involve a shape adjustment in several dimensions. However, they will presumably lead to the formation of structures that will move with almost constant shapes and slowly varying velocities, since these features only depend on the large ratio between $\mathcal{E}_{\text{high}}$ and $\mathcal{E}_{\text{low}}$. Finally, we expect that, in large enough boxes, the only stable configurations with average energy densities in the unstable region will be phase-separated configurations with a single high-energy domain.

It would be interesting to extend our analysis in several other directions, including a more systematic study of domain collisions or the inclusion of a conserved U(1) charge.

Acknowledgements

We thank Tomas Andrade, Jean-Paul Blaizot, Oscar Dias, Roberto Emparan, Jaume Garriga, Thanasis Giannakopoulos, Volker Koch, Guy Moore, Daniele Musso, Mikel Sanchez, Jorge Santos, Mikhail Stephanov and Benson Way for discussions. We thank the MareNostrum supercomputer at the BSC (project no. UB65) and the FinisTerra II at CESGA for significant computational resources. MA would like to thank the Helsinki Institute of Physics for their hospitality during late stages of this work. MA acknowledges support through H2020-MSCA-IF-2014 FastTh 658574, FPA2017-83814-P, Xunta de Galicia (Consellería de Educación), MDM-2016-0692 and HPC17YDH55 by HPC-EUROPA3 (INFRAIA-2016-1-730897). MZ acknowledges support through the FCT (Portugal) IF programme, IF/00729/2015. We are also supported by grants FPA2016-76005-C2-1-P, FPA2016-76005-C2-2-P, 2014-SGR-104, 2014-SGR-1474, SGR-2017-754 and MDM-2014-0369.

Appendix A Dispersion relation of sound modes

Let us assume that we start with a non-conformal fluid in equilibrium, such that its stress tensor is constant and characterised by some energy density ${\mathcal{E}}_{0}$ and pressure $P({\mathcal{E}}_{0})$. We now excite the fluid with small fluctuations in the sound channel, such that we can express the fluctuations in terms of $\delta{\mathcal{E}}$ and $v^{i}$, both small. Since the fluid is isotropic we can choose both the momentum $k$ and the velocity $v$ (in the sound channel) to lie along the $z$-direction. We can define the velocity and energy density fluctuations via the stress tensor components

[TABLE]

where

[TABLE]

is the enthalpy.

Using the general expression for the second-order stress tensor, the fluctuations of all other (relevant) stress tensor components is given by

[TABLE]

where we recall that $\Gamma$ is the sound attenuation constant (2.11) and $f_{L}$ is the transport coefficient in (5.7b) and (5.8a).

After performing a space Fourier transform we define

[TABLE]

In terms of these Fourier components the two conservation equations become

[TABLE]

Differentiating the first equation and substituting in the second one we obtain a wave-like equation for the energy fluctuations:

[TABLE]

This equation leads to a non-trivial dispersion relation. When $c^{2}_{s}<0$ there is an unstable mode with frequency given by

[TABLE]

In the small-frequency limit this yields

[TABLE]

Note also that the intercept of the unstable mode with zero, namely the edge of the unstable dome, becomes

[TABLE]

This implies that $f_{L}$ must be negative in order for the second-order hydrodynamic dispersion relation to become stable at some $k$. In addition to the unstable mode, there is a second, stable mode with frequency

[TABLE]

Equations (A.8) and (A.11) correspond to (2.10) in the main text.