$\Lambda_c \to N$ form factors from lattice QCD and phenomenology of $\Lambda_c \to n \ell^+ \nu_\ell$ and $\Lambda_c \to p \mu^+ \mu^-$ decays

A lattice QCD determination of the $\Lambda_c \to N$ vector, axial vector, and tensor form factors is reported. The calculation was performed with $2+1$ flavors of domain wall fermions at lattice spacings of $a\approx 0.11\:{\rm fm},\:0.085\:{\rm fm}$ and pion masses in the range $230\:{\rm MeV} \lesssim m_\pi \lesssim 350$ MeV. The form factors are extrapolated to the continuum limit and the physical pion mass using modified $z$ expansions. The rates of the charged-current decays $\Lambda_c \to n\, e^+ \nu_e$ and $\Lambda_c \to n\, \mu^+ \nu_\mu$ are predicted to be $\left( 0.405 \pm 0.016_{\,\rm stat} \pm 0.020_{\,\rm syst} \right)|V_{cd}|^2 \:{\rm ps}^{-1}$ and $\left( 0.396 \pm 0.016_{\,\rm stat} \pm 0.020_{\,\rm syst} \right)|V_{cd}|^2 \:{\rm ps}^{-1}$, respectively. The phenomenology of the rare charm decay $\Lambda_c \to p\, \mu^+ \mu^-$ is also studied. The differential branching fraction, the fraction of longitudinally polarized dimuons, and the forward-backward asymmetry are calculated in the Standard Model and in an illustrative new-physics scenario.


A. Lattice parameters and correlation functions
This calculation uses the Iwasaki action [16] for the gluons, the Shamir-type domain-wall action [17][18][19] for the u, d, and s quarks, and an anisotropic clover action for the c quark, with parameters tuned in Ref. [20] to reproduce the correct charmonium mass and relativistic dispersion relation. The gauge field ensembles were generated by the RBC and UKQCD Collaborations and are described in detail in Ref. [21]. This work is based on the same six sets of light-quark domain-wall propagators as Ref. [22]; the parameters of these sets and the resulting pion masses are listed in Table I. The parameters of the charm-quark action are given in the first four rows of Table II. The Λ c and nucleon two-point functions and the Λ c → N three-point functions were computed analogously to Ref. [22], with the bottom quark replaced by the charm quark. The c → u currents were renormalized using the "mostly nonperturbative" method [24,25], in which the bulk of the renormalization is absorbed by an overall factor of Z  [22], with residual matching factors and O(a)-improvement coefficients computed to one loop by Christoph Lehner [26,27] and listed in Table II. The full one-loop O(a) improvement was performed for the vector and axial vector currents for all source-sink separations in the three-point functions (instead of just a subset of separations as in Ref. [22]). For the c → u tensor currents, the one-loop calculation was not available, and the perturbative coefficients were evaluated at tree-level as in Ref. [28]. That is, the tensor currents were written as  [21] and light quark propagators. The C14/C24 and F23 data sets are based on the same gauge field ensembles as the C54 and F43 data sets, respectively, and differ only in the valence quark mass used for the propagators. The lattice spacing values given here were determined using the Υ(2S) − Υ(1S) splitting in Ref. [23].  (26) TABLE II. Charm-quark action parameters [20] and c → u current matching and O(a) improvement parameters [26,27]. The nonperturbative factors Z were determined in Refs. [30] and [22], respectively. The uncertainties of the residual matching factors and O(a)-improvement coefficients were estimated as ( where d 1 is the mean-field-improved heavy-quark "field rotation" coefficient [29], whose values are also listed in Table  II. The missing one-loop corrections result in larger systematic uncertainties for the tensor form factors, as discussed in Sec. III B. The three-point functions were computed for all source-sink separations in the ranges t/a = 4...15 (C14, C24, C54 data sets), t/a = 5...15 (F43 data set), and t/a = 5...17 (F63 data set). The Λ c momentum, p, was set to zero, and all nucleon momenta p with 1 2π L 2 ≤ |p | 2 ≤ 5 2π L 2 were used. From the three-point and twopoint correlation functions, the "ratios" , and R h+ (|p |, t), defined as in Eqs. (52-54), (58-60) of Ref. [22] and Eqs. (27)(28)(29)(30) of Ref. [28] (with the appropriate replacements of initial and final baryons), were computed. These ratios are equal to the form factors f ⊥ , f + , ... for the given momentum and lattice parameters, up to excited-state contributions that vanish exponentially as t → ∞. The ratios also depend on the baryon masses, which were obtained by fitting the two-point functions from the same data sets and are listed in Table III (the table also contains the D meson masses, which are needed at a later stage in the analysis).
The ground-state form factors f (|p |) were then extracted by performing correlated fits of the ratios of the form which include the leading excited-state contributions. Examples of the fits are shown in Fig. 1. At a given momentum, the fits were performed jointly for all six data sets, and jointly for the form factors associated with a given type of current, with constraints as explained in Ref. [22]. The constraints limit the variation of the energy gap parameters across data sets to physically reasonable values, and enforce that the relations between the form factors in the helicity basis and the "Weinberg basis" are preserved by the extrapolation [22,28].
The values of the start time slices t min were chosen to achieve The data shown here are from the C54 set at |p | 2 = 2 2π L 2 . At a given momentum, the fits were performed jointly for all six data sets with constraints as explained in Ref. [22]. Some data points at the shortest source-sink separations, plotted with open symbols, were excluded from the fits to achieve χ 2 /d.o.f. < ∼ 1. The curves going through the data points show the fit functions of the form  To estimate the remaining systematic uncertainties associated with the choices of fit ranges, additional fits were performed in which all t min values were increased simultaneously by one unit. As in Refs. [22,28,31], the systematic uncertainty in f (|p |) for a given momentum |p | and data set was estimated as the larger of the following two: i) the shift in f (|p |) at the given momentum, and ii) the average of the shifts in f (|p |) over all momenta. These systematic uncertainties were added to the statistical uncertainties in quadrature. The results for f (|p |) with the combined uncertainties are listed in Table VIII in the Appendix, and are also shown as the horizontal bands in Fig. 1 and as the data points in Figs. 2-4.

