First measurement of $e^+ e^- \to pK^{0}_{S}\bar{n}K^{-} + c.c.$ above open charm threshold

The process $e^+e^-\rightarrow pK^{0}_{S}\bar{n}K^{-} + c.c.$ and its intermediate processes are studied for the first time, using data samples collected with the BESIII detector at BEPCII at center-of-mass energies of 3.773, 4.008, 4.226, 4.258, 4.358, 4.416, and 4.600 GeV, with a total integrated luminosity of 7.4 fb$^{-1}$. The Born cross section of $e^+e^- \to p K^{0}_S\bar{n}K^- + c.c.$ is measured at each center-of-mass energy, but no significant resonant structure in the measured cross-section line shape between 3.773 and 4.600 GeV is observed. No evident structure is detected in the $pK^-$, $nK^{0}_S$, $pK^0_{S}$, $nK^+$, $p\bar{n}$, or $K^{0}_S K^-$ invariant mass distributions except for $\Lambda(1520)$. The Born cross sections of $e^+e^-\rightarrow\Lambda(1520)\bar{n}K^{0}_{S} + c.c.$ and $e^+e^-\rightarrow \Lambda(1520)\bar{p}K^{+} + c.c.$ are measured, and the 90\% confidence level upper limits on the Born cross sections of $e^+e^-\rightarrow\Lambda(1520)\bar{\Lambda}(1520)$ are determined at the seven center-of-mass energies.

The process e + e − → pK 0 Sn K − + c.c. and its intermediate processes are studied for the first time, using data samples collected with the BESIII detector at BEPCII at center-of-mass energies of 3.773, 4.008, 4.226, 4.258, 4.358, 4.416, and 4.600 GeV, with a total integrated luminosity of 7.4 fb −1 . The Born cross section of e + e − → pK 0 Sn K − + c.c. is measured at each center-of-mass energy, but no significant resonant structure in the measured cross-section line shape between 3.773 and 4.600 GeV is observed. No evident structure is detected in the pK − , nK 0 S , pK 0 S , nK + , pn, or K 0 S K − invariant mass distributions except for Λ(1520). The Born cross sections of e + e − → Λ(1520)nK 0 S + c.c. and e + e − → Λ(1520)pK + + c.c. are measured, and the 90% confidence level upper limits on the Born cross sections of e + e − → Λ(1520)Λ(1520) are determined at the seven center-of-mass energies. There is an evident difference in line shape and magnitude of the measured cross sections between e + e − → Λ(1520)(→ pK − )nK 0 S and e + e − → pK −Λ (1520)(→nK 0 S ).

I. INTRODUCTION
The experimental discovery of unexpected resonances has brought new opportunities to the study of quantum chromodynamics in the charmonium and bottomonium energy regions [1][2][3]. The state Y (4260) was discovered by the BaBar collaboration [4,5] in the initial state radiation (ISR) process e + e − → γ ISR π + π − J/ψ and confirmed by the CLEO [6] and Belle [7] collaborations in the same process. This state was further confirmed by BESIII [8] and again by Belle [9]. Located above the DD mass threshold, the Y (4260) with J P C = 1 −− anomalously couples to the hidden-charm final state ππJ/ψ [10]. The same phenomenon had been observed in other Y states, such as the Y (4360) and Y (4660) [1]. Just recently, BESIII first reported that two structures around 4.22 GeV and 4.39 GeV in e + e − line shape strongly couple to π + D 0 D * − [11]. These interesting but little-known phenomenons have prompted researchers to focus on this charmonium-like spectroscopy [1][2][3].
Light hadron decays of Y states have not been found above 4 GeV, nor have any such decays of charmonium resonances. The continued search for light hadron decays helps further the understanding of the nature of undefined states and charmonium resonances.
In this paper, we report the cross sections of e + e − → pK 0 Sn K − +c.c. and search for possible structures, such as Y states or higher charmonia , in the e + e − → pK 0 Sn K − + c.c. cross section line shape, using data samples with a total integrated luminosity of 7.4 fb −1 collected with the BESIII detector at center-of-mass (c.m.) energies between 3.773 and 4.6 GeV. All possible intermediate states in the pK − , nK 0 S , pK 0 S , nK + , pn, and K 0 S K − invariant mass spectra and charge conjugated modes, such as Λ * , Σ * , a 0 (980) + , and other excited or exotic states, are searched for in the e + e − → pK 0 Sn K − + c.c.
process. However, no significant structures, except for Λ(1520), are seen in any of the studied mass spectra. The Born cross sections of e + e − → Λ(1520)nK 0 S and Λ(1520)pK + + c.c. are measured. In the following analysis, charged conjugated modes are included unless otherwise indicated. The process e + e − → pK 0 Sn K − with all of the potential intermediate states included is denoted as the "pK 0 Sn K − mode" hereinafter. Similarly, the processes e + e − → Λ(1520)nK 0 S , Λ(1520)pK + and Λ(1520)Λ(1520) are denoted as "Λ(1520) modes".

