Radiative neutrino mass model from a mass dimension-11 $\Delta L =2 $ effective operator

We present the first detailed phenomenological analysis of a radiative Majorana neutrino mass model constructed from opening up a $\Delta L = 2$ mass-dimension-11 effective operator constructed out of standard model fields. While three such operators are generated, only one dominates neutrino mass generation, namely $O_{47} = \overline{L^C} L \overline{Q^C} Q \overline{Q} Q^C H H$, where $L$ denotes lepton doublet, $Q$ quark doublet and $H$ Higgs doublet. The underlying renormalisable theory contains the scalars $S_1 \sim (\bar{3},1,1/3)$ coupling as a diquark, $S_3 \sim (\bar{3},3,1/3)$ coupling as a leptoquark, and $\Phi_3 \sim (3,3,2/3)$, which has no Yukawa couplings but does couple to $S_1$ and $S_3$ in addition to the gauge fields. Neutrino masses and mixings are generated at two-loop order. A feature of this model that is different from many other radiative models is the lack of proportionality to any quark and charged-lepton masses of the neutrino mass matrix. One consequence is that the scale of new physics can be as high as $10^7$ TeV, despite the operator having a high mass dimension. This raises the prospect that $\Delta L = 2$ effective operators at even higher mass dimensions may, when opened up, produce phenomenologically-viable radiative neutrino mass models. The parameter space of the model is explored through benchmark slices that are subject to experimental constraints from charged lepton flavour-violating decays, rare meson decays and neutral-meson mixing. The acceptable parameter space can accommodate the anomalies in $R_{K^{(*)}}$ and the anomalous magnetic moment of the muon.


Introduction
The minimal Standard Model (SM) features massless neutrinos. However, the experimental observation of neutrino oscillations has established that at least two of the three known neutrinos are massive [1][2][3][4][5][6][7][8][9][10][11][12]. These experiments have measured the squared-mass differences ∆m 2 21 ≡ m 2 2 − m 2 1 and |∆m 2 32 | ≡ |m 2 3 − m 2 2 |, but are unable to probe the absolute neutrino mass scale. However, cosmological constraints derived from large-scale structure and cosmic microwave background measurements provide a strong upper bound on the sum of the neutrino masses of about 0.2 eV [13]. Independently of cosmology, β-decay endpoint measurements constrain the absolute mass scale to be at most about 1 eV [14][15][16]. With or without the cosmological constraint, it is clear that the neutrino mass eigenvalues are at least six orders of magnitude smaller than that of the lightest charged fermion, the electron. The neutrino mass problem is the determination of the dynamical mechanism by which neutrino masses are generated and why those masses are so small. All mechanisms require the introduction of as-yet undiscovered fields, and thus constitute physics beyond the Standard Model (BSM). (We will use BSM and "new physics" (NP) interchangeably.) A pivotal question for neutrino mass models is whether or not neutrinos are their own antiparticles. Being electrically neutral, neutrinos are the only Majorana fermion candidates in the SM. Thus, neutrino mass models fall into two categories: Dirac and Majorana. Dirac mass can be generated by introducing right-handed neutrino fields into the low-energy spectrum of the SM. Neutrino mass would then be generated through the same mechanism responsible for all SM fermion masses; however, the smallness of the neutrino masses would simply be due to unusually small Yukawa couplings -an unsatisfying resolution.
Majorana neutrino mass models can provide a more natural explanation. 1 The argument is as follows. All Majorana neutrino mass terms must take the form ν c L m ν ν L + ν L m * ν (ν L ) c , where ν L is a SM left-handed neutrino field and ν c L is the CP conjugate which is equivalent to a right-handed antineutrino field. Since both ν L and ν c L carry a lepton number of +1, these mass terms violate total lepton number by two units (∆L = 2), as necessary when neutrinos and antineutrinos are identical. As we review below, this feature helps us explore neutrino mass models in a systematic way. Now, recall that the quantum numbers for the left-handed lepton doublet, L ∼ (1, 2, −1/2), which contains the left-handed neutrinos, are such that a Majorana mass term breaks SU (2) L × U (1) Y symmetry. This issue can be resolved by introducing exotic fields that exist at an energy scale above the electroweak scale. These heavy exotic fields couple to SM particles in a gauge invariant and renormalisable way, and generate self-energy Feynman diagrams for the left-handed neutrinos at tree-level or loop-level. At energy scales below the electroweak scale, neutrino mass manifests (can be understood) through ∆L = 2 effective operators, obtained by integrating out the exotic heavy fields. The mass terms are then suppressed by the scale of new physics leading to a natural explanation for the smallness of neutrino masses. The three seesaw models [17][18][19][20][21][22][23][24][25], for example, are all UV-completions of the same mass-dimension 5 effective operator, called the Weinberg operator. The seesaw models are wonderfully minimal, however the high BSM scale typically invoked makes them challenging to experimentally test [26][27][28]. We should, therefore, consider other possibilities, for this reason, as well as for the sake of completeness.
Classifying Majorana neutrino mass models using ∆L = 2 effective operators, each of which can be "opened up" (UV-completed at tree level) to produce neutrino self-energy diagrams, is a systematic way to approach the neutrino mass problem. Babu and Leung [29] have published a near-complete list of ∆L = 2 effective operators which may be opened up using exotic fields such as massive scalars, vector-like fermions and massive Majorana fermions [30]. The resulting models generate Majorana neutrino mass either at tree-level or loop-level with most of the operators leading to models that produce the latter. An alternative and complementary approach to neutrino-mass model classification can be structured around loop-level completions of the ∆L = 2 Weinberg-like operators L C LHH(H † H) n [31][32][33][34]. The mass-dimension of an operator is necessarily odd when (∆B − ∆L)/2 is odd [35], where ∆B is the change in baryon number. All effective operators classified by Babu and Leung conserve baryon number and break lepton number by two, thus they all have odd mass-dimension.
Radiative neutrino mass models, which have mass generated at loop level, introduce additional suppression factors alongside the suppression that comes from the masses of the heavy exotic particles. Radiative neutrino mass models are attractive because they naturally produce small neutrino masses for three reasons: i a suppression of 1 (16π 2 ) l , where l is the number of loops in the neutrino self-energy diagram, from the numerical factor which automatically comes with each loop integration, ii a product of couplings which are potentially all smaller than 1 representing the interaction strengths of the exotic particles, and iii a suppression by v Λ p , where v is the vacuum expectation value (vev) of the Higgs field, Λ is the mass scale of the exotica coming from the exotic propagators introduced during the UV-completion, and the exponent p > 0 is model-dependent.
The trend is for the higher dimensional effective operators listed in [29] to include more suppression in the form of i. and ii., thus accordingly decreasing the scale of new physics (NP) needed to produce small neutrino masses. A combination of points i. and iii. is the reason Babu and Leung do not include effective operators of dimension 13 and greater in their list. It was believed that any exotic particles used to complete these models would have to be detectable at an energy scale that has already been probed [29] and therefore, dimension-11 operators that produce neutrino masses in agreement with current data at two-loop level or more would lie in a sweet spot -bringing the scale of BSM physics to a few TeV, an energy scale that is being directly probed at the Large Hadron Collider (LHC) and indirectly at precision-or luminosity-frontier experiments, and would be fully accessible at a future 100 TeV collider. However, in this paper, we present a radiative Majorana neutrino mass model derived from a mass-dimension 11 effective operator with a scale of NP Λ that can be as high as about 10 7 TeV. Our findings suggest that dimension-13, and possibly even dimension-15 effective operators should not be overlooked in the search for viable Majorana neutrino mass models. For a comprehensive review of radiative neutrino mass models and the effective operator method see [36].
In this paper, we present the first detailed radiative Majorana neutrino mass model derived from a mass-dimension 11 effective operator. In Section 2, we define our Model and explain how neutrino masses are generated. Then, in Section 3, we investigate the constraints imposed by experimental results from rare processes involving charged leptons and flavour physics and discuss the results in Section 4. In Section 5, we offer our conclusions.

