Generalizing the Scotogenic model

The Scotogenic model is an economical setup that induces Majorana neutrino masses at the 1-loop level and includes a dark matter candidate. We discuss a generalization of the original Scotogenic model with arbitrary numbers of generations of singlet fermion and inert doublet scalar fields. First, the full form of the light neutrino mass matrix is presented, with some comments on its derivation and with special attention to some particular cases. The behavior of the theory at high energies is explored by solving the Renormalization Group Equations.


Introduction
The experimental observation of neutrino flavor oscillations constitutes a milestone in particle physics and proves that the Standard Model (SM) is an incomplete theory. Although many questions remain open, such as the Majorana or Dirac nature of neutrinos or the possible violation of CP in the leptonic sector, the SM must certainly be extended to include a mechanism that accounts for non-zero neutrino masses and mixings.
Many neutrino mass models have been proposed along the years. Among them, radiative models are particularly appealing. After the pioneer models in the 80's [1][2][3][4], countless radiative models have been proposed and studied [5]. The suppression introduced by the loop factors allows one to accommodate the observed solar and atmospheric mass scales with sizable couplings and relatively light (TeV scale) mediators. This typically leads to a richer phenomenology compared to the usual tree-level scenarios and, in fact, the new mediators may even be accessible to current colliders. Furthermore, in some radiative models one can easily address a completely independent problem: the nature of the dark matter (DM) of the Universe. Discrete symmetries, connected to the radiative origin of neutrino masses, may be used to stabilize viable DM candidates, resulting in very economical scenarios [6].

The general Scotogenic model
The Scotogenic model [7] is a simple extension of the SM that induces radiative neutrino masses and provides a potential dark matter candidate. Here we consider a generalization of the model. The SM particle content is extended by an unspecified number, n N , of singlet fermions N , and also an arbitrary number, n η , of inert scalar doublets η. Particular cases of this particle spectrum can be labeled by their (n N , n η ) values. In addition, the symmetry group of the SM is enlarged with a dark Z 2 parity, under which all the new fields are odd, while the SM particles are even. The scalar and fermion particle content of the model, as well as their representations under the gauge group SU(3) c × SU(2) L × U(1) Y and the Z 2 parity of the model are given in Tab. 1.
The relevant Yukawa and bare mass terms for our discussion are where n = 1, . . . , n N , a = 1, . . . , n η and α = 1, 2, 3 are generation indices and y is a general complex n N × n η × 3 object. Besides, M N is a symmetric n N × n N Majorana mass matrix that has been chosen diagonal without loss of generality. Furthermore, one can also write Here all the indices are η generation indices. Therefore, m 2 η and λ 3,4,5 are n η × n η matrices while λ 2 is an n η × n η × n η × n η object. Note that λ 5 must be symmetric whereas λ 3,4 must be Hermitian. Again, m 2 η will be assumed to be diagonal without loss of generality. Finally, we highlight the presence of the scalar potential quartic couplings λ ab 5 , which play a major role in the neutrino mass generation mechanism, as shown in Sec. 3.
We will assume that the minimization of the scalar potential in Eq. (2) leads to the vacuum configuration with a = 1, . . . , n η . Therefore, only the neutral component of H acquires a non-zero vacuum expectation value (VEV), which breaks the electroweak symmetry in the standard way, while the η a scalars are inert doublets with vanishing VEVs. In this way, the Z 2 symmetry remains unbroken and the stability of the lightest Z 2 -charged particle is guaranteed. We will come back to the possibility of Z 2 breaking due to Renormalization Group Equations (RGEs) effects later. We now decompose the neutral component of the η a multiplets, η 0 a , as In the following we will assume that all the parameters in the scalar potential are real, hence conserving CP in the scalar sector. In this case, the real and imaginary components of η 0 a do not mix. After electroweak symmetry breaking, the n η × n η mass matrices for the real and imaginary components are given by and respectively. We note that M 2 R = M 2 I in the limit λ 5 → 0, in which all the elements of λ 5 vanish. This will be crucial in the calculation of neutrino masses, as shown below. Both mass matrices can be brought into diagonal form by means of a change of basis. The gauge eigenstates, η Aa , are related to the mass eigenstates,η A b , where A = R, I, by Here η A andη A are n η -component vectors. In general, the n η × n η matrices V A are unitary, where I nη is the n η × n η identity matrix. However, in the simplified scenario of CP conservation in the scalar sector, M 2 R and M 2 I are real symmetric matrices, and then the V A matrices are orthogonal, such that With these transformations, the diagonal mass matrices are given by The resulting analytical expressions for the mass eigenvalues m 2 Aa and mixing matrices V A involve complicated combinations of the scalar potencial parameters. However, under the assumptions 2 λ aa one can find simple expressions. The m 2 Aa mass eigenvalues are given by We note that the mass splitting m 2 Ra − m 2 Ia = λ aa 5 v 2 vanishes in the limit λ 5 → 0. In what concerns the V A orthogonal matrices, each of them can be expressed as a product of n η (n η − 1)/2 rotation matrices, with the scalar mixing angles given by where the κ 2 A sign (κ 2 R = +1 and κ 2 I = −1) has been introduced. The generation of neutrino masses takes place at the 1-loop levelà la scotogenic [7]. In the presence of the terms given in Eqs. (1) and (2), lepton number is explicitly broken in two units, hence inducing Majorana neutrino masses. Assuming that the potential is such that the η a scalars do not get VEVs, see Eq.

