Domain Collisions

Yago Bea,1 Jorge Casalderrey-Solana,2 Thanasis Giannakopoulos,3
David Mateos,2, 4 Mikel Sanchez-Garitaonandia,2 and Miguel Zilhão3

1 Department of Physics and Helsinki Institute of Physics, PL 64, FI-00014
University of Helsinki, Finland.

2 Departament de Física Quántica i Astrofísica and Institut de Ciències del
Cosmos (ICC),

Universitat de Barcelona, Martí i Franquès 1, ES-08028, Barcelona, Spain.

3 Centro de Astrofísica e Gravitação (CENTRA),

Departamento de Física, Instituto Superior Técnico, Universidade de Lisboa,

Av. Rovisco Pais 1, 1049-001 Lisboa, Portugal.

4 Institució Catalana de Recerca i Estudis Avançats (ICREA),

Passeig Lluís Companys 23, ES-08010, Barcelona, Spain.

###### Contents

  1. 1 Introduction
  2. 2 The model
  3. 3 Domains in motion
  4. 4 Domain collisions
    1. 4.1 Low velocity
    2. 4.2 Intermediate velocity
    3. 4.3 High velocity
  5. 5 Discussion

##  1 Introduction

Holography has provided numerous insights into the out-of-equilibrium
properties of hot, strongly-coupled, non-Abelian plasmas [1, 2, 3, 4, 5, 6, 7,
8, 9, 10, 11, 12, 13, 14, 15, 16, 21, 18, 19, 22, 17, 23, 24, 25, 26, 27, 20,
28, 29, 30, 31, 32] (see e.g. [33] for a review). Here we will apply it to a
four-dimensional gauge theory with a first-order, thermal phase transition.
Imagine placing 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. For
this reason, Refs. [18, 19, 20] used holography, i.e. they followed the
evolution by solving the time-dependent Einstein’s equations on the gravity
side. The evolution consists of four generic stages [20] visible in Fig. 2: 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 collide and merge; and a fourth stage in which the
system finally relaxes to a static, phase-separated configuration. In any
given evolution of this type, the velocities of the different domains are a
complicated function of the perturbation of the initial, homogeneous state.
This makes it difficult to perform a systematic study of the physics of the
collision as a function of the domain velocities. In this paper we will
overcome this difficulty by directly preparing initial states consisting of
domains moving towards each other at a wide range of velocities.

Figure 1:  Energy as a function of temperature for the theory considered in
this paper, with intrinsic scale $\Lambda$. Stable, metastable and unstable
states are represented by solid, dashed and dot-dashed curves, respectively.
The vertical grey line corresponds to the phase transition temperature.

$10^{2}$$\mathcal{E}/\Lambda^{4}$

Figure 2:  Evolution of the energy density as a function of time and space for
an initial homogeneous state in the spinodal region that is slightly perturbed
at $t=0$. Figure taken from [20].

On the gravity side the above dynamics translates into the physics of a time-
dependent, inhomogeneous horizon. Understanding this physics provides one
motivation for our work. For example, we will examine in detail how the area
of the horizon changes as two domains merge. Another motivation comes from
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 [34, 35, 36]. If this scenario is
realised then the real-time dynamics of the spinodal instability may play an
important role.

One technical aspect that will allow us to go significantly beyond [18, 20] is
our choice of scalar potential on the gravity side. The choice in those
references gave rise to a phase diagram in the gauge theory with a large
hierarchy between $\mathcal{E}_{high}$ and $\mathcal{E}_{low}$. This large
ratio resulted in a run-time for our code on the scale of months. In contrast,
here we will use the potential introduced in [37] and subsequently used in
[29, 31]. With this choice the ratio $\mathcal{E}_{high}/\mathcal{E}_{low}$ is
of order one and the run-time decreases to minutes.

We emphasize from the start that the domains we consider here should not be
confused with the bubbles created via a nucleation process in first-order
phase transitions. The expansion of those bubbles is driven by a pressure
difference between the inside and the outside of the bubble. Instead, in our
case the domains move simply because they are created with some initial
velocity. In fact, if a number of bubbles are nucleated in a large – but
finite-size – box, and if the average energy density in the box is between
$\mathcal{E}_{high}$ and $\mathcal{E}_{low}$, then the endpoint of the
evolution will be a phase-separated state with two domains with energies
$\mathcal{E}_{high}$ and $\mathcal{E}_{low}$.

##  2 The model

Our gravity model is described by the Einstein-scalar action

$S=\frac{2}{\kappa_{5}^{2}}\int d^{5}x\sqrt{-g}\left[\frac{1}{4}{\cal
R}-\frac{1}{2}\left(\nabla\phi\right)^{2}-V(\phi)\right]\,.$ (2.1)

Exclusively for simplicity we assume that the scalar potential $V(\phi)$ can
be derived from a superpotential $W(\phi)$ through the usual relation

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

Different choices of (super)potential correspond to different dual four-
dimensional gauge theories. As in [37, 29, 31] we choose

$\ell
W(\phi)=-\frac{3}{2}-\frac{\phi^{2}}{2}-\frac{\phi^{4}}{4\phi_{M}^{2}}+\frac{\phi^{6}}{\phi_{Q}}~{},$
(2.3)

where $\ell$ is the asymptotic curvature radius of the corresponding AdS
geometry and $\phi_{M}$, $\phi_{Q}$ are constants, which in this paper we will
set to $\phi_{M}=1.0$ and $\phi_{Q}=10$. The dual gauge theory is a conformal
field theory (CFT) deformed with a dimension-three scalar operator with source
$\Lambda$. This intrinsic scale in the gauge theory will set the
characteristic scale of much of the physics of interest. On the gravity side,
$\Lambda$ appears as a boundary condition for the scalar $\phi$. In the limit
$\phi_{Q}\to\infty$ the sextic term is absent and the model reduces to that in
[18, 20, 27]. The motivation for the choice here is simplicity: the dual gauge
theory is a non-conformal theory with a first-order phase transition with a
ratio $\mathcal{E}_{high}/\mathcal{E}_{low}$ of order one, as shown in Fig. 1,
and the dual gravity solution is completely regular even at zero temperature.
This model has $T_{critical}/\Lambda=0.3959$,
$\mathcal{P}_{critical}/\Lambda^{4}=0.1124$,
$\mathcal{E}_{high}/\Lambda^{4}=1.3736$ and
$\mathcal{E}_{low}/\Lambda^{4}=0.4526$.

