Gravitational waves from a scotogenic two-loop neutrino mass model

We propose a framework to account for neutrino masses at the two-loop level. This mechanism introduces new scalars and Majorana fermions to the Standard Model. It is assumed the existence of a global U(1) × Z 2 symmetry which after partial breaking provides the stability of the dark matter candidates of the theory. The rich structure of the potential allows for the possibility of first-order phase transitions (FOPTs) in the early Universe which can lead to the generation of primordial gravitational waves. Taking into account relevant constraints from lepton flavour violation, neutrino physics as well as the trilinear Higgs couplings at next-to-leading order accuracy, we have found a wide range of possible FOPTs which are strong enough to be probed at the proposed gravitational-wave interferometer experiments such as LISA.


I. INTRODUCTION
The Standard Model (SM) is a highly successful theory that describes the electromagnetic, strong and weak interactions whose predictions have been experimentally verified at the LHC with the highest degree of accuracy.However it has several unaddressed issues such, for example the current pattern of SM fermion masses and mixing angles, the number of SM fermion families, the measured amount of dark matter relic density and baryon asymmmetry observed in the Universe, among others.Experiments with solar, atmospheric, and reactor neutrinos have brought evidence of neutrino oscillations caused by nonzero masses.Several extensions of the SM have been constructed in order to explain the tiny masses of the active neutrinos; see e.g.Ref. [1] for a review and a non-extensive list [2][3][4][5] of some comprehensive studies of one and two loop radiative neutrino mass models.The most economical way to generate the tiny masses of the active neutrinos considering the SM gauge symmetry, is by adding two right handed Majorana neutrinos that mix with the light active neutrinos thus triggering a type I seesaw mechanism [6][7][8][9][10][11][12], where either the right handed Majorana neutrinos have to be extremelly heavy with masses of the order of the Grand Unification scale or they can be around the TeV scale thus implying very tiny Dirac Yukawa couplings, in order to successfully reproduce the neutrino data.In both scenarios, the mixing between the active and sterile neutrinos is very tiny thus leading to strongly suppressed charged lepton flavor (clfv) violating signatures, several orders of magnitude below the experimental sensitvity, thus making this scenario untestable via clfv decays.This makes models with treelevel type-I seesaw realizations difficult to test via charged-lepton flavor-violating decays.Alternatively, radiative seesaw models are examples of interesting and testable extensions of the SM explaining tiny neutrino masses.In most radiative seesaw models, the tiny neutrino masses arise at a one-loop level, thus implying that in order to successfully reproduce the experimental neutrino data, one has to rely either on very small neutrino Yukawa couplings (of the order of the electron Yukawa coupling) or on an unnaturally small mass splitting between the CP-even and CP-odd components of the neutral scalar mediators.Two-loop neutrino mass models have been proposed in the literature [13][14][15][16][17][18][19][20][21] to provide a more natural explanation for the tiny active neutrino masses than those ones based on one-loop level radiative seesaw mechanisms.In this work we propose a minimally extended IDM theory where the scalar sector is enlarged by the inclusion of two electrically neutral gauge singlet scalars and the fermion sector is augmented by adding 4 right handed Majorana neutrinos.The gauge symmetry of the SM model is extended by including a spontaneously broken U (1) X global lepton number symmetry as well as a preserved Z 2 discrete symmetry.The scalar sector of our model is similar to the one of Ref. [22] however the fermion sector is more economical since in our model there are four right handed Majorana neutrinos whereas the model of [22] includes nine right handed Majorana neutrinos in the fermionic spectrum.Additionally the model of Ref. [22] has a U (1) B−L gauge symmetry which is not present in our model, thus leading to a different phenomenology.On the other hand, unlike the model of Ref. [21], our model does not rely in doubly charged scalar fields to implement the two loop level realization of active neutrino masses.

