CPHT-RR051.062024

HIP-2024-17/TH

# Spinodal slowing down and scaling in a holographic model

Alessio Caddeo<EMAIL_ADDRESS>Department of Physics and
Instituto de Ciencias y Tecnologías Espaciales de Asturias (ICTEA),
Universidad de Oviedo, c/ Leopoldo Calvo Sotelo 18, ES-33007 Oviedo, Spain
Oscar Henriksson<EMAIL_ADDRESS>Faculty of Science and Engineering,
Åbo Akademi University, Henrikinkatu 2, 20500 Turku, Finland  Helsinki
Institute of Physics, P.O. Box 64, 00014 University of Helsinki, Finland
Carlos Hoyos<EMAIL_ADDRESS>Department of Physics and
Instituto de Ciencias y Tecnologías Espaciales de Asturias (ICTEA),
Universidad de Oviedo, c/ Leopoldo Calvo Sotelo 18, ES-33007 Oviedo, Spain
Mikel Sanchez-Garitaonandia<EMAIL_ADDRESS>CPHT, CNRS, École
polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau, France

###### Abstract

The dynamics of first-order phase transitions in strongly coupled systems are
relevant in a variety of systems, from heavy ion collisions to the early
universe. Holographic theories can be used to model these systems, with
fluctuations usually suppressed. In this case the system can come close to a
spinodal point where theory and experiments indicate that the the behaviour
should be similar to a critical point of a second-order phase transition. We
study this question using a simple holographic model and confirm that there is
critical slowing down and scaling behaviour close to the spinodal point, with
precise quantitative estimates. In addition, we determine the start of the
scaling regime for the breakdown of quasistatic evolution when the temperature
of a thermal bath is slowly decreased across the transition. We also extend
the analysis to the dynamics of second-order phase transitions and strong
crossovers.

###### Contents

  1. I Introduction
  2. II The holographic model and its equilibrium behavior
    1. II.1 Breaking of conformal invariance
    2. II.2 Holographic renormalization and alternative quantization
    3. II.3 Effective potential and multitrace boundary conditions
    4. II.4 Static solutions and phase diagram
  3. III Relaxation time and correlation length
    1. III.1 Relaxation time and critical slowing down
    2. III.2 Correlation length and dynamical exponent
  4. IV Breakdown of quasistatic time evolution and scaling
    1. IV.1 Quasistatic evolution
    2. IV.2 Critical and spinodal scaling
    3. IV.3 Full time evolution
      1. IV.3.1 Scaling in first and second-order phase transitions
      2. IV.3.2 Scaling in crossovers and weak first-order phase transitions
  5. V Discussion and outlook
  6. A Stationary points of the effective potential
  7. B Analytic estimate of the decay rate and correlation length
    1. B.1 Decay rate
    2. B.2 Correlation length
  8. C Quasinormal modes
  9. D Numerical methods for the time evolution

##  I Introduction

Matter going through a first-order phase transition as it cools down is a
commonly observed occurrence in daily life, as when we turn water into ice
cubes for our drinks. At much smaller subatomic scales, extremely hot matter
produced in heavy ion collisions may undergo a transition as a droplet of
quark-gluon plasma expands and cools down, turning into hadronic matter [1,
2]. In the other extreme, at cosmological scales, an electroweak or dark
matter phase transition in the early universe may leave an observable imprint
in the gravitational wave spectrum [3, 4]. If observed, the former would
validate phenomenological predictions for the QCD phase diagram, while the
latter would be a clear signal of physics beyond the Standard Model.

Holographic models have been used to describe and estimate the effects of such
phase transitions, for both heavy ion collisions [5, 6, 7] and cosmology [8,
9, 10, 11, 12, 13, 14, 15]. Such models use a classical gravitational theory
to compute non-trivial quantities of a dual strongly coupled field theory. In
particular, one can describe the out-of-equilibrium evolution of the field
theory as it goes through a phase transition by solving a few classical
differential equations. This makes the holographic duality a powerful tool to
study various aspects of such transitions, including the nucleation of bubbles
and the properties of domain walls [16, 17, 18, 19, 20, 21, 22, 23, 24, 25,
26], the formation and evolution of mixed phases [27, 28, 29, 30, 31, 32, 33,
34], and quenches [35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47].

In this paper, we continue this program by using a simple holographic model to
study slow evolution through first-order phase transitions, comparing the
results to second-order transitions and crossovers.

Holographic duality is typically used in the classical gravity limit, which,
in the dual field theory, is a large-$N$ or mean field limit. In the canonical
case $N$ is the rank of the gauge group [48], but more generally we can think
of it as a measure of the number of degrees of freedom. Restricting to a
classical gravitational description is therefore akin to completely removing
fluctuations around the mean field theory. Even if $N$ is not strictly
infinite, in order for the holographic approach to be tractable, fluctuations
must be strongly suppressed. This can have a significant effect on the
evolution of the system through a phase transition.

Consider the case of a system with a first-order phase transition. The order
parameter as a function of the temperature would typically behave as shown in
Figure 1. Within some interval of temperatures around the critical value, the
high- and low-temperature phases coexist; the one with lower free energy is
the true stable phase, while the other is metastable. Precisely at the
critical temperature the free energy of the two phases is the same; the phase
which is stable at larger temperatures is metastable at smaller temperatures,
and vice versa. Each of the stable phases ends at a _spinodal point_ , marked
as green dots in the figure, where they merge with an unstable third
branch.111Note that this is a mean field, or coarse grained, depiction of the
phase diagram, since the true free energy must be a convex function of
thermodynamical variables. The true equilibrium states include mixtures of the
two phases.

Figure 1: Order parameter as a function of the temperature for the case in the
first row of Table 1. The blue and red lines are the (meta)stable low and high
temperature phases respectively. The dashed purple line is the unstable phase.
The dashed black vertical line marks the critical temperature, and the green
dots the position of the spinodal points where metastable and unstable phases
merge.

We will study the case where the system of interest is put in contact with a
large thermal bath that determines the temperature, which we then change with
time. If the change is slow enough, the system could in principle remain
arbitrarily close to equilibrium, undergoing a quasistatic evolution of
equilibrium quantities.

Starting at high temperatures and cooling down, the system will typically
undercool, reaching temperatures below the critical one before fluctuations
drive it to the low temperature phase. Usually, the fluctuations take the form
of bubbles of the low temperature phase nucleating inside the high temperature
phase, which expand due to their higher inner pressure (lower free energy). If
the nucleating barrier is very low, the size of the bubble can become of the
same order of the width of the interface and in this case one should also take
into account the formation of ramified domains of the low temperature phase
[49]. As argued in [20, 11], bubble nucleation is exponentially suppressed
with $N$ in holographic models, so, for a large enough value of $N$, the time
evolution can be slow enough as to remain quasistatic, but fast enough to
completely neglect nucleation before reaching the vicinity of the spinodal
point. In other words, the large-$N$ limit can uncover the physics of the
spinodal points, which would otherwise be hidden by the nucleation process.

Although not widely discussed, spinodal points share some features with
critical points of second-order phase transitions, for example the scaling
behaviour of thermodynamic variables and correlators, as well as the critical
slowing down of perturbations [50, 51, 52]. These features may be relevant
even in the presence of thermal fluctuations and disorder, as long as they are
not too strong [53, 54, 55, 56], and they have been observed experimentally in
some strongly correlated materials [57, 58, 59, 60, 61, 62, 63, 64]. The
critical slowing down means that the time it takes for an out-of-equilibrium
perturbation to decay diverges as one approaches the critical or spinodal
points. Thus, even if the evolution is very slow, eventually the system is
driven out of equilibrium before reaching the critical or spinodal
temperature.

In this paper we provide a simple theoretical realization of a strongly
correlated system where the properties of transitions through spinodal and
critical points can be studied, confirming the above picture. By tuning the
parameters of our holographic model, we can realize phase transitions of first
and second order, as well as smooth crossovers. In the case of first-order
transitions, we show that the system can never reach the spinodal point; no
matter how slow the time evolution is, the system will fall out of equilibrium
before reaching it, even in the absence of bubble nucleation. We will also
show evidence for the scaling behaviour at the spinodal point in the
holographic model and compare with scaling at critical points and strong
crossovers. Although our results are obtained for a particular class of
models, we expect the qualitative features to be valid more generally for
phase transitions in strongly coupled theories with a holographic dual.

The rest of the paper is organized as follows: In section II we introduce the
holographic model we work with, and discuss its equilibrium configurations and
phase diagram. In section III, we introduce small perturbations around
equilibrium and extract the relaxation time and correlation length. We argue
analytically, and show numerically, that the model exhibits critical slowing
down at the spinodal points. In section IV, we determine the breakdown of
quasistatic evolution and the scaling close to critical and spinodal points,
both through analytic approximations and numerically. We conclude and discuss
possible future directions in section V. Several appendices go into more
detail regarding the phase diagram, the computation of decay rates and
correlation lengths, quasinormal modes, and the numerical methods used.

##  II The holographic model and its equilibrium behavior

We will work with a bottom-up gravitational model dual to a $2+1$-dimensional
strongly coupled CFT with a large number $\sim N_{c}^{2}$ of degrees of
freedom. Conformal invariance will be broken in a subsector involving $\sim
N_{c}$ degrees of freedom, akin to fundamental fields spanning the flavor
sector of a gauge theory with a gauge group of rank $N_{c}$. We will work in
the quenched approximation where the effect of this ‘flavor’ subsector over
the rest of the theory is neglected.

The gravity dual to this setup consists of Einstein gravity with a
cosmological constant plus a scalar field. The classical action is

$S=\frac{1}{2\kappa_{4}^{2}}\int
d^{4}x\sqrt{|g|}\left[R+\frac{6}{L^{2}}-\frac{1}{N_{c}}\left(g^{MN}\partial_{M}\phi\partial_{N}\phi+2P(\phi)\right)\right]\
,$ (1)

where $g_{MN}$ is the metric, $R$ the associated scalar curvature and $\phi$
the scalar field with a potential $P(\phi)$. The metric $g_{MN}$ corresponds
to the ‘glue’ sector with $\sim N_{c}^{2}$ degrees of freedom, while the
scalar $\phi$ will capture the physics of the flavor sector. Following the
quenched approximation, we will treat the scalar as a probe, neglecting the
backreaction over the metric. The metric at zero temperature will then be AdS4
with radius $L$.

When turning on the temperature the glue sector will act as a thermal bath for
the flavor sector. Since there are no phase transitions in the glue sector,
for a time evolution slow enough (compared to the scale set by the
temperature), we can assume that the evolution is quasistatic, i.e. the system
is at thermal equilibrium at each instant of time. We will then fix the metric
to be AdS4 black brane with a position of the horizon $z_{h}$ that may depend
on time

$ds^{2}=\frac{L^{2}}{z^{2}}\left(f(z)dt^{2}+\frac{dz^{2}}{f(z)}+dx^{i}dx^{i}\right)\
,\quad\qquad f(z)=1-\frac{z^{3}}{z_{h}^{3}}\ .$ (2)

The temperature in the dual field theory equals the Hawking temperature of the
black brane geometry,

$T=\frac{3}{4\pi z_{h}}\ .$ (3)

If $z_{h}$ depends on time this metric is not a solution of Einstein’s
equations, but can serve as a good approximation as long as the the time
derivatives are small enough.

The flavor sector will be at the temperature set by the bath. Normally it
should also follow a quasistatic evolution as the temperature is changed, but
when it comes close to a phase transition we expect this will no longer be
true. In order to elucidate the evolution through the phase transition we will
solve the equations of motion of the dual scalar field in the background
metric with a time-dependent horizon. We will assume homogeneous
configurations in the spatial directions, so there will be dependence only in
time and the holographic radial direction. However, before this can be done,
it will be necessary to set up the model such that it shows a phase
transition. This can only happen if conformal invariance is broken in some
way.

###  II.1 Breaking of conformal invariance

We will now specify the potential for the scalar field. For small values of
the amplitude it typically starts with a quadratic term

$P(\phi)=\frac{m^{2}}{2}\phi^{2}+\cdots\ .$ (4)

The asymptotic boundary expansion of the scalar field in the AdS4 metric is
determined by the value of the mass in units of the AdS radius. One finds two
independent solutions:

$\phi(z)\sim\phi_{-}z^{\frac{3}{2}-\nu}+\phi_{+}z^{\frac{3}{2}+\nu}\
,\qquad\nu=\sqrt{m^{2}L^{2}+\frac{9}{4}}\ ,$ (5)

with $\phi_{-}$ and $\phi_{+}$ the coefficients of the leading and subleading
solutions. For $\nu\geq 3/2$ the leading solution maps to a source, or
coupling, for the dual operator and the subleading solution determines the
expectation value. However, for $3/2>\nu>0$ one may choose an ‘alternative
quantization’ where the roles of the two solutions are reversed [65, 66]. In
the case of ordinary quantization the conformal dimension of the dual operator
is $\Delta_{+}=\frac{3}{2}+\nu$, while for alternative quantization it is
$\Delta_{-}=\frac{3}{2}-\nu$.

The simplest way to break conformal invariance is to turn on a source for a
relevant operator, which will amount to solving the equations of motion for
the scalar field fixing $\phi_{-}$ or $\phi_{+}$. It can be shown that this is
not enough to induce a phase transition with only a quadratic potential. A
possible solution is to design a more complicated potential $P(\phi)$, for
instance by adding higher powers of the scalar field with appropriate
coefficients as in [67]. We will instead follow the approach of [20], where
the bulk potential $P(\phi)$ remains quadratic, but a non-trivial effective
potential is obtained by using the alternative quantization for the scalar
field and tuning its boundary conditions. This corresponds to deforming the
dual field theory by various powers of the operator $\mathcal{O}$ dual to
$\phi$: $S_{CFT}\to S_{CFT}+g_{n}\mathcal{O}^{n}$ [68]. In a large-$N$ theory
such a deformation simply results in the addition of a similar term to the
effective potential depending on the expectation value of the operator,
$V(\langle\mathcal{O}\rangle)\to
V(\langle\mathcal{O}\rangle)+g_{n}\langle\mathcal{O}\rangle^{n}$ [69].

