{'fraction_non_alphanumeric': 0.03967692243640451, 'fraction_numerical': 0.014155460260627004, 'mean_word_length': 4.429475587703436, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
alpha^2 corrections to parapositronium decay: a detailed description

We present details of our recent calculation of alpha^2 corrections to the parapositronium decay into two photons. These corrections are rather small and our final result for the parapositronium lifetime agrees well with the most recent measurement. Implications for orthopositronium decays are briefly discussed.


I. INTRODUCTION
Precision measurements with nonrelativistic bound states provide a sensitive test of bound state theory based on Quantum Electrodynamics (QED). Theoretical predictions generally agree well with experimental results. However, there are a few observables where the agreement is not satisfactory. The best known example is the lifetime of the orthopositronium (o-Ps), where some experimental results [1,2] differ from the theory by several standard deviations. The theoretical prediction is not yet complete, because two-loop QED corrections to the o-Ps decay have not been evaluated yet. However, those corrections would have to be very large to reconcile the theory with experiment.
Recently, we have reported a calculation of corrections to the parapositronium (p-Ps) decay into two photons [3] in the second order in the fine structure constant α. We argued there that the smallness of the O(α 2 ) corrections to the p-Ps lifetime makes it unlikely that the analogous effects in o-Ps could alone explain the discrepancy between theory and experiment.
The purpose of the present paper is to present a detailed description of the calculation reported in [3]. We begin with a brief discussion of the current theoretical and experimental status of positronium (Ps) decays.
Theoretical predictions for p-Ps and o-Ps decay rates into 2 and 3 photons, respectively, can be expressed as series in α: where are the lowest order decay widths of the p-Ps and o-Ps, respectively, and the ellipses denote unknown higher order terms which we will neglect in our analysis. Corrections of O(α) were calculated in [4] for p-Ps. For o-Ps the most accurate result was obtained in [5], where references to earlier works can be found. The logarithmic two-loop correction was found in [6] for o-Ps and in [7] for p-Ps. The leading logarithmic correction at three loops was computed in [8]. Some partial results on the O(α 2 ) corrections for both p-Ps and o-Ps can be found in [5,[9][10][11][12]. A complete calculation of B p has been reported recently in [3] but the value of B o has not been obtained so far. Eqs. (1,2,3) give the following predictions for the lifetimes: On the experimental side, the Ann Arbor group [1,2] which, for B o = 0, differ from (5) by 9.4σ and 6.3σ respectively. This apparent disagreement of experiment with theory has been known as the "orthopositronium lifetime puzzle" for a long time. More recently, an independent measurement by the Tokyo group found [13] Γ exp o−Ps (SiO 2 measurement) = 7.0398(29) µs −1 , which agrees with the theory if B o is not too large. 1 Thus, the present experimental situation is rather unclear. If future experimental efforts confirm the Ann Arbor results (6), the orthopositronium lifetime puzzle could be solved if B o turns out to be unusually large, e.g. ∼ 250 for the vacuum measurement. Alternatively, one might speculate that some "New Physics" effects such as o-Ps decays involving axions, millicharged particles, etc., cause the excess of the measured decay rate over the QED predictions. Some of those exotic scenarios seem to have already been excluded by dedicated experimental studies. (For a review and references to original papers see e.g. [15].) p-Ps can decay into 2 photons and its lifetime is too short to be measured directly. For a long time its precise value remained unknown. However, it was realized in [16] that such a measurement would be very useful, since the calculation of the coefficient B p is much easier than that of B o . Therefore, p-Ps offers an easier precision test of the bound-state QED. This observation motivated the measurement of the p-Ps lifetime [16] which is 6.5 times more accurate than the best previous results.
To fully utilize that experimental result and enable the rigorous test of bound-state QED envisioned in [16], we undertook a complete calculation of the O(α 2 ) corrections to p-Ps rate [3]. We found that the non-logarithmic part of those corrections is small; the coefficient B p is and the theoretical prediction for the p-Ps lifetime becomes (half of the logarithmic α 3 term in (1) is taken as an error estimate) Γ theory p−Ps = 7989.64(2) µs −1 .
Comparing this number with the most recent experimental result [16] Γ exp p−Ps = 7990.9(1.7) µs −1 (10) we find excellent agreement between theory and experiment.
In the next Section we explain our approach to Ps decays with the example of the leading order calculation. In Sections III and IV we discuss, respectively, the so-called soft and hard O (α 2 ) contributions. The final result is presented in Section V. 1 Another measurement of the o-Ps lifetime was performed earlier in Mainz [14] with the result Γ exp o−Ps = 7.031(7) µs −1 . This result is consistent with the theoretical predictions but uncertainties are clearly larger. For a summary of older experimental studies see [15].

II. FRAMEWORK OF THE CALCULATION
Positronium is a bound state of an electron and a positron. Its energy levels and lifetimes can be well understood within the framework of nonrelativistic expansion in QED. The decay width of p-Ps into two photons can be written as: where P = (M p−Ps , 0) is the four-momentum of the positronium, k 1 and k 2 are photons' momenta, φ(p) is the positronium nonrelativistic wave function, A(λ, p) is the amplitude of the process e + e − → 2γ, λ denotes polarizations of the photons, and is the spinor part of the p-Ps wave function. The wave function of the p-Ps is normalized to unity, hence the usual normalization factor, 1/2M p−Ps , is absent in Eq. (11). From Eq. (11) one sees that Γ p−Ps is determined by two quantities: the annihilation amplitude of the electronpositron pair into a pair of photons, and the positronium bound state wave function. Since the typical velocity of electron and positron in Ps is O (α), the annihilation amplitude can be expanded both in α and in the relative momentum p of the electron and positron. Corrections to the positronium wave function are computed using the standard time-independent perturbation theory and Breit correction to the nonrelativistic Hamiltonian. We begin with the calculation of the decay rate to leading order. Since the electron and positron in Ps are nonrelativistic, we employ non-covariant perturbation theory. Although it is not the only way, the non-covariant technique permits an easier evaluation of the soft corrections to the p-Ps decay rate.
The on-shell amplitude of the process e + e − → γγ reads Here E p = √ m 2 + p 2 is the energy of the electron (or positron), Λ ± (p) are the projectors on the mass shell, and k is the three-momentum of the photon in the final state.
To the leading order one can neglect the dependence of the annihilation amplitude on the momentum |p| ∼ mα, which is small in comparison with m and |k| ∼ m. One obtains the following leading order amplitude: From Eqs. (16,12) it follows that the photon polarization vectors are perpendicular to each other, i.e. e 1 e 2 = 0. Computing the trace in Eq. (11) and integrating over p one obtains which gives the following leading order decay width: The calculation of higher order corrections to positronium lifetime is performed in the framework of the Nonrelativistic Quantum Electrodynamics (NRQED) [6] which we regularize dimensionally [17]. In this framework, all contributions are divided into soft and hard corrections. The soft contributions come from the momenta region of the order of k ∼ mα and thus are sensitive to the details of the bound state dynamics. On the other hand, the hard corrections arise as contributions of the relativistic momenta k ∼ m; their effect can be described by adding δ(r)-like terms to the nonrelativistic Hamiltonian.
Special care in the present calculation is required because of a single γ 5 matrix in the positronium wave function. Since a consistent treatment of single γ 5 is known to be a problem in dimensional regularization, below we describe how we have dealt with it.
We use the following fact: if the p-Ps decays into two photons with polarizations e 1 and e 2 , then the photons are polarized so that e 1 e 2 = 0. Let the first photon be polarized along the x and the second photon along the y axis; the three-momentum of the first photon k being along the z-axis. If we introduce two auxiliary vectors then the standard four-dimensional representation of the γ 5 matrix, γ 5 = iγ 0 γ x γ y γ z , can be written as: We use this equation to define γ 5 in d dimensions.
With an additional trick it is possible to avoid any reference to photon polarization vectors. Let us consider Eq. (11) and use there the explicit representation for the γ 5 matrix, Eq. (19). Then the following trace has to be computed: Taking into account that e 1 e 2 = e 1,2 k = e 1,2 P = 0, one concludes that the final result can only depend on e 2 1 and e 2 2 . In this situation, one can average Eq. (20) over directions of e 1 and e 2 , provided one respects the above constraints. We arrive at the following formula: where According to the discussion of different momentum regions contributing to the second order corrections, we divide the coefficient B p in Eq. (4) into three parts: where B squared p is the contribution of the one-loop amplitude squared and B hard, soft p are the hard and soft two-loop contributions. The square of the one-loop amplitude is easily obtained from the one-loop result: In contrast to the first order correction which arises from hard photon exchange only, the second order correction is more difficult to compute, because of the appearance of the soft scale effects. In addition, hard and soft corrections are not finite separately. Below we discuss the calculation of these corrections.

III. SOFT SCALE CONTRIBUTIONS
As follows from Eq. (11), in order to obtain B soft p one has to compute relativistic corrections to the annihilation amplitude e + e − → γγ and relativistic corrections to the positronium wave function induced by the Breit Hamiltonian. Accordingly, the soft contribution is separated into two pieces: A. Relativistic corrections to the amplitude The calculation of relativistic corrections to the amplitude is straightforward. One starts with the on-shell amplitude Eq. (14) and expands it up to relative order O(p 2 /m 2 ). The calculation of these corrections can be performed in three dimensions. To demonstrate this, let us write the correction to the amplitude in the form: 2 To calculate the correction to the p-Ps lifetime induced by δA, we have to compute the following integral: To this end we use the Schrödinger equation in the momentum space: Rewriting shifting the integration momenta in the first term p → p + k and using the fact that the scaleless integrals vanish in dimensional regularization, we arrive at We see that the amplitude A (2) ij is needed only in the ǫ → 0 limit, where it can be easily calculated. We obtain: which induces the following O(mα 2 ) correction to the p-Ps lifetime: Recently there has been some discussion in the literature [11,10,18] concerning the linearly divergent integral in Eq. (29). Our approach to the linear divergence is based on dimensional regularization which permits a consistent treatment of hard and soft corrections simultaneously. We dealt similarly with linearly divergent integrals in our recent calculations of the O(mα 6 ) corrections to positronium S-wave energy spectrum [19,20] and found agreement with earlier results obtained in a different regularization scheme [21]. This gives us confidence in the result given in Eq. (29).

B. Relativistic corrections to the wave function
Relativistic corrections to the positronium wave function can be computed using Breit Hamiltonian. Since we regularize all divergences dimensionally, we need that Hamiltonian in d dimensions, as it has been derived in [20].
Breit Hamiltonian projected on the S-states reads Here the Pauli matrices σ i and σ ′ i act on the two-component spinors of the nonrelativistic electron and positron respectively, and is the d-dimensional generalization of the Coulomb potential. The operator U(r, p) induces the following correction to the p-Ps width, where the first term is due to the δ(r)-part of the operator U and the second term is due to the remaining terms of that operator. We note also, that Γ LO in the above formula stands for the leading order decay width computed in d-dimensions, in contrast to the threedimensional result Eq. (18). All necessary formulas can be extracted from the calculation described after Eq. (36) of Ref. [20]. One obtains: Thus the wave function correction contribution to the positronium lifetime becomes On the LHS of the above equation we have separated the logarithm of the fine structure constant to be consistent with division of corrections introduced in Eq. (1).

C. Final result for the soft contributions
The sum of the corrections to the annihilation amplitude (31) and to the wave function (36) gives the final result for the soft contributions, Vacuum polarization insertions into the one-loop graphs (an example is shown in Fig. 2(VP)) were computed in [22,23], (39)

A. Light-by-light scattering contributions
One class of the second order corrections to the p-Ps lifetime arises from the photonphoton scattering, shown in Fig. 3. These contributions are relatively small; it is interesting, however, that the positronium lifetime measurement of increased precision may become sensitive to effects of non-linear QED.
For the diagrams shown in Fig. 3, we find that the planar ones are equal (a=b), and so are the non-planar ones (c=d). Further, both classes remain unchanged when we cross the internal photons. Therefore we only need to compute two diagrams, one planar and one non-planar. We will call them L 1 and L 2 . The contribution of those diagrams to the final expression for the amplitude is 4L 1 + 2L 2 . The factor 4 arises because we have two types (a) and (b) which differ by the orientation of the fermion loop, and both can have parallel or crossed internal photon lines. For the type L 2 we only have a factor 2 for the orientation of the fermion loop, because crossing of the internal photon lines simply exchanges (c) and (d).
The light-by-light diagrams in Fig. 3 have an imaginary part because of the two-photon cut. This complicates a numerical integration over Feynman parameters and we used a different method to evaluate them. The idea is to (formally) assign a large mass M to the internal fermion line. The diagram may then be expanded in ratio m/M using the so-called Large Mass Expansion [24][25][26]. This reduces the task to the calculation of two-loop vacuum or massless propagator integrals. The price to be paid is that the result is a series in m/M, while we are interested in the value of the series at m/M = 1.
Fortunately, using symbolic manipulation programs one can compute many terms of the series; in our calculation about twenty terms were computed for each diagram. The resulting series converge well, especially for L 2 , and the number of computed terms is sufficient to obtain an accurate estimate of this contribution at x = 1. The behavior of the series for L 1 is improved if we make a change of variables m/M = z/(2 − z). The n−th term of the series in variable z decreases roughly like ln 2 n/n 2 . Finally, we find L 1 = 0.494 (30), and the contribution of these diagrams to the coefficient B p is: B hard p (LL) = 4L 1 + 2L 2 = 1.28 (13). (41)

B. Two-photon corrections
Another class of corrections is generated by the two-photon diagrams, shown in Fig. 4. These are the most difficult diagrams we have to compute, since in general they diverge and a regularization is required. There arise two types of divergences: First, there are the ultraviolet (UV) divergences; they are removed by an appropriate renormalization. Second, there are infrared (threshold) singularities. They will remain in the final result for the hard scale contributions and vanish only in the sum with the soft contributions.
The evaluation of the hard corrections is made possible by a combination of the analytical and numerical methods. The idea is to construct an infrared safe expression from divergent Feynman amplitudes by subtracting appropriate counterterms, which can be computed analytically. The construction of the counterterms is based on the following observation: in a given Feynman diagram the infrared singularities appear when the loop momenta are small. In this situation the propagator of the virtual electron in the t-channel can be contracted to a point. As the result, the infrared behavior of those Feynman diagrams is identical to that of the two-loop three point functions considered earlier in [27,20].
We proceed in the following way: after constructing an infrared finite expression, we combine propagators using Feynman parameters; perform momentum integrations analytically; extract ultraviolet (UV) divergences and integrate numerically over typically 5 (in some cases 6) Feynman parameters in the finite expressions. For the numerical integration we use the adaptive Monte Carlo routine VEGAS [28].
Let us illustrate the basic steps of the calculation by considering as an example the nonplanar box diagram D 1 shown in Fig. 4. A power counting shows that this diagram is UV finite but IR divergent. To demonstrate how the IR counterterm is constructed we consider a symbolic expression for this diagram (after taking the trace over Dirac matrices): .
Here p = (m, 0) is the four momentum of the incoming electron or positron, l 3 = l 1 + l 2 is the sum of the loop momenta l 1,2 , and p 1 = p − q, where q is the four-momentum of the outgoing photon. The quantities f (i) in the numerator denote the uniform functions of the loop momenta: Only terms with f (0) and f (1) diverge in IR. We use the following identity: where the counterterm D (i=0,1) 1 ct is obtained by expanding the propagator of the electron in the t-channel in Taylor series in small loop momenta, Examining Eq. (44) one recognizes that the first two terms in that equation are finite, both in the UV and the IR, and hence can be evaluated numerically. The last term, D (i=0,1) 1 ct , is divergent. Since the t-channel propagator has been contracted to a point, this term corresponds to a three-point, rather than four-point Feynman amplitude. Such integrals were computed in a previous study [27] (see also [20]). Using those results one can obtain the counterterm D (i=0,1) 1 ct analytically: Similar procedure was applied to evaluate the remaining Feynman diagrams D i . In some cases, like e.g. for the planar box diagram D 2 , the overall subtraction is not sufficient and a more sophisticated approach is required. The results of the calculation are summarized in Table III. Finally, we have to consider the one-loop diagrams. We need the results including terms O(ǫ), because they will be multiplied by divergent renormalization constants. The results are summarized in Table I.
Here B i are the tree-level diagrams (Fig. 1), S i are one-loop diagrams (Table I), C i are the one-loop diagrams with mass insertions (Table II), and D i are the two-loop diagrams (Table  III). The final result for the second order non-logarithmic correction to the p-Ps decay rate into two photons is obtained as a sum of the soft (37) and hard (50) pieces, and the square of one-loop corrections (24). Adding them one finds and the theoretical prediction for the p-Ps lifetimes becomes Γ theory p−Ps = 7989.64(2) µs −1 .
In this equation we have not included the contribution of the decay mode p-Ps→ 4γ. This decay channel increases Γ p−Ps by approximately 0.01 µs −1 [30,31].

VI. CONCLUSION
We have described the numerical and analytical methods employed in our calculation of the second order QED corrections to the p-Ps lifetime. Our final result, the non-logarithmic correction O(α 2 /π 2 ) with the coefficient B p = 5.1(3) increases the decay rate by approximately 0.1 µs −1 . The resulting theoretical prediction is in excellent agreement with experiment.
The framework of this calculation is the Nonrelativistic QED with dimensional regularization. The dimensional regularization facilitates the separation of scales, the cornerstone of the effective theory. Its main technical advantage is that no additional scales, such as a photon mass, are introduced by the regularization, so that only single scale integrals need to be computed. Unfortunately, these integrals are still too complicated to be evaluated analytically. We computed them numerically by subtracting IR counterterms. Those were constructed using simpler integrals which are known analytically, as we explained in Section IV. The multidimensional finite integrals were computed using adaptive Monte Carlo integration routine VEGAS [28].
The approach described in this paper can be also applied to the calculation of the O(α 2 ) corrections to the o-Ps decay into three photons. In particular, the IR counterterms can be constructed in a similar manner. However, the larger number of diagrams and two additional integrations over the three-photon phase space make this problem significantly more difficult.
Somewhat surprising is the smallness of the second order corrections found in this paper. It resulted from a strong cancellation between soft and hard pieces computed in dimensional regularization. We would like to stress that the soft and hard pieces are not separately finite and depend on the regularization (in this sense they are "scheme-dependent"). For this reason, large constants accompanying divergent pieces in Eqs. (50,37) may have no direct physical meaning. Unambiguous information is provided by the scheme-independent results, Eqs. (24,31,39,41), which are all of order [several units] × (α/π) 2 .
In the absence of a complete result on O(α 2 ) to Γ o−Ps it is interesting to discuss what our result for p-Ps might imply for the o-Ps lifetime puzzle. Although nothing can be said rigorously, we believe that our result indicates that no dramatic enhancement in the O(α 2 ) effects in o-Ps decay is possible. In Table IV we have summarized available results on the second order effects for both p-Ps and o-Ps decays (we have not included the partial results for the soft corrections in o-Ps since they are scheme-dependent). One can see that, with the exception of B squared , all radiative corrections are comparable for o-Ps and p-Ps decays. On the other hand, significantly larger value of B squared for o-Ps can be traced back to a larger value of the one-loop correction to o-Ps→ 3γ rate. The relation of the one-loop corrections for p-Ps and o-Ps decay, however, is rather natural since the number of diagrams is approximately three times larger for o-Ps decay. Thus, apart from the factor related to the number of Feynman diagrams, there seems to be no significant difference in the structure of radiative corrections to o-Ps and p-Ps decays.
Therefore, it is difficult to imagine that a complete calculation of the O(α 2 ) correction to o-Ps decay will result in a dramatically large number, necessary to resolve the o-Ps lifetime puzzle. We believe that this puzzle will be solved by continuing experimental studies and we look forward to learning their results.

VII. ACKNOWLEDGMENTS
We are grateful to G. Adkins, R.
FIG. 4. Two-loop photonic diagrams for the p-Ps decay.     For the various coefficients B we use the notation introduced in the text. B nγ describes multiphoton processes, p-Ps → 4γ and o-Ps → 5γ, and B hard (LL) is due to effects of non-linear QED: γ * γ * → γγ for p-Ps and γ * → 3γ for o-Ps. Γ LO denotes the lowest order decay width of a given state.  [30,31]