Studying ∆ L = 2 Lepton Flavor Violation with Muons

Flavor violating processes in the lepton sector have highly suppressed branching ratios in the standard model. Thus, observation of lepton flavor violation (LFV) constitutes a clear indication of physics beyond the standard model (BSM). We review new physics searches in the processes that violate the conservation of lepton (muon) flavor by two units with muonia and muonium–antimuonium oscillations.


Introduction
Flavor-changing neutral current (FCNC) interactions serve as a powerful probe of physics beyond the standard model (BSM). Since no local operators generate FCNCs in the standard model (SM) at tree level, new physics (NP) degrees of freedom can effectively compete with the SM particles running in the loop graphs, making their discovery possible. This is, of course, only true provided the BSM models include flavor-violating interactions. If the new physics particles are heavier than the muon mass, their effect on muon transitions can be parameterized in terms of local operators of increasing dimension. In fact, any new physics scenario which involves lepton flavor violating interactions can be matched to an effective Lagrangian of the Standard Model Effective Field Theory (SM EFT), where the Wilson coefficients (WC) c i of L SMEFT in Equation (1) are determined by the UV physics that becomes active at some scale Λ [1][2][3]. The effective operators Q i 's defined in Equation (1) reflect degrees of freedom relevant at the scale at which a given process takes place. For heavy new physics scenarios, those operators should be written in terms of muon, electron, neutrino, or photon (gluon) fields. Then to perform low energy computations, it is often convenient to match the SM EFT Lagrangian to the low energy Lagrangian L eff only containing these degrees of freedom. In addition, it is convenient to classify the operators of L eff by their lepton quantum numbers L (with = e, µ, τ), We shall employ the operators that only contain the fermion fields, so the operators in Equation (2) start at dimension six.
The first term in Equation (2) contains both the standard model and the new physics contributions. The SM piece is given by the usual Fermi Lagrangian, Since Λ m W , which is required for the consistency of SM EFT, we can neglect the NP piece of the ∆L µ = 0 Lagrangian for studies of FCNC transitions. This implies that L ∆L µ =0 eff is suppressed by powers of G F ∼ M −2 W , the Fermi constant. We should emphasize that only the operators that are local at the scale of the muonium mass are retained in Equation (2).
The second term in Equation (2) contains ∆L µ = 1 operators. The most general low energy Lagrangian representing muon flavor change by one unit, L ∆L µ =1 eff , is conventionnaly written as [4][5][6] L ∆L µ =1 eff = − 1 where µ and e are the fermion fields with (µ, e) L,R = P L,R (µ, e), where P R,L = 1 2 1 ± γ 5 are the projection operators, and f represents other fermions that are not integrated out from the low energy effective theory. We introduced indices for the Wilson coefficients C f IK to indicate the Lorentz structure of effective operators: vector, axial-vector, scalar, pseudo-scalar, and tensor. The Wilson coefficients are in general different for different fermions f . The operators additionally suppressed by factors of m i m j G F usually arise from SM EFT operators of dimension higher than six [4,5].
The last term in the low energy EFT described in Equation (2), L ∆L µ =2 eff , represents a collection of effective operators changing the lepton quantum number by two units. The most general effective Lagrangian for such transition that is applicable to the physics we will discuss in this review is given by with the operators written entirely in terms of the muon and electron degrees of freedom, Other Lorentz structures of the operators could be related to the ones in Equation (6) via Fierz relations. A ∆L = 2 interaction described by this Lagrangian can lead to muon decays such as µ + → 3e. It can also change a µ + e − bound state into a µ − e + bound state, leading to muonium-antimuonium oscillations. We will discuss those below.
In addition, the operators built out of muon and electron fields, we can construct other ∆L µ = 2 local operators that contain muon, electron, and neutrino fields [7], where we only included operators that contain left-handed neutrinos [3,8]. The operators in Equation (7) could lead to both muon decays µ → eν µνe (note the "reversed" neutrino flavors), and muonium-antimuonium oscillations.
The goal of experimental studies of ∆L = 2 transitions would be to find the observables that are sensitive to various combinations of Wilson coefficients of the effective Lagrangian in Equation (2). Different processes are in general sensitive to different combinations of the Wilson coefficients, which raises hope that a sufficient number of measurements would allow placing constraints on individual WCs without additional assumptions such as single operator dominance [5,6]. In that respect, studies of both unbound muon and muonium decays are needed.