For simplicity, in this paper we will allow for dynamics only along one of the
three spatial directions of the gauge theory, which we call $z$. In other
words, we impose translational symmetry along $x,y$. We will also refer to $z$
as the longitudinal coordinate and to $x,y$ as the transverse directions. We
will decompose the gauge theory stress tensor accordingly and work with the
rescaled quantities

$({\cal
E},{{P}}_{T},{{P}}_{T},{{P}}_{L},\mathcal{J})=\frac{\kappa_{5}^{2}}{2\ell^{3}}(-T^{t}_{t},T^{x}_{x},T^{y}_{y},T^{z}_{z},T^{t}_{z})\,,$
(2.4)

where ${\cal E}$ is the energy density, ${{P}}_{L}$ and ${{P}}_{T}$ are the
longitudinal and transverse pressures, and $\mathcal{J}$ is the momentum
density in the $z$-direction. For an $SU(N_{\textrm{\tiny c}})$ gauge theory
the prefactor on the right-hand side typically scales as $N_{\textrm{\tiny
c}}^{-2}$. In addition, we compactify the $z$-direction on a circle of length
$L$. This infrared cut-off is technically convenient since it reduces the
number of unstable modes of homogeneous states in the spinodal region to a
finite number. In addition, compactifying the $z$-direction on a circle also
brings about interesting new effects even for a single, boosted domain. We
will start by analyzing this case in Sec. 3 and then we will move to
collisions of two domains in Sec. 4. Note that we are slightly abusing the
word “domain” in the following sense. At the end of the evolution in Fig. 2
there are two domains with energy densities $\mathcal{E}_{high}$ and
$\mathcal{E}_{low}$. Nevertheless, we will often use the word “domain” to
refer specifically to the high-density phase. In this sense, we will sometimes
think of the end state of Fig. 2 as a domain surrounded by a low-density bath.

##  3 Domains in motion

In a domain collision, each of the participants is initially moving
independently of the other, surrounded by a low-density bath. In this section
we will study the dynamics of such a single domain. We will construct an
initial, out-of-equilibrium state with a non-uniform fluid velocity in the
$z$-direction and watch it evolve into a steady-state with uniform velocity.
In addition to its intrinsic interest, this will also help us understand the
construction of initial states in the case of domain collisions.

We start with a single, static domain obtained as the end-state of a
simulation like that in Fig. 2. We then modify the near-boundary fall-offs of
the metric functions on the initial time slice on the gravity
side111Specifically the subleading terms, which determine the expectation
values. so as to add a finite momentum density along the $z$-direction,
$\mathcal{J}(z)$. We keep the total energy fixed as we inject this momentum;
as we will see below this results in a decrease in the entropy of the system.
The specific $z$-dependence is chosen to be

$\mathcal{J}(z)=j\,\Big{(}\mathcal{E}(z)-\mathcal{E}_{low}(T_{c})\Big{)}\,,$
(3.1)

where $j$ is an appropriately chosen constant. At this initial step the
momentum density has vanishing support on the low-density region. In order to
be able to initialise our numerical code the initial value of $j$ cannot be
too large. However, choosing a small value, letting the system evolve for a
few time steps, increasing the value of $j$, and iterating this procedure,
allows us to eventually reach high velocities. The reason for the iterations
is that, after each change in the near-boundary fall-offs, a few time steps of
evolution allow the bulk functions to adapt to the new near-boundary fall-
offs. In particular, this allows our code to find the new horizon within our
numerical grid. We will define as the initial time, $t\Lambda=0$, the time at
the end of this iterative process. At each step in the iteration process a
fraction of the total momentum gets transferred to the low-density bath. As a
consequence, at the end of the process the bath will not be completely at
rest, as illustrated by the blue curve in Fig. 3(top). Nevertheless, at this
time most of the momentum density is still in the domain, where the fluid flow
velocity is larger than in the bath. It is precisely this feature that tells
us that the resulting state at $t\Lambda=0$ is out of equilibrium.

Figure 3:  Time evolution of the fluid velocity (top) and of the energy
density profile (bottom) for a domain in motion. To facilitate the comparison,
each profile has been shifted by a constant amount to match the domain
midpoint with $z\Lambda=0$. The maximum of the fluid velocity at $t\Lambda=0$
is $v_{max}=0.66$.

After $t\Lambda=0$ we do not inject any more momentum into the system. The
subsequent evolution is characterized by a significant energy and momentum
transfer between the domain and the bath. The initial state is so far from
equilibrium that it splits into several fragments that eventually merge back
to form a single domain, as illustrated in Fig. 4.

Figure 4:  Plots of the energy density for the evolution of an initial,
single-domain state that splits into several fragments that subsequently merge
back into a single domain.

The stationary state is reached around $t\Lambda\sim 400$. At this and at
later times the fluid velocity field is constant along the $z$-direction with
value $v=0.364$, as seen in Fig. 3(top). In other words, the “system moves as
a whole”. This is illustrated in Fig. 5, where we see that the energy density
profile remains constant in time once the steady state is reached to a
precision better than 1%.

---

Figure 5:  (Left) Snapshots of the energy density at late times for the same
initial state as in Fig. 3. (Right) The same profiles shifted by an
appropriate constant. (Bottom) Relative difference between the profiles.

This steady-state energy profile is narrower than that of the out-of-
equilibrium state at $t\Lambda=0$, as can be seen in Fig. 3(bottom), and it is
also narrower than that of the static, equilibrium state with the same total
energy, as illustrated in Fig. 6 for two different velocities.

