Principles and prospects of Bose-Einstein correlation study at BESIII

The hadron production process in electron-positron collisions involves $ e^+ $ and $ e^- $ annihilating into a virtual photon ($ \gamma^\ast $), which splits into an initial quark-antiquark pair ($ q\bar{q} $). The hadronization process consists of three phases: (1) perturbative Quantum Chromodynamics (pQCD) cascade evolutions of quarks and gluons, (2) hadron source formation and preliminary hadron emission, and (3) the decay of unstable hadrons.

In the energy region of the Beijing Electron Positron Collider II (BEPCII), which ranges from 2 to 5 GeV, the preliminary hadron states include continuous states and resonant states characterized by the quantum numbers $J^{PC}=1^{--}$. Because the hadronization process cannot be predicted comprehensively and quantitatively by pQCD, phenomenological models are used to describe this phase. One of the models commonly used to accomplish this task is the semi-classical Lund string fragmentation model [1, 2]. The Lund model describes the hadron source as a string, which fragments into hadrons.

The production of vector resonant states, such as $ \rho $, $ \omega $, $ \phi $, $ J/\psi $, and $ \psi(3686) $, can be summarized as follows: the $ \gamma^\ast $ evolves into a $ q\bar{q} $ pair, which constitutes the vector meson. The meson then decays into final states, a process described by the vector meson dominance model [3].

In quantum mechanics, the propagation of hadrons can be described using wave functions. The wave function of a system of identical bosons is symmetric under the interchange of the space-time coordinates of any two bosons. This symmetry results in a higher probability of finding identical bosons in a small phase-space element than finding different bosons, as illustrated in Fig. 1. This phenomenon is known as Bose-Einstein correlation (BEC), which is independent of specific interactions. Measurements of the BEC effect provide valuable space-time information about the hadron source [4, 5].

Figure 1.  Identical bosons in hadron final states are more likely to be reassigned in closer phase space.

Theoretical studies of the BEC effect have been conducted for many years, leading to various models of the space-time distribution of the hadron source, including Gaussian- and ellipsoid-type distributions [4, 5], as well as the more generalized L$ \acute{e} $vy distribution [6−8].

Numerous collaborations have conducted experimental measurements of the BEC effect involving $ e^+e^- $, $ p\bar{p} $ collisions, and ep and Kp scattering to select various signal samples at energies ranging from tens of GeV to TeV. A selection of these studies are listed in Table 1, and additional details on other experiments are documented in Ref. [9]. All these experiments demonstrate the universal characteristics of the BEC effect.

Collaboration	Beam	$\sqrt{s}$ (GeV)	Boson pairs
MARKII [10]	$e^+e^-$	29	$\pi^\pm\pi^\pm$
AMY [11]	$e^+e^-$	58	$\pi^\pm\pi^\pm$
OPAL [12]	$e^+e^-$	91	$\pi^\pm\pi^\pm$, $\pi^0\pi^0$
L3 [13, 14]	$e^+e^-$	91	$\pi^0\pi^0$
NA22 [15]	$\pi^+p$,$K^+p$	21.7	$\pi^\pm\pi^\pm$
ZEUS [16]	ep	10.5	$\pi^\pm\pi^\pm$,$K^\pm K^\pm$
CMS [17]	pp	900, 2360	$\pi^\pm\pi^\pm$
ATLAS [18]	pp	13000	$\pi^\pm\pi^\pm$
ALICE [19]	pp	900, 7000	$\pi^+\pi^+$
PHENIX [7]	AuAu	200	$\pi^\pm\pi^\pm$

Table 1.  Status of measurements of the BEC effect.

However, in the BEPCII energy region, measurements of the BEC effect are still lacking. The BESIII detector at BEPCII has collected data samples with large-sample statistics at certain energies, providing the opportunity to measure the BEC effect. In this study, we first introduce the basic concepts of the BEC effect. Then we derive the emission amplitude of the hadron source based on the Lund model, as well as the wave functions of a two-pion system. Finally, we perform feasibility studies using the Monte Carlo (MC) method via the Lund model generator to explore the properties relevant for measuring the BEC effect.

The hadron source in the Lund model is not static, and it evolves in space-time according to the Lund area law [2]. Hadrons are emitted from convex vertices at various space-time points via string fragmentation. The mass and momentum of a hadron originate from the potential energy stored in the string. Figure 2 shows that meson i is produced at space-time point $ i=(x_i,t_i) $, and is composed of constituent quarks $ q_A $ and $ \bar{q}_B $ located at adjacent concave vertices A = $ (x_A,t_A) $ and B = $ (x_B,t_B) $. The mass and energy-momentum of meson i satisfy the relativistic mass-energy relation [1], which is expressed as

