aainstitutetext: Department of Physics and Helsinki Institute of Physics, PL
64, FI-00014 University of Helsinki, Finland.bbinstitutetext: Departament de
Física Quàntica i Astrofísica and Institut de Ciències del Cosmos (ICC),
Universitat de Barcelona, Martí i Franquès 1, ES-08028, Barcelona,
Spain.ccinstitutetext: Centro de Astrofísica e Gravitação – CENTRA,
Departamento de Física, Instituto Superior Técnico – IST, Universidade de
Lisboa – UL, Av. Rovisco Pais 1, 1049-001 Lisboa, Portugal ddinstitutetext:
Institució Catalana de Recerca i Estudis Avançats (ICREA), Passeig Lluís
Companys 23, ES-08010, Barcelona, Spain.eeinstitutetext: Departamento de
Matemática da Universidade de Aveiro and Centre for Research and Development
in Mathematics and Applications (CIDMA), Campus de Santiago, 3810-183 Aveiro,
Portugal

# Holographic Bubbles with Jecco: Expanding, Collapsing and Critical

Yago Bea b Jorge Casalderrey-Solana c Thanasis Giannakopoulos b Aron Jansen
b,d David Mateos b Mikel Sanchez-Garitaonandia e,c and Miguel Zilhão
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Cosmological phase transitions can proceed via the nucleation of bubbles that
subsequently expand and collide. The resulting gravitational wave spectrum
depends crucially on the properties of these bubbles. We extend our previous
holographic work on planar bubbles to cylindrical bubbles in a strongly-
coupled, non-Abelian, four-dimensional gauge theory. This extension brings
about two new physical properties. First, the existence of a critical bubble,
which we determine. Second, the bubble profile at late times exhibits a richer
self-similar structure, which we verify. These results require a new 3+1
evolution code called Jecco that solves the Einstein equations in the
characteristic formulation in asymptotically AdS spaces. Jecco is written in
the Julia programming language and is freely available. We present an outline
of the code and the tests performed to assess its robustness and performance.

##  1 Introduction

A first-order, thermal phase transition in the Early Universe would produce
gravitational waves that could be detected in current or future experiments.
Since the Standard Model of particle physics possesses no first-order
transitions Aoki:2006we ; Kajantie:1996mn ; Laine:1998vn ; Rummukainen:1998as
, the discovery of gravitational waves originating from a cosmological phase
transition would amount to the discovery of new physics beyond the Standard
Model. The transition may proceed via bubble nucleation (see e.g.
Hindmarsh:2020hop  for a review) or via the spinodal instability Bea:2021zol .
In this paper we will focus on the first case.

Maximising the discovery potential requires an accurate understanding of the
bubble properties. These range from the action of the critical bubble that
gets nucleated to the terminal velocity of expanding bubbles. The former
controls the nucleation rate, whereas the latter controls the characteristic
frequency of the produced gravitational waves. Computing these parameters from
first principles is challenging even in weakly coupled theories. The former
requires knowledge of the effective potential at finite temperature
Laine:2016hma ; Gould:2021ccf , whereas the latter requires an understanding
of out-of-equilibrium physics Moore:1995ua ; Bodeker:2017cim ; Hoche:2020ysm .

In Bea:2021zsu  we performed the first holographic calculation of the bubble
wall velocity in a strongly-coupled, non-Abelian, four-dimensional gauge
theory. Because of technical limitations, in this reference we focused on
planar bubbles, namely we imposed translational invariance along two of the
spatial directions, in such a way that the dynamics was effectively 1+1
dimensional in the gauge theory and 2+1 dimensional on the gravity side. In
this paper we will extend our analysis by imposing translational invariance
along only one of the spatial directions. Thus, the effective dynamics will be
2+1 dimensional in the gauge theory and 3+1 dimensional on the gravity side.
This will allow us to study bubbles that have the topology of a cylinder.
Since we impose translation invariance along the axis of the cylinder, we will
only plot the dependence of physical quantities on the two spatial directions
transverse to this axis. We emphasize that we will not impose any symmetries
on these directions, meaning that the dynamics on the plane transverse to the
cylinder axis will be completely general.

The extension from planar to cylindrical bubbles brings about two new physical
aspects. The first one is that the surface tension now plays a role. In
particular, we will be able to identify a critical bubble in which the inward-
pointing force due to the surface tension exactly balances the outward-
pointing force coming from the pressure difference between the inside and the
outside of the bubble. The second one is that the asymptotic, self-similar
profile of an expanding bubble possesses a richer structure than in the planar
case. We will verify this by plotting our holographic result for the gauge
theory stress tensor at late times as a function of the appropriate scaling
variable. We will also compare the holographic result with the hydrodynamic
approximation. As expected, we will find that hydrodynamics provides a good
approximation everywhere except at the bubble wall.

To obtain these results we have developed a new 3+1 evolution code called
Jecco that solves the Einstein equations in the characteristic formulation in
asymptotically anti de Sitter (AdS) spaces. The characteristic approach to
solving Einstein’s equations has a long history. It dates back to the Bondi-
Sachs formalism Bondi:1960jsa ; Sachs:1962wk , crucial to the modern
understanding of gravitational waves. For numerical applications, these
formulations provide advantages over more standard spacelike foliations in a
number of situations. In the context of extracting gravitational-wave
information, for instance, this approach exploits the fact that null
hypersurfaces reach future null infinity, thereby avoiding systematic errors
from extrapolation techniques. Further advantages of such formulations
include: the initial data are free (i.e. there is no need to solve elliptic
equations for the initial data); no second time derivatives (resulting in
fewer evolution variables); the field equations are conveniently cast as a set
of nested ordinary differential equations (ODEs) which can be efficiently
solved.

Though questions remain about the well-posedness of these formulations
Giannakopoulos:2020dih ; Giannakopoulos:2021pnh , characteristic codes have
shown remarkable stability. Indeed, the first ever long-term stable evolutions
of moving black holes was accomplished with a characteristic scheme
Gomez:1997pd . Applications of this approach include the Cauchy-Characteristic
extraction method for the computation of gravitational waveforms at future
null infinity, which has been numerically implemented in doi:10.1063/1.525904
; Bishop:1996gt ; Bishop:1997ik ; Gomez:1996ge ; Reisswig:2009rx ;
Handmer:2014qha . There is an extensive literature on this and related
subjects — see Winicour’s Living Review Winicour:2012znc  for an overview.

Despite all the successes and advantages of this approach, one serious
drawback that it faces is the possible formation of caustics, which typically
spoil the numerical simulation. This is particularly severe when evolving
binary black holes, and for this reason the characteristic approach in solving
Einstein’s equations lost some ground in favour of more traditional Cauchy
evolution schemes. More recently, though, the characteristic approach has
shown to be particularly well-adapted for evolutions in the Poincaré patch of
AdS spaces. Crucially for these simulations is the presence of a (non-compact)
planar horizon embedded in the asymptotically AdS space, effectively acting as
an infrared cut-off, which removes caustic formation from the computational
domain.

Through holography, this approach has facilitated the study of far-from-
equilibrium dynamics of strongly-coupled gauge theories, allowing for studies
of isotropization Chesler:2008hg ; Heller:2013oxa ; Gursoy:2016tgf ,
collisions of gravitational shockwaves (used as models for heavy-ion
collisions) Chesler:2010bi ; Casalderrey-Solana:2013sxa ; Chesler:2015wra ,
momentum relaxation Balasubramanian:2013yqa , turbulence Adams:2013vsa ,
collisions in non-conformal theories Attems:2016tby ; Attems:2017zam , phase
transitions and dynamics of phase separation Attems:2017ezz ; Janik:2017ykj ;
Attems:2019yqn ; Bellantuono:2019wbn ; Bea:2020ees ; Janik:2021jbq ;
Bea:2021ieq , collisions in theories with phase transitions Attems:2018gou ,
dynamical instabilities Gursoy:2016ggq , and even applications to
gravitational-wave physics Ahmadvand:2017xrw ; Ahmadvand:2017tue ;
Bigazzi:2020avc ; Ares:2020lbt ; Ares:2021nap  and bubble dynamics
Bigazzi:2020phm ; Bea:2021zsu ; Bigazzi:2021ucw ; Ares:2021ntv . See
Chesler:2013lia  for more references and a comprehensive overview of the
techniques involved, and also Bantilan:2012vu ; Bantilan:2020pay ;
Bantilan:2020xas  for equivalent approaches using Cauchy evolutions.

Here we present a new 3+1 code called Jecco (Julia Einstein Characteristic
Code) that solves Einstein’s equations in the characteristic formulation in
asymptotically AdS spaces. Jecco is written in the Julia programming language
and comes with several tools (such as arbitrary-order finite-difference
operators as well as Chebyshev and Fourier differentiation matrices) useful
for generic numerical evolutions. The evolution part of the code would allow
for the study of any of the problems mentioned in the previous paragraph;
herein, as mentioned in the beginning, we will focus on the study of bubble
dynamics. The code is publicly available and can be obtained from github
https://github.com/mzilhao/Jecco.jl and Zenodo jecco-2022 . To the best of our
knowledge, this is the first such freely available code (see however the
PittNull code Bishop:1998uk ; Babiuc:2010ze  for characteristic evolutions in
asymptotically flat spaces, freely available and distributed as part of the
Einstein Toolkit EinsteinToolkit:2020_05 ).

This paper is organized as follows. In Sec. 2.1 we introduce the class of
models to which our code can be applied, as well as the corresponding
equations of motion. In Sec. 2.2 we discuss the implementation of these
equations in the code and the numerical methods that we use. In Sec. 3 we
discuss our new results for cylindrical bubbles. In Sec. 4 we conclude with
some final remarks. The tests of our code are collected in Appendix A. We use
$G=c=\hbar=1$ units throughout.

##  2 Jecco: a new characteristic code for numerical holography

###  2.1 Equations

In this section we outline the theoretical background and equations that are
implemented in Jecco. Our approach is similar to that of Chesler:2013lia  and
generalises the code presented in Attems:2017zam  to the 3+1 dimensional case.
See also Winicour:2012znc  for an overview of the approaches and codes used in
the asymptotically flat setting.

####  2.1.1 Equations of motion and characteristic formulation

We consider a five-dimensional action consisting of gravity coupled to a
scalar field $\phi$ with a non-trivial potential $V(\phi)$. The action for
this Einstein-scalar model is

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

where $\kappa_{5}^{2}$ is the 5D gravitational coupling constant, which in our
units takes the value $\kappa_{5}^{2}=8\pi$. The resulting dynamical equations
of motion read

$\displaystyle E_{\mu\nu}$ $\displaystyle\equiv
R_{\mu\nu}-\frac{R}{2}g_{\mu\nu}-8\pi T_{\mu\nu}=0,$ (2) $\displaystyle\Phi$
$\displaystyle\equiv\square\phi-\partial_{\phi}V(\phi)=0,$

where

$8\pi T_{\mu\nu}=2\partial_{\mu}\phi\partial_{\nu}\phi-
g_{\mu\nu}\left(g^{\alpha\beta}\partial_{\alpha}\phi\partial_{\beta}\phi+2V(\phi)\right).$

Our potential $V(\phi)$ comes from a superpotential $W(\phi)$ with the form

$LW(\phi)=-\frac{3}{2}-\frac{\phi^{2}}{2}+\lambda_{4}\,\phi^{4}+\lambda_{6}\,\phi^{6}\,,$
(3)

and its explicit expression can be derived via

$V=-\frac{4}{3}W^{2}+\frac{1}{2}W^{\prime}{}^{2},$

resulting in

$\displaystyle L^{2}V(\phi)$
$\displaystyle=-3-\frac{3}{2}\phi^{2}-\frac{1}{3}\phi^{4}+\left(\frac{4\lambda_{4}}{3}+8\lambda_{4}^{2}-2\lambda_{6}\right)\phi^{6}+\left(-\frac{4\lambda_{4}^{2}}{3}+\frac{4}{3}\lambda_{6}+24\lambda_{4}\lambda_{6}\right)\phi^{8}$
$\displaystyle\quad{}+\left(18\lambda_{6}^{2}-\frac{8}{3}\lambda_{4}\lambda_{6}\right)\phi^{10}-\frac{4}{3}\lambda_{6}^{2}\,\phi^{12}\,.$
(4)

In these equations $\lambda_{4}$ and $\lambda_{6}$ are freely specifiable
dimensionless parameters related to the parameters $\phi_{M}$ and $\phi_{Q}$
used in e.g. Bea:2018whf ; Bea:2020ees  through

$\lambda_{4}=-\frac{1}{4\phi_{M}^{2}}\,,\qquad\lambda_{6}=\frac{1}{\phi_{Q}}\,.$
(5)

This potential has a maximum at $\phi=0$, where it admits an exact AdS
solution of radius $L$. For numerical purposes we set $L=1$. The holographic
dual field theory corresponds to a 3+1 dimensional conformal field theory
which is deformed by a source $\Lambda$ for the dimension-three scalar
operator $\mathcal{O}_{\phi}$ dual to the scalar field $\phi$. The
thermodynamical and near-equilibrium properties of this model were presented
in Attems:2016ugt ; Attems:2016tby ; Attems:2017ezz  for $\lambda_{6}=0$ and
in Bea:2018whf ; Bea:2020ees  for $\lambda_{6}\neq 0$.

Let us point out that even if here we will always make use of the particular
potential (4), the code implementation is such that more generic potentials
can be used provided that, for low values of the scalar field, they behave as

$L^{2}V(\phi)=-3-\frac{3}{2}\phi^{2}-\frac{\phi^{4}}{3}+\mathcal{O}\left(\phi^{6}\right).$
(6)

The constant term is fixed by the 4+1 dimensional AdS asymptotics and the
quadratic one is in correspondence with the scaling dimension of the dual
scalar operator $\mathcal{O}_{\phi}$. The quartic term, determined by the
other two in our case, ensures the absence of a conformal anomaly, which would
give rise to logarithms in the asymptotic expansions. Thence, a change in this
near boundary behaviour of the potential would alter the hard-coded asymptotic
expansions and variable redefinitions to be introduced in Secs. 2.1.2 and
2.1.3.

We now write the following 5D ansatz for the metric in Eddington-Finkelstein
(EF) coordinates

$\displaystyle ds^{2}=g_{\mu\nu}dx^{\mu}dx^{\nu}$
$\displaystyle=-Adt^{2}+2dt\left(dr+F_{x}dx+F_{y}dy\right)+S^{2}\Big{[}e^{-B_{1}-B_{2}}\cosh(G)dx^{2}$
(7)
$\displaystyle\quad{}+e^{B_{1}-B_{2}}\cosh(G)dy^{2}+2e^{-B_{2}}\sinh(G)dxdy+e^{2B_{2}}dz^{2}\Big{]},$

where all functions depend on the radial coordinate $r$, time $t$ and
transverse directions $x$ and $y$. Nothing depends on the coordinate $z$, so
this is effectively a 3+1 system. Physically, this means that in the gauge
theory we impose translation invariance along the $z$-direction together with
$z\rightarrow-z$ symmetry. Along the remaining $(t,x,y)$-directions general
dynamics is permitted. Note that we denote by $t$ the (ingoing) null bulk
coordinate usually labeled $v$ in EF coordinates. At the boundary, $t$ becomes
the usual time coordinate. The spatial part of the metric is written such that
$S$ encodes the area of constant $t$ and $r$ slices,

$\sqrt{g|_{dt,dr=0}}=S^{3}.$

We can recover the 2+1 system of Attems:2017zam  by setting

$\displaystyle F_{y}=G=0,$ $\displaystyle B_{1}$ $\displaystyle=\frac{3}{2}B,$
$\displaystyle B_{2}$ $\displaystyle=\frac{1}{2}B,\qquad\mathrm{or}$ (8)
$\displaystyle F_{x}=G=0,$ $\displaystyle B_{1}$
$\displaystyle=-\frac{3}{2}B,$ $\displaystyle B_{2}$
$\displaystyle=\frac{1}{2}B,$

for non-trivial dependence only along the $x$ or $y$ direction respectively.
The metric (7) is invariant under

$\displaystyle r$ $\displaystyle\rightarrow\bar{r}=r+\xi(t,x,y)\,,$ (9)
$\displaystyle S$ $\displaystyle\rightarrow\bar{S}=S\,,$ $\displaystyle B_{1}$
$\displaystyle\rightarrow\bar{B}_{1}=B_{1}\,,$ $\displaystyle B_{2}$
$\displaystyle\rightarrow\bar{B}_{2}=B_{2}\,,$ $\displaystyle A$
$\displaystyle\rightarrow\bar{A}=A+2\partial_{t}\xi(t,x,y)\,,$ $\displaystyle
F_{x}$ $\displaystyle\rightarrow\bar{F}_{x}=F_{x}-\partial_{x}\xi(t,x,y)\,,$
$\displaystyle F_{y}$
$\displaystyle\rightarrow\bar{F}_{y}=F_{y}-\partial_{y}\xi(t,x,y)\,.$

Plugging the ansatz (7) into (2) results in a nested system of ODEs in the
radial (holographic) direction $r$ at each constant $t$ that can be solved
sequentially. We illustrate this system in Table 1. Each row in the table
represents an equation, obtained from the particular combination of the
equations of motion (2) as indicated, that takes the form

