Multistep Strongly First Order Phase Transitions from New Fermions at the TeV Scale

TL;DR

This paper explores how a simple fermionic extension of the Standard Model with TeV-scale vector-like leptons can induce multistep strongly first order electroweak phase transitions, producing detectable gravitational wave signals and potential collider signatures.

Contribution

It introduces a minimal fermionic model that generates multistep first order phase transitions near the electroweak scale, linking particle physics with gravitational wave phenomenology.

Findings

Multistep strongly first order phase transitions occur at TeV-scale temperatures.

Distinct gravitational wave signatures are predicted for future space-based detectors.

Potential LHC signatures of vector-like leptons are discussed.

Abstract

In spite of the vast literature on the subject of first order Electroweak Phase Transitions (EWPTs), which can provide the necessary conditions for generating the Baryon Asymmetry in the Universe, fermion-induced EWPTs still remain a rather uncharted territory. In this paper, we consider a simple fermionic extension of the Standard Model involving one $SU(2)_L$ doublet and two $SU(2)_L$ singlet Vector-Like Leptons, strongly coupled to the Higgs boson and with masses close to the TeV scale. We show how such a simple scenario can give rise to a non-trivial thermal history of the Universe, involving strongly first order multistep phase transitions occurring at temperatures close to the electroweak scale. Finally, we investigate the distinct Gravitational Wave (GW) signatures of these phase transitions at future space--based GW detectors, such as LISA, DECIGO, and BBO, and briefly discuss the possible LHC signatures of the Vector-Like Leptons.

Figures

Tables

Table 1. Table 1: Benchmark Points.

	BM1	BM2	BM3
$y_{N_L}$	3.40	3.47	3.47
$y_{N_R}$	3.49	3.45	3.36
$y_{E_L}$	3.34	3.33	2.55
$y_{E_R}$	3.46	3.41	3.28
$m_L$ (TeV)	1.06	1.42	1.43
$m_N$ (TeV)	0.94	0.75	0.83
$m_E$ (TeV)	1.34	1.25	0.72
${\xi }_1$	2.34	2	1.56
${\xi }_2$	1.54	1.35	1.38
${\mu }_{\gamma \gamma }$	1.28	1.20	1.28
$\Delta {\chi }^2(S,T)$	1.33	3.60	4.57
$m_{N_1}$ (GeV)	400	401	466
$m_{E_1}$ (GeV)	592	740	460

Table 2. Table 2: Relevant parameters to the GW spectrum.

	BM1	BM2	BM3
$T_{P{T}_2}$ (GeV)	165	186	164
${\alpha }_2$	0.012	0.010	0.005
${\left(\frac{\beta }{H_{\mathrm{PT}}}\right)}_2$	6480	10880	9690
$T_{P{T}_1}$ (GeV)	82.8	97.7	118
${\alpha }_1$	0.074	0.060	0.016
${\left(\frac{\beta }{H_{\mathrm{PT}}}\right)}_1$	1834	3740	7710

Table 3. Table 3: Predictions for the masses, decay widths, branching ratios, and 13 TeV LHC production cross sections of the lighter VLL eigenstates for the three benchmark points. Here, f NP subscript 𝑓 NP f_{\rm NP} denotes the new fermions N 1 subscript 𝑁 1 N_{1} or E 1 subscript 𝐸 1 E_{1} , while f SM subscript 𝑓 SM f_{\rm SM} stands for a τ 𝜏 \tau or ν 𝜈 \nu .

	BM1	BM2	BM3
$m_{N_1}$ (GeV)	400	401	466
${\Gamma }_{N_1}$  (MeV)	1.32	0.38	0.96
$m_{E_1}$ (GeV)	592	740	460
${\Gamma }_{E_2}$  (GeV)	7.46	11	0.032
$\mathrm{BR}\left(E_1\to N_{1}W\right)$	0.995	0.996	-
$\mathrm{BR}\left(E_1\to \tau h\right)$	$3.5\times {10}^{-3}$	$2.5\times {10}^{-3}$	0.606
$\mathrm{BR}\left(E_1\to \nu W\right)$	$7.6\times {10}^{-4}$	$9.6\times {10}^{-4}$	0.29
$\mathrm{BR}\left(E_1\to \tau Z\right)$	$4.6\times {10}^{-4}$	$4.1\times {10}^{-4}$	0.104
$\sigma \left(pp\to E_1{E}_1\right)$ (fb)	0.32	0.13	0.41
$\sigma \left(pp\to E_1{N}_1\right)\left(fb\right)$	0.36	0.09	0.09
$\sigma \left(pp\to N_1{N}_1\right)\left(fb\right)$	0.31	0.11	0.08
$\sigma \left(pp\to f_{\mathrm{NP}}{f}_{\mathrm{SM}}\right)$ (fb)	$𝒪\left({10}^{-4}\right)$

Equations

$$L_{L,R}=\left(\!\!\!\begin{array}[]{c}N\\ E\end{array}\!\!\!\right)_{L,R}\sim(1,2)_{-1/2},\quad N^{\prime}_{L,R}\sim(1,1)_{0},\quad E^{\prime}_{L,R}\sim(1,1)_{-1},$$

$$-{\mathcal{L}}_{VLL}$$

$$+m_L{\overline{L}}_L{L}_R+m_N{\overline{N}}_L^{\prime }{N}_R^{\prime }+m_E{\overline{E}}_L^{\prime }{E}_R^{\prime }+h.c.,$$

$${\mathcal{M}}_N=\left(\begin{array}{cc}{m}_L& v_h{y}_{N_L}\\ v_h{y}_{N_R}& m_N\end{array}\right),{\mathcal{M}}_E=\left(\begin{array}{cc}{m}_L& v_h{y}_{E_L}\\ v_h{y}_{E_R}& m_E\end{array}\right),$$

$${\omega }_{W,Z}(\varphi )=\frac{m_{W,Z}^2}{v^2}{\varphi }^2,{\omega }_t(\varphi )=\frac{m_t^2}{v^2}{\varphi }^2,$$

$${\omega }_h(\varphi )=\frac{m_h^2}{2v^2}(3{\varphi }^2-v^2),{\omega }_{\chi }(\varphi )=\frac{m_h^2}{2v^2}({\varphi }^2-v^2).$$

$$m_W=80.4\mathrm{GeV},m_Z=91.2\mathrm{GeV},m_t=174\mathrm{GeV},m_h=125\mathrm{GeV}.$$

$$\displaystyle\omega_{X_{1,2}}(\phi)=\frac{1}{2}\Bigg{(}m_{L}^{2}+m_{X}^{2}+\frac{y_{X_{L}}^{2}+y_{X_{R}}^{2}}{2}\phi^{2}$$

$$\displaystyle\mp\sqrt{\left(m_{L}^{2}+m_{X}^{2}+\frac{y_{X_{L}}^{2}+y_{X_{R}}^{2}}{2}\phi^{2}\right)^{2}-\left(2m_{L}m_{X}-y_{X_{L}}y_{X_{R}}\phi^{2}\right)^{2}}\Bigg{)},$$

$$V_0(\varphi )=\frac{m_h^2}{8v^2}({\varphi }^2-v^2{)}^2.$$

$$\frac{{\rm d}\Delta V_{1}^{0}}{{\rm d}\phi}\bigg{|}_{\phi=v}=\frac{{\rm d}^{2}\Delta V_{1}^{0}}{{\rm d}\phi^{2}}\bigg{|}_{\phi=v}=0,$$

$$\Delta V_1^0(\varphi )=\Delta V_{1,\mathrm{SM}}^0(\varphi )+\Delta V_{1,\mathrm{VLL}}^0(\varphi ),$$

$$\Delta V_{1,\mathrm{SM}}^0(\varphi )$$

$$\Delta V_{1,\mathrm{VLL}}^0(\varphi )$$

$$a_i$$

$$b_i$$

$$\omega_{i}\equiv\omega_{i}(v),\qquad\omega_{i}^{\prime}\equiv\frac{{\rm d}\omega(\phi)}{{\rm d}\phi}\bigg{|}_{\phi=v},\quad{\rm and}\quad\omega_{i}^{\prime\prime}\equiv\frac{{\rm d}^{2}\omega(\phi)}{{\rm d}\phi^{2}}\bigg{|}_{\phi=v}.$$