II. THEORETICAL STRUCTURE OF THE MODEL A. Particle content and charge assignments
In this work, we augment the SM gauge groups with an extended global abelian symmetry U(1)×Z 2 , where the U(1) associated with the global lepton number is spontaneously broken while the Z 2 is preserved.Besides the standard particle content of the SM, additional fields are present in the spectrum, including an inert SU(2) doublet scalar field, η, and two scalar singlets, φ and σ.In our setup, σ acquires a nonzero vacuum expectation value (VEV) while φ does not, such that the potential DM candidates can arise from the neutral components of η or from φ.Additional Majorana fermions are also included into this scheme, namely, two N R fields and two Ψ R .While the electroweak (EW) gauge group does not distinguish N R and Ψ R , the extended symmetry group treats both differently.The charge assignments of the fields are shown in Tab.I for the scalar fields and in Tab.II for matter fields.As shown in Tab.I, the SU(2) scalar doublet η and the scalar singlet φ have non-trivial Z 2 charges.Since these two scalars do no acquire VEVs and the charge assignments do not let neutrinos to couple to the SM Higgs, light active neutrinos do not acquire masses at tree level.The specific U(1) × Z 2 charge assignments of the N k,R Majorana neutrinos and σ scalar field enables the radiative generation of neutrino mass and the corresponding mixing only at two-loop level that can be noticed from a typical topology in Fig. 1.Namely, this is possible through the virtual exchange of the neutral components of the inert SU(2) scalar doublet η fields as well as of the real and imaginary parts of the inert singlet scalar φ, with the loop closing through the σ VEV and the VEV of the active Higgs doublet Φ.For simplicity of the analysis, we assume that the charged lepton sector is purely diagonal, such that the only contribution to the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) mixing matrix comes exclusively from the loop topology shown in Fig. 1.
and the supplemental global symmetry U(1) X × Z 2 for the new scalar fields (σ, η and φ) and the Higgs doublet (Φ).

B. Yukawa sector and scalar potential
Based on the charge assignments of Tabs.I and II, the most general and renormalisable Yukawa Lagrangian reads as ) Table II: Charge assignments under the SM gauge group (SU(3) C × SU(2) L × U(1) Y ) and the supplemental global symmetry U(1) X × Z 2 for the new matter fields (N R , Ψ R ) and the SM-like lepton doublet (ℓ i,L ) and singlet (ℓ i,R ).Here, i = 1, 2, 3 and k = 1, 2.
where we have defined η = iσ 2 η † and the superscript c in the Ψ field indicates charge conjugation.We adopt the Einstein summation convention where repeated indices indicate sum over them, with i, j = 1, 2, 3 and m, n, k = 1, 2. Here, y ℓ is a 3 × 3 matrix, y Ω and m ψ are 2 × 2 matrices and y N is a 2 × 3 matrix.The neutrino mass matrix can then be induced by the combination of the y N and y Ω Yukawa matrices.The complete gauge invariant scalar potential in this model is given as, where we take all parameters of the potential to be real, i.e.CP symmetry is conserved in the considered scenario.Notice that the quartic couplings λ 7 and λ 13 are crucial for generating neutrino masses through out the two-loop level mechanism.The nonzero scalar VEVs are denoted as ⟨Φ⟩ = v ∼ 246 GeV and ⟨σ⟩ = v σ .Also, the terms in the last line of the potential that breaks the U(1) X symmetry softly are incorporated in order to make the CP-odd scalar associated with σ massive.To simplify our analysis, we take the couplings κ i to be extremely small fom here onwards.Once the scalars Φ and σ acquire VEVs and the EW symmetry is broken, the model will have nine physical scalars -two charged scalars (η ± ), four neutral CP-even scalars (η R , φ R , h 1 and h 2 ) and three neutral CP-odd scalars (η I , φ I and σ I ).It is to be noted that the scalars associated with η and φ mix only among themselves and so do the scalar fields coming from Φ and σ 1 .

C. Scalar mass spectrum
After employing the minimization conditions along v and v σ , the scalar mass matrices can be written as 2) where the zero-eigenvalues of the charged and CP-odd matrices are absorbed by the longitudinal components of the W ± and Z 0 gauge bosons.To avoid constraints coming from the Higgs sector, we work closely to the alignment limit, such that vv σ λ 8 ≪ 1.The above 1 Let us note that similar scalar potentials can be found in scenarios that explain Dirac neutrino masses at the one-loop level [23].A detailed investigation of such scenarios, following the lines of the present work, will be dedicated in a subsequent paper.
matrices can be analytically diagonalised, which results in the mass eigenvalues, Charged scalars : CP − even scalars : CP − odd scalars : (2.3) Here, the lightest scalar m h 1 is taken to be the SM Higgs boson of mass ∼ 125 GeV.To simplify the above expressions, we have defined λ H ± = λ 5 +λ 6 +λ 7 ±λ 9 , Λ H ± = ±λ 10 +λ 12 +λ 13 in the CP-even sector and λ A ± = λ 5 + λ 6 − λ 7 ± λ 9 , Λ A ± = λ 10 ± λ 12 ∓ λ 13 in the CP-odd sector.To make the numerical analysis simpler, we work in the limit of small mixings between the two heavy CP-even and between the two CP-odd states (effectively we take λ 14 and λ 15 to be small).Under this assumption, the above expressions can be inverted to express 7 out of the 13 quartic couplings in terms of the physical scalar masses.Thus, the physical masses can be taken as the input parameters in the numerical analysis.The relations that are used in this work are, (2.4) Expressions for λ 1 and λ 4 are exact, whereas all others are valid in the approximation λ 14 , λ 15 → 0.

