Production locality and spatial diffusion of heavy flavour quarks
at high energy densities

Collisions of heavy ions at relativistic energies provide insights into fascinating features of nuclear matter at high energy densities. This includes the creation and properties of the Quark-Gluon Plasma (QGP) [1] - a state of matter with quark and gluon degrees of freedom expected to exist in the Universe’s first moments. Moreover, there is a possibility of discovering the critical point of strongly interacting matter; for example, see Refs. [2, 3] and references therein. Impressive progress has been made in experimental and theoretical studies in the last decades. Still, many properties of high-density matter and particle creation in the medium remain to be uncovered. In particular, we propose testing the popular assumption of the local production of particles and antiparticles. This should be experimentally possible by studying heavy-quark (charm and bottom) correlations at properly selected collision energies.

Measurements of correlations between the charm meson and its anti-particle have been proposed to test the equilibration of charm [4, 5] in momentum space. In a semi-classical picture, the initial back-to-back momentum correlations between the $c$ and $\bar{c}$ quarks are reduced by the interactions with the medium and hadronization of the quarks (see, for instance, Ref. [6] and references therein). Thus, the charm hadron correlations provide means for quantifying transport properties of the strongly interacting medium, complementary to measurements of collective effects (via elliptic flow $\mathrm{v_{2}}$) [7] and modification of momentum spectra via nuclear modification factor $\mathrm{R_{AA}}$ [8].

In this paper, we present and discuss a different effect related to the creation mechanism of heavy quarks in a deconfined, partonic system. Using azimuthal correlations of charmed hadrons at low transverse momentum $p_{\textrm{T}}$, one can study the spatial distribution of charm quarks. Specifically, we show that the observed correlations at low $p_{\textrm{T}}$ provide direct insights into whether heavy quarks ($Q$ and $\bar{Q}$) are produced close to each other in the coordinate space (in other words, if there are produced locally), or if the production is non-local, i.e. production points of $Q$ and $\bar{Q}$ are distant. Clearly, the information on the creation locality is of fundamental importance. Moreover, the creation mechanism is an input assumption in models of charm and bottom quark interaction with the quark-gluon plasma, and all the modern experiments at SPS, RHIC and the LHC conduct programs that aim to quantify the QGP parameters using heavy quarks. The interpretation of the experimental observations (for instance, long-range correlations measured with a large difference in rapidity $\Delta y$) could drastically differ depending on whether one assumes the heavy quarks are initially created close or distant in the coordinate space. Thus, experimental verification of the usually assumed local production of a $Q\bar{Q}$ pair is crucial for understanding particle production in strong interactions. The study requires selecting heavy-ion collisions having at most a single $Q\bar{Q}$ pair created. This condition is fulfilled for charm quarks at the SPS energy and bottom quarks at RHIC at sufficiently low collision energies.

The idea utilizes the collective flow of charm hadrons measured in heavy-ion collisions at high energies [9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19]. We assume that final-state momenta of charm hadrons are given by the superposition of the flow and a random (thermal) contribution due to statistical hadronization. Hadronic rescattering and final-state interactions are neglected, supported by the recent measurements of interaction parameters of $D$-mesons with hadrons [20, 21]. The flow contribution depends only on the emission point in the freeze-out hypersurface, whereas the thermal contribution is a random effect, uncorrelated for different hadrons. This paper is motivated by the fact that, even for a locally thermalized and expanding medium, the momenta of charm and anti-charm hadrons originating from the same $c$ and $\bar{c}$ pair are correlated. Depending on the charm’s creation points and spatial diffusion properties in the medium, the charm hadron and the anti-hadron emission points can be either close or distant. In a locally thermalized and expanding system, the charm hadrons have an average momentum dependent on the fluid cell’s drift speed (flow). If the emission points of the hadrons are close, they will have a similar drift, and thus their momenta will be correlated. In the case of distant emission points the correlation should vanish.

