Observation of a near-threshold enhancement in the $\Lambda\bar{\Lambda}$ mass spectrum from $e^+e^-\to\phi\Lambda\bar{\Lambda}$ at $\sqrt{s}$ from 3.51 to 4.60 GeV

The process $e^+ e^- \to \phi \Lambda \bar{\Lambda}$ is studied using data samples collected with the BESIII detector at the BEPCII collider at center-of-mass energies $\sqrt{s}$ ranging from $3.51$ to $4.60~{\rm GeV}$ . An intermediate resonance structure is observed near the threshold of $\Lambda \bar{\Lambda}$. It has a mass of $(2262 \pm 4 \pm 28)~{\rm{MeV}}/c^{2}$ and a width of $(72 \pm 5 \pm 43)~\rm{MeV}$, where the quoted uncertainties are statistical and systematic, respectively. The $J^{PC}$ quantum numbers of $0^{-+}$ and $0^{++}$ are rejected, while other $J^{PC}$ hypotheses are possible, according to the helicity angle study. The energy-dependent cross section of the $e^+ e^- \to \phi \Lambda \bar{\Lambda}$ process is measured for the first time in this energy region, and contributions from excited $\psi$ states and vector charmonium-like $Y$-states are investigated.


I. INTRODUCTION
Early in this century, a number of exotic states have been discovered [1] in final states with a quarkonium and one or two light hadrons, or with heavy-flavor mesons. Among these states, there are vector states with J P C = 1 −− which are usually called Y states, such as the Y (4260) [2], Y (4360) [3,4], and Y (4660) [3]. The Y (4260) state is observed for the first time by the BABAR experiment with a mass of (4259 ± 8 +2 −6 ) MeV/c 2 using the initial state radiation (ISR) events e + e − → γ ISR π + π − J/ψ [2]. The observation was latter confirmed by the CLEO [5] and Belle experiments [6]. In 2017, a dedicated analysis performed by the BESIII experiment revealed that the so-called Y (4260) state is not simply one Breit-Wigner (BW) resonance and can be a combination of two states [7]. The first one has a lower mass and a much narrower width than the Y (4260), but is consistent with the Y (4220) state observed in e + e − → π + π − h c events [8,9], and the second one at around 4.32 GeV/c 2 is observed for the first time with a significance greater than 7.6σ. The lower-mass resonance was also observed in e + e − → ωχ c0 [10] and πDD * + c.c. events [11].
Until now, the Y (4260) and other vector charmonium-like states were only reported in final states containing a cc pair: either charmonium states or charmed mesons. Several analyses have been performed by the BESIII collaboration to search for light hadron decays of these states, for example, Y (4260) → π 0 (η)pp [12], K + K − π 0 [13], ΞΞ [14], etc. Although there is no significant contribution from the vector charmonium or charmonium-like states identified, the cross section lineshapes of these processes do suggest contributions from amplitudes beyond simple continuum production. In Ref. [15], the Y (4260) state is interpreted as a diquark-antidiquark state ([cs] [cs]). This interpretation implies that the Y (4260) state decays easily into final states containing a pair of ss. One of the dominant contributions to e + e − → π + π − J/ψ comes from Y (4260) → f 0 (980)J/ψ decays [16], and the f 0 (980) meson is known to have a large ss component. If the cc quarks in Y (4260) annihilate while the ss pair survives in the final state, we expect Y (4260) decays into strange mesons or baryons, such as φΛΛ. Such a signal will manifest itself as a lineshape distortion due to the interference between the amplitudes of the Y (4260) decay and the e + e − → φΛΛ continuum production.
The η(2225) and φ(2170) states [17] are interpreted as loosely bound states of ΛΛ in Ref. [18]. This suggests that the η(2225) couples to ΛΛ strongly above the threshold and it can be produced in e + e − → φΛΛ processes. Together with the strong enhancement of the e + e − → ΛΛ production cross section close to threshold [19], this can help establish if there is any connection between these two hadron molecule candidates. On the other hand, near-threshold enhancements are observed in several processes involving baryon-anti-baryon pairs such as J/ψ → γpp [20], B → Kpp [21], and B 0 → KΛΛ [22]. There are a few interpretations for these phenomena, including states near threshold as found in a model by Nambu and Jona-Lasinio [23], J P C = 0 ±+ isoscalar states coupled to a pair of gluons [24], and low-mass enhancements favored by the fragmentation process [24]. The isoscalar C = + ΛΛ threshold enhancement can be searched for in e + e − → φΛΛ and its spin-parity can be determined by studying the angular distribution if the data sample is large enough.
In this paper, we report the first observation of the process e + e − → φΛΛ analysing data samples taken at center-of-mass (CM) energies √ s ranging from 3.51 to 4.60 GeV. The vector charmonium-(like) states are studied based on the energy-dependent cross sections, and an intermediate state in the ΛΛ system is investigated to extract information on light mesons.