Neutrino masses
where m A ν an αβ is the contribution to (m ν ) αβ generated by the N n − η Aa loop, given by where D = 4 − ε is the number of space-time dimensions, the external neutrinos are taken at rest and k is the momentum running in the loop. We note that the term proportional to / k does not contribute because it is odd in the loop momentum. C A naα is the N n − η Aa − ν α L coupling, given by with κ R = 1 and κ I = i. Since we assume real parameters in the scalar sector, complex conjugation in V A will be dropped in the following. Replacing Eq. (15) into Eq. (14) and introducing the standard Passarino-Veltman loop function B 0 [80], where ∆ ε diverges in the limit ε → 0, Eq. (13) becomes Eq. (17) constitutes our central result for the 1-loop neutrino mass matrix in the model. It is important to note that the divergent pieces cancel exactly. Indeed, the κ 2 A factor implies that the term proportional to ∆ ε in Eq. (17) involves the combination which vanishes due to the orthogonality of the V A matrices, ensuring the cancellation of the divergent part of the B 0 functions. This was expected since the neutrino mass matrix is physical and therefore finite. While Eq. (17) provides a simple analytical expression for the neutrino mass matrix, the dependence on the fundamental parameters of the model is not explicit. The neutrino mass matrix involves a product of V A matrices and B 0 functions, both in general depending on the scalar potential parameters in a non-trivial way. In order to identify more clearly the role of the scalar potential parameters, we will work under the assumptions in Eq. (9) and derive an approximate form for the neutrino mass matrix, valid for small λ ab 5 couplings and small mixing angles in the scalar sector. First, it is convenient to make an expansion in powers of λ ab where the superindex (i) , with i = 0, 1, denotes the order in λ ab 5 . We highlight that the expansion begins at 1st order in λ 5 . This was indeed expected, since λ 5 = 0 would imply the restoration of lepton number and massless neutrinos. With this in mind, the origin of the two terms in Eq. (19) is easy to understand. In the first term, the λ ab 5 couplings are neglected in the V A matrices but kept at leading order in the B 0 functions. This term is proportional to the B 0 (0, m 2 Ra , M Nn ) − B 0 (0, m 2 Ia , M Nn ) difference, which would vanish for λ aa 5 = 0, see Eqs. (10) and (11). The mass matrices for the real and imaginary components of η 0 are equal at 0th order in λ 5 , M In the second term, the λ ab 5 couplings are neglected in the B 0 functions but kept at leading order in the V A mixing matrices. Since m We note that this term will only be non-zero when the λ 5 matrix contains non-vanishing off-diagonal entries, since this is the only way the ( vanish at 1st order in λ 5 . Next, we find approximate expressions for the V A mixing matrices. This is only feasible by assuming small scalar mixing angles, in agreement with Eq. (9). In this case one can expand V not only in powers of λ 5 , but also in powers of the small parameter which is defined for a = b and corresponds to sin θ ab R or sin θ ab I at 0th order in λ 5 , see Eq. (12). With this definition, one finds the general expression (V ) ab = δ ab + (1 − δ ab ) s ab + O (s 2 ). Analogous expressions are found for V R and V I replacing s by sin θ R and sin θ I , respectively. With all these ingredients, Eq. (19) can be written as where we have defined the dimensionless quantity and the loop functions Eq. (22) involves the quantity Γ abn , which we have written in Eq. (23) as the sum of two terms. The first term in Γ abn contributes only for a = b and involves only diagonal elements of λ 5 . The second term, which involves diagonal as well as off-diagonal elements of λ 5 , only contributes for a = b. We also note that g abn = −g ban . Eq. (22) is the main analytical result of our work. Under the assumptions of Eq. (9), it reproduces the neutrino mass matrix in very good approximation. It is valid for any n N and n η values. We will now show how in some particular cases it reduces to well-known expressions in the literature.
3.1 Particular case 1: (n N , n η ) = (3, 1) The first example we consider is the standard Scotogenic model originally introduced in [7] and obtained for (n N , n η ) = (3, 1). In this case, only one inert doublet η is introduced. Therefore all the matrices in the scalar sector become just scalar parameters: Besides, the Yukawa couplings become 3 × 3 matrices: y naα ≡ y n1α ≡ y nα . Similarly, f an ≡ f 1n ≡ f n , and the second term in Eq. (23) does not contribute. With these simplifications, the general Γ abn reduces to Γ (3,1) n , given by Replacing this into Eq. (22), one obtains the well-known neutrino mass matrix with This expression agrees with [7] up to a factor of 1/2 that was missing in the original reference. 3 3.2 Particular case 2: (n N , n η ) = (1, 2) A version of the Scotogenic model with one singlet fermion and two inert doublets, (n N , n η ) = (1, 2), has been considered in [77,78]. Since the model contains only one singlet fermion N , M Nn ≡ M N is just a parameter. The Yukawa couplings become 2 × 3 matrices: y naα ≡ y 1aα ≡ y aα . Finally, f na ≡ f 1a ≡ f a and g abn ≡ g ab1 ≡ g ab . Both references work in the basis in which the m 2 η matrix is diagonal. However, they take different simplifying assumptions about the scalar potential parameters.
In [77] the matrix λ 3 + λ 4 was assumed to be diagonal. In this case, which we denote as scenario (1, 2) I , (1 − δ ab )s ab = 0 and the general Γ abn reduces to Replacing this expression into Eq. (22) and arranging the different pieces properly, one obtains which agrees with the result in [77] up to a global factor of 1/4. On the other hand, a diagonal λ 5 matrix was taken in [78]. We denote this as scenario (1, 2) II . Again, this simplifies Γ abn , which becomes With this result, one can easily use Eq. (22) to derive with which agrees with the expression given in [78] if terms of order s 2 12 are neglected.

High-energy behavior
The conservation of the Z 2 parity is crucial for the Scotogenic setup to be consistent. In the absence of this symmetry, neutrinos would acquire masses at tree-level and the DM candidate would no longer be stable. This motivates the study of the conservation of Z 2 at high energies, a line of work initiated in [81]. As pointed out in this reference, the RGE flow in the Scotogenic model might alter the shape of the scalar potential at high energies and lead to the breaking of Z 2 . This issue was fully explored in subsequent works [83,84], which show that the breaking of the Z 2 parity actually takes place in large regions of the parameter space. A similar discussion for a variation of the Scotogenic model including scalar and fermion triplets was presented in [46]. Some general features of the high-energy behavior of the model, and in particular of the possible breaking of the Z 2 symmetry, can be understood by inspecting the 1-loop β function for the m 2 η parameter, shown in Appendix A. Eq. (39) generalizes the result previously derived in [81] and gives the 1-loop β function for the m 2 η matrix, valid for any values of (n N , n η ). In order to study the possible breaking of Z 2 , one must consider the sign (positive or negative) of the individual contributions to the running of m 2 η . In this regard, the negative contribution of the term proportional to Tr y † a M * N M N y b turns out to be crucial. In the following, we will refer to this term as the trace term. As first pointed out in [81] for the standard Scotogenic model, in case of large Yukawa couplings (equivalent to λ 5 1) and M 2 N m 2 η , the trace term dominates the m 2 η running and drives it towards negative values. Eventually, this leads to the breaking of the Z 2 symmetry at high energies, once m 2 η < 0 induces a minimum of the scalar potential with η = 0. The same behavior is expected in the general Scotogenic model. Other terms in Eq. (39) may counteract this effect. In particular, the terms proportional to the quartic scalar couplings may do so if their signs are properly chosen. The contribution to the m 2 η running will be positive for λ 2 > 0 and λ 3,4 < 0 (since m 2 H < 0), while their effect will reinforce that of the trace term otherwise. We will now explore the scalar potential of the model at high energies by solving the full set of RGEs numerically. In order to do that we will concentrate on two specific (but representative) versions of the general Scotogenic model: • The (3, 1) model, with three singlet fermions and one inert doublet. This is the original Scotogenic model [7].
We set all model parameters at the electroweak scale, which we take to be the Z-boson mass, m Z . Therefore, in the following all values for the input parameters must be understood to hold at µ = m Z . We compute m 2 H by solving the tadpole equations of the model and set the λ 1 value to reproduce the measured Higgs boson mass. The remaining scalar potential parameters are chosen freely, but always to values that guarantee that the potential is bounded from below (BFB) at the electroweak scale. This is a non-trivial requirement due to the complexity of the scalar potential of the general Scotogenic model. We refer to Appendix B for a detailed discussion on how we check boundedness from below. Finally, we must accommodate the neutrino squared mass differences and the leptonic mixing angles measured in neutrino oscillation experiments by properly fixing the Yukawa couplings of the model. In the two variants of the general Scotogenic model considered the Yukawa couplings become 3 × 3 matrices, and then they can be obtained by means of a Casas-Ibarra parametrization [85], adapted to the Scotogenic model as explained in [86][87][88][89]. This allows us to write the Yukawa matrices in full generality as Here U is a 3 × 3 unitary matrix, defined by the Takagi decomposition of the neutrino mass matrix with m i the three physical neutrino masses. R is a general 3 × 3 orthogonal matrix and we have defined  33) ensures compatibility with neutrino oscillation data. We consider neutrino normal mass ordering and the 1 σ ranges for the oscillation parameters obtained in the global fit [90], including the CP-violating phase δ, hence allowing for complex Yukawa couplings. For simplicity, we take m 1 = 0 and R = I, with I the 3 × 3 identity matrix. 4 First, we have rediscovered the parity problem in the standard (3, 1) Scotogenic model. This is shown on the left-hand side of Fig. 2, which displays the RGE evolution of the CPeven scalar mass m R with the energy scale µ. This is the most convenient parameter to study the breaking of the Z 2 symmetry. When m 2 R becomes negative, the lightest CP-even scalar becomes tachyonic, a clear sign that η = 0 is not the minimum of the potential. We have checked that the scalar potential is BFB at all energy scales in this figure. We note that due to our parameter choices the lightest singlet fermion, N 1 , has vanishing Yukawa couplings. For the same reason, y 2α y 3α and the effect is driven predominantly by N 3 . This explains the drastic change in the evolution of m R at µ = 2 TeV, when N 3 becomes active. Below this scale, N 3 effectively decouples and does not contribute to the RGE running. We point out that a much less pronounced change takes also place at µ = 1.5 TeV, when N 2 becomes active, but this is not visible on the figure. The Z 2 parity gets broken at µ 60 TeV, after which the η R state becomes tachyonic. These results agree well with those found in [81] and confirm the possible breaking of Z 2 in the original Scotogenic model. A very similar behavior is found for the (1, 3) model, which only has one singlet fermion, as shown on the right-hand side of Fig. 2. In this case, the three CP-even scalar masses m Ra are displayed. Again, we have checked that the scalar potential is BFB at all energy scales in this figure. As in the case of the standard Scotogenic model, when one of the CP-even scalar masses reaches zero the Z 2 symmetry gets broken. We see in this figure that this happens at µ 15 TeV, where   In the following we concentrate on the (1, 3) model. As already discussed, the singlet fermion mass M N drives the scalar masses towards negative values via the trace term, hence breaking the Z 2 parity at high energies. Fig. 3 shows the Z 2 breaking scale as a function of M N for several scalar parameter sets. The blue and red lines correspond to moderate values for the quartic couplings, λ aaaa 2 = λ aa 3 = λ aa 4 = 0.1, while the green line has increased (and additional) quartics, λ aaaa 2 = λ aabb 2 = λ aa 3 = λ aa 4 = 0.3. The λ 5 matrix is taken to be diagonal, with λ aa 5 = 10 −9 . We have explicitly checked that the scalar potential is BFB at the electroweak scale in all scenarios. 5 As expected, the Z 2 breaking scale decreases for larger M N since the effect of the trace term becomes stronger. While different scalar potential couplings may alter the outcome, this generic behavior is found in large portions of the parameter space. One should notice, however, that the green curve begins at M N 2 TeV. For this specific scenario, lower values of M N do not break the Z 2 symmetry, as we now proceed to discuss. Fig. 4 shows the evolution of the lightest scalar mass m R 1 as a function of the energy for the parameter values corresponding to the green curve in Fig. 3. The results have been obtained for several values of M N . It is important to note that this figure shows the mass of the lightest scalar at each energy scale, and not the mass of a single mass eigenstate at 5 We have allowed for (possible) non-BFB potentials at high energies, where some of the quartic couplings become negative due to running effects. We note that our algorithm gives us only sufficient (and not necessary) boundedness from below conditions, and in principle some of the possibly non-BFB potentials might actually be BFB. Morevoer, even non-BFB potentials may be realistic if the electroweak vacuum is metastable and has a large enough lifetime. This issue is already present in the SM and is clearly beyond the scope of our analysis, which focuses on the possible breaking of the Z 2 symmetry.  (39). The high multiplicity of λ 2 couplings reinforces the effect. Actually, we note that this Landau pole is present at very high energies, well above the Z 2 breaking scale, for many choices of the parameters at the electroweak scale.
We conclude our exploration of the high-energy behavior of the (1, 3) model with Fig. 5. In this case we plot the Z 2 breaking scale as a function of one of the λ 2 parameters, namely λ 2233 2 . This is done for three scenarios: the blue curve corresponds to λ aaaa TeV. In all cases we have checked that the scalar potential is BFB at the electroweak scale. For the blue and green lines, the impact of λ 2233 2 is relatively mild. This is because the high values of M N (5 and 9 TeV, respectively) make the trace term completely dominant and break the Z 2 symmetry at relatively low energies. In contrast, the Z 2 breaking scale has a much stronger dependence on λ 2233 2 in the red scenario, which has a lower M N = 1.25 TeV. For λ 2233