Here, we stress that the arguments are generally valid, but to directly measure the wanted correlations, one should have no more than one $c\bar{c}$ pair produced per collision. Otherwise, the measured two-particle correlation function includes pairs of $c$- and $\bar{c}$-hadrons coming from different and likely independent charm production processes. The magnitude of this unwanted contribution to the momentum-correlation results strongly depends on the multiplicity distribution of heavy-flavour pairs. This effect is especially important in the heavy-ion collisions at RHIC and the LHC. On average, one expects $\simeq 3$ $c\bar{c}$ pairs in the 10% most central Au+Au collisions at $\mathrm{\sqrt{s_{NN}}=200}$ GeV at RHIC [22, 23], and a few tens at the LHC (for example, $\simeq 30$ $c\bar{c}$ pairs in the 10% most central Pb+Pb reactions at $\mathrm{\sqrt{s_{NN}}=5.02}$ TeV [24, 25]). The multiplicity distribution of heavy flavour is difficult to access experimentally, and the wanted correlations at very high energies cannot be extracted in a model-independent way. Thus, the measurements should be performed at sufficiently low collision energies to minimize the bias due to unwanted correlations, where the mean multiplicity of $c\bar{c}$ pairs is below one. For this reason, we consider an example of central Pb+Pb collisions at the CERN SPS energies. This example can be straightforwardly extended to bottom hadron production at RHIC or at the LHC.

The heavy-flavour production and azimuthal correlations in heavy-ion collisions at very high energies were addressed theoretically in the past; for review, see Ref. [26]. In particular, they were considered as a tool for uncovering a mechanism behind the jet suppression [27, 28] and the study of charm energy-loss mechanism [29, 30, 31, 32]. The heavy-quark spatial diffusion in QCD matter was discussed recently in Refs. [33, 34, 35], see also references therein. The ATLAS experiment measured the azimuthal-angle correlations of muon pairs originating from heavy-flavour decays in 5.02 TeV Pb+Pb collisions [36]. One notes that the measured muon pairs come from jet-like correlations of high transverse-momentum heavy-flavour hadrons. However, a possibility to distinguish between local and non-local $Q\bar{Q}$ production is addressed for the first time in this work.

The paper is organized as follows. First, we briefly discuss theoretical challenges in predicting ab initio hadron correlations in heavy-ion collisions, Sec. II. Then, the qualitative idea of testing production locality is quantified using simple modelling presented in Sec. III. The physics meaning of different assumptions on charm-hadron spatial correlations in the emission volume is also discussed in this section. The feasibility of the corresponding measurements is estimated in Sec. IV, and the results are summarized in Sec. V.

Quantum Chromodynamics (QCD) is the commonly accepted theory of strong interactions. However, attempts to derive precise quantitative predictions for multi-particle production in high-energy collisions from the QCD have been unsuccessful. Predictions of QCD–inspired models suffer from large uncertainties. Here, we discuss them in aspects relevant to this work.

The most popular QCD-inspired approaches to predict hadron production in heavy-ion collisions are based on classical approximations. For instance, these are the relativistic kinetic theory and hydrodynamic models [37]. In the context of heavy-ion collisions, the quantum effects are expected to be large or even dominant with respect to the classical predictions, at least regarding the flow of energy and momentum [38]. This is because the typical action scale of the system, a few hundred (at most) of MeV’s of temperature and spatial changes in a fraction of a femtometer, is smaller than $\hbar c\simeq 200$ MeV$\cdot$fm. The surprising success of hydrodynamical models in describing nuclear reactions can probably be traced back either to the attractor dynamics [39] or the generalized off-shell hydrodynamic expansion [40]. In any case, these arguments hold only for the hydrodynamic variables and not, for instance, for the two-particle correlations. Operators’ expectation values, such as the energy density in hydrodynamics, are generally considered. Still, their fluctuations (e.g., variance and higher-order moments) and related correlations are more difficult to deal with.

