Scouring meson decays for true muonium

Rare meson decay experiments promise to measure branching ratios as small as $10^{-13}$. This presents an opportunity to discover the $\mu^+\mu^-$ bound state true muonium. We consider a set of possible channels, all with branching ratios above $\sim10^{-11}$. For the electromagnetic decays $\eta/\eta'\rightarrow(\mu^+\mu^-)\gamma$, theoretical and phenomenological form factors $F_{\eta/\eta'\gamma\gamma^*}(Q^2)$ allow predictions of $\mathcal{BR}(\eta'\rightarrow (\mu^+\mu^-) \gamma)\sim 4.8\times10^{-10}$ and $\mathcal{BR}(\eta'\rightarrow (\mu^+\mu^-) \gamma)\sim 3.7\times10^{-11}$ at the 5% level. Discussion of experimental prospects and potential backgrounds are made.

In this work, we present predictions for BRðη=η 0 → ðμ þ μ − ÞγÞ where true muonium is accompanied by a monochromatic photon (which in the η=η 0 rest frame are 233.2 and 455.6 MeV, respectively) including OðαÞ radiative corrections. These calculations improve upon previous constituent quark model calculations which estimated BRðη → ðμ þ μ − ÞγÞ ≈ 10 −9 [35,37]. Numerous studies of F η=η 0 γγÃ ðQ 2 Þ are available to estimate potential systematics for this process. Other discovery channels involving hadronic final states with BR ≥ 10 −12 are discussed and comments on backgrounds are made.
II. BRðX → Yðl + l − ÞÞ Before proceeding to the specific decays of interest, it is useful to consider where the scaling comes from. The order of magnitude of BRðX → Yðl þ l − ÞÞ can be estimated by multiplying BRðX → YγÞ by α 4 , which arises from jMðγ Ã →l þÃ l −Ã ÞMðl þÃ l −Ã →ðl þ l − ÞÞj 2 ∝ αjψð0Þj 2 . This implies that BRðX → Yðμ þ μ − ÞÞ ⪅ 10 −9 . In Table I, we have included a list of channels which led to branching ratios to true muonium of ⪆10 −12 . There are two processes with BR > 10%: η 0 → ργ and the electromagnetic η → γγ, the latter alongside the percent level η 0 → γγ we will discuss in more detail. The large number of J=ψ events being collected in the near future also presents the possibility of subpercent channels like J=ψ → η 0 γ, as well as the percent level but more complicated inclusive J=ψ → Y h γ [48]. Predictions of Y h þ ðμ þ μ − Þ decays require knowledge of hadronic transition form factors and are more limited in precision compared to processes like η=η 0 → ðμ þ μ − Þγ where the better studied electromagnetic form factors are available. These electromagnetic decays are the main focus of the rest of the paper.
Following previous calculations for the electromagnetic decay of mesons to atoms [35][36][37][38][39][40], the leading-order (LO) decay rate is where ζð3Þ ¼ P n 1=n 3 arising from the sum over all Including the next-toleading-order (NLO) contributions, the squared matrix element is The final jM NLO j 2 is Oðα 2 Þ and therefore can be neglected for our purpose. The remaining OðαÞ radiative corrections can be decomposed into six contributions: vacuum polarizations for the three lepton flavors and hadrons previously computed [35][36][37]40], the leptonic vertex correction [36], and the initial-meson-dependent double-virtual-photon contribution of Fig. 1. With the exception of the doublevirtual-photon contribution, all contributions only depend on jF η=η 0 γγÃ ðzÞj 2 and therefore modify Eq. (1) by a multiplicative factor 1 þ C α π . In contrast, the doublevirtual-photon contribution gives a contribution where N i represent the remaining matrix elements after the form factor has been factored out. The convolution depends upon the initial state and involves the double-virtual-photon form factor. Formally, this prevents writing the radiative corrections into the form of a coefficient C multiplying Eq. (1). Following previous works [50][51][52] in three-and four-body electromagnetic decays of mesons, we approximate the double-virtual-photon form factor as a constant F γγ Ã ð0; z TM Þ. This allows it to be moved out of the integral, and therefore this contribution also provides only a multiplicative correction to the leading-order term. After introducing a counterterm in the MS scheme to cancel out the 1 ϵ divergence of the integral, we are able to find the doublevirtual-photon contribution. This approximation has been compared in three-and four-body pseudoscalar meson decays (e.g., π → e þ e − γ and η → e þ e − μ þ μ − ) where the virtual momentum is not fixed [50] and Eq. (3) is evaluated directly using various models of the double-virtual-photon form factors. From these studies, the approximation introduces systematic errors less than 10% [51,52] where the largest (∼10%) was for η → 4μ, which is very distinct from η=η 0 → ðμ þ μ − Þγ and the smallest (<1%) for the π → e þ e − γ. We take 4% as a conservative estimate for our particular process, given that unlike the threeand four-body processes we do not integrate over F η=η 0 γγ Ã afterward to arrive at the final decay rate due to our twobody nature. Within this approximation the branching ratios (Γ=Γ total ) normalized to the two-photon decay can be expressed as FIG. 1. One of the Feynman diagrams of η=η 0 → γ Ã γ Ã → ðμ þ μ − Þγ which contribute to the branching ratio at Oðα 5 Þ and is proportional to F γ Ã γ Ã ðz 1 ; z 2 Þ.  [49] considered in this work. The first two are electromagnetic decays, while the others are strong decays. Branching ratios to true muonium can be estimated by multiplying by α 4 ≈ 2.8 × 10 −9 .
Channel BR where C η ¼ −0.35ð6Þ and C η 0 ¼ 0.35ð3Þ are the total radiative corrections, which we discuss in detail in the next section.