Muonium: The Simplest Bound State
Even though the hydrogen atom is a quintessential quantum-mechanical bound state and is usually presented as the easiest QED problem, it is not the simplest. Precision studies of hydrogen reveal effects related to the structure of the proton: its finite size and composite nature. In that respect, the simplest hydrogen-like system is built entirely out of leptons. Such a system that is especially clean for studies of BSM effects in the lepton sector is muonium M µ . The muonium is a QED bound state of a positively-charged muon and a negatively-charged electron, |M µ ≡ |µ + e − .
Since muon is unstable, M µ is unstable. The main decay channel for the muonium is determined by the weak decay of the muon, M µ → e + e −ν µ ν e , so the average lifetime of a muonium state τ M µ is expected to be the same as that of the muon, with τ µ = (2.1969811 ± 0.0000022) × 10 −6 s [9], apart from the tiny effect due to time dilation [10]. Note that Equation (8) represents the leading-order result. The results including subleading corrections are available [10]. Like a hydrogen atom, muonium could be formed in two spin configurations. A spin-one triplet state M V µ is called ortho-muonium, while a spin-zero singlet state M P µ called paramuonium. In what follows, we will drop the superscript and employ the notation M µ if the spin of the muonium state is not important for the discussion.

Muonium-Antimuonium Oscillations
Since ∆L = 2 interaction can change the muonium state into the anti-muonium one, the possibility to study muonium-anti-muonium oscillations arises. Theoretical analyses of conversion probability for muonium into antimuonium have been performed, both in particular new physics models [11][12][13][14][15][16], and using the framework of effective theory [7], where all possible BSM models are encoded in a few Wilson coefficients of effective operators. Observation of muonium converting into anti-muonium provides clean probes of new physics in the leptonic sector [17,18].

Phenomenology of Muonium Oscillations
In order to determine experimental observables related to M µ − M µ oscillations, we recall that the treatment of the two-level system that represents muonium and antimuonium is similar to that of meson-antimeson oscillations [1,19,20]. There are, however, several important differences. First, both ortho-and para-muonium can oscillate. Second, the SM oscillation probability is tiny, as it is related to a function of neutrino masses, so any experimental indication of oscillation would represent a sign of new physics.
In the presence of the interactions coupling M µ and M µ , the time development of a muonium and anti-muonium states would be coupled, so it would be appropriate to consider their combined evolution, The time evolution of |ψ(t) evolution is governed by a Schrödinger-like equation, where is a 2 × 2 Hamiltonian (mass matrix) with non-zero off-diagonal terms originating from the ∆L = 2 interactions. CPT-invariance dictates that the masses and widths of the muonium and anti-muonium are the same, so m 11 = m 22 , Γ 11 = Γ 22 . In what follows, we assume CP-invariance of the ∆L µ = 2 interaction 1 . Then, The off-diagonal matrix elements in Equation (11) can be related to the matrix elements of the effective operators introduced in Section 1, as discussed in [1,19], To find the propagating states, the mass matrix needs to be diagonalized. The basis in which the mass matrix is diagonal is represented by the mass eigenstates |M µ 1,2 , which are related to the flavor eigenstates M µ and M µ as where we employed a convention where CP|M µ ± = ∓|M µ ± . The mass and the width differences of the mass eigenstates are Here, M i (Γ i ) are the masses (widths) of the physical mass eigenstates |M µ 1,2 .
It is interesting to see how the Equation (12) defines the mass and the lifetime differences. Since the first term in Equation (12) is defined by a local operator, its matrix element does not develop an absorptive part, so it contributes to m 12 , i.e., the mass difference. The second term contains bi-local contributions connected by physical intermediate states. This term has both real and imaginary parts and thus contributes to both m 12 and Γ 12 .
Noting that Γ is defined by the standard model decay rate of the muon, and x and y are driven by the lepton-flavor violating interactions, we should expect that both x, y 1. The time evolution of flavor eigenstates follows from Equation (10) [7,19,20], where the coefficients g ± (t) are defined as As x, y 1 we can expand Equation (17) in power series in x and y to obtain The most natural way to detect M µ − M µ oscillations experimentally is by producing M µ state and looking for the decay products of the CP-conjugated state M µ . Denoting an amplitude for the M µ decay into a final state f as A f = f |H|M µ and an amplitude for its decay into a CP-conjugated final state f as Af = f |H|M µ , we can write the time-dependent decay rate of M µ into the f , where N f is a phase-space factor and we defined the oscillation rate R M (x, y) as Integrating over time and normalizing to Γ(M µ → f ) we get the probability of M µ decaying as M µ at some time t > 0, The equation Equation (21) [7] generalizes oscillation probability found in the papers [12,14] by allowing for a non-zero lifetime difference in M µ − M µ oscillations. We will review how x and y are related to the fundamental parameters of the Lagrangian below [7].