It is important to stress that, despite the success in predicting some observables, the approximations used in quantitative models are known to be inadequate for the full description of the data. For instance, hydrodynamics cannot properly address the high $p_{T}$ part of the spectra; relativistic kinetic theory assumes molecular chaos and removes two and three-particle correlations already at the classical level. Most physically proven quantum effects (diffraction, entanglement, etc.) are neglected and cannot be addressed by the current models. It is not simple to estimate the size of the quantum effects, lacking quantitative models including them. Thus, it is unclear whether the correlations produced by a classical treatment like relativistic kinematics are enough to describe the experimental results. For further discussion see Appendix A.

Considering the above, we follow a model-agnostic approach in this work. We do not consider questionable ab initio calculations of $Q$ and $\bar{Q}$ spatial and momentum correlations. Instead, we test whether extreme assumptions on heavy-quark emission locality lead to experimentally distinguishable predictions.

The following assumptions are made to quantify the intuitive expectations for the considered momentum correlation:

1. (i)

   The production of charm hadrons in head-on Pb+Pb collisions is considered. The collision energy is assumed to be adjusted to have a mean charm multiplicity below one, allowing to neglect production of more than one $c$- and $\bar{c}$ hadron pair in a single collision. This likely corresponds to the top CERN SPS energy ($\sqrt{s_{NN}}\approx 17$ GeV) [41].

2. (ii)

   The charm hadrons are emitted from the freeze-out hyper-surface of a spherical fireball undergoing a Hubble-like expansion. That is, the three velocity reads $\vec{v}=\vec{r}/t$, with $\vec{r}={x,y,z}$ being the distance from the centre of the fireball, and four-velocity $u^{\mu}=x^{\mu}/\tau=x^{\mu}/\sqrt{t^{2}-r^{2}}$. It was recently demonstrated that Hubble-like expansion is an appropriate approximation of velocity fields in heavy-ion collisions in the energy range of our interest [42].

3. (iii)

   The freeze-out hyper-surface is set by the freeze-out time $\tau=\tau_{fo}$ and the maximal radius $r\leq R_{max}$. They are assumed to be of the order of the Pb nucleus radius: $\tau_{fo}=9$ fm/$c$ and $R_{max}=6$ fm.

4. (iv)

   Emission probability of charm hadrons is independent of the fluid cell on the freeze-out hyper-surface. This is consistent with the method of predicting the spectra within the relativistic hydrodynamics approach. Note that the considered correlations are given by the conditional probability of the charm hadrons to be emitted from the same cell or another one with respect to the anti-charm one.

5. (v)

   In the rest frame of the flow, the charm hadron momentum $p$ distribution at the freeze-out hyper-surface is assumed to be the statistical one:

   $$\frac{d^{3}N}{dp\leavevmode\nobreak\ d^{2}\Omega}\leavevmode\nobreak\ \propto\leavevmode\nobreak\ p^{2}\exp{\left(\frac{-\sqrt{m^{2}+p^{2}}}{T_{\rm fo}}\right)}\leavevmode\nobreak\ ,$$

   (1)

   where $m=1.869$ $\mbox{Ge\kern-1.00006ptV}\!/\!c$ is the charm hadron mass assumed to be equal to the $D^{0}$ meson mass, and the temperature parameter is $T_{\rm fo}$ = 150 MeV. The statistical momenta of charm hadrons are drawn independently.

6. (vi)

   To calculate the hadron momentum in the collision rest frame, the obtained statistical four-momentum is boosted with the flow velocity, $\vec{v}=\vec{u}\leavevmode\nobreak\ /\leavevmode\nobreak\ \sqrt{1+u^{2}}$.

Note that for simplicity, we do not consider correlations between momenta $c$ and $\bar{c}$-hadrons resulting from the energy-momentum conservation and dynamics of the pair creation process. The change of these correlations during the system evolution was discussed in Refs. [4, 5] for heavy-ion collision at the top RHIC energy ($\sqrt{s_{NN}}=200$ GeV) and at the LHC. Given that we focus on producing charmed mesons with low $p_{\textrm{T}}$ in low-energy collisions, we expect the back-to-back correlation will not play a significant role in the measurement we consider in this work.