Figure 2.  Meson i with constituent quarks $q_A$ and $\bar{q}_B$ lying on hyperbolas of space-time.

where $ \kappa $ is the attractive tension constant between quarks and antiquarks.

Equation (1) indicates that points A and B lie on two hyperbolas, and that the space-time interval between A and B should be large enough to ensure that the produced hadron has the required mass. Consequently, the hadron source in the Lund model is discontinuous in space-time.

The Lund area law provides a strict solution for the Lund model. The probability of producing n hadrons is given by [1, 2]

where $ p_j $ is the momentum component of the hadron parallel to the string, and $ m_j $ is the transverse mass. The matrix element $ {\cal{M}}_n $ is expressed as

where $ {\cal{A}}_n $ is the light-cone area, and b is a dynamics constant.

The space-time $ \{x_i\} $ and energy-momentum $ \{p_i\} $ are related through Eq. (1), implying that these sets of quantities are equivalent. Thus, $ {\cal{M}}_n $ can be viewed as a function of either $ \{x_i\} $ or $ \{p_i\} $ equivalently, which can be expressed as

where the quantities in the first set of braces are the variables of $ {\cal{M}}_n $, and those in the second set of braces are the conjugate parameters. In principle, using the correspondence between $ \{x_i\} $ and $ \{p_i\} $ and the Jacobian transformation, $ d{\cal{P}}_n(\{p_j\}) $ in Eq. (2) can be rewritten as the conjugate expression $ d{\cal{P}}_n(\{x_j\}) $ in space-time coordinates, which is the space-time distribution of the hadron emission source.

To match experimental data, the Lund model incorporates phenomenological parameters that represent hadron production ratios at each vertex [2, 20]. Consequently, each matrix element $ {\cal{M}}_n $ is associated with a set of additional parameters that depend on the hadron type, which is denoted as the factor $ g_n $. The exclusive hadron emitting amplitude $ F_n(\{x_i\};\{p_i\}) $ can be expressed as

The emitting amplitude of a single particle produced at $ x_i $ with momentum $ p_i $ in an event of multiplicity n can be expressed as

which is a complex number and depends on the hadron category.

The quantity $ f_n(x_i;p_i) $ shown in Eq. (5) is impossible to calculate analytically. The single-particle emitting probability amplitude is numerically obtained from MC samples with a fixed multiplicity n via

where $N_{\rm evt}$ is the number of events in the MC sample. The inclusive single-hadron emitting probability amplitude in MC simulations corresponds to

In relativistic terms, a free and stable particle can be described by a plane wave function (strictly speaking, however, the motion of particles should be described by wave packets):

The lifetimes of certain hadrons, such as $ \pi^\pm $ and $ K^\pm $, are sufficiently long to be considered stable in experiments. For an unstable particle with a short life-time (resonance), the wave function can be written as

where m is the resonance peak position (the mass of the resonance), s is the square of the center-of-mass energy, and $ \Gamma $ is the decay width (which indicates the lifetime of the particle via $ \tau=\hbar/\Gamma $). For simplicity, we neglect the phase factor $ \exp(i\theta) $, as it does not affect the analysis of the BEC effect.

However, if a particle is influenced by Coulomb and strong interactions, the corresponding wave function becomes significantly more complex [21]. This study does not discuss the situation of final state interaction.

For the BESIII detector, the spatial coordinate origin O was defined as the center of the detector, which also served as the collision point of $ e^+e^- $, and the time origin was synchronized with the event trigger clock of the data acquisition system.

As shown in Fig. 3, two pions ($ \pi^+\pi^- $ and $ \pi^\pm\pi^\pm $) are emitted from the hadron source points at $ x_i=({\bf r_i},t_i) $ and $ x_j=({\bf r_j},t_j) $, respectively, and are detected by detector Da at $ x_a=({\bf r_a},t_a) $ and by detector Db at $ x_b=({\bf r_b},t_b) $. Here, $ {\it\bf r}_a $ and $ {\it\bf r}_b $ are the central positions of Da and Db relative to O, respectively, and $ t_a $ and $ t_b $ represent the times required for particles to travel from O to Da and Db, respectively.