D. Implications for collider physics and Dark Matter phenomenology
Finally, to close this section we wish to provide a brief discussion of the key phenomenological implications of the multi-scalar sector in the considered model.It is expected that the presence of electrically charged scalars in the inert doublet will provide an extra contribution to the Higgs diphoton decay rate h → γγ, which will yield a deviation of the Higgs diphoton signal strength from the SM expectation.However, that deviation is expected to occur in the experimentally allowed range in a large region of the model parameter space.Besides, given that singlet scalar field acquires a VEV at the TeV scale, its mixing with the CP-even part of the active doublet is suppressed implying that the couplings of the 125 GeV SM like Higgs boson will be very close to the SM expectation realising the so-called Higgs alignment limit.Thus, such a scenario appears to be consistent with the SM in such a limit, at least, at the leading order (LO) level.
It is generally expected that the triple Higgs coupling can be affected compared to its SM value in typical multi-scalar extensions of the SM.In our model, the trilinear Higgs coupling, λ hhh , deviates from the SM prediction through a mixing with the σ field.In particular, at LO a simple analytical expression can be derived where α h is the CP-even mixing angle between h 1 (Higgs) and h 2 .It is not hard to notice that in the limit of α h → 0, the tree-level SM Higgs coupling λ SM hhh = 6vλ 1 is successfully reproduced as we effectively work very close to the alignment limit.Hence, the potential corrections to λ LO hhh are very suppressed and we expect no significant deviations from the SM expectation.On the other hand, previous works [24][25][26] have suggested that one-loop corrections to the λ hhh can actually be quite sizable in some BSM extensions, resulting, in some cases, in an increase of order O(50%) for certain combinations coupling/mass parameters.In our case, next-to-leading order (NLO) contributions will be dominant and represent the bulk of the new physics contribution to λ hhh .So, for a correct understanding of the impact on λ hhh , one must take these contributions into account.In this regard, we have computed all relevant one-loop contributions to the triple Higgs coupling using the formalism of Ref. [27].Note that the formulas derived here are valid in the limit of zero external momentum, and therefore, only the contributions coming from the heavier states are relevant.In particular, we have considered all possible contributions coming from all physical scalars of the model (η ± , h 1 , h 2 , η R , φ R , σ I , η I and φ I ), as well as the loop diagram involving the virtual exchange of the top quark.The loops containing the gauge bosons and the lighter scalar bosons were neglected, as they are expected to be subdominant.The full one-loop formulas are listed in Appendix A. In what follows, in our numerical analysis we will take into account the existing constraints on λ hhh available from the LHC [28,29].Future collider measurements [30,31] can further probe possible deviations of λ hhh from the SM value.Another interesting implication of our model is for understanding the structure of DM.The Majoron σ I , which is the pseudo-Nambu-Goldstone boson associated with the spontaneous breaking of the global lepton number symmetry, can be identified as the WIMP DM candidate.Its relic density can be generated by the standard freeze-out mechanism.Assuming for simplicity that the cosmological DM is dominated by such a scalar DM candidate, it would annihilate mainly into W W , ZZ, tt, bb and h 1 h 1 channels in the early Universe via a Higgs portal scalar interaction.These interactions will contribute to the DM relic density, which can be accommodated for appropriate values of the scalar DM mass, which in most cases are at the TeV scale.Such a scalar DM candidate would also scatter off a nuclear target in a detector via Higgs boson exchange in the t-channel, giving rise to a constraint on the Higgs portal scalar interaction coupling.In the low mass region, the invisible Higgs decay constraints should be taken into account and would imply a lower bound for the scalar DM candidate mass of about 60 GeV.Besides that, in a small window around half of the SM Higgs boson mass, the DM relic density constraints will be successfully accounted for.In Fig. 3, we have plotted the parameter space that gives the correct relic density in the m σ I − λ 8 plane, where λ 8 is the Higgs portal coupling.In calculating the relic density, we have assumed for simplicity that the dominant annihilation channels of the DM are the ones into W W , ZZ, h 1 h 1 , tt, bb, and the components of the doublet η.The corresponding Feyn-man diagrams are given in Fig. 2. The cross sections for the relevant annihilation channels are given as [32], ) ) where √ s is the centre-of-mass energy, N c = 3 stands for the color factor, m h 1 = 125.7 GeV and Γ h 1 is the total decay width of the SM Higgs boson which is equal to 4.1 MeV.From the above scattering cross sections, the current DM relic abundance in the Universe can be obtained as (c.f.Ref. [33,34]), where ⟨σv⟩ is the thermally averaged annihilation cross section, A is the total annihilation rate per unit volume at temperature T and n eq is the equilibrium value of the particle density, which are given as [33], with K 1 and K 2 being the modified Bessel functions of the second kind of order 1 and 2, respectively.We have taken T = m mσ I /20 following Ref.[33].The DM relic density thus determined should match the required value [35], Ω DM h 2 = 0.1200 ± 0.0012. (2.13) In generating Fig. 3, we have varied λ 8 from 0 to 4π and the masses of η R , η I and η ± in the range 0.2 -5 TeV.In this figure, the pink band is disfavored by perturbativity bounds.The gray region is excluded by the constraints from XENON1T [36] whereas the blue line corresponds to the sensitivity of DARWIN [37].From this figure, one can see that most of the allowed parameter space corresponds to a dark matter mass in the TeV range though there are a few allowed points in the sub-TeV region as well.The gray region is excluded by XENON1T [36] while the blue line corresponds to the sensitivity of DARWIN [37].

