Discovering True Muonium in $K_{L}\rightarrow(\mu^{+}\mu^{-})\gamma$

Lepton universality predicts differences in electron and muon observables should occur only due to their mass difference. Measurements of $(g-2)_{\ell}$ Bennett et al. (2006), nuclear charge radii Antognini et al. (2013); Pohl et al. (2016), and rare meson decays Aaij et al. (2014); *Aaij:2015yra have shown hints of violations to this universality. The bound state of $(\mu^{+}\mu^{-})$, true muonium, presents a unique opportunity to study lepton universality in and beyond the Standard Model Tucker-Smith and Yavin (2011); *Lamm:2015gka; *Lamm:2016jim. To facilitate these studies, efforts are on-going to improve theoretical predictions Jentschura et al. (1997a); *Jentschura:1997ma; *PhysRevD.91.073008; *Lamm:2016vtf; *PhysRevA.94.032507; *Lamm:2017lib. Alas, true muonium remains undetected today.

Since the late 60’s, two broad categories of $(\mu^{+}\mu^{-})$ production methods have been discussed: particle collisions (fixed-target and collider) Bilenky et al. (1969); *Hughes:1971; *Moffat:1975uw; *Holvik:1986ty; *Ginzburg:1998df; *ArteagaRomero:2000yh; *Brodsky:2009gx; *Chen:2012ci, or through rare decays of mesons Nemenov (1972); *Vysotsky:1979nv; *Kozlov:1987ey; Malenfant (1987). Until recently, none have been attempted due to the low production rate ($\propto\alpha^{4}$). Currently, the Heavy Photon Search (HPS) Celentano (2014) experiment is searching for true muonium Banburski and Schuster (2012) via $e^{-}Z\rightarrow(\mu^{+}\mu^{-})X$. Another fixed-target experiment, but with a proton beam, DImeson Relativistic Atom Complex (DIRAC) Benelli (2012) studies the $(\pi^{+}\pi^{-})$ bound state and could look for $(\mu^{+}\mu^{-})$ in a upgraded run Chliapnikov (2014).

In recent years, a strong focus on rare kaon decays has developed in the search for new physics. The existing KOTO experiment at J-PARC Ahn et al. (2017) and proposed NA62-KLEVER at CERN Moulson (2017) hope to achieve sensitivities of $\mathcal{BR}\sim 10^{-13}$ allowing a 1% measurement of $\mathcal{BR}(K_{L}\rightarrow\pi^{0}\nu\nu)\sim 10^{-11}$. Malenfant was the first to propose $K_{L}$ as a source of $(\mu^{+}\mu^{-})$ Malenfant (1987). He estimated $\mathcal{BR}(K_{L}\rightarrow(\mu^{+}\mu^{-})\gamma)\sim 5\times 10^{-13}$ by approximating $F_{K_{L}\gamma\gamma*}(Q^{2}=4M_{\mu}^{2})\sim F_{K_{L}\gamma\gamma*}(0)$ where $Q^{2}$ is the off-shell photon invariant mass squared. This two-body decay is the reach of rare kaon decay searches and is an attractive process for discovering $(\mu^{+}\mu^{-})$. The decay has simple kinematics with a single, monochromatic photon (of $E_{\gamma}=203.6$ MeV if the $K_{L}$ is at rest) plus $(\mu^{+}\mu^{-})$ which could undergo a two-body dissociate or decay into two electrons (with $M_{\ell\ell}^{2}\sim 4M_{\mu}^{2}$).

Another motivation for the search for this rare decay is its unique dependence on the form factor. Previous extractions of the form factor relied upon radiative Dalitz decays, $K_{L}\rightarrow\ell^{+}\ell^{-}\gamma$, the most recent being from the KTEV collaboration Abouzaid et al. (2007); Alavi-Harati et al. (2001). In these analyses, the phenomenological form factor is integrated over bins in $Q^{2}$, and fit to differential cross section data. In contrast, the $(\mu^{+}\mu^{-})$ branching ratio gives the form factor at one $Q^{2}$ and fixes one-parameter form factors. Further, a measurement of $(\mu^{+}\mu^{-})$ would help to better understand the kaon form factor through a completely different set of systematic and statistical uncertainties to the existing measurements.

In this letter, we present the $\mathcal{BR}(K_{L}\rightarrow(\mu^{+}\mu^{-})\gamma)$ including full $\mathcal{O}(\alpha)$ radiative corrections and four different treatments of the form factor $F_{K_{L}\gamma\gamma*}(Q^{2})$, thereby avoiding Malenfant’s approximation. It is shown that the approximation underestimates the branching ratio by a model-dependent 15-60%. Possible discovery channels are discussed and brief comments on important backgrounds are made.

Following previous calculations for atomic decays of mesons Nemenov (1972); Vysotsky (1979); Kozlov (1988); Malenfant (1987), the branching ratio can be computed

