A near-minimal leptoquark model for reconciling flavour anomalies and generating radiative neutrino masses

We introduce two scalar leptoquarks, the SU$(2)_L$ isosinglet denoted $\phi\sim(\mathbf{3}, \mathbf{1}, -1/3)$ and the isotriplet $\varphi\sim(\mathbf{3}, \mathbf{3}, -1/3)$, to explain observed deviations from the standard model in semi-leptonic $B$-meson decays. We explore the regions of parameter space in which this model accommodates the persistent tensions in the decay observables $R_{D^{(*)}}$, $R_{K^{(*)}}$, and angular observables in $b\to s \mu\mu$ transitions. Additionally, we exploit the role of these exotics in existing models for one-loop neutrino mass generation derived from $\Delta L=2$ effective operators. Introducing the vector-like quark $\chi \sim (\mathbf{3}, \mathbf{2}, -5/6)$ necessary for lepton-number violation, we consider the contribution of both leptoquarks to the generation of radiative neutrino mass. We find that constraints permit simultaneously accommodating the flavour anomalies while also explaining the relative smallness of neutrino mass without the need for cancellation between leptoquark contributions. A characteristic prediction of our model is a rate of muon--electron conversion in nuclei fixed by the anomalies in $b \to s \mu \mu$ and neutrino mass; the COMET experiment will thus test and potentially falsify our scenario. The model also predicts signatures that will be tested at the LHC and Belle II.


Introduction
The detection of neutrino oscillations establishes that neutrinos have small, but nonzero, masses and that the flavour and mass eigenstates do not coincide. However, the dynamical origin of the tiny scale of neutrino masses remains a mystery, as does any potential impact on the flavour structure of the Standard Model (SM). A plausible explanation for the lightness of neutrinos is to not explicitly introduce a tree-level mass term, but rather to engineer the generation of mass at loop-level: radiative models. We restrict our attention to models that induce a Majorana mass term and violate lepton-number by two units (∆L = 2). The magnitudes of these masses are naturally loop-suppressed 1 . As such, this mechanism gives a neat explanation for the disparity between the sizes of neutrino masses relative to those of other SM fermions.
Extending the SM to permit neutrino-flavour violation invites a more thorough investigation of the flavour sector. The connection between neutrino physics and flavour physics is of interest given the capacity of precision measurements to explore parameter space for such beyond-the-SM (BSM) models. Indeed, a variety of measurements have hinted at violation of Lepton Flavour Universality (LFU) in precision observables: of particular interest are angular parameters [2,3] and branching ratios in the b → s transition, and in the ratios R D ( * ) [4][5][6][7][8][9][10] and R K ( * ) [11][12][13]. Such observables are often very sensitive to the virtual effects of exotics, such as those that are typically introduced in radiative neutrino mass models.
In this paper we address both of these problems by developing a radiative neutrino mass model that features BSM contributions to flavour observables. We begin by outlining the anomalies we aim to address, before detailing a motivation for this model and a summary of related previous work.

Neutrino masses and mixing
Neutrino oscillations in vacuo depend on non-degenerate neutrino masses and lepton flavour mixing. The Pontecorvo-Maki-Nakagawa-Sakata (PMNS) leptonic mixing matrix U relates the neutrinos in their flavour and mass eigenstates. In the Majorana case, U can be written as the product of three rotations, the second of which depends on a phase, and a diagonal matrix of phases P:  Table 1: Summary of the average neutrino oscillation parameters from the NuFIT collaboration [14], with the assumption of normal ordering. The angles, θ ij , can be taken without loss of generality to be within the first quadrant. We implement these central values throughout this work. Majorana phases α 1,2 are entirely unconstrained.
An overview of the most recent global fit to these oscillation parameters is given in table 1, from the NuFIT collaboration [14]. These fits show a preference for the so-called normal hierarchy of neutrino mass. This mass-ordering mimics the generational indices: that is to say, m 1 < m 2 < m 3 , where m i represent mass-eigenvalues of the associated combinations of flavour states. An inverted ordering is less preferred by fits to data, but still represents a viable alternative regime. We note that the Dirac phase δ CP is poorly constrained by current data, and the Majorana phases α 1,2 are entirely unconstrained. We exploit this freedom later in our analysis.

Anomalies in flavour
A number of flavour observables appear to indicate flavour-dependent coupling of BSM physics to SM particles. A strong hint of LFU violation exists in experimental measurements [15,16] of semileptonic B-meson decays. A summary of the most recent experimental and theoretical averages for these observables may be found in table 2.

Neutral Current Processes
The anomalies include those in the rare b → s flavour-changing neutral current (FCNC) transition. Recent measurements [11][12][13] have reinforced deviations from the SM in the ratios R K * and R K : Other branching ratios of exclusive decays in b → sµµ have also been measured to be in tension with the SM prediction [17,18]. Discrepancies are also present in measurements of angular observables in B → K * µ + µ − , and the most significant of these is seen in the quantity P 5 [2,19,20]. The CMS measurements of these angular observables are consistent with the SM values [21]. Recent fits to the anomalous b → s data present large preferences for new physics in operators contributing to b → sµµ [22][23][24][25][26][27][28][29].

Charged Current Processes
Persistent tension in the semi-leptonic transition b → cτ ν has been observed independently by many experimental collaborations: through the decay B → Dτ ν by Belle and Babar [6][7][8][9][10]; and B → D * τ ν by Belle, Babar and LHCb [4][5][6][7][8][9][10]. Each experiment has measured deviations from the SM in the following quantities: Together these measurements amount to a deviation 3σ [32,35,36] from SM predictions. Although the ratios R D ( * ) are our primary concern in this work, we also introduce a number of other observables relevant to the charged current process that form the basis of predictions of our model. Specifically, in section 5 we present the predicted values for the observables R J/ψ , f D * L , and various tau polarisation asymmetries. The first of these is the ratio of the tauonic mode to the muonic mode for B → J/ψ ν, measured recently by LHCb to be R J/ψ = 0.71 ± 0.17 ± 0.18 [37]. Although the ratio is also measured to be enhanced with respect to the SM prediction R SM J/ψ ≈ 0.25-0.29 [38][39][40][41][42][43][44][45][46][47][48][49][50], the central value of the measurement shows a very large effect that cannot be wellaccommodated with BSM contributions [35]. The observable f D * L , the longitudinal polarisation of the D * in B → D * τ ν, also differs from the SM expectation by ∼ 1.6σ: as measured by the Belle collaboration [51], and has been shown to have good discriminating power for BSM explanations of R D ( * ) . The third class of observables we consider are tau polarisation asymmetries (see ref. [52] for a detailed discussion in the context of explaining R D ( * ) ). The polarisation asymmetry in the longitudinal direction of the τ in the D * mode has also recently been measured by Belle [10]: Although the errors are large, the projected Belle II sensitivity at 50 ab −1 for the same observable in the D mode is estimated at about 3% [53], and we expect the P * τ to be measured even more precisely at Belle II.
1.3 Muon anomalous magnetic moment, (g − 2) µ Precise measurements of the deviation in the semi-classical value of the muon gyromagnetic ratio, g µ = 2, have demonstrated an inconsistency. This is parameterised by the quantity There is a persistent deviation between the SM prediction and the experimentally measured value [54,55], corresponding to a 3.6σ anomaly. The error values refer to the experimental and theoretical prediction errors, respectively. Similarly, recent experimental results have indicated a deviation from the SM for the electron anomalous magnetic moment, of 2.5σ significance [56]. The leading candidates to explain these anomalies involve flavour-dependent, loop-level, BSM effects [55].