The Mass Difference x
The physical mixing parameters x and y can be obtained from Equation (15). The mass difference x comes from the dispersive part of the correlator.
Neglecting the SM contribution to the local |M µ , so the dominant contribution is only suppressed by Λ 2 . We can neglect the contribution of the second term in Equation (22), as it is suppressed by either Λ 4 or by M 2 W Λ 2 . In the former case this suppression comes from the double insertion of the L ∆L µ =1 eff term, each of which is suppressed by Λ 2 , as follows from Equation (4). In the later, the suppression comes from the insertion of the SM L ∆L µ =0 eff and L ∆L µ =2 eff terms with the operators given in Equation (7).
Computation of the mass and the lifetime differences involves evaluating the matrix elements between the muonium states for both the spin-0 singlet and the spin-1 triplet configurations. Since M µ is a QED bound state, such matrix elements can be written in terms of the value of the muonium wave function at the origin.
In the non-relativistic approximation, which is applicable for the muonium, it is given by a Coulombic bound state the wave function of the ground state, where a M µ = (αm red ) −1 is the muonium Bohr radius, α is the fine structure constant, and m red = m e m µ /(m e + m µ ) is the reduced mass. The absolute value |ϕ(0)| at the origin can be written as where we substituted the value of the Bohr radius. The applicability of a non-relativistic approximation to muonium can be established from a simple scaling argument. The typical momentum in the muonium state is where, for a moment, we reinstated c andh. We can see from Equation (25) that v ∼ αc, justifying the non-relativistic approximation. However, it might be easier to apply a factorization approach familiar from the description of meson flavor oscillations. In this approach the matrix elements in Equation (22) are obtained by inserting a vacuum state to turn matrix elements of four-fermion operators M µ ... M µ into products M µ ...|0 0|... M µ of matrix elements of current operators. Such matrix elements can be further parameterized as where f M is the muonium decay constant [6,7], p α is muonium's four-momentum, and α (p) is the ortho-muonium's polarization vector. Note that f P = f V = f T = f M in the non-relativistic limit. The decay constant f M can be expressed in terms of the bound-state wave function using the QED version of Van Royen-Weisskopf formula, This factorization gives the exact result for the QED matrix elements of the six-fermion operators in the non-relativistic limit, as can be explicitly verified [7]. Note that muonium mass and lifetime differences, and decay probabilities are thus suppressed by |ϕ(0)| 2 ∼ m 3 e . Para-muonium. The matrix elements of the spin-singlet states can be obtained from Equation (6) using the definitions of Equation (26), Combining the contributions from the different operators and using the definitions from Equations (24) and (27), we obtain an expression for x P for the para-muonium state, Ortho-muonium. Computing the relevant matrix elements for the vector ortho-muonium state, we obtain the matrix elements Again, combining the contributions from the different operators, we obtain an expression for x V for the ortho-muonium state, The results in Equations (29) and (31) are universal and hold true for any new physics model that can be matched into a set of local ∆L = 2 interactions.

