Event generators for $\eta/\eta^{\prime}$ decays at BESIII

The light unflavoured meson $\eta/\eta^{\prime}$ 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 $\eta/\eta^{\prime}\to\gamma l^{+}l^{-}$, $\eta/\eta^{\prime}\to\gamma \pi^{+}\pi^{-}$, $\eta^{\prime}\to\omega e^{+}e^{-}$, $\eta\to\pi^{+}\pi^{-}\pi^{0}$, $\eta/\eta^{\prime}\to\pi^{0}\pi^{0}\pi^{0}$, $\eta^{\prime}\to\eta\pi\pi$ and $\eta^{\prime}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}/\pi^{+}\pi^{-}\pi^{0}\pi^{0}$, which have been developed for investigating $\eta/\eta^\prime$ 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 $\eta/\eta^\prime$ physics.


Introduction
Light unflavoured meson spectroscopy is an effective tool to study the microscopic structure of matter and non-perturbation quantum chromodynamics (QCD), test Standard Model and probe the new physics beyond Standard Model [1]. The mixed unflavoured pseudoscalar mesons η/η ′ plays a pivotal role in studying the interactions between light quarks and the interactions between quarks and gluons [2]. The η/η ′ physics also offers us an unique place to test fundamental symmetries in QCD at low energy region, determine the difference of light quark masses [3]. η/η ′ mesons involve the mixing of a flavour octet state and a singlet state, with a mixed angle about −20 • [4]. Different with the prediction of flavour symmetry breaking, a larger probability to be singlet state brings a larger mass of η ′ , because the axial anomaly only contributes to the singlet mass. The uncertain gluon content size may also dedicates to the larger η ′ mass, which is related to vacuum topology and the U(1) anomaly [5]. The chiral perturbation theory (ChPT) and vector meson dominant models (VMD) are used to describe the different decay modes of η/η ′ .
The dominant decay mode of η is diphoton decay which is understood well and usually used in the determination of the η flux in many experiments. With one of the photons being virtual and decaying into lepton pair, the Dalitz decays occur. The second dominant η decay mode is 3π 0 , which violates isospin invariance and is not possible in the strong interaction. The strong decays to 2π or 4π are also forbidden in lowest order by P and CP invariance. The most important radiated decay of η is γπ + π − , while γπ 0 , γ2π 0 , γ3π 0 are forbidden by C invariance [6]. The dominant decay modes of η ′ can be classified as two kinds: one of which is the modes with 3 pseudoscalar mesons, another of which is the radiate decays with vector meson, and with the photon being virtual and decaying into lepton pair, the V e + e − decay modes are obtained.
In more than 50 years since η/η ′ mesons are discovered [7], the branching fractions of some of the dominant decay modes are measured [8], while vast of the decay modes are still waiting to be observed. Quite a number of up limits results on η/η ′ decays show the vast of interest on this sector and several η/η ′ factories are running or under construction. Low energy η/η ′ facili-ties include e + e − colliders such as KLOE-2 at DAφNE and BESIII at BEPCII, fixed-target experiments such as WASA at COSY and photoproduction experiments such as Crystal Ball at MAMI and CBELSA/TAPS at ELSA. High energy η/η ′ facilities also intend for the stage. Two programs, GlueX in Hall D and CLAS12 in Hall B, at the 12 GeV upgraded JLab will also dedicate their researches on the light meson spectroscopy with photons. The first result at GlueX on π 0 and η photoproduction has come out and more analysis are under study [9]. The K ± π ∓ nonleptonic weak decay is searched at BESIII [10] and the π + π − CP-violating strong decays are studied at LHCb [11]. A new Fermilab project REDTOP based on η meson is on the stage to search for new physics and symmetry violations [12].
The world's largest sample of 1.3 × 10 9 J/ψ events produced at rest and collected with the BESIII detector therefore offers a good opportunity to study the physics issues involved in the decays of η/η ′ . Recently, BE-SIII has reported several studies of the η/η ′ decays [13][14][15][16][17]. Experimentally, to measure a decay or a scattering process, Monte Carlo simulation should be run firstly, of which generator is one of the essential components. It is essential to build dynamic generator for the rare decays of η/η ′ in order to simulate the process as close as possible to the physics, improve the detection efficiencies and suppress backgrounds. Reference [18] has offer reliable event generators for η/η ′ → π + π − e + e − /π + π − µ + µ − /e + e − µ + µ − , however, more event generators are required.