II. EXPERIMENTAL DATA AND MONTE CARLO SIMULATION
The Beijing Electron Positron Collider II (BEPCII) [12], a double-ring electron-positron collider with a peak luminosity of 10 33 cm −2 s −1 at the c.m. energy of 3.770 GeV at a beam current of 0.93 A, operates in the center-of-mass energy region from 2 to 4.6 GeV. The Beijing Spectrometer III (BESIII) [12], which operates at the BEPCII storage ring, covers 93% of 4π solid angle. A small-cell helium-gas-based main drift chamber (MDC) provides charged particle momentum and ionization loss (dE/dx) measurements for charged particle identification (PID). The momentum resolution is better than 0.5% at 1 GeV/c in a magnetic field of 1 T. A plastic scintillator time-of-flight (TOF) system with a time resolution of 80 ps (110 ps) for the barrel (end-caps) is utilized for additional charged particle identification. A CsI(Tl) crystal electromagnetic calorimeter obtains a photon energy resolution of 2.5% (5%) at 1 GeV in the barrel (end-caps). A resistive-plate-counter muon system provides a position resolution of 2 cm and can detect muon tracks with momenta greater than 0.5 GeV/c.
All data samples used in this analysis are listed in Ta-ble I. The c.m. energies are measured using the dimuon process e + e − → (γ ISR/FSR )µ + µ − with an uncertainty of 0.8 MeV [13], and the integrated luminosities are measured with large-angle Bhabha scattering events with an uncertainty of 1.0% [14][15][16], where FSR denotes final state radiation. The data sets with c.m. energy above 4 GeV are named "XY Z data". The geant4-based [17] Monte Carlo (MC) simulation framework boost [18], consisting of event generators, the detector geometry, and detector response, is used to evaluate the detector efficiency, estimate ISR correction, optimize event selection criteria, and analyze backgrounds. The effects of FSR are simulated by the photos [19] package. Exclusive phase space (PHSP) MC samples for signal modes are modeled with kkmc [20][21][22] and besevtgen [23,24], which consider ISR effects and beam energy spreads. Seven 'inclusive' MC simulated data samples at the different c.m. energies, equivalent to the respective integrated luminosity of each data set, are produced to investigate potential backgrounds. The main known processes and decay modes are generated by besevtgen with cross sections or branching fractions obtained from the Particle Data Group (PDG) [25], and the remaining unmeasured events associated with charmonium decays or open charm processes are generated with lundcharm [23,26], while continuum light hadronic events are generated with pythia [27].

