Event generators for η / η ′ decays at BESIII

The light unflavoured meson η / η ′ decays are valuable for testing non-perturbative quantum chromodynamics and exploring new physics beyond the Standard Model. This paper describes a series of event generators, including η / η ′ → γ l + l − , η / η ′ → γ π + π − , η ′ → ω e + e − , η → π + π − π 0 , η / η ′ → π 0 π 0 π 0 , η ′ → η π π and η ′ → π + π − π + π − / π + π − π 0 π 0 , which have been developed for investigating η / η ′ decay dynamics. For most of these generators, their usability has been validated in BESIII analyses for determining the detection efficiency, and background studies. The consistency between data and Monte Carlo shows that these generators work well in the BESIII simulation, and will also be useful for ongoing BESIII analyses and other experiments for studying η / η ′ physics.


Introduction
As the ground states of pseudoscalar nonets, the η and η mesons have been firmly established and their main decay modes are fairly well known [1]. However, they still attract theoretical and experimental attention due to their special role in understanding low energy quantum chromodynamics (QCD), even though they were discovered about half a century ago [2]. The decays of η/η are of interest as probes of some aspects of strong interactions, and also as sources of information on physics beyond the Standard Model (SM).
Due to flavor symmetry breaking, η/η mesons involve the mixing of an octet state and a singlet state, with a mixing angle of about −20 • [3]. A larger probability to be a singlet state gives a larger mass to the η , because the axial anomaly only contributes to the singlet mass. The gluonic content in the η , which is related to vacuum topology and the U (1) anomaly [4], may also contribute to its large mass.
The unflavored η/η mesons not only play an important role in studying the interactions between light quarks and the interactions between quarks and glu-ons [5], but also offer a unique place to test fundamental symmetries in QCD in the low energy region. In addition, their hadronic decays could be used to determine the difference of light quark masses [6]. Therefore, η/η physics is listed in the programs of many experiments, such as BESIII, KLOE-2, MAMI, GlueX, and CLAS12. Most recently, a new facility, REDTOP [7], is also proposed to study η/η decays.
Due to the large production rate of η/η mesons in J/ψ hadronic and radiative decays, the world's largest sample of 1.31×10 9 J/ψ events [8], collected with the BESIII detector, offers a good opportunity to study the decays of η/η . In recent years BESIII has achieved much progress on η/η decays [9][10][11][12][13][14]. Experimentally, a good description of the amplitude for each decay mode in the Monte Carlo (MC) simulation plays an important role in the optimization of the event selection criteria, determination of the detection efficiency, and background suppression. In the near future, about 10 billion J/ψ events will have been accumulated at the BESIII detector, which will allow us to search for the rare decays of η/η at an unprecedented level. In that case, the established η/η decays become the dominant background sources and proper event generators are hence essential to the background estimation.
To meet the challenge of precision measurement of η/η physics, in this paper we present a series of event generators for their decays, including η/η → γl + l − , η/η →γπ + π − , η →ωe + e − , η→π + π − π 0 , η/η →π 0 π 0 π 0 , η → ηππ and η → π + π − π + π − /π + π − π 0 π 0 , within the framework of chiral perturbation theory (ChPT) and the vector meson dominance (VMD) model. In Section 2, a brief introduction is given to describe the software framework for event generators. In Section 3, the decay amplitudes for these decays and the corresponding parameters used for generating events are provided. Meanwhile, validations for some of these generators are also performed to ensure that they work well in the BESIII MC simulation package. A short summary is presented in the last section.

Software framework for event generators
At the BESIII experiment, the geant4-based simulation software boost includes the geometric and material description of the BESIII detector, the detector response, and the digitization models, as well as the detector running conditions and performance [15]. The charmonium state, e.g., J/ψ, is simulated with the MC event generator kkmc [16,17], while the decays are generated by BesEvtGen [18] for known decay modes, with branching fractions set to the world average values in the Particle Data Group (PDG) [1], and by lundcharm [19] for the remaining unknown decays. The event generators are based on the framework of the BESIII offline software system (BOSS) [20].
In general, the MC simulation in BOSS is performed by passing the Lorentz-vector of all the particles produced by the event generator into a simulation package of the BESIII detector after taking into account the detector construction, the detector response, the interaction between the particles, and the material [21]. In this paper, the event generators were developed in accordance with the corresponding η/η decay amplitudes, which are described in detail in Section 3, and then were implemented in BesEvtGen.

