Ditauonium spectroscopy

We examine the properties of ditauonium, an exotic atom consisting of a pair of opposite-sign $\tau$ leptons bound together by the quantum electrodynamics (QED) interaction in a hydrogen-like state. The energy levels, decay modes and associated partial widths, as well as total widths and lifetimes of the ortho- and para-ditauonium states are calculated. Higher-order QED effects -- including Lamb shifts, hyperfine splitting structure, and partial decay widths corrections -- are incorporated up to approximately next-to-next-to-leading-order (NNLO) accuracy. Beyond the dominant diphoton and difermion decays, the rates of rare decay channels -- including Dalitz, radiative, triple-photon, double-Dalitz, four-fermion, and neutrinos final states -- are determined.

(LO) expressions for all quantities, followed by their more complete and accurate results including higher-order QED corrections. The wavefunctions of the bound states at the origin are the physical quantities of interest to perturbatively compute QED onium spectroscopy properties. The LO wavefunctions are determined by solving the nonrelativistic Schrödinger equation with the Coulomb potential. At LO, the square of the nS wavefunctions at the radial origin (r = 0) depend on the lepton mass and QED coupling, as 1 |ϕ nS (r = 0)| 2 = (αm τ ) 3 8πn 3 .
As we will see below, the zeroth-order energy levels and decay rates are proportional to α 2 m τ and α 5 m τ , respectively. Virtual next-to-leading-order (NLO) and next-to-next-to-leading-order (NNLO) corrections given by one-loop Feynman diagrams corrections to the Coulomb photon, i.e., the nonrelativistic Uehling's potential, have been theoretically calculated for the light leptonium systems. Because of its larger mass, the ditauonium system has accessible many more decay channels that are kinematically forbidden in the positronium and dimuonium cases, i.e., ditauonium has multiple higher-order real corrections that are nonexistent for its lighter siblings. In this work, we compute for the first time the full spectroscopic properties of ditauonium, including real and virtual NLO plus (leading) NNLO corrections to the energy levels (Lamb shifts of order α 3 m τ , and hyperfine splittings proportional to α 4 m τ , respectively), as well as to the decay rates (organising them as functions of zeroth-order widths multiplied by (α/π) and (α/π) 2 terms, respectively), for the para and ortho states separately. The paper is closed with a summary of the main results in Section IV.

II. ENERGY LEVELS
As particle-antiparticle systems, leptonia are intrinsically unstable against annihilation, and this feature makes them markedly different compared to normal atomic systems. The annihilation rates of leptonia are dependent on the overlap of the lepton and antilepton wavefunctions and, therefore, vary for different states. For leptonia in the ground state (where L = 0), the singlet (S = 0) and triplet (S = 1) states can only annihilate into even and odd numbers of photons respectively, because of the Landau-Yang selection rule [15,16] and C symmetry. Thus, para-ditauonium annihilates into two real photons via a t-channel process (Fig. 1, left and center), whereas ortho-ditauonium does it into a pair of charged fermions via an intermediate s-channel virtual photon (Fig. 1, right), or into three real photons via the t-channel process (see Fig. 6 bottom right, later on). As shown below, higher orbital states prefer to radiatively decay to the ditauonium ground state(s) rather than do it through τ + τ − annihilation. The following two subsections present, respectively, the basic LO results and their higher-order QED and relativistic corrections, of the energy levels of the ditauonium spectrum.