III. DATA ANALYSIS
In this analysis, the n andn candidates are not reconstructed, and a partial reconstruction technique is employed to select the signal events of interest, i.e., we reconstruct p, K 0 S (→ π + π − ), and K − only, whilen is treated as a missing particle. The presence of an is inferred using the mass recoiling against the pK 0 S K − system, M rec (pK 0 where E is the c.m. energy, E i is the energy of the ith track, and − → P i is the vector sum of the track momenta. For signal candidate events, the distribution of M rec peaks at then nominal mass [25].
Events with at least two positively charged tracks and two negatively charged tracks are selected. All charged tracks must be well reconstructed in the MDC with |cosθ| < 0.93, where θ is the polar angle between the charged track and the positron beam direction.
The K 0 S candidate is reconstructed with a pair of oppositely charged pions, where the point of the closest approach to the e + e − interaction point is required to be within ±20 cm in the beam direction. A charged track is identified as a pion by using the combined TOF and dE/dx information. To suppress random combinatorial backgrounds, we require that the π + π − pair satisfies a secondary vertex fit [28], and the decay length, which is the distance between production and decay vertexes, is required to be greater than twice its resolution. If there is more than one π + π − combination in an event, the one with the smallest χ 2 of the secondary vertex fit is retained. A K 0 S signal is required to have the π + π − invariant mass within |M π + π − −m K 0 S | ≤ 10 MeV/c 2 , where m K 0 S is the K 0 S nominal mass [25]. After the selection of the two π tracks from K 0 S decays, the remaining charged tracks are assumed to be p or K, and the point of closest approach of these tracks to the e + e − interaction point must be within ±10 cm in the beam direction and within 1 cm in the plane perpendicular to the beam. The net charge of the p and K combination must be zero, and multiple combinations of p and K are permitted.
If multiple pK − combinations are found, then candidate whose M rec (pK 0 S K − ) is closest to the world average anti-neutron mass value [25] is selected. After the above event selection criteria applied, a clearn signal is observed in the M rec (pK 0 S K − ) at each c.m. energy, as shown in Fig. 1.
To investigate non-K 0 S backgrounds, the K 0 S mass sideband regions are selected as 0.4676 < M π + π − < 0.4776 GeV/c 2 or 0.5176 < M π + π − < 0.5276 GeV/c 2 . According to the analysis of the inclusive MC samples in the K 0 S mass side-band regions, the main backgrounds are from many processes with the pK − π + π −n final state with one weakly decaying hyperon like Λ or Σ involved, where a small peak exists, as shown in Fig. 1 in the shaded histograms. Other background events, which form a smooth distribution in the M rec (pK − K 0 S ) spectra around 0.94 GeV/c 2 , are from numerous other processes, but none of them is dominant.
At each c.m. energy, an unbinned maximum likelihood fit to the M rec (pK 0 S K − ) spectra is performed to determine the signal and background yields in the selected candidates within [0.80, 1.10] GeV/c 2 and the normalized K 0 S mass side-band events. Then signal shape is obtained through the MC simulation at each c.m. energy smeared with a Gaussian function to account for the difference in the resolution between the data and the MC simulation. The samen line shape is used for K 0 S mass side-band events, and the other, non-peaking background contribution is described by a first-order polynomial function. Another first-order polynomial function associated with the function of the fit to K 0 S sidebands represents the remaining background contribution. The parameters of the Gaussian function and the two first-order polynomials are left free. The fits of the M rec (pK 0 S K − ) spectra at the seven c.m. energies are shown in Fig. 1, and the signal yields along with other numerical results are summarized in Table I. The dots with uncertainty bars are the signal candidate events in data, and the green-shaded histograms are shown as the normalized K 0 S mass side-band events in data. The red solid curves show the total fits, the blue dashed lines are the total background components of the fits, and the violet long dashed curves are the fits to K 0 S mass side-band events. Table I. The c.m. energy ( √ s), integrated luminosity (L), detection efficiency (ǫ), vacuum polarization ( 1 |1−Π| 2 ) and radiative correction factor (1 + δ), number of fittedn signal events (Nsig) excluding the peaking backgrounds, and Born cross section (σB) of e + e − → pK 0 Sn K − at each energy point. The first uncertainties are statistical and the second systematic. In addition to the common selection criteria for p or K candidates, the combined TOF and dE/dx information is used to calculate χ 2 PID (i) (i = p, K) for each hadron (i) hypothesis, and a one-constraint (1C) kinematic fit is performed with the p, K 0 S , K − , andn combination by constraining the missing mass of the undetectedn to its nominal mass [25]. The combination with the minimum in an event is selected, and we require χ 2 sum < 20, where χ 2 1C is the χ 2 of the kinematic fit.
After the event selection requirements have been applied, clear Λ(1520) signals are found for the processes e + e − → Λ(1520)nK 0 S and e + e − → Λ(1520)pK + in data at the c.m. energy of 3.773 GeV and in the full XY Z data, as shown in Fig. 2.
According to the analysis of the inclusive MC samples, the dominant background events are from processes with pK 0 Sn K − final states without a Λ(1520). The remaining backgrounds are from many processes with only a few events in each mode. No peaking background is found in the pK − or nK 0 S invariant mass spectra at around 1520 MeV/c 2 .
An unbinned maximum likelihood fit is performed to the pK − and nK 0 S invariant mass distributions to determine the Λ(1520) yields individually. The Λ(1520) signal is described by a D-wave relativistic Breit-Wigner (BW) function with an energy-dependent width Γ pK convolved with a Gaussian function. The Λ(1520) mass and width are fixed to the world average values [25]. The pK − and nK 0 S mass resolutions of the signal Gaussian functions are the same and fixed to 3 MeV/c 2 as determined from fits to the individual MC distributions. The background shape is parameterized by an ARGUS function [29].
To avoid double counting from e + e − → Λ(1520)Λ(1520) in calculating the Born cross sections of e + e − → Λ(1520)nK 0 S and Λ(1520)pK + , Λ(1520)Λ(1520) pair events are subtracted from the three-body decays, since it is difficult to perform a two-dimensional (2D) fit in the full range. The final number of the signal events and the corresponding statistical significance of the e + e − → Λ(1520)nK 0 S and Λ(1520)pK + signal at each c.m. energy are listed in Table II. The statistical significance is calculated using −2 ln(L 0 /L max ), where L max and L 0 are the likelihoods of the fits with and without the Λ(1520) signal included, respectively.
We extend the unbinned maximum likelihood fit described above into a 2D fit to the pK − versusnK 0 S mass spectra to determine the yield of the process e + e − → Λ(1520)Λ(1520) → pK −n K 0 S . We assume that the two discriminating variables M (pK − ) and M (nK 0 S ) are uncorrelated, and the 2D probability density function (PDF) is the product of two one-dimensional (1D) PDFs for the two variables. The total PDFs include four components: Λ(1520)Λ(1520), Λ(1520)nK 0 S , pK −Λ (1520), and non−Λ(1520). The signal shapes of Λ(1520)(→ pK − ) andΛ(1520)(→nK 0 S ) in the 2D fit are the same as those in e + e − → Λ(1520)nK 0 S and Λ(1520)pK + , respectively. All of the backgrounds are parameterized by an ARGUS function [29]. The projections of the 2D fits in data at the c.m. energy of 3.773 GeV and in the full XY Z data are shown in Figs. 3(a-d).
Since only a few Λ(1520)Λ(1520) pair signal events are observed and the statistical significance is less than 3.0 standard deviations (σ) at each c.m. energy, the upper limit on the number of signal events, N up sig , is determined at the 90% confidence level (C.L.) by solving where x is the number of fitted signal events, and L(x) is the likelihood function in the fit to data.
To account for the systematic uncertainties, the likelihood distribution is convolved with a Gaussian function G(x; 0, σ) with a standard deviation of σ = x × ∆, where µ is the expected number of signal events, L ′ (µ) indicates the expected likelihood distribution, and ∆ refers to the total relative systematic uncertainty discussed in Section V. The upper limit on the number of Λ(1520)Λ(1520) pair events and statistical significance at each energy are listed in Table II.
To investigate other two-body invariant mass distributions, we apply the further requirement |M (pK − /nK 0 S )− 1.5195| > 0.025 GeV/c 2 to veto the Λ(1520) resonance. The pK 0 S , nK + , pn and K 0 S K − invariant mass spectra in the full data are shown in Figs. 4 (a-d). No significant structures are visible.

IV. CROSS SECTION MEASUREMENT
The Born cross section is calculated using: where N sig is the number of signal events, L int is the integrated luminosity, 1 + δ is the radiative correction factor obtained from a QED calculation with 1% accuracy [31], 1 |1−Π| 2 is the vacuum polarization factor [32,33], ǫ is the detection efficiency from the PHSP MC simulation, B is the product of intermediate branching fractions, i.e. B(K 0 S → π + π − ) for e + e − → pK 0 for e + e − → Λ(1520)Λ(1520). The branching fractions B(K 0 S → π + π − ), B[Λ(1520) → pK − ], and B[Λ(1520) → nK 0 S ] are 0.692, 0.225, and 0.1125 [25], respectively. All calculated Born cross sections or the 90% C.L. upper limits on the Born cross sections are summarized in Table I for the pK 0 sn K − mode and Table II for the Λ(1520) modes.
The Born cross sections of e + e − → pK 0 Sn K − are shown in Fig. 5 at c.m. energies between 3.773 and 4.6 GeV. We fit the 1/s k dependence of the cross sections, as shown in Fig. 5 with a dashed line. The fit gives k = 1.9 ± 0.1 ± 0.2 with the goodness of the fit χ 2 /ndf = 3.8/5, where the first uncertainty is statistical and the second systematic.  Fig. 6. A fit with the functional form σ 0 /s k to each measured Born cross section is performed, as shown in Fig. 6 with the dashed lines, where σ 0 is a constant and k is a free parameter. The fits yield k = 2.8 ± 0.9 ± 0.2 (χ 2 /ndf = 2.0/5) and 3.5 ± 0.6 ± 0.2 (χ 2 /ndf = 5.6/5) for the Born cross sections of e + e − → Λ(1520)nK 0 S and e + e − → Λ(1520)pK + , respectively, where the first uncertainties are statistical and the second systematic. , integrated luminosity (L), detection efficiency (ǫ), vacuum polarization ( 1 |1−Π| 2 ), radiative correction factor (1 + δ), number of observed signal events (Nsig), statistical signal significance (S), and calculated (90% C.L. upper limit of) Born cross section (σB) are listed for the studied Λ(1520) modes at each energy point. The first uncertainties are statistical and the second systematic. For the 90% C.L. upper limits, the systematic uncertainties have been included.   Table III, and the total systematic uncertainty is obtained by adding all contributions in quadrature assuming that each source is independent. Detailed descriptions of the estimates of the systematic uncertainties are listed in the following subsections.

A. Integrated luminosity, tracking, and PID
The luminosity is measured using large-angle Bhabha scattering events with a total uncertainty of less than 1.0% [14][15][16], which is taken as its systematic uncertainty at each c.m. energy.
Using the control samples of e + e − → ppπ + π − at √ s > 4 GeV and J/ψ → K 0 S K ± π ∓ events, the tracking efficiency difference between MC simulation and data is found to be 2.0% for each proton and 1.0% for each kaon. Not counting the two charged pions from K 0 S decays, there are two charged tracks (p, K), and the uncertainty in the tracking efficiency is 3.0%.
Based on the measurements of the particle identification efficiencies of protons from e + e − → ppπ + π − events and kaons from e + e − → K + K − π + π − events, the difference between data and MC simulation yields uncertainties of 1.0% for each proton and 2.0% for each kaon. Thus, a total uncertainty associated with the PID of 3.0% is assigned for the Λ(1520) modes. 6.9 7.0 6.9 7.1 7.0 7.1 Λ(1520)pK + 6.5 6.7 6.7 6.7 6.9 6.8 6.9 Λ(1520)Λ(1520) 9.1 8.9 8.9 8.9 9.0 9.0 9.0 B. K 0 S reconstruction and K 0 S mass window The K 0 S reconstruction efficiency is studied using two control samples: J/ψ → K * (892) ± K ∓ → K 0 S π ± K ∓ and J/ψ → φK 0 S K ± π ∓ . The difference in the K 0 S reconstruction efficiency between the MC simulation and the data is 1.2% [34]. Considering the additional PID requirements for two opposite-charged pions from K 0 S decays in our analysis, the systematic uncertainty for K 0 S reconstruction is conservatively taken as 2.3%, where the PID efficiency difference of 1% for each pion between MC and data is included [35].
The uncertainty attributed to the K 0 S mass window requirement, which originates from the mass resolution difference between the data and the MC simulation, is estimated using |ε data − ε MC |/ε data , where ε data is the efficiency of applying the K 0 S mass window requirement by extracting K 0 S signal in the π + π − invariant mass spectrum of the data at each c.m. energy, and ε MC is the analogous efficiency from the MC simulation. The difference between the data and the MC simulation is considered as the systematic uncertainty at each c.m. energy.

C. Kinematic fit and MC generator
A correction is applied to the track helix parameters in the MC simulation to make the χ 2 distribution of the 1C kinematic fit from the MC simulation agree better with data [36]. The difference between the efficiencies with and without the correction is taken as the systematic uncertainty. In this analysis, the detection efficien-cies from the MC samples with the corrected track helix parameters are taken as nominal results.
The detection efficiencies are obtained from the PHSP MC samples for e + e − → Λ(1520)nK 0 S and Λ(1520)pK + . To estimate the uncertainty attributed to the MC generator, we determine the efficiency correction factor by comparing the Dalitz plots between the data and the MC simulation at the c.m. energy of 3.773 GeV, and the correction factor is assigned to the other c.m. energies due to the limited statistics. The relative differences in the efficiency with and without correction are 3.3% and 1.5% for e + e − → Λ(1520)nK 0 S and e + e − → Λ(1520)pK + , respectively, which are taken as the systematic uncertainties due to the MC generator at all of c.m. energies.
For e + e − → Λ(1520)Λ(1520), different MC samples with angular distributions of 1 + cos 2 θ and 1 − cos 2 θ are generated, where θ is the polar angle of Λ(1520) in e + e − c.m. frame. The largest difference of 5.5%, compared to the PHSP MC efficiency, is taken as the systematic uncertainty attributed to the MC generator.

D. Fit procedure
Signal yields are determined from the fits to the M rec (pK 0 S K − ) spectra for the pK 0 Sn K − mode and the M (pK − ) or M (nK 0 S ) spectra for Λ(1520) modes. Alternative signal and background shapes as well as multiple fit ranges are used to estimate the systematic uncertainty in the fit procedure. We generated simulated pseudoexperiments out of the fit to the data with alternative shapes and fitted them back using the nominal model. Any deviation of the pull distributions from the normal gives the systematic effect.

1.
pK 0 Sn K − mode 1) Signal shape: In the nominal fit, then signal shape is obtained from the MC simulation directly convolved with a Gaussian function. Alternatively, the incoherent sum of a Gaussian function and a Novosibirsk function [37] is taken as then signal shape.
2) Background shape: The background shape without K 0 S mass side-band events is described by a firstorder polynomial function in the nominal fit, and a second-order polynomial function is used to estimate the systematic uncertainty due to background shape. Assuming that all of the above sources are independent, the systematic uncertainties associated with the fit procedure are the quadrature sum of above three sources.

Λ(1520) modes
Since only a few Λ(1520) signal events are observed in each Λ(1520) mode at each c.m. energy in the XY Z data, we use the full XY Z data to estimate the uncertainty associated with the fit procedure for each XY Z data sample. 1) Signal shape: In the nominal fit, the Λ(1520) signal is described by a D-wave relativistic BW function convolved with a Gaussian function. Alternatively, the Λ(1520) signal shape is obtained from the signal MC simulated shape convolved with a Gaussian function.
2) Background shape: In the nominal fit, the background shape is described by an ARGUS function [29]. To estimate the uncertainty due to background shape, we use the alternative parameterized exponential function Assuming that all of the above sources are independent and adding them in quadrature, we obtain the systematic uncertainties associated with the fit procedure for e + e − → Λ(1520)nK 0 S , Λ(1520)pK + , and Λ(1520)Λ(1520), respectively.

E. Radiative correction, vacuum polarization and intermediate decays
The line shape of the cross section for e + e − → pK 0 Sn K − affects the radiative correction factor (1 + δ) and detection efficiency (ǫ). For our nominal results, the dependence of the Born cross sections on the c.m. energy are σ B ∝ 1/s 1.9 , 1/s 2.8 and 1/s 3.5 for e + e − → pK 0 Sn K − , Λ(1520)nK − , and Λ(1520)pK 0 S , respectively. The dependence is assumed to be 1/s 3 for e + e − → Λ(1520)Λ(1520) due to the limited statistics. We change the energy dependence to 1/s for the above processes, and the largest difference in (1 + δ)ǫ among the modes is conservatively taken as the systematic uncertainty.
The vacuum polarization factor is calculated with an uncertainty of less than 0.1% [32], which is negligible compared with other sources of uncertainties.

