Dark matter in the Scotogenic model with spontaneous lepton number violation

Scotogenic models constitute an appealing solution to the generation of neutrino masses and to the dark matter mystery. In this work we consider a version of the Scotogenic model that breaks lepton number spontaneously. At this scope, we extend the particle content of the Scotogenic model with an additional singlet scalar which acquires a non-zero vacuum expectation value and breaks a global lepton number symmetry. As a consequence, a massless Goldstone boson, the majoron, appears in the particle spectrum. We discuss how the presence of the majoron modifies the phenomenology, both in flavor and dark matter observables. We focus on the fermionic dark matter candidate and analyze its relic abundance and prospects for both direct and indirect detection.


Introduction
The origin of neutrino masses and the nature of the dark matter (DM) component of the Universe are two of the most relevant open questions in current physics.Regarding the former, neutrino oscillation experiments have robustly established the existence of non-zero neutrino masses and lepton mixings.In fact, some of the oscillation parameters have been already determined with great accuracy [1].Since the Standard Model (SM) of particle physics does not include a mechanism for the generation of neutrino masses, an extension is called for.Similarly, the Planck collaboration has determined that about 27% of the energymatter content of the Universe is in the form of DM [2].It is often assumed that the DM is made of particles, but no state in the SM spectrum has the required properties to play such a role.Again, this motivates the exploration of scenarios beyond the SM.
There are many neutrino mass models.Among them, radiative models (models that induce neutrino masses at the loop level) are particularly well motivated, since they naturally explain the smallness of neutrino masses due to the loop suppression.Pioneer work on radiative models can be found in [3][4][5][6], while for a recent review we refer to [7].Furthermore, tree-level contributions to neutrino masses are often forbidden by a conserved Z 2 symmetry which, in addition, stabilizes the lightest Z 2 -odd state.Provided it has the correct quantum numbers and can be produced in the early Universe in the correct amount, this state is a valid DM candidate.Therefore, radiative models offer a good solution to simultaneously address the origin of neutrino masses and the DM problem.A prime example of such class of models is the Scotogenic model [8].This model introduces an additional SU(2) L doublet, η, and three generations of fermion singlets, N , all charged under a Z 2 parity.These ingredients suffice to generate neutrino masses at the 1-loop level and provide a viable DM candidate.
In most neutrino mass models, neutrinos are Majorana fermions.This is precisely the case of the Scotogenic model.In this class of models, U(1) L -where L stands for lepton number -is broken in two units.The breaking can be explicit, due to the presence of lepton number violating parameters in the Lagrangian, or spontaneous, if the minimum of the scalar potential of the model does not preserve the symmetry.In the standard Scotogenic model [8] the breaking is explicit.In contrast, in this paper we consider a version of the Scotogenic model that breaks lepton number spontaneously.This is achieved by extending the particle content of the model with an additional singlet scalar, denoted as σ, which acquires a non-zero vacuum expectation value (VEV) and breaks the global U(1) L symmetry.As a consequence, the spectrum of the theory contains a massless Goldstone boson, the majoron, J [9][10][11][12][13].This state leads to novel phenomenological predictions, both in flavor observables (due to the existence of new channels such as ℓ α → ℓ β J) and in the DM sector (due to the existence of new processes in the early Universe).
Several works combining spontaneous lepton number breaking with the Scotogenic generation of fermion masses can be found in the literature.We highlight [14], which also studies the DM phenomenology of the Scotogenic model with spontaneous lepton number violation.We build upon this previous work and go beyond it in several ways.First of all, our analysis takes into account a wide variety of lepton flavor violating (LFV) constraints, including processes that involve the majoron either virtually or as a particle in the final state.We confirm the results of [14], but also discuss in further detail some aspects of the DM phenomenology of the model.A high-energy extension of the Scotogenic model featuring a massless majoron was also introduced in [15], while Ref. [16] proposes a model with spontaneous lepton number violation that induces a small 1-loop mass for a dark Majorana fermion à la Scotogenic.The authors of [17] studied electroweak baryogenesis in an extended Scotogenic scenario including a majoron, whereas the possible Scotogenic origin of the small lepton number violation of the inverse seesaw was discussed in [18].Finally, the spontaneous breaking of a gauged version of lepton number in a Scotogenic scenario was considered in [19].
The rest of the manuscript is organized as follows.We present the model in Sec. 2, where we define its basic ingredients, discuss its scalar sector and the generation of Majorana neutrino masses and briefly comment on the possible DM candidates.The most important experimental bounds that constrain our scenario are discussed in Sec. 3, while the results of our numerical study are presented in Sec. 4. Finally, we summarize and draw our conclusions in Sec. 5.