Then, the results are calculated for three different space correlations of the $c$- and $\bar{c}$-hadrons. These are

1. (a)

   The $c$- and $\bar{c}$-hadrons are emitted from the same fluid cell. Thus, the average of their momenta is set by the drift velocity of the cell. Their actual momenta are different because of the independence of their momenta in the fluid rest frame. This ansatz is labelled the local emission.

2. (b)

   The emission points of charm hadrons are independent of each other. Hence, they don’t have a common drift velocity, This ansatz is labelled the independent emission.

3. (c)

   The intermediate case is modelled, assuming the correlation function of the emission points to be the 3D Gaussian with $\sigma=\sigma_{x}=\sigma_{y}=\sigma_{z}=2$ fm. Note that the points are required to be within the fireball volume. The flow components of $c$- and $\bar{c}$-hadrons are different but correlated, leading to the correlation of their hadron momenta. Clearly in the limits of $\sigma\rightarrow 0$ and $\sigma\rightarrow\infty$ one recovers the local and independent emissions, respectively. This ansatz is labelled the correlated emission.

Figure 1 shows the distribution of $c$-$\bar{c}$ hadron pairs in the difference of azimuthal angles $\Delta\phi$ (left) and transverse momenta $\Delta p_{T}$ (right) for local, independent and correlated emission. The results are obtained using the Monte Carlo technique with $10^{7}$ events generated. The distributions of the pairs in $\Delta\phi$ significantly differ for local, independent and correlated emissions. The differences are smaller in the case of the transverse momentum difference. The uniform distribution in $\Delta\phi$ for the independent emission is independent of the flow and random momentum contributions modelling. The distributions in $\Delta\phi$ decrease monotonically for the local and correlated emission from $\Delta\phi=0$ to $\Delta\phi=\pi$, but the quantitative properties of this qualitative behaviour depend on model details. Nonetheless, the effect of correlation at $\Delta\phi\approx 0$ is remarkably different compared to the back-to-back correlations expected for charm pair production in hard parton scatterings [4, 5]. Thus, we expect experimental data will allow discrimination between these two different types of correlations.

It is clear that in the case of the $\Delta\phi$ distribution rather limited data statistics (see the next section) should allow us to distinguish between predictions obtained assuming different space correlations between the emitted charm hadrons and different production mechanisms of charm quarks. Encouraged by this conclusion, we turn to the standard approach to heavy-ion collisions [37] and, within it, discuss the implications of different possible outcomes of the experimental measurements.

The approach pictures heavy-ion collisions at high energies as a time sequence of the following stages:

1. (1)

   Initial stage - a high-density quark-gluon plasma is created. QCD is assumed to be a valid theory. Charm-anti-charm quark pairs are produced locally and in a limited number because of the high energy threshold.

2. (2)

   Expansion stage - the plasma expands [43], reaching the hadronisation temperature $T_{H}\approx 150$ MeV. The (anti-)charm quarks thermalize by interactions with the medium and acquire the medium flow.

3. (3)

   Hadronization stage - the plasma, including the $c$ and $\bar{c}$ quarks, is converted to hadrons and resonances following the statistical rules [44, 45] applied in the rest frame of a plasma fluid element. Thus, the flow and hadronization (local statistical process) contributions give the momenta of charm hadrons.

4. (4)

   Free-streaming stage - resonances decay, and non-interacting hadrons freely stream in the vacuum to a detector.

Many additional details, conceptual and quantitative, can be added [37], about the hydrodynamic evolution or the rescattering after hadronization. The stages listed above are the most relevant to this paper.

Within the standard heavy-ion approach,