B. Chiral and continuum extrapolations of the form factors
Following the extraction of the form factor values for each lattice data set and for each discrete momentum, global fits of the form factor shape and of the dependence on the lattice spacing and quark masses were performed using modified BCL z-expansions [32]. In the physical limit (a = 0, m π = m π,phys ), the fit functions reduce to the form where the expansion variable is defined as This definition maps the complex q 2 plane, cut along the real axis for q 2 ≥ t + , onto the disk |z| < 1. Here, t + is set equal to the threshold of D π two-particle states, The parameter t 0 determines which value of q 2 gets mapped to z = 0; in this work, is used. Furthermore, in Eq. (10), the lowest poles are factored out before the z expansion. The quantum numbers and masses of the D mesons producing these poles in the different form factors are given in Table IV.
To fit the lattice data, Eq. (10) was augmented with additional parameters to describe the dependence on the lattice spacing and on m 2 π (which serves as a proxy for the the u/d quark mass). As in Refs. [22,28], two independent fits were performed: a "nominal" fit, which provides the central values and statistical uncertainties of the form factors, and a "higher-order (HO)" fit that is used to estimate systematic uncertainties. In this work, the functions used for the nominal fit were with free parameters Here, the scales Λ χ = 4πf π with f π = 132 MeV and Λ had = 300 MeV were introduced to make all parameters dimensionless. The momentum transfers in lattice units, a 2 q 2 , were evaluated using the lattice results for the baryon masses from each data set, and their uncertainties and correlations were taken into account. In addition, the pole masses were rewritten as am f pole = am D + a∆ f , where ∆ f are the mass differences relative to the pseudoscalar D meson. These mass differences were fixed to their physical values according to Table  IV, while the lattice results from each data set were used to evaluate am D .
The systematic uncertainties associated with the choices of t min in the extractions of the lattice form factors from the ratios of correlation functions (cf. Sec. III A) were considered to be part of the statistical uncertainties in the z-expansion fits discussed here, and were therefore included in the data covariance matrix for both the nominal and higher-order fits. Given that the procedure for estimating these uncertainties was based on the magnitudes in the shifts, and conservatively used the larger of two choices, there is no obviously correct way of evaluating their correlations. These systematic uncertainties were therefore added to the diagonal elements of the covariance matrices only. As a result, the z-expansion fits have rather low values of χ 2 /d.o.f. (0.20 for the nominal fit and 0.19 for the higher-order fit, with d.o.f. = 235).
The functions used for the higher-order fit were The quantum numbers and masses of the D mesons producing poles in the different form factors [34]. To evaluate t+, the values mD = 1.870 GeV and mπ = 135 MeV should be used.
The nominal fit already provides a good description of the lattice data, and the additional terms in the higher-order fit are used only to estimate systematic uncertainties in a Bayesian approach. To this end, the additional parameters were constrained with Gaussian priors to be no larger than natural-sized. The priors for the parameters , and j f were chosen as in Ref. [28], while the a f 3 were constrained to be 0 ± 30 as in Ref. [31]. In the higher-order fit, additional sources of systematic uncertainties were simultaneously incorporated as follows: 1. When generating the bootstrap samples for the ratios , the residual matching factors and O(a)-improvement coefficients were drawn from Gaussian random distributions with central values and widths according to Table II. 2. To incorporate the systematic uncertainty in the tensor form factors due to the use of the tree-level values ρ T µν = 1 for the residual matching factors, nuisance parameters ρ T µν multiplying these form factors with Gaussian priors 1 ± σ ρ T µν were introduced in the fit. For the b → s currents in Ref. [28], σ ρ T µν was estimated to be equal to 2 times the maximum value of |ρ V µ − 1|, |ρ A µ − 1|, which was 0.05316, comparable in magnitude to actual one-loop results for |ρ T µν − 1| obtained for staggered light quarks in Ref. [33]. For the c → u currents in this work, ρ V µ and ρ A µ are much closer to 1 (see Table II), and the analogous procedure would yield σ ρ T µν = 0.0105. Given that the tensor current is scale-dependent and may exhibit qualitatively different behavior in a matching calculation (compared to the vector and axial vector currents), this appears to be too aggressive. The uncertainty estimate was therefore increased to σ ρ T µν = 0.05, which is approximately 10 times the maximum value of |ρ V µ − 1|, |ρ A µ − 1|. This estimate of the uncertainty (and the numerical values of the tensor form factors) should be understood as corresponding to the scale µ = m c .
3. The finite-volume errors in the Λ c → N form factors are expected to be approximately equal to those in the Λ b → N form factors, which were estimated to be 3% for the parameters of this calculation [22]. The missing isospin breaking effects are expected to be of order O((m d − m u )/Λ QCD ) ≈ 0.5% and O(α e.m. ) ≈ 0.7%. The uncertainties from these sources were added to the data correlation matrix used in the fit.
4. The lattice spacings and pion masses of the different data sets were promoted to fit parameters, with Gaussian priors chosen according to their known central values and uncertainties.
In the physical limit, the nominal and higher-order fits reduce to the form given in Eq. (10), with n max = 2 and n max = 3, respectively. The solid curves in Figs. 2, 3, and 4 with shaded error bands show the form factors in the physical limit. The results for the relevant parameters are given in Table V. Files containing the parameter values and the full covariance matrices are provided as supplemental material [35]. The systematic uncertainty of any quantity depending on the form factors is estimated as where O and O HO are the central values obtained from the nominal and higher-order parameterizations, and σ O and σ O,HO are the uncertainties propagated using the covariance matrices given in the supplemental material for the nominal and higher-order fit parameters. See also Ref. [28] for further explanations of this procedure. A breakdown of the form factor systematic uncertainties into individual sources is shown in Fig. 5. This was obtained by performing additional fits to the lattice data where each one of the above modifications to the fit functions or data covariance matrix was done individually, and then comparing each of these fits to the nominal fit as in Eq. (16).  Table I), along with the modified z-expansion fits evaluated at the lattice parameters (dashed lines at a ≈ 0.11 fm and dotted lines at a ≈ 0.085 fm) and in the physical limit (solid lines, with statistical and total uncertainties indicated by the inner and outer bands).  are common to the form factors g ⊥ , g+ and h ⊥ , h+, respectively. Files containing the parameter values with more digits and the full covariance matrices are provided as supplemental material. IV. PHENOMENOLOGY OF Λc → n + ν In terms of the helicity form factors, the Λ c → n + ν differential decay rate in the Standard Model reads Evaluating this expression for = e and = µ using the form factor results described in the previous section gives the results shown in Fig. 6. The q 2 -integrated rates are Γ(Λ c → n e + ν e ) |V cd | 2 = (0.405 ± 0.016 stat ± 0.020 syst ) ps −1 , with statistical and systematic uncertainties from the form factors. Using |V cd | = 0.22497(67) from UTFit [36] and the Λ c lifetime τ Λc = 0.200(6) ps from the Particle Data Group [34] yields the branching fractions where the uncertainties labeled "LQCD" are the total uncertainties resulting from the form factors. Previous predictions for Γ(Λ c → n e + ν e )/|V cd | 2 and Γ(Λ c → n µ + ν µ )/|V cd | 2 are summarized in Table VI. All references included there predicted decay rates somewhat lower 1 than the lattice QCD results (18) and (19). While Refs. [1,2,5,7] and [8] estimated the Λ c → n form factors using quark models and light-cone sum rules, respectively, Ref.
Λ c → n µ + ν µ FIG. 6. Predictions for the Λc → n e + νe and Λc → n µ + νµ differential decay rates in the Standard Model, without the factor of |V cd | 2 . The inner and outer bands show the statistical and total uncertainty originating from the form factors.  [2] correspond to nonrelativistic and semirelativistic kinetic terms in the model. Reference [6] predicted the branching fraction, which was converted here using |V cd | = 0.22497(67) [36] and τΛ c = 0.200(6) ps [34]. Reference [8] gave a value for Γ, which was converted using the value of |V cd | given there.