III. COMPUTING C η=η 0
To reiterate, the six contributions at OðαÞ are vacuum polarizations for the three lepton flavors and hadrons, the vertex correction, and the double-virtual-photon contribution. Together their sum is where the bracketed terms are independent of the initial meson and have been computed in the literature [35][36][37]40]. C iVP indicate vacuum polarization contributions from i ¼ e, μ, and hadrons (where as pointed out in [19] the electron vacuum polarization contributions to true muonium are substantial), C ver is the vertex correction term of [36], while C γ Ã γ Ã is the contribution from the diagram in Fig. 1 and alike. A similar calculation for positronium, where the τ, μ and hadronic vacuum polarization and the double-virtual-photon corrections are negligible, finds the α π coefficient to be C 0 ¼ C VP þ C ver ¼ −52=9 [36]. C iVP are found by computing from the spectral functions ImΠðtÞ. This standard QED vacuum polarization function is known to leading order analytically for the leptons and is derived from the experiment for the hadronic distribution.
To derive C γ Ã γ Ã , we computed the double-virtual-photon matrix element using FEYNCALC [53,54] and obtained its interference with the leading-order matrix element. This reproduces the result of [50] where it is called δ 1γIR;loop ðx; yÞ [Eq. (3.16)]. We then evaluate it for our process at x ¼ z TM and y ¼ 0 to obtain C γ Ã γÃ by subtracting the 1 ϵ pole in the MS scheme.