1. (A)

   The experimental data consistent with the local emission would imply a small space separation of $c$- and $\bar{c}$-quarks during the expansion stage. This should be confronted with the charm-quark spatial diffusion calculated using the QCD-inspired approaches; for recent examples, see Refs. [33, 34, 35].

2. (B)

   The experimental data consistent with the independent emission would imply a large spatial diffusion of the charm quarks in the plasma. Ultimately, for more accurate models, the experimental results may be inconsistent with a local creation of heavy quarks and anti-quarks because of the limited expansion time and maximum velocity (speed of light) in a semi-classical transport. This might imply non-local creation or charm transport with velocities larger than the speed of light in a high-energy-density medium. The idea is qualitatively discussed in Appendix A.

   For the historical record, this paper was motivated by the non-local, indeterministic toy model [46] requiring the teleportation transitions in its most symmetric version.

3. (C)

   The data consistent with the correlated emission would give a sensitive input for restricting the charm-quark spatial diffusion in the plasma.

4. (D)

   It is always wise to leave a door open to the unexpected. The approximations used to compute the spectra might fail to describe the correlations, and the experimental results could qualitatively disagree with the expectations.

This section discusses the requirements for the experimental measurements of correlations between charm and anti-charm hadrons produced in head-on heavy-ion collisions. The important physics condition is a mean multiplicity of charm being small enough to neglect the production of two or more pairs of charm hadrons. This requirement implies the measurements at relatively small collision energies, probably close to the top SPS energy of $\sqrt{s_{NN}}\approx 20$ GeV. It also suggests collecting data in the fixed target mode, which, due to the Lorentz boost of the centre-of-mass allows for high detection acceptance and efficiency. For now, we only consider measurements of the most abundant open-charm hadrons, $D^{0}$ and $\bar{D^{0}}$ mesons. The required statistics of recorded Pb+Pb central collisions can be derived from the average number of reconstructed $D^{0}\bar{D^{0}}$-pairs, $\langle D^{0}\bar{D^{0}}\rangle_{\textnormal{rec}}$. In Appendix B, we estimate that modern experiments at the CERN SPS should be able to record sufficient data to measure 1000 or more $D^{0}\bar{D^{0}}$-pairs. Figure 2 demonstrates the statistical precision of a signal from 1000 $D^{0}\bar{D^{0}}$-pairs, assuming that the statistical fluctuations of background pairs can be neglected. We conclude that this is sufficient to distinguish between the independent and local emission.

The discussed measurement of correlations between $c$ and $\bar{c}$ is only meaningful for the quarks produced as a pair. However, if multiple pairs of $c\bar{c}$ quarks were produced within the same event, we will observe an unavoidable background due to combining $c$- and $\bar{c}$-hadrons originating from different pairs. Quantifying this background suffices taking a ratio between a multiplicity of produced $c\bar{c}$ pairs in the event to a number of all possible combinations of $c$- and $\bar{c}$-hadrons that could be observed in the event. To compute this ratio for different values of $\langle c\bar{c}\rangle$, it was assumed that $c\bar{c}$-multiplicity follows a Poisson distribution, parameterized by the given $\langle c\bar{c}\rangle$. This yields a probability of having more than one $c\bar{c}$ pair within the same event, in which case one has unwanted combinations of $c$- and $\bar{c}$-hadrons. Obtained values of the ratio pairs to combinations are given in Appendix B, Table 1. Values around $50\%$ indicate that about half of the combinations of $c$- and $\bar{c}$-hadrons are unwanted. For the real-world analysis, background due to the misidentification of open-charm hadrons would likely be significant. Moreover, lower values of $\langle c\bar{c}\rangle$ also imply higher requirements for the event statistics. It suggests that $\langle c\bar{c}\rangle$ should be selected to maximize the results’ significance. Realistic estimates of $\langle c\bar{c}\rangle$ for head-on Pb+Pb collisions at the top CERN SPS energies range between $0.1$ and $1$ [41]. This further supports the conclusion that the measurements at the CERN SPS should considered.