Thermal effects and the fate of the Z 2 symmetry
To determine the cosmological impact of Z 2 breaking one needs to take into account thermal corrections. This is because the interaction with the hot, primordial plasma induces an effective potential for the scalar fields. This effective potential, at 1-loop order, is given by Here V CW is the standard Coleman-Weinberg potential for η at zero temperature while J B (m 2 N (η/T 2 )) and J F (m f (η) 2 /T 2 ) are the bosonic and fermionic functions, respectively. These functions admit a high-T expansion (see [91] for a review) which allows to write the scalar mass as The coefficient c depends on the details of the theory, such as the quartic, gauge and Yukawa couplings. 6 At any given time in the early Universe, as it can be seen in Eq. (36), the effect of the temperature is usually to restore the symmetry with the subsequent dilution of the effects of the running that we have discussed in the previous section. It is therefore mandatory to study if temperature has any impact on the fate of the Z 2 symmetry and, therefore, on the stability of DM. During inflation, the η field is expected to have large quantum fluctuations, comparable to the Hubble parameter in this period, H I . These fluctuations can be much larger than the scalar mass at zero temperature and, acting as a sort of random walk, might bring the field to a vacuum where the Z 2 is broken. Right after reheating, when the temperature is potentially very large, the thermal mass of the scalar field may be large enough to overcome all breaking effects. The reason is that, assuming the decay of the inflaton is fast enough (instantaneous reheating, Γ Φ ∼ H I ), the reheating temperature is roughly given by [94] T RH ∼ 10 −1 H I M P , where M P is the Planck mass. Note that this temperature is generically much larger than H I . If the number of e-folds is not exceedingly large, T RH is expected to be larger than any field excursion and we expect m η (T RH ) 2 > 0. In addition, this also implies that m η (T RH ) 2 ∼ c T 2 RH H(T RH ) 2 , meaning that the field will fastly roll down to the minimum, at zero value η = 0. 7 As the temperature decreases, it may happen that RGE effects make the Z 2 breaking to occur at some high-energy scale. However, the η field will be already at η = 0, meaning that it cannot experience such a breaking. As the temperature continues decreasing, we reach the freeze-out temperature. From this point on, any breaking of the dark parity would be a disaster for the DM stability. Note however, that since the η field is at its local minimum, η = 0, it cannot notice this high-energy RGE induced symmetry breaking as it will only feel the local properties of the vacuum around η = 0.
Of course, this does not mean that RGE effects are completely harmless for the Scotogenic model. In fact, the RGE-induced breaking could induce the appearance of deeper minima in the potential, implying that the stability of DM is just a local property of our vacuum, which could be a false vacuum, and not a global feature of the potential.