Neutrino mass and the one-leptoquark solution
Effective ∆L = 2 interactions involving SM fields were systematically studied by Babu and Leung [57] up to mass-dimension (D) eleven. By opening-up such operators at tree-level, and looping-off external fields, neutrino mass is generated at loop-level: radiative neutrino mass generation. Ref. [58] investigated the D = 7 operators in further detail, assessing the viability of minimal UV-completions for yielding neutrino masses consistent with the observed values. They identified the particle content of such completions, and explored the explicit phenomenology of one particular model: a completion of O 3b = (L i Q j )(L k d c )H l ij kl , involving the introduction of a scalar leptoquark (LQ) field, φ, and an exotic vector-like quark, χ ('Model 2' in table 3).
It is important to note that generating neutrino mass in 'Model 2' relies explicitly on mixing of the vector-like exotic with the SM b-quark. A direct consequence of this is that the mixing parameters are heavily constrained by measurements of b couplings and associated observables [59,60]. A similar radiative neutrino mass model ('Model 1' in table 3) containing the isotriplet, instead of the isosinglet, was also identified by ref. [58], although the implications of this model were not thoroughly explored.
The isosinglet leptoquark, φ, also features in a study by Bauer and Neubert [61] as a simple explanation for the flavour anomalies R D , R D * , and the anomalous b → s data. This Table 3: Particle content of two key radiative models identified by ref. [58]. The tuple entries refer to transformation properties under SU (3) C ⊗ SU (2) L ⊗ U (1) Y , and here we adopt the hypercharge convention Q = I 3 + Y . model was further studied in refs. [62][63][64][65], where its viability as a combined explanation for the flavour anomalies was evaluated in more detail. An idiosyncrasy of the model is that the charged current R D ( * ) anomalies are mediated at tree-level, while the b → s anomalies are explained by box diagrams with the leptoquark in the loop. We have noted the role of this LQ in models of radiative neutrino mass (see, e.g. [58,66,67]), and the connection between radiative neutrino mass and the flavour anomalies has also been explored more broadly in the literature [63,[68][69][70][71][72][73][74]. The authors of ref. [63] studied the overlap between the flavour anomalies and a two-loop radiative neutrino mass model containing the leptoquark φ. In that model mild tensions exist in explaining both R D ( * ) and R K ( * ) ; at best, it was found that this model could reconcile these anomalies together to within a 2σ region, if the model was restricted to the minimal particle content. If, however, the isosinglet leptoquark did not contribute to the b → s anomalies, it was found that a combined explanation of neutrino mass, (g − 2) µ and R D ( * ) was possible in the minimal model.

A motivation for near-minimality
Leptoquarks as BSM candidates have experienced a resurgence of interest in recent years. While use of scalar and vector LQs in constructing models of LFU-violation has a long history, the new measurements, particularly of anomalies, have provided additional motivation for these exotics. Additionally, they can play a pivotal role in models for explaining other SM problems, as detailed in ref. [75]. It is also interesting that many LQ models are motivated by unification, as they offer a direct portal between the quark and lepton sectors.
The limited success of the one-leptoquark solution of Bauer and Neubert motivates us to explore next-to-minimal models to explain these anomalies. Noting the relevance for flavour observables of introducing b-quark mixing to the models presented in table 3, vector-like fermion extensions to the SM are particularly intriguing. Further, it is known that the interactions of the isotriplet 2 ϕ [76] contribute to b → sµµ transitions at treelevel, whereas the isosinglet scalar leptoquark only generates loop level contributions. We therefore propose a merger of 'Model 1' and 'Model 2', to capture the beneficial features of both.
Neutrino oscillations imply a violation of family lepton number, whereas the flavour measurements imply a violation of LFU. This work will aim to explore the connection between these two phenomena, and extend upon earlier work in the field to construct a nonminimal model with a broader scope for explaining deviations from the SM of experimental observations.
The remainder of this work will be structured as follows. Section 2 will outline the mathematical structure of the model as a completion of a dimension-7 effective operator for radiative neutrino mass generation. Section 3 will develop the calculation and framework of contributions to the aforementioned flavour anomalies. Section 4 will proceed to investigate additional relevant flavour constraints and observables, establishing the phenomenology of this model. Section 5 contains the results from an investigation of parameter space, implementing aforementioned constraints. Section 6 contains a discussion of implications and prospects of this work.

The Model
Combining 'Model 1' and 'Model 2' we arrive at the BSM field content: In the SM gauge basis, the complete set of Yukawa interactions between SM and BSM fields are described by the following Lagrangian: with the BSM portion of the scalar potential, V, given by: where the notation '[ ] i ' means that bracketed fields are combined through a tensor product to produce the i-dimensional representation of SU(2) L . Additionally, the set of SM-BSM Yukawa interactions is given by the interaction lagrangian, L int : We may introduce indices {i, j} on Yukawa couplings, to reference the relevant fermion flavour 3 , for example: For clarity we often omit flavour indices, however when present we adopt the convention that for terms with two flavour indices, the first 4 refers to the lepton and the second to the quark. 3 Note that these generational indices will occasionally be replaced with their associated particle symbols (e.g. y23 → y µb ). 4 Although this is not often the convention in the literature, we choose this so that Lepto-Quark can act as a mnemonic, maintaining convention with ref. [58].
Note that we have imposed global U (1) B baryon-number conservation by switching off possible di-quark couplings for both leptoquarks. This circumvents bounds from proton stability [75,77].

Vector-like bottom partner and associated mixing
An appealing feature of this model, particularly for resolving anomalies in B-physics, is mixing of the vector-like fermion χ with the SM down-type quark. To avoid obvious constraints that would arise from mixing with the lighter flavours, we restrict it to the third generation b-quark. This mixing is a result of the term: where we have introduced the notation 'ˆ' to represent the gauge-basis fermion eigenstates.
A feature of the one-loop neutrino mass models derived from dimension-7 operators is that the neutrino mass matrix is proportional to the mass matrix of the SM fermion in the loop [58]. In this case, the coupling of the vector-like quark to the b quark dominates, and for this reason we restrict the consideration of theχ 1 mixing with the down-type quarks to the b quark. This mixing is chiral and in general quantified by two mixing angles, θ L and θ R , for the left-and right-chiral sectors, respectively: The associated mass matric is defined through where y b is the b-quark Higgs Yukawa coupling. Through singular-value decomposition, the unitary rotation matrices, V L and V R , rotate the fields into the mass basis: We bring the matrices to diagonal form by alternately left-and right-multiplying V † L MV R with its hermitian conjugate: The diagonalisation produces the following expressions for the mass and mixing parameters, in the limit that m χ m bB : (2.14) (2.15)

Structure of eigenbasis mapping
To calculate flavour observables, we first transform into the charged-fermion mass eigenbasis. Beginning with the SM fields, where {i, j, l} represent flavour indices, and defining L a , R b as unitary rotations 5 between the gauge and mass bases, we have that (2.16) Above, the 'ˆ' on the down-type quarks and exotics indicate that they are yet to be fully rotated into the mass basis -i.e. the mixing is yet to be incorporated. Hereν represents the neutrino weak-eigenstate, which is related to the neutrino mass-eigenstate basis via the usual PMNS matrix. We parameterise the PMNS matrix by the central values quoted in the NuFit collaboration 2018 global fit (table 1) [14], and in our numerical analysis we scan over values for the Majorana phases.
Redefining the coupling constants to absorb this transformation, the interaction Lagrangian becomes: (2.17) The correspondences between the couplings are given by the following redefinitions: (2. 18) in terms of unphysical mixing matrices. The matrices y Lϕ and y Lφ are related by the physical CKM matrix to x Lϕ and x Lφ , respectively, so they are not independent: The following holds for the mass-mixing of SM with BSM fields: Note that we have denoted the mass-eigenstate by χ 1 rather than the oft-used 'B' to avoid later confusion when discussing B-meson decays. At this point, it is important to explicitly note that some quartic field couplings in the scalar potential, V (eq. (2.4)), will be set to zero, for simplicity, in subsequent calculations -in particular, those which generate ϕ-φ mixing and mass-splitting between triplet components 6 , λ m and λ Hϕ 2 in 2.4. For the remainder of this work we will take that each of the isotriplet components is degenerate in mass, m ϕ .

Massive neutrinos: a tale of two regimes
The leading-order contributions to neutrino masses are from the one-loop diagrams in figure 1. After EWSB, we calculate the radiatively generated neutrino mass matrix, in the limit that m b m χ , to be: where the relative factor of two arises from the Clebsch-Gordon coefficients in eq. 2.17. Upon inspection of eq. 2.21, we note that there are two copies of similarly structured contributions: one from ϕ, and one from φ. In our analysis we consider, for simplicity, two distinct phenomenological regimes: Regime 1: The contribution to the neutrino masses comes solely from the isotriplet leptoquark ϕ (corresponding to Model 1).

Regime 2:
The contribution to the neutrino masses comes solely from the isosinglet leptoquark φ (corresponding to Model 2).
We leave a full exploration of the overlap of these phenomenological regimes, and the implications of LQ mixing, to future work.

Implementing a Casas-Ibarra-like parameterisation
Deriving the physical neutrino masses {m 1 , m 2 , m 3 } requires that we diagonalise eq.(2.21), via the standard PMNS matrix Expanding this relationship, we can express the neutrino mass matrix as a rank-2 matrix in terms of the low-energy parameters U ij and m i : The value of k depends on the hierarchy assumption for SM neutrinos. In each case, the lightest neutrino is taken to be massless in order to obtain the rank-two flavour structure of each term in eq. (2.21). For definiteness, we adopt normal ordering and set k = 3 in what follows. We use the notation η ∈ {φ, ϕ} to denote the LQ that contributes to neutrino mass generation: where δ = 0 for η ≡ ϕ, and δ = 1 for η ≡ φ. The matrices x Lη and y χη are column matrices of Yukawa couplings. Under the specified assumptions, couplings of the form are mutually constrained by neutrino oscillation measurements. We adopt a Casas-Ibarra-like procedure [78] and parameterise our ignorance of the coupling constants by introducing a parameter ζ ∈ C through re-expressing eq. (2.24) as where, upon expansion, ζ cancels out. This means that the neutrino oscillation parameters can be fitted for any value of ζ. Matching eq. (2.26) to eq. (2.23), we obtain expressions for the coupling matrices in terms of ζ, m i and u i : Parameterising the Yukawa couplings in terms of ζ enables us to efficiently scan over the parameter space that agrees with the oscillation measurements. The coupling values are inputted at the high-energy scale -consistent with the energy-scale of the LQ mediator. Upon assuming normal ordering as detailed above, the measurements of physical masssquared values can be translated into measurements of m 2 i : To demonstrate that our model fits neutrino mass, we simply need to show that we can assign reasonable values of ζ and m 0 , within appropriate limits defined by flavour-violating processes and perturbativity of the generated couplings. Imposing a perturbativity bound p on the couplings constrains ζ and m 0 such that ∀j: We impose these constraints with p = √ 4π in our analysis. The remaining couplings in this model remain free parameters to be assigned values in accordance with constraints in subsequent chapters.

Ameliorating Anomalies
In section 2.3 we identified the couplings in each mass regime which are fixed by the Casas-Ibarra parameterisation. To address the remaining goals of this model, it remains to calculate the corrections to the anomalous processes outlined in section 1: R K ( * ) , R D ( * ) and (g−2) µ . To parameterise the BSM contributions to these processes we frame our constraints in terms of effective operators, O i , weighted by the Wilson coefficients C i , such that the effective lagrangian at a particular energy scale is given by: The set {O i } represents an operator basis that encompasses the interactions of this model at low energy, usually corresponding to below the mass scale of the BSM mediator. There are a number of commonly used EFT bases, so we will be careful to specify the basis of interest and coefficient normalisation as we proceed with our discussion of constraints. Where these are relevant for application in computational procedures, they will be referenced in accordance with the Wilson Coefficient exchange format (WCxf) [79].
The leading-order contribution from our model to the b → sµµ transition is given by the isotriplet, ϕ, via the tree-level diagram shown in figure 2. The isosinglet φ also contributes a one-loop box contribution to this process, as detailed in refs. [61,63]. Using identities summarised in appendix B, the diagram in figure 2 corresponds to a BSM contribution to the effective operator O µµ LL : Typically, fits to the available experimental data on the b → s decays [80] involve the chiral-basis containing O µµ LL , and the related operator O µµ LR : These find a good fit to the data so long as the following expressions are satisfied [28] 7 : where the multiplier corresponds to normalisation for this particular EFT, quoted in this way to ensure consistency with eq. in R K ( * ) and significantly improves the SM fit to all of the b → s data, with a total pull from the SM of 6.5σ [28].
As proposed in ref. [61], the isosinglet-only model provides a one-loop contribution to these operators, and in general generates non-zero values for both C µµ LL and C µµ LR . In this case, fitting C µµ LR ≈ 0 requires positing a suppression of the appropriate right-handed Yukawas y Rφ . In contrast, the leading-order contribution generated by the SU (2)-triplet ϕ is at tree-level, and has C µµ LR ≈ 0 at the high scale. For consistency with the literature, we translate the quoted operators to the so-called 'flavio' basis for the Weak Effective Theory (WET), as outlined in ref. [79]. Implementing the Fierz transforms as quoted in Appendix B, the operator mapping for O µµ LL is: with the operators O µµ 9 and O µµ 10 defined as: We label the corresponding BSM contributions to these operators C µµ 9 and C µµ 10 . They are generated at tree-level by the ϕ LQ with coefficients Analogous to the fits quoted above, we take the central value for C µµ 9 = −C µµ 10 as (3.9) Note that the couplings above are derived from purely real values for the Wilson coefficients, and most fits assume real values for the Wilson coefficients. To ensure a C µµ 9 = −C µµ 10 consistent with these, we fix the value of C µµ 9 such that Im C µµ 9 ≈ 0. In Regime 2, this corresponds to scanning over real-only values for the two free-parameters x Lϕ * 22 and x Lϕ 23 , whereas in Regime 1 this corresponds to constraints on the Casas-Ibarra parameterisation. These will be discussed in section 4.4.
The contribution from this model to the b → cτ ν transition is described by three diagrams, as illustrated in figure 3. The dominant BSM effects manifest in contributions to the following 8 effective operators, expressed in the flavio basis for the Weak Effective Theory (WET): The neutrino index, here denoted j, can run over all three flavours -the leptoquark interactions need not conserve lepton flavour, and final-state neutrino flavour is rarely a direct observable, particularly in collider studies. Consequentially, the BSM contributions to the operator coefficients, calculated at the LQ mass scale, are as follows: We will often drop the neutrino-flavour index j to avoid unnecessary clutter when j = 3, the case that leads to constructive interference with the SM contribution.
In the limit of small mixing between χ and the SM b-quark, these Wilson coefficients are generated solely from the leftmost diagram in figure 3. The other diagrams are subject to suppression 9 by sin θ R ≈ 0. (See section 4.3 for a discussion of the constraints on mixing parameters).
Whilst the QCD Ward identity implies that the vector coupling C V L ,j does not run with energy scale, the relationship (multiplier) between the tensor and scalar couplings will change dramatically with energy-scale running. This will be accounted for in subsequent calculations. 8 The Fierz transformations used to derive these operators from the diagram in figure 3 are given in appendix B. 9 In both cases the topmost vertex originates from a nonzero coupling y χϕ i , therefore the effective coupling to the b-quark after SM-BSM mixing is proportional to sin θR. Figure 3: Dominant model contributions to the charged-current b → cτ ν process, not assuming a global conservation of lepton flavour, where j ∈ {e, µ, τ }.

Fit
Best fit 1σ region 2σ region -- Table 4: The results of our fit to R D and R D * including the new Belle combined measurement [32]. The first row shows the best fit point and σ-regions fitting to C V L with all other operator coefficients vanishing. The second row shows the same for C S L (Λ) = −4C T (Λ) for Λ = 2 TeV and all other coefficients set to zero. The third row shows the best fit point for a 2D fit to Re The values required for these Wilson coefficients to give a good fit to R D and R D * have been studied in the literature, typically under the assumption of lepton-flavour conservation. Existing fits [58] suggest a good match to data can be attained with contributions to the vector operator C V L . Contributions in the direction C S L = −4C T [63,65,81] can also provide a good fit to data, and this approach is subject to fewer constraints (See section 4).
We incorporate the new Belle combined measurement [32] into a fit of all measurements of R D and R D * using the fitting software flavio [82] 10 . The fit contours are shown in figure 4, with the fit excluding the new Belle measurement shown with dashed contours to indicate its effect. We find the best-fit point for the 2D fit to ReC V L and ReC S L (Λ) = −4C T (Λ) at Λ = 2 TeV. We also fit to C V L with C S L (Λ) = 0 and vice versa. These results are summarised in table 4. We comment here that in this model, for the isosinglet LQ φ contributing to the direction C S L (Λ) = −4C T (Λ), the vector operator will also always be non-zero. This follows from the relation in eq. (2.19).
The leading contribution where only x Lφ 33 is non-zero is suppressed by |V ts | ≈ 0.04, but this contribution can still be sizeable if x Lφ 33 is chosen to be large. A short discussion of this scenario is included in our phenomenological analysis (see section 5.2).   Table 4 for central values and the text for more details. Figure 5: Dominant contributions to the muon magnetic moment from BSM content in this model. Fermion lines are omitted as there are multiple valid assignments possible for these topologies, each of which will be considered in calculations.

Muon magnetic moment:
There are two viable leading-order diagrams for this correction ( figure 5). Since contributions from these diagrams are well-established in the literature [61,75] we simply quote and 10 We note that our fit does not include the measurements of f D * L and R J/ψ , since errors here are still large. Instead, we take the central values from our fits and discuss predictions for these observables in section 5. interpret the results. In the limit that m 2 ϕ , m 2 φ m 2 t , the contributions are given by: As m µ m t , the second, like-chirality, terms are small -leading to the requirement of non-vanishing right-chiral couplings to obtain an adequate fit to the anomaly. We are left with the contribution generated by φ in this limit: In both neutrino mass generation regimes at least one of the two necessary Yukawa couplings is a free parameter unconstrained by contributions to the other anomalies outlined in this section.

Constraints
Below we discuss the constraints relevant to our model and the limits we require in our subsequent analysis. We restrict our main discussion to what we consider to be the minimal scenario to explain the B anomalies and (g − 2) µ . Here, the isotriplet LQ ϕ explains the neutral current anomalies, while the SU (2)  The leptoquark that participates in the neutrino mass generation must have a non-zero Yukawa coupling to the electron, and this is the most import phenomenological consequence for the constraints we consider. This, together with the relation in eq. (2.19), means that constraints from processes involving the first generation of SM fermions cannot be avoided completely. In fact, the hierarchy present in the leptoquark couplings to charged leptons is fixed by measured PMNS matrix elements, while the couplings to light quarks are suppressed by CKM matrix elements. Explicitly Of course, the Lagrangian in eq. (2.17) contains many more parameters than these. For simplicity, we turn off any couplings not immediately related to the anomalies or neutrino mass. In reality these need only be small enough to respect any limits placed on them by experiment.
The choice of the minimal set of Yukawa couplings depends on the choice of neutrinomass regime. In this section, expressions are given in generality assuming x Lφ 13 = 0 and x Lϕ 13 = 0. If one chooses a particular regime, then the Yukawa coupling to the electron of the leptoquark that does not participate in the neutrino mass can be switched off, and we make this choice according to the principle of minimality discussed above. Thus only one leptoquark will have a Yukawa coupling to the electron at a time. Of course an additional phenomenological consequence of choosing a neutrino-mass regime is the absence of ∆L = 2 interactions for one of the leptoquarks. This can be achieved by turning off the associated couplings y χη i . In the absence of LQ mixing, the term is not generated at any order since the interactions of η now conserve lepton number.
Below we summarise these comments with concrete Yukawa-coupling textures. The constraints presented in this section assume the following set of non-zero Yukawa couplings: However, if the restriction is made to regime 1, it is understood that x Lφ 13 , x Lφ 23 = 0 and y χφ i = 0, while regime 2 implies x Lϕ 13 , x Lϕ 33 = 0 and y χϕ i = 0. We do discuss other Yukawacoupling textures throughout this section where appropriate. Notably, we comment briefly on explaining R D ( * ) with contributions only to the vector operator C V L , and the constraints associated with this scenario are presented in this section as well.
This parameter space is explored in the context of the constraints implied by fits to the flavour anomalies and neutrino mass. We use a suite of computational machinery for most of the calculations, and this setup is discussed in section 4.1. Where appropriate we explicitly write out the dominant contributions to observables where we consider that this provides useful insight. Some observables are also calculated separate to these methods, and these are also discussed in detail below.

Calculation pipeline
Using SARAH [83,84] we construct the model from the Lagrangian upward, using inbuilt machinery to encode the algebraic structure of the fields, associated global symmetries and mixing. SARAH generates an output module for use with SPheno [84], which can calculate the Wilson coefficients, decay rates and a subset of flavour observables, defined by FlavorKit [83], for a particular assignment of model parameters. A full discussion of the underlying machinery and symbioses of these programs can be found in ref. [85].
In addition to the above, Flavio [82] was utilised to process manually calculated Wilson coefficient dictionaries where appropriate, or to take as input the Wilson coefficient .json files outputted by SPheno. This enabled us access to a broader class of flavour observable calculations than would have been otherwise possible. The running of Wilson coefficients in Flavio is implemented using the Wilson package [86].
By implementing this combination of computational machinery, it was possible to construct an efficient parameter scan over a large number of dimensions. In doing so, it was possible to determine the regions of parameter space that establish this model as viable for its desired purpose: reconciling the anomalies, generating radiative neutrino mass, and satisfying the most compelling experimental constraints, as we discuss in this section.

Validation and benchmark regions
To begin, we consider this model solely as a radiative neutrino mass model, only switching on the couplings generated by the Casas-Ibarra parameter. Truncating parameter space so as to only include those important for radiative neutrino mass generation, we may validate our calculations against the scans performed in ref. [58] for the φ-χ model (Model 2 in table 3). In figures 6 and 7, we depict the two parameter regimes fixed by the Casas-Ibarra procedure: neutrino mass generated by ϕ (Regime 1), or by φ (Regime 2). This is proof-of-principle to validate this computational setup, and more generally to determine benchmark regions of parameter space for subsequent scans. Throughout this section we have implemented the perturbativity bound that the magnitude of each physical Yukawa coupling is less than √ 4π. Ref. [67] identify the most important flavour-changing processes as the LFV decays: µ → eγ, µ → eee, as well as µ − e conversion in nuclei, all of which depend on the LQ couplings fixed by ζ. We have scanned over the magnitude of the Casas-Ibarra parameter, |ζ|. The observables in question are proportional to the product of two LQ couplings generated by the Casas-Ibarra parameterisation, and therefore, ∝ |ζ| 2 . Any phase information is irrelevant for these calculations.
For Regime 2, where φ generates neutrino mass, these plots represent replication of those in ref. [58], with updated measurements and using the framework outlined above. The minimal contrast between these two results acts as initial validation of the calculations  Figure 6: Regime 1 Constraints on ζ for varied isotriplet ϕ LQ mass (a), and vector-like quark mass (b); in each case the alternate mass is fixed at 2 TeV and the isosinglet LQ couplings are switched off. Allowed points for each constraint lie between the two samecoloured lines. The 'dip' in both graphs is due to an accidental cancellation by virtue of parameter choice.   in SPheno/SARAH for these LFV processes. Although we have not imposed any further constraints in these scans, nevertheless from figure 6 and figure 7 we can begin refining parameter space by noting that |ζ| ∈ (10 −2 , 10 2 ) roughly defines the region which is capable of accommodating these LFV constraints with perturbative couplings. This is consistent with the bounds prescribed earlier in eq. (2.29).

Collider bounds
Many collider searches already exist for third-generation quark 'partners', as they are referred to in many SUSY models. Here, we focus on those presented in the context of each vector-like quark component as a part of the χ isodoublet. Contributions to the decay branching fraction limit were found to be sensitive only to the BSM-SM quark mixing and the mass of the exotic, and the following analysis considers a mixing Yukawa Y b,3 = 1.
The branching fraction of the exotic bottom-partner χ 1 can be expanded as per the following; Br(χ 1 → Zb) + Br(χ 1 → Hb) ∼ 1. (4. 3) The third possible decay channel, χ 1 → W t, is highly suppressed by the small mixing between the exotic and SM quark sectors assumed for this channel. Taking the constraints from the ATLAS collaboration, and calculating the the corresponding branching fraction using SPheno, the resultant limit on the vector-like quark mass can be read off the graph in figure 8: This represents a very low bound from decays of χ 1 , and decays of χ 2 impose a slightly stronger constraint. The ATLAS collaboration imply a lower bound on the mass of χ 2 (or 'Y ' as it is referred to) of: The masses of each component of the isodoublet χ are of comparable orders of magnitude, as we will see in section 4.3. Therefore, we adopt the more conservative bound for m χ = m χ 2 ∼ m χ 1 from χ 2 decay in subsequent discussion.
Similarly to the vector-like quarks, approximate bounds on the leptoquark masses can be inferred from collider searches. Most searches constrain the masses of the isosinglet and/or the isotriplet to be above ∼ 1 TeV, subject to model-dependent assumptions. Masses of scalar leptoquarks η with a reasonable coupling to bν (as is required by our model to satisfy b → cτ ν), are constrained by the CMS collaboration [88]to be: m η 1100 GeV. (4.5) High-p T dilepton production through the ϕ leptoquark has been shown to provide interesting constraints and signatures for the leptoquarks in our model [65]. The leptoquark contributes to the processes pp → through tree-level t-channel graphs whose effects can alter the tail of the differential cross-sections for pp → . We take the limits from ref. [65] for the muon and tau modes derived from 36 fb −1 of ATLAS data at 13 TeV [89,90]. We derive bounds on bb → ee and extract the 3000 fb −1 ATLAS sensitivity for the electron and muon modes from ref. [91]. These bounds are shown in figure 9. We find that the limits on cc → ee and uu → ee give less stringent bounds on |ζ|, and thus we do not include them in our numerical scans. ij . The limits are from LHC searches in pp → high-p T tails at 13 TeV from ATLAS [92], derived from ref. [91]. The Yukawa coupling being constrained depends on the process. For example, ss → µµ will constrain x Lϕ 22 .

Limits on the b-quark/vector-like quark mixing
Of central importance in this model is the mixing generated by the terms Y b,3bR Hχ L + h.c. between the b quark and the vector-like quark χ. This mixing is a necessary ingredient for the violation of lepton number by ϕ, and plays a governing role in the overall scale of the neutrino mass m 0 according to eq. (2.24). Its size also dictates the extent to which ∆L = 2 neutrino final states are important to consider, for example in B → K ( * ) νν. The mixing of the b with χ leads to new contributions to the oblique electroweak parameters S and T . These have been measured to high precision by LEP [93]. The mixing also leads to an alteration of the Zbb coupling at tree-level, for which global electroweak fits have suggested a small deviation from the SM value, e.g. [94]. The dominant contributions to these effects are encapsulated in the effective dimension-6 Lagrangian generated by the heavy χ at the scales probed by experiment: The first operator modifies the electroweak precision observables discussed above and the second affects Higgs measurements and is currently poorly constrained. We take the 2σ bounds on the operator coefficient C 33 Hd from the electroweak fit performed in ref. [  This agrees with ref. [59] which studied the effects of the doublet χ and other vector-like quarks. The relation θ L ≈ m b mχ θ R from eq. (2.15) implies the cos θ L factors appearing in eq. (3.7) and eq. (3.13)-(3.15) do not suppress the contributions to the anomalous observables 11 . Restricting this mixing to be small consequentially reduces the mass-splitting between the components of the exotic doublet χ, such that m χ ∼ m χ 1 ∼ m χ 2 remains a valid approximation.

Leptonic decays
The leptoquarks contribute to the LFV decays i → j γ and i → j k k at loop level through diagrams like those depicted in figure 11, as well as H penguins and box diagrams for the latter. The strongest of these constraints on the couplings are from first-and secondgeneration LFV processes. The experimental bounds on these processes can be found in table 5.
The related process of muon-electron conversion in nuclei is mediated at tree-level by the leptoquark participating in the neutrino mass generation (figure 10), although there are also loop-level topologies similar to those for i → j k k with quarks in place of the sameflavour lepton-antilepton pair. We presume that a coherent conversion process dominates, i.e the final state nucleon is the same as the initial state, therefore the coupling to the sea-quarks 12 is negligible [97]. The tree-level contributions to muon-electron conversion depicted by figure 10 will lead to vector, scalar and tensor effective operator contributions, via Fiertz transformation (appendix B). These will be discussed for the two neutrino mass regimes below.   Figure 11: Contributions to i → j γ processes, where additional splitting of the boson line gives contributions to i → j k k . The photon/Z can be attached to any of the four lines.

Ruling-out Regime 2
For neutrino mass as generated in Regime 2 (φ couplings fixed by Casas-Ibarra), contributions to muon-electron conversion are dominated by tree-level processes via the LQ mediator φ. Assuming the right-handed couplings to first-and second-generation quarks and leptons are switched-off, the dominant contribution is from the effective interaction: LQ couplings between charged leptons and up-type quarks are unavoidably generated by CKM mixing: Ref. [98] contains a study of the effective Lagrangian approach to constraining these contributions, and provides model-independent limits which we will interpret here. The strongest constraints from muon-electron conversion in nuclei presently come from measurements involving gold (Au) nuclei. In the absence of accidental cancellation, the experimental constraint in table 5, results in the following bound on the dominant contribution to this measurement: (4.10) Since both Yukawa couplings involved are fixed by the Casas-Ibarra procedure, the bound can be interpreted in terms of |ζ|. We find 13 To contrast with other constraints, we begin by parameterising the Casas-Ibarra parameter via ζ = (1 + iκ)Reζ, giving |ζ| = |Reζ| √ 1 + κ 2 , where κ ∈ R. The only viable scenario for this leptoquark to explain the charged current anomalies alone is for it to contribute to the scalar and tensor operators [63,65,99]. As discussed in section 3.2, there are always unavoidable contributions to the vector operator C V L present as well once x 33 is non-zero. Our fit results suggest Re C S L ≈ 0.14, Im C S L ≈ 0, (4.12) which implies, assuming c θ L ∼ 1 and y Rφ 32 ∈ R, and inputting numerical values for U ij : Combining this with the above bound on Re(ζ) 2 (eq.(4.11)), and saturating the perturbativity bound with |y Rφ * 32 | = √ 4π, we obtain (4.14) Re-expressing the upper perturbativity constraint given in eq.(2.29) in terms of numerical inputs gives which, when combined with eq.(4.13), reduces to: These two upper bounds are shown in figure 12 as a function of κ. In order to ensure that Im C S L ≈ 0, we can rephrase this requirement as needing a suppression of Im C S L relative to Re C S L , that is to say a minimal requirement of: To satisfy all three of these constraints already reduces the available parameter space significantly (the grey-shaded region is allowed). These are necessary, but not sufficient, conditions for these parameters to satisfy in order to ensure that these constraints are met. We may note from figure 12 that these constraints alone indicate a very small allowed mass region of 1 TeV m ϕ 2 TeV. 13 Throughout this discussion we have assumed m0 is positive-definite. Figure 12: Constraints on neutrino mass Regime 2 (φ generating neutrino mass). Coloured lines represent upper bounds on the LQ mass-squared, and the shaded region represents the allowed points when also considering suppression of imaginary component of respective Wilson coefficient C S L . The red-dashed line in each indicates the approximate lower massbound from collider constraints.

Constraints on Regime 1
Following the calculation procedure above, we find for Regime 2 that the effective interaction contributing to muon-electron conversion is given by: Implementing the bound determined in ref. [98], as above, constraints from muon-electron conversion in Gold nuclei give the following: which implies, via the Casas-Ibarra parameterisation, and using the notation consistent with the previous section: Requiring agreement with earlier discussed value for ReC 9 = −ReC 10 fit eq. (3.8) leads to the following central value: Also, re-expressing the perturbativity constraint, combining eq. (4.15) and eq. (4.21), gives: These provide a necessary, but not sufficient, bound on these parameters to satisfy the specified constraints. Simply from contrasting the size of the pre-factors between equations 4.22 and 4.23,4.14 and 4.16, the unconstrained parameter space for mass-squared in Regime 1 is significantly greater than for Regime 2. This is justification for concentrating solely on Regime 1 for the remainder of this work.

Z decays
The leptoquarks φ and ϕ will modify the Z coupling to leptons through one-loop diagrams involving SM quarks and the vector-like quark χ. For the contributions involving leptoquarks and SM fermions we use the results of ref. [100], which include corrections due to the external momenta of the Z. The additional diagrams with the vector-like quark in the loop are shown in figure 13. We find the contributions to the leptonic Z couplings from these to be The couplings y χϕ i are inversely proportional to ζ, and thus we expect these contributions to be suppressed when ζ and the χ-b mixing parameter Y b,3 are sizeable.

Charm meson decays
Since couplings to up-type quarks and charged leptons cannot be avoided for the leptoquark that couples to χ, the physics of operators of the form O ijkl ∼ (u i Γu j )( k Γ l ) is important to study. Here, we consider the leptonic decays of the D 0 meson, since a sizeable coupling to the charm quark can assist in the explanation of the large effects seen in the charged current anomalies. The isosinglet LQ φ generates the entire spectrum of operators which can in principle contribute to the leptonic decays of the D 0 , since it interacts with both left-and right-chiral SM fermions. Concretely, the dimension-6 Lagrangian  Figure 13: The leading contributions to Z → and Z → νν in our model. Fermion arrows omitted for brevity, such that each diagram can be associated with multiple flow assignments. is generated with tree-level contributions from both leptoquarks: . We impose the experimental upper limit Br(D 0 → µµ) < 7.6 · 10 −9 [102]. Contributions to the electronic mode from the vector operators are helicity suppressed and we ensure |y Rφ 1i | 1 in all numerical scans to avoid contributions to the electronic scalar and tensor operators.

Bottom meson decays
The b → sνν transition provides one of the theoretically cleanest FCNC processes. Predictions for this transition are devoid of hadronic uncertainty, beyond the form-factors, unlike the b → s transition. This makes measurements of B → Kνν highly useful for constraining BSM in the flavour sector.
In this model, there are two types of contribution to this process, depicted in figure 14: one with two neutrinos in the final state, and one with a neutrino and its charge-conjugate. Assuming the χ-quark mixing angle is small, then the dominant contribution is to: parameterised by the BSM coefficients: We do not presume that the final-state leptons are of the same flavour -all combinations of i and j can be incorporated into a fit using the Flavio software. However, as was noted earlier, of particular importance are the contributions from the lepton-number conserving processes which can interfere with the SM contributions. Additionally, the recovery of lepton-number as a global symmetry in the large-mass limit for exotics motivates this parameter choice.
Contrasting the contributions to w νν with the structure of the BSM contributions to other anomalies, we observe significant overlap in the relevant observables for fitting measurements of B → Kνν and in explaining the R D and R D * anomalies. Consequentially, we expect b → sνν measurements to provide some of the most stringent constraints on free parameters relevant for BSM in the b → cτ ν transition.

Meson mixing
The process of B s -B s mixing provides a complementary constraint on the same couplings involved in the b → sνν processes discussed in section 4.7. It was found in refs. [61,63] that the process leads to a weaker constraint than B → Kνν and B → K * νν for low φ masses, but becomes relevant for masses larger than a few TeV. In our model, we have contributions from both the isosinglet and isotriplet LQs through box diagrams with neutrinos and charged leptons in the loop. These contribute to the operator C bs 1 , where colour indices are contracted within parentheses. The combination C Bs exp 2iφ Bs = ∆m exp s /∆m SM s is calculated using SPheno [84,85,104], and we impose the UTfit collaboration's result [105] C Bs = 1.110 ± 0.090 (4.40) in our numerical scans. We will work with small imaginary parts for the couplings fixed by the neutrino mass and we maintain φ Bs ≈ 0 for the phase, consistent with UTfit's result.

Summary of constraints
In tables 5 and 6 we present summaries of the constraints of this section. The tables contain the observables we consider in our later phenomenological analysis as well as the limits we require.

Results and Discussion
Below we explore the extent to which this model can accommodate the charged-and neutralcurrent anomalies, the anomalous magnetic moment of the muon and neutrino mass in light of the constraints presented in the previous section. First, we review the minimal setup introduced in section 4 and present the results of our Monte Carlo analysis. We comment briefly on non-minimal scenarios in section 5.2.

Quantity
Requirement

Monte Carlo analysis
In the minimal scenario the deviations in R D ( * ) are explained by the isosinglet leptoquark φ with contributions in the direction C S L (m φ ) = −4C T (m φ ), implying O(1) values for the couplings x Lφ 33 and y Rφ 32 [63,65,110]. Contributions to the vector operator are more heavily constrained since they necessarily imply large effects in B → K ( * ) νν and B s -B s mixing, in the absence of any kind of cancellation (see section 5.2). The φ particle also explains the anomalous magnetic moment of the muon with the values of y Rφ 23 and y Lφ 23 = x Lφ 23 fixed according to eq. (3.18). The limits derived in section 4.4.1 suggest the extent to which φ can contribute to the generation of neutrino masses is small. Since we consider suppressed LQ mixing, this means there is no connection between the neutrino mass mechanism and the anomalies in R D ( * ) and (g − 2) µ in this model. For this reason, we fix m φ and the couplings involved in eq.  [52,110,111]. Here we expand briefly on some of these.
The fit we present in section 3.2 does not include the less-precisely measured observables R J/ψ , f D * L and P * τ , introduced in section. 1.2.2. We instead use the preferred values from our fit to make predictions for these observables, concentrating on the scalar-tensor solution, since this is the easiest to accommodate with the φ LQ. We note that this solution gives negligible efficiency variation from the SM for the measurement in the D * mode [9], predicts Br(B c → τ ν) 30% [36,63,65] and displays a q 2 spectrum that agrees well with experiment [112].
In figure 15 we project the 2σ preferred region for C S L (see Table 4) onto combinations of b → c related observables to illustrate the ability of combined measurements to close in on this scenario. Were possible, we have also shown Belle II 50 ab −1 sensitivity [53] in grey centred around the SM prediction in black. Current measurements are shown in red with their 1σ errors in orange. With contributions in the scalar-tensor direction, the φ leptoquark's contributions to f D * L are in the opposite direction to current measurements, although still within the 2σ region. If the central value of f D * L stays close to where it is, or moves down slightly, the model would then predict P τ ≈ 0.4, which compromises the potential mild improvement in R J/ψ the model can offer. This scenario leads to a SM-like P * τ , but potentially large deviations in the P Figure 15: A grid plot of the various b → c related observables in addition to R D and R D * considered in our analysis. Solid black lines represent the SM predictions around which the grey shaded regions are the Belle II 50 ab −1 sensitivities [53], bordered by the black dashed lines. Red lines are current measurements and orange regions are their 1σ errors. Where the Belle II sensitivity is unavailable we present only the SM prediction without a shaded region. The blue points explain R D ( * ) to 2σ.
ImC µµ 9 /ReC µµ 9 = tan arg C µµ 9 as a measure of the relative size of imaginary part of C µµ 9 , we find tan arg C µµ 9 ≈ − cot arg ζ + as contours with varying arg ζ and α 2 . This illustrates the behaviour discussed above but also investigates the effect on the other leptonic couplings. It is evident that the choice arg ζ ≈ π/2 also leaves the tau coupling mostly real, although this cannot be said for the electron coupling where the dependence on α 2 is significant. We nevertheless proceed with the choice arg ζ = π/2 in our numerical analysis and account for the possibility of a large imaginary part in the coupling x Lϕ 13 . To explore the extent to which this scenario can explain the flavour anomalies and neutrino mass, we perform a random scan over the 5 free parameters of the setup: |ζ|, x Lϕ 22 , m ϕ , m χ , α 2 . Random values are drawn uniformly over the intervals defined for these parameters in table 7. Notably, the Yukawa coupling Y b,3 is fixed to 0.25m χ /TeV, the lower edge of the 2σ region from eq. (4.7) needed to explain the small discrepancy in Z → bb. We generate 2 · 10 6 points which are filtered through all of the constraints presented in section. 4.
As discussed at length in section 4.4, the leptoquark ϕ mediates highly-constraining processes of muon-electron conversion in nuclei at tree-level, and the couplings involved are  directly related to those that explain the neutrino masses and the b → s anomalies.
We find that only about 12% of the points in our numerical scan are rejected on the basis of a constraint, but from among these almost all are disallowed because they violate the muon-electron conversion bound given in table 5. In figure 17 we present the results of our random scan with slices through the parameter space and various projections. We find that our model requires muon-electron conversion in Gold nuclei at a rate no less than 2 · 10 −13 to accommodate the preferred value of C µµ 9 . The COMET experiment [114][115][116][117], has a projected sensitivity of Br(µAl → eAl) < 10 −16 at 90% confidence. This will provide an improvement on the current limit by four orders of magnitude, and will test and potentially falsify this scenario. Interestingly, our model cannot simultaneously avoid the muon-electron conversion bound and explain the b → s anomalies with a vanishing Majorana phase α 2 , a result made clear in figure 18a. There, it is also apparent that the constraint can be avoided for 100 • α 2 300 • , a region that overlaps with that shown in figure 16, needed for a small imaginary part for the electron couplings x Lϕ 13 . We find that an additional two constraints cut into the parameter space significantly: bounds from D 0 → µµ and the ATLAS measurement of ss → µµ. Our model predicts the D 0 → µµ rate to be an order of magnitude larger than estimates of the SM contribution Br(D 0 → µµ) SM ∼ 3 · 10 −13 [118] (see figure 18b), while the ATLAS 3000 fb −1 projected limit from ss → µµ indicates that a non-observation would almost entirely rule out the model for low LQ masses (see figure 18c).

Comments on explaining R D ( * ) with the vector operator
In our analysis above we consider only contributions in the scalar-tensor direction to explain the charged-current anomalies in R D and R D * . This choice is made to avoid the dangerous contributions to B → K ( * ) νν, which necessarily exist in the presence of a large C V L . We explored two ways these constraints could be avoided in the context of our model: 1. As discussed previously in the literature [63,120], one way to avoid the constraints from B → K ( * ) νν and B s -B s mixing is to explain the R D ( * ) anomalies with a large x Lφ 33 while ensuring x Lφ 32 ≈ 0. The coupling y Lφ 32 required to explain R D ( * ) is generated through eq. 2.19, while keeping the strange-quark coupling to the neutrinos zero. Combining eq. 2.19 and eq. 3.14 with x Lφ 33 x Lϕ 33 gives which implies 1.7 |x Lφ 33 |/(m φ TeV) 7.2 for cos θ L ≈ 1 to explain R D ( * ) according to our fit to C V L (see table 4). We note here that even saturating the lower 2σ bound Figure 17: The results of the random scan projected onto Br(µAu → eAu) and C µµ 9 . All constraints except muon-electron conversion in Gold have been applied. The solid black line represents the central value of the fit of ref. [28] to the anomalous b → s data, and the heavier and lighter shaded regions are the 1 and 2σ regions. The dashed red line corresponds to the current most stringent limit on Br(µAu → eAu) from SINDRUM II [119] and the black dot-dashed line is the projected sensitivity of COMET [114][115][116][117]. The plot on the right is an enlarged look at the interesting region of the plot on the left. The colour axis represents the value of the Majorana phase α 2 .
on C V L leads to contributions to Z → τ τ that disagree with experiment, despite the reduction of the global average driven by the latest Belle result.
2. Ref. [121] proposed that the R D ( * ) anomalies could be explained through the vector operator by considering a cancellation between the φ and ϕ particles in this model to the processes B → K ( * ) νν. This was further studied in ref. [64] and ref. [122]. We have investigated this suggestion in considerable detail for this model, and could not find any parameter space that could simultaneously resolve the R D ( * ) anomalies and be consistent with constraints from B s -B s mixing. Our findings are in agreement with ref. [122]. We note that ref. [122] proposed some lines of investigation, such as the use of complex-valued couplings constants, that could potentially alter this conclusion, but an investigation of such a scenario is beyond the scope of this paper.

Conclusions
In this work, we have established a near-minimal scalar leptoquark model capable of producing significant flavour-specific BSM effects, and radiatively generating Majorana neutrino masses. Combining two existing completions of D7 ∆L = 2 effective operators, this model consists of two scalar LQs, ϕ ∼ (3, 3, −1/3) and φ ∼ (3, 1, −1/3), and the vector-like quark doublet χ ∼ (3, 2, −5/6). We developed the structure of this model, including mixing between the vector-like exotic χ and the SM b-quark, and explored the significance of this effect for one-loop generation of radiative neutrino masses in two distinct phenomenological regimes: Regime 1 (ϕ contribution) and Regime 2 (φ contribution). The tree-level contri- The points shown pass all of the constraints considered in our analysis. The colour axis represents the imaginary part of C µµ 9 . The plot shows the preference away from a vanishing α 2 , driven by the constraint Br(µAu → eAu), and the consistency of the available parameter space with a small imaginary part of C µµ 9 . (b) The results of the random scan projected onto Br(D 0 → µµ) and C µµ 9 . Points shown pass all constraints. Our model predicts Br(D 0 → µµ) 10 −12 , about an order of magnitude larger than the SM estimate from ref. [118]. We note that our calculation is not valid below the dashed orange line since it only represents the new-physics contribution. (c) The plot shows the influence of the ATLAS ss → µµ limits on the parameter space of our model. Coloured points lie in the 2σ region of the b → s fit we use. Dark blue points cannot explain the b → s data. The dashed orange line corresponds to the current ATLAS limit, while the solid red line is the 3000 fb −1 projection. The abrupt absence of points in the bottom left of the plot is due to the constraint D 0 → µµ.
butions of each of these scalar LQs to the anomalous processes b → cτ ν and b → sµµ were established, within the context of an EFT framework.
We then discussed the experimental constraints on this model for each mass-generation regime. Regime 2 was found to be significantly constrained by a combination of µ − e conversion in nuclei and the required contribution to ameliorate anomalies in b → cτ ν processes. Meanwhile, it was found that Regime 1 is capable of explaining neutrino masses, charged-and neutral-current anomalies, (g − 2) µ and the minor deviation in Z → bb, whilst also avoiding other notable constraints. In order to avoid strong constraints from leptonflavor violating processes, whilst also accommodating the neutral-current anomalies, we needed to float the Majorana phase α 2 . To avoid the constraint from µ − e conversion in nuclei, this model showed an interesting preference for the region in which the generated Yukawa couplings were mostly real-valued, with 100 • α 2 300 • .
The established scenario was found to be highly predictive and extremely testable. For µ − e conversion in nuclei, we predict a rate of no less than 2 × 10 −13 , a regime testable by the projected sensitivity of the COMET experiment. We also predict rates of D 0 → µµ an order-of-magnitude larger than current estimates of the SM contribution. This model is also testable at the LHC, where dimuon searches at high p T provide strong limits. Additionally, future measurements of the asymmetry observable P τ by Belle II can test this model, which prefers P τ ≈ 0.4, assuming the central values of f D * L and R D ( * ) remain constant. and insert it into the third term in (A.1) to obtain After EWSB, the term proportional to |H 0 | 2 in eq. (A.3) generates mass-mixing between φ and ϕ 2 . We may parameterise this through: where θ m is the scalar mixing angle, R is the rotation matrix, and η 1 , η 2 are the mass eigenstates. In order to derive their masses, we begin by re-expressing the relevant Lagrangian terms in matrix form: Requiring RM 2 s R † to be diagonal gives the following mixing parameters (here we have assumed, for simplicity, that the parameter κ is real-valued): Nevertheless, although this mixing can produce additional contribution to observable effects, the mass-insertion is likely to render this sub-dominant to contributions that are independent of φ-ϕ mixing.

B Fierz transforms and the C-operator
The following field rearrangements will be useful in calculations throughout this work. The general Fierz transform for anti-commuting fields [124] is: where here each bracket represents the respective fields that sandwich the operators: {Γ A } = {P R , P L , P L γ µ , P R γ µ , 1 2 σ µν }, {Γ A } = {P R , P L , P R γ µ , P L γ µ , σ µν }. We will also employ the actions of charge conjugation on certain operators, as outlined below: