Three-Loop Inverse Scotogenic Seesaw Models

We propose a class of models providing an explanation of the origin of light neutrino masses, the baryon asymmetry of the Universe via leptogenesis and offering viable dark matter candidates. In these models the Majorana masses of the active neutrino are generated by the inverse seesaw mechanism with the lepton number violating right-handed Majorana neutrino masses µ arising at three loops. The latter is ensured by the preserved discrete symmetries, which also guarantee the stability of the dark matter candidate. We focus on one of these models and perform a detailed analysis of the phenomenology of its leptonic sector. The model can successfully accommodate baryogenesis through leptogenesis in both weak and strong washout regimes. The lightest heavy fermion turns out to be a viable dark matter


Introduction
Several New Physics models beyond the Standard Model (SM) of particle physics have been proposed to accommodate the neutrino oscillation phenomenon and thus to generate masses and mixings for the active neutrinos.The simplest and most elegant extensions of the SM rely on adding new particles while keeping the same gauge symmetry to generate neutrino masses and mixings at tree level: that is, the seesaw mechanism.Among the ingredients, there are righthanded (RH) Majorana neutrinos for the type-I seesaw, scalar isospin doublets for the type-II, and a mixture between the latter two for the type-III seesaw, which necessitates the inclusion of fermion triplets.In most cases, to comply with neutrino data, the additional states are either too heavy to be detected or, if lighter, they couple to the SM via tiny Yukawa couplings.In both cases, the possibilities of testing tree-level neutrino mass generation are very limited, unless the model is enlarged with some extra symmetries (like in the case of minimal flavor violation) and further fields.In addition, it is difficult to comply with the total relic dark matter (DM) density of the Universe with these tree-level scenarios.
Alternatively, radiative seesaw models are viable and testable extensions of the SM explaining the tiny neutrino masses and their mixings, while the seesaw mediators play an important role in successfully accommodating the observed amount of DM relic density.In most radiative seesaw models, light neutrino masses are generated at the one-loop level.Additionally, to comply with neutrino data, one needs to assume very small neutrino Yukawa couplings, of the order of O(10 −6 ), or unnaturally small mass splitting between the CP-even and CP-odd components of the neutral scalar mediators; see e.g.Ref. [1] for a review and Refs.[2,3] for comprehensive studies of one-and two-loop radiative neutrino mass models.Two-loop neutrino mass models have been proposed in the literature [4][5][6][7][8][9][10][11][12] to provide a more natural explanation for the tiny active neutrino masses than those based on the one-loop radiative seesaw.
In this work, we consider neutrino masses generated at the three-loop level  to provide a more natural explanation for the smallness of active neutrino masses than those relying on oneor two-loop seesaw realizations.We recall that in the latter constructions, a significant number of new particles has to be considered, small Yukawa couplings are required, and most of the time one type of DM candidate is available, usually fermionic.In our previous work [31,33], we proposed an extended inert doublet model, in which light-active neutrino masses arise from a three-looplevel seesaw mechanism, realized by enlarging the SM group with a spontaneously broken global symmetry U (1) ′ and a preserved Z 2 parity that forbids the generation of neutrino masses at one-and two-loop orders.The scalar sector was also enlarged by including four neutral gauge singlet scalars, whereas the fermion content was augmented with two RH Majorana neutrinos.
We showed that this model is consistent with neutrino oscillation data and can successfully accommodate the DM relic abundance while being consistent with bounds arising from charged lepton flavor violation (cLFV) and electroweak precision observables (oblique parameters S, T , and U , in addition to being consistent with the recently observed W -mass anomaly [34]).
In the present work, we propose to further explain the origin of the RH Majorana neutrino masses at three-loop order, which, together with a dynamical origin of the lepton number violation (LNV) at the keV-MeV scale, leads to an inverse seesaw (ISS) realization.For this purpose, we investigate which class of three-loop seesaw realization can generate the ISS neutrino texture and consider two topologies of three-loop diagrams, based on which we propose three well-motivated models where Majorana mass submatrices are generated at the three-loop level.
These three models are characteristic examples of a class of models with different symmetries and particle field content that are capable of dynamically generating an ISS mechanism for light and active neutrino masses, with a LNV (∆L = 2) source parameter µ generated at three-loop order at the keV-MeV scale.This was inspired by the general formulation of the ISS [35], where the smallness of µ was attributed to the supersymmetry breaking effects in a E 6 scenario inspired by superstrings.In the context of a non-supersymmetric SO(10) model, which contains remnants of a larger E 6 symmetry, the other relevant terms of the neutrino mass matrix are generated at two loops, while µ is generated at higher loops, justifying its smallness [36].
Here, we focus on the potential of one of these three models and conduct a thorough study on the scalar sector, the dynamical generation of the ISS, the viable DM candidates, and their possible phenomenological impact.Regarding the other two example models, we briefly summarize them, leaving their more detailed investigation for future work.
The paper is organized as follows: After a detailed description of Model 1 in Section 2, in Section 3 we study its scalar sector and the neutrino mass generation.Section 4 is devoted to the phenomenological consequences of Model 1, particularly in the violation of the charged-lepton flavor and in the prospects for neutrinoless double-beta decay.Solutions to baryon asymmetry of the Universe (BAU) through leptogenesis and DM problems are discussed in Sections 5 and 6, respectively.The interplay between LNV, charged-lepton flavor violation, and DM is discussed in Section 7. Models 2 and 3 are briefly discussed in Section 8. Finally, the main findings of this work are collected in Section 9.