$$n_{\{W,Z,h,\chi ,t,\mathrm{VLL}\}}=\{6,3,1,3,-12,-4\},$$

$$\Delta V_1^T(\varphi ,T)=\frac{T^4}{2{\pi }^2}\left[\sum _{i=W,Z,h,\chi }{n}_i{J}_b\left(\frac{{\omega }_i(\varphi )}{T^2}\right)+\sum _{i=t,N_{1,2},E_{1,2}}{n}_i{J}_f\left(\frac{{\omega }_i(\varphi )}{T^2}\right)\right].$$

$$J_{b,f}(x^2)={\int }_0^{\infty }\mathrm{d}k{k}^2\log \left[1\mp {\mathrm{e}}^{-\sqrt{k^2+x^2}}\right],$$

$$\Delta V_{\mathrm{D}}(\varphi ,T)=\frac{T}{12\pi }\sum _{i=h,\chi ,W,Z,\gamma }{\bar{n}}_i\left[{\omega }_i^{3/2}(\varphi )-{\left({\omega }_i(\varphi )+{\Pi }_i(T)\right)}^{3/2}\right],$$

$${\Pi }_{h,\chi }(T)$$

$${\Pi }_W(T)$$

$${\Pi }_Z(T)$$

$${\Pi }_{\gamma }(T)$$

$$V(\varphi ,T)=V_0(\varphi )+\Delta V_{1,\mathrm{SM}}^0(\varphi )+\Delta V_{1,\mathrm{VLL}}^0(\varphi )+\Delta V_1^T(\varphi ,T)+\Delta V_D(\varphi ,T),$$

$$m_L,m_N,m_E\in [500,1500]\text{leavevmode}\mathrm{GeV},y_{N_{L,R}},y_{E_{L,R}}\in [2,\sqrt{4\pi }].$$

$$V({\varphi }_c,T_c)=0\mathrm{and}{V}^{\prime }({\varphi }_c,T_c)=0\mathrm{with}{\varphi }_c\ne 0.$$

$$\alpha \sim 0.01-0.1,\beta /H_{\mathrm{PT}}\sim 10^3-10^4.$$

$$S_3(T)=4\pi \int dr{r}^2\left[\frac{1}{2}{\left(\frac{d{\varphi }_b}{dr}\right)}^2+V({\varphi }_b(r),T)\right],$$

$$\frac{d^{2}\phi_{b}}{dr^{2}}+\frac{2}{r}\frac{d\phi_{b}}{d_{r}}=V^{\prime}(\phi,T),\quad{\rm with}\quad\phi_{b}(r=\infty)=\phi_{\rm true}\quad{\rm and}\quad\frac{d\phi_{b}}{dr}\bigg{|}_{r=0}=0,$$

$$\frac{S_3(T_{\mathrm{PT}})}{T_{\mathrm{PT}}}=142.$$

$$\alpha =\frac{\left[V({\varphi }_{\mathrm{false}},T_{\mathrm{PT}})-V({\varphi }_{\mathrm{true}},T_{\mathrm{PT}})\right]+T_{\mathrm{PT}}{\left[{\partial }_{T}V({\varphi }_{\mathrm{true}},T)-{\partial }_{T}V({\varphi }_{\mathrm{false}},T)\right]}_{T_{\mathrm{PT}}}}{{\rho }_{\mathrm{rad},\mathrm{true}}(T_{\mathrm{PT}})},$$

$$\displaystyle\frac{\beta}{H_{\rm PT}}=T_{\rm PT}\frac{d}{dT}\left(\frac{S_{3}(T)}{T}\right)\bigg{|}_{T_{\rm PT}},$$

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

Full text

1. 1 Introduction

2. 2 New Dirac Fermions and the Phase Structure of the Universe

3. 2.1 A Minimal Vector-Like Lepton Model for Strong Phase Transitions

4. 2.2 Phase Transition Calculation

5. 2.3 General Predictions of the Model

6. 3 Gravitational wave and collider signatures

7. 3.1 GW signature

8. 3.2 LHC signatures

9. 4 Summary and conclusions

10. A Calculation of Gravitational Wave Spectrum

11. B SM–VLL Mixing for Collider Phenomenology

12. C Theoretical Constraints

1 Introduction

The origin of the baryon asymmetry in the Universe (BAU), which is also known as the matter-antimatter asymmetry, is one of the most outstanding problems in modern cosmology. Baryogenesis is an appealing scenario of generating the matter-antimatter asymmetry dynamically. In 1967, Sakharov enunciated three conditions for a successful baryogenesis: baryon number violation, $C$ and $CP$ violation, and a departure from the thermal equilibrium [1].

Of particular interest for this work is the third condition, departure from thermal equilibrium, which can only be met if the Universe underwent a strongly first order phase transition (SFOPT) in its early stages. While there is no indication about the energy scale at which it happened, an attractive possibility would be that such a phase transition (PT) occurred around the Electroweak (EW) scale. Indeed, it is beyond any doubt that EW symmetry was broken at some point in the history of the Universe, and a straightforward possibility would be that the PT responsible for baryogenesis took place when EW symmetry was broken.

As any cosmological PT, an Electroweak phase transition (EWPT) would generate a stochastic GW background sourced in the collision of nucleating bubbles and the plasma motion induced by bubble collisions. In the specific case of a strong EWPT, it is expected that, due to redshift, the GW signal would nowadays peak at frequencies around the mHz-dHz range [2]. Interestingly, this frequency range overlaps with the frequency ranges which future space–based interferometers, such as LISA [2], DECIGO [3], and BBO [4], will be most sensitive to. Discovering such a gravitational wave signal would establish the EWPT as a new milestone in our understanding of the early universe.

In the Standard Model (SM), electroweak symmetry breaking (EWSB) would proceed via a smooth crossover unless the Higgs mass is below $\sim 70$ GeV [5, 6]. Therefore, the discovery of the SM Higgs boson with a mass $m_{h}=125$ GeV [7, 8] meant that the SM alone cannot satisfy the third Sakharov condition, i.e. departure from thermal equilibrium. 111Moreover, the amount of CP violation in the SM is too small for satisfying the second Sakharov condition, which means that new sources of CP violation should be considered. Addressing this problem is, however, beyond the scope of this work.

Consequently, the problem of strongly first order EWPTs was considered in several beyond the Standard Model (BSM) scenarios, most of them containing extra scalar states with respect to the SM. Examples include scalar singlet extensions [9, 10, 11, 12, 13, 14, 15, 16, 17, 18], Two Higgs Doublet Models (2HDMs) [19, 20, 21, 22, 23, 24], the (Next to) Minimal Supersymmetric Standard Model 222In the MSSM, the 1st order EWPT is excluded by the measured Higgs mass, $m_{h}=125$ GeV, and by light stop searches [25]. (MSSM/NMSSM) [26, 27], or Composite Higgs Models [28, 29, 30, 31, 32, 33]. Scalars are of particular interest as, when integrating out heavy scalars, a $(H^{\dagger}H)^{3}$ term can be generated and therefore induce a barrier in the Higgs potential [11, 34]. Also, scalars contribute to a negative cubic term in the Higgs effective potential in the high temperature expansion and therefore can induce a barrier via thermal effects [35]. There has been much less focus on fermionic extensions in the context of EW baryogenesis [36, 37, 38, 39, 40, 41] as fermions do not contribute as much in the high temperature expansion. However, when the critical temperature is low compared to the new fermion masses, fermions contribute to the Higgs effective potential in the same way as scalars, and therefore can lead to non-trivial effects in the thermal history of the early universe. For example, Ref. [38] considered two fermionic multiplets: a Majorana singlet and a vector–like (Dirac) $SU(2)_{L}$ doublet with the same quantum numbers as the SM leptonic doublet. In this setup, a strong first order EWPT can be induced by the neutral fermions, which couple strongly to the Higgs.

