Measurement of the $e^{+}e^{-}\rightarrow\Lambda\bar{\Lambda}$ cross section from threshold to 3.00 GeV using events with initial-state radiation

Electromagnetic form factors (EMFFs), which parametrize the inner structure of hadrons, are fundamental observables for understanding the strong interaction. In the timelike region, EMFFs are extensively studied in electron-positron collisions by measuring hadron pair production cross sections. For a spin-$1/2$ baryon ($B$), the cross section in the Born approximation of the one-photon-exchange process $e^{+}e^{-}\to B\bar{B}$ is parameterized in terms of electric and magnetic form factors $G_{E}$ and $G_{M}$ by Cabibbo and Gatto (1961):

$$\sigma^{B}(s)=\frac{4\pi\alpha^{2}C\beta}{3s}\left[\left|G_{M}(s)\right|^{2}+\frac{2m^{2}_{B}c^{2}}{s}\left|G_{E}(s)\right|^{2}\right],$$

(1)

where $\alpha$ is the fine-structure constant, $C$ is the Coulomb correction factor Brodsky and Lebed (2009), $\beta=\sqrt{1-4m^{2}_{B}c^{4}/s}$ is a phase-space (PHSP) factor, $s$ is the square of the center-of-mass (c.m.) energy, $m_{B}$ is the mass of the baryon, and $c$ is the speed of light. $C$ accounts for the electromagnetic interaction of the fermions in the final state, and in the point-like approximation, it is 1 for neutral baryons and $y/(1-e^{-y})$ with $y=\pi\alpha\sqrt{1-\beta^{2}}/\beta$ for charged baryons. Therefore, for charged baryon-pairs, the factor of $\beta$ due to PHSP is canceled by the Coulomb factor, which results in a non-zero cross section at the threshold when $\beta=0$. However, there is no cancelation in the neutral baryon-pair case, so the cross section is zero.

There have been many experimental studies on the charged and neutral baryon-pair production cross sections in the past decades, such as $e^{+}e^{-}\rightarrow p\bar{p}$ Ablikim et al. (2020a); Lees et al. (2013), $e^{+}e^{-}\rightarrow n\bar{n}$ Ablikim et al. (2021a), $e^{+}e^{-}\rightarrow\Lambda\bar{\Lambda}$ Aubert et al. (2007); Ablikim et al. (2018a, 2019a); Bisello et al. (1990), $e^{+}e^{-}\rightarrow\Sigma\bar{\Sigma}$ Ablikim et al. (2021b, 2022a), $e^{+}e^{-}\rightarrow\Xi\bar{\Xi}$ Ablikim et al. (2021c, d), and $e^{+}e^{-}\rightarrow\Lambda_{c}^{+}\bar{\Lambda}_{c}^{-}$ Ablikim et al. (2018b). Although the conclusions for some channels are questionable due to large uncertainties, there is a general tendency in the production cross sections for these baryon pairs to have a step near the threshold, which then decreases with the increase of the c.m. energy of the baryon pair Huang and Ferroli (2021).

The cross section of the process $e^{+}e^{-}\rightarrow\Lambda\bar{\Lambda}$ very close to the threshold has been measured in both the BaBar and the BESIII experiments. In the BaBar experiment, the cross section from the $\Lambda\bar{\Lambda}$ production threshold up to $M_{\Lambda\bar{\Lambda}}=$ 2.27 GeV$/c^{2}$ was measured as $204_{-60}^{+62}\pm 22$ pb Aubert et al. (2007). This result indicates a possible non-zero cross section at threshold which is in conflict with Eq. (1). However, due to the wide $\Lambda\bar{\Lambda}$ mass interval and large uncertainties, a solid conclusion cannot be drawn. The BESIII experiment also measured the cross section at the c.m. energy ($\sqrt{s}$) of $2.2324$ GeV, which is only $1.0$ MeV above the $\Lambda\bar{\Lambda}$ production threshold, to be $305\pm 45_{-36}^{+66}$ pb Ablikim et al. (2018a). This indicates a threshold enhancement phenomenon in the process $e^{+}e^{-}\rightarrow\Lambda\bar{\Lambda}$. Interestingly, in both the BaBar and BESIII experiments, a jump was observed in the process $e^{+}e^{-}\rightarrow K^{+}K^{-}K^{+}K^{-}$ near the $\Lambda\bar{\Lambda}$ production threshold Lees et al. (2012); Ablikim et al. (2019b).

To explain the near threshold enhancement, some theoretical studies have been performed, in which the effects of final-state radiation Haidenbauer and Meißner (2016) and vector-meson resonances Cao et al. (2018); Yang et al. (2019) have been taken into account. The enhancement in the case of neutral baryons may also be explained by an electromagnetic interaction occurring at the quark level Baldini et al. (2009). However, experimentally, the cross section measurements of $e^{+}e^{-}\rightarrow\Lambda\bar{\Lambda}$ near threshold are still limited and more measurements are needed to further understand this phenomenon.