The Model
We introduce three exotic colour-triplet scalar fields to the particle content of the SM: an SU (2) singlet, S 1 , and two SU (2) triplets, S 3 and φ 3 , with quantum numbers given by where the subscripts indicate the transformation property of the scalars under SU (2) L and the superscripts indicate the electric charge of each component of the exotic scalars. The first entry in the triples specifies the colour multiplet, the second the weak-isospin allocation, and the third the hypercharge, Y , normalised such that electric charge Q = I 3 +Y .
The three exotic scalars listed above generate three separate ∆L = 2, dimension-11 effective operators at tree level, and give rise to radiative Majorana neutrino masses at two-loops. It is important to note that these three scalars do not give rise to any lower dimension ∆L = 2 effective operators at tree-level. Thus, the neutrino self-energy diagrams generated in the UV-completion of these dimension-11 operators with our three exotic scalars will be the leading order contribution to the neutrino mass [37]. In the notation used by Babu and Leung in [29], the operators, depicted in Figure 1, are The scalars S 1 and S 3 can Yukawa-couple as leptoquarks, diquarks, as one of each, or as both. As leptoquarks, they appear together in models which tackle the flavour anomalies in the R K ( * ) and R D ( * ) observables [38][39][40][41][42]. S 1 coupling as a leptoquark is able to explain R D ( * ) (see, for example [43,44]) while S 3 coupling as a leptoquark is able explain R K ( * ) [45][46][47][48][49][50][51]. However, in order to generate neutrino mass in this model, the fermion content of the chosen effective operators forces us to have one of either S 1 or S 3 coupling as a leptoquark, and the other as a diquark. As we discuss later, the choice that leads to neutrino mass generation at an acceptable loop order has S 3 coupling as a leptoquark, and S 1 coupling as a diquark, both with flavour dependent couplings. Consequently, our model can only adequately explain the anomalies resolved by the leptoquark S 3 -specifically those in R K ( * ) . The scalar φ 3 only couples to other scalars and gauge bosons.

The Lagrangian
The general, gauge invariant, renormalisable Lagrangian produced when introducing the three scalars mentioned above can be found in Appendix A. The full Lagrangian has both leptoquark and diquark couplings for S 1 and S 3 , thus explicitly violating baryon number conservation. This is, of course, phenomenologically unacceptable unless the couplings that lead to proton decay are extremely small. In our analysis, we simply impose exact U (1) B symmetry so that baryon number conservation is exact.
Two neutrino mass models emerge once baryon number conservation is imposed on the general, gauge invariant, renormalisable Lagrangian, written in full in Appendix A. Two models are produced by the SM particles together with the three exotic scalar fields with baryon number assignments Model 2, in which S 3 is a diquark and S 1 is a leptoquark, leads to unacceptably small neutrino masses, as will be detailed in Section 2.6. Model 1, in which S 3 is a leptoquark and S 1 is a diquark produces non-vanishing neutrino mass associated with the UV-completion of the dimension-11 operators O 25 , O 47 , and O 55 . Although our model generates all three operators, we now show that, once dressed to become self-energy diagrams, the graph associated with operator O 47 dominates by several orders of magnitude compared to those associated with O 25 and O 55 .
Scale of new physics from effective operators-The generic type of neutrino mass diagrams generated from the UV-completion of each operator can be found in Figure 2. Let us start by analysing O 25 , whose neutrino mass diagram is depicted in Figure 2a. 2 Due to the chirality structure of O 25 , its UV-completions include two mass insertions. Consequently, contributions to the neutrino mass originating from this operator will depend on the mass of the up and down-type quarks in the loop. The neutrino mass will be where f, g, h, and λ are coupling constants, v is the Higgs vev, Y u and Y d are the Higgs Yukawa couplings for the up-type and down-type quarks, and Λ is the scale of new physics. The cubic scalar coupling, µ, we assume to be of the scale of new physics and thus it cancels with a factor of Λ in the denominator.
Operator O 47 is an interesting operator in that it produces neutrino mass contributions that are not constrained by the masses of SM particles. The neutrino mass generated from tree level completions of operator O 47 will look like 2 When the scale of new physics is less than or equal to 2 TeV, neutrino masses generated through the UV-completion of operator O25 will also include an extra contribution from a three-loop diagram obtained by closing the neutral Higgs bosons into a loop.
Finally, operator O 55 produces neutrino mass where g W is the weak coupling, and Y e is the SM charged lepton Yukawa coupling. As O 55 only contains one lepton field, a W boson is required to obtain a second neutrino, generating the self-energy diagram in Figure 2c.
That is, the contribution to neutrino masses coming from the insertions of O 25 is suppressed by at least 10 −3 , and those from O 55 by at least 3 10 −5 compared to O 47 . With this knowledge, we can sensibly approximate the contribution to neutrino mass generated by the UV-completions of O 47 to be dominant, ignoring the contributions associated with the other two operators.
From Equation 2.2 we expect that for indicative couplings f = g = h = λ = 1, the scale of NP is Λ ∼ 10 7 TeV, for couplings f = g = h = λ = 0.1, the scale of NP is Λ ∼ 10 3 TeV, and for couplings f = g = h = λ = 0.01, the scale of NP is Λ ∼ 10 −1 TeV. Thus, to an order of magnitude precision, we can expect that our neutrino mass model is viable with reasonably-valued exotic couplings.
This analysis leads us to conclude that, to a good approximation we can, and do, choose to consider only the neutrino mass diagrams associated with the direct closure of operator O 47 into neutrino self-energies. These are two-loop diagrams, with only exotic scalars and left-handed SM fermions running through the loops, as depicted in Figure 2b.
Even after imposing B-conservation with the Model 1 assignments, the Lagrangian retains a large number of parameters. In order to make exploring that parameter space tractable, we also make the simplifying assumption that all couplings that play no role in neutrino mass generation are zero.
where the coupling between the gauge bosons and exotic scalars, L gauge−S , is defined in equation A.3, S-F represents couplings between the exotic scalar and SM fermions and 4SB, 3SB and 2SB represent scalar only interactions between four, three and two scalar bosons, respectively. In the fermion sector, we define and in the scalar sector, we have (2.7) There are no Yukawa interactions allowed by the SM gauge symmetries between the scalar φ 3 and SM fermions. The τ k , k = 1, 2, 3, are the Pauli matrices; i, j = 1, 2, 3 are generation indices; a, b = 1, 2 are SU (2) flavour indices; ab = (iτ 2 ) ab ; and S k 3 are components of S 3 in SU (2) space. Colour indices are not explicitly shown. The diquark coupling to S 1 , z LL 1 is symmetric due to a combination of the antisymmetry of SU (2) structure and the colour structure of the fermion bilinear. We take the leptoquark coupling, y LL 3 , to be real. In the expansion of the SU (2) structure of the Lagrangian, we have also rotated into the mass eigenbasis of the quarks, using the convention that u i L → (V † CKM ) ij u j L and d i L → d i L . Ultimately, this simply amounts to a definition of the relevant coupling matrices: y LL 3 and z LL 1 . In the scalar sector, the notation [. . .] i indicates that the scalars enclosed couple to form an SU (2) singlet for i = 1 or triplet for i = 3.

Scalar Boson mixing
After electroweak symmetry breaking, Equation 2.7 produces mass mixing between likecharge components of S 3 and φ 3 , generating the mass matrices = cos θ 34 sin θ 34 −sin θ 34 cos θ 34 (2.9) The mixing angles θ 12 and θ 34 are related to the squared mass parameters µ 2 S 3 , µ 2 φ 3 and off-diagonal parameters in the mass matrix through There are seven BSM physical scalar states in our theory (one diquark, five leptoquarks and one coloured, electrically-charged scalar that does not couple to SM fermions) with squared masses given by (2.11) The squared masses of all physical particles are required to be positive, placing a bound on λ φ 3 S 3 H such that Thus, the parameters µ S 3 , µ φ 3 and λ S 3 φ 3 H determine the masses of the five leptoquarks and one charged scalar, and the mixing angles. We can also derive the relationship between the mixing angle and the squared mass difference of the leptoquarks involved in mixing from Equations 2.10 and 2.11, The UV-completion of operator O 47 with the introduction of three exotic scalars; S 1 , Closing the loops by joining the quarks leads to a neutrino mass diagram.
The neutrino self-energy, after electroweak symmetry breaking and rotating into the mass basis of the exotic scalars. Figure 3: The UV-completion of operator O 47 by the exotic scalars S 3 (leptoquark), S 1 (diquark) and φ 3 and its closure forming the neutrino self-energy Feynman diagrams for our model.

Neutrino mass generation
As discussed in Section 2.1, the dominant contributions to the neutrino mass for Model 1 are two-loop neutrino self-energy graphs, with exotic scalars and left-handed SM fermions running through the loops, generated by the completion of operator O 47 4 . The UV-completion of O 47 by the diquark S 1 , the leptoquark S 3 , and the scalar φ 3 , is depicted in the tree-level diagram of Figure 3a. Joining the quark Q to Q lines, and the second Q to Q C , gives two-loop self-energy diagrams, which generate the neutrino mass matrix. There is no mass insertion necessary in the quark lines for a chirality flip, thus the neutrino mass matrix generated by this model is not proportional to the mass of any SM fermion; this interesting feature characterises this model. In terms of the physical mass eigenstates, S 1 , r 1 , r 2 , r 3 , and r 4 , there are eight diagrams; half of them are obtained by reversing the flow of charge arrows in both loops of Figure 3c. Individually, each diagram is divergent. However, due to the absence of a bare neutrino mass in the Lagrangian, the divergences are guaranteed to cancel.
The neutrino mass matrix is obtained from the flavour sum of the self-energy diagrams with the freedom to choose the external momentum to be zero: (2.14) The PMNS matrix can then be used to obtain the physical masses of the neutrinos, the factor of 9 is a QCD colour factor, and the loop integral, I kl , is .
The tensor structure of the numerator arises from the chiral projection operators at the vertices, and the lack of proportionality to the SM fermion masses is a good sanity check when cross-referenced with the lack of the relevant mass insertions in the self-energy diagram in Figure 3. Although m r 1 = m r 3 and m r 2 = m r 4 at our level of approximation, we denote these individually so that the correspondence between terms in the loop integral and diagrams is manifest.
When evaluating the integral it is convenient to work in terms of the dimensionless parameters The numerator of Equation 2.15 can be rewritten as: where the antisymmetric commutator term vanishes due to the fact the integral is µ, ν symmetric, and I 4 is the 4 × 4 identity matrix in Lorentzspinor space. Factoring out m 2 S 1 and rescaling the momenta to dimensionless quantities allows us to write the integral as .

(2.17)
This is a sum of four integrals, each of which is evaluated in Appendix B, both in full generality and in the sensible limit that the quark masses are much smaller than the scale of new physics, m S 1 , i.e. in the limit s k , s l → 0. In this limit, the integral is independent of k, l, and we thus obtain I kl = I ∀k, l, which we calculate to be whereĝ(t α , t β ) andĝ(t α , 0) are defined in Appendix B, specifically in Equations B.5 and B.9 respectively. The behaviour of this integral for leptoquark mass parameter µ S 3 ranging from 1.1 − 100 TeV is shown in Figure 4. A combination of the loop suppression factor and suppression coming from the mass of the heavy exotic scalars allow the integral to give neutrino mass a substantial suppression. For this plot, the other exotic mass parameters and the quartic and cubic scalar coupling values have been fixed: Figure 5 shows an example of the calculated sum of neutrino masses for leptoquark mass parameter µ S 3 ranging from 1.1−100 TeV. The other parameters are set as for Figure 4, with leptoquark and diquark Yukawa couplings set to y LL 3ij = 0.001 and z LL 1ij = 0.01 respectively. Note that neutrino mass goes to zero as µ S 3 → ∞, as expected. It should be understood that Figure 5 shows only the typical scale of the neutrino mass, and our model has enough freedom to allow for more precise fitting to the experimental results, including the correct mass differences between the neutrino mass states -which can be achieved by enforcing a relationship between the leptoquark and diquark coupling matrices, as described below in Section 2.4.

Casas-Ibarra Parametrisation
After these simplifications the neutrino mass, in the flavour basis, is where we have absorbed all constants into m 0 = 18m S 1 S 3 φ 3 I and, for convenience, we define the dimensionless matrix κ ≡ mν m 0 . For a given κ, we thus see that the two coupling matrices, y LL 3 and z LL 1 , must be related. Their relationship can be obtained using the parametrisation method originally described by Casas and Ibarra [52]. Recalling that the diquark couplings z LL 1 must be symmetric, we can use Takagi's factorisation method to diagonalise z LL  that (2.20) Multiplying both sides of the equality by D −1 κ on the left and the right, we get where 1 3×3 is the identity matrix. This implies that where R is an orthogonal matrix (in general with complex entries). Thus, to produce the measured light neutrino masses contained in D κ , with mixing parameters contained in U PMNS , the most general leptoquark coupling is given by (2.23) Figure 5: Plot of the sum of neutrino masses, m ν , in eV, for leptoquark mass parameter µ S 3 ranging from 1 − 100 TeV. The other exotic mass parameters and the quartic and cubic scalar coupling values have been fixed: All entries in the leptoquark and diquark coupling matrices have been set to y LL 3ij = 0.001 and z LL 1ij = 0.01 respectively. The solid horizontal orange line indicates the experimental lower bound of 0.06 eV, while the green line shows the approximate cosmological upper limit of 0.2 eV.
Based on this equation, we see that y LL 3 depends on the known low-energy parameters contained in D κ and U PMNS , as well as the following free parameters: six real parameters from the symmetric diquark coupling, z LL 1 , and three, generally complex, parameters in R. Alternatively, since we initially place constraints on y LL 3 , we can rearrange Equation 2.23 to find z LL 1 as a function of y LL 3 , such that Notice that, due to its symmetric nature, z LL 1 is independent of the orthogonal matrix R. This makes sense since we still have nine free parameters, now all contained in y LL 3 .