$\left[A_{f}(t,u,x,y)\,\partial^{2}_{u}+B_{f}(t,u,x,y)\,\partial_{u}+C_{f}(t,u,x,y)\right]f(t,u,x,y)=-S_{f}(t,u,x,y),$
(10)

where $u\equiv 1/r$, $f$ is the corresponding function to be solved for and
the coefficients $A_{f}$, $B_{f}$, $C_{f}$ and $S_{f}$ are fully determined
once the preceding equations have been solved. Dotted functions denote an
operation defined as

$\dot{f}\equiv\left(\partial_{t}+\frac{A}{2}\partial_{r}\right)f,$ (11)

which are necessary to obtain this nested structure.

Table 1: Nested structure of the equations of motion. Function(s) | Combination
---|---
$S$ | $E_{rr}$
$F_{x}$ | $E_{rx}-g_{tx}E_{rr}$
$F_{y}$ | $E_{ry}-g_{ty}E_{rr}$
$\dot{S}$ | $E_{tr}-\frac{1}{2}g_{tt}E_{rr}$
$\dot{\phi}$ | $\Phi$
$A$ | $\frac{E_{zz}}{g_{zz}}+\left(g^{ry}g_{ty}+g^{rx}g_{tx}\right)E_{rr}+2g^{rx}\left(E_{rx}-g_{tx}E_{rr}\right)+2g^{ry}\left(E_{ry}-g_{ty}E_{rr}\right)$
|
$-4\left(E_{tr}-\frac{1}{2}g_{tt}E_{rr}\right)+2\frac{E_{xy}}{g_{xy}}+g_{xx}g^{xx}\left(\frac{E_{yy}}{g_{yy}}+\frac{E_{xx}}{g_{xx}}-2\frac{E_{xy}}{g_{xy}}\right)$
$\dot{B_{2}}$ | $E_{zz}$
$\dot{G}$ | $E_{xy}$
$\dot{B_{1}}$ | $E_{yy}$
$\ddot{S}$ | $E_{tt}-\frac{1}{2}g_{tt}E_{tr}-\frac{1}{2}g_{tt}\left(E_{tr}-\frac{1}{2}g_{tt}E_{rr}\right)$
$\dot{F_{x}}$ | $E_{tx}-\frac{1}{2}g_{tt}E_{rx}-g_{tx}\left(E_{tr}-\frac{1}{2}g_{tt}E_{rr}\right)$
$\dot{F_{y}}$ | $E_{ty}-\frac{1}{2}g_{tt}E_{ry}-g_{ty}\left(E_{tr}-\frac{1}{2}g_{tt}E_{rr}\right)$

There are three sets of (two) coupled equations, indicated in the table by the
absence of a separating line. These still take the form of (10), but now $f$
should be thought of as a vector of the two functions involved, as is the
source term $S_{f}$, while $A_{f}$, $B_{f}$ and $C_{f}$ become $2\times 2$
matrices. The equations themselves are lengthy and given in (72-80).

These equations need to be supplemented with boundary conditions specified at
the AdS boundary $u\equiv 1/r=0$, see Sec. 2.1.3. In addition, the functions
$B_{1}(t_{0},u,x,y)$, $B_{2}(t_{0},u,x,y)$, $G(t_{0},u,x,y)$ and
$\phi(t_{0},u,x,y)$ should be thought of as initial data which, unlike for
Cauchy-based approaches of solving Einstein’s equations, can be freely
specified provided they are consistent with AdS asymptotics.

####  2.1.2 Asymptotic expansions

The study of the near-boundary behaviour ($u\rightarrow 0$) of the functions
is relevant for two reasons. The first one is that, as usually for
asymptotically AdS (AAdS) spacetimes, some metric components diverge as one
approaches the boundary, and their expansion in powers of $u$ is useful to
redefine the variables in terms of new, finite ones. The second reason is that
it allows us to understand which boundary conditions to impose on the ODEs
(10).

For this purpose, we start with an ansatz that is compatible with the AAdS
condition,111This ansatz must be modified if (6) does not hold

$\displaystyle A(t,u,x,y)$
$\displaystyle=\frac{1}{u^{2}}+\sum_{n=-1}^{\infty}a_{(2n)}(t,x,y)u^{n}\,,$
$\displaystyle B_{1}(t,u,x,y)$
$\displaystyle=\sum_{n=1}^{\infty}b_{1n}(t,x,y)u^{n}\,,$ (12) $\displaystyle
B_{2}(t,u,x,y)$ $\displaystyle=\sum_{n=1}^{\infty}b_{2n}(t,x,y)u^{n}\,,$
$\displaystyle G(t,u,x,y)$
$\displaystyle=\sum_{n=1}^{\infty}g_{n}(t,x,y)u^{n}\,,$ $\displaystyle
S(t,u,x,y)$
$\displaystyle=\frac{1}{u}+\sum_{n=0}^{\infty}s_{n}(t,x,y)u^{n}\,,$
$\displaystyle F_{x}(t,u,x,y)$
$\displaystyle=\sum_{n=0}^{\infty}f_{xn}(t,x,y)u^{n}\,,$ $\displaystyle
F_{y}(t,u,x,y)$ $\displaystyle=\sum_{n=0}^{\infty}f_{yn}(t,x,y)u^{n}\,,$
$\displaystyle\phi(t,u,x,y)$
$\displaystyle=\sum_{n=1}^{\infty}\phi_{n-1}(t,x,y)u^{n}\,.$

Substituting into equations (72-80) and solving order by order, we obtain

$\displaystyle A(t,u,x,y)$
$\displaystyle=\frac{1}{u^{2}}+\frac{2}{u}\xi+\xi^{2}-2\partial_{t}\xi-\frac{2\phi_{0}^{2}}{3}+u^{2}a_{4}$
$\displaystyle\quad{}-\frac{2}{3}u^{3}\left(3\xi
a_{4}+\partial_{x}f_{x2}+\partial_{y}f_{y2}+\phi_{0}\partial_{t}\phi_{2}\right)+\mathcal{O}\left(u^{4}\right),$
(13a) $\displaystyle B_{1}(t,u,x,y)$
$\displaystyle=u^{4}b_{14}+\mathcal{O}\left(u^{5}\right),$ (13b)
$\displaystyle B_{2}(t,u,x,y)$
$\displaystyle=u^{4}b_{24}+\mathcal{O}\left(u^{5}\right),$ (13c)
$\displaystyle G(t,u,x,y)$
$\displaystyle=u^{4}g_{4}+\mathcal{O}\left(u^{5}\right),$ (13d) $\displaystyle
S(t,u,x,y)$
$\displaystyle=\frac{1}{u}+\xi-\frac{\phi_{0}^{2}}{3}u+\frac{1}{3}\xi\phi_{0}^{2}u^{2}+\frac{1}{54}u^{3}\left(-18\xi^{2}\phi_{0}^{2}+\phi_{0}^{4}-18\phi_{0}\phi_{2}\right)$
$\displaystyle\quad{}+\frac{\phi_{0}}{90}u^{4}\left(30\xi^{3}\phi_{0}-5\xi\phi_{0}^{3}+90\xi\phi_{2}-24\partial_{t}\phi_{2}\right)+\mathcal{O}\left(u^{5}\right),$
(13e) $\displaystyle F_{x}(t,u,x,y)$
$\displaystyle=\partial_{x}\xi+u^{2}f_{x2}$
$\displaystyle\quad{}-\frac{2}{15}u^{3}\left(15\xi
f_{x2}+6\partial_{x}b_{14}+6\partial_{x}b_{24}-\partial_{y}g_{4}-2\phi_{0}\partial_{x}\phi_{2}\right)+\mathcal{O}\left(u^{4}\right),$
(13f) $\displaystyle F_{y}(t,u,x,y)$
$\displaystyle=\partial_{y}\xi+u^{2}f_{y2}$
$\displaystyle\quad{}-\frac{2}{15}u^{3}\left(15\xi
f_{y2}-6\partial_{y}b_{14}+6\partial_{y}b_{24}-\partial_{x}g_{4}-2\phi_{0}\partial_{y}\phi_{2}\right)+\mathcal{O}\left(u^{4}\right),$
(13g) $\displaystyle\phi(t,u,x,y)$
$\displaystyle=\phi_{0}u-\xi\phi_{0}u^{2}+u^{3}\left(\xi^{2}\phi_{0}+\phi_{2}\right)+u^{4}\left(\partial_{t}\phi_{2}-3\xi\phi_{2}-\xi^{3}\phi_{0}\right)+\mathcal{O}\left(u^{5}\right),$
(13h)

where $\phi_{2}$ is _not_ the one in (12), but redefined as

$\phi_{2}(t,x,y)\rightarrow\phi_{2}(t,x,y)+\xi^{2}(t,x,y)\phi_{0}.$ (14)

Note that $\phi_{0}$ is a constant, while the remaining variables in this
expansion are functions of $(t,x,y)$. In reality, the near boundary expansions
depend on $s_{0}$ instead of $\xi$. The fact that the former is simply shifted
by $\xi$ under (9) means that we can identify $s_{0}$ with $\xi$ and exchange
them everywhere.

We also need the expansions of “dotted” variables, defined in (71), which take
the form

$\displaystyle\dot{B_{1}}(t,u,x,y)$
$\displaystyle=-2b_{14}u^{3}+\mathcal{O}\left(u^{4}\right),$ (15a)
$\displaystyle\dot{B_{2}}(t,u,x,y)$
$\displaystyle=-2b_{24}u^{3}+\mathcal{O}\left(u^{4}\right),$ (15b)
$\displaystyle\dot{G}(t,u,x,y)$
$\displaystyle=-2g_{4}u^{3}+\mathcal{O}\left(u^{4}\right),$ (15c)
$\displaystyle\dot{S}(t,u,x,y)$
$\displaystyle=\frac{1}{2u^{2}}+\frac{\xi}{u}+\frac{\xi^{2}}{2}-\frac{\phi_{0}^{2}}{6}+\frac{1}{36}u^{2}\left(10a_{4}-5\phi_{0}^{4}+18\phi_{0}\phi_{2}\right)+\mathcal{O}\left(u^{3}\right),$
(15d) $\displaystyle\dot{F_{x}}(t,u,x,y)$
$\displaystyle=\partial_{t}\partial_{x}\xi-
uf_{x2}+\mathcal{O}\left(u^{2}\right),$ (15e)
$\displaystyle\dot{F_{y}}(t,u,x,y)$
$\displaystyle=\partial_{t}\partial_{y}\xi-
uf_{y2}+\mathcal{O}\left(u^{2}\right),$ (15f)
$\displaystyle\dot{\phi}(t,u,x,y)$
$\displaystyle=-\frac{\phi_{0}}{2}+u^{2}\left(\frac{\phi_{0}^{3}}{3}-\frac{3}{2}\phi_{2}\right)+\mathcal{O}\left(u^{3}\right).$
(15g)

The function $\xi(t,x,y)$ encodes our residual gauge freedom, and the
functions $a_{4}(t,x,y)$, $f_{x2}(t,x,y)$, $f_{y2}(t,x,y)$ are further
constrained to obey

$\displaystyle\partial_{t}a_{4}$
$\displaystyle=-\frac{4}{3}\left(\partial_{x}f_{x2}+\partial_{y}f_{y2}+\phi_{0}\partial_{t}\phi_{2}\right),$
(16a) $\displaystyle\partial_{t}f_{x2}$
$\displaystyle=-\frac{1}{4}\partial_{x}a_{4}-\partial_{x}b_{14}-\partial_{x}b_{24}+\partial_{y}g_{4}+\frac{1}{3}\phi_{0}\partial_{x}\phi_{2},$
(16b) $\displaystyle\partial_{t}f_{y2}$
$\displaystyle=-\frac{1}{4}\partial_{y}a_{4}+\partial_{y}b_{14}-\partial_{y}b_{24}+\partial_{x}g_{4}+\frac{1}{3}\phi_{0}\partial_{y}\phi_{2},$
(16c)

where $b_{14}(t,x,y)$, $b_{24}(t,x,y)$, $g_{4}(t,x,y)$, $\phi_{2}(t,x,y)$, and
$\partial_{t}\phi_{2}(t,x,y)$ are understood to be read off from the
asymptotic behaviour of $B_{1}(t,r,x,y)$, $B_{2}(t,r,x,y)$, $G(t,r,x,y)$, and
$\phi(t,r,x,y)$ in equations (13b), (13c), (13d) and (13h). The functions
$a_{4}(t_{0},x,y)$, $f_{x2}(t_{0},x,y)$, $f_{y2}(t_{0},x,y)$, and
$\xi(t_{0},x,y)$ should also be thought of as initial data, which can be
freely specified. $\phi_{0}$ is a parameter that must also be specified and
corresponds to the energy scale $\Lambda$ of the dual boundary theory.

####  2.1.3 Field redefinitions and boundary conditions

For the numerical implementation we find it useful to split the numerical grid
into two parts: the outer grid region (deep bulk) and the inner grid region
(close to the AdS boundary, where boundary conditions are imposed and the
gauge-theory variables are read off). As mentioned earlier, some of the metric
functions diverge at the AdS boundary while others vanish, being convenient to
make some redefinitions inspired by the asymptotic behaviour of these
functions so that the variables employed in the inner grid remain of order
unity therein. For the outer grid we choose to make simpler redefinitions,
which is helpful for the equation used to fix the gauge variable $\xi$.
Denoting with the $g1$ ($g2$) subscript the variables defined in the inner
(outer) grid, the redefinitions that we choose to make are then

$\displaystyle A(t,u,x,y)$
$\displaystyle=\frac{1}{u^{2}}+\frac{2}{u}\xi(t,x,y)+\xi^{2}(t,x,y)-2\partial_{t}\xi(t,x,y)-\frac{2\phi_{0}^{2}}{3}+u^{2}A_{g1}(t,u,x,y)$
$\displaystyle=-2\partial_{t}\xi(t,x,y)+A_{g2}(t,u,x,y)$ $\displaystyle
B_{1}(t,u,x,y)$ $\displaystyle=u^{4}B_{1g1}(t,u,x,y)$
$\displaystyle=B_{1g2}(t,u,x,y),$ $\displaystyle B_{2}(t,u,x,y)$
$\displaystyle=u^{4}B_{2g1}(t,u,x,y)$ $\displaystyle=B_{2g2}(t,u,x,y),$
$\displaystyle G(t,u,x,y)$ $\displaystyle=u^{4}G_{g1}(t,u,x,y)$
$\displaystyle=G_{g2}(t,u,x,y),$ $\displaystyle S(t,u,x,y)$
$\displaystyle=\frac{1}{u}+\xi(t,x,y)-\frac{\phi_{0}^{2}}{3}u+\frac{1}{3}\xi\phi_{0}^{2}u^{2}+u^{3}S_{g1}(t,u,x,y)$
$\displaystyle=S_{g2}(t,u,x,y),$ $\displaystyle F_{x}(t,u,x,y)$
$\displaystyle=\partial_{x}\xi(t,x,y)+u^{2}F_{xg1}(t,u,x,y)$
$\displaystyle=\partial_{x}\xi(t,x,y)+F_{xg2}(t,u,x,y),$ $\displaystyle
F_{y}(t,u,x,y)$ $\displaystyle=\partial_{y}\xi(t,x,y)+u^{2}F_{yg1}(t,u,x,y)$
$\displaystyle=\partial_{y}\xi(t,x,y)+F_{yg2}(t,u,x,y),$
$\displaystyle\phi(t,u,x,y)$
$\displaystyle=\phi_{0}u-\xi(t,x,y)\phi_{0}u^{2}+u^{3}\phi_{0}^{3}\phi_{g1}(t,u,x,y)$
$\displaystyle=\phi_{g2}(t,u,x,y),$ $\displaystyle\dot{B_{1}}(t,u,x,y)$
$\displaystyle=u^{3}\dot{B}_{1g1}(t,u,x,y)$
$\displaystyle=\dot{B}_{1g2}(t,u,x,y),$ $\displaystyle\dot{B_{2}}(t,u,x,y)$
$\displaystyle=u^{3}\dot{B}_{2g1}(t,u,x,y)$
$\displaystyle=\dot{B}_{2g2}(t,u,x,y),$ $\displaystyle\dot{G}(t,u,x,y)$
$\displaystyle=u^{3}\dot{G}_{g1}(t,u,x,y)$
$\displaystyle=\dot{G}_{g2}(t,u,x,y),$ $\displaystyle\dot{S}(t,u,x,y)$
$\displaystyle=\frac{1}{2u^{2}}+\frac{\xi(t,x,y)}{u}+\frac{\xi^{2}(t,x,y)}{2}-\frac{\phi_{0}^{2}}{6}+u^{2}\dot{S}_{g1}(t,u,x,y)$
$\displaystyle=\dot{S}_{g2}(t,u,x,y),$ $\displaystyle\dot{\phi}(t,u,x,y)$
$\displaystyle=-\frac{\phi_{0}}{2}+u^{2}\phi_{0}^{3}\dot{\phi}_{g1}(t,u,x,y)$
$\displaystyle=\dot{\phi}_{g2}(t,u,x,y).$

Substituting these redefined variables into the system of equations (72-80),
we are left with two new versions of this system, one for the near boundary
region (inner grid), and the other one for the bulk region (outer grid). The
corresponding ODEs can then be integrated in the inner grid ($g1$) by imposing
the following boundary conditions