The model
We consider a variant of the original Scotogenic model.The SM particle content is extended by adding the SU(2) L scalar doublet η, the scalar singlet σ and three generations of fermion singlets N .The scalar doublets of the model can be decomposed into SU(2) L components as Here H is the usual SM Higgs doublet.We impose the conservation of a global U(1) L symmetry which can be identified with lepton number.Finally, we also introduce the usual Z 2 parity of the Scotogenic model, under which N and η are odd while the rest of the fields are even. 1 The particle content of the model and the representations under the gauge and global symmetries are summarized in Table 1.
The most general Yukawa Lagrangian, involving the new particles compatible with all symmetries, can be written as where y and κ are 3 × 3 matrices.In the following, we take κ to be diagonal without loss of generality.The most general scalar potential is given by where m2 H , m 2 η and m 2 σ are parameters with dimension of mass 2 and the rest of the parameters are dimensionless.

Symmetry breaking and scalar sector
We will assume that the scalar potential parameters are such that a minimum is found for the configuration Here v ≈ 246 GeV is the usual electroweak VEV.This vacuum preserves the Z 2 parity, which remains a conserved symmetry.In contrast, lepton number is spontaneously broken and a Table 1: Particle content of the model and their representations under the gauge and global symmetries.q L , ℓ L , u R , d R , e R and H are the usual SM fields.
Majorana mass term for the N singlets is induced, with The tadpole equations obtained by minimizing the scalar potential are given by Assuming the conservation of CP in the scalar sector, one can split the neutral scalar fields in terms of their real and imaginary components as The η R and η I fields do not mix with the rest of scalars due to the Z 2 parity.In this case, the scalar potential contains the piece where z = {H 0 , σ} and M 2 R and M 2 I are the 2 × 2 CP-even and CP-odd squared mass matrices, respectively.One finds and One can now use the tadpole equations in Eqs. ( 6)- (7) to evaluate these matrices at the minimum of the scalar potential.We obtain while the CP-odd mass matrix becomes identically zero as expected, since it has to provide two massless states: the unphysical Goldstone boson z that becomes the longitudinal component of the Z boson and a physical massless Goldstone boson associated to the spontaneous breaking of the lepton number, the majoron (J).Therefore, since σ is a gauge singlet field one can make the identification The CP-even states {S H , S σ } mix leading to two massive states, h 1 and h 2 as follows: where O is the 2 × 2 orthogonal matrix which diagonalizes the CP-even mass matrix, such that and the mass eigenvalues are given by One of the two scalar masses has to be associated with the ∼ 125 GeV SM Higgs boson and an additional CP-even state is present in the spectrum.The angle α is the doublet-singlet mixing angle and is given by We focus now on the Z 2 -odd scalars.The masses of the CP-even and CP-odd components of η 0 are given by thus as in the usual Scotogenic model the mass difference between η R and η I is controlled by the λ 5 coupling.Finally, the mass of the charged scalar fields η ± turns out to be

Neutrino masses
Neutrino masses are induced at the 1-loop level, in the same way as in the standard Scotogenic model, as shown in Fig. 1.One finds the 3 × 3 neutrino mass matrix where m η R and m η I are the η R and η I masses, respectively, and m 2 N b are the diagonal elements of M N .We note that neutrino masses vanish for m η R = m η I .This is consistent with the fact that m η R − m η I ∝ λ 5 , and in the limit λ 5 → 0 a conserved lepton number can be defined.This allows one to assume λ 5 ≪ 1 in a natural way [20].

Dark matter
The lightest Z 2 -odd state is completely stable and can, in principle, be a good DM candidate.In this model, as in the standard Scotogenic model, this role can be played either by the lightest N state or by a neutral η field (η R or η I , depending on the sign of λ 5 ).In this work we will concentrate on the fermion DM and thus consider N 1 , the lightest singlet fermion, to be our DM candidate.

