Probing Trans-Electroweak First Order Phase Transitions from Gravitational Waves

Abstract

We propose direct tests of very high energy first-order phase transitions, which are elusive to collider physics, deploying the gravitational waves’ measurements. We show that first-order phase transitions lying in a large window of critical temperatures, which is considerably larger than the electroweak energy scale, can be tested from advanced LIGO (aLIGO) and the Einstein Telescope. This provides the possibility to probe several inflationary mechanisms ending with the inflaton in a false minimum and high-energy first order phase transitions that are due to new scalar bosons, beyond the Standard Model of particle physics. As an important example, we consider the axion monodromy inflationary scenario and analyze the potential for its experimental verification, deploying the gravitational wave interferometers.

1. Introduction

The first direct detections of gravitational waves (GW) from the merging of black holes and neutron stars recently measured by the LIGO/VIRGO collaboration raised an urgent question: can we detect any hint of new physics beyond the Standard Model (SM) of particles and cosmology within the same frequency sensitivity range of LIGO/VIRGO?

In this letter, we show that first order phase transitions (FOPT) at very high energy, beyond the electroweak scale, can be tested in the next generation of experiments beyond LIGO and VIRGO, including the aLIGO and Einstein Telescope (ET). According to the standard scenario of cosmological FOPTs, our Universe was initially set in a false vacuum state. After a certain time scale, a tunneling toward the true vacuum state happened. This induced the nucleation of bubbles, relativistically expanding with constant acceleration. The violent acceleration was sourced by a difference of pressure between the exterior false vacuum and the interior true vacuum. Such bubbles nucleating in the early Universe have produced the stochastic GW background through three distinct processes: bubble-bubble scattering, acoustic shock waves, and magnetohydrodynamic (MHD) turbulence.

From the experimental perspective, the temperature of the FOPT determines the characteristic frequency window of GWs. For instance, the phase transitions with critical temperatures around $100\mathrm{GeV}÷1\mathrm{TeV}$ produce the GW signals that are peaked within the frequency range of $1÷10\mathrm{MHz}$ [1,2,3,4,5,6,7,8]. Intriguingly, these frequencies can be probed by the next generation of interferometers, including LISA, BBO, and U-DECIGO [9,10,11] (see also a more recent analysis in [12,13,14,15,16,17] and [18] for a complete review on recent developments in the subject). Several theoretical models associated with an FOPT at the energy scale of about $100\mathrm{GeV}÷1\mathrm{TeV}$ were hitherto studied [9,10,11]. Nonetheless, FOPTs with a critical temperature much larger than the electroweak scale cannot be observed by LISA, U-DECIGO, and BBO. This seems to preclude a possibility of testing the very high-energy FOPTs due to new physics beyond the SM.

In this work, we show that FOPTs at energy scales much higher than the electroweak energy scale lie in the frequency window of LIGO/VIRGO, KAGRA, aLIGO, and ET. The FOPT scenarios can be very elusive when the current data from LIGO and VIRGO are considered. Specifically, aLIGO and ET may achieve very interesting sensitivity ranges, allowing testing FOPTs at energies much higher than the electroweak scale.

Although these considerations may have very strong implications for probing the physics of the early Universe, these were never discussed in the literature. By testing the cosmological FOPTs, we can infer valuable information on such stages of the Universe’s evolution as inflation, Grand Unified Theories (GUTs), baryogenesis, Primordial Black Holes’ (PBHs) production, etc. For instance, many models of inflation, with the inflationary regime ending in a false minimum after reheating, can be tested and eventually falsified [19,20,21,22]. Several GUT models, including the $SO\left(10\right)$ Pati–Salam, predict a high energy pattern of FOPTs sourced from a rich Higgs sector, for example as the spontaneous symmetry breaking pattern $SO\left(10\right)\to$ $SU\left(4\right)\times SU{\left(2\right)}_L\times SU{\left(2\right)}_R\to$ $S{U}_c\left(3\right)\times U_{B-L}\left(1\right)\times SU{\left(2\right)}_L\times SU{\left(2\right)}_R\to$ $S{U}_c\left(3\right)\times S{U}_L\left(2\right)\times U_Y\left(1\right)$ [23] realizes. On the other hand, PBHs may compose a part of dark matter and be efficiently produced from FOPTs, through bubble-bubble scatterings [24].