IV. FORM FACTOR DEPENDENCE
Evaluating Eq. (4) with our computed radiative corrections and the experimentally measured BRðη=η 0 → ðμ þ μ − ÞγÞ [49] BRðη → γγÞ ¼ 39.41ð20Þ%; we arrive at predictions for the branching ratios that depend only on the value of the η=η 0 electromagnetic form factors BRðη → ðμ þ μ − ÞγÞ ¼ 4.14ð3Þ × 10 −10 jfðz TM Þj 2 ; where the dominant error is from BRðη=η 0 → γγÞ. The radiative corrections correspond to 0.07% corrections to the branching ratios, and therefore the statistical uncertainty dominating BRðη=η 0 → γγÞ prevents their measurement via this normalization. This error could be reduced, yielding an ∼0.07% precision level measurement, from a 10 4 times larger η=η 0 sample. This is within reach of near-term experimental proposals. Results for the form factors F η=η 0 γγ Ã ðQ 2 Þ can be broadly classified into three groups: theoretical predictions [55][56][57][58], experimental extractions [59][60][61][62][63][64], and dispersion analyses [65][66][67][68]. The standard parametrization for fðzÞ is a series expansion in z which for all but the dispersion analyses of [67,68] truncate at first order. The theoretical predictions make vastly different assumptions about the coupling between η=η 0 and the photon, as well as different modeling of the mixing between the two mesons. The experimental results rely upon integrating the functional form of Eq. (9) in Q 2 bins, and then nonlinearly fitting b η=η 0 . The dispersive analyses rely upon connecting experimental data for multiple processes to the virtual-photon form factors through analyticity and crossing symmetry. We have tabulated all of the form factors considered in this work in Table II. With the exception of the Brodsky-Lepage interpolation predictions of η [56], the branching ratios predicted for true muonium agree within the uncertainties.

V. LEPTON UNIVERSALITY TESTS
If this discovery channel for true muonium is utilized, one should also be able to measure the positronium (e þ e − ) branching ratio ∼10 −10 . This would allow for the construction of a lepton universality test ratio, where K η ¼ 0.62469ð8Þ and K η 0 ¼ 0.87340ð5Þ. The dominant error in K η=η 0 is from the approximation of the doublevirtual-photon contribution to C η=η 0 and would require an unreasonably precise measurement of the ratio to probe. R could be used in two ways. First, if form factors are taken as inputs, R tests lepton universality at the 5% level due to the double-virtual-photon approximation and Oðα 2 Þ contributions and would be competitive with other lepton universality tests in the light meson sector. The other use would be to constrain the form factors themselves within the uncertainties via R's fixed Q 2 value dependence allowed by leptonic atom production, unlike other decays that require integration over Q 2 . This would partially break the degeneracy between coefficients in Eq. (9) because it does not require integrating the differential cross section over Q 2 bins.