$\displaystyle S_{g1}|_{u=0}$
$\displaystyle=\frac{1}{54}\left(-18\xi^{2}\phi_{0}^{2}+\phi_{0}^{4}-18\phi_{0}\phi_{2}\right),$
(17a) $\displaystyle\partial_{u}S_{g1}|_{u=0}$
$\displaystyle=\frac{\phi_{0}}{90}\left(30\xi^{3}\phi_{0}-5\xi\phi_{0}^{3}+90\xi\phi_{2}-24\partial_{t}\phi_{2}\right),$
(17b) $\displaystyle F_{xg1}|_{u=0}$ $\displaystyle=f_{x2},$ (17c)
$\displaystyle\partial_{u}F_{xg1}|_{u=0}$
$\displaystyle=-\frac{2}{15}\left(15\xi
f_{x2}+6\partial_{x}b_{14}+6\partial_{x}b_{24}-\partial_{y}g_{4}-2\phi_{0}\partial_{x}\phi_{2}\right),$
(17d) $\displaystyle F_{yg1}|_{u=0}$ $\displaystyle=f_{y2},$ (17e)
$\displaystyle\partial_{u}F_{yg1}|_{u=0}$
$\displaystyle=-\frac{2}{15}\left(15\xi
f_{y2}+6\partial_{y}b_{14}+6\partial_{y}b_{24}-\partial_{x}g_{4}-2\phi_{0}\partial_{y}\phi_{2}\right),$
(17f) $\displaystyle\dot{S}_{g1}|_{u=0}$
$\displaystyle=\frac{1}{36}\left(18a_{4}-5\phi_{0}^{4}+18\phi_{0}\phi_{2}\right),$
(17g) $\displaystyle\dot{B}_{1g1}|_{u=0}$ $\displaystyle=-2b_{14},$ (17h)
$\displaystyle\dot{B}_{2g1}|_{u=0}$ $\displaystyle=-2b_{24},$ (17i)
$\displaystyle\dot{G}_{g1}|_{u=0}$ $\displaystyle=-2g_{4},$ (17j)
$\displaystyle\dot{\phi}_{g1}|_{u=0}$
$\displaystyle=\frac{1}{3}-\frac{3\phi_{2}}{2\phi_{0}^{3}},$ (17k)
$\displaystyle A_{g1}|_{u=0}$ $\displaystyle=a_{4},$ (17l)
$\displaystyle\partial_{u}A_{g1}|_{u=0}$ $\displaystyle=-\frac{2}{3}\left(3\xi
a_{4}+\partial_{x}f_{x2}+\partial_{y}f_{y2}+\phi_{0}\partial_{t}\phi_{2}\right).$
(17m)

Once again we note that functions $B_{1}$, $B_{2}$, $G$, $\phi$, $a_{4}$,
$f_{x2}$, $f_{y2}$ and $\xi$ encode the freely-specifiable data. Once the
inner grid ODEs have been solved, we evaluate each function at the interface
with the outer grid to obtain the boundary conditions for the $g2$ variables
and integrate the corresponding equations.

####  2.1.4 Gauge fixing

To fully close our system we still need to fix the residual gauge freedom (9).
It is advantageous for the numerical implementation to have the Apparent
Horizon (AH) lie at constant radial slice $r=r_{H}$ at all times, so it will
be convenient to fix a gauge that enforces this throughout the numerical
evolution. We thus want to guarantee that $\Theta|_{r=r_{H}}=0$ at all times,
where $\Theta$ is the expansion of outgoing null rays. Its explicit expression
for the metric (7) is shown in Appendix C.

A simple way to enforce $\Theta|_{r=r_{H}}=0$ at all times during the
numerical evolution is to impose a diffusion-like equation of the form

$\left(\partial_{t}\Theta+\kappa\Theta\right)|_{u=u_{H}}=0$ (18)

with $\kappa>0$, ensuring that the expansion $\Theta$ is driven towards the
fix point $\Theta|_{u=u_{H}}=0$ as the time evolution runs, pushing the AH
surface to $u=u_{H}=\mathrm{constant}$.

The way to proceed is the following. We expand equation (18) using (89) and
also the equations of motion for both $\ddot{S}$ and $\dot{F}_{x,y}$. Then we
substitute all the variables by the outer grid redefinitions, $g2$, and
evaluate them at $u=u_{H}$. We obtain a linear PDE for $\partial_{t}\xi$ of
the type

$\left(A^{(\xi)}_{xx}\partial^{2}_{x}+A^{(\xi)}_{xy}\partial_{x}\partial_{y}+A^{(\xi)}_{yy}\partial^{2}_{y}+B^{(\xi)}_{x}\partial_{x}+B^{(\xi)}_{y}\partial_{y}+C^{(\xi)}\right)\partial_{t}\xi(t,x,y)=-S^{(\xi)}\,,$
(19)

which can be readily integrated with periodic boundary conditions in $x$ and
$y$.

####  2.1.5 Evolution algorithm

Having solved equations (72-80), we use the definition of the “dot” operator,
cf. equation (71), to write

$\displaystyle\partial_{t}B_{1}(t,u,x,y)=\dot{B}_{1}(t,u,x,y)+\frac{u^{2}}{2}A(t,u,x,y)\partial_{u}B_{1}(t,u,x,y)\,,$
(20)

and analogously for $B_{2}$, $G$ and $\phi$. This tells us how to march these
quantities forward in time.222In practice we write explicitly the evolution
equations in terms of the redefined $g1$ and $g2$ functions.

As outlined in the previous subsections, we decompose our computational grid
(in the $u$-direction) into two domains: an inner (near boundary) domain and
an outer (bulk) domain. The outer domain can further be split into subdomains.
We therefore need to match the evolution variables across these domains. The
procedure is outlined in Appendix A of Attems:2017zam  which, for
completeness, we here summarize.

The evolution equation for $B_{1}$ (the case for the remaining evolution
variables is analogous) has the generic form

$\partial_{t}B_{1}(t,u,x,y)=c(t,u,x,y)\partial_{u}B_{1}(t,u,x,y)+F_{B_{1}}(t,u,x,y)\,,$
(21)

with

$c(t,u,x,y)=\frac{u^{2}}{2}A(t,u,x,y)\,.$ (22)

$c(t,u,x,y)$ is locally the propagation speed, and in the vicinity of some
$u=u_{0}$ lying at the interface between two domains $i$ and $i+1$ we can
formally write the solution of this equation (ignoring from now on the $x,y$
dependence) as

$B_{1}(t,u_{0})\simeq h(u_{0}+c\,t)+\int F_{B_{1}}$

for any given function $h$.

Therefore, for $c>0$ ($c<0$), information is propagating from domain $i+1$ to
domain $i$ (domain $i$ to domain $i+1$). In order to consistently solve this
system, the procedure we employ is to use equation (21) (and corresponding
ones for the remaining domains) on all interior points; at the junction point
$u=u_{0}$ we check the propagation speed $c$ at each $x,y$ point and copy the
values according to the propagation direction at the interface junction:

  * •

$c>0$

$\displaystyle\partial_{t}B_{1}^{(i+1)}|_{u=u_{0}}$
$\displaystyle=c(u_{0})\partial_{u}B_{1}^{(i+1)}|_{u=u_{0}}+F_{B_{1}}(u_{0})\,,$
(23) $\displaystyle\partial_{t}B_{1}^{(i)}|_{u=u_{0}}$
$\displaystyle=\partial_{t}B_{1}^{(i+1)}|_{u=u_{0}}\,,$

i.e., we copy the modes leaving domain $i+1$ to domain $i$.

  * •

$c<0$

$\displaystyle\partial_{t}B_{1}^{(i)}|_{u=u_{0}}$
$\displaystyle=c(u_{0})\partial_{u}B_{1}^{(i)}|_{u=u_{0}}+F_{B_{1}}(u_{0})\,,$
(24) $\displaystyle\partial_{t}B_{1}^{(i+1)}|_{u=u_{0}}$
$\displaystyle=\partial_{t}B_{1}^{(i)}|_{u=u_{0}}\,,$

i.e., we copy the modes leaving domain $i$ to domain $i+1$.

We can now schematically outline the evolution algorithm, which is as follows.

  1. 1.

Initial conditions $B_{1}(t_{0},u,x,y)$, $B_{2}(t_{0},u,x,y)$,
$G(t_{0},u,x,y)$, $\phi(t_{0},u,x,y)$, $a_{4}(t_{0},x,y)$,
$f_{x2}(t_{0},x,y)$, $f_{y2}(t_{0},x,y)$ and $\xi(t_{0},x,y)$ are provided for
some initial time $t_{0}$.

  2. 2.

Equations (72-80) are solved in succession for the redefined variables in the
inner grid $g1$, imposing the boundary conditions (17), and then the same
equations are solved for the outer grids $g2$, forcing the variables to match
their values at grid interfaces.

  3. 3.

Equation (19) is solved to find $\partial_{t}\xi(t_{0},x,y)$. Expression (20)
is then used to evaluate $\partial_{t}B_{1}(t_{0},u,x,y)$,
$\partial_{t}B_{2}(t_{0},u,x,y)$, $\partial_{t}G(t_{0},u,x,y)$,
$\partial_{t}\phi(t_{0},u,x,y)$. These variables communicate at domain
interfaces through equations (23) and (24).

  4. 4.

Obtain $\partial_{t}a_{4}(t_{0},x,y)$, $\partial_{t}f_{x2}(t_{0},x,y)$ and
$\partial_{t}f_{y2}(t_{0},x,y)$ through (16).

  5. 5.

Advance $B_{1}$, $B_{2}$, $G$, $\phi$, $a_{4}$, $f_{x2}$, $f_{y2}$ and $\xi$
to time $t_{1}$.

See Fig. 1 for a cartoon picture of the coordinates used and the evolution
scheme (at constant $x,y$).

Figure 1: Penrose diagram of the evolution procedure, at constant $x,y$
slices. The shaded region represents the region covered by the computational
domain.

####  2.1.6 Gauge theory expectation values

The gauge theory expectation values can be obtained from the asymptotic
behaviour of the bulk variables in a way similar to Attems:2017zam . The
result is:

$\displaystyle\mathcal{E}$
$\displaystyle=\tfrac{\kappa_{5}^{2}}{2L^{3}}\,\langle T^{tt}\rangle$
$\displaystyle=-\frac{3}{4}a_{4}-\phi_{0}\phi_{2}+\left(\frac{7}{36}-\lambda_{4}\right)\phi_{0}^{4},$
(25) $\displaystyle\mathcal{P}_{x}$
$\displaystyle=\tfrac{\kappa_{5}^{2}}{2L^{3}}\,\langle T^{xx}\rangle$
$\displaystyle=-\frac{a_{4}}{4}-b_{14}-b_{24}+\frac{\phi_{0}\phi_{2}}{3}+\left(\frac{-5}{108}+\lambda_{4}\right)\phi_{0}^{4},$
$\displaystyle\mathcal{P}_{xy}$
$\displaystyle=\tfrac{\kappa_{5}^{2}}{2L^{3}}\,\langle T^{xy}\rangle$
$\displaystyle=-g_{4},$ $\displaystyle\mathcal{P}_{y}$
$\displaystyle=\tfrac{\kappa_{5}^{2}}{2L^{3}}\,\langle T^{yy}\rangle$
$\displaystyle=-\frac{a_{4}}{4}+b_{14}-b_{24}+\frac{\phi_{0}\phi_{2}}{3}+\left(\frac{-5}{108}+\lambda_{4}\right)\phi_{0}^{4},$
$\displaystyle\mathcal{P}_{z}$
$\displaystyle=\tfrac{\kappa_{5}^{2}}{2L^{3}}\,\langle T^{zz}\rangle$
$\displaystyle=-\frac{a_{4}}{4}+2b_{24}+\frac{\phi_{0}\phi_{2}}{3}+\left(\frac{-5}{108}+\lambda_{4}\right)\phi_{0}^{4},$
$\displaystyle\mathcal{J}_{x}$
$\displaystyle=-\tfrac{\kappa_{5}^{2}}{2L^{3}}\,\langle T^{tx}\rangle$
$\displaystyle=f_{x2},$ $\displaystyle\mathcal{J}_{y}$
$\displaystyle=-\tfrac{\kappa_{5}^{2}}{2L^{3}}\,\langle T^{ty}\rangle$
$\displaystyle=f_{y2},$ $\displaystyle\mathcal{V}$
$\displaystyle=\tfrac{\kappa_{5}^{2}}{2L^{3}}\,\langle\mathcal{O}_{\phi}\rangle$
$\displaystyle=-2\phi_{2}+\left(\frac{1}{3}-4\lambda_{4}\right)\phi_{0}^{3}.$

For an $SU(N)$ gauge theory the prefactor $\kappa_{5}^{2}/2L^{3}$ in these
equations typically scales as $N^{-2}$, whereas the stress tensor scales as
$N^{2}$. The rescaled quantities are therefore finite in the large-$N$ limit.
The stress tensor and the expectation of the scalar operator are related
through the Ward identity

$\langle T^{\mu}_{\mu}\rangle=-\Lambda\langle\mathcal{O}\rangle\,.$ (26)

###  2.2 Implementation

As already mentioned, we have implemented the algorithm of Sec. 2.1.5 in a new
numerical code called Jecco jecco-2022 , written in Julia Julia-2017 . Julia
is a dynamically-typed language with good support for interactive use and with
runtime performance approaching that of statically-typed languages such as C
or Fortran. Even though a relative newcomer to the field of scientific
computing, its popularity has been steadily growing in the last few years. It
boasts a friendly community of users and developers and a rapidly growing
package ecosystem.

Jecco was developed as a Julia module and is freely available at
https://github.com/mzilhao/Jecco.jl. This code is a generalization of the 2+1
C code introduced in Attems:2017zam , and completely written from scratch. The
codebase is neatly divided into generic infrastructure, such as general
derivative operators, filters, and input/output routines (which are defined in
the main Jecco module) and physics, such as initial data, evolution equations,
and diagnostic routines (which are defined in submodules).

In Jecco we have implemented finite-difference operators of arbitrary order
through the Fornberg algorithm fornberg1998classroom  as well as Chebyshev and
Fourier differentiation matrices. These methods are completely general and can
be used with any Julia multidimensional array. We have also implemented output
methods that roughly follow the openPMD standard openPMD-2018  for writing
data.

####  2.2.1 Discretization

For our numerical implementation of the algorithm in Sec. 2.1.5 we have
discretized the $x$ and $y$ directions on uniform grids where periodic
boundary conditions are imposed, while along the $u$ direction we break the
computational domain into several (touching) subdomains with $N_{u}$ points.
In each subdomain a Lobatto-Chebyshev grid is used where the collocation
points, given by

$X_{i+1}=-\cos\left(\frac{\pi\,i}{N_{u}}\right)\qquad(i=0,1,\ldots,N_{u}-1)\,,$
(27)

are defined in the range $[-1:+1]$, and can be mapped to the physical grid by

$u_{i}=\frac{u_{R}+u_{L}}{2}+\frac{u_{R}-u_{L}}{2}X_{i}\qquad(i=1,\ldots,N_{u})\,,$
(28)

where $u_{L}$ and $u_{R}$ are the limits of each of subdomain. For the
subdomain that includes the AdS boundary ($u=0$), the inner grid variables of
Sec. 2.1.3 are used; all remaining subdomains use the outer grid variables.

Derivatives along the $x$ and $y$ directions are approximated by (central)
finite differences. Although in Jecco operators of arbitrary order are
available, we have mostly made use of fourth-order accurate ones for our
applications. In the radial direction $u$, the use of the Chebyshev-Lobatto
grid allow us to use pseudo-spectral collocation methods Boyd2001 . These
methods are based in approximating solutions in a basis of Chebyshev
polynomials $T_{n}(X)$ but, in addition to the spectral basis, we have an
additional physical representation – the values that functions take on each
grid point – and therefore we can perform operations in one basis or the other
depending on our needs. Discretization using the pseudo-spectral method
consists in the exact imposition of our equations at the collocation points of
the Chebyshev-Lobatto grid.

The radial equations that determine our grid functions have the schematic form
of equation (10), where $f$ represents the metric coefficients and scalar
field $\phi$. Once our coordinate $u$ is discretized, the differential
operator becomes an algebraic one acting over the values of the functions in
the collocation points taking the form (at every point in the transverse
directions $x,y$)

$\sum_{j=1}^{N_{u}}\left[A_{f}^{i}(t,x,y)\mathcal{D}_{uu}^{ij}+B_{f}^{i}(t,x,y)\mathcal{D}_{u}^{ij}+C_{f}^{i}(t,x,y)\mathbb{I}^{ij}\right]f^{j}(t,x,y)=-S_{f}^{i}(t,x,y)$
(29)

(no sum in $i$), where $\mathcal{D}_{uu}$, $\mathcal{D}_{u}$ represent the
derivative operators for a Chebyshev-Lobatto grid in the physical
representation (see for instance trefethen2000spectral  for the explicit
expression) and $i$, $j$ indices in the $u$ coordinate. Boundary conditions
are imposed by replacing full rows in this operator by the values we need to
fix: at the inner grid $g1$, we impose the boundary conditions in (17); at the
outer grids these are read off from the obtained values in the previous
subdomain.