In this work, we propose to test the locality of heavy-quark creation by studying the momentum correlation of $c$- and $\bar{c}$-hadrons produced in heavy-ion collisions at collision energies with the mean multiplicity of $c\bar{c}$ pairs below one.

We show that the azimuthal correlations of charm hadrons observed in an experiment are sensitive to whether $c$ and $\bar{c}$ quarks are emitted locally, close to each other in the coordinate space, or they have distant, ultimately uncorrelated, emission points. Beyond its fundamental significance, the locality test will provide important input for modelling charm quark interaction with a nuclear medium.

Since the emission of multiple, uncorrelated, pairs of $c$- and $\bar{c}$-hadrons from a single collision would spoil the wanted correlations, it is mandatory to perform the measurements at sufficiently low collision energies, granting a low production probability of multiple-charm pairs. The proposed method can also be used for bottom hadrons.

As a quantitative example, we consider charm hadron measurements in head-on Pb+Pb collisions at the CERN SPS. Assuming typical values of data-taking parameters for the NA61/SHINE experiment at SPS, we show that the required measurements would need a data-taking rate of 10k Hz or more. These rates are easily allowed by the current detector technologies. Thus, the corresponding measurements should be possible by the upgraded NA61/SHINE and the new NA60++ experiments after the CERN LS3 upgrade period.

Appendices

A. Classical vs quantum-mechanical approach to heavy-quark production. Here, we discuss why the classical approximation for charm production in the limit of a single pair is incorrect. The simplest model to address is hydrodynamics. The main equation is the local four-momentum conservation or, more precisely, the expectation value of it

with respect to the density matrix $\hat{\rho}$ of the system. Additional equations, the equation of state and the treatment of the non-ideal part (transport coefficients) allow us to solve the system for the expectation values $T^{\mu\nu}$. Statistical hadronization is then used to calculate predictions for particle production. Scattering after the hadronization is usually considered with a separate transport phase. In some cases, the baryon number conservation equation

is added to the hydrodynamics equations. In principle, the electric current and the other conserved charges should also be considered when calculating the charge densities. Moreover, the link between tensors in space-time and particles in phase space, necessary at the hadronization stage, is through the relativistic Wigner distribution $W(x,p)$ [47]. The latter is the generalization of the classical distribution function $W(x,p)\xrightarrow{\rm classical\;limit}\propto\delta(p^{2}-m^{2})f(x,{\bf p})$. It depends on the bilinearity of the fields and the one-particle reduced density matrix. It does not depend on the two-particle ones and higher orders. By construction, regardless of the ansatz (local equilibrium, viscous corrections), the hadronization formula is for the one-particle observables only. All of the content about particle correlations must come from somewhere else.
Despite being a very different model, similar considerations hold for the relativistic Boltzmann equation because it stems from $W(x,p)$ too. Relativistic kinetic theory is a limit of the evolution of the Wigner distribution. As explained in [47], the approximations needed to use the Relativistic Boltzmann equation instead of the more general equations for the evolution of $W(x,p)$ include both arbitrarily small gradients and arbitrarily weak interaction. Then, one can neglect the coupling with the two-particle reduced density matrix, and the only “quantum leftover” is the cross-section, which must be evaluated in the framework of axiomatic field theory. These two conditions are enough to question whether relativistic kinetic theory can be used for the QGP. Strong interactions and large gradients are needed to fit the experimental data on top of a realistic (non-ideal) state equation that already requires phenomenological modifications to the simple relativistic Boltzmann equation. All the phenomenological modifications used in the state-of-the-art models (temperature-dependent masses, off-shell cross sections, etc.) do not insert any contribution from the $n$-particles reduced density matrix. This sector of the microscopic theory is systematically neglected. The spectra can be deduced from $W(x,p)$ alone, a one-particle object. If the evolution of $W(x,p)$ couples mostly to itself, one can argue that these extensions of the relativistic kinetic theory have a good chance to reproduce the spectra (and $v_{2}$ and other one particle objects). Still, the same cannot be said about correlations.
If one prefers a more intuitive approach to quantum fields, some considerations must be made from first principles. Because of the Heisenberg uncertainty, one cannot have an arbitrarily sharp wave function in both position and momentum at the same time. The more $c\bar{c}$ pair is well-defined in momentum, the more it must be delocalized. If the quarks are assumed to be produced as close as possible to momentum eigenstates, to forget about the details of the wavefunction in momentum space, one has to consider them substantially delocalized in space. They cannot be considered in one cell, and the wave function in the configuration space gives a weight regarding which part of the medium is ‘’felt more” by the heavy flavours. In any case, neither hydrodynamics nor kinetic theory are equipped to treat such wave functions dynamically.