VI. EXPERIMENTAL OUTLOOK
Compared to this theoretical precision, the current and near-future experimental outlooks are less optimistic. While K L experiments are reaching sensitivities of 10 −13 , the sensitivity of J=ψ and η=η 0 searches is worse. At present, the BESIII experiments have the largest J=ψ and η=η 0 datasets. From the 1.3 × 10 9 J=ψ events collected and using the two largest branching ratios, γη=η 0 and ϕη=η 0 , one anticipates 2.0 × 10 6 η events and 7.7 × 10 6 η 0 events with a factor to 10 increase expected in the next decade [70]. This 10 10 J=ψ production capability is right at the edge of what is necessary for an inclusive search. Another similar sized dataset exists from the A2 collaboration in Mainz where η is produced through γp → ηp and have 6.2 × 10 6 η events [71]. The recently approved JLab Eta Factory experiment anticipates collecting 1.3 × 10 8 η and 9.8 × 10 7 η 0 events with 200 days of beam time [72], and beyond that can be run in parallel with the GlueX experiment if the latter is extended beyond 2023 [72]. Further into the future, proposals like REDTOP at FermiLab discussed methods to accrue 10 13 η and 10 11 η 0 events [73].
Clearly the present 10 6 η=η 0 BESIII and A2 data are insufficient for observing true muonium. From these, one would expect to place upper limits on the branching ratios BR ≲ 10 −5 which is 10 4 -10 5 times larger than the predicted rates. This should be compared with the situation for K L → ðμ þ μ − Þγ where upper bounds of BR ≲ 10 −9 are possible at KTEV [74,75], which are superior limits but still 10 3 times as large as predicted [40]. In the next decade, BESIII's larger dataset is still inadequate for even single-event detection through the η=η 0 processes, although the inclusive J=ψ channel is potentially viable. The JLab Eta Factory experiment would be competitive with the possible bounds from KTEV based on the Standard Model predictions. What is required is a proposal like REDTOP at FermiLab which would be sufficient for a discovery of true muonium through the decay of η with N s > 100 events and potentially an observation of the η 0 decay.
The most promising signal channel for discovering true muonium in mesonic decays is e þ e − γ, with a large background from the free decays η=η 0 → e þ e − γ. This background can be computed by integrating the differential cross section in an invariant mass bin M bin , centered around the (μ þ μ − ) peak defined as ½2m μ − M bin =2; 2m μ þ M bin =2   [67,68] BRðη 0 → e þ e − γÞ bin ¼ 3 × 10 −7 M bin where M bin is in MeV. This large raw background, which corresponds to N b ∼ 10 4 × N s , must be reduced or even a REDTOP-like proposal is insufficient. If we estimate the significance of an observation by counting statistics, σ ≈ N s = ffiffiffiffiffiffi N b p , then at least N s > 10 4 is required for an observation. Luckily, several distinct features of the three-body decay compared to true muonium decays can plausibly reduce the background by at least an order of magnitude, as we will explain in what follows.
The two two-body decay topology suggests cuts on momenta and angular distribution would be powerful in background suppression, but exact values of the suppression will be highly detector efficiency dependent. As an example, for radiative Dalitz decay the angle θ e between the e þ e − can be arbitrary, but from the true muonium decay θ e ∼ m TM =E TM ∼ 50°× GeV E η=η 0 . In BESIII, where the typical η=η 0 is produced from the decay of J=Ψ, one finds θ e ∼ 30°. This correlation can be more precisely measured than the invariant mass and can yield a factor of 10 in background suppression. Full reconstruction of the η=η 0 allows for cuts on the γ energy, where the bin resolution is OðM bin Þ. The antiparallel correlation between the γ and the true muonium from the η=η 0 decay yields further background suppression. To estimate this suppression, we compute the three-body decay rate where the angle between e þ and γ is restricted to be greater than 120°(a wider range than the 165°suggested by θ e ). These cuts result in an additional factor of 3 reduction in background. The highly collimated nature of true muonium production in both angle and invariant mass means these cuts will have minimum effect on signal events. With these, we find N b ∼ 300 × N s , sufficient for an observation with a REDTOP-like experiment where hundreds of signal events occur. Further reduction of background is possible with more aggressive or detector-dependent cuts. If the vertex resolution can be made much better than 0.5 mm, cuts can be made using the proper lifetime of true muonium ground state cτ ¼ 0.5 mm which can strongly remove backgrounds but exponentially reduces the signal as well. Otherwise nearly all the e þ e − will be insufficiently separated from the primary η=η 0 → e þ e − γ vertex to distinguish. The advent of picosecond timing in experiments also presents opportunities for clean discrimination from e þ e − backgrounds based on time of flight, but without severe degradation of signal.

VII. CONCLUSIONS
In this paper, we have considered a set of channels for discovery of true muonium, all with branching ratios above ∼10 −11 . For the electromagnetic decays η=η 0 → ðμ þ μ − Þγ, theoretical and phenomenological form factors F η=η 0 γγ Ã ðQ 2 Þ allow predictions of BRðη 0 → ðμ þ μ − ÞγÞ ∼ 4.8 × 10 −10 and BRðη 0 → ðμ þ μ − ÞγÞ ∼ 3.7 × 10 −11 at the 5% level. While present experimental luminosities are inadequate, there is reason for hope in the near future. We have shown how some of the potential backgrounds can be reduced given the special topology of the decay, but future work is necessary to rigorously model the realistic detector backgrounds. Given that rare meson decays with BR ∼ 10 −11 are presently being measured, there seems no fundamental obstacle to discovering true muonium in this way, although it will be a similar challenge.