Model setup
In this section, we construct a model with scotogenic ISS at three-loop level, which is an extension of the SM with gauge singlet fields: three scalars φ 1,2 , σ, two left-handed Majorana neutrinos Ω 1,2 and two vector-like neutral leptons Ψ 1,2 .The SM gauge symmetry is extended with the global symmetry U (1) ′ ⊗ Z 2 .As will be discussed in the next section, the model scalar potential develops tree-level instability forming vacuum expectation values (VEVs) of ⟨ϕ⟩ = v ϕ , ⟨σ⟩ = v σ that break the symmetry according to Field where ϕ corresponds to the SM Higgs doublet.The global U (1) ′ symmetry is spontaneously broken at the TeV scale by the VEV of ⟨σ⟩ down to a residual preserved Z 4 symmetry.This results in the new particles acquiring masses on the TeV scale and hence inducing interesting phenomenological consequences, as discussed in the next sections.Other new scalars with nontrivial Z 2 charges do not acquire VEVs to maintain this symmetry unbroken.The assignment of charge for the leptonic fields for the extended symmetry in Eq. ( 1) is shown in Table 1.We do not show the quark field assignments since in the present work we focus on the lepton sector.
Furthermore, it is worth mentioning that spontaneous breaking of the global U (1) ′ symmetry will give rise to a Goldstone boson corresponding to the imaginary part of the gauge singlet scalar field σ.Assuming that the global U (1) ′ symmetry is also softly broken, such a Goldstone boson will acquire a mass via a U (1) ′ soft-breaking mass term µ 2 sb σ 2 + H.c. in the scalar potential.Note also that this term does not generate other soft breaking terms radiatively and, therefore, its introduction does not change other aspects of the model except for engendering a mass to the U (1) ′ Goldstone boson.Alternatively, the Goldstone boson can also acquire a mass from non-perturbative QCD effects associated with the mixed (SU (3) C ) 2 U (1) ′ anomaly, in an extended version of our model where the quark fields will have non-trivial U (1) ′ assignments.
As discussed in detail in Ref. [37], the Goldstone boson arising from the spontaneous breaking of the global U (1) ′ symmetry can be a scalar DM candidate for appropriate values of its mass and couplings and can play a crucial role for Electroweak Phase Transition.Here, we do not consider a detailed analysis of that scenario, which deserves another study, which is left beyond the scope of our work.The Yukawa interactions relevant to the neutrino mass generation in our model are given by Scalar fields φ 1 and φ 2 do not acquire vacuum expectation values, but together with the heavy neutral leptons Ψ kR , Ψ kL and Ω kL (with k = 1, 2) induce the LNV Majorana mass term µ nk N nR N C kR at the three-loop level according to the diagram shown in Fig. 1.Note that we assign lepton numbers L to the fields so that the accidental U (1) L lepton number symmetry is softly broken by the mass term m Ω Ω L Ω C L .Therefore, the tiny active neutrino masses are protected in our model by this accidental U (1) L symmetry and are technically natural.
It is worth mentioning that the case of only one generation of neutral leptons ν R , N R , Ψ L,R , Ω L would imply that the entries of the light active neutrino mass matrix arising from the inverse seesaw will be of the form M ν ij = z i z j (i, j = 1, 2, 3).Such a low-energy neutrino mass matrix will only have one nonvanishing eigenvalue, thus implying that only one active neutrino will acquire a mass, which is in contradiction with the experimental data on neutrino oscillations.
As requested by our strategy, the U (1) ′ ⊗ Z 2 symmetry, as well as the Z 4 symmetry arising from the spontaneous breaking of the U (1) ′ , forbid tree, one-loop and two-loop-level mass generation for active neutrinos, while allowing the three-loop contribution in Fig. 1.
In the next sections, the phenomenology of this model will be studied in detail, in particular the scalar potential, the generation of neutrino masses, the DM phenomenology and the dynamical generation of the BAU through leptogenesis.

