General relativistic bubble growth in cosmological phase transitions

We use a full general relativistic framework to study the self-similar expansion of bubbles of the stable phase into a flat Friedmann-Lema\^itre-Robertson-Walker Universe in a first order phase transition in the early Universe. With a simple linear barotropic equation of state in both phases, and in the limit of a phase boundary of negligible width, we find that self-similar solutions exist, which are qualitatively similar to the analogous solutions in Minkowski space, but with distinguishing features. Rarefaction waves extend to the centre of the bubble, while spatial sections near the centre of the bubble have negative curvature. Gravitational effects redistribute the kinetic energy of the fluid around the bubble, and can change the kinetic energy fraction significantly. The kinetic energy fraction of the gravitating solution can be enhanced over the analogous Minkowski solution by as much as $\mathcal{O}(1)$, and suppressed by a factor as larger as $\mathcal{O}(10)$ in case of fast detonations. The amount of negative spatial curvature at the centre of the bubble is of the same order of magnitude of the naive expectation based on considerations of the energy density perturbation in Minkowski solutions, with gravitating deflagrations less negatively curved, and detonations more. We infer that general relativistic effects might have a significant impact on accurate calculations of the gravitational wave power spectrum when the bubble size becomes comparable to the cosmological Hubble radius, affecting the primary generation from the fluid shear stress, and inducing secondary generation by scalar perturbations.


