Reconciling $B$-decay anomalies with neutrino masses, dark matter and constraints from flavour violation

Motivated by an explanation of the $R_{K^{(*)}}$ anomalies, we propose a Standard Model extension via two scalar SU(2)$_L$ triplet leptoquarks and three generations of triplet Majorana fermions. The gauge group is reinforced by a $Z_2$ symmetry, ensuring the stability of the lightest $Z_2$-odd particle, which is a potentially viable dark matter candidate. Neutrino mass generation occurs radiatively (at the three-loop level), and leads to important constraints on the leptoquark couplings to leptons. We consider very generic textures for the flavour structure of the $h_1$ leptoquark Yukawa couplings, identifying classes of textures which succeed in saturating the $R_{K^{(*)}}$ anomalies. We subsequently carry a comprehensive analysis of the model's contributions to numerous high-intensity observables such as meson oscillations and decays, as well as charged lepton flavour violating processes, which put severe constraints on the flavour structure of these leptoquark extensions. Our findings suggest that the most constraining observables are $K^+ \to \pi^+ \nu \bar \nu$ decays, and charged lepton flavour violating $\mu -e$ conversion in nuclei (among others). Nevertheless, for several classes of flavour textures and for wide mass regimes of the new mediators (within collider reach), this Standard Model extension successfully addresses neutrino mass generation, explains the current $R_{K^{(*)}}$ tensions, and offers a viable dark matter candidate.