Summary and discussion
The Scotogenic model is a well-known radiative scenario for the generation of neutrino masses. The introduction of three singlet fermions and one inert scalar doublet, all charged under a new Z 2 parity, leads to 1-loop Majorana neutrino masses and, as a bonus, provides a viable weakly-interacting dark matter candidate. In this work we have considered a generalization of this setup to any numbers of generations of singlet fermions and inert doublets. After computing the 1-loop neutrino mass matrix in the general version of the model, we have studied its high-energy behavior, focusing on two specific variants: the original Scotogenic model with (n N , n η ) = (3, 1) and a new multi-scalar variant with (n N , n η ) = (1, 3). Our main conclusion is that all the features of the original model are kept in the multi-scalar version, with some particularities due to the presence of a more involved scalar sector.
Our generalization of the Scotogenic model offer several novel possibilities. For instance, flavor model building could benefit from an interesting feature of multi-scalar versions of the model. In the (n N , n η ) = (1, 3) model, one obtains three massive neutrinos and leptonic mixing can be fully explained even if the Yukawa matrices are completely diagonal. In this case the leptonic mixing matrix would be generated by mixing in the scalar sector. This could be relevant in some flavor models. For example, it may be a crucial ingredient to rescue models where lepton mixing is predicted to be similar to quark mixing. Novel phenomenological signatures might exist as well. The η doublets can be produced at the Large Hadron Collider due to their couplings to the SM gauge bosons. Exotic signatures might be possible in models with many η generations, such as the (n N , n η ) = (1, 3) model. Cascade decays initiated by the production of the heaviest η doublets would lead to striking multilepton signatures, including missing energy due to the production of the lightest Z 2 -odd state at the end of the decay chain. Finally, the dark matter production rates in the early Universe might be affected as well by the presence of additional scalar degrees of freedom. These interesting questions certainly deserve further study.