Introduction
Phase transitions are characteristic phenomena predicted by many gauge field theories.They naturally arise as a consequence of a spontaneous breaking of a gauge symmetry [1][2][3].The standard Hot Big Bang model predicts at least two transitions [4,5], one at the electroweak scale (∼ 100 GeV) and one at the scale of quark confinement (∼ 100 MeV).In the Standard Model of particle physics, these transitions complete with a crossover [6,7].First order phase transitions are possible in many extensions of the Standard Model [8].The characteristic feature of a first order phase transition (FOPT) is the presence of a metastable quasi-equilibrium state just below the critical temperature.Quantum and thermal fluctuations push the field across a potential barrier and initiate the production of bubbles of the stable equilibrium phase, which eventually will expand and merge converting larger and larger portions of the Universe to the stable phase [9][10][11][12][13].
FOPTs enjoy great interest as a necessary component of electroweak baryogenesis [14][15][16], and the possibility of observing the stochastic gravitational wave background [17][18][19][20] they produce.The scenario for phase transitions becomes particularly compelling if a FOPT took place around the electroweak scale, because the gravitational wave signal that this transition is expected to generate might lie in the frequency window of future gravitational wave detectors as LISA (Laser Interferometer Space Antenna) [21].
The large energy density perturbations created by a strong phase transition might also seed a population of primordial black holes [22][23][24][25][26][27], which in turn would produce gravitational waves when they become gravitationally bound (see e.g.[28] for a review, also [29][30][31]).Gravitational wave production by primordial black holes has been under intense investigation since the first direct detection of gravitational waves from binary black holes merging by the Virgo-LIGO collaboration [32][33][34] in 2015.
Predictions of the gravitational wave power spectrum directly produced at a FOPT [8,[35][36][37] indicate that the estimated energy density in the gravitational waves Ω GW is an increasing function of (R ⋆ H e ), with R ⋆ the mean bubble spacing and H e the Hubble constant at some convenient reference time during the transition which we take to be the time when a fraction 1/e of the Universe remains in the metastable phase.In a conventional thermal phase transition, this is also the time of peak bubble nucleation rate, sometimes written H n [38].One is naturally led to consider the possible scenarios that maximize this factor, as these would be the transitions most likely to be observable.In this paper we focus on the study of the hydrodynamics of large bubbles, by which we mean bubbles of the stable equilibrium phase whose radial size R b is comparable to the Hubble radius, R b ∼ H −1 .This possibility is allowed if the nucleation rate of bubble of the stable phase is sufficiently small to allow the bubbles to expand for enough time before merging.We will see that gravitation affects the fluid velocity and energy density profiles around the bubbles, which will lead to changes in the gravitational wave power spectrum [36].
Furthermore, work on gravitational wave production from FOPTs has until now considered only primary gravitational waves, sourced by the transverse and traceless part of the energy momentum tensor [8,[35][36][37].Secondary gravitational waves sourced by scalar perturbations have been neglected, because they are expected to provide a negligible contribution as long as the bubbles are smaller than the cosmological horizon.However, in slow transitions, if the bubbles' size become comparable with the Hubble radius, we estimate that scalar induced (secondary) gravitational waves can play an important role.In light of these considerations, we are motivated to study the limit of large bubbles (R b H ∼ O(1)), that require a full general relativistic treatment.
Right after the nucleation, a pressure difference is established in the plasma filling the Universe at the interface between the symmetric and broken phases.Once the bubble is set in motion by the pressure difference at the interface, the friction between the bubble and the particles of the plasma induces a negative acceleration that slows down the expansion of the bubble [13,39].The competition between the friction and the driving force induced by the pressure difference determines whether, after a sufficient period of time, the bubble reaches a steady state with a constant terminal velocity, or continue to accelerate [40][41][42][43].Steadystate solutions are also possible without entropy generation [38,40,44].Phase transition fronts that expand at subsonic speed have been shown to be unstable to perturbations with wavelength λ c ≲ λ ≲ R b , where λ c is a critical length scale set by the surface tension and the kinetic energy difference across them [45].The effect of non-spherical perturbations is to increase the bubble phase boundary and thereby accelerate the bubble expansion [46] to supersonic speeds.We leave investigation of instabilities for future work.
On the macroscopic scale the cosmic plasma around the bubble admits a fluid description.Numerical analysis of the fluid profile of the bubble [36,47] as well as numerical simulations of FOPTs [48][49][50] have been carried on a flat background, and the effect of the expansion of the Universe has only been included as a rescaling of the energy-momentum tensor of the radiation fluid.This approximation relies of the fact that, if the transition completes in a period of time that is much shorter than the Hubble time, then the bubble size will remain much smaller then the Hubble radius, and, in a first approximation, the resulting gravitational waves behave as if in Minkowski space, after a conformal rescaling.As we said, this approximation breaks in slow phase transitions when the bubble size become comparable to the Hubble radius.The hydrodynamics of bubble in an expanding Universe described by a flat Friedmann-Lemaître-Robertson-Walker (FLRW) solution were studied in [51].However, the back-reaction of the bubble expansion on the bulk FLRW metric was neglected.
In this work we focus on the expansion of a single bubble during its steady state evolution and we wish to go beyond the approximation of Ref. [51] considering back-reaction between the spacetime metric and the bubble.Full general relativity has already been applied to study the collision between two vacuum bubbles [52].Here we wish to gain a full general relativistic understanding of the fluid dynamics around a single bubble, which will be important for accurate predictions of both primary and secondary gravitational wave production in slow phase transitions.
We find convenient to use geometrized units in this paper, where the Newton constant G and the speed of light c are set to unity, G = c = 1.

Phase transitions in the Early Universe
The content of the early Universe before the phase transition can be described on large scales as a perfect fluid in thermal equilibrium with energy momentum with the enthalpy density w = e + p, where e and p are the energy density and pressure respectively; u µ is the four-velocity of the fluid, and g µν the background metric.The simplest way to model a system with a first order phase transition is to add a scalar field ϕ coupled to the plasma which acts as an order parameter for the transition.The energy momentum tensor of the full scalar-fluid system is then [53] T where w, e and p depend on both the temperature T and the scalar field ϕ.We suppose that ϕ = 0 in the high-temperature phase.

Free energy and effective potential
All the thermodynamic properties of the scalar-fluid system are encoded in the pressure p(ϕ, T ), or equivalently in the free energy density F(ϕ, T ) = −p(ϕ, T ) [36,54,55].Given the free energy, one can identify the potential of the zero temperature theory V 0 (ϕ) and the effective potential V eff (ϕ, T ), arising from thermal corrections and interactions between the scalar field ϕ and the particles in the plasma as Both the energy and enthalpy density derive from the free energy according to [51] At any temperature T , the equilibrium states of the system are found as the states of minimal free energy In the general scenario [36], at high temperature there is only one global minimum of the free energy, which uniquely defines the equilibrium state of the system.As the Universe expands, the temperature drops, and the free energy features a new local minimum.At the critical temperature T c , the two minima become degenerate in the free energy density F. Below the critical temperature the system is initially in a metastable quasi-equilibrium state (sometimes called false vacuum).An effective potential barrier separates states with homogeneous ϕ and prevents the system from classically rolling to the new stable equilibrium state (sometimes called true vacuum), which is now the global minimum for the free energy density.Nonetheless, quantum and thermal fluctuations allow the scalar field ϕ to jump through the barrier via inhomogeneous configurations and to reach the new equilibrium state.Since the transition between the two equilibrium state is sometimes accompanied by a (spontaneous) symmetry breaking, it is common to refer to the system in the false vacuum as symmetric phase, and the system in the true vacuum as broken phase.We will adopt here the same nomenclature, even though we do not require any symmetry breaking to take place.After nucleation, the bubbles whose radius exceeds a critical value R c [3], start growing until they eventually merge and and the stable phase percolates.
The bubble wall separates two distinct phases with different equations of state.In the interior of the bubble the scalar-fluid system is in the broken phase, while outside the bubble the system is in the old symmetric metastable phase.Between the two phases there is a wall, that is a region where the order parameter ϕ smoothly interpolates between the values in the two vacuum states, over a microscopic distance of order T −1 c .The wall itself carries a surface tension σ that is determined by the effective potential, and is of order T 3 c .We mentioned in the Introduction 1 that our focus is on large bubbles, that is bubble whose broken phase has a radial size R b comparable with the Hubble radius H −1 .As the bubble grows large, we are less and less interested in resolving the microphysics at the wall.Denoting with W the proper width of the wall, the difference in scale W/R b ≪ 1 motivates us to adopt the thin wall approximation, that is to consider the wall as an infinitesimally thin surface, and to treat the order parameter ϕ as a step function with a discontinuity at the wall.

Equation of state
After nucleation, the pressure difference at the interface between the two phases forces the bubble to expand.To describe this process we need an equation of state (EOS) that links the thermodynamic variables p(ϕ, T ), e(ϕ, T ).For the purpose of this paper, we will consider a simple scenario where the scalar-fluid system behaves as a barotropic perfect fluid in both phases, with equation of state where ω is the EOS parameter, taken to be constant.For simplicity we will fix the EOS parameter in the broken phase to ω − = 1/3, and consider different values of the EOS parameter ω + in the symmetric phase.We further assume that ω ± evolves slowly with the temperature T throughout the phase transition.For numerical calculations, this allows us to consider, in a first approximation, the EOS parameter ω ± to be constant during the expansion of the bubble.In order for the bubble to expand we need ω + < ω − ; this condition reflects the fact that some particle species acquire a mass in the broken phase and become non-relativistic.Qualitatively speaking, the number of particle species that become non-relativistic in the broken phase determines the strength of the phase transition.To quantify this concept, it is useful to introduce the trace anomaly θ θ ≡ 1 4 (e − 3p). (2.8) We can then define the boundary strength parameter of the phase transition α + as [36] α where the plus (minus) sign indicates that the quantity must be evaluated just ahead (behind) the interface.We also introduce the related idea of the global strength parameter where all quantities are evaluated at the same temperature, typically that at which the bubble nucleates, and the subscripts s or b indicates that the quantity is evaluated in the symmetric or broken phase.The two quantities are not equal in general, as the temperature in front of the wall can be higher than the nucleation temperature, and is different from the temperature behind the wall.However, in our simplified equation of state, the quantities θ and w scale the same way with temperature, and α + = α.

Mean bubble separation
Once the temperature of the Universe drops below the critical temperature T c , spherical bubbles of the new phase start nucleating from thermal and quantum fluctuations of the field ϕ.Due to the random nature of quantum and thermal fluctuations, the appearance of bubbles of the new broken phase is a stochastic process that can happen everywhere in the Universe.After nucleation, the bubbles start growing until they eventually merge and percolate, and the Universe converts to the new phase.We introduce a new length scale R ⋆ which defines the mean bubble centre spacing, so that where n b is the number density of bubbles.The length scale R ⋆ is important for observations, since it sets the peak frequency of the gravitational waves generated by the phase transition [8,35,36,[54][55][56].The interest of this work is in the limit of large bubbles (R ⋆ H ∼ O(1)), which requires a full general relativistic treatment.

Bubble growth in Minkowski spacetime
In this section we briefly review the results on bubble dynamics on a Minkowski background.This might help the reader to establish a comparison with the general relativistic treatment of the following sections.The bubble wall is the interface between the symmetric and broken phases, and below the critical temperature there is a pressure difference across the wall, which drives the expansion of the bubble.The microscopic description of the expansion is given by a solution of the effective field equation for the scalar field ϕ [53] □ϕ where η(ϕ, T ) encapsulates the microscopic out-of-equilibrium dynamics leading to friction between the expanding bubble and the particles of the plasma [40,41,57].The friction function η(ϕ, T ) is difficult to calculate and is often treated as a free constant parameter [53].Its main function is to determine the asymptotic speed of the wall as the bubble grows large.
The steady flow solution is usually reached after a first period of accelerated expansion; once the wall reaches its terminal velocity, the bubble continues to expand at constant speed until colliding with other bubbles in the Universe.In this paper we will focus on describing bubbles during their steady flow state.Equation (2.12) is non trivial only at the wall, where the field ϕ interpolates between its values in the broken and symmetric phase.By the equivalence principle, there must be a system of coordinates that is comoving with the wall and in which the metric is locally Minkowski; we call this system the Gaussian normal frame [58].Therefore, in Gaussian normal coordinates the coupled scalar field -fluid equation (2.12) is the same in the general relativistic case and in the flat Minkowski approximation, and so should be relationship between the friction η(ϕ, T ) and the fluid speed in the wall frame should also be the same [53].
Let us consider a flat static Universe.In polar coordinates x µ = (t, r, θ, ϕ), the metric is and the fluid four-velocity with u ≡ dr/dt the outward fluid three-velocity.Once the bubble reaches its steady state, the wall moves with a constant velocity u w with respect to the Universe rest frame (Universe frame).In the rest frame of the wall (wall frame) the conservation of the energy-momentum at the wall implies where v is inward fluid velocity in the wall frame2 and γ = (1 − v 2 ) −1/2 ; for convenience we introduced the bracket notation (2.16b) The subscript positive (negative) sign convention, consistently with equation (2.9), denotes quantities evaluated just outside (inside) the interface.Equations (2.15) admit an analytical solution v ± (α + , v ∓ ) for a given strength parameter α + (2.9).Physical solutions of (2.17a) demands the positive sign when v − > 1/3 and the negative sign when v − ≤ 1/3.Analogously, physical solutions of (2.17b) demand the positive sign for v + < 1/3 and negative sign for v + ≥ 1/3.This means that both velocities inside and outside the wall can be either subsonic (deflagration) or supersonic (detonation).The strength parameter α + for subsonic solutions has an upper bound, α + < 1/3, given by the requirement that v + must be positive.The choice of sign in (2.17) depends on the wall speed in the Universe rest frame and the sound speed c s of the broken phase where ϕ(T ) is the equilibrium value of the field.As stated in Section 2.2, we suppose that the effective potential is such that the scalar-fluid system is a barotropic perfect fluid in thermal equilibrium, so that the pressure is a function of the energy density only, p = p(e), and we compute the speed of sound as We will comment more about the sign choices in (2.17) later.The profiles of the fluid velocity and enthalpy density as functions of r and t are found by projecting the conservation equations onto the fluid four-velocity u µ and its spacelike orthogonal vector ūµ ≡ γ u (−u, 1, 0, 0): (2.20) Let us now choose our time coordinate t such that the bubble nucleates at t = 0. Since the bubble is expanding on a flat background, there is no physical length scale involved in the problem.We then expect the solution of conservation equations to manifest a self-similar behaviour.This means that the physical variables of the system, namely the fluid threevelocity and enthalpy density, depend on the radial and time coordinates only through the combination ξ ≡ r/t.For the purpose of numerical calculation, it is convenient to write the system of equations of motion in a parametric form, adding a new parameter τ such that ) ) where we have defined the outward fluid velocity in a frame that is moving outward at speed ξ.The system (2.21) has three fixed points in the phase space (ξ, u); one at the origin (0, 0), one on the light cone at (1, 1) and one at (c s , 0).Technically speaking, the point (c s , 0) is an unstable improper node [47], and all the trajectories approach it tangentially to the axis (ξ, u = 0).For a given wall speed u w , the system (2.21) uniquely provide the bubble profile once we impose boundary conditions at the origin of coordinates, at ξ = u w , and at infinity.At ξ = 0 we expect u = 0 because of the spherically symmetric configuration.On the other hand, ξ = 1 represent the position of the light cone.Since information cannot travel faster than light, we also expect the velocity field to reach zero before the light cone, that is u(ξ = 1) = 0.In between we use the junction conditions (2.15) to join interior and exterior solutions.
Depending on the way we meet the boundary conditions, we distinguish three cases: • Subsonic deflagration: The bubble expands in a deflagration mode when the inward fluid velocity in the wall frame behind the wall is slower than the sound speed of the broken phase, v − < c s− .The fluid is at rest inside the wall in the Universe frame, hence u w = v − .We must choose the negative sign in (2.17a) in order to ensure that the fluid flows into the bubble, i.e. u + = µ(v + , ξ w ) > 0. In front of the wall, the fluid velocity profile shows a compression wave that decays until the fluid reaches, via a shock, the configuration at rest with the Universe frame (u = 0).
• Detonation: A detonation is characterized by a supersonic inward fluid exit speed in the wall frame v − > c s− , while the fluid is at rest with the Universe frame everywhere outside the wall u + = µ(v + , ξ w ) = 0.The velocity profile shows a rarefaction wave behind the wall that reaches the improper node at (c s− , 0).The condition v − > c s− is satisfied choosing the positive sign in equation (2.17b), and it implies that there is a minimum possible value for v + given by v CJ ≡ v + (α + , c s− ).This is the Chapman-Jouguet speed; in our case3 we have c s− = 1/ √ 3, for which the Chapman-Jouguet speed is Figure 1: Profiles of outward fluid velocity v (top panels) and enthalpy w normalized to its value outside the shock w 1 (bottom panels) against the self-similar coordinate ξ = r/t with the constant sound speed equation of state in the case of a bubble expanding on a flat Minkowski background for deflagration (left panels), hybrid (central panels) and detonation (right panels).In the constant sound speed model ω − (ω + ) is the speed of sound squared in the broken (symmetric) phase, and α is the transition strength parameter (2.10), which is temperature-independent in this model.Therefore in a detonation u w > v CJ .
• Supersonic deflagration (hybrid): A hybrid is a solution at wall speeds in between the two aforementioned cases; more precisely, in a hybrid solution the inward fluid exit speed in the wall frame is precisely v − = c s− , and the wall speed in the Universe frame is bounded by c s− < u w < v CJ .The velocity profile shows both a rarefaction and a compression wave behind and outside the wall respectively.The rarefaction wave ends at the fixed point (c s− , 0) as in a detonation, while the compression wave ends with a shock as in a deflagration.At the wall, we choose the negative sign in (2.17a) in order to ensure u + > 0. Outside the shock the fluid is at rest with the Universe frame u sh + = 0. Hybrids may not be stable solutions for some equations of state [57].
The solutions for the three different cases have been studied and discussed in [47,53,55].In Figure 1 we show an example of fluid velocity and enthalpy profiles where the sound speed c 2 s ≡ dp/de is constant in the two phases (see Eq. (3.9) and (3.10)) for a mathematical definition of the EOS).Plots were produced using a modification of the code snippet provided in Ref. [59].In the following we will use this class of Minkowski spacetime solution with constant sound speed EOS as a reference to compare to the gravitating bubble solutions.
3 General relativistic description of bubble growth