The cross section and EMFFs of the $\Lambda$ hyperon have been measured via the annihilation channel $e^{+}e^{-}\rightarrow\Lambda\bar{\Lambda}$ using the energy scan technique Ablikim et al. (2018a, 2019a); Bisello et al. (1990), in which the c.m. energy of the collider is varied according to the experimental plan and the cross section is measured at each c.m. energy. In addition, the radiative return channel $e^{+}e^{-}\rightarrow\gamma\Lambda\bar{\Lambda}$ as illustrated in Fig. 1, where $\gamma$ is a hard photon from the initial-state radiation (ISR) process, offers a technique complementary to the energy scan technique for the $\Lambda$ hyperon cross section measurement. This technique has been used in the BaBar experiment to measure the cross section and effective form factor of the $\Lambda$ hyperon Aubert et al. (2007).

The differential Born cross section for the $e^{+}e^{-}\rightarrow\gamma\Lambda\bar{\Lambda}$ process, integrated over the $\Lambda(\bar{\Lambda})$ momenta and the photon polar angle, is written as Druzhinin et al. (2011):

$$\frac{d\sigma_{e^{+}e^{-}\rightarrow\gamma\Lambda\bar{\Lambda}}\left(q^{2}\right)}{dq^{2}}=\frac{1}{s}W(s,x)\sigma_{\Lambda\bar{\Lambda}}\left(q^{2}\right),$$

(2)

where $\sigma_{\Lambda\bar{\Lambda}}(q^{2})$ is the cross section for the $e^{+}e^{-}\rightarrow\Lambda\bar{\Lambda}$ process, $q$ is the momentum transfer of the virtual photon whose squared value represents the invariant mass squared of $\Lambda\bar{\Lambda}$, $x=\frac{2E_{\gamma}^{*}}{\sqrt{s}}=1-\frac{q^{2}}{s}$, and $E_{\gamma}^{*}$ is the energy of the ISR photon in the $e^{+}e^{-}$ c.m. system. The function Kuraev and Fadin (1985)

	$\displaystyle W(s,x)=kx^{k-1}\left[1+\frac{\alpha}{\pi}\left(\frac{\pi^{2}}{3}-\frac{1}{2}\right)+\frac{3}{4}k+k^{2}\left(\frac{37}{96}-\frac{\pi^{2}}{12}-\frac{1}{72}\ln\frac{s}{m_{e}^{2}}\right)\right]-k\left(1-\frac{1}{2}x\right)$		(3)
	$\displaystyle+\frac{1}{8}k^{2}\left[4\left(2-x\right)\ln\frac{1}{x}-\frac{1+3\left(1-x\right)^{2}}{x}\ln\left(1-x\right)-6+x\right],k=\frac{2\alpha}{\pi}\left[\ln\frac{s}{m_{e}^{2}}-1\right],$

describes the probability for the emission of an ISR photon with energy fraction $x$, and $m_{e}$ is the electron mass.

In this analysis, we present the measurement of the $e^{+}e^{-}\rightarrow\Lambda\bar{\Lambda}$ cross section from the production threshold up to 3.00 GeV$/c^{2}$ using the ISR process $e^{+}e^{-}\rightarrow\gamma\Lambda\bar{\Lambda}$. The used data sets, corresponding to a total integrated luminosity of 11.957 fb^-1, are collected at twelve c.m. energies between 3.773 and 4.258 GeV with the BESIII detector Ablikim et al. (2010) at the BEPCII Collider  Wang (2006).

The BESIII detector Ablikim et al. (2010) records symmetric $e^{+}e^{-}$ collisions provided by the BEPCII storage ring Wang (2006) in the c.m. energy range from 2.00 up to 4.95 GeV, with a peak luminosity of $1\times 10^{33}$ cm^-2s^-1 achieved at $\sqrt{s}=3.77\;\text{GeV}$. BESIII has collected large data samples in this energy region Ablikim et al. (2020b). 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 magnetic field Huang et al. (2022). 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 d$E/$d$x$ resolution is 6% for 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 in the TOF barrel region is 68 ps, while that in the end cap region used to be 110 ps. The end cap TOF system was updated in 2015 using multi-gap resistive plate chamber technology, providing a time resolution of 60 ps Li et al. ; Guo et al. ; Cao et al. (2020).