The resulting operators are then factorized through an LU decomposition and
the linear systems (29) are subsequently solved using Julia’s left division
(ldiv!) operation. Recall that we need to solve one such radial equation per
grid point in the $x,y$ transverse directions. Since these equations are
independent of each other, we can trivially parallelize the procedure using
Julia’s Threads.@threads macro.

Equation (19) for $\partial_{t}\xi$ is a linear PDE in $x,y$. To solve it,
after discretizing in a $N_{x}\times N_{y}$ grid, we flatten the solution
vector using lexicographic ordering

$\bm{g}\equiv\begin{pmatrix}\partial_{t}\xi(t,x_{1},y_{1})\\\
\partial_{t}\xi(t,x_{2},y_{1})\\\ \vdots\\\
\partial_{t}\xi(t,x_{N_{x}},y_{1})\\\ \partial_{t}\xi(t,x_{1},y_{2})\\\
\vdots\\\ \partial_{t}\xi(t,x_{N_{x}},y_{N_{y}})\end{pmatrix}$

and introduce enlarged differentiation matrices, which can be conveniently
built as Kronecker products

$\displaystyle\hat{\mathcal{D}}_{x}$ $\displaystyle=\mathbb{I}_{N_{y}\times
N_{y}}\otimes\mathcal{D}_{x},$ $\displaystyle\qquad\hat{\mathcal{D}}_{y}$
$\displaystyle=\mathcal{D}_{y}\otimes\mathbb{I}_{N_{x}\times N_{x}},$ (30)
$\displaystyle\hat{\mathcal{D}}_{xx}$ $\displaystyle=\mathbb{I}_{N_{y}\times
N_{y}}\otimes\mathcal{D}_{xx},$ $\displaystyle\qquad\hat{\mathcal{D}}_{yy}$
$\displaystyle=\mathcal{D}_{yy}\otimes\mathbb{I}_{N_{x}\times N_{x}},$

where $\mathcal{D}_{x}$, $\mathcal{D}_{y}$, $\mathcal{D}_{xx}$,
$\mathcal{D}_{yy}$ are the first and second derivative finite-difference
operators. The cross derivative operator is built as a matrix product,
$\hat{\mathcal{D}}_{xy}=\hat{\mathcal{D}}_{x}\hat{\mathcal{D}}_{y}$. The PDE
(19) then takes the algebraic form

$\sum_{J=1}^{N_{x}\times
N_{y}}\left[A_{xx}^{I}\hat{\mathcal{D}}_{xx}^{IJ}+A_{xy}^{I}\hat{\mathcal{D}}_{xy}^{IJ}+A_{yy}^{I}\hat{\mathcal{D}}_{yy}^{IJ}+B_{x}^{I}\hat{\mathcal{D}}_{x}^{IJ}+B_{y}^{I}\hat{\mathcal{D}}_{y}^{IJ}+C^{I}\mathbb{I}^{IJ}\right]\bm{g}^{J}=-S_{\bm{g}}^{I}$
(31)

(no sum in $I$), where $I,J=1,\ldots,N_{x}\times N_{y}$. The $x$ and $y$
directions are periodic, so no boundary conditions need to be imposed. See for
example Krikun2018  for a pedagogical overview of these techniques.

As before, the operator defined inside the square brackets is factorized
through an LU decomposition and the linear system (31) is then solved with the
left division operation. Since all the matrices are sparse, we store them in
the Compressed Sparse Column format using the type SparseMatrixCSC.

####  2.2.2 Time evolution

For the time evolution we use a method of lines procedure, where we find it
convenient to pack all evolved variables (across all subdomains) into one
single state vector. This state vector is then marched forwarded in time with
the procedure of Sec. 2.1.5 using the ODEProblem interface from the
DifferentialEquations.jl Julia package DifferentialEquations.jl-2017 . This
package provides a very long and complete list of integration methods. For our
applications, since evaluating the time derivative of our state vector is an
expensive operation, we find it convenient for reasons of speed and accuracy
to use the Adams-Bashforth and Adams-Moulton family of multistep methods.
Depending on the application, we find that the (third order) fixed step method
AB3 and the adaptive step size ones VCAB3 and VCABM3 seem to work particularly
well. The integration package automatically takes care of the starting values
by using a lower-order method initially.

We use Kreiss-Oliger dissipation Kreiss1973  to remove spurious high-frequency
noise common to finite-difference schemes. In particular, when using finite-
difference operators of order $p-1$, we add Kreiss-Oliger dissipation of order
$p$ to all evolved quantities $f$ as

$f\leftarrow
f+\sigma\frac{(-1)^{(p+3)/2}}{2^{p+1}}\left(h_{x}^{p+1}\frac{\partial^{(p+1)}}{\partial
x^{(p+1)}}+h_{y}^{p+1}\frac{\partial^{(p+1)}}{\partial y^{(p+1)}}\right)f$
(32)

after each time step, where $h_{x}$ and $h_{y}$ are the grid spacings and
$\sigma$ is a tuneable dissipation parameter which we typically set to 0.2
unless explicitly stated otherwise. This procedure effectively works as a low-
pass filter.

Along the $u$-direction we can damp high order modes directly in the spectral
representation. After each time step, we apply an exponential filter to the
spectral coefficients of our $u$-dependent evolved quantities $f$ (see for
instance KANEVSKY200641 ). The complete scheme is

$\left\\{f_{i}\ \right\\}\stackrel{{\scriptstyle\rm
FFT}}{{\longrightarrow}}\left\\{\hat{f}_{k}\right\\}\stackrel{{\scriptstyle}}{{\longrightarrow}}\left\\{\hat{f}_{k}\;e^{-\alpha(k/M)^{\gamma
M}}\right\\}\stackrel{{\scriptstyle\rm
FFT}}{{\longrightarrow}}\left\\{f_{i}\right\\}$ (33)

where $M\equiv N_{u}-1$, $k=0,\ldots,M$, $\alpha=\log\epsilon$ where
$\epsilon$ is the machine epsilon (for the standard choice of
$\epsilon=2^{-52}$, $\alpha=36.0437$) and $\gamma$ is a tuneable parameter
which we typically fix to $\gamma=8$. This effectively dampens the
coefficients of the higher-order Chebyshev polynomials.

We performed a thorough set of tests on this implementation, which is detailed
in Appendix A.

##  3 Bubble dynamics

The Jecco code described in the previous section was first applied to the
study of gravitational waves produced by the spinodal instability in a
cosmological first-order phase transition Bea:2021zol . We now turn to a new
application, namely the dynamics of bubbles in a strongly-coupled, four-
dimensional gauge theory. For this purpose we will focus on a holographic
model of the type described by equations (1) and (4) with the same value of
the parameters (5) as in Bea:2021zsu , namely

$\phi_{M}=0.85\,,\qquad\phi_{Q}=10\,.\iff\lambda_{4}=-0.346021\,,\qquad\lambda_{6}=0.1\,.$
(34)

The motivation for the general class of models under consideration is that
they provide simple examples of non-conformal theories with first-order phase
transitions (for appropriate values of $\phi_{M}$ and $\phi_{Q}$) whose dual
gravity solutions are completely regular even at zero temperature. The
motivation for the choice (34) is that it leads to a sizeable bubble wall
velocity, as we will see in Sec. 3.4.

###  3.1 Thermodynamics

The thermodynamics of the gauge theory can be extracted from the homogeneous
black brane solutions on the gravity side (see e.g. Gubser:2008ny ). Figure 2
shows the result for the energy density as a function of temperature, where we
see the usual multivaluedness associated to a first-order phase transition.

Figure 2: Energy density as a function of temperature for the gauge theory
dual to the holographic model (1)-(4) with parameters (34). The squares
$B_{c}$ and $A_{c}$ correspond, respectively, to the states inside and outside
of the closest-to-critical bubble studied in Sec. 3.3. The dots $B$ and $C$
correspond to the initial states inside and outside the expanding bubble
studied in Sec. 3.4, respectively. At late times, the state $B$ inside the
bubble evolves into $C$, and a heated region is created in front of the bubble
that can be characterized in terms of the point $D$ in the phase diagram.

At high and low temperatures there is only one phase available to the system.
Each of these phases is represented by a solid, blue curve. At the critical
temperature

$T_{c}=0.418\Lambda$

the state that minimizes the free energy moves from one branch to the other.
The first-order nature of the transition is encoded in the non-zero latent
heat, namely in the discontinuous jump in the energy density given by

$\mathcal{E}_{\mathrm{latent}}=\mathcal{E}_{\mathrm{high}}-\mathcal{E}_{\mathrm{low}}\,,\qquad\mathcal{E}_{\mathrm{low}}=0.225\Lambda^{4}\,,\qquad\mathcal{E}_{\mathrm{high}}=2.123\Lambda^{4}\,.$
(35)

Note that the phase transition is a transition between two deconfined plasma
phases, since both phases have energy densities of order $N^{2}$ and they are
both represented by a black brane geometry with a horizon on the gravity side.

In a region

$T_{s}=0.3879\Lambda<T<T_{s}^{\prime}=0.4057\Lambda$ (36)

around the critical temperature there are three different states available to
the system for a given temperature. The thermodynamically preferred one is the
state that minimizes the free energy, namely a state on one of the blue
curves. The states on the dashed, brown curves are not globally preferred but
they are locally thermodynamically stable, i.e. they are metastable. This
follows from the fact that specific heat

$c_{v}\equiv\frac{d\mathcal{E}}{dT}$ (37)

is positive on these branches. At the temperatures $T_{s}$ and
$T_{s}^{\prime}$ the metastable curves meet the dotted-dashed, red curve,
known as the “spinodal branch”. States on this branch are locally unstable
since their specific heat is negative and have energies comprised between

$\mathcal{E}^{\prime}_{s}=0.717\Lambda^{4}\,,\qquad\mathcal{E}_{s}=1.141\Lambda^{4}\,.$
(38)

Note that the characteristic scale for all the quantities above is set by the
microscopic scale in the gauge theory, $\Lambda$, given holographically by
$\Lambda=\phi_{0}$ in terms of the leading term in the near-boundary fall-off
of the scalar field in (13h).

###  3.2 Initial data

As any other thermal system with a first order phase transition, the gauge
theory can be overcooled past the critical temperature $T_{c}$. The
homogeneous, overcooled state, represented by a point on the upper, brown
branch in Fig. 2, is stable against small fluctuations, including thermal
ones, but not against sufficiently large fluctuations. A particular class of
large fluctuations are bubbles, namely inhomogeneous configurations in which
the energy density of a certain region of space within the overcooled
homogeneous phase is reduced. For sufficiently large bubbles, the energy
density in the centre of this region lies in the stable branch of the phase
diagram, represented by the lower, blue curve in Fig. 2, and the bubble
smoothly interpolates between the stable and the metastable phases.

In a homogeneous and isotropic thermal system it is expected that the
nucleated bubbles are spherical. However, given our symmetry restrictions we
will study cylindrical bubbles. This is enough to bring about two new physical
aspects compared to our previous work Bea:2021zsu  for planar configurations.
The first one is that the surface tension now plays a role. In particular, we
will be able to identify a critical bubble in which the inward-pointing force
due to the surface tension exactly balances the outward-pointing force coming
from the pressure difference between the inside and the outside of the bubble.
The second one is that the asymptotic profile of an expanding bubble possesses
more structure than in the planar case.

Our first task is to construct initial data corresponding to a bubble. By
definition, this is a configuration consisting of a cylindrical region filled
with the stable phase (the inside of the bubble) connected to an asymptotic
region filled with the metastable phase (the outside of the bubble) through an
appropriate interface. The stable and metastable phases correspond to the
points labelled $B$ and $A$ in Fig. 2, respectively, and both have $T<T_{c}$.
As we will now explain, our strategy to construct these bubbles will be to
start with a phase-separated state, which has $T=T_{c}$, and to rescale it
appropriately.

Figure 3: Phase-separated configurations in a box of size
$L_{x}\Lambda=L_{y}\Lambda=20$ with average energy densities
$\bar{\mathcal{E}}/\Lambda^{4}=1.0$ (top), $\bar{\mathcal{E}}/\Lambda^{4}=1.6$
(middle) and $\bar{\mathcal{E}}/\Lambda^{4}=1.8$ (bottom).

Phase-separated states are configurations in which the two homogeneous phases
with energy densities $\mathcal{E}_{\mathrm{high}}$ and
$\mathcal{E}_{\mathrm{low}}$ coexist in equilibrium at $T=T_{c}$. This is
possible because at this temperature the free energy densities, and hence the
pressures, are equal in the two phases. Three examples of such configurations
in a box of constant size are shown in Fig. 3. The difference between the
three cases is the relative fraction of the total volume occupied by each
phase. For a box of fixed size, changing this relative fraction is equivalent
to changing the average energy density in the box, $\bar{\mathcal{E}}$. The
larger the average energy density, the larger the size of the high-energy
region, and vice-versa. We will use this fact to our advantage when we search
for the critical bubble below.

Strictly speaking, phase-separated states only exist in infinite volume, since
only in that case the two coexisting phases become arbitrarily close to being
homogeneous sufficiently far away from the interface. The middle and bottom
panels of Fig. 3 correspond to states that are fairly close to this limit, but
deviations can still be seen with the naked eye. For example, the energy
density in the region outside the bubbles is slightly below $2\Lambda^{4}$,
whereas the energy density in the high-energy phase at $T=T_{c}$ has
$\mathcal{E}_{\mathrm{high}}$ above $2\Lambda^{4}$, as given in (35). The
state in Fig. 3(top) is even more affected by finite box-size effects because
the size of the low-energy region is comparable to the size of the box. In any
case, these deviations will have no implications for our purposes, since we
are not interested in phase-separated states per se but only in using them to
construct initial data for bubble configurations.

The value of $\bar{\mathcal{E}}$ in a box of fixed size is conserved upon time
evolution. Therefore, phase-separated states with an average energy density in
the region $\mathcal{E}_{s}^{\prime}\leq\bar{\mathcal{E}}\leq\mathcal{E}_{s}$
can be generated by starting with a homogeneous state in the spinodal region
of Fig. 2, perturbing it slightly, and letting evolve until it settles down to
a phase-separated configuration Attems:2019yqn ; Bea:2021zol . To initialize
the code we specify some $\phi_{2}$ that is not too far away from the value of
the thermal state and generate a simple bulk profile for the scalar,
$\phi(t=0,u)$, given by the truncated series in (13h) to third order. This is
not the geometry associated to the black brane of such energy density, but it
would relax fast to the true static solution. The value for $a_{4}$ is
obtained by using the energy expression in (25). On top of it we add a
sinusoidal perturbation, so the final $a_{4}$ reads

$a_{4}(t=0,x,y)=\bar{a}_{4}\left[1+\delta
a_{4}\left(\cos\left(\frac{2\pi}{L_{x}}(x-x_{\rm
mid})\right)+\cos\left(\frac{2\pi}{L_{y}}(y-y_{\rm
mid})\right)\right)\right],$ (39)

where $\bar{a}_{4}$ is the value we obtained above, $L_{x}$ and $L_{y}$ are
the lengths of the box, $x_{\rm mid}$ and $y_{\rm mid}$ correspond to the
central point and $\delta a_{4}$ represents the amplitude of the perturbation,
equal for both $x$\- and $y$-directions. The fastest way to arrive at a phase-
separated configuration is to assign the largest possible value to $\delta
a_{4}$ compatible with keeping the apparent horizon within our grid. We have
found that $\delta a_{4}\sim 10^{-3}$ is a convenient choice. The state in
Fig. 3(top) was generated following this method with
$\phi_{2}=0.3\Lambda^{3}$. After a time $t\Lambda=300$ the system has settled
down to the configuration shown in the figure.

Phase-separated configurations with average energy densities in the regions

$\mathcal{E}_{s}\leq\bar{\mathcal{E}}<\mathcal{E}_{\mathrm{high}}\qquad\mbox{and}\qquad\mathcal{E}_{\mathrm{low}}<\bar{\mathcal{E}}\leq\mathcal{E}_{s}^{\prime}$
(40)

also exist, but they cannot be found directly via time evolution of an initial
state in the spinodal region. Instead, to obtain them we follow Bea:2020ees .
We take initial data corresponding to a phase-separated state with
$\bar{\mathcal{E}}$ in the spinodal region, and we modify it by increasing or
decreasing the value of $\bar{a}_{4}$ so that the new $\bar{\mathcal{E}}$
takes the desired value. We then let the system evolve. In a time around
$t=100/\Lambda$ the system relaxes to a new inhomogeneous, static
configuration. The phase-separated configurations in the middle and bottom
panels of Fig. 3 have
$\mathcal{E}_{s}\leq\bar{\mathcal{E}}\leq\mathcal{E}_{\mathrm{high}}$ and were
obtained with this procedure.