A. Leading-order results
At LO, the energy levels of the exotic ditauonium atom can be described by the Schrödinger equation with the Coulomb potential, yielding the same Bohr expression for a hydrogen atom with reduced mass m red = m τ /2, where the last numerical expression is obtained using the values of the tau lepton mass m τ and QED coupling at zero momentum α listed in Table I. The binding energy of the ground state (n = 1) of true tauonium (T ) is, thus, E n=1 = −23.655 keV, and its mass is where the uncertainty is dominated by the current precision of the tau lepton mass [17]. Note that the ∼23.65 keV binding energy of the ditauonium ground states is about ten times smaller than the current uncertainty of the central value of the ditauonium mass itself. Values of the masses of the leptons and (approximate) constituent quarks, tau lifetime, QED coupling, and hadronic photon vacuum polarization self-energy for N f = 3 quark flavours, ∆α (3) had (m 2 T ), and R had (m 2 T ) ratio in e + e − collisions, both evaluated at the T mass scale, used in this work [17]. The quoted value of ∆α (3) had (m 2 T ) is computed using alphaQED19 [18,19]. The Bohr radius of the ditauonium ground-state is a 0 = 2/(αm τ ) = 30.4 fm, and its Rydberg constant amounts to R ∞ = α/(4πa 0 ) = 3.76 keV. Namely, ditauonium is the smallest of all leptonium atoms, and has the largest "photon ionization" energy among them, i.e., it is the most strongly bound of all leptonia. The velocity of each tau in the n-th Bohr orbit is β = 1/(n m τ a 0 ) = α/(2n), which justifies the use of nonrelativistic bound-state perturbation theory (NRQED) [20] to calculate its properties as commonly done for the lighter positronium and dimuonium systems.
FIG. 2: Leading-order energy levels and lifetimes of the three lowest (n = 1, 2, 3) para-(n 1 S 0 ) and ortho-(n 3 S 1 ) ditauonium states decaying into a pair of photons and of lighter charged fermions ( f f = e + e − , µ + µ − , qq), respectively. Figure 2 shows the LO energy levels, determined from Eq. (2), and the LO decay lifetimes determined as explained in Section III, for the three lowest ditauonium states (n = 1, 2, 3). The excited spectrum is obtained considering that a n 2S +1 S 1 ditauonium state can decay via an electric dipole transition, which conserves the spin quantum number, to a n 2S +1 P J state with the emission of a photon with energy E n ∝ α 2 m τ . The radiative transitions from the 3 3 S 1 state to the 2 3 P state, as well as the transition from the latter to the 1 3 S 1 state, have energies of Namely, the Lyman-α photon line of a ditauonium atom transitioning between the first excited (n = 2) and the ground (n = 1) states has an energy of 17.74 keV.

B. Higher-order corrections
Equation (2) predicts that all ditauonium states sharing the same principal quantum number n are degenerate in energy. However, two types of higher-order corrections are known to break such a degeneracy for QED bound states: • The Lamb shift, due to loop quantum corrections of the Coulomb potential, leads to energy shifts among different L states (S, P, D,...) with the same principal quantum number n.
• The hyperfine splittings (hfs) -due to relativistic corrections, the (spin-spin) interaction of the magnetic moments of the components of the bound system, and spin-orbit interactions-lead to separations between singlet and triplet states.

Lamb shift
The Lamb shift is an energy splitting of order O(α 3 ), i.e., it is an NLO correction to the leading-order energy spectrum given by Eq. (2), and can be computed via the expression where ϕ( r) is the ditauonium wavefunction, and V U (r) is the nonrelativistic Uehling's potential characterizing difermion insertions into the Coulomb's photon propagator, which is discussed in detail later around Eq. (23). Numerically, accounting for the three lepton and quark loop contributions, the energy corrections to the first levels read For the lowest 1S ditauonium states, we see that the Lamb effect leads to a ∼115 eV downwards shift of their energy levels. Namely, the ditauonium ground-state mass changes from m T = 3553.6963 ± 0.2400 MeV (adding one extra significant digit to Eq. (3) for visualization of the effect) to m T = 3553.6962 ± 0.2400 MeV, an effect that is one thousand times smaller than the current uncertainty of the ditauonium ground-level mass (driven, as aforementioned, by the current tau lepton mass precision). The uncertainties of the numerical values listed in Eq. (6) from missing NLO corrections to the Lamb shift (two-loop contribution to the Uelhing potential) are of order O(α 4 m τ ), and thus a factor of 10-100 smaller than the values quoted. Such NNLO corrections to the energy levels only impact the accuracy of the energy shifts between states of different L, but not the splitting of different J, S states. In the next section, we consider NNLO corrections to the ditauonium energy spectrum that lead to hyperfine splittings among energy levels with different J, S quantum numbers.
From these results, one can see, first, that the energy splittings of the S-wave family are positive (negative) for the n 3 S (n 1 S) states, whereas the mass hierarchy in the P waves start with n 3 P 2 , followed by n 1 P 1 and n 3 P 1 , with the lightest state being n 3 P 0 . Figure 3 displays the ditauonium mass spectrum including the higher-order effects due to Lamb shifts, Eq. (6), and leading contributions to the hfs structure, Eq. (9), computed here. Calculations exist for the higher-order hfs c i j coefficients of Eq. (7) for (i, j) (0, 0), in the case of positronium and dimuonium [23][24][25]. Such higher-order c i j coefficients are of order unity but suppressed by increasing powers of (α) i ln j (1/α), and therefore they are numerically in the sub-eV range. Given the short lifetimes of the ditauonium systems and the difficulty to produce them in very large quantities for experimental study [13,14], it is unlikely that one will ever be able to carry out such accurate spectroscopic studies of their higher-order hyperfine splitting structure. For this reason, we just retain the first term of Eq. (7), ∆E hfs = α 4 m τ · c 00 , that indicates an ortho-minus para-state energy difference of ∆E 3S−1S hfs = 1.286 + 1.653 ≈ 2.94 eV for the ditauonium ground states.