where $\zeta(3)=\sum_{n}1/n^{3}$ arising from the sum over all allowed $(\mu^{+}\mu^{-})$ states, $z_{TM}=M_{TM}^{2}/M_{K}^{2}\approx 4M_{\mu}^{2}/M^{2}_{K}$, and $f(z)=F_{K_{L}\gamma\gamma^{*}}(z)/F_{K_{L}\gamma\gamma^{*}}(0)$. Previous computation of radiative corrections considered only the vacuum polarization from the flavor found in the final state Vysotsky (1979). We have computed the full $(\mu^{+}\mu^{-})$ results including the electronic, muonic, and hadronic vacuum polarization Ji and Lamm (2016) as well as the QED process $K_{L}\rightarrow\gamma^{*}(k)+\gamma^{*}(P_{K_{L}}-k)\rightarrow\gamma+\text{TM}$ demonstrated by Fig. 1 where $P_{K_{L}}$ is the four-momentum of the $K_{L}$. In this contribution, one should take the convolution of the QED amplitude with double-virtual-photon form factor $F_{K_{L}\gamma^{*}\gamma^{*}}(k^{2}/M_{K}^{2},(P_{I}-k)^{2}/M_{K}^{2})$. For our purpose, however, taking the form factor to be $F_{\gamma\gamma^{*}}(0,z_{TM})$ is a sufficient approximation as shown in Kampf et al. (2006). A similar calculation for positronium, where other lepton flavors and hadronic loop corrections are negligible, finds the $\displaystyle\frac{\alpha}{\pi}$ coefficient is $-52/9$.

$F_{K_{L}\gamma\gamma^{*}}(0)$ is fixed to the experimental value of $\mathcal{BR}(K_{L}\rightarrow\gamma\gamma)=5.47(4)\times 10^{-4}$ Patrignani et al. (2016). Evaluating Eq. (Discovering True Muonium in $K_{L}\rightarrow(\mu^{+}\mu^{-})\gamma$), we find $\mathcal{BR}(K_{L}\rightarrow(\mu^{+}\mu^{-})\gamma)=5.13(4)\times 10^{-13}|f\left(z_{TM}\right)|^{2}$, where the dominant error is from $\mathcal{BR}(K_{L}\rightarrow(\mu^{+}\mu^{-})\gamma)$, preventing the measurement of these radiative corrections from this ratio. An improved value of $\mathcal{BR}(K_{L}\rightarrow(\mu^{+}\mu^{-})\gamma)$ or constructing a different ratio, as we do below, can allow sensitivity to these corrections.

The theoretical predictions for $f(z)$ are computed as a series expansion to first order in $z$ with slope $b$. It is typically decomposed into $b=b_{V}+b_{D}$. $b_{V}$ arises from a weak transition from $K_{L}\rightarrow P$ followed by a strong-interaction vector interchange $P\rightarrow V\gamma$ and concluding with the vector meson mixing with the off-shell photon. Here, we denote with $P$ the pseudoscalars ($\pi^{0},\eta,\eta^{\prime})$ and with $V$ the vector mesons ($\rho,\omega,\phi$). The second term, $b_{D}$, arises from the direct weak vertex $K_{L}\rightarrow V\gamma$ which then mixes with $\gamma+\gamma^{*}$ which requires modeling. Following D’Ambrosio and Portoles (1997), the predictions of $b_{V}$ and $b_{D}$ are divided into whether nonet or octet symmetry in the light mesons is assumed.

To compute $b_{V}$, one integrates out the vector mesons from the $P\rightarrow V\gamma$ vertex and assuming a particular pseudoscalar symmetry, the effective Lagrangian is derived and low energy constants can be used. $b_{V}^{octet}=0$ at leading order due to the cancellation between $\pi^{0}$ and $\eta$ in the Gell-Mann-Okubo relation Gell-Mann (1961); Okubo (1962). In the nonet realization, a nonzero contribution coming from $\eta^{\prime}$ yields $b_{V}^{nonet}=r_{V}M^{2}_{K}/M^{2}_{\rho}\sim 0.46$ Ecker (1990), where $r_{V}$ is a model-independent parameter depending on the couplings of each decomposed meson fields in the effective Lagrangian and are ultimately determined by experimental data.

For $b_{D}$, the derivation is more complicated and relies on models. In the naive factorization model (FM) Pich and de Rafael (1991); *Ecker:1992de; *Ecker:1993cq, the dominant contribution to the weak vertex is assumed to be factorized current$\times$current operators which neglect the chiral structure of QCD. A free parameter, $k_{F}$, is introduced that is related to goodness of the factorized current approximation. If this factorization was exact, $k_{F}=1$. In this scheme, $b^{nonet}_{D}=2b^{octet}_{D}=1.41k_{F}$. This model predicts the process $K_{L}\rightarrow\pi^{0}\gamma\gamma$ as well, and we use the unweighted average of the two most recent measurements of this process to fix $k_{F}=0.55(6)$ Lai et al. (2002); Abouzaid et al. (2008).

In the Bergström-Massó-Singer (BMS) model Bergstrom et al. (1983); *Bergstrom:1990uh, the direct transition is instead assumed to be dominated by a weak vector-vector interaction ($K_{L}\rightarrow\gamma+K^{*}\rightarrow\gamma+\rho,\omega,\phi\rightarrow\gamma+\gamma^{*}$). BMS further assumes that no $\Delta I=\frac{1}{2}$ enhancement occurs. This model produces a complete form factor:

The two terms correspond to the vector interchange and direct transition, respectively. Expanding this expression in powers of $z$, we find the BMS model predicts

Under the model assumptions, $-\alpha_{K^{*}}$ is theoretically estimated to be $\sim 0.2-0.3$ Bergstrom et al. (1983); *Bergstrom:1990uh. $C=2.7(4)$ depends on a number of other mesonic decay rates Ohl et al. (1990); *AlaviHarati:2001wd, and we used the modern values Patrignani et al. (2016). The error comes from the experimental uncertainty which is dominated by the two $K^{*}$ measurements. $\mathcal{BR}(K^{*}\rightarrow K^{0}\gamma)$ contributes $\Delta C\sim 13\%$ and $\Gamma_{K^{*},tot}$ contributes $\Delta C\sim 4\%$ due to a disagreement between decay modes. This choice of $C$ and $\alpha_{K^{*}}$ is consistent with the measured rates for $K_{L}\rightarrow\ell^{+}\ell^{-}\gamma$.

D’Ambrosio et. al. advocates the view that $b_{D,BMS}$ is one of a series of contributions to $b_{D}$, which should be summed together with the model-independent $b_{V}$ D’Ambrosio and Portoles (1997). They construct another contribution by factorizing the vector coupling (FMV) similar to FM but first restricting the Lagrangian to left-handed currents. For the different symmetry realizations, $b_{D}^{nonet}=3.14\eta\sim 0.66$ and $b_{D}^{octet}=2.42\eta\sim 0.51$ where $\eta$ is a coefficient multiplying the naive weak coupling $G_{8}$ and like $k_{F}$ is related to the quality of the factorization assumption. We use their value of $\eta=g_{8}^{Wilson}/|g_{8}|_{K\rightarrow\pi\pi,LO}=0.21$. Our theoretical results are compiled in Table 1. These values disagree outside their error, and a 10% precision measurement would be able to discriminate between them. This is in contrast to the radiative Dalitz decays, where the theoretical values are consistent.

The BMS form factor also has been used to phenomenologically fit $K_{L}\rightarrow\ell^{+}\ell^{-}\gamma$ for both $\ell=e,\mu$, and $C\,\alpha_{K^{*}}$ is derived from the differential cross sections of these processes; yielding $(C\alpha_{K^{*}})_{e}=-0.517(30)_{stat}(22)_{sys}$ Abouzaid et al. (2007) and $(C\,\alpha_{K^{*}})_{\mu}=-0.37(7)$ Alavi-Harati et al. (2001), which are each input into our prediction for $(\mu^{+}\mu^{-})$.

We also consider the D’Ambrosio-Isidori-Portolés (DIP) phenomenological $F_{\gamma^{*}\gamma^{*}}(z_{1},z_{2})$ D’Ambrosio et al. (1998):

where $z_{1}=z_{TM},z_{2}=0$ for $(\mu^{+}\mu^{-})$ production. To set $\alpha_{DIP}$, we take the values from $K_{L}\rightarrow e^{+}e^{-}\gamma$, $\alpha_{DIP,e}=-1.729(43)_{stat}(28)_{sys}$ Abouzaid et al. (2007), and from $K_{L}\rightarrow\mu^{+}\mu^{-}\gamma$, $\alpha_{DIP,\mu}=-1.54(10)$ Alavi-Harati et al. (2001). Our phenomenological results are compiled in Table 2. Comparing the phenomenological form factors, they are indistinguishable within uncertainty in $(\mu^{+}\mu^{-})$ production. This is perhaps unsurprising because they arise from the same underlying data, but the difference in functional forms could be discriminated by higher precision data.

Due to the small value of $z_{Ps}\approx 4M_{e}^{2}/M^{2}_{K}$, the branching ratio to positronium, $\mathcal{BR}(K_{L}\rightarrow(e^{+}e^{-})\gamma)=9.31(5)\times 10^{-13}$, is independent of the form factor within the error of $\mathcal{BR}(K_{L}\rightarrow\gamma\gamma)$ and slightly larger than $(\mu^{+}\mu^{-})$. While this branching ratio also has not been measured, one can construct a ratio