Figure 3.  (color online) $\pi^+\pi^-$ and $\pi^+\pi^+$ are emitted from the hadron source and detected by detectors Da and Db.

For a real detector, Da is not a point but has a finite range, and time and momentum are measured with uncertainties. Thus, the detected interval $ x_{ia} $ and $ p_a $ can be expressed as

where $ \delta x_{ia} $ and $ \delta p_a $ denote the detected space-time and energy-momentum uncertainties. The wave function of $ \pi^+ $ arriving at $ x_a $ is given by

where $ \delta_{ia}=p_a\delta x_{ia}+\delta p_ax_{ia} $ is a random phase angle, possibly with observable effects on the measurement of the correlation functions.

Similarly, the wave function of $ \pi^- $ emitted from $ x_j $ with momentum $ p_j $ and detected by Db at $ x_b $ with momentum $ p_b $ is given by

where $ x_{jb}=x_b-x_j+\delta x_{jb} $, $ p_b=p_j+\delta p_b $, and $ \delta_{jb}= p_b\delta x_{jb}+ \delta p_bx_{jb} $.

The joint wave function of $ \pi^+\pi^- $ arriving at Da and Db is expressed as

This expression for $ \psi_d $ is also applicable to $ \pi^{\pm}\pi^0 $.

As depicted in Fig. 3(B), $ \pi^+ $ and $ \pi^+ $ are detected by Da and Db with momenta $ p_a $ and $ p_b $, respectively. The two identical pions are indistinguishable according to measurements, meaning Da and Db cannot determine whether a $ \pi^+ $ was emitted from $ x_i $ or $ x_j $. Consequently, the joint wave function of $ \pi^+\pi^+ $ has a symmetrized form for interchanging source points i and j:

The expression for $ \psi_s $ is also applicable to $ \pi^-\pi^- $ and $ \pi^0\pi^0 $.

If the preliminary state contains unstable particles that decay into stable final states, owing to the fact that the lifetimes of most unstable particles are only a few femtometers (much smaller than the spatial resolution of the detector), the decay vertex is considered as part of the experimental uncertainty in the emitting position.

The inclusive detection probability amplitudes for $ \pi^+ $ and $ \pi^- $, emitted from points $ x_i $ and $ x_j $ and detected by detectors Da and Db with momenta $ p_a $ and $ p_b $, respectively, are expressed as

In these calculations, we assume that $ p_a=p_i $ and $ p_b=p_j $ (incorporating measurement uncertainties into the random phase angles $ \delta_{ia} $ and $ \delta_{jb} $), and that they can be absorbed into the functions $ f(x_i;p_a) $ and $ f(x_j;p_b) $, respectively. Therefore, the detection probability amplitudes not only encode information about the hadron sources but are also influenced by measurement uncertainties.

The joint probability amplitude that $ \pi^+ $ and $ \pi^- $ are detected by Da and Db is given by

For the case in which two identical $ \pi^+\pi^+ $ are detected by Da and Db, four probability amplitudes exist, as shown in Fig. 3(B). Besides the two already given by Eqs. (16) and (17), the additional amplitudes are

The joint probability amplitude for detecting $ \pi^+\pi^+ $ is then

which is symmetric with respect to the interchange of source points i and j.

The expression for the correlation function for a two-pion system varies depending on whether the measurement is coherent or incoherent.

One of the key considerations in BEC effect measurements is determining whether the measurement is coherent or incoherent. These two types of measurements are represented by different formulas, leading to different results. In other words, the results of the measurement are influenced not only by the detected particles but also by the measurement method and the status of the detector. A more detailed explanation of how to distinguish between coherent and incoherent measurements can be found in Ref. [4].

The principles determining whether a measurement of the BEC effect is coherent or incoherent are analogous to those of the well-known electron double-slit experiments. If the experiment does not detect the slit through which the electron passes, the screen displays interference (coherent) patterns. Conversely, if the electron's path is detected, the screen shows non-interferential (incoherent) patterns. Similar results have been observed in double-slit experiments using neutrons and atoms, indicating the universal quantum characteristics of microscopic particles.

Based on the principles discovered in double-slit experiments, we can define the criteria for a BEC measurement to be coherent or incoherent as follows:

● Incoherent measurement: If the detector can determine the emission position and the trajectory of a particle, such as those found by the BESIII detector shown in Fig. 4(A), using a vertex detector (VD) or track reconstruction in the main drift chamber (MDC), the measurement is considered incoherent. In this scenario, the detection of particle paths eliminates interference effects.

