Gravitational tests of electroweak relaxation

We consider a scenario in which the electroweak scale is stabilized via the relaxion mechanism during inflation, focussing on the case in which the back-reaction potential is generated by the confinement of new strongly interacting vector-like fermions. If the reheating temperature is sufficiently high to cause the deconfinement of the new strong interactions, the back-reaction barrier then disappears and the Universe undergoes a second relaxation phase. This phase stops when the temperature drops sufficiently for the back-reaction to form again. We identify the regions of parameter space in which the second relaxation phase does not spoil the successful stabilization of the electroweak scale. In addition, the generation of the back-reaction potential that ends the second relaxation phase can be associated to a strong first order phase transition. We then study when such transition can generate a gravitational wave signal in the range of detectability of future interferometer experiments.


Introduction
The 2010s decade has been marked by two scientific milestones: the Higgs discovery at the Large Hadron Collider (LHC) in 2012 by the ATLAS and CMS collaborations [1,2] and the first direct detection of gravitational waves (GW) on Earth in 2015 by the LIGO and VIRGO collaborations [3,4]. While the latter has given access to previously unaccessible phenomena, like the merging of black holes binary systems, the first discovery has exacerbated the hierarchy problem, i.e. the question of how the electroweak (EW) scale can be so much smaller than the Standard Model (SM) cutoff without the need for a large degree of fine tuning. Traditional symmetry based solutions like supersymmetry and composite dynamics are nowadays pushed in quite tuned regions of parameter space by the null LHC searches.
This has motivated the scientific community to consider alternative solutions to the problem of the instability of the EW scale. A compelling possibility is the one where the Higgs mass is driven to a value much smaller than the SM cutoff by a dynamical evolution in the early Universe. This mechanism has been firstly proposed in [5] and goes under the name of cosmological relaxation. The basic idea is as follows: the Higgs squared mass parameter H † H is made dynamical by its coupling with a new scalar degree of freedom, the relaxion φ, generally assumed to be a pseudo Nambu Goldstone boson (pNGB). The evolution of the relaxion field during the early Universe evolution, governed by an opportune potential V (φ), scans the Higgs mass parameter, making it evolving from large positive values up to the critical value in which electroweak symmetry breaking (EWSB) is triggered. Once the Higgs develops a vacuum expectation value (VEV), a back-reaction potential turns on and stops the relaxion evolution, dynamically selecting the measured value for the EW scale.
For the mechanism to work, two ingredients are essential: i ) a friction mechanism that slows down the relaxion evolution and avoids the overshooting of the back-reaction barrier and consequently of the correct EW scale, and ii ) a mechanism to generate the back-reaction itself. In much of the explicit realizations of the relaxion mechanism, the friction is provided by the Hubble expansion during inflation [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24] 1 , so that the relaxion field slow rolls during the cosmological relaxation phase. Alternatives are however possible: the friction can be generated by particle production [26,[29][30][31][32], the relaxion can fast-roll during inflation [33], it can be stopped by a potential instability [34], or it can fragment [35,36]. Cosmological relaxation after inflation is also possible [37]. As for the back-reaction, we can essentially distinguish between familon models [10] and models in which the potential barrier is generated by the confinement of a new strongly interacting dynamics with new vector fermions, as already discussed in the original paper [5].
In this work, we focus on the second scenario. This opens up an interesting possibility: after a relaxation phase during inflation in which the EW scale is dynamically selected, the Universe may be reheated to temperatures above the critical temperature of the new confining interactions. If this happens, the back-reaction barrier disappears and the Universe undergoes a second relaxation phase. When the temperature of the Universe drops again below the confinement scale of the new strong dynamics, the barrier is once again generated and the relaxion stops again its evolution. Crucially, depending on the gauge group of the new confining dynamics, the number of new fermions and their representations under the gauge group, this phase transition can be of first order [38] and can thus give rise to a stochastic GW background [39]. The signal might be then detected at present and future interferometer experiments [40][41][42][43][44][45][46][47][48], thus connecting in this way the two milestones discovery of the 2010s. A sketch of the situation we are considering is shown in Fig. 1.
The paper is organized as follows: in Section 2 we review the cosmological relaxation mechanism, where the back-reaction potential is generated by a new strongly interacting dynamics. In Section 3 we analyze the relaxion evolution after reheating.  Figure 1: Sketch of the framework we are considering. The thick black line represents an approximate sketch of the temperature evolution. During inflation, a first relaxation phase selects the correct EW vacuum. When reheating happens at t RH , the temperature increases and the radiation domination phase begins. For simplicity, we will consider reheating to be an instantaneous process. If the reheating temperature T RH is larger than the nucleation temperature T n , the back-reaction potential disappears and a second relaxation phase occurs. After a while, the temperature falls below the nucleation temperature of the new strongly interacting group, and a phase transition occurs which may generate a GW signal.
We show how the equation of motion of the relaxion during this second relaxation phase can be analytically solved for a certain range of temperatures, and what are the bounds on the parameter space that we obtain by requiring the additional relaxation phase not to spoil the dynamical selection of the EW scale achieved during the first relaxation phase. In Section 4 we discuss the GW spectrum generated in the first order phase transition associated with the generation of the back-reaction potential after reheating and its detectability at future experiments. We thus conclude in Section 5. We also add two Appendices where we collect more technical material for the interested reader. In Appendix A we discuss in detail various models of strongly interacting vector-like fermions and review how to describe their low energy dynamics and in particular their vacuum energy, which is relevant for the relaxion mechanism. In Appendix B we instead explicitly show how the constraints on the relaxion parameter space used throughout the paper are derived.

Relaxation with strongly interacting fermions
In this section we briefly summarize the relaxation mechanism, highlighting its main features. Following [5] we consider a scalar potential where Λ is the cutoff of the theory, λ, and r are positive dimensionless parameters and V BR (H, φ) is the back-reaction potential. The dependence of the first term of Eq. (2.1) on φ makes the Higgs mass squared parameter a dynamical quantity. As for the back-reaction potential, V BR (H, φ), we consider it to be of the form where Λ 4 BR ( H ) is a (model dependent) monotonically increasing function of the Higgs VEV and F is the scale at which the NGB φ appears. Notice that the backreaction potential respects a discrete shift symmetry φ → φ + 2πF . This symmetry is explicitly broken by the spurion , which can thus be taken small in a technically natural way.
We assume that the cosmological relaxation mechanism takes place during inflation, with φ evolving in a slow-roll regime. We furthermore take the initial value of the field φ to be sufficiently small to guarantee the condition Λ 2 − Λφ in > 0 to be satisfied at the beginning of the relaxation phase, thus implying an unbroken EW symmetry.
As φ evolves, increasing its value due to the −r Λ 3 φ term in the potential, the Higgs squared mass parameter becomes smaller and smaller, until it crosses zero and causes the breaking of the EW symmetry. Once H = 0, the back-reaction potential is switched on. The amplitude of the oscillating term grows with the Higgs VEV, until it's large enough to stop the evolution of φ, dynamically selecting the value of H . This is the essence of the relaxion mechanism (see Fig 7).
A possible origin for the back-reaction potential, already considered in the original work [5], involves fermions charged under new strong interactions G dark as well as under EW interactions. Provided that the scale of inflation is smaller that the confinement scale of G dark , which we will denote as g ρ f , and that the Higgs boson interacts with the new fermions, the back-reaction potential forms and, as described in App. A, Λ BR takes the form where the constants µ 2 B and Λ 4 0 depend on the specific model considered. The form of Eq. (2.3) is however generic. We conclude this section by listing the conditions that are needed to guarantee the successful realization of the relaxion mechanism, see Appendix B for more details.
Validity of the EFT: We are assuming the field φ to be an angular degree of freedom, with a decoupled radial mode. Since the mass of the radial mode is O(F ) we obtain the condition F Λ. We will also require F M Pl , where M Pl 10 18 GeV is the Planck mass.

Conditions on inflation:
We ask the dynamics of inflation to be decoupled from the one of the relaxion. To achieve this we require the relaxion energy density ρ φ ∼ Λ 4 to be smaller than the inflaton energy density, ρ infl ∼ H 2 I M 2 Pl , where H I is the Hubble scale during inflation. This implies the lower bound Since we are assuming that the back-reaction is generated by some new strong force, we need to guarantee that the barriers can form during inflation. This requires where Λ d denotes the scale of confinement of the new strong interactions.
We also demand the classical relaxion evolution not to be spoiled by quantum fluctuations. To achieve this condition we require the classical spread of the relaxion field in one Hubble time, ∆φ cl ∼φ/H I ∼ V /H 2 I , to be larger than the quantum spread ∆φ quantum ∼ H I , where V indicates the relaxion potential and the derivative is with respect to the field φ. This implies a second upper bound on the scale of inflation that reads Conditions on the back-reaction: The term Λ 4 0 is present in the back-reaction potential even before EWSB. To guarantee that it does not stop the relaxion evolution we require In addition, to ensure that after EWSB there will be a period of evolution in which the height of the barrier grows, we need The EW scale is an output of the relaxation: An approximate expression for the EW VEV in terms of the parameters of the model is As shown in Appendix B, the cutoff satisfy Λ 6.7 × 10 6 GeV µ B 10 GeV favoring thus very small values of . This result can be problematic in two aspects: i ) to completely solve the hierarchy problem an additional protection mechanism must be present for scales between Λ and M Pl , and ii) a successful relaxation requires the relaxion excursion to be at least ∆φ ∼ Λ/ in order for the Higgs mass parameter to change sign. Given the very small needed for the mechanism to work, the resulting excursion is transplanckian. The first issue can be solved assuming supersymmetry [49][50][51], a composite dynamics [52] or a warped dimension [53] to be present above Λ. The second problem requires more model building effort, but can be solved in the context of clockwork models [54][55][56][57]. In the following we assume that one of these mechanisms is present to stabilize the EW scale all the way up to the Planck scale, and focus only on the effective field theory defined by Eq. (2.1). Let us conclude this section with a comment on the compatibility between the well-known upper bound on the reheating temperature, T RH √ H I M P L , and the conditions on the Hubble parameter during inflation needed for the relaxation mechanism to work. If we take H I close to the lower bound of Eq. (2.4) we see that T RH Λ, i.e. the phase following reheating can be described inside the validity of the EFT considered. Since the value of the Hubble constant during reheating has no effect on our subsequent considerations, we will take H I ∼ Λ 2 /M P L from now on.