Introduction
Despite the many successes of the Standard Model (SM) in interpreting experimental data and predicting new phenomena, three independent observations signal the need to consider New Physics (NP) scenarios: neutrino oscillations, the baryon asymmetry of the Universe (BAU), and the lack of a viable dark matter (DM) candidate. Although no direct evidence for NP states has been unveiled at the LHC, certain experimental measurements have revealed non-negligible tensions with respect to the SM predictions. In addition to the anomalous magnetic moment of the muon, recent data hints to several discrepancies with respect to the SM in some B-meson decay modes, potentially suggesting the violation of lepton flavour universality (LFUV). In the SM's original formulation, gauge interactions (both charged and neutral currents) are strictly lepton flavour universal; precise measurements of related electroweak observables -for instance Z → decays [1] have so far been in agreement with the SM's predictions.
The so-called R K ( * ) observables are built from the comparison of the branching ratios (BR) of B into di-muon and di-electron plus K ( * ) final states, and are parametrised as In the above ratio of BRs, the hadronic uncertainties cancel out to a very good approximation, and consequently these observables are sensitive probes of NP contributions [2]. First results on the measurement of R K by LHCb were reported in 2014 [3], having been obtained (as denoted in subscript) for the dilepton invariant mass squared bin q 2 ∈ [1, 6] GeV 2 . This corresponds to a 2.6σ deviation from the SM prediction of R SM K = 1.00±0.01 [4,5]. The corresponding measurements for the decays into K * were reported in 2017 [6], respectively corresponding to 2.3σ and 2.6σ deviations from the expected SM values, R SM K * [0.045,1.1] ∼ 0.92 ± 0.02 and R SM K * [1.1,6] ∼ 1.00 ± 0.01 [4,5]. Interestingly, deviations from SM expectations have also been observed in B → D ( * ) ν decays, in particular in the ratio of tau to final states composed of light leptons, The measured value for R D = 0.407 ± 0.039 ± 0.024 [7], reported by several experiments [8][9][10][11][12][13][14], already deviates from the SM prediction R SM D = 0.299±0.003 [15] by about 2.3 σ. The experimental value of R D * = 0.304 ± 0.013 ± 0.007 [7,13,14,16] is also larger than the SM expectation (R SM D * = 0.260 ± 0.008 [17]), exhibiting a 2/6σ deviation. When combined, the latter experimental results point towards a deviation of 4.1 σ from the SM prediction [7, 16,18,19]. Other anomalies in B meson decays have emerged concerning the angular observable P 5 in B → K * + − processes. While LHCb's results for P 5 in B → K * µ + µ − decays manifest a slight discrepancy with respect to the SM (either due to NP contributions, or possibly a result of SM QCD effects [20]), the Belle Collaboration [21] reported that when compared to the muon case, P 5 results for electrons show a better agreement with respect to the SM prediction. The P 5 results could thus be interpreted as suggestive of the fact that NP effects may be dominant for the second generation of leptons.
The leptoquark hypothesis -both in its scalar and vector field realisations -has been extensively explored in recent years, as have its implications for both quark and flavour dynamics. In particular, having leptoquark couplings which are necessarily flavour non-universal in the lepton sector has multiple implications concerning lepton observables, ranging from the anomalous moment of the muon, to charged lepton flavour violating (cLFV) decays and transitions. Likewise, efforts have been made to connect the neutrino mass generation problem (itself calling upon a modified lepton sector) with an explanation of the flavour tensions via leptoquarks; many such models lead to radiative generation of the light neutrino masses 1 , for instance through mixings with the standard model Higgs boson (in SM extensions via vector leptoquarks [70]), using scalar leptoquarks and color-octet Majorana fermion [55], or calling upon a "coloured Zee-Babu model" [71,72] (with the addition of scalar leptoquarks and diquarks), which leads to two-loop radiative neutrino mass generation [73]. Extensions of SM via both leptoquarks and additional Majorana fermions aimed at connecting B-decay anomalies to neutrino masses (radiatively generated at the three-loop level), and to a solution of the dark matter problem [74,75], relying on the so-called "coloured KNT models" [76]. These studies also evaluated the impact of the BSM construction to cLFV decays, focusing on the rôle of radiative charged lepton decays, → γ.
Building upon the previous analysis, in this work we consider a scalar leptoquark model, which aims at simultaneously explaining the B meson decay anomalies, accounting for neutrino oscillation data and putting forward a viable dark matter candidate. The SM is extended via two scalar leptoquarks h 1,2 and three generations of triplet neutrinos Σ i R . The SM symmetry group is enlarged by a discrete Z 2 symmetry under which only h 2 and Σ i R are odd. While effectively forbidding the realisation of a tree-level type III seesaw [77], the Z 2 symmetry is instrumental to ensure the stability of the lightest Z 2 -odd particle (the neutral component of the lightest triplet), which is found to be a viable dark matter candidate. Neutrino masses can be radiatively generated and, as we argue here, complying with oscillation data turns out to severely constrain the leptoquark Yukawa couplings. Focusing on saturating the R K ( * ) anomalies, we carry a thorough analysis of this phenomenological model. Our study relies in assuming generic perturbative textures for the leptoquark Yukawa couplings: in particular, and contrary to previous analyses, we do not forbid couplings of the leptoquarks to the first generation of quark and leptons. Moreover we take into account a comprehensive set of flavour observables (meson and lepton rare decays and transitions); this allows to identify several classes of textures for the leptoquark Yukawa couplings in agreement with observation. Our findings suggest that the most severe constraints arise from K → πνν decays and neutrinoless µ−e conversion in nuclei -and not from radiative muon decays, as suggested by other studies; furthermore, the joint interplay of these high-intensity observables also disfavours several ansätze for the leptoquark textures previously considered (see [74,75]).
The paper is organised as follows: after presenting the building blocks of the model in Section 2, and discussing neutrino mass generation in Section 3, Section 4 is devoted to establishing first constraints on the model from the requirement of having a viable DM candidate. The B meson anomalies are presented in Section 5, and the discrepancy between SM prediction and observation is parametrised in terms of the leptoquark couplings. Sections 6 and 7 are dedicated to the constraints potentially arising from rare meson processes and from cLFV decays. Finally, our results (both in what concerns identifying viable textures for the leptoquark Yukawa couplings, as well as numerical studies of the model's parameter space) are collected in Section 8. We summarise the most important points, as well as our final remarks, in the Conclusions.
As highlighted in the Introduction, the primary goals of this model are to simultaneously address the problem of neutrino mass generation, and provide a viable DM candidate, while explaining the observed anomalies in B meson decays. If sufficiently long-lived or stable, the neutral component of the lightest triplet Σ 1 R gives a potential cold dark matter (CDM) candidate: the quantum corrections generate a mass splitting such that the neutral component is indeed the lightest one; its stability can be ensured by reinforcing the SM gauge group by a discrete Z 2 symmetry under which both h 2 and Σ i R are odd, while all the SM fields and h 1 are even. Since -and as mentioned before -Σ 1 R is the lightest state, the final DM relic abundance is solely governed by the relevant electroweak (EW) gauge interactions and m Σ 1 R , independent of its Yukawa interactions. This is in contrast with scenarios in which a SU(2) L singlet fermion is considered as a dark matter candidate, subject only to Yukawa interactions.
It is important to notice that a consequence of the Z 2 symmetry is that it forbids a conventional type III seesaw mechanism; however, neutrino masses can still be radiatively generated at higher orders (as discussed in the following section), from diagrams involving the new exotic states and down-type quarks.
The complete particle spectrum is presented in Table 1. The Lagrangian of the present SM extension can be cast as in which L h,Σ int and V H,h scalar respectively denote the interactions of h 1 , h 2 and Σ i R with matter 2 , and the scalar potential, while L Σ mass encodes the Majorana mass term for the fermion triplets. The new interaction and Majorana mass terms are given by In the above i, j = 1 . . . 3 denote generation indices, while a, b = 1, 2 are SU(2) indices; τ c are the Pauli matrices (c = 1, 2, 3), and we have further defined ab = (iτ 2 ) ab . Finally, C denotes charge conjugation. The scalar potential (including SM terms) can be written as As can be inferred from inspection of Eq. (6), the simultaneous presence of the first two terms violates baryon number, and can thus lead to B−L conserving dimension-6 contributions to proton decay. In the following analysis, we will assume that these interactions are absent (i.e., z ij = 0), an hypothesis that can naturally arise from the embedding of the model into an ultraviolet (UV) complete framework, as discussed in [79,80]. The absence of the diquark couplings then allows to unambiguously assign baryon and lepton number to the scalar leptoquarks h 1,2 .
To cast the interaction Lagrangian in a more explicit way, it is convenient to work in the U(1) em basis: the physical states, respectively with electric charges 4/3, −2/3, 1/3 can be written in terms of the SU(2) components as follows Using the above redefinitions, the interaction Lagrangian of Eq. (6) can be rewritten as Following the decomposition of Eq. (9), the most relevant interaction term for neutrino mass diagrams can now be written as In what follows (and for simplicity), we will further assume the couplings λ Hh 1,2 to be negligible; this leads to having degenerate physical masses for the components of the scalar triplets h 1,2 (m 2 h 1(2) = µ 2 h 1(2) + λ Hh 1(2) v 2 /2, with v the SM Higgs vacuum expectation value), and thus allows to comply with EW precision constraints on oblique parameters.

Radiative neutrino mass generation and leptonic mixings
As mentioned in the previous section, the Z 2 symmetry precludes any coupling between the fermion triplets and neutral leptons, which effectively dismisses a type III seesaw explanation of neutrino mass generation (at the tree level). The absence of right-handed neutrinos further excludes the possibility of Dirac-type masses for the neutral leptons. Interestingly, the particle content of the model does allow a natural explanation to the smallness of neutrino masses: these are radiatively generated, from higher order contributions.
The first non-vanishing contributions to m ν arise at the three-loop level, and are induced by the diagrams displayed in Fig. 1, from the exchange of leptoquarks h 1,2 and neutral (charged) fermion triplets, calling upon chirality flips in the internal down-type quark lines (proportional to the down quark masses). Despite the different particle content, the diagrams are akin to those originally proposed in [76], which were mediated via colourless scalars and a Z 2 -odd singlet Majorana fermion. (In the latter case, the computation of the contributions to m ν has been carried in [81].) The computation of the different diagrams of Fig. 1 leads to the following contributions to neutrino masses (which are for simplicity cast in the weak interaction basis): in which Q 1,2 denotes the momenta of the internal down-type quarks, m D = diag(m d , m s , m b ) is the diagonal down-quark mass matrix 3 , and P L,R = (1 ∓ γ 5 )/2. For clarity, we have explicitly highlighted the contributions of the four-scalar leptoquark vertices of diagrams (a) and (b) (as well as the Σ − counterpart of (b)), writing them as products of symmetry and colour factors, respectively s and κ c . Diagram (a) is associated with s a = 4, while (b)-like diagrams lead to s b = 2; for the colour factor one has κ c = 15, due to the possible distinct contractions of the four-scalar leptoquark vertex (in agreement with [75]). Using the appropriate loop integrals 4 , and after a Wick rotation we obtain where m D is again the diagonal down type mass matrix (in the computation of the loop integrals, the down-quark masses are neglected when compared to the heavier h 1,2 and Σ masses in the loop), and we recall that y (ỹ) denotes the Yukawa couplings of the Z 2 -even leptoquark h 1 to matter (Z 2 -odd leptoquark h 2 and triplet fermion Σ R to down quarks); finally G(a, b) is defined as and in the case in which (for simplicity and without loss of generality) m Σ is assumed to be diagonal, G(m 2 Σ j /m 2 h 2 , m 2 h 1 /m 2 h 2 ) will also be diagonal. The neutrino mass eigenstates can be obtained using the transformation 3 For convenience, and without loss of generality, we chose to work in a basis in which the down-quark Yukawa couplings are taken to be diagonal (Y d ij δij = mD i /v), while parametrising the up-quark Yukawa couplings as Y u ij = V † ij mU j /v, with V the Cabibbo-Kobayashi-Maskawa (CKM) mixing matrix. 4 In particular, one makes use of the identity where B0 is the Passarino-Veltman function defined (in terms of the renormalisation scale µ) as [82] B0 k 2 , m 2 1 , m 2 where U is the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) unitary mixing matrix, which we parametrise as follows (18) In the above, s ij ≡ sin(θ ij ) and c ij ≡ cos(θ ij ); δ D is the CP violating Dirac phase while α and β are Majorana phases.
It is important to notice that as can be seen from Eqs. (15) and (17), neutrino masses (and leptonic mixings) do indeed depend on both of the Yukawa couplings involving the leptoquarks, y andỹ. In particular, the former will be at the source of a number of flavour transitions, including rare meson decays, neutral meson-antimeson oscillations as well as charged lepton flavour violating processes. This implies that a strong connection between neutrino phenomena and flavour nonuniversal processes is established via the flavour structure of the Yukawa matrix y.
In the absence of a complete framework proving a full theory of flavour (which would suggest a structure for y andỹ), one can nevertheless parametrise one of the Yukawa couplings -for example, y -using a modified Casas-Ibarra parametrisation [83], which further allows to accommodate neutrino oscillation data.
In order to construct the modified Casas-Ibarra parametrisation, we first notice that from Eqs. (15) and (17) one can write the diagonal neutrino mass matrix as in which we have omitted flavour (generation) indices for simplicity, and where As noted before, for a diagonal m Σ , G and thus F are also diagonal matrices in generation space, which allows to write the identity with R an arbitrary complex orthogonal matrix (R T R = 1) which can be parametrised in terms of three complex angles, θ i (i = 1 − 3). Finally, from Eq. (21) one obtains the modified Casas-Ibarra parametrisation, which allows to write the Yukawa couplings of the h 2 leptoquark in terms of observable quantities (light neutrino masses, leptonic mixings, down quark masses, triplet and leptoquark masses), and of two unknown quantities -the h 1 leptoquark Yukawa couplings and a complex orthogonal matrix -asỹ The above parametrisation (which is in agreement to a similar approach carried in [75]), allows to write the couplingsỹ in terms of y up to a complex orthogonal matrix. As will be discussed in detail in Section 8.1, once the approximate flavour texture of y is inferred from various experimental constraints, that ofỹ can be derived (up to the mixings due to R).