III. NEUTRINO MASS AND LEPTON FLAVOR VIOLATION
Once Φ and σ develop VEVs, a non-trivial contribution to the neutrino mass matrix is generated at two-loop level, whose diagram is illustrated in Fig. 1.The active light neutrino mass matrix in this case is given as, where n, k, r = 1, 2. Here, the loop integral I can be written as [22] G(x 2 , y 2 , z 2 )= Note that the lightest active neutrino is massless which implies that there is only one Majorana phase in the PMNS matrix.A convenient way to incorporate the bounds from the neutrino oscillation data is to express the Yukawa coupling matrix using the Casas-Ibarra parametrization.For the model considered in this paper, one can express y N as, where U * PMNS is the 3 × 3 PMNS neutrino mixing matrix, M diag ν = diag(m 1 , m 2 , m 3 ) is the diagonal neutrino mass matrix, R is a 2 × 3 matrix that is given as, R = 0 cosz −sinz 0 sinz cosz for normal hierarchy (NH) of light neutrino masses (m 1 = 0) and R = cosz −sinz 0 sinz cosz 0 for inverted hierarchy (IH) of light neutrino masses (m 3 = 0).
The matrix A is given by Under this parameterization, the experimentally observed neutrino mass differences and the PMNS mixing matrix can be given as input in the numerical scan, with the Yukawa coupling being found as output.Since we do not have control over the magnitude of the Yukawa couplings, we only allow solutions which respect perturbativity, i.e. |Y N |, |Y Ω | < √ 4π.For the purpose of numerical scans, we only consider the normal ordering scenario for the neutrino masses (that is, m 1 = 0).The presence of the Majorana fermions can induce LFV decays, such as µ → eγ, which are strongly constrained by experiment.In our model, such decays are mediated at one-loop level via virtual exchanges of the neutral fermions and the charged scalars.The branching fraction for the two-body decay process ℓ i → ℓ j γ, where i = e, µ, τ is given as [38][39][40][41] with s = 1, 2. Here, α em = 1/137 is the fine-structure constant, x , with V being the left-handed charged lepton mixing matrix, G F = 1.166364 × 10 −5 GeV −2 is the Fermi constant, which in our case is the identity matrix, m φ ± are the masses of the charged scalar components of the SU(2) L inert doublet φ, and m N sR (s = 1, 2) correspond to the masses of the right-handed Majorana neutrinos N sR , which due to the charge assignment, are generated at one-loop level.The loop function F is written as, The masses of the right handed neutrinos N sR are given as, where we have taken y Ω and m ψ to be diagonal and degenerate for simplicity.
The most stringent bounds for LFV come from muon decay measurements, namely, µ → eγ.
Current experimental results put an upper bound on the branching ratio, which reads as BR (µ → eγ) < 4.2 × 10 −13 [42].The model can also give rise to µ → e conversion in the atomic nucleus.But the existing bounds are not as restrictive as the ones from BR (µ → eγ).However, among various planned experiments, the µ → e conversion in the Al nuclei has the best projected sensitivity of ∼ 10 −17 [43].The branching ratio for this process can approximately be expressed as [41], CR(µ Al → e Al) ≈ 1 350 BR (µ → eγ) .For all the points shown in Fig. 4, neutrino oscillation data is respected, with the relevant experimental parameters (namely, neutrino mass differences and the PMNS mixing angles) being given as input and allowed to vary within their respective 3σ uncertainties [44][45][46].The Majorana phase is also varied in the range [0, π] and y Ω is taken to be diagonal where the entries are allowed to vary in the range [0, 1].Here, z is defined as x 1 + y 1 i with x 1 and y 1 being varied in the range [−2, +2], the fermions Ψ kR are taken to be degenerate, with their masses varied in the range [0.5,5] TeV.For the scalar potential parameters, we employ the inverted equations in (2.4) to allow for the physical masses to be given as input.In this regard, we have considered the following parameter ranges: • λ 2,5,8,9,10,11,12 → [0, 1].The remaining λ's are treated as outputs in the inversion procedure.Notice that we have inverted the expressions for λ 5 and λ 9 , which can be solved such that µ 2 η and µ 2 φ are instead given as output parameters, where we have taken them to be both greater than 0. Also, in the parameter space that has been considered, we have ensured that all the quartic couplings, λ i > 0, even though this is only a sufficient condition but not necessary.

