The ScotoSinglet Model: A Scalar Singlet Extension of the Scotogenic Model

The Scotogenic Model is one of the most minimal models to account for both neutrino masses and dark matter (DM). In this model, neutrino masses are generated at the one-loop level, and in principle, both the lightest fermion singlet and the lightest neutral component of the scalar doublet can be viable DM candidates. However, the correct DM relic abundance can only be obtained in somewhat small regions of the parameter space, as there are strong constraints stemming from lepton flavour violation, neutrino masses, electroweak precision tests and direct detection. For the case of scalar DM, a sufficiently large lepton-number-violating coupling is required, whereas for fermionic DM, coannihilations are typically necessary. In this work, we study how the new scalar singlet modifies the phenomenology of the Scotogenic Model, particularly in the case of scalar DM. We find that the new singlet modifies both the phenomenology of neutrino masses and scalar DM, and opens up a large portion of the parameter space of the original model.


Introduction
ScSM, and comparing our results against the usual ScM. We also study how the parameter space expands with respect to the case of pure singlet or pure doublet DM via turning on/off the relevant couplings. The ScSM has some interesting features (mostly due to the presence of a scalar singletdoublet mixing, and a trilinear coupling between ϕ, Φ and the Higgs doublet H): new contributions to neutrino masses, 3 potential DM candidates (one of which is a mixture of singlet and doublet), and the possibility to maintain the Z 2 symmetry up to high energy scales. In light of these features, we perform, for the first time, a convergent global fit of the ScSM. As we will show, the presence of the singlet significantly opens up the allowed parameter space of the CP-even scalar DM, which now has a non-negligible singlet component. Due to the presence of the trilinear coupling with the singlet, there exists a splitting between the CP-even components. This naturally translates into a significant mass splitting between the lightest CP-even and CP-odd scalars, which allows to naturally evade the stringent DD limits.
The rest of the paper is organised as follows. In section 2, we introduce the ScSM. 3 The phenomenology of the ScSM and various theoretical/observational constraints that we impose are described in section 3. Sections 4 and 5 are devoted to our numerical analysis and results, respectively. Our conclusions are presented in section 6. A list of appendices provide supplementary information for understanding various expressions in the paper.

The ScotoSinglet Model (ScSM)
The new particle fields of the ScSM and their quantum numbers are presented in table 1. The Lagrangian for the Majorana fermion fields Ψ ≡ (ψ 1 , . . . , ψ N ) T is given by where M Ψ is an N × N diagonal mass matrix with real and positive values, and y Ψ is an N × 3 complex matrix of Yukawa couplings. Without loss of generality, we take m ψ 1 ≥ m ψ 2 ≥ . . . ≥ m ψ N . To reproduce the neutrino masses, N ≥ 2. Here we study the minimal case (N = 2). 4 The most general form of a Z 2 symmetric scalar potential is 3 (H † Φ) 2 + H.c. Here H ≡ 0, (v + h)/ √ 2 T is the SM Higgs doublet after electroweak symmetry breaking (EWSB). 5 Without loss of generality, λ HΦ,3 can be made real by performing a rotation of Φ. Although κ, in general, can be complex, we also take it as real in our study. We assume that the Z 2 symmetry is exactly preserved such that in the vacuum state, Φ = ϕ = 0. The electrically neutral and charged field components of the weak scalar doublet Φ are After EWSB, the physical masses of charged scalar φ + and CP-odd pseudoscalar A are In the (φ R , ϕ) basis, the squared mass matrix is non-diagonal, namely Notice from eq. (2.2) that κ controls the mixing between φ R and ϕ, and λ HΦ,3 the mass splitting between φ R and A. To diagonalise the mass matrix in eq. (2.5), we perform a rotation into the physical mass basis (η 1 , η 2 ) by  The mixing angle θ ∈ [0, π] and the quadrant is determined by the signs of a − b and c. For instance, a − b > 0 and c > 0 implies θ ∈ (0, π/4), whereas a − b < 0 and c > 0 implies θ ∈ (π/4, π/2). The physical scalar η 1,2 masses are given by Note that by convention, m η 1 ≥ m η 2 . The DM candidate can either be the CP-even scalar η 2 , the CP-odd pseudoscalar A or the lightest Majorana fermion ψ 2 (as m ψ 1 ≥ m ψ 2 by convention).

Phenomenology
After EWSB and rotation into the physical mass basis, the ScSM contains 12 free model parameters, namely 5 masses (1 mixing angle) : The remaining parameters in eq. (2.2) can be expressed as (see appendix A) We also have free parameters within the complex Yukawa matrices, for which we use a Casas-Ibarra parametrisation [38]. For N = 2 fermionic singlets, there are N = 2 real angles from the complex orthogonal Casas-Ibarra matrix R, and N − 1 = 1 Majorana phase in the PMNS matrix U (as one of the neutrinos is massless); see appendix B for more details. Thus, we can express the 6 complex Yukawa couplings in terms of the lowenergy neutrino oscillation data (3 masses, 3 angles, and 2 phases), the 2 heavy Majorana masses parameters, and the 2 real angles of the R matrix. These are considered as nuisance parameters in our study. In the following subsections, we discuss various theoretical/observational constraints that are imposed on the allowed model parameter space.

Perturbativity
From perturbativity arguments, we require all quartic couplings in the potential to satisfy In addition, we also impose the co-positivity conditions discussed in appendix C. As for the Yukawa couplings (derived parameters in our scan), we require where i = 1, 2 and α = e, µ, τ .

Naturalness
Using eq. (2.8) and expressing everything in terms of the physical scalar masses, the trilinear coupling κ can be bounded as where we have used m 2 η 1 m 2 η 2 in the last step. However, as the heavy neutral scalar masses are unknown, this upper bound is not very useful.
We can make use of the fact that the trilinear coupling gives a correction to the Higgs boson mass at the one-loop level, where the new scalars run in the loop. Up to factors of 2 and logs, we estimate For the fine-tuning to be under control, i.e., δm h /m h < , we demand that where we have assumed 1. Similar considerations were made in the Zee model [39]. If the fine-tuning condition is relaxed (i.e., 1), the upper limit on κ can also be relaxed.

Z 2 symmetry breaking at tree-level
If |κ| is much larger than the scalar masses, it can lead to a deeper minimum than the SM one, thereby breaking the Z 2 symmetry. Looking at different field directions, and using eq.
The above requirement, although more robust than the naturality one in eq. (3.7), turns out to be weaker, particularly for m Φ and m ϕ larger than the EW scale. In order for the model to be valid above the EW scale (such that the Z 2 symmetry is preserved), it is important to check that the RGE evolution does not only preserve m 2 Φ > 0 and m 2 ϕ > 0, but also that |κ| is not too large compared to the rest of the scalar masses. The latter requirement, however, is expected to be easily satisfied, as κ renormalises multiplicatively (see section 3.9).

Neutrino masses
Neutrino masses are generated at the one-loop level with the neutral fields running in the loop, see figure 1. From the Yukawa Lagrangian, eq. (2.1), and the scalar potential, eq. (2.2), we can see that the lepton number is violated by 2 units in the presence of y Ψ , M Ψ and either: 1. The quartic coupling λ HΦ, 3 . This contribution is same as in the ScM; see the left diagram in figure 1.
2. The square of the dimensionful trilinear coupling, κ 2 . This is the new contribution in the ScSM; see the right diagram in figure 1.
These are the parameter combinations that enter in the expression for the neutrino masses, namely 7 where the loop function is It is instructive to expand in the limit of small λ HΦ,3 and κ. For a > b and small κ, the mixing is θ ≈ 0, thus η 1 (η 2 ) is mostly doublet (singlet). In this case, the mass splitting between η 1 and A reads m 2 (3.11) 7 We believe there are two typos in the expressions of ref. [23]: i) a factor of π is missing in the denominator; and ii) an extra contribution to the scalar masses after EWSB is missing -this is due to the λHΦ,1 term in the potential (λ4 in the notation of ref. [23]).
For m η 1 m ψ k m η 2 , we get Alternatively, for b > a, θ ≈ π/2 and thus η 1 (η 2 ) is mostly singlet (doublet). The mass splitting between η 2 and A is m 2 whereas for m η 1 m ψ k m η 2 , we get Note how the Scotogenic-like contribution (proportional to λ HΦ,3 v 2 ) is always suppressed by the fermion or doublet mass, depending on which one is the heaviest state in the spectrum.

Integrating-out the heavy scalar singlet
If the singlet ϕ is the heaviest particle in the spectrum, e.g., m 2 , it can be integrated out (see ref. [40] for a similar study, and ref. [41] for the example of a charged scalar singlet). For scalar DM and θ π/2, the DM candidate η 2 is mainly doublet and this is a good approximation. Assuming that the weak and mass eigenstates are similar, i.e., for small mixing (small θ), after integrating out the singlet at tree-level before EWSB, we obtain the following expression for the scalar potential at dimension-6: This agrees with our expectation from considerations of lepton number violation. We see that neutrino masses can be suppressed either by small λ HΦ,3 and κ, or by cancellations among these terms. Notice that the last term in eq. (3.16) gives a contribution to neutrino masses via a dimension-7 Weinberg-like operator, which is expected to be suppressed compared to the usual dimension-5 one (the rest of the terms). In our numerical scan, we check that the combination in eq. (3.16) is indeed fixed (with the scale set by the neutrino masses). Let us elaborate a bit more on the threshold corrections to λ HΦ,2 and λ HΦ,3 at the scale of the singlet ϕ. Interestingly, the threshold effects increase their values for high energies. The stability conditions, however, must to be applied differently above and below m ϕ . This effect has been used to keep the Higgs potential stable up to high energy scales in ref. [42]. Notice also that the signs of λ HΦ,3 and κ are preserved under the RGE flow due to their multiplicative renormalization. The behaviour in eq. (3.16) can also be understood from the scalar mass matrix, see eqs. (2.5) and (2.6). In the symmetric Higgs phase (scales above the EW scale) and for m 2 ϕ m 2 Φ , the neutral scalar masses are see-saw like, i.e., from eqs. (2.9), we have We see that the lightest state (by convention η 2 ) is mainly doublet (Φ), while the heaviest state (η 1 ) is mainly singlet (ϕ), corresponding to a mixing angle of θ ≈ π/2. This repulsion of mass eigenvalues is the same effect that we see in the quartic couplings λ HΦ,2 and λ HΦ,3 , see eq. (3.15).

Electroweak precision tests
The new particles in the ScSM contribute to the W ± and Z boson self-energies. These contributions are parametrised by the oblique S, T and U parameters [43,44]. The strongest constraint comes from the T parameter, which bounds the mass splitting of the neutral and charged scalars. It is given by [45] where the loop function (symmetric in x and y) is The S parameter is given by where The Passarino-Veltman function B 22 [46] (symmetric in the last two arguments) is in d space-time dimensions and γ E 0.577 is the Euler-Mascheroni constant. We use the compact analytic expressions from appendix B of ref. [39]. Our expressions agree with the ones for the Inert Doublet Model (IDM) [47] in the appropriate limits. The oblique parameters are constrained from the global electroweak fit [48], along with the following correlation matrix:

Higgs decay into di-photons
The coupling between the SM Higgs h and charged scalar φ + in eq. (2.2) modifies the decay rate of the h → γγ process. The ratio with respect to the SM value is [49][50][51] where j = φ + , t and W , and f = arcsin 2 ( √ x) for m h < 2m i .