In this paper, we study a SM extension containing only Dirac fermions, and investigate the impact of these fermions on the thermal history of the Universe. We choose a BSM spectrum containing only vector–like leptons (VLLs), as current limits from the LHC push vector–like quarks (VLQs) at masses above $1$ TeV, making them too heavy to considerably influence the EWPT. We show that such a model can indeed accommodate a SFOPT capable of generating the BAU as long as the new fermions couple strongly to the Higgs. Interestingly, we find that, although having only one scalar field – the SM Higgs, our model predicts a three–step phase transition, consisting of one smooth crossover at high temperatures and two SFOPTs at temperatures of the order of the EW scale. Moreover, we calculate the GW signal and collider impact of our model, comparing the signals with present and/or future collider (LHC) and GW searches.

This paper is structured as follows. In Sec. 2, we construct a model containing new Vector-Like Leptons (VLLs), show how it can accommodate SFOPTs, and discuss the general predictions of this model. Sec. 3 is dedicated to a more detailed analysis of three benchmark points, for which we calculate the GW signature and several collider observables, such as the VLL production cross sections and branching ratios. Finally, we summarize and present our conclusions in Sec. 4.

2 New Dirac Fermions and the Phase Structure of the Universe

In this section, we build a minimal Dirac VLL model that can produce SFOPTs in the Early Universe. To calculate the PT strength in our model, we construct its 1-loop effective potential, comprising both zero– and finite–temperature pieces, and then include the so-called Daisy resummation contribution. We then discuss the thermal behavior and current constraints of this model.

2.1 A Minimal Vector-Like Lepton Model for Strong Phase Transitions

Let us start by outlining the logical steps followed when building our model. Firstly, we require that the new fermions couple strongly to the Higgs field. In the limit of null Yukawa couplings for the new fermions, adding them would leave the 1-loop effective potential unchanged with respect to the SM-only case. In order to dramatically change the thermal evolution of the 1-loop effective potential, which exhibits only a crossover in the SM, one should therefore consider strong Yukawas for the VLLs.

The simplest VLL model would correspond to adding a single VLL multiplet to the SM. In such a case, the only possible non-SM Yukawa terms would couple a SM lepton and the VLL multiplet to the Higgs doublet. However, such a term would mix the VLL and the SM lepton (which we assume to be the $\tau$ lepton, to avoid stronger constraints from electrons and muons). As explained in the previous paragraph, the new Yukawa coupling has to be large, which would result in a strong mixing between the VLL and $\tau$ and hence a significant departure from the SM prediction of the $\tau$ couplings. This forces us to discard such a scenario, as the $\tau$ couplings are tightly constrained to be SM-like by experimental measurements such as $Z\to\tau\tau$ decays at LEP  [42] or $h\to\tau\tau$ decays at LHC [43, 44]. Therefore, throughout this section, we neglect the mixing between VLLs and the SM fermions 333We will come back to this point in Sec. 3, as the corresponding Yukawa couplings would have an insignificant effect on the phase structure of the Universe.

The next logical choice would be to augment the SM with one VLL doublet and one VLL singlet, since this configuration would allow for a Yukawa term coupling the two VLL multiplets to the Higgs doublet. However, such a model with strong Yukawas would badly violate custodial symmetry, giving rise to unacceptable contributions to the $T$ parameter [45, 46, 47]. As we have checked, one cannot accommodate SFOPTs in this scenario without dramatically exceeding the experimental bounds on the $T$ parameter.

Therefore, the minimal solution is to add one VLL doublet and two VLL singlets, since such a configuration can accommodate an (approximate) custodial symmetry, which allows for large Yukawas while avoiding significant contributions to the $T$ parameter. We choose the new leptons to have similar $SU(3)_{c}\times SU(2)_{L}\times U(1)_{Y}$ quantum numbers as their SM counterparts:

[TABLE]

where $L,R$ stand for the VLL chiralities. The new fermions $N^{(\prime)}_{L,R}$ and $E^{(\prime)}_{L,R}$ have zero and $-1$ electric charge, respectively, so we shall refer to them as neutral VLLs or VL neutrinos and charged VLLs or VL electrons. For denoting the multiplets in the equation above, we use the standard $(SU(3)_{c},SU(2)_{L})_{Y}$ notation, where the hypercharge $Y$ is given by the difference between the electric charge and the third isospin componen, i.e. $Y=Q-T_{3}$.

As stated previously, we neglect for the time being the mixing between the SM leptons and the VLLs. Therefore, the most general renormalizable VLL Yukawa Lagrangian, consistent with $SU(3)_{c}\times SU(2)_{L}\times U(1)_{Y}$ gauge symmetry, reads

[TABLE]

where $H$ represents the SM Higgs doublet and $\tilde{H}$ its $SU(2)$ conjugate, $y$’s the dimensionless Yukawa couplings, and $m$’s the vector-like masses. Upon Electroweak Symmetry Breaking (EWSB), one can write the neutral and charged mass matrices, respectively, as

[TABLE]

with $v_{h}\equiv v/\sqrt{2}\simeq 174$ GeV, where $v\simeq 246$ GeV is the Higgs vacuum expectation value (vev). The physical masses, which we denote and order as $m_{N_{1}}<m_{N_{2}}$ and $m_{E_{1}}<m_{E_{2}}$, are obtained as usual by bi-diagonalizing the mass matrices from the above equation. The physical couplings of the VLL eigenstates can be calculated from the corresponding rotation matrices that bi-diagonalize the mass matrices.

2.2 Phase Transition Calculation

To study the thermal history of the Universe, we first write down the 1-loop effective scalar potential, taking into account the effect of SM particles strongly coupled to the Higgs ($W$ and $Z$ bosons, $t$ quark, $h$ boson, and Goldstone bosons, $\chi$) and of the VLLs. We denote the background field-dependent squared mass as $\omega_{i}(\phi)$, where $i$ labels the particles and $\phi$ is the background field value. For the SM fields coupling strongly to the Higgs, the various $\omega$’s are

[TABLE]

Throughout this work, we use the following values for the masses of the SM particles [48]:

[TABLE]

The $\phi$-dependent VLL squared eigenmasses are obtained by diagonalizing $\mathcal{M}_{X}^{\dagger}\mathcal{M}_{X}$, with $X=N,E$, and $\mathcal{M}_{X}$ defined in Eq. (3). We thus have:

[TABLE]

where “$-$” corresponds to $X_{1}$ (the lighter eigenstate) and “$+$” to $X_{2}$ (the heavier eigenstate). To make the connection with previous notations, we remind the reader that $\omega_{X_{j}}(v)=m_{X_{j}}^{2}$, with $j=1,2$.

The Effective Potential

We now proceed to calculating the 1-loop effective potential, comprised of the SM tree level part and the zero and finite temperature 1–loop contributions (for both the SM particles and the VLLs), to which we add the so-called Daisy contribution. The Daisy (or ring) contribution is a finite temperature effect coming from the resummation of higher–loop IR-divergent diagrams of a certain topology [49], whose sum amounts to a finite result.

The SM tree level contribution is given by

[TABLE]

The 1-loop zero-temperature contribution. We work in the on-shell renormalization scheme for the 1-loop contribution,

[TABLE]

where $\Delta V_{1}^{0}(\phi)$ is the (finite) Coleman-Weinberg potential, which includes the counter-terms and the zero-temperature piece of the 1-loop correction to the potential. We further split $\Delta V_{1}^{0}(\phi)$ into two pieces,

[TABLE]

$\Delta V_{\rm 1,\,SM}^{0}$ comprising the effect of SM particles and $\Delta V_{\rm 1,\,VLL}^{0}$ capturing the VLL contribution (at $T=0$). In the Landau gauge, which we are going to use throughout this work, these two contributions read 444Note that, since $\omega_{\chi}(v)=0$, the Goldstone contribution at 1–loop diverges due to the logarithmic term. This IR divergence can be cured by imposing that the second derivative of the total potential (tree plus 1–loop contributions) at $\phi=v$ is equal to the 1–loop Higgs mass evaluated at $p^{2}=m_{h}^{2}$ (and not $p^{2}=0$), where $p$ is the external momentum in the Higgs self-energy diagram. This amounts to replacing $\omega_{\chi}(v)$ with $\omega_{h}(v)$ [50] in the logarithmic term corresponding to $i=\chi$ from Eq. (10). For a more detailed discussion of the matter, see Refs. [51, 52].

