Holographic Bubbles with Jecco: Expanding, Collapsing and Critical

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 circular 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.

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.

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 [1][2][3][4], 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. [5] for a review) or via the spinodal instability [6]. 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 [7,8], whereas the latter requires an understanding of out-of-equilibrium physics [9][10][11].
In [12] 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 [13,14], 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 [15,16], characteristic codes have shown remarkable stability. Indeed, the first ever long-term stable evolutions of moving black holes was accomplished with a characteristic scheme [17]. 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 [18][19][20][21][22][23]. There is an extensive literature on this and related subjectssee Winicour's Living Review [24] 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.
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 [56]. To the best of our knowledge, this is the first such freely available code (see however the PittNull code [57,58] for characteristic evolutions in asymptotically flat spaces, freely available and distributed as part of the Einstein Toolkit [59]). 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 =h = 1 units throughout.
2 Jecco: a new characteristic code for numerical holography

Equations
In this section we outline the theoretical background and equations that are implemented in Jecco. Our approach is similar to that of [52] and generalises the code presented in [34] to the 3+1 dimensional case. See also [24] for an overview of the approaches and codes used in the asymptotically flat setting.

Equations of motion and characteristic formulation
We consider a five-dimensional action consisting of gravity coupled to a scalar field φ with a non-trivial potential V (φ). The action for this Einstein-scalar model is where κ 2 5 is the 5D gravitational coupling constant, which in our units takes the value κ 2 5 = 8π. The resulting dynamical equations of motion read Our potential V (φ) comes from a superpotential W (φ) with the form and its explicit expression can be derived via In these equations λ 4 and λ 6 are freely specifiable dimensionless parameters related to the parameters φ M and φ Q used in e.g. [60,39] through This potential has a maximum at φ = 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 Λ for the dimension-three scalar operator O φ dual to the scalar field φ. The thermodynamical and near-equilibrium properties of this model were presented in [61,33,35] for λ 6 = 0 and in [60,39] for λ 6 = 0. Let us point out that even if here we will always make use of the particular potential (2.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 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 O φ . 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 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 → −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, We can recover the 2+1 system of [34] by setting (2.8) Table 1. Nested structure of the equations of motion.

Function(s) Combination
for non-trivial dependence only along the x or y direction respectively. The metric (2.7) is invariant under r →r = r + ξ(t, x, y) , S →S = S , x, y) . (2.9) Plugging the ansatz (2.7) into (2.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.2) as indicated, that takes the form where u ≡ 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 aṡ 11) which are necessary to obtain this nested structure. 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 (2.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 × 2 matrices. The equations themselves are lengthy and given in (B.2-B.8).
These equations need to be supplemented with boundary conditions specified at the AdS boundary u ≡ 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 φ(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.

Asymptotic expansions
The study of the near-boundary behaviour (u → 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 (2.10).
For this purpose, we start with an ansatz that is compatible with the AAdS condition, 1 (2.12) Substituting into equations (B.2-B.8) and solving order by order, we obtain where φ 2 is not the one in (2.12), but redefined as (2.14) Note that φ 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 ξ. The fact that the former is simply shifted by ξ under (2.9) means that we can identify s 0 with ξ and exchange them everywhere. We also need the expansions of "dotted" variables, defined in (B.1), which take the formḂ The function ξ(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 where b 14 (t, x, y), b 24 (t, x, y), g 4 (t, x, y), φ 2 (t, x, y), and ∂ t φ 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 φ(t, r, x, y) in equations (2.13b), (2.13c), (2.13d) and (2.13h). The functions a 4 (t 0 , x, y), f x2 (t 0 , x, y), f y2 (t 0 , x, y), and ξ(t 0 , x, y) should also be thought of as initial data, which can be freely specified. φ 0 is a parameter that must also be specified and corresponds to the energy scale Λ of the dual boundary theory.

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 ξ. Denoting with the g1 (g2) subscript the variables defined in the inner (outer) grid, the redefinitions that we choose to make are then Substituting these redefined variables into the system of equations (B.2-B.8), 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 Once again we note that functions B 1 , B 2 , G, φ, a 4 , f x2 , f y2 and ξ encode the freelyspecifiable 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.

Gauge fixing
To fully close our system we still need to fix the residual gauge freedom (2.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 Θ| r=r H = 0 at all times, where Θ is the expansion of outgoing null rays. Its explicit expression for the metric (2.7) is shown in Appendix C.
A simple way to enforce Θ| r=r H = 0 at all times during the numerical evolution is to impose a diffusion-like equation of the form with κ > 0, ensuring that the expansion Θ is driven towards the fix point Θ| u=u H = 0 as the time evolution runs, pushing the AH surface to u = u H = constant. The way to proceed is the following. We expand equation (2.18) using (C.6) and also the equations of motion for bothS andḞ 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 ∂ t ξ of the type which can be readily integrated with periodic boundary conditions in x and y.

Evolution algorithm
Having solved equations (B.2-B.8), we use the definition of the "dot" operator, cf. equation (B.1), to write and analogously for B 2 , G and φ. This tells us how to march these quantities forward in time. 2 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 [34] 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 ∂ t B 1 (t, u, x, y) = c(t, u, x, y)∂ u B 1 (t, u, x, y) + F B 1 (t, u, x, y) , (2.21) with c(t, u, x, y) = u 2 2 A(t, u, x, y) . (2.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 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 (2.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: i.e., we copy the modes leaving domain i + 1 to domain i.
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.
4. Obtain ∂ t a 4 (t 0 , x, y), ∂ t f x2 (t 0 , x, y) and ∂ t f y2 (t 0 , x, y) through (2.16). See Fig. 1 for a cartoon picture of the coordinates used and the evolution scheme (at constant x, y).

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 [34]. The result is: (2.25) For an SU (N ) gauge theory the prefactor κ 2 5 /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 T µ µ = −Λ O . (2.26)

Implementation
As already mentioned, we have implemented the algorithm of Sec. 2.1.5 in a new numerical code called Jecco [56], written in Julia [62]. 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 [34], 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 [63] 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 [64] for writing data.

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 are defined in the range [−1 : +1], and can be mapped to the physical grid by 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 [65]. 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 (2.10), where f represents the metric coefficients and scalar field φ. 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) where D uu , D u represent the derivative operators for a Chebyshev-Lobatto grid in the physical representation (see for instance [66] 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 (2.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 (2.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 (2.19) for ∂ t ξ is a linear PDE in x, y. To solve it, after discretizing in a N x × N y grid, we flatten the solution vector using lexicographic ordering and introduce enlarged differentiation matrices, which can be conveniently built as Kronecker productsD where D x , D y , D xx , D yy are the first and second derivative finite-difference operators. The cross derivative operator is built as a matrix product,D xy =D xDy . The PDE (2.19) then takes the algebraic form  [67] 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 (2.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.

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 [68]. 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 [69] 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 after each time step, where h x and h y are the grid spacings and σ 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 [70]). The complete scheme is where is the machine epsilon (for the standard choice of = 2 −52 , α = 36.0437) and γ is a tuneable parameter which we typically fix to γ = 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.

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 [6]. 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 (2.1) and (2.4) with the same value of the parameters (2.5) as in [12], namely 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 φ M and φ Q ) whose dual gravity solutions are completely regular even at zero temperature. The motivation for the choice (3.1) is that it leads to a sizeable bubble wall velocity, as we will see in Sec. 3.4. 2.5 Figure 2. Energy density as a function of temperature for the gauge theory dual to the holographic model (2.1)-(2.4) with parameters (3.1). 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.

Thermodynamics
The thermodynamics of the gauge theory can be extracted from the homogeneous black brane solutions on the gravity side (see e.g. [71]). 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. 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 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 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 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 is positive on these branches. At the temperatures T s and T s 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 Note that the characteristic scale for all the quantities above is set by the microscopic scale in the gauge theory, Λ, given holographically by Λ = φ 0 in terms of the leading term in the near-boundary fall-off of the scalar field in (2.13h).

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 [12] 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.
Phase-separated states are configurations in which the two homogeneous phases with energy densities E high and E 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,Ē. 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Λ 4 , whereas the energy density in the high-energy phase at T = T c has E high above 2Λ 4 , as given in (3.2). 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Ē in a box of fixed size is conserved upon time evolution. Therefore, phase-separated states with an average energy density in the region E s ≤Ē ≤ 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 [37,6]. To initialize the code we specify some φ 2 that is not too far away from the value of the thermal state and generate a simple bulk profile for the scalar, φ(t = 0, u), given by the truncated series in (2.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 (2.25). On top of it we add a sinusoidal whereā 4 is the value we obtained above, L x and L y are the lengths of the box, x mid and y mid correspond to the central point and δ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 δa 4 compatible with keeping the apparent horizon within our grid. We have found that δa 4 ∼ 10 −3 is a convenient choice. The state in Fig. 3(top) was generated following this method with φ 2 = 0.3Λ 3 . After a time tΛ = 300 the system has settled down to the configuration shown in the figure.
Phase-separated configurations with average energy densities in the regions 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 [39]. We take initial data corresponding to a phase-separated state withĒ in the spinodal region, and we modify it by increasing or decreasing the value ofā 4 so that the newĒ takes the desired value. We then let the system evolve. In a time around t = 100/Λ the system relaxes to a new inhomogeneous, static configuration. The phase-separated configurations in the middle and bottom panels of Fig. 3 have E s ≤Ē ≤ E high and were obtained with this procedure. The phase-separated states interpolate between the energy densities E low and E high . To construct initial data for bubble configurations that interpolate between two energy densities E B and E A we proceed as follows. Let f 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 PS = f PS (u, x, y), or one of the boundary functions such as a 4 , in which case f PS = f PS (x, y). We assume that the centre of the region with energy density E 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 E 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 .
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 lowenergy region and 1 at the edge of the box. As a consequence, f 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. 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 [12], we expect that the subsequent time evolution will quickly select a dynamically preferred state C = B inside the bubble. We could also multiply the bulk metric functions B 1 and B 2 in (2.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).

Critical bubbles
Consider a cylindrical bubble of radius ρ 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 Ac . 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 This follows from approximating the interface by a zero-width surface with free energy density γ, assigning a well defined pressure P Bc , and hence a free energy density −P Bc , 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 Ac is close to the phase transition temperature T c , namely for T Ac 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 Ac is sufficiently close to the turning point at T = T s , namely when T s T Ac . 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.  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  runtime was around 250h. We see that bubbles with initial radius Λρ c ≥ 3.75 eventually expand, whereas bubbles with radius Λρ c ≤ 3.69 eventually collapse. This means that the critical radius must be in between these two values. Substituting into (3.9) we then obtain an estimate for the surface tension γ. Thus, 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 -25 -   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 ρ 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 Λt = 20 so that the result is not contaminated by the fast-decaying, transient oscillations present around Λt = 0 in Fig. 5.
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. [7,8]). 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 1.5 2. 2.5 3. Figure 8. Relative difference between the energy density profiles at tΛ ∼ 20 of bubbles with different initial radii. We take as a reference the profile for a bubble with initial radius Λρ = 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 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 [39].

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 ρ/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.

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.
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 [12], this state is dynamically determined. While the initial configuration at Λ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 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, with E loc the energy density of the fluid in the local rest frame. The flow velocity v = u ρ /u 0 , with u ρ 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 Λ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. 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 0.31 can be extracted from Fig. 9 via a linear fit to the wall position of the form ρ wall (t) = ρ wall,0 + v wall t .   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 latetime form. Despite this, Fig. 11 shows that the wall has a fixed size set by the microscopic scale of the theory, Λ. In particular, the size of the wall does not grow with time, in contrast with the overheated region in front of the bubble wall.
In the case of planar bubbles, Ref. [12] 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 Fig. 6 with Λρ 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 E R as the value of the energy density at the maximum located right in front of the wall.

.3 and with the walls
of phase-separated planar and cylindrical configurations. Following [12], 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 E(ρ) but the combination (E(ρ) − E L )/(E R − E L ), with E L and 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 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%.

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 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 ρ/t, as described in e.g. [72]. 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 shift and ρ shift that we will define below. In other words, we define 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. 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 ρ = ρ 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 ρ wall (t) is given by (3.11) and ρ interface (t) = ρ interface,0 + v interface t . We then impose that, as soon as this regime starts, the values of ξ 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 shift and ρ shift so that the following two conditions are satisfied: In Fig. 14 we show the energy density and fluid velocity profiles for different simulation times as a function of ξ. In both plots we see two regions of fast change that separate three smooth regions. The first region of fast change occurs around ξ = v 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 ξ-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 ρ-coordinate decreases with time. However, going to the ξ-coordinate enhances this slope, since at late times dE/dξ ∼ t dE/dρ. 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 ξ-coordinate in this limit. The second abrupt region occurs at ξ 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 ξ-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 [72]. At very late times, when the bubble profile depends only on the scaling variable ξ, the ideal hydrodynamic equations lead to the following equation for the energy density and the velocity field of a cylindrical bubble where γ = 1/ √ 1 − v 2 is the Lorentz factor, c s is the speed of sound, E loc is the energy density in the local rest frame of the fluid, is the enthalpy density, and It is well known that the ideal hydrodynamic equation (3.15) 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 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 [12]. 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 [12], and the state in the overheated region was dubbed D 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 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 The state D is represented by a black dot in Fig. 2.

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 [34] to 3+1 dimensional settings and further allows, for instance, the usage of other choices for the scalar potential V (φ). The code is written in the Julia programming language [62] and is freely available at github https://github.com/mzilhao/Jecco.jl and Zenodo [56].
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 [6]. 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 [12] 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 [39]. 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 [73]. 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 [74,75] or to condensed matter systems [76][77][78], we expect that the first-principle nature of the holographic approach will shed new light on this problem too.

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 [34]. 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 =h = L = 1.
We note that we solve the equations of motion of our Einstein-scalar model (2.1) using the ingoing Eddington-Finkelstein gauge (equation (2.7)), which is a Bondi-like gauge, and the resulting PDE system is expected to be only weakly hyperbolic [15]. We thus restrict our tests to smooth data, where the effect of weak hyperbolicity is not expected to be manifested [15].
As mentioned in the main text, for the moment we have only implemented sharedmemory 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 φ 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: using the field redefinitions of Sec. 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 where δa 4 = 5 · 10 −4 , and the remaining free data functions (B 2 , G, φ, 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 [34], for the same setup. We use one inner radial domain spanning the region u ∈ [0, 0.1] discretized with 12 grid points, and another (outer) domain spanning the region u ∈ [0.1, 1.01] with 48 grid points. The transverse direction x spans x ∈ [−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 ξ 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.

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) [79]. 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 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 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 ξ than the one specified in Sec. 2.1.4. In particular, we impose the advection equation 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 ξ, which we initialize to where L x ≡ x max − x min . For such a configuration, the solution to equation (A.6) is and the exact solution of the metric function A is given by (A.1), where ξ is now provided by (A.8).
For the tests presented herein we have fixed For the numerical discretization we have employed one inner radial domain with 12 grid points (spanning the region u ∈ [0, 0.1]) and three equal-sized outer domains for the region u ∈ [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 (2.32) with σ = 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. 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 (3.6) withᾱ 4 = 1 and δα = 10 −3 , as well as φ 0 = Λ = 1, φ 2 = 0.3, andĒ = 1. The size of the box is L x Λ = L y Λ = 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 ∈ [0.1, 1.05]. We use second order finite difference operators and set σ = 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Ē of the above configuration. Notice that the runs performed for these tests reach tΛ = 21.69, whereas the top phase-separated profile of Fig. 3 corresponds to tΛ = 300 of the setup. Since we are using a low value for the dissipation parameter (σ = 10 −5 ), it is not possible to perform such long runs. The reason for this choice is that high values of σ seem to non-trivially affect the convergence properties of these configurations. However, we have checked that when performing the same runs with σ = 0.2, which is sufficient for long runs, the difference when comparing to the setups with σ = 10 −5 drops fast with increasing resolution, as illustrated in Fig. 19.    Figure 19. The difference inĒ for runs with low and high Kreiss-Oliger dissipation, with σ = 10 −5 and σ = 0.2, respectively. The norm of this difference is illustrated until Λ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 λ 4 = −0.25 and λ 6 = 0.1 for the model given by equation (2.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 ξ as follows where 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 φ 2 = 0.29819 and φ 0 = 1.0(= Λ). 3 Motivated by its near boundary behaviour, we initialize the scalar field to (A.10) 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 E.
We performed a total of 16 runs with energies evenly distributed in the interval E/Λ 4 ∈ [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 λ 4 and λ 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 [60,39]. Notice that the off-diagonal pressure, P xy , and the energy fluxes, J x and J y , are not shown as they are vanishing for these solutions.
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 [34], and was performed using a single core in a 16GB memory machine running for a few minutes.
We fix φ 0 = 1, λ 4 = 0.0025 and λ 6 = 0. We set as initial conditions G = 0 together with f x2 = f y2 = 0. a 4 is obtained from equation (A.9) with an average energy density E = 0.379686 and φ 2 = 0.0868357, which corresponds to the equilibrium value of the nonconformal uniform black brane with that same energy density. For this test ξ = 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: 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 [61] for different values of λ 4 . In particular, the boundary variables φ 2 , b 14 and b 24 will evolve in time according to with the mode 1, f (1) (t), being the longest lived (smallest ω i ).
We can obtain the different parameters from the data for f − f 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 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 lowestfrequency 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 [81,82] 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. [83]) truncated at first order in the hydrodynamic gradient expansion take the form where E loc is the energy density in the local rest frame of the fluid, u µ is the local fluid velocity, ∆ µν :=η µν + u µ u ν is the projector, P (E loc ) is the equation of state, and η and ζ are the shear and bulk viscosities, respectively. Moreover, we define σ µν := 2∇ µ u ν , where · 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 [39,61]). We consider as initial state a homogeneous black brane solution with small sinusoidal perturbations along x and y. The scalar bulk profile is again chosen to be given by equation (A.10), with φ 0 = 1.0(= Λ) and φ 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 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.
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 φ(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.  (Bottom-left) Difference between the T xy obtained from the numerical evolution and T 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.
T 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 con-figuration.

C Apparent horizon finder
In order to find the AH we need to compute the expansion of the outgoing null rays. We can construct the tangent vector to such outgoing rays using the ingoing null rays, n, together with the form perpendicular to the AH, s, s = N s (−∂ t σdt − ∂ y σdy − ∂ y σdy + dr) , n = −N n ∂ r , (C.1) from where we can compute the vector s by simply raising the indices. The normalisation factors, N s and N n , can be computed by imposing s 2 = 1 and s · n = −1/ √ 2. Combining these two vectors we can construct another vector tangent to outgoing trajectories, so that it is null, l 2 = 0, and properly normalised, l · n = −1. The expansion of these rays can be computed as where h µν = g µν + l µ n ν + l ν n µ (C.4) is the induced metric over hypersurfaces normal to both in-and out-going null rays. The AH location is given by the condition θ l = 0. Imposing it at a generic surface, r = σ(x, y), we obtain the following equation, (C.5) 2e B 2 (F y + ∂ y σ) S e B 1 cosh(G) G + G (F x + ∂ x σ) + e B 1 sinh(G) B 2 + B 2 (F x + ∂ x σ) + cosh(G) B 1 (F y + ∂ y σ) +B 1 − cosh(G) B 2 (F y + ∂ y σ) +B 2 − sinh(G) G (F y + ∂ y σ) +Ĝ + e B 1 sinh(G) S − 2S (F x + ∂ x σ) − cosh(G) S (F y + ∂ y σ) +Ŝ − 2e B 1 +B 2 (F x +∂ x σ) S e B 1 cosh(G) B 1 +B 1 (F x +∂ x σ) +cosh(G) B 2 +B 2 (F x +∂ x σ) +sinh(G) G +G (F x +∂ x σ) − sinh(G) B 2 (F y + ∂ y σ) +B 2 − cosh(G) G (F y + ∂ y σ) +Ĝ + e B 1 cosh(G) S + S (F x + ∂ x σ) − sinh(G) S (F y + ∂ y σ) +Ŝ where every function is evaluated at the r = σ(x, y) surface defining the AH. When the AH is located at constant radial surfaces, i.e. σ(t, x, y) = r = constant -which is what we impose to find the evolution equation for the gauge function ξ -equation (C.5) reduces to (C.6) Θ ≡ −2e B 1 +B 2 F x S e B 1 cosh(G) B 1 +B 1 F x +cosh(G) B 2 +B 2 F x +sinh(G) G +F x G − sinh(G) B 2 F y +B 2 − cosh(G) F y G +Ĝ + e B 1 cosh(G) S + F x S − sinh(G) F y S +Ŝ + 2e B 2 F y S e B 1 cosh(G) G + F x G + e B 1 sinh(G) B 2 + B 2 F x + cosh(G) B 1 F y +B 1 − cosh(G) B 2 F y +B 2 − sinh(G) F y G +Ĝ + e B 1 sinh(G) S − 2F x S − cosh(G) F y S +Ŝ x cosh(G)S + 3e B 2 F 2 y cosh(G)S = 0 One can check that when going to the 2+1 case, by imposing conditions (2.8), equation (3.17) of [34] is recovered.