The phase-separated states interpolate between the energy densities
$\mathcal{E}_{\mathrm{low}}$ and $\mathcal{E}_{\mathrm{high}}$. To construct
initial data for bubble configurations that interpolate between two energy
densities $\mathcal{E}_{B}$ and $\mathcal{E}_{A}$ we proceed as follows. Let
$f_{\textrm{PS}}$ be any of the functions specifying the initial data of a
phase-separated state. This could be one of the metric components in the bulk
or the scalar field, in which case $f_{\textrm{PS}}=f_{\textrm{PS}}(u,x,y)$,
or one of the boundary functions such as $a_{4}$, in which case
$f_{\textrm{PS}}=f_{\textrm{PS}}(x,y)$. We assume that the centre of the
region with energy density $\mathcal{E}_{\mathrm{low}}$ is at $x=y=0$, and
that the point at the edge of the box $x=y=L/2$ lies in the region with energy
density $\mathcal{E}_{\mathrm{high}}$. Let $f_{A}$ and $f_{B}$ be the
corresponding functions for the states $A$ and $B$. Since these states are
homogeneous, $f_{A}$ and $f_{B}$ depend on $u$ for a bulk function and are
just constants for a boundary function. We then define the corresponding
initial data for a bubble through the rescaling

$f_{\textrm{bubble}}(u,x,y)=f_{B}(u)+\Big{(}f_{A}(u)-f_{B}(u)\Big{)}\left[\frac{f_{\textrm{PS}}(u,x,y)-f_{\textrm{PS}}(u,0,0)}{f_{\textrm{PS}}(u,L/2,L/2)-f_{\textrm{PS}}(u,0,0)}\right]\,.$
(41)

If $f$ is a boundary function then there the dependence on $u$ is absent. At
any fixed value of $u$, the term in square brackets interpolates smoothly
between 0 at the centre of the low-energy region and 1 at the edge of the box.
As a consequence, $f_{\textrm{bubble}}(u,x,y)$ interpolates smoothly between
$f_{B}$ and $f_{A}$, as desired. A state generated with this procedure is
shown in Fig. 4. If the subsequent time evolution leads to an expansion of the
bubble, it is convenient to further enlarge the size of the box before
starting the evolution, in order to prevent the bubble from reaching the
boundary of the box before it has reached an asymptotic state. This can be
done simply by “adding” more metastable bath outside the initial box.

Figure 4: Initial energy density profile of a bubble.

Variations of an initial bubble state can be obtained in a simple way. For
example, we can choose different states $B$ for a fixed $A$. As in Bea:2021zsu
, we expect that the subsequent time evolution will quickly select a
dynamically preferred state $C\neq B$ inside the bubble. We could also
multiply the bulk metric functions $B_{1}$ and $B_{2}$ in (7) by some factor,
thus changing the pressure distribution (the anisotropy) along the wall but
not the energy profile. We could further consider initial bubbles whose cross
sections are not perfectly circularly symmetric by starting with an initial
phase-separated state whose low-energy region is comparable to the size of the
box, as in Fig. 3(top).

###  3.3 Critical bubbles

Consider a cylindrical bubble of radius $\rho$ such that the states inside and
outside the bubble correspond to the points marked as $B_{c}$ and $A_{c}$ in
Fig. 2, respectively. The pressure difference between these states generates
an outward-pointing force on the bubble wall. In turn, the surface tension of
the bubble wall results in an inward-pointing force on the wall. A critical
bubble is one for which these two forces exactly balance each other. Since
these bubbles are static, they correspond to equilibrium states. As a
consequence, the temperature must be constant across the entire system and, in
particular, it must be equal to $T_{A_{c}}$. It follows that the state $B_{c}$
is determined by $A_{c}$. If the radius of the bubble is large compared to the
width of the interface between $A_{c}$ and $B_{c}$, then the radius of the
critical bubble takes the form

$\rho_{c}=\frac{\gamma}{{\mathcal{P}}_{B_{c}}-{\mathcal{P}}_{A_{c}}}\,.$ (42)

This follows from approximating the interface by a zero-width surface with
free energy density $\gamma$, assigning a well defined pressure
${\mathcal{P}}_{B_{c}}$, and hence a free energy density
$-{\mathcal{P}}_{B_{c}}$, to the interior of the bubble, and requiring that
the critical bubble locally extremizes the free energy. The fact that this
extremum is a maximum means that the critical bubble is in unstable
equilibrium. This expression for the critical radius is only valid for large
critical bubbles, which are realized when $T_{A_{c}}$ is close to the phase
transition temperature $T_{c}$, namely for $T_{A_{c}}\lesssim T_{c}$. This is
the reason for our choice of the point ${A_{c}}$ in Fig. 2. If the bubble is
not large enough then the phase inside the bubble is not approximately
homogeneous and it cannot be clearly separated from the interface. In this
case one cannot assign a meaningful surface tension to the interface or a well
defined pressure to the interior of the bubble. This situation is realized
when $T_{A_{c}}$ is sufficiently close to the turning point at $T=T_{s}$,
namely when $T_{s}\lesssim T_{A_{c}}$. In this paper we will only discuss
large critical bubbles; small bubbles will be analysed elsewhere.

The fact that critical bubbles are unstable means that supercritical bubbles
expand, whereas undercritical bubbles collapse. Critical bubbles are therefore
the static configurations that separate these two sets of large,
inhomogeneous, cylindrically-symmetric fluctuations of the plasma. This is
precisely the feature that will allow us to identify the critical bubbles with
Jecco.

Figure 5: Time evolution of the wall position for several different initial
bubbles. The critical bubble radius has to be $3.69<\Lambda\rho_{c}<3.75$.

Following the procedure outlined in Sec. 3.2, we generate a family of initial
cylindrical bubbles with different radii and we numerically evolve them with
Jecco. As expected from the discussion above, large bubbles expand and small
bubbles collapse. This is illustrated in Fig. 5, where we plot the radius of
each bubble, defined as the position of the inflection point of the energy
density profile, as a function of time. Each simulation presented in this
section were performed in MareNostrum 4 using 1 node with 48 cores. The
typical runtime was around 250h. We see that bubbles with initial radius
$\Lambda\rho_{c}\geq 3.75$ eventually expand, whereas bubbles with radius
$\Lambda\rho_{c}\leq 3.69$ eventually collapse. This means that the critical
radius must be in between these two values. Substituting into (42) we then
obtain an estimate for the surface tension $\gamma$. Thus,

$3.69<\Lambda\rho_{c}<3.75\,,\qquad 0.116<\frac{\gamma}{\Lambda^{3}}<0.118\,.$

Figure 6: Time evolution of bubbles with different initial radii. The black
curves represent the position of the wall, defined as the inflection point in
the energy density profile. The radius of the critical bubble lies in the
interval $3.69<\Lambda\rho_{c}<3.75$. The bubbles on the left column are
supercritical and they expand. The bubbles on the right column are subcritical
and they collapse. The bubbles in the bottom row are closer to the critical
bubble than those on the top row and hence they evolve more slowly. Videos of
each of the evolutions can be found at
https://youtube.com/playlist?list=PL6eUQq2UUQ4JJTD_pRfJt-ShPGKxGBRW9

Figure 7: Snapshots of the energy density profile of bubbles with different
initial radii. The radius of the critical bubble lies in the interval
$3.69<\Lambda\rho_{c}<3.75$. The bubbles on the left column are supercritical
and they expand. The bubbles on the right column are subcritical and they
collapse. The bubbles in the bottom row are closer to the critical bubble than
those on the top row and hence they evolve more slowly.

As we approach the critical bubble, the dynamics becomes slower and slower.
This feature can be seen in the contour plots of Fig. 6 and in the energy
density snapshots of Fig. 7. In these figures the bubbles in the bottom row
evolve more slowly than those in the top row because their initial radii are
closer to $\rho_{c}$. By fine-tuning the radius of the initial bubble we can
get closer and closer to the critical bubble. Fig. 8 shows that, as we
approach this limit both from above and from below, the bubble profile
converges to a single profile. In this figure we evaluate the profiles at
$\Lambda t=20$ so that the result is not contaminated by the fast-decaying,
transient oscillations present around $\Lambda t=0$ in Fig. 5.

Figure 8: Relative difference between the energy density profiles at
$t\Lambda\sim 20$ of bubbles with different initial radii. We take as a
reference the profile for a bubble with initial radius $\Lambda\rho=3.75$,
which is close to the critical radius. We see that, as this value is
approached both from above and from below, the profiles converge to a single
profile.

The fact that we can approach the critical bubble by fine-tuning a single
parameter is consistent with the fact that the critical bubble should possess
a single unstable mode (see e.g. Laine:2016hma ; Gould:2021ccf ). Indeed, the
latter property means that, in the infinite-dimensional space of
configurations around the critical bubble, the hypersurface of stable
perturbations has codimension one. As we change a single parameter in our
initial data, we trace a curve in the space of configurations that will
generically intersect this hypersurface. If we were to start the time
evolution exactly on this hypersurface, we would remain within it and we would
be attracted to the exact, static critical bubble solution. By tuning the
radius of the bubble in our initial data we come close to this situation and
therefore the dynamics becomes slower and slower.

Since the critical bubble is a static solution, an alternative method to
determine it would be to solve an elliptic problem in two dimensions in AdS,
along the lines of Bea:2020ees .

###  3.4 Expanding bubbles

We now turn to the analysis of expanding bubbles, which play an important role
in the dynamics of first order phase transitions. At sufficiently late times,
the wall of these bubbles is expected to move with a constant velocity, which
results from the balance between the friction that the plasma exerts on the
wall and the pressure difference between the inside and the outside of the
bubble. Moreover, the energy density profile should approach a characteristic
and time-independent shape when plotted as a function of $\rho/t$. In this
section we will use holography to determine both the bubble wall velocity and
the asymptotic profile.

The simulation presented in this section was performed in MareNostrum 4 using
1 node with 48 cores. The typical runtime was around 800h.

####  3.4.1 Wall profile, wall velocity and hydrodynamics

For computational reasons, it is easier to identify the late-time limit for
bubbles that expand at high velocity, since for these configurations the
evolution is faster and we need to run our code for a shorter time to reach
the late-time, asymptotic limit. Based on the mechanical picture we described
above we expect that, as the pressure difference between the inside and the
outside of the bubble grows, the wall velocity will grow too. Therefore, we
will focus on bubbles formed in the large overcooling limit, when the
metastable phase is close to the limit of local stability and the pressure
difference between the inside and outside of the bubble is the largest. For
this reason we will choose the state $A$ outside the bubble as indicated in
Fig. 2, whereas for the state inside we choose the one indicated as $B$.
Following Sec. 3.2, we then construct a bubble that interpolates monotonically
between the states $B$ inside and $A$ outside, as in Fig. 4. This is our
initial state at $t=0$.

Figure 9: (Top) Snapshots of the energy density profile (top) and of the fluid
velocity (bottom) for an expanding bubble. The bubble at $t=0$ interpolates
monotonically between the states $B$ inside and $A$ outside, as in Fig. 4. At
late times the state inside the bubble evolves dynamically to $C$. The states
$A,B$ and $C$ are indicated by black dots in Fig. 2. A full time evolution
video of the energy density can be found at https://youtu.be/wFLp0FSeO8Q

In Fig. 9(top) we show snapshots of the subsequent evolution of the energy
density of the bubble and in https://youtu.be/wFLp0FSeO8Q we show a video of
the full time evolution. As time progresses, the energy density in the
interior of the bubble evolves until it reaches the value corresponding to the
state $C$ in Fig. 2. This means that, as in Bea:2021zsu , this state is
dynamically determined. While the initial configuration at $\Lambda t=0$
interpolates monotonically between the stable and meta-stable branches of the
phase diagram, the expanding bubbles quickly develop a non-monotonic energy
density profile. As illustrated in Fig. 9, the propagation of the bubble leads
to an overheating of the region in front of the bubble that gradually
decreases back to $\mathcal{E}_{A}$ sufficiently far away from the bubble
front. This overheated region possesses non-vanishing energy and momentum
fluxes, which allows us to define a flow velocity via the Landau matching
condition,

$T^{\nu\mu}u_{\mu}=-\mathcal{E}_{\mathrm{loc}}u^{\nu}\,,$ (43)

with $\mathcal{E}_{\mathrm{loc}}$ the energy density of the fluid in the local
rest frame. The flow velocity $v=u^{\rho}/u^{0}$, with $u^{\rho}$ the radial
component of the flow field, for this configuration is shown in Fig.
9(bottom). As we can see in these figures, the region between the bubble wall
and the asymptotic metastable state grows linearly with time as the bubble
expands. As a consequence, we expect that, at late times, the gradients of the
bubble profile decrease and most of the dynamics is captured by hydrodynamics.
We can test this expectation by checking the validity of the hydrodynamic
constitutive relations for the stress tensor in the Landau frame. After
extracting the rest frame energy density and the fluid velocity from the
holographic stress tensor, we can predict the rest of the components of the
stress tensor via the constitutive relations with or without viscous
corrections. The result of this comparison at $\Lambda t=110$ is shown in Fig.
10. We see that hydrodynamics becomes a very good approximation for the
dynamics of the entire system except for the bubble wall, where the failure of
hydrodynamics is expected on general grounds.

Figure 10: Comparison of the holographic stress tensor with the ideal and
viscous hydrodynamic approximations based on the constitutive relations at a
time $\Lambda t=110$ at which the bubble wall is located at
$\Lambda\rho=40.7$. $\mathcal{P}_{z},\mathcal{P}_{\rho},\mathcal{P}_{\varphi}$
and $\mathcal{P}_{\rho\varphi}$ are the stress tensor components in the
$zz-,\rho\rho-,\varphi\varphi-$ and $\rho\varphi-$directions, respectively.
Figure 11: Comparison of the bubble wall profiles at several different times
for the expanding bubble of Fig. 9. To facilitate the comparison, we shift the
position of each curve such that the inflexion point of the different walls at
different times coincide with one another.

Despite its non-hydrodynamic nature, the dynamics of the bubble wall becomes
remarkably simple at sufficiently late times: it moves almost rigidly at
constant velocity. The velocity $v\simeq 0.31$ can be extracted from Fig. 9
via a linear fit to the wall position of the form

$\rho_{\mathrm{wall}}(t)=\rho_{\mathrm{wall},0}+v_{\mathrm{wall}}t\,.$ (44)

To illustrate the rigidity, in Fig. 11 we compare the bubble wall profiles at
several different times. To facilitate the comparison, we shift the position
of each curve such that the inflexion point of the different walls at
different times coincide with one another. We see that the way that the wall
deviates from the inner region $C$ is identical for all sufficiently late
times. In contrast, the maximum value of the energy density at the end of the
wall grows slowly with time. As we will explain in the next section, this
growth indicates that, in the times covered by our simulation, the bubble has
not yet reached the asymptotic late-time form. Despite this, Fig. 11 shows
that the wall has a fixed size set by the microscopic scale of the theory,
$\Lambda$. In particular, the size of the wall does not grow with time, in
contrast with the overheated region in front of the bubble wall.

Figure 12: Comparison of wall profiles for several configurations. “Critical”
refers to the bubble of Fig. 6 with $\Lambda\rho_{0}=3.75$. “Expanding” refers
to the bubble of Fig. 9. “Phase sep.” refers to phase-separated
configurations, be they planar or cylindrical. Each profile has been shifted
and rescaled so that it interpolates between 0 on the left of the wall and 1
on the right. In the case of the expanding bubble, we define $\mathcal{E}_{R}$
as the value of the energy density at the maximum located right in front of
the wall.

In the case of planar bubbles, Ref. Bea:2021zsu  showed that the late-time
wall profile only depends on the asymptotic metastable state $A$. In other
words, the profile is independent of the initial conditions used to generate
the bubble in the first place, as long as they lead to an expanding bubble. We
expect the same conclusion to hold for the cylindrical bubbles considered
here, but it would be interesting to verify it explicitly. Assuming this, it
is interesting to check how the wall profile of an expanding bubble compares
to those of (almost) static walls. For this purpose, in Fig. 12 we compare the
profile of the expanding wall of Fig. 9 with that of the critical bubble of
Sec. 3.3 and with the walls of phase-separated planar and cylindrical
configurations. Following Bea:2021zsu , to facilitate the comparison we shift
and rescale each profile appropriately so that it interpolates between 0 on
the left of the wall and 1 on the right. We achieve this by plotting not just
the energy density $\mathcal{E}(\rho)$ but the combination
$(\mathcal{E}(\rho)-\mathcal{E}_{L})/(\mathcal{E}_{R}-\mathcal{E}_{L})$, with
$\mathcal{E}_{L}$ and $\mathcal{E}_{R}$ the values of the energy density on
the left and on the right of the wall, respectively. In the case of the
expanding bubble, we define $\mathcal{E}_{R}$ as the value of the energy
density at the maximum located right in front of the wall. We see from the
figure that, while all profiles are fairly similar, differences can be seen
with the naked eye. These are more pronounced in the regions where the second
derivative is larger, where they are of the order of 9%.

####  3.4.2 Late-time self-similar solution

As we have seen, for sufficiently late time the bubble wall becomes rigid and
moves at a constant velocity $v_{\mathrm{wall}}$. This implies that the radius
of the region inside the bubble grows linearly with time. Since the energy
density in this region is lower than that in the asymptotic, metastable phase,
this linear growth of the bubble radius must be compensated by a linear growth
in the size of the overheated region in front of the bubble. At very late
times, when all the microscopic scales become irrelevant, this behaviour leads
to a self-similar solution for the bubble that only depends on the ratio
$\rho/t$, as described in e.g. Espinosa:2010hh . In this section we study how
our numerical solutions approach this late-time self-similar solution. For
this purpose, we shift the time and radial coordinates by appropriate amounts
$t_{\rm{shift}}$ and $\rho_{\rm{shift}}$ that we will define below. In other
words, we define