Coordinate frame and equations of motion
We will consider that, during the expansion of the bubble, both the inner and outer regions are described by a spherically symmetric spacetime.We parametrize the expanding region as a foliation of concentric spherical shells labeled with a radial coordinate r comoving with the fluid [60].The metric is with a, b and R functions of t and r only.This metric was previously studied in the context of gravitational collapse [61], droplets in the quark-hadron transition [62], and primordial black hole formation [60].As the radial coordinate is comoving with the fluid, the radial component of the fluid 4-velocity is zero, as discussed in [63].In particular The time coordinate can be fixed by specifying the metric coefficient a on a timelike curve [63].
We chose to set a(0, t) = 1, so that the cosmic time t refers to the time measured by an observer at the centre of the bubble.We discussed in Section 2.4 that steady state solutions with constant wall speed always exist in Minkowski spacetime.Defining a wall speed in general relativity can be more difficult, and we rather characterize steady state solutions as bubbles whose wall is located at R w = ξ w t, with ξ w a constant.The position of the wall in comoving coordinates r w (t) is determined implicitly from R(r w (t), t) = ξ w t.We will see that the solutions we find have constant fluid speeds at the wall, similar to those in the Minkowski solutions.
Using the notation of [61], we define the proper time and radial derivatives where ūµ is a radial spacelike vector orthogonal to the fluid four-velocity u µ , and we introduce the quantities The former represents the radial component of the fluid four-velocity in a (non-comoving) Eulerian frame where R is the radial coordinate [64] (hence R = A/4π, with A the constantr spatial surface, is sometimes called the Schwarzschild radial coordinate).The quantity Γ can be interpreted as a generalization of the fluid Lorentz factor.Let us consider the invariant quantity ∆, defined as which in turn defines a coordinate invariant quantity M that we will interpret as mass.
Using the Einstein equations we find which gives to M the interpretation of the mass-energy contained inside a shell of radius R.
Finally we introduce two new dimensionless quantities, the surface energy Ω on a shell of radius R, and the gravitational potential Φ at R as In terms of U, Ω, Φ the system of Einstein and hydrodynamic equations is [60,64] To close the system, we require an equation of state p = p(e) .Remembering the discussion in Section 2.2 we consider the simple EOS with ω(t, r) a constant, but different, value in the two phases.The time dependence in ω(t, r) is caused by the expansion of bubble wall into the symmetric phase.We can write the EOS parameter ω as a distribution defined on the whole spacetime as with Θ(r) the Heaviside distribution and r w (t) the position of the wall in comoving radial coordinate at time t.Analogously, the speed of sound can be expressed as This equation of state is a special case of the constant sound speed model [59], with zero vacuum energy in both phases.For simplicity we fix ω − = 1/3 in the broken phase, while the EOS in the symmetric phase ω + is a free constant parameter.Fixing ω − = 1/3 implies that the boundary transition strength parameter α + is independent of the temperature, and given by The temperature-independence, and the fact that θ − = 0, also means that the global transition strength parameter α = α + .Therefore we do not need to make a distinction between the two in this special equation of state.