IV. PHASE TRANSITIONS AND PRIMORDIAL GRAVITATIONAL WAVES
One of the more interesting phenomenological consequences of extended scalar sectors compared to the SM framework is the possibility for the existence of FOPTs.Indeed, such a phenomenon is not present in the SM (both the EW and QCD phase transitions are of second order), making it pure beyond-SM driven (BSM) physics.Such transitions can lead to the generation of PGWs.With the detection of GWs at LIGO [47], a new era of multi-messenger astronomy has arisen, where such events may give us unique perspectives into new physics beyond the SM.It has been recently proposed that the cosmological phase transitions between different vacua at finite temperatures (such as those associated with symmetry breaking in BSM models) may give rise to gravitational imprints that should be observable in future ground-and space-based experiments [48,49].The PGWs signals can also provide hints into the nature of the neutrino masses in the context of distinct seesaw mechanisms [50].As was mentioned above, PGWs can arise from FOPTs, which produce a stochastic GW background in the early Universe which is formed out of equilibrium due to fast-expanding vacuum bubbles [51].These bubbles eventually collide and merge to give rise to the GW echoes [52].The characteristics of the phase transition such as its strength and inverse duration are determined by the structure of the effective potential at finite temperatures and are highly sensitive, in particular, to the properties of the potential barrier between the phases and to the critical temperature of the transition (see e.g.Refs.[48,[53][54][55][56][57] and references therein).The phenomenon of strong FOPTs is quite generic for multi-scalar extensions of the SM (in particular, for those that undergo several symmetry breaking steps) such as the one we address in this article.Therefore, PGWs emerging from cosmological FOPTs are often considered as a potentially useful source of phenomenological information about multi-scalar extensions of the SM which may be complementary to constraints coming from collider searches.In order to fully classify and study the PGW spectra emergent in our model, there are a total of five main ingredients.The first one is inherently dependent on the considered model structure -the finite temperature effective potential [58][59][60].Generically, it can be written as [61,62] where the first term, V (0) , is the tree-level scalar potential given in Eq. (2.2).The second term is the Coleman-Weinberg one-loop correction, where F is 0 for bosons and 1 for fermions, m a (ϕ b ) is the ϕ a -field dependent mass of the particle a, n a is the number of degrees of freedom for each particle a and µ is a renormalisation scale.The particle degrees of freedom can be computed as (−1) 2s QN (2s + 1), where s is the spin of the particle, Q = 1(2) for neutral (charged) particles, and N is the number of colours.In this work, we consider the MS renormalisation scheme such that we have C a = 3/2 for scalars, fermions and longitudinally polarized gauge bosons, and C a = 1/2 for transverse gauge bosons.Additionally, the renormalisation scale is fixed as where i runs over all BSM scalars in the model.Since we allow to vary the masses of the scalar fields, µ will be different for each sampled point.The third term encodes thermal corrections to the potential, where T is the temperature, J 0 (for bosons) and J 1 (for fermions) are the standard thermal integrals found e.g. in Ref. [61].At the leading (m/T ) 2 order in thermal expansions of J a , Eq. ( 4.3) can be approximated as where in the last sum, all the fermions of the model are included.Here, M is the fielddependent Hessian matrix in the classical field approximation where all fields are fixed to their classical (background) components, that is and the degrees of freedom, n i , for different states are as follows [61]: for quarks n a = 12, for charged leptons n a = 4, for neutral leptons n a = 2, for the gauge bosons n W = 6, n Z = 3, n γ = 2 and for the longitudinal component of the photon n γ L = 1.Additional higher-order corrections must also be included in the form of Daisy (ring) diagrams [58,[63][64][65].These can be parameterized as thermal corrections to the quadratic mass terms of the bare potential and read as For the considered model, we have where g (g Y ) are the SU(2) (U(1)) gauge couplings and y i are the Yukawa couplings for each of the SM-like third-generation fermions.The gauge boson masses also gain thermal corrections and, since there are no new gauge bosons in the model, the formula is identical to what is found in the previous literature (see e.g.Ref. [61]).
The last relevant part of the effective potential in Eq. (4.1) is the counter-term potential V CT .The counter-terms are defined such that the minima of the tree-level potential remains fixed at zero temperature (i.e. the tree-level masses remain the same at one-loop level), which is guaranteed by imposing that the first and second derivatives of the counter-term potential matches the first and second derivatives of the Coleman-Weinberg one-loop potential (for a pedagogical discussion, see for instance Ref. [55] and references therein).In our case, this results in the following conditions while for other parameters, we set δλ i = δµ i = 0.The remaining four ingredients determine the FOPT dynamics, which include [52] • The inverse time-scale of the phase transition, β.This parameter essentially determines the bubble nucleation rate.Normalising it conventionally to the Hubble expansion rate, H, one writes where S(V eff ) is the Euclidean action which depends on the structure of the effective potential in Eq. (4.1).
• Nucleation temperature, T n .It represents the temperature of the Universe directly after the phase transition, and can be computed by requiring that the probability for a single bubble nucleation per horizon volume is equal to one [66].
• Phase transition strength α, which is related to the released latent heat during the phase transition.It is found in terms of the difference between the effective potential values before and after the transition via the following relation [67,68] being also normalised to the radiation energy density of the Universe at the nucleation temperature.Here, g * is the effective number of relativistic degrees of freedom in the cosmic plasma.The effective potential is evaluated at initial (metastable) phase, V eff,i , and at the final (stable) phase, V eff,f .
• Bubble wall velocity, v b .This quantity represents the speed of the new phase at the interface of the nucleating bubble.In the current we will simply assume the case of supersonic detonations with v b = 0.95.An improved analysis using the recently developed methods for the determination of the speed of sound [69] and wall velocity [70] is left for future work.
From this set of FOPT characteristics, the power spectrum associated with the production of PGWs can arise from three distinct sources, Here, h 2 Ω BC is the power for GW spectrum produced via collision of vacuum bubbles and is a particularly important effect when such bubbles "run-away", that is, the friction created by the thermal plasma at the boundaries is not enough to stop the acceleration of the wall [52].A study has shown that this effect is subleading compared to the rest, and can usually be neglected [71].The second term, h 2 Ω SW is the power of GWs emerging from sound-waves triggered in the plasma once the bubble wall sweeps through the plasma at relativistic speeds.In this case, a vast portion of the released energy is converted to the waves' propagation in the plasma that surrounds the walls as sound waves [52].The last term is due to magnetohydrodynamic turbulence effects that arise due to non-linear effects on the sound waves in the plasma that are largely uncertain but expected to be less relevant for the considered FOPTs than the sound waves.For analytical expressions that connect each term of Eq. (4.11) with the various parameters of the FOPTs introduced above, see [52,72,73] and references therein.

A. Numerical results
We have implemented the thermal effective potential of the considered scotogenic two-loop neutrino mass model into the CosmoTransitions package [74] which enable us to compute numerically the bounce action and hence the probability for bubble nucleation at any given temperature.The key parameters that control the dynamics of the phase transitions, β/H and α, depend on derivative of the Euclidean effective action, hence, they can exhibit numerical instabilities if the action is not well-behaved.To minimize the corresponding numerical uncertainties, we perform an interpolation of the action on a point-by-point basis.The technique for the smoothing of the Euclidean action is described in detail in some of the author's previous work [55].Furthermore, for the scan, we have extended the ranges of various parameters to be more inclusive.Namely, we now perform a logarithmic scan over the quartic couplings in the ranges of [10 −8 , 4π] and the masses of the BSM scalars to be above 200 GeV.Constraints from neutrino physics and LFVs are also imposed in the data points.
First, we show in Fig. 5 the spectral GW peak signal, h 2 Ω peak GW , as a function of the peak frequency, f peak in Hz; both axes are shown in logarithmic scale.Focusing first on the panels (a), (b) and (c) where we display in the colour axis the relevant GW observables for the dynamics of the signal, namely, the phase transition strength (a), the inverse time-scale of the phase transition (b), and the nucleation temperature (c).Let us note that there are several planned GW interferometers that could explore the PGW signals predicted by the model.They include the Laser Interferometer Space Antenna (LISA) [73,75,76], the Big  5 TeV whereas the h 2 mass is below 3 TeV.As for the last panel (f), we do not find any correlation with µ → eγ process.Indeed, we find that various distinct possibilities for its branching ratio can be accommodated across different values of the frequency and the signal strength.
In general, finding certain complementarity between the physical parameters of the model and PGW observables is relevant for exploring the model's parameter space.In the context of PGWs studies, it has been noted in Refs.[80][81][82][83] that potentially large corrections to the trilinear Higgs coupling compared to the SM expectation can significantly affect the type and the strength of the cosmological phase transitions.Indeed, strong FOPTs are typically driven by large trilinear scalar couplings as they strongly influence the potential barrier between different vacua.We show in Fig. 7 the scatter plots for the BSM-to-SM ratio of the trilinear Higgs coupling computed for our BSM scenario versus its SM value, i.e. κ ≡ λ BSM hhh /λ SM hhh , as a function of the sine of the CP-even mixing angle, with the logarithm of the GW signal strength in colour.Here, λ BSM hhh = λ LO hhh +λ NLO hhh is the full NLO result for the trilinear Higgs coupling in our model and λ SM hhh corresponds to the SM prediction.Most of the generated points satisfy sin α h < 0.07.Indeed, for such low values of the scalar mixing angle, Eq. (2.5) asymptotically tends towards the SM prediction.However, NLO contributions result in a sizeable contribution to λ hhh .We can also observe from panel (b) that the points with larger GW signal strength tend to accumulate as one approaches the Figure 7: The BSM-to-SM ratio of trilinear Higgs couplings κ ≡ λ BSM hhh /λ SM hhh as a function of the sine of the CP-even mixing angle α h .On the colour axis we show the logarithm of the GW intensity.On the panel (a), the results from the full collected dataset are showcased, whereas on (b) a cut on h 2 Ω peak peak was imposed, such that only points which are visible at future GW experiments (LISA, BBO and DECIGO) are displayed.We note that all constraints from neutrino physics and LFV are applied here.Regions marked by gray bands are excluded at 95% confidence level (CL) based on the latest results from the LHC, the blue band region corresponds to the expected sensitivity at the high-luminosity phase of the LHC [30] while the red band region corresponds to the expected sensitivity at the high-energy phase of the LHC [31].
alignment limit.However, this does not represent a solid trend as we find signal strengths well below the sensitivity ranges of LISA, BBO and DECIGO also for low values of the CP-even mixing angle, as can be readily noted in panel (a).Besides, one observes that there is a part of the parameter space that can already be excluded by the current data (gray bands) and it can be further reduced by future measurements both at the high-luminosity and high-energy phases of the LHC.However, the message to retain is that even within those regions there still exist parameter space points that feature FOPTs potentially detectable at future GW facilities.In Fig. 8 we show a selection of plots showing the impact of different quartic couplings and physical masses in the energy density amplitude of the GW spectra (colour axis).For panels (a), (b) and (c), we find various islands with distinct ranges of quartic couplings which lead to a strong PGW signal having both O(1) couplings, such as e.g.λ 2 and λ 10 and much smaller couplings such as λ 8 .Besides, we note that λ 8 is correlated with the size of the mixing for the CP-even states, and its smallness is related to the fact that we work very close to the alignment limit.Additionally, taking into account the results shown in Fig. 3, the constraints on λ 8 imposed by DM direct detection (both DARWIN and XENON1T) do not impact the presence of phase transitions and the corresponding GW spectra, since much lower values of λ 8 allow for the generation of phase transitions.Much tighter limits on the Higgs portal coupling would be needed if one would need to confront both DM and GW physics.In panels (d), (e) and (f) we notice that a wide range of physical scalar masses can be successfully accommodated with the observable GW signals.Let us stress that we did not find neither correlations nor further restrictions on the quartic couplings λ 7 and λ 13 involved in the neutrino mass generation.These parameters can take either small or large values without modifying the conclusions of this section.It can also be noted that no states is small, we have performed a parameter scan of the model taking into account relevant constraints from LFV processes (such as µ → eγ decay) and those from neutrino physics.
We have explored the model's parameter space in detail and identified the domains that can successfully accommodate the existing constraints and giving rise to strong FOPTs that lead to observable PGWs at LISA.In addition, we have taken into account the constraints on the trilinear Higgs coupling at the NLO, and found a number of points with strong PGW signals that can be probed in both the collider and GW measurements in the future.We have found that the current constraints from the Higgs portal coupling does not impose severe constraints to the presence of observable GWs coming from phase transitions.where m t ≈ 172.76 GeV is the top mass.The NLO contribution coming from the scalar fields reads where for simplicity of presentation we have k = 1 . . .8 ≡ η ± , h 1 , h 2 , η R , φ R , σ I , η I , φ I .The superscript "(0)" indicates tree-level accuracy.The loop function F is defined as [27] F(m a 1 . . .m a N ) =

Figure 1 :
Figure 1: Two-loop Feynman diagram which is responsible for the generation of the neutrino masses and mixing.Here, ν are the neutrinos, η R,I are the real and imaginary parts of the neutral component of the η doublet, φ R,I are the real and imaginary parts of the φ singlet, v and v σ are the VEVs of the Φ doublet the σ singlet fields, respectively.

Figure 2 :
Figure 2: Dominant Feynman diagrams contributing to the DM annihilation for the scenario considered.

8 Figure 3 :
Figure3: Parameter space that gives the correct relic density in the m σ I − λ 8 plane.The pink band is disfavored by perturbativity.The gray region is excluded by XENON1T[36] while the blue line corresponds to the sensitivity of DARWIN[37].

(3. 7 )Figure 4 :
Figure 4: The predictions for BR (µ → eγ) (green points) and CR(µ Al → e Al) (red points) as functions of Tr[y † N y N ].All the points shown are within the existing experimental bounds marked by the horizontal dashed line.The cyan points correspond to the region within the projected sensitivity for µ Al → e Al conversion.In this regards, in Fig. 4 shows the predictions for BR (µ → eγ) (green points) and CR(µ Al → e Al) (red points) as functions of Tr[y † N y N ].Only the values that are within the existing bounds are shown, whose limit is indicated by a dashed horizontal line in the plot.The cyan points correspond to the region within the projected sensitivity for µ Al → e Al conversion.For all the points shown in Fig.4, neutrino oscillation data is respected, with the relevant experimental parameters (namely, neutrino mass differences and the PMNS mixing angles) being given as input and allowed to vary within their respective 3σ uncertainties[44][45][46]. The Majorana phase is also varied in the range [0, π] and y Ω is taken to be diagonal where the entries are allowed to vary in the range [0, 1].Here, z is defined as x 1 + y 1 i with x 1 and y 1 being varied in the range [−2, +2], the fermions Ψ kR are taken to be degenerate, with their masses varied in the range [0.5,5] TeV.For the scalar potential parameters, we employ the inverted equations in (2.4) to allow for the physical masses to be given as input.In this regard, we have considered the following parameter ranges:

Figure 6 :
Figure 6: Signal-to-noise ratio (SNR) plot for a LISA mission profile of 3 years.In the x-axis we show the strength of the phase transition (α) and in the y-axis the inverse time duration in Hubble units.In the colour axis, we indicate the magnitude of the phase transition h 2 Ω peak GW .The coloured isolines (corresponding to the values 1, 5, 20, 50, 100 and 200) are representative of the expected values for the SNR.As a general rule of thumb, points with an SNR greater than 5 can be potentially observed.The grey area corresponds to the region where the sound-wave component is dominant.All points shown in the plot correspond to the points potentially detectable by the LISA experiment (that is, above the red LISA curves of Fig. 5).

Table I :
(3arge assignments under the SM gauge group (SU(3