Software framework for event generators
Monte Carlo (MC) simulations are used to determine the mass resolutions and detection efficiencies. 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 [19]. The charmonium state, like J/ψ, is simulated with the MC event generator kkmc [20,21], while the decays are generated by BesEvtGen [22] for known decay modes with branching fractions set to the Particle Data Group (PDG) world average values [8], and by lundcharm [23] for the remaining unknown decays. The event generators are based on the framework of the BESIII offline software system (BOSS) [24].
In general, the MC simulation on the 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 [18]. In this paper, the event generators are simplified by generating the phase space events firstly, and then a accept-reject algorithm is used to simulate the decays of η/η ′ with the Doing It Yourself (DIY) model in BesEvtGen [22]. The accept-reject algorithm is determined by the square of the amplitudes described in the next section.

Theoretical formulas and simulations
3.1 η/η ′ → γl + l − Electromagnetic Dalitz decays of light pseudoscalar mesons, such as, P → γ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. The decay, which a pseudoscalar meson P (η or η ′ ) decaying into a photon γ and a lepton-lepton pair, is the so-called single off-shell decay or Dalitz decay. In this case, the P decays into two photons first, then one of the photons is off-shell (γ * ) decays into the l + l − pair. The four-momenta for the process P → γ(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 [25] 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 xxx-2 where α = 1/137, f π = 92.4MeV, f 0 = 1.04f π , f 8 = 1.3f π and θ mix = −20 • . The vector meson dominance factor is where Γ is the width of the vector meson which satisfies and Θ function is where s is defined as s = p 2 = (p 1 +p 2 ) 2 . In the DIY generators, the parameters are set to be m ρ = 775.49MeV, m π 0 = 134.98MeV, and Γ ρ = 149.1 MeV. The different vector meson dominance models can be switched by inserting different values of c 3 .
The DIY event generators are performed in the hidden gauge case (c 3 = 1). η → γe + e − , η ′ → γe + e − ,η → γµ + µ − and η ′ → γµ + µ − decays are simulated by Be-sEvtGen package at the truth level, when the detector response, the detector construction, and the interaction between the particles and the materials are ignored. The invariant masses of the dileptons obtained from the MC simulations are shown in Fig. 1. The phase space model from η → γe + e − is taken as an example in Fig. 1 (a), comparing with the DIY model, and the discrepancy is very obvious, which indicates that reliable event generators are very important for the study of η/η ′ physics. The typical strong peaks around the small masses are emerged clearly in Figs 2) and (d) in the invariant mass spectra from the decays η ′ → γe + e − and η ′ → γµ + µ − , indicating the contribution from the vector meson, which are consistent with the theoretical results [25]. Especially for the η ′ → γe + e − channel, the shape of M (e + e − ) is also consistent with the BESIII result [14], which proves that the DIY event generator is performing very well.
The decay η/η ′ → γπ + π − is governed by the so called box anomaly proceed [25]. For the η(P ) → γ(k)π + (p + )π − (p − ) decay, the unpolarized squared decay amplitude in the DIY generator is based on [ where the variables in the formula are s ππ = (p + + p − ) 2 , the p π + p π − rest frame with respect to the direction of the flight of the p π + p π − in the pseudoscalar rest frame, and λ(m 2 P , s ππ , with and where Γ(s ππ ) is similar with Eq. (5), when m l ± is replaced by m π ± . The electric form factor E G is (11) where the form factor in item 2 is used in the generator with e = √ 4πα, F (s ππ ) ∼ F (0) = 0.19, and G = 1. The decay η ′ → γπ + π − is described by η ′ → γρ with the subsequently decay of ρ → π + π − in early years. The squared decay amplitude of η ′ → γπ + π − is similar with Eq. (7). However, the ρ mass extracted from dipion mass spectrum by different experiments was about 20 MeV larger than that from e + e − annihilation [26]. This effect is accounted for the higher term of Wess-Zumino-Witten (WZW) ChPT Lagrangian describing the non-resonant part of the coupling [27]. A simple Breit-Wigner function for ρ in the form factor is not enough in Eq. (10). As Ref [13] shown, a DIY event generator for η ′ → γπ + π − is introduced, taking the ρ − ω interference and the box anomaly into account. Deduced from the ones used in Refs [28], the decay rate formula can be written as where k γ is the photon energy and q π (m) is the momentum of pion in the π + π − rest frame. |δ| represents the contribution from ω resonance and the complex phase of δ represents the interference between ω and ρ(770) resonance. m ρ is the mass of the ρ(770) resonance. β is the constant ratio, which represents the non-resonant contribution. BW ω represents a simple Breit-Wigner function for ω, BW GS ρ is the Breit-Wigner distribution in GS parametrization [29] BW GS ρ = where and π is the momentum of pion in the π + π − rest frame with m π = 139.57MeV, q π (m 2 ρ ) = m 2 ρ /4 − m 2 π is the momentum of pion in the π + π − rest frame with m = m ρ , and The MC samples are generated to check the DIY event generators at truth level. The distributions from the MC samples are shown in Fig. 2. In the generator, η → γπ + π − is simulated in the hidden gauge case (c 3 = 1). The shapes of the invariant masses are consistent with theoretical results. For the simulation of η → γπ + π − decay, the ρ − ω interference and the box anomaly are taken into account and the shape of M (π + π − ) is in tune with the BESIII result [13].
It is interesting to study the decay η ′ → V e + e − (V represents 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 η ′ meson and momentum dependence of the transition form factor. In 2013, BESIII reported the measurement of η ′ → π + π − e + e − [13], which is found to be dominant by η ′ → ρe + e − , in agreement with the theoretical predictions [30,31]. More recently, the branching fraction measurement of B(η ′ → ωe + e − ) = [1.97±0.34±0.17]×10 −4 is reported for the first time via J/ψ radiative decays [17], which is consistent with theoretical predictions [30,32].
In the DIY event generator for η ′ → ωe + e − , the square of the amplitude can be written as [30] |A(P → V e + e − )| 2 where p 2 = (p l + + p l − ) 2 , p l ± is the four-momentum of the dilepton; Γ P →γV is the branching ratio of η ′ → γω; m P , m V are the mass of the pseudoscalar and vector meson in the process of P → γV ; β p = 1 − 4m 2 l ± p 2 ; θ p is the polar angle of p l ± in the p l + p l − rest frame with respect to the direction of the flight of the p l + p l − in the pseudoscalar rest frame. The vector meson dominance form factor can be written as The MC sample is generated to check the DIY event generator at truth level. The distribution from the MC samples is shown in Fig. 3. The shape of the invariant mass spectrum is consistent with theoretical result and the BESIII result [17].

η
The decays η/η ′ → 3π violate G parity and are induced dominantly by the strong interaction via u − d quark mass difference. It offers an ideal laboratory for testing ChPT and provides validation of models for π−π final state interaction [33,34]. BESIII has measured branching fractions of η/η ′ → 3π both charged and neutral channels [35], and a Dalitz plot analysis is also reported [15]. The internal dynamics of charged decay channel (η → π + π − π 0 ) can be described by two independent Dalitz plot variables which can be defined as [36] xxx-5 where T π denotes the kinetic energy of a given pion in the η rest frame, Q = m η −m π + −m π − −m π 0 is the excess energy of the reaction, and m η/π are the nominal masses from the PDG [8]. 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. The parameters in the DIY event generators are set for the case of charge conjugation invariance, where the Dalitz plot matrix elements for η → π + π − π 0 are taken from BESIII result [15] a = −1.128, where the parameters c and e are set to be zero, which indicates that no charge symmetry breaking here [15].
The neutral decay channels (η/η ′ → π 0 π 0 π 0 ) can be described by just one free parameter due to the symmetry of three neutral pions [37] and the parametrization of decay amplitude reads where T i denotes the kinetic energies of each π 0 in η/η ′ rest frame, Q = m η/η ′ − 3m π 0 , N is the normalized factor, and α is the slope parameter of the Dalitz plot. A nonzero α indicates final-state interactions. The value of the Dalitz plot slope parameter α is also taken from BESIII result [15] α = −0.055, for η; −0.640, for η ′ .
The MC samples are generated to check the DIY event generators at truth level. The distributions from the MC samples are shown in Fig. 4. The shapes of the distributions are consistent with the BESIII results [15]. (d) Fig. 4. (a) is the X distribution from η → π + π − π 0 , (b) is the Y distribution from η → π + π − π 0 , (c) is the Z distribution from η → π 0 π 0 π 0 and (d) is the Z distribution from η ′ → π 0 π 0 π 0 . xxx-6

η ′ → ηππ
The related matrix elements of η ′ → ηππ have been measured by many experiments [38]. The Dalitz plot distribution for the charged channel η ′ → ηπ + π − is described by the following two variables For the neutral channel η ′ → ηπ 0 π 0 , due to the symmetry of the two π 0 , the variable X is replaced by where T π , T η are the kinetic energies of the mesons in η ′ rest frame and Q = m η ′ − m η − 2m π is the excess energy of the decay, m η/π are the nominal masses in PDG.
The squared absolute value of the decay amplitude in the Dalitz plot is "general parameterized", defining as (29) where N is a normalization factor, and a, b, c and d are real parameters.
An alternative way to parameterize the squared amplitude is call the "linear parameterized" one where the complex parameter α, could be compared with the general parameterized one 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 if b > a 2 /4. The MC samples are generated to check the DIY event generators at truth level. The parameters are set for the "general parameterized" case, when the elements for η ′ → ηπ + π − are taken from BESIII result [38] a = −0.047, and these values are also used in the simulation for the neutral channel η ′ → ηπ 0 π 0 . The distributions from the MC samples are shown in Fig. 5. The shapes of the Y variables in Fig. 5 (b) and (d) are similar. The X variables in Fig. 5 (a) and (c) are different because of the different formulas for X variables.  5. (a) is the X distribution from η ′ → ηπ + π − , (b) is the Y distribution from η ′ → ηπ + π − , (c) is the X distribution from η ′ → ηπ 0 π 0 and (d) is the Y distribution from η ′ → ηπ 0 π 0 . xxx-7 3.6 η ′ → π + π − π + π − /π + π − π 0 π 0 In the ChPT, η ′ → 4π is believed to be governed by the WZW term via chiral anomalies. Recently Guo, Kubis, and Wirzba calculated the branching fractions of the η ′ decays into four pions based on a combination of chiral perturbation theory and vector-meson dominance model [39], and the theoretical results are obtained to be B(η ′ → π + π − π + π − ) = (1.0 ± 0.3) × 10 −4 , B(η ′ → π + π − π 0 π 0 ) = (2.4 ± 0.7) × 10 −4 . In 2014, BESIII reported the first observation of η ′ → π + π − π + π − and η ′ → π + π − π 0 π 0 decays coming from J/ψ → γη ′ radiative decay events [16]. The measured branching fractions are B(η ′ → π + π − π + π − ) = [8.53 ± 0.69 ± 0.64] × 10 −5 and B(η ′ → π + π − π 0 π 0 ) = [1.82±0.35±0.18]×10 −4 , which are consistent with theoretical predictions. The decay amplitudes in the DIY event generators are based on Ref. [39]. The four-momenta are defined as The amplitudes can be described in terms of the invariant variables s ij = (p i + p j ) 2 , i, j = 1, ..., 4, which are subject to the constraint (in the isospin limit of equal pion masses) The five-meson vertices of the WZW term can be deduced from the Lagrangian where ϕ∂ µ ϕ∂ ν ϕ∂ α ϕ∂ β ϕ represents the trace in flavour space, f π = 92.4MeV is the pion decay constant, and N c = 3 is the number of colors. Meanwhile, The η − η ′ mixing is described as The decay amplitudes are given by D ρ (s) is the inverse ρ propagator. Γ ρ is the width of ρ, which is set to be 149.1 MeV in the generator. In this DIY model, the coupling constants are assigned to be [39] The MC samples are generated to check the DIY event generators at truth level. The distributions from MC samples are shown in Fig. 6. The shapes of the invariant mass spectra are consistent with theoretical results and the BESIII result [16]. (c) Fig. 6. The invariant mass spectra obtained from the MC simulations. (a) is from η ′ → π + π − π + π − (there are 4 entries per event); (b) and (c) are from η ′ → π + π − π 0 π 0 . xxx-8
The amplitudes of η/η ′ decays calculated by ChPT and VMD models have been examined by the consistence between our simulation distributions and the BESIII results so that they will be useful in other η/η ′ factories, such as GlueX and CLAS12 experiments. It is no doubt that the event generators will provide powerful tools to explore η/η ′ physics. These generators will also be improved by new experiment results. The updating measurements of η and η ′ decays will give deep insights into the nature of pseudoscalar mesons and non-perturbation QCD, aspects of Standard Model and new physics beyond Standard Model.