Studying ΔL = 2 Lepton Flavor Violation with Muons

Abstract

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.

1. 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),

$${\mathcal{L}}_{\mathrm{SMEFT}}=-\frac{1}{{\Lambda }^2}\sum _i{c}_i\left(\mu \right)Q_i,$$

(1)

where the Wilson coefficients (WC) $c_i$ of ${\mathcal{L}}_{\mathrm{SMEFT}}$ in Equation (1) are determined by the UV physics that becomes active at some scale $\Lambda$ [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 ${\mathcal{L}}_{\mathrm{eff}}$ only containing these degrees of freedom. In addition, it is convenient to classify the operators of ${\mathcal{L}}_{\mathrm{eff}}$ by their lepton quantum numbers $L_{\ell }$ (with $\ell =e,\mu ,\tau$),

$${\mathcal{L}}_{\mathrm{eff}}={\mathcal{L}}_{\mathrm{eff}}^{\Delta L_{\mu }=0}+{\mathcal{L}}_{\mathrm{eff}}^{\Delta L_{\mu }=1}+{\mathcal{L}}_{\mathrm{eff}}^{\Delta L_{\mu }=2}$$

(2)

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,

$${\mathcal{L}}_{\mathrm{eff}}^{\Delta L_{\mu }=0}=-\frac{4G_F}{\sqrt{2}}\left({\overline{\mu }}_L{\gamma }_{\alpha }{e}_L\right)\left({\overline{{\nu }_e}}_L{\gamma }^{\alpha }{{\nu }_{\mu }}_L\right).$$

(3)

Since $\Lambda \gg m_W$, which is required for the consistency of SM EFT, we can neglect the NP piece of the $\Delta L_{\mu }=0$ Lagrangian for studies of FCNC transitions. This implies that ${\mathcal{L}}_{\mathrm{eff}}^{\Delta L_{\mu }=0}$ is suppressed by powers of $G_F\sim M_W^{-2}$, 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 $\Delta L_{\mu }=1$ operators. The most general low energy Lagrangian representing muon flavor change by one unit, ${\mathcal{L}}_{\mathrm{eff}}^{\Delta L_{\mu }=1}$, is conventionnaly written as [4,5,6]

$$\begin{array}{ccc}\hfill {\mathcal{L}}_{\mathrm{eff}}^{\Delta L_{\mu }=1}=& -& \frac{1}{{\Lambda }^2}\sum _f[\left(C_{VR}^f{\overline{\mu }}_R{\gamma }^{\alpha }{e}_R+C_{VL}^f{\overline{\mu }}_L{\gamma }^{\alpha }{e}_L\right)\overline{f}{\gamma }_{\alpha }f\hfill \\ & +& \left(C_{AR}^f{\overline{\mu }}_R{\gamma }^{\alpha }{e}_R+C_{AL}^q{\overline{\mu }}_L{\gamma }^{\alpha }{e}_L\right)\overline{f}{\gamma }_{\alpha }{\gamma }_{5}f\hfill \\ & +& m_e{m}_f{G}_F\left(C_{SR}^f{\overline{\mu }}_R{e}_L+C_{SL}^f{\overline{\mu }}_L{e}_R\right)\overline{f}f\hfill \\ & +& m_e{m}_f{G}_F\left(C_{PR}^f{\overline{\mu }}_R{e}_L+C_{PL}^f{\overline{\mu }}_L{e}_R\right)\overline{f}{\gamma }_{5}f\hfill \\ & +& m_e{m}_f{G}_F\left(C_{TR}^f{\overline{\mu }}_R{\sigma }^{\alpha \beta }{e}_L+C_{TL}^f{\overline{\mu }}_L{\sigma }^{\alpha \beta }{e}_R\right)\overline{f}{\sigma }_{\alpha \beta }f+h.c.],\hfill \end{array}$$

(4)

where $\mu$ and e are the fermion fields with ${(\mu ,e)}_{L,R}=P_{L,R}(\mu ,e)$, where $P_{R,L}=\frac{1}{2}\left(1\pm {\gamma }^5\right)$ 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_{IK}^f$ 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), ${\mathcal{L}}_{\mathrm{eff}}^{\Delta L_{\mu }=2}$, 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

$${\mathcal{L}}_{\mathrm{eff}}^{\Delta L_{\mu }=2}=-\frac{1}{{\Lambda }^2}\sum _i{C}_i^{\Delta L=2}\left(\mu \right)Q_i\left(\mu \right).$$

(5)

with the operators written entirely in terms of the muon and electron degrees of freedom,