[TABLE]

where the coefficients $a_{i}$ and $b_{i}$, with $i=N_{1},N_{2},E_{1},E_{2}$, are given by [36, 38]

[TABLE]

In Eqs. (11)–(13), $\mu$ denotes an arbitrary energy scale which is introduced to make the logarithm arguments adimensional. Moreover, we used the following notations:

[TABLE]

Finally, the number of degrees of freedom for the fields running in the loops are

[TABLE]

with $n_{\rm VLL}=n_{N_{1,2}}=n_{E_{1,2}}=-4$, since the VLLs we introduce are Dirac fermions.

The 1-loop finite-temperature contribution is obtained, in the imaginary time formalism, by performing a Wick rotation and compactifying the resulting Euclidean time dimension on a circle of radius $R=(2\pi T)^{-1}$, with $T$ denoting the temperature. The fields are thus Fourier expanded along the periodic time dimension, with eigenfrequencies given by the Matsubara frequencies, which are given by $2n\pi T$ for bosons (periodic on the time circle) and $(2n+1)\pi T$ for fermions (anti–periodic on the time circle). Performing the infinite sum on these frequencies gives rise to the finite–$T$ contribution to the potential, which reads

[TABLE]

Here, $J_{b}$ and $J_{f}$ are adimensional integrals accounting for the thermal contribution of boson and fermions, respectively, and are given by (see e.g. Ref. [53])

[TABLE]

where the “$-$” (“$+$”) sign applies to bosons (fermions).

The Daisy contribution. At finite temperature, the perturbative expansion in terms of a small coupling breaks down due to IR divergences coming from thermal loops involving bosonic 0–modes, which have vanishing Matsubara frequencies. This problem can be fixed by performing a resummation of the a certain class of diagrams, the so–called Ring or Daisy diagrams, which are $N$–loop diagrams in which $N-1$ loops are each attached to the main one through one and only one 4–point vertex. While all these diagrams are IR–divergent if taken one by one, their sum adds up to a finite result, which corresponds to the following contribution to the effective potential:

[TABLE]

with $\bar{n}_{\{h,\chi,W,Z,\gamma\}}=\{1,3,2,1,1\}$ represent either scalar ($h,\chi$) or gauge boson longitudinal degrees of freedom, since only these degrees of freedom acquire a thermal mass (transverse gauge modes are protected by gauge symmetry). Note that fermionic diagrams do not produce any IR divergences, as their 0–modes have non-vanishing Matsubara frequencies. In the equation above, the $\Pi_{i}$’s represent the thermal or Debye masses of the fields appearing in the sum  [52],

[TABLE]

In the expressions above, we use the high–temperature approximation for the contribution of the SM particles, which is justified by the fact that the Daisy diagrams give a negligible contribution for $\omega_{i}(\phi)\gtrsim T^{2}$. Moreover, because the strong phase transitions in our model occur at temperatures well below the VLL masses, we neglect the VLL contribution to the $\Pi_{i}$’s in Eqs. (19)–(22).

Adding all the pieces together, the effective potential we use for calculating the strength of the phase transition is given by

[TABLE]

with all the contributions detailed in Eqs. (7), (10), (11), (16), and (18). Finally, since only potential differences have a physical impact in our analysis, we shift the potential by a constant such that $V(\phi=0,T)=0$ for every $T$.

Scan for the PT strength calculation

We now calculate the strength of the phase transition in our model. We scan over the following range of parameter values:

[TABLE]

In our initial scans, we allowed for wider ranges, and found out that parameter values inside the ranges shown above are more likely to lead to strong phase transitions. Also, as noted in Ref. [38], having $y_{N_{R}}y_{N_{L}}>0$ and $y_{E_{R}}y_{E_{L}}>0$ favors SFOPTs, which is why we chose all the Yukawas to be positive.

After each point in parameter space is generated, we check whether the said point is in agreement with experimental constraints. Firstly, as the VLLs under consideration have $SU(2)_{L}$ quantum numbers, they affect the electroweak gauge boson self-energies, so one constraint is the contribution to the oblique parameters $S$ and $T$ [46, 45], for which we use the $2\sigma$ values quoted in Ref. [54]. Secondly, the charged VLLs change the loop-induced $h\gamma\gamma$ coupling with respect to its SM value. This coupling is probed at the LHC through the diphoton Higgs signal strength, $\mu_{\gamma\gamma}\equiv\Gamma(h\to\gamma\gamma)/\Gamma(h\to\gamma\gamma)_{\rm SM}$. As the experimental bound for this observable, we use the 2$\sigma$ interval quoted by the ATLAS collaboration, $0.71<\mu_{\gamma\gamma}<1.29$ [55]. 555A subsequent CMS measurement shows a slightly higher value for the $h\to\gamma\gamma$ signal strength, $\mu_{\gamma\gamma}=1.18^{+0.17}_{-0.12}$ at 1$\sigma$ [56]. Thirdly, we impose a lower bound on the masses of the lighter eigenstates, $m_{E_{1}}>100$ GeV and $m_{N_{1}}>90$ GeV [42].

From the theoretical point of view, we also impose a lower limit on the depth of the EW minimum at the lower minimum, $\left|V(\phi=v,T=0)\right|$ (we remind the reader that, by convention, $V(0,T)=0$). As illustrated in Ref. [57], the lower the depth of the present day EW minimum, the more delayed and thus the stronger the phase transition is. For our analysis, we choose $\left|V(\phi=v,T=0)\right|<8.5\times 10^{7}$ GeV4, a value for which we checked explicitly that most of the surviving points exhibit SFOPTs.

For the points surviving the constraints listed above, we calculate the phase transition strength (or order parameter), which is defined as $\xi\equiv\phi_{c}/T_{c}$, where $\phi_{c}$ and $T_{c}$ are the critical field value and critical temperature, respectively. $T_{c}$ is defined as the temperature at which the values of the potential at the minima located at $\phi=0$ (“symmetric minimum”) at $\phi\neq 0$ (“broken minimum”) become degenerate. $\phi_{c}$ is defined as the position along the field axis of the broken minimum for $T=T_{c}$:

[TABLE]

In order to find the critical field value and temperature for a given point in the parameter space spanned by the VLL Yukawas and diagonal mass terms, we numerically solve the system of equations Eq. (25), whose solution is given by $\phi_{c}$ and $T_{c}$ from which we calculate the phase transition strength, $\xi=\phi_{c}/T_{c}$. In the above equations, the prime symbol denotes differentiation with respect to $\phi$.

To speed up computations, without spoiling the reliability of our results, we calculate $\phi_{c}$ and $T_{c}$ for the 1-loop effective potential without the subleading scalar (Higgs+Goldstone) and Daisy contributions. We checked that adding taking into account these contributions amounts to a $\mathcal{O}(5-10\%)$ correction to the values of $\xi$, which justifies a posteriori our approximation. 666Nevertheless, when studying the benchmark points in Sec. 3, we are going to take into account the scalar and Daisy contributions too. Therefore, in the following, unless mentioned otherwise, all the scatter plots resulting from our scan will not feature the scalar and Daisy contributions.

2.3 General Predictions of the Model

