Simulation of electron-positron annihilation into hadrons with the event generator PHOKHARA

The precise determination of the cross section for electron-positron annihilation into hadrons is one of the central tasks of ongoing experiments at low energy colliders. These measurements have to be complemented by Monte Carlo generators which simulate a large number of final states and include higher order radiative corrections. With this motivation in mind the generator PHOKHARA is extended to version 8.0, thus allowing for the simulation of final states with zero, one or two real photons. At the same time corrections from the emission of one or two virtual photons are included, such that a full next-to-next-to leading order generator is available. The stability and consistency of the program is tested. The results (for muon-pair final states) are compared to the programs KKMC and MCGPJ and implications for the analysis of various hadronic final states are investigated.


Introduction
The importance of a precise measurement of the total cross section for electron-positron annihilation into hadrons has been emphasised on many occasions (For recent reviews see e.g. [1][2][3][4][5]). Its low energy behaviour governs the running of the electromagnetic coupling from the Thompson limit to higher energies and is, therefore, a decisive input for all precision analyses of electroweak interactions. Indeed, it may well be one of the limiting factors for the interpretation of the Standard Model, in particular the test of the relation between masses of the top quark, the W and the Higgs boson [6,7]. It is, furthermore, one of the limiting ingredients for the theory prediction of the anomalous magnetic moment of the muon. Last not least, this cross section gives crucial input to dedicated analyses based on perturbative QCD (pQCD), be it the determination of the strong coupling constant α s (for a review see [8,9]), quark mass determinations [10] or low energy quantities like the pion or nucleon form factors. At high and intermediate energies a purely perturbative treatment of the total cross section is assumed to provide sufficiently precise predictions, however at low energies and in the charm-and bottom-quark threshold region no ab-initio prediction based on perturbative QCD (pQCD) can be made.
To determine the total cross section, one may either perform an inclusive measurement, or identify individually all the different multi-hadron final states. From the theoretical side the inclusive cross section can be well predicted on the basis of perturbative QCD (at least for properly chosen energy regions), exclusive channels carry more detailed information on form factors, resonances, isospin symmetry and breaking. In view of the importance of these measurements at low energies, where the number of different exclusive modes is still limited, and considering the fact, that this region plays the dominant role for the aforementioned theoretical investigations, it seems useful to develop a Monte Carlo generator which is tailored to the simulation of exclusive hadronic final states.
The measurement of the annihilation cross section proceeds in two conceptually different ways. The traditional and most obvious method is based on a variable center of JHEP08(2013)110 mass energy of the electron-positron collider which allows the experiment to scan through an energy range which is trivially dictated by the available beam energy (scanning mode). As a second alternative one may use the "radiative return", exploiting the fact that initial state radiation leads to a variable "effective" energy and invariant mass of the hadronic system. In view of the large luminosity of current electron-positron colliders, the loss in cross section, resulting from the factor of α to be paid for the photon emission, is compensated by the advantage of running at a stable fixed beam energy.
Various generators have been developed to simulate events with radiative return into a variety of hadronic and leptonic final states in leading order (LO) [11,12] and NLO [13][14][15]. Also for the scanning mode a number of precision generators have been developed which include electromagnetic corrections from initial state radiation in next-to-leading (NLO) and partly even next-to-next-to-leading (NNLO) order [16][17][18][19][20]. Note, that the counting of orders is different for the scanning mode and the radiative return, since one photon emission constitutes the leading order process in the latter, while this same process contributes to the NLO corrections in the scanning mode. It is the aim of this work to construct a new NNLO generator for the scanning mode, which is based on the already existing PHOKHARA 7.0, 1 as far as radiative corrections from one and two photon emission are concerned. Evidently, the only missing ingredient, the two loop virtual corrections, are available in the literature since long [21] and can be implemented into our generator in a straightforward way.
At present, all generators simulating the scanning mode have been constructed for leptonic final states only and to a limited extent for two body (ππ and KK) hadronic states. In contrast, the Monte Carlo event generator PHOKHARA has been designed from the beginning to simulate exclusive hadronic final states, using specific models for the form factors. Indeed many two-, three-and four-body final states have been implemented by now (For the detailed listing of the final states and the underlying assumptions about the form factors see section 2.) As said above, it is a fairly straightforward task to extend PHOKHARA 7.0 such that it is applicable for the scanning mode, thus providing a generator which includes the full NNLO corrections, at least as far as initial state radiation is concerned. This program will be complementary to the two main generators currently in use, KKMC 4.13 [16,19] and MCGPJ [17,18,20]. The former is based on YFS exponentiation, thus includes the emission of an arbitrary number of soft photons and hence a certain class of corrections which may be important for the measurement of a rapidly varying cross section, in particular close to a narrow resonance. This feature is not implemented in PHOKHARA. On the other hand, in contrast to PHOKHARA 8.0 the currently available version 4.13 of KKMC does not include the full NNLO corrections in the region of hard real photon emission (Although these corrections, originally evaluated in [13] and recalculated in [22,23], were implemented in a private version of the program and used for analysis in [24].). Furthermore, KKMC 4.13 is restricted to leptonic final states and thus cannot serve for an analysis of the multitude of hadronic states mentioned above.