Figure 4.  (color online) Illustrations of two kinds of detectors. Detector (A) can determine from where the particles were emitted, but detector (B) cannot.

● Coherent measurement: If the detector cannot determine the emission position and trajectory, as reflected by the setup shown in Fig. 4(B), and the VD and MDC are turned off and only the electromagnetic calorimeter (EMC) is operational, the measurement is considered coherent. In this case, the inability to track the particle paths allows interference patterns to emerge.

For a real detector that falls somewhere between the status of Fig. 4(A) and Fig. 4(B), the detected correlation functions are a combination of both incoherent and coherent terms. This mixed measurement results from partial information about particle trajectories, leading to partial interference effects.

When measurements are performed using data collected with an ideal and perfect detector, which is shown in Fig. 4(A), the measurement is guaranteed to be incoherent. This incoherence is a consequence of the inherent quantum characteristics of the system.

● Distribution of a single pion: When a hadron source emits a pion from position $ x_i $ with momentum $ p_i $, and it is detected by detector Da with momentum $ p_a $, the inclusive probability distribution is given by

which depends on the space-time distribution of only the hadron source.

● Distribution of two different pions: For a pair of $ \pi^+\pi^- $ emitted from position $ x_i $ with momentum $ p_i $ and position $ x_j $ with momentum $ p_j $, and detected by Da with momentum $ p_a $ and by Db with momentum $ p_b $, the incoherent correlation distribution is expressed as

This equation shows that the incoherent correlation distribution for different pion types is the product of individual single-pion probability distributions. The same expression also applies to pairs of $ \pi^\pm\pi^0 $.

● Distribution of two identical pions: For a pair of $ \pi^+\pi^+ $, emitted from position $ x_i $ with momentum $ p_i $ and position $ x_j $ with momentum $ p_j $, and detected by Da with momentum $ p_a $ and by Db with momentum $ p_b $, the incoherent correlation distribution is

where the additional term is defined as:

In Eq. (24), the first term represents the genuine incoherent correlation distribution, which is similar to $ P_{2d}^I(p_a,p_b) $ in Eq. (23). The second term, $ P_s^{I}(p_a,p_b) $, represents the interference, which accounts for phase averaging between the detection amplitudes of the two $ \pi^+ $ in $ {\cal{A}}_s $ given by Eq. (21). The interference term $ P_s^{I}(p_a,p_b) $ in Eq. (25) disappears if the phases $ \{\delta_{ia}\} $ and $ \{\delta_{jb}\} $ in the wave functions described by Eqs. (12) and (13) fluctuate randomly. The same expressions also apply to pairs of $ \pi^-\pi^- $ and $ \pi^0\pi^0 $.

When data samples are collected using the detector shown in Fig. 4(B), the measurements are coherent.

● Distribution of a single particle: The inclusive distribution of a pion emitted from $ x_i $ with momentum $ p_i $ and detected by Da with momentum $ p_a $ is given by

In this expression, the first term $ P_1^I(p_a) $, given by Eq. (22), represents the incoherent contribution, which depends on the intensity of individual sources, and the second term represents the coherent contribution, which reflects the interference between different source points. In the summation of these cross terms, each amplitude has a corresponding conjugate term, and their imaginary parts cancel each other out, leaving only the real parts.

● Correlation between two different pions: For a pair of $ \pi^+ $ and $ \pi^- $ detected by Da with momentum $ p_a $ and by Db with momentum $ p_b $, the coherent correlation distribution is

where

In Eq. (27), the first term is the incoherent correlation function given by Eq. (23), and the second term, $ P_d^{C}(p_a,p_b) $, is the true coherent one.

● Correlation between two identical pions: For a pair of identical $ \pi^+\pi^+ $ detected by Da with momentum $ p_a $ and by Db with momentum $ p_b $, the coherent correlation distribution is

where $ P_s^I(p_a,p_b) $ is given by Eq. (25), and

is the true coherent term of the symmetric wave function. The sign (positive or negative) of the summation of the cross terms is determined by the average of the random phases in Eqs. (12) and (13). This inherently quantum property produces experiment-dependent measurement results owing to the random phase angles caused by the experimental uncertainties.

Many factors contribute to correlation effects in the hadron final state. For instance, the $ \delta $ function in Eq. (2) imposes strong restrictions on the hadron distribution, leading to energy-momentum correlations among all final state particles. The correlation effect caused by energy-momentum conservation depends on the topology of the event multiplicity. To avoid bias, both the signal and reference samples should be selected from events with the same multiplicity.