A viable dark matter candidate
Reinforcing the SM gauge group via a discrete Z 2 symmetry ensures the stability of the lightest state which is odd under Z 2 . If the lightest Z 2 -odd particle (LZoP) is neutral, and the strength of its interactions is such that its relic abundance is in agreement with observational data, then it can be indeed a viable dark matter candidate. In our analysis we assume that the spectrum of the new states is such that one has m Σ 1 R < m h 2 . At the tree level, all the components of a generation Σ i R have the same mass m Σ i,tree = m Σ i,± = m Σ i,0 . The degeneracy between the components is broken by EW radiative corrections, which render the charged states heavier than the neutral one. Dropping for simplicity the generation indices (in this section our discussion is focused on the components of the lighest triplet, Σ 1 R ), the splitting between the neutral and charged components is given by [84] and in which α 2 = α e / sin 2 θ w , with α e the fine structure constant and θ w the weak mixing angle.
In the limit m Σ M W,Z (justified by negative searches at LHC and EW precision measurements), the EW radiative corrections are found to be of order m Σ ± − m Σ 0 ∼ 166 MeV. While the latter mass difference is enough to ensure that the neutral component of the lightest triplet, Σ 0 , is indeed the LZoP, it is sufficiently small to be neglected in the subsequent (numerical) analysis.
The relic abundance of the LZoP Σ 0 is determined by its interactions and by the annihilation and coannihilation channels open in view of the particle mass spectrum. In addition to the Yukawa interactions with h 2 , the Σ R triplets are subject to SU(2) L gauge interactions. Since, as previously highlighted, we assume that m Σ 1 R < m h 2 , the relic density of Σ 0 is, to first order approximation 5 , solely determined by its gauge interactions, which govern the distinct annihilation and coannihilation channels (involving also the charged components Σ ± ). In particular, the annihilation and coannihilation processes involve the following channels: Other processes involving the charged components of the triplet fermion must be also taken into account in the Boltzmann equations leading to the computation of the relic density. These include The relevant cross sections for the above mentioned processes are given by [85] where only the coefficients a ij are kept in the relative velocity (v) expansion of the cross section, i.e. σ ij |v| = a ij + b ijv 2 . The computation of the relic abundance follows closely the method of [86], in which the freezeout temperature of the LZoP Σ 0 (x f ≡ m Σ /T f ) is obtained in terms of the thermally averaged effective cross section σ eff |v| . The relevant channels above referred to contribute to the thermally averaged effective cross section σ eff |v| as and the freeze-out temperature is then recursively given by in which g * ∼ 106.75 is the total number of effective relativistic degrees of freedom the freeze-out, and M Pl = 1.22 × 10 19 GeV is the Planck scale. The quantity g eff is related to the degrees of freedom of the triplet components, g 0 = 2 and g ± = 2, respectively for Σ 0 and Σ ± and to ∆ m Σ (the mass splitting between the charged and neutral components of Σ R ), and can be written as Finally, the relic abundance is given by with I a the annihilation integral, which is defined as In the above, one has used the approximation a eff ∼ σ eff |v| (we recall that we have not taken into account the second and higher order terms in the relative velocity expansion of the effective cross section). For our purpose of constraining the parameter space of the model, we will rely on a simple iterative solution of Eq. (29), which allows to infer first limits on the values of m Σ 1 R leading to a relic density in agreement with the most recent observational data [87], The results of this (approximative) numerical analysis are displayed in Fig. 2, which reveals that having an LZoP mass in the range 2.425 TeV < m Σ < 2.465 TeV leads to a dark matter relic abundance in agreement with the latest data. In what follows, we will use m Σ ∼ 2.45TeV as an illustrative benchmark value for the lightest fermion triplet mass.
Although already mentioned, we nevertheless stress again that several approximations were done in our computation of the relic abundance; moreover, one should also explicitly solve the Figure 2: Estimated relic abundance (see text for details) of the dark matter candidate, Σ 0 , as a function of its mass. The rose-coloured band corresponds to the latest observational data Ωh 2 = 0.1198 ± 0.0015 [87].
coupled Boltzmann equations to numerically obtain the allowed mass range of the LZoP (the lightest Σ 0 ). Other than complying with the dark matter relic abundance, the potential candidate is also subject to the increasingly strong constraints from direct and indirect searches (see, e.g. [84,88]); a detailed discussion of the latter (and the dedicated facilities) lies beyond the scope of this work.