Theoretical formulas and simulations
3.1 η/η →γl + l − Electromagnetic Dalitz decays of η/η → γl + l − (l ± stands for µ ± or e ± ) play an important role in revealing the structure of hadrons and the interaction mechanism between photons and hadrons. This is the so-called single off-shell decay in which the l + l − pair originates from the off-shell photon (γ * ). The four-momenta for the pro-cess η/η → γ(k)γ * (p) → γ(k)l + (p 1 )l − (p 2 ) are defined as P =k+p=k+p 1 +p 2 . The square of the amplitude can be written as [22] |A(P →l + l − γ)| 2 where p=p 1 +p 2 , β p = 1− 4m 2 l ± p 2 and θ p is the helicity angle. M P (p 2 ,k 2 =0) is the form factor, which is described as where where α = 1/137, f π = 92.4 MeV, f 0 = 1.04f π , f 8 = 1.3f π and θ mix =−20 • [23]. The η/η →γl + l − form factor measurements support the theoretical prediction from the VMD model that the dominant contribution is from the ρ [10,24]. The VMD form factor for η → γl + l − in the generator is written as where s is defined as s=p 2 =(p 1 +p 2 ) 2 . Γ (s) is the width of the vector meson [25] and the Θ function is For η →γl + l − , the phase space allows production of ρ and ω. Even although the contribution from ω is small, it cannot be directly ignored, since its interference with ρ may lead to a sizable contribution. By combining with ρ and ω, the VMD form factor for η → γl + l − can be described by where BW ρ (BW ω ) represents a simple Breit-Wigner function for ρ(ω).
In the generator, the parameters in the above formula are set to be m ρ =775.49 MeV, m π =139.57 MeV, 013001-2 and Γ ρ = 149.1 MeV [1]. The weight factors w ρ and w ω are subjected to w ρ : w ω = 3 : 1 [26]. The different vector meson dominance models can be switched by inserting different values of c 3 [27,28].
By implementing the above amplitudes in the BESIII simulation package, the η/η → γl + l − events are generated. The mass spectra of the leptonic pairs (the solid histograms) at the truth level are shown in Fig. 1, where the VMD form factor used for generating η → γl + l − events is that for the case of the hidden gauge (c 3 = 1). For a comparison, we also generate the events with V M D(s)=1, as indicated by the hatched histograms in Fig. 1. The discrepancies between the events generated with and without VDM form factors are quite obvious, in particular for η/η → γµ + µ − , which also shows that reliable dynamic generators are very important for the study of η/η physics.