JHEP08(2013)110
MCGPJ, on the other hand, is built on the usage of structure function method, an approach which does not allow to simulate the proper kinematics in the case of hard, non-collinear photon emission. Furthermore, among the hadronic final states only π + π − , Let us conclude this section with a brief comment concerning final state radiation. Inclusion of one-photon real and virtual corrections is straightforward, as far as lepton final states are concerned. For two body hadronic states an ansatz based on point-like pions, kaons and protons has been implemented in PHOKHARA from the very beginning. This ansatz has been successfully compared to data [26] as long as relatively soft photons of several hundred MeV are concerned. No further comparison between data and model has been performed to our knowledge. For multi-hadron final states use of PHOTOS ( [27] and updates) is recommended.

The NNLO corrections
As stated above, we would like to develop an event generator for the processes e + e − → hadrons and e + e − → µ + µ − , which relies on fixed order formulae and simulates exact kinematics. The starting point for this development is the event generator PHOKHARA [14], and in fact its latest version PHOKHARA7.0 [28] has been used for this purpose. The generator profits from the dedicated and continuously updated hadronic currents and the well tested implementation of QED radiative corrections. The master formula for the fixed order calculations, which is now implemented in the code reads dσ(e + e − → hadrons + photons) = dσ(e + e − → hadrons) +dσ(e + e − → hadrons + one hard photon) +dσ(e + e − → hadrons + two hard photons) . (2.1) The present version 8.0 of the program is limited to initial state emission and only the photon emission from electron and positron is taken into account. Emission of real and virtual lepton pairs will be treated in a later version. The radiative corrections in dσ(e + e − → hadrons + one hard photon) are described in [13,29] and the real emission of two photons in [14]. The implementations of various hadronic modes with up-to-date hadronic currents modelling are discussed in [28] for π + π − , K + K − andK 0 K 0 final states, in [30] forpp andnn final states, in [31] for π + π − π 0 final state, in [32] for π + π − 2π 0 and 2π + 2π − final states and in [33] for Λ(→ π − p)Λ(→ π +p ) final state. Since the implementation of the ηπ + π − channel was never documented in the literature we add a description of this hadronic current in appendix (A). For the readers convenience the most important formulae are repeated in the following.
(i) Single photon emission in LO and NLO: the differential rate in dσ(e + e − → hadrons + one hard photon) is [14] cast into the product of a leptonic and a hadronic tensor and the corresponding factorised phase space: dσ(e + e − → hadrons + one hard photon) where dΦ n (Q; q 1 , ·, q n ) denotes the hadronic n-body phase space including all statistical factors and Q 2 is the invariant mass of the hadronic system. The physics of the hadronic system, whose description is model-dependent, enters only through the hadronic tensor where the hadronic current has to be parametrised through form factors, which depend on the four momenta of the final state hadrons. The form of the leptonic tensor, including the one-loop correction and emission of a second soft photon as well as the vacuum polarisation corrections, can be found in [14] (eq. (5)).
(ii) Two photon emission: for the evaluation of dσ(e + e − → hadrons+two hard photons) the helicity amplitude method was used as described in section 3 of [14]. This part contains as well the vacuum polarisation corrections.
(iii) Zero photon emission in LO, NLO and NNLO: the new part, dσ(e + e − → hadrons), added to the code is written as with the same hadronic tensor, but calculated at a different from eq. (2.2) kinematic point. The leptonic tensor contains the virtual and soft radiative corrections up to the second order [21] L 0 µν = 4 p 1µ p 2ν − g µν where ∆ V P (s) is the vacuum polarisation correction, ∆ = ∆ virt,1ph + ∆ sof t,1ph + ∆ virt,2ph + ∆ sof t,2ph + ∆ virt,sof t,1ph (2.6)  where: We use a bit different from [21] division of the two-photon phase space (see figure 1) into soft-soft, soft-hard and hard-hard parts. Hence the soft photon contribution coming from two soft photons and calculated analytically is also different. This division is more suitable for the implementation into our generator as we generate directly the photon energies.
To generate the phase space of the hadrons + two photons, the program generates first the invariant mass of the hadronic system q 2 and the angles of the photons in the centerof-mass-frame of the initial fermions. Then the energy of one of the photons is generated in this frame and the second is calculated from the relation where E 1 , E 2 are the energies of the photons and θ 12 is the angle between their momenta. The curve from eq. (2.14) gives the boundary of the allowed energies in figure 1.

JHEP08(2013)110
For the sum of all contributions we have tested the independence of the integrated cross section from the separation parameter between soft and hard parts. The results are summarised in figure 2. The recommended cut to be used is w = 10 −4 , while for w = 10 −5 we start to observe negative weights and the result is obtained using weighted events. Perfect agreement between the results for √ s = 3.65 GeV (in fact we have performed tests for several energies with similar results) demonstrates that the generation in the soft photon region is implemented properly. The small disagreement at √ s = 1.02 GeV is due to the presence of a relatively narrow resonance φ and it disappears when one further decreases the w cut. This difference together with the appearing of the negative weights also indicates that the missing multi-photon emission is important in this region. If one runs with w = 10 −4 the difference to w = 10 −5 is an estimate of additional error of the generator (the error estimate will be given at the end of this section).
Apart of many technical tests of the one-and two-photon parts we would like to recall a comparison presented in [24], where it was shown that the muon pair invariant mass distributions obtained with PHOKHARA 2.0 on the one hand and with KKMC upgraded for this purpose on the other hand are in excellent agreement. This is valid except in the kinematical region where the muon pair invariant mass is very close to √ s. In this case multi-photon emission is most important and leads to differences up to a few percent.
For any scan experiment one measures the photon-inclusive cross section with a loose cut on the invariant mass of the hadronic/muonic system, say q 2 min > s/2 . The agreement between the two codes KKMC 4.13 and PHOKHARA 8.0 is even better for this type of comparison. This is demonstrated in figure 3, where we show the ratio of the integrated cross sections   The difference between the predictions of PHOKHARA 8.0 and KKMC 4.13 for the inclusive cross sections for muon production depends only slightly on the center-of-mass energy as can be observed in figure 4. The difference is even smaller for MCGPJ. However, below 2 GeV the version of the MCGPJ code [34] we were using was not stable numerically when trying to improve the accuracy of the result. The indicated errors are the smallest we were able to get. For many experimental applications of an event generator, for example studies of acceptance corrections, the fully differential cross section is needed. A typical modification of an angular distribution by the radiative corrections is shown for muon polar angle distribution in figure 5. The NNLO corrections, even if they are small (figure 6), are relevant in the era of precision hadronic physics [1].
For some particular cases, when for example at a nominal energy of the experiments the form factor is small, the radiative corrections might be bigger than the LO result. It is shown in figure 7 for neutral kaons for the center of mass energy 1.2 GeV. The large correction reflects the large number of events from the radiative return to the Φ with JHEP08(2013)110    subsequent decay into neutral kaons, i.e. from e + e − → γΦ(→ K 0K0 ). For comparison we also show the corresponding plot for charged kaons. The charged kaon form factor at 1.2 GeV is significantly larger than the one for the neutral kaon and the relative impact of γΦ(→ K + K − ) is correspondingly smaller.
As another cross check of our new generator the differential cross section with respect to the missing transverse momentum (carried away by one or two photons) generated with PHOKHARA 8.0 and with KKMC 4.13 is shown in figure 8. The missing momentum is defined as the difference between incoming and outgoing fermion momenta. The small difference between these distributions can be attributed to the missing multi-photon corrections in PHOKHARA 8.0 generator. A comparison of the angular distributions given by KKMC 4.13 and PHOKHARA 8.0 codes is shown in figure 9. The substantial difference between the two codes seen in the left plot comes only from the region of invariant masses of the muon pairs close to the threshold and drops to permille level when the threshold region is excluded (right plot of figure 9). We attribute this difference to the missing mass terms in the version of KKMC 4.13 we use and expect they will disappear when the missing mass corrections are included as in [24].
From the studies and comparisons presented in this section one can conclude that the accuracy of the code, as far as ISR corrections are concerned, is from 0.3 permille for cumulative observables up to 2-3 percent for some specific differential distributions in regions of small event rates.

Conclusions
In its original form up to version 7.0 the generator PHOKHARA was constructed to simulate events with at least one photon from initial state radiation. In combination with the input based on the two-loop form factor [21] derived long time ago, the complete NNLO corrections from ISR are available and can be used to construct the corresponding generator, now version 8.0, for electron-positron annihilation into muon pairs or hadrons. This JHEP08(2013)110 generator is complementary to the two other Monte Carlo computer codes currently in use, KKMC and MCGPJ. In particular PHOKHARA 8.0 does not include multi-photon (beyond two) emission, however, it can be used to simulate a multitude of exclusive hadronic final states and, in contrast to structure function methods, the present approach includes exact kinematics in the one-and two photon emission. Furthermore already in its present version it includes a large number of exclusive hadronic final states.
In the paper we demonstrate the independence of the results from the soft photon cutoff and, for muon-pairs, compare a few selected distributions with the results from other programs. Furthermore we study the impact of NNLO compared to NLO corrections on the predicted cross section and its dependence on the cut on the mass of the hadronic system and find a strong dependence on the choice of the final state and the center-of-mass energy of the experiment.
Since our choice for the amplitude for the production of the three-meson state ηπ + π − has not yet been documented in the literature, a brief description is presented in the appendix.

Acknowledgments
We would like to thank Staszek Jadach for clarification of topics concerning KKMC-PHOKHARA comparisons in [24]. We would like to thank our experimental friends for a constant push towards finishing this project, especially Achim Denig. H.C. would like to thank Simon Eidelman for discussions concerning the ηπ + π − hadronic current. This work is a part of the activity of the "Working Group on Radiative Corrections and Monte Carlo Generators for Low Energies" [http://www.lnf.infn.it/wg/sighad/]. This work is supported by BMBF project 05H12VKE, SFB/TR-9 "Computational Particle Physics" and Polish Ministry of Science and High Education from budget for science for 2010-2013 under grant number N N202 102638. M. Gunia was supported byŚwider PhD program.
The hadronic current J µ em responsible for this reaction can be related to the corresponding charged current J µ cc governing the τ decay τ + →νπ + π 0 [35] Following [36,37] the amplitude is normalised to its chiral limit, which is predicted by the Wess-Zumino-Witten term [38,39]. Its behaviour away from this point is governed by a form factor F , which includes the dominant ρ resonance and its radial excitations and is normalised to one for s = (p η + p π − + p π + ) 2 = 0 and s 1 = (p π − + p π + ) 2 = 0. This leads to the following ansatz ε µ νρσ p ν π + p ρ π − p σ η F, (A.5) with F = 1 1 + e iφ 1 c 1 + e iφ 2 c 2 BW ρ 0 (s 1 ) (A.6) × BW ρ 0 (s) + BW ρ 1 (s)e iφ 1 c 1 + BW ρ 2 (s)e iφ 2 c 2 where s 1 is the π + π − system invariant mass, We performed an 8 parameter fit of the form factor (A.7) of the hadronic current (A.5) to the experimental data from BaBar [41], CMD-2 [42], DM1 [43], DM2 [44], ND [45] and SND. Table 1 contains values of the fit's parameters together with its χ 2 , the result of the fit is compared to the data in figure 10.  Open Access. This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.