Muonium-antimuonium oscillations in effective field theory

Flavor violating processes in the lepton sector have highly suppressed branching ratios in the standard model mainly due to the tiny neutrino mass. This means that observing lepton flavor violation (LFV) in the next round of experiments would constitute a clear indication of physics beyond the standard model (BSM). We revisit a discussion of one possible way to search for LFV, muonium-antimuonium oscillations. This process violates muon lepton number by two units and could be sensitive to the types of BSM physics that are not probed by other types of LFV processes. Using techniques of effective field theory, we calculate the mass and width differences of the mass eigenstates of muonium. We argue that its invisible decays give the parametrically leading contribution to the lifetime difference and put constraints on the scales of new physics probed by effective operators in muonium oscillations.


I. 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. An especially clean system to study BSM effects in lepton sector is muonium M µ , a QED bound state of a positively-charged muon and a negatively-charged electron, |M µ ≡ |µ + e − .
The main decay channel for both states is driven by the weak decay of the muon. The average lifetime of a muonium state τ Mµ is expected to be the same as that of the muon, τ µ = (2.1969811 ± 0.0000022) × 10 −6 s [1], apart from the tiny effect due to time dilation, (τ Mµ − τ µ )/τ µ = α 2 m 2 e /(2m 2 µ ) = 6 × 10 −10 [2]. Just like a positronium or a Hydrogen atom, muonium could be produced in two spin configurations, a spin-one triplet state called ortho-muonium, and a spin-zero singlet state called para-muonium. We shall denote the para-muonium state as M P µ and the ortho-muonium state as M V µ . If the spin of the state does not matter, we shall employ the notation |M µ .
So far, we have not yet observed FCNC in the charged lepton sector. This is because in the standard model with massive neutrinos the charged lepton flavor violating (CLFV) transitions are suppressed by the powers of m 2 ν /m 2 W , which renders the predictions for their transition rates vanishingly small, e.g. B(µ → eγ) νSM ∼ 10 −54 [3,4]. Yet, experimental analyses constantly push the bounds on the CLFV transitions. It might be that in some models of NP, such as a model with the doubly-charged Higgs particles [5][6][7][8], the effective ∆L = 2 transitions could occur at a rate that is not far below the sensitivity of currentlyoperating experiments. Alternatively, it might be that no term that changes the lepton flavor by two units is present in a BSM Lagrangian. But even in this case, a subsequent application of two ∆L = 1 interactions would also generate an effective ∆L = 2 interaction.
Such a ∆L = 2 interaction would then change the muonium state into the anti-muonium one, leading to the possibility of muonium-anti-muonium oscillations. As a variety of wellestablished new physics models contain ∆L = 2 interaction terms [3], observation of muonium converting into anti-muonium could then provide especially clean probes of new physics in the leptonic sector [4,9]. Theoretical analyses of conversion probability for such transi-tions have been actively studied, mainly using the framework of particular models [10][11][12][13]. It would be useful to perform a model-independent computation of the oscillation parameters using techniques of effective theory that includes all possible BSM models encoded in a few Wilson coefficients of effective operators. We do so in this paper, computing all relevant QED matrix elements. Finally, employing similar effective field theory (EFT) techniques for computation of the contributions that are non-local at the muon mass scale, we present them in terms of the series of local operators expanded in inverse powers of m µ [16,17].
In this paper we discuss the most general analysis of M µ −M µ oscillations in the framework of effective field theory. We review phenomenology of muonium oscillations in Sec. II, taking into account both mass and lifetime differences in the muonium system. We compute the mass and width differences in Sec. III. In Sec. IV we constrain the BSM scale Λ using experimental muonium-anti-muonium oscillation parameters. We conclude in Sec. V.
Appendix VI contains some details of calculations.