We now present some general predictions of our model. The most striking one concerns the non-trivial thermal evolution of the potential from Eq. (23). For the scan points that pass the constraints, three distinct phase transitions are predicted in the Early Universe: one crossover and two SFOPTs. As an illustration of this fact, we show in Fig. 1 the behavior of the potential for the benchmark point BM1 presented in Sec. 3. In Fig. 1, we plot the potential as a function of the background field value, $\phi$, for six different temperatures. In its early stages, the Universe starts in the symmetric phase, i.e. at $\phi=0$, and then, as it expands and cools down, has a crossover to the broken phase, $\phi\neq 0$, which results in EWSB. In the following, we shall refer to the broken phase also as the EW or broken minimum. As the temperature lowers, the EW minimum becomes less and less deep, and a potential barrier starts developing between the symmetric and broken phases. At a critical temperature, which we denote by $T_{c_{2}}$, the two minima become degenerate, and the Universe starts tunneling back to the symmetric phase, thus undergoing a SFOPT and restoring EW symmetry. We denote the critical field value for this first SFOPT as $\phi_{c_{2}}$, and the corresponding PT strength as $\xi_{2}=\phi_{c_{2}}/T_{c_{2}}$. After this SFOPT, as the Universe cools down, the EW minimum becomes a local minimum, and for this particular point in parameter space even disappears altogether. However, this trend reverses at even lower temperatures, and the EW minimum becomes again degenerate with the symmetric one, and EW symmetry breaking is triggered once again. This happens for a critical temperature $T_{c_{1}}$ and a corresponding critical field value $\phi_{c_{1}}$, with the PT strength given by $\xi_{1}=\phi_{c_{1}}/T_{c_{1}}$. For temperatures below $T_{c_{1}}$, the Universe remains in the EW minimum.

We now discuss the non–trivial behavior of the potential, as shown in Fig. 1, focusing on the contributions determining the three PTs – the crossover and the two SFOPT. Firstly, the symmetry–breaking crossover occurring at high $T\gtrsim 500-600$ GeV is induced through thermal effects: the light VLL ($N_{1}$ and $E_{1}$) thermal contribution, which favors the broken phase, competes with the SM and heavy VLL ($N_{2}$ and $E_{2}$) thermal pieces, which tend to keep the Universe in the symmetric phase. The crossover takes place when the first contribution overtakes the other one. Secondly, the earlier, symmetry–restoring SFOPT, which typically takes place at $T\sim 150-200$ GeV, is sourced by both finite and zero temperature effects. In this case, the SM thermal part and the $T=0$ VLL contribution tend to restore EW symmetry, whereas the light VLL thermal part and the SM $T=0$ contribution work towards keeping the Universe in the broken minimum. The net effect is a barrier between the two minima, and the symmetry–restoring PT is thus strongly first–order. Concerning the heavy VLLs, they are decoupled from the thermal bath at this stage and play a negligible role in the dynamics of the earlier SFOPT. Finally, the later SFOPT, which occurs at temperatures around $100$ GeV and is responsible with generating the BAU, is again the result of an interplay between zero– and finite–T effects. On one side, the SM thermal piece and the $T=0$ VLL contribution favor the Universe to remain in the symmetric phase, while the SM $T=0$ part favors EWSB. Adding together these competing effects, a barrier is again generated between the two minima, and EW symmetry is broken via a SFOPT. The thermal impact of the new fermions is negligible for this PT, as both light and heavy VLLs are decoupled from the plasma when this PT occurs.

Among the two SFOPTs in this model, only the later one is responsible for generating the BAU. During the later one, the universe undergoes a phase transition from the symmetric phase to the broken phase. This leads to the formation and expansion of bubbles of the true vacuum in the false, symmetric vacuum. In the presence of the $CP$ violation, particle interactions with the expanding bubbles can lead to the creation of an excess of baryons inside the bubbles through baryon number violating processes induced by sphalerons [58]. The sphaleron process, when in equilibrium, wipes off the created excess of baryons. The sphaleron rate in the broken phase is proportional to ${\rm e}^{-\sqrt{\omega_{W}(\phi)}/T}$, and is suppressed if the phase transition is of strong first order, $\xi$, is greater than $1.3$ [59]. In our case, this condition translates to $\xi_{1}\gtrsim 1.3$, which, as we are going to see in the following, is satisfied by most of the surviving points from our scan. The earlier phase transition, however, can not generate baryon asymmetry. The sphaleron rate is unsuppressed and behaves as $T^{4}$ [59] in the symmetric phase, so any matter anti–matter asymmetry thus generated would be wiped out by the high sphaleron rate in the symmetric phase, to which the Universe tunnels during the earlier SFOPT.

We now display some scatter plots resulting from our scan. As mentioned previously, the points shown in these scatter plots correspond to a computation done without incorporating the scalar and Daisy contributions to the effective potential.

Firstly, we display in Fig. 2 the correlation between $\xi_{1}$ and $\xi_{2}$, the strengths of the two SFOPTS. A feature implied by this plot is that, in our model, the later PT is in general stronger than the earlier one, i.e. $\xi_{1}>\xi_{2}$. In simple terms, this can be explained by noting that the critical field values are rather close, $\phi_{c_{1}}\simeq\phi_{c_{2}}\simeq\mathcal{O}(v)$ (see the example in Fig. 1), whereas there is a larger separation between the two critical temperatures, $T_{c_{2}}>T_{c_{1}}$. Furthermore, as we shall see in Sec. 3, the later SFOPT generally produces a stronger Gravitational Wave (GW) signal than the earlier one, as $\xi_{1}>\xi_{2}$.

Secondly, we plot in Fig. 3 the masses of the lighter VLL eigenstates, $m_{N_{1}}$ and $m_{E_{1}}$, versus the corresponding values of $\xi_{1}$. We find that, when applying the constraints mentioned in Sec. 2.2, the lightest new particles predicted by our model lie in the hundreds of GeV range, with masses typically comprised between $250$ GeV and $900$ GeV. We also notice a correlation between the masses of $N_{1}$, $E_{1}$ and the strength of the later PT: the larger is $\xi_{1}$, the lower are the upper bounds of $m_{N_{1}}$ and $m_{E_{1}}$. This correlation can be understood by using a decoupling argument: the more massive the new fermions become, the more we expect their influence on EW scale physics to diminish. It is also worth noting that, in the right panel of Fig. 3, as $\xi_{1}$ increases, the lower bound on the mass of $E_{1}$ increases. This behavior is induced by the $\mu_{\gamma\gamma}$ constraint: larger values of $\xi_{1}$ need larger values of Yukawa couplings for the charged VLLs, which, in order to satisfy the limit on $\mu_{\gamma\gamma}$, need to be compensated by a larger $m_{E_{1}}$.

Finally, we show in Fig. 4 the correlation between the Higgs-diphoton signal strength, $\mu_{\gamma\gamma}$, and $\xi_{1}$. The diphoton signal strength is influenced in our model by the charged VLLs running in the $h\to\gamma\gamma$ loop. We observe a clear trend in the figure: the stronger the later SFOPT is, the more $\mu_{\gamma\gamma}$ deviates from its SM value, $\mu_{\gamma\gamma}^{\rm SM}=1$. Moreover, the deviation of the diphoton signal strength is always positive, our model predicting an increase of order 15–30$\%$ in $\mu_{\gamma\gamma}$. This correlation can be understood as follows: in order to obtain strong PTs, the new VLLs – both neutral and charged – need to be not too heavy and to have large Yukawa couplings to the SM Higgs. In turn, relatively light charged VLLs which couple strongly to the Higgs give a sizeable contribution to the $h\gamma\gamma$ loop, hence the deviation in $\mu_{\gamma\gamma}$. The positive interference between the charged VLL and SM contributions comes from the fact that strong SFOPTs favor same–sign Yukawa couplings in the charged mass matrix written in Eq. (3). We also mention that $\xi_{1}$ and the $S,T$ parameters are uncorrelated.

Before moving to the next section, we would like to stress the importance of $\mu_{\gamma\gamma}$ in testing our model. As explained in the previous paragraph, for our model, achieving a FOPT strong enough to accommodate the generation of the BAU entails a deviation of more than 10$\%$ of the loop-induced $h\gamma\gamma$ coupling from its SM value. The $\mu_{\gamma\gamma}$ observable is currently measured by the two LHC collaborations, ATLAS and CMS, and the present error is of order $\sim 15\%$ at $1\sigma$ [55, 60]. However, in the future high–luminosity phase of LHC, HL–LHC, the error $\mu_{\gamma\gamma}$ is expected to reduce to $\sim 5\%$ [61, 62, 63]. Therefore, we conclude that the HL–LHC will be able to fully test our model of VLL–induced EW baryogenesis.

3 Gravitational wave and collider signatures