The experimental data sets used in this analysis are listed in Table 1. To optimize the event selection criteria, Monte Carlo (MC) simulations are performed with Geant4-based Agostinelli et al. (2003) software, which includes the description of geometry and material, the detector response and the digitization model, as well as a database for the detector running conditions and performances. In this analysis, the event generator ConExc Ping (2014) is used to generate the signal process $e^{+}e^{-}\rightarrow\gamma\Lambda\bar{\Lambda}$ ($\Lambda{\rightarrow}p\pi^{-}$ $\bar{\Lambda}\rightarrow\bar{p}\pi^{+}$) with 1 million events at the different c.m. energies up to ISR leading order (LO), i.e. with only one ISR photon, and vacuum polarization (VP) is included. The selection efficiencies are estimated by the signal MC samples. An alternative event generator, PHOKHARA10.0 pho , is used to study the systematic uncertainty of the MC model. The cross section lineshape used for the generation of the signal MC samples is from Ref. Li et al. (2022). Inclusive MC samples at $\sqrt{s}=3.773$ and $4.178$ GeV are used to investigate possible background contamination. They consist of inclusive hadronic processes ($e^{+}e^{-}\rightarrow q\bar{q}$, $q=u,d,s$) modeled with the LUARLW Ablikim et al. (2022b) at $\sqrt{s}=3.773$ GeV and KKMC Jadach et al. (2001, 2000) at $\sqrt{s}=4.178$ GeV, and the ISR production of vector charmonium states ($e^{+}e^{-}\rightarrow\gamma J/\psi$, $\gamma\psi(2S)$, $\gamma\psi(3773)$) generated with BesEvtGen Ping (2008) using the VECTORISR model Bonneau and Martin (1971); Lange (2001). In addition, several exclusive MC samples are generated to study the background, with different event generators and models.

The complete process we study is $e^{+}e^{-}\rightarrow\gamma\Lambda\bar{\Lambda}\rightarrow\gamma(p\pi^{-})(\bar{p}\pi^{+})$, with the final state $\gamma p\pi^{-}\bar{p}\pi^{+}$, where $\gamma$ is the ISR photon. To provide a clean sample in the threshold region, the ISR photon is detected (tagged). However, the differential cross section of the ISR reaction (such as $e^{+}e^{-}\rightarrow\gamma\Lambda\bar{\Lambda}$) as a function of the ISR photon polar angle reaches its highest value when the photon is emitted at a small angle relative to the direction of the electron (or positron) beam Druzhinin et al. (2011). Since this is out of the angular acceptance of the EMC, photons falling in this region cannot be detected, resulting in a reduction of signal efficiency. Moreover, the detection efficiency is further reduced by the low momenta of the pions, which, according to the study of the signal MC samples, are mostly less than 0.2 GeV$/c$. We categorize the reconstruction of signal candidates into two modes: mode I corresponds to fully reconstructed events, i.e. all particles in the final state are identified; in mode II, a partial reconstruction method with a missing pion is used to increase the efficiency.

Charged tracks detected in the MDC are required to be within $\left|\cos\theta\right|<0.93$, where $\theta$ is the polar angle with respect to the $z$ axis, which is the symmetry axis of the MDC. The distance of closest approach of each charged track to the interaction point must be less than 30 cm along the $z$ direction and less than 10 cm in the transverse plane. For each signal candidate, at least three charged tracks are required.

The combined information of d$E/$d$x$ and TOF is used to calculate particle identification (PID) probabilities for the pion, kaon, and proton hypotheses, and the particle type with the highest probability is assigned to the track.

A secondary vertex fit is performed to obtain the decay vertex of the $\Lambda(\bar{\Lambda})$ candidate, and the $\Lambda(\bar{\Lambda})$ candidate is reconstructed by fitting the $p\pi^{-}(\bar{p}\pi^{+})$ tracks to a common decay vertex. If there is more than one $\Lambda(\bar{\Lambda})$ candidate, the one with the minimum chi-square value of the secondary vertex fit is selected. The reconstructed mass of $\Lambda(\bar{\Lambda})$ candidate $\left(M_{\Lambda(\bar{\Lambda})}\right)$ is required to be within 6.4 MeV$/c^{2}$ of the nominal $\Lambda$ mass ($m_{\Lambda}$) Workman et al. (2022), as shown in Fig. 2. There is no requirement on the decay length of $\Lambda(\bar{\Lambda})$. Both a $\Lambda$ and a $\bar{\Lambda}$ are required in mode I, while either a $\Lambda$ or a $\bar{\Lambda}$ is required in mode II.

Information on the electromagnetic showers in the EMC is used to select the photon candidates. It is required that the shower time is within 700 ns of the event’s start time to suppress electronic noise and energy deposits unrelated to the events. A photon candidate is selected if its deposited energy is greater than 0.4 GeV. For each candidate signal event, at least one photon is required which is considered as the ISR photon.