We now focus exclusively on the BEC effect. The BEC function (BECF) in experiments is defined as the ratio of the correlation functions of identical bosons (the signal sample) to the correlation functions of different bosons (the reference sample). This ratio isolates the BEC effect by counteracting other correlation effects.

In an ideal scenario, if data samples of $ \pi^+\pi^+ $ and $ \pi^+\pi^- $ are collected using a perfect detector, as depicted in Fig. 4(A), this is termed an ideal measurement because the precise information of the hadron source is detected. The ratio of the corresponding incoherent correlation functions, based on Eqs. (23) and (24), is given by

where

In a more realistic scenario, if data samples of $ \pi^+\pi^+ $ and $ \pi^+\pi^- $ are collected using an imperfect detector, as shown in Fig. 4(B), this is termed a non-ideal measurement because no information about the hadron emission source is detected. In this case, the measurement is coherent, and the ratio of their correlation functions, based on Eqs. (27) and (29), is given by

where

Ideally, the data collected by the detector shown in Fig. 4(A) would be used to measure the BECF. However, real detectors are imperfect, often because of factors such as electromagnetic noise in the MDC or an incomplete reconstruction of the particle tracks, which indicates that certain emitting points of the hadron source are not detected. Thus, the measured BECF is expressed as

where $ P_{2s}(p_a,p_b) $ and $ P_{2d}(p_a,p_b) $ are combinations of incoherent and coherent correlation functions of the signal and reference samples, respectively:

The quantities $ \alpha $ and $ \beta $ reflect the degree of imperfection in the detector, which causes some of the hadron emissions to be measured in a coherent-like manner. Typically, $ \alpha $ and $ \beta $ are small, and the first-order approximation is given by

where $ R_2^I(p_a,p_b) $ is the ideal two-pion correlation function defined in Eq. (32), and

with

where $ \lambda^{CI} $ represents the deviation of the measured BECF from the ideal BECF owing to the imperfections of the real detector. The coefficient $ \lambda^{CI} $ is distinct from the incoherent parameter $ \lambda $ in phenomenological hadron source models and is adopted to fit the experimental data, which arises from different factors [4].

Currently, two main approaches for describing the origin of the BEC effect exist: the wave function approach and the more general field theory approach. The wave function approach is intuitive for explaining the BEC effect, as it assumes that the BECF can appear coherent or incoherent depending solely on the measurement method and experimental equipment, regardless of the coherence of the hadron source. In contrast, the field theory approach assumes that the hadron source is divided into incoherent and coherent field components [4−6].

The BEC effect is a natural consequence of the Lund model [1, 22]. In this model, multiple string configurations can fragment into the same hadron state. For example, as shown in Fig. 14.1 of Ref. [1], there are two configurations with light-cone areas $ {\cal{A}}_{12} $ and $ {\cal{A}}_{21} $, both of which fragment into the same hadron final state but with different momenta. In these configurations, two identical bosons are produced at convex vertices 1 and 2. The difference between these two configurations is simply the exchange of the two identical bosons. The total matrix element for this hadron state can be expressed as

leading to

where the term $ {\cal{H}} $ is defined as

where $ \Delta{\cal{A}} $ is the difference between $ {\cal{A}}_{12} $ and $ {\cal{A}}_{21} $. This shows that the interference between the two configurations $ {\cal{A}}_{12} $ and $ {\cal{A}}_{21} $ results in an enhancement factor of $ 1+{\cal{H}} $. This factor increases as $ \Delta{\cal{A}} $ decreases, which occurs when the momenta $ p_1 $ and $ p_2 $ of the two identical bosons become more similar in magnitude.

In quantum field theory, the distribution functions are related to the differential cross sections. For the inclusive process $ a+b\to c+X $ (where X represents all other particles except c), the differential cross section is defined as

where ${\rm d}\phi$ is the phase-space element of particle c, $ <n> $ is the average multiplicity, and $ \rho_1(p) $ is the single-particle density distribution. Similarly, for the two-particle differential cross section in the process $ a+b\to c_1+d_2+X $, we have

where $ \rho_2(p_1,p_2) $ is the two-particle density distribution.

The relationships between the dimensionless probability distributions and the density distributions are given by

In theoretical and modeling studies, the correlation function often uses the ratio

which is a six-dimensional representation. However, it is more convenient to use the one-dimensional two-particle correlation function, which is expressed as