II. DETECTOR AND DATA SAMPLES
The BESIII detector is a magnetic spectrometer [25] located at the Beijing Electron Positron Collider (BEPCII) [26]. The cylindrical core of the BESIII detector 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 magnetic field. The solenoid is supported by an octagonal flux-return yoke with resistive plate counter muon identifier modules interleaved with steel. The acceptance of charged particles and photons is 93% over 4π solid angle. The charged-particle momentum resolution at 1 GeV/c is 0.5%, and the dE/dx resolution is 6% for the electrons from Bhabha scattering. The EMC measures photon energies with a resolution of 2.5% (5%) at 1 GeV in the barrel (end cap) region. The time resolution of the TOF barrel part is 68 ps, while that of the end cap part is 110 ps. The end cap TOF system is upgraded in 2015 with multi-gap resistive plate chamber technology, providing a time resolution of 60 ps [27].
The experimental data used in this analysis were taken at the CM energies ranging from 3.51 to 4.60 GeV as shown in Table I. Simulated samples produced with the GEANT4-based [28] Monte Carlo (MC) package which includes the geometric description of the BESIII detector and the detector response, are used to determine the detection efficiency and to estimate the background contributions. The simulation includes the beam energy spread and ISR in the e + e − annihilations modeled with the generator KKMC [29]. Inclusive MC simulation samples generated at √ s = 4.178 GeV are used to analyze the possible background contributions. In total, these samples are 40 times larger than the data sample. They consist of open charm production processes, ISR production of vector charmonium(-like) states, and continuum processes (e + e − → qq , q = u, d, s). The open charm production processes are generated using CONEXC, and the ISR production is incorporated in KKMC [29]. The known decay states are modeled with EVTGEN [30] using branching fractions taken from the Particle Data Group (PDG) [31], and the remaining unknown decays from the charmonium states with LUNDCHARM [32]. The final state radiations (FSR) from charged final state particles is incorporated with the PHOTOS package [33]. It should be pointed out that only e + e − → K + K − ΛΛ is simulated in the inclusive MC samples.
The signal MC samples of e + e − → φΛΛ are generated with EVTGEN [30] along with KKMC [29] to handle the e + e − annihilations and ISR production. The signal events are generated with three-body phase space (PHSP) model where the φΛΛ distributed uniformly in the phase space. The data samples used in this analysis have been collected by BESIII at 28 CM energies between 3.51 and 4.60 GeV, as listed in Table I, along with the CM energy and corresponding integrated luminosity. The total integrated luminosity is 19.5 fb −1 .