Constraints
Several experimental and theoretical constraints will be considered in our numerical analysis.

Boundedness from below
We demand the scalar potential to be bounded from below, which implies the following set of conditions [21]:

Higgs boson production and decays
In our model, all Higgs boson production cross-sections at the LHC are suppressed with respect to the SM by c 2 α , where c α = cos α and α is the mixing angle in the CP-even scalar sector.In addition, Higgs decays are also affected in two ways.First, the rates of all visible Higgs decay channels are universally reduced due the abovementioned mixing.And second, new decay channels are available.The Higgs boson may decay invisibly to a pair of majorons or to a pair of DM particles, h → J J and h → N 1 N 1 .The former will always be kinematically available, since the majoron is massless, whereas the latter requires m N 1 ≤ m h /2.The CMS collaboration has searched for invisible Higgs boson decays at the LHC [22], assuming a completely SM-like Higgs boson production through vector boson fusion.Therefore the limit derived in [22] translates into c 2 α BR(h → invisible) < 0.19 at 95% C.L..
A proper phenomenological analysis must take into account both Higgs production and decays, including visible and invisible ones.In fact, the recent analysis [23] has clearly shown that the strongest constraints on the parameter space of our model are obtained by combining the bounds from visible and invisible Higgs decays.In particular, Figure 9 of this reference shows the limits obtained for our scenario.These are the constraints that will be considered in our numerical analysis.

Electroweak precision data
Bounds from electroweak precision data can also be used to contrain the parameter space of our model.In particular, the oblique parameters S, T and U [24] are known to capture the effect of heavy new fields affecting the gauge boson propagators.Their current determination is in good agreement with the SM expectations, although there is some room for new physics.In our analysis we considered the bounds [25] S = −0.01 ± 0.10 , (26) T = 0.03 ± 0.12 , (27) U = 0.02 ± 0.11 . (28)

Neutrino oscillation data
All the parameter points considered in our analysis comply with the constraints from neutrino oscillation experiments.This is guaranteed by means of a modified Casas-Ibarra parametrization [26], properly adapted to the Scotogenic model [27][28][29], which allows us to express the y Yukawa matrix as Here Λ is a matrix defined as Λ = diag(Λ b ), with while R is an orthogonal matrix (R T R = RR T = I), generally parametrized by three complex angles.Finally, U is the unitary matrix that brings m ν to diagonal form as ) the neutrino physical masses.The entries of the unitary matrix U as well as the neutrino squared mass differences are measured in neutrino oscillation experiments.Our analysis will use the results of the global fit [1].

Majoron diagonal couplings to charged leptons
The interaction Lagrangian of majorons with charged leptons can be written as [30] where ℓ α,β are the standard light charged leptons and P L,R are the usual chiral projectors.The S L,R couplings are induced at the 1-loop level in our model, as shown in [15].The diagonal S ββ = S ββ L + S ββ * R couplings are purely imaginary, due to the fact that majorons are pseudoscalar states, and are strongly constrained due to their potential impact on astrophysical observations.Large couplings to electrons or muons are excluded since they would lead to an abundant production of majorons in dense astrophysical media and an efficient cooling mechanism [31][32][33][34][35].The authors of [33] used data from white dwarfs to set the bound while the supernova SN1987A was considered in [35] to establish the limit2 Finally, there are also laboratory bounds on the majoron diagonal couplings to charged leptons.In [30], the results of the OSQAR experiment [36], a light-shining-through-a-wall experiment, were used to find the approximate bounds S ee ≲ 10 −7 and S µµ ≲ 10 −5 .We note that these are clearly less stringent than the bounds obtained from astrophysical observations.

Lepton flavor violation
As in most neutrino mass models, LFV is a powerful constraint that strongly restricts the allowed parameter space of our model.Several processes will be considered in our analysis: • The radiative decays ℓ α → ℓ β γ, which turn out to be the most constraining ones in most neutrino mass models.In particular, the MEG experiment restricts the µ → eγ branching ratio to be smaller than 4.2 × 10 −13 [37].We also consider the analogous limits on τ LFV decays [25], but they are less stringent.
• The 3-body decays ℓ α → ℓ β ℓ γ ℓ γ , with β = γ and β ̸ = γ.In this case we follow [38] and include the usual photon penguin contributions as well as other usually less relevant contributions, such as box diagrams.Majoron mediated contributions are also included, using the results derived in [30].
• The decays ℓ α → ℓ β J with the majoron in the final state are also considered, as they constrain the off-diagonal S L,R couplings directly.For instance, the null results obtained in the search for the decay µ → eJ at TRIUMF [39] can be translated into the bound |S eµ | < 5.3 × 10 −11 [40].We use the analytical expressions for the majoron off-diagonal couplings to charged leptons in the Scotogenic model found in [15].
• µ − e conversion in nuclei, again following the analytical results in [38].