$$\begin{array}{ccc}\hfill Q_1& =& \left({\overline{\mu }}_L{\gamma }_{\alpha }{e}_L\right)\left({\overline{\mu }}_L{\gamma }^{\alpha }{e}_L\right),Q_2=\left({\overline{\mu }}_R{\gamma }_{\alpha }{e}_R\right)\left({\overline{\mu }}_R{\gamma }^{\alpha }{e}_R\right),\hfill \\ \hfill Q_3& =& \left({\overline{\mu }}_L{\gamma }_{\alpha }{e}_L\right)\left({\overline{\mu }}_R{\gamma }^{\alpha }{e}_R\right),Q_4=\left({\overline{\mu }}_L{e}_R\right)\left({\overline{\mu }}_L{e}_R\right),\hfill \\ \hfill Q_5& =& \left({\overline{\mu }}_R{e}_L\right)\left({\overline{\mu }}_R{e}_L\right).\hfill \end{array}$$

(6)

Other Lorentz structures of the operators could be related to the ones in Equation (6) via Fierz relations. A $\Delta L=2$ interaction described by this Lagrangian can lead to muon decays such as ${\mu }^{+}\to 3e$. It can also change a ${\mu }^{+}{e}^{-}$ bound state into a ${\mu }^{-}{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 $\Delta L_{\mu }=2$ local operators that contain muon, electron, and neutrino fields [7],

$$Q_6=\left({\overline{\mu }}_L{\gamma }_{\alpha }{e}_L\right)\left({\overline{{\nu }_{\mu }}}_L{\gamma }^{\alpha }{{\nu }_e}_L\right),Q_7=\left({\overline{\mu }}_R{\gamma }_{\alpha }{e}_R\right)\left({\overline{{\nu }_{\mu }}}_L{\gamma }^{\alpha }{{\nu }_e}_L\right),$$

(7)

where we only included operators that contain left-handed neutrinos [3,8]. The operators in Equation (7) could lead to both muon decays $\mu \to e{\nu }_{\mu }{\overline{\nu }}_e$ (note the “reversed” neutrino flavors), and muonium–antimuonium oscillations.

The goal of experimental studies of $\Delta 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.

2. 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_{\mu }$. The muonium is a QED bound state of a positively-charged muon and a negatively-charged electron, $|M_{\mu }\rangle \equiv |{\mu }^{+}{e}^{-}\rangle$.

Since muon is unstable, $M_{\mu }$ is unstable. The main decay channel for the muonium is determined by the weak decay of the muon, $M_{\mu }\to e^{+}{e}^{-}{\overline{\nu }}_{\mu }{\nu }_e$, so the average lifetime of a muonium state ${\tau }_{M_{\mu }}$ is expected to be the same as that of the muon,

$$1/{\tau }_{M_{\mu }}=\Gamma \left(M_{\mu }\to e^{+}{e}^{-}{\nu }_e{\overline{\nu }}_{\mu }\right)\approx \Gamma \left(\mu \to e^{+}{\nu }_e{\overline{\nu }}_{\mu }\right)=\frac{G_F^2{m}_{\mu }^5}{192{\pi }^3}=1/{\tau }_{\mu },$$

(8)

with ${\tau }_{\mu }=(2.1969811\pm 0.0000022)\times {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 $\left|M_{\mu }^V\right.〉$ is called ortho-muonium, while a spin-zero singlet state $\left|M_{\mu }^P\right.〉$ called para-muonium. In what follows, we will drop the superscript and employ the notation $\left|M_{\mu }\right.〉$ if the spin of the muonium state is not important for the discussion.

3. Muonium–Antimuonium Oscillations

Since $\Delta 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].

3.1. Phenomenology of Muonium Oscillations

In order to determine experimental observables related to $M_{\mu }-{\overline{M}}_{\mu }$ 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_{\mu }$ and ${\overline{M}}_{\mu }$, the time development of a muonium and anti-muonium states would be coupled, so it would be appropriate to consider their combined evolution,

$$|\psi \left(t\right)\rangle =\left(\begin{array}{c}a\left(t\right)\\ b\left(t\right)\end{array}\right)=a\left(t\right)|M_{\mu }\rangle +b\left(t\right)|{\overline{M}}_{\mu }\rangle .$$

(9)

The time evolution of $\left|\psi \right(t)\rangle$ evolution is governed by a Schrödinger-like equation,

$$i\frac{d}{dt}\left(\begin{array}{c}\left|M_{\mu }\left(t\right)\right.〉\\ \left|{\overline{M}}_{\mu }\left(t\right)\right.〉\end{array}\right)=\left(m-i\frac{\Gamma }{2}\right)\left(\begin{array}{c}\left|M_{\mu }\left(t\right)\right.〉\\ \left|{\overline{M}}_{\mu }\left(t\right)\right.〉\end{array}\right).$$

(10)

where ${\left(m-i\frac{\Gamma }{2}\right)}_{ik}$ is a $2\times 2$ Hamiltonian (mass matrix) with non-zero off-diagonal terms originating from the $\Delta 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}$, ${\Gamma }_{11}={\Gamma }_{22}$. In what follows, we assume CP-invariance of the $\Delta L_{\mu }=2$ interaction1. Then,