Self-similar assumption
Motivated by the results in Minkowski space (see Section 2.4), we seek self-similar solutions of the physical dimensionless quantities U (ξ), Ω(ξ), and Φ(ξ).Self-similar collapsing solutions are already known in the study of primordial black hole formation [60].We here adopt a similar formulation and introduce the self-similar variable 4ξ ≡ R(t, r) t .
(3.13) which greatly simplifies the set of Einstein equations to a system of three differential equations where Γ(ξ) is given by We notice that the system (3.14) is autonomous and enjoys the shift symmetry ln ξ → ln ξ +λ, with λ an arbitrary constant.Furthermore, introducing the self-similar variable and combining equations (3.8b) and (3.8c), we find that the metric coefficient a(t, r) also has a self-similar behavior given by The junction conditions in Minkowski are stated in the frame at rest with the wall (see Section 2.4).Analogously we want to define here the notion of a fluid three-velocity for some observer moving outward with the wall.As the wall is assumed to be stationary at a fixed position ξ w , we find convenient to introduce a set of constant-ξ observers, that is observers that are moving outwards with respect to the fluid rest frame with constant ξ, defined by the condition dξ/dt = 0.The worldline of these observers in the (t, r) coordinates is described by a 4-velocity (3.17) We can re-express the definition (3.17) of v using the expression (3.15) for a as This is a more convenient expression to use to evaluate junction conditions and to find critical points.
The relative speed v rel between two observers moving with 4-velocities u µ 1 and u µ 2 can be computed as shown in [65] with n µ 2 the radial orthonormal vector to u µ 2 , and a • b = g µν a µ a ν .One can verify that with u µ the fluid 4-velocity (3.2) and V µ ξ the 4-velocity of a constant-ξ observer (3.16).Hence we interpret the 3-velocity v as the inward fluid velocity for a constant-ξ observer.At ξ = ξ w , constant-ξ observers are comoving with the wall, and v(ξ w ) is the inward fluid speed in the wall rest frame.
The simplification from the system (3.8) of partial differential equations to the system (3.14) of three ordinary differential equations comes with some assumptions.First of all, we implicitly assumed that the bubble nucleates at t ′ = 0. Otherwise, one should define the self-similar variable as ξ = R/(t−t ′ ), and explicit time dependent factors t ′ /t would appear in the system (3.14),spoiling the self-similar behavior.The EOS parameter ω defined in (3.10) is consistent with H ∝ t −1 .Steady state solutions R w = ξ w t imply R w H ∼ ξ w .Therefore, self-similar steady state solutions of (3.14) represent bubbles whose radial size is a constant fraction, approximately given by ξ w , of the Hubble radius H −1 .In other words, self-similarity is compatible with a steady flow.However, in a real phase transition the EOS parameter is expected to decreases as the Universe cools down.Therefore, the two pictures are compatible if we assume that ω evolves slowly in time and the time-dependent solutions evolve slowly with ω.

Kinematic and geometric properties of an expanding bubble
The kinematic properties of fluid worldlines encode important information about the behavior of the congruence of the fluid four-velocity u µ .These properties can be read from the irreducible tensorial decomposition of the covariant derivative of the fluid four-velocity [66,67] with a µ ≡ u ν ∇ ν u µ the fluid four-acceleration and hµν = g µν + u µ u ν the projector onto the hypersurface orthogonal to u µ , and are the vorticity tensor, shear tensor, and local expansion scalar respectively.Vorticity describes a rotation of the congruences around a fixed axis; a non vanishing shear introduces distortion on the congruences without changing the volume of spacetime they cover; the expansion scalar instead describes the volume changing of the congruences.We compute these quantities in the case of a spherically symmetric expanding bubble.Vorticity is always zero around a single bubble, while for the shear and expansion scalar we have respectively In a FLRW Universe, U = RH and d ln U/d ln ξ = 1, so that σ = 0 and θ = 3H.
The curvature properties of the spacetime can be described by the mean of coordinate invariant contractions of the Riemann and Weyl tensor [68].For the purposes of this paper we will be interested only in the Ricci curvature scalar R (4) and Kretschmann scalar K (4) R (4) ≡ R ρ µν σ g µ ρ g νσ , (3.25a) with R ρ µν σ the Riemann tensor.Using the Einstein equations (3.8), after some algebra we find The Ricci scalar contains informations about the energy content of the spacetime.In particular, in our framework R (4) vanishes if the Universe is empty or if its energy content is described by an EOS parameter ω = 1/3.The Kretschmann scalar instead is a useful probe to use to investigate spacetime singularities.Constant-time slices of the spacetime provide an easier visual understanding of the geometry of the Universe.We compute the Ricci curvature of the spatial sections as where i, j, ℓ, k ∈ {1, 2, 3}.In Section 7.4 we will show how the Ricci curvature of fixed-time slices is connected to the energy under-density inside the bubble.

Matching conditions
For a given wall position, the system (3.14) is sufficient to solve the evolution of the bubble in both phases.At the wall, the order parameter ⟨ϕ⟩ moves from the symmetric to the broken phase in a continuous way, obeying (2.12) [53].However, as the size of the bubble grows larger and larger, we are not interested in resolving the width of the wall, and we can replace the continuous transition of the order parameter with a discontinuous jump at the wall.The interface between the the symmetric and broken phase is then a singular surface layer described by a δ-function singularity in the stress-energy tensor T µν .The set of rules that connects quantities on the two sides of the layer are given by the Israel junction conditions [58].

Junction conditions
In this section we give a brief summary of the Israel junction conditions, while a deeper discussion can be found in the Appendix A. Let y µ = (t, r, θ, ϕ) be the coordinate system of the global space-time, and x a = (τ, θ, ϕ) a coordinate system on the interface of the transition (namely the surface defining the wall).Such surface is defined by a constraint C(y µ ) = 0, and defines a tetrad of orthogonal vectors (n µ , e µ τ , e µ θ , e µ ϕ ) foliating the global space-time, with n µ the orthogonal vector to the hypersurface, and e µ a the tangent vectors to it.The two Israel junction conditions are [58,69] [h ab ] ± = 0, (4.1a) with D the covariant derivative induced by the intrinsic metric.Solving equations (4.1) and (4.5) simultaneously is a difficult problem, and in general requires an involved numerical analysis.For the current purposes, we can apply some simplifications given the difference in scales between the size of the bubble R b and the width of the wall W , and expand the junction conditions at zero order in the expansion parameter ϵ = W/R b ≪ 1 First, we suppose that the energy density e is everywhere the same order of magnitude, including the region in the domain wall.Combining the Friedmann equation in geometrized units H 2 ∼ e and the requirement that the bubble remains sub-horizon R b H ≲ 1, we estimate the energy density scale of the bubble as e ≲ R −2 b .By dimensional analysis we estimate the surface stress energy tensor on the hypersurface as S ∼ eW ≲ ϵR −1 W ≪ R −1 b .On the other hand, since R b is the only length scale that the intrinsic metric is sensitive to, we estimate D, K ∼ R −1 b .Therefore the right hand side of equation (4.1b) can be neglected compared to K.Moreover, the left hand sides of equations (4.5) can be estimated in this way as DS ∼ KS ∼ eϵ ≪ e.On the other hand the non-zero components of the energy momentum tensor T ∼ e, and, as before, at zeroth order in ϵ we neglect the right hand side of equations (4.5).
The full system of junction conditions that we consider is then [T n a ] ± = 0, (4.6c) The last two equations in (4.6) are analogous to the junction condition (2.15) in Minkowski, while the first two equations (4.6), trivially satisfied in Minkowski, are needed to quantify the discontinuity in the additional degrees of freedom in the metric g µν .Considering the constraint C(x µ ) = R − ξt = 0 we find e µ θ = (0, 0, 1, 0) (4.7c) e µ ϕ = (0, 0, 0, 1) (4.7d) with v the inward fluid velocity in the constant ξ frame (3.17) and γ = 1 − v 2 −1/2 .At ξ = ξ w , the constant ξ observer is at rest with the wall, so that v(ξ w ) is understood as the inward component of fluid velocity in the wall rest frame.Inserting (4.7) and the metric (3.1) in the junction conditions (4.6) we find (see Appendix A for a detailed derivation of (4.8a) and (4.8b)) [R] ± = 0, (4.8a) [Φ] ± = 0, (4.8b) ) We can understand the continuity of R from the fact that 4πR 2 is the proper area of a spherical shell at fixed time t, and two observers on the different sides of the interface must agree on the measurement of the proper area.The continuity of Φ instead derives from the fact that the mass M is computed through (3.6) as the integral of the energy density Ω from the origin of coordinates to R. A finite jump in Ω at the interface would then produce a discontinuity in the first derivative of Φ, without spoiling the continuity of the function itself.At first order in the limit ϵ = W/R b ≪ 1, the junction conditions (4.8c) and (4.8d) share the same functional structure of the ones in a flat Minkowski space-time [54,55].Hence, solutions (2.17a) and (2.17b) apply for v(ξ w ) with strength parameter α Now, from (4.8c), we can extract the energy density outside the interface Finally, inverting (3.17), we get a quadratic equation for U + in terms of the other variables.The positive root is while the negative root for U + can be neglected since it represents an inflowing solution in the non-comoving Eulerian frame.