LFV process Upper limit (90% CL)
Ref. Table 2: Upper limits on the branching ratios/conversion rates of various lepton flavour violation (LFV) processes at 90% CL. Here µ → e refers to conversion rate in gold (Au) and titanium (Ti).

Lepton flavour violation
We use the usual expressions for the lepton flavour violation (LFV) processes (including µ − e conversion rates in various elements) that are applicable for the ScM [11,22]. For instance, radiative decays are given by where The various LFV processes that we include in our study are summarised in table 2.

Relic abundance
The ScSM permits both scalar and fermionic DM candidates, but we focus on the case of scalar DM for two reasons: 1. The fermion DM case is very similar to the usual ScM where strong constraints from LFV exist on the Yukawa couplings. It requires special textures, and/or coannihilations and/or fermion triplets instead of singlets to saturate the observed DM abundance. In our scan, we indeed find the need for coannihilations (see section 5).
2. It is interesting to study the rich phenomenology of CP-even scalar DM (η 2 ) as an admixture of singlet (pure singlet case is θ 0) and doublet (pure doublet case is θ π/2) components, see eq. (2.7). The doublet case (both CP-even and CP-odd) has the phenomenology of the ScM, while the singlet has some extra terms from the scalar potential with respect to just a pure singlet scalar. As the CP-odd scalar A can be a DM candidate, there are a few possible mass hierarchies: In general, several (co-)annihilation channels are possible in the ScSM. All of these are included in micrOMEGAs v5.2.0 [54] which we use to compute the DM relic abundance. The DM abundance is required to be equal to or smaller than the Planck (2018) measured abundance [55]: Thus, the Planck measurement provides an upper limit on the DM relic abundance.

Direct detection
Direct detection (DD) typically imposes strong constraints on scalar DM candidates with a non-zero hypercharge due to the presence of t-channel Z/h-mediated diagrams [8]. The gauge interactions stem from the kinetic term for the doublet, namely Notice that the Z-mediated interactions are inelastic, and always involve a CP-even scalar and pseudoscalar A. A small enough mixing can suppress DD limits via the last term in eq. (3.32), as in this case, η 2 is mainly singlet and does not directly couple to the Zboson. In addition, a large enough mass splitting ( MeV) can kinematically forbid the scatterings between η 2 and A. This implies that λ HΦ,3 10 −6 . Two cases are possible: 1. For b a (and b c), η 2 is mainly doublet (with a correction proportional to κ) and the mass splitting with A is given by λ HΦ,3 (which can be made naturally smaller than in the ScM); 2. For b a (and a c), η 2 is mainly singlet and m A is given in terms of other parameters, so there is no reason for the mass splitting to be small, and thus inelastic scatterings are expected to be forbidden.
Interactions mediated by the SM Higgs h lead to the usual elastic (and also inelastic) spin-independent (SI) scattering; the resulting limits are also quite severe [56,57], but they can be suppressed by small scalar couplings unlike in the case of Z-mediated interaction [58]. For non-zero mixing, the dimensionless h-η 2 -η 2 coupling is given by [28] where λ 123 ≡ λ HΦ,1 + λ HΦ,2 + λ HΦ,3 . As λ eff decreases, the h-mediated direct detection cross section also decreases. This coupling will vanish (i.e., no overall h-η 2 -η 2 coupling) when the mixing angle satisfies the following relation: We confirm that when the coupling λ Hϕ (λ 123 ) vanishes, the Higgs-DM coupling is zero only for pure singlet (doublet) DM. In the event that the couplings are exactly equal and dominate over the bare masses, λ 123 = λ Hϕ |m 2 Φ − m 2 ϕ |/v 2 , the effective Higgs coupling λ eff vanishes for maximal mixing θ = π/4.
Using micrOMEGAs, we compute the effective SI DM-proton scattering cross section, assuming that the local DM energy density scales proportional to the global one: where Ω X is the X DM abundance (X = η 2 , A or ψ 2 ) and Ω DM is the Planck measured abundance in eq. (3.31). We then recast the observed exclusion limit from XENON1T [59] within micrOMEGAs [60] for a fixed DM mass and compare it against the effective SI cross section in eq. (3.35). There are a number of planned xenon-based experiments that will increase the sensitivity significantly, e.g., XENONnT [61], PandaX [62], LZ [63] and DARWIN [64]. In our plots, we will illustrate the expected reach of LZ only.