We want to construct an effective potential which is quartic in
$\langle\mathcal{O}\rangle$, as this will let us realize first- and second-
order thermal phase transitions, as well as smooth crossovers. In order not to
have irrelevant deformations of the CFT, this means that the conformal
dimension of the operator should be $\Delta\leq 3/4$. We find it convenient to
saturate this condition, meaning that $m^{2}L^{2}=-27/16$. Implementing both
alternative quantization and the boundary conditions requires adding
additional boundary terms to the gravitational on-shell action, after it has
been properly renormalized, as we now review.

###  II.2 Holographic renormalization and alternative quantization

Since we are working in the probe approximation, the Einstein and scalar
actions are renormalized independently. For us only the scalar action is
relevant. Following the usual procedure of holographic renormalization, we
introduce a cutoff in the holographic radial coordinate $z=z_{UV}$, and add a
boundary counterterm of the form

$S_{c.t.}=-\frac{3}{8L\kappa_{4}^{2}N_{c}}\int
d^{3}x\sqrt{|h|}\phi^{2}\Big{|}_{z=z_{UV}}\ ,$ (6)

where $h_{\mu\nu}=g_{\mu\nu}\Big{|}_{z=z_{UV}}$ is the induced metric on a
radial slice.

The variation of the regularized scalar on-shell action is

$\delta S_{\phi}=\delta S_{\text{on-shell}}+\delta
S_{c.t.}=\frac{1}{\kappa_{4}^{2}N_{c}}\int
d^{3}x\left(\sqrt{|g|}\,g^{zz}\partial_{z}\phi-\frac{3}{4L}\sqrt{|h|}\phi\right)\delta\phi\Big{|}_{z=z_{UV}}\
.$ (7)

Taking the $z_{UV}\to 0$ limit, we obtain the variation of the generating
functional of the dual field theory in ordinary quantization

$\delta{\cal W}_{\Delta_{+}}[\phi_{-}]=\frac{3L^{2}}{2\kappa_{4}^{2}N_{c}}\int
d^{3}x\,\phi_{+}\delta\phi_{-}\ .$ (8)

In order to move to alternative quantization we add the boundary term [69]

$S_{a.q.}=-\frac{3L^{2}}{2\kappa_{4}^{2}N_{c}}\int d^{3}x\,\phi_{+}\phi_{-}\
.$ (9)

Combining this with the variation of the renormalized on-shell action gives
the variation of the generating functional of the dual field theory when the
operator has dimension $\Delta_{-}=3/4$:

$\delta{\cal W}_{\Delta_{-}}[\phi_{+}]=\delta S_{\phi}\Big{|}_{z_{UV}\to
0}+\delta S_{a.q.}=-\frac{3L^{2}}{2\kappa_{4}^{2}N_{c}}\int
d^{3}x\,\phi_{-}\delta\phi_{+}\ .$ (10)

The generating functional itself (rather than the variation) can also be
explicitly computed in this case. First, the on-shell action for the scalar is

$S_{\phi}=\frac{1}{2\kappa_{4}^{2}N_{c}}\int
d^{3}x\,\sqrt{|g|}g^{zz}\partial_{z}\phi\,\phi\,\Big{|}_{z=z_{UV}}^{z=z_{h}}\
.$ (11)

For regular solutions the contribution at the horizon vanishes. Then, adding
all the boundary terms, the generating functional is

${\cal W}_{\Delta_{-}}[\phi_{+}]=-\frac{3L^{2}}{4\kappa_{4}^{2}N_{c}}\int
d^{3}x\,\phi_{+}\phi_{-}\ .$ (12)

###  II.3 Effective potential and multitrace boundary conditions

The effective potential can be obtained from the generating functional via a
Legendre transform that trades the dependence on the source by a dependence on
the expectation value of the operator. First, we split the coefficient of the
subleading term $\phi_{+}$ into the coupling breaking conformal invariance,
$\Lambda$, plus an external source, $J$, i.e. $\phi_{+}=-\Lambda+J$. The
generating functional is

${\cal W}_{\Delta_{-}}[\Lambda,J]=-\frac{3L^{2}}{4\kappa_{4}^{2}N_{c}}\int
d^{3}x\,(J-\Lambda)\phi_{-}[J]\ .$ (13)

The expectation value of $\mathcal{O}$ is obtained by varying the generating
functional with respect to the source:

$\left\langle\cal O\right\rangle=\frac{\delta{\cal W}_{\Delta_{-}}}{\delta
J}=-\frac{3L^{2}}{2\kappa_{4}^{2}N_{c}}\phi_{-}\equiv\frac{L^{2}}{2\kappa_{4}^{2}N_{c}}\psi\
.$ (14)

The Legendre transform of the generating functional results in the effective
action

$\Gamma_{\Delta_{-}}[\psi]={\cal W}_{\Delta_{-}}[\Lambda,J]-J\left\langle\cal
O\right\rangle=-\frac{L^{2}}{4\kappa_{4}^{2}N_{c}}\int
d^{3}x\,(J[\psi]+\Lambda)\psi\ .$ (15)

Using that $J=\phi_{+}+\Lambda$, up to an overall factor the effective
potential is

$V(\psi)=\frac{1}{2}\phi_{+}[\psi]\psi+\Lambda\psi\ .$ (16)

For a quadratic scalar action like ours, the equations of motion for the
scalar field are linear, implying that $\phi_{+}$ will be linear in $\psi$.
The effective potential cannot have several minima in this case, so there
cannot be a phase transition.

We can modify the effective potential by adding additional boundary terms,
that we write using a function $W(\psi)$ whose meaning will be clear
momentarily:

$S_{W}=\frac{L^{2}}{2\kappa_{4}^{2}N_{c}}\int
d^{3}x\left(\psi\,\partial_{\psi}W(\psi)-W(\psi)\right)\ .$ (17)

With the new terms, the variation of the renormalized on-shell action in
alternative quantization becomes

$\delta{\cal W}_{\Delta_{-}}[\phi_{+}]+\delta
S_{W}=\frac{L^{2}}{2\kappa_{4}^{2}N_{c}}\int
d^{3}x\,\psi\,\delta\left(\phi_{+}+\partial_{\psi}W(\psi)\right)\equiv\frac{L^{2}}{2\kappa_{4}^{2}N_{c}}\int
d^{3}x\,\psi\,\delta J\ .$ (18)

We see that the new terms modify how the source should be identified, thus
implementing the multitrace boundary conditions [68, 69] where

$J=\phi_{+}+\Lambda+\partial_{\psi}W(\psi)\ .$ (19)

After performing the Legendre transform to get the effective action, one finds
that the effective potential is now modified to

$V(\psi)=\frac{1}{2}\phi_{+}[\psi]\psi+\Lambda\psi+W(\psi)\ .$ (20)

Note that for a linear relation between the coefficients of leading and
subleading terms $\phi_{+}[\psi]\propto\psi$ 222This will be true only in the
static case, otherwise $\phi_{+}$ also depends on time derivatives of $\psi$.,
the source is equal to the first derivative of the potential

$J=\partial_{\psi}V(\psi)\ .$ (21)

Thus, configurations with no source coincide with the critical points of the
effective potential. We will have to impose the condition $J=0$ in order to
stay within a theory with a fixed coupling $\Lambda$.

Having derived the effective potential (20), we will pick the multitrace
deformation to be a quartic polynomial

$W(\psi)=\frac{a}{2}\psi^{2}+\frac{b}{3}\psi^{3}+\frac{c}{4}\psi^{4}\ .$ (22)

We can then tune the parameters $a,b,c,\Lambda$ such that, for example, the
effective potential has several critical points and the theory exhibits a
phase transition as the temperature changes.

###  II.4 Static solutions and phase diagram

To fully specify the effective potential, we now determine the function
$\phi_{+}[\psi]$. This is done by finding the static solutions for the scalar
field $\phi=\phi(z)$ in the background AdS4 black brane (2). The solutions
should satisfy the no-source boundary condition and be regular at the black
brane horizon.

From the scalar action (1) with a quadratic potential, the equations of motion
are

$\partial_{z}\left(\frac{f(z)}{z^{2}}\partial_{z}\phi\right)-\frac{m^{2}L^{2}}{z^{4}}\phi=0\
.$ (23)

Taking into account that $m^{2}L^{2}=-27/16$, there are two independent
solutions that we identify according to the asymptotic expansion at the AdS
boundary

$\phi(z)=\phi_{-}z_{h}^{3/4}\eta_{3/4}\left(\frac{z}{z_{h}}\right)+\phi_{+}z_{h}^{9/4}\eta_{9/4}\left(\frac{z}{z_{h}}\right)\
,$ (24)

where

$\eta_{3/4}(Z)=Z^{3/4}\,{}_{2}F_{1}\left(\frac{1}{4},\frac{1}{4},\frac{1}{2},Z^{3}\right)\
,\qquad\eta_{9/4}(Z)=Z^{9/4}\,{}_{2}F_{1}\left(\frac{3}{4},\frac{3}{4},\frac{3}{2},Z^{3}\right)\
.$ (25)

Expanding close to the horizon, we find

$\phi(z)\sim C_{h}z_{h}^{-9/4}\left(\phi_{+}+\gamma
z_{h}^{-3/2}\phi_{-}\right)\log(z_{h}-z)+\cdots\ ,$ (26)

where

$C_{h}=-\frac{4\pi^{3/2}}{\Gamma\left(\frac{1}{4}\right)^{2}}\
,\qquad\gamma=\frac{2\Gamma\left(\frac{3}{4}\right)^{2}}{\Gamma\left(\frac{1}{4}\right)^{2}}\
.$ (27)

Regularity of the solution requires the coefficient of the logarithmic term to
vanish. The regular solution is the linear combination

$\eta_{R}(Z)=\eta_{3/4}(Z)-\gamma\,\eta_{9/4}(Z)\ .$ (28)

This fixes the relation between $\phi_{+}$ and $\phi_{-}$ in equilibrium
solutions

$\phi_{+}=-\gamma z_{h}^{-3/2}\phi_{-}=\frac{\gamma}{3}z_{h}^{-3/2}\psi\ .$
(29)

Combining this result with the expression for the temperature (3) and the
effective potential (20), (22) we arrive at

$V(\psi)=\Lambda\psi+\frac{1}{2}a_{T}\psi^{2}+\frac{b}{3}\psi^{3}+\frac{c}{4}\psi^{4}\
,$ (30)

where the total coefficient of the quadratic term is

$a_{T}=a+\frac{\gamma}{3}\left(\frac{4\pi T}{3}\right)^{3/2}\ .$ (31)

Finally, the value of $\psi$ is fixed by the no-source condition
$J=\partial_{\psi}V(\psi)=0$:

$\Lambda+a_{T}\psi+b\psi^{2}+c\psi^{3}=0\ .$ (32)

The real solutions of (32) determine the stationary points of the effective
potential (30) and thus the equilibrium states of the dual field theory.
Minima correspond to stable or metastable equilibrium states, while a maximum
corresponds to an unstable equilibrium state in the dual field theory. The
theory features a first-order phase transition if the potential has a single
minimum at high and low temperatures and three stationary points in some
intermediate range $T_{1}>T>T_{2}$. In appendix A we analyze in detail the
conditions the parameters in the effective potential need to satisfy in order
to have a first-order transition, a second-order transition, or a strong
crossover. In particular, we define a parameter $A$ which distinguishes
between theories with first-order transitions ($A>0$), second-order
transitions ($A=0$), and crossovers ($A<0$). By tuning $A$ we can then study
the interplay between the scaling associated to the second-order transition
and to the first-order spinodal point, as well as their effect on strong
crossovers near the second-order critical point.

We will focus on the set of parameters shown in Table 1, for which we present
numerical results in section IV.3. In Figure 2, we display the order parameter
as a function of temperature for some of these parameters. In Table 1 and
elsewhere, we write the parameters in units of a reference mass scale $M$ that
can be chosen arbitrarily as part of the definition of the theory.

Order | $A$ | $\Lambda M^{-9/4}$ | $aM^{-3/2}$ | $bM^{-3/4}$ |  $c$ | $T_{c}M^{-1}$ | $T_{1}/T_{c}$ | $T_{2}/T_{c}$ | Fig. 2
---|---|---|---|---|---|---|---|---|---
1st | $0.2014$ | $-1$ | $-1$ | $-5$ | $1$ | $4.934$ | $1.046$ | $0.8135$ | -
2nd | 0 | $-4.630$ | $0.6765$ | $-5$ | $1$ | $5.162$ | - | - | black
$\infty$ | $-5\times 10^{-3}$ | $-4.699$ | $0.7181$ | $-5$ | $1$ | $\sim 5.16$ | - | - | -
$\infty$ | $-5\times 10^{-4}$ | $-4.637$ | $0.6806$ | $-5$ | $1$ | $\sim 5.16$ | - | - | blue
$\infty$ | $-5\times 10^{-5}$ | $-4.630$ | $0.6768$ | $-5$ | $1$ | $\sim 5.16$ | - | - | purple
1st | $+5\times 10^{-3}$ | $-4.557$ | $0.6348$ | $-5$ | $1$ | $5.1609$ | $1+2\cdot 10^{-4}$ | $1-2\cdot 10^{-4}$ | -
1st | $+5\times 10^{-4}$ | $-4.623$ | $0.6723$ | $-5$ | $1$ | $5.1618$ | $1+5\cdot 10^{-6}$ | $1-6\cdot 10^{-6}$ | red
1st | $+5\times 10^{-5}$ | $-4.629$ | $0.6761$ | $-5$ | $1$ | 5.1618 | $1+2\cdot 10^{-7}$ | $1-2\cdot 10^{-7}$ | orange
Table 1: Parameters of the effective potential and critical and spinodal
temperatures for the phase transitions and crossovers studied in this paper.
Some of these entries are also displayed in Figure 2. For the crossovers
$T_{c}$ refers to the approximate temperature where the expectation value of
the scalar operator changes rapidly. Figure 2: Equilibrium expectation value
of the scalar operator as a function of the temperature. The curves represent:
first-order phase transition with $A=5\cdot 10^{-4}$ (red), $A=5\cdot 10^{-5}$
(orange); second-order phase transition (black); crossover with $A=-5\cdot
10^{-4}$ (blue), and crossover with $A=-5\cdot 10^{-5}$ (purple). Here,
$T_{\text{snd}}\approx 5.162M$ is the critical temperature of the second-order
phase transition.