where $ Q^2 $ is defined as the invariant squared difference of the four-momenta of the two particles:

The correlation functions are obtained via

where the pseudo-correlation in $ \rho_2(p_1,p_2) $ is deducted.

In experimental studies, BECFs are often parameterized as

where $ \lambda $ is known as the incoherent or chaotic parameter, and R is often interpreted as the source radius [4, 5]. The line shapes of Eqs. (53) and (54) are shown in Fig. 5. In the figure, $ \lambda=0.75 $ and $R = 0.4 $ fm are used, but the actual values are to be determined by fitting the experimental data. According to Eqs. (53) and (54), the correlation function has two important features: it goes to 1 at large Q values, and it goes to $ 1+\lambda $ at $Q=0 $.

Figure 5.  (color online) Line shapes of the exponential and Gaussian forms of the BECF.

Because $ Q^2 $ is defined as a four-dimensional form in Eq. (50), the parameter R in Eqs. (53) and (54) does not necessarily reflect the true spatial scale of the hadron source, as it also depends on the lifetime of the hadron source. Various models attempt to provide a physical explanation for the parameter $ \lambda $ [4−6, 23]. While $ \lambda $ may consist of multiple physical factors, the term expressed in Eq. (39) is certainly a component of them.

The shape of the BECF, as measured by the experimental data, can be fit using the correlation density, which is expressed as

where $ \Delta n(Q_i) $ is the number of particle pairs falling within the interval $ (Q_i,Q_i+\Delta Q_i) $, and K is a normalization factor.

The observed density correlation functions are obtained by selecting signal and reference particle pairs from the experimental data. Some preliminary hadrons are unstable and decay into final-state hadrons before reaching the detector. The correlation functions obtained from the decayed final-state hadrons are given by

The final state BECF is defined as

Ideally, the BECF should be measured using preliminary-state hadrons. If the decayed final states can be traced back to the preliminary states, the preliminary-state density correlation functions are given by

and the observed preliminary BECF is

Both the BECFs expressed by Eqs. (58) and (61) depend on the data analysis method, and the physical results account for the efficiency correction.

To derive the physical BECF, $ C^{\rm{phy}}(Q) $, from the observed $ C_{\rm{pre}}^{\rm{datobs}}(Q) $ and $ C_{\rm{fin}}^{\rm{datobs}}(Q) $, two sets of corrections are required: one that accounts for the efficiency in the analysis-dependent quantities, and another that converts the final-state $ C_{\rm{fin}}(Q) $ into a preliminary-state one $ C_{\rm{pre}}(Q) $.

If MC simulations align well with the experimental data, then the relationship

holds, where

are the BECFs obtained by MC simulations at the generator and detector levels, respectively.

For both preliminary- and final-state hadrons simulated at the generator and detector levels, we have:

The experimental physics results of the BECFs can be derived from the observed values by applying the efficiency corrections

where the efficiencies $ \epsilon_{\rm{pre}} $ and $ \epsilon_{\rm{fin}} $ are obtained via MC simulations:

The relationship between the experimental raw data and physics result is

where the recover factor is

The MC generators for the Lund model are LUARLW (for the low- to medium-energy region) [2] and JETSET (for the high-energy region) [20].

The LUBOEI subroutine in LUARLW and JETSET employs algorithms based on mean-field potential attraction between identical bosons to simulate the BEC effect. This is achieved by shifting the final-state momenta of identical bosons closer together while maintaining energy-momentum conservation within the event [24−26]. LUBOEI offers two options for simulating the BEC effect: the exponential form (given in Eq. (53)) and the Gaussian form (given in Eq. (54)). An analysis of the BESIII data indicated that the Gaussian form agreed well with the experimental data; thus, it was used in the MC simulations.

In this section, we present various BECF line shapes simulated at the generator level. The signal samples typically include $ \pi^+\pi^+ $, $ \pi^-\pi^- $, and $ \pi^0\pi^0 $, while the reference samples include $ \pi^+\pi^- $ and $ \pi^\pm\pi^0 $. We provide comparisons between the preliminary and final states for different signal and reference samples at the generator level.