Majoron coupling to neutrinos
Constraints on the couplings of the majoron to neutrinos can be derived from astrophysics, since the majoron can have a significant impact on the explosion and cooling of supernovae (see e.g.[41]) and also from cosmic microwave background data [42].Laboratory experiments searching for neutrinoless double beta decays (e.g.[43]) and possible effects on meson and lepton decays [44] also set limits on the magnitude of neutrino-majoron couplings.Among these, the most stringent ones are those derived from astrophysics, with constraints in the ∼ 10 −7 ballpark.In our model, the interaction of majorons with neutrinos arises at the 1loop level, similarly to the interaction with charged leptons, and the corresponding couplings are thus expected to be of the same order.Since the constraints on the couplings to charged leptons are orders of magnitude more stringent, we can safely ignore the constraints on the couplings to neutrinos in our analysis.

Numerical results
We now proceed to discuss the results of our analysis.To perform the numerical scan, we have first implemented the model in SARAH (version 4.11.0)[45], a Mathematica package for the analytical evaluation of all the information about the model. 3With this tool, we have created a source code for SPheno (version 4.0.2) [47,48], thus allowing for an efficient numerical evaluation of all the analytical expressions derived with SARAH.We have also Table 2: Values of the main input parameters for the numerical scan.
computed several observables of interest in our model, including the lepton flavor violating ones, both analytically and with the help of FlavorKit [49], for an in-depth cross-check of their expressions.Finally, we have used micrOmegas (version 5.0.9)[50] to obtain the main DM observables, namely the DM relic density and direct and indirect detection predictions.
As already mentioned, while both the scalar η R,I and the fermions N i are, in principle, viable DM candidates in this model, in our analysis we focus on the lightest Majorana fermion N 1 as the main component of the DM.We summarize our choice of parameters for the numerical scan in Tab. 2.Moreover, λ 1 is fixed by the condition of requiring m h 1 = 125 GeV.In some numerical scans we have fixed the value of m h 2 = 500 GeV or the value of m 2 η such that m η R,I ,η ± −m N 1 ≲ 20 GeV, as we will discuss in more detail below.With our choice of parameters, all the parameter points considered in our numerical scans easily pass the bounds from the S, T, U parameters.We note that with our choice of λ 4 and λ 5 values the mass splitting between the neutral and charged components of the η doublet is small (see Eqs. ( 18) and ( 19)).We have explicitly checked that only for λ 4 ≳ 2 the electroweak precision bounds become relevant, but we did not explore this region of parameter space in our scans.Since we want to focus on N 1 as the DM candidate, we further require that m N 1 < m N 2,3 , m η R,I .We have chosen normal hierarchy for the neutrino spectrum and considered the best-fit values for the neutrino oscillation parameters found by the global fit [1].Finally, the three angles in the orthogonal R matrix are assumed to be real and taken randomly in our numerical scans. 4e first show in Fig. 2 the relic abundance of N 1 , as a function of its mass.For this specific scan, we have fixed m h 2 = 500 GeV to highlight the s-channel annihilation of N 1 via h 2 .In this figure, grey points denote solutions either leading to overabundant DM or excluded by any of the constraints listed in Sec. 3, or where the spin-independent N 1 -nucleon elastic scattering cross section is excluded by the most recent data from the LUX-ZEPLIN experiment [52].Red points denote solutions which can reproduce the observed cold DM relic density, as they fall within the 3σ range obtained by the Planck satellite data [2], Ω N 1 h 2 = 0.120 ± 0.0036 (blue thin band).Solutions leading to underabundant DM (which would then require another DM candidate to explain the totality of the observed cold DM relic density) are depicted in blue.As can be seen from the plot, most of the solutions lead to overabundant DM, except for points falling in the following regions: (i) a resonant region where m N 1 ∼ m h 1 /2 ∼ 60 GeV, (ii) a second resonant region where m N 1 ∼ m h 2 /2 ∼ 250 GeV and (iii) a region of coannihilations at higher m N 1 .To explore in more detail the third, high-mass region, we performed a second numerical scan in which we have varied the mass difference ∆ = m η R − m N 1 in the [0, 20] GeV range.In such a way, we have enforced N 1 to be in the ∼ [100 − 3000] GeV region, where coannihilations with η R,I and η ± are very relevant, thus reducing the relic abundance of N 1 .Figure 3 shows this region in parameter space, in which the DM relic density is set by coannihilations.The color code is the same as in Fig. 2. Compared to the result of the previous scan (Fig. 2), we can see that if coannihilations are relevant, more viable solutions can be found in the m N 1 ∼ [100, 3000] GeV region.We clarify that ∆ < 20 GeV is not motivated by any symmetry argument, but just a convenient parameter choice to focus our numerical analysis on a region in which coannihilations are more effective.Finally, we should also make a comment about the region with light η states (m η R,I , m η ± ≲ 250 GeV).These states can be pair-produced at the LHC via Drell-Yan processes.However, our choice of ∆ implies a compressed spectrum with N , η R,I and η ± in a narrow window of just 20 GeV, thus implying soft leptons in the η → N ℓ final decay.Searches for this type of signal exist, although not dedicated to our specific scenario.For instance, the ATLAS collaboration looked for direct slepton production with a compressed spectrum in [53].This analysis assumes mass-degenerate 1st and 2nd generation sleptons decaying to flavor conserving final states, whereas our scenario contains only a copy of the η doublet (and not two), with both flavorconserving and flavor-violating decays.Therefore, the obtained limits are not applicable.Nevertheless, we note that some points with ∆ ∼ 10 GeV, where the experimental searches are more efficient, must be excluded.A detailed analysis including this constraint is clearly beyond the scope of our work and would not have any impact on our conclusions.Next we discuss the results for N 1 direct detection.In order to maximize the number of viable solutions, we focus again on the coannihilation region, and we show in Fig. 4 the spin-independent N 1 -nucleon elastic scattering cross section, σ SI , as a function of the DM mass, m N 1 .The cross section shown in this figure is weighted by the relative abundance ξ, defined as where Ω DM,Planck h 2 = 0.120 [2].We apply the same color code as in Fig. 2, that is red points indicate solutions explaining the totality of observed DM, while blue points denote underabundant DM.The plain green line and dashed area indicate the current most stringent limit from the LUX-ZEPLIN experiment (LZ-2022) [52], while the black dashed line denotes the constraint from XENON1T (XENON1T-2018) [54].Other (less stringent) constraints on σ SI apply from the liquid xenon experiment PandaX-II [55] and from liquid argon experiments like DarkSide-50 [56] and DEAP-3600 [57], although they are not shown here.Future facilities including XENONnT [58], DarkSide-20k [59], ARGO [59] and DARWIN [60,61] (see [62] for an overview) will be able to further inspect the parameter space of this model.
As for general reference, we further illustrate the expected discovery limit corresponding to the so-called "ν-floor" from coherent elastic neutrino-nucleus scattering (CEνNS) for a Ge target [63] (dashed orange line).The green area is already excluded by the LUX-ZEPLIN experiment (LZ-2022) [52], while the black dashed line denotes the constraint from XENON1T (XENON1T-2018) [54].The dashed orange curve indicates the expected discovery limit corresponding to the ν-floor from CEνNS of solar and atmospheric neutrinos for a Ge target [63].
Finally, we have explored the predictions for the velocity-averaged cross section of N 1 annihilation into gamma rays.These are among the most suitable messengers to probe DM via indirect detection.We focus once more on the coannihilation region, i.e. on the high-mass range m N 1 ∼ 0.1−2 TeV, where the annihilation channels The hadronization of the final-state gauge bosons and Higgs bosons will produce neutral pions, which in turn can decay into photons thus giving rise to a gammaray flux with a continuum spectrum which may be within reach of DM indirect detection experiments.While a detailed calculation of the gamma-ray energy spectra produced by the annihilation of two N 1 particles in this specific model should be performed, in order to correctly compute exclusion bounds from existing gamma-ray data, this is out of the scope of this work.However, we can notice that the main annihilation channels in this high-mass range include Higgs bosons in the final state.The gamma-ray energy spectrum from DM DM → h 1 h 1 annihilation channel is very similar to that from DM DM → W + W − at m N 1 ∼ 1 TeV (see for instance Fig. 15 of [67]).In the following, for the sake of simplicity, we will compare our predictions with bounds obtained assuming W + W − as the main annihilation channel, to get an overall idea of how current data can constrain the parameter space of this model.Charged cosmic rays can also be used to look for N 1 annihilations, even though their detection is more challenging due to uncertainties in the treatment of their propagation.For instance, AMS-02 data on the antiproton flux and the Boron to Carbon (B/C) ratio can be used to constrain the N 1 annihilation cross section [68][69][70].With some caveats concerning the astrophysical uncertainties on the p production, propagation and on solar modulation (see e.g.[71][72][73]), these bounds turn out to be stronger than gamma-ray limits from dwarf spheroidal satellite galaxies in some mass ranges.Following the same considerations as before, i.e. that the antiproton energy spectrum from DM DM → h 1 h 1 annihilation channel is very similar to that from DM DM → W + W − at m N 1 ∼ 1 TeV, we will compare our predictions to current limits on the N 1 annihilation cross section set by combination of p and B/C data of AMS-02 [68,69] assuming W + W − as the dominant annihilation channel.We show in Fig. 5 the N 1 total annihilation cross section -weighted by ξ 2 -versus its mass.The color code follow the same scheme as in Figs. 3, 4. We also depict the 95% C.L. upper limits currently set by the Fermi-LAT with gamma-ray observations of Milky Way dSphs (6 years, Pass 8 event-level analysis) [74] (red solid curve and shaded area) and from a combination of p and B/C data of AMS-02 [68,69] (green), both assuming N 1 N 1 → W + W − as main annihilation channel due to the considerations made before.We see that few solutions already fall within the region currently excluded by AMS-02 data.As already highlighted, while a dedicated analysis should be performed for this specific model, we can conclude that current p and B/C data may be already excluding a relevant part of the parameter space.Forthcoming data will allow to further probe N 1 as a DM candidate via its multi-messenger signals.
As in many scenarios for neutrino mass generation, LFV processes strongly restrict the available parameter space of the model.In addition to µ → eγ, very commonly considered in phenomenological studies, our model also leads to signatures with the majoron in the final state, like µ → eJ. Figure 6 shows BR(µ → eJ) as a function of BR(µ → eγ).Again, we have focused on the coannihilation region.We first notice that some parameter points are already excluded by the current experimental limits on these LFV branching ratios.However, one can also see that our numerical scan also finds many valid parameter points leading to very low values of both BR(µ → eγ) and BR(µ → eJ), clearly below the discovery reach of planned experiments.This is not surprising, since we take random R matrices in our numerical scans, hence accidentally finding parameter points with suppressed µ − e flavor violation.While a slight correlation among these two observables can be observed in Fig. 6, µ → eγ receives contributions from additional loop diagrams that do not involve the majoron.The two observables are hence independent.Interestingly, we find that BR(µ → eJ) is generally more constraining than BR(µ → eγ), although the difference is not very significant.