The Lifetime Difference y
The lifetime difference in the muonium system, defined in Equation (15), is obtained from the absorptive part of Equation (12) and comes from the on-shell intermediate states common to both M µ and M µ [21], where ρ n is a phase space function for a given intermediate state. in Equation (32). According to Equation (4), it appears that, quite generally, this contribution is suppressed by Λ 4 , i.e., will be much smaller than x. The decays of muonia to the γγ intermediate states are generated by higher-dimensional operators and therefore are suppressed by even higher powers of Λ or the QED coupling α than the contributions considered here.
The only other possible contribution to y comes from the on-shell νν intermediate state. This intermediate state can be reached by the standard model tree level decay M µ → ν µ ν e and the ∆L µ = 2 decay M µ → ν µ ν e , i.e., it is common for both M µ and M µ . This contribution is only suppressed by Λ 2 M 2 W and represents the parametrically leading contribution to y [7]. Writing y in terms of the absorptive part of the correlation function, where the H is given by the ordinary standard model Lagrangian of Equation (3), and H ∆L µ =2 eff only contributes through the operators Q 6 and Q 7 of Equation (7). Since the decaying muon injects a large momentum into the two-neutrino intermediate state, the integral in Equation (33) is dominated by small distance contributions, compared to the scale set by 1/m µ . We can compute the correlation function in Equation (33) by employing a short distance operator product expansion, systematically expanding it in powers of 1/m µ [7].
Using Cutkoski rules to compute the discontinuity (imaginary part) of the transition amplitude (see Figure 1), calculating the relevant phase space integrals, and taking the matrix elements for the spin-singlet and the spin-triplet states of the muonium we arrive at the lifetime differences for the two spin states [7]. Para-muonium. The relevant matrix elements of the spin-singlet state can be read off the Equation (28). Recalling the definitions in Equations (24) and (27), we obtain an expression for the lifetime difference y P for the para-muonium state, One should note that if C ∆L=2

6
= C ∆L=2 7 current conservation assures that no lifetime difference is generated at this order in 1/Λ for the para-muonium.

Ortho-muonium. Employing the matrix elements for the spin-triplet state computed in Equation (30), the lifetime difference for the vector muonia is
We emphasize that Equations (34) and (35) represent parametrically leading contributions to muonium lifetime difference, as they are only suppressed by two powers of Λ. Nevertheless, additional suppression by G F makes the observation of the lifetime difference in the muonium system extremely challenging.

