Probing $\mu e \gamma \gamma$ contact interactions with $\mu \to e$ conversion

Contact interactions of a muon, an electron and two photons can contribute to the decay $\mu \to e \gamma \gamma$, but also to the conversion of a muon into an electron in the electric field of a nucleus. We calculate the $\mu \to e$ conversion rate, and show that for the coefficients of operators involving the combination $FF \propto |\vec{E}|^2$ (as opposed to $F\tilde{F} \propto \vec{E} \cdot \vec{B}$), the current bound on $\mu \to e$ conversion is more sensitive than the bound on $\mu \to e \gamma \gamma$. Constraining these $\bar{e}\mu FF$ operators gives the best sensitivity to some dimension six operators.


Introduction
The observed neutrino masses imply the existence of contact interactions where charged leptons change flavour. This is referred to as (Charged) Lepton Flavour Violation (CLFV) and is reviewed for muon decays in, eg [1]. Current constraints on several µ ↔ e flavour changing processes are restrictive, and experiments under construction [2,3,4] aim to reach BR ∼ 10 −16 . Some bounds and future sensitivities are given in table 1.
If CLFV is discovered, experimental bounds on, or observations of, a multitude of independent processes would assist in discriminating among models. This motivates our interest in the less commonly considered contact interactions involving a muon, an electron and two photons. Such interactions could mediate various processes, such as µ → eγγ and µ → e conversion in the electric field of a nucleus. The rate for µ → eγγ was calculated by Bowman, Cheng, Li and Matis (BCLM) [5], whose results are reviewed in section 2, and an experimental search with the Crystal Box detector obtained BR(µ → eγγ) ≤ 7.2 × 10 −11 [6].
We will parametrize CLFV interactions via contact interactions involving Standard Model (SM) particles. This would be appropriate if the new particles involved in CLFV are heavy, but may not be generic for µ → eγγ. This decay could be mediated by µ → ea [7] followed by a → γγ, where a is a light (pseudo) scalar such as an Axion-Like Particle [8]. Recently, the MEG experiment searched for the collinear photons from this process [9]. They found that the branching ratio of µ + → e + a, a → γγ is smaller than O 10 −11 when the mediator a has a mass of 20-45MeV and a lifetime below 40 ps.
In this manuscript, we calculate the µ → e conversion rate induced by contact interactions of µ, e and two photons, and discuss the contribution of such operators to constraining and distinguishing among New Physics models. Section 2 introduces the basis of operators (previously given by BCLM [5]), and gives their contribution to µ → eγγ. The operators are of dimension seven and eight, so we recall that some of the dimension seven operators can be generated from dimension six scalar operators. Our calculation of µ → e conversion mediated by the eµF F operator is presented in Section 3, where we find a significant dip in the rate for targets of intermediate Z. The final discussion section integrates our results in the usual expression for the spin independent Branching Ratio of µ → e conversion, and discusses the resulting sensitivities to various operator coefficients. An Appendix A gives some details of our claim that the eµF F operators give the best sensitivity to a some dimension six operators with a scalar contraction of heavy fermions. This Appendix employs both the QED×QCD-invariant EFT suitable below the weak scale, and the SMEFT (where we also mention an accidental cancellation in the contribution of the LFV Higgs coupling operator O EH to µ → e conversion). process current sensitivity future µ → eγ < 4.2 × 10 −13 (MEG [10]) ∼ 10 −14 (MEG II [11]) µ → eγγ < 7.2 × 10 −11 (Crystal Box [6]) µ → eēe < 1.0 × 10 −12 (SINDRUM [12]) ∼ 10 −16 (Mu3e [4]) µA→ eA < 7 × 10 −13 (SINDRUM II [13]) ∼ 10 −16 (COMET [2], Mu2e [3]) ∼ 10 −18 (PRISM/PRIME [14]) Table 1: Current bounds on the branching ratios for various CLFV processes, and the expected reach of upcoming experiments.