The large variety of possible instantiations of cosmological FOPTs at energies beyond the electroweak scale calls for a model-independent analysis, in the framework of effective field theories. We investigated a space of parameters that corresponds to a large class of models related to high-energy FOPTs. As an inverse scattering approach, our analysis can be useful in order to infer new bounds on the inflaton potentials, as well as to constrain extensions of the Higgs sector beyond the SM, deploying the GWs data from LIGO and VIRGO, and future data from other experiments, probing the same frequency window, but with more sensitivity in the GWs energy-density. Elusive to any possible laboratory experiment, this possibility is certainly intriguing and worth a thorough discussion.

2. GW Signals from FOPT

The FOPTs in the early Universe could produce vacuum bubbles that expanded at relativistic velocities. Such a mechanism produces a stochastic background of gravitational radiation. GWs are sourced via bubble-bubble collisions, sound waves, and turbulence generated by the bubbles’ expansion in the plasma.

The three contributions to the GW spectrum are related to the effective potential of the scalar field, $V_{\mathrm{eff}}\left(T\right)$, which is responsible for the bubbles’ coalescence. In the early Universe plasma, the scalar potential received extra contributions that depended on the plasma temperature and the coupling constants to the other fields.

The bubble nucleation has a rate that formally reads as $\mathrm{\Gamma }\left(t\right)=A\left(t\right)e^{-S\left(t\right)}$, where S is the Euclidean action describing the bubble dynamics. From this expression, one can define a critical temperature $T_{*}$ for the FOPT, at which the nucleation rate becomes large enough to induce at least the nucleation of one bubble per horizon volume with a probability of order one. From the $\mathrm{\Gamma }\left(t\right)$ rate, one therefore defines:

$$\begin{array}{c}\hfill \beta =-{\left[\frac{dS}{dt}\right]}_{t=t_{*}}\simeq \frac{\dot{\mathrm{\Gamma }}}{\mathrm{\Gamma }},\end{array}$$

(1)

as an inverse characteristic time-scale of the FOPT. The hierarchy among the $\beta$ and the Hubble rate at the FOPT temperature $H_{*}\equiv H\left(T_{*}\right)$ is crucial to ensure an efficient production of GWs. The ratio can be cast as:

$$\frac{\beta }{H_{*}}=T_{*}{\left[\frac{dS}{dT}\right]}_{T=T_{*}}.$$

Another important parameter is:

$$\begin{array}{c}\hfill \alpha =\frac{{\rho }_{\mathrm{vac}}}{{\rho }_{\mathrm{rad}}\left(T_{*}\right)},\end{array}$$

(2)

which measures the vacuum energy density released in the FOPT compared to the energy density of radiation at $T_{*}$; the latter is ${\rho }_{\mathrm{rad}}\left(T_{*}\right)=g_{*}{\pi }^2{T}_{*}^4/30$, with $g_{*}$ denoting the number of cosmological relativistic degrees of freedom at $T_{*}$. Other useful parameters characterizing the GW spectrum are:

$$\begin{array}{c}\hfill {\kappa }_{\varphi }=\frac{{\rho }_{\varphi }}{{\rho }_{\mathrm{vac}}},{\kappa }_v=\frac{{\rho }_v}{{\rho }_{\mathrm{vac}}},\end{array}$$

(3)

measuring the inverse fractions of the vacuum energy converted into gradient energy of the scalar field and into the fluid motion, respectively. Finally, a relevant parameter characterizing the GW spectrum is the bubble wall velocity, $v_w$, in the rest frame of the fluid that is found asymptotically far from the bubble.

It is worth noting that all these quantities have direct relations to the thermally-corrected scalar potential, appearing in the Euclidean action and in the released vacuum energy density, i.e.,

$$\begin{array}{c}{\rho }_{\mathrm{vac}}\left(T_{*}\right)=-\left[T\frac{d{V}_{\mathrm{eff}}}{dT}-V_{\mathrm{eff}}\left(T\right)\right]{|}_{T=T_{*}},\hfill \\ S\left(T\right)\simeq \frac{S_3\left(T\right)}{T},S_3\equiv \int d^{3}r[{\partial }_i{\sigma }^{†}{\partial }_i\sigma +V_{\mathrm{eff}}(\sigma ,T)],\hfill \\ {\mathrm{\Gamma }}_0\left(T\right)\sim T^4,\mathrm{\Gamma }={\mathrm{\Gamma }}_0\left(T\right)exp[-S\left(T\right)],\hfill \end{array}$$