$\xi=\frac{\rho-\rho_{\rm{shift}}}{t-t_{\rm{shift}}}\,.$ (45)

These shifts are motivated by the fact that our initial configuration has a
finite size, and that it takes a certain amount of time for the configuration
to become sufficiently close to the late-time asymptotic solution. While at
asymptotic times these shifts become irrelevant, we find that this procedure
accelerates the convergence to the self-similar regime in our finite-time
simulations.

Figure 13: Expanding bubble profile at late times and inflection points,
marked with vertical grey lines, used to define the size of the overheated
region.

The shifts in question are defined as follows. Consider the overheated region
in front of the bubble wall. This region is connected with the asymptotic
region $A$ by an interface. We begin by locating the inflection point on this
interface, indicated by a vertical line at $\rho=\rho_{\mathrm{interface}}$ in
Fig. 13. We then consider sufficiently late times such that both the wall and
the interface positions move with constant velocity. In this regime
$\rho_{\mathrm{wall}}(t)$ is given by (44) and

$\rho_{\mathrm{interface}}(t)=\rho_{\mathrm{interface},0}+v_{\mathrm{interface}}t\,.$
(46)

We then impose that, as soon as this regime starts, the values of $\xi$ at the
positions of the wall and of the interface immediately agree their late-time
limits. In other words, we adjust the two parameters $t_{\rm{shift}}$ and
$\rho_{\rm{shift}}$ so that the following two conditions are satisfied:

$\xi_{\rm{wall}}=v_{\rm{wall}}\,,\qquad\xi_{\rm{interface}}=v_{\rm{interface}}\,.$
(47)

Figure 14: Energy density (top) and fluid velocity (bottom) profiles for
different simulation times as a function of the scaling variable (45). The
black solid curves correspond to the ideal hydrodynamic prediction.

In Fig. 14 we show the energy density and fluid velocity profiles for
different simulation times as a function of $\xi$. In both plots we see two
regions of fast change that separate three smooth regions. The first region of
fast change occurs around $\xi=v_{\mathrm{wall}}$ and connects the state $C$
in the interior of the bubble, at rest and with a fixed energy density, with
the overheated boosted region in front of the bubble. This abrupt behaviour is
associated to the presence of the bubble wall. Since the size of the wall
remains approximately constant in time, its width in the $\xi$-coordinate
decreases with time. As a consequence, the wall becomes a discontinuity at
asymptotically late times. The shape of the overheated region in front of the
wall is not constant in time. In particular, its slope in the
$\rho$-coordinate decreases with time. However, going to the $\xi$-coordinate
enhances this slope, since at late times $d\mathcal{E}/d\xi\sim
t\,d\mathcal{E}/d\rho$. The curves in Fig. 14 indicate that these two effects
exactly cancel each other at asymptotically late times, resulting in a
constant, non-zero value of the slope in the $\xi$-coordinate in this limit.
The second abrupt region occurs at $\xi\simeq 0.52$ and corresponds to the
interface between the overheated region and the asymptotic metastable region
$A$. In the times covered by our simulations, the width of this interface
grows with time, but this growth is slower than linear. However, it is
possible that, at sufficiently late times, the width of this interface
approaches a constant value. It would be interesting to verify this in the
future through longer simulations. In any case, this interface also approaches
a discontinuity in the $\xi$-coordinate at late times. Despite this, both the
interface and the overheated region are well described by hydrodynamics at
late times, as we saw in Fig. 10.

This discussion suggests that, at asymptotically late times, the bubble
profile should consist of a static inner region $C$ and an outer static region
$A$ connected through discontinuities with an intermediate overheated region
with non-zero fluid velocity. This behaviour agrees with hydrodynamic analysis
of large bubbles, as performed for example in Espinosa:2010hh . At very late
times, when the bubble profile depends only on the scaling variable $\xi$, the
ideal hydrodynamic equations lead to the following equation for the energy
density and the velocity field of a cylindrical bubble

$\displaystyle\gamma^{2}\Big{[}1-\xi\,v(\xi)\Big{]}\Big{[}c_{s}^{2}\mu^{2}-1\Big{]}v^{\prime}(\xi)-\frac{v(\xi)}{\xi}$
$\displaystyle=0\,,$ (48)
$\displaystyle\frac{c_{s}^{2}}{\mathcal{W}}\,\Big{[}1-\xi
v(\xi)\Big{]}\mathcal{E}_{\rm
loc}^{\prime}(\xi)-\gamma^{2}\Big{[}\xi-v(\xi)\Big{]}v^{\prime}(\xi)$
$\displaystyle=0\,,$ (49)

where $\gamma=1/\sqrt{1-v^{2}}$ is the Lorentz factor, $c_{s}$ is the speed of
sound, $\mathcal{E}_{\rm loc}$ is the energy density in the local rest frame
of the fluid,

$\mathcal{W}=\mathcal{E}_{\rm loc}+\mathcal{P}_{\rm eq}(\mathcal{E}_{\rm
loc})$ (50)

is the enthalpy density, and

$\mu=\frac{\xi-v}{1-\xi v}\,.$ (51)

It is well known that the ideal hydrodynamic equation (48) for the fluid
velocity does not posses non-trivial continuous solutions with zero velocity
in the interior and exterior of the bubble. Therefore, in this approximation
the description of an expanding bubble requires the introduction of
discontinuities in the hydrodynamic fields. These discontinuities are
constrained by energy-momentum conservation: although the local energy density
or the fluid velocity may be discontinuous, the energy-momentum flux across
the discontinuity must be continuous. For each value of the wall velocity,
these “junction conditions” at the discontinuities, together with the
hydrodynamic equations elsewhere, determine the entire bubble profile in terms
of the energy density in $A$. This is the reason why a microscopic model is
needed in order to determine the wall velocity. In our case, this model is
provided by holography. Using the holographic prediction for
$v_{\mathrm{wall}}$ as an input, we have solved the hydrodynamic equations
plus the junction conditions and we have determined the profiles represented
by the black solid lines in Fig. 14. The result is consistent with the
holographic profiles at late times in the sense that the holographic curves
approach the black curves more and more as time progresses.

Incidentally, these results allow us to define an analogue of “the state in
front of the bubble wall” for planar bubbles. In the planar case the entire
overheated region in front of the bubble has constant energy density and moves
with constant fluid velocity $v_{D}$ Bea:2021zsu . Using this velocity one can
boost the overheated region to its rest frame and thus define a state in the
phase diagram of Fig. 2. This state was dubbed $D$ in Bea:2021zsu , and the
state in the overheated region was dubbed $D_{\mathrm{boosted}}$. The
difference between $A$ and $D$ gives an intuitive idea of the intensity of the
overheating in front of the wall, since in the absence of it we would have
$A=D$. In the cylindrical case we can obtain a similar idea by defining the
state $D_{\mathrm{boosted}}$ in terms of the maximum values of the black solid
curves in Fig. 14 as we approach the bubble wall discontinuity from the right.
The values we obtain are

${\mathcal{E}}_{D_{\mathrm{boosted}}}=2.26\Lambda^{4}\,,\qquad
v_{D}=0.292\,,\qquad{\mathcal{E}}_{D}=2.06\Lambda^{4}\,.$ (52)

The state $D$ is represented by a black dot in Fig. 2.

##  4 Final remarks

We have presented a new code called Jecco (Julia Einstein Characteristic
Code), which is able to evolve Einstein’s equations coupled to a scalar field
in asymptotically AdS spacetimes using a characteristic formulation. This
implementation generalises the one presented in Attems:2017zam  to 3+1
dimensional settings and further allows, for instance, the usage of other
choices for the scalar potential $V(\phi)$. The code is written in the Julia
programming language Julia-2017  and is freely available at github
https://github.com/mzilhao/Jecco.jl and Zenodo jecco-2022 .

Jecco is written in a modular way, making it an interesting tool to attack
other physical setups. Different problems can be implemented as separate Julia
modules (containing, for example, evolution equations, initial data, and
diagnostic tools) which could be tackled by taking advantage of the general
infrastructure in Jecco (such as finite-difference and pseudo-spectral
derivative operators, filtering tools, and input/output routines).

In the main body of this paper we have presented the formulation, equations of
motion, numerical methods, and the corresponding implementation currently
present in the code. Moreover, in Appendix A we show several tests of this
implementation in various setups, including convergence tests, comparisons
with analytical solutions and an independent numerical implementation,
recovering thermodynamical and quasi-normal mode properties of known
solutions, and checking the constitutive relations of hydrodynamics through
the fluid/gravity prescription. We obtained very good results in all the tests
performed, which reassures us that the implementation is working as intended.

The first new physical application of Jecco was the calculation of the
gravitational wave spectrum produced by a first-order phase transition that
takes place via the instability of the spinodal branch of the phase diagram of
Fig. 2 Bea:2021zol . In this paper we have presented a second application to
the dynamics of bubbles in a strongly-coupled four-dimensional gauge theory.
This extends our previous work on planar bubbles Bea:2021zsu  to cylindrical
bubbles and brings about two new physical aspects. The first one is that the
surface tension now plays a role, and therefore a critical bubble exists in
which the inward-pointing force due to the surface tension exactly balances
the outward-pointing force coming from the pressure difference between the
inside and the outside of the bubble. We have shown that our numerical code
allows us to construct configurations that are arbitrarily close to this
critical bubble. The fact that we can do this with a time evolution code by
fine-tuning a single parameter (which we chose to be the radius of the bubble)
is compatible with the fact that the space of perturbations of a critical
bubble has only one unstable direction. Nevertheless, since the critical
bubble is static, it would be interesting to find it by solving an elliptic 2D
problem in AdS along the lines of Bea:2020ees . This would allow for an
efficient exploration of the bubble properties for the entire range of
temperatures on the metastable branch.

The second new physical aspect brought about by cylindrical bubbles is that
the asymptotic, self-similar profile of an expanding bubble possesses a richer
structure than in the planar case. We have verified this by plotting our
holographic result for the gauge theory stress tensor at late times as a
function of the appropriate scaling variable. We have also compared the
holographic result with the hydrodynamic approximation. As expected, we have
found that hydrodynamics provides a good approximation everywhere except at
the bubble wall.

An immediate extension of this work is to consider multiple expanding bubbles
bubbles . This is an extremely interesting problem because the resulting
bubble collisions will generate gravitational waves. As in previous
applications of holography to the quark-gluon plasma Casalderrey-
Solana:2011dxg ; Busza:2018rrf  or to condensed matter systems Zaanen:2015oix
; Hartnoll:2016apf ; Nastase:2017cxp , we expect that the first-principle
nature of the holographic approach will shed new light on this problem too.

###### Acknowledgements.

It is a pleasure to thank Bartomeu Fiol, Oscar Henriksson, David Hilditch,
Mark Hindmarsh, Carlos Hoyos, Christiana Pantelidou and Oriol Pujolàs for
discussions. YB acknowledges support from the European Research Council Grant
No. ERC-2014-StG 639022-NewNGR and the Academy of Finland grant no. 333609. TG
acknowledges financial support from FCT/Portugal Grant No. PD/BD/135425/2017
in the framework of the Doctoral Programme IDPASC-Portugal. AJ acknowledges
support from the European Research Council Grant No. ERC-2016-AvG
692951-GravBHs. MSG acknowledges financial support from the APIF program,
fellowship APIF_18_19/226. JCS, DM and MSG are also supported by grants
SGR-2017-754, PID2019-105614GB-C21, PID2019-105614GB-C22 and the “Unit of
Excellence MdM 2020-2023” award to the Institute of Cosmos Sciences
(CEX2019-000918-M). MZ acknowledges financial support provided by FCT/Portugal
through the IF programme grant IF/00729/2015 and CERN project CERN/FIS-
PAR/0023/2019, as well as the support by the Center for Research and
Development in Mathematics and Applications (CIDMA) through FCT/Portugal,
references UIDB/04106/2020, UIDP/04106/2020 and the projects PTDC/FIS-
AST/3041/2020 and CERN/FIS-PAR/0024/2021. We further acknowledge support from
the European Union’s Horizon 2020 research and innovation (RISE) program
H2020-MSCA-RISE-2017 Grant No. FunFiCO-777740. The authors thankfully
acknowledge the computer resources, technical expertise and assistance
provided by CENTRA/IST. Computations were performed in part at the cluster
“Baltasar-Sete-Sóis” and supported by the H2020 ERC Consolidator Grant “Matter
and strong field gravity: New frontiers in Einstein’s theory” grant agreement
No. MaGRaTh-646597. We also thank the MareNostrum supercomputer at the BSC
(activity Id FI-2020-3-0010, FI-2021-1-0008 and FI-2021-3-0010) for
significant computational resources.

##  Appendix A Tests of Jecco

To gauge the performance, accuracy and reliability of Jecco we conduct a
number of tests. These tests include comparing the data from numerical
simulations against known analytical results, as well those from the 2+1 SWEC
code introduced in Attems:2017zam . We also perform convergence tests and
contrast obtained results against expected physical quantities and properties
of our model systems, such as the black brane entropy density and the
frequencies of its quasi-normal modes. Unless specifically mentioned, results
will be presented in “code units”, where $G=c=\hbar=L=1$.

We note that we solve the equations of motion of our Einstein-scalar model (1)
using the ingoing Eddington-Finkelstein gauge (equation (7)), which is a
Bondi-like gauge, and the resulting PDE system is expected to be only weakly
hyperbolic Giannakopoulos:2020dih . We thus restrict our tests to smooth data,
where the effect of weak hyperbolicity is not expected to be manifested
Giannakopoulos:2020dih .

As mentioned in the main text, for the moment we have only implemented shared-
memory parallelism using Julia’s Threads.@threads macro. We have performed
some simple scaling tests with an AMD Ryzen 9 5950X 16-Core Processor and we
see a speedup factor of 2.7 when running with 4 threads, 3.5 with 8 threads,
and 4.5 with 16 threads. The bottleneck comes from an operation within the
DifferentialEquations.jl package which does not seem to be parallelized. We
plan to investigate this further in the near future.

###  A.1 Analytical black brane

In these tests the code is initiated in a homogeneous black brane
configuration, which is a static exact solution of the equations of motion
with $\phi_{0}=0$ (conformal case). The functions specified in the initial
data vanish and the only non-vanishing boundary data are $a_{4}=-4/3$. For
most of these tests, we do not perform a time evolution but instead we just
solve the whole nested system at $t=0$ and compare the last bulk function to
be computed, that is $A$, against its analytic form:

$\displaystyle
A=\frac{1}{u^{2}}+\frac{2\xi}{u}+\xi^{2}+\frac{a_{4}\,u^{2}}{1+2\,\xi\,u+\xi^{2}\,u^{2}}\,,$
(53)

using the field redefinitions of Sec. 2.1.3 appropriately. From (53) we see
that the gauge fixing can be performed via

$\displaystyle\xi=\left(-a_{4}\right)^{-1/4}-1/u_{H}\,,$ (54)

with $u_{H}=1$ the gauge fixed position of the apparent horizon for the tested
configuration. Since Jecco provides us with the possibility of multiple outer
spectral domains, we wish to understand to what extent faster configurations
compromise the accuracy of the numerical solution. We vary the number of nodes
in the $u$-domains, as well as the number of outer $u$-domains, to examine the
accuracy of the code for different configurations of the spectral grid. The
inner $u$-domain discretizes the region $[0,0.1]$ and the outer one the region
$[0.1,1.0]$. The domain of both the transverse directions $x$ and $y$ is
$[-5,5)$ and is discretized uniformly with 128 nodes in each case.

Figure 15: The maximum relative errors for the bulk function $A$, in the outer
radial domains, for different configurations of the test against the
analytical homogeneous black brane static solution. The same accuracy for this
test is achieved e.g. by three outer radial domains with 32 nodes per domain,
and a single domain with 56 nodes. The former configuration is faster.

The maximum relative error of $A$ for the inner spectral domain remains below
$O(10^{-10})$ for a range of nodes between 12 and 36\. The respective error
for different configurations of outer spectral domains is shown in Fig. 15. A
maximum relative error below $O(10^{-5})$ in the outer region can be achieved
with one or multiple domains, where the latter typically provides faster
configurations. The difference in orders of magnitude between the maximum
relative error of the inner and outer domains is due to the near boundary
field redefinition. This redefinition factors out the near boundary radial
dependence of the field and allows for a more accurate numerical solution. For
completeness, we perform a time evolution for one of the aforementioned
configurations, even if the evolution is expected to be trivial since we are
investigating a static setup. For a configuration with 12 nodes in the inner
domain and 28 nodes on each of the three outer domains we have verified that
the maximum error maintains its expected value even after 550 timesteps, which
corresponds to $t_{f}=2$ in code units. For the time integration the third
order Adams-Moulton method with adaptive step is used.

For a generic physical setup we find that some experimentation may be required
to find the optimal numerical parameters, like the number of outer domains and
nodes per domain, the choice of time integrator, etc. For instance, if
accuracy of temporal derivatives of the solution is important one might
consider chosing a fixed timestep integrator with a small timestep instead of
an adaptive one. If the main focus is the late-time behaviour of the solution,
perhaps an adaptive step integrator is preferable.