Summary and discussion
Most SM extensions aiming at an explanation of neutrino oscillation data consider Majorana neutrinos.This option breaks the accidental U(1) L lepton number symmetry of the SM in two units.If the breaking of lepton number is spontaneous, a Goldstone boson appears in the particle spectrum of the theory, the majoron.In this work we have analyzed the dark matter phenomenology of this scenario in the context of the popular Scotogenic model.Focusing on the fermionic DM candidate N 1 , we have found that it can explain the observed DM abundance in three regions of parameter space: (i) a resonant region where it annihilates via h 1 , with m N 1 ∼ 60 GeV, (ii) a second resonant region where s−channel annihilations via h 2 are relevant and (iii) a region of coannihilations at m N 1 ∼ 1 TeV.In particular, if coannihilations are relevant, more allowed solutions are found, either explaining the totality of DM or at least a sizeable part of it.While some of these solutions are already excluded by the recent LUX-ZEPLIN result, most of them are within the reach of near-future direct detection experiments.Interestingly, indirect detection searches seem to constitute another promising tool to further probe N 1 as a DM candidate via its multi-messenger signals, mainly gamma rays and antiprotons.All in all, the presence of the majoron and of a second Higgs open up the allowed parameter space of N 1 as DM, compared to the standard Scotogenic model.The majoron has an impact also on the phenomenology of LFV observables, as it leads to new interesting signatures, where it appears in the final state.Among these, we found that BR(µ → eJ) is generally more constraining than the most common BR(µ → eγ).Moreover, let us comment that the presence of a massless majoron may have relevant implications on the early-Universe cosmology.In particular, it can affect cosmological and astrophysical environments, and can contribute to ∆N eff .In principle, these bounds could be relaxed if the majoron acquires a small mass (for instance from quantum gravity considerations).In such a case, the majoron would decay before Big Bang nucleosynthesis and would not affect cosmological observations.However, let us notice that, in order not to alter the phenomenological analysis presented in this paper, the majoron mass should be smaller than the electron one and hence its only available decay channel would be into active neutrinos.On the other hand, the majoron can be produced from the Higgs decay, or the annihilation of N 1 or even via freeze-in through its small coupling with the active neutrinos.If it is massless and thermalizes, in order to avoid constraints from ∆N eff , one should require the majoron to decouple before T ∼ 0.5 GeV (see e.g.[75]) to avoid the current constraint from Planck.This can be easily obtained if λ Hσ 3 is set small enough (≲ 10 −5 ).In such a case, all majoron production channels through SM particles would be suppressed, and it could only be produced via interactions with N 1 , through the (sizeable) coupling κ.If this is the case, the majoron would freeze out at around the same time as N 1 , that is at T F ∼ m N 1 /20, thus not substantially contributing to ∆N eff .We have checked that by imposing λ Hσ 3 ≲ 10 −5 our results remain almost unchanged, with the only exception of the second resonance shown in Fig. 2.This region would disappear, due to the fact that a tiny λ Hσ 3 suppresses all vertices involving a (heavy or light) Higgs.
Finally, another interesting scenario consists in N 1 having very tiny couplings, so that it does not reach thermal equilibrium in the early Universe and it is instead produced via freeze-in.Such a production mechanism, yet together with the presence of the majoron, should also lead to some interesting phenomenology.We leave such analysis for a follow-up of this paper.