VI. SYSTEMATIC UNCERTAINTIES FOR FITS TO THE CROSS-SECTION DISTRIBUTIONS
The systematic uncertainty in the measured cross section is divided into two categories: the uncorrelated part among the different c.m. energies, which comes from the fit to then or Λ(1520) mass spectrum to determine the signal yields, and the correlated part, which includes all other uncertainties common for the whole data set. By including the uncorrelated uncertainty in the fit to the cross-section distributions, the systematic uncertainties in the parameter k are estimated to be 8.6%, 6.5%, and 4.3% for the decays of e + e − → pK 0 Sn K − , Λ(1520)nK 0 S , and Λ(1520)pK + , respectively.

VII. SUMMARY AND DISCUSSION
In summary, we study for the first time the processes e + e − → pK 0 Sn K − , Λ(1520)nK 0 S , Λ(1520)pK + , and Λ(1520)Λ(1520) using data samples with a total integrated luminosity of 7.4 fb −1 collected with the BESIII detector at c.m. energies of 3.773, 4.008, 4.226, 4.258, 4.358, 4.416, and 4.600 GeV. The Born cross sections of e + e − → pK 0 Sn K − are measured, and no structure in the cross section line shape between 3.773 and 4.60 GeV is visible.
Furthermore, Λ(1520) signals are observed in the pK − and nK 0 S invariant mass spectra for the processes e + e − → Λ(1520)nK 0 S and Λ(1520)pK + with statistical significances equal to or greater than 3.0σ, and the corresponding Born cross sections are measured. For e + e − → Λ(1520)Λ(1520), the statistical significances are less than 3.0σ, and the 90% C.L. upper limits on the Born cross sections are determined. No other significant structure is found in the pK − , nK 0 S , pK 0 S , nK + , pn or K 0 S K − invariant mass spectra in any of the data samples. As a consequence, no light hadron decay modes of Y states or conventional charmonium resonances are observed in our analysis. However, we note that there is an evident difference in line shape and magnitude of the measured cross sections between e + e − → Λ(1520)nK 0 S and e + e − → Λ(1520)pK + (the statistical significance of the cross-section difference is 3.1σ at the c.m. energy of 3.770 GeV). Such an isospin violating effect may be due to the interference between I = 1 and I = 0 final states. The final states Λ(1520)nK 0 S and Λ(1520)pK + can be produced from pK − andnK 0 S systems either in I = 1 or I = 0 states, namely, excited Σ * or Λ * states. These two states can decay into both pK − andnK 0 S final states, but with a sign difference from Clebsch-Gordan coefficients. Another possible approach is e + e − → K * K with highly excited K * decays into Λ * p or Λ * n . The K * K system can be produced from I = 1 (excited ρ) or I = 0 (excited ω or φ) states, where the interference effect can occur. If the final state is pN * or nN * , a similar pattern could be observed. More experimental data are desirable to confirm these interpretations and speculations in the future.