Figure 6:  Steady-state energy profiles of domains with different velocities.
The total energy in the box is the same in all cases. The average energy is
$\bar{\mathcal{E}}/\Lambda^{4}=0.817$ and the box has size $L\Lambda=42.0719$.
Figure 7:  Comparison of the energy profile obtained by Lorentz-transforming
the curve in Fig. 6 (dashed, red curve) with a configuration obtained directly
as an inhomogeneous, static configuration with the same average energy (3.5)
in a box of size $L^{\prime}=\gamma L$ (solid, blue curve).

We see in the figure that the domains become narrower as the velocity
increases, and that the energy densities in the domain and in the bath are
higher than $\mathcal{E}_{high}(T_{c})$ and $\mathcal{E}_{low}(T_{c})$,
respectively. These features are consistent with the fact that, by
construction, all these configurations have the same total energy.

In the steady state, a comoving observer traveling with velocity $v$ sees a
fluid configuration at rest. This observer uses coordinates
$t^{\prime},z^{\prime}$ related to $t,z$ through the usual Lorentz
transformation

$t^{\prime}=\gamma(t-vz)\,,\qquad
z^{\prime}=\gamma(z-vt)\,,\qquad\gamma=(1-v^{2})^{-1/2}\,.$ (3.2)

The original identification $(t,z)\sim(t,z+L)$, translates into a non-trivial
identification for the comoving observer

$(t^{\prime},z^{\prime})\sim(t^{\prime}-\gamma vL,z^{\prime}+\gamma L)\,.$
(3.3)

In other words, for the co-moving observer the identification is not purely
along the spatial direction but it involves time as well. This imposes a non-
trivial constraint on any physical configuration seen by this observer, since
any such configuration must be invariant under (3.3). However, for
configurations that are static with respect to the comoving observer this
condition reduces to periodicity along the $z^{\prime}$-direction with period
$L^{\prime}=\gamma L$. This is the case for the fluid configuration obtained
by applying the Lorentz transformation (3.2) to the configurations of Fig. 6.
In particular, the energy density in the comoving frame takes the form

$\mathcal{E}^{\prime}=\gamma^{2}\left(\mathcal{E}+v^{2}P_{z}-2v\mathcal{J}\right)\,,$
(3.4)

and the average energy in this frame is

$\bar{\mathcal{E}}^{\prime}=\bar{\mathcal{E}}-v\bar{\mathcal{J}}\,.$ (3.5)

The result of applying Eq. (3.4) to the $v=0.187$ case of Fig. 6 is shown by
the dashed, red curve in Fig. 7. As expected, this matches the solid blue
curve, corresponding to a profile obtained directly as an inhomogeneous,
static configuration with average energy (3.5) in a box of size $L^{\prime}$.

We see that, in the case of a single domain, we could have obtained the steady
states of Fig. 6 by applying the inverse Lorentz transformation to the
corresponding static configurations. In contrast, in the case of two domains
moving towards each other that we will consider in the next section, it will
be crucial to be able to construct the moving domains via the method described
around Eq. (3.1), because in that case we will need to match two baths with
opposite velocities.

##  4 Domain collisions

We now turn to the construction of initial states whose evolution will result
in the collision of two domains. We start with a configuration consisting of
two static domains separated by some distance from one another, like the one
described by the dashed, red curve in Fig. 8(left). This can be constructed,
for example, by gluing together two individual domains along the low-energy
bath. We place the two domains at antipodal points of the $z$-circle, in which
case the resulting configuration is unstable but static [20]. If this
configuration is slightly perturbed then, upon time evolution, the domains
will approach each other and eventually collide. However, this collision will
typically happen at low relative velocity [20]. Therefore, in order to explore
the collision physics over a significant range of relative velocities, we will
modify the initial configuration following the procedure described around Eq.
(3.1). The crucial difference is that now we inject momentum densities of
equal magnitudes but opposite signs in each domain. As in the case of a single
domain, at the end of the iteration procedure some of the momentum is
contained in the bath but most of it is in the domains, as illustrated in Fig.
8(right). Also, the shape of the domains, described by the solid, blue curve
in Fig. 8(left), is modified with respect to the static situation. Note that,
by symmetry, the fluid velocity must vanish at the two middle points in
between the domains, i.e. at $z\Lambda=0$ and $z\Lambda=L/2$. This is the
reason why this initial state cannot be obtained by gluing together two
individual steady states of the previous section moving towards each other,
since at the middle points the baths would be moving with non-zero, opposite
velocities, and hence the velocity field would be discontinuous. One may think
that an alternative method would be to use superposition instead of gluing.
Specifically, we recall that, in the characteristic formulation that we use on
the gravity side (see [21] for details), the data on the initial time slice is
free. We can therefore superpose the initial data corresponding to two
individual domains moving towards each other. Although this procedure is
conceptually sensible, it is technically challenging to find the apparent
horizon in the resulting geometry. For this reason we followed the procedure
described around Eq. (3.1).

---
Figure 8:  (Left) Initial configuration (dashed red) used to generate two
domains in motion (solid blue) that will collide with each other. (Right)
Fluid velocity for the domains in motion on the left.

We performed a series of identical simulations except for the fact that we
varied the maximum velocity of the fluid in the initial state in the range
$v_{\mathrm{max}}\in(0,0.73)$. We use this velocity as one possible
characterization of the initial state. In all cases we fixed the distance
between the mid points of the domains to $d\Lambda\simeq 42$. We found three
qualitatively different dynamical regimes. For low values
$v_{\mathrm{max}}\lesssim v_{1}\simeq 0.09$ the domains initially slow down,
enter a quasi-static regime, and eventually collide and merge. For
intermediate values $v_{1}\lesssim v_{\mathrm{max}}\lesssim v_{2}\simeq 0.2$
the quasi-static phase is absent. For high values $v_{\mathrm{max}}\gtrsim
v_{2}$ the domains collide but do not merge. Instead they break apart and
merge only in subsequent collisions.

