1 Introduction

We are currently witnessing the dawn of a new era in astrophysics and cosmology, started by the LIGO/Virgo observations of gravitational waves (GWs) from black hole mergers [1, 2]. Many experiments are planned to further explore GWs in a broad frequency range in the coming decades [3,4,5,6,7,8,9,10]. In addition to transient GW signals, such as those from black hole mergers, these experiments are able to probe stochastic GW backgrounds. In fact, recent results from NANOGrav pulsar timing observations [11] may already indicate the first observation of a stochastic GW background [12,13,14,15,16,17,18,19,20,21].

Observations of stochastic GW backgrounds could allow us a glimpse of the very early Universe as many high-energy 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 [23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71].

Fig. 1
figure 1

The strength of the transition \(\alpha \), and the dimensionless parameters \(\tilde{\lambda }\) and \({\tilde{g}}\) (see Eq. (3A5)) 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 of \({\tilde{g}}\)

Full size image

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 \(\propto \omega ^{2.3}\) at low frequencies with a \(\propto \omega ^{-2.9}\) fall above the peak.

2 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 neglected [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

$$\begin{aligned} {\mathcal {L}} = -\frac{1}{4} (F_{\mu \nu })^2 + |D_\mu \phi |^2 - V(|\phi |) \,, \end{aligned}$$
(1)

where \(F_{\mu \nu } = \partial _\mu A_\nu - \partial _\nu A_\mu \) and \(D_\mu = \partial _\mu + ig A_\mu \) are the electromagnetic field strength tensor and the gauge covariant derivative. In classically conformal models the tree-level scalar potential is quartic, \(V(|\phi |) = \lambda |\phi |/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

$$\begin{aligned} V(|\phi |) = \frac{g^2}{2} T^2 |\phi |^2 + \frac{3 g^4}{4\pi ^2} |\phi |^4 \left[ \ln \frac{|\phi |^2}{v^2} - \frac{1}{2} \right] \,, \end{aligned}$$
(2)

where T denotes temperature of the plasma and v the vacuum expectation value of \(|\phi |\) 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 \(\alpha \equiv \Delta V(T=0)/\rho _{\mathrm{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, \(\Delta V(T=0)\) and \(\rho _{\mathrm{rad}}(T)\), contribute to the expansion rate, and approximate the effective number of relativistic degrees of freedom by its Standard Model value [107]. If \(\alpha >1\) the transition finishes only after a vacuum energy dominated period. By strong supercooling we refer to \(\alpha \gg 1\).

For the following analysis we define dimensionless parameters \({\tilde{g}}\) and \(\tilde{\lambda }\) as

$$\begin{aligned} {\tilde{g}} = \frac{g v^2}{\sqrt{\Delta V}} \,, \qquad \tilde{\lambda } = \frac{g^2 v^2 T^2}{2 \Delta V} \,, \end{aligned}$$
(3)

such that \(\tilde{\lambda }\) determines the shape of the scalar potential and \({\tilde{g}}\) 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 transitions \({\tilde{g}}^2 \approx 8\pi /(3 g^2)\) and \(\tilde{\lambda }^2 \approx 60/(g_*(T_n)\alpha )\).

Fig. 2
figure 2

Evolution of the GW source in collision of two bubbles averaged over simulations with different initial complex phase differences. The collision happens at \(t=t_c\). Solid curves correspond to different values of \({\tilde{g}}\) and the dashed curve to \(\Delta \varphi =0\). The dotted lines show \(\propto t^{-2}\) and \(\propto t^{-3}\) power-laws

Full size image

3 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 direction \({\hat{k}}\) at an angular frequency \(\omega = |\vec {k}|\) of the GWs emitted in the phase transition is given by [108]

$$\begin{aligned} \frac{{\mathrm{d}}E}{{\mathrm{d}}\Omega _k{\mathrm{d}}\omega } = 2G\omega ^2 \Lambda _{ij,lm}({\hat{k}}) T_{ij}^*(\vec {k}) T_{lm}(\vec {k}) \,, \end{aligned}$$
(4)

where \(\Lambda _{ij,lm}\) is the transverse-traceless projection tensor. As \(\Lambda _{ij,lm}\delta _{ij} = 0\), the part of the stress energy tensor that is proportional to the metric tensor \(g_{\mu \nu }\) does not contribute to formation of GWs. We therefore define \(T_{\mu \nu }\) as (see Appendix A for an explicit form)

$$\begin{aligned} T_{\mu \nu } \equiv \left( \frac{\partial {\mathcal {L}}}{\partial (\partial ^\mu \phi )} \partial _\nu \phi + \mathrm{c.c}\right) + \frac{\partial {\mathcal {L}}}{\partial (\partial ^\mu A_\alpha )} \partial _\nu A_\alpha \,. \end{aligned}$$
(5)

The evolution of the system is governed by the equations of motion, given in the Lorentz gauge (\(\partial _\mu A^\mu = 0\))Footnote 1 by

$$\begin{aligned} \begin{aligned}&\Box A_\mu = i g (\phi ^* \partial _\mu \phi - \phi \partial _\mu \phi ^*) - 2 g^2 A_\mu |\phi |^2 \,, \\&\Box \phi + \frac{{\mathrm{d}}V}{{\mathrm{d}}\phi ^*} = -i2g A_\mu \partial ^\mu \phi + g^2 A^2 \phi \,, \end{aligned}\end{aligned}$$
(6)

which we solve on a lattice starting from a configuration where \(A_\mu = 0\)Footnote 2 and two O(4) symmetric scalar field bubblesFootnote 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,

$$\begin{aligned} \begin{aligned} T_{zz}&= \,2|\partial _z \phi | - (\partial _z A_t - \partial _t A_z) \partial _z A_t \\&\quad - i g A_z (\phi ^*\partial _z\phi - \phi \partial _z\phi ^*) \,. \end{aligned} \end{aligned}$$
(7)

The bubble nucleation breaks the U(1) symmetry inside the bubble, as the complex phase of the scalar field, which we denote by \(\varphi \) (i.e. \(\phi = |\phi | e^{i\varphi }\)), takes a value in the range \(0\le \varphi < 2\pi \).Footnote 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.

Table 1 Fitted values for the parametrization of the spectral shape (12)
Full size table

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 = \pm d/2 \mp 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 \(\tilde{\lambda }\) 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 \(\tilde{\lambda }\).Footnote 5 Second, the larger \({\tilde{g}}\) is, the closer the behaviour of the GW source is to the case where the complex phases inside the colliding bubbles are equal, \(\Delta \varphi =0\). Moreover, the smaller \(\tilde{\lambda }\) is, the faster the scaling reaches the \(\Delta \varphi =0\) case as a function of \({\tilde{g}}\). From Fig. 1 we see that \(\tilde{\lambda }\) is very small, \(\tilde{\lambda }\ll 1\), and \({\tilde{g}}\) is large, \({\tilde{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 \(\propto t^{-3}\) behaviour after the collision.Footnote 6 Instead, for example in the case of breaking of a global U(1) symmetry, corresponding to \({\tilde{g}}=0\), \(\propto t^{-2}\) scaling can be realized.

4 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, \(\Gamma \propto e^{\beta t}\), and write the abundance of GWs produced in bubble collisions in a logarithmic frequency interval as [80]

$$\begin{aligned} \Omega _{\mathrm{GW}}(\omega ) \equiv \frac{1}{E_{\mathrm{tot}}}\frac{{\mathrm{d}}E}{{\mathrm{d}}\ln \omega } = \left( \frac{H}{\beta }\right) ^2\left( \frac{\alpha }{1+\alpha }\right) ^2 S(\omega ) \,, \end{aligned}$$
(8)

where

$$\begin{aligned} S(\omega ) \!=\! \left( \frac{\omega }{\beta }\right) ^3 \frac{3\beta ^5}{8\pi V_s} \int \!{\mathrm{d}}\Omega _k \left[ |C_+(\omega )|^2 + |C_\times (\omega )|^2 \right] \end{aligned}$$
(9)

gives the spectral shape of the GW background. The volume over which \(\Omega _{\mathrm{GW}}\) is averaged is denoted by \(V_s\). The functions \(C_{+,\times }\) are for \({\hat{k}} = (0,0,1)\), in the thin-wall limit given by (see Appendix B for the derivation)

$$\begin{aligned} \begin{aligned} C_{+,\times }(\omega )&\approx \frac{1}{6\pi } \sum _n \int _{t_n} {\mathrm{d}}t\, {\mathrm{d}}\Omega _x\, \sin ^2\theta _x\, g_{+,\times }(\phi _x) \\&\quad \times R_n^3 f(R_n) \,e^{i\omega (t - z_n - R_n\cos \theta _x)} \,, \end{aligned}\end{aligned}$$
(10)

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_{+,\times }\) are defined as \(g_+(\phi _x) = \cos (2\phi _x)\) and \(g_\times (\phi _x) = \sin (2\phi _x)\). The function \(f(R_n)\) accounts for the scaling of the GW source,

$$\begin{aligned} f(R_n) = \min \left[ 1,\left( R_{n,c}/R_n\right) ^{\xi +1}\right] \,, \end{aligned}$$
(11)

following the results of our lattice simulations, which showed that the maximum of \(T_{zz}\) scales as \(R_n^{-\xi }\) after the collision. The bubble radius at the collision moment, \(t=t_c\), is denoted by \(R_{n,c}\).

Fig. 3
figure 3

The spectral shape of GWs (see Eq. (8B13)) from vacuum bubble collisions. The curves show broken power-law fits to the simulation results for different decay-laws of the GW source after collisions and in the envelope approximation. The solid curve is realised in the case of breaking of a gauge U(1) symmetry. The corresponding parameters of the fit, and their errors, are given in Table 1

Full size image

We calculate S by performing thin-wall simulations where we nucleate bubbles according to the rate \(\Gamma \propto e^{\beta t}\) in a cubic box of size \(12/\beta \) with periodic boundary conditions (see Appendix B for the details of the thin-wall simulations). Our results are calculated from 40 simulations. We parametrize the results as a broken power-law,

$$\begin{aligned} S_{\mathrm{fit}}(\omega ) =\frac{{\bar{S}}\,(a+b)^c}{\left[ b \left( \frac{\omega }{\bar{\omega }}\right) ^{-a/c} + a \left( \frac{\omega }{\bar{\omega }}\right) ^{b/c}\right] ^c} \,, \end{aligned}$$
(12)

where \({\bar{S}}\) and \(\bar{\omega }\) are the peak amplitude and angular frequency of the spectrum, \(a,b>0\) are the low- and high-frequency 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 1 for \(\xi =2,3,4\) and the envelope approximation,Footnote 7 and illustrate these fits in Fig. 3. The spectrum today can be obtained from Eq. (8B13) by red-shifting [80, 81]. At super-horizon scales the spectrum scales as \(\omega ^3\) as the source is diluted by the Hubble expansion [109, 110].

For \(\xi =3\), corresponding to the case of breaking of a gauge U(1) symmetry, the low- and high-frequency tails of the spectrum are \(\propto \omega ^{2.3}\) and \(\propto \omega ^{-2.9}\). Instead, for \(\xi =2\), which can be realized for example in the case of breaking of a global U(1) symmetry, they are \(\propto \omega ^{1.0}\) and \(\propto \omega ^{-2.6}\). The spectrum peaks in both cases slightly below \(\omega = \beta \) with an amplitude \(S\approx 0.03\). We find that increasing \(\xi \) brings the low-frequency power-law quickly closer to envelope result, \(a=3.0\), as shown by the \(\xi =4\) case. The high-frequency power-law instead seems to decrease very slowly for \(\xi >3\).Footnote 8 Getting a high-frequency tail that agrees with the envelope approximation, \(\propto \omega ^{-0.9}\), requires an extremely rapid decay of the GW source after the bubble collisions.

5 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 \(\propto R^{-3}\). The resulting spectrum, shown by the green solid curve in Fig. 3, follows \(\Omega _{\mathrm{GW}}\propto \omega ^{2.3}\) at low frequencies and \(\Omega _{\mathrm{GW}} \propto \omega ^{-2.9}\) at high frequencies. By calculating the transition temperature in classically conformal scalar electrodynamics we have shown that this limit with \(\tilde{\lambda } \ll 1\) and \({\tilde{g}} \gg 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 to the underlying particle physics model. Our result shows that probing a signal that falls roughly as \(\omega ^{-3}\) at high frequencies would point to a very strong phase transition. Sources associated to plasma dynamics, that dominate GW production in weaker transitions, in general predict different high-frequency slopes [111]: \(\omega ^{-4}\) from sound waves and \(\omega ^{-2/3}\) from turbulence. Moreover, compared to the envelope approximation, our result indicates that observing the bubble collisions signal in multiple GW detectors at different frequencies will be less likely as the spectrum is narrower than previously thought.

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}\propto R^{-2}\) could be realised in models where the U(1) symmetry is global (i.e. \({\tilde{g}} = 0\)) or modified transition dynamics allows \({\tilde{g}}\ll 1\). In this case the resulting spectrum follows \(\Omega _{\mathrm{GW}}\propto \omega \) at low and \(\Omega _{\mathrm{GW}}\propto \omega ^{2.6}\) at high frequencies. Moreover, we have shown that the envelope result is unlikely to be able to describe realistic spectra especially at high frequencies.