Figure 1 :
Figure 1: Generation of neutrino masses at the 1-loop level.In this diagram, η 0 denotes the real and imaginary components of the neutral component of the η doublet.

Figure 2 :
Figure 2: Relic abundance of N 1 as a function of m N 1 .Red points depict solutions in agreement with the cold DM measurement obtained from Planck data [2] (the blue thin band shows the 3σ interval) while blue points depict solutions leading to underabundant DM.Gray points are excluded by any of the constraints listed in Sec. 3 or due to an overabundant DM relic density.

Figure 3 :
Figure 3: Relic abundance of N 1 as a function of m N 1 in the coannihilation region, where ∆ ∈ [0, 20] GeV.Same color code as in Fig. 2.

44 XE NO N1 T -2 01 8 LZ -2 02 2 νFigure 4 :
Figure4: Spin-independent N 1 -nucleon elastic scattering cross section -weighted by the relative abundance -as a function of m N 1 .The green area is already excluded by the LUX-ZEPLIN experiment (LZ-2022)[52], while the black dashed line denotes the constraint from XENON1T (XENON1T-2018)[54].The dashed orange curve indicates the expected discovery limit corresponding to the ν-floor from CEνNS of solar and atmospheric neutrinos for a Ge target[63].

Figure 5 :
Figure5: N 1 total annihilation cross section as a function of m N 1 .The red and green lines refer to the corresponding 95% C.L. upper limits currently set by Fermi-LAT gamma-ray data from dSphs[74] and from the antiproton and B/C data of AMS-02[69], respectively.

Figure 6 :
Figure 6: BR(µ → eJ) as a function of BR(µ → eγ) in the coannihilation region, where ∆ ∈ [0, 20] GeV.Same color code as in Fig. 2. The horizontal and vertical lines correspond to the current experimental limits, discussed in Sec. 3.