Partial wave analysis of $J/\psi\rightarrow\gamma\eta\eta'$

Based on a sample of (10.09$\pm$0.04)$\times$10$^{9}$ $J/\psi$ events collected with the BESIII detector operating at the BEPCII storage ring, a partial wave analysis of the decay $J/\psi \rightarrow \gamma\eta\eta'$ is performed. An isoscalar state with exotic quantum numbers $J^{PC}=1^{-+}$, denoted as $\eta_1(1855)$, has been observed for the first time with statistical significance larger than 19$\sigma$. Its mass and width are measured to be (1855$\pm$9$_{-1}^{+6}$)~MeV/$c^{2}$ and (188$\pm$18$_{-8}^{+3}$)~MeV, respectively. The product branching fraction ${\cal B}(J/\psi$$\rightarrow$$ \gamma\eta_1(1855)$$\rightarrow$$\gamma\eta\eta')$ is measured to be (2.70$\pm 0.41 _{-0.35}^{+0.16}) \times$10$^{-6}$. The first uncertainties are statistical and the second are systematic. In addition, an upper limit on the branching ratio ${\cal B}(f_0(1710)$$\rightarrow$$\eta\eta')$/${\cal B}(f_0(1710)$$\rightarrow$$\pi\pi)$ is determined to be $1.61 \times 10^{-3}$ at 90\% confidence level, which lends support to the hypothesis that the $f_0(1710)$ has a large glueball component.

a Also at the Moscow Institute of Physics and Technology, Moscow 141700, Russia b Also at the Novosibirsk State University, Novosibirsk, 630090, Russia c Also at the NRC "Kurchatov Institute", PNPI, 188300, Gatchina, Russia Confinement is a unique property of quantum chromodynamics (QCD), and can be probed via the spectrum of mesons.While the quark model describes a conventional meson as a bound state of a quark and an antiquark, lattice QCD (LQCD) and QCD-motivated models predict a more rich spectrum of mesons that includes bound states with gluonic degrees of freedom, such as glueballs and hybrids.Radiative decays of the J/ψ meson provide a gluon-rich environment and are therefore regarded as one of the most promising hunting grounds for gluonic excitations [1][2][3][4].
In this paper, based on a sample of (10.09±0.04)×10 9J/ψ events collected with the BESIII detector [37], we present a partial wave analysis of J/ψ → γηη ′ to search for 1 −+ and investigate the decay property of f 0 (1710).The η is reconstructed via the decay channel γγ, and the η ′ is reconstructed via the decay channels η ′ → γπ + π − and η ′ → ηπ + π − .This paper is accompanied by a letter submitted to Physical Review Letters [38].

II. BESIII DETECTOR AND MONTE CARLO SIMULATION
The BESIII detector [39] records symmetric e + e − collisions provided by the BEPCII storage ring, which operates with a peak luminosity of 1 × 10 33 cm −2 s −1 at the centerof-mass energy 3.89 GeV.BESIII has collected large data samples in the energy region between 2.0 and 4.9 GeV [40].The cylindrical core of the BESIII detector covers 93% of the full solid angle and consists of a helium-based multilayer drift chamber (MDC), a plastic scintillator time-of-flight system (TOF), and a CsI(Tl) electromagnetic calorimeter (EMC), which are all enclosed in a superconducting solenoidal magnet providing a 1.0 T (0.9 T in 2012) magnetic field.The solenoid is supported by an octagonal flux-return yoke with resistive plate counter muon identification modules interleaved with steel.The charged-particle momentum resolution at 1 GeV/c is 0.5%, and the dE/dx resolution is 6% for electrons The EMC measures photon energies with a resolution of 2.5% (5%) at 1 GeV in the barrel (end cap) region.The time resolution in the TOF barrel region is 68 ps, while that in the end cap region is 110 ps.The end cap TOF system was upgraded in 2015 using multigap resistive plate chamber technology, providing a time resolution of 60 ps [41][42][43].
Simulated data samples produced with a GEANT4based [44] Monte Carlo (MC) package, which includes the geometric description of the BESIII detector and the detector response, are used to optimize the event selection criteria, to determine detection efficiencies, and to estimate backgrounds.Signal MC samples for the process J/ψ → γηη ′ with the subsequent decays η → γγ and η ′ → ηπ + π − are generated uniformly in phase space.The decay η ′ → γπ + π − is simulated by taking into account both ρ − ω interference and the box anomaly [45].
An inclusive MC sample with 10.01×10 9 J/ψ decays is used to study backgrounds.The known decay modes are modeled with EVTGEN [46] by incorporating branching fractions taken from the Particle Data Group [14], and the remaining unknown decays are generated using the LUNDCHARM [47] generator.The simulation includes the beam energy spread and initial state radiation (ISR) in the e + e − annihilations modeled with the generator KKMC [48].Final state radiation (FSR) from charged particles is incorporated with the PHOTOS package [49].

III. EVENT SELECTION
Charged tracks are reconstructed from hits in the MDC and are required to have |cosθ| < 0.93, where θ is the polar angle defined with respect to the symmetry axis of the MDC in the laboratory frame.Tracks must approach within 10 cm of the interaction point in the beam direction and 1 cm in the plane perpendicular to the beam where the distances are defined in the laboratory frame.Each track is assumed to be a pion, and no particle identification is applied.
Photon candidates are required to have energy deposition above 25 MeV in the barrel region (|cosθ|<0.80)or 50 MeV in the end caps (0.86<|cosθ|<0.92).To exclude spurious photons caused by hadronic interactions and final state radiation, photon candidates must be at least 10 o away from any charged tracks when extrapolated to the EMC.To suppress spurious photons due to electronic noise or energy deposits unrelated to the event, candidate showers are required to occur within 700 ns of the event start time.
For the J/ψ → γηη ′ , η ′ → ηπ + π − channel, events are reconstructed with two oppositely charged tracks and at least five candidate photons.A six-constraint (6C) kinematic fit under the hypothesis J/ψ → γηηπ + π − is performed by constraining energy-momentum conservation and the masses of two pairs of photons to m η .If there is more than one combination, the combination with the minimum χ 2 6C is retained.The resulting χ 2 6C is required to be less than 45.To suppress backgrounds with four or six photons in the final state, 4C kinematic fits are performed by constraining energy-momentum conservation under the hypotheses J/ψ→4γπ + π − , J/ψ→5γπ + π − , and J/ψ→6γπ + π − .The χ 2 4C (5γπ + π − ) is required to be less than all possible χ 2 4C (4γπ + π − ) and χ 2 4C (6γπ + π − ).The ηπ + π − combination with the minimum |M (ηπ The invariant mass distributions of ηη ′ and the Dalitz plot for the selected γηη ′ candidate events from the two decay channels are shown in Fig. 1(c, d) and Fig. 1(e, f), respectively.Clear structures in the ηη ′ invariant mass spectrum are observed.
Potential backgrounds are studied using the inclusive MC sample of 10.01×10 9 J/ψ decays (as described in Sec.II).No significant peaking background is observed in the invariant mass distribution of the η ′ .Backgrounds are estimated by the η ′ sidebands in the data.The sideband regions of η ′ → ηπ + π − and η ′ → γπ + π − are defined as M (ηπ + π − ) ∈ [0.900, 0.920] [0.995, 1.015] GeV/c 2 and M (γπ + π − ) ∈ [0.903, 0.915] [1.000, 1.012] GeV/c 2 , respectively.The normalization factors for events in the two sideband regions are obtained by a fit to data of the invariant mass spectrum of M (ηπ + π − ) and M (γπ + π − ).The η ′ signal shapes are determined from the shapes of signal MC samples described with RooHistPdf [50], and the backgrounds are described by 2nd degree polynomial functions.The definition of the sidebands and the fit results are shown in Fig. 1(a) and Fig. 1(b).

A. Analysis method
Using the GPUPWA framework [51], a combined PWA fit is performed to the selected samples of J/ψ → γηη ′ , η ′ → γπ + π − and J/ψ → γηη ′ , η ′ → ηπ + π − .Quasi twobody amplitudes in the sequential radiative decay processes J/ψ → γX, X → ηη ′ and hadronic decay process J/ψ → ηX, X → γη ′ and J/ψ → η ′ X, X → γη are constructed using the covariant tensor amplitudes described in Ref. [52].Let A X be the amplitude for a J/ψ decay process including intermediate resonance X.For J/ψ radiative decays, the general form of A X is where the summation is over the number of independent amplitudes and for J/ψ hadronic decays, the general form of A X is where ψ µ (m 1 ) is the polarization four-vector for the J/ψ; e * ν (m 2 ) is the polarization four-vector for the photon; m 1 and m 2 are the spin projections of the J/ψ and photon, respectively; U µν k is the kth independent partial wave amplitude of J/ψ radiative decays to intermediate resonance X with coupling strength determined by a complex parameter Λ k ; and U µ k is the kth independent partial wave amplitude of J/ψ hadronic decays to intermediate resonance X with coupling strength determined by a complex parameter Λ k .The partial wave amplitudes U µν k and U µ k are constructed using the four-momenta of the reconstructed γ, η, η ′ .
The amplitudes for the J/ψ radiative decay processes J/ψ → γf 0 , γf 2 , γf 4 are given in Ref. [52].For J/ψ → γη 1 , where the η 1 is an isoscalar state with exotic quantum numbers J P C = 1 −+ , the η 1 can decay into ηη ′ in a Pwave [34,35] with two amplitudes: where g µν is the metric tensor, p ψ is the four-momentum of the J/ψ, q is the four-momentum of the radiative photon, and f (η1) is the propagator for the process η 1 → ηη ′ .Blatt-Weisskopf barrier factors [53][54][55] are included in the orbital angular momentum covariant tensors t.Due to the special properties (massless and gauge invariance) of the photon, the number of independent partial wave amplitudes for a J/ψ radiative decay is smaller than for the corresponding decay to a massive vector meson, the details are given in Ref. [52].
For the J/ψ hadronic decay processes J/ψ → Vη ′ , V → γη, where V is vector meson that has quantum numbers J P C = 1 −− , such as φ, ρ, ω and their excitations, the corresponding amplitude is where, ǫ µ ′ ν ′ αβ is the totally antisymmetric tensor, p 1 is the four-momentum of the η ′ from the J/ψ decay, T β Vη ′ is the orbital angular momentum covariant tensor of the process J/ψ → Vη ′ , and f γη is the propagator for the process V → γη.The subscript P indicates J/ψ → Vη ′ is in a P wave.The amplitude of J/ψ → Vη, V → γη ′ is analogously to Eq. 5.
For the process J/ψ → h 1 η ′ , h 1 → γη, the corresponding two independent amplitudes are where gµν , p h1 is the four-momentum of the h 1 , and f (h1) is the propagator for the process h 1 → γη.The superscript (2) on T indicates the orbital momentum between h 1 and η ′ is 2. The subscripts S and D indicate J/ψ → h 1 η ′ is in an S wave and D wave, respectively.The amplitudes for the process J/ψ → h 1 η, h 1 → γη ′ are analogously to Eq. 6 and Eq. 7.
In this analysis, resonance decays are described by a Breit-Wigner (BW) function, parametrized by a constant-width, relativistic BW propagator, where M and Γ are the mass and width of the intermediate resonance X, and √ s is the invariant mass of the ηη ′ , γη, or γη ′ system.
The complex coefficients of the amplitudes (relative magnitudes and phases) and resonance parameters (masses and widths) are determined by an unbinned maximum likelihood The events between the blue arrows are the selected signal events, and events between the red arrows on the same side are selected sideband events.The black points with error bars are data, the red solid line is the fit result, the blue dashed line is the η ′ signal shape from signal MC, and the green dashed line is the 2nd degree polynomial function for the background.(c,d) Invariant mass distributions of the ηη ′ for the selected γηη ′ candidates.The points with error bars are data and the red line shows the background events estimated from the η ′ sideband.(e,f) The corresponding Dalitz plots for the selected γηη ′ candidates.
fit to the data.The likelihood is constructed following a method similar to that used in Ref. [56].
The probability to observe the ith event characterized by the measurement ξ i , i.e., the measured four-momenta of the particles in the final state, is where ǫ(ξ i ) is the detection efficiency, Φ(ξ i ) is the standard element of phase space, and M (ξ i ) = X A X (ξ i ) is the ma-trix element describing the decay processes from the J/ψ to the final state γηη ′ .A X (ξ i ) is the amplitude corresponding to intermediate resonance X as defined in Eq. 1 and Eq 2.
The joint probability for observing N events in the data sample is For technical reasons, rather than maximizing L, −lnL is minimized, with for a given dataset.The third term is a constant and has no impact on the determination of the parameters of the amplitudes or on the relative changes of −lnL values.In the fitting, the third term will not be considered.
The free parameters are optimized by MINUIT [57].The normalization integral σ ′ is evaluated using MC with importance sampling [58,59].An MC sample of N gen is generated with signal events distributed uniformly in phase space.These events are put through the detector simulation, subjected to the selection criteria and yield a sample of N acc accepted events.The normalization integral is computed as: where the constant value of the phase space integral Φ(ξ)dξ is ignored.
Instead of modeling the background, the likelihood is defined by the signal PDF [Eq.11] and the contribution to the negative log-likelihood from background events in the signal region is removed by subtracting out the negative loglikelihood of events in the η ′ sideband region in proper proportion [60], i.e., where −lnL signal is the likelihood for the signal, lnL data is the likelihood calculated by Eq. 11 using the data sample, lnL i background is the likelihood calculated by Eq. 11 using the events of ith sideband, and w i is the normalization factor for background events in the ith sideband region, which is determined from the fit results of Fig. 1(a) and Fig. 1(b).
The number of fitted events N X for an intermediate resonance X is defined as: where N is the number of selected events after background subtraction, and is calculated with the same MC sample as the normalization integral σ ′ .The detection efficiency ǫ X for an intermediate resonance X is obtained by the partial wave amplitude weighted MC sample, A combined unbinned maximum likelihood fit is performed for the two decay channels by adding the negative loglikelihood of signal, −lnL signal , for J/ψ → γηη ′ , η ′ → ηπ + π − and that for J/ψ → γηη ′ , η ′ → γπ + π − together.In the combined fit, the two decay modes share the same set of masses, widths, relative magnitudes, and phases.
The product branching fraction of J/ψ → ηX, X → γη ′ and J/ψ → η ′ X, X → γη can be calculated in a similar way to Eq. 17.

B. PWA results
To construct a set of two-body amplitudes to use in the PWA fit, a "PDG-optimized" set of amplitudes is first determined.To describe the ηη ′ and γη (′) spectra, all kinematically allowed resonances with J P C = 0 ++ , 2 ++ , and 4 ++ (for the ηη ′ system) and J P C = 1 +− and 1 −− (for the γη (′) systems) listed in the PDG [14] are considered.Within the allowed phase space (PHSP) of the ηη ′ system, four additional states [the f 0 (2102), f 0 (2330), f 2 (2240), and f 4 (2283)] reported in Ref. [61] and an additional scalar state (the f 0 (1810)) reported in Ref. [62] are also considered.Table I shows the complete set of resonances considered from the PDG, Ref. [61], and Ref. [62].All possible sets of amplitudes corresponding to resonances listed in Table I are evaluated.The statistical significance for each resonance is determined by examining the probability of the change in negative log-likelihood values when this resonance is included or excluded in the fits, where the probability is calculated under the χ 2 distribution hypothesis taking into account the change in the number of degrees of freedom.The masses and widths of the resonances near ηη ′ mass threshold [f 0 (1500), f 2 (1525), f 2 (1565), and f 2 (1640)] as well as those with small fit fractions (<3%) are always fixed to the PDG [14] values.The mass and width of the f 0 (2330), which corresponds to a clear structure around 2.3 GeV/c 2 in the ηη ′ mass spectrum, are free parameters.All other masses and widths are also free parameters in the fit.The final PDG-optimized set of amplitudes is the combination where each included resonance has a statistical significance larger than 5σ.Results from the PWA fit using the PDG-optimized set of amplitudes, including the masses, the widths, and the statistical significances of each component, are shown in Table II, where the uncertainties are statistical only.
In the next step, a search is performed for additional resonances with , and 1 −− γη (′) , where the subscript labels the composition of the resonance, by individually adding each possibility to the PDG-optimized solution and scanning over its mass and width.The significance of each additional resonance at each mass and width is evaluated.The result indicates that a significant 1 −+ contribution (>7σ) is needed around 1.9 GeV/c 2 in the ηη ′ system.The significances for all other additional contributions are less than 5σ.Therefore, an η 1 state is included in the PWA.
In the final step, a baseline set of amplitudes is determined by adding the η 1 state, with its mass and width as free parameters, to the PDG-optimized set of amplitudes.The statistical significances of all resonances in the PDG-optimized set are then reevaluated in the presence of the η 1 state.Contributions from the f 0 (2100), h 1 (1595) γη ′ , ρ(1700) γη ′ , φ(2170) γη , f 2 (1810), and f 2 (2340) in the PDG-optimized set of amplitudes become insignificant (< 3σ) and are thus omitted from the baseline set of amplitudes, where the subscript labels the composition of the resonance.The statistical significance of the f 4 (2050) is reduced from 5.6σ to 4.6σ, but is still retained.By introducing the η 1 , the mass and width of the f 0 (2020) becomes more consistent with the average values in the PDG [14].In addition, a nonresonant contribution modeled by a 0 ++ ηη ′ system uniformly distributed in the phase space, is included with a significance of 15.7 σ.After this amplitude selection process, the baseline set of amplitudes includes eleven components.The isoscalar state with exotic quantum numbers J P C = 1 −+ , the η 1 , has a mass of (1855±9 stat ) MeV/c 2 and a width of (188±18 stat ) MeV with a statistical significance of 21.4σ.It is denoted as η 1 (1855).
The results of the PWA with the baseline set of amplitudes, including the masses and widths of the resonances, the product branching fractions B(J/ψ → γX)B(X → ηη ′ ) and B(J/ψ → η ′ X)B(X → γη), and the statistical significances, are summarized in Table III.The fit fractions for each component and their interference fractions are shown in Table IV.The measured masses and widths of the f 0 (2020) and f 2 (2010) are consistent with the PDG [14] average values.The measured mass of the f 0 (2330), which is unestablished in the PDG [14], is consistent with the results of Ref. [61], but our measured width is 79 MeV smaller (3.4σ).
All other resonances considered have statistical significance less than 3σ when added to the baseline set of amplitudes, as shown in Table V.To investigate additional possible contributions, resonances with different , and 1 −− γη (′) .)and with different masses and widths are added to the baseline set of amplitudes.No significant contributions from additional resonances with conventional quantum numbers are found.The most significant additional contribution (4.4σ) comes from an exotic 1 −+ component around 2.2 GeV.Changing the J P C assignment of the 0 ++ PHSP component in the baseline set of amplitudes to , and 1 −− γη (′) , results in a worse negative log-likelihood by at least 57.Furthermore, additional nonresonant contributions with all other J P C assignments are found to be insignificant.II.The mass (M ), width (Γ), PDG mass (MPDG), PDG width (ΓPDG), significance (Sig.) and the product branching fractions B(J/ψ → γX)B(X → ηη ′ ), B(J/ψ → η ′ X)B(X → γη) and B(J/ψ → ηX)B(X → γη ′ ) (B.F.) for each component in the PWA fit using the PDG-optimized set of amplitudes.The uncertainties are statistical.

Decay mode
Resonance M (MeV/c ground subtracted) and the PWA fit projections, respectively.Figure 2 (d) shows the cosθ η distribution, where θ η is the angle of the η momentum in the ηη ′ (Jacob and Wick) helicity frame (in which the ηη ′ system is at rest and the z-axis is defined by the momentum of the photon) [63].This angle carries information about the spin of the particle decaying to ηη ′ .The χ 2 /n bin value is displayed on each figure to demonstrate the goodness of fit, where n bin is the number of bins in each histogram, and χ 2 is defined as: where n i and v i are the number of events for the data and the fit projections with the baseline set of amplitudes in the ith bin of each figure, respectively.In comparison, the χ 2 /n bin values for the PDG-optimized set of amplitudes for the distributions of M (ηη ′ ), M (γη), M (γη ′ ), and cosθ η are 0.26, 0.43, 0.12 and 0.30 worse than baseline set of amplitudes, respectively.Figure 3 shows the Dalitz plots for the PWA fit projection from the baseline set of amplitudes, the selected data, and the background estimated from the η ′ sideband.Figure 2 and Figure 3 indicates that the data and the PWA fit result (baseline set of amplitudes) are in good agreements.Compared with the PDG-optimized set, the negative log-likelihood value of the baseline set is improved by 32 units and the number of free parameters is reduced by 16.
C. Further checks on the η1(1855) Various checks are performed to validate the existence of the η 1 (1855).The fits are carried out by assigning all other possible J P C to the η 1 (1855), and the negative loglikelihoods are worse by at least 235 units (21.8σ).To probe the significance of the BW phase motion, the BW parameterization of the η 1 (1855) in the baseline PWA is replaced with an amplitude whose magnitude matches that of a BW function but with constant phase (independent of s).This alternative fit has a negative log-likelihood 43 units (9.2σ) worse than the baseline fit, which indicates that a resonant structure is favorable.In the scenario η 1 (1855) is removed from the baseline set of amplitudes, the significance of η 1 (1855) with different masses and widths are evaluated.The changes of negative loglikelihood value are shown in Fig. 4. The result shows that a significant 1 −+ contribution is needed around 1.85 GeV/c 2 .
To visualize the agreement between the PWA fit results and data, angular moments as a function of M (ηη ′ ) can be calculated for data (with background subtracted) and the PWA model.For events within a given region of M (ηη ′ ), the cosθ η distribution can be expressed as an expansion in terms of Legendre polynomials.The coefficients, which are called the unnormalized moments of the expansion, characterize the spin of the contributing ηη ′ resonances.The moment for the TABLE III.The mass (M ), width (Γ), PDG mass (MPDG), PDG width (ΓPDG), significance (Sig.) and the product branching fractions B(J/ψ → γX)B(X → ηη ′ ) and B(J/ψ → η ′ X)B(X → γη) (B.F.) of each component in the PWA fit using the baseline set of amplitudes.
The first uncertainties are statistical and the second are systematic.

Decay mode
Resonance M (MeV/c   kth bin of M (ηη ′ ) is For data, N k is the number of observed events in the kth bin of M (ηη ′ ) and W i is a weight used to implement background subtraction.For the PWA model, N k is the number of events in a PHSP MC sample, which is generated with signal events distributed uniformly in phase space, and W i is the intensity for each event calculated in the PWA model.
Neglecting ηη ′ amplitudes with spin greater than 2, and ignoring the effects of symmetrization and the presence of resonance contributions in the γη and γη ′ subsystems, the moments are related to the spin-0 (S), spin-1 (P ) and spin-2 (D) amplitudes by [64,65]: where φ P and φ D are the phases of the P wave and D wave relative to the S wave. Figure 5 shows the moments computed for the data and the PWA model, using Eq. 19, where good data/PWA consistency can be seen.The need for the η 1 (1855) P-wave component is apparent in the Y 0 1 moment [Fig.5(b)].
Figure 6 shows a comparison of the data and the PWA projection of cosθ η in different M (ηη ′ ) regions ([1.5,1.7],[1.7,2.0] and [2.0,3.2]GeV/c 2 ).There is a clear asymmetry in the cosθ η distribution in the region [1.7,2.0]GeV/c 2 large-ly due to the η 1 (1855) signal, and the χ 2 /n bin of this region indicates a good agreement between data and the fitting results.
D. Discussion of the f0(1500) and f0(1710) The dominant contributions in the baseline PWA are from the scalar resonances.A significant signal for the f 0 (1500) is observed with a large product branching fraction (B(J/ψ → γf 0 (1500)) B(f 0 (1500) →ηη ′ ) = (1.81±0.11stat )×10 −5 ). 2. Background-subtracted data (black points) and the PWA fit projections (lines) for (a,b,c) the invariant mass distributions of (a) ηη ′ , (b) γη, and (c) γη ′ , and (d) the distribution of cosθη, where θη is the angle of the η momentum in the ηη ′ helicity coordinate system.The red lines are the total fit projections from the baseline PWA.The blue lines are the total fit projections from a fit excluding the η1 component.The dashed lines for the 1 −+ , 0 ++ , 2 ++ , 4 ++ and 1 +− contributions are the coherent sums of amplitudes for each J P C .Since the mass of the f 0 (1500) is close to the ηη ′ mass threshold and the f 0 (1500) has other decay modes, we parameterize the f 0 (1500) with a Flatté-like form with its mass and width as free parameters.The Flatté-like propagator is where the Γ in the first term of Γ(s) is an effective parameter corresponds to the decay mode f 0 (1500) → ηη ′ , l is orbital angular momentum of ηη ′ system, g is B(f 0 (1500) → ηη ′ ), which is estimated to be 0.02 from the PDG [14], the second term corresponds to all other decay modes of the f 0 (1500), and Γ 0 is a constant which represents the total width of the f 0 (1500) listed in the PDG [14].ρ(s) is the momentum magnitude of η (′) in the ηη ′ rest frame : The impact of using the Flatté-like parametrization for the f 0 (1500) is assigned as a systematic uncertainty, which is discussed further in Sec.V.
If the f 0 (1810) is replaced by the f 0 (1710) with mass and width taken from the PDG, the product branching fraction B(J/ψ→γf 0 (1710))B(f 0 (1710)→ηη ′ ) becomes (7.16±2.21stat )×10 −7 and the negative log-likelihood is worse by 29 units (7.3σ).There is therefore no significant evidence for J/ψ→γf 0 (1710) → γηη ′ .To determine the upper limits on B(J/ψ→γf 0 (1710))B(f 0 (1710)→ηη ′ ) for different scenarios, the same approach as that in Ref. [66] is used.For each alternative fit, the upper limit is determined at 90% of the integral of a Gaussian distribution with mean and width equal to the fitting yield and the statistical uncertainty.The maximum value is taken as the upper limit.

V. SYSTEMATIC UNCERTAINTIES
The sources of systematic uncertainty are divided into two categories.The first category concerns the systematic uncertainties related to event selection, which are applicable to mea-surements of the branching fractions.These sources of systematic uncertainties are described below.The second category of systematic uncertainties concerns the PWA and will be treated later.
(i) Pion tracking.The MDC tracking efficiency of charged pions is investigated using a clean control sample of J/ψ → ppπ + π − [67].The difference in tracking efficiency between data and MC simulation is 1% for each charged pion.
(ii) Photon detection efficiency.The photon detection efficiency is studied with a clean sample of J/ψ →ρ 0 π 0 [68].The result shows that the data-MC efficiency difference is 1% per photon.
(iii) Kinematic fit.To investigate the systematic uncertainty associate with the kinematic fit, the track helix parameter correction method [69] is used.The difference in the detection efficiency between using and not using the helix correction is taken as the systematic uncertainty.
(iv) η ′ mass resolution.The difference in the mass resolution between data and MC simulation leads to uncertainties related to the η ′ mass window requirement.This is investigated by smearing the MC simulation to improve the consistency between data and MC simulation.The difference of the detection efficiency before and after smearing is assigned as the systematic uncertainty for the η ′ mass window requirement.
For the two η ′ decay modes, the systematic uncertainties from pion tracking, four photon detection, number of J/ψ events, B(η → γγ) are common systematic uncertainties and the other systematic uncertainties are independent systematic uncertainties.The combination of common and independent systematic uncertainties for these two decay modes are calculated with the weighted least squares method [70], and the total systematic uncertainty is determined to be 4.8%.A summary of all systematic uncertainties related to event selection is shown in Table VI.Systematic uncertainties from the PWA impact the branching fractions and resonance parameters.These are studied below, and the statistical significance of the η 1 (1855) is recalculated in every variation.
Uncertainty from the BW parametrization is estimated by the changes in the fit result caused by replacing the constant width Γ 0 of the BW for the threshold state f 0 (1500) with a mass dependent width as described in Sec.IV D. The statistical significance of the η 1 (1855) in this case is 21.8σ.(ii) Uncertainty from resonance parameters.In the baseline fit, the resonance parameters of the f 0 (1500), f 2 (1565), f 4 (2050), h 1 (1415) γη , and h 1 (1595) γη are fixed to PDG [14] average values, and the resonance parameters of f 0 (1810) are fixed to the previous measurment [62].An alternative fit is performed where resonance parameters are allowed to vary within one standard deviation of the PDG values [14] and Ref. [62], and the changes in the results are taken as systematic uncertainties.The statistical significance of the η 1 (1855) in this case is 20.6σ.
(iii) Background uncertainty.To estimate the uncertainty due to the background estimation, alternative fits are performed using different η ′ sideband regions and different background normalization factors.In detail, the background normalization factors are varied by one standard deviation, which is determined from the fit results of Fig. 1(a) and Fig.For each alternative fit performed to estimate the systematic uncertainties in the PWA fit procedure, the changes of the measurements are taken as the one-sided systematic uncertainties.For each measurement, the individual uncertainties are assumed to be independent and are added in quadrature to obtain the total systematic uncertainty on the negative and positive sides, respectively.The sources of systematic uncertainties affecting the measurements of masses and widths of the f 0 (2020), f 0 (2330), η 1 (1855), f 2 (2010) and their contributions are summarized in Table VII.The relative systematic uncertainties relevant to the branching fraction measurements are summarized in Table VIII.
To include the systematic uncertainties in the upper limit B(f 0 (1710)→ηη ′ )/B(f 0 (1710)→ππ), the additive systematic uncertainties, i.e., the systematic uncertainties on the upper limit associated with the PWA, are considered by performing alternative fits and taking the maximum value as the upper limit.The multiplicative systematic uncertainties, i.e. the other systematic uncertainties, are taken into account by dividing by the factor (1 − δ combined ), where δ combined is the systematic uncertainties associated with the event selection and the uncertainty of branching fractions and the uncertainty of J/ψ→γf 0 (1710)→γππ [14].The upper limit on B(f 0 (1710)→ηη ′ )/B(f 0 (1710)→ππ) at 90% C.L. is determined to be 2.87×10 −3 .
FIG.2.Background-subtracted data (black points) and the PWA fit projections (lines) for (a,b,c) the invariant mass distributions of (a) ηη ′ , (b) γη, and (c) γη ′ , and (d) the distribution of cosθη, where θη is the angle of the η momentum in the ηη ′ helicity coordinate system.The red lines are the total fit projections from the baseline PWA.The blue lines are the total fit projections from a fit excluding the η1 component.The dashed lines for the 1 −+ , 0 ++ , 2 ++ , 4 ++ and 1 +− contributions are the coherent sums of amplitudes for each J P C .

5 FIG. 4 .
FIG.4.The change in negative log-likelihood values for a range of η1 resonance parameters in the baseline set of amplitudes.The red dashed line indicates that the statistical significances of the hypotheses under the red dashed line are higher than 5σ.

FIG. 5 .
FIG. 5.The distributions of the unnormalized moments Y 0L (L = 0, 1, 2, 3 and 4) for J/ψ → γηη ′ as functions of the ηη ′ mass.Black dots with error bars represent the background-subtracted data weighted with angular moments; the red solid lines represent the baseline fit projections; and the blue dotted lines represent the projections from a fit excluding the η1 component.

FIG. 6 .
FIG.6.Background-subtracted data (black points) and the PWA fit projections (lines) for the cosθη distribution when the ηη ′ mass is restricted to the regions: (a) [1.5,1.7],(b) [1.7,2.0], and (c) [2.0,3.2]GeV/c 2 .The red lines are the total fit projections from the baseline PWA.The blue lines are the total fit projections from a fit excluding the η1 component.

d
Also at Goethe University Frankfurt, 60323 Frankfurt am Main, Germany e Also at Key Laboratory for Particle Physics, Astrophysics and Cosmology, Ministry of Education; Shanghai Key Laboratory for Particle Physics and Cosmology; Institute of Nuclear and Particle Physics, Shanghai 200240, People's Republic of China f Also at Key Laboratory of Nuclear Physics and Ion-beam Application (MOE) and Institute of Modern Physics, Fudan University, Shanghai 200443, People's Republic of China g Also at State Key Laboratory of Nuclear Physics and Technology, Peking University, Beijing 100871, People's Republic of China h Also at School of Physics and Electronics, Hunan University, Changsha 410082, China i Also at Guangdong Provincial Key Laboratory of Nuclear Science, Institute of Quantum Matter, South China Normal University, Guangzhou 510006, China j Also at Frontiers Science Center for Rare Isotopes, Lanzhou University, Lanzhou 730000, People's Republic of China k Also at Lanzhou Center for Theoretical Physics, Lanzhou University, Lanzhou 730000, People's Republic of China l Also at the Department of Mathematical Sciences, IBA, Karachi , Pakistan (Dated: March 8, 2023) I. INTRODUCTION

TABLE I .
The set of all intermediate resonances considered when the PDG-optimized set of amplitudes is developed.States with quantum numbers J P C = 0 ++ , 2 ++ , and 4 ++ in the ηη ′ spectrum and states with quantum numbers J P C = 1 +− and 1 −− in the γη (′) spectra are considered.

TABLE IV .
The fit fractions for each component and the interference fractions between two components(%) in the PWA fit with the baseline set of amplitudes.The uncertainties are statistical only.

TABLE V .
Additional resonances considered, their J P C , the change in negative log-likelihood (∆lnL) when each is added to the baseline set of amplitudes, the change in the number of free parameters (∆dof), and the resulting statistical significance (Sig.).