Reference
Method symmetry. The low rate obtained in that way can likely be attributed in part to the assumption m n = m Λ , which results in an underestimated phase space.
The authors of Ref. [4] also calculated the Λ c → N form factors using QCD light-cone sum rules, but only at q 2 = 0. As shown in Table VII, the results are consistent with the lattice QCD determination within 1 to 2 σ.
The theory of exclusive c → u + − transitions in the Standard Model is complicated by the dominance of nonlocal hadronic matrix elements that cannot easily be calculated in lattice QCD. The form factors computed in this work only describe local matrix elements of the form p|ūΓc|Λ c . In the following, two different approaches for expressing the Λ c → p µ + µ − observables in terms of these form factors will be considered: a perturbative calculation of effective Wilson coefficients at next-to-next-to-leading order (Sec. V A), and a phenomenological Breit-Wigner model for the contributions of intermediate φ, ω, and ρ 0 resonances (Sec. V B). The resulting predictions for the Λ c → p µ + µ − differential branching fraction and angular distribution are given in Sec. V C. In the effective-field-theory description, possible heavy new physics beyond the Standard Model contributes only via the local matrix elements given by the form factors. In the case of Λ c → p µ + µ − (but not, for example, in lepton-flavor-violating modes such as Λ c → p e + µ − ), such contributions still interfere with Standard-Model contributions. Nevertheless, the Λ c → p µ + µ − forward-backward asymmetry is nonzero only in the presence of new physics, and therefore provides a clean test of the Standard Model.