In this section, we perform an in-depth analysis of three benchmark points selected from our model. The first benchmark point corresponds to the strongest later first-order PT in the parameter space. The second benchmark point has a relatively lower value of $\mu_{\gamma\gamma}$, but at the same time accommodating two SFOPTs. The third benchmark point corresponds to the lowest value (among the points retained from our scan) of the mass of the lightest charged eigenstate, $m_{E_{1}}$. For all three benchmark points, the later PT, which is the one responsible for generating the BAU, is strongly first order, $\xi_{1}>1.3$. Moreover, the Higgs diphoton signal strength deviates from unity, as expected from the high Yukawas, but is still within the $2\sigma$ interval quoted by the most recent ATLAS measurement [55].

Furthermore, the masses of the lightest VLLs are in the $400$–$800$ GeV range, and the (correlated) values of the $S$ and $T$ parameters are well within the $2\sigma$ range. The benchmark scenarios are listed in Table 1. In addition to the observables discussed in Sec. 2.3, we calculate the gravitational wave (GW) spectrum and several collider observables (VLL branching ratios and production cross sections) corresponding to these benchmarks and comment on the detectability of such signals at future GW experiments and the LHC. We remind that in this section we take into account the scalar and Daisy contributions as well when calculating the effective potential.

3.1 GW signature

A first order EWPT is expected to generate a stochastic background of gravitational waves [64, 65]. In a radiation–dominated Universe, there are three sources of GW production at a SFOPT: bubble collisions, in which the localized energy density generates a quadrupole contribution to the stress-energy tensor, which in turn gives rise to gravitational waves, plus sound waves in the plasma and magnetohydrodynamic turbulence. The latter two are generated after the bubbles have collided.

There are two key parameters that determine the spectrum of the stochastic GW background generated during a SFOPT. The first one, denoted by $\alpha$, represents the ratio between the latent heat released during the PT and the radiation energy density at the temperature at which the PT occurs. The $\alpha$ parameter plays a role in setting the strength of the GW signal: the higher is $\alpha$, the stronger is the predicted GW stochastic background. The second one, commonly referred to as $\beta$, sets the inverse timescale associated with the PT duration. $\beta$ influences both the strength of the GW signal, as well as the frequency at which the GW peaks. A higher value of $\beta$ implies a weaker GW signal and a shifting towards higher values of the peak frequency of the signal. For details regarding our computation of the GW spectrum we refer the reader to the Appendix.

For the points of our scan that accommodate a strong enough first–order later PT, the typical values of these parameters are

[TABLE]

By comparison, the more popular singlet scalar models feature higher values of $\alpha$ and lower values of $\beta$ [66, 67], which is why our VLL model predicts weaker GW signals than the singlet extension of the SM. Moreover, due to the high values of $\beta/H_{\rm PT}$, the typical GW signal of our VLL model peaks at frequencies in the $0.01-1$ Hz range. This can be understood immediately from Eq. (A.10), which states that the peak frequency depends linearly on $\beta/H_{\rm PT}$. An interesting prediction of the VLL model under study is the multi–peaked GW signature [68, 69, 70], which is generated by the two SFOPTs featured for the points selected by our scans. We find that the GW signature coming from the earlier SFOPT is weaker and peaks at higher frequencies than the one resulting from the later SFOPT. These features are a consequence of having $\alpha_{2}<\alpha_{1}$ and $\left(\beta/H_{\rm PT}\right)_{2}>\left(\beta/H_{\rm PT}\right)_{1}$.

The parameters relevant to the GW spectrum are listed in Table. 2. We show in Fig. 5 - Fig. 7 the spectrum $h^{2}\Omega_{\rm GW}$, as a function of the frequency $f$, of the GWs released during the two SFOPTs predicted for the three benchmarks listed in Table 1. The solid black line corresponds to the later PT, while the dashed one stands for the earlier PT. Owing to the fact that $\beta_{1}<\beta_{2}$ and $\alpha_{1}>\alpha_{2}$, the GW signal from the later (and stronger) PT is generally 2–3 orders of magnitude stronger than the signal from the earlier PT. The GW1 signal, corresponding to the later PT, peaks at frequencies of $\sim 0.05$ Hz, $\sim 0.1$ Hz, and $\sim 0.3$ Hz for the three benchmark points respectively, while the other signal’s peak is located at $\sim 0.4$ Hz, $\sim 0.7$  Hz, and $\sim 0.5$  Hz for the three benchmark points respectively. The colored regions represent the future sensitivity of the LISA experiment for four possible configurations, C1–C4 [2], and of the DECIGO [3] and BBO [4] experiments. The sensitivity curves are taken from Ref. [2] for LISA and from Ref.  [71] for DECIGO and BBO. We observe that, even if the GW1 signal has a strength comparable to the sensitivity of the C1 configuration of LISA, it peaks at a frequency which does not correspond to the maximum LISA sensitivity, which lies in the mHz range. Therefore, we conclude that this GW signal would not be detectable by the LISA experiment. However, the stronger GW1 signals, such as the ones from benchmarks BM1 and BM2, can be detected by the DECIGO and BBO experiments, whose projected sensitivities are maximized close to the peak frequency of GW1. On the other hand, the weaker GW2 signal would escape detection at all three GW experiments we consider in this work.

3.2 LHC signatures

The main VLL production mechanism at the LHC is pair production. The total VLL pair production cross sections for our benchmarks are of $\mathcal{O}(0.1)$ - $\mathcal{O}(0.4)$ fb at the LHC, whereas the SM fermion + VLL ones are severely suppressed by the small VLL–$\tau$ mixing (see Appendix). For the first two benchmarks, because of $m_{E_{1}}>m_{N_{1}}+m_{W}$, the $E_{1}\to N_{1}W$ decay channel is open. As it is not suppressed by SM–VLL mixing, it is the dominant decay mode of $E_{1}$. This also explains the sizeable difference between the widths of $E_{1}$ and $N_{1}$. For the third benchmark, the approximate degeneracy of $E_{1}$ and $N_{1}$ drastically reduces the $N_{1}\to E_{1}W$ branching ratio. As a result, $E_{1}$ has a width of $\sim 30$ MeV (instead of $\mathcal{O}(10)$ GeV in the previous two examples), its most probable decay channel being $E_{1}\to\tau h$, with subleading $\nu W$ and $\tau Z$ branching ratios.

These findings, alongside with other collider predictions, are collected in Table 3 for all three benchmarks. In summary, the cross section for VLL pair production at a 13 TeV proton–proton collider are below the fb level, rendering direct searches at the LHC challenging. Moreover, single VLL production (in association with a SM lepton) is below the ab level, which means that such a process is undetectable at the LHC. However, as explained at the end of Sec. 2.3, a more promising search avenue is measuring the Higgs to diphoton signal strength, $\mu_{\gamma\gamma}$, for which our model predicts a rather significant enhancement (but still within experimental limits) with respect to its SM value.

4 Summary and conclusions

In this work, we have studied the impact of a minimal Dirac VLL model on the thermal history of the Universe. We have shown that, indeed, TeV-scale VLLs can induce strongly first order EW phase transitions, which would generate favorable conditions for a dynamical origin of the baryon asymmetry in the Universe. We have also discussed the collider and GW predictions of such a scenario, and assessed how it can be tested at the LHC and at future GW experiments, such as LISA.

Remarkably, such a simple setup predicts a complex phase structure of the Universe, involving three PTs: a crossover in the very early Universe ($T\gtrsim 500$ GeV), which results in EWSB, followed by two SFOPTs, both at EW–scale temperatures, the first one restoring EW symmetry, and the last one breaking it again. This non-trivial succession of PTs can be traced back to competing thermal and non–thermal effects coming from the VLLs and the SM contribution. To the best of our knowledge, our model provides a first example of single–field multistep strongly first order EW phase transitions.

Since our model exhibits two SFOPTs, it also predicts a GW signature featuring two peaks. We noticed that, for nearly all the points from our scan, the later SFOPT produces a stronger GW signal and peaks at a lower frequency than the earlier one. Generally, the signal peaks are located at frequencies in the $0.01-1$ Hz range. Meanwhile, the maximum LISA sensitivity (corresponding to the C1 scenario) is achieved for $f\sim 4$ mHz, whereas DECIGO and BBO are most sensitive to $f\sim 0.1-0.3$ Hz. Thus, even if for some points of parameter space our model feature a later PT with a GW signal strength comparable to the reach of LISA, the offset between the GW peak frequency and the frequency of maximum LISA sensitivity leads us to conclude that LISA is unlikely to detect the GW signature predicted by our model. However, the detection prospects at DECIGO and BBO are more optimistic, as their typical maximum sensitivity is achieved at frequencies close to the peak frequencies of the GW1 signal.

Finally, on the collider front, direct production of the VLLs is not a promising way of testing our model at the LHC, as the predicted cross–sections of VLL pair–production are around the $0.1$ fb level. Even without a dedicated collider analysis, we estimate that such low rates would make pair–produced VLLs difficult to study at the LHC. Furthermore, production of a VLL in association with a SM lepton has a cross section of $\mathcal{O}(10^{-4})$ fb, which is beyond doubt out of the reach of LHC. Instead, the measurement of the diphoton Higgs signal strength, $\mu_{\gamma\gamma}$, is a powerful collider probe of our scenario. The relatively light charged VLLs which, as required for inducing a SFOPT in the Early Universe, couple strongly to the Higgs, enhance $\mu_{\gamma\gamma}$ by at least $15\%$ with respect to its SM value of 1. While the current $\mu_{\gamma\gamma}$ searches still allow for such departures from the SM, the high–luminosity (HL) option of LHC will constrain $\mu_{\gamma\gamma}$ at the level of $5\%$ [61, 62, 63]. Therefore, we conclude that, through $\mu_{\gamma\gamma}$, HL–LHC will be able to fully test our scenario of VLL-induced SFOPTs.

Acknowledgments

We thank K.S. Babu, J. Braathen, D. Egana-Ugrinovic, I. Lewis, A. Morais, K. Sinha, D. Teresi, and C. Wagner for helpful discussions. This work is supported by University of Nebraska-Lincoln, National Science Foundation under grant number PHY-1820891, and the NSF Nebraska EPSCoR under grant number OIA-1557417.

Appendix A Calculation of Gravitational Wave Spectrum

Our calculation of the GW spectrum relies on the results collected in Ref. [2], while the notation follows the one from Ref [66].

We start by calculating the $\alpha$ and $\beta$ parameters, which have been already discussed in the main text. In order to determine their values, several computational steps need to be taken. First, the “bounce” action, $S_{3}(T)$, has to be evaluated (see e.g. Ref. [72]):

[TABLE]

with $\phi_{b}(r)$ being the $SO(3)$–symmetric bounce solution, which describes the field value profile of an expanding spherically symmetric bubble, with $r$ measuring the distance from the center of the bubble. For a given temperature $T$, the bounce solution $\phi_{b}(r)$ satisfies the following differential equation:

[TABLE]

where $\phi_{\rm true}$ represents the field coordinate of the true vacuum (the global minimum of the potential at a given temperature) and the prime symbol denotes differentiation with respect to $\phi$.

At finite temperature, the nucleation rate of true vacuum bubbles behaves as $\Gamma_{n}\sim A(T)\exp(-\frac{S_{3}(T)}{T})$ (see e.g. Ref. [73]). Denoting by $H$ the Hubble expansion rate, the phase transition begins when $\Gamma_{n}\simeq H$, this condition being equivalent to $\frac{S_{3}(T_{n})}{T_{n}}\simeq 142$, which serves as a definition of the nucleation temperature, $T_{n}$. The PT continues until a fraction of order unity of the universe is filled by true vacuum bubbles, which occurs at a temperature denoted by $T_{p}$ (percolation temperature). Then, as the bubbles collide, the energy stored in their walls is transferred to the plasma as heat, causing the plasma to reheat to a temperature $T_{\rm reh}>T_{p}$. Following Ref. [66], we make the simplifying approximation that $T_{n}\simeq T_{p}\simeq T_{\rm reh}\equiv T_{\rm PT}$, which is justified a posteriori by the high values of $\beta/H_{\rm PT}$ that imply a short duration of the PT. We thus define the temperature at which the PT occurs as

[TABLE]

The temperature at which the PT occurs, $T_{\rm PT}$, is an essential ingredient for computing the $\alpha$ and $\beta$ parameters, whose physical meaning was described at the beginning of this appendix. Mathematically, the two parameters are defined as [74]

[TABLE]

with $H_{\rm PT}$ being the Hubble rate corresponding to $T_{\rm PT}$, and $\rho_{\rm rad,true}(T_{\rm PT})$ the radiation energy density in the true vacuum phase, $\rho_{\rm rad,t}(T_{\rm PT})=(\pi^{2}/30)g_{\rm eff,PT}T_{\rm PT}^{4}$, where $g_{\rm eff,PT}$ denotes the number of relativistic degrees of freedom (again in the true vacuum phase) at temperature $T_{\rm PT}$. $\phi_{\rm false}$ is the $\phi$ coordinate of the false vacuum, while $\phi_{\rm true}$ has the same meaning as in Eq. (A.2). For the SFOPTs present in our model, either $\phi_{\rm false}$ or $\phi_{\rm true}$ are equal to 0.

Having defined $\alpha$ and $\beta$, we now concentrate on the calculation of the GW spectrum. As discussed in Sec. 3.1, there are three sources of GWs produced at a SFOPT: bubble collisions, sound waves in the plasma, and magnetohydrodynamic (MHD) turbulence. As in Ref. [2], we suppose that the three contributions to the stochastic GW background combine linearly, giving the total GW signal:

[TABLE]

The spectrum of the GWs produced by bubble collisions reads

[TABLE]

where $\kappa_{\rm col}$ is the fraction of latent heat converted into bubble wall kinetic energy, and $v_{w}$ the bubble wall speed. The bubble collision spectrum has a peak frequency given by

[TABLE]

The sound wave contribution is given by

[TABLE]

and peaks at a frequency

[TABLE]

In Eq. (A.9), $\kappa_{v}$ stands for the fraction of latent heat transferred into the bulk motion of the plasma. Finally, the MHD turbulence decay contributes to the GW spectrum as

[TABLE]

where $\kappa_{\rm turb}$ is the fraction of latent heat transferred to turbulent plasma motion, and

[TABLE]

The question of choosing the efficiency factors $\kappa_{\rm col,sw,turb}$ and the bubble wall velocity $v_{w}$ is model-dependent and involves certain calculations and assumptions regarding the dynamics of the bubble walls. Therefore, such a task is beyond the scope of the current work. Instead, we resort to a much simpler approach. Using the results from Ref. [75], in which the authors numerically express $\kappa_{v}$ as a function of $v_{w}$ 777In Ref. [75], $\kappa_{v}$ is denoted simply as $\kappa$ and $v_{w}$ as $\xi_{w}$. for different values of $\alpha$, we choose the bubble wall velocity $v_{w}$ such that it corresponds to the maximum value of $\kappa_{v}$ for a given value of $\alpha$. In our model, the strongest PTs typically have $\alpha\lesssim 0.1$, and we choose $v_{w}=0.6$, from which it follows that $\kappa_{v}\simeq 0.4$ [75]. Concerning the turbulence efficiency factor, it is given by $\kappa_{\rm turb}=\epsilon\kappa_{v}$, with the choice $\epsilon=0.05$ [66]. Finally, we take for definiteness $\kappa_{\rm col}=0.5$, but nevertheless mention that the choice for $\kappa_{\rm col}$ plays little role in our analysis, as the sound wave contribution is the one dominating by far the GW spectrum predicted by our scenario.

Appendix B SM–VLL Mixing for Collider Phenomenology

In absence of mixing with the SM fermions, the lightest VLL of our model would be stable. For the range of $N_{1}$ and $E_{1}$ masses predicted by our scenario, this would not be a viable option. On the one hand, if $m_{E_{1}}>m_{N_{1}}$, then a stable $N_{1}$ would not be a suitable Dark Matter (DM) candidate: the large Yukawas necessary for a SFOPT would induce a strong $SU(2)_{L}$–doublet component in $N_{1}$. This would imply a sizeable $ZN_{1}N_{1}$ coupling, which would be in conflict with null results from DM direct detection experiments [76]. On the other hand, for $m_{N_{1}}>m_{E_{1}}$, $E_{1}$ would be a stable charged particle, but we choose to not pursue this possibility. Therefore, in order to avoid stable VLLs, our model has to feature a mixing between the VLL sector and the SM fermions.