We emphasize that there is no sharp transition between these regimes and hence
the differences are qualitative in nature. For the specific simulations with
the domains of Fig. 8 we found the quoted values $v_{1}\simeq 0.09,v_{2}\simeq
0.2$. Repeating these simulations for different initial domains we found that
the values of $v_{1},v_{2}$ depend on the choice of initial domains. However,
the existence of three qualitatively different regimes persisted, suggesting
that it is a robust property of the collision dynamics.

###  4.1 Low velocity

Fig. 9 shows the evolution of the energy density and the fluid velocity for a
simulation with $v_{max}=0.08$. We can distinguish several stages of the
evolution: the initial slowing down of the domains, the quasi-static regime,
the merging into a single domain and the subsequent relaxation to equilibrium.

In the first stage, $t\Lambda\lesssim 100$, the domains slow down as they
transfer momentum to the low-energy bath in between them. This results in an
increase of the energy density and of the longitudinal pressure at the center,
as shown at early times in Fig. 10. We use the area density of the apparent
horizon on the gravity side as a proxy for the entropy density of the system.
As shown in Fig. 11, the dissipation in the initial phase causes a
considerable increase in the entropy of the system, meaning that it is
comparable to the entropy increase in the final merger.

Figure 9: Evolution of the energy density (top) and of the fluid velocity
(bottom) for a simulation with $v_{\mathrm{max}}=0.08$.

In the period $100\lesssim t\Lambda\lesssim 500$ the energy density (the
pressure) at the centre increases (decreases) slowly. The quasi-static nature
of this regime can be clearly seen in the entropy, which displays a plateau
that extends up to $t\Lambda\sim 500$. Note that the entropy of the plateau
coincides with that of the initial two-domain configuration (dashed, red curve
in Fig. 8(left)), which is static. Presumably, this agreement is due to the
fact that the velocity of the two domains in the quasi-static regime is
extremely low, as can be seen in Fig. 9(bottom). This is consistent with the
fact that the domains move as rigid bodies whose profiles are very well
approximated by the equilibrium profile, as shown in Fig. 12(left). The right
plot in this figure also shows that the dynamics is driven by a small pressure
deficit in the region in between the domains.

The fact that the entropy of the system at times $t\Lambda\lesssim 100$ is
lower than the entropy of two static domains means that the procedure by which
we inject momentum into the system at constant total energy decreases the
entropy. This can be easily understood on the gravity side, where this
procedure increases the total momentum of the dual black brane while keeping
its mass fixed. If one performs a Kaluza-Klain reduction along the
$z$-direction, this corresponds to increasing the electric charge of the
resulting, lower-dimensional black brane while keeping its mass fixed. Just in
a Reissner-Nordström black hole, this procedure is expected to reduce the
entropy.

---
Figure 10:  Energy density (left) and longitudinal pressure (right) as a
function of time at $z\Lambda=0$ for the collision of Fig. 9. Figure 11:  Time
evolution of the total entropy (per unit transverse area) for the collision of
Fig. 9. $S_{2}$ is the entropy of a static configuration with two domains (the
dashed, red curve in Fig. 8). $S_{1}$ is the entropy of a static configuration
with the same total energy but a single domain.
---
Figure 12:  (Left) Energy profiles (shifted by an appropriate amount) of the
domains during the quasi-static phase for the collision of Fig. 9. (Right)
Snapshot of the longitudinal pressure distribution at $t\Lambda=350$.

The collision takes place around $t\Lambda\sim 500$. The two domains merge
into a single, out-of-equilibrium, larger domain that will eventually reach
equilibrium through damped oscillations. At sufficiently late times we expect
to be able to describe the dynamics in terms of an equilibrium profile plus
small deviations, $\delta\mathcal{E}/\mathcal{E}\ll 1$, so that we can apply a
linear analysis. Following [20], we model the system at late times as two
phases of energies $\mathcal{E}_{low}$ and $\mathcal{E}_{high}$ separated by
two interfaces or walls that can move as a whole but whose shapes are rigid.
Since the system is symmetric under $z\to-z$ we will often focus only on $z>0$
and speak of “the” wall. We define the positions of the wall as the point at
which the energy density is the average between $\mathcal{E}_{high}$ and
$\mathcal{E}_{low}$, namely we require that
$\mathcal{E}(t,z_{wall})=(\mathcal{E}_{high}+\mathcal{E}_{low})/2$. The size
of the high-energy phase is then $\ell_{1}=2z_{wall}$ and the size of the low-
energy phase is $\ell_{2}=L-\ell_{1}$. Rigidity implies that the perturbations
of each phase must vanish at the wall, meaning that the only allowed
perturbations are odd cosine harmonics of wavelength $\lambda_{i}=2\ell_{i}$.
In general there could also be sine harmonic of wavelength
$\lambda_{i}=\ell_{i}$, but these are excluded in our particular case because
of the reflection symmetry. As both subsystems are coupled and the wall is
rigid we expect the two phases to oscillate in time as a whole, with identical
frequencies. Since translation invariance is broken, all the modes in the
Fourier decomposition of
$\delta\mathcal{E}=\mathcal{E}(t,z)-\mathcal{E}_{eq}(z)$, with
$\mathcal{E}_{eq}(z)$ the equilibrium profile at asymptotic times, will mix
with one another. Thus, at late times they should all oscillate with the
frequency of the longest-lived mode (see Fig. 25 in [29] for a related
discussion). We verified this for modes up to $n=14$, i.e. for modes with
momentum up to $k=14\times 2\pi/L$. Fig. 13 illustrates the result for the
four lowest modes.

Figure 13:  Time evolution of the spatial Fourier modes of
$\delta\mathcal{E}=\mathcal{E}(t,z)-\mathcal{E}_{eq}(z)$ after the merging.
All modes evolve as
$\delta\mathcal{E}_{n}\sim\exp(\omega_{1}t)\cos(\omega_{2}t)$ with the same
$(\omega_{1},\omega_{2})=(-0.012,0.063)$.

This result motivates the following ansatz for the perturbations of each phase
at late times:

$\delta\mathcal{E}_{i}(t,z)=\begin{cases}e^{\omega_{1}(t-t_{0})}\sum_{n=0}^{\infty}a_{n}^{i}\cos\Big{[}\omega_{2}(t-t_{0})+\gamma_{n}^{i}\Big{]}\cos\Big{[}\frac{(2n+1)\pi}{\ell_{i}}\left(z-z_{0}^{i}\right)\Big{]},\quad\left|z-z_{0}^{i}\right|\leq\frac{\ell_{i}}{2}\\\
0,\quad\mathrm{otherwise},\end{cases}$ (4.1)

where the index $i$ refers to each of the two phases and $z_{0}^{i}$ refers to
their midpoints. The form of the energy density at late times is then

$\mathcal{E}(t,z)=\mathcal{E}_{eq}\left(z\mp\frac{\Delta\ell_{1}(t)}{2}\right)+\delta\mathcal{E}_{high}(t,z)+\delta\mathcal{E}_{low}(t,z),$
(4.2)

where $\mp$ means that we use $-$ ($+$) if $z>0$ ($z<0$). The second and third
terms on the right-hand side describe the oscillation in time and space of the
two phases. The first term is due to the rigid motion of the walls as the
phase domains oscillate. If we call $\Delta\ell_{1}(t)=-\Delta\ell_{2}(t)$ the
changes in the sizes of the two phases caused by the motion of the walls, then
mathematically this means that the oscillations happen on top of phases with
sizes $\ell_{i}+\Delta\ell_{i}(t)$, whose energy density can be written
approximately as the first term in (4.2). The form of this term is just a
simple way of “stretching” (for positive $\Delta\ell_{i}$) or “compressing”
(for negative $\Delta\ell_{i}$) the domain profile by “gluing in” or “cutting
out” a small piece at the centre of each phase, taking advantage of the fact
that the energy density is almost exactly constant there. In order to
determine $\Delta\ell_{1}$ as a function of time, we simply impose
conservation of energy, namely that the energy change associated to the rigid
shift of the walls is exactly compensated by the energy change associated to
the oscillations of the phases:

$2\left(\mathcal{E}_{high}-\mathcal{E}_{low}\right)\Delta\ell_{1}(t)+\int_{-L/2}^{L/2}dz\Big{[}\delta\mathcal{E}_{high}(t,z)+\delta\mathcal{E}_{low}(t,z)\Big{]}=0\,.$
(4.3)

The approximation (4.2) becomes applicable soon after the merger, as
illustrated by Fig. 14.

---
Figure 14:  Oscillations of the high- and low-energy phases during the
relaxation period after the merger.

To describe correctly the perturbations of the high-energy phase it is enough
to truncate the series (4.1) at $n=1$, whereas in order to obtain a good
result for the low-energy phase we included the second mode, $n=2$. The values
that we found for the fit parameters in (4.1) are

$\begin{array}[]{c c
c}\omega_{1}/\Lambda=-0.0048,&\omega_{2}/\Lambda=0.032,&t_{0}\Lambda=800,\\\
L_{high}\Lambda=33.33,&L_{low}\Lambda=50.82,&L\Lambda=84.14,\\\
z_{0}^{high}=0,&z_{0}^{low}=\pm L/2,&a_{1}^{high}/\Lambda^{4}=0.0375,\\\
\gamma_{1}^{high}=0.227,&a_{1}^{low}/\Lambda^{4}=0.0142,&\gamma_{1}^{low}=2.346,\\\
a_{2}^{low}/\Lambda^{4}=-0.0022,&\gamma_{2}^{low}=-0.244.\end{array}\\\ $
(4.4)

###  4.2 Intermediate velocity

A simulation in this regime with $v_{max}=0.11$ is shown in Fig. 15. In this
case there is no quasi-static regime, as seen by the absence of an
intermediate plateau in the entropy density in Fig. 16. The overall increase
in the entropy is greater than in the low-velocity case, confirming the
intuition that this is a more violent collision.

Figure 15: Evolution of the energy density (top) and of the fluid velocity
(bottom) for a simulation with $v_{\mathrm{max}}=0.11$. Figure 16:  Time
evolution of the total entropy (per unit transverse area) for the evolution in
Fig. 15. $S_{2}$ is the entropy of a static configuration with two domains
(the dashed, red curve in Fig. 8). $S_{1}$ is the entropy of a static
configuration with the same total energy but a single domain.
---
Figure 17:  (Left) Energy density $\mathcal{E}/\Lambda^{4}$ and (Right)
longitudinal pressure $P_{l}/\Lambda^{4}$ as a function of time for constant
coordinate $z\Lambda=0$ for the evolution in Fig. 15.

The rest of the evolution is qualitatively similar to the low-velocity case:
there is an initial excess of pressure between the domains that slows them
down and they merge to form a domain that relaxes to equilibrium through the
same kind of oscillations, see Fig. 17.

---
Figure 18:  Oscillations of the high- and low-energy phases during the
relaxation period after the merger for the evolution in Fig. 15.

Applying (4.2) leads to a good description of the oscillations at times soon
after the collision, see Fig. 18. The parameters in this case are

$\begin{array}[]{c c
c}\omega_{1}/\Lambda=-0.0047,&\omega_{2}/\Lambda=0.032,&t_{0}\Lambda=200,\\\
L_{high}\Lambda=33.33,&L_{low}\Lambda=50.82,&L\Lambda=84.14,\\\
z_{0}^{high}=0,&z_{0}^{low}=\pm L/2,&a_{1}^{high}/\Lambda^{4}=0.1065,\\\
\gamma_{1}^{high}=2.09,&a_{1}^{low}/\Lambda^{4}=-0.0397,&\gamma_{1}^{low}=1.306,\\\
a_{2}^{low}/\Lambda^{4}=-0.0082,&\gamma_{2}^{low}=1.26\,.\end{array}\\\ $
(4.5)

