Gravitational waves from a supercooled electroweak phase transition and   their detection with pulsar timing arrays

Archil Kobakhidze; Cyril Lagger; Adrian Manning; Jason Yue

arXiv:1703.06552·hep-ph·November 9, 2017

Gravitational waves from a supercooled electroweak phase transition and their detection with pulsar timing arrays

Archil Kobakhidze, Cyril Lagger, Adrian Manning, Jason Yue

PDF

TL;DR

This paper explores the gravitational wave background generated by a supercooled electroweak phase transition, highlighting the potential for detection by pulsar timing arrays and the impact of extreme supercooling on the spectrum.

Contribution

It introduces a detailed analysis of gravitational waves from a supercooled electroweak phase transition, including a model with a cubic Higgs coupling and a comprehensive treatment of the transition dynamics.

Findings

01

Gravitational wave spectrum peaks at very low frequencies (~10^{-9}-10^{-7} Hz).

02

Certain parameter ranges make the spectrum detectable by future pulsar timing arrays.

03

The scenario can occur in models with a non-linear realization of the electroweak gauge group.

Abstract

We investigate the properties of a stochastic gravitational wave background produced by a first-order electroweak phase transition in the regime of extreme supercooling. We study a scenario whereby the percolation temperature that signifies the completion of the transition, $T_{p}$ , can be as low as a few MeV (nucleosynthesis temperature), while most of the true vacuum bubbles are formed much earlier at the nucleation temperature, $T_{n} \sim 50$ GeV. This implies that the gravitational wave spectrum is mainly produced by the collisions of large bubbles and characterised by a large amplitude and a peak frequency as low as $f \sim 1 0^{- 9} - 1 0^{- 7}$ Hz. We show that such a scenario can occur in (but not limited to) a model based on a non-linear realisation of the electroweak gauge group, such that the Higgs vacuum configuration is altered by a cubic coupling. In order to carefully quantify the…

Tables1

Table 1. Table 1 : Key parameters describing the phase transition for κ = [ − 1.92 , − 1.87 ] m h 2 v 𝜅 1.92 1.87 superscript subscript 𝑚 ℎ 2 𝑣 \kappa=[-1.92,-1.87]\frac{m_{h}^{2}}{v} . This range of κ 𝜅 \kappa values has been selected in prevision of its relevance for GW production.

$κ$ $[m_{h}^{2} / \| v \|]$	$T_{c}$ GeV	$T_{n}$ GeV	$T_{p}$ GeV	${(\bar{R} H_{p})}^{- 1}$	$ϵ_{kin} / ϵ_{rad}$
$- 1.87$	98.3	48.9	43.4	8.79	0.57
$- 1.88$	98.0	48.9	31.2	2.76	1.88
$- 1.89$	97.7	49.0	14.4	1.41	37.8
$- 1.9$	97.4	48.7	4.21	1.09	$5.09 \cdot 10^{3}$
$- 1.91$	97.1	48.6	0.977	1.02	$1.73 \cdot 10^{6}$
$- 1.92$	96.8	48.5	0.205	1.00	$8.80 \cdot 10^{8}$

Equations75

H (x) = \frac{ρ ( x )}{2} e^{\frac{i}{2} π^{i} (x) T^{i}} (01), i \in {1, 2, 3}

H (x) = \frac{ρ ( x )}{2} e^{\frac{i}{2} π^{i} (x) T^{i}} (01), i \in {1, 2, 3}

V^{(0)} (ρ) = - \frac{μ ^{2}}{2} ρ^{2} + \frac{κ}{3} ρ^{3} + \frac{λ}{4} ρ^{4} .

V^{(0)} (ρ) = - \frac{μ ^{2}}{2} ρ^{2} + \frac{κ}{3} ρ^{3} + \frac{λ}{4} ρ^{4} .

μ^{2}

μ^{2}

λ

V (ρ, T) = V^{(0)} (ρ) + V_{C W}^{(1)} (ρ) + V^{(1)} (ρ, T) + V_{D ai sy} (ρ, T),

V (ρ, T) = V^{(0)} (ρ) + V_{C W}^{(1)} (ρ) + V^{(1)} (ρ, T) + V_{D ai sy} (ρ, T),

Γ (T) \approx A (T) e^{- S (T)}

Γ (T) \approx A (T) e^{- S (T)}

S [ρ, T] = 4 π \int_{0}^{β} d τ \int_{0}^{\infty} d r r^{2} [\frac{1}{2} (\frac{d ρ}{d τ})^{2} + \frac{1}{2} (\frac{d ρ}{d r})^{2} + \tilde{F} (ρ, T)],

S [ρ, T] = 4 π \int_{0}^{β} d τ \int_{0}^{\infty} d r r^{2} [\frac{1}{2} (\frac{d ρ}{d τ})^{2} + \frac{1}{2} (\frac{d ρ}{d r})^{2} + \tilde{F} (ρ, T)],

\frac{\partial ^{2} ρ}{\partial τ ^{2}} + \frac{\partial ^{2} ρ}{\partial r ^{2}} + \frac{2}{r} \frac{\partial ρ}{\partial r} - \frac{\partial F ~}{\partial ρ} (ρ, T) = 0,

\frac{\partial ^{2} ρ}{\partial τ ^{2}} + \frac{\partial ^{2} ρ}{\partial r ^{2}} + \frac{2}{r} \frac{\partial ρ}{\partial r} - \frac{\partial F ~}{\partial ρ} (ρ, T) = 0,

\frac{\partial ρ}{\partial τ}_{τ = 0, \pm β /2} = 0, \frac{\partial ρ}{\partial r}_{r = 0} = 0, r \to \infty lim ρ (r) = v_{T}^{(+)} .

\frac{\partial ρ}{\partial τ}_{τ = 0, \pm β /2} = 0, \frac{\partial ρ}{\partial r}_{r = 0} = 0, r \to \infty lim ρ (r) = v_{T}^{(+)} .