###  A.2 Comparison with SWEC

For this test the code is initialized with an $x$-dependent perturbation on
top of a homogeneous black brane configuration. The initial data are

$\displaystyle B_{1}(0,u,x,y)$ $\displaystyle=0.01u^{4}\,,$ (55)
$\displaystyle a_{4}(0,x,y)$ $\displaystyle=-\frac{3}{4}\left[1+\delta
a_{4}\cos\left(2\pi
k_{x}\frac{x-x_{\textrm{mid}}}{x_{\textrm{max}}-x_{\textrm{min}}}\right)\right]\,,$
$\displaystyle\xi(0,x,y)$ $\displaystyle=\left(\frac{4}{3}\right)^{1/4}-1\,,$

where $\delta a_{4}=5\cdot 10^{-4}$, and the remaining free data functions
($B_{2}$, $G$, $\phi$, $f_{x2}$, $f_{y2}$) are set to zero. We compare the
error of the numerical solution provided by Jecco against that of the SWEC
code used in Attems:2017zam , for the same setup.

We use one inner radial domain spanning the region $u\in[0,0.1]$ discretized
with 12 grid points, and another (outer) domain spanning the region
$u\in[0.1,1.01]$ with 48 grid points. The transverse direction $x$ spans
$x\in[-10,10)$, which is discretized with 128 grid points, while the $y$ has
trivial dynamics for this setup (and 6 grid points are used so that the finite
difference operator fits in the domain). The time evolution is performed using
the fourth-order accurate Adams-Bashforth method. The evolution is performed
for a total of 2000 time steps. The choice of a single outer radial domain in
Jecco is made for a more explicit comparison against SWEC, since the latter
does not offer the possibility of multiple outer radial domains. It is worth
noticing, however, that there are still differences between the setups in the
two codes. For instance, the inner and outer domains of Jecco share only one
common radial point, whereas in SWEC there is an overlapping $u$-region
between them.

We show relative differences between the $a_{4}$ and $\xi$ functions obtained
in the two codes in Fig. 16. The pattern observed was similar for the metric
function $B_{1}$. To compare the output of the two codes exactly on the same
grid points we perform cubic spline interpolation on the data and use the
values of the interpolated functions for the comparison. It is reassuring that
the results from the two codes agree so well.

Figure 16: Relative errors for the $a_{4}$ and $\xi$ functions at the end of
the evolution. Results obtained with the SWEC code are used as benchmark.

###  A.3 Convergence tests

We now show convergence tests using numerical solutions obtained only from
Jecco. For this, we solve the same physical setup with increasing resolution
and inspect the rate at which the numerical solution tends to the exact one.
The rate at which numerical error tends to zero with increasing resolution is
determined by the approximation accuracy. The latter is the degree to which a
discretized version of a PDE system approximates the correct continuum PDE
system, and such a discretized version is called consistent. If its numerical
solution is bounded at some arbitrary finite time by the given data of the
problem in a discretized version of a suitable norm, it is furthermore called
stable. The Lax equivalence theorem states that consistency of the finite
difference scheme and stability with respect to a specific norm guarantee
convergence for linear problems (and the converse) LaxRic56 .

For our present case, since the spatial discretization is performed with a
mixture of finite-difference and pseudo-spectral techniques, we fix the number
of grid points along the spectral direction and vary only the number of grid
points in the uniform grid along the transverse directions $x,y$. The finite-
difference operators dominate the numerical error, so the expected convergence
rate is controlled by the rate at which we increase the resolution in the
uniform grid, as well as the approximation order of the operators.

Let us denote by $f$ the solution to the continuum PDE problem and by $f_{h}$
its numerical approximation. We have

$\displaystyle f=f_{h}+O(h^{n})\,,$ (56)

where $h$ is the grid spacing and $n$ the accuracy of the finite-difference
operators.

Consider performing numerical evolutions with coarse, medium and fine
resolutions $h_{c}$, $h_{m}$ and $h_{f}$ respectively. Then one can construct
the quantity

$\displaystyle
Q\equiv\frac{h_{c}^{n}-h_{m}^{n}}{h_{m}^{n}-h_{f}^{n}}=\frac{f_{h_{c}}-f_{h_{m}}}{f_{h_{m}}-f_{h_{f}}}\,,$
(57)

often called the convergence factor, which informs us about the rate at which
the numerical error induced by the finite-difference scheme converges to zero.
Comparison of grid functions corresponding to different resolutions is to be
understood by the use of the common grid points among the different
resolutions.

Using a physical setup with known exact solution provides a clear benchmark to
compare with, and we can prepare such a setup by evolving a homogeneous black
brane with only gauge dynamics. This can be achieved by using a different
choice for the evolution of the gauge function $\xi$ than the one specified in
Sec. 2.1.4. In particular, we impose the advection equation

$\displaystyle\partial_{t}\xi(t,x,y)=-v_{x}\,\partial_{x}\xi(t,x,y)\,,$ (58)

which introduces non-trivial dynamics to the numerical evolution.

The only non-vanishing initial data for this setup is the boundary function
$a_{4}$, which we set to $a_{4}(t,x,y)=-1$, and the gauge function $\xi$,
which we initialize to

$\displaystyle\xi(0,x,y)=\xi_{0}+A_{x}\sin\left(\frac{2\pi\,n_{x}}{L_{x}}\left(x_{\textrm{max}}-x\right)\right)\,,$
(59)

where $L_{x}\equiv x_{\textrm{max}}-x_{\textrm{min}}$. For such a
configuration, the solution to equation (58) is

$\displaystyle\xi(t,x,y)=\xi_{0}+A_{x}\sin\left(\frac{2\pi\,n_{x}}{L_{x}}\left(x_{\textrm{max}}-x+v_{x}t\right)\right)\,,$
(60)

and the exact solution of the metric function $A$ is given by (53), where
$\xi$ is now provided by (60).

For the tests presented herein we have fixed

$\displaystyle\xi_{0}=0\,,\quad A_{x}=0.1\,,\quad n_{x}=1\,,\quad
x_{\textrm{max}}=5\,\quad x_{\textrm{min}}=-5\,.$

For the numerical discretization we have employed one inner radial domain with
12 grid points (spanning the region $u\in[0,0.1]$) and three equal-sized outer
domains for the region $u\in[0.1,1.2]$ with 28 grid points each. For the
transverse directions we use 16, 32, and 64 grid points for coarse, medium and
fine resolution respectively. The time integration is done with the third-
order accurate Adams-Moulton method, with adaptive timestep. We have performed
these tests with both second- and fourth-order accurate (periodic) finite
difference operators, where Kreiss-Oliger dissipation is used with the
prescription of equation (32) with $\sigma=0.01$. The tests were run on a
laptop with an Intel Core i7-10510U at 1.80GHz CPU. For the fourth-order
accurate finite difference case, the coarse resolution ran with a single
thread and was completed within 36 minutes. The corresponding medium and high
resolution cases were performed with two threads and were completed within 66
and 271 minutes, respectively.

Figure 17: (Left) Pointwise convergence of the metric function $A$ along the
$x$ direction, at $t=9.98$ (code units), $u=0.83$ and $y=0.625$, for the
medium and fine resolutions. (Right) Convergence rate for the metric function
$A$ in terms of rescaled norms. Perfect overlap of curves should be understood
as perfect convergence. Top rows show the second-order finite difference
approximation case and bottom rows show the fourth-order one. The expected
convergence factor for the former is $Q=4$ and the latter $Q=16$ for our
specific setup.

Convergence tests for the $A$ metric function can be seen in Fig. 17. As
mentioned above, the comparison of the grid functions against the exact
solution is performed only on grid points that are common to all three
resolutions. The expected convergence factor for this setup is $Q=4$ for
second-order finite difference operators and $Q=16$ for fourth-order ones,
which is indeed what we observe in the left column. The same convergence rate
is expected when we perform a norm comparison. The discretized version of the
$L^{2}$-norm that we employ here is simply the square root of the sum of the
squared grid function under consideration (over all domains). In the right
column of the figure we again see very good agreement for the norm convergence
rate.

We also perform convergence tests for the setup that results in the top phase-
separated configuration of Fig. 3. In this case, the initial data comprises of
the sinusoidal perturbation (39) with $\bar{\alpha}_{4}=1$ and
$\delta\alpha=10^{-3}$, as well as $\phi_{0}=\Lambda=1$, $\phi_{2}=0.3$, and
$\bar{\mathcal{E}}=1$. The size of the box is $L_{x}\Lambda=L_{y}\Lambda=10$.
The discretization of the transverse and holographic domains, as well as the
time integrator are the same as for the previous convergence test, with the
only difference that the outer holographic domain here resolves the region
$u\in[0.1,1.05]$. We use second order finite difference operators and set
$\sigma=10^{-5}$. Since we do not have an exact solution, we perform self
convergence tests using only numerical results. The comparisons are performed
again using the common points of the coarse grid.

In Fig. 18 we present pointwise and norm convergence tests for the boundary
energy density $\bar{\mathcal{E}}$ of the above configuration. Notice that the
runs performed for these tests reach $t\Lambda=21.69$, whereas the top phase-
separated profile of Fig. 3 corresponds to $t\Lambda=300$ of the setup. Since
we are using a low value for the dissipation parameter ($\sigma=10^{-5}$), it
is not possible to perform such long runs. The reason for this choice is that
high values of $\sigma$ seem to non-trivially affect the convergence
properties of these configurations. However, we have checked that when
performing the same runs with $\sigma=0.2$, which is sufficient for long runs,
the difference when comparing to the setups with $\sigma=10^{-5}$ drops fast
with increasing resolution, as illustrated in Fig. 19.

Figure 18:  (Left) Pointwise convergence of the boundary energy density
$\bar{\mathcal{E}}$ along the $x$ direction for fixed $y=-1.25$, at $\Lambda
t=21.69$ (code units). (Right) Convergence rate for the boundary energy
density in terms of rescaled norms, until $\Lambda t=21.69$. Perfect overlap
of curves should be understood as perfect convergence. We see good second
order convergence throughout the evolution.  Figure 19:  The difference in
$\bar{\mathcal{E}}$ for runs with low and high Kreiss-Oliger dissipation, with
$\sigma=10^{-5}$ and $\sigma=0.2$, respectively. The norm of this difference
is illustrated until $\Lambda t=21.69$. We observe that this difference
decreases with increasing resolution.

###  A.4 Thermodynamics tests

Let us now explore how well the code can recover known properties of non-
conformal homogeneous black branes. For concreteness we will focus on cases
with $\lambda_{4}=-0.25$ and $\lambda_{6}=0.1$ for the model given by equation
(4).

We initialize the code to some homogeneous (along $x$ and $y$) and isotropic
state, setting $B_{1}=B_{2}=G=0$ and, as we are not interested in non-zero
momenta, $f_{x2}=f_{y2}=0$. We set $a_{4}$ and the (initial) gauge parameter
$\xi$ as follows

$a_{4}=-\frac{4}{3}\left(\mathcal{E}+\phi_{0}\phi_{2}+\left(\lambda_{4}-\frac{7}{36}\right)\phi_{0}^{4}\right),\qquad\xi=\left(-a_{4}\right)^{-1/4}-\frac{1}{u_{H}},$
(61)

where $\mathcal{E}$ is the energy density of the black brane, $u_{H}=1.0$ is
the location at which the apparent horizon will be placed, and we choose
$\phi_{2}=0.29819$ and $\phi_{0}=1.0(=\Lambda)$.333The value of $\phi_{2}$
should not be too far away from the equilibrium value of the non-conformal
black brane with energy density $\mathcal{E}$. Otherwise, the initial choice
for the gauge function $\xi$ may not ensure that the apparent horizon lies
inside the numerical domain. Motivated by its near boundary behaviour, we
initialize the scalar field to

$\phi(u)=\phi_{0}u-\xi_{0}\phi_{0}u^{2}+u^{3}\left(\xi^{2}\phi_{0}+\phi_{2}\right).$
(62)

We then let the code evolve. Since this scalar profile is not an equilibrium
configuration, the system will relax in a few time units to the non-conformal
uniform black brane with the given energy density $\mathcal{E}$.

Figure 20: (Up-left) Pressure, (Up-right) temperature, (Bottom-left) entropy
density and (Bottom-right) the scalar VEV, computed solving the static
Einstein’s equations (blue line) compared against the results obtained from
Jecco after an evolution until $t\Lambda=60$ (red dots).

We performed a total of 16 runs with energies evenly distributed in the
interval $\mathcal{E}/\Lambda^{4}\in[0.4,2.0]$ and compared the results with
those obtained from directly integrating the static solution of Einstein’s
equations for the same physical configuration. Each run was performed using a
single core in a 16GB memory machine, the runtime being a few minutes. The
chosen range of energies is of most relevance because it completely contains
the first order thermal phase transition exhibited for this value of
$\lambda_{4}$ and $\lambda_{6}$ (see Fig. 20). For higher and lower energies
the theory tends to the conformal case, which was explored previously. Fig. 20
shows that the results obtained by both methods lie on top of each other. The
pressures along the 3 boundary directions are equal to each other (there are
no anisotropies), and the behaviour of the energy density as a function of the
temperature (up-right panel) shows the typical behaviour of a theory with a
first order phase transition, see Bea:2018whf ; Bea:2020ees . Notice that the
off-diagonal pressure, $\mathcal{P}_{xy}$, and the energy fluxes,
$\mathcal{J}_{x}$ and $\mathcal{J}_{y}$, are not shown as they are vanishing
for these solutions.

Figure 21: (Left) Absolute error of the vanishing quantities and (Right)
relative error of the non vanishing ones. We took the value obtained by
solving the static Einstein’s equations as the reference for the relative
error.

In Fig. 21 we show the differences between the quantities obtained with these
two methods; we plot absolute differences for those quantities that vanish and
relative differences for the non-trivial ones. As can be seen, the off-
diagonal pressure and the energy fluxes have vanishingly small values; for the
non-trivial quantities, the pressure presents the largest relative error,
which is smaller than $0.03\%$. We thus see that Jecco is returning the
expected properties of these solutions with very good accuracy.

###  A.5 Quasi-normal mode tests

We now show results for a time-dependent test, where we recover expected
quasi-normal mode frequencies. This test is a replica of the one performed for
SWEC in Attems:2017zam , and was performed using a single core in a 16GB
memory machine running for a few minutes.

We fix $\phi_{0}=1$, $\lambda_{4}=0.0025$ and $\lambda_{6}=0$. We set as
initial conditions $G=0$ together with $f_{x2}=f_{y2}=0$. $a_{4}$ is obtained
from equation (61) with an average energy density $\mathcal{E}=0.379686$ and
$\phi_{2}=0.0868357$, which corresponds to the equilibrium value of the non-
conformal uniform black brane with that same energy density. For this test
$\xi=0$ (initially) is good enough a choice. The vanishing wave number ($k=0$)
perturbation is inserted by activating the anisotropy together with a non
equilibrium scalar field profile:

$\displaystyle\phi(u)$ $\displaystyle=\phi_{0}u+\phi_{2}u^{3}\,,$ (63)
$\displaystyle B_{1}(u)$ $\displaystyle=0.15u^{8}\,,$ $\displaystyle B_{2}(u)$
$\displaystyle=0.05u^{8}\,.$

Notice that, since the perturbation is independent of both $x$ and $y$,
$B_{1}$ and $B_{2}$ will behave identically up to an overall constant factor.

The system will relax to the equilibrium state through damped oscillations
whose parameters were extracted in Attems:2016ugt  for different values of
$\lambda_{4}$. In particular, the boundary variables $\phi_{2}$, $b_{14}$ and
$b_{24}$ will evolve in time according to

$\displaystyle f(t)$ $\displaystyle=f_{\rm eq}+f^{(1)}(t)+f^{(2)}(t)+\cdots,$
(64) $\displaystyle f^{(n)}(t)$
$\displaystyle=A_{n}e^{-\omega_{i}^{(n)}t}\cos(\omega_{r}^{(n)}t+\delta_{n}),$
$\displaystyle\omega^{(n)}$
$\displaystyle=\omega_{r}^{(n)}+i\omega_{i}^{(n)},$

with the mode 1, $f^{(1)}(t)$, being the longest lived (smallest
$\omega_{i}$).

Figure 22:  (Left column) Time evolution of the boundary variables. Red dotted
lines correspond to the longest lived mode, while the black dotted ones are
the full combination given by equation (64). (Right column) Log plots of
$f-f_{\rm eq}$ for the boundary variables. Red dotted lines correspond to the
exponential decay that best fits the data at late times.

We can obtain the different parameters from the data for $f-f_{\rm eq}$, whose
log-plots are shown in Fig. 22. For late times, the longest-lived mode
dominates and the data clearly behaves as a damped oscillation. We use this
fact to fit the $f^{(1)}(t)$ to the data. Once we have the parameters for mode
1, the shorter lived mode can be obtained by fitting $f(t)$ at early times,
where its presence is still important. As a consequence we get an improvement
in the description, specially at early times, as can be seen in the left
column of the figure.

We find the following values for the frequencies