Evolution after reheating
As already mentioned, we are considering a situation in which the relaxation of the EW scale happens during inflation, i.e. the correct EW VEV is selected at the end of this phase. This leaves open the question of what happens after reheating? One possibility is for reheating to leave the Universe in a bath with temperature T RH Λ d . In this case there is no further dynamical evolution in the relaxion field direction, since φ remains stuck in its minimum. This is what is implicitly assumed in the original work [5]. If T RH > T EW , with T EW the scale of EW phase transition, we expect thermal corrections in the Higgs field direction to recover the EW symmetry. As the Universe expands and cools down, EWSB is again triggered as usual 2 .
A second possibility, on which we focus in the following, is for reheating to happen at T RH > max(T d , T n ), where T d and T n are the temperatures associated with the deconfinement and confinement of the additional strong interaction G dark . If this is the case, after reheating the back-reaction disappears and the Universe undergoes a second period of relaxation. This possibility has received less attention from the literature, see e.g. [8,25,58,59]. As the Universe expands and cools down, it will reach a temperature T ∼ T n at which the new strong sector again confines, thus producing again the back-reaction barrier and ultimately stopping the relaxion evolution 3 . As we are going to see at the end of this Section, the relaxion evolution stops soon after the barrier is again formed. If the transition producing the backreaction barrier is strongly first order, it will produce a GW signal that might be observable at interferometer experiments, as we will study in Sec. 4. This allows to open a new window on this type of solutions of the hierarchy problem.
Let us now discuss more in detail the relaxion evolution when T RH > max(T d , T n ). We want to understand whether this additional relaxation phase can spoil the solution of the hierarchy problem. We achieve this by imposing additional conditions on the parameters in order to achieve the correct EW minimum today. We will also make the following simplifying assumptions: (i) reheating is instantaneous and (ii) the energy density of the relaxion after reheating is subdominant with respect to the energy density in radiation. All our assumptions aim at avoiding additional model building related to the reheating or the deconfinement phases 4 . We also ignore a potential gravity wave signal generated during the deconfining phase, since the spectrum would depend on the details of the process. While interesting, we leave the study of these problems to future work, focussing only on what happens after reheating. With our assumption that the barrier disappears when the universe enters the radiation domination phase after inflation we need to consider three situations, depending on the relative hierarchy between T RH , T n and T EW : • If T n < T RH < T EW the Universe is reheated to a phase in which the scalar potential has already a non-trivial minimum in the Higgs direction, while there are no barriers in the relaxion direction. Using the potential in Eq. (2.1) in the equation of motionθ + 3H(t)φ + ∂V /∂φ = 0 we obtain is the Higgs VEV in the absence of the back-reaction barrier. We have used H = 1/(2t), as appropriate for a radiation dominated Universe. The solution for temperatures T n < T < T RH can be written analytically in terms of modified Bessel functions I n and K n of first and second kind, respectively 5 . To write a compact expression we define a new time variable τ = Λt/ √ 2λ, the functions and the constant term The solution can be written as where the primes denote the derivative with respect to τ , and the subscript 0 indicates that the quantity must be computed at the initial time, i.e. at reheating. Since we are assuming reheating to be an instantaneous process, we can identify θ 0 with the value of the relaxion field at the end of the inflationary period and we can assume that the relaxion field does not acquire a relevant dynamics during reheating, leading toθ 0 0. We also notice that for large values of τ we have the asymptotic behavior The dependence on hyperbolic trigonometric functions implies that the relaxion field will evolve very quickly, if enough time is allowed to pass. We will comment later on the dynamics taking place for T < T n ; • If T n < T EW < T RH there are two phases in the evolution. The first one applies for T EW < T < T RH while the second one for T n < T < T EW . In the first phase the equation of motion to be solved is Eq (3.1) with v 2 (θ) = 0. The solution can again be written analytically: Notice that the evolution is a power law in this regime, and the relaxion field evolves much less than what happens when v 2 (θ) = 0. Using the relation between time and temperature in a radiation domination Universe where g denotes the number of relativistic degrees of freedom we can rewrite the solution as (3.9) Once EWSB is triggered the equation of motion to be solved is again Eq. (3.1) with v 2 (θ) = 0, whose solution is given in Eq. (3.5) with τ 0 ≡ τ EW ; • The last case is the one in which T EW < T n < T RH . For T n < T < T RH the solution is given by Eq. (3.7), while we will discuss in the next paragraph what happens for T < T n .
To compute the solution for T < T n we need to turn on not only the Higgs VEV, for T < T EW , but the back-reaction barrier as well. To the best of our knowledge, no analytical solution can be found in this case. We can however qualitatively expect that, once the barriers form, the relaxion will find itself trapped in a position that is displaced from the minimum in the relaxion direction. It will then oscillate around this minimum losing energy. It is in this phase that it can behave like dark matter, as studied in [59]. Whether or not the relaxion stops its evolution when encountering the first barrier will depend on the velocity it has acquired during the second relaxation phase, and cannot be inferred analytically. As shown at the end of App. B however, the change in the Higgs VEV between two subsequent minima is relatively small, so that it is likely that the second stopping phase will not modify dramatically the value of the Higgs VEV at the end of the relaxion evolution. What can modify dramatically this value is however the relaxion evolution before the barrier is again formed. To understand whether this is the case, we show in Fig. 2 and Fig. 3 in blue the regions in which In these regions the Higgs VEV at the end of the second relaxation phase differs by more than a factor of 2 with respect to the observed VEV, and the solution of the hierarchy problem is spoiled. In Fig. 2 we consider the case in which T n < T RH < T EW , while in Fig. 3 we show the case in which T n < T EW < T RH . We take T EW 160 GeV. We have fixed the parameters of the model to the same representative values chose in Eq. (2.10), with the exceptions of and F , whose values are reported in the plots. For simplicity we have taken a fixed value of g = 100, although by varying it the overall picture does not change. Furthermore, we have supposed that the initial relaxion value at reheating gives v EW and has vanishing velocity. The quantity v[θ(T n )] is computed using Eqs. (3.5) and (3.9). We also impose a lower bound T n , T RH 4 MeV, represented by the red regions, where the limit on T RH is taken from [61] while the limit on T n is imposed to be safely far from the Big Bang nucleosynthesis epoch T ∼ 1 MeV. By inspecting Fig. 2 we see that  .10) is larger than unity (blue) and therefore the solution of the hierarchy problem is spoiled and in which the nucleation or the reheating temperatures are smaller than 4 MeV (red). We suppose the initial VEV to be v EW . All the plots show the T n < T RH < T EW case. See the text for more details.
for reheating temperatures below the EW phase transition one, there are regions in parameter space in which the solution of the hierarchy problem is completely spoiled, depending on the choice of parameters. This is due to the exponential evolution of the relaxion field in this regime. On the contrary, when T n < T EW < T RH most of the relaxion evolution happens when Eq. (3.9) is valid. Since the relaxion field does not evolve much in this regime, there are large regions of parameter space in which the modifications of the VEV are small. Finally, when T EW < T n < T RH the evolution  .10) is larger than unity (blue) and therefore the solution of the hierarchy problem is spoiled and in which the nucleation or the reheating temperatures are smaller than 4 MeV (red). We suppose the initial VEV to be v EW . All the plots show the T n < T EW < T RH case. See the text for more details.
follows again Eq. (3.9) and is not enough to drastically modify the value of the VEV. For this reason we do not show any plot for this case. Let us finish this Section noticing that we can imagine a situation in which the VEV at the end of the first relaxation phase is smaller than the observed Higgs VEV. If this is the case, the second relaxation phase could provide the additional evolution needed to be compatible with experiments. In this paper we will focus on a situation in which the correct Higgs VEV is selected at the end of the first relaxation phase, leaving for future work the analysis of what happens in the opposite situation.