Asymptotic FLRW solutions
At the centre of the bubble r = 0 and we expect that the Schwarzschild circumference coordinate R and Ω tend to zero.Moreover, the spherical symmetry of the system imposes the vanishing of the fluid velocity U at the centre of the bubble, while the physical interpretation of M as a mass inside a shell of radius R suggests that M (and then Φ) tends to zero in this limit.This motivates us to parametrize U = jξ γ , Φ = kξ α , Ω = ℓξ β in a sufficiently small neighborhood of the origin ξ = 0, with α, β, γ, j, k, ℓ unknown constants.The equations of motion (3.14) with the gauge condition a(t, r = 0) = 1 allow to restrict the solutions to a one parameter family ) Since R ≥ 0, the weak (WEC) and strong (SEC) energy conditions [68,69] demand k > 0.
The value of k controls the amount of curvature of the spatial sections near the origin.Using the asymptotic solution (4.12) in the expression (3.27), we find This equation highlights the presence of a unique value k F = 2/9(1 + ω − ) 2 for which the spatial sections around the origin are flat.For k > k F the spatial curvature is positive, while for k < k F the spatial curvature is negative.As the Minkowski solutions manifest a deficit in the energy density in the interior of the bubbles, following this analogy we are led to consider only the values of k smaller than k F , that is Once the boundary conditions at the origin are fixed, the system (3.14) can be integrated out to the wall, where the junction conditions (4.8) provide the matching between the two phases and the initial conditions for the solutions outside the wall.Outside the wall the system can be integrated until it is matched via a shock to a spatially flat Friedmann-Lemaître-Robertson-Walker (FLRW) Universe.The outermost part of the solution is described by the metric (3.1) with Here, the subscript F denotes quantities evaluated on a FLRW background.Then the Friedmann equations imply and given the solution we can directly compute

.20)
Outer constant-ξ observers live in a FLRW Universe and are described by the 4-velocity (3.16) with the additional conditions (4.18).The outer observers can signal inwards to the inner constant-ξ observer inside the bubble about the Universe they live in; in particular they can communicate the boost that inner observers should apply in order for them to be comoving with what would be the cosmic fluid in an unperturbed FLRW Universe.There then exists a class of observers U µ that are comoving with what would be cosmic fluid in a FLRW Universe defined by where we denoted the boost transformation and v F the value of v at ξ in a FLRW Universe with equation of state parameter After some simple algebra we find the relative velocity between a constant-ξ observer V µ ξ and what would be a constant-ξ in an unperturbed FLRW Universe By Lorentz symmetry, this coincides with the relative fluid speed that the observers V ξ and V F ξ would measure.Hence, u represents the local departure from a FLRW Universe that outer constant-ξ observers can communicate to the inner constant-ξ observers.By construction, u vanishes outside the shock, where the Universe is described by the flat FLRW solution.
5 Properties of self-similar solutions

Fixed points
In the case of a Minkowski background, see Section 2.4, the equations of motion (2.21) for the fluid velocity and energy density show three fixed points [51,54,55].In analogy with the Minkowski case, in order to look for fixed points, it is convenient to introduce an auxiliary parameter τ such that ) ) (5.1d) A fixed point is found whenever dξ dτ = dU dτ = 0, as at this point the right hand side of the full system vanishes identically.We find three different solutions to the fixed point conditions.One fixed point is the trivial solution (ξ, U, Ω, Φ) = (0, 0, 0, 0).A second solution is ) where the star denotes variables that satisfy the fixed point conditions dξ dτ = dU dτ = 0 and c s− is the sound speed of the broken phase (3.11).In light of the expression (3.18) of the inward fluid speed for a constant-ξ observer, we can replace the first condition (5.2a) with the more convenient expression v ⋆ = c s− . (5. 3) The third fixed point allowed by the system (5.1)satisfies the fixed point conditions with v ⋆ = c s− and Ω ⋆ = Φ ⋆ .This class of solutions is compatible either with an EOS parameter ω − = −1 or with a fluid flow U that is positive around the origin, crosses U ⋆ = 0 and becomes negative after the fixed point.The latter situation occurs when k > k s , that is for solutions with positive spatial curvature R (3) around the origin.We interpret this solution as the time reversal of the ones studied in [60].However, since this fixed point does not play any role in our solutions, we will not discuss this further.Due to the ln ξ shift symmetry of the equations of motion (3.14), we can constrain the study of fixed points to the phase-space (U, Ω, Φ).Indeed the conditions dξ dτ = dU dτ = 0 are met regardless of the value of ξ.Conditions (5.2) define a codimension-2 hypersurface in the phase space (U, Ω, Φ), i.e. a line, that we will denote as fixed line, or line of fixed points.⃗ Y ⋆ ≡ (U ⋆ , Ω ⋆ , Φ ⋆ ) is then a line of fixed points for any ξ, and (ξ, ⃗ Y ⋆ ) is a surface of fixed point in the four-dimensional phase-space that we call a fixed surface.Similarly, we will refer to the solution of the Einstein equation ⃗ Y (τ ) ≡ (U, Ω, Φ) as the trajectory of the solution in the phase space.Figure 2 shows that, integrating outward from the origin with boundary conditions (4.12) for a given parameter k, the trajectory of the solution ⃗ Y k (τ ) always ends up on a fixed point.We denote the value of ξ at the end point of ⃗ Y k (τ ) with ξ ⋆ (k).The ln ξ shift symmetry of the equations of motion is broken by the initial conditions (4.12) with the choice a(ξ → 0) = 1.As mentioned above (see Sections 3.1 and 4.2), this choice fixes the time coordinate t and the scale of ξ = R/t.Therefore the scale of ξ ⋆ (k) is set by the condition a(ξ → 0) = 1, and the set of points γ ⋆ (k) ≡ (ξ ⋆ (k), ⃗ Y ⋆ (k)) is reduced to a curve in the four dimensional phase space.As one can see from Figure 3 and 2, varying the parameter k the end point of the trajectory of the solution ⃗ Y moves along the line of fixed points γ ⋆ (k).In other words, k parametrizes the line of fixed points γ ⋆ (k).

Subsonic deflagrations
In a deflagration the fluid speed behind the wall in the constant ξ frame is always smaller than the sound speed in the broken phase, i.e. v − < c s− .Hence the trajectory of the solution never reaches the line of fixed points.At the wall, in order to match on to an inflowing solution U + > 0, the negative sign in (2.17a) must be chosen.The fluid outside the wall then moves slower than inside, forming a compression wave in front of the wall.In order for the bubble to expand over the symmetric phase, we need to verify that the inward fluid speed remains positive at the wall v + > 0. According to (2.17a), this condition puts an upper bound on the strength parameter for deflagration solutions (5.4) The compression wave must revert to a Friedmann solution within the boundary imposed by causality in order to ensure that patches of the Universe causally disconnected with the bubble remain unperturbed.Causal connection is established in the region v < v lc where v lc is the value of the fluid speed for constant-ξ observes v at the light cone.The light cone is found by looking at the region where the vector (3.16) has zero norm, i.e. g µν V µ ξ V ν ξ = 0, which occurs at v lc = 1. (5.5) The compression wave matches on to a FLRW Universe through a shock, i.e. with a new discontinuity that must satisfy the junction conditions (4.8) while matching to the FLRW solutions (4.18) outside the shock, which satisfy (5.6)

Detonation
In a detonation, the fluid speed v − just behind the wall is always larger than the sound speed of the broken phase v − > c s− , and the fluid moves as in a FLRW Universe for constant-ξ observers outside the wall, u = 0.This provides a natural initial condition for the inward integration.As in the Minkowski case, the positive sign in (2.17b) must be chosen in order to ensure that the fluid moves with a supersonic velocity v − > c s− behind the wall.The condition v + > 0 in (2.17a) is always satisfied with the choice of the positive sign, so that no constraint on α + comes from the junction conditions.However, according to Eq. (3.12), α + > 1/3 means that the sound speed in the symmetric phase has negative square: c s 2 + < 0. Although these solutions seem to be allowed by the equations of motion (5.1) and the junction conditions (2.17a), a negative sound speed squared would indicate an instability, and we therefore look for solution within the bounds (5.4).
Equipped with the initial conditions (U − , Ω − , Φ − ), we can integrate inwards towards the line of fixed points (U * (k), Ω * (k), Φ * (k)), obtaining the profile of a rarefaction wave behind the wall.The asymptote of the solution at the line determines the value of the parameter k, and how to match to a trajectory smoothly reaching the origin.
The condition v − > c s− implies that there is a minimum fluid velocity outside the wall given by v + (α + , c s− ); this is the Chapman-Jouguet speed v CJ .As explained in Section 2.2 and after equation (3.9), in our model the sound speed in the broken phase is fixed at c s− = 1/ √ 3, for which, from equation (2.17a), one finds Hence, in a detonation the fluid velocity in the constant ξ frame just outside the wall v + is bounded from above and below by v CJ ≤ v + < 1. (5.8)