A kinematic fit is applied to further suppress background. For mode I, a four-constraint (4C) kinematic fit requiring energy-momentum conservation under the hypothesis of a $\gamma\Lambda\bar{\Lambda}$ final state is applied to the signal candidates. If there is more than one photon candidate, the combination with the minimum $\chi^{2}_{\rm 4C}$ is selected. To suppress the background with one more photon than the signal process, we require $\chi^{2}_{\rm 4C}\leq\chi^{2}_{\rm 4C,\gamma\gamma}$, where $\chi^{2}_{\rm 4C}$ and $\chi^{2}_{\rm 4C,\gamma\gamma}$ are the chi-square values under the hypotheses of $\gamma\Lambda\bar{\Lambda}$ and $\gamma\gamma\Lambda\bar{\Lambda}$ final states. For mode II, a one-constraint (1C) kinematic fit with a missing $\pi^{+}(\pi^{-})$ under the hypothesis of a $\gamma\Lambda\bar{p}\pi^{+}(\gamma\bar{\Lambda}p\pi^{-})$ final state is applied to the signal candidates. Combining all $\gamma\bar{p}(\gamma p)$ pairs with the reconstructed $\Lambda(\bar{\Lambda})$, 1C kinematic fits are applied with the invariant mass of $\bar{p}\pi^{+}(p\pi^{-})$ being constrained to the nominal $\Lambda$ mass Workman et al. (2022) and the mass of $\pi^{+}(\pi^{-})$ being unconstrained. The $\gamma\bar{p}(\gamma p)$ combination with the minimum $\chi^{2}_{\rm 1C}$ is selected, where $\chi_{\rm 1C}^{2}$ is the chi-square of the 1C kinematic fit. A requirement of $\chi^{2}_{\rm 4C}\leq 50$ ($\chi^{2}_{\rm 1C}\leq 5$) is optimized for the signal candidates for mode I (mode II).

For the candidates of mode II, the distribution of the mass squared of the missing $\pi$ ($M_{\pi}^{2}$), obtained from energy-momentum conservation, is shown in Fig. 3. To suppress background, a requirement of $0.012\leq{M_{\pi}^{2}}\leq 0.025$ GeV${}^{2}/c^{4}$ is applied.

The distribution of the selection efficiencies obtained from signal MC samples as a function of invariant mass of $\Lambda\bar{\Lambda}$ ($M_{\Lambda\bar{\Lambda}}$) is shown in Fig. 4, where the efficiencies at the c.m. energies between 4.128 and 4.258 GeV are combined and weighted according to the effective luminosity of the ISR process. It should be noted that to improve the mass resolution of $M_{\Lambda\bar{\Lambda}}$, we correct $M_{\Lambda\bar{\Lambda}}$ to $\left(M_{\Lambda\bar{\Lambda}}-M_{\Lambda}-M_{\bar{\Lambda}}+2\times m_{\Lambda}\right)$. The mass resolution is given by the root-mean-square deviation of $\left(M_{\Lambda\bar{\Lambda}}-M_{\Lambda\bar{\Lambda}}^{\rm truth}\right)$ of the signal MC sample, where $M_{\Lambda\bar{\Lambda}}^{\rm truth}$ is the set value of the invariant mass of $\Lambda\bar{\Lambda}$ when generating the MC events. In this paper, the correction of the $M_{\Lambda\bar{\Lambda}}$ is implied unless specified. The $M_{\Lambda\bar{\Lambda}}$ spectrum of the accepted candidates from all data sets is shown in Fig. 5, in which 817 events are retained. The contributions from $J/\psi\rightarrow\Lambda\bar{\Lambda}$ and $\psi(2S)\rightarrow\Lambda\bar{\Lambda}$ decays are clearly seen. About 60% of the signal candidates have $M_{\Lambda\bar{\Lambda}}$ below 3.00 GeV/$c^{2}$, and the number of signal candidates ($N_{\rm obs}$) in each $M_{\Lambda\bar{\Lambda}}$ interval is listed in the first column of Table 2.

Potential background channels are investigated in the inclusive MC samples with a topology analysis Zhou et al. (2021); they consist of channels containing $\Lambda\bar{\Lambda}$ and channels without $\Lambda\bar{\Lambda}$. The background channels containing $\Lambda\bar{\Lambda}$, such as the processes of $e^{+}e^{-}\rightarrow\pi^{0}\Lambda\bar{\Lambda}$, $e^{+}e^{-}\rightarrow\gamma(\Lambda\bar{\Sigma}^{0}+c.c.)$, and $e^{+}e^{-}\rightarrow\gamma J/\psi(\psi(2S))$ with $J/\psi(\psi(2S))$ decaying to $\gamma\Lambda\bar{\Lambda}$, are studied individually, while the non-$\Lambda\bar{\Lambda}$ background is estimated with the sideband method.