Notation and Review
A set of QED-invariant operators that could mediate the decay µ → eγγ was given by Bowman, Cheng, Li and Matis (BCLM) [5]: where two changes have been made to their notation: the New Physics scale in the denominator is taken to be the Higgs vacuum expectation value v ≃ m t , with 2 √ 2G F = 1/v 2 (BCLM took m µ ), and we use chiral fermions, because this facilitates matching onto the full SM at the weak scale, and because the out-going electrons are relativistic so ≈ chiral.
This basis of operators is constructed to include all possible Lorentz contractions that give the desired external particles (there is no tensor, because there is no two-index antisymmetric combination of F F or FF to contract with eσµ), so corresponds to a general parametrization of the interaction at lowest order in a momentum expansion. The resulting operators are of dimension seven and eight.
Curiously, all the operators of eqn (1) induce a matrix-element-squared for µ(P µ ) → e(p e ) + γ(k) + γ(q) that is proportional to [5] |M| 2 ∝ P µ · p e (k · q) 2 , giving a Branching Ratio where Equation (2) is 1/2 the result given in BCLM, because they label their γ 1 as more energetic than γ 2 , whereas we divide by 2 for identical final state particles. The experimental bound from Crystal Box ( [6]) given in table 1 therefore corresponds to In many heavy New Physics models, LFV arises at dimension six. So for a sufficiently high New Physics scale Λ N P , it is reasonable to neglect the dimension ≥ 7 operators that could be generated at Λ N P , because their contributions to observables will be suppressed by additional factors of E expt /Λ N P . However, some of the operators of eqn (1) can arise at O(1/Λ 2 N P ) in a CLFV New Physics model, with a SM mass scale providing the additional dimensions in the denominator. For example, if the scalar operators are present in the Lagrangian as δL = 1 Λ 2 N P (C ψψ S,XX O ψψ S,XX + C ψψ S,XY O ψψ S,XY ), then at a heavy fermion mass scale m ψ , they match onto the two-photon operator O F F,X via the diagram of figure 1, with coefficient (6) where N c = 3 for heavy quarks, and is one otherwise. This result is related to the conformal anomaly [15], and is the QED version of the matching of scalar heavy quark operators onto gluons, performed by Shifman, Vainshtein and Zakharov [16], and included in µ → e conversion by [17]. The dimension eight two-photon operators O V F F,X =ēγ σ P X µF αβ ∂ β F ασ appear more difficult to obtain at O(1/Λ 2 N P ). Furry's Theorem says that dimension six vector operators, such as (ēγ α P X µ)(ψγ α ψ), do not match onto O V F F,X via the diagram of figure 1, because an odd number of vector current insertions appear on the fermion loop. Writing the loop of figure 1 with external legs amputated and an axial heavy fermion current (ēγ α P X µ)(ψγ α γ 5 ψ) in the grey blob, gives a vacuum matrix element that is even under Charge conjugation, but odd under CP. Analogously to Furry's theorem, it should vanish in a CP invariant theory, so we do not calculate this diagram in the approximation of CP invariance.
In the next section, we will calculate the contribution of the scalar O F F,X operators to coherent µ → e conversion. Since µ → e conversion can be induced by dimension six LFV operators, we do not consider the contribution of the dimension eight O V F F,X operators, on the assumption that they are relatively suppressed by O(E 2 expt /Λ 2 N P ).