III. EVENT SELECTION
The selection of charged tracks is based on the following criteria. For each charged track, the polar angle in the MDC must satisfy | cos θ| < 0.93, and the point of closest approach to the e + e − interaction point (IP) must be within ±20 cm in the beam direction and within 10 cm in the plane perpendicular to the beam direction. The particle identification (PID) of kaons, pions, and protons is based on the dE/dx and TOF information. Assumption of a given particle identification is based on the largest of the all PID hypotheses probabilities. The φ meson is reconstructed using candidate K + K − pairs. OneΛ (Λ) baryon is assumed to be missing in order to improve the reconstruction efficiency. Thus we require that there should be at least one proton and one pion with opposite charge, and one K + K − pair in the final state.
Since the Λ baryon has a relatively long lifetime, it travels a certain distance before it decays. A vertex fit is applied to its decay products pπ − (pπ + ) to ensure that their tracks are pointing back to the same vertex. The Λ (Λ) baryon is reconstructed combining the pπ − (pπ + ) final state passed the vertex fit. Then, to verify that the selected K + K − Λ(Λ) candidates originate from the IP, another vertex fit is performed. Only events with a good quality vertex fit are retained. The flight distance between the IP and the Λ decay vertex is required to be greater than two times of its resolution. The momenta corrected by the vertex fit are used for kinematic fit.
To improve the track momentum resolution and to reduce the background, a kinematic fit is applied to the K + K − Λ(Λ) candidates constraining the missing mass to the nominal mass of Λ. The fraction of events containing more than one Λ baryon is about 31%. For events with multiple candidates, we choose the combination with the smallest χ 2 combining the two vertex fits, and the kinematic fit. The distributions of the combined χ 2 versus the invariant mass of pπ are shown in Fig. 1. To make a better comparison, the one dimensional χ 2 distribution compared between data and MC is shown in Fig. 1 (right). The possible difference are considered as systematic uncertainty coming from the kinematic fit. The sum of the χ 2 values of the vertex and kinematic fits is required to be less than 30. The invariant mass of selected pπ final state should be within the interval [1.112, 1.120] GeV/c 2 , which covers about ±3σ of the Λ signal region.  Here, L int is the integrated luminosity, N sig is the number of signal events from the fit to M (K + K − ) distributions with statistical uncertainty only, ε is the efficiency, (1 + δ) is the radiative correction, and σ dress (pb −1 ) is the cross section quoted with a statistical and systematic uncertainty, respectively. Studies of the inclusive MC simulation indicate that the main background contribution comes from the process e + e − → (γ)K + K − ΛΛ, which does not peak around the φ signal area. It should be pointed out that at √ s = 4.6 GeV, the CM energy is above the threshold of e + e − → Λ + cΛ − c , and there is a background contribution from Λ ± c → ΛK ± decays. However, according to the Born cross section reported in Ref. [34], the contribution of this background in the whole fitting range (the invariant mass of K + K − M(K + K − ) ∈ [0.98, 1.20] GeV/c 2 ) is estimated to be only 8.8±0.1 events. Other sources of background are considered are found to be negligible.
To obtain the signal yields, an unbinned maximum likelihood fit is performed to the invariant mass spectrum of the K + K − pair for each CM energy point. The signal distribution is described by a MC-simulated shape, and the background shape is described by an inverted ARGUS [35] function whose threshold is fixed to 2m K ± , where m K ± is the nominal kaon mass [31]. The fit result for √ s = 4.178 GeV is shown in Fig. 2 as an example, and the numbers of signal events (N sig ) at 28 energy points are listed in Table I. We perform a study to investigate possible intermediate structures to better estimate the reconstruction efficiency. The Dalitz plot distribution of the φΛΛ candidates at √ s = 4.178 GeV is shown in Fig. 3, after requiring that M(K + K − ) ∈ [1.01, 1.03] GeV/c 2 . It is clear that most of the events in the data are deposited near the ΛΛ threshold, which is different from the PHSP MC sample generated with a uniform distribution. Signal MC samples are generated at 28 energy points to study the reconstruction efficiency and resolutions. The efficiency and resolution curves are shown in Fig. 4 for MC samples at √ s = 4.178 GeV. We can see that the reconstruction efficiency is quite smooth near the threshold and the resolution is relatively small.   The invariant mass distribution of the ΛΛ candidates for the full data sample is shown in Fig. 9. There are serval dynamics to generate such an enhancement, including final state interaction (FSI), a tail of a lower mass resonance, and so on [20,24,31]. To describe the lineshape of this enhancement, an extended unbinned maximum likelihood fit is performed to all the data samples simultaneously. We firstly perform a fit using a Breit-Wigner function (BW) to describe the signal. Three components are considered in the fit: a near-threshold enhancement, a component distributed uniformly in PHSP, and a non-φ background component. The interference between the resonant signal and non-resonant signal is ignored here. The following formula is used to describe the lineshape of the enhancement [36]: where ε(M ΛΛ ) is the mass-dependent efficiency obtained from MC simulation. Here the MC sample is generated with the nonuniform angular distributions measured in data (to be described latter).M ΛΛ is the invariant mass of ΛΛ system, is the momentum of ΛΛ system in the e + e − rest frame, where m φ is the nominal mass of φ [31], is the momentum of the Λ baryon in ΛΛ system rest frame, l is the orbital angular momentum between φ and ΛΛ system, L d is the orbital angular momentum between Λ andΛ, f L is the Blatt-Weisskopf barrier factor, with f 2 0 (z) = 1, f 2 1 (z) = 1/(1 + z), and f 2 2 (z) = 1/(9 + 3z + z 2 ). The relativistic Breit-Wigner function with a mass-dependent width used here is defined as where 2 , m and Γ 0 are the mass and width of the BW function, respectively, and q 0 is equal to q * for M(ΛΛ) = m. In the fit, the mass and width are shared parameters between all the data samples, and are left free, as well as the signal yields. The orbital angular momentum between φ and the ΛΛ system is l = 0, and the orbital angular momentum between the Λ andΛ baryons is L d = 1, assuming this is a 1 ++ or 2 ++ state. Please note that even we use a BW function here to describe the near-threhold enhancement, we are not suggesting that this enhancement is a resonant-(like) structure. The resolution effect is ignored here because it is relatively small compared with such a broader distribution.
The shape for PHSP signal is obtained from MC simulation. The shape of the nonφ background is obtained from the φ sideband region (M(K + K − ) ∈ [0.99, 1.005] or [1.075, 1.090] GeV/c 2 ), and is parameterized with a Landau function. The number of background events is extrapolated from the sidebands to the φ signal region using the inverted ARGUS background function. The fit result using all data samples [37] is show in Fig. 9 (left). We also zoom in on the lower mass side to have a closer look on the rise of the enhancement, as shown in Fig. 9 (right). The mass and width of the BW formula are fitted as (2262 ± 4) MeV/c 2 and (72 ± 5) MeV, respectively.
Alternatively, we perform a fit to estimate the rise rate near the threshold with the formula: where P 3 is a third-order polynomial whose parameters are free, p 0 is a free parameter, and ∆M ΛΛ ≡ M ΛΛ − 2m Λ . The fit result is shown in Fig. 9 (right), with p 0 = 33 ± 11 MeV/c 2 . Compared with the lineshapes of PHSP events weighted by angular distribution and cross section from each energy point, the rising rate in data is much faster.
To further understand the nature of this enhancement, the helicity angles of the φ and Λ candidates are studied. The helicity angle is defined as the angle between the momentum of the φ or Λ in its parent's rest frame and the momentum of φ or Λ's parent in its grandparent's rest frame. The helicity angular distributions for events in the φ signal region after efficiency correction are shown in Fig. 6, combining all data samples. The unbinned maximum likelihood fit is performed simultaneously to the angular distributions, considering the same components as the ones contributing to the fit of the M(ΛΛ) distribution. Fraction of each component is fixed to that obtained from the ΛΛ mass spectra fit. Possible interference is not considered in the fit. The shapes of resonant signal are described with the formula constructed according to Ref. [38]. The details of the formula we used are provided in the supplementary material. The shapes of the background and PHSP signal distributions are assumed to be flat. The number of events in each component is fixed to the fit result of M(ΛΛ). The fit favors the hypotheses of J P C = 1 ++ , 2 ++ , or 2 −+ , where the results lead to the similar fit quality. The hypothesis of this enhancement having spin zero is rejected with significance greater than 7σ compared with other hypotheses. The fit result with different J P C hypotheses are shown in Fig. 6.
Data driven reconstruction efficiencies are obtained by re-weighting the signal MC samples in the generator level. The contribution of the near-ΛΛ-threshold enhancement and non-uniform angular distributions measured in data are considered. The energy-dependent reconstruction efficiencies for the PHSP and the re-weighted models are shown in Fig. 7. The fine structures observed in the efficiency curve is due to the deformation in the cross section line shape which is considered in MC generation to obtain the correct ISR factor and efficiency.