##  III Relaxation time and correlation length

In this section, we fix the temperature and study small perturbations around
thermal equilibrium, in order to show that critical slowing down is a feature
of both spinodal points and second-order critical points in our model. We will
also obtain the dynamical exponent $z$ relating the correlation length $\xi$
with the relaxation time $\tau$, $\xi\sim\tau^{1/z}$ by considering
fluctuations with non-zero spatial momentum.

Starting from the action (1) and the background metric (2), we first put the
metric in Eddington-Finkelstein form with dimensionless coordinates. In order
to do this, we rescale the coordinates by the horizon position,

$z=z_{h}Z\ ,\quad\quad t=z_{h}X^{0}\ ,\quad\quad x^{i}=z_{h}X^{i}\ .$ (33)

After this, we perform the change of variables

$X^{0}=v+r(Z)\ ,\quad\quad r^{\prime}(Z)=\frac{1}{F(Z)}\ ,\quad\quad
F(Z)=1-Z^{3}\ .$ (34)

The solution for $r$ is

$r(Z)=r_{0}+\frac{1}{6}\left[\log\left(Z^{2}+Z+1\right)-2\log(1-Z)+2\sqrt{3}\tan^{-1}\left(\frac{2Z+1}{\sqrt{3}}\right)\right]\
.$ (35)

The resulting metric is

$ds^{2}=\frac{L^{2}}{Z^{2}}\left(-F(Z)dv^{2}-2dvdZ+dX^{i}dX^{i}\right)\ .$
(36)

The equation of motion we need to solve becomes

$0=\partial_{Z}^{2}\phi+\left(\frac{F^{\prime}(Z)}{F(Z)}-\frac{2}{Z}\right)\partial_{Z}\phi+\frac{27/16}{Z^{2}F(Z)}\phi+\frac{2}{ZF(Z)}\partial_{v}\phi-\frac{2}{F(Z)}\partial_{Z}\partial_{v}\phi+\frac{1}{F(Z)}\partial_{i}^{2}\phi\
.$ (37)

According to the discussion in section II.3, we need to solve (37) imposing
regularity at the horizon $Z=1$ and the no-source, multitrace, boundary
condition $J=0$, with $J$ given by (19).

Consider the regular static solution for the scalar field $\eta_{R}(Z)$ found
in (28). The equilibrium solution satisfying regularity is then

$\phi_{0}(Z)=-\frac{\overline{\psi}_{0}}{3}\eta_{R}(Z),\ \
\psi_{0}=\left(\frac{4\pi T}{3}\right)^{3/4}\overline{\psi}_{0}\ ,$ (38)

where $\overline{\psi}_{0}$ is a constant that should be chosen to satisfy the
no-source condition $J=\partial_{\psi}V(\psi_{0})=0$.

We now introduce a small perturbation around this equilibrium solution:

$\phi(v,Z)=\phi_{0}(Z)+\delta\phi(v,Z,X)\ .$ (39)

Inserting this in (37), the equation of motion for the perturbation can be
written as

$0=\partial_{Z}\left(\frac{F(Z)}{Z^{2}}\partial_{Z}\delta\phi\right)+\frac{27/16}{Z^{4}}\delta\phi+\frac{2}{Z^{3}}\partial_{v}\delta\phi-\frac{2}{Z^{2}}\partial_{Z}\partial_{v}\delta\phi+\frac{1}{Z^{2}}\partial_{i}^{2}\delta\phi\
.$ (40)

We can identify the leading and subleading coefficients in the boundary
expansion of the perturbation,

$\delta\phi(v,Z,X)\sim\delta\phi_{-}(v,X)Z^{3/4}+\delta\phi_{+}(v,X)Z^{9/4}+\cdots\,.\qquad\delta\phi_{-}=-\frac{1}{3}\delta\psi\
.$ (41)

The no-source boundary condition expanded to linear order determines the
boundary condition for the perturbation:

$\delta J=\delta\phi_{+}(v,X)+\partial_{\psi}^{2}W(\psi_{0})\delta\psi(v,X)=0\
.$ (42)

###  III.1 Relaxation time and critical slowing down

We can restrict the analysis to homogeneous perturbations
$\delta\phi(v,Z,X)=\delta\phi(v,Z)$. The no-source condition (42) and the
equation of motion (40) are consistent with an exponential decay of the
perturbation:

$\delta\phi(Z,v)=e^{-\Gamma v}\delta\phi_{0}(Z)\
,\qquad\delta\psi(v)=e^{-\Gamma v}\delta\psi_{0}\ .$ (43)

In appendix B, we obtain the following analytical estimate of the
dimensionless decay rate $\Gamma$, under the assumption that it is small:

$\Gamma\approx\frac{9}{2(1-2I_{1})}\left(\frac{4\pi
T}{3}\right)^{-3/2}V^{\prime\prime}(\psi_{0})\ ,$ (44)

where $I_{1}\approx 0.141115$. Restoring units, the relaxation time is the
inverse of the decay rate

$\tau^{-1}=\Gamma_{T}=\frac{\Gamma}{z_{h}}=\left(\frac{4\pi
T}{3}\right)\Gamma\approx\frac{9}{2(1-2I_{1})}\left(\frac{4\pi
T}{3}\right)^{-1/2}V^{\prime\prime}(\psi_{0})\ .$ (45)

If we approach the critical point of a second-order phase transition or the
spinodal point of a first-order phase transition, the curvature of the
effective potential vanishes, $V^{\prime\prime}(\psi_{0})\to 0$. Then, we have
now shown that the decay rate also vanishes $\Gamma_{T}\to 0$, so that
perturbations take longer and longer times to decay. This corresponds to the
anticipated critical slowing down.

The decay rate $\Gamma_{T}$ also corresponds to the purely imaginary frequency
$\omega=-i\Gamma_{T}$ of the lowest quasinormal mode (QNM) in the spectrum of
perturbations of the scalar field in the black brane geometry. We have
computed numerically the first few QNM as detailed in Appendix C. In Figure 3,
we show that the numerical results agree with the analytic formula in the
vicinity of spinodal points.

Figure 3: Comparison between the (dimensionless) decay rate computed from the
analytical formula (44) (red) and the numerical result obtained as described
in appendix C (blue) for temperatures close to the spinodal temperature
$T_{2}$.

###  III.2 Correlation length and dynamical exponent

In order to compute the correlation length and dynamical exponent, we expand
the perturbation in plane waves,

$\delta\phi(v,Z,X)=\int\frac{d\bar{\omega}d^{3}\bar{k}}{(2\pi)^{3}}\delta\phi_{\omega,k}(Z)e^{-i\bar{\omega}v+i\bar{k}\cdot
X}\ ,$ (46)

where $\bar{\omega}=3\omega/(4\pi T)$ and $\bar{k}=3k/(4\pi T)$.

The two-point function for the order parameter has a pole corresponding to the
lowest QNM; we find that the associated dispersion relation reads, for small
momenta,

$\omega=-i\Gamma_{T}-iD_{T}k^{2}+\cdots\ ,$ (47)

where $D_{T}$ depends only on $T$, and $\Gamma_{T}\approx\tau^{-1}$
corresponds to the decay rate previously introduced and shown in Figure 3. At
zero frequency, the Fourier transform of the two-point function has the form

$\left\langle\psi(k)\psi(-k)\right\rangle\approx\frac{Z}{k^{2}+\xi^{-2}}\ ,$
(48)

where the pole at imaginary momentum $k=\pm i\xi^{-1}$ determines the
correlation length, implying that

$\xi^{-2}\approx\frac{\Gamma_{T}}{D_{T}}\ .$ (49)

Close to the critical and spinodal points the correlation length diverges and,
as shown in appendix B, we obtain the follwing approximate analytic
expression:

$\xi^{-2}\approx\frac{9}{2I_{2}}\left(\frac{4\pi
T}{3}\right)^{1/2}V^{\prime\prime}(\psi_{0})\sim\Gamma\sim\tau^{-1}\ ,\qquad
I_{2}\approx 1.829\ .$ (50)

This shows that the exponent of the momentum in the dispersion relation
coincides with the dynamical exponent $z=2$. Thus, it can be obtained by
measuring the small momentum dependence of the lowest QNM frequency. We also
obtain the dimensionless ratio

$D\equiv\frac{4\pi T}{3}D_{T}=\frac{4\pi
T}{3}\frac{\xi^{2}}{\tau}\approx\frac{I_{2}}{1-2I_{1}}\approx 2.548\ .$ (51)

Our numerical results are shown in Figure 4, where we show the momentum
dependence of the frequency of the lightest QNM frequency for the first- and
second-order phase transitions in the first two lines of Table 1. We have
added a $k^{2}$ curve (gray line), showing that the dynamical exponent is
$z=2$ for both kind of transitions. The numerical calculation also confirms
that the dimensionless coefficient $D$ extracted from (47) is very close in
both kinds of phase transitions, and agrees with the analytical estimate (51)
$D\approx 2.548$.

Figure 4: Momentum dependence of the lightest QNM frequency for (on the left)
the first-order transition in the first row of Table 1 near the spinodal
temperature, and (on the right) the second-order transition in the second row
of Table 1 near the critical temperature.

##  IV Breakdown of quasistatic time evolution and scaling

In this section we depart from equilibrium and allow the temperature to change
slowly in time. The background geometry is approximately (2) with a time-
dependent position of the horizon $z_{h}$, related to the temperature through
(3). Applying the change of variables introduced at the beginning of section
III, we arrived at (37), where all explicit dependence on the horizon position
has been removed from the equation of motion of the scalar. Thus, the time
evolution is completely encoded in the boundary conditions of the scalar
field. In the following we will restrict to homogeneous configurations and
drop the terms with derivatives along the spatial field theory directions
$x^{i}$.

###  IV.1 Quasistatic evolution

If the change of the temperature is slow enough, the solution for the scalar
field will in principle remain close to the static solution (38), with small
corrections proportional to the time derivative of the expectation value
and/or temperature. The solution for the scalar field can be split in
“equilibrium” and “out-of-equilibrium” terms,

$\phi(v,Z)=\phi_{0}(v,Z)+\delta\phi(v,Z)\ ,$ (52)

where the equilibrium solution is taken to be

$\phi_{0}(v,Z)=-\frac{\overline{\psi}(v)}{3}\eta_{R}(Z)\
,\qquad\overline{\psi}=\left(\frac{4\pi T}{3}\right)^{-3/4}\psi\ .$ (53)

The out-of-equilibrium term can be expanded in terms proportional to time
derivatives of the equilibrum solution

$\delta\phi(v,Z)\sim\partial_{v}\overline{\psi}\,\delta\phi^{(1)}(Z)+\partial_{v}^{2}\overline{\psi}\,\delta\phi^{(2)}(Z)+\cdots,$
(54)

that can be solved iteratively oder by order in equation (37). Then, at the
leading order in this expansion, the equation of motion for the out-of-
equilibrium term is

$\partial_{Z}\left(\frac{F(Z)}{Z^{2}}\partial_{Z}\delta\phi\right)+\frac{27/16}{Z^{4}}\delta\phi=\frac{1}{3}\partial_{v}\overline{\psi}(v)\left[\frac{2}{Z^{3}}\eta_{R}(Z)-\frac{2}{Z^{2}}\eta_{R}^{\prime}(Z)\right]+O(\partial_{v}^{2}\overline{\psi})\
.$ (55)

The solution is of the same form as the one used to estimate the decay rate
(100), replacing
$-\delta\overline{\psi}_{0}\Gamma\to\partial_{v}\overline{\psi}$. This allows
for an immediate translation of the results in appendix B to this case. The
leading and subleading terms in the boundary expansion of the solution are

$\phi_{-}(v)=-\frac{1}{3}\overline{\psi}\
,\qquad\phi_{+}(v)\approx\frac{\gamma}{3}\overline{\psi}+\frac{4}{9}\left(\frac{1}{2}-I_{1}\right)\partial_{v}\overline{\psi}+O(\partial_{v}^{2}\overline{\psi})\
.$ (56)

Expanded to lowest order in derivatives, the no source condition $J=0$ becomes

$\frac{2}{9}\left(1-2I_{1}\right)\left(\frac{4\pi
T}{3}\right)^{9/4}\partial_{v}\left[\left(\frac{4\pi
T}{3}\right)^{-3/4}\psi\right]+\partial_{\psi}V(\psi)\approx 0\ .$ (57)

The expectation value is a small perturbation around the equilibrium value
$\psi(v)=\psi_{0}[T(v)]+\delta\psi(v)$, satisfying
$\partial_{\psi}V(\psi_{0})=0$ and

$\delta\psi\approx\frac{\partial_{v}T}{\Gamma[\psi_{0}]}\left(\frac{3}{4T}\psi_{0}-\partial_{T}\psi_{0}\right)\
,$ (58)

where $\Gamma[\psi_{0}]$ is given in (44). Even if the temperature changes
very slowly, $\partial_{v}T/T\ll 1$, when it approaches a critical or spinodal
value, the assumption that $\delta\psi$ is small breaks down, since $\Gamma\to
0$. This shows that critical slowing down results in the breakdown of the
quasistatic approximation.

###  IV.2 Critical and spinodal scaling

Critical points of second-order phase transitions exhibit characteristic
scaling properties that also affect the equilibrium states close to them. The
correlation length and relaxation time diverge with particular exponents of
the temperature difference close to the critical point, while the deviation
from the critical values goes to zero:

$\xi\sim|T-T_{c}|^{-\nu}\ ,\qquad\tau\sim|T-T_{c}|^{-z\nu}\
,\qquad\Delta\psi\sim|T-T_{c}|^{\beta}\ .$ (59)

Let us characterize the time evolution close to the critical point with a
linear function

$\frac{T_{c}-T}{T_{c}}=\frac{t}{\tau_{{}_{Q}}}\ ,$ (60)

where $\tau_{{}_{Q}}$ determines the rate of change of the temperature. The
larger $\tau_{{}_{Q}}$ is, the slower the time evolution. We have selected
$t=0$ as the point where the critical temperature is reached when one
approaches it from from higher temperatures. At some point in the time
evolution, the time remaining before reaching the critical point, $\Delta
t=|t|$, will be the same as the relaxation time, $\Delta t=\tau(t)\equiv
t_{*}$, so the system will not have time to equilibrate before undergoing the
phase transition. This is known as the freeze-out time [70]; note that by our
definition, $t_{*}$ is positive even though the freeze-out occurs at negative
times. We can use it, with the corresponding temperature $T_{*}=T(-t_{*})$, to
define when the quasistatic approximation breaks down and the system goes out
of equilibrium. The relaxation time then satisfies

$\tau=\tau_{c}\left(\frac{T_{*}-T_{c}}{T_{c}}\right)^{-z\nu}=\tau_{c}\left(\frac{t_{*}}{\tau_{{}_{Q}}}\right)^{-z\nu}=t_{*}\
.$ (61)

From this, we obtain the scalings

$t_{*}\sim\tau_{{}_{Q}}^{\frac{z\nu}{1+z\nu}},\ \
T_{*}\sim\tau_{{}_{Q}}^{-\frac{1}{1+z\nu}}\ .$ (62)

In the case of a closed system with a phase transition producing spontaneous
symmetry breaking, this type of scaling determines the size of domains in the
broken phase through the Kibble-Zurek mechanism [71, 70].

Our case differs in several ways from the type of transitions to which the
Kibble-Zurek mechanism applies. In the first place, the two phases have the
same symmetry realization, there is no spontaneous breaking of symmetry.
Secondly, our initial state is homogeneous and fluctuations are completely
suppressed, so there cannot be formation of domains. And lastly, we are
describing an open system, since the sector undergoing the phase transition is
coupled to a thermal bath.

Although the freeze-out argument is appealing, we need a more direct
characterization of the departure from equilibrium. A simple choice is to
compare the expectation value of the operator or the effective potential with
their equilibrium values. This, however, leads to a different scaling
exponent. Let us declare the quasistatic evolution to be valid as long as the
deviation from the equilibrium value in (58) remains bounded and relatively
small compared to the equilibrium value,

$\left|\frac{\delta\psi}{\psi_{0}}\right|<\epsilon<1\ ,$ (63)

for some fixed $\epsilon$.

Consider a first-order transition, and let $T_{a}$, $a=1,2$ be the two
spinodal temperatures. We can expand close to the spinodal values,
$\psi_{0}=\psi_{a}+\Delta\psi_{a}$, $T=T_{a}+\Delta T$. The following
conditions are satisfied:

$\partial_{\psi}V(\psi_{0},T)=0=\partial_{\psi}V(\psi_{a},T_{a})\
,\qquad\partial_{\psi}^{2}V(\psi_{a},T_{a})=0\ .$ (64)

From the condition that $\psi_{0}$ is an equilibrium value at temperature $T$
we obtain

$\Delta\psi_{a}\approx\pm(-\sigma_{a}\Delta T)^{1/2}\
,\qquad\sigma_{a}=\frac{2\partial_{T}\partial_{\psi}V(\psi_{a},T_{a})}{\partial_{\psi}^{3}V(\psi_{a},T_{a})}\
,$ (65)

fixing the value of the “critical” exponent $\beta=1/2$. Note that we can
approach a given spinodal point only from higher or lower temperatures, but
not from both. This is reflected in that $\Delta\psi$ is real only if $\Delta
T$ has the correct sign. The dimensionless rate (44) is approximately

$\Gamma_{a}\approx\frac{9}{2(1-2I_{1})}\left(\frac{4\pi
T_{a}}{3}\right)^{-3/2}\partial_{\psi}^{3}V(\psi_{a},T_{a})\Delta\psi_{a}\ .$
(66)

In order to approach the point from a (meta)stable branch, we should select
the sign such that $\partial_{\psi}^{3}V(\psi_{a},T_{a})\Delta\psi_{a}>0$.
Since $\Gamma_{a}$ is the inverse of the relaxation time, this gives
“critical” exponent $z\nu=1/2$ for the spinodal point.

Considering instead a second-order phase transition, close to the critical
temperature $T=T_{c}+\Delta T$, $\psi_{0}=\psi_{c}+\Delta\psi_{c}$, the
conditions on the derivatives of the effective potential are

$\partial_{\psi}V(\psi_{0},T)=0=\partial_{\psi}V(\psi_{c},T_{c})\
,\qquad\partial_{\psi}^{2}V(\psi_{c},T_{c})=0=\partial_{\psi}^{3}V(\psi_{c},T_{c})\
.$ (67)

Now the temperature dependence appears with a modified exponent

$\Delta\psi_{c}\approx(-\sigma_{c}\Delta T)^{1/3}\
,\qquad\sigma_{c}=\frac{6\partial_{T}\partial_{\psi}V(\psi_{c},T_{c})}{\partial_{\psi}^{4}V(\psi_{c},T_{c})}\
,$ (68)

fixing the critical exponent $\beta=1/3$. In this case we can always find a
real root for $\Delta\psi_{c}$. The dimensionless rate (44) can be
approximated by

$\Gamma_{c}\approx\frac{9}{4(1-2I_{1})}\left(\frac{4\pi
T_{c}}{3}\right)^{-3/2}\partial_{\psi}^{4}V(\psi_{c},T_{c})(\Delta\psi_{c})^{2}\
.$ (69)

Identifying $\Gamma_{c}$ with the inverse of the relaxation time results in
the exponents $z\nu=2/3$ for the critical point of the second-order phase
transition. Although the values for the exponents might look unusual, they are
mean field values for the effective potential (30). The results for the
scaling exponents of both types of transitions in our model are summarized in
Table 2.

Order |  $z$ |  $\beta$ |  $\nu$
---|---|---|---
1st | 2 | 1/2 | 1/4
2nd | 2 | 1/3 | 1/3
Table 2: Scaling exponents as defined in (59) for first and second-order phase
transitions in the holographic model.

From (58) and (63) we arrive at the following condition for the spinodal
points:

$\epsilon>\bar{\tau}_{a}\left|\frac{\partial_{v}T}{T-T_{a}}\right|\ .$ (70)

For a critical point we have instead

$\epsilon>\bar{\tau}_{c}\left|\frac{T_{c}^{1/3}\partial_{v}T}{(T-T_{c})^{4/3}}\right|\
.$ (71)

In the above equations, we have defined

$\bar{\tau}_{a}=\left(\frac{4\pi
T_{a}}{3}\right)^{3/4}\left|\frac{1-2I_{1}}{9\partial_{\psi}^{3}V(\psi_{a},T_{a})}\right|\
,\qquad\bar{\tau}_{c}=\left(\frac{4\pi
T_{c}}{3}\right)^{-3/4}\left|\frac{4(1-2I_{1})}{27(T_{c}\sigma_{c})^{1/3}\partial_{\psi}^{4}V(\psi_{c},T_{c})}\right|\
.$ (72)

Then, assuming a linear dependence on time close to the spinodal/critical
point $\Delta T_{a,c}=T_{a,c}t/\tau_{{}_{Q}}$, the time $t=-t_{*}$ when the
system goes out of equilibrium before reaching the spinodal/critical point is

$\displaystyle(t_{*})_{a}$ $\displaystyle\approx$ $\displaystyle\frac{3}{4\pi
T_{a}}\frac{\bar{\tau}_{a}}{\epsilon}\sim\tau_{{}_{Q}}^{0}\ ,$ (73)
$\displaystyle(t_{*})_{c}$ $\displaystyle\approx$
$\displaystyle\left(\frac{3}{4\pi
T_{c}}\right)^{3/4}\left(\frac{\bar{\tau}_{c}^{3}}{\epsilon^{3}}\tau_{{}_{Q}}\right)^{1/4}\sim\tau_{{}_{Q}}^{1/4}\
.$ (74)

The temperature when this happens is

$\displaystyle\left|\frac{T_{*}-T_{a}}{T_{a}}\right|$ $\displaystyle\approx$
$\displaystyle\left(\frac{3}{4\pi
T_{a}}\right)\frac{\bar{\tau}_{a}/\epsilon}{\tau_{{}_{Q}}}\sim\tau_{{}_{Q}}^{-1}\
,$ (75) $\displaystyle\left|\frac{T_{*}-T_{c}}{T_{c}}\right|$
$\displaystyle\approx$ $\displaystyle\left(\frac{3}{4\pi
T_{c}}\right)^{3/4}\left(\frac{\bar{\tau}_{c}/\epsilon}{\tau_{{}_{Q}}}\right)^{3/4}\sim\tau_{{}_{Q}}^{-3/4}\
.$ (76)

Although critical slowing down happens for both spinodal and critical points,
the critical scaling found close to each type of point is different. In both
cases it is different from (62). The general formulas for the scaling in terms
of the critical exponents are

$t_{*}\sim\tau_{{}_{Q}}^{\frac{z\nu-\beta}{1+z\nu-\beta}}\
,\qquad\left|\frac{T_{*}-T_{c}}{T_{c}}\right|\sim\tau_{{}_{Q}}^{-\frac{1}{1+z\nu-\beta}}\
,\qquad(\beta<1)\ .$ (77)

###  IV.3 Full time evolution

So far we have used the properties of the effective potential and near-
equilibrium solutions close to the critical/spinodal points to estimate the
approximate behaviour of the solutions as these points are slowly approached.
However, this does not give us precise quantitative information about how slow
the evolution should be in order to enter the scaling regime before going out
of equilibrium and how close to the spinodal/critical point this will happen.

In order to answer these questions we will now perform a full-fledged
numerical calculation of the time evolution and compare the results against
our analytical estimates of the scalings (73)-(76).

In order to proceed we need to specify the time dependence of the temperature
$T(v)$. We will pick as initial time of the evolution $v=0$, where the scalar
solution will be equal to the equilibrium solution (38) at the initial
temperature. In order for the equation of motion of the scalar to be satisfied
at the initial time we also impose $\partial_{v}T(0)=0$. We also want a
function that changes slowly to remain in the quasistatic regime and has a
linear behaviour close to the time $v_{a,c}$ where the spinodal/critical
temperature $T_{a,c}$ is reached, with a slope that can be tuned.

A possible function satisfying these requirements is

$\frac{T}{T_{\text{first}}}=\frac{T_{a,c}}{T_{\text{first}}}-\mathcal{B}\tanh\left[\frac{\sigma}{\mathcal{B}}(v-v_{a,c})+\frac{\sigma}{2v_{a,c}\mathcal{B}}(v-v_{a,c})^{2}\right]\
,$ (78)

with the parameter $\sigma$ determining the slope close to the
spinodal/critical point and $\mathcal{B}$ determining the initial temperature

$\tau_{{}_{Q}}=\frac{3}{4\pi
T_{\text{first}}\sigma}\,,\quad\frac{T(0)}{T_{\text{first}}}=\frac{T_{a,c}}{T_{\text{first}}}+\mathcal{B}\tanh\left(\frac{\sigma
v_{a,c}}{2\mathcal{B}}\right)\ .$ (79)

In these formulas we use as a reference scale $T_{\text{first}}\approx
4.934M$, which is the critical temperature of the first-order transition
corresponding to the first row in Table 1.

In order to avoid starting the evolution too close to the spinodal/critical
point, we will require that as the evolution becomes slower $\sigma\ll 1$, the
time to reach the spinodal/critical point is larger $v_{a,c}\gg 1$. In
addition, we will also require starting at a higher temperature so that, if
the system wants to move away from equilibrium before reaching the
spinodal/critical time, it has enough time to do so. This last condition
requires $\sigma v_{a,c}/\mathcal{B}\gtrsim 1$ and $\mathcal{B}\gg 1$. These
conditions can be achieved by imposing the following relations among the
parameters

$\sigma=\frac{2}{v_{a,c}}\left(4\log_{10}v_{a,c}-\frac{T_{a,c}}{T_{\text{first}}}\right)\
,\qquad\frac{T(0)}{T_{\text{first}}}=2\log_{10}v_{a,c}\ .$ (80)

This leaves $v_{a,c}$ as the only parameter of the time evolution.

The strategy we followed is, given a theory determined by the parameters in
the effective potential, choose a set of increasing values for $v_{a,c}$,
which corresponds to another set of increasing values of $\tau_{{}_{Q}}$. For
each value, we numerically solve the equations as explained in appendix D and
define the time $t_{*}$ at which the system falls out-of-equilibrium by (63)
for several sufficiently small values of $\epsilon$ that let us explore the
scaling regime.

####  IV.3.1 Scaling in first and second-order phase transitions

For the case of the first-order phase transition with parameters as in the
first row of Table 1, we obtained the results shown in Figure 5 for a few
choices of $\epsilon$. In gray we added the line corresponding to the scaling
that we estimated analytically for a first-order phase transition (73), (75).
We observe that all curves show a scaling compatible with the analytical
estimate, being its associated range in $\tau_{{}_{Q}}$ larger for smaller
values of $\epsilon$. As $\epsilon$ decreases, the scaling regime will start
at larger values of $\tau_{{}_{Q}}$ while deviations from it will take place
at even larger values. Additionally, the time to fall out-of-equilibrium
$t_{*}$ increases, letting the system probe the regime closer to the spinodal
point. Furthermore, in Figure 6, we show that the scaling of the temperature
and the time is compatible with the prediction at fixed $\tau_{{}_{Q}}$ and as
a function of $\epsilon$. The deviations observed at large values of
$\tau_{{}_{Q}}$ are related to the finite resolution of the spinodal point
determined by $\epsilon$. Namely, we cannot observe deviations from
equilibrium if the evolution is so slow that they happen at separations from
the spinoidal point smaller than $\epsilon$.