Z 2 symmetry at high energies
The ScSM also has interesting features in light of the model viability up to high-energy scales [21,[65][66][67]. The trilinear coupling κ, due to the presence of the scalar singlet, gives positive contributions to the RGE evolution of the new scalar mass-squared parameters if it is real, thereby helping prevent the breaking of the model symmetries even in the absence of finite temperature effects. The compatibility with the evolution of the bare Higgs mass (and thus of EWSB) requires a separate study, which we leave for a future work. In principle, one could incorporate the evolution of Renormalisation Group Equations (RGEs) of model parameters into our numerical scan to verify that the Z 2 symmetry remains unbroken at high-energy scales, but the practical implementation remains computationally difficult. Thus, we do not include RGE effects in our numerical analysis.

Numerical analysis
To efficiently sample the allowed parameter space of the ScSM, we use the Importance Nested Sampling algorithm implemented in MultiNest v3.10.0 [68] with 25,000 live points (nlive) and a stopping tolerance (tol) of 10 −3 . 8 The composite log-likelihood used is ln L total (θ) = ln L κ (θ) + ln L EWPT (θ) + ln L Rγγ (θ) where θ are the free parameters of the ScSM. The individual log-likelihood contributions are described below: 1. ln L κ (θ): log-likelihood for the trilinear coupling κ. It is Gaussian in nature, centered at 0 TeV with a standard deviation of 1.5 TeV, see eq. (3.7).
2. ln L EWPT (θ): log-likelihood for the electroweak precision tests (EWPT) (see subsection 3.4). It is given by [69] ln where ∆O i are the central values for the shifts in eq. (3.25), Σ 2 ij ≡ σ i ρ ij σ j is the covariance matrix, ρ ij is the correlation matrix in eq. (3.26) and σ i are the associated errors in eq. (3.25).

ln L
It is a Gaussian likelihood function, centered at the PDG measured value of 1.1 with a standard deviation of 0.1 [53]. The SM expectation is R γγ = 1.0.
4. ln L LFV (θ): log-likelihood for the LFV processes (see subsection 3.6). It is given by Each of the individual likelihood functions are Gaussian, and centered at a branching ratio/conversion rate of 0 with a standard deviation equal to the respective upper limit shown in column 2 of table 2.
5. ln L Ωh 2 (θ): log-likelihood for the DM relic density Ωh 2 (see subsection 3.7). It is a onesided Gaussian, i.e., a flat likelihood on Ω X ≤ Ω DM and Gaussian for Ω X > Ω DM . The Planck measured uncertainty is also combined in quadrature with a 5% theoretical uncertainty (stemming from our assumed uncertainty on the relic density calculation in micrOMEGAs).
6. ln L DD (θ): log-likelihood for the XENON1T experiment (see subsection 3.8). It is a simple step-function-like likelihood, i.e., parameter points are allowed (rejected) if the effective SI cross section in eq. (3.35) is below (above) the official XENON1T exclusion limit [59] for a given DM mass.
The ranges and priors for the 21 (12 free + 9 nuisance) model parameters in normal ordering (NO) and inverted ordering (NO) are summarised in table 3. 9 Due to the presence of coannihilations, we find it efficient to scan over δ 1 and δ A (instead of m η 1 and m A ) where (4.4) 9 We keep the CP-even (mη 1,2 ) and CP-odd (mA) scalar masses 100 GeV to avoid constraints from h/Z invisible decays and collider limits.   By convention, δ 1 ≥ 0 as m η 1 ≥ m η 2 . Meanwhile, δ A > 0 for η 2 as DM, δ A < 0 for A as DM, and either for ψ 2 as DM.
In the next section, we show various two-dimensional (2D) plots of the profile likelihood ratio (PLR) [70] in the relevant parameter planes or key observables of interest. Model parameters that are not shown in those plots are profiled over, i.e., the composite loglikelihood function in eq. (4.1) is maximised with respect to those parameters. Using Wilks' theorem [71], the PLR can be used as a test statistic to approximately construct the 1σ (∼ 68.3%) and 2σ (∼ 95.4%) CL contours [72].

Results
We start by showing results for the scalar DM in the case of no mixing between the singlet and doublet (e.g., θ = 0, π/2). For scalar doublet DM, we reproduce the standard results [73,74]; we consider the CP-even scalar, but the results for the CP-odd scalar are similar. For singlet DM, we also reproduce the results from the literature [56,57,75,76]. 10 Next, we turn on the mixing angle θ and consider separately the case of real scalar η 2 , real pseudoscalar A and Majorana fermion ψ 2 as DM candidates. We pay special attention on studying how the parameter space opens up in each cases with respect to the usual ScM.

No mixing case: Scotogenic model + scalar singlet
From eq. (2.7), the physical state η 2 in the case of no mixing is In both cases, the trilinear scalar coupling κ = 0, and the model reduces to the usual Scotogenic model plus a scalar singlet; we indeed recover the results from the literature.
In the left (right) panel of figure 2, we plot the η 2 relic density (effective SI η 2 scattering cross section with protons) versus the singlet DM η 2 mass for θ = 0. We see two allowed disconnected regions, at masses around 150 GeV and above ∼ TeV. Notice that only in the latter region, the singlet can constitute 100% of the observed DM abundance [56]. As we see from the plot in the right panel, next-generation DD experiments will be able to test this model, given the fact that the quartic couplings are large enough to reproduce the abundance.
Similarly in figure 3, we plot the η 2 relic density (effective SI η 2 -p cross section) versus the CP-even doublet DM η 2 mass for θ = π/2. 11 The upper left corner in the left plot implies a too-large annihilation cross-section to reproduce the observed DM abundance. Indeed, for the doublet DM to saturate the abundance, its mass should be 500 GeV. In this case, next-generation DD experiments will still leave a large portion of the parameter space unexplored. This is due to the fact that the annihilation cross section, driven by gauge interactions, is decoupled from the DD cross section (unless the mass splitting with the CP-odd scalar is MeV). 10 There is an extra parameter with respect to both models separately, λΦϕ, but as expected, we see that it does not affect the model phenomenology. 11 Notice the under-sampling of the profile likelihood surface at small cross sections. This is expected due to the nature of our XENON1T likelihood (a step-function-like), i.e., parameter points have same likelihood if they are all compatible with the official XENON1T limit. A full coverage of this region requires a dedicated scan over extremely small h-η2-η2 couplings, which is not the main goal of our study. similar structure is seen in the case of the trilinear coupling κ (top-right plot) for the same reasons, although in this case, the region is much less pronounced. In the bottomleft plot, we see how the Yukawa couplings are correlated amongst each other, with the heaviest sterile (y 1 ) being larger than the lightest one (y 2 ); this correlation becomes even more pronounced at large couplings. In the bottom-right plot, we observe how the ScSM demands a relationship between the trilinear and quartic couplings, such that neutrino masses are reproduced. This can be understood analytically in the limit of heavy scalar singlet masses, see eq. (3.16), as represented by a solid brown line in the plot.
In the top-left plot in figure 5, the η 2 abundance is much smaller for m η 2 ∼ m h due to direct annihilation process η 2 η 2 → h h. Final states with gauge bosons (W + W − /ZZ) are always open. Indeed, the dominant annihilation channels that determine the η 2 abundance involve gauge bosons, and less often the Higgs bosons. Annihilations into tt, leptons or photons are sometimes present. When the scalar masses are degenerate enough, coannihilations can be important. An upper limit of m η 2 5 TeV is obtained at the 1σ CL, but is somewhat below the upper limit from the prior of 10 TeV. In the top-right plot, we show the effective SI η 2 -proton scattering cross section. The rise in upper limit with DM mass is expected from the DM number density for heavier masses. The small cross sections arise from a cancellation in eq. (3.34); in any case, they are well below the sensitivity of next-generation DD experiments.
In the bottom-left plot in figure 5, we observe that the Higgs to di-photon rate with respect to the SM is enhanced (suppressed) for λ HΦ,1 < 0 (λ HΦ,1 > 0). It is evident that the currently allowed value, shown by horizontal brown lines, demands λ HΦ,1 < 0. Finally, in the bottom-right plot, we show the effective neutrino mass parameter m ee ≡ U 2 ei m i , which enters in the expression for the lifetime of neutrinoless double beta (0νββ) decay [78]. As expected, we reproduce the expected result for NO and IO. In particular, IO results imply that 0.012 eV ≤ m ee ≤ 0.05 eV. These values can potentially be tested in coming years, see refs. [79,80] for recent reviews.
In the left panel of figure 6, we plot an LFV radiative decay versus the Yukawa of the heaviest fermion singlet; here we observe a V -shaped region. The behaviour for large Yukawas goes as ∼ |y| 4 , as expected. For |y 1 | 10 −4 , the contribution from other neutrino (∝ |y 2 |) dominates. The upper-left region corresponds to |y 2 | |y 1 |; we see that it is not allowed, as the mass of lightest fermion singlet is too light to suppress enough LFV. In the right panel of figure 6, we see how the next-generation DD experiments (e.g., LZ projected sensitivity for 50 GeV DM [77] -dashed orange line) test complementary parts of the parameter space to those of LFV experiments (e.g., expected sensitivity of µ − e conversion rate (an improvement by 4 orders of magnitude) -solid brown line).
We have checked that different LFV observables (see table 2) show a clear correlation among themselves. This is expected from the fact that Yukawa couplings are smaller than one, so that the box (dipole) contributions are suppressed (dominant): It is interesting to highlight that the sensitivity to BR(µ → 3e) is expected to improve by up to 4 orders of magnitude [81]. These conclusions have already been obtained in the literature, although for somewhat different versions of the model [11,13,22]. In addition, we have checked that there are no significant differences between the two mass orderings (NO vs IO).
In the left panel of figure 7, we show how the η 2 relic abundance changes with the dimensionless h-η 2 -η 2 coupling λ eff , see eq. (3.33). We observe how the smallest relic abundance is obtained for λ eff values of order one, in which the annihilations proceed via a Higgs-mediated s-channel diagram. In the right panel, we see that the SI DD cross section scales linearly with λ eff , meaning that Z-mediated processes do not contribute significantly to the DD cross section.

DM candidate: real pseudoscalar, A
The allowed parameter space and phenomenology of the ScSM in this case is similar to the real (CP-even) scalar coming from the doublet (see figure 3). The relic abundance of A is also bounded from above for low masses, as shown in the left panel of figure 8. This translates into a lower limit of m A 300 GeV for pseudoscalar A to saturate the observed DM abundance, somewhat smaller than for the CP-even candidate. This difference comes from the coupling λ HΦ,3 , which enters differently into the scalar masses.

DM candidate: Majorana fermion, ψ 2
We do not investigate the case of fermion DM in detail, as it has already been explored in a model with similar phenomenology [22]. Here we simply confirm that we arrive at the same conclusions, namely that the fermion DM case requires coannihilations for its abundance to match the Planck measured value. In figure 9, we see that viable parameter space requires m ψ 2 to be degenerate with m η 2 (left panel ) and m φ + (right panel ). This is so because the scalar masses are set by similar combinations of Lagrangian parameters, and EWPT demands them to be close in mass, although there is somewhat a wider region at small masses for the charged scalar. Thus, the case of fermion DM in the ScSM introduces a fine-tuning that is not present in the case of scalar DM. Apart from this difference, the Yukawa couplings and LFV processes are similar to the case of η 2 and A DM. In regards to DD in the Scotogenic model with a fermion singlet, it occurs at one-loop level, and is typically suppressed [9,22,82,83].

Conclusions
We have proposed a simple variation of the original Scotogenic Model (ScM), namely with an extra real scalar singlet. The model, termed the ScotoSinglet Model (ScSM), is arguably the simplest extension of the popular ScM, with a very rich phenomenology and several interesting features: 1. It allows for DM to be scalar (CP even or odd) with a naturally-suppressed direct detection rate, either due to a typically large mass splitting with the opposite-CP scalar (not only dependent on λ HΦ,3 as in ScM), or due to a small mixing, θ 0. In this case, DM is mainly singlet with a small doublet component, and Z-boson mediated interactions are suppressed; 2. There are two contributions to neutrino masses: the usual Scotogenic one (∝ λ HΦ, 3 ) and the new singlet one (∝ κ 2 ). In principle, lepton number violation (e.g., the smallness of neutrino masses) demands both couplings to be small, which is technically natural. However, in the limit of large singlet mass, a large λ HΦ,3 can be cancelled with the trilinear coupling term (∝ κ 2 /m 2 ϕ ), see eq. (3.16). This allows for the cou-plings to be larger than usual and contribute to the DM phenomenology, e.g., to DM annihilations and scatterings.
3. The presence of the singlet improves the stability of the Z 2 symmetry up to highenergy scales when the trilinear coupling κ is real, as it contributes positively to the evolution of Renormalisation Group Equations (RGEs) for m 2 ϕ and m 2 Φ . The above features significantly open up the parameter space with respect to the ScM. Extensions to three sterile fermions are not expected to change our results significantly. Other variations, such a proper RGE study, or a Generalised ScotoSinglet Model [22], is left out for a future work.
The origin of neutrino masses, the nature of DM, and their possible connection remains an open question that may possibly take several decades to fully understand. While we wait eagerly for a positive signal, the study of simplified models (such as the ScSM) allows us to gain insight into the big puzzles, and search for new correlations among different observables that can help us in distinguishing models among a plethora of possibilities.

A Mass eigenstate basis
For the weak eigenstates A = (φ R , ϕ) T , the mass-term is given by is a (non-diagonal) squared mass matrix. To diagonalise M 2 , we define the mass eigenstates (η 1 , η 2 ) as   η 1 where θ is the mixing angle. Thus,  In the mass eigenstate basis, a squared mass matrix satisfies the following relation: Following the analysis in appendix B of ref. [85], we find that where the last equality can also be expressed as In matrix notation, the above expressions read as  By taking θ → −θ, we can express (a, b, c) in terms of (m 2 η 1 , m 2 η 2 , θ) as Using the relations for (a, b, c) from section 2, we get

B Parameterisation of the Yukawa couplings
Following the working in the original Casas-Ibarra paper [38] (see also ref. [86] for a oneloop parametrization, and ref. [87] for a general parametrization), we write the neutrino mass matrix in terms of our high energy parameters, equivalent to eq.   [88,89].
It is easy to check that this satisfies the low energy definition by substituting back into the high energy expression. We parametrise the R matrix as (B.10) where c ij ≡ cos θ ij and s ij ≡ sin θ ij (θ 12 , θ 13 , and θ 23 are the 3 lepton mixing angles), α (δ CP ) is the Majorana (Dirac) phase. As the lightest neutrino is massless with just two fermionic singlets, there is only one physical Majorana phase. For the neutrino oscillation parameters, we use the results based on a global fit from the Nu-FIT collaboration [88,89] with SK atmospheric data (see table 4).

D Renormalisation Group Equations (RGEs)
Here we provide the Renormalisation Group Equations (RGEs) for the ScSM at one-loop level, as computed using the SARAH package [92].