C. Cross section measurement
The cross section at a certain CM energy is calculated as where N sig is the number of φΛΛ signal events obtained from the fit to the M(K + K − ) distribution, L int is the integrated luminosity, ε is a weighted value of the efficiencies from the process e + e − → φΛΛ where M(ΛΛ) following the lineshape of the near-threshold enhancement as well as the angular distributions, and the process e + e − → φΛΛ uniformly distributed in phase space, B is the product of the branching fraction of the intermediate decays φ → K + K − and Λ → pπ, which are taken from Ref. [31], (1 + δ) is the ISR correction factor.
To obtain the proper ISR correction factor, an iterative procedure is used. First, a series of signal MC samples are generated for all energy points with a constant cross section using KKMC. The cross sections are calculated based on the reconstruction efficiencies and ISR correction factors obtained from the signal MC simulation. We use the Lowess [39] method to smoothen the lineshape of the measured cross sections, then use the method introduced in Ref. [40] to get the ISR correction factors and efficiencies with the new lineshape. A new series of cross sections could be obtained, and after several iterations, the cross section results become stable.
However, when the iteration is performed, the cross section results at each energy point are correlated. To take the correlation into consideration, we use pseudoexperiments. First, a large pseudodata sample is generated by sampling a Gaussian distribution, the mean value of which is the nominal cross section result, and its width is the statistical error from the fit to the data. Then, the iteration described in the previous paragraph with this new lineshape is performed. The resulting cross section distributions at each energy point are fitted with Gaussian functions. Their mean and width values are taken as the final results for the cross sections and their corresponding uncertainties, respectively. The final results are shown in Fig. 8, and listed in Table I, including the statistical and systematic uncertainties.