Experimental Studies of Muonium Oscillations
Both x and y are the observable parameters, so experiments studying M µ − M µ oscillations could probe them by producing the M µ state and looking for the decay products of the M µ state. Such experiments must overcome several considerable challenges.
First, muonium states must be produced inside the target and moved to the decay volume of the detector. The efficiency of such a process strongly depends on the target material. For example, the Muonium-Antimuonium Conversion Spectrometer (MACS) at PSI [18] used SiO 2 powder target with efficiency of 4-5%. A recently proposed Muonium-to-Antimuonium Conversion Experiment (MACE) [22] at the China Spallation Neutron Source (CSNS) will use an aerogel target with laser-drilled channels. The efficiency of producing muonia on such target could reach 40%.
Second, the strategy of looking for the "wrong-sign" final state, while working well in the studies of B 0 or D 0 oscillations, faces new challenges when applied to muonium oscillation searchers. The main reason is that the final state in the muonium decay M µ → e + e −ν µ ν e and that in the antimuonium decay M µ → e + e −ν e ν µ , differ only by the flavor combination of the neutrino states. Since the neutrinos are not detected, a method based on the kinematics of the decay must be employed: the decay products of the M µ involve fast electron with the momentum of about 53 MeV and slow positron with the momentum of about 13.5 eV. This follows from the fact that the main decay channel of the state is the decay of a muon.
Third, the experiments usually involve a setup where produced muonia propagate in the magnetic field B 0 . This magnetic field suppresses oscillations by removing degeneracy between M µ and M µ [15,23]. It also has a different effect on different spin configurations of the muonium state and the Lorentz structure of the operators that generate mixing [24,25]. Fortunately, these effects can be taken into account. MACS experiment corrected the oscillation probability by introducing a factor S B (B 0 ) [18], The values of S B (B 0 ), presented in Table II of [18], are different for different values of magnetic field and chiral structure of the operators governing the M µ − M µ oscillations.
We can now use the derived expressions for x and y to place constraints on the BSM scale Λ (or the Wilson coefficients C i ) from the experimental constraints on muonium-anti-muoium oscillation parameters. Since both spin-0 and spin-1 muonium states were produced in the experiment [18], we should average the oscillation probability over the number of polarization degrees of freedom, where P(M µ → M µ ) exp is the experimental oscillation probability from Equation (36). We shall use the values of S B (B 0 ) for B 0 = 2.8 µT from the Table II of [18], as it will provide us the best experimental constraints on the BSM scale Λ. We set the corresponding Wilson coefficient C i = 1. We report those constraints in Table 1. Table 1. Constraints on the energy scales probed by different ∆L = 2 operators of Equations (6) and (7) in MACS experiment (from [7]).

Operator Interaction Type
As one can see from Equations (29), (31), (34) and (35), each observable depends on the combination of the operators. Assuming that only one operator at a time gives a dominant contribution, i.e., employing the single operator dominance hypothesis, it is possible to constrain the Wilson coefficients of each operator. While this ansatz is not necessarily realized in many particular UV completions of the LFV EFTs, as cancellations among contributions of different operators are possible, it is, however, a useful tool in constraining parameters of L eff .
We must emphasize that the constraints presented in Table 1 use the data obtained more than twenty years ago! New results from new experiments are therefore highly desired. The MACE experiment [22] is expected to improve the sensitivity to C ∆L µ =2 i /Λ 2 by at least two orders of magnitude. A new experiment at J-PARC is also expected to improve the constraints on the muonium-antimuonium oscillation parameters [26].

Constraints on Explicit Models of New Physics
It might be instructive to consider explicit models of new physics which could be probed by M µ − M µ oscillations. This approach lacks the universality of the EFT approach described above. However, what it lacks in the universality it compensates for in the applicability: new degrees of freedom introduced in explicit models can contribute to other processes, some of which might not even include FCNCs or muons. Even restricting our attention to the sector of those frameworks containing ∆L µ = 1 and/or ∆L µ = 2 interactions reveals such a multitude of models that reviewing each and every one of them here would be impractical. We will take another approach.
We will review two specific models, one with heavy NP particles (a doubly-charged Higgs model), and one with light NP states (a model with flavor-violating axion-like particle (ALPs)) as examples. Then we discuss classes of NP interactions that can be probed by muonium oscillations, and refer the readers to further examples.
One way to classify the models of NP that contribute to M µ − M µ mixing is by the masses of NP particles m. The models with heavy (m m M µ ) new physics degrees of freedom can be matched to the effective Lagrangian of Equation (5). With that, experimental constraints on the parameters of this Lagrangian discussed in Section 3.4 and Table 1 lead to constraints on model parameters. The constraints on models with light (m ≤ m M µ ) new physics degrees of freedom can be implemented within the framework of either concrete or simplified models.