Supersonic deflagrations
Supersonic deflagration (or hybrid) solutions are found between the two aforementioned cases when the fluid exit speed in the constant ξ frame is v − = c s− , while the fluid is at rest with the FLRW cosmic fluid for constant-ξ observers, i.e. u = 0, outside the bubble.The outer FLRW Universe provides a natural initial condition for the inward integration from the shock.The fluid speed ahead of the shock for constant-ξ observers is bounded by (5.9) In the region v sh + > v CJ the solution is a detonation, while in the region v sh + < c s− the solution is a deflagration.The junction conditions are the same for subsonic and supersonic deflagrations: at the shock and we integrate inward along the compression wave until the wall, where the junction condition (2.17b) has to be taken with the negative sign in order to ensure that U + > 0. As in a subsonic deflagration, the bounds (5.4) apply.The junction conditions at the wall provide the initial conditions (U − , Ω − , Φ − ) for a second inward integration along the rarefaction wave towards the line of fixed points γ ⋆ (k).The intersection between the trajectory of the solution and the fixed line determines the value of the parameter k and the matching to a trajectory that smoothly reaches the origin, as done for the detonation case.Integrating outward from the shock, the we stop the solutions at the light cone where v = 1.

Numerical methods
In our model, the EOS inside the bubble is fixed with ω − = 1/3.We wish to find solutions as a function of the EOS parameter in the symmetric phase ω + , which determines the strength of the transition at the wall (3.12), and the wall position ξ w that discriminates between the three possible cases: subsonic deflagration, detonation, and supersonic deflagration.In order to find the solutions, we must determine the parameters k and a F , which characterize the asymptotic solutions at the origin and at the light cone respectively, and ξ sh , which determines the position of the shock in a deflagration.In this section we provide some numerical methods that we used to find our solutions and comment on the results.In all the cases of study of this paper, we performed the integration of the system (3.14) using the fourth-order Runge-Kutta scheme encoded in the Python SciPy solver integrate.solve_ivp.We adopted different shooting algorithms to obtain the optimal values of the parameters k, a F and ξ sh using the bisection routine carried out by the SciPy method optimize.bisect.
In all the discussed cases, the initial displacement in ξ from the origin that we used for numerical integration is 10 −3 .To maximize the precision of the solutions we use the minimum relative (rtol = 2.3e-14) and absolute (atol =1E-18) tolerance levels allowed by the solver integrate.solve_ivp.We will explain in the following paragraphs that the position of the shock (wall) in a subsonic (supersonic) deflagration as well as the position of the fixed line in detonation and hybrid solutions are dynamically found when the solver integrate.solve_ivpmeets certain requirements.The precision with which we can fix the interfaces and the fixed line is controlled by the sampling spacing ∆ξ in the ξ coordinate.A smaller sampling spacing provides a better resolution, at the price of a lower efficiency of the code.Depending on the case, we aim to find the best compromise, keeping as an upper bound ∆ξ < 10 −4 .

Subsonic deflagration
Integrating outward from the origin with fixed ω + and ξ w , we shoot with the parameter k to determine the shock position ξ sh .The parameter a F can be read afterwards from the junction conditions at the shock, but it does not play any role in the shooting procedure for deflagrations.Keeping in mind that the solution must match on to a FLRW Universe (4.18) outside the shock, the shooting algorithm for optimal values of k and ξ sh proceeds as follow: 1. for a given value of k in the range (4.14) we integrate forward from the origin up to the wall.
2. The junction conditions (4.8) at the wall give the boundary conditions for the trajectory of the solutions ⃗ Y just outside the wall.We integrate forward until the first shock condition (4.6d) [T n n ] ± sh = 0 is met, where and we used the fact that the junction conditions (2.15) imply Φ + sh = Φ − sh and outside the shock the relations (4.18) hold.

At this same point we evaluate the jump [T n
0 ] ± sh (4.6c), which turns out to be approximately a linear function of k with a unique zero.
4. We use a bisection algorithm to find the unique value of k for which [T n 0 ] ± sh = 0 at the same value of ξ where [T n n ] ± sh = 0, knowing that k is bracketed by (4.14).The value of ξ at which this condition is realized defines the shock position ξ sh .

The trajectory of solutions ⃗
Y + sh just outside the shock given by the solution of the junction conditions (2.15) determines the value of a F through (4.20).
All the unknown parameters, k, ξ sh and a F , are thus set.

Detonation
In Section 5.1 we have shown that all the solutions starting from the origin with initial conditions (4.12) end on the line of fixed points γ ⋆ (k), and k is a continuous parameter of the curve.The result of this analysis is shown in Figure 3.The strategy is to split the shooting algorithm in two sections.First we integrate inward from the wall aiming to reach the fixed line γ ⋆ (k), and then integrate outward from the origin aiming for the same point on the fixed line.As there is no shock in a detonation, the only two parameters to set are k and a F .In details, the algorithm works as follow: 1.For a fixed value of ξ w we select an arbitrary value of a F in the range (5.8) that defines the trajectory of solutions ⃗ Y outside the wall .The junction conditions (2.15) give the boundary conditions for the trajectory ⃗ Y just behind the wall.
2. We integrate inward from the wall until the condition v = c s− is met.The end point of ⃗ Y is now on a fixed point5 .
3. By direct inspection, we find that the value of ξ at the end point of the trajectory of the solution, that we will denote as ξ end , is approximately a linear function of a F , while the values of the other three variables ⃗ Y (ξ end ) at the end point are approximately independent of a F . 4. There is a unique value of a F such that ξ end = ξ ⋆ , that is such that (ξ end , ⃗ Y (ξ end )) reaches the fixed point line γ ⋆ (k).We find the optimal parameter a F with a bisection method, using the fact that a F is bracketed by (5.8).
5. The intersection between the trajectory of the solution and the fixed point line identifies a point (ξ end ⋆ , ⃗ Y end ⋆ ) which allows us to extract the parameter k by inverting the parametrization of the curve The unknown parameters k and a F are thus determined.

Supersonic deflagrations
It was mentioned in Section 5.2 that, differently from the subsonic deflagration case, the asymptotic FLRW solution outside the bubble provides a natural choice of initial conditions for the inward integration from the shock.The most convenient choice of input parameters is then ω + and shock position ξ sh , while through the shooting algorithm we aim to determine the metric coefficient a F , the wall position ξ w and the critical density parameter at the origin k.As in the detonation case, we use the analysis of fixed points in Figure 3 to split the shooting algorithm in two different sections, the first for a F and ξ w , and the second for k.In the following we explain the algorithm for supersonic deflagrations step by step: 1.With fixed ω + and ξ sh , we choose a value of a F that respects the bounds (5.9).This parameter uniquely defines the trajectory of solutions ⃗ Y outside the shock.
2. The junction conditions (5.10) at the shock provide the initial condition for the inward integration.
3. Integrate inward until the wall position ξ w defined by v − = c s− .According to (2.17b), this condition is realized when At each step of the integration we verify whether the above condition is met.The point in ξ where this happen defines the wall position, and there we stop the integration.
4. The junction conditions (4.8) at the wall ensure that v − = c s− and provide the initial condition for a second inward integration toward the fixed line.However, we discussed in Section 5.1 that v = c s− is one of the two conditions that define a fixed point.To leave the fixed point we introduce a small displacement v − = c s− + δ with δ = 10 −7 .This provides safe initial conditions for the inward integration.
5. Integrate inward from the wall toward the fixed line until the condition v = c s− is met.As in the detonation case, we verify that the end point of the trajectory of the solution ⃗ Y is a fixed point.From here, the algorithm follow the same strategy of last three bullets drawn for detonation solutions.
After this shooting routine, all the necessary parameters are fixed.