Note that, up to fitting errors, the frequencies in this case and in the low-
velocity case are the same. The reason is that the final equilibrium state is
the same in both cases, and therefore the frequencies of the linear
perturbations around it are identical. The only difference between the
simulations is therefore in the amplitudes of the oscillations. The following
table collects the ratios of the coefficients for the $n=1,2$ modes,
$a_{2}^{i}/a_{1}^{i}$, for both phases for various simulations:

$\begin{array}[]{|c| c| c|}\hline\cr
v_{max}&\mathrm{low}&\mathrm{intermediate}\\\ \hline\cr\mathrm{low-
energy\,\,phase}&-0.158&0.2058\\\ \hline\cr\mathrm{high-
energy\,\,phase}&-0.043&-0.044\\\ \hline\cr\end{array}$ (4.6)

We see that they are of the same order for different simulations, and that the
ratio is smaller in the high-energy phase, consistently with our neglect of
the $n=2$ (and higher) modes in this phase.

###  4.3 High velocity

For velocities $v_{max}\geq v_{2}=0.2$ the domains collide but do not merge in
the first collision. Instead, the excited state that results from the
collision breaks apart into pieces that will subsequently collide again until
they finally merge. Fig. 19 shows two illustrative cases.

Figure 19:  Evolution of the energy density for simulations with
$v_{\mathrm{max}}=0.25$ (top) and $v_{\mathrm{max}}=0.73$ (bottom).

As shown by Fig. 20, the peak energy density right after the first collision
increases significantly with $v_{max}$, meaning that the resulting
configuration is a large deviation from an equilibrium domain. The maximum
relative deformations are $\mathcal{E}_{peak}/\mathcal{E}_{high}\simeq
1.071,1.096,1.762,3.966$, respectively, for the increasing velocities shown in
Fig. 20. As a consequence, non-linear dynamics becomes important and the
excited state breaks apart into smaller components rather than relaxing to a
single equilibrium domain as in previous sections.

Figure 20:  Snapshots of the energy density at the times at which the maximum
energy is reached for collisions with different $v_{max}$. Figure 21:
Snapshots of the energy density for the collision of Fig. 19(top) with
$v_{max}=0.25$. The blue curve corresponds to a time soon after the
fragmentation of the excited state. Note that the horizontal axis in this
figure has been shifted to show the new domain at the center of the figure and
the central fragments on the sides of it.

Figure 22:  Snapshots of the energy density for the evolution in Fig.
19(bottom) with $v_{max}=0.73$. The top (bottom) plot corresponds to a time
shortly after the second (fourth) collision. The blue curve corresponds to a
time soon after the fragmentation of the excited state.

In both cases shown in Fig. 19, the initial excited state first “emits”
perturbations that travel away from the central blob, and subsequently it
splits into two fragments. The perturbations travel away from the central
blob, colliding with each other on the other side of the box, at
$z\Lambda\simeq 42$. We will refer to this as the second collision. This
collision creates an excitation that quickly becomes a new domain smaller than
the central fragments, as illustrated in Fig. 21. Note from the horizontal
axis in this figure that the new domain is shown in the center of the figure
and the central fragments on the sides of it. At later times the two central
fragments reach the location of the new domain and a third collision, this
time a three-body one, takes place. In cases with $v_{max}\gtrsim v_{2}$, such
as the one in Fig. 19(top), this third collision results in the formation of a
single domain that then relaxes to equilibrium as described in previous
sections.

In contrast, in cases with large enough $v_{max}$, such as the one of Fig.
19(bottom), the excitations emitted by the initial excited state, as well as
the fragments to which it decays, travel faster to the opposite side of the
box and do not have time to relax to domains, as shown in Fig. 22(top). As
above, the perturbations collide in what we call the second collision.
However, in this case the third collision, i.e. the collision between the
resulting new structure and the fragments from the first excited state, is
again a high-velocity collision. As a consequence, it results in the emission
of new perturbations plus fragments that in this case travel from
$z\Lambda\simeq 42$ towards $z\Lambda=0$. These products now have time to
relax to approximate domains, as illustrated in Fig. 22(bottom). The
perturbations merge in a fourth collision at $z\Lambda=0$, and finally a fifth
collision occurs, again a three-body collision, which results in the formation
of a single domain that relaxes to equilibrium as in previous sections.

The entropy production in Fig. 23 reflects the qualitative differences
described above between the two evolutions of Fig. 19. On the left plot we can
identify three different regimes corresponding to the three different
collisions. In contrast, no such clear distinction is apparent on the right
plot. In both cases, however, the relative amount of entropy increase is much
larger than in the low- and intermediate-velocity collisions. For comparison,
the relative entropy production in the evolutions with $v_{max}=0.11$,
$v_{max}=0.25$ and $v_{max}=0.73$ are 2, 9 and 45 times larger than in the
$v_{max}=0.08$ case, respectively.

---
Figure 23:  Time evolution of the total entropy (per unit transverse area) for
the evolutions of Fig. 19 with $v_{max}=0.25$ (left) and $v_{max}=0.73$
(right). $S_{1}$ is the entropy of a static configuration with the same total
energy but a single domain. The vertical dashed lines represent the
approximate time of each collision.

##  5 Discussion

The time evolution of the spinodal instability in a theory with a first-order
phase transition typically results in the creation of phase domains that
subsequently collide with one another [20]. The velocities of the different
domains are a complicated function of the perturbation of the initial,
homogeneous, unstable state. This makes it difficult to perform a systematic
study of the physics of the collision as a function of the domain velocities.
In this paper we have overcome this difficulty by directly preparing initial
states consisting of domains moving towards each other at velocities $0\leq
v_{max}\leq 0.73$, where $v_{max}$ is the maximum fluid velocity in the
initial state. Going to higher velocities becomes challenging due to numerical
issues. Nevertheless, the above range sufficed to uncover three qualitatively
different dynamical regimes.