Addressing the B B B meson decay anomalies
As can be seen from the interaction Lagrangian of Eq. (10), the scalar leptoquark h 1 couples to both down-type quarks and charged leptons, the couplings having a priori a non-trivial structure in flavour space. This will open the door to new contributions to numerous rare flavour changing transitions and decays, and -most importantly -can potentially lead to lepton flavour nonuniversal effects, such as those currently suggested by the reported LHCb anomalies.
In this section we will explore to which extent the current model succeeds in addressing the B meson decay anomalies, R K ( * ) and R D ( * ) , respectively associated with neutral and charged current transitions. The abundant constraints arising from negative searches for NP effects in meson oscillation and decays, as well as in charged lepton flavour violating observables will be discussed in Sections 6 and 7.
5.1 Neutral current anomalies: R K and R K * As mentioned in the Introduction, recent measurements [3,6] of the ratios of branching ratios of B → K(K * ) decays into pairs of muons over those into di-electrons exhibit non-negligible deviations when compared to the SM predictions [4,5]; as already stated one has where the dilepton invariant mass squared bin (in GeV 2 ) is identified by the subscript. The comparison of SM predictions with observation respectively reveals deviations of 2.6σ, 2.4σ and 2.5σ. For the leptoquark mass regime considered here (multi-TeV), the neutral current effects induced by the heavy degrees of freedom (SM and NP contributions) in the quark level transitions d j → d i + − can be described by the following effective Hamiltonian [89,90] in which G F is the Fermi constant and V is the CKM mixing matrix. The effective operators present in the above equation can be defined as (for simplicity, hereafter we drop the ij superscripts, which are set to i, j = s, b for the process b → s + − ): In the above, e is the electric charge and σ µν = i[γ µ , γ ν ]/2. The set of primed operators (X = 7 , 9 , 10 , S , P ) comprises those of opposite chirality, and can be obtained by replacing P L ↔ P R in the quark currents. The contribution of the right-handed current operators is negligible in the SM. Flavour universality of lepton-gauge interactions in the SM implies that the Wilson coefficients of operators O i are universal for all lepton flavours ( = e, µ, τ ), and the strict conservation of individual lepton flavour further precludes cLFV Wilson coefficients C i ( = ).
In the present scalar leptoquark model, once the heavy degrees of freedom have been integrated out (under the assumption that M 2 t,W,Z m 2 h 1 ), the NP effective Hamiltonian is given by [91] H comparing the above with the operator basis of Eq. (33), it is possible to infer the following contributions to the Wilson coefficients The deviations from the SM lepton flavour universality imply that the modifications to the Wilson coefficients are necessarily non-universal for the muon and electron entries; the modelindependent fit of [91] suggests the following corrections at the 1σ range: (notice that the second constraint is compatible with zero, and can be fulfilled by setting C µµ,ee 9 ,10 = 0). In the present NP construction, leptoquark couplings to both muons and electrons are present, and are of left-handed nature. Given that C 9 = −C 10 (cf. Eq. (36)), the best fit to the C µµ,ee 9,10 NP Wilson coefficients of Eq. (36) can be recast as Global fits to a large number of observables probing lepton flavour universality in relation to b → sµ ± µ ∓ , b → se ± e ∓ and b → sγ processes also suggest NP scenarios consistent with the above fit to the LFUV R K ( * ) observables. A common conclusion that can be generically drawn is that the NP responsible for the observed discrepancies in b → s data appears to predominantly couple to muons, and is strongly manifest in vector operators (as O µµ 9,10 ). Recent studies and fits by a number of authors (see, e.g. [5,[92][93][94][95]) advocate NP contribution to C µµ 9 (∼ −1) only, or then SU(2) L invariant scenarios (C µµ 9 = −C µµ 10 ∼ −0.6) as preferred NP solutions to alleviate the tensions with the SM.
In terms of the h 1 leptoquark mass and couplings, the expression of Eq. (39) translates into the following condition [91] In Fig. 3, for an illustrative value of m h 1 = 1.5 TeV, we display the (y b y s ) parameter space consistent with the observed values of R K and R K * at the 1σ level, in agreement with Eq. (40). The inset plot shows the regimes of y bµ and y sµ compatible with R K ( * ) (also for m h 1 = 1.5 TeV, and for fixed values of the couplings to the electron, y be y se = 2 × 10 −5 ).

Anomalies in
Several experimental collaborations have also reported deviations from lepton flavour universality in association with B → D * ν decays (charged current b → c ν i transitions). A scalar charged leptoquark with non-trivial couplings to quarks and leptons can mediate d k → u j ν i transitions at the tree level, via the exchange of a charged h 1/3 1 . The SM effective Hamiltonian governing these transitions is thus modified as follows: Figure 3: On the main plot, 1σ region in the (y be y se −y bµ y sµ ) parameter space consistent with R K and R K * data, cf. Eq. (40), for an example with m h 1 = 1.5 TeV. On the inset plot we display the R K ( * ) allowed region in the (y bµ − y sµ ) plane, for fixed values of the leptoquark electron Yukawa couplings y be y se = 2 × 10 −5 (again for m h 1 = 1.5 TeV).
The experimentally measured decay probability is an incoherent sum over the (untagged) neutrino flavour i; one thus finds where x j = (V y * ) j v 2 /4V cb m 2 h 1 , and in which one has used the unitarity of the PMNS matrix. The SM width for the decay B → D * ν (i.e., for d k = b, u j = c) will thus be corrected by an overall factor Pure leptoquark contributions are suppressed by an additional v 2 /m 2 h 1 factor with respect to the SM-NP interference term, and can be hence neglected as a first approximation.
Defining R D ( * ) as the ratio of the decay widths of tau and muon modes -that is, (equal to one in the absence of NP). After combining current experimental world averages with the SM predictions, the current anomalous data can be parametrised as in which the statistical and systematical errors have been added in quadrature. Similar ratios comparing distinct final state lepton flavours can be built to test possible LFUV in the corresponding sectors. For example, the Belle Collaboration has reported measurements of the ratios R After having detailed the leptoquark contributions to the meson observables currently exhibiting a significant deviation from the SM expectation, we now address other processes which -being in agreement with SM predictions (negative searches or compatible measurements) can constrain the masses of the new states and their couplings.

Constraints from rare meson decays and oscillations
Various observables involving mesons lead to important bounds on leptoquark couplings; here we discuss the (leptoquark) NP contributions to leptonic and semi-leptonic meson decays (occurring at the tree level), and to meson oscillations and rare radiative decays, both at the loop level. The SM predictions and current experimental bounds for the processes here discussed are summarised in Table 2. The stringent bounds on NP contributions arising from the now observed decay B s → µ + µ − are not discussed here, as they have been implicitly taken into account in defining the allowed ranges for C µµ 10,NP and the y 22,23 Yukawa couplings in the previous section. Likewise, we do not include the constraints arising from semileptonic K-decays into charged dileptons (K → π ), as theoretical (SM) predictions are plagued by important uncertainties, and are not yet up to par with the precision of the experimental results (see for example [105,106]). 6.1 Rare K K K and B B B meson decays: The s → dνν and b → sνν transitions provide some of the most important constraints on NP scenarios aiming at addressing the anomalies in R K ( * ) and R D ( * ) data.
Following the convention of [107], at the quark level the |∆S| = 1 rare decays K + (K L ) → π + (π 0 ) ν ν and B → K ( * ) ν ν can be described by the following short-distance effective Hamiltonian for d j → d i ν ν transitions [108,109]   in which i, j denote the down-type flavour content of the final and initial state meson, respectively. In the SM, only the lepton flavour conserving left-handed operator is present; the associated Wilson coefficient C L,ij corresponding to a d j → d i transition is given by In the above, X SM t and X c ( = e, µ, τ ) are the loop functions associated with the (lepton flavour conserving) contributions from the top and charm quarks, respectively. In the SM, computations of the top loop function (including Next-to-Leading Order QCD corrections [110][111][112] and the full two-loop electroweak corrections [113]) have led to the result X SM t = 1.481(9) [98]; the charm loop functions X c have been computed at NLO [112,114] and at Next-to-NLO [115,116], and their numerical value is given by 1 3 X c /|V us | 4 = 0.365 ± 0.012. Once the heavy leptoquark degrees of freedom have been integrated out, the general effective Hamiltonian describing the NP contribution to d j → d i ν ν processes induced by h 1 is given by and comparison with Eq. (47) leads to the following NP Wilson coefficient, in which we notice the presence of the PMNS matrix U , and the possibility that the tree level NP contribution, via the exchange of h 1/3 1 , may lead to different lepton flavours in the final state, i.e. = .
The branching fraction for K + → π + ν ν (corresponding to setting i = d, j = s in the previous discussion) is given by [108,109] in which κ + = 5.173(25) × 10 −11 (|V us |/0.225) 8 , ∆ em = −0.003 is the electromagnetic correction from photon exchanges and ∆ Pc,u = 0.04(2) denotes the long-distance contribution from light quark loops, computed in [117]. The loop function associated with NP exchanges, X t , is given by Notice from both Eqs. (51,52) that the SM top and charm contributions are lepton flavour conserving. Experimentally, the measurement of the charged kaon decay mode [99], is expected to be improved in the near future by the results of the NA62 collaboration. The NA62 recent measurement [100] (one event) BR(K + → π + νν) exp = 28 +44 −23 × 10 −11 still suffers from large statistical uncertainties, but these are expected to improve by an order of magnitude in coming months. The experimental bounds will allow to better constrain X t , and hence the leptoquark contributions, encoded in C , LQ The branching fraction for the neutral mode K L → π 0 ν ν can be written as [108,109] with κ L = 2.231(13) × 10 −10 (|V us |/0.225) 8 in the framework of the SM. Concerning the experimental status, only a 90% C.L. bound has been reported for the decay [101] BR which -as occuring for the charged decay modes -can also be used to constrain X t and C , LQ L,ij (albeit leading to weaker constraints than those inferred from the K + → π + νν mode).
The latest experimental data from Belle [102] allows to infer the following bound, R νν K ( * ) < 3.9(2.7) at 90% C.L., which can be used to constrain combinations of leptoquark couplings (lepton flavour conserving and violating) as given in Eq. (57).

Neutral meson mixings and oscillations
where [57] with {i, j} respectively denoting {b, s}, {b, d} or {d, s} for P = B 0 s , B 0 d or K 0 mesons. Let us begin by discussing B 0 s mixing, which is potentially the process most sensitive to the couplings y b and y s which are at the origin of the R K ( * ) anomalies. Following [103,118,119], one can define the ratio of the total contribution (NP and SM) to the SM one In the leptoquark model we consider, the relative NP contribution p s can be cast as  [119]. Both analyses obtain their tightest 1σ constraints for imaginary NP contributions (i.e., arg p s = −π/2): |p s | < 0.016. In our study, and for simplicity, we will assume real Yukawa couplings y q and a real p s . In this case one has p s ≥ 0, and both analyses lead to p s < p max s ≈ 0.17, which translates into the following 1σ upper bound for the (real) h 1 leptoquark couplings: Following the same global analyses, a similar bound can be inferred from data on B 0 d mixing (still in the case of real couplings), which leads to p d < p max d ≈ 0.13. One thus has For the case of K 0 −K 0 mixing, the SM effective coupling can be expressed as where F (x c , x t ) denotes the contribution of the distinct Inami-Lim functions, and is defined as 120,121]; the coefficients η cc = 1.87 (76), η tt = 0.5765(65) and η ct = 0.496(47) encode NNLO QCD corrections [104,[122][123][124]. The last two terms in Eq. (65) can be safely neglected (as they are CKM-suppressed); the first (and dominant) term is associated to important theoretical uncertainties, O(40%). From Eq. (59), the leptoquark contribution can be written as in whichr ≈ 0.95 allows to take into account the difference between the relevant scales (M t and m h 1 ) [125]. Real couplings y q only affect ∆m K = Re K 0 |H K |K 0 /m K . Taking into account only the (dominant) first term of Eq. (65), the SM short distance contribution (cf. Eq. (64)) gives [104,126] (∆m K ) SM = (3.1 ± 1.2) × 10 −15 GeV = (0.89 ± 0.34) (∆m K ) exp Comparing (∆m K ) exp with the theoretical prediction ((∆m K ) SM + (∆m K ) NP + (∆m K ) LD ), thus allows to obtain a conservative upper bound, p K < p max K ≈ 0.45/0.55 = 0.81. Leading to the latter limit, one has taken the lowest values of both (∆m K ) SM (∼ 55%) and the long distance (LD) contributions (∆m K ) LD , which are hard to evaluate but are expected to be positive like p k [126]. This translates into an upper bound on the leptoquark couplings:

Leptonic decays of pseudoscalar mesons
The leptonic decays of pseudoscalar mesons are known to provide stringent constraints on models of NP with modified lepton (and/or quark) sectors; leptoquark models are no exception, and in what follows we discuss some relevant modes of K and B mesons.
Leptonic B s meson decays B 0 (s) → ± ∓ are well predicted in the SM (the only hadronic uncertainty coming from the decay constant f Bs ). At present, only the B s → µ + µ − decay mode has been observed, and it is in agreement with the SM. The LHC collaborations have reported values for its branching fraction of (2.8 +0.7 −0.6 ) × 10 −6 [127,128]. As mentioned before, these bounds have been taken into account upon saturation of the B decay anomalies.
For both K and B mesons, the cLFV leptonic decays have been shown to lead to important constraints on NP models: the cLFV B decay modes are particularly relevant for leptoquark SM extensions [129]; although the hard to quantify long-distance QCD corrections render non-trivial an estimation of the leptonic K L decays [130]), the upper bounds on the cLFV mode K L → µ ± e ∓ prove to be one of the most stringent constraints on the couplings of leptoquarks to the first two generations of leptons and down-type quarks.
Following the effective Hamiltonian conventions adopted in Eq. (33), the decay width of P → ± ∓ is governed by the O ij; 9 and O ij; 10 operators. In terms of the corresponding Wilson coefficients C ij; 9,10 , the decay width can be written [129] in which m i,j denotes the mass of the meson valence quarks, the index q refers to up-type quarks (the sum being dominated by the top quark contribution), and one defines