$$\begin{array}{c}\hfill m_{12}=m_{21}^{*},{\Gamma }_{12}={\Gamma }_{21}^{*}.\end{array}$$

(11)

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],

$${\left(m-\frac{i}{2}\Gamma \right)}_{12}=\frac{1}{2M_M}〈{\overline{M}}_{\mu }\left|{\mathcal{H}}_{\mathrm{eff}}\right|M_{\mu }〉+\frac{1}{2M_M}\sum _n\frac{〈{\overline{M}}_{\mu }\left|{\mathcal{H}}_{\mathrm{eff}}\right|n〉〈n\left|{\mathcal{H}}_{\mathrm{eff}}\right|M_{\mu }〉}{M_M-E_n+iϵ}.$$

(12)

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_{\mu }}_{1,2}\rangle$, which are related to the flavor eigenstates $M_{\mu }$ and ${\overline{M}}_{\mu }$ as

$$|{M_{\mu }}_{1,2}\rangle =\frac{1}{\sqrt{2}}\left[|M_{\mu }\rangle \mp |{\overline{M}}_{\mu }\rangle \right],$$

(13)

where we employed a convention where $CP|{M_{\mu }}_{\pm }\rangle =\mp |{M_{\mu }}_{\pm }\rangle$. The mass and the width differences of the mass eigenstates are

$$\begin{array}{c}\hfill \Delta m\equiv M_1-M_2,\Delta \Gamma \equiv {\Gamma }_2-{\Gamma }_1.\end{array}$$

(14)

Here, $M_i$ (${\Gamma }_i$) are the masses (widths) of the physical mass eigenstates $|{M_{\mu }}_{1,2}\rangle$.

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 ${\Gamma }_{12}$.

It is often convenient to introduce dimensionless quantities,

$$x=\frac{\Delta m}{\Gamma },y=\frac{\Delta \Gamma }{2\Gamma },$$

(15)

where the average lifetime $\Gamma =({\Gamma }_1+{\Gamma }_2)/2$, and $M_M=\left(M_1+M_2\right)/2$ is the muonium mass. Noting that $\Gamma$ 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\ll 1$.

The time evolution of flavor eigenstates follows from Equation (10) [7,19,20],

$$\begin{array}{ccc}\hfill \left|M\left(t\right)\right.〉& =& g_{+}\left(t\right)\left|M_{\mu }\right.〉+g_{-}\left(t\right)\left|{\overline{M}}_{\mu }\right.〉,\hfill \\ \hfill \left|\overline{M}\left(t\right)\right.〉& =& g_{-}\left(t\right)\left|M_{\mu }\right.〉+g_{+}\left(t\right)\left|{\overline{M}}_{\mu }\right.〉,\hfill \end{array}$$

(16)

where the coefficients $g_{\pm }\left(t\right)$ are defined as

$$g_{\pm }\left(t\right)=\frac{1}{2}{e}^{-{\Gamma }_{1}t/2}{e}^{-i{M}_{1}t}\left[1\pm e^{\Delta \Gamma t/2}{e}^{i\Delta mt}\right].$$

(17)

As $x,y\ll 1$ we can expand Equation (17) in power series in x and y to obtain

$$\begin{array}{ccc}\hfill g_{+}\left(t\right)& =& e^{-{\Gamma }_{1}t/2}{e}^{-i{M}_{1}t}\left[1+\frac{1}{8}{\left(y-ix\right)}^2{\left(\Gamma t\right)}^2\right],\hfill \\ \hfill g_{-}\left(t\right)& =& \frac{1}{2}{e}^{-{\Gamma }_{1}t/2}{e}^{-i{M}_{1}t}\left(y-ix\right)\left(\Gamma t\right).\hfill \end{array}$$

(18)

The most natural way to detect $M_{\mu }-{\overline{M}}_{\mu }$ oscillations experimentally is by producing $M_{\mu }$ state and looking for the decay products of the CP-conjugated state ${\overline{M}}_{\mu }$. Denoting an amplitude for the $M_{\mu }$ decay into a final state f as $A_f=\langle f\left|\mathcal{H}\right|M_{\mu }\rangle$ and an amplitude for its decay into a CP-conjugated final state $\overline{f}$ as $A_{\overline{f}}=\langle \overline{f}|\mathcal{H}|M_{\mu }\rangle$, we can write the time-dependent decay rate of $M_{\mu }$ into the $\overline{f}$,

$$\Gamma (M_{\mu }\to \overline{f})\left(t\right)=\frac{1}{2}{N}_f{\left|A_f\right|}^2{e}^{-\Gamma t}{\left(\Gamma t\right)}^2{R}_M(x,y),$$