Gravitational waves signal
In relaxion models with strongly interacting fermions there are two main sources of GW: i) the confinement of G dark if it proceeds through a first order phase transition [62,63] and ii) the possible penetration of the barrier by the relaxion field before stopping at the minimum [64][65][66]. 6 The last process is much more difficult to analyze we will focus here on the first signal, leaving the second one for future work. Since the confinement dynamics of G dark is a non-perturbative process, it is difficult to make quantitative predictions on the GW spectra that can be produced. The phase transition ultimately depends on the phase diagram of the theory, whose determination requires computationally expensive lattice computations. We can however have an order of magnitude estimate of the expected signal using effective models to parametrize the confining dynamics. As shown in [62], in QCD-like theories different effective models give GW spectra whose peak amplitude might differ even by two orders of magnitudes between each other. Nevertheless they can give a useful indication of what kind of theories can produce a GW spectrum that is close enough to current and future experimental sensitivities, and that thus deserve further dedicated theoretical studies to allow for a more precise study of the phase transition and the associated GW signals. In this work we decide to focus on the linear sigma model description of a confining strong dynamics, which is defined in Eq. (A.7) of Appendix A.1. We nevertheless stress that our results should be taken as a preliminary indication for what the real GW spectrum could be, as recently also emphasized in [62,63]. Furthermore, it is important to stress that there are many sources of theoretical uncertainties in the computation of the GW spectrum. Among others, we can list the renormalization scale dependence, higher loop corrections, different ways to deal with the high temperature approximation, gauge dependence, and non-perturbative and nucleation corrections (see [67]). It is also worthwhile to say that not all the uncertainties can be accurately quantified. For simplicity, we do not consider these additional uncertainties in our estimate of the GW spectrum since, as we have seen, we are describing non-perturbative processes that are, in any case, difficult to deal with.
As already mentioned in the Introduction, a well-known argument by Pisarski and Wilczek [38] implies that SU (N d ) gauge theories with N d ≥ 3 feature a first order phase transition if N F ≥ 3 light 7 flavors are present. Motivated by this result, we 6 As noticed above, we ignore a possible signal coming from the deconfining of the strongly interacting group. We have estimated that, in ample regions of the parameter space shown in Figs. 2 and 3, the time passed between deconfinement and confinement is sufficiently long to avoid interference between the GW signals generated in the two processes. 7 That is with a mass smaller than the new gauge theory confinement scale. study here SU (N d ) gauge theories with N F = 3 and N F = 4 flavors, considering both the situations in which the new confining phase transition happens at a scale below and above the EW critical temperature T EW . A general discussion of these models in the relaxion framework can be found in App. A.2. We focus on the regime in which the anomaly term dominates over the mass term. Since the only dependence of the Sigma model potential on the Higgs field comes through the mass term, we expect the Higgs dynamics to have only a very small impact on the phase transition. For simplicity, we will neglect such impact in our computation. We keep our discussion general, but to make it more concrete it is useful to take as an example the variations of the L + N model shown in Fig. 6, where L and N are new vector-like fermions in the fundamental of SU (N d ) with the quantum number of the SM lepton doublet and of a total singlet, respectively. We see that the physics of the N F = 3 case is captured by the minimal L + N model with dark confinement scale above the EW one (model A), as well as the physics of models in which we add 3 singlets to the spectrum independently on the confinement scale (models B and D). Notice that the EW charged fermion L cannot be too light to evade current experimental direct searches bounds, and we then assume it to have a mass larger than Λ EW . The physics of the N F = 4 dark flavors, on the contrary, is captured by the minimal model with two copies of L fermions and a decoupled singlet (model C), and of models in which L is decoupled and 4 light dark N flavors are added to the spectrum (models B and D with n N = 4). Experimental bounds on the minimal L + N model have been considered in [68][69][70]. Although important, these bounds allow for rather different spectra at the level of sigma model without spoiling the successful relaxation of the EW scale.