Heavy new physics.
Let us consider a model which contains a doubly-charged Higgs boson [27][28][29]. Such states often appear in the context of left-right models [30,31], where an additional Higgs triplet is introduced to introduce neutrino masses A coupling of the doubly charged Higgs field ∆ ++ to the lepton fields can be written as where c = C T is the charge-conjugated lepton state. Integrating out the ∆ −− field, this Lagrangian leads to the following effective Hamiltonian [27,31] below the scales associated with the doubly-charged Higgs field's mass M ∆ . Examining Equation (40) we see that this Hamiltonian matches onto our operator Q 2 (see Equation (6)) with the scale Λ = M ∆ and the corresponding Wilson coefficient C ∆L=2 2 = g ee g * µµ /8. The constraints on the masses and coupling constants depend on a particular model and other assumptions, such as whether the hierarchy of the neutrino masses is direct or inverse. It is however claimed that future experiments, such as MACE, could provide constraints on the doubly-charged Higgs state mass of m ∆ < 3 TeV [29].
It is interesting to point out that such a doubly-charged scalar state would also contribute to the anomalous magnetic moment of the muon, both at one-loop [32] and at two-loops via the Barr-Zee type of mechanism [33]. It was shown that ∆ ++ with the mass of a few hundred GeV could explain the discrepancy between the theoretical prediction and experimental measurement of (g − 2) of the muon [34], particularly due to the enhancement from the two units of electric charge of the ∆ ++ .
Light new physics. Let us consider a model with light axion-like particles that couple derivatively to the lepton current [16,35]. For the flavor off-diagonal interactions of the ALP a, the Lagrangian would contain a term where C V,A i j are Hermitian matrices of coupling constants, and f a is a decay constant related to the scale of the symmetry breaking associated with the ALP. Since the mass of the ALP, m a , is a free parameter, it is possible that, while a is light, to have m a > m M µ . In this case the ALP can be integrated out and the constraints from the Table 1 imply that [35] 1 3.8 TeV In the opposite case m a < m M µ the constraint can be obtained by taking a limit m µ /m a → 1 [35]. Note that in both cases one can only constrain the combination C V,A µe / f a . These results have direct implications for ALP contributions to muon (g − 2) [16] (c.f., [36]), the cross-section of e + e − → eeµµ [35], and other quantities.
Both decays contribute to the "invisible" width of the muonium, where ellipses denote decays into multineutrino states, which are further suppressed by the phase space factors. Experimental studies of the invisible width of muonia are very challenging. One possibility to constrain the invisible width indirectly [40] involves comparing the measurements of the positively-charged muon lifetime inside the material, where it forms the muonia and that of the free muon. Since the lifetime difference between the muonium in the target and free muon in vacuum was estimated to be negligible [10], the possible difference can be ascribed to the invisible width of the muonium. Using the lifetime measurements performed by the MuLan Collaboration [41], the limit B M µ → invisible < 5.7 × 10 −6 (47) can be established at 90% C.L. [40] with an assumption that the fraction of triplet muonium state in the MuLan's quartz target is 3/4. This bound does not yet provide a meaningful probe of either an SM process of Equation (44) or the ∆L µ = 2 process of Equation (43). A more sensitive probe of the ∆L µ = 2 process does not require a muonium bound state and can be obtained by comparing the lifetime of the muon in a vacuum to theoretical prediction in the framework of the standard model. The results will be reported elsewhere.

Conclusions
Muonium is the simplest QED bound state. Yet, it holds the potential to probe new physics in ways that are unique and complementary to other NP searches with muons. In this review, we discussed how lepton-flavor violating ∆L µ = 2 transitions could be probed in the muonium system, both in muonium decays and in muonium-antimuonium oscillations. We discussed the phenomenology of such oscillations, arguing that the presence of ∆L µ = 2 interactions would lead to both mass and lifetime splittings in the muonium-antimuonium system. We showed how to compute those oscillation parameters in an effective field theory approach, showing that both parameters scale as 1/Λ 2 with respect to the new physics scale Λ.
We discussed the current status and future perspectives of experiments that can constrain muonium-antimuonium oscillation parameters, and how such measurements can be used to place bounds on the parameters of the effective Lagrangians governing the ∆L µ = 2 transitions. These bounds can be further translated into the constraints on the masses and couplings of new particles in specific models of new physics, and, in case of the observation of the oscillation phenomena would help us to identify the proper extension of the standard model.