This is consistent with the scaling region becoming larger at smaller
$\epsilon$. As a consequence of this limitation in resolution, for larger
values of $\epsilon\gtrsim 10^{-3}$ the scaling region ceases to be
observable. Taking this into account, we see that the cooling rate should be
rather small in order for the scaling to be observed
$1/\tau_{{}_{Q}}\lesssim(10^{-6}-10^{-8})T_{2}$, depending on $\epsilon$. On
the other hand, the temperature where the scaling regime is entered does not
seem to be as sensitive to the value of $\epsilon$ as the cooling rate.
Indeed, declaring that the scaling region starts when the deviation of the
temperature scaling from the theoretical value is $10\%$, for all the cases
studied in Figure 5, we find that the scaling regime starts when
$(T_{*}-T_{2})/T_{2}\approx 0.07-0.08$.

Figure 5: Scaling of $T_{*}$ and $t_{*}$ with the characteristic cooling rate
$\tau_{{}_{Q}}$ for the first-order phase transition with parameters displayed
in the first row of Table 1. The different curves were obtained for values of
epsilon equal to $3\cdot 10^{-4}$ (blue), $4\cdot 10^{-5}$ (orange), $7\cdot
10^{-6}$ (green), and $10^{-6}$ (red). The gray lines correspond to the
theoretical expectation for the scaling of both time and temperature for a
first-order phase transition.

Figure 6: First-order phase transition (first row of Table 1). Scaling of
$T_{*}$ and $t_{*}$ with $\epsilon$ for time evolution for $\tau_{Q}\approx
1.16\cdot 10^{8}T_{2}^{-1}$. The red dots correspond to the data extracted
from the numerical evolutions. The blue lines correspond to the theoretical
curves.

In the case of the second-order phase transition, whose parameters are given
in the second row of Table 1, we expect to find a different scaling, given in
(74), (76). We show the results for the same choices of $\epsilon$ as earlier
in Figure 7, in which the gray lines represent the scaling estimated
analytically. We observe that our numerical results exhibit scaling regimes
and agree with the analytical estimate. Once again, as $\epsilon$ is
decreased, the system requires greater $\tau_{{}_{Q}}$ to enter the scaling
region. However, this time there are no clear signs of deviations for the
three smaller choices of $\epsilon$ for the largest $\tau_{{}_{Q}}$ shown
here. Similarly to the previous case, we find a compatible scaling with
$\epsilon$ at fixed $\tau_{{}_{Q}}$, as shown in Figure 8.

Compared to the first-order phase transition, the necessary cooling rate to
observe the scaling region is still comparably small
$1/\tau_{Q}\lesssim(10^{-7}-10^{-9})T_{c}$, although the scaling region is
still relatively large for the larger value of $\epsilon$ shown in Figure 7.
The temperature seems also to have a weak dependence on $\epsilon$ with the
onset of the scaling regime closer to the critical temperature than in the
first-order phase transition. Using the same criterion as for the first-order
phase transition, we find that for all the cases shown in Figure 7 the scaling
region starts at $(T_{*}-T_{c})/T_{c}\approx 0.01$.

Figure 7: Scaling of $T_{*}$ and $t_{*}$ with the characteristic cooling rate
$\tau_{{}_{Q}}$ for the second-order phase transition with parameters
displayed in the second row of Table 1. The different curves were obtained for
values of epsilon equal to $3\cdot 10^{-4}$ (blue), $4\cdot 10^{-5}$ (orange),
$7\cdot 10^{-6}$ (green), and $10^{-6}$ (red). The gray lines correspond to
the theoretical expectation for the scaling of both time and temperature for a
second-order phase transition.

Figure 8: Second-order phase transition (second row of Table 1). Scaling of
$T_{*}$ and $t_{*}$ with $\epsilon$ for time evolution with
$\tau_{Q}=1.16\cdot 10^{11}T_{c}^{-1}$. The red dots correspond to the data
extracted from the numerical evolutions. The blue lines correspond to the
theoretical curves.

####  IV.3.2 Scaling in crossovers and weak first-order phase transitions

As the second-order phase transition is transformed into a crossover, we
expect to stop observing a scaling regime. However, as long as the crossover
is strong enough, one could in principle observe some approximate scaling in a
limited region of cooling rates. This is displayed in Figure 9, where we show
time and temperature in which the system falls out-of-equilibrium as a
function $\tau_{{}_{Q}}$ for two different crossovers, given by the parameters
in the third and fourth rows of Table 1, and for a single choice of
$\epsilon=10^{-3}$.

We observe that the stronger crossover still exhibits a region at large values
of $\tau_{{}_{Q}}$ in which there is an approximate scaling not far from the
one expected for a second-order phase transition. The behaviour of $T_{*}$ and
$t_{*}$ for the weaker crossover closely follows the stronger one, but
noticeable departures from it take place at around $\tau_{{}_{Q}}\sim
10^{7}T_{c}^{-1}$, not exhibiting a clear scaling regime.

As $\tau_{{}_{Q}}$ increases, a big drop in the values of $t_{*}$ and $T_{*}$
happens and deviation from equilibrium are not observed any further. The
reason behind it is that, as the second derivative of the potential does not
vanish at $T_{c}$, we can always remain closer and closer to equilibrium in an
evolution by increasing more and more $\tau_{{}_{Q}}$. With the specified
criterion of $\epsilon=10^{-3}$, the system remains always in equilibrium for
larger values of $\tau_{{}_{Q}}$ than the ones for which we displayed data in
Figure 9. We checked that decreasing $\epsilon$ increases the range of cooling
rate for which we find data, but this range always remains finite. As it is
expected, we find data for a larger range of $\tau_{{}_{Q}}$ for the stronger
crossover as it is closer to a second-order phase transition.

Figure 9: Scaling of $T_{*}$ and $t_{*}$ with the characteristic cooling rate
$\tau_{{}_{Q}}$ for the crossovers with $A=5\cdot 10^{-3}$ (red), $A=5\cdot
10^{-4}$ (orange), and $A=5\cdot 10^{-5}$ (blue). Here, $\epsilon=10^{-3}$.
The gray line corresponds to the theoretical expectation for the scaling of
both time and temperature for a second-order phase transition.

A related interesting case arises when considering weak first-order phase
transitions, meaning those that are close to being a second-order one. In such
case one would expect to find a first-order phase transition scaling regime
for slow enough evolution of the temperatures. However, at intermediate
cooling rates, when the system falls out-of-equilibrium not too close to the
spinodal point, one may expect to find a scaling close to that of a second-
order phase transition if the transition is weak enough. The three examples we
will consider to explore this possibility are defined by the parameters given
in the last three rows of Table 1.

Figure 10 shows the dependence of the time and temperature in which the field
falls out of equilibrium for the choice of $\epsilon=10^{-5}$ and for the
three mentioned theories. We have added for visual aid a solid and dashed gray
lines representing the second and first-order scalings respectively. These
lines were forced to pass through the curve for the weakest and strongest
transitions respectively.

The results show that for fast coolings, all three curves agree and seem to
start following the second-order scaling (solid gray line). The weaker the
transition, the closer it gets to exhibit such scaling, showing deviations for
larger values of $\tau_{{}_{Q}}$. These deviations at large $\tau_{{}_{Q}}$
seem to be getting the curves towards the first-order scaling (dashed gray
line). For the strongest transition we studied (green dots), the first-order
scaling does clearly exhibit itself in the studied range of $\tau_{{}_{Q}}$.
The mid strength transition results (orange dots) show the scaling at the very
last values of cooling rate, being this a region that is likely prolonged
toward higher values of $\tau_{{}_{Q}}$.

Our results are compatible then with the idea that weak first-order phase
transition exhibit a second-order scaling before showing the first-order one
for very slow temperature evolution.

Figure 10: Scaling of $T_{*}$ and $t_{*}$ with the characteristic cooling rate
$\tau_{{}_{Q}}$ for the weak phase transition with $A=5\cdot 10^{-3}$ (green),
$A=5\cdot 10^{-4}$ (orange), and $A=5\cdot 10^{-5}$ (blue). For all the
curves, $\epsilon=10^{-5}$. The gray line correspond to the theoretical
expectation for the scalings of a first-order (dashed) and second-order phase
(solid) transitions.

Finally, Figure 11 compares $T_{*}$ and $t_{*}$ for different values of
$\epsilon$ for the first-order phase transition with $A=5\cdot 10^{-4}$ (see
Table 1). In order for the curves to overlap, we have made a shift in both
axes so that the intermediate region falls in the same range.

We observe that all curves lie on top of each other. This implies that the
intermediate region of the curves, before the first-order phase scaling
arises, follows an equivalent behaviour and changing $\epsilon$ does not help
finding closer scalings to a second-order phase transition. The only change
that the choice of $\epsilon$ produces is the displacement of the region in
which the first-order scaling emerges to larger values of $\tau_{{}_{Q}}$.

In conclusion, our results show that the proximity to a second-order phase
transition is the main factor determining whether an intermediate second-order
scaling region manifests itself.

Figure 11: Comparison of $T_{*}$ and $t_{*}$ as a function of the cooling rate
$\tau_{Q}$ for the weak first-order phase transition with $A=5\cdot 10^{-4}$
for different choices of epsilon: $10^{-4}$ (blue), $10^{-5}$ (orange),
$10^{-6}$ (green), and $10^{-7}$ (red). We have performed a shift in both axes
for the curves to overlap.

##  V Discussion and outlook

Our results highlight similarities in the behaviour of a system when driven
slowly through a first- or a second-order phase transition, if fluctuations
are suppressed: In both cases there is critical slowing down of fluctuations
and the departure from equilibrium obeys a scaling law related to the
equilibrium “critical” exponents of critical and spinodal points. In addition,
we have obtained specific approximate values for the relative temperature when
the scaling regimes are entered, and it will be interesting to explore how
much these are model dependent or whether they hold more generally than the
cases we have studied.

In practice, whether the scaling regime can be reached will depend on the
relative suppression of fluctuations, which depend on the properties of the
particular system undergoing the transition. However, the fact that spinodal
scaling has been observed in experiments shows that understanding this regime
is not a purely academic exercise, but that can also be relevant for real
systems.

It should be noted that there can be different behaviours to the one presented
here for a first-order phase transition when it is the volume of a closed
system rather than the temperature of an open system the parameter that is
evolving. In the first case the energy density can be equivalently used as
time-evolving parameter so that, in contrast to the open system, there is a
single homogeneous equilibrium state at each time and in principle the
unstable state lying between the spinodal points could be reached.

There are several interesting directions one can pursue using the simple
holographic model described here:

  * •

Although we have studied only the case of undercooling, one can similarly
study the complementary overheating case, although we expect no qualitative
differences.

  * •

We have focused on the approach towards the spinodal/critical points. Evolving
further in time, one could study possible freeze-out behavior, as well as the
approach to the new equilibrium state.

  * •

Our model could easily be adapted to phase transitions involving spontaneous
symmetry breaking. A crucial difference with the transitions we have studied
here is that the symmetric high temperature state is still present as an
unstable state at low temperatures. The non-trivial evolution to the ground
state then happens through (small, random) perturbations. Adding spatial
dependence, one could realize the Kibble-Zurek mechanism of domain formation.

  * •

Adding spatial dependence to our field, there would be other interesting
inhomogeneous configurations whose time-evolution could be studied, including
the expansion of bubbles [17] and spinodal decomposition [11].

  * •

One could simulate fast, instead of slow, quenches333Quenches that are fast
compared to the relaxation rate of the scalar field but that satisfy
$T^{\prime}(t)/T^{2}\ll 1$ can be captured with the formalism presented here.
Quenches not far from the spinodal point are good candidates to work in this
regime as $\Gamma_{T}\ll T$. Alternatively, one can realize fast quenches by
fixing the temperature and making the couplings of the theory time-dependent..
In particular, one could ask if the system can transition before reaching the
spinodal point, if the evolution is fast enough to give it sufficient energy
to jump the potential barrier. That is, could the phase transition be
completed while the false vacuum is still metastable, not by bubble nucleation
but using energy pumped in during the quench?

  * •

The model could also be used to study dynamical phase transitions [72] and
Floquet states [73] 444See [74, 75, 76, 77] for holographic realizations of
Floquet states. by turning on an oscillating source. It could be particularly
interesting when the oscillation happens around the critical or spinodal
points.

  * •

In previous work, some of us have discussed how to use holography to compute
the field theory’s effective action in a derivative expansion [20]. This could
easily be done in the current model; it would then be interesting to study how
well the numerical evolution reported here could be captured by solving the
equation of motion coming from the effective action (which is an ordinary
differential equation instead of a partial one).

  * •

Instead of a black hole background geometry one could consider an AdS soliton
geometry or a global AdS geometry. In this case the time evolution of the
confinement scale could produce quantum phase transitions.

In addition to the points listed above, there are several possible extensions
and generalizations of the holographic model that could be interesting to
consider and for which similar questions can also be addressed:

Phase transitions involving a single scalar could also be realized by
introducing a non-trivial potential in the gravitational action, rather than
by modifying boundary conditions as done in this work, see e.g. [67].

Treating the background evolution non-adiabatically would allow to consider
more general evolutions than the ones mentioned earlier, with which to probe
the far-from-equilibrium regime. Widely used examples of time-evolving
analytic backgrounds are the dual to a boost invariant plasma [78, 79] and
Vaidya geometries used to describe thermalization [80].

Finally, going beyond the probe approximation and taking into account the
backreaction of the scalar over the metric will in general modify the
effective potential and even could result in the introduction of non-analytic
terms, allowing for a generalization to higher dimensions as in e.g.[20].
Adding other types of fields, such as charged scalars and gauge fields will
extend the phase diagram and allow for transitions produced by changes in the
chemical potential rather than the temperature.

## Acknowledgments

