Next-to-leading order scalar contributions to $\mu\rightarrow e$ conversion

Within a class of models in which lepton flavor violation is induced dominantly by scalar particle exchanges, we estimate the $\mu \to e$ conversion rate in several nuclei. We include next-to-leading order (NLO) terms in the one- and two-nucleon interactions in chiral effective theory, rectifying some incorrect results in the previous literature. We provide an uncertainty budget for the conversion rates and we find that NLO contributions affect the amplitudes at the level of $10\%$, which could be larger than the uncertainty on the leading order couplings, dominated by the strange and non-strange nucleon sigma terms. We study the implications of our results for testing Higgs-mediated CLFV in the future by combining results from various experimental searches, such as $\mu \to e$ conversion in multiple target nuclei and $\mu \to e \gamma$.


I. INTRODUCTION
Lepton flavor violating processes involving charged leptons (CLFV) are among the theoretically cleanest probes of physics beyond the Standard Model (BSM). The minimal model extending the Standard Model with neutrino masses predicts CLFV amplitudes proportional to ∆m 2 ν /m 2 W [1][2][3][4], leading to branching ratios forty orders of magnitude below current experimental sensitivity and hence a huge discovery window. Moreover, CLFV processes also test the origin of flavor breaking in the lepton sector, ultimately related to the structure of neutrino mass matrices. So far, the strongest experimental bounds on CLFV have been given on µ → e transitions. The current limit on the branching ratio of µ → eγ is B µ→eγ < 4.2 × 10 −13 at 90% C.L. [5], while that of µ → e conversion in gold is B µ→e (Au) < 7 × 10 −13 at 90% C.L. [6]. The next-generation searches aim to achieve higher sensitivities. For example, the MEG II experiment at the Paul Scherrer Institute (PSI) is expected to reach B µ→eγ < 6 × 10 −14 [7,8]. The Mu2e experiment at Fermilab and the COherent Muon to Electron Transition (COMET) experiment at Japan Proton Research Complex (J-PARC) plan to increase their sensitivity reaches by four orders of magnitude, i.e. to the level of B µ→e < O(10 −17 ) [9][10][11].
Compared to the non-hadronic decay µ → eγ, the process for µ → e conversion occurs in nuclei, wherein a muon is trapped to form a muonic atom. The µ → e conversion process has the potential to discern signs of various BSM physics (for early studies see Ref. [12]), since it can be generated by not only photonic dipole operators but also non-photonic contact four-fermion interactions involving two leptons and two quarks, with various Lorentz structures. Thus, estimations of the conversion process involve multiple scales (hadronic, nuclear, atomic) and require careful treatments. In light of the expected improvements in experimental sensitivity, in order to maximize the constraining power (in case of null signal) or the model-diagnosing power (in case of discovery), accurate theoretical predictions within various classes of models are desirable.
The calculation of µ → e conversion rates has a long history, starting with the pioneering work of Weinberg and Feinberg [13], in which the electron wavefunction was taken as a plane wave and muon wavefunction was approximated to a constant. It was later realized that relativistic effects can be important in medium and heavy nuclei [14], and detailed calculation of the conversion rate including relativistic lepton wavefunctions was performed in Refs. [15,16]. The conversion rates in various nuclei and for dipole, scalar, and vector operators have been calculated in [16], where the uncertainty induced by neutron and proton densities in the ground state is also estimated. The discriminating power of various interactions in µ → eγ and µ → e conversion was further assessed in [17], where the uncertainty from scalar form factors was also discussed. The effect of nucleon spin-dependent operators was first discussed in Refs. [18,19] and effective field theory studies including renormalization group evolution, along with their phenomenological implications, have appeared [20][21][22]. Very recently, a nuclear-level effective theory for µ → e conversion was developed [23].
Among the spin-independent operators, the largest theoretical uncertainties arise in the scalar sector. An improved treatment of scalar matrix elements based on arXiv:2203.09547v2 [hep-ph] 24 May 2022 SU(2) chiral perturbation theory (ChPT) was introduced in Ref. [24]. More recently, in the framework of SU (2) ChPT, the impact of next-to-leading order (NLO) nucleon interactions induced by quark-level scalar densities was discussed in Ref. [25]. The NLO contributions involve both single-nucleon scalar form factors and twonucleon interactions, which were reduced in Ref. [25] to an effective one-nucleon interaction by performing an average of the interaction over a Fermi gas model. It was found that, when scalar operators are the dominant sources for the conversion process, the NLO interactions could bring destructive contributions relative to LO contributions and significantly reduce the branching ratio.
Inspired by those previous studies, we revisit the calculation of branching ratios of µ → e conversion in several nuclei including the NLO nucleon interactions induced by scalar quark densities. Our analysis is performed in a model-independent way, i.e. the Standard Model Effective Field Theory (SM-EFT) [26][27][28][29], focusing on a particular class of SM-EFT operators: photonic dipole and scalar four-fermion operators. This setup is well motivated by BSM models with heavy scalar particles such as Two-Higgs Doublet Model and Leptoquark Models. As a byproduct of our analysis, we correct two errors appearing in [25]: (1) we eliminate two un-physical contributions to the overlap integrals (denoted by τ (2,3) ); (2) we provide a corrected expression and corresponding numerical result for the effective one-nucleon interaction resulting from the average over the Fermi gas model. We also develop a new method to include the momentumdependence of the scalar form factor in the overlap integrals. With these results at hand, we assess the NLO contributions to the conversion process and discuss the current hadronic and nuclear uncertainties. We discuss the prediction of this scalar-dominance model for the ratio of µ → eγ over µ → e conversion in 27 Al and for the ratio of conversion rates in different target nuclei. Finally, we apply the analysis to a simple model in which CLFV is mediated by CLFV Yukawa couplings of the SM Higgs, which realizes the current setup.
The paper is organized as follows. In Section II we set up the EFT framework, starting from quark-level interactions and matching to nucleon-level interactions. In Section III we present the results for the overlap integrals and the µ → e conversion rate, including NLO chiral effects. In Section IV we discuss the implication for scalarmediated CLFV first in the EFT setup and subsequently in a model with 'minimal' Higgs-mediated CLFV, i.e. an extension of the Standard Model in which the only new interactions are CLFV Yukawa couplings of the Higgs to leptons. We present our conclusions in Section V and relegate some technical details to the Appendices.