Charged lepton flavour violating processes
Due to the presence of new states with non-negligible couplings to neutral and charged leptons, which are a source of LFUV, one expects that the model under consideration will give rise to important contributions to cLFV observables. While most cLFV decays correspond to higher order (loop) processes, it is important to notice that transitions occurring in the presence of matter, such as neutrinoless muon-electron conversion in nuclei can now occur at the tree level. In the following, we address the contributions of the model to several cLFV observables 6 , whose experimental status (current bounds and future sensitivities) is summarised in Table 3. The bounds on leptonic observables will give rise to stringent constraints on the parameter space of the model: other than neutrino oscillation data, bounds on several lepton flavour violating observables will play a significant rôle in identifying the regimes which can successfully lead to an explanation of the B meson anomalies. 6 In what concerns three-body decays we focus our discussion on the case of same-flavour final lepton state, i.e., → 3 .
7 For a relevant discussion concerning the future perspectives of the searches for µ → eγ process see [133].

Radiative decays → γ
In the present model, radiative charged lepton decays are induced at the loop level (h 1 leptoquarks and quarks running in the loop). The effective Lagrangian for → γ decays can be written as [144] with F µν the electromagnetic field strength. The flavour violating coefficients can be cast as the new states and interactions give rise to the following effective coefficients σ L(R) , In the above, X q ≡ − √ 2y q and Q S = 4/3 for down-type quarks (q = d), while X q ≡ −V * qq y q and Q S = 1/3 for up-type quarks (q = u). The loop functions, cast in terms of x q = m 2 q /m 2 h1 , are given by Finally, the cLFV radiative decay width is given by

Three body decays →
The photonic interactions at the source of the radiative decays (parametrised by the couplings C L,R , see Eq. (71)) will also induce the three-body cLFV decays; moreover, direct four-fermion interactions are responsible for additional contributions. Following [145,146], the low-energy effective Lagrangian including the four-fermion (contact) operators responsible for → decays can be written as In the model under study, there are several types of diagrams contributing to the 3-body cLFV decays: photon penguins (dipole and off-shell "anapole"), Z penguins and box diagrams, all due to flavour violating interactions involving the scalar leptoquark h 1 and quarks. Neglecting Higgsmediated exchanges, the distinct diagrams will give rise to non-vanishing contributions to the dipole operators (C L,R ), as well as to g 3 , g 4 , g 5 and g 6 . The box and the Z penguin diagrams contribute to g 4 and g 6 as follows [57] g box,Z in which x = M 2 Z /m 2 h 1 . The off-shell γ-penguin diagrams induce non-vanishing contributions to g 3,5 and g 4,6 , which are given by [145] in which the form factors can be cast as and are defined in such a way thatf E0 (q 2 ) andf M 0 (q 2 ) are finite at q 2 → 0 (q being the fourmomentum transfer). The cLFV loops involving h 1 leptoquarks and up (or down) quarks contribute to the off-shell penguin form factors as follows We recall that in the above X i ≡ − √ 2y q for i = d, and X i ≡ −V * ij y j for i = u; the loop function f γ (r i ) is given by with r i = m 2 i /(−q 2 ), where i denotes the quark in the loop. As an example, for the case of µ → 3e decays, one is led to the following branching ratio [145,146] Analogous expressions can be easily inferred for the other cLFV 3-body decay channels.