For low velocities the domains initially slow down, enter a period of quasi-
static evolution, and finally collide and merge into a single domain that then
relaxes to equilibrium through damped oscillations. The quasi-static period is
clearly visible as a plateau in the entropy of the system, which we computed
from the area of the dual horizon on the gravity side. For intermediate
velocities the evolution is qualitatively identical except for the fact that
no quasi-static period is present. For high velocities the domains can collide
several times before they eventually merge into a final, single domain that
then relaxes to equilibrium. Our investigations suggest that, while the
precise values of the velocities that distinguish these three regimes depend
on the size of the domains in question, the existence of these regimes is a
robust property of the collision dynamics.

## Acknowledgements

YB acknowledges support from the European Research Council Grant No.
ERC-2014-StG 639022-NewNGR and the Academy of Finland grant no. 333609. TG
acknowledges financial support from FCT/Portugal Grant No. PD/BD/135425/2017
in the framework of the Doctoral Programme IDPASC-Portugal. MSG acknowledges
financial support from the APIF program, fellowship APIF_18_19/226. JCS, DM
and MSG are also supported by grants SGR-2017-754, PID2019-105614GB-C21,
PID2019-105614GB-C22 and the “Unit of Excellence MdM 2020-2023” award to the
Institute of Cosmos Sciences (CEX2019-000918-M). MZ acknowledges financial
support provided by FCT/Portugal through the IF programme grant IF/00729/2015,
and CERN/FIS-PAR/0023/2019. The authors thankfully acknowledge the computer
resources, technical expertise and assistance provided by CENTRA/IST.
Computations were performed in part at the cluster “Baltasar-Sete-Sóis” and
supported by the H2020 ERC Consolidator Grant “Matter and strong field
gravity: New frontiers in Einstein’s theory” grant agreement No.
MaGRaTh-646597. We also thank the MareNostrum supercomputer at the BSC
(activity Id FI-2021-1-0008) for significant computational resources.