III. DECAY WIDTHS, BRANCHING FRACTIONS, AND LIFETIMES
Two main differences separate ditauonium from its lightest relatives positronium and dimuonium: the much larger mass and much shorter lifetime of its constituents leptons. First, the significantly more massive ditauonium system has accessible decays to lighter fermion pairs that are kinematically forbidden in the positronium and dimuonium cases. Indeed, ortho-ditauonium can decay into five different pairs of standard model (SM) charged fermions 2 , (τ + τ − ) 1 → f f (with f = e, µ, u, d, s) via the annihilation of both constituent τ leptons into an intermediate single photon (Fig. 1,  right), whereas ortho-positronium can decay to none (except to νν with negligible branching ratio [26,27], if one considers all fermions and not just the charged ones, see Section III C 2), and ortho-dimuonium can only decay to e + e − pairs. Similarly, para-ditauonium has higher-order (τ + τ − ) 0 → γ f f Dalitz decays that are not possible in the positronium and dimuonium cases. This opens up multiple different decay partial widths for the ditauonium states, which increase by a few percent their total width compared to just considering their leading diphoton and difermion decays as we show below. Second, whereas positronium is constituted by a pair of e + e − that are each one of them individually stable, and whereas the µ ± in dimuonium have very large lifetimes (τ = 2.197 µs) for particle physics standards, the τ leptons composing ditauonium are short-lived objects (τ = 290.3 fs) that can individually decay weakly before their bound system itself does. When any of the two tau leptons decays through the weak interaction, the bound state disappears with an effective decay width 3 of As we show below, this leaves only the two lowest ditauonium states as those that can form actual bound states before their inner leptons decay freely. The subsections below present the leading-and higher-order QED corrections for the para-and ortho-ditauonium decays. Electroweak corrections are comparatively suppressed by an O(α 7 m τ ) factor [28], and neglected here. Higher quantum-chromodynamics (QCD) corrections impact decay modes involving quark-antiquark pairs, and they are either phenomenologically incorporated into the computed ortho-ditauonium decays through the experimentally measured R had ratio of hadronic to dimuon cross sections in e + e − collisions, or suppressed by at least a factor of αα s with respect to the leading decay width in the para-ditauonium case, and also neglected hereafter.

