Isotope dependence of muon decay in orbit

The decay $\mu^-\to e^- \bar{\nu}_e\nu_\mu$ of a muon that is bound to a nucleus poses an unavoidable background for experiments such as Mu2e, COMET, and DeeMe searching for lepton-flavor-violating $\mu^-\to e^-$ conversion and thus requires a precise understanding. We calculate the electron-energy spectra of muon decays in orbit near the electron energy endpoint for all stable isotopes and estimate uncertainties due to the nuclear charge distribution. Our results enable background studies of mixed or enriched target materials that are necessary for the next generation of $\mu^-\to e^-$ conversion experiments.


I. INTRODUCTION
The search for lepton flavor violation (LFV) is a sensitive probe of physics beyond the Standard Model [1]. Rare muon decays are a particularly interesting avenue due to their long lifetime and the relative ease to produce them in large numbers. Particularly clean processes are the decays µ → eγ and µ → 3e as well as the coherent µ − → e − conversion in the field of a nucleus, µ − +(A, Z) → e − +(A, Z), the latter being especially sensitive to effective LFV operators involving quarks [2,3]. An immense jump in sensitivity to µ − → e − conversion in nuclei is expected in the near future with experiments such as DeeMe [4,5], COMET [6,7], and Mu2e [8,9]. Since the µ − → e − conversion rate depends on the target nucleus (notably its charge Z), it is possible to differentiate between different LFV effective operators using complementary target nuclei in the event of a positive signal observation [3,[10][11][12][13][14]. DeeMe will use graphite or silicon carbide targets while COMET and Mu2e study µ − → e − conversion in aluminium. Previous experiments set limits using copper [15], sulfur [16], lead [17], titanium [18], and gold [19] targets. µ − → e − conversion of bound muons produces approximately monoenergetic electrons with energy where E b α 2 Z 2 m µ /2 and E recoil = (m µ − E b ) 2 /(2m N ) are small corrections due to the muon's binding and nuclear recoil, respectively, in a nucleus with charge Z and mass m N . While electron energies ∼ m µ are seemingly far away from the typical electron energies of the competing µ → eνν decay in orbit (DIO), the presence of a heavy nucleus ensures that this Standard-Model decay distribution has a tail up to E end , providing an irreducible background. Precise calculations are necessary to predict and understand this background DIO spectrum near the electron-energy endpoint in order to chose the optimal signal window. Such calculations exist [20][21][22][23][24][25], most importantly for aluminium, but are lacking sufficient precision for other target nuclei.
Here we set out to provide the DIO spectrum near the electron endpoint -as relevant for µ − → e − conversion experiments -for all stable nuclei, including isotopes. We use the most up-to-date data for the necessary nuclear charge distributions and include the dominant radiative corrections due to the soft photon bremsstrahlung.

II. MUON DECAY IN ORBIT
Calculating the µ → eνν decay of a bound muon requires solving Dirac equation for both the bound muon and the outgoing electron in the electric field of the nucleus [26][27][28][29][30][31]. For spherically symmetric charge distribution ρ, the electric potential V is given by with the Heaviside θ function. The charge distribution is normalized as follows: For a given charge distribution the potential V can be calculated and the Dirac equation for muon and electron be solved numerically. The relevant formulae are given in Ref. [21], which we follow closely. The Dirac equation provides us with the muon binding energy E b and the wavefunctions of bound muon and free electron in the Coulomb background, which ultimately de-arXiv:2110.14667v2 [hep-ph] 18 Mar 2022 termine the decay rate. Ref. [21] also provided an expansion of the DIO spectrum near the endpoint of the form is the leading-order free-muon decay rate and B is a coefficient of mass dimension −6 that depends on several overlap integrals over products of wavefunctions.
Here, we improve the endpoint expansion of Ref. [21] by implementing the leading radiative corrections: with δ = α(2 log[2m µ /m e ] − 2)/π 0.023 [23] coming from soft photon radiation and a shifted endpoint energy [25] due to the vacuum polarization effects. Here, E end is given by Eq. (1). Notice that E end coincides with the electron energy in µ − → e − conversion. In the above formulas, α = 7.297 352 5693(11) × 10 −3 [32] is the fine structure constant. For a target that consists of several isotopes, the above spectrum dΓ/dE e is to be summed over the isotopes, weighted by their abundance. This is necessary because both the coefficient B and the binding energy E b depend to a small degree on the number of neutrons in the nucleus. In the following we will determine how large this dependence is.

III. NUCLEAR CHARGE DISTRIBUTIONS
The electric-charge distribution of nuclei will be approximated as spherically symmetric in the following and denoted as ρ(r), normalized according to Eq. (3). Experimental information about ρ can be obtained via spectroscopy in (muonic) atoms and through scattering. For example, electron-nucleus scattering cross sections with momentum transfer q in the Born approximation are proportional to the squared form factor |F (q)| 2 , which is related to ρ(r) by Fourier transformation: We have normalized F (q) so that it gives 1 for a point charge. Since scattering experiments can only be performed over a finite range of q, it is not possible to obtain ρ model-independently; rather, once a parametrization for ρ(r) is chosen, F (q) can be calculated and fitted to data. Popular parametrizations for spherically symmetric charge distributions with varying degrees of complexity are listed below: 1. Three-parameter Fermi model (3pF) [33,34]: The two-parameter Fermi model (2pF) can be obtained as the special case w = 0; the one-parameter Fermi function (1pF) is defined here through w = 0 and z = 0.52 fm (which corresponds to a constant surface thickness of 2.3 fm).
5. Sum of Gaussians (SOG) [37]: with n k Q k = 1. The normalization ρ 0 can be calculated analytically in all cases and the potential V in Eq. (2) as well, except for the 3pG model, and are given in appendix A.
The simplest parametrization, 1pF, takes as input only the nuclear charge radius, listed for all elements and isotopes of interest in Ref. [38]. For many isotopes this is the only available parametrization, making it extremely valuable despite its lack of substructure. For the multiparameter parametrizations we use the available data tables from Refs. [39][40][41][42][43], generally choosing the newest possible data set for a given isotope.
For several parametrizations the errors on the fit parameters are either not given in the literature (e.g. for FB) or not representative of the uncertainties in the charge distribution (e.g. the errors on the charge radius relevant for 1pF are tiny). To obtain an estimate for the error on our DIO spectrum from the uncertainty of the nuclear charge distribution we will compare the results of several parametrizations for a given isotope.
We restrict our survey to the 237 stable isotopes with natural abundance above 1%. Only half of these isotopes have a charge parametrization with two or more parameters, for the other half we have to rely on the simplistic 1pF model.

IV. RESULTS
For a given charge distribution we solve the Dirac equations numerically using a fourth-order Runge-Kutta method to extract the muon binding energy E b and evaluate the coefficient B, which can then be put into Eq. (4). Our results agree with Ref. [21] for the six elements presented there, using their charge parametrizations. Our results for the binding energy also match the point-charge [44,45]. Table I provides our results and will be discussed in the following.

A. Binding energy
Our results for the muon binding energy E b as a function of Z are shown in Fig. 1 (top) for all elements, isotopes, and available parametrizations. E b (Z) is very well behaved and agrees with the point-nucleus approxima- For Z > 50, the Z dependence becomes linear and is approximately given by Since the differences between the different parametrizations cannot be seen in Fig. 1 (top), we show in Fig. 1 (middle) the difference between the 1pF parametrization result with all the other available parametrizations. The relative difference is below 1% for all Z and markedly smaller for small Z. We take this as a sign that the endpoint energy does not depend strongly on the details of the charge distribution and is hence given by the 1pF results for all isotopes with at least 1% accuracy, and even far better for smaller Z (say Z < 20). Using the 1pF results we can illustrate the isotope dependence of E b by taking the difference between E b (Z, N ) and E b (Z, N 0 ), N 0 being the neutron number of the lightest stable isotope. This is shown in Fig. 1 (bottom). The relative changes between isotopes are below 1%, once again smaller for smaller Z. Overall they are of similar size as the parametrization uncertainty.
Ultimately, the quantity of interest for µ − → e − conversion experiments is not the binding energy E b but the related endpoint energy E end , given by Eq. (5). Since E b is itself only a small correction to the leading-order estimate E end ∼ m µ , small uncertainties in E b are even less relevant in E end . In fact, the relative parametrization uncertainty as well as the isotope dependence of E end are below one per mille for all Z. Therefore, the endpoint energies (shown in Fig. 2) are well approximated by the 1pF results, available for all isotopes of interest.
Notice that we have approximated the nuclear-recoil energy in Eq. (1) at leading order in 1/m N , which is an excellent approximation for heavy nuclei; for light nuclei, e.g. 6 3 Li, the subleading recoil terms can be of order m 3 µ /m 2 N ∼ 0.04 MeV and pose the largest remaining uncertainty of the endpoint energy.  Endpoint energy E end from Eq. (5) vs Z for all isotopes in all available charge-distribution parametrizations. The relative uncertainty as well as isotope dependence is below one per mille.

B. B coefficient
Having determined the electron endpoint for all isotopes, we can move on to the more difficult task of determining the B coefficient in Eq. (4), which provides the normalization of the spectrum and hence the number of electrons within the signal window of a given µ − → e − conversion experiment. Our results are shown in Fig. 3, which illustrates that the dependence of B on the charge distribution, and by extension the isotope, is much larger than for the binding or endpoint energies, at least for large Z.
To be more quantitative, we show in Fig. 4 (top) the differences in B for charge parametrizations relative to 1pF. The differences for small Z are typically below 5% and grow to 10% for larger Z. This serves as an estimate for the uncertainty on B due to ρ. In Fig. 4 (bottom) we also see that the isotope differences are often of the same order or even larger than these uncertainties and should not be neglected.

C. Connection to form factor
As mentioned above, the finite size of the nucleus in electron-nucleus scattering can, to first order, be incorporated by multiplying the point-nucleus cross section by the form factor squared at the relevant momentum transfer q: |F (q)| 2 . The same argument can be made for DIO, where we expect B ∝ |F (m µ )| 2 at leading order, with q = m µ the relevant momentum-transfer scale [23]. Since the nuclear charge distribution complicates the numerical solution of the Dirac equation it would be convenient if this effect could indeed be factored out. To test the scaling B ∝ |F (m µ )| 2 , we calculated the ratio of B/|F (m µ )| 2 for the 1pF charge distribution and compared it to the available multi-parameter charge distri- butions (FB, SOG, 3pG, 2pF/3pF). The result is shown in Fig. 5 and illustrates that the scaling B ∝ |F (m µ )| 2 holds to better than 2% for all Z, at least in the region of interest for us. This makes it possible to simply rescale our B results obtained with ρ 1pF by the ratio of form factors to obtain fairly accurate estimates of B with new data, without the need to numerically solve Dirac equations again. In addition, the scaling B ∝ |F (m µ )| 2 opens up a different way to estimate the uncertainty on B due to the charge distribution. Rather than trying to estimate the B uncertainty from the uncertainty of ρ(r) -which is itself difficult to estimate as ρ is not a measured observable -we can look directly at electron-nucleus scattering data to obtain |F (q = m µ )| 2 including uncertainties. Despite the ∼ 2% uncertainty in the scaling B ∝ |F (m µ )| 2 (Fig. 5) this can be beneficial compared to the usual chain of measuring F (q), fitting ρ(r), calculating V (r), and solving Dirac equation to get B, which is quite tedious to propagate errors. Notice that electron-nucleus scattering data around q ∼ m µ 0.53 fm −1 is not available for all isotopes, nor is the actual form factor plus error bars given explicitly in many cases, rendering this  approach unsuitable in general. As a concrete example for the above procedure we consider 6 3 Li. The B coefficients using 1pF and FB [43] are B = 1.20 × 10 −19 MeV −6 and B = 1.29 × 10 −19 MeV −6 , respectively. Neither come with any reliable error bars so the best we can do is interpret the difference between the two parametrizations as an estimate of the uncertainty, which is roughly 7%. 1 Using instead the form factor data from Ref. [43] underlying the FB parametrization, we have the interpolated value |F (m µ )| 2 = 0.57 with 2.5% (systematic) uncertainty; since the experimental form factor value is larger than the FB fit at q ∼ m µ [43], B should be larger as well, and can be estimated as with an error of 3%, generously adding the experimental uncertainty and the B ∝ |F (m µ )| 2 scaling uncertainty in quadrature. This gives us a more reliable result for B and its error, which in this case happens to be even smaller than our naive estimate.

V. CONCLUSIONS
The search for lepton flavor violation is one of our most sensitive probes of physics beyond the Standard Model. Experiments searching for µ − → e − conversion such as COMET, DeeMe, and Mu2e, promise to improve existing limits by several orders of magnitude. An observation would make it possible to track the nucleus dependence of the µ − → e − conversion rate by looking at different target materials, which would then help to discriminate the possible underlying new-physics models and operators. Muon decay in orbit is an irreducible background in all µ − → e − conversion searches and thus needs to be calculated to sufficient precision. In this article, we have performed a broad study of the DIO spectrum near the endpoint for all stable isotopes that could conceivably be used experimentally, extending previous studies. Our results are sufficient for background studies involving broad range of different or enriched target materials. Once the experiments choose alternative target materials, the next step in our investigation will be a detailed study of the DIO spectrum over the complete energy range.
Appendix A: Analytic expressions for the potentials For completeness we give the normalization constants ρ 0 (defined via Eq. (3)) and the electric potentials V (defined via Eq. (2)) for all charge parametrizations of interest.

Three-parameter Fermi model (3pF)
The normalization is analytically given by where Li s is the sth-order polylogarithm function. The analytic formula of the nuclear potential is where χ 0 = exp (−c/z) and χ = χ 0 exp (r/z).
Here the integral I(r) is defined by

Modified-harmonic oscillator model (MHO)
The normalization is and the potential takes the form where erf is the error function.

Fourier-Bessel expansion (FB)
The parameters a k (1 ≤ k ≤ n) are normalized by The potential is

Sum of Gaussians (SOG)
The potential is