Events of $e^{+}e^{-}\rightarrow\pi^{0}\Lambda\bar{\Lambda}$ are easily mistaken as signal events if a soft photon from the high-energy $\pi^{0}$ is missing. A data-driven method is used to estimate their contribution. A sample of $\pi^{0}\Lambda\bar{\Lambda}$ events is selected from data, and its background is estimated with the sideband method. The sideband regions are chosen in the distribution of the invariant mass of $\gamma\gamma$ ($M_{\gamma\gamma}$). The number of events of this sample is calculated by $N_{\pi^{0}}^{\rm data}=N_{\pi^{0}}^{\rm SigReg}-N_{\pi^{0}}^{\rm Side}/2$, where $N_{\pi^{0}}^{\rm SigReg}$ and $N_{\pi^{0}}^{\rm Side}$ are the numbers of events from the signal and the sideband regions of the $\pi^{0}\Lambda\bar{\Lambda}$ sample, respectively. Next, the contribution from the remaining $\pi^{0}\Lambda\bar{\Lambda}$ background ($N_{\pi^{0}}^{\rm bkg}$) in the signal candidates is determined by:

$$N_{\pi^{0}}^{\rm bkg}=N_{\pi^{0}}^{\rm data}\times\frac{N_{\rm ISR}^{\rm MC}}{N_{\pi^{0}}^{\rm MC}},$$

(4)

where $N_{\rm ISR}^{\rm MC}$ and $N_{\pi^{0}}^{\rm MC}$ are the numbers of the events selected by the signal and $\pi^{0}\Lambda\bar{\Lambda}$ selection criteria from the $\pi^{0}\Lambda\bar{\Lambda}$ MC sample. The $\pi^{0}\Lambda\bar{\Lambda}$ MC sample is generated with the ConExc Ping (2014) event generator up to ISR LO, and the lineshape is obtained with the data sets collected at c.m. energies from 2.644 GeV to 3.080 GeV by BESIII.

In the reaction $e^{+}e^{-}\rightarrow\gamma(\Lambda\bar{\Sigma}^{0}+c.c.)$, the $\Sigma^{0}(\bar{\Sigma}^{0})$ decays to $\gamma\Lambda(\bar{\Lambda})$ with a branching ratio of 100% Workman et al. (2022), where the $\gamma$ has low energy. Therefore, if the photon from the $\Sigma^{0}(\bar{\Sigma}^{0})$ decay is missing, this event can be misidentified as signal. To estimate the background from this reaction, a MC sample with a total of 2 million events is generated with the ConExc Ping (2014) event generator up to ISR LO, and the lineshape used to generate the MC events is determined with the data sets collected at c.m. energies from 2.309 GeV to 3.080 GeV by BESIII. After applying the signal ($\gamma\Lambda\bar{\Lambda}$) selection criteria to this sample, we obtain the number of the surviving $\gamma(\Lambda\bar{\Sigma}^{0}+c.c.)$ events ($N_{\Lambda\Sigma}^{\rm MC}$). A scaling factor is obtained by $f=N_{\rm exp}/N_{\rm gen}$, where $N_{\rm exp}$ is the expected number of the $\gamma(\Lambda\bar{\Sigma}^{0}+c.c.)$ events estimated with the $(\Lambda\bar{\Sigma}^{0}+c.c.)$ cross section lineshape, and $N_{\rm gen}$ is the number of MC simulated events. Finally, the number of $\gamma(\Lambda\bar{\Sigma}^{0}+c.c.)$ background events ($N_{\Lambda\Sigma}^{\rm bkg}$) is estimated by $N_{\Lambda\Sigma}^{\rm bkg}=f\times N_{\Lambda\Sigma}^{\rm MC}$. Some other background channels, such as the processes $e^{+}e^{-}\rightarrow\eta\Lambda\bar{\Lambda}$ and $e^{+}e^{-}\rightarrow\gamma J/\psi(\psi(2S))$, are negligible.

Next, the sideband method is used to study the non-$\Lambda\bar{\Lambda}$ background. For mode I, two-dimensional (2D) sideband regions of $M_{\Lambda}$ versus $M_{\bar{\Lambda}}$ are adopted, and for mode II, one-dimensional (1D) sideband regions in the distribution of $M_{\pi}^{2}$ are used. The distributions of $M_{\Lambda(\bar{\Lambda})}$ and $M_{\pi}^{2}$ of inclusive MC samples after removing the channels containing the $\Lambda\bar{\Lambda}$ pair are nearly flat, so it is reasonable to use the sideband method. The 2D sideband regions (shown in Fig. 2) are chosen as: B1: $1.0901\leq M_{\Lambda}\leq 1.1029$ GeV$/c^{2}$ and $1.1285\leq M_{\bar{\Lambda}}\leq 1.1413$ GeV$/c^{2}$; B2: $1.1285\leq M_{\Lambda}\leq 1.1413$ GeV$/c^{2}$ and $1.1285\leq M_{\bar{\Lambda}}\leq 1.1413$ GeV$/c^{2}$; B3: $1.0901\leq M_{\Lambda}\leq 1.1029$ GeV$/c^{2}$ and $1.0901\leq M_{\bar{\Lambda}}\leq 1.1029$ GeV$/c^{2}$; and B4: $1.1285\leq M_{\Lambda}\leq 1.1413$ GeV$/c^{2}$ and $1.0901\leq M_{\bar{\Lambda}}\leq 1.1029$ GeV$/c^{2}$. The 1D sideband regions (shown in Fig. 3) are chosen as $-0.024\leq M_{\pi}^{2}\leq 0$ GeV${}^{2}/c^{4}$ and $0.029\leq M_{\pi}^{2}\leq 0.031$ GeV${}^{2}/c^{4}$. The numbers of events from sideband regions of data ($N_{\rm non-\Lambda\bar{\Lambda}}^{\rm data}$) are calculated by:

$$N_{\rm non-\Lambda\bar{\Lambda}}^{\rm data}=\frac{1}{4}\times N_{\rm 2D}+\frac{1}{2}\times N_{\rm 1D},$$

(5)

where $N_{\rm 2D}$ and $N_{\rm 1D}$ are the numbers of the events from the 2D and 1D sideband regions of data, respectively. The same sideband regions are used for the $\pi^{0}\Lambda\bar{\Lambda}$ and $\gamma(\Lambda\bar{\Sigma}^{0}+c.c.)$ MC samples, and the numbers of events from sideband regions of these MC samples ($N_{\rm non-\Lambda\bar{\Lambda}}^{\rm MC}$) are obtained with Eq. (5). The number of non-$\Lambda\bar{\Lambda}$ background events ($N_{\rm non-\Lambda\bar{\Lambda}}^{\rm bkg}$) is estimated by:

$$N_{\rm non-\Lambda\bar{\Lambda}}^{\rm bkg}=N_{\rm non-\Lambda\bar{\Lambda}}^{\rm data}-N_{\rm non-\Lambda\bar{\Lambda}}^{\rm MC}.$$

(6)

The numbers of events for the three main background channels above ($N_{\pi^{0}}^{\rm bkg}$, $N_{\Lambda\Sigma}^{\rm bkg}$, $N_{\rm non-\Lambda\bar{\Lambda}}^{\rm bkg}$) are calculated in each $M_{\Lambda\bar{\Lambda}}$ interval when measuring the Born cross section.

The distributions of $M_{\Lambda\bar{\Lambda}}$ of the main background events from all data sets are shown in Fig. 6, and the numbers of background events over all data sets for the three main background channels in each $M_{\Lambda\bar{\Lambda}}$ interval are listed in Table 2.

Several sources of systematic uncertainties are considered in the cross section measurement. The combined results of different reconstructed methods and different data sets are summarized in Tables 3 and 4 for the correlated and uncorrelated parts, respectively. The correlated and uncorrelated parts are summed in quadrature to determine the total uncertainty.

The integrated luminosity is measured with an uncertainty of 0.5% at $\sqrt{s}=3.773$ GeV and an uncertainty of 1.0% at other c.m. energies Ablikim et al. (2016b, 2015a, 2022c). In this analysis, the effective luminosity of the ISR process is calculated based on Eq. (3), and a 0.5% uncertainty is estimated Noh (2021). Thus, the total systematic uncertainty on the luminosity is 0.8% at $\sqrt{s}=3.773$ GeV and 1.2% at other energy points.

The uncertainties from the reconstruction of $\Lambda$ and $\bar{\Lambda}$ are studied by a control sample of $J/\psi\rightarrow\ pK^{-}\bar{\Lambda}+c.c.$, and determined to be 2.8% and 3.8% at $\sqrt{s}=3.773$ GeV, and 2.6% and 3.4% at other energy points, respectively. A 1.0% uncertainty is taken for the ISR photon detection Ablikim et al. (2019c).

For mode II, the uncertainties due to the $p(\bar{p})$ tracking and PID are 1.0% for each Ablikim et al. (2015b). The uncertainty due to the $M^{2}_{\pi^{-}}(M^{2}_{\pi^{+}})$ window is also studied by the control sample of $J/\psi\rightarrow\ pK^{-}\bar{\Lambda}+c.c.$, and estimated as 1.4% (0.8%) at $\sqrt{s}=3.773$ GeV, and 1.5% (0.9%) at other energy points. The uncertainty due to the branching fraction of $\Lambda\left(\bar{\Lambda}\right)\rightarrow p\pi^{-}\left(\bar{p}\pi^{+}\right)$, $\mathcal{B}\left(\Lambda\rightarrow p\pi\right)$, is obtained from the PDG Workman et al. (2022) to be 1.6%.