We would like to thank Jorge Casalderrey Solana, Mark Hindmarsh, Javier Mas,
David Mateos and Alfonso Ramallo for valuable discussions and useful comments.
The work of O.H. is supported by the Research Council of Finland (grant number
330346) and the Waldemar von Frenckell foundation. The work of A.C. is
supported by Ministerio de Ciencia e Innovación de España under the program
Juan de la Cierva-formación. A.C. and C.H. are partially supported by the AEI
and the MCIU through the Spanish grant PID2021-123021NB-I00. The work of
M.S.G. is supported by the European Research Council (ERC) under the European
Union’s Horizon 2020 research and innovation program (grant agreement
No758759).

##  Appendix A Stationary points of the effective potential

Here we study in greater detail the equilibrium states of our system; these
can be extracted from the effective potential (30), which we reproduce here
for convenience:

$V(\psi)=\Lambda\psi+\frac{1}{2}a_{T}\psi^{2}+\frac{b}{3}\psi^{3}+\frac{c}{4}\psi^{4}\qquad\text{with}\qquad
a_{T}=a+\frac{\gamma}{3}\left(\frac{4\pi T}{3}\right)^{3/2}\ .$ (81)

Firstly, we take all coefficients $\\{\Lambda,a,b,c\\}$ real, and $c>0$ in
order to avoid a runaway behaviour of the potential. The no-source boundary
condition (32) has at least one and at most three solutions for $\psi$ real.
If there is only one solution, the effective potential has one minimum, while
if there are three real solutions the potential has two minima and a maximum.
Two solutions correspond to limiting cases where the maximum coincides with a
minimum and the potential has an inflection point.

It will be convenient to write all dimensionful quantities in units of a
reference mass scale $M$:

$a_{T}=\bar{a}_{T}M^{3/2}\ ,\qquad b=\bar{b}M^{3/4}\
,\qquad\Lambda=\bar{\Lambda}M^{9/4}\ ,\qquad\psi=\bar{\psi}M^{3/4}\ .$ (82)

The cubic equation (32) has three real roots when the discriminant is
positive. In terms of the dimensionless parameters the discriminant becomes

$\text{Discr}=\frac{1}{27c^{3}}\left[4(\bar{b}^{2}-3\bar{a}_{T}c)^{3}-(2\bar{b}^{3}-9\bar{a}_{T}\bar{b}c+27c^{2}\bar{\Lambda})^{2}\right]\
.$ (83)

We see that for $T/M\to\infty$, $\bar{a}_{T}\sim(T/M)^{3/2}\to\infty$ and
$\text{Discr}\sim-\bar{a}_{T}^{3}$, so at high temperatures the discriminant
is always negative, meaning that there is a single equilibrium state. When we
lower the temperature there is a possibility that the discriminant becomes
positive. This requires

$\bar{b}^{2}-3\bar{a}_{T}c>0\
,\qquad(\bar{b}^{2}-3\bar{a}_{T}c)^{3}\geq\frac{1}{4}(2\bar{b}^{3}-9\bar{a}_{T}\bar{b}c+27c^{2}\bar{\Lambda})^{2}\
.$ (84)

Let us identify $T_{1}$ as the temperature where a second real root first
appears. We introduce the parametrization

$a=-\frac{\gamma}{3}\left(\frac{4\pi
T_{1}}{3}\right)^{3/2}+\frac{b^{2}}{3c}(1-A)\
,\qquad\Lambda=\frac{b^{3}}{27c^{2}}\left[1-A(B+2)\right]\ .$ (85)

Then, at $T=T_{1}$,

$\bar{a}_{T_{1}}=\frac{\bar{b}^{2}}{3c}(1-A)\ ,$ (86)

which gives the condition

$A^{2}\left[4A-(B-1)^{2}\right]=0\ .$ (87)

So we should impose $A=(B-1)^{2}/4$. Let us now identify $T=T_{2}<T_{1}$ as
the temperature where two of the real roots coincide before moving to the
complex plane and leaving a single real root. In this case

$\bar{a}_{T_{2}}=\frac{\bar{b}^{2}}{3c}(1-A-\Delta A)\ ,$ (88)

where we have introduced a parametrization

$\frac{\gamma}{3}\left(\frac{4\pi}{3}\right)^{3/2}\frac{T_{1}^{3/2}-T_{2}^{3/2}}{M^{3/2}}=\frac{\bar{b}^{2}}{3c}\Delta
A>0\ .$ (89)

The second condition in (84) is saturated when $\Delta A=0$, and for

$\Delta
A_{\pm}=\frac{3}{4}+\frac{3}{8}B(2-B)\pm\frac{B+2}{8}\left[3(4-B^{2})\right]^{1/2}\
.$ (90)

The range where there are two positive solutions $\Delta A_{\pm}>0$ is
$1<B<2$. The temperature $T_{2}$ corresponds to $\Delta A_{-}$, and we must
demand that

$\frac{\bar{b}^{2}}{3c}\Delta
A_{+}>\frac{\gamma}{3}\left(\frac{4\pi}{3}\right)^{3/2}\frac{T_{1}^{3/2}}{M^{3/2}},\quad\Rightarrow\quad\frac{\bar{b}^{2}}{3c}(\Delta
A_{+}-\Delta
A_{-})>\frac{\gamma}{3}\left(\frac{4\pi}{3}\right)^{3/2}\frac{T_{2}^{3/2}}{M^{3/2}}\
,$ (91)

so that there is a single real root in the interval $T_{2}>T\geq 0$.

In the case where there is a second-order phase transition the potential
develops an inflection point but there is always a single minimum. This
corresponds to the discriminant having a double root, which happens at the
extreme of the allowed interval for $B$: $B=1$, $A=\Delta A_{-}=0$. The
critical temperature can be identified as $T_{c}=T_{1}$, and should satisfy
the constraint

$\frac{9}{4}\frac{\bar{b}^{2}}{3c}>\frac{\gamma}{3}\left(\frac{4\pi}{3}\right)^{3/2}\frac{T_{c}^{3/2}}{M^{3/2}}\
.$ (92)

A crossover can be obtained by taking $B=1$ and $A<0$ but small, such that the
discriminant remains negative at the would-be double root. The smaller $A$,
the more pronounced the crossover would be. In this case $T_{\rm cross}=T_{1}$
would serve as a crossover temperature, with a constraint similar to (92) save
for a small correction proportional to $A$.

##  Appendix B Analytic estimate of the decay rate and correlation length

In this appendix we estimate the decay rate in amplitude of the perturbation
and the correlation length. Let us introduce the differential operator acting
on static solutions

${\cal
D}_{Z}\,\phi\equiv\partial_{Z}\left(\frac{F(Z)}{Z^{2}}\partial_{Z}\,\phi\right)+\frac{27/16}{Z^{4}}\phi\
.$ (93)

The associated Green’s function is

${\cal D}_{Z}G(Z,Z^{\prime})=\delta(Z-Z^{\prime}),\ \
G(Z,Z^{\prime})=-\frac{2}{3}\left\\{\begin{array}[]{lc}\eta_{R}(Z)\eta_{9/4}(Z^{\prime}),&Z^{\prime}<Z\\\
\eta_{9/4}(Z)\eta_{R}(Z^{\prime}),&Z<Z^{\prime}\end{array}\right.$ (94)

Where the solutions $\eta_{9/4}(Z),\eta_{R}(Z)$ were defined in (25) and (28).

###  B.1 Decay rate

Using the Green’s function the equation for a homogeneous perturbation can be
recast as the integral equation

$\delta\phi=-\frac{1}{3}\delta\overline{\psi}(v)\eta_{R}(Z)-2\int_{z_{{}_{UV}}}^{1}dZ^{\prime}\,G(Z,Z^{\prime})\left(\frac{1}{{Z^{\prime}}^{3}}\partial_{v}\delta\phi-\frac{1}{{Z^{\prime}}^{2}}\partial_{Z^{\prime}}\partial_{v}\delta\phi\right)\
.$ (95)

We have introduced a cutoff $0<z_{{}_{UV}}\ll 1$ in the integration for
convenience; we will take $z_{{}_{UV}}\to 0$ at the end of the calculation.

We can assume an exponential decay of the perturbation

$\delta\phi(Z,v)=e^{-\Gamma v}\delta\phi_{0}(Z)\
,\qquad\delta\overline{\psi}(v)=e^{-\Gamma v}\delta\overline{\psi}_{0}\ .$
(96)

In this case, the integral equation simplifies to

$\delta\phi_{0}(Z)=-\frac{1}{3}\delta\overline{\psi}_{0}\eta_{R}(Z)+2\Gamma\int_{z_{{}_{UV}}}^{1}dZ^{\prime}\,G(Z,Z^{\prime})\left(\frac{1}{{Z^{\prime}}^{3}}\delta\phi_{0}-\frac{1}{{Z^{\prime}}^{2}}\partial_{Z^{\prime}}\delta\phi_{0}\right)\
.$ (97)

In principle a solution could be found using Neumann’s series

$\delta\phi_{0}(Z)=\sum_{n=0}^{\infty}\Gamma^{n}\chi_{n}(Z)\ .$ (98)

The first two orders in the expansion give

$\chi_{0}(Z)=-\frac{1}{3}\delta\overline{\psi}_{0}\eta_{R}(Z)\ ,$ (99)

and

$\chi_{1}(Z)=-\frac{2}{3}\delta\overline{\psi}_{0}\Gamma\int_{z_{{}_{UV}}}^{1}dZ^{\prime}\,G(Z,Z^{\prime})\left[\frac{1}{{Z^{\prime}}^{3}}\eta_{R}(Z^{\prime})-\frac{1}{{Z^{\prime}}^{2}}\eta_{R}^{\prime}(Z^{\prime})\right]\
.$ (100)

We will need the $z_{{}_{UV}}\to 0$ expansion of the integral,

$\int_{z_{{}_{UV}}}^{1}dZ^{\prime}\eta_{R}(Z^{\prime})\left[\frac{1}{{Z^{\prime}}^{3}}\eta_{R}(Z^{\prime})-\frac{1}{{Z^{\prime}}^{2}}\eta_{R}^{\prime}(Z^{\prime})\right]\simeq-\frac{1}{2}(1-z_{{}_{UV}}^{-1/2})+I_{1}+O(z_{{}_{UV}})\
,$ (101)

where $I_{1}\approx 0.141115$.

We now expand the solution for $Z=z_{{}_{UV}}\to 0$ and identify the
coefficients of the leading and subleading solutions

$\begin{split}\chi_{0}(z_{{}_{UV}})&\simeq-\frac{1}{3}\delta\overline{\psi}_{0}\left(z_{{}_{UV}}^{3/4}-\gamma
z_{{}_{UV}}^{9/4}\right)+\cdots\ ,\\\
\chi_{1}(z_{{}_{UV}})&\simeq\frac{4}{9}\delta\overline{\psi}_{0}\Gamma\left[\frac{z_{{}_{UV}}^{7/4}}{2}-\left(\frac{1}{2}-I_{1}\right)z_{{}_{UV}}^{9/4}\right]+\cdots\
.\end{split}$ (102)

From this we can read the coefficient of the subleading solution

$\delta\overline{\phi}_{+}(v)\simeq\left[]\frac{\gamma}{3}-\frac{4\Gamma}{9}\left(\frac{1}{2}-I_{1}\right)\right]\delta\overline{\psi}(v)\
.$ (103)

Then, imposing $\delta J=0$ in (42), we arrive at

$\left[]\frac{\gamma}{3}-\frac{4\Gamma}{9}\left(\frac{1}{2}-I_{1}\right)\right]\left(\frac{4\pi
T}{3}\right)^{3/2}+W^{\prime\prime}(\psi_{0})=0\ .$ (104)

Let us note now that, at equilibrium,

$0=J_{0}=V^{\prime}(\psi_{0})=\Lambda+W^{\prime}(\psi_{0})+\frac{\gamma}{3}\left(\frac{4\pi
T}{3}\right)^{3/2}\psi_{0}\ .$ (105)

Then, we can arrange terms in (104) such that the condition boils down to

$-\frac{4\Gamma}{9}\left(\frac{1}{2}-I_{1}\right)\left(\frac{4\pi
T}{3}\right)^{3/2}+V^{\prime\prime}(\psi_{0})=0\ .$ (106)

Hence, the decay rate is determined by

$\Gamma=\frac{9}{2(1-2I_{1})}\left(\frac{4\pi
T}{3}\right)^{-3/2}V^{\prime\prime}(\psi_{0})\ .$ (107)

###  B.2 Correlation length

We can restrict to time-independent solutions. For small spatial momentum the
zero frequency plane wave perturbation (46) can be expanded as

$\delta\phi_{k}(Z)=\delta\phi_{k}^{(0)}(Z)+\bar{k}^{2}\delta\phi_{k}^{(1)}(Z)+O(k^{4})\
.$ (108)

The equations for the first two terms in the expansion are

$\displaystyle{\cal D}_{Z}\delta\phi_{k}^{(0)}$ $\displaystyle=$
$\displaystyle 0\ ,$ (109) $\displaystyle{\cal D}_{Z}\delta\phi_{k}^{(1)}$
$\displaystyle=$ $\displaystyle\frac{1}{Z^{2}}\delta\phi_{k}^{(0)}\ ,$ (110)

with the solutions

$\displaystyle\delta\phi_{k}^{(0)}(Z)$ $\displaystyle=$
$\displaystyle-\frac{\delta\bar{\psi}_{k}}{3}\eta_{R}(Z)\ ,$ (111)
$\displaystyle\delta\phi_{k}^{(1)}(Z)$ $\displaystyle=$
$\displaystyle\int_{0}^{1}\frac{dZ^{\prime}}{{Z^{\prime}}^{2}}G(Z,Z^{\prime})\delta\phi_{k}^{(0)}(Z^{\prime})\
.$ (112)

In this case a radial cutoff is not necessary since the integral with the
Green’s function remains finite. The near-boundary expansion of the solution
to this order is