In the MC simulations, we could generate either preliminary hadron states only or allow unstable hadrons to decay into final states. The inclusive line shapes of the normalized density distributions $ \rho_{\rm{presig}}^{\rm{simgen}}(Q) $ and $ \rho_{\rm{preref}}^{\rm{simgen}}(Q) $ in Eq. (65) and those of $ \rho_{\rm{finsig}}^{\rm{simgen}}(Q) $ and $ \rho_{\rm{finref}}^{\rm{simgen}}(Q) $ in Eq. (66) are shown in Fig. 6. The line shapes of $ \rho_{\rm{presig}}^{\rm{simgen}}(Q) $ and $ \rho_{\rm{finsig}}^{\rm{simgen}}(Q) $ tend towards the low Q end, which is a characteristic feature of the BEC effect.

Figure 6.  (color online) Density correlation function $\rho(Q)$ of the preliminary states (A) and final states (B) of the signal and reference samples at 3 GeV.

In the reference sample, some $ \pi^+\pi^- $ pairs are the decay final states of resonances such as $ \rho $, $ \eta $, $ \eta' $, $ \omega $, and $ \phi $. In Fig. 6(B), the convexity of $ \rho_{\rm{finref}}^{\rm{simgen}}(Q) $ is caused by resonant effects. We can exclude these $ \pi^+\pi^- $ pairs by requiring their squared invariant mass $ q^2=(p_{\pi^+}+p_{\pi^-})^2 $ to lie outside the resonance peak $[(M_R-\Gamma)^2,\;(M_R+\Gamma)^2]$. However, this requirement inevitably removes some $ \pi^+\pi^- $ pairs from non-resonant states.

Various schemes for selecting reference samples have been investigated by different collaborations [18]. However, these schemes may violate energy-momentum conservation. As noted in Ref. [4], the correlations arising from energy-momentum conservation are stronger than those arising from the BEC effect. It is expected that the biases introduced by these schemes can be corrected using high-quality MC simulations and the correction factor given in Eq. (70).

Figure 7 shows the density correlation functions $ \rho(Q) $ of the preliminary and final states for $ \pi^+\pi^- $, $ \pi^{\pm}\pi^0 $, and $ K^+K^- $. If $ \rho(Q) $ consists of only $ \pi^+\pi^- $ pairs or all of the $ \pi^+\pi^- $, $ \pi^{\pm}\pi^0 $, and $ K^+K^- $ pairs, the line shape of $ C(Q) $ differs.

Figure 7.  (color online) Density correlation functions of the preliminary and final states for the reference samples $\pi^+\pi^-$, $\pi^{\pm}\pi^0$, and $K^+K^-$ at 3 GeV.

Figure 8 illustrates the BECF of the inclusive preliminary and final states, along with the corresponding ratio as presented in Eq. (70). The figure indicates that the line shape of $ C(Q) $ is energy-dependent within the BEPCII energy range. Figure 9 shows the dependence of the BECF on the preliminary and final states for the semi-exclusive continuous states at 3 GeV.

Figure 8.  (color online) BECF of the inclusive preliminary and final states and the corresponding ratio at 2 GeV (top) and 3 GeV (bottom).

Figure 9.  (color online) BECF of the semi-exclusive continuous preliminary production and final states.

The $ J/\psi $ particle is produced via $ e^+e^-\to\gamma^\ast\to c\bar{c}\to J/\psi $, and it decays to preliminary states and/or stable final states. It is interesting to examine whether the BEC effect manifests in the preliminary and secondary final states of $ J/\psi $ decay. MC simulations indicate that both preliminary and secondary decay final states exhibit the BEC effect, as shown in Fig. 10, which shows that the line shapes of $ C(Q) $ and $ J/\psi $ are similar to those of the continuous states.

Figure 10.  (color online) BECF of the semi-exclusive preliminary and secondary decay final states of the three-gluon decay mode of $J/\psi$.

The Lund model predicts that the hadron emission amplitude with a fixed multiplicity given by Eq. (7) differs from that of the inclusive case given by Eq. (8). Figure 11 shows the BECF for preliminary hadron states with $ n= $ 4, 5, and 6, as well as for inclusive states.

Figure 11.  (color online) BECF of inclusive hadron states with different preliminary multiplicities.

Figure 12 demonstrates that the line shape of the BECF depends on the choice of signal and reference samples. This implies that the BECF of the hadron source should be measured exclusively, as inclusive measurements can only capture the average features of the hadron source.

Figure 12.  (color online) BECF of a selection of semi-exclusive channels. The reference sample includes $\pi^+\pi^-$ and $\pi^\pm\pi^0$.