A. Standard-Model Wilson coefficients in perturbation theory
The c → u + − effective weak Lagrangian at µ < m b , after integrating out the b quark, has the form where P 1,...,6 are four-quark operators, P 7 and P 8 are electromagnetic and gluonic dipole operators, and P 9 , P 10 are semileptonic operators [9,38]. Following Refs. [9,39], the matrix elements are written as with rescaled electromagnetic dipole and semileptonic operators and q 2 -dependent effective Wilson coefficients which include the perturbative matrix elements of the four-quark and gluonic dipole operators. The Wilson coefficient C 10 is zero in the Standard Model due to CKM unitarity, since all down-type quarks are treated as massless at µ = m W and C 10 does not mix under renormalization [38,39]. At next-to-next-to-leading order, including the recently derived two-loop QCD matrix elements of P for arbitrary momentum transfer [39,40], the effective Wilson coefficients are given by and C eff(q) 7 2,q (m 2 c , q 2 ) + F where the Wilson coefficientsC i are expanded in the strong coupling as and the notationC The functions L(m 2 , q 2 ), L(0, q 2 ), F 8 (m 2 c , q 2 ), and F 8 (m 2 c , q 2 ) can be found in Appendix B of Ref. [9]. The functions F 1,q (m 2 c , q 2 ), F 2,q (m 2 c , q 2 ), and F 1,q (m 2 c , q 2 ) = − 1 6 F 2,q (m 2 c , q 2 ) are given in Ref. [40], and were evaluated here using the files fit F*.dat provided as supplemental material in Ref. [40]. The values ofC i , and αs 4π 2C (2) i at µ = m c were taken from Table 2.2 of Ref. [39]. For the purpose of estimating the perturbative uncertainties in the observables, the values of these coefficients were provided by Stefan de Boer additionally for µ = √ 2 m c and µ = m c / √ 2 [41]. The low-energy value of α e was used, and the strong coupling α s (µ) for four flavors was evaluated using the RunDec package [42]. The quark masses were set to m MS c = 1.27 GeV, m pole c = 1.47 GeV, m pole s = 0.13 GeV as in Ref. [40]. The CKM matrix elements were taken from UTFit [36].