$\delta\phi_{k}(Z)\sim-\frac{\delta\bar{\psi}_{k}}{3}\left[Z^{3/4}-\left(\gamma+\frac{2I_{2}}{3}\right)Z^{9/4}\right]\
,$ (113)

where

$I_{2}=\int_{0}^{1}dZ^{\prime}\left(\frac{\eta_{R}(Z^{\prime})}{Z^{\prime}}\right)^{2}\approx
1.829\ .$ (114)

Reading the asymptotic coefficients from the expression above, the no source
condition (42) becomes

$\left(V^{\prime\prime}(\psi_{0})+\frac{2I_{2}}{9}\left(\frac{3}{4\pi
T}\right)^{1/2}k^{2}\right)\delta\psi_{k}=0\ .$ (115)

This is only satisfied for a purely imaginary value of the momentum that
determines the correlation length

$-k^{2}=\xi^{-2}=\frac{9}{2I_{2}}\left(\frac{4\pi
T}{3}\right)^{1/2}V^{\prime\prime}(\psi_{0})\ .$ (116)

##  Appendix C Quasinormal modes

In this appendix, we provide some details on the computation of the
quasinormal modes of our model. These modes describe the decay of the scalar
field fluctuations close to the spinodal points studied in section III.1.

We did this for both homogeneous, and non-homogeneous fluctuations satisfying
the scalar equation (40) and use the plane wave expansion (46). We impose
ingoing conditions at the horizon $Z=1$ and the no-source boundary condition
$\delta J=0$ at $Z=0$. In the following, it will be convenient to adopt the
coordinate $u=Z^{1/2}$ and to introduce the rescaled scalar field
$\Phi(u)=u^{-3/2}\delta\phi_{\omega,k}(u)$. The equation of motion to be
solved then reads

$4u\left(u^{6}-1\right)\Phi^{\prime\prime}(u)+8\left(2u^{6}-2iu^{2}\bar{\omega}+1\right)\Phi^{\prime}(u)+\left(9u^{5}+8iu\bar{\omega}+16u^{3}\bar{k}^{2}\right)\Phi(u)=0\
.$ (117)

Given the boundary conditions, only a discrete set of $\bar{\omega}$ give non-
trivial solutions for each $\bar{k}$.

To numerically determine these values, we adopt two paths. For the case of
vanishing momenta, $\bar{k}=0$, we adopt a standard shooting method. We take
the expansion of the field at the boundary and at the horizon,

$\Phi_{\text{b}}(u)=\sum_{n=0}^{\infty}j_{n}u^{n}\
,\quad\quad\quad\Phi_{\text{h}}(u)=\sum_{n=0}^{\infty}h_{n}(1-u)^{n}\ .$ (118)

Notice that by expanding $\Phi_{\text{h}}$ in this manner, we are selecting
the ingoing solution, as the outgoing one is not regular at the horizon if one
uses Eddington-Finklestein coordinates. Furthermore, the no-source boundary
condition $\delta J=0$ establishes a relationship between $j_{3}$ and $j_{0}$
according to

$j_{3}=3\left(\frac{4\pi}{3}\right)^{-3/2}\left[\frac{a}{T^{3/2}}+\frac{2b}{T^{3/4}}\left(\frac{\psi_{0}}{T^{3/4}}\right)+3c\left(\frac{\psi_{0}}{T^{3/4}}\right)^{2}\right]j_{0}\
,$ (119)

where $\psi_{0}$ is the expectation value of the equilibrium solution as
defined in (14). Solving the equations close to the boundary and the horizon,
we find the coefficients $j_{n}$ and $h_{n}$. The first ones read

$\displaystyle j_{1}=0\ ,\quad\quad\quad j_{2}=-i\bar{\omega}j_{0}\ ,$ (120)

and

$\displaystyle h_{1}=-\frac{8\bar{\omega}-9i}{8(2\bar{\omega}+3i)}h_{0}\
,\quad\quad
h_{2}=-\frac{64\bar{\omega}^{2}+1152i\bar{\omega}-135}{256\left(2\bar{\omega}^{2}+9i\bar{\omega}-9\right)}h_{0}\
.$ (121)

We then solve equation (117) numerically twice, utilizing the Taylor
expansions described above as boundary conditions. This yields two numerical
solutions, $\Phi_{\text{h}}(u)$ depending on the free coefficient $h_{0}$, and
$\Phi_{\text{b}}(u)$ depending on the free coefficient $j_{0}$. These two
solutions smoothly match at an arbitrary location $u=u_{0}$ along the radial
direction only for a discrete set of values of $\bar{\omega}$ that render the
Wronskian of the solutions zero. Therefore, by imposing

$\Phi_{\text{h}}(u_{0})\Phi_{\text{b}}^{\prime}(u_{0})-\Phi_{\text{b}}(u_{0})\Phi_{\text{h}}^{\prime}(u_{0})=0\
,$ (122)

we numerically find the quasinormal mode values of

$\bar{\omega}=\bar{\omega}_{R}-i\Gamma\ .$ (123)

For the case of finite momenta, $\bar{k}\neq 0$, we found the quasinormal mode
frequencies by solving (117) as a generalized eigenvalue problem (see e.g.
[81]). In this procedure, we discretize (117) by using a Chebyshev grid,
leading to the matrix problem

$(M_{0}+\bar{\omega}M_{1})\Phi=0\ ,$ (124)

with the matrices being

$M_{0}=4u(u^{6}-1)D^{2}_{u}+8(2u^{6}+1)D_{u}+(9u^{5}+16u^{3}\bar{k}^{2})\mathds{1},\quad
M_{1}=-16iu^{2}D_{u}+8iu\mathds{1}\ .$ (125)

Here, $D_{u}$ and $\mathds{1}$ refer to the derivative matrix on the Chebyshev
grid and the identity one respectively.

Once the problem has been written in this way, we solved the system (124)
using the LinearSolve routine in Mathematica. We did this for two different
number of grid points, what allows one to distinguish those modes that are
physical from those that are spurious. In all our cases we used $N_{u}=40$ and
$60$ grid points, with 30 digits of precision, keeping those eigenmodes whose
relative difference was below $10^{-10}$.

An example of the result of the QNM computation is shown in Figure 12, for two
different values of momenta, $\bar{k}=\\{0,0.161\\}$, and first-order phase
transition in the first row of Table 1. The computation was done for a
temperature very close to the spinodal point, $T=T_{2}+10^{-8}T_{c}$, so that
the relaxation time of the scalar field is practically zero.

Figure 12: Lowest quasinormal mode frequencies for (left) $\bar{k}=0$ and
(right) $\bar{k}=0.161$ and first-order phase transition in the first line of
Table 1. The chosen temperature is $T=T_{2}+10^{-8}T_{c}$.

##  Appendix D Numerical methods for the time evolution

In this appendix, we provide some details on the numerical approach that we
have taken in order to study the time evolution of the scalar field in our
holographic model.

Assuming only time and radial dependence for the scalar field,
$\phi=\phi(v,Z)$, the equation of motion we need to solve (37) becomes

$0=\partial_{Z}^{2}\phi+\left(\frac{F^{\prime}(Z)}{F(Z)}-\frac{2}{Z}\right)\partial_{Z}\phi+\frac{27/16}{Z^{2}F(Z)}\phi+\frac{2}{ZF(Z)}\partial_{v}\phi-\frac{2}{F(Z)}\partial_{Z}\partial_{v}\phi\
.$ (126)

We introduce $\varphi_{-}$ and $\varphi_{+}$ as the dimensionless leading and
subleading terms in the boundary expansion of the scalar field,555The relation
between $\varphi_{\pm}$ and $\phi_{\pm}$ used in section II is
$\phi_{\pm}=z_{h}^{-\Delta_{\pm}}\varphi_{\pm}=\left(\frac{4\pi
T}{3}\right)^{\Delta_{\pm}}\varphi_{\pm}\ .$ (127)

$\phi(v,Z)=Z^{3/4}\left[\varphi_{-}(v)+\partial_{v}\varphi_{-}(v)Z+\varphi_{+}(v)Z^{3/2}+\cdots\right]\
.$ (128)

Following the discussion in section II.3, we need to solve (37) imposing
regularity at the horizon $Z=1$ and the no-source, multitrace, boundary
condition $J=0$, with $J$ given by (19). This can be written in terms of
dimensionless quantities as follows

$\varphi_{+}+\widetilde{\Lambda}\left(\frac{T}{T_{\text{first}}}\right)^{-9/4}-3\widetilde{a}\left(\frac{T}{T_{\text{first}}}\right)^{-3/2}\varphi_{-}+9\widetilde{b}\left(\frac{T}{T_{\text{first}}}\right)^{-3/4}\varphi_{-}^{2}-27c\varphi_{-}^{3}=0\
.$ (129)

Here we have introduced the dimensionless coefficients

$\widetilde{\Lambda}=\left(\frac{4\pi
T_{\text{first}}}{3}\right)^{-9/4}\Lambda,\quad\widetilde{a}=\left(\frac{4\pi
T_{\text{first}}}{3}\right)^{-3/2}a,\quad\widetilde{b}=\left(\frac{4\pi
T_{\text{first}}}{3}\right)^{-3/4}b,$ (130)

with the parameters $\\{a,b,c,\Lambda\\}$ chosen as in Table 1 depending on
the case at hand. $T$ is the temperature evolving in time according to (78)
and $T_{\text{first}}\approx 4.934M$ is the critical temperature of the first-
order phase transition in the first row of the table.

From the numerical point of view, it is more convenient to work with
expansions displaying integer powers, so we change the holographic variable
according to $Z=u^{2}$. Moreover, since we need to impose the multitrace
boundary condition (129) that relates the subleading and the leading terms in
the field expansion, it turns out to be useful to formulate the problem
introducing the field $\chi(v,u)$ through

$\phi(v,u)=u^{3/2}\left[\varphi_{-}(v)+u^{2}\partial_{v}\varphi_{-}(v)+u^{3}\chi(v,u)\right]\
.$ (131)

The final equation of motion is,

$\displaystyle
81u^{5}\chi-4(4-10u^{6})\partial_{u}\chi+4u\left[4u\partial_{u}\partial_{v}\chi-(1-u^{6})\partial_{u}^{2}\chi+10\partial_{v}\chi\right]+$
$\displaystyle+\
49u^{4}\partial_{v}\varphi_{-}+9u^{2}\varphi_{-}+24\partial_{v}^{2}\varphi_{-}=0\
.$ (132)

whereas the boundary condition to be imposed is $\chi(0,v)=\varphi_{+}(v)$
with $\varphi_{+}(v)$ given by (129).

We decided to numerically solve the problem by using the “method of lines”.
That is, we discretized the holographic direction $u$, including the horizon
$u=1$ and the boundary $u=0$ as points of the discretization grid,
transforming (D) into a set of coupled ordinary differential equations for the
function values at each grid point, $\chi_{i}(v)=\chi(u_{i},v)$, with
$i=1,2...N$, where $N$ is the number of grid points. By replacing
$\partial_{u}$ with the numerical derivative matrix in (D), we obtained the
aforementioned system of equations.

We solved this system of ordinary differential equations for the variables
$\varphi_{-}(v)$ and $\chi_{i}(v)$ utilizing the NDSolve routine of
Mathematica, choosing two different differentiation matrices: finite-
differences and Chebyshev. We found that the results coming from the two
choices agree.

