Gravitational wave spectra from strongly supercooled phase transitions

We study gravitational wave (GW) production in strongly supercooled cosmological phase transitions, taking particular care of models featuring a complex scalar field with a U$(1)$ symmetric potential. We perform lattice simulations of two-bubble collisions to properly model the scalar field gradients, and compute the GW spectrum sourced by them using the thin-wall approximation in many-bubble simulations. We find that in the U$(1)$ symmetric case the low-frequency spectrum is $\propto\omega$ whereas for a real scalar field it is $\propto\omega^3$. In both cases the spectrum decays as $\omega^{-2}$ at high frequencies.


I. INTRODUCTION
The direct detection of gravitational waves (GWs) from a binary black hole merger by LIGO [1] marked the dawn of a new era in astrophysics and cosmology. In the next decades various experiments will probe GWs in a wide range of frequencies [2][3][4][5][6][7][8][9]. In addition to the astrophysical GW sources, such as compact object binaries, these experiments will probe also cosmological GW backgrounds providing a unique probe of the early Universe, as, unlike electromagnetic signals, GWs can propagate freely from the very beginning of the Universe.
Cosmological first-order phase transitions constitute one possible source of GWs from the early Universe [10], which can be probed by the upcoming GW experiments [11,12]. In a first-order phase transition the false vacuum is separated from the true vacuum by a barrier as the transition proceeds. As a result the unstable vacuum decays through nucleation of bubbles, corresponding to the field trapped in the false vacuum tunnelling through the barrier [13][14][15]. After nucleation these bubbles grow until they collide, eventually converting the whole Hubble volume into the new phase.
The bubbles typically grow by many orders of magnitude between their nucleation and collisions, releasing a lot of energy. This energy goes into gradient and kinetic energy of the bubble walls and into motion of the plasma as the particles in the plasma interact with the bubble wall. GWs from a phase transition are sourced by the scalar field gradients [16] and motions in the plasma [17]. For very strongly supercooled phase transitions the plasma friction can be negligible. In this case the bubble walls reach velocities near the speed of light before they collide [18,19] and the GW signal is dominated by the scalar field gradients [20].
The GW signal from scalar field gradients was first calculated in the envelope approximation in Ref. [16]. In this approximation the bubble walls are treated as thin shells that disappear in the collisions. The resulting GW spectrum is a broken power-law that at low frequencies grows as ω 3 and at high frequencies decays as ω −1 [21][22][23]. In Refs. [23,24] the envelope approximation was extended in order to model colliding fluid shells. In this so called bulk flow model the bubble wall energy was assumed to decay as R −2 after the collision as a function of the bubble radius R. The bulk flow approximation results in a GW spectrum that turns from ω 1 behaviour to ω −2 at around the same frequency at which the spectrum in the envelope approximation peaks.
Recently the GW spectrum was calculated in 3D lattice simulations [25][26][27]. These simulations are very difficult because of the large separation between the characteristic length scales in the problem, that is the size of the growing bubble and its thinning wall. Due to numerical limitations in such simulations it is impossible account for realistically large bubble wall velocities. However, the GW spectrum can still be computed. At high frequencies it was found that the spectrum lies somewhere between the envelope and bulk flow approximations. The low-frequency behaviour of the spectrum is especially difficult to resolve and a ω 3 behaviour is typically assumed.
In this paper we approximate the GW spectrum from strongly supercooled phase transitions by first studying the scaling of the gradient energy in two-bubble collisions by lattice simulations and then calculating the GW spectrum by performing many-bubble simulations in thinwall approximation. In this way we can efficiently perform large simulations with realistic behaviour of the GW source. In contrast to the earlier works, we study also the case of a complex scalar field with a U(1) symmetric potential which is often realized in particle physics models 1 . We find that in this case the gradient energy quickly reaches an R −2 scaling after the collision, whereas for a real scalar field we find that the decay is much faster. In the former case we find that the GW spectrum is near the bulk flow result, and in the latter case we find a GW spectrum that grows as ω 3 at low frequencies and decays as ω −2 at high frequencies.