which is independent of the $\mathcal{BR}(K_{L}\rightarrow\gamma\gamma)$ uncertainty and directly measures lepton universality without an uncertainty due to $Q^{2}$ binning. By taking the largest and smallest theoretical values of $b$ to give a gross range, we predict $R=0.76(14)$. Applying the same procedure to the phenomenological form factors yields $R=0.707(9)$.

We now focus upon the experimental situation. Throughout, we assume a 10% acceptance. The largest previous experimental data set that could be used to study $\mathcal{BR}(K_{L}\rightarrow(\mu^{+}\mu^{-})\gamma)$ is KTEV. We estimate from the number of events reported for $\mathcal{BR}(K_{L}\rightarrow\ell^{+}\ell^{-}\gamma)$ Alavi-Harati et al. (2001); Abouzaid et al. (2007) that at least 1000 times the luminosity would be required for just one $(\mu^{+}\mu^{-})$ event. From the existing data, one might expect to place a limit on the order of $\mathcal{BR}(K_{L}\rightarrow(\mu^{+}\mu^{-})\gamma)\lesssim 10^{-9}$.

The KOTO experiment at J-PARC has reported $3.560(0.013)\times 10^{7}$ $K_{L}$ per 2$\times 10^{14}$ protons on target (POT) Shiomi et al. (2012). Their 2013 physics run accumulated 1.6$\times 10^{18}$ POT Ahn et al. (2017) which would correspond to 0.015 $(\mu^{+}\mu^{-})$ events. Through their 2015 physics run, 20 times the $K_{L}$ decays have been recorded Ahn et al. (2017), indicating 0.3 produced $(\mu^{+}\mu^{-})$ events and a limit of $\lesssim 10^{-11}$. Unfortunately, the KOTO experiment is designed to detect only photons, and detecting purely photon decay products of $(\mu^{+}\mu^{-})$ would be difficult. The J-PARC kaon beam hopes to run into the 2020s with an additional flux upgrade so a discovery is quite possible in an experiment with lepton identification. The NA62-KLEVER proposal Moulson (2017) for a rare $K_{L}$ beam at CERN hopes to start by 2026 and accumulate 3$\times 10^{13}$ $K_{L}$ over 5 years, which would also be nearly sufficient for single-event sensitivity.

A few channels are available to measure the branching ratio of true muonium: dissociated $\mu^{+}\mu^{-}$ with or without $\gamma$, decayed $e^{+}e^{-}$ with or without $\gamma$, or $\ell^{\pm}\gamma$ similar to SUSY searches with invisible decays Hinchliffe et al. (1997); *Allanach:2000kt. The decay to $\pi^{0}\gamma$ is suppressed by $10^{-5}$ but KOTO can search for it without modification Czarnecki and Karshenboim (2017).

For each channel, different backgrounds matter. The dominant backgrounds will arise from the free decays $K_{L}\rightarrow\ell^{+}\ell^{-}\gamma$. We compute the branching ratio for this by integrating the differential cross section in an invariant mass bin, $M_{bin}$, around the $(\mu^{+}\mu^{-})$ peak to obtain a background estimate. In the case of electrons, the bin is centered around the $(\mu^{+}\mu^{-})$ peak; for muon final states it is defined as $[2m_{\mu},2m_{\mu}+M_{bin}]$. This difference in binning reflects that the muons are above threshold. For bin size similar to KTEV, the values are $\mathcal{BR}(K_{L}\rightarrow e^{+}e^{-}\gamma)_{bin}=1.2\times 10^{-8}M_{bin}$, and $\mathcal{BR}(K_{L}\rightarrow\mu^{+}\mu^{-}\gamma)_{bin}=5.0\times 10^{-9}M_{bin}$ where $M_{bin}$ is in MeV. This large raw background ($\sim 10^{5}\times$ the signal) will have to be reduced, but it has distinct features compared to true muonium decays which can be leveraged.

The smoothness of the background differential cross section around the $(\mu^{+}\mu^{-})$ peak should allow accurate modeling from the sidebands. Reconstruction of the $K_{L}$ allows the energy of the $K_{L}$ to be used to cut on the $\gamma$ and leptonic energies. The two two-body decay topology suggests cuts on momenta and angular distribution would be powerful in background suppression. As an example, for radiative Dalitz decay the angle $\theta_{e}$ between the electrons can be arbitrary, but from the true muonium decay $e$ will have $\theta_{e}\sim m_{TM}/E_{TM}\sim 50^{o}\times\frac{\text{GeV}}{E_{K_{L}}}$. This suggests the higher energy of the proposed CERN beamline would be desirable. Additionally, vertex cuts can be made using the proper lifetime of true muonium $c\tau=0.5n^{3}$ mm, where $n$ is the principal quantum number. A more rigorous study of backgrounds is planned for the future.