A. Uncertainties on Cross Sections
The systematic uncertainties include contributions from luminosity, tracking and PID efficiencies of the kaons, Λ reconstruction, radiative correction factor associated with the efficiency, background, and branching fractions of the intermediate states.
The integrated luminosity is measured using Bhabha scattering events, with uncertainty smaller than 1.0% [41]. The uncertainty related to the tracking efficiency of kaons is estimated to be 1.0%, and the uncertainty arising from the kaon PID efficiency is determined to be 1.0% using a control sample e + e − → K + K − π + π − . With the control sample, the tracking or PID requirement efficiency is separately measured in the MC simulation sample and in the data sample. The difference between the efficiencies from MC simulation and data samples is taken as the systematic uncertainty.
The systematic uncertainty due to the Λ reconstruction efficiency including the tracking and PID of its decay products pπ, as well as the decay length requirement, is studied with the control sample Λ c → Λ + X decays. The resulting systematic uncertainty is 1.1% [42].
The systematic uncertainty due to the Λ(Λ) and φ mass window selection criteria accounts for the mass resolution discrepancy between the MC simulation and experimental data. The φ and Λ mass distributions from signal MC sample and data are fitted with double-Gaussian functions and compared with each other. The difference between the fit results is negligible.
For the uncertainty due to the ISR correction factor, we change the lineshape with a power law function 1/s n . The difference between the nominal result and the alternative parameterization is taken as the systematic uncertainty. To estimate the systematic uncertainty related to the background model, we vary the fit range of the M(K + K − ) distribution. We also use a secondorder polynomial as an alternative background model. The largest value among all variations is taken as systematic uncertainty for this source. Due to limited sample sizes at most CM energy points, the uncertainty from the data sample collected at √ s = 4.178 GeV is used for all the data sets.
The uncertainty of the kinematic fit is estimated by comparing the reconstruction efficiency before and after the helix parameter correction using the method described in Ref. [43]. It should be pointed out that the reconstruction efficiency after the helix parameter correction is used as the nominal result. The uncertainties of the branching fractions are taken from the PDG [31].
The summary of the systematic uncertainties at √ s = 4.178 GeV is presented in Table II.

B. Uncertainties on M (ΛΛ) Lineshape
The systematic uncertainties for the mass and width of the BW formed lineshape include those from the mass calibration, efficiency curve, signal parameterization, and background estimation.
To calibrate the mass component, a maximum likelihood fit of the K + K − invariant mass distributions is performed for all the data samples. The difference between the fitted mass and the known mass of the φ meson [31] is 0.4 MeV/c 2 . According to the conservation of energy and momentum, the difference on the ΛΛ side should be 0.3 MeV/c 2 . This value is taken as the systematic uncertainty.
To evaluate the systematic uncertainty from the efficiency curve estimation, we use unweighted PHSP MC sample instead of the nominal one to extract the efficiency curve. The changes on the mass and width, 0.6 MeV/c 2 and 8.0 MeV, respectively, are taken as the systematic uncertainties.
To account for the systematic uncertainty from the signal model, we change the parameterization form from a Landau to a BW function. The differences between the two parameterizations, 22.9 MeV/c 2 and 13.5 MeV, are taken as the systematic uncertainties on the mass and width, respectively. Another source of uncertainty in the signal parameterization is the quantum number assumption. We change the assignment of l/L d = 0/1 to 1/2. The differences on the mass and width, 1.7 MeV/c 2 and 36.0 MeV, are taken as the systematic uncertainties. For the background estimation, we get different yields by varying the fit range of the K + K − invariant mass distributions and repeating the fit. The differences on the mass and width, 3.1 MeV/c 2 and 1.2 MeV, respectively, are taken as the systematic uncertainties. Replacing the background parameterization with a BW function leads to changes in the measurement of the mass and width, 15.9 MeV/c 2 and 18.0 MeV, which are are taken as the systematic uncertainties from the background model. Table III summarizes these three sources of uncertainties. The sum of all the above uncertainties in quadrature, 28.1 MeV/c 2 and 43.2 MeV for the mass and width, respectively, are taken as the total uncertainties. The systematic uncertainty for the reversed exponential parameter come from similar sources. We broaden or narrow the fit range by 0.02 GeV, and the change of the exponential parameter 4.8 is taken as systematic uncertainty. Efficiency curve is flat near the threshold, which will not affect the exponential parameter. The mass calibration will not change the lineshape so this source is also ignored. We vary the background estimation and background shape, the largest changes of the exponential parameter is 0.4. The uncertainties are also summarized in Table III. The sum of all the above uncertainties in quadrature, 4.8, is taken as the total systematic uncertainty.