B. Example estimate of event statistics and data-taking time. Here, we present a simple estimate of the event statistics and data-taking time assuming detector setup and performance similar to the NA61/SHINE experiment at CERN [48] recording Pb+Pb collisions at $\sqrt{s_{NN}}=17.3$ GeV. Assuming that processes that impact the reconstruction of a $D^{0}$ and a $\bar{D^{0}}$ mesons within an event are approximately uncorrelated, we estimate the average number of reconstructed pairs as

where $\langle c\bar{c}\rangle$ is the average number of $c\bar{c}$-pairs per event. The $P(c\rightarrow D^{0})$ = 0.31 is a probability for $c$-quark to hadronize into the $D^{0}$ meson evaluated within the PHSD model [49], BR$(D^{0}\rightarrow K^{+}\pi^{-})$ = 3.98% is a branching ratio of decay channel used in the measurements [50], $P(\textnormal{acc})$ = 0.5 is a probability for $D^{0}$ to be within an acceptance region of the detector, $P(\textnormal{sel})$ = 0.2 is a probability for $D^{0}$ to pass background-suppressing selection of charm meson candidates, and $P(\textnormal{rec})$ = 0.9 is a probability of reconstructing the meson. The value of $P(\textnormal{acc})$ was evaluated using the Geant4 simulation with the detector setup for November 2022, $P(\textnormal{sel})$ is taken from the pilot analysis of $D^{0}$ and $\bar{D^{0}}$ production [51], and $P(\textnormal{rec})$ was obtained from a Geant4 simulation with the setup for November 2022 and reconstruction software used for previous open charm analysis using 2017 and 2018 data  [52, 51].

Finally, given $\langle D^{0}\bar{D^{0}}\rangle_{rec}$, an estimate of the required event statistics can be obtained via

The $\langle c\bar{c}\rangle$ value is neither reliably predicted by models nor measured by experiments. However, considering available estimates [41], we expect that the value of $\langle c\bar{c}\rangle$ for head-on Pb+Pb at $\sqrt{s_{NN}}\approx 17$ GeV should range from 0.1 up to 1.

Putting all together, estimates on the run time needed to collect 1000 $D^{0}\bar{D^{0}}$-pairs for different event rates of the upgraded NA61/SHINE experiment and for different values of $\langle c\bar{c}\rangle$ are given in Table 1.

A typical ion beam period at CERN is about four weeks. Entries in Table 1 with a data-taking time of 100 days or more correspond to scenarios where the measurement may take longer than a period between the CERN accelerators’ long shutdowns. Moreover, at the moment, the event rate of 100 kHz would require a significant upgrade of the NA61/SHINE detector and its beamline. However, a setup corresponding to 10 kHz should be achievable within the nearest years. Thus, we find that for $\langle c\bar{c}\rangle>0.2$, it should be possible to perform the measurements of $c\bar{c}$-correlations by NA61/SHINE in the CERN Run 4 period (2028-2032). The additional possibility for the experimental study would be constructing a new experiment optimized for charm measurements. The corresponding letter of intent was recently submitted to the CERN SPSC [53].