where $\sigma$ is the scalar field undergoing the FOPT.

3. Bubble-Bubble Collision

The GW energy spectrum produced during the collision of two bubbles depends on the false vacuum energy and the bubble size at collision. This implies that the GW spectrum depends only on the ratio between the rate of bubbles’ nucleation and the Hubble rate of the Universe, at the temperature of the FOPT, as well as on the latent heat energy and the bubble wall velocity. The redshift of GWs introduces a dependence on the number of cosmological degrees of freedom. During the FOPT, the GW intensity peak produced by the bubbles’ collisions reads [4,5]:

$$\begin{array}{c}\hfill \frac{f_{*}}{\beta }=\left(\frac{0.62}{1.8-0.1v_w+v_w^2}\right),\end{array}$$

(4)

where $f_{*}$ is determined by the characteristic time-scale of the FOPT, i.e., $1/\beta$. Using the inverse Hubble time at the GW production, redshifted today, from the well-known relation:

$$h_{*}=16.5\times {10}^{-3}\mathrm{MHz}\left(\frac{f_{*}}{\beta }\right)\left(\frac{\beta }{H_{*}}\right)\left(\frac{T_{*}}{100\mathrm{GeV}}\right){\left(\frac{g_{*}}{100}\right)}^{\frac{1}{6}},$$

one obtains the frequency peak:

$$\begin{array}{c}{f}_{\mathrm{coll}}=16.5\times {10}^{-3}\mathrm{MHz}\left(\frac{T_{*}}{100\mathrm{GeV}}\right){\left(\frac{g_{*}}{100}\right)}^{\frac{1}{6}}.\hfill \end{array}$$

(5)

This estimate was obtained assuming that the FOPT had happened in a radiation-dominated epoch. However, in some particular instantiations beyond the SM, this assumption may be relaxed. For a number of degrees of freedom that is close to the SM one, the GWs’ frequency ranges in the domain $1÷100\mathrm{Hz}$, which can be achieved for a range of temperatures $T_{*}\simeq {10}^7÷{10}^9\mathrm{GeV}$. This motivates our proposal in terms of order of magnitude.

At the same time, one should achieve an energy-density of GWs that is high enough to reach the sensitivity curves of the next generation of experiments probing the $1÷100\mathrm{Hz}$ frequency range. The GW energy-density peak corresponds to [4,5]:

$$\begin{array}{ccc}\hfill h^2{\mathrm{\Omega }}_{\mathrm{coll}}\left(f\right)& =& 1.67\times {10}^{-5}{\left(\frac{H_{*}}{\beta }\right)}^2{\left(\frac{\kappa \alpha }{1+\alpha }\right)}^2\hfill \\ & \times & {\left(\frac{100}{g_{*}}\right)}^{\frac{1}{3}}\left(\frac{0.11v_w^3}{0.42+v_w^2}\right),\hfill \end{array}$$

(6)

where $\kappa$ is the fraction of the latent heat deposited, in the envelope approximation, on the front of the FOPT. Assuming $\kappa \sim 10\%$, $H_{*}/\beta =1/100$, $v_w\sim 1$, and $g_{*}\sim g_{*}^{SM}\sim 100$, we can immediately estimate that $h^2{\mathrm{\Omega }}_{coll}\left(f\right)$ reaches the range ${10}^{-9}÷{10}^{-10}$. Although very far from the LIGO/VIRGO current sensitivity, this signal can be reached by the ET. More favorable cases are the ones with $\beta /H_{*}=1÷10$, for which $h^2{\mathrm{\Omega }}_{coll}\left(f\right)$ falls in the range ${10}^{-5}÷{10}^{-7}$, all the other parameters being taken as above, and these can be then tested by aLIGO.

4. Sound Waves

The estimates of the GW peak from sound waves is affected by $O\left(1\right)$ numerical uncertainties [25,26,27]. The order of magnitude can be inferred anyway, considering the characteristic scale set by the average bubble separation, i.e., $R_{*}={\left(8\pi \right)}^{1/3}{v}_w/\beta$. From a conservative estimate, the frequency peak reads:

$$f_{\mathrm{sw}}=\frac{2}{\sqrt{3}}\frac{\beta }{v_w},$$

which, redshifted until the Universe today, reads:

$$\begin{array}{c}\hfill f_{\mathrm{sw}}=1.9\times {10}^{-2}\mathrm{MHz}\frac{1}{v_w}\left(\frac{\beta }{H_{*}}\right)\left(\frac{T_{*}}{100\mathrm{GeV}}\right){\left(\frac{g_{*}}{100}\right)}^{\frac{1}{6}}.\end{array}$$

(7)

Contrary to Equation (5), in the latter expression, there is also a dependence of the frequency peak on $\beta /H_{*}$. In the reasonable range $\beta /H_{*}=10÷1000$, the possible frequencies are $1÷100\mathrm{Hz}$. For $v_w\simeq 1$, ${\kappa }_v\sim 1÷10\%$, this latter corresponds to the range of temperatures $T_{*}={10}^6÷{10}^{10}\mathrm{GeV}$. Nonetheless, one should also consider that the energy-density scales are the inverse of $\beta /H_{*}$. In particular, the peak can be estimated to be:

$$\begin{array}{c}\hfill h^2{\mathrm{\Omega }}_{\mathrm{sw}}=2.65\times {10}^{-6}\left(\frac{H_{*}}{\beta }\right)\\ \hfill \times {\left(\frac{{\kappa }_v\alpha }{1+\alpha }\right)}^2{\left(\frac{100}{g_{*}}\right)}^{\frac{1}{3}}{v}_w.\end{array}$$

(8)

Detection from LIGO/VIRGO seems also in this case impossible, and very elusive even for aLIGO. Nonetheless, the ET can test this scenario for ${\kappa }_v\sim 10\%$, $v_w\sim 1$, $H_{*}/\beta \sim 1÷10$, $T_{*}\sim {10}^7÷{10}^8\mathrm{GeV}$ and $\alpha \sim 1$, since the corresponding amplitudes would peak in the range $h^2{\mathrm{\Omega }}_{sw}\sim {10}^{-7}÷{10}^{-8}$.

5. MHD Turbulence

Finally, the frequency peak for MHD turbulence was recovered to be [9]:

$$\begin{array}{c}\hfill f_{\mathrm{turb}}=2.7\times {10}^{-2}\mathrm{MHz}\frac{1}{v_w}\left(\frac{\beta }{H_{*}}\right)\left(\frac{T_{*}}{100\mathrm{GeV}}\right){\left(\frac{g_{*}}{100}\right)}^{\frac{1}{6}}.\end{array}$$

In the reasonable domain $\beta /H_{*}=10÷1000$, one obtains possible frequencies in the $1÷100\mathrm{Hz}$ range. Specifically, for $v_w\simeq 1$, the temperatures’ range spans the region $T_{*}={10}^4÷{10}^8\mathrm{GeV}$. The GW energy-density peak is expressed by:

$$\begin{array}{c}\hfill h^2{\mathrm{\Omega }}_{\mathrm{turb}}=3.35\times {10}^{-4}\left(\frac{H_{*}}{\beta }\right){\left(\frac{{\kappa }_t\alpha }{1+\alpha }\right)}^{\frac{3}{2}}{v}_w{\left(\frac{100}{g_{*}}\right)}^{\frac{1}{3}},\end{array}$$

(9)

where ${\kappa }_t$ denotes the fraction of latent heat transformed into turbulence of the plasma. Interesting situations similar to the ones already discussed for the case of sound waves can be recovered also for turbulence. Of course, approximations $O\left(1\right)$ are implied as well while considering the turbulence effects, and the competition of turbulence and sound waves is very much sensitive to the initial conditions. In both cases, these are issues not relevant to our discussion, since one cannot claim these predictions to be in the domain of high precision physics.

6. Dynamics of Bubbles

The GW signals from FOPTs can radically differ from one another, depending on different regimes considered, in which either bubbles’ collision, or MHD turbulence, or sound waves’ contributions may dominate. The dynamical regime is selected by the model-dependent component of the effective parameters introduced above. In our discussion, which is model independent, we will consider the most likely scenarios. These can be divided into two classes: (i) non-runaway (in the plasma) bubbles; (ii) runaway bubbles.

(i) Non-runaway bubbles in the plasma: Bubbles expanding in a cosmic plasma can reach a relativistic velocity. In this situation, the energy stored in the scalar field can be neglected, since it would only be scaled with the surface of the bubble despite its volume. The most relevant contributions to the signal are expected from the cosmic fluid motion induced from the bubbles’ expansion, i.e., from sound waves and/or MHD turbulence. The energy-density spectrum can be approximated to be:

$$h^2{\mathrm{\Omega }}_{GW}\simeq h^2{\mathrm{\Omega }}_{turb}+h^2{\mathrm{\Omega }}_{sw}.$$

(10)

For this estimate, it is relevant to include ${\kappa }_v$ as the efficiency factor of the conversion of the latent heat into the bubbles’ motion through the bulk. The efficiency factor of non-runaway bubbles can be modeled with a good approximation and is found to be in the two opposite velocity regimes:

$$\begin{array}{c}{\kappa }_v\simeq \alpha {(0.73+0.083\sqrt{\alpha }+\alpha )}^{-1},v_w\sim 1,\hfill \\ {\kappa }_v\simeq v_w^{6/5}6.9\alpha {(1.36-0.037\sqrt{\alpha }+\alpha )}^{-1},v_w\le 1.\hfill \end{array}$$

The efficiency of the energy transfer to the turbulence is then ${\kappa }_t\simeq ϵ{\kappa }_v$, such that $ϵ$ denotes the efficiency for the transfer of the bulk motion into the plasma turbulence; from numerical simulations, one finds that $ϵ\sim 5÷10\%$ [9,28]. In Figure 1, we display several simulations derived within this regime. We show that in some viable cases, sensitivity curves of the ET can be reached.

(ii) Runaway bubbles in the plasma: Another possibility is that the bubble wall experiences unbounded acceleration, running away with $v_w\to 1$. This case corresponds to a very high energy density stored in the scalar field profile, i.e., for $\alpha \gg 1$:

$$\begin{array}{c}\hfill h^2{\mathrm{\Omega }}_{\mathrm{GW}}\simeq h^2{\mathrm{\Omega }}_{\varphi }+h^2{\mathrm{\Omega }}_{\mathrm{sw}}+h^2{\mathrm{\Omega }}_{\mathrm{turb}},\end{array}$$

where ${\mathrm{\Omega }}_{\varphi }$ is sourced by the gradients of the scalar field, modeled in the envelope approximation. The variable ${\alpha }_{\mathrm{\infty }}$ is defined as the minimum threshold value of $\alpha$ such that the bubbles start to run away. The efficiency of the energy transfer is parametrized by:

$$\begin{array}{c}{\kappa }_{\varphi }=(\alpha -{\alpha }_{\mathrm{\infty }})/\alpha \ge 0,{\kappa }_v=\frac{{\alpha }_{\mathrm{\infty }}}{\alpha }{\kappa }_{\mathrm{\infty }},\hfill \\ {\kappa }_t=(1-{\alpha }_{\mathrm{\infty }})\frac{{\alpha }_{\mathrm{\infty }}}{\alpha },{\kappa }_{\mathrm{\infty }}=\frac{{\alpha }_{\mathrm{\infty }}}{0.73+0.083\sqrt{{\alpha }_{\mathrm{\infty }}}+{\alpha }_{\mathrm{\infty }}}.\hfill \end{array}$$

The ${\alpha }_{\mathrm{\infty }}$ parameter is model dependent and reads:

$${\alpha }_{\mathrm{\infty }}\simeq \frac{30}{24{\pi }^2}\frac{{\sum }_i{c}_i\mathrm{\Delta }{m}_i^2\left({\varphi }_{*}\right)}{g_{*}{T}_{*}^2},$$

(11)

where ${\varphi }_{*}$ is the field value acquired inside the bubble immediately after the tunneling process, i runs all over the particles, $\mathrm{\Delta }{m}_i^2\left({\varphi }_{*}\right)$ are the squared mass differences in the two phases, and $c_i$ is equal to $N_i$ for bosons and $N_i/2$ for fermions in terms of the numbers of species $N_i$. In Figure 2, we display several cases of runaway bubbles.

7. Trans-Electroweak Phase Transitions in Axion Monodromy Inflation

We now focus on the well-known case of axion monodromy inflation, where the inflaton has a potential of the form:

$$V\left(a\right)=\frac{1}{2}{m}^2{a}^2+{\mathrm{\Lambda }}^{4}cos\left(\frac{a}{f}+\overline{\varphi }\right),$$

(12)

f and m denoting the axion decay constant and mass, $\mathrm{\Lambda }$ standing for the typical non-perturbative axion scale; this is analogous to the QCD dimensional transmutation energy scale of the QCD axion.