$\displaystyle\phi_{2}:$ $\displaystyle\omega^{(1)}$
$\displaystyle=2.313004+i1.264244,\qquad$ $\displaystyle\omega^{(2)}$
$\displaystyle=4.091354+i2.944895,$ (65) $\displaystyle b_{14}:$
$\displaystyle\omega^{(1)}$ $\displaystyle=3.039641+i2.120511,\qquad$
$\displaystyle\omega^{(2)}$ $\displaystyle=4.953161+i3.737359,$ $\displaystyle
b_{24}:$ $\displaystyle\omega^{(1)}$ $\displaystyle=3.039649+i2.120516,\qquad$
$\displaystyle\omega^{(2)}$ $\displaystyle=4.952340+i3.741593,$

The frequencies can also be obtained directly from the equilibrium solution by
solving linear perturbation equations. Taking the final equilibrium solution
of this evolution and computing the linearized fluctuations around this
background using the Mathematica package QNMspectral Jansen:2017oag , we
obtain

$\displaystyle\phi_{2}:$ $\displaystyle\omega^{(1)}$
$\displaystyle=2.313080+i1.264337,\qquad$ $\displaystyle\omega^{(2)}$
$\displaystyle=4.108219+i2.931352,$ (66) $\displaystyle b_{14}:$
$\displaystyle\omega^{(1)}$ $\displaystyle=3.039399+i2.120359,\qquad$
$\displaystyle\omega^{(2)}$ $\displaystyle=4.934072+i3.739264,$

and for $b_{24}$ identical results to those for $b_{14}$ as they are both
tensor fluctuations.

The agreement among the values obtained by both methods is excellent for the
lowest-frequency modes, the easiest to extract, with a relative error under
$0.01\%$ in all cases. For the shorter-lived modes, as expected, the relative
error is higher, but the agreement is still very good, always below $0.5\%$.

###  A.6 Fluid/gravity tests

The fluid/gravity duality establishes a precise map between the equations of
relativistic hydrodynamics in $d$ dimensions and the Einstein equations with
negative cosmological constant in $d+1$ dimensions in a specific regime (see
Hubeny:2011hd ; Rangamani:2009xk  for comprehensive reviews). This is a map
between non-linear equations, and solutions on one side map to solutions on
the other side. Even if originally derived from holography, fluid/gravity is
an independent statement and constitutes a duality between two classical
theories. This represents a complementary test to the one of Sec. A.5.

We will now use this fluid/gravity mapping to test the code. The idea of this
test is to consider a microscopic holographic evolution, which by construction
is in the regime of hydrodynamics, and then compare this microscopic evolution
against the constitutive relations of hydrodynamics at every spacetime point.

The constitutive relations of hydrodynamics (see e.g. Romatschke:2009kr )
truncated at first order in the hydrodynamic gradient expansion take the form

$T^{\rm hydro}_{\mu\nu}=T^{\rm ideal}_{\mu\nu}+T^{\rm
1st}_{\mu\nu}+\cdots~{},$ (67)

with

$\displaystyle T^{\rm ideal}_{\mu\nu}$ $\displaystyle=\mathcal{E}_{\rm
loc}u_{\mu}u_{\nu}+P(\mathcal{E}_{\rm loc})\Delta_{\mu\nu}~{},$ (68)
$\displaystyle T^{\rm 1st}_{\mu\nu}$ $\displaystyle=-\eta(\mathcal{E}_{\rm
loc})\sigma_{\mu\nu}-\zeta(\mathcal{E}_{\rm loc})(\nabla\cdot
u)\Delta_{\mu\nu}~{},$ (69)

where $\mathcal{E}_{\rm loc}$ is the energy density in the local rest frame of
the fluid, $u_{\mu}$ is the local fluid velocity,
$\Delta_{\mu\nu}:=\bar{\eta}_{\mu\nu}+u_{\mu}u_{\nu}$ is the projector,
$P(\mathcal{E}_{\rm loc})$ is the equation of state, and $\eta$ and $\zeta$
are the shear and bulk viscosities, respectively. Moreover, we define
$\sigma_{\mu\nu}:=2\nabla_{\langle\mu}u_{\nu\rangle}$, where
$\langle\cdot\rangle$ indicates symmetrization, tracelessness and
orthogonality to the velocity. The equation of state and the viscosities are
determined by the specific microscopic theory under consideration, and in our
case we obtain them by constructing the set of homogeneous black branes and by
using Kubo formulas (see Bea:2020ees ; Attems:2016ugt ).

We consider as initial state a homogeneous black brane solution with small
sinusoidal perturbations along $x$ and $y$. For this test we choose
$\lambda_{4}=-0.25$, $\lambda_{6}=0.1$, $x,y\in[-100,100)$ with
$N_{x}=N_{y}=100$ grid points in the $x$ and $y$ directions. The average
energy density is fixed to $\bar{\mathcal{E}}=10\Lambda^{4}$ and the $a_{4}$
function is chosen to be

$a_{4}(0,x,y)=\bar{a}_{4}\left[1+0.001\left(\cos\left(\frac{\pi
x}{100}\right)+\cos\left(\frac{\pi y}{100}\right)\right)\right],$ (70)

where the average value, $\bar{a}_{4}$, is determined by equation (61) (left)
with our chosen value of $\bar{\mathcal{E}}$. The initial gauge parameter
$\xi$ is obtained from equation (61) (right) replacing $a_{4}$ by
$\bar{a}_{4}$. The scalar bulk profile is again chosen to be given by equation
(62), with $\phi_{0}=1.0(=\Lambda)$ and $\phi_{2}=0.29819$. If the momentum of
the perturbation $k$ is small compared to the temperature of the black brane
$T$, the system will be within the regime of hydrodynamics. For this
simulation the ratio is $k/T\simeq 0.051$.

We evolve this initial configuration with Jecco and compare the obtained
boundary stress tensor as a function of time with the constitutive relations
of hydrodynamics. As we will see below, we find very good agreement for all
components of the stress tensor.

Figure 23: Real-time evolution from Jecco and comparison with the constitutive
relations of hydrodynamics. (Top-left) Initial energy density configuration: a
homogeneous state plus small sinusoidal perturbations along $x$ and $y$. (Top-
right) $T_{xy}$ component of the stress tensor as a function of time and $y$
(extracted at $x=50/\Lambda$) obtained from the Jecco. (Bottom-left)
Difference between the $T_{xy}$ obtained from the numerical evolution and
$T^{\rm Id}_{xy}$ given by the constitutive relations of ideal hydrodynamics.
This difference is very small compared to $T_{xy}$, indicating that ideal
hydrodynamics provides a good description. (Bottom-right) We further include
the first order terms of hydrodynamics in the previous subtraction, obtaining
an even better description.

See Fig. 23 (top-left) for the initial energy density configuration. The
system has vanishing initial velocity. The fact that we are not initializing
Jecco with the equilibrium $\phi(u)$ does not affect the results, as the time
that this takes to decay (through quasi-normal modes) to the equilibrium
profile for the specified average energy density is much shorter than the time
scale of the dynamics triggered by the sinusoidal perturbation.

In the following, we present the results for the $T_{xy}$ component of the
stress tensor. For all other components the results are similar. The component
$T_{xy}$ is particularly interesting because it allows to test proper dynamics
in 2+1 dimensions on the boundary. Moreover, the constitutive relations of
hydrodynamics for this component are purely non-linear (the linearized
expression vanishes), so this provides also a truly non-linear test.

Figure 23 (top-right) shows the values of $T_{xy}$ obtained from Jecco, at
constant $x=50$ as a function of time, and the bottom-left panel shows the
difference between $T_{xy}$ and $T^{\rm Ideal}_{xy}$ – given by the
constitutive relations of ideal hydrodynamics. This difference is very small
compared to $T_{xy}$, indicating that ideal hydrodynamics provides a very good
description (within $0.01\%$). In Fig. 23 (bottom-right) we further include
the first order terms of hydrodynamics in the previous subtraction, obtaining
an even better description (within $0.0001\%$). Presumably, this difference
would be well described by second order hydrodynamics, but we lack the
corresponding coefficients to do this check.

We conclude that hydrodynamics provides a very good description of the system,
in consonance with the fluid/gravity mapping. In particular, we observe that
first order hydrodynamics further improves the ideal description, as expected
from the hydrodynamic gradient expansion. We emphasise that this test
constitutes a truly non-linear precision test of both the code and the
fluid/gravity correspondence in a real-time dynamical configuration.

##  Appendix B Radial equations

For completeness, here we list the radial equations obtained from ansatz (7).
It is useful to introduce the following operators to make all expressions more
compact

$\displaystyle f^{\prime}$ $\displaystyle=\partial_{r}f\,,$ (71)
$\displaystyle\dot{f}$
$\displaystyle=\left(\partial_{t}+\frac{A}{2}\partial_{r}\right)f\,,$
$\displaystyle\tilde{f}$
$\displaystyle=\left(\partial_{x}-F_{x}\partial_{r}\right)f\,,$
$\displaystyle\hat{f}$
$\displaystyle=\left(\partial_{y}-F_{y}\partial_{r}\right)f\,,$
$\displaystyle\bar{f}$
$\displaystyle=\left(\partial^{2}_{x}-2F_{x}\partial_{r}\partial_{x}+F_{x}^{2}\partial^{2}_{r}\right)f\,,$
$\displaystyle f^{\star}$
$\displaystyle=\left(\partial^{2}_{y}-2F_{y}\partial_{r}\partial_{y}+F_{y}^{2}\partial^{2}_{r}\right)f\,,$
$\displaystyle f^{\times}$
$\displaystyle=\left(\partial_{x}\partial_{y}-F_{x}\partial_{r}\partial_{y}-F_{y}\partial_{r}\partial_{x}+F_{x}F_{y}\partial^{2}_{r}\right)f\,.$

As shown in Table 1, by combining Einstein’s equations (2) in a particular way
we obtain a nested system of radial ODEs where one can sequentially solve for
the different variables. This is common to these characteristic approaches,
see e.g. Winicour:2012znc ; Chesler:2013lia . In this case, some of these
equations will be coupled. The full set of equations, in the order to be
solved, is:

$6S^{\prime\prime}+S\left(\cosh^{2}(G)\left(B_{1}^{\prime}\right)^{2}+3\left(B_{2}^{\prime}\right)^{2}+\left(G^{\prime}\right)^{2}+\left(\phi^{\prime}\right)^{2}\right)=0$
(72)

$2e^{B_{1}}S^{2}F_{x}^{\prime\prime}+e^{B_{1}}\left(S^{2}\left(-2\left(\cosh^{2}(G)\left(\tilde{B_{1}}^{\prime}-B_{1}^{\prime}F_{x}^{\prime}\right)+B_{2}^{\prime}\left(3\tilde{B_{2}}-F_{x}^{\prime}\right)+\tilde{G}\left(B_{1}^{\prime}\sinh(2G)+G^{\prime}\right)+\tilde{B_{2}}^{\prime}+4\tilde{\phi}\phi^{\prime}\right)-2\tilde{B_{1}}B_{1}^{\prime}\cosh^{2}(G)\right)+S\left(-6\tilde{S}\left(B_{1}^{\prime}\cosh^{2}(G)+B_{2}^{\prime}\right)-8\tilde{S}^{\prime}+2S^{\prime}F_{x}^{\prime}\right)+8\tilde{S}S^{\prime}\right)+S^{2}\left(-2G^{\prime}\left(\hat{B_{1}}+F_{y}^{\prime}\right)+\sinh(2G)\left(\hat{B_{1}}^{\prime}-B_{1}^{\prime}\left(\hat{B_{1}}+F_{y}^{\prime}\right)\right)+2\hat{G}B_{1}^{\prime}\cosh(2G)+2\hat{G}^{\prime}\right)+3\hat{S}S\left(B_{1}^{\prime}\sinh(2G)+2G^{\prime}\right)=0$
(73)
$2S^{2}F_{y}^{\prime\prime}+e^{B_{1}}\left(S^{2}\left(2\left(G^{\prime}\left(\tilde{B_{1}}-F_{x}^{\prime}\right)+\tilde{G}^{\prime}\right)-\sinh(2G)\left(B_{1}^{\prime}\left(\tilde{B_{1}}-F_{x}^{\prime}\right)+\tilde{B_{1}}^{\prime}\right)-2\tilde{G}B_{1}^{\prime}\cosh(2G)\right)-3S\tilde{S}\left(B_{1}^{\prime}\sinh(2G)-2G^{\prime}\right)\right)+2S^{2}\left(\cosh^{2}(G)\left(\hat{B_{1}}^{\prime}-B_{1}^{\prime}F_{y}^{\prime}\right)+B_{2}^{\prime}\left(F_{y}^{\prime}-3\hat{B_{2}}\right)+\hat{G}\left(B_{1}^{\prime}\sinh(2G)-G^{\prime}\right)-\hat{B_{1}}B_{1}^{\prime}\cosh^{2}(G)-\hat{B_{2}}^{\prime}-4\hat{\phi}\phi^{\prime}\right)+S\left(6\hat{S}\left(B_{1}^{\prime}\cosh^{2}(G)-B_{2}^{\prime}\right)+2S^{\prime}F_{y}^{\prime}-8\hat{S}^{\prime}\right)+8\hat{S}S^{\prime}=0$
(74)

$12e^{B_{1}}S^{3}\dot{S}^{\prime}+e^{B_{1}+B_{2}}\left(S^{2}\left(2\cosh(G)\left(-\hat{G}\left(\tilde{B_{1}}+\tilde{B_{2}}-F_{x}^{\prime}\right)+\tilde{G}\left(\hat{B_{1}}-\hat{B_{2}}+F_{y}^{\prime}\right)+G^{\prime}\left(\tilde{F_{y}}+\hat{F_{x}}\right)-2G^{\times}\right)+2\sinh(G)\left(B_{2}^{\prime}\left(\tilde{F_{y}}+\hat{F_{x}}\right)+F_{y}^{\prime}\left(\tilde{B_{2}}-F_{x}^{\prime}\right)+\hat{B_{2}}\left(F_{x}^{\prime}-4\tilde{B_{2}}\right)-2{B_{2}}^{\times}+\tilde{F_{y}}^{\prime}-2\hat{G}\tilde{G}-4\hat{\phi}\tilde{\phi}+\hat{F_{x}}^{\prime}\right)\right)+S\left(2\sinh(G)\left(\hat{S}\left(F_{x}^{\prime}-4\tilde{B_{2}}\right)+\tilde{S}\left(F_{y}^{\prime}-4\hat{B_{2}}\right)+4S^{\prime}\left(\tilde{F_{y}}+\hat{F_{x}}\right)-8S^{\times}\right)-8\cosh(G)\left(\hat{S}\tilde{G}+\hat{G}\tilde{S}\right)\right)+8\hat{S}\tilde{S}\sinh(G)\right)+e^{2B_{1}+B_{2}}\left(S^{2}\left(2\sinh(G)\left(\tilde{G}\left(2\tilde{B_{1}}+\tilde{B_{2}}-F_{x}^{\prime}\right)-G^{\prime}\tilde{F_{x}}+\bar{G}\right)+\cosh(G)\left(2\left(-\left(B_{1}^{\prime}+B_{2}^{\prime}\right)\tilde{F_{x}}+\bar{B_{1}}+\bar{B_{2}}-\tilde{F_{x}}^{\prime}+\tilde{G}^{2}+2\tilde{\phi}^{2}\right)-2\left(\tilde{B_{1}}+\tilde{B_{2}}\right)F_{x}^{\prime}+2\left(\tilde{B_{1}}{}^{2}+\tilde{B_{2}}\tilde{B_{1}}+2\tilde{B_{2}}{}^{2}\right)+\left(F_{x}^{\prime}\right){}^{2}\right)\right)+S\left(2\cosh(G)\left(\tilde{S}\left(4\left(\tilde{B_{1}}+\tilde{B_{2}}\right)-F_{x}^{\prime}\right)+4\left(\bar{S}-S^{\prime}\tilde{F_{x}}\right)\right)+8\tilde{G}\tilde{S}\sinh(G)\right)-4\tilde{S}^{2}\cosh(G)\right)+e^{B_{2}}\left(S^{2}\left(2\sinh(G)\left(\hat{G}\left(-2\hat{B_{1}}+\hat{B_{2}}-F_{y}^{\prime}\right)-\hat{F_{y}}G^{\prime}+G^{\star}\right)+\cosh(G)\left(2\left(\left(B_{1}^{\prime}-B_{2}^{\prime}\right)\hat{F_{y}}-B_{1}^{\star}+B_{2}^{\star}-\hat{F_{y}}^{\prime}+\hat{G}^{2}+2\hat{\phi}^{2}\right)+2\left(\hat{B_{1}}-\hat{B_{2}}\right)F_{y}^{\prime}+2\left(\hat{B_{1}}{}^{2}-\hat{B_{2}}\hat{B_{1}}+2\hat{B_{2}}{}^{2}\right)+\left(F_{y}^{\prime}\right){}^{2}\right)\right)+S\left(8\hat{G}\hat{S}\sinh(G)-2\cosh(G)\left(\hat{S}\left(4\hat{B_{1}}-4\hat{B_{2}}+F_{y}^{\prime}\right)+4\hat{F_{y}}S^{\prime}-4S^{\star}\right)\right)-4\hat{S}^{2}\cosh(G)\right)+e^{B_{1}}\left(8S^{4}V(\phi)+24\dot{S}S^{2}S^{\prime}\right)=0$
(75)