II. BUBBLE COLLISIONS
We begin by studying collisions of two complex scalar field bubbles. In order for the scalar field gradients to be the dominant source of GWs the phase transition has to be severely supercooled [20]. Typically this can not be realized in models based on polynomial potentials [36], but in models that are classically scale invariant a prolonged period of supercooling is possible [29,30,33,[37][38][39][40][41][42][43][44][45][46][47][48][49][50]. In these models the symmetry breaking originates from radiative corrections [51] and finite temperature effects give raise to a potential energy barrier between the symmetric and the symmetry-breaking minima. The one-loop effective potential is of the form where B and C are dimensionless constants that depend on the couplings of the scalar field φ (see e.g. [42]), v is the vacuum expectation value of |φ| and T denotes the temperature of the plasma. Motivated by this, we consider logarithmic potential where κ is a dimensionless parameter. This potential is U(1) symmetric and its global minimum The parameters B and C of Eq. (1) are related to κ and ∆V via The radial initial profile of the modulus |φ| for an O(4) symmetric bubble is obtained as the solution of with boundary conditions ∂ r φ = 0 at r = 0 and φ → 0 at r → ∞. Due to the U (1) symmetry of the potential every bubble will be nucleated with a complex phase ϕ of the field φ chosen from the range ϕ ∈ [0, 2π[ with equal probability. We assume that the phase transition finishes within a Hubble time and therefore neglect the background expansion. Collision of two initially O(4) symmetric scalar field bubbles is O(1, 2) symmetric, and it is convenient to define new coordinates [52] where r 2 = x 2 + y 2 . The bubbles lie at the z axis. The Klein-Gordon equation for the real (X = R) and imaginary (X = I) parts of the field in these coordinates simplify to where + and − signs correspond to the regions t ≥ r and t < r, respectively. The collision of bubbles occurs in the region t ≥ r where we solve the above equation numerically. In the region t < r the evolution is given by analytical continuation of the initial bubble solution where z j denotes the position of the bubble j.
In Fig. 1 we show the result from two-bubble collision in three values of the phase difference between the colliding bubbles: ∆ϕ = 0, ∆ϕ = π/2 and ∆ϕ ≈ π 2 . The numerical lattice calculation is performed in dimensionless variables, obtained by scaling φ → φ/v and x µ → √ ∆V x µ /v, and the results in Fig. 1 are shown in the simulation units. We see that the bubble walls quickly accelerate after nucleation, approaching velocities near the speed of light before their collision. 3 After the collision we see that a sharp phase wall continues to propagate with a constant velocity near the speed of light. This can also be seen in the lower panels, where we show the evolution of the gradient energy density of the scalar field, ρ grad = |∂ z φ| 2 /2. It is also clear from these plots that the energy loss of the gradients after the collision is much faster in the case where the bubbles have equal complex phases effectively corresponding to the case with a real scalar field.
In Fig. 2 we show by the blue solid line the time evolution of the maximum of gradient energy density averaged over various values of the phase difference ∆ϕ. We obtain the case of a real scalar field from our analysis by taking only the case ∆ϕ = 0. This is shown by the red solid line. The collision happens at t = t c . Before the collision the total released energy scales as E rel ∝ t 3 , the surface area as A ∝ t 2 and the wall thickness as L ∝ 1/t due to increase in the Lorentz factor of the wall. Therefore the energy density at the wall scales as ρ grad ∝ E rel /(AL) ∝ t 2 , which is what we see also in Fig. 2. After the collision, if the bubble wall velocity is constant and the energy remains localized at the bubble wall, its gradient energy density scales as ρ grad ∝ t −2 as E rel and L remain constant. From Fig. 2 we see that this is not a perfect description. In reality some of the energy is spread into the collided volume and therefore the maximal gradient energy density decays faster. In the case of a real scalar field we see that the scaling of gradient energy density reaches ∝ t −3 behaviour after the collision, while in the U(1) symmetric case it scales significantly slower, ∝ t −2 . This is a consequence of the phase difference between the colliding bubbles that after the collision for a long time continues to propagate as a sharp phase domain wall.  The results shown in Figs. 1 and 2 are from simulations with κ = 0.2, but we have checked that the scaling of the gradient energy density after the collision is not sensitive to the value of κ. We have also checked that the same scaling results hold in the case of a polynomial potential 4 . In addition, we have run simulations with different values for the initial bubble separation d, that in Figs. 1 and 2 was set to d = 20 (in the simulation units). Our results from simulations with d ∈ [8,200] indicate that the bubble separation dictating that the γ factor of the wall at collision and its energy does not change the scaling behaviour.

III. PRODUCTION OF GRAVITATIONAL WAVES
Next we will study the production of GWs from scalar field gradients using the scaling results of the bubble wall energy obtained in the previous section. Following Ref. [55], the total energy spectrum in a directionk at frequency ω of the GWs emitted in the phase transition is given by where Λ ijlm is the projection tensor, 4 In the class of potentials we focus on false vacuum trapping in the collision typically does not occur. If it did that would change the result, as most of the gradient energy would remain around the collision point [53,54]. and T ij is the traceless part of the stress energy tensor, In the thin-wall limit, the gradient energy carried by an uncollided element of the bubble wall at solid angle dΩ x can be approximated as [16] where x n denotes the position vector of the bubble center,x is a unit vector that points from the centre of the bubble in the direction dΩ and R n ≈ t − t n is the radius of the bubble that nucleated at time t n . After the element of the bubble wall at dΩ x has collided with another bubble, its energy starts to decrease. Assuming that the velocity of the wall element does not change in the collision, the scaling of the energy can be accounted by multiplying Eq. (11) by a function f (t) which depends on the time t n,c = t n,c (x) when the wall element collides with another bubble. Before the collision f (t < t n,c ) = 1, and the envelope approximation corresponds to taking f (t > t n,c ) = 0. Assuming instead that the bubble wall loses energy ∝ R −2 we get the bulk flow approximation [23,24] where On the basis of the results shown in Fig. 2, we find that the decay of the maximum of the gradient energy density can be approximated as a broken power-law, that changes from an ∝ R −ξ1 behaviour to ∝ R −ξ2 . In the thin wall approximation, the bubble wall energy is simply dE wall ∝ dΩ x R 2 L max[ρ grad ], where the bubble wall width L after the collision is constant, and therefore For the U(1) symmetric case we use b 1 = 0.6, b 2 = 0.4, ξ 1 = 8, ξ 2 = 2, and in the case of a real scalar field we use b 1 = 0.93, b 2 = 0.07, ξ 1 = 8, ξ 2 = 3. These approximations are shown in Fig. 2 with the yellow dotdashed and dashed lines, respectively. The contribution from N bubbles on the traceless part of the stress energy tensor (10) can now be written as We can rotate the coordinate system such that a given k =k(φ k , θ k ) after the rotation points to z direction, k →k = (0, 0, 1). Then, the projection in Eq. (8) simplifies to [21] dE where with g + (φ x ) = cos(2φ x ) and g × (φ x ) = sin(2φ x ). The spatial angles in the rotated coordinate system are We consider the exponential bubble nucleation rate per unit volume, Γ ∝ e βt . The parameter β, with mass dimension 1, sets the time and length scale of the transition. The abundance of GWs produced in bubble collisions in a logarithmic frequency interval is then (17) where α = ∆V /(ρ tot − ∆V ) characterizes the strength of the transition, H 2 = 8πGρ tot /3 is the Hubble rate, and gives the spectral shape of the GW background. The volume over which Ω GW is averaged is denoted by V s . We note that dΩ k |C + | 2 +|C × | 2 ∝ V s /β 5 . Next we calculate the S(ω) function numerically.

IV. GRAVITATIONAL WAVE SPECTRUM
In order to determine the spectral shape of the GW signal we simulate the phase transition by nucleating bubbles according to the rate Γ ∝ e βt inside a cubic simulation volume with periodic boundary conditions. We neglect the initial bubble sizes, assume that their wall is infinitesimally thin and evolve the bubble radii as R n = t − t n . We generate points on the surface of a bubble and find the time t n,c when each of these points collides with another bubble surface by the bisection method [56]. As an example, in Fig. 3 we show the surface of the first bubble that nucleated with the color coding indicating the collision time t n,c . Once the collision times are known, we can simply integrate the functions C +,× , and finally compute the GW spectrum (18).
The time integral in Eq. (15) can be evaluated analytically if f (t, t n,c ) is a (broken) power-law. We perform the remaining integrals overk andx directions numerically. As the simulation volume is not spherically symmetric we calculate the spectrum only for 6k directions that correspond to the normal vectors of the cube faces. In order to accurately determine the GW spectrum, we calculate the average spectrum over multiple simulations. 5 In Fig. 4 the red and blue solid curves show, respectively, our result for the spectral shape of the GW spectrum in the U(1) symmetric case, and the case where all bubbles have equal complex phase, ∆ϕ = 0, which corresponds to having a real scalar field. The green and yellow dashed curves, for comparison, the result in the envelope and bulk flow approximations. The results are obtained by taking geometric mean over 40 simulations with around 200 bubbles each, while the error bands show the corresponding geometric standard deviations. 5 Averaging many smaller simulations is equivalent to running one much larger simulation but much more easily parallelised. We calculate a broken power-law fit to the spectrum parametrized as where A andω correspond to the peak amplitude and frequency, c determines the width of the peak, and a, b > 0 are the low-and high-frequency slopes near the peak of the spectrum. The first bracket parametrizes the change in the low-frequency slope which in the U(1) symmetric case resembles the envelope result shortly after collision with emission from a slowly decaying gradient dominating later on. In the other cases we take d = a for which the first bracket in (19) gives 1. 6 The parameter values of the fit are shown in Table I.
Finally, the present day GW spectrum can be obtained by simple red-shifting as [17] is the inverse Hubble time at the transition redshifted to today 7 . The transition temperature is denoted by T * and the effective number of relativistic degrees of freedom at the temperature T * by g * . At scales larger than the horizon scale at the time of the transition the spectrum scales as ω 3 because the source is diluted by the Hubble expansion [59,60]. At the present time this corresponds to frequencies ω < h * /(2π).

V. CONCLUSIONS
We studied the GW spectrum produced in a strongly supercooled phase transition. We started from a lattice simulation of two-bubble collisions in order to model the evolution of the scalar field gradients which source GWs. We used a logarithmic potential typical for very strong phase transitions, and considered the impact of a complex scalar potential possessing an U (1) symmetry. We then simulated the production of GWs assuming that after collisions the gradients will continue moving at velocities near the speed of light losing energy according to the lattice results.
We found that the collision fronts disappear much more slowly in collisions of bubbles with different complex phases than in the case where the phases are equal. Therefore in the U (1) symmetric case the GW source is decaying significantly slower and the resulting GW spectrum is less steep at low frequencies than in the case of a real scalar field. Our results for the GW spectra are shown in Fig. 4, from which we see that the low-frequency behaviours are in these cases very different: In the case of a real scalar field the low-frequency spectrum resembles the envelope approximation result ∝ ω 3 , whereas in the U(1) symmetric case it is closer to the bulk flow result ∝ ω. In both cases the high-frequency power-law is ω −2 . We provided a simple broken power-law fits to these spectra, convenient for phenomenological studies.
Our treatment accounts only for the gradients at the bubble wall. Since the gradient energy is concentrated where the bubble wall would propagate also after the collisions, accounting for the field gradients in detail should result only in minor changes in the GW spectrum. However, to check this detailed lattice simulations are needed.