II. PHENOMENOLOGY OF MUONIUM OSCILLATIONS
Phenomenology of M µ − M µ oscillations is very similar to phenomenology of mesonantimeson oscillations [18,19]. There are, however, several important differences that we will emphasize below. One major difference is related to the fact that both ortho-and para-muonium can, in principle, oscillate. While most studies only considered muonium oscillations due to the BSM heavy states, below we also discuss the possibility of oscillations via the light states. Since such states can go on mass shell, these contributions would lead to the possibility of a lifetime difference in the M µ − M µ system.
If the new physics Lagrangian includes lepton-flavor violating interactions, 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 equation, CPT-invariance dictates that the masses and widths of muonium and anti-muonium are the same, m 11 = m 22 , Γ 11 = Γ 22 , while CP-invariance of the ∆L µ = 2 interaction, which we assume for simplicity, dictates that The presence of off-diagonal pieces in the mass matrix signals that it needs to be diagonalized.
The mass eigenstates |M µ 1,2 can be defined as where we neglected CP-violation and employed a convention where CP |M µ ± = ∓|M µ ± .
The mass and the width differences of the mass eigenstates are where M i (Γ i ) are the masses (widths) of the mass eigenstates |M µ 1,2 . We defined ∆m and ∆Γ to be either positive or negative, which is to be determined by experiment. It is often convenient to introduce dimensionless quantities, where the average lifetime Γ = (Γ 1 +Γ 2 )/2. It is important to note that while Γ is defined by the standard model decay rate of the muon, x and y are driven by the lepton-flavor violating interactions. It is then expected that both x, y 1.
The time evolution of flavor eigenstates follows from Eq. (2) [18,19], where the coefficients g ± (t) are defined as As x, y 1 we can expand Eq. (8) to get Denoting an amplitude for the muonium 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 R M (x, y) is the oscillation rate, Integrating over time and normalizing to Γ(M µ → f ) we get the probability of M µ decaying as M µ at some time t > 0, This equation generalizes oscillation probability computed in the classic papers [11,13] by accounting for the lifetime difference in the muonium system, making it dependent on both the normalized mass x and the lifetime y differences. We will compute those in the next section.
We shall use the data from the most recent experiment [9] in order to place constraints on the oscillation parameters. To do so, we have to account for the fact that the setup described in [9] had muonia propagating in a magnetic field B 0 . This magnetic field suppresses oscillations by removing degeneracy between M µ and M µ . It also has a different effect on different spin configurations of the muonium state and the Lorentz structure of the operators that generate mixing [14,15]. Experimentally these effects were accounted for by introducing a factor S B (B 0 ). The oscillation probability is then [9], We shall use different values of S B (B 0 ), presented in Table II of [9] when placing constraints on the Wilson coefficients of effective operators in the next section.