(19)

where $N_f$ is a phase-space factor and we defined the oscillation rate $R_M(x,y)$ as

$$R_M(x,y)=\frac{1}{2}\left(x^2+y^2\right).$$

(20)

Integrating over time and normalizing to $\Gamma (M_{\mu }\to f)$ we get the probability of $M_{\mu }$ decaying as ${\overline{M}}_{\mu }$ at some time $t>0$,

$$P(M_{\mu }\to {\overline{M}}_{\mu })=\frac{\Gamma (M_{\mu }\to \overline{f})}{\Gamma (M_{\mu }\to f)}=R_M(x,y).$$

(21)

The equation Equation (21) [7] generalizes oscillation probability found in the papers [12,14] by allowing for a non-zero lifetime difference in $M_{\mu }-{\overline{M}}_{\mu }$ oscillations. We will review how x and y are related to the fundamental parameters of the Lagrangian below [7].

3.2. 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.

$$x=\frac{1}{2M_M\Gamma }\mathrm{Re}\left[2\langle {\overline{M}}_{\mu }\left|{\mathcal{H}}_{\mathrm{eff}}\right|M_{\mu }\rangle +\langle {\overline{M}}_{\mu }\left|i\int d^{4}x\mathrm{T}\left[{\mathcal{H}}_{\mathrm{eff}}\left(x\right){\mathcal{H}}_{\mathrm{eff}}\left(0\right)\right]\right|M_{\mu }\rangle \right].$$

(22)

Neglecting the SM contribution to the local $\Delta L_{\mu }=2$ Hamiltonian, $\langle {\overline{M}}_{\mu }|{\mathcal{H}}_{\mathrm{eff}}|M_{\mu }\rangle =\langle {\overline{M}}_{\mu }|{\mathcal{H}}_{\mathrm{eff}}^{\Delta L_{\mu }=2}|M_{\mu }\rangle$, so the dominant contribution is only suppressed by ${\Lambda }^2$. We can neglect the contribution of the second term in Equation (22), as it is suppressed by either ${\Lambda }^4$ or by $M_W^2{\Lambda }^2$. In the former case this suppression comes from the double insertion of the ${\mathcal{L}}_{\mathrm{eff}}^{\Delta L_{\mu }=1}$ term, each of which is suppressed by ${\Lambda }^2$, as follows from Equation (4). In the later, the suppression comes from the insertion of the SM ${\mathcal{L}}_{\mathrm{eff}}^{\Delta L_{\mu }=0}$ and ${\mathcal{L}}_{\mathrm{eff}}^{\Delta L_{\mu }=2}$ 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_{\mu }$ 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,

$$\phi \left(r\right)=\frac{1}{\sqrt{\pi a_{M_{\mu }}^3}}{e}^{-\frac{r}{a_{M_{\mu }}}},$$

(23)

where $a_{M_{\mu }}={\left(\alpha m_{\mathrm{red}}\right)}^{-1}$ is the muonium Bohr radius, $\alpha$ is the fine structure constant, and $m_{red}=m_e{m}_{\mu }/(m_e+m_{\mu })$ is the reduced mass. The absolute value $\left|\phi \right(0\left)\right|$ at the origin can be written as

$${\left|\phi \left(0\right)\right|}^2=\frac{{\left(m_{red}\alpha \right)}^3}{\pi }=\frac{1}{\pi }{\left(m_{red}\alpha \right)}^3,$$

(24)

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

$$p\simeq \frac{\hslash }{a_{M_{\mu }}}=\frac{m_{\mathrm{red}}{e}^2}{\hslash }=m_{\mathrm{red}}\left(\frac{e^2}{\hslash c}\right)c\equiv m_{\mathrm{red}}v$$

(25)

where, for a moment, we reinstated c and ℏ. We can see from Equation (25) that $v\sim \alpha 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 $\langle {\overline{M}}_{\mu }|...|M_{\mu }\rangle$ into products $\langle {\overline{M}}_{\mu }|...\left|0\right.〉\left.〈0\right|...\left|M_{\mu }\right.〉$ of matrix elements of current operators. Such matrix elements can be further parameterized as

$$\begin{array}{ccc}\hfill \left.〈0\right|\overline{\mu }{\gamma }^{\alpha }{\gamma }^{5}e|M_{\mu }^P\rangle & =& i{f}_P{p}^{\alpha },\left.〈0\right|\overline{\mu }{\gamma }^{\alpha }e|M_{\mu }^V\rangle =f_V{M}_M{ϵ}^{\alpha }\left(p\right),\hfill \\ \hfill \left.〈0\right|\overline{\mu }{\sigma }^{\alpha \beta }e|M_{\mu }^V\rangle & =& i{f}_T\left({ϵ}^{\alpha }{p}^{\beta }-{ϵ}^{\beta }{p}^{\alpha }\right),\hfill \end{array}$$