Calculating the µ → e conversion rate
We consider the coherent µ → e conversion described by the first two terms of eqn (1). Assuming a particle picture in which the muon and the scattering electron are independently described by its own wave function in a Coulomb potential, the transition matrix is where ψ 1s µ and ψ e are the wave functions of a 1s bound muon and a scattering electron, respectively. Here, we omit spin indices for simplicity. |N denotes the ground state of a nucleus. For an ordinary nucleus, we can safely assume that the electric field E (r) is spherically symmetric and the magnetic field is negligible. Then we relate the hadronic matrix element with a classical field strength as follows, With the amplitude M, the conversion probability is given by where the summation includes spin averaging of the initial state, and E conv e is the energy of the signal electron, given by E conv Here B µ is the binding energy of initial muon in the muonic atom. The lepton wave functions ψ ℓ (ℓ = e, µ) obey the Dirac equation in a nuclear Coulomb potential; our formulation below follows [18,19].
For a spherically symmetric potential, one can represent the wave function of the bound muon as where χ is a two-component spherical spinor 1 . The differential equations for the radial wave functions G(r) and F (r) are obtained from the Dirac equation as follows, The nuclear Coulomb potential V C is calculated with a nuclear charge density ρ(r) as, For the nuclear density, we adopted two different models, the two-parameter-Fermi distribution (2pF) and threeparameter Gaussian distribution (3pG), given by The normalization has been used such that Ze = 4π ∞ 0 ρ(r)r 2 dr with the normalization factors ρ 0 for each type of distribution. The parameters, ω, c and z, are listed in Refs. [20,21]. For simplicity of formulation, we express the wave function of the outgoing electron of momentum p e using the partial wave expansion: where j κ and l κ are the total and orbital angular momentum, respectively. We introduced an integer quantum number κ that runs from −∞ → ∞ skipping 0, and determines j and l as j κ = |κ| − 1/2 and l κ = j κ + κ/2|κ|. Due to angular momentum conservation, only the waves with κ = ∓1 contributes to µ → e conversion. δ κ is a phase shift of the partial wave with κ, and the incoming boundary condition is taken from [19]. (l κ , m, 1/2, s e |j κ , ν) is the Clebsch-Gordan coefficient, and Y m lκ (p e ) is a spherical harmonic. The radial Dirac equations for each partial wave are The normalization of the wave functions is the same as Ref. [18]. Then the conversion probability is where the overlap integrals F − A and F + A for a target nucleus A are Neglecting the electron mass, we have Eventually we obtain the conversion probability, where we redefine Fig. 2. For instance, F A for aluminum (Z = 13) and gold (Z = 79) are 3.8 × 10 −4 and 6.1 × 10 −3 , respectively. F A for other targets are listed in Appendix B.
A few nuclei are modeled by both the 2pF and 3pG distributions, in which case we give the results with the latest distribution. Apart from the dip around Z = 38 (which will be discussed later), different distribution models lead to the same results within O(1)% accuracy. The squared electric field of heavy nuclei E(r) 2 ∝ Z 2 renders the overlap of lepton wave functions and the electric field large.
In Fig. 2, one sees a dip in the overlap integral in the range 30 Z 50. In order to interpret this cancellation, the integrand of F A is plotted as a function of radius in Fig. 3 for Z = 13, 38, and 79. The oscillations arise from the electron wave function g −1 , whose first node is at r ≃ π/m µ ≃ 5.8 fm. Since the electric field is maximized around the nuclear surface, there is a significant cancellation between the interior and exterior contributions to the integral when the first node of g −1 is close to the nuclear radius. The overlap integral is minimised at 35 Z 40, where the nuclear radius is about 5.5 fm. We define the branching ratio as which is shown in Fig. 4. Here Γ cap stands for the muon capture rate in a muonic atom [22]. For aluminum and gold, we have |C F F,L | 2 + |C F F,R | 2 = 6.7 × 10 −9 for 27 Al 9.1 × 10 −8 for 197 Au . Figure 4 shows the branching ratios for targets of atomic number Z normalized by that for aluminum target. The branching ratio via the dipole CLFV operator, O D,X =ēσ αβ P X µF αβ , is also shown to highlight the difference of Z dependence. (The µ → e conversion by the dipole operator has been studied in detail [18].) The conversion rate via theēµF F operator receives an additional overall Z factor compared with the dipole case, due to the extra F µν factor, and more rapidly increases with Z than the muon capture rate Γ cap ∝ Z 4 ef f [23]. The resultant branching ratio continues to increase at large Z, which differs from the high-Z falloff of the better known dipole, scalar or vector operators [18].
A noticeable dip occurs at 30 Z 50 in Fig. 4 due to the accidental cancellation in the overlap integral F A for the operatorēP L(R) µF F . We stress that this dependence could be used to discriminate this operator from other CLFV operators. However, the precise prediction of the dip are difficult, since the overlap integral at 30 Z 50 is very sensitive to the nuclear model, and some parameters in the muonic atom (such as the muon binding energy at O(1)% level). In order to reliably predict BRs for these targets, it would be necessary to model the nuclear distributions with considerable accuracy. Figure 4: The µ → e conversion branching ratios normalized by that for aluminum target. Solid and dashed lines represent the cases where the O F F,X =ēP X µF αβ F αβ and O D,X =ēσ αβ P X µF αβ dominate, respectively.

