Gravitational waves from colliding vacuum bubbles in gauge theories

We study production of gravitational waves (GWs) in strongly supercooled cosmological phase transitions in gauge theories. We extract from two-bubble lattice simulations the scaling of the GW source, and use it in many-bubble simulations in the thin-wall limit to estimate the resulting GW spectrum. We find that in presence of the gauge field the GW source decays with bubble radius as $\propto R^{-3}$ after collisions. This leads to a GW spectrum that follows $\Omega_{\rm GW} \propto \omega^{2.3}$ at low frequencies and $\Omega_{\rm GW} \propto \omega^{-2.4}$ at high frequencies, marking a significant deviation from the popular envelope approximation.

Observations of stochastic GW backgrounds could allow us a glimpse of the very early Universe as many highenergy processes are predicted to be potential sources of such backgrounds. In this paper we will focus on cosmological first-order phase transitions, which are one example of such a source [22]. Many beyond Standard Model scenarios predict first-order phase transitions and a significant amount of work has already been put into the possibility of exploring them through GWs .
In a first-order phase transition the Universe starts in a metastable false vacuum. The transition proceeds via nucleation and subsequent expansion of bubbles of the true vacuum [72][73][74]. Eventually these bubbles collide and convert the whole Hubble volume into the new phase. In this process GWs are sourced by the bubble collisions [75][76][77][78][79][80] and plasma motions generated by the interactions of the plasma with the bubble walls [81][82][83][84][85][86][87]. In strongly supercooled transitions the former source dominates [77,88].
For the calculation of the GWs from colliding vacuum bubbles the equations of motion of the fields sourcing GWs need to be solved, requiring, in principle, 3D lattice simulations [76,79,89]. These simulations are computationally very expensive as very large simulation volumes are needed in order to simulate multiple bubbles, and very dense lattices to resolve the thinning bubble walls. Therefore, it is practical to develop approximations that provide a realistic description of the phase transition dynamics and an accurate estimate of the resulting GW spectrum, but are computationally less expensive than full 3D lattice simulations.
For a long time the envelope approximation, introduced in Ref. [75] and studied further in Refs. [90][91][92], has been used to estimate the GW spectrum sourced by the bubble collisions. In this approximation the collided parts of the bubble walls are completely neglected and the GW spectrum is calculated in the thin-wall limit. Improved modeling was developed in Refs. [93][94][95][96] as an attempt to model the behaviour of the plasma after the transition. Following a similar approach in Ref. [80] we developed a new estimate for the GW spectrum from bubble collisions by accounting for the scaling of the GW source after the collisions. Our estimate lead to a spectrum significantly different from the envelope approximation.
In this paper we consider a class of realistic models where bubble collisions can give the dominant contribution to the GW production. Furthermore, we describe breaking of a gauge U(1) symmetry, and study with lattice simulations the evolution of the scalar and gauge fields in two-bubble collisions. We find that the gradients in the complex phase of the scalar field are quickly damped after the collision by the gauge field. As a result, in gauge theories the GW source after the collision scales similarly to the case of just a real scalar, and the resulting GW spectrum follows ∝ ω 2.3 at low frequencies with a ∝ ω −2.4 fall above the peak.