In this model, the inflaton is an axion-like field that cannot be recast from the solution of the strong CP problem, but is largely discussed in the string phenomenology literature [29,30,31,32,33,34,35]. Such a scenario is particularly interesting since it is related to a large production of primordial B-modes in the Cosmic Microwave Background radiation. This amounts to a large r-parameter, characterizing the ratio tensor over scalar perturbations, which renders the model falsifiable in the near future [29].

Here, we will show that such a model may be tested also in GWs interferometers, since it is connected to a trans-electroweak phase transition. The scenario represents a particularly important example where our model-independent analysis discussed above can be applied. The axion-monodromy potential has a pattern of periodic local false minima, with the global minima at the bottom of it (see Figure 3).

Finally, we emphasize that in the case of warm inflationary scenarios, where thermal corrections are higher than the Hubble scale, the gravitational waves’ signal may be even stronger, as suggested by [36].

The inflaton rolls down to the bottom of the potential through a series of tunneling transitions from the local false minima. The last tunneling from the last false minimum to the true minimum receives relevant thermal corrections, since it is closest to the reheating process. This can trigger an FOPT at energies around the typical reheating scale $T_R\sim {10}^7÷{10}^9\mathrm{GeV}$ [33] (see Figure 3). Thermal radiative corrections to the axion-like field are provided not only from the axion self-interactions, but also from typical couplings of the axion to the SM fields, including:

$$c_{1}a\overline{f}{\gamma }_{5}f\mathrm{and}{c}_{2}a{F}_{\mu \nu }{\tilde{F}}^{\mu \nu }.$$

Here, f are the SM fermions (quarks and leptons), and F are the field strengths of SM gauge bosons (electroweak and strong gauge bosons). The coupling constants are not necessarily related to the axion mass and the f-scale, like for QCD axions, and in principle, can be just free parameters of the model. Using the standard DESYcorrections, including thermal field theory one-loop corrections techniques (see, e.g., [8] for applications to electroweak phase transitions) adapted to the axion model case, several FOPT points, testable at the ET interferometer, can be found in a large part of the parameter space.

In Figure 4, we show the number distribution of FOPT events testable at the ET as a function of the $\alpha$ (x-axis) and the $\beta$ (y-axis) parameters. The scan for such events was realized considering the following conditions: axion self-interaction scales as ${\mathrm{\Lambda }}^4/f^2{m}^2\simeq 50$, $\overline{\theta }=0$, $m\simeq {10}^{-6}{M}_{Pl}$; $T_{\mathrm{Reheating}}=0.5\times {10}^9$, and couplings $c_{1f}=0.1÷1$, $c_2\ll c_{1f}$, which correspond to the case of the axion strongly coupled to the SM fermions, while all the others couplings with the gauge bosons are negligible. In principle, the axion may be coupled to the Higgs boson, but we do not consider this case, since this would turn this scenario into a technically much more complicated one. We found that most of the events were picked around the point $\beta =12.5$ and $\alpha =0.51$, while far from it, the viable points disappeared. The points found in Figure 4 correspond to the non-runaway bubbles’ case. Let us remark that our parameter space scan probably would not include all possible interesting FOPTs for technical reasons, but still represents well the whole picture and a potential for future experimental analysis of the considered scenario. It would be interesting in the future to further explore even more viable parameter-space islands, which correspond to testable FOPTs.

8. Conclusions and Remarks

We explored the possibility of testing trans-electroweak FOPTs in the early Universe, deploying radio astronomy constraints. The FOPT may occur around ${10}^6÷{10}^8\mathrm{GeV}$ scales and generate a characteristic stochastic potentially observable GW background. Depending on the specific subclasses of cases considered, aLIGO and the ET can measure the frequencies around $1÷100\mathrm{Hz}$, with the possibility of probing FOPT temperature at energies much higher than the electroweak scale, in a range around ${10}^8\mathrm{GeV}$ or so. We remark that our results imply the exciting possibility to test very high energy mechanisms, including inflaton FOPTs, as well as the GUT, baryogenesis, and PBH production models. As an important example, we have made new predictions for the trans-electroweak phase transitions in the so-called axion monodromy inflation and demonstrated the potential for testability in future GW measurements. Our main conclusion is that the new data from aLIGO and the ET will provide a crucial test for FOPTs lying at very high energy scales. This may provide a new important frontier for very early Universe physics.