The uncertainty from the kinematic fit is divided into two parts: the contribution of the ISR photon and the contribution of the remainder. The former is determined by a control sample of the radiative Bhabha process $e^{+}e^{-}\rightarrow\gamma e^{+}e^{-}$, and estimated as 0.4%, 0.2% and 1.1% for the cases of full reconstruction, missing $\pi^{-}$ and missing $\pi^{+}$, respectively. The later is studied by a control sample of $J/\psi\rightarrow\Lambda\bar{\Lambda}$, and is 0.2% (0.2%), 2.4% (2.2%) and 2.2% (2.2%) at $\sqrt{s}=3.773$ GeV (other energy points), for the cases of full reconstruction, missing $\pi^{-}$ and missing $\pi^{+}$, respectively. Thus, the uncertainty due to the kinematic fit at $\sqrt{s}=3.773$ GeV (other energy points) is 0.6% (0.6%), 2.6% (2.4%) and 3.3% (3.3%) for the cases of full reconstruction, missing $\pi^{-}$ and missing $\pi^{+}$, respectively.

The signal MC samples are generated with PHSP. The angular distribution of the $\Lambda\bar{\Lambda}$ pair, the spin correlation between $\Lambda$ and $\bar{\Lambda}$, and the polarization of $\Lambda(\bar{\Lambda})$ decay are not taken into account. To estimate the uncertainty due to these factors, signal MC samples with an angular amplitude including these effects are generated. The parameterization of the angular amplitude is the same as that in Ref. Ablikim et al. (2019a), and the corresponding parameters are cited from it when $M_{\Lambda\bar{\Lambda}}\leq 2.52$ GeV$/c^{2}$ and obtained with the data set at $\sqrt{s}=2.900$ GeV when $M_{\Lambda\bar{\Lambda}}\geq 2.52$ GeV$/c^{2}$. The relative difference of the detection efficiency to that of the PHSP mode is regarded as the uncertainty.

The uncertainty from the MC model is considered by changing the event generator from ConExc Ping (2014) to PHOKHARA10.0 pho . The relative difference of the detection efficiency of these two event generators is taken as the uncertainty.

For the channel of $e^{+}e^{-}\rightarrow\pi^{0}\Lambda\bar{\Lambda}$, the sideband regions on the $M_{\gamma\gamma}$ spectrum are used to estimate the background of the $\pi^{0}\Lambda\bar{\Lambda}$ sample. Here, the 2D sideband regions (sideband of $M_{\Lambda}$ and $M_{\bar{\Lambda}}$) and 3D sideband regions (sideband of $M_{\Lambda}$, $M_{\bar{\Lambda}}$ and $M_{\gamma\gamma}$) are also used. The values of $\left|\frac{N_{M_{\gamma\gamma}}-N_{\rm 2D}}{N_{\rm sig}}\right|$ and $\left|\frac{N_{M_{\gamma\gamma}}-N_{\rm 3D}}{N_{\rm sig}}\right|$ are obtained, where $N_{\rm sig}$ is the number of signal events, $N_{M_{\gamma\gamma}}$, $N_{\rm 2D}$ and $N_{\rm 3D}$ are the estimated numbers of $\pi^{0}\Lambda\bar{\Lambda}$ events based on $M_{\gamma\gamma}$, 2D and 3D sidebands, respectively. The larger of the two values is taken as the uncertainty of this channel.

For the channel of $e^{+}e^{-}\rightarrow\gamma(\Lambda\bar{\Sigma}^{0}+c.c.)$, one of the parameters of the lineshape is changed by adding and subtracting a standard deviation ($\pm 1\sigma$). Based on the different lineshapes, different estimated numbers of $\gamma(\Lambda\bar{\Sigma}^{0}+c.c.)$ events are obtained. Further, the same method as for the $e^{+}e^{-}\rightarrow\pi^{0}\Lambda\bar{\Lambda}$ channel is used here to obtain the uncertainty of this channel.

For the non-$\Lambda\bar{\Lambda}$ background, we move the sideband regions by 0.002 GeV/$c^{2}$ and 0.002 GeV${}^{2}/c^{4}$ towards the signal for the 2D and the 1D sidebands, respectively, and obtain the new estimated numbers of non-$\Lambda\bar{\Lambda}$ background events. The relative difference between the old and new results is regarded as the uncertainty. For the $M_{\Lambda\bar{\Lambda}}$ interval of 2.70-2.80 GeV$/c^{2}$, since $N_{\rm sig}$ is extremely small ($0.8\pm 2.3$) at $\sqrt{s}=3.773$ GeV, the estimation of this uncertainty at $\sqrt{s}=3.773$ GeV is significantly larger than that in other intervals. Except for the three main background sources mentioned above, several other background channels are neglected, and their contribution is considered as a systematic uncertainty, which is 2.2% at $\sqrt{s}=3.773$ GeV and 1.1% at other energy points.

In this analysis, twelve data sets are used and three reconstruction methods (full reconstruction, and partial reconstruction with missing $\pi^{-}$ or $\pi^{+}$) are applied. We divide the data sets into two groups, where the first group only includes the data set at $\sqrt{s}=$ 3.773 GeV and the second group includes the other data sets at c.m. energies from 4.128 to 4.258 GeV. The uncertainties of the second group are studied together or inherited from the result at $\sqrt{s}=4.178$ GeV. Thus, the systematic uncertainties are combined in two steps, where the first step combines the three reconstruction methods in each group and the second step combines the two groups. Uncertainties of the three reconstruction methods (two data set groups) are combined as the average value weighted by detection efficiencies (products of detection efficiency and effective luminosity). The weighted average formula is:

$$\sigma_{\rm tot}^{2}=\sum_{i=1}^{3(2)}\omega_{i}^{2}\sigma_{i}^{2}+\sum_{i,j=1;i\neq j}^{3(2)}\rho_{ij}\omega_{i}\omega_{j}\sigma_{i}\sigma_{j},$$

(7)

with

$$\omega_{i}=\frac{\varepsilon_{i}}{\sum_{i=1}^{3}\varepsilon_{i}}~{}\left(\omega_{i}=\frac{\varepsilon_{i}\mathcal{L}_{i}}{\sum_{i=1}^{2}\varepsilon_{i}\mathcal{L}_{i}}\right),$$

(8)

where $\omega_{i}$, $\sigma_{i}$ and $\varepsilon_{i}$ with $i=1,2,3$ ($i=1,2$) are the weight, systematic uncertainty and efficiency for the reconstruction method (data set group) $i$, and $\rho_{ij}$ is the correlation parameter for two different reconstruction methods (data set groups) $i$ and $j$, and $\mathcal{L}_{i}$ is the effective luminosity for the data set group $i$. For the systematic uncertainties arising from background the $\rho_{ij}$ values are set to 0, and for other systematic uncertainties the $\rho_{ij}$ are set to 1.

The cross section for $e^{+}e^{-}\rightarrow\Lambda\bar{\Lambda}$ is calculated from the $M_{\Lambda\bar{\Lambda}}$ spectrum by:

$$\sigma_{\Lambda\bar{\Lambda}}\left(M_{\Lambda\bar{\Lambda}}\right)=\frac{\left(dN_{\rm sig}/dM_{\Lambda\bar{\Lambda}}\right)}{\varepsilon\cdot\mathcal{B}^{2}\left(\Lambda\rightarrow p\pi\right)\cdot d\mathcal{L}_{\rm int}/dM_{\Lambda\bar{\Lambda}}},$$

(9)

where $\left(dN_{\rm sig}/dM_{\Lambda\bar{\Lambda}}\right)$ is the $M_{\Lambda\bar{\Lambda}}$ spectrum of data corrected for resolution effects after subtracting the background, $\varepsilon$ is the detection efficiency from MC simulation as a function of $M_{\Lambda\bar{\Lambda}}$, and $\mathcal{B}\left(\Lambda\rightarrow p\pi\right)=0.639\pm 0.005$ Workman et al. (2022). The effective ISR luminosity $d\mathcal{L}_{\rm int}/dM_{\Lambda\bar{\Lambda}}$ is calculated by $d\mathcal{L}_{\rm int}/dM_{\Lambda\bar{\Lambda}}=W(s,x)\cdot\mathcal{L}_{\rm int}$, where $W(s,x)$ is described by Eq. (3). This effective luminosity includes the first-order radiative correction but does not take into account VP, so the obtained cross section is the “dressed” cross section.

The dependence of the mass resolution on $M_{\Lambda\bar{\Lambda}}$ is determined, and accordingly the $M_{\Lambda\bar{\Lambda}}$ is divided into 16 intervals from the threshold up to 3.00 GeV$/c^{2}$. To reduce the impact of the mass resolution, the width of the $M_{\Lambda\bar{\Lambda}}$ bin is at least 5 times larger than the mass resolution, so we do not correct the mass spectrum for resolution effects. The measured cross sections for the process $e^{+}e^{-}\rightarrow\Lambda\bar{\Lambda}$ in these intervals are listed in Table 5. A comparison between the results of this work and those of previous ones Ablikim et al. (2018a, 2019a); Bisello et al. (1990); Aubert et al. (2007) is displayed in Fig. 7.

A search for a threshold effect is made by performing a least chi-square fit to the cross section from the production threshold up to 3.00 GeV with different assumed functions. The systematic uncertainty is included in the fit with the correlated and uncorrelated parts considered separately.

The first fit function is a perturbative QCD (pQCD) driven energy power function Pacetti et al. (2015)

$$\sigma(s)=\frac{c_{0}\cdot\beta(s)\cdot C}{\left(\sqrt{s}-c_{1}\right)^{10}},$$

(10)

where $c_{0}$ and $c_{1}$ are free parameters and the Coulomb correction factor is $C=1$ for neutral baryons. The fit result is shown as the blue dashed line in Fig. 7, with $c_{0}=(1.07\pm 0.74)\times 10^{3}$ pb$\cdot$GeV¹⁰, $c_{1}=1.27\pm 0.08$ GeV and the fit quality $\chi^{2}/d.o.f=19.06/14$.