We thus choose to introduce a small $\tau$ lepton–VLL mixing in our model, which is achieved by adding the following Yukawa terms to the Lagrangian in Eq. (2):

[TABLE]

where $L^{3}_{L}$ is the third generation SM lepton doublet. For simplicity, we suppose that the SM neutrinos do not mix with the neutral VLLs. When presenting our collider predictions for the benchmarks in Sec. 3, we choose $y_{1}=y_{2}=0.05$. We have explicitly checked that these values for $y_{1,2}$ predict deviations from the SM values in the $W\tau\nu$ and $Z\tau\tau$ couplings which are below the sensitivity achieved at the LEP experiment [42] or in $\tau$ lifetime measurements [48]. More precisely, the precision in the measurement of the $Z\tau\tau$ axial coupling and the $W\tau\nu$ couplings (from $\tau$ lifetime 888Since all the (tree level) $\tau$ decays are proportional to the $W\tau\nu$ coupling squared, a rescaling of the latter would change the $\tau$ lepton lifetime accordingly. measurements) is at the permille level, while the precision for the $Z\tau\tau$ vectorial coupling measurement is at the percent level (see e.g. PDG [48]). Meanwhile, in our model, we checked that $\tau$–VLL mixing amounts to shifts with respect to the SM which are typically one order of magnitude below the corresponding experimental sensitivities.

We now briefly discuss the deviation in the $h\tau\tau$ Yukawa coupling induced by $\tau$–VLL mixing. We start by noting that, in order to address this problem, we would in principle need to add another term in the Lagrangian, namely the SM–like term $y_{\tau}\overline{L}_{L}^{3}H\tau_{R}+{\rm h.c.}$, with $y_{\tau}$ a free parameter. However, once the values of all the other Yukawas are specified, $y_{\tau}$ can be chosen such that the $\tau$ mass central experimental value, $m_{\tau}\simeq 1.777$ GeV, is reproduced. Once this step is performed, the result is a deviation in the $h\tau\tau$ physical Yukawa coupling. We have checked that for our values of $y_{1,2}$ and the scanned values of the other parameters in Eq. (2), the $h\tau\tau$ coupling never deviates by more than $\sim$5% from its SM value, which is below the current sensitivity achieved by meaurements of the $h\to\tau\tau$ decay at the LHC [77, 60].

Appendix C Theoretical Constraints

In this Appendix, we briefly overview the theoretical constraints affecting our model. Any model containing new fermions that couple strongly to the Higgs suffers from theoretical inconsistencies, and can only be viewed as a low–energy effective description of a more fundamental UV completion. Indeed, the effect of large Yukawas is two–fold: on the one hand, they tend to push the one–loop effective potential towards negative values for large values of $\phi$; on the other hand, the Renormalization Group Equations (RGEs) drive the strong Yukawas towards Landau poles, which signal the breakdown at high energies of the theory under consideration. Moreover, too high values of the Yukawa couplings can lead to unitarity violation in various processes, such as VLL–VLL scattering. Therefore, in the following, we discuss the problems of vacuum stability, Yukawa Landau poles, and unitarity in our model. For concreteness, we analyze the benchmark point “BM1”, discussed in Sec. 3. Since BM1 exhibits the strongest later PT among the points in our scan, we expect that the above theoretical consideration will be most constraining for this point, as the other points typically present weaker Yukawa couplings.

We start by considering the issue of vacuum stability. Due to the large Yukawa couplings we consider, loop effects associated with the new fermions destabilize the effective potential at high field values, rendering it unbounded from below. This has catastrophic consequences for the EW vacuum, which becomes unstable. For instance, in the case of benchmark point BM1, the zero–temperature one–loop effective potential becomes negative at $\phi\simeq 2.5$ TeV, which is slightly above the mass scale of the heavier VLL eigenstates, $m_{N_{2}}\simeq 1.6$ TeV and $m_{E_{2}}\simeq 1.8$ TeV.

The simplest solution for stabilizing the EW vacuum is to add the dimension–6 effective operator $\left(H^{\dagger}H\right)^{3}/\Lambda^{2}$ [38], which can arise, for example, from a UV completion featuring compositeness/strong dynamics. Indeed, we find that, for $\Lambda\leq 3.2$ TeV, the potential increases monotonically for $\phi>v$ and thus is no longer unbounded from below. Therefore, the EW minimum becomes absolutely stable. For the slightly higher value of $\Lambda=5$ TeV, the zero–temperature potential develops a second minimum, deeper than the EW one, at $\phi\sim 6.4$ TeV. In this case, tunneling from the EW minimum to the deeper one becomes possible. We find that the $O(4)$–symmetric bounce action for this transition is $S_{4}\simeq 20$, and, using the formalism from Refs. [78, 79], we estimate the EW vacuum lifetime to be roughly 40 orders of magnitude lower than the age of the Universe. Clearly, such a situation is ruled out. For higher values of the cutoff scale, the EW vacuum becomes even more short–lived.

Before moving on, we mention that, as argued in Ref. [38], we expect the $\left(H^{\dagger}H\right)^{3}$ operator to have a negligible contribution (of the order $v^{2}/\Lambda^{2}$) on the dynamics of the PTs in our model.

We now turn our attention to the Landau poles appearing in our model. For our RGE analysis, we have adapted the general 1–loop Yukawa RGEs from Ref. [80] to our case, taking into account only the running of the VLL and top Yukawas and neglecting the subleading contribution of the gauge couplings. The resulting beta functions are given by

[TABLE]

with $\beta_{y}\equiv\frac{dy}{d\log\mu}$. The beta functions for $y_{N_{R}}$ and $y_{E_{L,R}}$ can be obtained from the second line of Eq. (C.1) by an appropriate interchange of the indexes, i.e. $N\leftrightarrow E$ and/or $L\leftrightarrow R$. As an example, interchanging $N$ and $E$ in the expression of $\beta_{y_{N_{L}}}$ gives $\beta_{y_{E_{L}}}$. For the initial conditions of the RGEs, we consider a starting energy scale $\mu_{0}$ and take $y_{t}(\mu_{0})=1$, whereas for the VLL Yukawas at $\mu_{0}$ we take the values listed in Table 1 in the BM1 column. Under these assumptions, we find that the Yukawa Landau poles occur at a scale around $\sim 20\,\mu_{0}$. Thus, taking $\mu_{0}=500$ GeV, which corresponds to a scale associated with the lighter VLL mass eigenstates, would lead to a Landau pole at $\mu\sim 10$ TeV. This finding is consistent with a UV–completion at $\Lambda\leq 5$ TeV that would stabilize the potential, as discussed in the previous paragraph.

Finally, we study the partial wave unitarity constraints on the possible values of Yukawa couplings. We focus on the $2\to 2$ process of VLL–VLL scattering, as its amplitude scales as the square of VLL Yukawa couplings. We neglect the subleading contributions coming from gauge boson exchanges, and, using the Feynman gauge, we consider only Feynman diagrams involving internal scalars (Higgs and Goldstones). We work in the high–energy limit where all the masses involved in the process can be neglected. More specifically, we consider the full transition matrix of $\psi_{i}\bar{\psi}_{j}\to\psi_{k}\bar{\psi}_{l}$ scattering amplitudes, with $i,j,k,l$ labeling the four possible VLL states, and restrain ourselves to the $++++$ helicity states (see e.g. Ref. [81] for a discussion). The $J=0$ (vanishing total angular momentum) partial wave unitarity criterion, applied to the highest eigenvalue of the $\psi_{i}\bar{\psi}_{j}\to\psi_{k}\bar{\psi}_{l}$ transition matrix, leads to the following constraint:

[TABLE]

which is automatically satisfied by imposing $|y_{\rm VLL}|<\sqrt{4\pi}$, as we did in our scan.