II. PHASE TRANSITION
In order for the bubble collisions to give the dominant GW source, the phase transition has to be strongly supercooled [77,88]. Such strong supercooling is not typically realized in models with a polynomial scalar potential [77,85]. Instead, in models featuring classical scale invariance [27,36,38,40,44,50,52,60,61,[97][98][99][100][101][102][103][104][105] the transition can be so strongly supercooled that the interactions of the bubble wall with the plasma can be ne- The strength of the transition α, and the dimensionless parametersλ andg (see Eq. (A5)) as a function of the gauge coupling g. Different curves correspond to different values of v. The dashed curve in the right panel shows T = 0 limit ofg.
glected [77,88]. Many such models also include a gauge U(1) symmetry under which the scalar field is charged, and the dominant contribution on the effective potential arises from the gauge field loops. The phase transition in these models is therefore similar to that in classically conformal scalar electrodynamics, which we choose as a benchmark model. Scalar electrodynamics is described by the gauge U(1) symmetric Lagrangian where F µν = ∂ µ A ν − ∂ ν A µ and D µ = ∂ µ + igA µ are the electromagnetic field strength tensor and the gauge covariant derivative. In classically conformal models the tree-level scalar potential is quartic, V (|φ|) = λ|φ|/4. A non-trivial minimum is revealed when the radiative corrections are taken into account [106], and finite temperature effects induce a potential energy barrier between the symmetric and the symmetry-breaking minima. The one-loop effective potential of classically conformal scalar electrodynamics is where T denotes temperature of the plasma and v the vacuum expectation value of |φ| at T = 0. The symmetric and broken vacua are degenerate at a critical temperature T = T c . The bubble nucleation temperature T n < T c is defined as the temperature at which the probability of nucleating at least one bubble in a horizon volume in a Hubble time approaches unity [74]. In the left of Fig. 1 we show the parameter α ≡ ∆V (T = 0)/ρ rad (T ), that characterizes the strength of the transition, as a function of g for different values of v. We assume that only vacuum and radiation energy densities, ∆V (T = 0) and ρ rad (T ), contribute to the expansion rate, and approximate the effective number of relativistic degrees of freedom by its Standard Model value [107]. If α > 1 the transition finishes only after a vacuum energy dominated period. By strong supercooling we refer to α 1.
For the following analysis we define dimensionless parametersg andλ as such thatλ determines the shape of the scalar potential andg the strength of the coupling between the gauge field and the scalar field. In the middle and right panels of Fig. 1 we show these parameters at T = T n . For strongly supercooled transitionsg 2 ≈ 8π/(3g 2 ) andλ 2 ≈ 60/(g * (T n )α).

III. GRAVITATIONAL WAVE SOURCE
Next we study two-bubble collisions in order to find how the GW source scales after the collision. The total energy spectrum in a directionk at an angular frequency ω = | k| of the GWs emitted in the phase transition is given by [108] where Λ ij,lm is the transverse-traceless projection tensor. As Λ ij,lm δ ij = 0, the part of the stress energy tensor that is proportional to the metric tensor g µν does not contribute to formation of GWs. We therefore define T µν as (see Appendix A for an explicit form) The evolution of the system is governed by the equations of motion, given in the Lorentz gauge (∂ µ A µ = 0) 1 by 1 Our results are independent of the gauge choice because Tµν is gauge invariant. which we solve on a lattice starting from a configuration where A µ = 0 2 and two O(4) symmetric scalar field bubbles 3 have nucleated simultaneously with their centers lying on z-axis (see Appendix A for details of the lattice simulation). Then, along the collision axis only the zz component of T ij is non-zero, The bubble nucleation breaks the U(1) symmetry inside the bubble, as the complex phase of the scalar field, which we denote by ϕ (i.e. φ = |φ|e iϕ ), takes a value in the range 0 ≤ ϕ < 2π. 4 Eventually, as the bubbles expand, they will collide with bubbles containing different complex phases. Therefore, to get the average scaling of the GW source, we average T zz over simulations with different initial complex phase differences.
Our lattice simulations show that a T zz deviates from zero in a very narrow region around the bubble wall and this feature continues propagating almost at the speed of light after the collision. In Fig. 2 we show by the solid curves the scaling of the maximal T zz as a function of time, which much after nucleation is obtained roughly at z = ±d/2 ∓ t, where d denotes the distance between the bubble centers. Two important remarks are in order: First, we see that the steep drop after the collision becomes shorter asλ decreases. This can be traced to false vacuum trapping (field bouncing back to the false vacuum in the collision region) which becomes increasingly likely for larger values ofλ. 5 Second, the largerg is, the closer the behaviour of the GW source is to the case where the complex phases inside the colliding bubbles are equal, ∆ϕ = 0. Moreover, the smallerλ is, the faster the scaling reaches the ∆ϕ = 0 case as a function ofg. From Fig. 1 we see thatλ is very small,λ 1, andg is large, g > 10, in the region where supercooling is strong and the bubble collision signal can be the dominant contribution. As can be seen from the right panel of Fig. 2 the scaling in this case quickly reaches ∝ t −3 behaviour after the collision. 6 Instead, for example in the case of breaking of a global U(1) symmetry, corresponding tog = 0, ∝ t −2 scaling can be realized.

IV. GRAVITATIONAL WAVE SPECTRUM
Next, following a similar approach as one would using the envelope approximation, we perform many-bubble simulations in the thin-wall limit. Whereas in the envelope approach the collided parts of the walls are neglected, we instead use the scaling obtained in the previous section.
We consider an exponential bubble nucleation rate per unit volume, Γ ∝ e βt , and write the abundance of GWs produced in bubble collisions in a logarithmic frequency interval as [80] 5 In Ref. [80] we used larger values ofλ and found scaling resembling the left panel of Fig. 2. Here we focus onλ 1 which is more relevant for very strong transitions. 6 We have also checked that subsequent collisions do not change the scaling.   Table I. where gives the spectral shape of the GW background. The volume over which Ω GW is averaged is denoted by V s . The functions C +,× are fork = (0, 0, 1), in the thin-wall limit given by (see Appendix B for the derivation) where t n , z n and R n = t − t n denote the nucleation time, the z coordinate of the bubble nucleation center and the radius of the bubble n. The functions g +,× are defined as g + (φ x ) = cos(2φ x ) and g × (φ x ) = sin(2φ x ). The function f (R n ) accounts for the scaling of the GW source, following the results of our lattice simulations, which showed that the maximum of T zz scales as R −ξ n after the collision. The bubble radius at the collision moment, t = t c , is denoted by R n,c .
We calculate S by performing thin-wall simulations where we nucleate bubbles according to the rate Γ ∝ e βt a cubic box of size (7/β) 3 with periodic boundary conditions (see Appendix B for the details of the thin-wall simulations). We perform the angular integrals over the bubble surfaces by discretising each of them with 10 6 evenly distributed points. Our results are calculated from 60 simulations. We parametrize the results as a broken power-law, whereS andω are the peak amplitude and angular frequency of the spectrum, a, b > 0 are the low-and highfrequency slopes of the spectrum respectively, and c determines the width of the peak. We show the parameter values and their errors resulting from fits to our simulation results in Table I for ξ = 2, 3, 4 and the envelope approximation, 7 and illustrate these fits in Fig. 3. The spectrum today can be obtained from Eq. (8) by redshifting [80,81]. At super-horizon scales the spectrum scales as ω 3 as the source is diluted by the Hubble expansion [109,110]. For ξ = 3, corresponding to the case of breaking of a gauge U(1) symmetry, the low-and high-frequency tails of the spectrum are ∝ ω 2.3 and ∝ ω −2.4 . Instead, for ξ = 2, which can be realized for example in the case of breaking of a global U(1) symmetry, they are ∝ ω 1.0 and ∝ ω −2.2 . The spectrum peaks in both cases slightly below ω = β with an amplitude S ≈ 0.04. We find that increasing ξ brings the low-frequency power-law quickly closer to envelope result, a = 3.0, as shown by the ξ = 4 case. The high-frequency power-law instead seems to change very mildly for ξ > 3 not obviously converging to a slope that agrees with the envelope approximation.

V. CONCLUSIONS
Vacuum bubble collisions give the dominant source of GWs in a cosmological first-order phase transition if the transition is sufficiently strongly supercooled. This can be realized in classically conformal models. The simplest realistic examples of such involve breaking of a U(1) gauge symmetry. Motivated by these observations, we have studied the formation of GWs in a first-order phase transition in classically conformal scalar electrodynamics.
We have estimated the GW spectrum by first studying the scaling of the GW source in two-bubble lattice simulations, and then using that scaling in many-bubble simulations in the thin-wall limit. We have found that the presence of the gauge field brings the results close to the simple real scalar case where the GW source decays with the bubble size as ∝ R −3 . The resulting spectrum, shown by the green solid curve in Fig. 3, follows Ω GW ∝ ω 2.3 at low frequencies and Ω GW ∝ ω −2.4 at high frequencies. By calculating the transition temperature in classically conformal scalar electrodynamics we have shown that this limit withλ 1 andg 1 is realised in most of the parameter space of interest where the bubble collision signal can give the dominant contribution to the GW spectrum.
Ascertaining the shape of the signal is crucial as it could shine light on the underlying particle physics model. Sources associated to plasma dynamics, that dominate GW production in weaker transitions, in general predict high-frequency slopes [111] ω −4 from sound waves and ω −2/3 from turbulence. Our result shows that probing a signal that falls between the power-laws ω −2.5 and ω −2.1 at high frequencies would point to a very strong phase transition.
While our result should describe both real and gauged scalar field transitions, we have explored also different energy decay laws which could be realised in other models. Most notably, T zz ∝ R −2 could be realised in models where the U(1) symmetry is global (i.e.g = 0) or modified transition dynamics allowsg 1. In this case the resulting spectrum follows Ω GW ∝ ω at low and Ω GW ∝ ω −2.2 at high frequencies. Moreover, we have shown that the envelope result is unlikely to be able to describe realistic spectra especially at high frequencies. The project is co-financed by the Polish National Agency for Academic Exchange within Polish Returns Programme under agreement PPN/PPO/2020/1/00013/U/00001. This work was also supported by the grants FPA2017-88915-P and SEV-2016-0588. IFAE is partially funded by the CERCA program of the Generalitat de Catalunya.

Appendix A: Lattice simulation
We perform lattice simulations of two-bubble collisions, starting from a configuration where A µ = 0 and two O(4) symmetric scalar field bubbles have nucleated simultaneously at (x, y, z) = (0, 0, ±d/2). The radial profile of an O(4) symmetric initial configuration is obtained as the solution of with boundary conditions ∂ r |φ| = 0 at r = 0 and |φ| → 0 at r → ∞. A system of two simultaneously nucleated O(4) symmetric bubbles is conveniently described in coordinates (s, z, ψ, θ) defined via tan θ = x/y, t = s cosh ψ and r = s sinh ψ where r 2 = x 2 + y 2 . We consider the region t > r as this is where the bubbles collide (see Ref. [78] for details). The d'Alembertian in these coordinates reads For the lattice implementation, we write the equations of motion of the scalar and gauge fields in dimensionless where φ R and φ I are the real and imaginary parts of φ, defined such that φ = (φ R + iφ I )/ √ 2. The dimensionless scalar potential V = V /∆V , where ∆V denotes the vacuum energy difference between the symmetric and broken vacua at temperature T , can be written as Here we have defined dimensionless parametersg andλ as We solve the equations of motion (A3) numerically on a diamond-shaped sz lattice as in Ref. [112]. To ascertain the numerical stability of the simulation, we have checked that the gauge condition, ∂ µ A µ = 0 remains satisfied througout the simulation. We have also performed the simulations with different lattice spacings finding that the results are unchanged unless the grid is much less dense than what we use in the following results (δs = δz = 0.005).
In Fig. 4 we show the result from a simulation withλ = 0.04,g = 10, initial bubble separation d = 20 (in the dimensionless units) and initial complex phase difference ∆ϕ = π/2. The left panel shows the evolution of the complex phase of the scalar field. As can be seen from the equations of motion (A3), gradients in ϕ source the gauge field. Therefore, it is expected that the gauge field deviates from zero where the gradients in ϕ are large. We see this in the middle and right panels of Fig. 4, which show the gauge field components A s and A z : A sharp feature in the gauge field propagates roughly at the speed of light after collision.
In Fig. 5 we show the zz component of the stress-energy tensor, by the color coding for three different two-bubble collisions. In the left panelg = 0, and in the right panel the complex phase of the scalar field inside the colliding bubbles is the same, ∆ϕ = 0. In these cases only the scalar field gradients contribute to T zz , and the result agrees with the ones shown in Ref. [80]: If there is a complex phase difference between the colliding bubbles, the scalar field gradients propagate much longer after the collision than in the case where the complex phases are equal. The middle panel of Fig. 5 shows the case where the complex phases are different andg > 0. We see that the result in that case roughly matches with the ∆ϕ = 0 case. This can be understood as decay of the gradients in the complex phase of the scalar field to gauge fields.
The second crucial piece of information we get from Fig. 5 is that the gradients are well localised in space not only as the walls accelerate and become thinner but also after the time of the collision. In fact the spatial localisation of gradients becomes even more narrow as bubbles grow bigger before colliding. In a realistic transition the bubbles would grow many orders of magnitude in size before colliding which means it is well justified to assume the thickness of the walls and gradients after the collision is negligible compared to size of the colliding bubbles. This is the well known thin-wall approximation we will utilise in Appendix B. Next we generalize the treatment of Ref. [80] to the case where the stress-energy tensor is not given solely by the scalar field gradients after the bubble collisions. The Fourier transform of the stress-energy tensor is given by In the thin wall limit, by breaking the spatial integral into regions around each bubble nucleation center and taking k = (0, 0, ω), we get where t n denotes the nucleation time of the bubble n, z n the z coordinate of the bubble nucleation center, and R n = t − t n the bubble radius at time t > t n . It is convenient to write T ij (r) in a coordinate system where the z axis points to the radial direction,ẑ =r. The coordinate transformation of T ij (r) is If only zz component of is non zero, we get dr r 2 T ij (r) =x ixj dr r 2 T zz (r) .
In the thin-wall limit we approximate where L denotes the bubble wall width. Before the wall element in the solid angle dΩ collides with another bubble T zz is given solely by the scalar field gradients, T zz = |∂ z φ| 2 . By energy conservation where R n,c denotes the bubble radius at the time of the collision. Therefore, before the collision From the definition of f and using f (R n = R n,c ) = 1 we get that the wall width at the collision moment is Before the collision the bubble wall gets thinner as the Lorentz factor of the bubble wall increases, but after the collision the wall element moves roughly at a constant velocity as no more energy is injected into it. Therefore we assume that after the collision the wall thickness L remains constant, L = L c . We can then write the f (R n ) function after the collision as Since T ij ( k) is symmetric, we can write the transverse-traceless projection fork = (0, 0, 1) as The functions C + ≡ T 11 − T 22 and C × ≡ 2T 12 are given by C +,× (ω) ≈ 1 6π n tn dt dΩ x sin 2 θ x g +,× (φ x )R 3 n f (R n ) e iω(t−zn−Rn cos θx) , where g + (φ x ) = cos(2φ x ) and g × (φ x ) = sin(2φ x ). The total energy spectrum in a directionk at an angular frequency ω = | k| of the GWs emitted in the phase transition is given by [108] dE dΩ k dω = 2Gω 2 Λ ij,lm (k)T * ij ( k)T lm ( k) .
Using Eq. (B10) and the definition α ≡ ∆V /ρ rad we write the abundance of GWs produced in bubble collisions in a logarithmic frequency interval as where E tot = V s (ρ rad + ∆V ) and S(ω) = ω β 3 3β 5 8πV s dΩ k |C + (ω)| 2 + |C × (ω)| 2 (B14) gives the spectral shape of the GW background. Here V s denotes the volume over which Ω GW (ω) is averaged. We consider exponential bubble nucleation rate, Γ ∝ e βt , which implies that dΩ k [|C + (ω)| 2 + |C × (ω)| 2 ] ∝ V s /β 5 . We simulate the phase transition by nucleating bubbles according to the exponential bubble nucleation rate inside a cubic box with periodic boundary conditions. Following the thin-wall approximation, we simulate the bubbles as spherical shells. We discretise the bubble surfaces and find the time when each of these bubble wall elements collides with another bubble wall by bisection method. The corresponding radius is denoted by R n = R n,c . Once we know R n,c for each bubble wall element of each bubble in the simulation, we integrate the functions C +,× (ω) for a given value of ω. We note that if f (R) is a (broken) power-law the temporal integral can be performed analytically. The spectral function S(ω) is then simply obtained by integrating over thek directions. In practice, since our simulation box is cubic, the integral overk directions is done by summing over 6 directions, corresponding to the normal vectors of the faces of the cube, with equal weights 2π/3. Finally, to reduce the errors, we calculate the final result by averaging the spectrum over many simulations with different randomly generated bubble nucleation histories.