The scalar potential and generation of neutrino masses
Hereinafter, the phenomenology of the described model will be studied in detail.The most general scalar potential invariant under the symmetry in Eq. ( 1) reads where the coefficients µ i have mass dimension, while the quartic couplings λ i are dimensionless.
The scalar fields of the model can be expanded as with k = 1, 2. Here, ϕ ± and ϕ 0 I are the SM would-be-Goldstone bosons, while σ I is a massless physical CP-odd Goldstone boson, sterile with respect to the SM interactions.Due to the preserved Z 2 symmetry, the neutral CP-even components ϕ 0 R and σ R do not mix with the remaining neutral CP-even scalar fields φ kR .The squared mass matrix for the neutral CP-even scalar fields that transform trivially under the preserved Z 2 ⊗Z 4 symmetry, on the basis ϕ 0 R , σ R can be written as where in the decoupling limit λ 5 → 0, ϕ 0 R corresponds to the SM-like 125 GeV Higgs boson.The squared masses for the inert CP-even scalar fields φ 1R , φ 2R and for the CP-odd scalars φ 1I , φ 2I are given by Finally, the conditions for minimization of the scalar potential yield the following relations

Stability of the scalar potential
In what follows, we derive tree-level stability conditions that the scalar potential must satisfy in order to be bounded from below.To this end, it is sufficient to analyze the quartic interactions because they dominate the behavior of the scalar potential in the region of very large values of the field components.We introduce the following Hermitian bilinear combinations of the scalar fields which allow to express the quartic couplings in the form or, equivalently, to Following the procedure used for analyzing the stability described in Refs.[38][39][40][41], we find that the scalar potential is stable when the following conditions are fulfilled: These conditions are imposed in the numerical analysis presented in the following sections.