Computation of the gravitational wave spectrum
We now quickly remind the reader how the computation of the spectrum of the stochastic GW background produced in a first order phase transition proceeds. Three contributions must be considered [71]: the one from true vacuum bubble collisions, Ω col h 2 , the one from the propagation of sound waves in the plasma, Ω sw h 2 , and the one from magnetohydrodynamic turbulence effects, Ω turb h 2 . A fourth contribution may be generated by quantum fluctuations [72], but since its impact is not well understood we will not consider it in the following. As mentioned, we describe the confining dynamics of the new strongly interacting sector through the linear sigma model described in Eq. (A.7). In this setup, the radial degree of freedom σ interacts with the light particles in the plasma. We thus expect the interactions between the scalar shell and the plasma to be important, causing the behavior of the bubble to be a non-runaway one [40], that is the bubble walls do not keep accelerating until the bubbles collide. Since in this case most of the latent heat of the phase transition is transferred to the plasma in the form of sound waves and turbulence, in the following we will focus on the Ω sw h 2 and Ω turb h 2 contributions only. Their expressions can be written in a compact way as [40,41,62,73] where i = sw or turb for the sound waves and turbulence contributions, respectively. The details of the phase transition enter in the parameters α, β/H and T n , which we will discuss in detail below. As for the other parameters, the numerical coefficients a i in the two cases are given by [40,41,62,73] a sw = 2.65 × 10 −6 , a turb = 3.35 × 10 −4 (4.2) and the exponents b i are b sw = 2 , b turb = 3/2 . Again, g denotes the number of relativistic degrees of freedom and H is the Hubble parameter, where both quantities are computed at the nucleation temperature. The quantity v wall is the wall velocity. It is known that scenarios with large v wall lead to stronger GW signals [40,62,63,74]. While in the non-runway behavior the bubble walls stop their acceleration and reach a terminal velocity, this velocity might still be relativistic. Calculating the specific value of v wall is beyond the scope of the present work. We will then concentrate on the case of highly relativistic bubbles, v wall ∼ 1, since it is the most interesting regime from an observational point of view. Decreasing v wall down to v wall ∼ 0.75 doesn't drastically modify the overall picture, see e.g. [62]. Finally, the coefficients κ i are the efficiencies of each process. The former is the efficiency to convert the latent heat into bulk motion, while the latter is the part that is converted into vorticity in the plasma. Numerical fits of κ sw suggests the form [40,74] κ sw α 0.73 + 0.083 for v wall ∼ 1 where the quantity α is related to the strength of the transition and will be defined below. κ turb is determined trough κ turb = κ sw . Previous studies suggest the conservative value of = 5 × 10 −2 [40,75,76]. As regarding the spectral shapes S i (f ), they read In the previous expressionh is the Hubble rate at t n redshifted to today, h = 1.65 × 10 −5 Hz T n 100 GeV g 100 while the peak frequencies in the two cases are given by similar expressions, with f 0 sw = 1.9 × 10 −5 Hz and f 0 turb = 2.7 × 10 −5 Hz. It has been suggested in [41,73,76] that when β/H 1, the sound wave and the turbulence contributions shown above overestimate and underestimate the signal, respectively. Following the suggestion in the same works, we modify Ω sw h 2 and Ω turb h 2 to where is related to the duration of sound waves and U f is the root-mean-square four velocity of the plasma.
Let us now discuss the parameters directly related to the phase transition under consideration: the nucleation temperature T n , the inverse time duration β/H and the strength of the transition α. The first two parameters are defined in terms of the nucleation rate [77,78] Γ where S 3 is the Euclidian action computed at its bounce and A (T ) is a non-perturbative factor. The nucleation temperature is defined as the temperature at which one of the nucleated bubbles reaches a size comparable to the Hubble radius at that time, Γ(T n ) H 4 (T n ). Considering the process to take place in the radiation domination epoch, the above statement implies [72,79] S 3 T n 164.56 − 2 log g * 100 − 4 log T n 1 GeV . (4.10) The exact numerical value of the first term on the right-hand side depends on the precise form of A (T ) , for simplicity we assume A (T ) T 4 . The inverse time duration β/H is instead defined as the rate of change of Γ at the nucleation time t nucl via [80] β The strength of the transition, α, is instead not univocally defined. In the literature many definitions can be found [62,63,73,74]. The main two are based on the latent heat and on the trance anomaly. In the first one α is identified as the ratio between the latent heat of the transition and the energy density of the radiation in the plasma ρ rad . In the second one it is defined through the trace of the the stress-energy tensor and ρ rad . From the practical point of view, both definitions can be summarized as  where ∆V = V eff (σ false , T ) − V eff (σ true , T ) is the difference of the free energy density between the two phases and n = 1 (1/4) when α is defined via the latent heat (trace of the stress-energy tensor).