(26)

where $f_M$ is the muonium decay constant [6,7], $p^{\alpha }$ is muonium’s four-momentum, and ${ϵ}^{\alpha }\left(p\right)$ 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,

$$\begin{array}{c}\hfill f_M^2=4\frac{{\left|\phi \left(0\right)\right|}^2}{M_M}.\end{array}$$

(27)

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 ${\left|\phi \left(0\right)\right|}^2\sim m_e^3$.

Para-muonium. The matrix elements of the spin-singlet states can be obtained from Equation (6) using the definitions of Equation (26),

$$\begin{array}{ccc}\hfill \langle {\overline{M}}_{\mu }^P|Q_1\left|M_{\mu }^P\right.〉& =& f_M^2{M}_M^2,\langle {\overline{M}}_{\mu }^P|Q_2\left|M_{\mu }^P\right.〉=f_M^2{M}_M^2,\hfill \\ \hfill \langle {\overline{M}}_{\mu }^P|Q_3\left|M_{\mu }^P\right.〉& =& -\frac{3}{2}{f}_M^2{M}_M^2,\langle {\overline{M}}_{\mu }^P|Q_4\left|M_{\mu }^P\right.〉=-\frac{1}{4}{f}_M^2{M}_M^2,\hfill \\ \hfill \langle {\overline{M}}_{\mu }^P|Q_5\left|M_{\mu }^P\right.〉& =& -\frac{1}{4}{f}_M^2{M}_M^2.\hfill \end{array}$$

(28)

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,

$$x_P=\frac{4{\left(m_{red}\alpha \right)}^3}{\pi {\Lambda }^2\Gamma }\left[C_1^{\Delta L=2}+C_2^{\Delta L=2}-\frac{3}{2}{C}_3^{\Delta L=2}-\frac{1}{4}\left(C_4^{\Delta L=2}+C_5^{\Delta L=2}\right)\right].$$

(29)

Ortho-muonium. Computing the relevant matrix elements for the vector ortho-muonium state, we obtain the matrix elements

$$\begin{array}{ccc}\hfill \langle {\overline{M}}_{\mu }^V|Q_1\left|M_{\mu }^V\right.〉& =& -3f_M^2{M}_M^2,\langle {\overline{M}}_{\mu }^V|Q_2\left|M_{\mu }^V\right.〉=-3f_M^2{M}_M^2,\hfill \\ \hfill \langle {\overline{M}}_{\mu }^V|Q_3\left|M_{\mu }^V\right.〉& =& -\frac{3}{2}{f}_M^2{M}_M^2,\langle {\overline{M}}_{\mu }^V|Q_4\left|M_{\mu }^V\right.〉=-\frac{3}{4}{f}_M^2{M}_M^2,\hfill \\ \hfill \langle {\overline{M}}_{\mu }^V|Q_5\left|M_{\mu }^V\right.〉& =& -\frac{3}{4}{f}_M^2{M}_M^2.\hfill \end{array}$$

(30)

Again, combining the contributions from the different operators, we obtain an expression for $x_V$ for the ortho-muonium state,

$$x_V=-\frac{12{\left(m_{red}\alpha \right)}^3}{\pi {\Lambda }^2\Gamma }\left[C_1^{\Delta L=2}+C_2^{\Delta L=2}+\frac{1}{2}{C}_3^{\Delta L=2}+\frac{1}{4}\left(C_4^{\Delta L=2}+C_5^{\Delta L=2}\right)\right].$$

(31)

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 $\Delta L=2$ interactions.

3.3. 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_{\mu }$ and ${\overline{M}}_{\mu }$ [21],

$$y=\frac{1}{\Gamma }\sum _n{\rho }_n〈{\overline{M}}_{\mu }\left|{\mathcal{H}}_{\mathrm{eff}}\right|n〉〈n\left|{\mathcal{H}}_{\mathrm{eff}}\right|M_{\mu }〉,$$

(32)

where ${\rho }_n$ is a phase space function for a given intermediate state. There are many possible intermediate states composed of electrons, photons, and neutrinos. However, only the intermediate state containing neutrinos gives the largest contribution. This follows from the following argument. Noting that the contributions of multibody intermediate states are suppressed by the phase space factors ${\rho }_n$ for $n>2$, only two body intermediate states need to be considered. All possible SM two-body intermediate states that can contribute to y are $e^{+}{e}^{-}$, $\gamma \gamma$, and $\nu \overline{\nu }$.