S[\rho,T]\approx\left\{\begin{array}[]{llrl}S_{4}[\rho,T]&=&\displaystyle 2\pi^{2}\int^{\infty}_{0}d\tilde{r}\ \tilde{r}^{3}\left[\frac{1}{2}\left(\frac{d\rho}{d\tilde{r}}\right)^{2}+\tilde{\mathcal{F}}(\rho,T)\right],&T\ll R_{0}^{-1}\\ \frac{1}{T}S_{3}[\rho,T]&=&\displaystyle\frac{4\pi}{T}\int^{\infty}_{0}dr\ r^{2}\left[\frac{1}{2}\left(\frac{d\rho}{dr}\right)^{2}+\tilde{\mathcal{F}}(\rho,T)\right],&T\gg R_{0}^{-1}\end{array}\right.

S[\rho,T]\approx\left\{\begin{array}[]{llrl}S_{4}[\rho,T]&=&\displaystyle 2\pi^{2}\int^{\infty}_{0}d\tilde{r}\ \tilde{r}^{3}\left[\frac{1}{2}\left(\frac{d\rho}{d\tilde{r}}\right)^{2}+\tilde{\mathcal{F}}(\rho,T)\right],&T\ll R_{0}^{-1}\\ \frac{1}{T}S_{3}[\rho,T]&=&\displaystyle\frac{4\pi}{T}\int^{\infty}_{0}dr\ r^{2}\left[\frac{1}{2}\left(\frac{d\rho}{dr}\right)^{2}+\tilde{\mathcal{F}}(\rho,T)\right],&T\gg R_{0}^{-1}\end{array}\right.

\frac{d ^{2} ρ}{d r ^{2}} + \frac{α}{r} \frac{d ρ}{d r} - \frac{\partial F ~}{\partial ρ} (ρ, T) = 0, \frac{d ρ}{d r}_{r = 0} = 0, r \to \infty lim ρ (r) = v_{T}^{(+)}

\frac{d ^{2} ρ}{d r ^{2}} + \frac{α}{r} \frac{d ρ}{d r} - \frac{\partial F ~}{\partial ρ} (ρ, T) = 0, \frac{d ρ}{d r}_{r = 0} = 0, r \to \infty lim ρ (r) = v_{T}^{(+)}

A(T)\approx\left\{\begin{array}[]{ll}\displaystyle\frac{1}{R_{0}^{4}}\left(\frac{S_{4}(T)}{2\pi}\right)^{2}&,T\ll R_{0}^{-1}\\ \displaystyle T^{4}\left(\frac{S_{3}(T)}{2\pi T}\right)^{3/2}&,T\gg R_{0}^{-1}\end{array}\right.

A(T)\approx\left\{\begin{array}[]{ll}\displaystyle\frac{1}{R_{0}^{4}}\left(\frac{S_{4}(T)}{2\pi}\right)^{2}&,T\ll R_{0}^{-1}\\ \displaystyle T^{4}\left(\frac{S_{3}(T)}{2\pi T}\right)^{3/2}&,T\gg R_{0}^{-1}\end{array}\right.

p (t) = exp - := I (t) \frac{4 π}{3} \int_{t_{c}}^{t} d t^{'} Γ (t^{'}) a^{3} (t^{'}) r^{3} (t, t^{'}),

p (t) = exp - := I (t) \frac{4 π}{3} \int_{t_{c}}^{t} d t^{'} Γ (t^{'}) a^{3} (t^{'}) r^{3} (t, t^{'}),

r (t, t^{'}) = \int_{t^{'}}^{t} d t^{''} \frac{v ( t ^{''} )}{a ( t ^{''} )},

r (t, t^{'}) = \int_{t^{'}}^{t} d t^{''} \frac{v ( t ^{''} )}{a ( t ^{''} )},

\frac{d N}{d R} (t, t_{R}) = Γ (t_{R}) (\frac{a ( t _{R} )}{a ( t )})^{4} \frac{p ( t _{R} )}{v ( t _{R} )} .

\frac{d N}{d R} (t, t_{R}) = Γ (t_{R}) (\frac{a ( t _{R} )}{a ( t )})^{4} \frac{p ( t _{R} )}{v ( t _{R} )} .

\frac{d}{d t _{R}} (\frac{d N}{d R} (t_{p}, t_{R}))_{t_{R} = t_{n}} = 0.

\frac{d}{d t _{R}} (\frac{d N}{d R} (t_{p}, t_{R}))_{t_{R} = t_{n}} = 0.

E_{kin} = 4 π \int_{t_{n}}^{t_{p}} d t \frac{d R}{d t} (t, t_{n}) R^{2} (t, t_{n}) \tilde{ϵ} (t)

E_{kin} = 4 π \int_{t_{n}}^{t_{p}} d t \frac{d R}{d t} (t, t_{n}) R^{2} (t, t_{n}) \tilde{ϵ} (t)

p = \frac{1}{3} ϵ, a (t) \propto t^{1/2}, t = := ξ (\frac{45 M _{p}^{2}}{16 π ^{3} g _{⋆}})^{1/2} \frac{1}{T ^{2}},

p = \frac{1}{3} ϵ, a (t) \propto t^{1/2}, t = := ξ (\frac{45 M _{p}^{2}}{16 π ^{3} g _{⋆}})^{1/2} \frac{1}{T ^{2}},

R (T, T^{'})

R (T, T^{'})

p (T)

E_{kin}

\frac{β}{H ( T ~ _{n} )} = \tilde{T}_{n} \frac{d}{d T} (\frac{S _{3} ( T )}{T})_{T = \tilde{T}_{n}} .

\frac{β}{H ( T ~ _{n} )} = \tilde{T}_{n} \frac{d}{d T} (\frac{S _{3} ( T )}{T})_{T = \tilde{T}_{n}} .

h^{2} Ω_{G W} (f) = \frac{h ^{2}}{ϵ _{c}} \frac{d ϵ _{G W}}{d ( ln f )},

h^{2} Ω_{G W} (f) = \frac{h ^{2}}{ϵ _{c}} \frac{d ϵ _{G W}}{d ( ln f )},

h^{2} Ω_{GW} ≃ h^{2} Ω_{co l} + h^{2} Ω_{s w} + h^{2} Ω_{MHD} .

h^{2} Ω_{GW} ≃ h^{2} Ω_{co l} + h^{2} Ω_{s w} + h^{2} Ω_{MHD} .

h^{2} Ω_{GW} ≃ h^{2} Ω_{co l} .

h^{2} Ω_{GW} ≃ h^{2} Ω_{co l} .

\begin{array}[]{cll}f_{0}&=\displaystyle f_{p}\frac{a(t_{p})}{a(t_{0})}&=\displaystyle 1.65\times 10^{-7}\ \text{Hz}\left(\frac{f_{p}}{H_{p}}\right)\left(\frac{T_{p}}{1\ \text{GeV}}\right)\left(\frac{g_{*}}{100}\right)^{1/6}\\ \Omega_{GW,0}&=\displaystyle\Omega_{col}\left(\frac{a(t_{p})}{a(t_{0})}\right)^{4}\left(\frac{H_{p}}{H_{0}}\right)^{2}&\displaystyle=1.67\times 10^{-5}h^{-2}\left(\frac{100}{g_{*}}\right)^{1/3}\Omega_{col}.\end{array}

\begin{array}[]{cll}f_{0}&=\displaystyle f_{p}\frac{a(t_{p})}{a(t_{0})}&=\displaystyle 1.65\times 10^{-7}\ \text{Hz}\left(\frac{f_{p}}{H_{p}}\right)\left(\frac{T_{p}}{1\ \text{GeV}}\right)\left(\frac{g_{*}}{100}\right)^{1/6}\\ \Omega_{GW,0}&=\displaystyle\Omega_{col}\left(\frac{a(t_{p})}{a(t_{0})}\right)^{4}\left(\frac{H_{p}}{H_{0}}\right)^{2}&\displaystyle=1.67\times 10^{-5}h^{-2}\left(\frac{100}{g_{*}}\right)^{1/3}\Omega_{col}.\end{array}

Ω_{co l} \sim (\overset{ˉ}{R} H_{p})^{2} \frac{ϵ _{kin}^{2}}{ϵ _{t o t}^{2}} \sim (\overset{ˉ}{R} H_{p})^{2} \frac{α ^{2}}{( 1 + α ) ^{2}} .

Ω_{co l} \sim (\overset{ˉ}{R} H_{p})^{2} \frac{ϵ _{kin}^{2}}{ϵ _{t o t}^{2}} \sim (\overset{ˉ}{R} H_{p})^{2} \frac{α ^{2}}{( 1 + α ) ^{2}} .

h^{2} Ω_{co l} (f) = 1.67 \times 1 0^{- 5} (\frac{100}{g _{*}})^{1/3} (\frac{β}{H _{p}})^{- 2} κ_{ρ}^{2} (\frac{α}{1 + α})^{2} (\frac{0.11 v ^{3}}{0.42 + v ^{2}}) S (f),

h^{2} Ω_{co l} (f) = 1.67 \times 1 0^{- 5} (\frac{100}{g _{*}})^{1/3} (\frac{β}{H _{p}})^{- 2} κ_{ρ}^{2} (\frac{α}{1 + α})^{2} (\frac{0.11 v ^{3}}{0.42 + v ^{2}}) S (f),

S (f)

S (f)

f_{0}

V_{C W}^{(1)} (ρ) = i = W, Z, t, h \sum n_{i} \frac{m _{i}^{4} ( ρ )}{64 π ^{2}} (ln (\frac{m _{i}^{2} ( ρ )}{v ^{2}}) - \frac{3}{2}) .

V_{C W}^{(1)} (ρ) = i = W, Z, t, h \sum n_{i} \frac{m _{i}^{4} ( ρ )}{64 π ^{2}} (ln (\frac{m _{i}^{2} ( ρ )}{v ^{2}}) - \frac{3}{2}) .

V^{(1)} (ρ, T)

V^{(1)} (ρ, T)

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

Gravitational waves from a supercooled electroweak phase transition and their detection with pulsar timing arrays

Archil Kobakhidze1, Cyril Lagger1, Adrian Manning1 and Jason Yue2

Abstract

We investigate the properties of a stochastic gravitational wave background produced by a first-order electroweak phase transition in the regime of extreme supercooling. We study a scenario whereby the percolation temperature that signifies the completion of the transition, $T_{p}$ , can be as low as a few MeV (nucleosynthesis temperature), while most of the true vacuum bubbles are formed much earlier at the nucleation temperature, $T_{n}\sim 50$ GeV. This implies that the gravitational wave spectrum is mainly produced by the collisions of large bubbles and characterised by a large amplitude and a peak frequency as low as $f\sim 10^{-9}-10^{-7}$ Hz. We show that such a scenario can occur in (but not limited to) a model based on a non-linear realisation of the electroweak gauge group, such that the Higgs vacuum configuration is altered by a cubic coupling. In order to carefully quantify the evolution of the phase transition of this model over such a wide temperature range, we go beyond the usual fast transition approximation, taking into account the expansion of the Universe as well as the behaviour of the nucleation probability at low temperatures. Our computation shows that there exists a range of parameters for which the gravitational wave spectrum lies at the edge between the exclusion limits of current pulsar timing array experiments and the detection band of the future Square Kilometre Array observatory.

1 Introduction

Cosmological phase transitions (PT) are predicted by many particle physics models with important consequences on the dynamics of the Universe. In particular, first-order PTs produce a stochastic background of gravitational waves (GWs) from the collision of true vacuum bubbles and their interaction with the surrounding hot plasma [1, 2, 3, 4, 5]. The observation of a stochastic GW spectrum would then provide an opportunity to obtain more information on the early Universe and potentially new physics. In particular, the Standard Model itself does not predict a first-order electroweak PT [6] and thus neither the production of the associated GWs. However, several extensions of the SM do accommodate such a transition, allowing the matter-antimatter asymmetry of the Universe to be explained via electroweak baryogenesis [7] (cf. e.g. [8, 9, 10, 11] for recent reviews on the topic).

The peak frequency of a stochastic GW background produced by a PT near the electroweak scale, $T_{\text{EW}}\sim 100$ GeV, is expected to lie in the millihertz range, which coincides with the projected sensitivity of the future eLISA space-based interferometer [12]. This motivated a series of investigations into the production of GWs in various BSM models, see e.g. [13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26]. The characteristic frequency and amplitude of the spectrum are derived from the dynamics of the PT and depend on a few key parameters: the duration of the transition, the size of colliding bubbles, the bubble-walls velocity and the fraction of vacuum energy transferred into the bubble-walls. In the aforementioned studies, these quantities are computed under the assumption that the PT occurs on a time scale much shorter than the Hubble time. The instant at which most of the bubbles are nucleated is thus very close to the time when they collide and cover a significant volume of the Universe. In this article, we shall however consider the case of a prolonged electroweak PT with a non-negligible amount of time between nucleation and collision. We expect the GW background to be predominantly produced by large bubbles colliding much later than in a typical electroweak PT previously discussed in the literature. A range of lower peak frequencies is then observed in such a scenario.

In order to illustrate this phenomenon, we consider a model based on a non-linearly realised electroweak gauge group [27, 28, 29, 26] where the Higgs potential admits a cubic term at tree-level. In other words, a barrier exists between the two different phases of the Higgs field from the electroweak scale down to zero temperature, allowing a significant amount of supercooling. Interestingly, this model exhibits a range of parameters for which the PT is long-lasting, meaning that most of the true vacuum bubbles are nucleated around $T\sim 50$ GeV but collide well below the electroweak scale, as low as $T\sim[0.1,10]$ GeV. Precise results depend on the exact equation of state of the Universe which is complicated to compute in this context. Naively, vacuum energy is expected to dominate over radiation energy below a given temperature, potentially leading to an inflationary stage. However, we shall argue that such a scenario is unlikely to happen as a significant amount of bubbles are produced early enough (during the radiation dominated era) and subsequently act as a source of inhomogeneity which prevents inflation to occur. Therefore, our model differs from previous studies of scale-invariant models [30, 31, 32] in which the nucleation of true vacuum bubbles occurs at very low temperatures, namely after inflation started. We should note that extreme supercooling [33, 34] is also a feature of another class of scale-invariant models of electroweak symmetry breaking with a very light scalar particle [35].

Once the details of the PT are known, rough estimates of the peak frequency and peak amplitude of the GW spectrum can be derived from dimensional arguments, although more precise predictions usually require numerical analysis. Since significant supercooling occurs in our case, bubble collisions become the dominant source of GWs, whilst the interactions of bubbles with the surrounding plasma can be neglected. We then rely on previous numerical simulations employing the so-called envelope approximation [2, 3, 4, 5, 36, 37]. It is important to notice that these simulations have been performed under the assumption of a rapid PT and it is not clear a priori whether they are applicable to longer transitions. However, we take a special care to include the effects of the expansion of the Universe on the related growth of the vacuum bubbles and to study the behaviour of the nucleation probability at low temperature. Together with the previous simulations, this approach is expected to yield good estimates of the GW characteristics. Our calculations show that frequencies of the GWs produced by supercooled phase transitions are in the range $10^{-9}-10^{-7}$ Hz and allows for detection by pulsar timing array (PTA) experiments111Note also that if the QCD transition were first-order, it could be probed by PTA detectors as well [38, 39]. , such as the future Square Kilometre Array (SKA) [40].

This article is organised as follows. In Sec. 2, we briefly present the model based on a non-linear realisation of the electroweak gauge group. In Sec. 3, we describe the dynamics of a supercooled and long-lasting phase transition. Then we apply this formalism to the aforementioned model. In Sec. 4, we estimate the GW spectrum produced by this phase transition and compare the results with current and future PTA detectors. Finally, we discuss our results and approximations in the conclusion.

2 A non-linearly realised electroweak gauge group

It has been recently shown that a model with a non-linearly realised electroweak gauge group can accommodate a first-order phase transition [26]. We emphasise that the aim of this article is not to focus on the specifics of this model. Rather, such a “toy-model” allows us to illustrate a realistic mechanism that is capable of producing GWs associated with a supercooled and long-lasting PT. Henceforth, a brief description of the key features of this model is given below (cf. e.g. [27, 28, 29] for further details).

In this approach, the coset group $\mathcal{G}_{\text{coset}}=SU(2)_{L}\times U(1)_{Y}/U(1)_{Q}$ is gauged and the Higgs boson appears as a singlet $\rho(x)\sim(\mathbf{1},\mathbf{1})_{0}$ under the SM gauge group. The SM model Higgs doublet can then be identified as:

[TABLE]

where the three would-be Goldstone bosons spanning the coset space $\mathcal{G}_{\text{coset}}$ are represented by the $\pi^{i}(x)$ fields. The broken generators associated with this coset group are correspondingly $T^{i}=\sigma^{i}-\delta^{i3}\mathbb{I}$ , with $\sigma^{i}$ denoting the Pauli matrices. The physical Higgs $h$ is then identified as the fluctuation of $\rho$ around the electroweak vacuum expectation value $v=246$ GeV such that $\rho=v+h$ .

Following [26], all the SM configurations are assumed, with the exception of the Higgs potential. Indeed, as $\rho$ is a singlet, an anomalous cubic term is allowed in the following way:

[TABLE]

This tree-level potential explicitly depends on the three parameters $\mu$ , $\kappa$ and $\lambda$ . However, the relations $\left.\frac{dV}{d\rho}\right|_{\rho=v}=0$ and $\left.\frac{d^{2}V}{d\rho^{2}}\right|_{\rho=v}=m_{h}^{2}\approx(125\text{ GeV})^{2}$ allow the model to be controlled by a single free parameter, which is chosen to be $\kappa$ . Taking the example at tree level, the above relations can be solved analytically giving:

[TABLE]

The same process can be used to express $\mu$ and $\lambda$ as a function of $\kappa$ consistently at each order of perturbation, at least numerically. In this article, we solve the relations at one-loop level.

In order to describe the behaviour of the Higgs field in the early Universe, we require the one-loop finite temperature potential. It is usually written as follows [41, 42, 43, 44]:

[TABLE]

where $V^{(0)}$ is the classical potential (2.2), $V_{CW}^{(1)}$ is the one-loop Coleman-Weinberg potential at $T=0$ , $V^{(1)}(\rho,T)$ is the finite temperature contribution and $V_{Daisy}$ are correction terms dealing with infrared divergences. The explicit expressions of these contributions are given in the appendix A.1. For each value of $\kappa$ , the potential (2.4) can be numerically computed and the thermal behaviour of the Higgs field can be analysed.

3 Prolonged electroweak phase transition

The cosmological behaviour of the Higgs field is mainly described by its free energy density $\mathcal{F}(\rho,T)=V(\rho,T)$ identified as the effective potential (2.4). We can summarise its dynamics as a first-order PT using a few key temperatures ( $\tilde{T}>T_{c}>T_{n}>T_{p}$ ) as follows. For $T>\tilde{T}$ , $\mathcal{F}(\rho,T)$ admits a single minimum at $\rho=v_{T}^{(+)}$ called the symmetric phase of the Higgs field. As the Universe cools and reaches $\tilde{T}$ , a second minimum, called the broken phase, forms at $\rho=v_{T}^{(-)}$ with a free energy density initially higher than that of the symmetric phase. This free energy density then decreases until the two vacua become degenerate at the critical temperature $T=T_{c}$ . For $T<T_{c}$ , the free energy density of the broken phase keeps decreasing, causing the symmetric phase to become metastable: in other words, the Higgs field may tunnel through the potential barrier between $v_{T}^{(+)}$ and $v_{T}^{(-)}$ . If the decay probability is high enough, bubbles of true vacuum nucleate and expand in the surrounding symmetric phase. The nucleation temperature, $T_{n}$ is then defined as the temperature at which most of the bubbles are produced. On the other hand, the percolation temperature, $T_{p}$ corresponds to the instant when a significant volume of the Universe (whose value would be specified later on) has been converted from the symmetric to the broken phase. We expect most of bubble-collisions to occur around $T_{p}$ and not $T_{n}$ , unlike typical high-temperature (short-lived) phase transitions.

We should mention that some models admit a temperature $T_{0}<T_{c}$ below which the barrier between the two vacua disappears. If a transition has not occurred by this time, the Higgs field will then roll down the potential without forming any bubbles. This situation does not occur in our model of interest since the potential (2.2) admits a cubic term at zero temperature. In other words, the barrier will never vanish and a first-order PT can occur a priori at arbitrarily low $T$ , unless the Higgs field stays trapped in its metastable state. Although an electroweak PT is usually assumed to be quick, with $T_{p}\sim T_{n}$ , we can in this case consider a longer transition with $T_{p}\ll T_{n}$ and a large amount of supercooling. We now show how we can explicitly compute these temperatures from both the nucleation probability and the bubble dynamics.

3.1 Decay probability

The tunnelling of the Higgs field between the two vacua is characterised by the decay probability $\Gamma$ per unit time per unit volume. Quantum fluctuations drive this process at zero temperature [45, 46] while thermal fluctuations dominate at finite $T$ [47]. Therefore, $\Gamma$ is expressed as a function of the temperature of the Universe and can be written in the semiclassical approximation as follows:

[TABLE]

where $A(T)$ is a prefactor of mass dimension $4$ and $S(T)$ is the Euclidean action $S[\rho,T]$ evaluated along the bounce trajectory $\rho_{B}$ . In full generality, the Euclidean action is the functional over the Higgs field $\rho$ defined as [48]:

[TABLE]

where $\tau=-it$ is the Euclidean time, $\beta=\frac{1}{T}$ and $\tilde{\mathcal{F}}(\rho,T):=V(\rho,T)-V\Big{(}v_{T}^{(+)},T\Big{)}$ is the free energy density normalised according to its value in the unbroken phase. The bounce trajectory $\rho_{B}(\tau,r)$ is the solution which minimises the Euclidean action and thus satisfies the following equation of motion:

[TABLE]

with the boundary conditions

[TABLE]

The specific shape of the bounce $\rho(\tau,r)$ depends on the temperature [47]. At zero or low temperature, it reduces to an $O(4)$ symmetric solution $\rho(\tilde{r})$ with $\tilde{r}=\sqrt{\tau^{2}+r^{2}}$ , while at high temperature it is given by an $O(3)$ -symmetric and time-independent solution $\rho(r)$ . The temperature scale that allows us to distinguish between these regimes is given by the mass scale of the problem or equivalently by the size $R_{0}$ of the $O(4)$ symmetric bubble at $T=0$ . In both the limits $T\ll R_{0}^{-1}$ and $T\gg R_{0}^{-1}$ , the action (3.2) simplifies as follows:

[TABLE]

In these limits, the equations of motion for the bounce become:

[TABLE]

with $\alpha=2$ for $T\gg R_{0}^{-1}$ and $\alpha=3$ with $r$ replaced by $\tilde{r}$ for $T\ll R_{0}^{-1}$ . The prefactor $A(T)$ in equation (3.1) also admit different forms in the low and high temperature limits:

[TABLE]

The difference in these expressions comes from the fact that the $O(4)$ -symmetric bounce has 4 zero-modes contributing a factor $\left[S/(2\pi)\right]^{1/2}$ each, while the $O(3)$ -symmetric solution only has 3 zero-modes.

In the case of a rapid phase transition occurring around the electroweak scale $T_{EW}\sim 100$ GeV, the high-temperature formula provides a good approximation. However, it is not clear a priori how $T_{p}$ and $R_{0}$ will scale if the transition occurs with a significant amount of supercooling. In particular if $T_{p}\lesssim R_{0}^{-1}$ , approximating $S$ by $S_{3}/T$ might not be accurate anymore, requiring the use of the exact expression (3.2) (or $S_{4}$ at even lower temperature). For this reason, we compare how each of the three different actions $S(T)$ , $S_{4}(T)$ and $S_{3}(T)$ behaves as a function of the temperature. To do this, the bounce equations of motion must be solved numerically. In the low and high temperature regime, (3.6) is an ODE and can be integrated using the shooting method222Note that there is no bounce solution when the two vacua are exactly degenerate at $T_{c}$ and that tunnelling occurs only for $T<T_{c}$ . Numerically, the shooting method provides solutions only for a wide enough energy separation between the vacua, namely for $T\leq T_{\star}<T_{c}$ . Although tunnelling solutions can exist for $T_{\star}<T<T_{c}$ and be estimated through the thin-wall approximation, they are negligible for the PT as $\Gamma$ is more and more suppressed as the vacua are more and more degenerate (see e.g. Sec. IV in [45]). . On the other hand, the spacetime dependent equation (3.3) is a PDE and thus more difficult to address. Following [49, 48], we discretise the spacetime over a lattice. The PDE and the boundary conditions reduce then to a set of non-linear algebraic equations located at each point of the lattice. This set of equations is solved according to the Newton’s method: starting from a guess solution we build a new solution which minimises the error and iterate until the error becomes small enough. For this method to converge, the choice of the guess is important. In our case, we use the zero-temperature $O(4)$ solution, found from the shooting method, as a guess to solve (3.3) at $T=0+\Delta T$ . This solution is then used to solve the problem at $T^{(n+1)}=T^{(n)}+\Delta T$ recursively. The numerical solutions will be presented in Sec. 3.4.

3.2 Bubble dynamics and energy

Given the nucleation probability $\Gamma(T)$ discussed in the previous section, we can now describe the dynamics of a first-order PT. We apply the general formalism provided in [50]. We consider an expanding Universe with scale factor $a(t)$ and Hubble rate $H=\dot{a}/a$ . The probability $p(t)$ for a given point of spacetime to be in the symmetric phase at time $t$ is then given by [50]:

[TABLE]

where $\mathcal{I}(t)$ corresponds to the volume occupied by the true vacuum bubbles333Note that the exponentiation of $\mathcal{I}$ in (3.8) corrects the fact that regions with overlapping bubbles have been counted twice in $\mathcal{I}$ .. Indeed, bubbles which have nucleated at $t^{\prime}<t$ with probability $\Gamma(t^{\prime})$ would have then grown until $t$ reaching a (coordinate) radius $r(t,t^{\prime})$ given by:

[TABLE]

with $v(t)$ being the bubble wall velocity. In the previous equation, we have neglected the initial radius of the bubble which rapidly becomes negligible compared to the expanding size.

The condition that the phase transition completes can be translated to the condition that $p(t)\to 0$ for $t>t_{c}$ . As we are ultimately interested in the production of gravitational waves from bubble collisions, we are looking for the transition time corresponding to the period of maximum bubble collisions. This period can be estimated by the percolation time $t_{p}$ [15, 16]. According to numerical simulations performed with spheres of equal size, percolation occurs when approximately $29\%$ of space is covered by bubbles [51]. As suggested by [15, 16], we thus define $t_{p}$ from the condition $p(t_{p})\approx 0.7$ .

Knowing the collision time, we can then look for the distribution of number of bubbles at that time as a function of their size. From (3.9), a bubble formed at time $t_{R}$ will have a physical size $R(t,t_{R})=a(t)r(t,t_{R})$ at time $t$ . The number density of such bubbles is then given by [50]:

[TABLE]

For $t=t_{p}$ , the peak of this distribution gives us the size $\bar{R}$ of the majority of the bubbles which are colliding. Equivalently, it also provides the time $t_{\bar{R}}$ when most of these bubbles have been produced. We call this moment the nucleation time $t_{n}$ (rather than $t_{\bar{R}}$ ) and it can be explicitly computed via:

[TABLE]

As we shall see in Sec. 4, $\bar{R}:=R(t_{p},t_{n})$ is the key parameter to determine the peak frequency of the GW spectrum produced by bubble collisions.

Another important parameter, related to the amplitude of the GW spectrum, is the kinetic energy stored in the bubble walls. This kinetic energy comes from the vacuum energy released during the transition from the unbroken phase to the broken phase of the scalar field $\rho$ . In order to derive this quantity, we briefly remind the basic thermodynamic properties of this field. As described above, its free energy is given by the effective potential, $\mathcal{F}(\rho,T)=V(\rho,T)$ , and this allows us to define the pressure $p=-\mathcal{F}$ and the energy density $\epsilon(\rho,T)=\mathcal{F}-T\frac{d\mathcal{F}}{dT}$ . The released vacuum energy density is associated with the following latent heat: $\tilde{\epsilon}(T)=\epsilon(v_{T}^{(+)},T)-\epsilon(v_{T}^{(-)},T)$ . During the transition, this latent heat is converted into the formation of the bubbles (surface energy and kinetic energy of the walls) and into the reheating and fluid motion of the plasma. Following the notation of [12], we write $\kappa_{\rho}$ the fraction of energy which goes into the kinetic energy of the bubbles (i.e. the scalar field $\rho$ ). In our case, we can assume that $\kappa_{\rho}\sim 1$ as we are considering a very strong phase transition (see the discussion below for more details). As a result, the kinetic energy of a bubble is given by $\tilde{\epsilon}$ and the portion of space it has converted. For bubbles produced at $t_{n}$ , their kinetic energy at the percolation time $t_{p}$ is then

[TABLE]

where we have taken into account the fact that the latent heat varies with time (namely temperature). In the case of a short phase transition or a slowly varying $\tilde{\epsilon}$ , the above equation reduces to $E_{\text{kin}}=\frac{4\pi}{3}\bar{R}^{3}\tilde{\epsilon}$ as we should expect.

In order to explicitly compute $\bar{R}$ and $E_{\text{kin}}$ , we need to determine the bubble growth which depends on the velocity $v(t)$ and the scale factor $a(t)$ according to Eq. (3.9). We discuss the details of the evolution of the background Universe in the next section. Regarding the velocity, it is usually a difficult task to calculate precisely $v(t)$ as its depends on the interaction between the bubble wall and the plasma. However, it has been shown that for phase transitions with a sufficient amount of supercooling the produced bubbles quickly reach the speed of light [12]. They are referred to as runaway bubbles. Indeed, the amount of converted vacuum energy is such that the energy deposited in the plasma saturates and the majority goes into accelerating the bubble wall. This also confirms the previous assumption that $\kappa_{\rho}\sim 1$ . Note that this statement has been rigorously verified in [26] (see their Sec. 4) for our model of interest given in Sec. 2.

3.3 Equation of state

In order to carefully describe the dynamics of a long-lasting phase transition, the expansion of the Universe cannot be neglected and this requires to determine the scale factor $a(t)$ . In the same way, it is also important to know how the temperature of the Universe, $T(t)$ , evolves during the process. Both these quantities depend on the equation of state (EOS) of the different components of the Universe and which of them dominate. In particular, if the Universe is dominated by a single component with EOS $p=w\epsilon$ , the scale factor is then given by $a(t)\propto t^{\gamma}$ with $\gamma=\frac{2}{3(w+1)}$ ( $w\neq-1$ ). For $w<-1/3$ ( $\gamma>1$ ), it follows that the Universe undergoes an accelerated expansion (power-law inflation). In the same way, the case $w=-1$ (vacuum domination) also leads to an accelerating phase with $a(t)\propto e^{Ht}$ (exponential inflation).

In the general scenario of electroweak PT, bubbles nucleate near the electroweak scale, $T_{EW}\sim 100$ GeV, and percolate rapidly. During such a process, the Universe is radiation dominated with

[TABLE]

with $M_{p}\sim 1.22\times 10^{19}$ GeV the Planck mass and $g_{\star}\sim 100$ the effective number of relativistic degrees of freedom in the symmetric phase. However, this equation of state might not be valid in the case of strong supercooling or for a prolonged transition. The reason comes from the fact that as the Universe cools down, the vacuum energy density of the scalar field which remains in the unbroken phase starts dominating over the radiation energy density, $\epsilon_{\text{rad}}=\pi^{2}g_{\star}T^{4}/30$ , possibly leading to the phase of inflation described above.

In order to have a general understanding of transitions which such a behaviour, we introduce the time $t_{e}$ of radiation-vacuum equality satisfying $\epsilon_{\text{vac}}(t_{e})=\epsilon_{\text{rad}}(t_{e})$ .444For simplicity, we assume here that the Universe is dominated by a single component at a time and that the transition is sharp between radiation and vacuum domination. In the standard case, $t_{n}\lesssim t_{p}\ll t_{e}$ , and the vacuum energy is released into the bubbles before it could dominate. Then, scenarios with an inflationary background have been considered for some classes of scale-invariant models in [30, 31, 32]. In such cases, most of the bubbles nucleate after radiation-vacuum equality, namely $t_{e}<t_{n}\lesssim t_{p}$ . On the other hand, the process we want to describe in this article (prolonged PT) is different from the two previous ones in the sense that $t_{n}<t_{e}<t_{p}$ , namely bubbles are produced before vacuum energy would dominate and percolation requires a long time to complete. We now address this type of transition in more details.

The large separation between nucleation and percolation comes from a decay probability $\Gamma$ weaker than in the standard case, such that less bubbles are produced per unit volume and more time is required for them to collide. Let us clarify this reasoning by assuming that all bubbles are nucleated at $t_{n}$ such that $\Gamma(t)=\bar{\Gamma}(t_{n})\delta(t-t_{n})$ . The exponent in Eq. (3.8) becomes $\mathcal{I}(t)=\frac{4\pi}{3}\bar{\Gamma}(t_{n})a^{3}(t_{n})r^{3}(t,t_{n})$ and this clearly shows how a larger radius (i.e. longer time) compensates for a weaker nucleation probability. However, this last expression also indicates than the transition might never complete if the Universe is in accelerating expansion, since in such a case $r(t,t_{n})$ is bounded555This can easily be seen explicitly. Assuming $a(t)\propto t^{1/2}$ for $t<t_{e}$ and $a(t)\propto t^{\gamma}$ for $t>t_{e}$ , Eq. (3.9) (with $v\sim 1$ ) gives $r(t,t_{n})\propto 2(t_{e}^{1/2}-t_{n}^{1/2})+t_{e}^{\gamma-1/2}(t^{1-\gamma}-t_{e}^{1-\gamma})/(1-\gamma)$ , such that when $t\to\infty$ $r(t,t_{n})\to+\infty$ if $\gamma<1$ and $r(t,t_{n})\to(2+1/(\gamma-1))t_{e}^{1/2}-2t_{n}^{1/2}$ if $\gamma>1$ . when $t\to\infty$ . In other words, there is the possibility for the bubbles to not grow fast enough in order to reach each other and to collide.

However, we now argue that such a scenario (with no percolation) is unlikely to occur as long as enough bubbles are produced during the radiation dominated period, namely before $t_{e}$ . Indeed, as bubbles nucleate, vacuum energy is converted into kinetic energy of the wall motion such that the energy budget at the time $t_{e}$ is not simply dominated by vacuum energy density even if the bubbles have not yet collided. Actually, the bubbles are acting as inhomogeneity in the background of the expanding space-time and this render difficult to naively estimate what would be the corresponding dynamics of the Universe. According to several studies including numerical simulations [52, 53, 54] (see [55] for a recent review), it has been shown that small-field inflation is very unlikely to proceed with inhomogeneous initial conditions. We shall then assume in the following part of the article that for a sufficient number of bubbles produced at $t_{n}$ , the Universe expansion will not accelerate around $t_{e}$ and that percolation does occur at a given time $t_{p}>t_{e}$ .

An exact description of the evolution (namely a precise value of $\gamma$ ) would require numerical simulations which are beyond the scope of this study. As we expect no acceleration because of the previous argument, we have $\gamma<1$ and so we assume for simplicity that the Universe remains radiation dominated during the entire process ( $\gamma=1/2$ ). Deviation of the value of $\gamma$ in the range $[0,1]$ would change the estimation of the parameters describing the transition, in particular $\bar{R}$ and $E_{\text{kin}}$ , but not the qualitative picture. Moreover, we expect such a deviation to be compensated by a shift in the value of the initial conditions describing the underlying particle physics model (the parameter $\kappa$ in Eq. (2.2) in our case).

Under the aforementioned assumptions ( $\gamma\sim 1/2$ , $v\sim 1$ , $\kappa_{\rho}\sim 1$ ) and using Eq. (3.13), we can simplify Eqs. (3.8), (3.9) and (3.12) and write them in terms of temperature rather than time. Regarding the evolution of the temperature, we also recall that for a strong transition the dominant part of the vacuum energy is transformed into kinetic energy of the bubble walls meaning that we can neglect heating of the plasma. In the same way, the kinetic energy is subsequently transformed into GW energy through bubble collisions such that again heating is negligible. We eventually obtain the following key equations:

[TABLE]

It is now possible to numerically evaluate the previous expressions and to derive the key parameters $T_{n}$ , $T_{p}$ , $\bar{R}$ and $E_{\text{kin}}$ defining the phase transition. We present the results in the next section.

It is worth mentioning how the above formalism simplifies in the case of a quick phase transition, which is the main situation investigated in the literature. In that case, the PT is assumed to proceed rapidly around the temperature $\tilde{T}_{n}$ when at least one bubble has been produced per Hubble volume, namely $\int_{\tilde{T}_{n}}^{T_{c}}dT\frac{\Gamma(T)}{H^{4}(T)T}\sim 1$ . In this context, $\tilde{T}_{n}$ is called the nucleation temperature and replaces our expression $T_{n}$ derived from Eq. (3.11). Then the decay probability can be expanded around that instant as $\Gamma(t)\approx\Gamma(\tilde{t}_{n})e^{\beta(t-\tilde{t}_{n})}$ , where $\beta^{-1}$ gives the time scale of the transition. As such a PT is not expected to proceed too far below the electroweak scale, we have $\Gamma(T)\approx A(T)e^{-S_{3}(T)/T}$ and hence:

[TABLE]

The characteristic size and energy of the bubbles are then expected to be $\tilde{R}=v\beta^{-1}$ and $\tilde{E}_{\text{kin}}=\frac{4\pi}{3}\tilde{R}^{3}\epsilon(\tilde{T}_{n})$ respectively.

3.4 Numerical solutions

We give the numerical results of the previous formalism applied to the model described by the potential (2.2) and (2.4). The range of the parameter $\kappa\in[-1.85,-1]\frac{m_{h}^{2}}{v}$ has already been investigated in [26]. In that case, the phase transition occurs quickly and can be described by a rapid phase transition (as explained in the last paragraph of the previous section). It results in the production of a GW spectrum potentially detectable by eLISA. However, for $\kappa<-1.85\frac{m_{h}^{2}}{v}$ , the transition lasts longer and we require the more general prescription outlined in this paper.

First, we show in Fig. 1 how the Euclidean action behaves as a function of temperature, for $\kappa=-1.9\frac{m_{h}^{2}}{v}$ . It appears that the action $S(T)$ given by Eq. (3.2) is not only well approximated by $S_{4}(T)$ at $T\ll R_{0}^{-1}\approx 6$ GeV and $S_{3}(T)/T$ at $T\gg R_{0}^{-1}$ , but is also close to $R_{0}^{-1}$ . In other words, the $\min\{S_{4},S_{3}/T\}$ provides a good approximation of the action over the entire range of temperatures considered, as also suggested in [47]. This observation is important from the computational point of view as it means we can avoid solving the time-dependent PDE (3.3) which is computationally expensive. Moreover, we observe that the action becomes, and stays, large at low temperatures ( $S\sim S_{4}\sim 930$ ), meaning that its effect on the PT dynamics is exponentially suppressed (see for example (3.1)).

Second, $p(T)$ is computed from $\Gamma(T)$ according to Eq. (3.14). The results for several key values of $\kappa$ are given in Fig. 2. As expected, the phase transition can be identified as a rapid change in $p(T)$ from $1$ to [math]. The corresponding nucleation and percolation temperatures are given in Table 1. We observe that $T_{n}\sim 49$ GeV for each $\kappa$ . This is due to the fact that most of the bubbles are produced when the action $S_{3}(T)/T$ reaches its minimum, whose location only slightly changes with $\kappa$ . On the other hand, we notice that $T_{p}$ varies through several order of magnitudes. This is because the number density of bubbles produced at the nucleation time changes as a function of $\kappa$ . This confirm the expectation that as the decay probability decreases, more time is needed for the transition to complete. We have verified that these results are consistent with our assumption that most bubbles are produced before vacuum energy dominates. Indeed, vacuum-radiation equality would occur at $T_{e}\sim 35.5$ GeV for this range of $\kappa$ , if no bubbles were produced earlier. This confirms $T_{n}>T_{e}$ .

The remaining task is to compute the characteristic bubble size $\bar{R}$ and the kinetic energy $E_{\text{kin}}$ of the bubbles at the percolation temperature. For later convenience, we rescale them compared to the Hubble radius and radiation energy density at this time. To this end, Table 1 uses the dimensionless parameters $(\bar{R}H_{p})^{-1}$ and $\alpha=\epsilon_{\text{kin}}/\epsilon_{\text{rad}}(T_{p})$ with $\epsilon_{\text{kin}}=E_{\text{kin}}/\bar{R}^{3}$ . In this way, one can easily compare our results with those mentioned in the literature (as $v\sim 1$ , $(\bar{R}H_{p})^{-1}$ takes the role of $\beta/H$ given by Eq. (3.15)).666Note that Eq. (3.15) could clearly not have been applied in this scenario as $\beta$ would have been negative for temperature below the minimum of $S_{3}/T$ .

We confirm from Table 1, that there exists some parameters for which the PT completes well below the electroweak scale. Lower values of percolation temperature are also possible for lower values of $\kappa$ , but they are constrained from various considerations. First, the obvious lower bound $T_{p}\gtrsim 1$ MeV is given by nucleosynthesis constraints. Second, we shall see in the next section that the GW spectrum for $\kappa\lesssim-1.92\frac{m_{h}^{2}}{v}$ is already excluded from current PTA surveys.

Note also that at temperatures around $T\sim 100$ MeV, the QCD phase transition, which is believed to be second-order, takes place. As discussed in [33], a $q\bar{q}$ condensate with non zero vacuum expectation value is thus expected to form and to contribute to the thermal potential (2.4) via a linear term in the Higgs field. This effect has no significant influence on our model since the cubic term in the tree-level potential (2.2) induces a large barrier which is not affected by such a linear term. It is however of importance for models in which the barrier becomes weaker at low temperature, see e.g. [33, 56, 34].

Another observation from Table 1 is that the two quantities $(\bar{R}H_{p})^{-1}$ and $\frac{\alpha}{\alpha+1}$ (which will be important in the next section) approach $1$ for lower and lower values of $\kappa$ . The fact that $(\bar{R}H_{p})^{-1}\to 1$ means that bubbles are almost of horizon size when they collide and that they never become of super-horizon size because $H^{-1}$ also increases linearly with time. Note also that the increase in $\alpha$ for lower $\kappa$ is mainly due to $\alpha\propto T_{p}^{-4}$ rather than to the change of kinetic energy stored in the bubbles.

4 Gravitational wave production

A stochastic background of gravitational waves is usually described in terms of its contribution to the energy density of the Universe per frequency interval:

[TABLE]

where $f$ is the frequency, $\epsilon_{GW}$ the gravitational wave energy density, and $\epsilon_{c}=3H_{0}^{2}/(8\pi G)$ is the critical energy density today. The production of GWs from a first-order phase transition originates from three sources: the collisions of bubbles walls [1, 2, 3, 4, 5, 37, 36], sound waves in the plasma formed after collision [57, 58, 59, 60] and magnetohydrodynamics turbulences in the plasma [61, 62, 63, 64, 65]. As these three contributions should approximately linearly combine [12], the total energy density can be written as

[TABLE]

For models with a significant amount of supercooling (with $v\sim 1$ and $\kappa_{\rho}\sim 1$ as described previously), the collision term $\Omega_{col}$ is dominant [12]. In this article, we can then assume that:

[TABLE]

We can thus expect the GW spectrum to be mainly produced around the percolation temperature $t_{p}$ of the phase transition with a characteristic frequency $f_{p}$ . The amplitude of this signal then decreases as $a^{-4}(t)$ up to today while its frequency redshifts as $a^{-1}(t)$ . In other words, the energy density stored in the GWs and the peak frequency today are given by [5]:

[TABLE]

We can estimate the properties of the above GW spectrum from dimensional analysis. We expect the peak frequency to scale with the inverse size of the bubbles at their collision, namely $f_{p}\sim(\bar{R})^{-1}$ . Regarding the amplitude, we first have $\Omega_{col}=\epsilon_{GW}/\epsilon_{tot}$ where $\epsilon_{tot}$ is the total energy density of the Universe at the percolation time. The energy density of gravitational waves is then given by [66]: $\epsilon_{GW}=\frac{1}{8\pi G}\langle\partial_{t}h_{\mu\nu}\partial_{t}h^{\mu\nu}\rangle$ , where $h_{\mu\nu}$ is the metric perturbation. It satisfies Einstein’s equations $\partial_{\alpha}\partial^{\alpha}h_{\mu\nu}\sim 8\pi GT_{\mu\nu}$ with $T_{\mu\nu}$ the energy momentum tensor describing the source. In our case, we expect that $\partial_{\alpha}\sim\bar{R}^{-1}$ (the characteristic size of bubbles) and $T_{\mu\nu}\sim\epsilon_{\text{kin}}$ (the kinetic energy density stored in the bubbles). This implies $\epsilon_{GW}\sim 8\pi G\bar{R}^{2}\epsilon^{2}_{\text{kin}}$ . Substituting $G$ from the Friedmann equation $H_{p}^{2}=\frac{8\pi G}{3}\epsilon_{tot}$ , we get:

[TABLE]

where in the last equality we have used the fact that $\epsilon_{tot}=\epsilon_{\text{kin}}+\epsilon_{\text{rad}}$ . As expected, the GW amplitude can then be estimated from the two parameters $\bar{R}H_{p}$ and $\alpha$ which have been computed in the previous section and are given in Table 1. Note that in general Eq. (4.5) would also depend on $v$ and $\kappa_{\rho}$ , which again have been assumed to be close to unity in our case.

Several studies provide more accurate expressions for the GW spectrum from bubble collisions, beyond the simple dimensional analysis. They usually rely on the envelope approximation and numerical simulations [4, 5, 37], although some analytical formula have also been suggested [36, 67]. Using the results of [37], the spectrum today can be described as follows:

[TABLE]

where:

[TABLE]

It is important to realise that these formula have been derived under the assumption of a short-lasting PT as described at the end of Sec. 3.2. This is why they use the coefficient $\beta$ , defined by Eq. (3.15), which corresponds to the time scale of the transition. In this article, we shall instead substitute $\beta\sim v\bar{R}^{-1}\sim\bar{R}^{-1}$ in Eq. (4.6-4.7) in agreement with the dimensional estimate given by Eq. (4.5). However, it is not clear if these fitted formula accurately describe a long-lasting transition. Obtaining a more precise result would require to further numerical simulations of bubble collisions without the assumption that the PT completes in a time less than the Hubble time. Such a task is beyond the scope of this article. Indeed, our main aim is to find in which broad range of frequency and amplitude lies the stochastic GW background produced by a supercooled and long-lasting transition and if this significantly deviates from the usual case. If so, a more precise description of the GW spectrum would be interesting to compute in future works.

It is now straightforward to evaluate Eqs. (4.6)-(4.7) for the parameters given in Table 1, describing the phase transition occurring in our model of interest. It appears that the peak frequency $f_{0}$ can be as low as $\sim 10^{-9}-10^{-7}$ Hz and thus lies in the detection range of PTA experiments. In Fig. 3, we compare the GW spectrum for several values of $\kappa$ and the current status of PTA detectors. Three collaborations have published limits on the amplitude of a stochastic GW background: the European Pulsar Timing Array (EPTA) [68], the Parkes Pulsar Timing Array (PPTA) [69] and the North American Nanohertz Observatory for Gravitational waves (NANOGrav) [70]. All these three limits are of similar amplitudes and thus we only display the EPTA results777Note that the exclusion line is computed as explained in [71].. The sensitivity area should be improved in the future by the Square Kilometre Array (SKA) [40] whose expected detection range is also given in Fig. 3 [71].

It is clear that our investigated model predicts GWs detectable by the aforementioned detectors and consequently we argue that this demonstrates a new method for probing first-order electroweak PTs. We might be tempted now to define the range of $\kappa$ values which is not yet excluded by EPTA and which could be tested by SKA. Although Fig. 3 gives us a rough estimate around $\kappa\sim-1.9\frac{m_{h}^{2}}{v}$ , this task would require a more precise computation of the GW spectrum taking into account the long duration of the transition in the numerical simulations. We also remind that these results rely on our assumption of a radiation dominated period. As explained earlier, the exact equation of state describing the Universe during the transition is expected to be more complicated with the potential effect of changing the value of the parameters entering in the GW spectrum (but not the general behaviour). However, we emphasize that the model we discussed in Sec. 2 was mainly introduced as a case of study to illustrate the dynamics of a supercooled and long-lasting PT. We expect other theories to present similar features, especially a class of scale-invariant models with very light scalar particles [33, 35, 34].

5 Conclusion

In this article, we investigated the production of gravitational waves during a strongly supercooled electroweak phase transition. Considering a particle physics model based on a nonlinear realisation of the electroweak gauge group, we carefully computed the dynamics of the Higgs field during such a transition. In particular, we found an interesting range of parameters for which the PT completes at a temperature significantly below the electroweak scale. This analysis required to take into account the expansion of the Universe and to study the behaviour of the nucleation probability of true vacuum bubbles at low temperature. Regarding the latter point, we compared the usual low and high temperature expressions of the Euclidean action to a more general formula based on a time-dependent bounce solution. We observed that the two approximate equations actually provide a good estimate of the action over the entire range of valid temperatures.

In this scenario, we argued that GWs are produced by the collisions of large bubbles (of the order the horizon size) with a kinetic energy density significantly higher than the radiation energy density of the Universe at the time of collision. This results in a large amplitude stochastic GW background in the frequency range $10^{-9}-10^{-7}$ Hz which can be probed by pulsar timing arrays. We derived this prediction from both dimensional arguments and the use of previous numerical simulations of colliding bubbles. Although it is clear from these analysis that our model of interest predicts a GW spectrum in the sensitivity band of PTA detectors, more refined simulations would be needed to improve the accuracy of our results. In particular, a better estimation of the exact equation of state and scale factor during the phase transition would increase the accuracy of the relation between GW predictions and specific values of the parameter $\kappa$ of the particle physics model we considered.

In summary, we showed that a sufficiently supercooled electroweak phase transition can be detected with pulsar timing arrays. This enlarges the way of probing first-order cosmological PT in addition to previous proposals with space-base interferometers such as eLISA. This also increases the prospects of new generation PTA detectors like SKA.

Acknowledgments

We would like to thank Csaba Balazs and Alexander Kusenko for useful discussions. This work was partially supported by the Australian Research Council.

Appendix A Appendix

A.1 Finite Temperature Potential

We give here the detailed expression of the finite temperature potential (2.4). The Coleman-Weinberg contribution $V_{CW}^{(1)}$ at $T=0$ is given by:

[TABLE]

The finite temperature part $V^{(1)}(\rho,T)$ is defined via the thermal function $J$ :

[TABLE]

where $s_{i}$ corresponds to the spin and $n_{i}$ to the number of degrees of freedom of the particle species $i$ . The field-dependent masses $m_{i}(\rho)$ are given by

[TABLE]

The correction from Daisy terms can be described by a shift in the longitudinal components of the respective boson masses by their Debye correction (cf. e.g. [44, 72]):

[TABLE]

where:

[TABLE]

The number of degrees of freedom is then:

[TABLE]

Bibliography72

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1[1] M. S. Turner and F. Wilczek, Relic gravitational waves and extended inflation , Phys. Rev. Lett. 65 (1990) 3080–3083 . · doi ↗
2[2] A. Kosowsky, M. S. Turner and R. Watkins, Gravitational radiation from colliding vacuum bubbles , Phys. Rev. D 45 (1992) 4514–4535 . · doi ↗
3[3] A. Kosowsky, M. S. Turner and R. Watkins, Gravitational waves from first order cosmological phase transitions , Phys. Rev. Lett. 69 (1992) 2026–2029 . · doi ↗
4[4] A. Kosowsky and M. S. Turner, Gravitational radiation from colliding vacuum bubbles: envelope approximation to many bubble collisions , Phys. Rev. D 47 (1993) 4372–4391 , [ astro-ph/9211004 ]. · doi ↗
5[5] M. Kamionkowski, A. Kosowsky and M. S. Turner, Gravitational radiation from first order phase transitions , Phys. Rev. D 49 (1994) 2837–2851 , [ astro-ph/9310044 ]. · doi ↗
6[6] K. Kajantie, M. Laine, K. Rummukainen and M. E. Shaposhnikov, Is there a hot electroweak phase transition at m(H) larger or equal to m(W)? , Phys. Rev. Lett. 77 (1996) 2887–2890 , [ hep-ph/9605288 ]. · doi ↗
7[7] V. A. Kuzmin, V. A. Rubakov and M. E. Shaposhnikov, On the Anomalous Electroweak Baryon Number Nonconservation in the Early Universe , Phys. Lett. B 155 (1985) 36 . · doi ↗
8[8] M. Trodden, Electroweak baryogenesis , Rev. Mod. Phys. 71 (1999) 1463–1500 , [ hep-ph/9803479 ]. · doi ↗