Gravitational wave spectrum: QCD-like case
We now compute the GW spectrum using the linear sigma model defined in App. A.1 with the formalism of Sec. 4.1. We refer the reader to the Appendix for a definition of the notation we use for the linear sigma model. We focus here on the situation in which the masses of the σ and η mesons satisfy m σ,η f π , where f π is the pion decay constant. Since this is the well-known spectrum of QCD, we use the term "QCD-like" to refer to this case. We will study in Sec. 4.3 a different physical spectrum. In our computation we consider the chiral limit, in which the mass of the lightest fermions is much smaller than the confinement scale Λ d and we focus on the N F = 3 and N F = 4 cases. We present our results for two specific values of T n : one where T n < T EW and one where T n > T EW , thus capturing the different possibilities discussed in Sec. 3. In particular we fix the values of the linear sigma model parameters for the two cases as reported in Tab. 1, where we also indicate the values of the chiral symmetry breaking VEV, the relevant mesons masses, the nucleation temperature T n and of the main parameters entering the GW spectrum computation. In choosing these benchmark point we have scanned the linear sigma models parameters trying to maximize the GW spectrum for both nucleation temperatures. In the lack of analytical expressions for α and β/H this represents a challenging numerical process. As already stressed in [62] it's interesting to notice that large values of β/H are found. This is to be contrasted with the naive estimate Since for the sound waves and turbulence contribution the GW amplitude decreases linearly with β/H, see Eq. (4.1), while the peak frequency increases linearly with it, see Eq. (4.7), the direct computation of β/H through an explicit, although effective, model has a strong consequence for the observability of the GW spectra. All together our results for the QCD-like case are shown in Fig. 4 for the N F = 3 and N F = 4 models in the upper and lower panels respectively. There, the left and right panels correspond to the T n < T EW and T n > T EW cases. In all plots the blue line shows the spectrum computed using Eq. (4.1), while the red line is computed using the modification presented in Eq. (4.8), which explicitly show the suppression factor due to the decrease of the sound wave contribution. The situation depicted in the left panels can be realized, in the relaxation framework we are considering, when T n < T RH < T EW or when T n < T EW < T RH . As shown in Sec. 3 in the first case large regions of parameter space end up with large variations on v EW from the second relaxation phase, leaving thus T EW < T RH as the preferred case 8 . The situation of the right panels refers instead to the case T EW < T n < T RH , for which we showed that large regions of parameter space are compatible with a small variation of the EW VEV during the second relaxation phase. In all the figures the colored regions represent the sensitivities of future interferometer experiments. In particular we show the projected reach from AstroD-GW [42,43], eLISA [42,81], BBO [42,44], DECIGO [42,44], B-DECIGO [42,44], AION [45], MAGIS [46], ET [47] and CE [48].
As we see from Fig. 4, both for N F = 3 and N F = 4 cases the GW signal from the dark phase transition is a few orders of magnitude below the region that can be probed in future experiments, in agreement with previous results obtained in similar frameworks [62,63]. We stress again, however, that these results must be interpreted as an order-of-magnitude estimate since, as shown in [62], different effective models for the strong sector confining dynamics can give results that might differ even by two orders of magnitude for the amplitude of the signal and that might then fall on the edge of detectability. Notice also that changing N F does not dramatically change the region in which the signal falls, making thus difficult the identification of the underlying model in case a signal is detected. We conclude by emphasizing once more that, according to the discussion in Appendix A.2, the strongly interacting models generating the back-reaction potential in the relaxion mechanism do not suffer from very strong experimental limits. In particular, it is not difficult to obtain in such models the spectrum used in this section.