## References

  * [1] A. Pandav, D. Mallick, and B. Mohanty, Search for the QCD critical point in high energy nuclear collisions, Prog. Part. Nucl. Phys. 125 (2022) 103960, [arXiv:2203.07817].
  * [2] L. Du, A. Sorensen, and M. Stephanov, The QCD phase diagram and Beam Energy Scan physics: a theory overview, 2, 2024.  arXiv:2402.10183.
  * [3] M. B. Hindmarsh, M. Lüben, J. Lumma, and M. Pauly, Phase transitions in the early universe, SciPost Phys. Lect. Notes 24 (2021) 1, [arXiv:2008.09136].
  * [4] P. Athron, C. Balázs, A. Fowlie, L. Morris, and L. Wu, Cosmological phase transitions: From perturbative particle physics to gravitational waves, Prog. Part. Nucl. Phys. 135 (2024) 104094, [arXiv:2305.02357].
  * [5] 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 (2019), no. 6 066004, [arXiv:1805.00882].
  * [6] 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 (2018), no. 26 261601, [arXiv:1807.05175].
  * [7] R. Rougemont and W. Barreto, Holographic entropy production in a Bjorken expanding hot and dense strongly coupled quantum fluid, Phys. Rev. D 106 (2022), no. 12 126023, [arXiv:2207.02411].
  * [8] F. R. Ares, M. Hindmarsh, C. Hoyos, and N. Jokela, Gravitational waves from a holographic phase transition, JHEP 21 (2020) 100, [arXiv:2011.12878].
  * [9] Z.-R. Zhu, J. Chen, and D. Hou, Gravitational waves from holographic QCD phase transition with gluon condensate, Eur. Phys. J. A 58 (2022), no. 6 104, [arXiv:2109.09933].
  * [10] F. R. Ares, O. Henriksson, M. Hindmarsh, C. Hoyos, and N. Jokela, Gravitational Waves at Strong Coupling from an Effective Action, Phys. Rev. Lett. 128 (2022), no. 13 131101, [arXiv:2110.14442].
  * [11] Y. Bea, J. Casalderrey-Solana, T. Giannakopoulos, A. Jansen, S. Krippendorf, D. Mateos, M. Sanchez-Garitaonandia, and M. Zilhão, Spinodal Gravitational Waves, arXiv:2112.15478.
  * [12] A. Escrivà and J. G. Subils, Primordial black hole formation during a strongly coupled crossover, Phys. Rev. D 107 (2023), no. 4 L041301, [arXiv:2211.15674].
  * [13] Y. Chen, D. Li, and M. Huang, Bubble nucleation and gravitational waves from holography in the probe approximation, JHEP 07 (2023) 225, [arXiv:2212.06591].
  * [14] F. Bigazzi, A. Caddeo, A. L. Cotrone, and A. Paredes, Fate of false vacua in holographic first-order phase transitions, JHEP 12 (2020) 200, [arXiv:2008.02579].
  * [15] F. Bigazzi, A. Caddeo, A. L. Cotrone, and A. Paredes, Dark Holograms and Gravitational Waves, JHEP 04 (2021) 094, [arXiv:2011.08757].
  * [16] F. Bigazzi, A. Caddeo, T. Canneti, and A. L. Cotrone, Bubble wall velocity at strong coupling, JHEP 08 (2021) 090, [arXiv:2104.12817].
  * [17] Y. Bea, J. Casalderrey-Solana, T. Giannakopoulos, D. Mateos, M. Sanchez-Garitaonandia, and M. Zilhão, Bubble wall velocity from holography, Phys. Rev. D 104 (2021), no. 12 L121903, [arXiv:2104.05708].
  * [18] R. A. Janik, M. Jarvinen, and J. Sonnenschein, A simple description of holographic domain walls in confining theories — extended hydrodynamics, JHEP 09 (2021) 129, [arXiv:2106.02642].
  * [19] O. Henriksson, Black brane evaporation through D-brane bubble nucleation, Phys. Rev. D 105 (2022), no. 4 L041901, [arXiv:2106.13254].
  * [20] F. R. Ares, O. Henriksson, M. Hindmarsh, C. Hoyos, and N. Jokela, Effective actions and bubble nucleation from holography, Phys. Rev. D 105 (2022), no. 6 066020, [arXiv:2109.13784].
  * [21] Y. Bea, J. Casalderrey-Solana, T. Giannakopoulos, D. Mateos, M. Sanchez-Garitaonandia, and M. Zilhão, Domain collisions, JHEP 06 (2022) 025, [arXiv:2111.03355].
  * [22] R. A. Janik, M. Jarvinen, H. Soltanpanahi, and J. Sonnenschein, Perfect Fluid Hydrodynamic Picture of Domain Wall Velocities at Strong Coupling, Phys. Rev. Lett. 129 (2022), no. 8 081601, [arXiv:2205.06274].
  * [23] Y. Bea, J. Casalderrey-Solana, T. Giannakopoulos, A. Jansen, D. Mateos, M. Sanchez-Garitaonandia, and M. Zilhão, Holographic bubbles with Jecco: expanding, collapsing and critical, JHEP 09 (2022) 008, [arXiv:2202.10503]. [Erratum: JHEP 03, 225 (2023)].
  * [24] J.-C. Wang, Z.-Y. Yuwen, Y.-S. Hao, and S.-J. Wang, General bubble expansion at strong coupling, arXiv:2311.07347.
  * [25] L. Li, S.-J. Wang, and Z.-Y. Yuwen, Bubble expansion at strong coupling, Phys. Rev. D 108 (2023), no. 9 096033, [arXiv:2302.10042].
  * [26] M. Sanchez-Garitaonandia and J. van de Vis, Prediction of the bubble wall velocity for a large jump in degrees of freedom, arXiv:2312.09964.
  * [27] M. Attems, Y. Bea, J. Casalderrey-Solana, D. Mateos, M. Triana, and M. Zilhao, Phase Transitions, Inhomogeneous Horizons and Second-Order Hydrodynamics, JHEP 06 (2017) 129, [arXiv:1703.02948].
  * [28] 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 (2017), no. 26 261601, [arXiv:1704.05387].
  * [29] 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].
  * [30] L. Bellantuono, R. A. Janik, J. Jankowski, and H. Soltanpanahi, Dynamics near a first order phase transition, JHEP 10 (2019) 146, [arXiv:1906.00061].
  * [31] W.-C. Yang, C.-Y. Xia, H.-B. Zeng, and H.-Q. Zhang, Phase Separation and Exotic Vortex Phases in a Two-Species Holographic Superfluid, Eur. Phys. J. C 81 (2021), no. 1 21, [arXiv:1907.01918].
  * [32] 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].
  * [33] M. Attems, Holographic approach of the spinodal instability to criticality, JHEP 08 (2021) 155, [arXiv:2012.15687].
  * [34] X. Zhao, Z.-Y. Nie, Z.-Q. Zhao, H.-B. Zeng, Y. Tian, and M. Baggioli, Dynamical evolution of spinodal decomposition in holographic superfluids, JHEP 02 (2024) 184, [arXiv:2311.08277].
  * [35] M. J. Bhaseen, J. P. Gauntlett, B. D. Simons, J. Sonner, and T. Wiseman, Holographic Superfluids and the Dynamics of Symmetry Breaking, Phys. Rev. Lett. 110 (2013), no. 1 015301, [arXiv:1207.4194].
  * [36] P. Basu, D. Das, S. R. Das, and K. Sengupta, Quantum Quench and Double Trace Couplings, JHEP 12 (2013) 070, [arXiv:1308.4061].
  * [37] A. M. García-García, H. B. Zeng, and H. Q. Zhang, A thermal quench induces spatial inhomogeneities in a holographic superconductor, JHEP 07 (2014) 096, [arXiv:1308.5398].
  * [38] J. Sonner, A. del Campo, and W. H. Zurek, Universal far-from-equilibrium Dynamics of a Holographic Superconductor, Nature Commun. 6 (2015) 7406, [arXiv:1406.2329].
  * [39] S. R. Das and T. Morita, Kibble-Zurek Scaling in Holographic Quantum Quench : Backreaction, JHEP 01 (2015) 084, [arXiv:1409.7361].
  * [40] P. M. Chesler, A. M. Garcia-Garcia, and H. Liu, Defect Formation beyond Kibble-Zurek Mechanism and Holography, Phys. Rev. X 5 (2015), no. 2 021015, [arXiv:1407.1862].
  * [41] T. Ishii, E. Kiritsis, and C. Rosen, Thermalization in a Holographic Confining Gauge Theory, JHEP 08 (2015) 008, [arXiv:1503.07766].
  * [42] M. Natsuume and T. Okamura, Kibble-Zurek scaling in holography, Phys. Rev. D 95 (2017), no. 10 106009, [arXiv:1703.00933].
  * [43] C.-Y. Xia, Z.-Y. Nie, H.-B. Zeng, and Y. Zhang, Uniform quenching processes in a holographic s + p model with reentrance, Eur. Phys. J. C 81 (2021), no. 10 882, [arXiv:2102.01083].
  * [44] M. Flory, S. Grieninger, and S. Morales-Tejera, Critical and near-critical relaxation of holographic superfluids, arXiv:2209.09251.
  * [45] Q. Chen, Y. Liu, Y. Tian, B. Wang, C.-Y. Zhang, and H. Zhang, Critical dynamics in holographic first-order phase transition, JHEP 01 (2023) 056, [arXiv:2209.12789].
  * [46] Q. Chen, Y. Liu, Y. Tian, X. Wu, and H. Zhang, Quench dynamics in holographic first-order phase transition, Phys. Rev. D 108 (2023), no. 10 106017, [arXiv:2211.11291].
  * [47] C.-Y. Xia, H.-B. Zeng, A. Grabarits, and A. del Campo, Kibble-Zurek Mechanism and Beyond: Lessons from a Holographic Superfluid Disk, arXiv:2406.09433.
  * [48] O. Aharony, S. S. Gubser, J. M. Maldacena, H. Ooguri, and Y. Oz, Large N field theories, string theory and gravity, Phys. Rept. 323 (2000) 183–386, [hep-th/9905111].
  * [49] K. Binder, Theory of first-order phase transitions, Reports on Progress in Physics 50 (July, 1987) 783–859.
  * [50] A. Compagner, On pseudocritical exponents at endpoints of metastable branches, Physica 72 (Feb., 1974) 115–122.
  * [51] M. E. Fisher and A. N. Berker, Scaling for first-order phase transitions in thermodynamic and finite systems, Phys. Rev. B 26 (Sep, 1982) 2507–2513.
  * [52] C. Sasaki, B. Friman, and K. Redlich, Density fluctuations in the presence of spinodal instabilities, Phys. Rev. Lett. 99 (Dec, 2007) 232301\.
  * [53] N. Liang and F. Zhong, Renormalization-group theory for temperature-driven first-order phase transitions in scalar models, arXiv e-prints (Nov., 2015) arXiv:1511.08258, [arXiv:1511.08258].
  * [54] X. An, D. Mesterházy, and M. A. Stephanov, On spinodal points and Lee-Yang edge singularities, J. Stat. Mech. 1803 (2018), no. 3 033207, [arXiv:1707.06447].
  * [55] F. Zhong, Universal scaling in first-order phase transitions mixed with nucleation and growth, arXiv e-prints (Oct., 2017) arXiv:1710.05198, [arXiv:1710.05198].
  * [56] A. Banerjee and T. Bar, Finite-dimensional signature of spinodal instability in an athermal hysteretic transition, Physical Review B 107 (Jan., 2023).
  * [57] T. Bar, S. K. Choudhary, M. A. Ashraf, K. S. Sujith, S. Puri, S. Raj, and B. Bansal, Kinetic Spinodal Instabilities in the Mott Transition in V2 O3 : Evidence from Hysteresis Scaling and Dissipative Phase Ordering, Phys. Rev. Lett.  121 (July, 2018) 045701, [arXiv:1808.00693].
  * [58] Y. Zhu, J. Hoffman, C. E. Rowland, H. Park, D. A. Walko, J. W. Freeland, P. J. Ryan, R. D. Schaller, A. Bhattacharya, and H. Wen, Unconventional slowing down of electronic recovery in photoexcited charge-ordered la1/3sr2/3feo3, Nature Communications 9 (May, 2018).
  * [59] T. Furukawa, K. Kobashi, Y. Kurosaki, K. Miyagawa, and K. Kanoda, Quasi-continuous transition from a fermi liquid to a spin liquid in kappa-(et)2cu2(cn)3, Nature Communications 9 (Jan., 2018).
  * [60] S. Kundu, T. Bar, R. K. Nayak, and B. Bansal, Critical slowing down at the abrupt mott transition: When the first-order phase transition becomes zeroth order and looks like second order, Physical Review Letters 124 (Mar., 2020) [arXiv:2301.08254].
  * [61] S. Pal, K. Kumar, and A. Banerjee, Universal scaling of charge-order melting in the magnetic field–pressure–temperature landscape, Physical Review B 102 (Nov., 2020).
  * [62] T. Bar, A. Ghosh, and A. Banerjee, Suppression of spinodal instability by disorder in an athermal system, Phys. Rev. B 104 (Oct., 2021) 144102, [arXiv:2102.11347].
  * [63] S. Kumbhakar, S. Islam, Z. Mao, Y. Wang, and A. Ghosh, Glassy electrons at the first-order Mott metal-insulator transition, Phys. Rev. B 106 (Nov., 2022) L201112.
  * [64] S. Kundu, R. K. Patel, S. Middey, and B. Bansal, Dynamic hysteresis at a noisy saddle node shows power-law scaling but nonuniversal exponent, Physical Review E 108 (Aug., 2023).
  * [65] P. Breitenlohner and D. Z. Freedman, Positive Energy in anti-De Sitter Backgrounds and Gauged Extended Supergravity, Phys. Lett. B 115 (1982) 197–201.
  * [66] I. R. Klebanov and E. Witten, AdS / CFT correspondence and symmetry breaking, Nucl. Phys. B 556 (1999) 89–114, [hep-th/9905104].
  * [67] Y. Bea and D. Mateos, Heating up Exotic RG Flows with Holography, JHEP 08 (2018) 034, [arXiv:1805.01806].
  * [68] E. Witten, Multitrace operators, boundary conditions, and AdS / CFT correspondence, hep-th/0112258.
  * [69] I. Papadimitriou, Multi-Trace Deformations in AdS/CFT: Exploring the Vacuum Structure of the Deformed CFT, JHEP 05 (2007) 075, [hep-th/0703152].
  * [70] W. H. Zurek, Cosmological Experiments in Superfluid Helium?, Nature 317 (1985) 505–508.
  * [71] T. W. B. Kibble, Topology of cosmic domains and strings, Journal of Physics A Mathematical General 9 (Aug., 1976) 1387–1398.
  * [72] M. Heyl, Dynamical quantum phase transitions: a review, Rept. Prog. Phys. 81 (2018), no. 5 054001, [arXiv:1709.07461].
  * [73] N. Tsuji, Floquet States, arXiv:2301.12676.
  * [74] A. Biasi, P. Carracedo, J. Mas, D. Musso, and A. Serantes, Floquet Scalar Dynamics in Global AdS, JHEP 04 (2018) 137, [arXiv:1712.07637].
  * [75] T. Ishii and K. Murata, Floquet superconductor in holography, Phys. Rev. D 98 (2018), no. 12 126005, [arXiv:1804.06785].
  * [76] A. Garbayo, J. Mas, and A. V. Ramallo, Holographic Floquet states in low dimensions, JHEP 10 (2020) 013, [arXiv:2007.12115].
  * [77] M. Berenguer, A. Garbayo, J. Mas, and A. V. Ramallo, Holographic Floquet states in low dimensions (II), JHEP 12 (2022) 020, [arXiv:2209.03884].
  * [78] R. A. Janik and R. B. Peschanski, Asymptotic perfect fluid dynamics as a consequence of Ads/CFT, Phys. Rev. D 73 (2006) 045013, [hep-th/0512162].
  * [79] M. P. Heller, P. Surowka, R. Loganayagam, M. Spalinski, and S. E. Vazquez, Consistent Holographic Description of Boost-Invariant Plasma, Phys. Rev. Lett. 102 (2009) 041601, [arXiv:0805.3774].
  * [80] S. Bhattacharyya and S. Minwalla, Weak Field Black Hole Formation in Asymptotically AdS Spacetimes, JHEP 09 (2009) 034, [arXiv:0904.0464].
  * [81] A. Jansen, Overdamped modes in Schwarzschild-de Sitter and a Mathematica package for the numerical computation of quasinormal modes, Eur. Phys. J. Plus 132 (2017), no. 12 546, [arXiv:1709.09178].