The $e^{+}{e}^{-}$ intermediate state corresponds to a $\Delta L_{\mu }=1$ decay $M_{\mu }\to e^{+}{e}^{-}$, which implies that ${\mathcal{H}}_{\mathrm{eff}}={\mathcal{H}}_{\mathrm{eff}}^{\Delta L_{\mu }=1}$ in Equation (32). According to Equation (4), it appears that, quite generally, this contribution is suppressed by ${\Lambda }^4$, i.e., will be much smaller than x. The decays of muonia to the $\gamma \gamma$ intermediate states are generated by higher-dimensional operators and therefore are suppressed by even higher powers of $\Lambda$ or the QED coupling $\alpha$ than the contributions considered here.

The only other possible contribution to y comes from the on-shell $\nu \overline{\nu }$ intermediate state. This intermediate state can be reached by the standard model tree level decay $M_{\mu }\to \overline{{\nu }_{\mu }}{\nu }_e$ and the $\Delta L_{\mu }=2$ decay ${\overline{M}}_{\mu }\to \overline{{\nu }_{\mu }}{\nu }_e$, i.e., it is common for both $M_{\mu }$ and ${\overline{M}}_{\mu }$. This contribution is only suppressed by ${\Lambda }^2{M}_W^2$ and represents the parametrically leading contribution to y [7].

Writing y in terms of the absorptive part of the correlation function,

$$\begin{array}{ccc}\hfill y& =& \frac{1}{2M_M\Gamma }\mathrm{Im}\left[\langle {\overline{M}}_{\mu }\left|i\int d^{4}x\mathrm{T}\left[{\mathcal{H}}_{\mathrm{eff}}\left(x\right){\mathcal{H}}_{\mathrm{eff}}\left(0\right)\right]\right|M_{\mu }\rangle \right]\hfill \\ & =& \frac{1}{M_M\Gamma }\mathrm{Im}\left[\langle {\overline{M}}_{\mu }\left|i\int d^{4}x\mathrm{T}\left[{\mathcal{H}}_{\mathrm{eff}}^{\Delta L_{\mu }=2}\left(x\right){\mathcal{H}}_{\mathrm{eff}}^{\Delta L_{\mu }=0}\left(0\right)\right]\right|M_{\mu }\rangle \right],\hfill \end{array}$$

(33)