Parameters and notation for analysis
Our model, given the simplifying assumptions, has 14 free parameters: four coming from the mass-dimension 1 couplings µ S 1 , µ S 3 , µ φ 3 , and m S 1 S 3 φ 3 , and another 10 coming from the dimensionless coupling constants z LL 1 , y LL 3 , and λ S 3 S 1 H , with z LL 1 and y LL 3 being related by Equation 2.24. In Section 3, we discuss several phenomenological constraints on the leptoquark couplings. From here on in, we will simplify our notation such that leptoquark couplings read (y LL 3 ) ij ≡ y ij , with the index i (j) representing the generation of the contributing quark (lepton). We will similarly denote diquark couplings by (z LL 1 ) ij ≡ z ij .
Given that 14 parameters is too large a space to sample properly, and the results would be difficult to visually present, we are forced to fix the majority of the parameters at benchmark values. We choose to scan over four leptoquark coupling constants, y 11 , y 12 , y 21 and y 22 , and one mass parameter, µ S 3 . The benchmark values allocated to µ S 1 , µ φ 3 , m S 1 S 3 φ 3 and λ S 3 S 1 H can be found in Table 1. In order to give a representative idea of our model's robustness as well as investigate a variety of possible conclusions drawn from future particle experiments, Table 1 also includes three benchmark textures for the leptoquark coupling matrix.
have been fixed to the values indicated. Three different leptoquark coupling matrix textures were used to scan the parameter space, each matrix texture having five of the nine leptoquark couplings fixed, as above, while the other four leptoquark couplings, specifically y 11 , y 12 , y 21 and y 22 , are scanned over. In texture A, A = 10 −5 , in texture B, B = 10 −3 , and in texture C, C = 10 −1 .

Model 2 and vanishing neutrino self-energies
Before investigating the phenomenology of Model 1, we will use this section to tie up the loose end of the discarded alternative completion of O 47 . Recall that in Model 2 S 3 couples as a diquark and S 1 couples as a leptoquark. When integrated out, these exotic fields give rise to an operator with SU (2) Working in two-component Weyl spinor notation, the specific Lorentz structure generated is (LL)(QQ)(Q † Q † )HH, where parentheses indicate contracted spinors. In Appendix C we present the calculation of the two-loop contribution to the neutrino mass in this model. Curiously, we find that the neutrino masses vanish due to the symmetry properties of the integral and the antisymmetry of a set of couplings. (This antisymmetry is enforced by Fermi-Dirac statistics.) This suggests that the neutrino mass for this model arises at some higher-loop order. Below we show that the leading-order contribution to the neutrino masses arising from this particle content vanishes. This leading-order argument contains essentially the same ingredients as those required to see the behaviour in the UV theory, and we point readers to the calculation in the appendix for more detail.
Writing all indices (SU (2), SU (3), Lorentz, and flavour) explicitly, Model 1 generates at tree level, where Greek letters (α, β, . . .) represent spinor indices, Latin letters from the middle of the alphabet (i, j, . . .) represent SU (2) indices, Latin letters from the end of the alphabet (r, s...) represent flavour indices and capital letters represent SU (3) indices. The neutrino masses arising from a single insertion of operators of this type will vanish since they depend on integrals with an odd number of loop momenta in the numerator [53]. We thus consider neutrino masses arising from insertions of the dimension-13 operator with a derivative acting on each Q † . This operator will also be generated by the particle content of the UV-completion of operator O k 47 . However, we now show that this contribution also vanishes.
We first show that the operator must be anti-symmetric under exchange of the v and w flavour indices: which confirms that the diquark coupling of S 3 is anti-symmetric in flavour as stated in [54]. The neutrino mass generated by this operator is then represented by is a Wilson coefficient obtained from the evaluated self-energy diagrams. The square brackets indicate the anti-symmetry under interchange of v and w discussed above.
It should be noted that O k 47 is missing from the list of ∆L = 2 effective operators listed in [53]. This may be due to an implicit assumption that the number of SM fermion generations is not more than one. In this case O k 47 itself vanishes since Q m Q n H m H n = 0.
One might worry about the validity of this claim in light of the "extended black box" theorem [55], which states that any non-vanishing ∆L = 2 effective operator leads to nonvanishing Majorana neutrino mass. This is remedied by the fact that we are only closing off the effective operator in the simplest way to generate neutrino masses. The theorem tells us that there must be non-zero contributions to neutrino mass coming from O k 47 since the operator itself does not vanish when flavour is considered. We therefore surmise that neutrino masses arise at higher loop order, and are probably too small to meet the lower bound of |∆m 2 32 | 0.05 eV with phenomenologically acceptable exotic particle masses.
There are several other ∆L = 2 effective operators which exhibit the same property, including, but not limited to O 11b , O 12a and O 48 . The two-loop contributions coming from completions of all these operators vanish, implying that there must be nonzero higher-loop contributions. Similar remarks about the 0.05 eV lower bound pertain. This observation could potentially be used to eliminate a sizable number of effective operators from the pool of neutrino-mass-model candidates.

Constraints from Rare Processes and Flavour Physics
In this section, we investigate the phenomenology of our Model 1, and place constraints on the values of the coupling constants responsible for generating neutrino mass. This investigation is conducted in three parts. First, the leptoquark couplings y ij are constrained via the model's BSM contribution to rare processes of charged leptons, including µ to e conversion in nuclei, the decays µ → eγ and µ → eee, and the anomalous magnetic moment of the muon. Second, the leptoquark couplings are constrained via BSM contributions to rare meson decays. Finally, the diquark couplings, z ij , are constrained via experimental results from neutral meson anti-meson mixing.

Rare processes of charged leptons
In the absence of neutrino flavour oscillations, lepton number is conserved in the SM. While lepton flavour has been shown to be violated by neutrino oscillations, it has as-yet not been observed in the charged lepton sector. The lepton flavour violating (LFV) terms in our Lagrangian are thus constrained by charged LFV processes. The most stringent upper bounds on LFV processes in leptoquark models come from µ → eee and µ → eγ decays and µ − e conversion in nuclei.

µ → e conversion in nuclei
The strongest bounds on the branching ratio Br(µ → e) come from µ → e conversion off titanium and gold nuclei. The current constraints, which were set by the SINDRUM Collaboration [56,57], are of order 10 −12 (Table 2), with future experimental sensitivities predicted to improve by several orders of magnitude. The most promising are the COMET [58] and Mu2e/COMET [59] experiments, aiming for sensitivities of order 10 −16 , and the PRISM/PRIME proposal [60], boasting a possible sensitivity of 10 −18 .  µ , the branching ratio at 90% confidence level and the total capture rates for different nuclei.
The most general interaction Lagrangian for this process, in the notation of [61], is where G F is the Fermi constant, m µ is the muon mass, and the A L,R and g's are all dimensionless coupling constants corresponding to the relevant operators.
The branching ratio is defined to be where ω conv is the µ to e conversion rate, and ω capt is the total muon capture rate. The conversion rate, ω conv is calculated from the effective Lagrangian in Equation 3.1 to be where S, D and V are overlap integrals, and n and p superscripts refer to processes interacting with a neutron or proton respectively. The coefficients G q,p and G q,n , associated with S (p,n) , are calculated in [62], but do not play a role in our model. This is due to the fact that µ → e conversion does not generate scalar operators in our model, thus g LS(q) = 0. Similarly, the coefficient associated with tensor operators, namely A R , vanishes. Accordingly, we only provide values relevant to our model: the V overlap integral values for titanium and gold can be found in Table 2, while other values can be found in [61].
In our model, the dominant contributions to µ − e conversion in nuclei come from diagrams with the leptoquark S 3 mediating interactions between the charged leptons and the  Figure 6: Tree level processes contributing to µ to e conversion in nuclei. Note the notation in the left diagram represents two diagrams, one mediated by r 1 and another by r 2 . The arrows here, as in all other diagrams in this paper. represent the chirality of the field. Arrows pointing towards the vertex represent left-handed fields.
three lightest quarks, as can be seen in Figure 6. The effective Lagrangian, calculated using Feynman rules for fermion number violating interactions found in [63], is and (3.7) The matrix element involving strange quarks vanishes as coherent conversion processes dominate and the vector coupling to sea quarks is zero. This leads to For fixed leptoquark masses, this process places the most stringent constraints on the product y 11 y 12 through ωconv ωcapt < 7.0 × 10 −13 for gold and < 4.3 × 10 −12 for titanium.