Gravitational wave spectrum: non QCD-like case
We now turn to the discussion of more exotic strongly interacting models, i.e. models in which, unlike the case of QCD, the theory spectrum features m σ f π . The behavior of gauge theories with different number of flavors has been studied on the lattice in a certain number of situations, with the surprising result that a non-QCD behavior in which the σ meson is lighter than expected might emerge. In particular, light composite σ scalars have been found in an SU(3) gauge theory with 8 flavors in the fundamental [82][83][84][85], with 2 flavors in the symmetric representation [86][87][88] Non QCD-like case N F = 3 3  and with 4 light and 8 heavy flavors [89]. They also appear in SU(2) theories with one adjoint flavor [90]. The behavior seems to be quite generic and is typically associated with gauge theories near to their conformal limit. In the case of N F = 3 and N F = 4, this is expected to happen when the number of colors is larger than 4 and the fermions transform in the antisymmetric representation [91]. In all the cases mentioned above, the σ meson is found to be roughly degenerate with the pions, at least in the limit of large chiral symmetry breaking. Such behavior can be captured by a sigma model, as shown in [92]. The limit of small chiral symmetry breaking is however more difficult to describe on the lattice. In the absence of conclusive data, we will suppose that the physics is correctly captured by a sigma model, at least in a first approximation. In [93] such sigma model has been extended to include the effects of the η mass. An interesting result emerging from the analysis is that the degeneracy between m σ and m π can be explained by an approximate cancellation between the VEV and m η . If this happens in the chiral limit, an immediate consequence is that typically m σ VEV, and thus m σ f π . This is the situation studied in this section. Our results are shown in Fig. 5, where the color codes are the same as in Fig. 4.
As we see, allowing for m σ f π allows to boost the GW signal amplitude, in agreement with the results of [63], where detectable GW spectra where found for large values of m η /m σ . Altogether we find that, while in the N F = 3 case also in the case of non QCD-like theory the predicted GW signal lies well below the reach of future experiments, in the N F = 4 case, the signal could be potentially detected     both for the T c < T EW and T c > T EW cases, in the frequency range 10 −3 Hz − 1 Hz. Numerically, we have also explicitly checked that other than a large m η /m σ , also the condition of a small m σ /f π should be satisfied to enhance the signal towards the reach of future experiments. If any one of these two conditions fails to be satisfied, the signal typically lies well below the region of future detectability.

Conclusions
In this paper we have considered the framework in which the EW scale is stabilized through the relaxation mechanism. We have assumed that this happens during infla-tion and that the back-reaction potential needed to stop the relaxion evolution is generated by new vector-like fermions charged under a new strongly interacting SU (N d ) gauge group. We have focused on a configuration where the reheating temperature is above the confinement scale of the new strong dynamics. This causes the restoration of a deconfined phase after inflation, the disappearance of the back-reaction potential and the presence of a second relaxation phase during the early Universe thermal evolution. This second relaxation phase can in principle completely spoil the solution of the hierarchy problem. Once the temperature of the plasma drops again below the confinement scale of the new strong dynamics, the barrier again forms and the relaxation finally stops its evolution. Crucially, the phase transition between the confined and unconfined dynamics might be strongly first order and can then produce a stochastic GW background, that can be detected at present and future interferometer experiments. We have studied the relaxation evolution during the second relaxation phase, finding analytical solutions for its equation of motions for various ranges of temperatures. In particular, we have shown that, depending on the model parameters and on the hierarchy between the nucleation, reheating and EW phase transition temperatures, there are ample regions in parameter space where the second relaxation phase does not spoil the solution to the hierarchy problem. Such regions are large when T n < T EW < T RH and T EW < T n < T RH , but are typically small or inexistent when T n < T RH < T EW .
We have then studied the GW signal that can be generated during the confining phase transition that ends the second relaxation stage, considering SU (N d ) gauge theories with 3 and 4 light flavors present in the spectrum. To quantitatively describe the strong dynamics we have employed a linear sigma model, considering both QCDlike spectra, in which the σ meson is heavier than the symmetry breaking scale, and non QCD-like spectra, in which the σ meson can be lighter than it. The latter behavior may emerge in theories close to their conformal window, although additional lattice studies are needed to establish whether this is the case or not. While in the first case we find that the predicted signals lie below the present and future experimental sensitivities, in the case of non-QCD like spectra signals close and within the experimental reach can be obtained for the N F = 4 case. We however observe that, even if a GW signal will be detected in the future, the reconstruction of the underlying model will in general be challenging. On the one hand, as we have shown, there is little difference in the signal shapes expected for N F = 3 and N F = 4 cases we have analyzed. On the other hand, many different models of strongly interacting vector-like fermions can give rise to the same relaxion backreaction potential and can be described through the same linear sigma model studied in this paper. We finally stress that all results obtained in this work by describing the dynamics of a strongly interacting theory through effective models suffer by large uncertainties, that can affect the peaks positions and heights of the predicted GW spectra [62]. Nevertheless we believe that is of paramount interest that BSM physics that can offer a solution to the hierarchy problem through the relaxion mechanism might generate a GW signal in the range of detectability of future experiments, and this makes even more important a thorough study of such theories through first principle calculations.