where the ${\mathcal{H}}_{\mathrm{eff}}^{\Delta L_{\mu }=0}=-{\mathcal{L}}_{\mathrm{eff}}^{\Delta L_{\mu }=0}$ is given by the ordinary standard model Lagrangian of Equation (3), and ${\mathcal{H}}_{\mathrm{eff}}^{\Delta L_{\mu }=2}$ 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_{\mu }$. 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_{\mu }$ [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,

$$y_P=\frac{G_F}{\sqrt{2}{\Lambda }^2}\frac{M_M^2}{{\pi }^2\Gamma }{\left(m_{red}\alpha \right)}^3\left(C_6^{\Delta L=2}-C_7^{\Delta L=2}\right).$$

(34)

One should note that if $C_6^{\Delta L=2}=C_7^{\Delta L=2}$ current conservation assures that no lifetime difference is generated at this order in $1/\Lambda$ 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

$$y_V=-\frac{G_F}{\sqrt{2}{\Lambda }^2}\frac{M_M^2}{{\pi }^2\Gamma }{\left(m_{red}\alpha \right)}^3\left(5C_6^{\Delta L=2}+C_7^{\Delta L=2}\right).$$

(35)

We emphasize that Equations (34) and (35) represent parametrically leading contributions to muonium lifetime difference, as they are only suppressed by two powers of $\Lambda$. Nevertheless, additional suppression by $G_F$ makes the observation of the lifetime difference in the muonium system extremely challenging.

3.4. Experimental Studies of Muonium Oscillations

Both x and y are the observable parameters, so experiments studying $M_{\mu }-{\overline{M}}_{\mu }$ oscillations could probe them by producing the $M_{\mu }$ state and looking for the decay products of the ${\overline{M}}_{\mu }$ 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_{\mu }\to e^{+}{e}^{-}{\overline{\nu }}_{\mu }{\nu }_e$ and that in the antimuonium decay ${\overline{M}}_{\mu }\to e^{+}{e}^{-}{\overline{\nu }}_e{\nu }_{\mu }$, 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 ${\overline{M}}_{\mu }$ 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_{\mu }$ and ${\overline{M}}_{\mu }$ [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\left(B_0\right)$ [18],

$$P(M_{\mu }\to {\overline{M}}_{\mu })\le 8.3\times {10}^{-11}/S_B\left(B_0\right).$$

(36)

The values of $S_B\left(B_0\right)$, presented in Table II of [18], are different for different values of magnetic field and chiral structure of the operators governing the $M_{\mu }-{\overline{M}}_{\mu }$ oscillations.

We can now use the derived expressions for x and y to place constraints on the BSM scale $\Lambda$ (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,

$$P{(M_{\mu }\to {\overline{M}}_{\mu })}_{\exp }=\sum _{i=P,V}\frac{1}{2S_i+1}P({M_{\mu }}^i\to {{\overline{M}}_{\mu }}^i),$$

(37)

where $P{(M_{\mu }\to {\overline{M}}_{\mu })}_{\exp }$ is the experimental oscillation probability from Equation (36). We shall use the values of $S_B\left(B_0\right)$ for $B_0=2.8$$\mu$T from the Table II of [18], as it will provide us the best experimental constraints on the BSM scale $\Lambda$. We set the corresponding Wilson coefficient $C_i=1$. We report those constraints in

Table 1. Constraints on the energy scales probed by different $\Delta L=2$ operators of Equations (6) and (7) in MACS experiment (from [7]).

Operator	Interaction Type	${\mathit{S}}_{\mathit{B}}\left({\mathit{B}}_0\right)$ (from [18])	Scale $\mathbf{\Lambda }$, TeV
$Q_1$	$(V-A)\times (V-A)$	0.75	$5.4$
$Q_2$	$(V+A)\times (V+A)$	0.75	$5.4$
$Q_3$	$(V-A)\times (V+A)$	0.95	$5.4$
$Q_4$	$(S+P)\times (S+P)$	0.75	$2.7$
$Q_5$	$(S-P)\times (S-P)$	0.75	$2.7$
$Q_6$	$(V-A)\times (V-A)$	0.75	$0.58\times {10}^{-3}$
$Q_7$	$(V+A)\times (V-A)$	0.95	$0.38\times {10}^{-3}$

.

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 ${\mathcal{L}}_{\mathrm{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_i^{\Delta L_{\mu }=2}/{\Lambda }^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].

3.5. Constraints on Explicit Models of New Physics

It might be instructive to consider explicit models of new physics which could be probed by $M_{\mu }-{\overline{M}}_{\mu }$ 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 $\Delta L_{\mu }=1$ and/or $\Delta L_{\mu }=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_{\mu }-{\overline{M}}_{\mu }$ mixing is by the masses of NP particles m. The models with heavy ($m\gg m_{M_{\mu }}$) 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\le m_{M_{\mu }}$) 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

$$\Delta =\left(\begin{array}{cc}{\Delta }^{+}/\sqrt{2}& {\Delta }^{++}\\ {\Delta }^0& -{\Delta }^{+}/\sqrt{2}\end{array}\right)$$

(38)

A coupling of the doubly charged Higgs field ${\Delta }^{++}$ to the lepton fields can be written as

$${\mathcal{L}}_R=-\frac{1}{2}{g}_{{\ell }_i{\ell }_j}{\overline{\ell }}_i^c\Delta {\ell }_j+H.c.,$$

(39)

where ${\ell }^c=C{\overline{\ell }}^T$ is the charge-conjugated lepton state. Integrating out the ${\Delta }^{--}$ field, this Lagrangian leads to the following effective Hamiltonian [27,31]

$${\mathcal{H}}_{\Delta }=-\frac{g_{ee}{g}_{\mu \mu }^{*}}{8M_{\Delta }^2}\left({\overline{\mu }}_L{\gamma }_{\alpha }{e}_L\right)\left({\overline{\mu }}_L{\gamma }^{\alpha }{e}_L\right)+H.c.,$$

(40)

below the scales associated with the doubly-charged Higgs field’s mass $M_{\Delta }$. Examining Equation (40) we see that this Hamiltonian matches onto our operator $Q_2$ (see Equation (6)) with the scale $\Lambda =M_{\Delta }$ and the corresponding Wilson coefficient $C_2^{\Delta L=2}=g_{ee}{g}_{\mu \mu }^{*}/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_{\Delta }<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 ${\Delta }^{++}$ 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 ${\Delta }^{++}$.

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

$${\mathcal{L}}_{\mathrm{ALP}}=\sum _{i\ne j}\frac{{\partial }_{\mu }a}{2f_a}{\overline{\ell }}_i{\gamma }^{\mu }\left(C_{{\ell }_i{\ell }_j}^V+C_{{\ell }_i{\ell }_j}^A{\gamma }_5\right){\ell }_j,$$

(41)

where $C_{{\ell }_i{\ell }_j}^{V,A}$ 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_{\mu }}$. In this case the ALP can be integrated out and the constraints from the Table 1 imply that [35]

$$\frac{1}{3.8\mathrm{TeV}}>\frac{1}{2f_a}{\left|{\left(C_{\mu e}^A\right)}^2-{\left(C_{\mu e}^V\right)}^2\right|}^{1/2}\left(\frac{m_{\mu }}{m_e}\right).$$

(42)

In the opposite case $m_a<m_{M_{\mu }}$ the constraint can be obtained by taking a limit $m_{\mu }/m_a\to 1$ [35]. Note that in both cases one can only constrain the combination $C_{\mu e}^{V,A}/f_a$. These results have direct implications for ALP contributions to muon $(g-2)$ [16] (c.f., [36]), the cross-section of $e^{+}{e}^{-}\to ee\mu \mu$ [35], and other quantities.

A useful classification of lepton-flavor-violating interactions probed by the muonium oscillations was given in [37]. Such interactions can be classified by the way they break the lepton flavor.

• $\Delta L_e=\Delta L_{\mu }=0$. In such models, the interactions do not violate lepton flavor quantum numbers. The muonium oscillations are induced at one-loop order by the mass terms of the fields that transform as singlets under the SM gauge group and have $L_{e,\mu }=\pm 2$. Examples of such models include constructions with Majorana neutrinos [13,37].
• $\Delta L_e=\pm 2, \Delta L_{\mu }=0$ and $\Delta L_e=0, \Delta L_{\mu }=\pm 2$. In such models, the lepton flavor quantum numbers are separately broken. The mediators of such interactions, sometimes called dilepton bosons, have electric charge $2e$ and lepton number equal to two. The mediators could be scalar or vector bosons. Examples of such models include the doubly-charged Higgs model considered above [24,25,27,28,29,32,33] and many other constructions [37]. Muonium oscillations can be generated at tree level in such models.
• $\Delta L_e=-\Delta L_{\mu }=\pm 1$. In such models, interaction terms violate both lepton flavor numbers. The mediators of such interactions are electrically neutral and could be both scalar and vector bosons. Examples of such models include models with flavor-violating Higgs boson [38], and many others [37]. Muonium oscillations can also be generated at the tree level in such models. These models can also be probed in muon conversion experiments unless the mediator is introduced such that it only interacts with the leptons.
• $\Delta L_e=\pm 1, \Delta L_{\mu }=0$ and $\Delta L_e=0, \Delta L_{\mu }=\pm 1$. In such models, effective operators mediating muonium oscillations are generated at one-loop order. These models can also be probed in other muon transitions, such as $\mu \to 3e$.

As was pointed out in [37] if a $Z_n$ discrete symmetry is imposed to suppress the $\Delta L_e-\Delta L_{\mu }=\pm 1$ interactions, while allowing for $\Delta L_e-\Delta L_{\mu }=\pm 2$ terms, the oscillation rates are not well-probed by other experiments and can be as large as allowed by the current experimental bound. Clearly, $M_{\mu }-{\overline{M}}_{\mu }$ oscillations probe a wide variety of NP models and serve as effective tools that are complementary to other searches.

4. Muonium Decays

Processes that change lepton flavor by two units can, in principle, also be studied in the decays of muons and muonia. The transitions into neutrino states, governed by the operators in Equation (7), fit the bill. The decay rate for such transition is

$$\Gamma \left(M_{\mu }^V\to {\overline{\nu }}_e{\nu }_{\mu }\right)=\frac{f_M^2{M}_M^3}{9\pi {\Lambda }^4}{\left|C_6^{\Delta L_{\mu }=2}+C_7^{\Delta L_{\mu }=2}\right|}^2.$$

(43)

Note that in the limit of massless neutrinos only the spin-one ortho-muonium state would have non-zero width. While the invisible decays of the spin-zero para-muonium would not be strictly forbidden in this limit, as they would be dominated by the four-neutrino final state [39], the decay rate of the para-muonium would be very small.

The decay rate of Equation (43) will add incoherently to the SM decay rate

$$\Gamma \left(M_{\mu }^V\to {\overline{\nu }}_{\mu }{\nu }_e\right)=\frac{G_F^2{f}_M^2{M}_M^3}{12\pi }.$$

(44)

Even in the SM, this is a very rare process. Numerically, the branching fraction can be computed to be

$$\mathcal{B}\left(M_{\mu }\to {\overline{\nu }}_{\mu }{\nu }_e\right)=6.6\times {10}^{-12}.$$

(45)

Both decays contribute to the “invisible” width of the muonium,

$$\Gamma \left(M_{\mu }^V\to \mathrm{invisible}\right)=\Gamma \left(M_{\mu }^V\to {\overline{\nu }}_{\mu }{\nu }_e\right)+ \Gamma \left(M_{\mu }^V\to {\overline{\nu }}_e{\nu }_{\mu }\right)+...,$$

(46)

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

$$\mathcal{B}\left(M_{\mu }\to \mathrm{invisible}\right)<5.7\times {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 $\Delta L_{\mu }=2$ process of Equation (43).

A more sensitive probe of the $\Delta L_{\mu }=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.

5. 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 $\Delta L_{\mu }=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 $\Delta L_{\mu }=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/{\Lambda }^2$ with respect to the new physics scale $\Lambda$.

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 $\Delta L_{\mu }=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.