µ → eγ
The most stringent constraints on this process are obtained from the non-observation of LFV muonic decays by the MEG experiment [64], with a measured branching ratio of Br(µ → eγ) = 4.2 × 10 −13 at 90% CL. Future prospects are looking to improve on this by an order of magnitude. Specifically, the MEG-II experiment [65,66] is predicted to start searching for µ → eγ decays this year, with a target sensitivity of 4 × 10 −14 .
The effective Lagrangian for µ → eγ is where l = µ and l = e, σ µν = i 2 [γ µ , γ ν ], F µν and σ L(R) are Wilson coefficients. The partial decay width for µ → eγ is (3.10) Our model will have contributions from the leptoquark mass states r 1 and r 2 with up-type quarks running in the loop and leptoquark S

4/3 3
with down-type quarks running in the loop, for each of the four diagrams in Figure 7. In total there are 36 contributing diagrams leading to Wilson coefficients  leptoquark and down-type quarks, or r 1 and r 2 (electric charge 1/3) and up-type quarks running through the loop.
We thus obtain the following constraint on the leptoquark coupling constants for first and second generation leptons: where Γ tot µ = 2.99 × 10 −19 GeV.

Anomalous magnetic moment of the muon
The SM predicts the anomalous magnetic moment of the muon to be a SM µ = 1.16591803(70)× 10 −3 [67,68], while the most precise experimental measurement, which comes from the E821 experiment, is a exp µ = 1.16592080(63) × 10 −3 . The difference between the SM prediction and the experimental measurement, ∆a µ ≡ a exp µ − a SM µ = (2.8 ± 0.9) × 10 −9 , suggests the possible presence of BSM contributions. In our model, the leptoquark couplings with the muon provide such a contribution, given by (3.14) The contribution lies inside the bounds of ∆a µ , ameliorating the anomaly, without placing a strong constraint on the leptoquark couplings involved. This is consistent with previous results found in literature [69][70][71]. However, when combined with other leptoquark solutions, the leptoquark S 3 has been shown to explain the discrepancy between theory and experiment in the anomalous magnetic moment of the muon [72].

µ → eee
To date, the strongest constraint on Br(µ → eee) remains the 1.0 × 10 −12 achieved by the SINDRUM collaboration in 1988 [73]. Looking ahead, the Mu3e collaboration [74] is promising to improve the current constraint by four orders of magnitude.
The interaction Lagrangian for this process involves interactions between the S 3 leptoquark and both the gauge sector and the quark sector. At one-loop level, µ → eee decays receive contributions from three types of Feynman diagrams: γpenguins, Z-penguins and box diagrams, as depicted in Figure 8. Thus, the µ → eee probability amplitude consists of three parts Photon-penguins -The µ → eee photon penguin diagrams closely resemble the µ → eγ decay diagrams, however this time the photon is internal, and thus not on-shell. The amplitude for the µ → eee photon penguin diagrams is [75][76][77] A with the Wilson coefficients as follows: (3.17) The u f and v f are the usual free-particle spinors. The variables in the loop functions are ratios of the squared masses of the quarks and the leptoquarks: and t ji = m 2 u i /m 2 r j and the loop functions F 1 (x) and F 2 (x) are leptoquark and down-type quarks, or r 1 and r 2 (electric charge 1/3) and up-type quarks running through the loop. The situation is similar for the box diagrams, except there will also be contributions where r 1 and r 2 are simultaneously present in the loop. must also consider all possible combinations of leptoquark mass states and quarks running through the loop.

Z-penguins -The amplitude of the Z-penguin diagrams is
with the Wilson coefficients: (3.20) Here θ W is the weak angle, and x i and t ji are as above. We also note that the gauge coupling between the Z-boson and leptoquark mass states includes flavour changing contributions. Consequently, the mass states which involve mixing, specifically r 1 and r 2 , must be treated together when considering the coupling between the Z-boson and the leptoquark S 3 in the mass basis. The loop functions in 3.20 are: (3.21b) In the limit of vanishing mixing angle and x 1 → x 2 we have the following simplification (

3.22)
Box diagrams -For our model, the non-vanishing amplitude from the contribution of box diagrams to the µ → eee decay is with (3.24) The loop function for the box diagrams is defined as for i = 1, 2, 3.
µ → eee Amplitude -Using the form factors defined above, we calculate the µ → eee decay rate to be with (3.27) Thus, we obtain a strong constraint on the leptoquark coupling constants to first and second generation leptons, given by where Γ tot µ = 2.99 × 10 −19 GeV.

Rare meson decays
In the SM, rare meson decays in the form of flavour-changing neutral currents (FCNCs) arise at loop level and are thus heavily suppressed, leading them to be highly sensitive to BSM contributions. A plethora of precision experiments have placed stringent bounds on these rare decays. These processes occur at tree-level in our model, thus the couplings involved are severely constrained.
Carpentier and Davidson [78] published a comprehensive list of (order of magnitude) constraints on two-lepton-two-quark (2l2q) operators. They work with the effective Lagrangian (3.29) where C ijkn /2m NP are the Wilson coefficients. The coefficient relevant to our model is that accompanying a dimension six, left-handed chiral vector effective operator (3.30) The bounds are set on dimensionless coefficients, ijkn (n)lq , related to the respective Wilson coefficients by Bounds are placed on ijkn by analysing the contribution of a relevant effective operator to the branching ratios of rare meson decays, one effective operator at a time. In doing so, there is a risk of overlooking possible destructive interference effects. Therefore, the constraints in this section are only order of magnitude estimates. The analysis discussed in this section allows us to place constraints on all nine leptoquark couplings, y ij , through simply calculating the contribution of leptonic rare meson decays and semi-leptonic neutral current decays to 2l2q effective operators. The results are summarised in Table 3.

Leptonic meson decays
Starting with the leptonic meson decays K 0 L → eµ and K 0 L → µµ we place bounds on products of the first two generations of leptoquark couplings. The Wilson coefficient associated with the operator O 1212 The relevant Feynman diagram, which contributes to the process K 0 L → eµ, is depicted in Figure 9. This process places constraints on the first and second generation diagonal leptoquark couplings. Since the constraints on ijkn are equivalent under the exchange (l i γ µ P L l j )(q k γ µ P L q n ) Br(K + →π + νν) Br(K + →π 0 eνe)

Constraint on
Tab. 2 and Tab. 12]. The left most column specifies the generation indices ijkn, the second column gives the best constraint on ijkn , obtained from the observable indicated in the third column and the experimental bounds given in the last column. The bounds also apply under permutation of lepton and/or quark indices. of quark or lepton indices, the same constraints also apply to the product y 12 y 21 . Similar bounds can be placed on y 12 y 22 from the process K 0 L → µµ, with Feynman diagram depicted in Figure 9.