III. EFFECTIVE THEORY OF OSCILLATIONS
Muonium-anti-muonium oscillations could be effective probes of flavor-violating new physics in leptons. One of the issues is that at this point we do not know which particular model of new physics will provide the correct ultraviolet (UV) extension for the standard model. However, since the muonium mass is most likely much smaller than the new particle masses, it is not necessary to know it. Any new physics scenario which involves lepton flavor violating interactions can be matched to an effective Lagrangian, L eff , whose Wilson coefficients would be determined by the UV physics that becomes active at some scale Λ [20,21], where the c i 's are the short distance Wilson coefficients. They encode all model-specific information. Q i 's are the effective operators which reflect degrees of freedom relevant at the scale at which a given process takes place. If we assume that no new light particles (such as "dark photons" or axions) exist in the low energy spectrum, those operators would be written entirely in terms of the SM degrees of freedom. In the case at hand, all SM particles with masses larger than that of the muon should also be integrated out, leaving only muon, electron, photon, and neutrino degrees of freedom.
It would be convenient for us to classify effective operators in Eq. (14) by their lepton quantum numbers. In particular, we can write the effective Lagrangian as The first term in this expansion contains both the standard model and the new physics contributions. It then follows that the leading term in L where G F ∼ M −2 W is the Fermi constant, µ and e are the fermion fields, (µ, e) L,R = P L,R (µ, e). P R,L = 1 2 (1 ± γ 5 ) are the projection operators, and f represents other fermions that are not integrated out at the the muonium scale. The subscripts on the Wilson coefficients are for the type of Lorentz structure: vector, axial-vector, scalar, pseudo-scalar, and tensor. The Wilson coefficients would in general be different for different fermions f . Note that the Lagrangian Eq. (16) also contains terms that do not follow from the dimension six in the standard model effective field theory (SMEFT), but could be generated by higher order operators. This is taken into account by introducing mass and G F factors emulating such suppression [23,24].
The last term in Eq. (15), L ∆Lµ=2 eff , represents the effective operators changing the lepton quantum number by two units. The leading contribution to muonium oscillations is given by the dimension six operators. The most general effective Lagrangian can be written with the operators written entirely in terms of the muon and electron degrees of freedom, We did not include operators that could be related to the presented ones via Fierz relations.
It is important to note that some of the operators in Eq. (18) are not invariant under the SM gauge group SU (2) L × U (1). This means that they receive additional suppression, as they may be generated from the higher-dimensional operators in SMEFT [21].
Other ∆L µ = 2 local operators that will be important later in this paper can be written as where we only included SMEFT operators that contain left-handed neutrinos [21,22]. In order to see how these operators (and thus new physics) contribute to the mixing parameters, it is instructive to consider off-diagonal terms in the mass matrix [18]  We can rewrite Eq. (20) to extract the physical mixing parameters x and y of Eq. (6).
For the mass difference, Assuming the LFV NP is present, the dominant local contribution to x comes from the last term in Eq. (15), In order to evaluate the mass difference contribution, we need to take the matrix elements.
As explained in the Introduction, we expect that both spin-0 singlet and spin-1 triplet muonium states would undergo oscillations. The oscillation parameters would in general be different, as the matrix elements would differ for those two cases.
Using factorization approach familiar from the meson flavor oscillation, the matrix elements can be easily written in terms of the muonium decay constant f M [25,26].
where p α is para-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 can be written in terms of the bound-state wave function, which is the QED's version of Van Royen-Weisskopf formula. For a Coulombic bound state the wave function of the ground state is 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. Then, In the non-relativistic limit factorization gives the exact result for the QED matrix elements of the six-fermion operators. Nevertheless, we explicitly verified that this is indeed the case (see Appendix VI).
Para-muonium. The matrix elements of the spin-singlet states can be obtained from Eq. (18) using the definitions of Eq. (23), Combining the contributions from the different operators and using the definitions from Eqs. (24) and (26), we obtain an expression for x P for the para-muonium state, This result is universal and holds true for any new physics model that can be matched into a set of local ∆L = 2 interactions.
Ortho-muonium. Using the same procedure, but 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, Again, this result is universal and holds true for any new physics model that can be matched into a set of local ∆L = 2 interactions.
It might be instructive to present an example of a BSM model that can be matched into the effective Lagrangian of Eq. (17) and can be constrained from Eqs. (27,29). Let us consider a model which contains a doubly-charged Higgs boson [5,6]. Such states often appear in the context of left-right models [7,8]. 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 [5,8] below the scales associated with the doubly-charged Higgs field's mass M ∆ . Examining Eq. (32) we see that this Hamiltonian matches onto our operator Q 2 (see Eq. (18) The lifetime difference in the muonium system can be obtained from Eq. (20) [27]. It comes from the physical intermediate states, which is signified by the imaginary part in Eq. (20) and reads, Writing y similarly to x in Eq. (21), i.e. in terms of the correlation function, we obtain where the H ∆Lµ=0 eff = −L ∆Lµ=0 eff is given by the ordinary standard model Lagrangian, and H ∆Lµ=2 eff only contributes through the operators Q 6 and Q 7 .
Since the decaying muon injects a large momentum into the two-neutrino intermediate The leading term is obtained by contracting the neutrino fields in Eq. (36) into propagators, where S F (x) represents the propagator in coordinate representation. In what follows we will consider neutrinos to be Dirac fields for simplicity.
Using Cutkoski rules to compute the discontinuity (imaginary part) of T and calculating the phase space integrals we get We can now compute the lifetime difference y by using Eq. (34) and take the relevant matrix elements for the spin singlet and the spin triplet states of the muonium.
Para-muonium. The matrix elements of the spin-singlet state have been computed above and presented in Eq. (27). Computing the matrix elements in Eq. (34) using their definitions from Eqs. (24) and (26), we obtain an expression for the lifetime difference y P for the para-muonium state, It is interesting to 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. Similarly, using the matrix elements for the spin-triplet state computed in Eq. (29), the expression fo Eq. (38) leads to the lifetime difference We emphasize that Eqs. (39) and (40) represent parametrically leading contributions to muonium lifetime difference, as they are only suppressed by two powers of Λ.