Summary
In this manuscript, we calculated the contribution ofēµγγ contact interactions to µ → e conversion on nuclei. We considered the first two operators of eqn (1), which are CP-even, of dimension seven, and involve F µν F µν ⊃ | E| 2 . Other possibilities are discussed in section 2.
Our calculation is outlined in section 3, and the resulting conversion rate is defined in eqn (21), and given in Appendix B. It relies on the overlap integrals in the nucleus, of the electron and muon wave functions with the electric-field-squared, which are given in eqns (19,20). Our results can be included in the well-known Branching Ratio for Spin Independent µ → e conversion [18] as where D A and S (N ) A are the overlap integrals inside the nucleus A, with respectively the electric field or the appropriate nucleon (N ∈ {n, p}) distribution, which can be found in [18]. Γ cap is the muon capture rate on nucleus A [22], C D,X is the dipole coefficient, {C qq S,X } are the coefficients of 2 √ 2G F (ēP X µ)(qq), and the "+..." represents the contributions of vector operators involving a light quark bilinear. This expression uses the quark densities in the nucleon of References 2 [24,25,26], the gluon density [16,17] and we included the overlap integral F A for the operators O F F,X , calculated in this work and given in figure 2 and Appendix B.
The SINDRUM II experiment searched for µ → e conversion on gold, and obtained the upper bound BR(µAu → eAu) ≤ 7 × 10 −13 [13]. This corresponds to a bound on the "gauge boson operator" coefficients: which gives, in the absence of C D,X and C GG,X , This is a better sensitivity than that given in eqn (4) Table 2: µēγγ operator coefficients bounded by µ → eγγ [6], and the sensitivity of µAu → eAu [13] obtained in this manuscript. The operators are given in eqn (1) The upcoming COMET and Mu2e experiments plan to start with an Aluminium target. Combining eqn (24) and F Al = 3.8 × 10 −4 , we obtain a future sensitivity of The sensitivity to C F F,X would be improved by an order of magnitude with an expected branching ratio of ∼ 10 −16 on the light target Aluminum. The contribution of the F F operator to µA → eA has an interesting dependence on Z (the electric charge of the target nucleus), illustrated in Fig. 4 and discussed at the end of section 3. There is a significant cancellation in the branching ratio at intermediate Z, beyond which the BR increases monotonically, (unlike more familiar operators, where there is a falloff at large Z).
Finally, we comment on the interest of the (ēP L,R µ)F F operators in identifying heavy New Physics in the lepton sector. These operators are of dimension seven in the QED×QCD-invariant EFT below m W , and dimension eight above. However, section 2 mentioned that they can be induced in matching out heavy fermion scalar operators of dimension six, as illustrated in figure 1, and given in eqn (5). Some such scalar operators, in particular (ēP L µ)(τ P R τ ) and (ēP R µ)(τ P L τ ) do not appear to contribute at one-loop to any other well-constrained observables. In the Appendix A, it is shown that µA → eA gives the best sensitivity to the coefficients of such operators at the weak scale: The "bound" is also quoted in the SMEFT operator basis {in curly brackets}. This is the best sensitivity to these coefficients that the authors are aware of.
In this appendix, we assume that New Physics generates dimension six CLFV operators at some scale Λ N P > m W , which is high enough that dimension ≥ 7 operators can be neglected. Then we explore whether including O F F,X in the calculation of µ → e conversion improves the sensitivity to any dimension six operators. Incidentally, we observe that in the SMEFT, the contributions of the LFV-Higgs coupling operator O HE to µ → e conversion, via the dipole and GG operators, are of comparable magnitude and opposite sign, weakening the sensitivity of µ → e conversion to O HE by a factor ∼ 1/2, so this is also discussed.
We first consider the QED×QCD invariant EFT below m W . As illustrated in figure 1, the scalar operators O ψψ S,XX , O ψψ S,XY (see eqn 5), match at m ψ onto O F F,X and O GG,X : where Q ∈ {c, b, t}. Here we take Λ N P = v. We focus on ψ ∈ {τ, c, b, t} a heavy fermion, because the operators with ψ ∈ {e, u, d, s} contribute at tree level to µ → e conversion or µ → eēe, and for ψ = µ, the operator contributes at one loop to µ → eγ. It is interesting to pursue the loop effects of these heavy-fermion scalars, because the two heavy fermions make the operators difficult to probe directly in experiment.
The O ψψ S,XX scalar operators (with the same chiral projector in both bilinears) contribute to the dipole operator via Barr-Zee diagrams. The log 2 -enhanced part is given by the one-loop RGEs of QED [28,29] as where C D,X is the dipole coefficient at the experimental scale, and the coefficients on the right are evaluated at m W . For the heavy quark ( Q ∈ {c, b}) scalar operators O QQ S,XX , O QQ S,XY , the contribution to µ → e conversion via the O GG,X operator is clearly larger than via O F F,X (see eqns 31,24). And it turns out, that despite the restrictive current experimental bound on µ → eγ (see table 1), the SINDRUMII search for µ → e conversion has better sensitivity to the scalar heavy quark operators via eqn (31). For scalar operators with a τ bilinear, the dipole has better sensitivity to the O τ τ S,XX operators, but only O F F,X allows to probe O τ τ S,XY . So indeed, including O F F,X allows to probe additional operators of dimension six.
The discussion so far has been in the context of QED×QCD invariant operators below the weak scale. However, since we assume Λ N P is large, it is relevant to translate to the SMEFT, where SU(2) invariance restricts the operator basis to three scalar four-fermion operators at dimension six: the XX scalar for u-type quarks, and the XY scalars for d-type quarks and charged leptons. There is also a flavour-changing Higgs coupling. Including also the dipoles, these operators appear in the SMEFT Lagrangian as δL SMEF T = 1 v 2 C µe EH H † Hl µ He + C eµ EW y β (ℓ e τ a Hσ µν e µ )W a µν + C eµ EB y β (ℓ e Hσ µν e µ )B µν (33) +C eτ τ µ LE (ℓ e γ µ ℓ τ )(e τ γ µ e µ ) + C τ µeτ LE (ℓ τ γ µ ℓ µ )(e e γ µ e τ ) +C eµnn LEQU (ℓ A e e µ )ε AB (q B n u n ) + C eµnn LEDQ (ℓ e e µ )(d n q n ) + h.c. , where the capitalized SU(2) indices are explicit when not contracted in the parentheses, ℓ and q are doublets, u, d, e are singlets, flavour indices are superscripts, n ∈ {c, t, b}, and the operator labels are according to [30]. The O eµ EW and O eµ EB will combine to the dipole, the O LE operators Fiertz to XY scalar operators with a τ bilinear, and in the quark sector, O LEQU is a Y Y -scalar operator (same chiral projector twice), whereas O LEDQ is XY .
Loop effects between Λ N P and the weak scale can be partially included via the Renormalization Group Equations of the SMEFT. Gauge boson loops can renormalize the coefficients, and mix the C eµnn LEQU coefficients into the u−type tensor operator, and then to the dipole (as occurs below m W for Y Y scalars). Higgs exchange can mix these scalars into vector four-fermion operators (to which there could be better experimental sensitivity), but for O LEDQ and O LE , this is negligible because suppressed by ∼ y µ y ψ /(16π 2 ) (ψ ∈ {τ, b}). We therefore suppose that the coefficients in eqn (33) are given at the weak scale m W , since the one-loop RGEs above m W do not appear to significantly mix the XY -scalars into more experimentally accessible operators.
The coefficients from eqn (33) can be matched at m W onto those of QED×QCD-invariant scalar four-fermion operators, relevant at low energy. All the scalar operators below m W are generated at tree level, just that some arise due to Higgs exchange with a flavour-changing coupling from the O HE operator, leading to correlations in the coefficients. One obtains [28] where s W = sin θ W , all the masses and couplings are running, and are evaluated at the weak scale. The two-loop Barr-Zee diagrams involving top and W loops were included in the matching to the dipole, and the top loop matching the scalar operators onto F F and GG was included for these operators. These coefficients then run down to the experimental scale with the Renormalization Group Equations of QED×QCD (see, eg, [29]). QCD effects are numerically significant, although they only renormalize the coefficients 3 . The scalar quark operators run like quark masses, and the operator O GG,X runs like the gluon kinetic term, which is accounted for by the wave function renormalization of the gluons. So the running parameters in the coefficient are evaluated at the matching scale.
Including the interesting QED effects (such as the mixing of LL and RR scalars to tensors and then into the dipole), and the matching at the heavy fermion mass scales given in eqn (31), gives the coefficients that will contribute to µ → e experiments as : where on the right appear SMEFT coefficients evaluated at m W , and the coefficients on the left can be input into the rate for µ → e conversion. Here C eµ eγ = c W C eµ EB − s W C eµ EW .
Combining with eqn (26), one can see that the largest contribution to µ → e conversion of the operators O eµcc LEQU , O eµtt LEQU , and O eµbb LEDQ is via O GG,X 4 , and, to the accuracy we calculate, the only contribution of the O LE to µ → e conversion arises via O F F,X . The contribution of the LFV Higgs interactions O eµ EH via O GG,X is of opposite sign and 1 3 the magnitude of the dipole contribution. The contribution via the light quark (u, d, s) scalar operators is slightly smaller than the GG contributions and of same sign, which worsens the sensitivity of µ → e conversion to O EH 5 . Including both effects, µAu → e + Au cannot see C eµ EH ≤ 4.9 × 10 −5 (50) whereas including only the dipole would give a sensitivity of < ∼ 1.6 × 10 −5 .