Generation of neutrino masses and lepton mixing
The neutrino Yukawa terms in Eq. ( 4) give rise to the following neutrino mass terms where (m While the submatrices m νD and M are generated at tree level, the entries of the 2 × 2 submatrix µ, which violate the total lepton number, are generated at three-loop level, as can be seen in Fig. 1.
The entries of the 3 × 2 Dirac submatrix m νD are given by with i = 1, 2, 3 and k = 1, 2. The entries of the 2 × 2 matrix µ induced at three-loop are with n, k, r, s, p = 1, 2, in term of the 3-loop function [13,42,43] As shown in detail in Ref. [44], the full rotation matrix that diagonalizes a neutrino mass matrix of the form of Eq. ( 19) is given by: where It is worth mentioning here that the whole 7 × 7 matrix U is unitary.In the above matrix, R ν , R ν , respectively.Regarding the latter, the obtained mass spectrum exhibits the typical spectrum pattern of the ISS mechanism: light masses corresponding to the active neutrino states and quasi-degenerate masses corresponding to the heavy, mostly sterile states forming pseudo-Dirac neutrino pairs, in which the mass gap is proportional to entries of the lepton number violating mass matrix µ.This can be seen in the following mass matrix eigenstates obtained from the diagonalization [45][46][47][48] where M ν corresponds to the (mostly) active neutrino mass matrix, while M (−) ν and M (+) ν are the exotic neutrino mass matrices.Note that the physical neutrino spectrum is composed of three light active neutrinos and four exotic neutrinos, forming two pairs of pseudo-Dirac neutrino states, with masses ∼ ± 1 2 M + M T and a small splitting µ in each pair.Also note that the rotation matrix R ν diagonalizing the light neutrino mass one, M ν , corresponds to the U PMNS lepton mixing matrix, which would deviate from unitarity, given the presence of heavy exotic neutrinos mixing with the active ones with the mixings in R 1 and R 2 .
Furthermore, using Eq. ( 23) we find that the neutrino fields T are related with the neutrino mass eigenstates by the where the total mixing matrix UN L is given by Eq. ( 23), ν jL (j = 1, 2, 3), N where being m 1 , m 2 and m 3 the masses of the light active neutrinos, R ν the U PMNS leptonic mixing matrix (without unitarity violations) and O a complex orthogonal rotation matrix.

Neutrinoless double beta decay
In our model, the only source of LNV is the Majorana neutrino mass terms.Consequently, neutrinoless double beta decay 0νββ is mediated by the well-known neutrino mass mechanism.
In this case, its amplitude is proportional to the effective mass parameter [53,54] written for N = 7 Majorana neutrino mass states existing in the model.Here, p 2 is usually interpreted as the mean square of the Fermi momentum of the nucleon for the decaying nucleus.
Its value depends on the nucleus and the nuclear structure model used to calculate it.Here, we take the average value for the experimentally interesting isotopes p 2 ≃ −(150 MeV) 2 [54,55].
The physical neutrino states in our model are light (active) neutrinos and pairs of heavy pseudo-Dirac neutrinos with opposite CP parity.Due to the last fact, the contribution of the pseudo-Dirac pairs to the effective mass in Eq. ( 30) almost completely annihilates.The remaining contribution comes from light neutrinos.Thus, we have which we use in our analysis.
As mentioned in the previous section, the sterile neutrino spectrum of Model 1 is composed of four TeV-scale neutrinos that are practically degenerate.These heavy sterile neutrinos mix with the active ones, with mixing angles of the order of 1 The admixture of the heavy sterile neutrinos in the left-handed charged current SU (2) L ⊗ U (1) Y weak interaction gives rise to the decay l i → l j γ at one-loop level, with a branching ratio given by [32,[56][57][58][59][60][61][62] Br where Γ µ ≃ 3×10 −19 GeV is the total muon decay width and R is given in Eq. (25).Furthermore, to successfully reproduce the experimental neutrino oscillation data, the Dirac submatrix m νD , in the basis of the diagonal SM charged lepton mass matrix, should have the form given in Eq. (28).

CR(µ − e, N)
The µ − e conversion occurs in a muonic atom formed when a muon is captured, falling into the first state of a target nucleus N .The conversion rate is defined as Box and penguin diagrams contribute as where Γ capt (Z) denotes the capture rate of a nucleus with atomic number Z [63], G F is the Fermi constant, m µ the muon mass, α ≡ e 2 /(4π), with s w corresponding to the sine of the weak mixing angle.The form factors F µe q (q = u, d) are given by where Q q denotes the quark electric charge The quantities F µe γ , F µe Z and F µeqq box correspond to the different form factors of the diagrams, and G µe γ corresponds to the dipole term; all expressions are collected in Appendix A. The relevant nuclear information (nuclear form factors and averages over the atomic electric field) is encoded in the form factors D, V (p) , and V (n) .In our analysis, we use the numerical values presented in Ref. [63].
The expression for BR(µ → eee) [64,65] we use is which contains the same form factors as those entering in CR(µ − e, N), although in different combinations, see Appendix A.

Leptogenesis
In this section, we will analyze the implications of our model for the baryon asymmetry of the Universe and the possibility of generating it via leptogenenesis.In the following, we will only consider its viability in the presence of out-of-equilibrium CP and lepton number violating processes.We will not solve the Boltzmann equations, which deserve a dedicated study beyond the scope of this work.Instead, we will discuss whether or not the regimes explaining mixings and light neutrino masses also fulfill the first necessary conditions for a viable leptogenesis.
To simplify our analysis, we assume that M is a diagonal matrix and we consider the case where We further assume that the gauge singlet neutral lepton Ω L as well as the dark scalar singlet φ 1 are heavier than the lightest pseudo-Dirac fermions N ± 1 = N ± , while for simplicity we work on the basis of a diagonal SM charged lepton mass matrix.In the scenario mentioned above, only the first generation of N ± k (k = 1, 2) can contribute to the BAU.Then, the lepton asymmetry parameter, which is induced by the decay process of N ± , is [66,67] with Neglecting the interference terms involving the two different sterile neutrinos N ± , the washout parameter K N + +K N − is huge, as mentioned in Ref. [47].However, small mass splitting between pseudo-Dirac neutrinos leads to destructive interference in the scattering process [68].The washout parameter including the interference term has the form where and H corresponds to the Hubble expansion rate of the universe.In the case of a standard cosmological scenario where the total energy density is dominated by SM radiation, where g ⋆ corresponds to the number of relativistic degrees of freedom in the SM bath, and GeV is the reduced Planck mass.
In the weak and strong washout regimes, the BAU is related to the lepton asymmetry [67] as follows In the analysis, we constrain the model in order to successfully reproduce the experimental value for the baryon asymmetry parameter [69] Y ∆B = (0.87 ± 0.01) × 10 −10 .
6 Dark matter relic density The residual symmetry protects the lightest odd state of the dark sector, rendering it a viable candidate for DM.Depending on the mass hierarchy, we could have scalar or fermionic DM candidates.To simplify our analysis, we assume without loss of generality that φ 2DM is sufficiently heavy to allow its decay, then implying that the lightest among the φ 1R , φ 1I , that is, φ 1DM is the scalar DM candidate.In addition to that, in this scenario, the DM candidate scalar φ 1DM will annihilate into a couple of SM particles and pairs of BSM scalars σ R σ R , σ I σ I through the Higgs portal scalar interactions λ 6 (ϕ † ϕ)(φ * 1 φ 1 ) and λ 8 (σ * σ)(φ * 1 φ 1 ), respectively.Furthermore, the scalar DM candidate φ 1DM can also annihilate into neutrino pairs ν i ν j (active-active), (where k, r = 1, 2 and Ψ (±) k are the physical states arising from m Ψ appearing in the second term of Eq. ( 18)) via the t channel exchange of the heavy neutral leptons Ψ Regarding the fermionic DM possibility, which we focus on from now on, Ψ (±) 2 .Here in this work, to simplify our analysis, we assume that Ω k (k = 1, 2) is sufficiently heavy to allow its decay into the φ 2 Ψ nL final state, then implying that the lightest stable candidates are Ψ (±) 1 having a common mass m Ψ 1 .This case is particularly interesting because the Yukawa coupling that mediates DM annihilation also participates in the calculation of the µ parameter, allowing us to correlate the cFLV, DM, and neutrino masses.In this case, for DM masses m Ψ 1 in the GeV-to-TeV ballpark and Yukawa couplings at the electroweak scale, DM could have been generated in the early Universe via the WIMP mechanism.The evolution of the DM number density n can be tracked with the help of the Boltzmann equation where ⟨σv⟩ is the total thermally-averaged annihilation cross section of a couple of DM particles into lighter states, n eq (T ) corresponds to the DM number density in equilibrium at a temperature T , given by in the nonrelativistic limit, and H corresponds to the Hubble expansion rate of the Universe assumed to be dominated by SM radiation. 1  In the early Universe, pairs of non-relativistic DM Ψ (±) 1 can annihilate mainly into right-handed neutrinos N nR through the t-channel exchange of φ 1 .Later, N nR decays into Higgs and lepton doublets.The corresponding squared amplitude for the annihilation is given by and hence the thermally-averaged annihilation cross section ⟨σv⟩ is where K 1 is the modified Bessel function and σ R (s) is the reduced cross section It is interesting to note that Ψ if their mass difference is typically smaller than ∼ 20% [71].In that case, an additional Boltzmann equation for the other Ψ (±) 2 must be taken into account.Here, however, we assume a larger splitting so that coannihilation processes can safely be ignored.
To match the entire observed DM relic density, it is required that where n 0 /s 0 is the asymptotic value of the DM yield at low temperatures, s 0 ≃ 2.69 × 10 3 cm −3 is the present entropy density [72], ρ c ≃ 1.05 × 10 −5 h 2 GeV/cm 3 is the critical energy density of the Universe, and Ωh 2 ≃ 0.12 is the observed DM relic abundance [69].

Interplay between DM, cLFV, LNV, 0νββ and leptogenesis
In this section, we discuss the implications of the considered model for DM, cLFV, LNV, 0νββ, and leptogenesis.To this end, we performed a random scan in the ranges m Ψ 1R ∈ [50, 500] GeV, 1 DM freeze out in nonstandard cosmological scenarios or during a low-temperature reheating era is also possible, but will not be discussed here [70].We show in Fig. 2 the allowed parameter space in the |y ν |-m N 1 plane consistent with the constraints imposed by the measured value of the DM relic abundance and by the experimental neutrino oscillation data.We consider both normal and inverted neutrino mass hierarchies, corresponding to the left and right panels of Fig. 2, respectively.The consistency with the experimental values of the neutrino mass squared differences and with the constraints arising from the charged-lepton violation requires Yukawa couplings y ν ik associated with the Dirac submatrix smaller than 10 −2 .On the other hand, the gray points in Fig. 2 correspond to large mixing angles between active and heavy neutrinos that give rise to unacceptably large rates for the lepton flavor-violating process, higher than their corresponding upper experimental limits, and thus they are excluded by cLFV constraints.should have the form [73,74] where M ν diag defined in Eq. ( 29), and U PMNS the PMNS leptonic mixing matrix (without unitarity violations).It proves convenient to employ the µ parameterization in this case to avoid the apparent non-decoupling behavior with the mass of the heavy neutrinos when using the Casas-Ibarra parameterization, which would lead to a constant value of the cLFV rates for large m N 1 [73].The intensity of the colors in Fig. 3 corresponds to different values of the baryon asymmetry parameter.As shown in Fig. 3, there are several points corresponding to the decay rate µ → eγ lower than its upper experimental limit of 4.2 × 10 −13 and that yields values for the baryon asymmetry parameter consistent with its experimental value.
Figure 4 shows the baryon asymmetry parameter as a function of the trace of the Majorana submatrix µ for scenarios of normal (left panel) and inverted (right panel) neutrino mass hierarchies.
The consistency with the neutrino oscillation data is ensured in this case by the Casas-Ibarra parameterization.The gray points in Fig. 4 are excluded by the constraints arising from cLFV.
The purple points correspond to the values of the baryon asymmetry parameter Y ∆B within the experimentally allowed range at 3σ CL.As shown in Fig. 4, this model is consistent with the experimental range of the baryon asymmetry parameter, provided that the entries of the Majorana submatrix µ are in the keV to MeV range.
Having studied the parameter space favored by DM and leptogenesis separately, we next combine this information to show that both cosmological and electroweak precision data can be successfully accommodated at the same time.Again, the experimental neutrino oscillation data are also guaranteed by the Casas-Ibarra parameterization.spond to rates for cLFV processes larger than their experimental upper limits and are therefore excluded.
The correlation among the cLFV observables for the scenarios of normal and inverted neutrino mass hierarchies is shown in the left and right panels of Fig. 6, respectively.As indicated in Fig. 6, an increase of the µ → eee decay rate produces larger values for the µ → e conversion rate as well as for the µ → eγ branching ratio.Current bounds are shown with full lines, while projections are shown with dashed lines.All the points shown in Fig. 6 agree with the experimental values of the DM relic abundance, BAU, and neutrino oscillation.Figure 6 shows a nice interplay between cosmological and particle physics observables.Although it shows that the cLFV constraints already exclude many otherwise viable points to account for DM and leptogenesis, it also shows that the upcoming cLFV experiments will be able to probe a large portion of the remaining ones.
It is worth mentioning that despite the large number of similarities in the different correlation plots among several observables, in the plot of Fig. 5, which displays the correlation between the trace of the Majorana submatrix µ and the mass m N 1 of the lightest pseudo-Dirac lepton, in the region of large m N 1 values, there are more gray points (points excluded by the cLFV constraints) for the case of inverted hierarchy than for the normal neutrino mass spectrum.

Other three-loop inverse seesaw models
For completeness, here we present two more versions of the model setup leading to three-loop ISS neutrino mass generation.We do not explore the phenomenological and cosmological implications of these models here and leave this task to future work.

Model 2
This second example is based on the same gauge group of Eq. ( 1), but with a different spontaneous symmetry breaking scheme The global U (1) ′ symmetry is spontaneously broken at the TeV scale by the VEV of ⟨σ⟩.Other new scalars with non-trivial charges Z 2 do not acquire VEVs to maintain this unbroken sym-Field metry.The charges for the fields under the symmetry of Eq. ( 55) are shown in Table 1.The neutrino Yukawa interations invariant with respect to this group are Note that the mass splitting between the real and imaginary parts of the scalar singlet φ is generated only at the two-loop level and therefore is small.With this at hand and assuming that the neutral leptons Ψ kR are heavy enough, we can dynamically generate the small Majorana mass terms (µ) nk N nR N C kR at the three-loop level, as shown in Fig. 7.The preserved Z 2 discrete symmetry allows for a stable scalar or fermionic DM candidate corresponding to the lightest electrically neutral particle odd under Z 2 .The fermionic DM candidate can be the lightest of the two Ψ kR states.Field are given by the Lagrangian density

Model 3
In our model the scalar fields σ, χ, ϕ develop non-zero VEVs ⟨σ⟩, ⟨χ⟩, ⟨ϕ⟩ breaking the model symmetry in two ways: for and for In both cases, the residual Z 6 originates from the spontaneously broken U (1) B−L and survives after the electroweak symmetry breaking.In the second case, v χ > v σ , in the first step of symmetry breaking by v χ = ⟨χ⟩ there appears another residual symmetry Z ′ 6 ⊂ U (1) ′ , which is, however, completely broken in the second step by v σ .
Note that the Z 2 parity is preserved at low energies, offering stable scalar or fermionic DM candidates.In this model, the scalar fields φ, ρ, and ζ do not acquire VEVs and as they have non-trivial Z 2 charges the lightest of their real and/or imaginary parts can be a viable DM candidate.
In both scenarios ( 58) and ( 59), the LNV required for the generation of Majorana neutrinos masses arises due to the spontaneous breaking of U (1) B−L by ⟨χ⟩ ̸ = 0. Thus, in this model, there is an SM singlet majoron ∼ Im χ.The LNV comes into play only through the Majorana masses at tree level ∼ y⟨χ⟩ of the seesaw mediators Ψ iR , charged under U (1) B−L .The LNV is then transferred to active neutrinos at the loop level.Similarly to Model 1, the U (1) ′ ⊗ Z 2 symmetry forbids active neutrino masses at one-and two-loop levels but allows them via ISS with the Majorana mass term (µ) nk N nR N C kR dynamically generated at the three-loop level, as shown in Fig. 8.

Conclusions
We have constructed three models where the tiny masses of the active neutrino are radiatively generated from an inverse seesaw mechanism at three-loop level, thanks to preserved discrete symmetries.In these models, the Majorana mass submatrix corresponding to the lepton number violating (LNV) µ parameter is radiatively generated via the three-loop level exchange of neutral leptons and gauge singlet scalars.Such a three-loop suppression for the LNV µ parameter allows values for the entries of the Majorana mass submatrix in the keV to MeV range.Focusing on one of these models, where the SM particle content is extended by the inclusion of several neutral leptons and gauge singlet scalar fields and the SM gauge symmetry is extended with the global symmetry U (1) ′ ⊗ Z 2 , we analyzed in detail the scalar sector and explored the implications of the model in neutrino masses and in the lepton sector phenomenology for both scenarios of normal and inverted neutrino mass hierarchies.In such a model, the Z 2 symmetry is preserved, and the U (1) ′ symmetry is spontaneously broken down to a residual conserved Z 4 symmetry.
These discrete symmetries are crucial for ensuring the radiative nature of the inverse seesaw mechanism as well as the stability of fermionic and scalar dark matter candidates, particles that mediate the three-loop level radiative seesaw mechanism that yields the LNV µ parameter.
We found that the considered model successfully complies with the constraints imposed by the neutrino oscillation experimental data, neutrinoless double beta decay, dark matter relic density, charged lepton flavor violation, electron-muon conversion, and provides essential means for efficient low-scale resonant leptogenesis.We have also shown that in the model considered for the analysis, the charged lepton flavor violating decays µ → eγ, µ → eee as well as the electron-muon conversion processes get sizable rates, which are within the reach of sensitivity of the forthcoming experiments.Finally, we found that consistency with the measured value of the baryon asymmetry of the Universe requires that the entries of the Majorana submatrix corresponding to the LNV µ parameter acquire values in the keV to MeV range.
A Form factors µ-e conversion and BR(µ → eee) In this appendix, we provide the analytical expressions for the form factors and the loop functions entering in the amplitudes of the muon-electron conversion process.
In the following, we provide the expressions for the radiative decay and ℓ β → 3ℓ α decays (the full expression for the most general case of the 3-body decay can be found in Ref. [64]), closely following the notation of Ref. [75].The rates for the radiative and three-body decays in the SM extended via heavy RH neutrinos, are given by [65] BR(ℓ β → ℓ α γ) = α 3 w s where M W is the W -boson mass, α w = g 2 w /4π denotes the weak coupling, s w the sine of the weak mixing angle, and m β (Γ β ) the mass (total width) of the decaying charged lepton of flavor β.The form factors G βα γ , F βα γ , and F β3α box are given by [64,65] In the above expressions, the sums are understood to be taken over all neutral mass eigenstates.
C ij is defined as The neutrinoless conversion rate is given by [65] CR(µ − e, in which Γ capt denotes the capture rate for the nucleus N, with D, V (p) and V (n) corresponding to nuclear form factors (see Ref. [63]), and e is the unit electric charge.The above form factors are given by [64,65] to which one must add Here, x q = m 2 q /M 2 W and V is the Cabibbo-Kobayashi-Maskawa quark mixing matrix.The different loop functions (with arguments defined as x i = m 2 i /M 2 W ) are summarized below.

Photon dipole and anapole functions
Z-penguin: two-and three-point functions

M
are rotation matrices that diagonalize the physical mass matrices, M ν , M kL(k = 1, 2) are the three active neutrinos and four exotic neutrinos, respectively.In order to successfully reproduce the neutrino oscillation experimental data, the Dirac submatrix m νD , in the basis of diagonal SM charged lepton mass matrix, should have the following form[45,[49][50][51][52] 1,2) and neutrinos ν i and N (±) i (i = 1, 2, 3), respectively.Moreover, φ 1DM can also annihilate in neutrino pairs Ψ r via the t channel exchange of Ω k and Ψ (±) k (k, r = 1, 2), respectively.These annihilation channels will contribute to the DM relic density, which can be successfully accommodated for appropriate values of the scalar DM mass and couplings of the aforementioned scalar portal interactions.

Figure 2 :
Figure 2: Allowed values in the |y ν |-m N1 plane consistent with the DM relic abundance and neutrino oscillation data assuming normal (inverted) ordering in the left (right) panel.The small m N1 leads to Yukawa couplings smaller than 10 −2 , to generate the correct active neutrino masses.Higher Yukawa couplings, although still allowed by oscillation data, are excluded by cLFV constraints (gray points).

Figure 3
Figure3shows the values obtained for the branching ratio for the decay of µ → eγ for random values of the entries of the Dirac submatrix m νD and different masses of the lightest pseudo-Dirac neutral lepton N 1 with m N 2 fixed to be equal to 2 TeV and 5 TeV, for the left and right panels, respectively.Consistency with the neutrino oscillation data in Fig.3is ensured by using the µ parameterization, which implies that to successfully reproduce the experimental values of the neutrino mass squared splittings, leptonic mixing parameters, and the leptonic Dirac CP phase, the Majorana submatrix µ, in the basis of diagonal SM charged lepton mass matrix,

Figure 3 :
Figure 3: BR(µ → eγ) as a function of m N1 for fixed values of m N2 at 2 TeV and 5 TeV.The points reproduce neutrino oscillation data (using the µ parameterization), and are colored according to the size of the baryon asymmetry parameter Y ∆B .Points that have Y ∆B of the order of the correct experimental value typically pass the constraint of µ → eγ, shown by the horizontal black line.

Figure 5
shows solutions in the plane Tr[µ]-m N 1 that reproduce the DM relic abundance and the BAU simultaneously.The colors correspond to the values of the y ν e1 Dirac Yukawa coupling.The gray points in Fig. 5 corre-

Figure 4 :
Figure 4: Baryon asymmetry parameter Y ∆B as a function of the trace of the Majorana submatrix µ.All points comply with the DM relic abundance and neutrino oscillation data assuming normal (inverted) ordering in the left (right) panel.The purple points have the correct Y ∆B within 3σ, while the gray points are excluded by cLFV.The values of µ in the keV-MeV range are favored by leptogenesis.

Figure 5 :
Figure 5: Correlation between the trace of the Majorana submatrix µ and m N1 .Points comply with the BAU, DM relic abundance and neutrino oscillation data assuming normal (inverted) ordering in the left (right) panel.The points are colored according to the size of y νe1 .Gray points are excluded by cLFV constraints.

Figure 6 :
Figure 6: Correlation among cLFV observables.Points comply with the BAU, DM relic abundance and neutrino oscillation data assuming normal (inverted) ordering in the left (right) panel.Current bounds are shown as full lines, while projections are shown as dashed lines.

Model 3
is similar to Model 2, with the important difference that an additional U (1) B−L gauge symmetry is introduced, together with U (1) ′ ⊗ Z 2 symmetry.This model, compared to Model 1, has one extra scalar and one RH Majorana neutrino, which means that in total it contains five electrically neutral scalar singlets φ, ρ, ζ, σ, χ, and six RH neutrinos ν kR , N kR , Ψ iR (with k = 1, 2 and i = 1, 2, 3).The increased number of fermions comes from the requirement of cancelation of chiral anomalies.The quark and lepton fields and their assignments are shown in Table3.Note that, in contrast to Models 1 and 2, the quark and lepton sectors are correlated here by the condition of anomaly cancelation.The corresponding neutrino Yukawa interactions