Neutrinoless µ-e conversion in Nuclei
One of the most important constraints on SM extensions with scalar leptoquarks 8 arises from the nuclear assisted µ − e conversion. Phenomenologically -and contrary to other cLFV transitions which remain loop-mediated processes -neutrinoless µ − e conversion can occur at the tree-level in the presence of lepton-quark-leptoquark interactions. Moreover, as summarised in Table 3, the experimental prospects for this cLFV observable are particularly promising: not only current bounds (obtained for Gold nuclei) are already O(10 −13 ), but in the near future several dedicated experiments should bring the sensitivity down to O(10 −17,−18 ) [142,143]. The conversion ratio is defined as in which Γ capture (Z) denotes the capture rate for a nucleus with atomic number Z, and Γ(µ − e, N) is the cLFV width, which can be generically cast as follows [144,148]: where the (tree-level) flavour violation is encoded in the following quantities [144] g (d) Other than the (dominant) tree-level exchanges, we have also included the photon-penguin contributions in the expression of the conversion width 9 ; these are associated with the C L(R) coefficients, which have been previously defined in Eq. (72). The relevant nuclear information (nuclear form factors and averages over the atomic electric field) are encoded in the D, V (p) and V (n) form factors. The latter overlap integrals have been numerically estimated for various nuclei [148]; Table 4 Nucleus D[m  µ ) and total capture rates for Gold and Aluminium [148].
summarises some of the above quantities for Gold and Aluminium nuclei (in units of m 5/2 µ ), as well as the corresponding capture widths.
Current bounds (from Gold nuclei) already allow to infer the following stringent constraints [144]: g

Further leptonic observables
Other tensions between SM predictions and observation have also fuelled the need for NP. One such case is the anomalous magnetic moment of the muon; we notice here that the present leptoquark construction does provide a non-vanishing contribution to (g − 2) µ , albeit with the "wrong" sign [57], so that it cannot ease the current discrepancy.
Being a priori complex, the leptoquark couplings can also induce contributions to the electric dipole moments of quarks and leptons, at the two loop level. Although these could possibly allow to further constrain the couplings (in particular, the CP violating phases), such a detailed analysis lies beyond the scope of the current work.

Accommodating B-anomalies, dark matter and neutrino data
In the previous sections we have discussed in detail the distinct observations and experimental tensions that the present leptoquark model is called upon to explain; moreover, we have also addressed a comprehensive set of observables (encompassing numerous quark and lepton flavour transitions) that are expected to lead to important constraints on specific realisations of the model.
In this section, we finally identify the different regimes thus allowing to: (i) accommodate the latest data on neutrino oscillation parameters; (ii) account for a correct relic abundance for the dark matter candidate; (iii) explain the R K ( * ) and R D ( * ) anomalies; (iv) be compatible with all available bounds on leptoquark couplings and masses arising from direct searches, as well as from the relevant leptonic and semi-leptonic meson decays and transitions (including neutral meson oscillations and rare meson decays) -both tree level and higher order processes -, and cLFV processes (radiative and three-body decays, and µ − e conversion in nuclei).
The results of the approximative numerical study of Section 4 suggested that a viable dark matter candidate could be obtained for a LZoP (the lightest Σ 1,0 ) mass in the range 2.425 TeV m Σ 2.465 TeV, as inferred from Fig. 2. As working benchmark values, we will thus set the masses of the three generations 10 of m Σ i as 2.45, 3.5 and 4.5 TeV. By construction, the other Z 2 -odd particle, h 2 , must be necessarily heavier than Σ 1 ; in order to comply with the hierarchy of the Z 2 -odd spectrum, we choose m h 2 ∼ 2.6 TeV. Notice that h 1 is not subject to any DM-related arguments; its mass is not related to that of the LZoP, nor to m h 2 , and can in principle vary in the TeV range. The chosen (illustrative) benchmark values of the scalar leptoquarks and fermion triplets are in agreement with the current limits established by negative collider searches; we refer to [57,75] for a detailed discussion. We nevertheless highlight here a few important points, and current experimental bounds. Both leptoquarks can be pair produced via strong interactions pp → h 1(2) h 1 (2) ; each of the Z 2 -even h 1 can subsequently decay into quark-lepton pairs (either neutrino or charged lepton). Searches for dilepton+dijet signals have been carried by the ATLAS and CMS collaborations: for the 13 TeV run data, and considering decays into ue and cµ, ATLAS has set lower bounds on leptoquark masses of 1100 GeV (900 GeV), respectively assuming 100% (50%) branching fractions [149]; mass limits on leptoquarks decaying to bτ have been established by CMS, which has reported bounds 850 GeV (550 GeV) assuming 100% (50%) branching fractions [150]. The Z 2 -odd h 2 can decay into Σ and a down type quark; the decay mode associated with the neutral component of the triplet leads to a dijet + / E T signal (common to several supersymmetry search channels). Preliminary bounds on the mass of h 2 can thus be inferred from current squark mass limits: about 1.3 TeV for first generation scalar down quarks [151] and 800 GeV for third generation sbottoms [152].
A survey of the previous sections dedicated to neutrino masses and flavour observables reveals that the Yukawa couplings of the leptoquark h 1 to matter are at the core of the distinct observables so far discussed: on the one hand, y is responsible for saturating the B-meson anomalies and accounting for ν-oscillation data (recall thatỹ can be inferred from y using a modified Casas-Ibarra parametrisation, see Eq. (22)); on the other, its different entries are severely constrained by the strong bounds arising either from negative searches or apparent SM-compatibility of a vast array of flavour observables. Our first goal will thus be to identify the most minimal flavour textures that can comply the with points listed above.

Towards a parametrisation of the scalar leptoquark Yukawa couplings
Similarly to what occurs with the quark and lepton Yukawa couplings in the SM, the new couplings are associated with numerous degrees of freedom (being complex, y contains 18 free parameters) which, in the absence of a full theory of flavour, can only be moderately constrained by data.
A possible approach to circumvent the latter problem relies in extending the symmetry group to include flavour symmetries, which can effectively reduce the number of free parameters. To this end, there have been attempts to obtain hierarchical leptoquark patterns by embedding the extended particle content in a Froggatt-Nielsen framework, which can also explain the fermion mass hierarchies as well as the CKM mixing pattern [153]. The Froggatt-Nielsen (FN) mechanism is usually implemented via a U(1) symmetry 11 and a singlet scalar, non-trivially charged under the U(1) FN . The singlet scalar then acquires a vacuum expectation value v FN at some high scale Λ FN , resulting in a suppression of the non-renormalisable Yukawa interactions by a factor (v FN /Λ FN ) n , where n is the sum of the fermion U(1) FN charges [70]. Alternatively, a weakly broken U(2) 5 flavour symmetry has also been proposed in the context of possible interpretations of the B-decay anomalies [158].
A systematic and comprehensive study of the allowed textures for the leptoquark couplings, relying on symmetry-inspired flavour constructions, clearly lies beyond the scope of the analysis; here, we will adopt a phenomenological approach, and identify possible textures for the new Yukawa couplings from the requirements of explaining the B-meson anomalies while complying with all available experimental bounds. As a starting point (inspired by generic FN-like flavour patterns), we consider generic parametrisations of y in terms of powers of a small parameter (taken to be positive and real), with each entry 12 weighed by an O(1) real coefficient a ij : with denoting that there is no summation implied over i, j.
As a first step, we set the individual coefficients a ij = 1, and use the requirement of saturating the R K ( * ) tensions to infer the size of the parameter : in Section 5 we have seen that at the leading order, the explanation of the R K ( * ) anomalies constrains combinations of the quark-"muon" couplings y 22 and y 32 (further depending on inverse powers of the h 1 leptoquark mass). For a benchmark value of m h 1 ∼ 1.5 TeV, one is led to the following relation in which we have elected n 22 + n 32 = 4 as a natural choice (so that ∼ O(1)).
Following the above, and having fixed ≈ 0.215 (notice that this value reflects the choice of m h 1 ), we express the most general texture written in terms of the parameter and positive integers n ij ≥ 1 with i, j = 1, 2, 3: subject to the constraint n 22 + n 32 = 4 to explain the R K ( * ) anomalies. The experimental bounds on rare meson and charged lepton decays can now be used to identify generic textures 13 for y which are compatible with observation. As can be inferred from the analytical expressions presented in Sections 6 and 7, the Yukawa couplings of h 1 to the first two generations of quarks are stringently constrained from rare meson decays; likewise, its couplings to the first generations of leptons are expected to be limited by cLFV transitions in the µ − e sector. A numerical scan of all possible textures (i.e., thorough tests on the viability of each y(n ij ) -cf. Eq. (88), for fixed values of and m h 1 and with n 22 + n 32 = 4) has shown that the most constraining observables turn out to be the rare decay K + → π + νν and, on the lepton sector, µ − e conversion in nuclei and the radiative µ → eγ decay.
The numerical study has further allowed to identify generic classes of representative textures which are in agreement with the B-meson anomalies as well as all leptonic and mesonic processes taken into account (for the above mentioned ( , m h 1 ) benchmark): µ → eγ, τ → eγ, τ → µγ, µ → 3e, τ → 3µ, µ − e conversion, and K + (K L ) → πνν, B → K * νν, B 0 s −B 0 s oscillations 14 , as well as the cLFV decays B s → µe and K L → µe.
For each of the classes identified, we have chosen an illustrative case (setting n ij in agreement with Table 5), and we have evaluated the associated contributions to the different leptonic and mesonic observables mentioned above. The information is summarised in Table 6. We do not include here bounds from the neutral meson-antimeson oscillations since they are considerably less constraining than the meson decay processes involving the same set of leptoquark couplings. Figure 4 graphically summarises the information given in Table 6: for each type of texture, we display the associated predictions for µ → eγ, τ → µγ, τ → eγ, µ → 3e, τ → 3µ, µ − e conversion, K + (K L ) → π + (π 0 )νν, R νν K ( * ) as well as B s → µe and K L → µe. For each process we include the current experimental bounds and future sensitivities, and when applicable, the SM predictions.
Before proceeding, we briefly comment on the other LFUV observables discussed in Section 5.2. The present leptoquark construction leads to SM-like predictions to the distinct R D ( * ) ratios (independently of the texture type and/or mass regime for the exchanged pseudoscalar). If on the one hand this means that, once the distinct experimental constraints have been taken into account, the muon to electron ratios R µ/e D ( * ) are consistent with experimental measurements, on the other hand it also implies that the current experimental measurement of R D ( * ) (tau to muon ratio, exhibiting a significant deviation from SM predictions) cannot be accounted for. Should the latter R D ( * ) discrepancy be confirmed in the future, then the present leptoquark construction will be ruled out, at least in this minimal version.

Constraining the leptoquark parameter space
The previously chosen textures, as well as the numerical results (both for the parameter and for the contributions to the distinct observables) were obtained for a benchmark value of the h 1   Table 5 (viable for the  To explain the R K ( * ) anomalies, the BSM construction must comply with the conditions given in Section 5.1, in particular with the interval for the C ee,µµ 9,NP couplings given in Eq. (39); this can also be written as a condition on the ratio of the relevant Yukawa couplings to the h 1 leptoquark mass, 0.64 × 10 −3 Re[y bµ y * sµ −y be y * se ] (m h 1 /1 TeV) 2 1.12 × 10 −3 . Varying the mass of the h 1 leptoquark over a wide interval -in agreement with LHC direct search bounds -leads to new ranges for the relevant entries of the Yukawa couplings y ij (and thus new values for ). Since for increasing values of m h 1 saturating the R K ( * ) anomalies calls for larger y 22,32 -and hence for larger -, bounds on other observables are expected to become more severe, leading to the exclusion of a given realisation. This is displayed on the distinct panels of Fig. 5: the coloured regions denote the contributions for a given observable arising from varying (i.e., y ij ) in the R K ( * ) favoured interval given above. Light (dark) regions correspond to allowed (excluded) regimes in view of current experimental bounds.
bounds obtained in association with the maximal and minimal values of ).
Despite being manifest in all panels, the "kinks" associated with the contributions to CR(µ−e, Au) are particularly apparent for type III textures 15 . The behaviour has been well identified in the literature (see, e.g. [159]), and stems from having a localised cancellation of opposite sign up-and down-type quark contributions to the conversion rate (due to different charge and weak isospin).
At this point, it is important to recall that the proposed parametrisation for y, as given in Eq. (86), allows each element to be weighed by a real coefficient a ij , O(1). In order to understand how generic perturbations of the unconstrained entries of y affect the phenomenological viability of the model, we have thus taken a type I texture for y, and varied one a ij at a time 16 in the range [0.4, 1.6] (with a step size of 0.1). Although not explicitly displayed here, the numerical studies revealed that µ → eγ and µ → 3e are predominantly sensitive to y 31 (a 31 ), with a mild secondary dependence on y 21 (a 21 ); likewise, τ → µγ and τ → 3µ are controlled by y 33 (a 33 ); τ → eγ and τ → 3e are sensitive to variations from both y 31 and y 33 . Finally, K + → π + νν exhibits a significant dependence on y 13 (a 13 ), y 21 (a 21 ) and y 23 (a 23 ). As a general qualitative statement, for all the radiative and 3-body decays mentioned above, the variation of the a ij coefficients in 10 -7 m h1 (GeV) Figure 5: Contributions to BR(µ → eγ), BR(µ → 3e), CR(µ − e, Au) and BR(K + → π + νν) as a function of the h 1 leptoquark mass m h 1 , for type I, II and III textures of y (respectively from left to right, top to bottom), complying with the current interval for R K ( * ) . Light (solid) surfaces denote currently allowed (excluded) regimes due to the violation of the associated experimental bound.
the interval [0.4, 1.6] leads to a variation of about one order of magnitude in the prediction for the observable. Occurring at the tree level, the neutrinoless µ − e conversion strongly depends on y 11 (a 11 ), y 12 (a 12 ), and y 21 (a 21 ) -the dominant element depending to a certain extent on the leptoquark mass regime. The same study can be carried for type II and III textures, with the results reflecting the relative n ij dependence.

Final constraints from neutrino oscillation data
The requirements of having a viable DM candidate, and of accounting for the R ( * ) K anomalies while complying with all available data on meson and lepton rare decays and transitions have allowed to identify viable flavour textures for the y leptoquark Yukawa couplings, as well as mass regimes for the new states.
As mentioned in Section 3, once the flavour structure of y has been fixed (be it from theoretical arguments or, as in the present case, from a comprehensive phenomenological analysis), the modified Casas-Ibarra parametrisation of Eq. (22) readily allows to determineỹ, while complying with current neutrino oscillation data. As a final step in our study, we thus consider the three textures already discussed in the previous section, and for each one we vary the n ij powers of in agreement with the ranges given in Table 5, as well as the associated a ij prefactors (in the range [0.4,1.6]). The h 1 leptoquark mass is, for simplicity, set to the benchmark value of 1.5 TeV (although the results here discussed qualitatively hold for other choices -in agreement with the discussion of the previous subsection).
Concerning neutrino data, we use the best-fit values from the global oscillation analysis of [160], taking a normal ordering for the neutrino spectrum, with the lightest neutrino mass taken in the range [10 −8 eV, 0.001 eV]. As already mentioned, we take the right-handed triplet masses to be m Σ =2.45, 3.5 and 4.5 TeV. The remaining degrees of freedom in the modified Casas-Ibarra parametrisation are randomly sampled from the following intervals: [0, 2π] for the phases, and [−4π, 4π] for the angles; one further has λ h 4π (see Eq. (8)).
Each of the thus obtained couplings (y andỹ) are again subject to the various flavour constraints previously discussed; moreover, each entry of the couplings must comply with perturbativity requirements, |y(ỹ)| 4π.
In order to illustrate our findings, we display in Fig. 6 the results of the scan, for the three types of textures. Since neutrinoless conversion in nuclei and the K + → π + νν decay are the most constraining observables, we display the corresponding predictions of the randomly sampled textures in the plane spanned by the latter two observables; the colour code distinguishes between perturbative and non-perturbative entries of the y andỹ couplings.
As is manifest from inspection of Fig. 6, accommodating ν-oscillation data from type I textures for the leptoquark y couplings, in agreement with experimental data, and for perturbativeỹ does not excessively constrain the remaining degrees of freedom. Even though perturbative regimes forỹ are more likely to be associated with large values of BR(K + → π + νν), one can easily find regimes which are phenomenologically allowed. Notice however that a near-future improvement in the associated experimental sensitivities (for instance CR(µ − e, Al) ∼ 10 −15 , and BR(K + → π + νν)∼ 10 −10 ) should allow to probe the present leptoquark construction (and possibly falsify it).
For type II textures, the associated panel of Fig. 6 reveals that perturbativeỹ couplings are far harder to accommodate, especially due to the excessive contributions to CR(µ − e, Au). Finally, notice that only a tiny subset of the sampled type III textures is in agreement with flavour observables, and no sub-region of the latter leads to perturbativeỹ. One thus concludes that type III textures (despite being marginally compatible with all the quark and lepton observables here discussed) do not lead to a satisfactory leptoquark construction.
We have also explored the possibility of having a distinct ordering (inverted) for the light neutrino spectrum: in what concerns type I textures, we found no significant changes, so that in fact both orderings can be easily accommodated; in the case of type II textures for y (which do allow to accommodate oscillation data for a normal ordering) we failed to find viable solutions for an inverted ordering; finally, type III textures remain unable to account for oscillation data with perturbative values of the couplings even in the case of a inverted ordering of the neutrino spectrum.
Before moving to our final remarks, it is worth mentioning that it would have been theoretically appealing to have FN-inspired textures for both y andỹ couplings; as can be indirectly inferred from the above discussion, we did not succeed in finding phenomenologically viableỹ couplings with textures mirroring those of y. Figure 6: Predictions for CR(µ−e, N) and BR(K + → π + νν) associated with randomly sampled y andỹ (see text), for type I, II and III textures (top to bottom, left to right). The horizontal/vertical lines denote the current experimental bounds, and the colour code identifies perturbative (blue) or non-perturbative (red) regimes ofỹ.

Concluding remarks
In this work we have carried a comprehensive phenomenological study of a SM extension via two scalar leptoquarks h 1,2 and three generations of triplet neutrinos Σ i R , further reinforcing the SM gauge group via a discrete Z 2 symmetry under which h 2 and Σ i R are odd (all other fields being even).
The present New Physics construction aims at simultaneously addressing two long-standing SM observational problems -neutrino mass generation, and a viable dark matter candidate -while further offering a solution to the currently reported anomalies in B meson decays, R K ( * ) .
The Z 2 symmetry ensures the stability of the LZoP, rendering the neutral component of the lightest Σ R a viable cold dark matter candidate for well defined intervals of its mass. In the absence of a full theory of flavour, we have identified several classes of flavour textures for the h 1 leptoquark Yukawa couplings which succeed in saturating the R K ( * ) anomalies. These textures (loosely based on Froggatt-Nielsen inspired ansätze) were subjected to a vast array of flavour conserving and flavour violating observables (including meson decays, neutral meson oscillations and cLFV decays), which allowed to infer stringent constraints on the h 1 leptoquark mass and couplings. Contrary to previous claims in the literature, our findings suggest that the strongest constraints on these leptoquark extensions do arise from cLFV µ−e conversion in nuclei, and from the rare K + → π + νν decays. Furthermore, we also verified that numerous ansätze (identified as promising ones for leptoquarks couplings, see e.g. [74,75]) were in fact phenomenologically disfavoured by several of the here considered observables.
It is important to emphasise that the constraints on leptoquark couplings arising from flavour observables are not intrinsic (nor peculiar) to the leptoquark realisation here considered; in fact these are valid for any SM extension via scalar triplet leptoquarks.
The present BSM realisation leads to a scenario in which neutrino masses are radiatively generated (at the three-loop level, from the exchange of leptoquarks, down-type quarks and lepton triplets, see Fig. 1). Neutrino oscillation data can be accounted for by means of a modified Casas-Ibarra parametrisation: avoiding non-perturbative regimes for the Yukawa couplings, y andỹ, establishes the final constraints on the parameter space of the model.
The inclusion of Majorana states opens the door to lepton number violating processes; the radiatively induced masses for the light (left-handed) neutrinos are one such example. The new interactions and couplings further allow for additional sources of CP violation. It is thus only natural to envisage the possibility of accounting for the baryon asymmetry of the universe. In the present realisation, one can have tree-level processes which are lepton number violating: one such example can be obtained from the neutrino mass (loop) diagrams -see Fig. 1 -, by "cutting" the inner fermion propagators. This would lead to tree-level LNV decays of the heavier neutral Σ 2,0 R into, for instance, Σ 1,0 R +d+d+dν L +dν L (which could have CP violating interferences with higher order diagrams). However, these appear to be heavily suppressed processes and, in the absence of a detailed evaluation, it remains unclear whether one could indeed generate a significant lepton asymmetry. In addition, Σ R decoupling would be required to occur above the EW phase transition to have an efficient conversion into a baryon asymmetry.
In summary, and following a thorough study of an extensive array of observables, we have proposed realisations of a SM scalar leptoquark extension capable of accommodating neutrino oscillation data, a viable DM candidate, and saturating the observed discrepancies for R K ( * ) . We notice that the present construction cannot account for the tensions in R D ( * ) , nor for the discrepancy between observation and SM prediction in what concerns the muon anomalous magnetic moment. Should the latter persist, then the candidate model here studied will have to be extended, or then embedded in a larger framework [62].
In the near future, a number of high-intensity experiments will put the present leptoquark construction to the test, in particular several cLFV-dedicated facilities (searches for radiative and three-body muon decays, in addition to neutrinoless conversion in nuclei) and a possible measurement of the rare decay K + → π + νν. Hopefully, positive signals or new stringent bounds emerging from negative searches, will allow to further constrain the model's parameter space, or possibly disfavour it as a candidate New Physics model.