Acknowledgments
We are indebted to Djuna Croon and Rachel Houtz for clarifications on their related works, and to Bithika Jain for collaboration in the early stages of the project.

A Strongly interacting models for the relaxation of the EW scale
We collect in this Appendix some useful formulas regarding strongly interacting vector-fermion models and their vacuum energy. We start with some general results and then specialize to the models used in the relaxion framework.

A.1 General setup
Consider N F vector-like fermions Ψ charged under a new confining group, which for simplicity we take to be SU (N d ). We will assume a situation similar to what happens in QCD, namely that there is a small explicit breaking of the chiral symmetry SU (N F ) L × SU (N F ) R due to non vanishing fermion masses. Moreover, we take the axial part of the global flavor group to be anomalous. The Lagrangian we consider is where M is the n × n mass matrix of the vector-like fermions Ψ and D µν is the field strength tensor of the SU (N d ) gauge group, withD µν its dual. Using a SU (N F ) L × SU (N F ) R transformation 9 the mass matrix can always be put in the form where ϕ M = arg det M and M i ≥ 0. We use this basis in the following. We can write the low energy theory by using the following transformations and spurions under an axial transformation We collect the low energy degrees of freedom in a matrix where σ is the radial degree of freedom, T a are the SU (N ) generator, S a are CPeven scalars and U is the matrix of the NBGs, including the dark η . Notice that we denote the VEV with v, as opposed to the EW VEV, which has been called v EW . We write U explicitly as Following [93] we define the pion decay constant in a theory with N F flavors as The most general potential invariant under SU ( where · denotes the trace and we have used an axial transformation to put the phase dependence in the anomaly term, θ = θ 0 + ϕ M . Let us start with the computation of the vacuum energy, since it is essential to write the relaxion potential in Eq. (2.1). In doing this we ignore the heavy degrees of freedom σ and S a in Eq. (A.7). Using the results of [94], working in the basis in which the fermion mass matrix is diagonal forces the Σ matrix in the vacuum to be diagonal. We write it as which gives a vacuum energy In the previous equation we have usedB = Bv 2/N F andμ Σ = 2µ Σ (v/ √ 2N F ) N F . The Dashen equations give the minimum conditions We can find approximate solutions to the Dashen equations when the anomaly term dominates over the mass term, which is equivalent to assume that the mass of the η is larger than the masses of the mesons. In this limit we have The remaining Dashen equations can now be written as In the case in which all the fermion masses are approximately of the same order the previous equation is solved by θ i θ k θ/N F . This means that in the limit considered, i.e. dominance of the anomaly term over the mass term and approximate degeneracy of the fermion masses, the vacuum matrix is given by If instead a hierarchy is present among the fermion masses the situation drastically changes. For definitiveness, let us consider the hierarchy M 1 M i , i = 2, . . . , N F when the anomaly term dominates. Eq. (A.12) is now approximately solved by θ 1 θ and θ i 0 for i = 2, . . . , N F . The conclusion is that, when a clear hierarchy is present among the fermion masses, only the lightest fermions contribute significantly to the vacuum energy. As last step, we remind that the inclusion of the relaxion field can be achieved simply promoting θ to a dynamical parameter, We now move on with the discussion of the potential of the σ particle, since it drives the phase transition discussed in Sec. 4. Having discussed the effect of the VEV of the light modes in the computation of the vacuum energy, see Eq. (A.9), we now simply focus on the potential driven by σ. To write the full potential we implement thermal effects using the Truncated Full Dressing (TFD) procedure [95] V eff (σ, The first term, V (σ), is the tree level potential as a function of the homogenous background field σ, and can be obtained from Eq. (A.7). Following refs. [63,96] we consider the Coleman-Weinberg contribution already included in the tree level term, since its inclusion just renormalizes the tree-level couplings. Each thermal contribution depends on the multiplicity n i , on the bosonic thermal integral and on the thermal masses M i (σ, T ). In a sigma model with N F flavors they read The terms proportional to µ Σ are present only for N F ≥ 3. We have already included the "Debye" masses Π(N F ) computed at one loop, the so called "hard thermal loop" [95]. In a theory with N F flavor we obtain

A.2 Explicit models
In the original paper [5] the back-reaction barrier is generated by the so-called L + N model. This is a theory in which vector-like fermions charged under a new confining group are introduced. These fermions have EW quantum numbers to allow for Yukawa interactions with the Higgs boson. More specifically, the model consists in a vector-like pair L, L c (where L has the same quantum numbers as the SM lepton doublet and L c has conjugated charges) and by a second vector-like pair N , N c of SM singlets. This is only one of the possibilities since, as we are now going to show, all models that allow for a Yukawa interaction with the Higgs doublet and in which there is a clear mass hierarchy can generate the required back-reaction. Before starting, let us remind the reader that in light of the Pisarski-Wilczek argument [38] already mentioned in Sec. 4, we require the presence of at least 3 light flavors below the confinement scale to produce a strong first order phase transition. The most general Lagrangian we consider is where the quantum numbers of the vector-like fermions are such that it is possible to write the Yukawa interactions. Moreover, M L , M N , Y and Y are matrices whose dimensions depend on the number of fermions. Suppose now there is a hierarchy M χ M ψ . Integrating out the heavy fermions we obtain the effective Lagrangian As explained in Sec. A.1, the computation of the vacuum energy can be done analytically when all the fermions are approximately degenerate or when there is a clear mass hierarchy, in which case only the lightest fermions contribute. This allows us to conclude that the heavy ψ states do not substantially contribute to the vacuum energy even if their masses happen to be below the confinement scale. The problem thus is to compute the eigenvalues of the χ mass matrix. To write approximate analytical formulas we take all mass matrices to be proportional to the identity, M χ = m χ 1 and M ψ = m ψ 1 and all Yukawa matrices to be real with equal entries, Y ( ) ij = y ( ) . In this limit the vacuum energy is equal to where n χ and n ψ are the number of χ and ψ fermions. This equation justifies the form of the back-reaction used in Eq. (2.3). Experimental limits on the L + N model have been studied in [68] for the situation in which only the N flavor confines, and in [69,70] for the situation in which both L and N confine. Although relevant, none of these bounds put significant restrictions on the parameters of the linear sigma model used in the computation of the GW spectra. Let us comment on two further points before concluding this section. There is a caveat to the above argument: when the fermions that form bound states carry EW quantum numbers, loop contributions to the effective Lagrangians can be important, see e.g. [69]. We have assumed so far that such loops are negligible, but this is not necessarily the case. When loops are important the argument leading to Eq. (A.8) fails, and the computation of the vacuum energy can in general be done only numerically. Finally, let us give some concrete example of models that can lead to strong first order phase transition. Focussing for simplicity only on variations of the original L + N model used in Ref. [5], we summarize the different possibilities in Fig. 1.
A B C D Figure 6: Summary of the models that can generate the back-reaction for the relaxion mechanism to work and, at the same time, satisfy the Pisarski-Wilczek condition. We show only variations of the L + N model defined in Sec. A.2, although more general possibilities are possible.

B Successful relaxation of the EW scale with strongly interacting fermions
We now describe in more detail how the conditions listed in Sec. 2 are obtained. Some of these results have been presented in the literature in approximate form, but we give here more complete expressions. To keep the discussion generic, we will use the form of the back-reaction potential shown in Eq. (2.3). Minimizing the potential in Eq. (2.1) we obtain where we have defined the dimensionless field θ = φ/F , v 2 (θ) is the θ-dependent Higgs minimum, and S = sign[µ 2 B v 2 (θ) − Λ 4 0 ]. As usual, the first equation applies when v 2 (θ) is a positive quantity, otherwise the Higgs VEV vanishes. From the minimum equations we see that for small θ we have v 2 (θ) = 0, and S = −1. From the equations of motion it follows that ∂V /∂φ = 0 corresponds to the stopping of the relaxion evolution. To avoid this, we need to make sure that the minimum equation has no solution while v 2 (θ) = 0. This amounts to require which corresponds to Eq. (2.7). As time passes and θ increases, the system reaches a critical value for the relaxion field in which ΛF θ c − µ 2 B cos θ c = 2λΛ 2 and EWSB is triggered. From Eq. (B.1) we see that v 2 (θ) starts growing essentially linearly with θ. Figure 7: Sketch of the electroweak relaxation mechanism. At the beginning of the relaxion evolution the Higgs squared mass parameter is positive and H = 0. As φ evolves it drives H = 0. Once the Higgs VEV is turned on, a back-reaction is generated, stopping the evolution. We also show the sign of the determinant of the matrix of second derivatives of the potential, important in deriving Eq. (B.8).
Looking at Eq. (B.2), and keeping into account Eq. (B.3), we conclude that the right hand side must grow to guarantee the existence of a solution after EWSB. This can happen only if the factor multiplying v 2 (θ) is positive. This requirement amounts to S = +1 (at least around the EW scale), i.e. µ 2 B v 2 EW > Λ 4 0 , and to µ 2 B > ΛF , which are Eq. (2.8).
We now discuss the computation of the EW scale in terms of the parameters of the model. Once the EW minimum is reached for a value θ 0 of the relaxion field we must have 2λv 2 EW = −Λ 2 + ΛF θ 0 + µ 2 B cos θ 0 , r Λ 3 F = µ 2 B sin θ 0 − ΛF v 2 EW − Λ 4 0 sin θ 0 .

(B.4)
Solving these equations for sin θ 0 and cos θ 0 , and using sin 2 θ 0 + cos θ 2 0 = 1 we can determine the value of θ 0 . The two solutions are (B.5) The positive sign gives cos θ 0 < 0, while the negative sign gives cos θ 0 > 0. We now analyze the minimum conditions. Requiring det V > 0, where V is the matrix of second derivative of the potential, see the sketch in Fig. 7, we obtain Since the right hand side is a positive quantity, we immediately conclude that the solution with cos θ 0 < 0 corresponds to a saddle point, while the solution with cos θ 0 > may be a local minimum in some region of parameter space. To determine this region we first translate Eq. (B.6) in a maximum equation for sin θ 0 . We then combine this maximum equation with Eq. (B.9), obtaining Noticing that the term inside the square brackets is very close to one we end up with F Λ 3 µ 2 B v 2 EW − Λ 4 0 as condition to guarantee the existence of a minimum. As for the saddle point solution, we notice that it corresponds to a minimum in the Higgs direction and a local maximum in the relaxion direction, see Fig. 7. By focussing in the φ direction, the difference between the potential in the maximum and in the minimum reads (B.8) Let us now go back to Eq. (B.4). The second equation can be used to compute Requiring | sin θ 0 | ≤ 1 we obtain which is the inequality of Eq. (2.10). In order to maximize the allowed value of Λ we must be close to θ 0 = π/2 + 2πκ, with κ an integer. This agrees with the results of Ref. [97]. Using this in Eq. (B.1) allows us to compute how much the Higgs VEV changes as a function of κ. We obtain In the case κ − κ = 1, i.e. for two subsequent minima, we obtain ∆v 2 π λ where we have used the condition on the parameter space required by the existence of a minimum. We then see that the change of the Higgs VEV between subsequent minima is very small.