A Renormalization Group Equations
At the 1-loop order, the RGEs of a model can be written as where t ≡ log µ, µ is the renormalization scale and β x is the 1-loop β function for the parameter x. In our analysis, the full 1-loop running in the Scotogenic model with arbitrary numbers of N and η generations has been considered. Analytical expressions for all the 1-loop β functions have been derived with the help of SARAH [95][96][97][98][99]. 8 These have been included in a code that solves the complete set of RGEs numerically.
We are mainly interested in the possible breaking of the Z 2 parity at high energies, and this is associated to the running of the m 2 η matrix. The corresponding 1-loop β functions are given by Here y a ≡ [y naα ] is a n N × 3 matrix, being the first index a singlet fermion family index and the third one a charged lepton family index. We have explicitly checked that for n N = 3 and n η = 1, Eq. (39) reduces to the m 2 η β function in the standard Scotogenic model [81].

B Boundedness from below
In order to ensure the existence of a stable vacuum, the scalar potential of the theory must be BFB. There exist several approaches to analyze boundedness from below. Ideally, one would like to have a BFB test that provides necessary and sufficient conditions. This way, one could not only guarantee that all potentials that pass the test are BFB (sufficient condition), but also discard potentials that fail it (necessary condition). In this regard, a major step forward was recently given in [100]. The algorithm proposed in this reference provides necessary and sufficient conditions for boundedness from below in a generic scalar potential using notions of spectral theory of tensors. However, applying this algorithm beyond a few simple cases turns out to be impractical due to the computational cost involved. For this reason, in phenomenological analyses one usually resorts to less ambitious approaches which only provide sufficient conditions, but not necessary. These methods are overconstraining, since one must reject potentials not passing the test, even though they might actually be BFB. Nevertheless, if the potential passes the test, one can fully trust that boundedness from below is guaranteed.
Here we will employ the copositivity criterion, which combined with a recently developed mathematical algorithm, never applied to a high-energy physics scenario, will give us sufficient (but not necessary) conditions. To the best of our knowledge, the first paper relating copositivity with boundedness from below was [101]. One must first express the quartic part of the scalar potential, V 4 , as a quadratic form of the n real fields ϕ a (a = 1, 2, . . . n) in the theory, The scalar potential is BFB if and only if the matrix of quartic couplings Λ ab is copositive. A real matrix A is said to be copositive if x T A x 0 for every non-negative vector x 0, that is, x i 0. If the inequality is strict, the matrix is strictly copositive. Therefore, checking for the copositivity of the matrix of quartic couplings would in principle provide sufficient and necessary boundedness from below conditions. However, in complicated models such as the general Scotogenic model, one cannot write V 4 as a quadratic form without introducing mixed bilinears (scalar field combinations involving two different fields). For this reason, this method only leads to sufficient conditions, as we now explain.
In order to write the quartic part of the scalar potential as a quadratic form we define with |ρ ij | ∈ [0, 1] by virtue of the Cauchy-Schwarz inequality. Thus, we can express the boundedness from below condition as with x = h 2 1 . . . h 2 i . . . h 2 ij . . . and the matrix V 4 is given by a combination of the quartic couplings, the λ's, as well as the ρ's and φ phases. 9 The reason why this method provides only sufficient conditions is the presence of the mixed bilinears. Notice that the direction given by h 2 ij is not independent of h 2 i and h 2 j . Therefore, imposing x T V 4 x 0 for every non-negative x vector is overconstraining, since unphysical directions would be included in the test. Nonetheless, when the test is positive, the potential is BFB. In summary, a scalar potential is BFB if the associated V 4 matrix is copositive. However, when the matrix is not copositive nothing can be said about the potential.
There is mathematical work showing that a symmetric matrix A of order n is (strictly) copositive if and only if every principal submatrix B of A has no eigenvector w > 0 with associated eigenvalue κ < 0 ( 0) [102]. However, these theorems are of little practical value when the matrix has a large order, since there will be 2 n − 1 principal submatrices. Luckily, we can make use of [103] instead. The authors of this work proposed an algorithm that leads to necessary and sufficient conditions for the copositivity of unit diagonal matrices (matrices with all diagonal elements equal to 1). Although the algorithm in [103] could only be applied for up to 7 × 7 matrices, incidentally the case in the (1, 3) Scotogenic model, more recent work by the same authors contains indications to extend it to higher orders [104].
After all these considerations, our procedure to check for copositivity is as follows: 1. We replace all the quartic couplings in V 4 by the numerical values in the scalar potential we want to test.
2. We transform each element of the matrix to the worst case scenario. This is achieved by treating the remaining ρ and φ parameters as independent variables and setting them to the configuration for which the term is minimal. 10 3. We check if the matrix has null entries in the diagonal. If it does, we remove the corresponding rows and columns. The original matrix will be copositive if the remaining one is and the removed elements are non-negative. 4. We need the matrix to have unit diagonal to be able to apply the algorithm in [103]. Therefore, we divide all its entries by the smallest element in the diagonal and we replace all the values greater than 1 by 1. The original matrix will be copositive if the new one is.
5. We finally check the copositivity of the resulting matrix with the algorithm in [103].
A final remark about our method is in order. The stability in charge-breaking directions is ignored in many analyses. However, since we are being overly restrictive treating all the ρ moduli and φ phases as independent variables in the different entries of V 4 , charge-breaking directions are included as well in our BFB test. In order to prove it, let us parametrize the scalar doublets of the model under consideration as This parametrization and an example of how to use it to explore boundedness from below is shown in [105]. Let us consider a contraction of scalar doublets φ † i φ j = √ r i r j sin α 1 sin α 2 + cos α 1 cos α 2 e −i(β 1 −β 2 ) , and take the modulus of the term in square brackets sin α i sin α j + cos α i cos α j e iβ 2 = sin 2 α i sin 2 α j + cos 2 α i cos 2 α j + sin α i sin α j cos α i cos α j e iβ + e −iβ = sin 2 α i sin 2 α j + cos 2 α i cos 2 α j + 2 sin α i sin α j cos α i cos α j cos β 1 .
As expected, the product is, at most, as large as the modulus of the fields, √ r i . Therefore, if we treat the factors that multiply √ r i r j as independent variables (that is, being overly restrictive as explained in footnote 10), ρ ij e iφ ij , and make all combinations minimal, our method will cover boundedness from below in charge-breaking directions as well.