## References

  * [1] P. M. Chesler and L. G. Yaffe, “Horizon formation and far-from-equilibrium isotropization in supersymmetric Yang-Mills plasma,” Phys. Rev. Lett. 102, 211601 (2009) doi:10.1103/PhysRevLett.102.211601 [arXiv:0812.2053 [hep-th]].
  * [2] P. M. Chesler and L. G. Yaffe, “Boost invariant flow, black hole formation, and far-from-equilibrium dynamics in N = 4 supersymmetric Yang-Mills theory,” Phys. Rev. D 82, 026006 (2010) doi:10.1103/PhysRevD.82.026006 [arXiv:0906.4426 [hep-th]].
  * [3] M. P. Heller, R. A. Janik and P. Witaszczyk, “The characteristics of thermalization of boost-invariant plasma from holography,” Phys. Rev. Lett. 108, 201602 (2012) doi:10.1103/PhysRevLett.108.201602 [arXiv:1103.3452 [hep-th]].
  * [4] P. M. Chesler and L. G. Yaffe, “Holography and colliding gravitational shock waves in asymptotically AdS5 spacetime,” Phys. Rev. Lett. 106, 021601 (2011) doi:10.1103/PhysRevLett.106.021601 [arXiv:1011.3562 [hep-th]].
  * [5] M. P. Heller, D. Mateos, W. van der Schee and D. Trancanelli, “Strong Coupling Isotropization of Non-Abelian Plasmas Simplified,” Phys. Rev. Lett. 108, 191601 (2012) doi:10.1103/PhysRevLett.108.191601 [arXiv:1202.0981 [hep-th]].
  * [6] M. P. Heller, D. Mateos, W. van der Schee and M. Triana, “Holographic isotropization linearized,” JHEP 1309, 026 (2013) doi:10.1007/JHEP09(2013)026 [arXiv:1304.5172 [hep-th]].
  * [7] J. Casalderrey-Solana, M. P. Heller, D. Mateos and W. van der Schee, “From full stopping to transparency in a holographic model of heavy ion collisions,” Phys. Rev. Lett. 111, 181601 (2013) [arXiv:1305.4919 [hep-th]].
  * [8] J. Casalderrey-Solana, M. P. Heller, D. Mateos and W. van der Schee, “Longitudinal Coherence in a Holographic Model of Asymmetric Collisions,” Phys. Rev. Lett. 112, no. 22, 221602 (2014) [arXiv:1312.2956 [hep-th]].
  * [9] P. M. Chesler and L. G. Yaffe, “Holography and off-center collisions of localized shock waves,” JHEP 1510, 070 (2015) [arXiv:1501.04644 [hep-th]].
  * [10] P. M. Chesler, “Colliding shock waves and hydrodynamics in small systems,” Phys. Rev. Lett. 115, no. 24, 241602 (2015) [arXiv:1506.02209 [hep-th]].
  * [11] P. M. Chesler and W. van der Schee, “Early thermalization, hydrodynamics and energy loss in AdS/CFT,” Int. J. Mod. Phys. E 24, no. 10, 1530011 (2015) [arXiv:1501.04952 [nucl-th]].
  * [12] A. Buchel, M. P. Heller and R. C. Myers, “Equilibration rates in a strongly coupled nonconformal quark-gluon plasma,” Phys. Rev. Lett. 114 (2015) no.25, 251601 [arXiv:1503.07114 [hep-th]].
  * [13] P. M. Chesler, “How big are the smallest drops of quark-gluon plasma?,” JHEP 1603, 146 (2016) [arXiv:1601.01583 [hep-th]].
  * [14] M. Attems, J. Casalderrey-Solana, D. Mateos, D. Santos-Oliván, C. F. Sopuerta, M. Triana and M. Zilhão, “Holographic Collisions in Non-conformal Theories,” JHEP 1701, 026 (2017) [arXiv:1604.06439 [hep-th]].
  * [15] J. Casalderrey-Solana, D. Mateos, W. van der Schee and M. Triana, “Holographic heavy ion collisions with baryon charge,” JHEP 1609, 108 (2016) doi:10.1007/JHEP09(2016)108 [arXiv:1607.05273 [hep-th]].
  * [16] M. Attems, J. Casalderrey-Solana, D. Mateos, I. Papadimitriou, D. Santos-Oliván, C. F. Sopuerta, M. Triana and M. Zilhão, “Thermodynamics, transport and relaxation in non-conformal theories,” JHEP 1610, 155 (2016) [arXiv:1603.01254 [hep-th]].
  * [17] U. Gürsoy, A. Jansen and W. van der Schee, “New dynamical instability in asymptotically anti-de Sitter spacetime,” Phys. Rev. D 94, no. 6, 061901 (2016) [arXiv:1603.07724 [hep-th]].
  * [18] M. Attems, Y. Bea, J. Casalderrey-Solana, D. Mateos, M. Triana and M. Zilhao, “Phase Transitions, Inhomogeneous Horizons and Second-Order Hydrodynamics,” JHEP 1706, 129 (2017) [arXiv:1703.02948 [hep-th]].
  * [19] R. A. Janik, J. Jankowski and H. Soltanpanahi, “Real-Time dynamics and phase separation in a holographic first order phase transition,” Phys. Rev. Lett. 119, no. 26, 261601 (2017) [arXiv:1704.05387 [hep-th]].
  * [20] M. Attems, Y. Bea, J. Casalderrey-Solana, D. Mateos and M. Zilhão, “Dynamics of Phase Separation from Holography,” JHEP 01 (2020), 106 [arXiv:1905.12544 [hep-th]].
  * [21] M. Attems, J. Casalderrey-Solana, D. Mateos, D. Santos-Oliván, C. F. Sopuerta, M. Triana and M. Zilhão, “Paths to equilibrium in non-conformal collisions,” JHEP 1706, 154 (2017) [arXiv:1703.09681 [hep-th]].
  * [22] R. Rougemont, R. Critelli, J. Noronha-Hostler, J. Noronha and C. Ratti, “Dynamical versus equilibrium properties of the QCD phase transition: A holographic perspective,” Phys. Rev. D 96 (2017) no.1, 014032 [arXiv:1704.05558 [hep-ph]].
  * [23] R. Critelli, R. Rougemont and J. Noronha, “Holographic Bjorken flow of a hot and dense fluid in the vicinity of a critical point,” Phys. Rev. D 99, no. 6, 066004 (2019) doi:10.1103/PhysRevD.99.066004 [arXiv:1805.00882 [hep-th]].
  * [24] A. Buchel, “Non-conformal holographic Gauss-Bonnet hydrodynamics,” JHEP 1803, 037 (2018) [arXiv:1801.06165 [hep-th]].
  * [25] A. Czajka, K. Dasgupta, C. Gale, S. Jeon, A. Misra, M. Richard and K. Sil, “Bulk Viscosity at Extreme Limits: From Kinetic Theory to Strings,” JHEP 07 (2019), 145 [arXiv:1807.04713 [hep-th]].
  * [26] A. Czajka, K. Dasgupta, C. Gale, S. Jeon, A. Misra, M. Richard and K. Sil, “On bulk viscosity at weak and strong ’t Hooft couplings,” Mod. Phys. Lett. A 35 (2020) no.27, 2030012 [arXiv:1807.07950 [hep-th]].
  * [27] M. Attems, Y. Bea, J. Casalderrey-Solana, D. Mateos, M. Triana and M. Zilhão, “Holographic Collisions across a Phase Transition,” Phys. Rev. Lett. 121, no. 26, 261601 (2018) doi:10.1103/PhysRevLett.121.261601 [arXiv:1807.05175 [hep-th]].
  * [28] H. Bantilan, P. Figueras and D. Mateos, “Real-time Dynamics of Plasma Balls from Holography,” Phys. Rev. Lett. 124 (2020) no.19, 191601 [arXiv:2001.05476 [hep-th]].
  * [29] Y. Bea, O. J. C. Dias, T. Giannakopoulos, D. Mateos, M. Sanchez-Garitaonandia, J. E. Santos and M. Zilhao, “Crossing a large-$N$ phase transition at finite volume,” JHEP 02 (2021), 061 [arXiv:2007.06467 [hep-th]].
  * [30] J. Casalderrey-Solana, C. Ecker, D. Mateos and W. Van Der Schee, “Strong-coupling dynamics and entanglement in de Sitter space,” JHEP 03 (2021), 181 [arXiv:2011.08194 [hep-th]].
  * [31] Y. Bea, J. Casalderrey-Solana, T. Giannakopoulos, D. Mateos, M. Sanchez-Garitaonandia and M. Zilhão, “Bubble Wall Velocity from Holography,” [arXiv:2104.05708 [hep-th]].
  * [32] C. Ecker, W. van der Schee, D. Mateos and J. Casalderrey-Solana, “Holographic Evolution with Dynamical Boundary Gravity,” [arXiv:2109.10355 [hep-th]].
  * [33] J. Casalderrey-Solana, H. Liu, D. Mateos, K. Rajagopal and U. A. Wiedemann, “Gauge/String Duality, Hot QCD and Heavy Ion Collisions,” book:Gauge/String Duality, Hot QCD and Heavy Ion Collisions. Cambridge, UK: Cambridge University Press, 2014 [arXiv:1101.0618 [hep-th]].
  * [34] M. A. Stephanov, K. Rajagopal and E. V. Shuryak, “Signatures of the tricritical point in QCD,” Phys. Rev. Lett. 81, 4816 (1998) [hep-ph/9806219].
  * [35] M. A. Stephanov, K. Rajagopal and E. V. Shuryak, “Event-by-event fluctuations in heavy ion collisions and the QCD critical point,” Phys. Rev. D 60, 114028 (1999) [hep-ph/9903292].
  * [36] M. Stephanov and Y. Yin, “Hydrodynamics with parametric slowing down and fluctuations near the critical point,” Phys. Rev. D 98 (2018) no.3, 036006 [‘arXiv:1712.10305 [nucl-th]].
  * [37] Y. Bea and D. Mateos, “Heating up Exotic RG Flows with Holography,” JHEP 1808 (2018) 034, arXiv:1805.01806 [hep-th].