IV. EXPERIMENTAL CONSTRAINTS
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 [9], we should average the oscillation probability over the number of polarization degres of freedom, Operator Interaction type S B (B 0 ) (from [9]) Constraints on the scale Λ, TeV  (18) and (19). We set the corresponding Wilson coefficient C i = 1.
where P (M µ → M µ ) exp is the experimental oscillation probability from Eq. (13). We shall use the values of S B (B 0 ) for B 0 = 2.8 µT from the Table II of [9], as it will provide us the best experimental constraints on the BSM scale Λ. We report those constraints in Table   I. As one can see from Eqs. Since it is the combination C i /Λ 2 that enters the theoretical predictions for x and y, one cannot separately measure C i and Λ. We choose to constrain the scale Λ that is probed by the corresponding operator and set the corresponding value of the Wilson coefficient C i to one.
The results are reported in Table I. As can be seen, the experimental data provide constraints on the scales comparable to those probed by the LHC program, except for Q 6 and Q 7 . The constraints on the lepton-flavor violating neutrino operators Q 6 and Q 7 are understandably weaker, as the lifetime difference is suppressed by a factor G F /Λ 2 , while the mass difference is only suppressed by a factor of 1/Λ 2 . We note that the best constraints on the oscillation parameters come from the data that is over 20 years old [9]. We urge our experimental colleagues to further study muonium-antimuonium oscillations.

V. CONCLUSIONS
Lepton flavor violating transitions provide a powerful engine for new physics searches.
In this work we revisited phenomenology of muonium-antimuonium oscillations. We argued that in generic models of new physics both mass and lifetime differences in the muonium system would contribute to the oscillation probability. We computed the normalized mass difference x in the muonium system with the most general set of effective operators for both spin-singlet and the spin-triplet muonium states. We setup a formalism for computing the lifetime differences and computed the parametrically leading contribution to y.

VI. APPENDIX
In this Appendix we show that the vacuum insertion approximation leads to the same answer as a direct computation of a four-fermion matrix element relevant for the muoniumanti-muonium oscillations. We shall show that by computing a matrix element of the Q 1 operator as an example. The matrix elements is defined as for both pseudoscalar and vector muonium states. In order to compute the matrix element in Eq. (42) we need to build the muonium states. We can employ the standard Bethe-Salpeter formalism. Since the muonium state is essentially a a nonrelativistic Coulomb bound state of a µ + and an e − , we can conventionally define it [28] This state is normalized as M µ (P)|M µ (P ) = 2E p (2π) 3 δ 3 (P − P ). The muonium state in Eq. (43) is projected from a two-particle state of a muon and an electron |p, p = |0 with the help of the Fourier transform of the spatial wave equation describing the bound state ϕ(p), We expand each electron and muon field in the operator of Eq. (42) as We will work in non-relativistic approximation and neglect the momentum dependence of the spinors, which are defined as Here ξ and η are the two-component spinors [28]. There are four ways to Wick contract the fields in the operator with those generating the state. Using anti-commutation relation where we indicated that the spinors still need to be projected onto the spin-triplet or or the spin-singlet states. This projection can be illustrated explicitly by considering the first term in Eq. (47), (uγ α P L v) (vγ α P L u), the rest can be computed in a complete analogy to that.
as the sum over polarizations is µ µ * = −3. Following the same procedure for the rest of the terms in Eq. (47) and using we get Q 1 for spin-0 and spin-1 which is identical to the definitions in Eq. (27) and (29), provided that the Van Royen-Weisskopf formula of Eq. (24) is used. The proof for the rest of the operators follows the same steps.