VI. SUMMARY AND DISCUSSION
In summary, we observe the e + e − → φΛΛ process for the first time with data samples at CM energies ranging from 3.51 to 4.60 GeV. The energy-dependent cross sections of e + e − → φΛΛ are measured. Due to the limited sample sizes, we cannot resolve the composition of the resonance structure, and the lineshape might not be simply described with a continuum process parameterized as 1/s n (n = 3.3 ± 0.3).
Moreover, a near-threshold enhancement is observed on ΛΛ with a significance greater than 25σ compared with the pure phase space distribution. By fitting the lineshape with a BW function, we obtain the mass and width as (2262±4±28) MeV/c 2 and (72±5±43) MeV, repectively, where the first uncertainties are statistical and the second ones systematic. By fitting the linehshape with an reversed exponential function, we obtain the rising rate (exponential parameter) as 33 ± 11 ± 6 MeV/c 2 .
According to the helicity angle study, the quantum numbers of ΛΛ system J P C = 0 ++/−+ is rejected with a significance of 7σ. The interpretation of the ΛΛ system originating from a decay η(2225) → ΛΛ is rejected. The J P C quantum numbers could be 2 ++ , 2 −+ , or 1 ++ , but they cannot be distinguished because of the limited data sample sizes. Another interpretation of a lower mass resonance is that this could be a f 2 (2300) meson. However, according to previous measurements using the decay modes J/ψ → γφφ and J/ψ → γΛΛ [44,45], the f 2 (2300) meson is more likely to decay into a φφ final state rather than to ΛΛ. The cross sections of the e + e − → φφφ process are measured [46] at BESIII with similar cross sections to those of the e + e − → φΛΛ process, but no structure around 2.23 GeV/c 2 is observed in the φφ mass spectrum. Therefore the interpretation of f 2 (2300) → ΛΛ is also rejected.
This enhancement does not match any known resonance [31] seen before, and could be the same thing observed in B → KΛΛ decays [22]. If so, the theoretical explanation, an isoscalar state with J P C = 0 ±+ coupled to a pair of gluons [24], could be discarded since the hypothesis implying quantum numbers 0 −+ is rejected based on the angular distribution study. Also, the author of Ref. [24] implies that the observed threshold enhancements in low-mass baryon-antibaryon systems might not be limited to the ordinary quantum numbers of the qq system. Further studies of the ΛΛ system would be helpful to understand the nature of this threshold enhancement. For example, a search for a threshold enhancement in the e + e − → ηΛΛ process could provide a crucial test because the states produced in this mode have the exact opposite C parity to the state in our analysis.

VII. ACKNOWLEDGEMENTS
The BESIII collaboration thanks the staff of BEPCII and the IHEP computing center for their strong support. This work is supported in part by National

Appendix B: Helicity Angular Distribution Evalution
The angular distribution parameterization is constructed according to Ref. [38]. The differential cross section where M = ±1 is the magnetic moment of e + e − system, µ, ν are the helicity eigenstates of φ, X(2260), µ ′ , ν ′ are the helicity eigenstates of Λ andΛ, and where a C 1 is constant parameter which we ignore, F i µ,nu and F j µ ′ ,ν ′ are the amplitudes where i = 1 is the spin of e + e − system, and j is the spin of ΛΛ system, and D i M,µ−ν is the Wigner D-matrix. After substituting the different hypothesis for J P C in the Wigner D-matrix, the formula is simplified and integrated over the the azimuthal angle, we obtain the angular distribution of the final state for different hypotheses.
For the hypotheses of 0 −+ and 0 ++ cases, the Wigner D-matrix D 0 M,µ−ν is a constant, which means the helicity angle distribution of Λ(Λ) should be uniform. The fitting results with J P C = 0 −+ hypothesis is shown in Fig. 10 as an example. Clearly the fitting result on Λ(Λ) helicity angular distribution is bad. For J P C = 1 ++ , we have where A = F 01 + F 10 + 2F 11 (B4) The fitting results with J P C = 1 ++ hypothesis is shown in Fig. 11.