A. Leading-order results
As shown in Fig. 1 (left and center), the dominant para-ditauonium decay is into two photons. At zeroth order, the diphoton decay rate of a spin-singlet bound state of principal quantum number n is given by [21,22] 0.018384 eV, for the n = 1 state 0.002298 eV, for the n = 2 state 0.000681 eV, for the n = 3 state .
Ortho-ditauonium decays preferentially into all pairs of SM fermions lighter than half the ditauonium mass, 1 3 S 1 → f f ( f = e, µ, u, d, s), through τ + τ − annihilation into an intermediate single photon (Fig. 1, right) with LO partial widths, where Q f = −1 for charged leptons and Q f = 2 3 , − 1 3 , − 1 3 for u, d, s quarks, respectively, and N c, f is the number of colours of the fermion, i.e., N c, f = 1 for leptons and N c, f = N c = 3 for quarks. Given the smallness of the electron and muon masses compared to the tau one, the product of the two last terms in Eq. (12) is basically unity for all numerical purposes, and the zeroth-order dilepton decay width of ortho-ditauonium reads 0.006128 eV, for the n = 1 state 0.000766 eV, for the n = 2 state 0.000227 eV, for the n = 3 state .
with decays into e + e − and µ + µ − basically happening at the same rate. The values of the inclusive quark-pair decay partial widths can be numerically derived from the same dilepton expression above multiplied by the experimentally However, such modes are suppressed by a O(αv 3 c ) factor (where v c is the relatively low velocity of the bound charm quarks) and a small phase space compared to the leading decay channels, and are therefore ignored here. 3 Numerical conversion of lifetimes and widths in natural units is done via: Γ (eV) = 0.658212/τ (fs). measured value of the ratio of hadronic to dimuon cross sections, R had ≡ σ(e + e − → γ * → qq)/σ(e + e − → γ * → µ + µ − ) ≈ q=u,d,s N c Q 2 q evaluated at the ditauonium mass scale. Namely, 0.013482 eV, for the n = 1 state 0.001685 eV, for the n = 2 state 0.000499 eV, for the n = 3 state , .
where the last approximate numerical equality uses the current experimental value of the R had ratio ( Table I).
The total LO width of ortho-ditauonium -adding dielectron, dimuon, and diquark decays-amounts thus to: 0.025738 eV, for the n = 1 state 0.003217 eV, for the n = 2 state 0.000953 eV, for the n = 3 state. .
For excited states with L > 0, the overlap of the wavefunctions is zero resulting in preferential radiative decays to lower 1S states rather than annihilation decays. The radiative decay widths of the excited (3S and 2P) states can be derived from the generic expressions for the ratios of decay widths of true leptonium states, The partial decays widths from higher to lower orbitals are thus given by The annihilation partial decay widths of P-wave states are at least of order O(α 7 ) and thus comparatively negligible. The zeroth-order results for the energy levels and leading decays of ditauonium states are shown in Fig. 2, where the lifetime has been derived from each width via τ = 1/Γ. From this plot, we see first that only the two lowest-energy levels, para-ditauonium 1 1 S 0 (with τ = 35.8 fs) and ortho-ditauonium 1 3 S 1 (with τ = 25.6 fs) have lifetimes smaller than half that of its constituent tau leptons (τ = 145.15 fs), and can thereby form a transient bound state before the weak decay of any of them. Second, one can observe that the lifetimes of para-and ortho-states are relatively similar for ditauonium, at variance with the positronium case where the ortho-state can only kinematically decay into the 3-photon mode (which is, at least, α times suppressed compared to the para-positronium 2-photon decay) and thereby features three orders-of-magnitude longer lifetime than the para-state.
B. Para-ditauonium decay 1. Virtual NLO and (partial) NNLO corrections Figure 4 shows the NLO virtual corrections to the diphoton para-ditauonium partial width, Γ(n 1 S 0 → γγ). They include virtual photon exchanges between the tau leptons, as well as modifications of the ditauonium wavefunction at the origin, |ϕ nS (0)|, due to vacuum polarization effects affecting the common Coulomb potential of the two tau leptons and that are effectively of order O(α/π). These latter effects are indicated in the rightmost panel of Fig. 4 by the fermion loop insertion to the Coulomb potential (dashed line) typical of NRQED calculations [29]. Calculations for all these higher-order corrections exist for the dimuonium case [7,28,[30][31][32][33], and are extended and applied to ditauonium in this section. The radiative virtual QED corrections are given by the set of four leftmost diagrams displayed in Fig. 4. These virtual NLO corrections were exactly calculated in [34,35] for para-positronium. In the nonrelativistic limit, the overall sum of such diagrams amounts to a relative width correction of size The universality of the correction in this limit justifies its use also for the para-dimuonium case [30,32], and hence also for para-ditauonium in the present study. The loop in the Coulomb's photon propagator of Fig. 4 (right) represents modifications of order O(α/π) of the wavefunction at the origin. For para-ditauonium, the corrections incorporate all fermion loops, though they are dominated by the contribution of the lightest electron loop. The first approach to obtain this correction, according to the expression was performed in [7] and further improved in [30,32]. This NLO coefficient for the nS para-ditauonium states can be computed via the integral where, first, g nS ( ) are the reduced Green's functions of the dimensionless radial variable = r/a 0 for the nonrelativistic Coulomb potential, given by the following expressions for the 1S and 2S states, respectively, with γ E ≈ 0.577216 the Euler's constant, and where the integration kernels read Secondly, the V U (r) function in Eq. (19) is the nonrelativistic Uehling's potential, characterizing polarization insertions into the Coulomb's photon propagator [30,32], given by where v f is the integration variable defined as v 2 f = 1 − 4m 2 f /s, with s the square of the invariant mass of the fermion pair in the loop, and where the effective Coulomb's photon mass λ f (with dimension L −1 ) takes the values, upon insertion of a fermion f loop. Such a potential is known fully analytically expressed in terms of modified Bessel functions and Bickley functions [36]. It is worth mentioning that the integral's kernel term proportional to does not depend on m f and so do not either the integration bounds. Hence, the mass dependence on the type of loop mainly arises from λ f in the argument of the exponential term. After cross-checking that with Eq. (19) we can reproduce the same values of the NLO coefficients derived for the para-dimuonium case in refs. [30,32], we determine the numerical coefficients for S-wave ditauonium by inserting in the loop all fermions f = {e, µ, τ, u, d, s} with the masses listed in Table I. These Coulomb corrections are clearly dominated by the electron loop (accounting for about 90% of the total value) and, therefore, the relatively large uncertainties of the constituent masses of the inserted light quarks are not numerically important.
In addition, virtual NNLO corrections related to the low energy logarithmic contribution arising from spin-spin contact interactions and pure orbit corrections in the Breit Hamiltonian [22], which were first derived for parapositronium [37,38], are found to amount to: for the n = 1 (valid also for n = 2) para-ditauonium case.

Real NLO and (partial) NNLO corrections
The higher-order real corrections to para-ditauonium decays include, first, the splitting of one of the decay photons into any kinematically allowed fermion-antifermion pair (dielectron, dimuon, and light diquarks). The corresponding diagrams are shown in the left and center panels of Fig. 5. Beyond this, the rightmost panel of Fig. 5 shows the NNLO double photon splitting into four fermions. The Dalitz γ f f decay width can be calculated at the lowest order as Such a perturbative formula applies well for the f = e, µ leptons, for which we find Γ NLO Dalitz (n 1 S 0 → γµ + µ − ) = 6.55457 For the quarks, we need in addition to sum over the three quark flavour charges as follows, In the equation above, we have replaced the purely perturbative expression with quark masses of the first line with the second-line expression that employs the hadronic quantities ∆α (3) had (m 2 T ) and R had (m 2 τ ) listed in Table I, where the hadronic running of the QED coupling, ∆α (3) had (m 2 T ), can be obtained from the experimental R had (s) ratio via dispersion relations.
Finally, we compute the NNLO double-real correction given by the "double Dalitz" diagram shown in Fig. 5 (right). The partial widths can be cast into the following generic form: with the values of C NNLO 2Dalitz ( f f f f ) calculated numerically with the phase space integrator of HELAC-Onia [39,40], and amounting to C NNLO 2Dalitz (e + e − e + e − ) = 23.40, C NNLO 2Dalitz (e + e − µ + µ − ) = 13.87, C NNLO 2Dalitz (e + e − qq) = 12.07, In the above equation, the notation "qqq q " indicates the sum of all four quark final states, including both same flavour and different flavour channels. It is interesting to notice that in the asymptotic heavy-m τ limit, the coefficients for the same-flavour 4-lepton final states approach the results of the expression with the double logarithm term matching that given in Refs. [32,37,38] for positronium and dimuonium. Similarly, the asymptotic formulas read, for the same-flavour 4-quark final states, and, for the generic different-flavour double-Dalitz case.

Combined higher-order corrections
The total annihilation decay width of para-ditauonium (we focus just on the n = 1 state hereafter) can be written in a compact form, by combining all previous results, as follows The LO total annihilation decay width of para-ditauonium has a single one-channel contribution, i.e., which is given by Eq. (11). The P NLO 1S,tot parameter accounts for the numerical NLO coefficients derived in Eqs. (17) where one can see that the largest higher-order contribution is from the γe + e − Dalitz decay, which accounts for about half of all NLO corrections. P NNLO 1S,tot is a numerical coefficient accounting for the partial virtual and real NNLO corrections given, respectively, by Eqs. (26), (31) and (32), The corresponding values of all individual coefficients are listed in Table II. To visualize the relative size of the real and virtual higher-order corrections, and where Γ NLO Dalitz and Γ NNLO 2Dalitz sum up, respectively, all individual Dalitz and double-Dalitz decays. From these results, a few quantitative facts can be highlighted: (i) altogether, the higher-order corrections augment the total 1 1 S 0 annihilation decay width by +5.0%, (ii) the virtual NLO + NNLO corrections increase the dominant diphoton decay by +0.8%, and (iii) the single-and double-Dalitz decays occur with a combined B = 3.25% rate (over the total width including single-τ weak decays, given by Eq. (41) below). Although the impact of the higher-order effects computed here may seem relatively small, the presence of final states with charged particles can facilitate the first experimental detection of para-ditauonium as discussed in Refs. [13,14].  Table IV summarizes all the properties of the 1 1 S 0 state including all individual partial widths, up to the highest NLO + NNLO accuracy computed here, and associated B X = Γ X /Γ tot branching fractions. The total ditauonium width is determined by adding all individual partial widths plus the effective width due to the weak tau decay, i.e., The last significant figures of all values listed in Table IV have been rounded off to approximately match the associated theoretical accuracy of each width. Theoretical uncertainties due to missing higher-order corrections are very small for the total width and main diphoton decay channel, as they have been computed here including up to the most important NNLO corrections. The uncertainty is therefore at the NNLO level, (α/π) 2 ≈ 10 −4 accounting for O(10) coefficient prefactors (this is a realistic order-of-magnitude estimate, although partial cancellations between the complete set of NNLO virtual and real corrections, which have not been fully computed here, are not excluded). For the Dalitz γ f f para-ditauonium decays, since they are LO for this mode, their relative uncertainty is O(α/π) ≈ 10 −2 . The propagated parametric uncertainty due to the tau mass precision is around 7 · 10 −5 for all quantities, which linearly depend on m τ through the LO widths. The uncertainty of the tau decay width due to the tau lifetime, Eq. (10), is around 2 · 10 −3 , which propagates into 4 · 10 −4 relative uncertainty of the para-ditauonium total width. As one can see, theoretical uncertainties (intrinsic and parametric) are very small and very likely well beyond the reach of any potential experimental precision. Let us deal first with the NLO virtual corrections to the 1 3 S 1 wavefunction at the origin (Fig. 6, top right), which can be calculated through fermion loop insertions in the Coulomb's photon propagator, as done for the para-ditauonium case in Eq. (18), Here, the C NLO,1S Coul = 5.805 coefficient has the same numerical value as computed before, Eq. (25). We address next the ortho-ditauonium radiative corrections. The two γ-emission diagrams in Fig. 6 (bottom left and center) produce divergent infrared double logarithms that need to be canceled out against similar terms produced by the virtual correction diagrams shown in the three first top panels of the figure. The real and virtual radiative corrections can be calculated with the standard techniques by defining x f = m 2 f /m 2 τ , and ignoring power corrections of the form O(αx f ), as follows: The first terms of this equation account for inclusive photon exchanges/emissions and leptonic vacuum polarization loops (the symbol indicates the real part of the expression in parenthesis), whereas the last term of the inner sum stands for the hadronic vacuum polarization contributions of the virtual photon at the ditauonium mass (i.e., for the N f = 3 hadronic loop contributions of the third top diagram of Fig. 6) quantified by the ∆α (3) had (m 2 T ) = 0.0077 term (Table I).
Equation (43) determines the NLO contribution to the ortho-ditauonium width from photon emissions and exchanges inclusively. Here, next, we evaluate the size of the contributions from 1 3 S 1 decays with photons explicitly tagged or measured in the final state, shown in the bottom (left and center) diagrams of Fig. 6, where a photon energy cutoff needs to be introduced to avoid infrared divergences. We restrict ourselves to the case of the dilepton decays 4 with one final-state photon, 1 3 S 1 → + − γ. Collinear divergences can be regularised by keeping the outgoing lepton masses nonzero, while a threshold on the photon energy in the ditauonium rest frame is required to remove soft divergences, which we take as E γ > E thr γ ≈ O(m ). By introducing two variables, x = m 2 /m 2 τ and = E thr γ /m τ , the partial width then reads, with the collinear ln x and soft ln logarithms appearing, and where we have ignored power corrections of order O( , x ). We have also computed the full expression with the and x dependencies, which is a bit lengthy and we refrain ourselves from writing it here. The numerical values for = m /m τ are 23.76 for = µ.
To estimate the size of the purely virtual NLO corrections of ortho-ditauonium, we can combine Eqs. (43) and (45) and determine the NLO partial width for 1 3 S 1 → + − without real photon emission (with an energy above m ), as follows, The leading NNLO virtual contribution from the wavefunction correction from the Breit Hamiltonian [22] is This (small) NNLO virtual correction is negative and will slightly decrease the partial Γ f f (γ) decay widths although, as we see next, there are extra real NLO + NNLO contributions (3-photons and 4-fermions final states) that contribute positively to the total ortho-ditauonium decay rate. Figure 6 (bottom right) shows the 3-photon decay channel of ortho-ditauonium, which is suppressed by an extra α factor compared to the diphoton para-ditauonium decay given by Eq. (11). At its lowest order, this partial decay width is The (infrared-finite) NLO virtual contribution of this channel can be calculated by the standard perturbative techniques, and amounts to ∆Γ NLO virt.exch (n 3 S 1 → γγγ) = −13.44 The NLO correction from the wavefunction at the origin for the 3-γ decay is Coul,3γ α π 2 Γ (0) (n 3 S 1 ) , with C NLO,nS Coul,3γ = C NLO,nS Coul · C 3γ = 1.6026 for n = 1.
One can see that the net effect of the NLO corrections to the 3-photon decay, combining Eqs. (49) and (50), leads to a small decrease (by about −1.8%) of this rate. Finally, the ortho-ditauonium can also have 4-fermion decay channels contributing at real NNLO accuracy to the total width (Fig. 7). The partial widths can be again cast into the following generic form with the numerical values of the C NNLO 4 f ( f f f f ) coefficients, derived with the help of HELAC-Onia [39,40], given by FIG. 7: Representative Feynman diagrams of (real NNLO) ortho-ditauonium decays into four fermions.

Rare weak decays
Of the two spin states of ditauonium, only 5 the ortho state can decay weakly into a pair of neutrinos. The two amplitudes that contribute to the decay T 1 → ν ν , are W exchange in the t channel and Z annihilation in the schannel (Fig. 8). The Z diagram is a tiny correction to the virtual photon annihilation (Fig. 1, right), and there is a destructive interference between the W and Z exchange amplitudes. The corresponding partial widths were found negligible for positronium and dimuonium -for positronium they are O(10 −18 ) and O(10 −21 ) for decays into like-and unlike-flavour between the e + e − and neutrinos respectively [27]-, but they are comparatively enhanced by powers of the tau to electron (or muon) masses for ditauonium, and it is worth to estimate their importance here. For ortho-ditauonium, the neutrino-pair decay widths read where m W,Z are the W and Z boson masses, and s w and c w are the sine and cosine of the Weinberg angle. Both channels turn out to be also very rare for ditauonium due to the strong suppression driven by the m 4 τ /m 4 W,Z factor.

Combined higher-order corrections
The total annihilation decay width for n = 1 ortho-ditauonium can be now obtained from all previous results and written in the compact form where Γ (0) (1 3 S 1 ) is given in Eq. (15). , respectively-, whereas the (small) NNLO corrections are larger for the para-than for the ortho-state as one can see comparing Eqs. (39) and (55). This implies that the total NLO + NNLO corrections increase by almost the same amount, about 5.0%, the LO annihilation rates of both states.  T state  Table VI displays the partial annihilation decay widths of ortho-ditauonium grouped at LO, NLO, and NNLO accuracies to visualize the relative size of the real and virtual higher-order corrections, where we define Table VII lists all the properties of the 1 3 S 1 state computed here. The total ortho-ditauonium width listed is determined by adding all individual partial widths plus the effective width due to the constituent tau weak decays, i.e., Our main findings about ortho-ditauonium decay rates can be summarized as follows: (i) the radiative real and virtual NLO corrections increase the difermion decays by +5.0%, (ii) the virtual NNLO corrections are tiny and negative but "compensated" in the total width by positive real NNLO corrections from new 4-fermion channels that open up at this level of accuracy (with combined branching fractions of B ≈ 0.04%), and (iii) the decays into a pair of neutrinos have O(10 −7 -10 −9 ) rates. The ortho-ditauonium branching fractions are thus dominated by decays into a pair of light diquarks (with or without γ emission), B qq(γ) = 44.82%, with the actual hadronic final states mostly consisting of a few charged and neutral pions and/or, to a less extent, kaons. The combined dilepton final states with or without photon emission, 1 3 S 1 → e + e − (γ), µ + µ − (γ), have a branching fraction of B + − (γ) = 2 × 20.37% = 40.74%. The presence in the ortho-ditauonium decays of final states with different charged particles can facilitate the measurement of this exotic atom, given that the experimental momentum and vertex resolutions are better for them than for photons [13,14].
As discussed for the para-ditauonium case, theoretical uncertainties due to missing higher-order corrections for the total width and main difermion decay channels are very small, since we have computed them including up to the most important NNLO corrections. The uncertainty is therefore at the NNLO level, α 2 ≈ 10 −4 (this is a realistic orderof-magnitude estimate, although partial cancellations between the complete set of NNLO virtual and real corrections, which have not been fully computed here, are not excluded). For the decay to hadrons, the perturbative uncertainty is also NNLO (i.e., 10 −4 ), but the uncertainty of ∆α (3) had (m 2 T ) and R had (m 2 τ ) can propagate into a larger/dominant value. Since we have provided the analytic result, Eq. (14), the numbers can be easily updated whenever the experimental value of R had (m 2 τ ) is refined. In all numerical evaluations, which linearly depend on m τ via the LO widths, the propagated parametric uncertainty due to the tau mass precision is around 7 · 10 −5 . The uncertainty in the tau decay width propagates into a 3 · 10 −4 relative uncertainty of the ortho-ditauonium total width. Therefore, theoretical uncertainties (both of intrinsic and parametric nature) are very small, and very likely well beyond the precision of any actual experimental measurement.

IV. SUMMARY
We have presented the first study of the spectroscopic structure of the purely leptonic system consisting of two τ leptons, bound by their mutual QED interaction, known as ditauonium. First, the basic zeroth-order expressions for its energy levels and dominant decay widths have been presented. The ground state (1S) has a binding energy of −23.655 keV, leading to a ditauonium mass of m T = 3553.696 ± 0.240 MeV (where the uncertainty is dominated by the current tau lepton mass precision). Ditauonium decays dominantly through annihilations into pairs of photons and of lighter charged fermions for para-(1 1 S 0 ) and ortho-(1 3 S 1 ) states, respectively. Secondly, QED corrections at NLO and (partially) NNLO accuracy have been calculated for the energy levels and for the rates of all decay modes kinematically accessible at each level of accuracy. For the ground state, we find that the Lamb shift decreases its binding energy by about 0.115 keV, whereas the hyperfine splitting separates the para-and ortho-states by O(3 eV).
A detailed study of all partial decay widths of para-and ortho-ditauonium has been carried out. Including all NLO and the most important NNLO corrections, the annihilation decay widths of 1 1 S 0 and 3 1 S 1 ditauonium states increase TABLE VII: Main properties (mass m X , J PC quantum numbers, total width Γ tot , lifetime, as well as partial decay widths Γ X and associated B X branching fractions) of the lowest-energy ortho-ditauonium bound state computed in this work. The two inclusive decay modes + − (γ) have been also broken-down into the exclusive modes + − and + − γ (where the latter require threshold photon energies E thr γ > m e,µ for the dielectron and dimuon channels, respectively, see text for details).
Ditauonium is the heaviest and most compact purely leptonic "atomic" system and remains experimentally unobserved to date. The results presented here are of usefulness to carry out potential experimental measurements of its production and study its properties at current and future colliders via multiple different final states whose decay rates have been quantified here for the first time. Concrete feasibility cases, beyond those described in [13], will be presented in an upcoming work [14]. Ditauonium studies can provide novel tests of bound-state QED that are sensitive to physics beyond the standard model at higher energies than those of its lighter siblings positronium and dimuonium.