Results and discussion
We present in this section our numerical solutions for expanding bubbles and comment on the results.The plots in Figures 4 to 9 are organized as follow: on the left we display quantities that are directly obtained from the integration of the equations of motions (3.14), U , the radial component of the fluid four-velocity for an Eulerian observer, Ω, the energy on a shell of radius R, Φ, the gravitational potential at R, Γ, the generalized Lorentz factor, and a, the time dilation with respect to the origin.On the right we display quantities that are derived from the previous ones.On top we put quantities related to the fluid: the fluid speed for constant-ξ observers v and the relative fluid speed u as measured by constant-ξ observers V ξ and V F ξ , the ratio Ω/Ω F between the energy on a shell of radius R and its value in a FLRW Universe with EOS parameter ω + , and the shear scalar σ normalized to the local expansion rate (θ/3).At the bottom we put quantities related to the curvature of the spacetime: the Ricci curvature scalar R (4) and the Kretchmann scalar K (4) normalized to the values that would have in the outer FLRW Universe, R (4) and K (4) respectively.We studied three different expansion modes, subsonic deflagration, detonation, and supersonic deflagration.
For each mode we perform two different analysis corresponding to two different figures for each case; in one we fix the EOS inside and outside the bubble and plot a set of solutions with different wall position ξ w , while in the other we fix the wall position (shock position for hybrid solutions) and plot a set of solutions for different values of the EOS parameter ω + .

Subsonic deflagrations
We show numerical results for deflagration solutions with fixed strength α + and different values of ξ w in Figure 4, while deflagrations with fixed ξ w and for different values of the strength parameter at the wall α + in Figure 5.
The fluid speed u contains information about the fluid velocity as seen by constant-ξ observers.Constant-ξ observers outside the bubble signal to the constant-ξ observers inside the bubble about the Universe they live in; inner constant-ξ observers then measure a nonzero relative velocity u with respect to a hypothetical FLRW constant-ξ observer at the same position, which defines a rest frame for the cosmic FLRW fluid.Inspecting Figure 4 (top right) we see that inner constant-ξ observers would measure a rarefaction wave behind the wall with respect to the cosmic FLRW fluid.Ahead of the wall, constant-ξ observers see a compression wave that ends at the shock.Outside the shock, constant-ξ observers measure zero relative velocity u, by definition, and they agree that they live in a flat FLRW Universe.In deflagration solutions in Minkowski spacetime (Figure 1), the fluid is completely at rest with the Universe frame behind the wall.
Information about curvature can be extracted by looking at constant-r observers; these observers are comoving with the fluid and see an anisotropic expansion rate inside the bubble, as shown by the departure of the shear scalar σ from zero in the interior of our solutions.This means that the Universe inside the bubble, is not isotropic, except the origin, and therefore F and K F .
cannot be reduced to any FLRW Universe.The sharp peaks in σ/(θ/3) at the shock come from the fact that the divisor, the local average expansion rate, becomes very small.We verify with the Kretschmann scalar K (4) that this does not lead to a curvature singularity.The energy ratio profile Ω/Ω F that constant-r observers measure can be compared to the enthalpy density profile of the Minkowski spacetime solutions; in both cases, constant-r observers measure a deficit of energy respect to a FLRW Universe with EOS parameter ω + inside the wall, and an excess of energy in the compression wave.A more quantitative comparison will be performed later in Figure 10.

Detonations
In Figure 6 we display detonation solutions for different values of the wall position ξ w , with fixed α + , and in Figure 7 we show detonation solutions for different values of ω + (or α + equivalently) with fixed ξ w .The insets in the first four plots of Figure 6 zoom in to the region where the solutions cross the fixed line, showing the small differences in the solutions from the centre to the fixed point.The solutions and their derivatives are continuous there.Detonation solutions meet the FLRW Universe just outside the wall, so that constant-ξ observers do not measure a compression wave ahead of the wall.Inside the wall, u differs from zero, and inner constantξ observers measure a departure from the FLRW universe signalled in by outer constantξ observers.While in the Minkowski spacetime solutions rarefaction waves in detonations extend inward only to the fixed point at the sound speed, in the gravitating solutions the rarefaction wave behind the wall extends to the centre of the bubble.We notice further that the profile of u always intersects the line u = 0 between the fixed point an the wall.This means that there exists a constant-ξ observer inside the bubble that sees the fluid moving inward with the same speed that an FLRW constant-ξ observer would measure.Although this observer measures a vanishing local departure from FLRW in the fluid speed, the observer can still realize that the local Universe is not FLRW by looking at the Kretschmann scalar K (4) .The fact that K (4) differs from the value that it has outside the wall means that the spacetime inside the wall is geometrically different from the spacetime outside the wall.
The statement that the the inner and outer sides of the bubble have different geometries is also supported by the fact that constant-r observers measure an anisotropic expansion σ inside the wall.As mentioned in the deflagration case, this points to the fact that the inner side of the bubble is not a FLRW Universe except in the origin.We further notice the the position of the fixed point in the four dimensional phase space (ξ, U, Ω, Φ) mostly depends of the EOS parameter ω + , while it is almost independent on ξ w .

Supersonic deflagrations
In Figure 8 we show some examples of hybrid solutions for different values of the shock position with fixed transition strength parameter α + , while in Figure 9 we display some example solutions with fixed shock position but different equation of state parameter in the symmetric phase ω + .In supersonic deflagrations the discontinuities at the wall and the shock are the strongest.Constant-ξ observers measure a compression wave ahead of the wall as in the deflagration case, and a rarefaction wave behind the wall as in detonation case; the rarefaction wave passes through the fixed line and extends to the origin.

Newtonian estimation of the curvature of spatial sections near the origin
In Section 4.2 we introduced the curvature of spatial sections near the origin (4.13) and its relation with the parameter k.In this section we want to provide a quasi-Newtonian calculation of the amount of spatial Ricci curvature that we would expect from the Minkowski spacetime solutions.
We consider linear scalar perturbations around a flat FLRW Universe [71] in the conformal Newtonian gauge [72]  where η is the conformal time, and the Bardeen potentials Φ B and Ψ B are two gauge invariant scalar perturbations.The curvature of spatial sections can be then computed at fist order in scalar perturbations as On the other hand, neglecting for the sake of simplicity the perturbations in the radial momentum density T 0i , the linearised Einstein equations [71,72] reduce to the Poisson equation for the Bardeen potential where in the last step we used the Friedman equation H 2 = 8πe/3, δe is the perturbation in the energy density of the fluid and δ = δe/e.Combining the above expressions, we give the quasi-Newtonian approximation of the spatial curvature induced by the expanding bubble as where δ is the energy density contrast in the Minkowski solutions, which we obtain by modifying the code snipped provided in Ref. [59].Given the barotropic EOS p = ωe, we can write δ = w/w 1 − 1, with w 1 the enthalpy outside the shock.
In comparing the gravitating and Minkowski solutions, we run into the problem of comparing the coordinate ratios ξ: the gravitating solutions extend to ξ > 1, which would be outside the light cone in Minkowski space.We choose instead to compare by the ratio of xi to its value on the light cone, ξ lc .In the Minkowski spacetime solutions ξ lc = 1 always, while in the gravitating solutions the light cone position is defined by the condition v lc = 1 (5.5).
In Figure 10 we plot the Ricci curvature of spatial sections near the origin R (3) (ξ → 0) against ξ w /ξ lc and compare it with the quasi-Newtonian estimate that we obtain from the solutions in Minkowski spacetime.
The two cases show an overall similar dependence on the wall position ξ/ξ lc , with a maximum in the hybrid regime.The amplitudes of the spatial Ricci curvature at the origin are around the same order of magnitude for the gravitating and quasi-Newtonian approximation cases, but gravitating deflagrations have relatively less negative curvature, and detonations more.
Gravitating solutions are more subject to numerical errors, because the system of Einstein and continuity equations (3.14) is highly non-linear.In particular supersonic deflagration seems to be the mode that mostly suffers from numerical error, as the compression wave can become very narrow and not well resolved (see the central panels of Figure 1).Hybrid solutions indeed have an additional source of error that come from the fact that, as explained in Section 6.3, the inward integration from the wall starts at v − = 1/3, and a small displacement is needed to leave the fixed point.4 and the text for a description of the quantities plotted.

Energy redistribution
As the transition proceeds and the bubble expands, part of the energy difference between the symmetric and broken phase is converted into kinetic energy of the fluid around the bubble [73].The fluid motion turns on and enhances shear stress perturbations, and becomes a source for gravitational radiation [17].The energy budget of the phase transition is crucial to quantify the power spectrum of gravitational waves [48,54,55].The total amount of energy that was contained in the volume occupied by the bubble before the transition happened is where e F is the homogeneous and isotropic energy density of the outer FLRW Universe.During the transition, the energy density of the fluid as measured by constant-ξ observers U µ (4.24), that is observers that move with speed u relative to the fluid, is T µν U µ U ν .We extract the kinetic part by subtracting the energy density in the rest frame of the fluid from the total energy density [36,51] e K ≡ T µν U µ U ν − T µν u µ u ν , (7.6) and we get e K = wu 2 γ 2 u .(7.7)The kinetic energy E K of the fluid is then calculated by integrating (7.7) over the volume of the bubble and the ratio measures the amount of energy contained in the bubble that has been converted into kinetic energy during the transition.We show the kinetic energy fraction obtained from the gravitating solutions in the top panel of Figure 11.We compare these results with the kinetic energy fraction that we obtain from a single bubble with the same EOS (3.9) (3.10) expanding on a Minkowski flat background (central panel).For the latter case, we made use of the publicly available code detonation sector shows a decrease of K with the wall position ξ w /ξ lc in both cases, but the curves are much steeper in the gravitating solutions.

Conclusions
We presented new self-similar solutions of the general relativistic hydrodynamic equations that describe the expansion of bubbles of the broken phase into the symmetric phase, in the limit that the ratio of the proper wall width and the bubble size go to zero.
To find these solutions, we assumed a simple equation of state compatible with selfsimilarity.The transition must be strongly supercooled, as the bubbles need a long time t ≫ t n = t(T n ) to reach sizes comparable to the Hubble radius.In the simplified equation of state model, the energy density of the metastable phse (2.8) decreases as T −4 , which is different from the standard scenario of supercooled phase transitions [74][75][76][77][78][79][80], where the energy density of the false vacuum is asymptotically constant in time.Nonetheless, we expect that the solutions given some guidance on the gravitational effects to be expected in the standard scenario, where large density perturbations δe/e ∼ O(1) are also produced.
Relaxing the self-similar assumption requires solving the equations of motion as a system of partial differential equations, and might be a challenging topic for further studies.The self-simular solutions will be a useful benchmark for such simulations.
Our solutions show qualitative and quantitative differences respect to the analogous Minkowski solutions.Constant-ξ observers see a local departure from a FLRW fluid u as a rarefaction wave behind the wall that extends to the origin.In the Minkowski solutions this rarefaction wave extends inward only as far as the sound speed in the detonation and supersonic deflagration case, and it is not present at all in deflagrations.As expected from the Minkowski solutions, the supersonic deflagration is the mode with strongest discontinuities in the fluid speed v and energy Ω at the wall.
Gravitational effects modify the distribution of the kinetic energy of the fluid around the bubble.We found that the kinetic energy fraction K (7.9) is enhanced by gravitational effects in the regime of supersonic deflagrations and fast subsonics deflagration, while detonations and slow deflagrations on a curved spacetime are less efficient than the same solutions on a Minkowski flat spacetime.This will have direct consequences in the estimation of the power spectrum of gravitational waves [48,54,55], and might significantly change its amplitude when the relative difference between the two cases is of order O(1).
In all cases, the bubble shows a deficit in the energy density inside the wall respect to the outside.This fact implies that the spatial sections inside the wall are negatively curved.We quantified the general relativistic effect on the energy deficit inside the wall by looking at the amount of spatial curvature at the centre of the bubble and comparing with a quasi-Newtonian estimate derived by the solutions in Minkowski spacetime.The amplitude of spatial Ricci curvature has a peak in the hybrid regime and is, but the latter case shows larger negative curvature in the subsonic deflagration cases, while the former has larger negative curvature in the detonation regime.
Finally we mention prospects for gravitational wave observations.A simple analysis of the Weyl tensor components shows that spherical symmetry protects the single bubble configuration from sourcing gravitational waves.Primary gravitational waves can instead arise when the spherical symmetry is broken, i.e. by collision of bubbles and sound waves, whose initial conditions are set by the velocity and energy density profiles around the bubble [54,55].
Gravitational effects change these profiles, which can be expected to lead to modifications in the gravitational wave power spectrum.
We also expect bubbles to produce metric perturbations which source secondary gravitational waves.The source of primary gravitational wave is approximately δT ij ∼ H 2 (δe/e), while the source of scalar induced gravitational waves can be estimated as ⋆ (HR ⋆ ) 4 (δe/e)

with h ab ≡ g µν e µ a e ν b ( 4 . 2 )
the induced metric on the hypersurface (the first fundamental form), andK ab ≡ e µ a e ν b ∇ µ n ν (4.3)the extrinsic curvature (the second fundamental form).Moreover, we have defined the surface stress-energy tensor S ab as the singular part of the energy momentum tensor at the interfaceT µν (y) = wu µ u ν + pg µν + S µν δ(C(y)).(4.4)Combining (4.1b) with the Gauss-Codazzi equations one finds[70]

Figure 4 :
Figure4: Deflagration solutions against self-similar variable ξ for fixed EOS parameter in the broken (ω − = 1/3) and symmetric phase (ω + = 0.2) and different values of the wall positions ξ w .On the left quantities used in the numerical integration U, Ω, Φ, Γ, a; on the right fluid speed v, u, energy Ω normalized to its value in FLRW Ω F and shear σ normalized to the local expansion θ/3 and curvature related quantities R(4) and K(4) normalized to their respective values in the outer flat FLRW Universe R

Figure 5 :
Figure 5: Deflagration solutions against self-similar variable ξ for fixed wall position ξ w and varying equation of state parameter of the symmetric phase ω + .See the caption to Fig. 4 and the text for a description of the quantities plotted.

Figure 6 :
Figure 6: Detonation solutions against self-similar variable ξ for fixed EOS parameter in the broken (ω − = 1/3) and symmetric phase (ω + = 0.2) and different values of the wall positions ξ w .See the caption to Fig. 4 and the text for a description of the quantities plotted.

Figure 7 :
Figure 7: Detonation solutions against self-similar variable ξ for fixed wall position ξ w and varying equation of state parameter of the symmetric phase ω + .See the caption to Fig. 4 and the text for a description of the quantities plotted.

Figure 8 :
Figure 8: Supersonic deflagration solutions against self-similar variable ξ for fixed EOS parameter in the broken (ω − = 1/3) and symmetric phase (ω + = 0.2) and different values of the wall positions ξ w .See the caption to Fig. 4 and the text for a description of the quantities plotted.

Figure 9 :
Figure9: Supersonic deflagration solutions against self-similar variable ξ for fixed shock position ξ sh and varying equation of state parameter of the symmetric phase ω + .See the caption to Fig.4and the text for a description of the quantities plotted.

Figure 10 :
Figure 10: Ricci curvature of spatial sections evaluated at the origin, normalized by the local expansion scalar θ/3, in the gravitating solutions (upper panel) and quasi-Newtionian estimate inferred by the Minkowski spacetime solutions (central panel), against the fractional distance of the wall from the light cone ξ w /ξ lc ; in the bottom panel relative difference between the two cases above.The light cone in the Minkowski solutions is always ξ lc = 1.In a FLRW background the local expansion θ/3 coincides with the global Hubble constant H.Some of the gravitating solutions in the hybrid regime belong to the deflagration regime in the Minkowski case for the same value of ξ w /ξ lc .For these solutions we use a grey colormap in the bottom plot.
[81]8.1)with ΦB the scalar Bardeen potential.Telative contribution of scalar induced gravitational waves∂ i Φ B ∂ j Φ B δT ij ∼ (HR ⋆ ) 2 δe e (8.2)becomes relevant when R ⋆ H ∼ O(1).A better estimate of the secondary gravitational waves generated by the expanding bubble in slow FOPTs can be achieved by performing a complete analysis of scalar cosmological perturbations[81].