η/η →γπ + π −
The decay η/η →γπ + π − receives a contribution from the box anomaly [22]. For the η→γ(k)π + (p 1 )π − (p 2 ) decay, the unpolarized squared decay amplitude is [ p 2 , and θ p is the polar angle of p π ± in the p π + p π − rest frame with respect to the direction of the flight of the p π + p π − in the rest frame of the pseudoscalar meson. The electric form factor E G is used to describe the CP violation in the decay, which is set to be zero in this paper. The Källén function with and where Γ (s) is the same as in Eq. (5). At first, the decay η → γπ + π − was believed to be dominated by η → γρ with the subsequent decay 013001-3 ρ → π + π − . In this case, the squared amplitude of η → γπ + π − would be similar to Eq. (8). However, the ρ mass extracted from the dipion mass spectrum by different experiments is about 20 MeV/c 2 larger than that from e + e − annihilation [29]. This effect is accounted for by the higher term of the Wess-Zumino-Witten (WZW) ChPT Lagrangian describing the non-resonant part of the coupling [30]. A simple Breit-Wigner function for ρ in the form factor in Eq. (12) is not enough to describe the data. The ρ−ω interference and the box anomaly should be taken into account for η → γπ + π − . Deduced from the ones used in Ref. [31], the decay rate [9] can be expressed by where m 2 = (p π + +p π − ) 2 , and p π + and p π − are the fourmomenta in the laboratory frame. k γ is the photon energy and q π (m) is the momentum of the pion in the π + π − rest frame. The parameter δ is a complex number, for which |δ| = 5.59×10 −4 represents the contribution from the ω resonance and the complex phase of δ (argδ = −3.78 rad) represents the interference between the ω and the ρ(770) resonance [9]. m ρ is the mass of the ρ(770) resonance. β = −19.33 is the box anomaly constant ratio, which represents the non-resonant contribution [9]. BW ω represents a simple Breit-Wigner function for ω. BW GS ρ is the Breit-Wigner distribution in the GS parametrization [32], where and when π is the momentum of the pion in the π + π − rest frame with m π = 139.57 MeV [1], π is the momentum of the pion in the π + π − rest frame with m = m ρ , and Γ ρ (m) = At the truth level, Figs. 2(a) and (b) show the π + π − mass spectra for η → γπ + π − and η → γπ + π − , respectively. At the detector level, the generator for η → γπ + π − has been validated in the BESIII measurement of η →4π [12] by providing a good description of the background events from η →π + π − η(η→γπ + π − ). The generator for η →γπ + π − has also been used in the determination of the detection efficiency of J/ψ→γη (η →γπ + π − ) at the BESIII experiment. It was found that the π + π − mass spectrum from the MC simulation is consistent with that of data. More model-independent approaches can be found in Refs. [33,34].

η →ωe + e −
It is interesting to study the decay η → V e + e − (V represents a vector meson), which is related to the twobody radiative decay into a vector meson and an off-shell photon. The e + e − invariant mass distribution will provide us useful information on the internal structure of the η meson and the momentum dependence of the transition form factor. In 2013, BESIII reported the measurement of η → π + π − e + e − [9], which is found to be dominated by η → ρe + e − , and the results are in agreement with the theoretical predictions [35,36]. More recently, the branching fraction B(η → ωe + e − ) = [1.97±0.34±0.17]×10 −4 was reported for the first time via J/ψ radiative decays [13], and is also consistent with the theoretical predictions [35,37]. Within the framework of effective meson theory [35], the square of the amplitude of η →ωe + e − can be written as where p 2 =(p l + +p l − ) 2 , p l ± is the four-momentum of the dilepton; Γ P →γV is the decay width of η →γω; m P and m V are the mass of the pseudoscalar and vector meson respectively in the process P →γV ; β p = 1− 4m 2 l ± p 2 ; and θ p is the polar angle of p l ± in the p l + p l − rest frame with respect to the direction of flight of the p l + p l − in the pseudoscalar rest frame. The vector meson dominance form factor can be written as where m ω (m φ ) and Γ ω (Γ φ ) are the mass and width, respectively, of ω(φ) in the PDG [1]. BW ρ (BW ω ) represents a simple Breit-Wigner function for ρ(ω). The weight factors w ω and w φ are subjected to w ω : w φ = 1 : 4 [38].  Figure 3 shows the e + e − invariant mass distribution at the truth level. In the observation of η →e + e − ω [13], the MC events generated with this generator provided a good description of the data.
3.4 η→π + π − π 0 ,η/η →π 0 π 0 π 0 The decays η/η → 3π violate G parity and are induced dominantly by the strong interaction via the u−d quark mass difference. It offers an ideal laboratory for testing ChPT and provides validation of models for the π-π final-state interaction [39]. BESIII has measured the branching fractions of η/η → 3π for both charged and neutral channels [40], and a Dalitz plot analysis is also reported [11]. The internal dynamics of charged decay channel (η→π + π − π 0 ) can be described by two independent Dalitz plot variables [41] where T π denotes the kinetic energy of a given pion in the η rest frame, Q=m η −m π + −m π − −m π 0 is the kinetic energy in the reaction, and m η/π are the nominal masses from the PDG [1]. The decay amplitude of η→π + π − π 0 can be parameterized as where N is a normalization factor and the coefficients a,b,c,... are called Dalitz parameters, when a non-zero value for c or e may imply the violation of charge conjugation.
Since no evidence of the charge-conjugation violation is seen in the previous measurement, we did not take this effect into account in the generator. The parameters for η → π + π − π 0 taken from the BESIII measurement [11] For convenience, we also provide an option by including the item X 2 Y in the generator, and the parameters are from the KLOE-2 measurement [42]:  For η/η → π 0 π 0 π 0 , the density distribution of the Dalitz plot has threefold symmetry, due to the three identical particles in the final states. Hence, the density distribution can be parameterized using the polar variable [43], and the parametrization of the decay amplitude is given by [44] |A(Z)| 2 =N (1+2αZ+2βZ 3/2 sin(3φ)...), where T i denotes the kinetic energies of each π 0 in the η/η rest frame, Q=m η/η −3m π 0 , φ=arctan(Y /X), and N is the normalized factor. α and β are the Dalitz plot parameters. A nonzero α indicates final-state interactions. β has not been measured yet, so it is set to zero in the generator, while the value of α is taken from the BESIII measurement [11], α= −0.055, for η→π 0 π 0 π 0 ; −0.640, for η →π 0 π 0 π 0 .
At the truth level, the distributions of the Dalitz plot variables from the MC simulation are shown in Fig. 4.

η →ηππ
The matrix elements of η → ηππ have been measured by many experiments [45]. The Dalitz plot for the charged channel η → ηπ + π − can be described by two variables, For the neutral channel η →ηπ 0 π 0 , due to the symmetry of the two π 0 , the variable X is replaced by where T π and T η are the kinetic energies of the mesons in the η rest frame and Q=m η −m η −2m π is the kinetic energy in the decay, with m η/π the nominal masses in the PDG [1]. For the general representation, the squared amplitude is parameterized as where N is a normalization factor, and a, b, c and d are real parameters. A non-zero c parameter indicates violation of C parity for η → ηπ + π − and violation of Bose symmetry for η →ηπ 0 π 0 . An alternative parameterization is the so-called "linear representation", which is written as follows: where the complex parameter α can be compared with the general parameterization with a = 2Re(α) and b = Re 2 (α)+Im 2 (α). The real component of the complex constant α is a linear function of the kinematic energy of the η. The two parameterizations are equivalent in the case of b>a 2 /4. In the case of the general representation, the parameters in the generator for η →ηπ + π − are taken from the BESIII measurement [45], which are: These values are also used in the simulation for the neutral channel η →ηπ 0 π 0 . With the above generator, the distributions from the MC simulation at the truth level are shown in Fig. 5. The shapes of the Y variables in Fig. 5(b) and Fig. 5(d) are similar as a result of isospin symmetry. There is a deviation between Fig. 5(a) and (c) due to different kinematics for X (no dynamics in the matrix elements). The physical cusp effect was not taken into account in the generator because no evidence has been observed yet [46].
The η−η mixing is described as The mixing angle is θ mix =−20 • [23]. The decay amplitudes are then given by where f π =92.4 MeV is the pion decay constant [23], and N c =3 is the number of colors.
At present, the world's largest sample of J/ψ events, 1.31 × 10 9 events collected with the BESIII detector, provides a unique opportunity to investigate η/η decay dynamics. Besides the achievements obtained from the J/ψ radiative or hadronic decays into η/η [9][10][11][12][13][14], many analyses on η/η physics are in progress, which offer an opportunity to further evaluate the usability of these generators by examining whether they can provide a good description of data. In addition to the BESIII experiment, these event generators could also be a useful tool to investigate η/η decays in other experiments, e.g., GlueX, CLAS12, and KLOE-2.
Nian Qin thanks Dr. Xiao-Lin Kang, Dr. Xin-Ying Song, Dr. Li-Qing Qin and Ms. Hui-Juan Li for their support in developing the generators and helpful discus-sions.