Figure 13 shows the BECF for the preliminary and final states of $ \pi^\pm\pi^\pm $, $ \pi^0\pi^0 $, and $ K^{\pm}K^{\pm} $, $ \pi^\pm\pi^\pm $, $ \pi^0\pi^0 $, and $ K^{\pm}K^{\pm} $, indicating that the BECF depends on the choice of signal samples.

Figure 13.  (color online) BECF of the preliminary and final states of $\pi^\pm\pi^\pm$, $\pi^0\pi^0$, and $K^{\pm}K^{\pm}$.

Eqs. (53) and (54) are simplified parametric forms that describe the main characteristics of the BECF. In the MC simulations, LUBOEI achieves the BEC effect via the mean-field method, and the actual effects also depend on the specific hadron states. Compared with Eqs. (53) and (54), the practical BECFs are multiplied by the corresponding normalization factors, as is done in most experimental data analyses [12, 15], and thus the actual distributions shown in Figs. 8−13 appear different to those shown in Fig. 5, but the overall trends are similar.

In previous experiments, signal samples were selected using cyclic combinations, in which any boson could be combined with other identical ones multiple times. This method results in non-independent Q values, introducing extra correlations into $ \rho(Q) $. Figure 14 shows that an event with four identical pions can have six possible cyclic combinations.

Figure 14.  (color online) Four identical $\pi$ productions detected by Da and Db.

Figure 15 illustrates the line shapes of $ C(Q) $ obtained for circular and noncircular combinations of pion pairs. The former introduces extra correlation at small Q values, which is not due to the BEC effect. The correct approach to selecting signal samples is the latter, but the number of identical signal pions is lower.

Figure 15.  (color online) Comparison between the BECF line shapes for circular and noncircular combinations of $\pi^\pm\pi^\pm$.

An alternative generator, HYBRID [27], was also used in this study, and similar results were obtained. It combines PHOKARA [28], CONEXC [29], and LUARLW, which were developed for the R value measurement at BESIII [30].

The semi-classical Lund string fragmentation model is currently the only phenomenological model that can capture an accurate physical picture of the hadron production source and correlate the spatiotemporal distribution of the hadron source with the energy-momentum distribution of the produced hadrons.

Starting with the Lund model, we explored basic concepts and formulas related to the BEC effect and performed MC simulations. Through these feasibility studies, several properties of the BEC effect were identified:

(1) The Lund model establishes the relationships between the space-time points of the hadron source and the momentum of the emitted hadron.

(2) In addition to the coherent or incoherent components of the hadron source, the status of the particle detector also affects the measured values of the incoherence parameter $ \lambda $ and the source scale R.

(3) Different choices of signal and reference samples result in different line shapes for the BECF.

(4) Experimental measurements of the BEC effect should be conducted exclusively.

(5) The space-time distribution characteristics of the hadron source are energy-dependent.

(6) We demonstrated how to convert the observed final state BECF into the preliminary state one using MC simulations.

Table 2 presents the integrated luminosity and number of hadron events of the data samples collected using large-sample statistics at BESIII, which enable precise measurements of the BECF at these energies.

$ E_{\rm cm} (\rm GeV) $	$ {\cal{L} }_{\rm int.} ({\rm{pb} }^{-1}) $	$N_{\rm had}\, (10^6)$
2.125	108.5	5.89
2.396	66.9	2.85
2.645	67.8	2.35
2.900	106.0	3.06
3.080	126.0	3.09
J/$\psi$		$10^4$
$\psi$(2S)		3$\times10^3$

Table 2.  Potential data samples for BEC study at BESIII.

The measurement of the BECF using the data collected at BESIII is promising. The pion pairs $ \pi^\pm\pi^\pm $ and $ \pi^+\pi^- $ were selected as the signal and reference samples, respectively. The generator LUARLW was used to simulate the hadron production, while GEANT4 [31] was used to simulate the propagation of the particles and their interaction with the detector. The primary systematic uncertainty originated from the hadron generators, for which the difference between LUARLW and HYBRID was taken as an estimation.

A preliminary analysis of the BESIII data indicated that several million same-sign charged pion pairs are available at the five continuum energies. It is expected that the measurement precision for the hadron source parameters, obtained by fitting the line shape of $ C(Q) $, will be better than 6% for $ \lambda $ and better than 3% for R. In the future, data collected for the $ J/\psi $ and $ \psi(2S) $ peaks will be utilized to measure the BEC effect, further enhancing our understanding of the hadron source.

The authors thank Mr. Wilson J. Huang for proofreading the manuscript.