B. Breit-Wigner model of resonant contributions
Similarly to Ref. [9], the contributions from intermediate φ, ω, and ρ 0 resonances were modeled using an effective Wilson coefficient C R 9 (q 2 ) given by and setting C 7 (q 2 ) = 0. Here, the relative magnitude and sign between the ρ 0 and ω amplitudes were fixed as for D + → π + + − in Ref. [43]. In the case of Λ c → p + − , the quark flow diagrams are different, but those diagrams that are expected to dominate yield the same relation. The resonance masses and widths were taken from the Particle Data Group [34]. To determine the couplings a ω and a φ , the Λ c → p µ + µ − branching fraction was computed using only C R 9 (q 2 ) and keeping only the ω or φ contribution, and demanding that where the right-hand side was evaluated using the following experimental inputs: = 0.23 ± 0.08 ± 0.03 [12].
The phases δ ω and δ φ were varied independently in the ranges 0 to 2π when calculating the Λ c → p µ + µ − observables presented in the next section.
C. Results for the Λc → p µ + µ − observables The two-fold differential decay rate of Λ c → p + − with unpolarized Λ c can be written as where θ is the angle of the + in the dilepton rest frame with respect to the direction of flight of the dilepton system in the Λ c rest frame, and the coefficients K 1ss , K 1cc , and K 1c are functions of q 2 . The q 2 -differential decay rate is obtained by integrating over cos θ , The fraction of longitudinally polarized dimuons and the forward-backward asymmetry are defined as and For the effective Lagrangian (23) with operators Q 7 , Q 9 , Q 10 , the expressions for K 1ss , K 1cc , and K 1c in terms of the form factors and the Wilson coefficients C i (q 2 ) can be obtained from Ref. [44] (in the approximation m = 0) or Ref. [45] (for m = 0; used here). These references consider the similar process Λ b → Λ + − ; to adapt the equations to Λ c → p + − , the factor of |V tb V * ts | 2 needs to be removed and the masses need to replaced appropriately. The predictions for dB/dq 2 = τ Λc dΓ/dq 2 , F L , and A F B for Λ c → p µ + µ − , using either the perturbative Wilson coefficients (28) and (29) or the resonant model (32), are shown in the left panels of Fig. 7. In the right panels, an example new-physics contribution of C NP 9 = −0.6, C NP 10 = 0.6 was added to the Wilson coefficients to illustrate the effect on the observables. A contribution of this magnitude is not yet excluded by experimental measurements of rare charm meson decays [9]. Also shown in Fig. 7 is the LHCb upper limit of B(Λ c → p µ + µ − ) < 7.7 × 10 −8 (at 90% confidence level) in the q 2 region excluding ±40 MeV intervals around m ω and m φ [12]. The LHCb upper limit on B(Λ c → p µ + µ − ) also does not yet exclude the new-physics scenario considered here, but comes close.
The error bands of the nonresonant SM predictions are dominated by the perturbative uncertainty, which was estimated by computing the changes in the observables when varying the renormalization scale from µ = m c to µ = √ 2 m c and µ = m c / √ 2. While doing this scale variation, the renormalization-group running of the tensor form factors was included for consistency, by multiplying these form factors with where γ (0) T = 2 C F = 8/3 is the anomalous dimension of the tensor current [46] and β 0 = (11 N c − 2 N f )/3 = 25/3 is the leading-order coefficient of the QCD beta function for 4 active flavors. The error bands of the predictions using the resonant model (32) are dominated by the phase uncertainty, which was estimated by independently varying δ ω and δ φ in the ranges 0 to 2π and showing the resulting ranges of the observables, and the uncertainties in the couplings a ω and a φ as given in Eqs. (37) and (38).
The branching fractions integrated over the entire q 2 range are found to be where, for the perturbative SM prediction, the first uncertainty given is the form factor uncertainty from the lattice calculation, and the second uncertainty is the perturbative uncertainty estimated by varying the renormalization scale. As can be seen in Fig. 7, the low-q 2 region gives most of the perturbative SM contribution. The value for B(Λ c → p µ + µ − ) Perturbative SM calculated here is approximately 10 3 times higher than that obtained in Refs. [47,48]. While [47] does not give a reference for the Wilson coefficients, [48] reportedly uses Wilson coefficients from Ref. [49]. The Wilson coefficients from Ref. [49] actually tend to give higher branching ratios than the Wilson coefficients employed here [38][39][40][41], and the very small branching fractions obtained in [47,48] are puzzling. Note that Refs. [47,48] write the matrix elements of the tensor current in terms of six form factors, of which only four are independent. The numerical parameterizations given there for the six tensor form factors violate the exact kinematical relations between these form factors.
The Λ c → p µ + µ − forward-backward asymmetry, shown at the bottom of Fig. 7, vanishes in the Standard Model because it contains an overall factor of C 10 . New physics giving a nonzero C 10 would produce a nonzero forwardbackward asymmetry, as shown in the bottom-right panel of Fig. 7. While the actual values of A F B strongly depend on the details of the resonant contributions to C 9 , this observable still provides a clean null test of the Standard Model. = −0.6, C NP 10 = 0.6 (right panels). Also indicated is the LHCb upper limit (at 90% confidence level) on B(Λc → p µ + µ − ) in the q 2 region excluding ±40 MeV intervals around mω and m φ [12]. The blue dashed curves use the perturbative results (28) and (29) for the Standard-Model contribution, with error bands including the perturbative scale uncertainties and the form factor uncertainties. The orange solid curves instead use the resonant Breit-Wigner model (32) for the Standard-Model contribution, with error bands including the changes in the observables under independent variations of the phases δω and δ φ in the ranges 0 to 2π, as well as the uncertainties in the couplings aω and a φ and the form factor uncertainties. The curves show the averages over all values of phases.

VI. CONCLUSIONS
In this paper, a precise lattice QCD determination of the Λ c → N (N = n, p) vector, axial vector, and tensor form factors was reported. The results provide Standard-Model predictions for the Λ c → n e + ν e and Λ c → n µ + ν µ decay rates with an uncertainty of 6.4%. The rates calculated here for these decays are higher than those predicted in Refs. [1,2,[5][6][7][8] using quark models, sum rules, or SU (3) symmetry, by factors ranging from 1.4 to 2.
The Λ c → p µ + µ − observables in the Standard Model are dominated by long-distance contributions from nonlocal matrix elements, whose treatment is still very unsatisfactory. However, the forward-backward asymmetry is nonzero only in the presence of new physics, and a measurement would provide a clean test of the Standard Model. More detailed phenomenological studies, including other observables such as CP asymmetries and lepton-flavor-violating decay modes, are warranted.  [50] was used for the lattice calculations.

NOTE ADDED
In this version, I corrected an error in the lepton-mass contributions to F L in Fig. 7. I thank Gudrun Hiller for bringing this to my attention.