II. EFFECTIVE INTERACTIONS: FROM QUARKS TO NUCLEONS
We assume in this work that the BSM physics responsible for CLFV originates at energies above the electroweak scale. In this case contributions from any BSM physics are captured by effective operators expressed in terms of SM fields, with appropriate couplings that contain information about the underlying model -this is the SM-EFT framework. We restrict our attention to a particular class of SM-EFT operators mediating CLFV transitions, namely the ones mediated by heavy scalar particles, including the SM Higgs itself. As our goal is to assess the uncertainties in scalar-mediated CLFV, we take as starting point below the electroweak scale the following effective Lagrangian (for a complete set of operators see Refs. [16,17]) where F µν is the field strength of the photon, P L,R = (1 ∓ γ 5 )/2 are the chirality projectors, Λ represents the new physics scale, and the Wilson coefficients C Dα , C (q) Sα are dimensionless. With this normalization, the chirality flip in lepton and quark bilinears is accompanied by a muon or a quark mass insertion, respectively. 1 This choice comes without loss of generality and simplifies many intermediate steps in the analysis. Finally, the factors of m q and α s multiplying the quark scalar bilinears and the gluonic operator ensure that the corresponding Wilson coefficients do not run under QCD renormalization. After integrating out the heavy quarks, at the GeV scale the effective Lagrangian takes the form of Eq. (1), with q = u, d, s and [30] ( The scalar quark operators in Eq. (1) induce singleand multi-nucleon momentum-dependent operators at low-energy, which eventually lead to nuclear transitions. The form of one-and two-nucleon operators has been derived to NLO in both SU(3) ChPT [31] and SU (2) ChPT [24,25,32,33], and we will work here in the SU(2) case. Making the following identifications (recall α ∈ {L, R} is a chirality label for the lepton bilinear appearing in the scalar operators) The isoscalar and isovector combinations of scalar Wilson coefficients are given by: The single-nucleon scalar form factor is given by We note that the first order Taylor expansion provides a representation of the full expression accurate to 2% for x ≤ 1 and to 5% for 1 < x ≤ 2. This means that it is quite safe to use the linear term in our nuclear analysis. In Appendix. A, we show how the momentum-transfer expansion corresponds to derivative operators acting on nucleon density functions.
The hadronic scalar current defined in this way carries uncertainties due to both input parameters and higher order terms in the chiral expansion. Higher order terms in the momentum-independent part of the current are effectively resummed by using the physical values of σ πN , σ s , and δm N . The hadronic inputs entering Eqs. (4) are defined as For the sigma term, we use as baseline input the analysis of the Roy-Steiner equations of Ref. [34], namely σ πN = 59.1 (3.5) MeV. This value is in tension with the (currently) more uncertain lattice QCD calculations (see [35] and references therein) which indicate σ πN = 65(13) MeV (with dynamical charm quark [36]) and σ πN = 40(4) MeV (no dynamical charm quark [37][38][39]). However, the recent lattice QCD analysis of Ref. [40] suggests that inclusion of excited state effects reconciles the tension, so we shall use the input from Ref. [34]. For the strange sigma term we use the lattice QCD average [35] σ s = 41(9)MeV (with dynamical charm [41]), while for the slope we will takeσ s = 0.3(2) GeV −1 [42]. For the strong-isospin contribution to the nucleon mass splitting we take the lattice QCD determination δm N = 2.32 (17) MeV from Ref. [43], consistent with the earlier lattice calculation of Ref. [44]. Finally, we take = 0.365(23) from the FLAG average [35]. Higher chiral orders in the momentum dependence of the single-nucleon form factor are expected to be sizable. In fact, a comparison of the NLO heavy baryon ChPT prediction [45] with a recent dispersive determination [42] indicates that the NLO result accounts for about 60% of the dispersive result.
Finally, note that the one-and two-nucleon amplitudes µN (k) → eN (k ) and µN (k 1 )N (k 2 ) → eN (k 1 )N (k 2 ) take the form where ēP α µ denotes the leptonic amplitude and the physically relevant combinations of hadronic scalar currents are which take the form If the gluonic operator is sourced only by integrating out the heavy quarks, one has the relation The NLO two-nucleon contribution can be reduced to an effective single-nucleon operator by averaging the twonucleon operator over a Fermi gas model of the target nucleus. In this approximation, the effect of J (2) α is captured by the shift where k F is the Fermi momentum in the target nucleus and f SI eff is the effective single-nucleon coupling resulting from averaging over the Fermi gas model. Although this procedure was carried out in [25], an error in that calculation resulted in incorrect results for the effective single-nucleon form factors f SI (q,k) and f SD (q,k) which propagated to all values obtained from these form factors. In Appendix C, we present the corrected expressions. The corrections significantly affect the values obtained after momentumaveraging, reducing the overall magnitudes by roughly a factor of two and leading to f SI eff = 0.43 +0.03 −0.22 and f SD eff = 0.43 +0.03 −0.22 . As discussed in Ref. [25], the effective coupling f SI eff obtained through the Fermi gas average is likely an overestimate of the underlying two-nucleon contribution. This expectation is based on a study of nuclear anapole moments [47] where, in addition to the Fermi gas average, nuclear shell model wave functions were used to directly evaluate a variety of two-nucleon currents -none of them the operator of present concern -in the nuclei 133 Cs and 205 Tl. Across all the currents tested, the Fermi gas average tended to overestimate the two-nucleon contribution by 2-3 times compared to the shell model.
To verify that this behavior persists in the present case, we evaluated the two-nucleon operator using shell model wave functions for two of our nuclei of interest, 27 Al and 48 Ti. Details of this calculation are presented in Appendix B. Fully correlated shell model wave functions were generated for 27 Al and 48 Ti using the configuration-interaction code BIGSTICK [48,49] and the USDB 1d [51] interactions, respectively. Harmonic oscillator bases with oscillator parameters b of 1.84 and 1.99 fm for Al and Ti, respectively, were employed in the calculations.
The effective one-body couplings which reproduce the shell model results are given in Table I. We find that the Fermi gas average estimate is ≈ 2 − 3 times larger than the shell model estimate, in good agreement with the anapole study. Although we did not carry out the shell model calculation for the heavy nuclei 197 Au and 208 Pb, given that the anapole study observed the overestimation in 205 Tl, it is likely that a similar result would be found for 197 Au and 208 Pb in our case.
The shell model results may still represent an overestimation of the two-nucleon contribution. The shell model 27 13 [46], Rp (Rn) and a are the parameters of the proton (neutron) density profile, f SI eff,FGA is the value of the spin-independent form factor obtained via the Fermi gas average, and f SI eff,NSM is the value implied by the nuclear shell model calculation (without additional correlation function).
The upper uncertainty of f SI eff,FGA is due to the error incurred in the momentum average overk and the uncertainty in the Fermi momentum (±5 MeV) whereas the lower uncertainty reflects the expectation that the FGA result tends to overestimate the strength of the operator by roughly a factor of 2. See discussion in Appendix B and C for more details. wave functions that we employ are constructed in very soft Hilbert spaces which lack the high-momentum modes necessary to properly resolve the strong repulsion of two nucleons at short distance. We find that the two-nucleon operator is sensitive to the short-range nucleon-nucleon physics. Introducing an ad hoc short-range correlation function [52] in the shell model calculation further reduces the estimated strength of the two-nucleon operator by ≈ 50%. For this reason, the value of f SI eff obtained from the shell model calculation can likely be considered as an upper limit on the strength of the two-nucleon contribution. A complete treatment involving the introduction of effective operators and wave function renormalization to account for the truncated shell model space is beyond the scope of this paper.

III. TRANSITION RATE INCLUDING NLO CORRECTIONS
The rate of the coherent µ → e conversion process depends on the behavior of the bound muon and outgoing electron. The lepton wave functions are obtained by solving the Dirac equation [15,16,53,54], with Here, the energy, potential and mass of the leptons are given by W, V (r) and m. σ are the Pauli matrices,r is a unit vector in the radial direction, l is the orbital angular momentum defined by l = −ir × ∇. We define the wave functions as where µ and κ represent the eigenvalues of the z component of the total angular momentum J z and K, respectively. The two-component spinors χ µ κ are the spinangular functions, with the properties The initial muon state corresponds to the ground state of the muonic atom, implying κ µ = −1 . On the other hand, the outgoing electron has two states of κ e = ±1. Normalization of the bound muon state is defined by Neglecting nuclear recoil, the final state electron carries energy W = m µ − B µ , where B µ is the binding energy of the 1s muonic atom. Its wave function is normalized as Inserting the expressions of the wave functions into the the spherical polar form of the Dirac equation, one can obtain d dr Utilizing the shoot-and-match procedure [55], we solve these coupled equations numerically.
For µ → e conversion, the branching ratio is defined by the conversion-to-capture ratio where A and Z are mass and atomic numbers, respectively. The standard muon capture rates Γ capt ≡ κ capt m 5 µ /v 4 for the nuclei of interest are listed in Table  II. Taking into account all the spin configurations of the initial muon and final electron states, one can express the branching ratio as in terms of dimensionless overlap integrals where the upper index indicates the quantum number κ e . The overlap integrals for the dipole operator are given by For the contributions from the scalar operator, we split the overlap integrals into two contributions, momentumtransfer independent (τ (±1) ρ ) and dependent (τ where C ρ N α and C f N α correspond to the constant and momentum-dependent parts of the nucleon scalar form factor, respectively: The nucleon-level couplings C ρ,f N α are dimensionless, and the minus (plus) sign in the second term of C ρ N α is for proton (neutron). The dimensionless overlap integrals where ρ N (r) is the nucleon density and the function  27 13 Al, 48 22 Ti, 197 79 Au and 208 82 Pb.
In our analysis, we employ the two-parameter Fermi model for the proton and neutron densities the parameters of which are given in [56] for Al, Au, Pb and [57] for Ti. We reproduce these parameters in Table I. The nucleon density profiles for Al, Au, and Pb include a separate determination of the neutron density from measurements of pionic atoms, whereas for Ti we have only the proton density, determined from electron scattering. In this case, we assume R n = R p . It should be noted that, assuming m e = 0, one can obtain relations between two electron states, g  Table III shows the results of the dimensionless overlap integrals τ 27 13 Al, 48 22 Ti, 197 79 Au and 208 82 Pb. We find that the values of the 2 We have numerically checked that the fourth derivative term is roughly two-orders of magnitude smaller than the second derivative. Therefore, we neglect higher-order derivatives in f N (r). 3 The overlap integrals defined in [16] are obtained by scaling τ is computed with the effective one-body operator strength obtained via Fermi gas average, f SI eff,FGA (see Table I.)  momentum-dependent overlap integrals τ f N for Al and Ti are roughly a factor of 2/3 smaller than the corresponding momentum-independent integrals τ ρ N . In heavier nuclei, τ f N is further suppressed relative to τ ρ N .
In the case of non-zero CLFV couplings to u and d quarks, we can compare the LO contributions with those from the NLO interactions by defining the following quantities: The resulting values for the nuclei of interest are presented in Table IV, where input parameters such as σ πN and f SI eff are fixed at their central values. As discussed in [25], the NLO nucleon interactions bring negative relative contributions compared to the LO value, thereby reducing the overall µ → e decay rate. The NLO loop contribution to the amplitude is less than 5% of the LO contribution, while the two-nucleon interactions roughly amount to 10%. If we assume C In the case of the strange quark, the NLO correction arises only from the momentum-dependent term, τ f N . Taking the central value ofσ s , we see that the NLO term could reduce the strange-quark contribution to the CLFV amplitude by 10% in Al and Ti, while the corrections are less than a few % in Au and Pb. Assuming equivalent Wilson coefficients C (0) Sα , the total (LO + NLO) contribution to the CLFV amplitude from the strange quark is reduced by ≈ 30% in Al and Ti and ≈ 20% in Au and Pb, relative to the contribution from u, d quarks.
In contrast to the light quarks, the NLO contributions to the gluonic coupling C Gα are less than 1% relative to the leading term which is dominated by the nucleon mass. Recalling the relation between the heavy quark (q = c, b, t) Wilson coefficients and the gluonic coupling, Eq. (2), we find that the prefactor of the heavy quark Wilson coefficients in Eq. (32a) is ≈ 20% larger than the prefactor of the isoscalar coupling C  Table  I. The lower uncertainty reflects the fact that the central value of f SI eff obtained from the Fermi gas average is likely too large by a factor of two. We find a corresponding increase in the value of τ Overall, the variation of the ratio in the lower plot is small compared to the upper plot. This is because, as discussed in [24], the case has two additional heavyquark contributions, leading to less impacts from the parameters that we currently focus on. We find that the scalar contribution is dominantly affected by the uncertainty in σ s , which roughly amounts to ±7% (3%) for the C d,s,b SR (C q=all SR ) = 0 case. Varying all of these parameters, we see that the deviation from the central value of τ (−1) S,c 4 Since the errors in , δm, and k F are negligible, we do not include them in our analysis.

FIG. 1. Uncertainty budget for the overlap integrals.
We show the dependence of τ in Al on each parameter: σπN , σs, f SI eff (denoted by 2N(f SI eff )) and one-body form factors for (u, d) and s (represented by 1N FF (u, d) and (s)). Only right-handed down-type operators are nonzero in the upper plot, while the lower plot takes all the right-handed operators to be non-vanishing. In both cases, the nonzero Wilson coefficients are assumed to be equal; that is, C is roughly ±10% and ±5% for the upper and lower case, respectively.
The relative importance of the NLO contributions can depend significantly on the underlying CLFV physics. In particular, when CLFV primarily arises from light quark scalar couplings the NLO contribution can in fact be larger than the LO uncertainty. We illustrate this with the following two examples: • If only the two lightest quarks contribute, C , we see that the negative NLO contributions are ≈ 25%, which is consistently larger than the LO uncertainties ≈ ±13%. This is because the two light-quark contributions dominate the conversion process if the scalar operators are generated with m q C (q) This analysis implies two important take-home messages: (1) the overall uncertainty is dominated by the LO amplitude, in particular by the sigma terms; (2) the central value of the NLO corrections could be larger than the uncertainty on the LO term if light quarks (q = u, d) have nonzero LFV couplings, making the analysis of NLO effects phenomenologically relevant. As we find in the next section, the relative significance of the NLO contributions can be diminished in scenarios with large contributions from either gluonic couplings generated by heavy quarks or dipole operators.
It should be noted that in our analysis we have adopted a central value and uncertainty for the two-nucleon contribution based on the Fermi gas model. In light of the nuclear shell model results and their strong dependence on short-range correlations, it is at present impossible to rigorously quantify the uncertainty in the two-nucleon sector. Nonetheless, it is reasonable to regard the values considered here as an upper limit on the relative strength of the two-body contribution (and in turn the NLO contribution).

IV. PHENOMENOLOGY IMPLICATIONS
We next discuss some phenomenological implications of the improved analysis of transition rates.

A. Dipole-scalar dominance model
First, we consider a restricted EFT setup in which only dipole and scalar operators are generated, assuming that the ratio C D /C S is characterized by a real parameter r. As in Ref. [17], we assume C DR = (r/8e)C SR with SR , while the rest of the operators are zero. This scenario may be explicitly realized in some regions of the R-parity conserving SUSY see-saw parameter space [58] (large tan β and relatively low "heavy" Higgs sector) and within R-parity violating SUSY [59][60][61][62]. The nonzero dipole operators generate the µ → eγ process as well, whose branching ratio is simply expressed as Figure 2 shows B µ→e (Al)/B µ→eγ in the upper plot and B µ→e (Ti)/B µ→e (Al) in the lower plot. In this parametrization, the dominant contribution to the µ → e conversion switches from the scalar operator to the dipole one around r ≈ 10 −6 . The band in the upper plot is obtained by taking the 1σ range of the input parameters (σ πN , σ s ,σ s , , δm), the uncertainty in f SI eff from Table I, and a 50% error in the one-body form factor. For example, at r = 10 −7 where the scalar operator dominates B µ→e (Al), the uncertainty in B µ→e (Al)/B µ→eγ corresponds to ±20% which can be understood from the analyses of τ S,c in the previous section. On the other hand, taking the ratio of the conversion process between two nuclei, one can see that the uncertainty becomes negligible as seen in the lower plot. Here, we take the central values of f SI eff for Al and Ti, since we expect that the uncertainty in the effective one-body coupling is correlated across all isotopes. We see a few % differences from the results in [17].
However, it should be noted that the ratio B µ→e (Ti)/B µ→e (Al) is affected by the uncertainty in the overlap integrals due to the neutron densities. As given in Table I, the uncertainty in the neutron density parameter R n determined from pionic atom experiments is roughly 6% in 27 Al. This propagates to a 5% error in the neutron overlap integrals for Al. For Ti, in absence of a direct measurement of the neutron density, we employ the same density profile as the measured proton density in Ti. Therefore one would expect a significantly larger uncertainty on the Ti overlap integrals stemming from the neutron density; for example, in the nearby nucleus 56 Fe -for which a measurement of the neutron density is available -the difference between using a neutron density measured in pionic atoms compared to assuming identical density profiles for protons and neutrons results in ≈ 7% change in the neutron overlap integrals. Accounting for this discrepancy, as well as the uncertainty in the neutron profile parameters, we assume an 8% error on the Ti neutron overlap integrals below in Eq. (46). None of these overlap integral uncertainties are reflected in Figure 2, although we expect these errors to be relevant in the scalar-dominated region where r 10 −5 .

B. CLFV Yukawa couplings
Finally, we apply our analysis to the Higgs-mediated CLFV model, where the following Yukawa interactions generate µ → e transitions These Yukawa interactions induce µ → eγ at one-and two-loop level as discussed in [63]. 5 On the other hand, the scalar operators arise from a tree-level process mediated by the Higgs particle. The Wilson coefficients are given by where the Higgs vacuum expectation value v = 246 GeV and the Higgs mass m h = 125 GeV. Note that C (q) Sα in this model becomes independent of the label q. Having the parametrization of C D = r/(8e)C S that we employ in the previous section, we obtain r = 3.4 × 10 −6 in this model. Figure 3 depicts bounds on the CLFV Yukawa couplings Y eµ and Y µe . The gray line represents the upper limit on the two couplings from B µ→e (Au)< 7 × 10 −13 . The bound originating from µ → eγ is presented by the orange region, which corresponds to Y eµ , Y µe 10 −6 . As shown by the black dashed line, the next-generation µ → e experiments will provide a sensitivity to Y eµ and Y µe that is stronger than MEG II [7] and ten times stronger than current limits. 6 The uncertainty resulting from the hadronic input parameters and NLO interactions is not visible on the scale of the plot.
The plurality of probes, namely µ → eγ and µ → e conversion in possibly more than one target nucleus, provides an opportunity to test underlying new physics CLFV mechanisms. The minimal Higgs-mediated CLFV scenario considered here produces at low-energy a specific combination of scalar and dipole operators, which leads to the following pattern of branching ratios: 7 Here, we assign a 5(8)% error to neutron overlap integrals τ ρn and τ fn in Al (Ti). While the dominant uncertainty in B µ→e (Al)/B µ→eγ arises from hadronic input parameters, the ratio B µ→e (Ti)/B µ→e (Al) is primarily affected by the uncertainties in neutron densities of Al and Ti whereas those of the hadronic parameters are negligible. Our predictions with quantified uncertainties offer a clean path towards testing the Higgs-mediated CLFV scenario, in case of discovery in the next generation experiments.

V. CONCLUSIONS
In this paper we have studied µ → e conversion in nuclei within a class of models in which lepton flavor violation is induced dominantly by the exchange of scalar particles, such as the SM Higgs. We have estimated the µ → e conversion rate in several nuclei of experimental interest, including NLO effects in the one-and two-nucleon interactions in chiral effective theory. The one-nucleon effects involve the momentum dependent scalar form factor, and we have developed an efficient method to take into account the momentum dependence in the overlap integrals. The two-nucleon terms are evaluated both by reducing them to one-body terms via an average over a Fermi gas model (as done in Ref. [25]) and within the nuclear shell model. The two approaches provide a way to estimate the uncertainty associated with the two-body NLO contribution. In the process, we correct the result of Ref. [25], finding a smaller effective one-body interaction and hence a smaller impact on the µ → e conversion rates.
For the light-quark scalar operators, the NLO corrections interfere destructively with the LO terms. At the amplitude level, the NLO form-factor contribution is at the 5% level, while the two-nucleon interactions amount to about 10%, with a total NLO impact on the decay rates at the 20-30% level. We provide an uncertainty budget for the conversion amplitude at LO and NLO and find that the overall uncertainty is dominated by the LO amplitude, in particular by the sigma terms. Importantly, we find that the central value of the NLO corrections could be larger than the uncertainty on the LO term, making the analysis of NLO effects phenomenologically relevant.
We studied the implications of our results for testing scalar-mediated CLFV processes. For the phenomenologically interesting case of Higgs-mediated CLFV, in which both scalar and dipole operators appear at low-energy, we have shown that: (i) in the next generation experiments µ → e conversion will have stronger sensitivity than µ → eγ to the CLFV Yukawa couplings Y µe,eµ ; (ii) ratios of branching ratios such as B µ→e (Al)/B µ→eγ and B µ→e (Al)/B µ→e (Ti) can be predicted with quantified uncertainties. In particular, B µ→e (Al)/B µ→eγ is affected by NLO corrections at the same level as LO uncertainties. These ratios offer a clean path towards testing the Higgs-mediated CLFV scenario in case of discovery in the next generation experiments. Let us consider µ → e conversion mediated by a generic scalar operator of the form: In the more commonly studied case of lepton scattering (electron or neutrino scattering), one typically uses the mode expansion of the leptonic fields in terms of plane waves. Later on, we will consider this case as a consistency check on the formalism we use. However, it turns out that using plane waves is not the most convenient choice for the problem of µ → e conversion. In this case it is more convenient to perform the mode expansion of the leptonic field in terms of solutions to the Dirac equation in the Coulomb field of the nucleus. The muon initial state corresponds to the 1s hydrogen-like state, while the electron final state corresponds to an outgoing wave in the continuum, with energy given by the muon mass minus the binding energy B µ . The scalar Hamiltonian at low energy takes the form where C L,R ≡ 1/Λ 2 L,R are dimensionful Wilson coefficients andŜ(x) is the hadronic realization of the quark scalar density. In chiral EFT this operator contains terms likeN N , π + π − , etc. Its matrix element between free nucleon states takes the form where g S (q 2 ) is the nucleon scalar form factor (for simplicity we do not display isospin indices here and throughout the discussion). In the regime we are working, we could further Taylor expand the scalar form factor as follows: g S (q 2 ) = g 0 + g 2 q 2 + g 4 (q 2 ) 2 + ... .

(A4)
Ultimately, we wish to have a non-relativistic realization of the scalar density operator, to be inserted between nuclear many-body wavefunctions. We denote this object byŜ(x) and we determine its form by requiring that the matrix elements within free one-nucleon states are the same in the two representations. Explicitly, Using Eq. (A3) and the Taylor expansion of the form factor, one can verify that the desired coordinate-space representation of the scalar density involves the delta function and its derivativeŝ .. , (A6) where the summation runs over nucleons.
With this result at hand, the conversion amplitude takes the form where |Ψ 0 is the nuclear many-body ground state wavefunction. The leptonic matrix element is expressed in terms of the solutions of the Dirac equation ψ (e),(µ) (x): The electron wavefunction has an oscillatory behavior with frequency set by m µ − B µ , but it is not a plane wave. For the nuclear part of the matrix element one has and using the definition of one-body density 8 8 Note that with this definition the proton and neutron densities are normalized as ∞ 0 dr 4πr 2 ρp(r) = Z and ∞ 0 dr 4πr 2 ρn(r) = A − Z.
Eq. (A12) involves derivatives acting on the nucleon densities in the nuclear ground state and differs from the standard results encountered in the analysis of leptonnucleus scattering. In the case of lepton scattering, the external leptonic wavefunctions are plane waves, and the analogue of Eq. (A12) becomes (here q = p f − p i is the leptonic momentum transfer) After integration by parts one recovers the familiar factorized form The factorized form of (A14) does not directly apply to µ → e conversion because the leptonic wavefuctions ψ (µ) (x) and ψ (e) (x) are not plane waves. However, one can perform a Fourier decomposition of ψ (µ) (x) and ψ (e) (x) and apply (A14) to each term in the double Fourier expansion in p e and p µ , leading to e f |Ĥ L,R |µ (1s) = C L,R dp e dp µū (p e )P L,R u(p µ ) × g S (p e − p µ ) 2 d 3 x e i(pe−pµ)·x ρ (1) (x) .
This form is equivalent to (A12) but in practice, as long as g S (q 2 ) can be approximated by a few terms in its series expansion around q 2 = 0, Eq. (A12) is computationally more convenient than (A15).

(B8)
The resulting operator J=0,M =0 (q T ) τ 1 · τ 2 C These functions depend not only on the magnitude of the dimensionless momentum transferq and average momentum k, but on their relative angle. Fortunately, for the physically relevant values of these momenta, f SI and f SD do not vary significantly over the range of possible angular values. Therefore we may replace each function by its angular average. The angle-averaged functions then depend only on the magnitude of the momentum transferq and the average momentumk. For a given nucleus, the three-momentum transfer in coherent µ → e conversion satisfies where M T is the mass of the target nucleus and E bind µ is the (positive) binding energy of the captured muon. Fixing the value of q T for a given nucleus, f SI and f SD become functions of the dimensionless average nucleon momentumk, as shown in Fig. 4. In order to recover a local one-body effective operator, we replace these slowly-varying functions ofk by a constant, weighting our average by the single-nucleon momentum probability distribution. The resulting momentum-averaged values f SI eff and f SD eff are shown in Fig. 4. Table I reports the values of f SI eff for the four nuclei of interest as well as the physical parameter values employed in each calculation.