Semi-leptonic meson decays
The semi-leptonic meson decays, the most tightly constrained being K + → π + νν and B + → K + νν, place bounds on third generation leptoquark couplings, as well as additional bounds on the couplings already discussed. The process K + → π + νν, depicted in Figure 10a, induces the following Wilson coefficient, associated with the operator O ij12 Figure 9: Figure 10: Dominant contributions to (a) K + → π + νν and (b) B + → K + νν.
where we sum over the neutrino flavour. The bounds on ijkn , found in Table 3, are calculated one flavour at a time, with all other contributions set to zero. After this analysis, the only unconstrained leptoquark parameters are those involving third generation quarks.
The leptoquark couplings to third generation quarks can be constrained via the process The Feynman diagram associated with this operator is pictured in Figure 10b, and the Wilson coefficients C ij23 m 2

LQ
, are identical to Equation 3.33 apart from the quark indices.
The constraints applied to each set of parameters are summarised in Table 4. Note that we do not include constraints from the anomalous decay B s → µµ, which will be discussed in Section 3.4.

Neutral meson anti-meson mixing
Mixing of neutral mesons occurs in the SM through box diagrams with W -bosons and top quarks as the propagators. Since neither S 3 nor S 1 have any restrictions with respect to the generation of the SM fermions they couple to, meson mixing gets contributions from diagrams with leptoquark couplings as well as diquark couplings. Another consequence of unrestrained flavour couplings is that there are contributions to the mixing of all neutralmeson species. We focus on neutral kaon mixing, B s − B s and B d − B d mixing. The most general effective Hamiltonian for neutral meson mixing is with the effective operator relevant to meson mixing both in the SM and in our model being The SM Wilson coefficient for meson mixing is where m t is the mass of the top quark, m W is the mass of the W boson and S 0 (x) is the Inami-Lim function [80], (3.37) As can be seen in Figure 11, which depicts the NP Feynman diagrams for kaon mixing, our model contributes to meson mixing through box diagrams with both diquark and leptoquark propagators. When considering the contributions from the diquark, S 1 , we only include contributions which include the top quark as a propagator, just as in the SM calculation. This is because the CKM matrix elements involving the top quark dominate over the others.  Table 4: Constraints applied to expression involving leptoquark couplings and masses, derived from leptonic and semi-leptonic rare meson decays. Figure 11: Dominant contributions to K − K mixing.
Leptoquark contributions occur through box diagrams with either neutrinos and r 1,2 or charged leptons and S (3.38) The effective operators and Wilson coefficients depend on the renormalisation scheme and scale. However, since we are only interested in an order of magnitude estimate, we neglect the running from the scale of new physics to the top quark mass and take the ratio of the new physics contribution to the SM contribution. This ratio is independent of QCD running and is simply a ratio of the respective Wilson coefficients, The UTfit collaboration has published model-independent constraints on ∆F = 2 oper-ators [79], with the results from the latest fit published in [81]: (3.40b) . (3.40c) The current best fit values for these parameters are given by [81]

Solving the R K ( * ) anomaly
While not the focus of this paper, another key feature of our model is its ability to explain the R K( * ) flavour anomalies due to the presence of the S 3 scalar leptoquark. The R K ( * ) flavour anomalies are a set of deviations from SM predictions in the decays of B-mesons. The anomalous quantities are the ratios of branching fractions The SM prediction of R K ( * ) = 1.0003 ± 0.0001 [82] is close to unity due to lepton flavour universality, with the only difference in the measured branching fractions coming from their dependence on the masses of the final state leptons. Most recently, the LHCb collaboration has updated the measurement of R K by combining Run-1 data with 2 fb −1 of Run-2 data [83]. They found where q 2 is the dilepton invariant mass squared, the first uncertainty is statistical and the second is systematic. This result continues to be in tension with theory at 2.5σ. The R K * results have also been updated, with preliminary measurements by Belle [84] of  Interestingly, the large uncertainties in the Belle measurement allow it to be in agreement with both the SM and the LHCb measurement, which deviates from the SM prediction by ∼ 2.5σ.
Our model contributes to R K ( * ) via semi-leptonic B-meson decays mediated by S 3 . The b → sll transition can be described by the effective Hamiltonian where α e is the fine structure constant. The relevant effective operator for our model, presented here in the chiral basis, is Reference [86] argues that contributions to the effective operator O µµ LL with Wilson coefficient ameliorates the R K ( * ) anomalies, with the uncertainties giving the 1(2)σ bounds. In our model, this corresponds to a central value of (3.52)

Results
We will now present the results of our random parameter scans, and discuss the predictions and limitations of the model. The discussion is broken up into three sections, each discussing one of the three leptoquark coupling matrix textures introduced in Table 1. The strongest constraint on the model comes from µ → e conversion in gold nuclei, which is mediated by leptoquark S 3 at tree level. Consequently, we will also explore the potential consequences for the model if future µ → e experiments fail to measure a signal at the promised prospective sensitivities.
The scans were performed over leptoquark couplings y 11 , y 12 , y 21 and y 22 , and the leptoquark mass parameter µ S 3 , with all other couplings fixed. In 2018 the CMS experiment at the LHC set an exclusion on diquark masses below 7.2 TeV at 95% confidence interval [87]. While the exclusion was calculated for a particular diquark which features in superstring inspired E 6 models [88], similar bounds are expected to apply to other diquarks. The CMS experiment has also set limits on leptoquarks masses, with masses below 1.1 TeV being excluded at 95% confidence interval for third generation leptoquarks decaying to bν [89].
We are not aware of current limits on exotic scalars that only couple to other scalars and gauge bosons, such as φ 3 . With these considerations in mind, we stick to conservative lower bounds that have the potential to be directly probed at the LHC and indirectly probed at precision-or luminosity-frontier experiments. We fix the diquark mass parameter in our model to a benchmark value of µ S 1 = 7.5 TeV, the mass parameter for the scalar φ 3 to µ φ 3 = 1.5 TeV, and scan over the leptoquark mass parameter such that 1.1 ≤ µ S 3 ≤ 10 TeV. The limits on leptoquark masses are set assuming sufficiently large leptoquark couplings of y ≥ 10 −7 to guarantee prompt decay of leptoquarks in the detector. We also found that due to the relationship between y LL 3 and z LL 1 displayed in Equation 2.24, leptoquark couplings below 10 −5 had a high likelihood of requiring large, non-perturbative diquark couplings in order to guarantee the desired neutrino masses. Thus, the free leptoquark couplings: y 11 , y 12 , y 21 and y 22 are allowed to vary between 10 −5 and 1. The other five leptoquark couplings are set to benchmark values. We investigate three leptoquark matrix textures, as indicated in Table 1 and detailed below. . The orange region of the parameter space is constrained by Br(µ → e) Au , the teal region is allowed; specifically it passes all constraints discussed in Section 3 as well as perturbativity constraints on the diquark coupling constants z LL 1ij . The light blue region is excluded by the perturbativity constraints on diquark couplings. The black (pink) region solves the R K ( * ) anomalies to 1(2)σ.

Texture A
The leptoquark matrix texture investigated first is Texture A y =    y 11 y 12 A y 21 y 22 with A = 10 −5 . If leptoquark couplings to third generation quarks and leptons are very weak, as they are for Texture A, ample parameter space is available when accounting for current constraints, including regions which allow the model to explain the R K ( * ) flavour anomalies, as summarised in Figures 12, 13 and 14. In fact, with leptoquark couplings of the order of ∼ 10 −5 , this neutrino mass model is able to explain R K ( * ) at a 1σ level even if NP is not discovered by future experiments, such as Mu2e/COMET and PRISM. This is evident in Figure 12, where viable parameter space is plotted in teal, points that solve the R K ( * ) anomalies to 1(2)σ are black (pink) and the dotted orange (yellow) lines represent prospective constraints for Br(µ → e) Au from the Mu2e/COMET (PRISM) experiments. The orange region of the parameter space is constrained by Br(µ → e) Au , the teal region is allowed; specifically it passes all constraints discussed in Section 3 as well as perturbativity constraints on the diquark coupling constants z LL 1ij . The light blue region is excluded by perturbativity constraints on diquark couplings. Figure 14: Indicative plot of the allowed parameter space in the y 22 vs. y 11 plane, for leptoquark coupling matrix Texture A. The yellow section is ruled out by Br(K → µe), the orange section is ruled out by Br(µ → e) Au , while the teal region is allowed parameter space; specifically it passes all constraints discussed in Section 3 as well as perturbativity constraints on the diquark coupling constants z LL 1ij . The light blue region is excluded by perturbativity constraints on diquark couplings. Figures 12 and 13 both contain a curious feature: a region of parameter space (in the form of a band in µ S 3 ) which is excluded due to perturbativity constraints placed on the diquark coupling constants. This feature can be understood as follows. Small leptoquark couplings lead to a higher probability of non-perturbative diquark couplings. This is a consequence of Equation 2.24, which parametrises the diquark couplings in terms of the leptoquark couplings to ensure the desired neutrino masses are computed. The value of z LL 1 is also inversely proportional to the value of the integral. The curious feature of the band of excluded parameter space in the µ S 3 ∼ 2 − 4 TeV, displayed in Figures 12 and 13, corresponds to the sign change in the integral. This can be verified in Figure 4. Simply put, there exits a region of µ S 3 for which the value of the leptoquark coupling and the integral are sufficiently small, making the diquark couplings z LL 1ij too large to be perturbative.
The R K ( * ) anomalies can be solved for a sub-region of the allowed parameter space, for y 11 10 −2 and 10 −3 y 22 10 −1 . There is a band around y 22 10 −2 for which the model cannot explain the R K ( * ) anomalies, which corresponds to the non-perturbative band in Figure 13. Figure 14 also shows that while µ → e conversion is found to be the most constraining process for this model in general as it appears at tree level, there are other important signals. Since the neutrino mass of this model does not depend on quark masses, there is no reason to set couplings to first generation quarks to zero, meaning our model is also sensitive to probes involving first generation quarks, such as K 0 L →ēµ and K + → π +ν ν. The process K 0 L →ēµ is the most constraining for our model, with Figure 14 showing that the leptoquark couplings y 11 and y 22 cannot be simultaneously close to unity. There exists a trade off due largely to constraints coming from Br(K → µe). It should be noted that µ → e conversion, shown in orange in Figures 14 and 13, is also strongly constraining in the regions excluded by perturbativity and Br(µ → e) Au constraints.

Texture B
The leptoquark matrix texture considered next is Texture B where B = 10 −3 . For this leptoquark matrix coupling texture, constraints due to the perturbativity of the diquark couplings are no longer a major concern. Texture B gives similar results to Texture A in that it is strongly constrained by Br(µ → e) Au and Br(K → µe), yet still able to solve the R K ( * ) anomalies. In fact, with leptoquark couplings to third generation quarks of order 10 −3 , there exists parameter space which is able to solve the R K ( * ) anomalies to 1σ(2σ) even if NP is not discovered by future experiments such as Mu2e/COMET and PRISM, as can be seen in Figure 15 in black (pink). Figure 15 also shows, in light green, that we now have a constraint from Br(µ → eγ) on the leptoquark The light green region shows parameter space that is excluded by Br(µ → eγ). The orange region of the parameter space is constrained by Br(µ → e) Au , the teal region is allowed; specifically it passes all constraints discussed in Section 3 as well as perturbativity constraints on the diquark coupling constants z LL 1ij . The black (pink) region solves the R K ( * ) anomalies to 1(2)σ.
As for Texture A, leptoquark couplings of y 11 10 −2 and 10 −3 y 22 10 −1 can currently solve the R K ( * ) anomalies, as depicted in Figure 17. In can also be seen that y 22 must be less than unity, or the model is not viable at all due to constraints from µ → e conversion. Constraints from Br(K → µe) enforces a trade off between large y 11 and y 22 couplings for Texture B. additionally, Figure 16 shows that y 22 10 −1 is necessary for µ S 3 1.2 TeV. The orange region of the parameter space is constrained by Br(µ → e) Au , the teal region is allowed; specifically it passes all constraints discussed in Section 3 as well as perturbativity constraints on the diquark coupling constants z LL 1ij . The light green region, which in this plot has been superimposed on the orange, shows parameter space that is excluded by Br(µ → eγ). Figure 17: Indicative plot of the allowed parameter space in the y 22 vs. y 11 plane, for leptoquark coupling matrix Texture B. The yellow section is ruled out by Br(K → µe), the orange section is ruled out by Br(µ → e) Au , while the teal region is allowed parameter space; specifically it passes all constraints discussed in Section 3 as well as perturbativity constraints on the diquark coupling constants z LL 1ij . The light green region shows parameter space that is excluded by Br(µ → eγ). The orange region of the parameter space is constrained by Br(µ → e) Au , the teal region is allowed; specifically it passes all constraints discussed in Section 3 as well as perturbativity constraints on the diquark coupling constants z LL 1ij . The black (pink) region solves the R K ( * ) anomalies to 1(2)σ. where C = 10 −1 . When couplings between electrons or muons and third generation quarks are ≥ 10 −2 there is a strong constraint on the leptoquark mass parameter coming from Br(µ → eγ), which excludes parameter space less than µ S 3 ≈ 4 TeV, as seen in Figure 18. A similar bound in the y 22 versus µ S 3 slice of parameter space comes from the Br(µ → e) Au constraint. This can be seen in Figures 19 and 20. Constraints from Br(K → µe) still enforce that y 11 and y 22 cannot be simultaneously large, but it is the Br(µ → e) Au constraint that places a bound on the leptoquark couplings such that y 22 0.5.
The model can solve the R K ( * ) anomalies to 1(2)σ for µ S 3 4 TeV, y 11 10 −2 and 10 −2 y 11 10 −1 , as can be seen in black (pink) in Figures 20, 18 and 21. Unsurprisingly, regardless of the leptoquark coupling matrix texture, constraints from processes involving first generation quarks and leptons, expecially Br(µ → e) provide the strongest bounds on the parameter space of our model and would also be the most promising signals.  : Indicative plot of the allowed parameter space in the y 22 vs. y 11 plane, for leptoquark coupling matrix Texture C. The yellow section is ruled out by Br(K → µe), the orange section is ruled out by Br(µ → e) Au , while the teal region is allowed parameter space; specifically it passes all constraints discussed in Section 3 as well as perturbativity constraints on the diquark coupling constants z LL 1ij . Figure 21: Indicative plot of the allowed parameter space in the y 11 vs. m S 3 plane, for leptoquark coupling matrix Texture C. The teal region shows allowed parameter space, while the orange region indicates parameter space excluded by constraints on Br(µ → e) Au , while the yellow region show parameter space excluded by Br(K → µe). The black (pink) points indicate bounds on the parameters y 11 and µ S 3 needed to explain the R K ( * ) anomalies to within 1(2)σ.

Conclusion
There are a large number of candidate radiative Majorana neutrino mass models. The ∆L = 2 interactions responsible may be classified according to the dominant low-energy effective operators generated at tree-level from integrating out the massive exotic fields. The baryon-number conserving ∆L = 2 operators occur at odd mass dimension, and an extensive list has been compiled up to mass dimension 11 [37]. An interesting question is: Beyond what mass dimension are phenomenologically viable models no longer possible? One generally expects that the higher the mass dimension, the higher also will be the number of vertices and loops in the neutrino self-energy graphs. Each additional vertex and loop contributes to additional suppression of the scale of neutrino mass, provided that the coupling constants at the vertices are small enough. There is also the prospect of suppression from powers of the ratio of electroweak scale to the new physics scale. At some point, the net suppression should become so strong that the 0.06 eV lower bound on the neutrino mass scale cannot be generated for phenomenologically-acceptable exotic particle masses. Indeed, it has been argued that models constructed from opening up mass dimension 13 and higher operators are unlikely to be viable. Our findings in this paper, arrived at by a detailed examination of a specific model constructed from a dimension-11 operator, cast doubt on this tentative conclusion.
The basic reason is evident from Equation 2.2: While there is a product of a few coupling constants in the numerator, and there is the (1/16π 2 ) 2 two-loop suppression factor, the mass suppression is only v 2 /Λ, so identical with that of the usual seesaw models. Formally, this is due to the neutrino mass diagram generating the same dimension-5 Weinberg operator as underpins the seesaw models, with the main difference being that it is generated at loop-rather than tree-level. Additional insertions of v are often produced in radiative models from the need to use quark or charged-lepton mass insertions, but there is no such necessity in this model: the dominantly induced Weinberg-type operator is the standard one at dimension-5, rather than a higher-dimension generalisation obtained by multiplying by powers of H † H. At the level of the underlying renormalisable theory, one also observes that one of the contributing vertices in the numerator is a trilinear scalar coupling, thus having the dimension of mass, and most naturally set at the scale of the new physics. One possible source of suppression is thus absent. With order one dimensionless couplings, and the scalar trilinear coupling set at the new physics scale, the masses of the exotic scalars can be pushed as high as 10 7 TeV. With the dimensionless couplings at 0.01, so of fine-structure constant magnitude, the scale of new physics drops to the hundreds of GeV level. Hence, while the existence of new particles at the 1-100 TeV scale explorable at current and proposed colliders is consistent with this model, it is not inevitable. From this perspective, models which require some of the couplings in the neutrino mass diagram to be standard-model Yukawa couplings (other than that of the top quark) are more experimentally relevant. From the model-building perspective, however, our analysis raises the prospect that some models based on tree-level UV completions of mass dimension 13 (and possibly even higher) operators might be viable.
The specific model analysed in this paper, consisting of an isotriplet scalar leptoquark, an isosinglet diquark and a third exotic scalar multiplet that has no Yukawa interactions, successfully generates neutrino masses and mixings at two-loop level consistent with experimental bounds from a variety of processes, of which µ → e conversion on nuclei proved to be the most stringent. It can also ameliorate the discrepancies between measurements and standard-model predictions for R K ( * ) and the anomalous magnetic moment of the muon.

A The full Lagrangian
The most general SU (3) c ⊗ SU (2) L ⊗ U (1) Y gauge-invariant, renormalisable Lagrangian will consist of the SM Lagrangian, together with additional NP terms, that involve interactions between the following: the exotic scalars, denoted S ( ) , and the SM fermions, F ; the exotic scalars and the Higgs boson, H; and the exotic scalars themselves. It also includes the gauge interactions of the exotic scalars. The accidental lepton and baryon number symmetries which are conserved by the SM Lagrangian have not been imposed here.

A.1 The gauge sector
The gauge covariant derivative acting on the scalar fields is where Y is the hypercharge of the scalar, I k are the SU (2) representation matrices (I k = 0 for SU(2) singlets and −i k for SU (2) triplets), λ A for A = 1, ..., 8 are the Gell-Mann matrices, g 1 , g 2 and g 3 are the respective coupling constants and after electroweak symmetry breaking we have e = g 2 sin θ W > 0 and g 1 /g 2 = tan θ W [54]. with the appropriate selections of hypercharge and SU (2) representation matrices.

A.2 The fermion sector
The part of the Lagrangian involving couplings between the fermions and the scalar bosons is There are no Yukawa interactions allowed by the SM gauge symmetries between the scalar φ 3 and SM fermions. The τ k , k = 1, 2, 3, are the Pauli matrices; i, j = 1, 2, 3 are generation indices; a, b = 1, 2 are SU (2) flavour indices; ab = (iτ 2 ) ab ; and S k 3 are components of S 3 in SU (2) space. The Levi-Civita tensor is needed in order to conserve charge for SU (2) triplets and doublets. Superscript C stands for the charge conjugation operation. For a fermion field ψ: ψ R,L = P R,L ψ, ψ = ψ † γ 0 , and ψ C = Cψ T , where P R,L = (1 ± γ 5 )/2 and C = iγ 2 γ 0 . The diquark coupling to S 1 , z LL 1 is symmetric due to a combination of the antisymmetry of SU (2) structure and of the colour structure of the fermion bilinear, while the diquark coupling to S 3 , z LL 3 , must be antisymmetric for similar reasons. In the expansion of the SU (2) structure of the Lagrangian, we have also rotated into the mass eigenbasis of the quarks, using the convention that u i L → (V † CKM ) ij u j L and d i L → d j L .

A.3 The scalar sector
The scalar part of the interaction Lagrangian will contain quartic interactions with dimensionless couplings, and cubic and quadratic interactions with dimensionful couplings. The full scalar Lagrangian, without baryon number conservation imposed, is x ± = 1 2 (−1 + s + t ± w) w = 1 + s 2 t 2 − 2(s + t + st). (B.7) We need not worry about the non-finite parts of the integrals as we are guaranteed that they will cancel.
Finally, we give the evaluation of I kl in the limiting case where the ratios of SM fermion squared masses on m 2 S 1 go to zero, specifically s k , s l → 0. In this limit lim s k →0,s l →0 (I klαβ ) = I 00αβ = π 4 2 ln t α ln t β +ĝ (t α , t β )(1 + t α + t β ) −ĝ(t α , 0)(1 + t α ) t α t β , where taking the limitĝ(s, t)  . In a similar fashion to that detailed for Model 1 in Section 2.3 neutrino mass is then obtained from the flavour sum of the self-energy diagrams with the freedom to set external momentum to zero: where y LL 1 is the leptoquark coupling matrix, z LL 3 is the diquark coupling matrix and m S 1 S 3 φ 3 represents the cubic exotic scalar coupling. The integral I klαβ is defined as in Equation B.2, with the k 1 · k 2 being suppressed to reduce clutter. The indices k, l run over the three possible quark flavours and are summed over, while α, β run over the two possible leptoquark mass states obtained from the mixing of S * 1 and φ −1/3 3 . We also define R αβ = − cos θ cos θ − sin θ cos θ − sin θ sin θ , where θ is the mixing angle between S * 1 and φ −1/3 3 .
The argument showing that Equation C.1 leads to vanishing neutrino mass goes as follows : 1. The diquark coupling to S 3 , z LL 3 is antisymmetric in flavour space, i.e. z LL 3ij = −z LL 3ji .
2. The integral I klαβ is symmetric under simultaneous relabelling of k ↔ l and α ↔ β, and the matrix R αβ is symmetric under α ↔ β.

The neutrino mass then evaluates to
where in the last line we used points 1. and 2., specifically that z LL 3kl = −z LL 3lk and I klαβ = I lkβα . The neutrino mass arising from the UV-completion of O 47 with S 3 coupling as a diquark, S 1 coupling as a leptoquark and φ 3 not coupling to any SM fermions vanishes.