First study of the two-body scattering involving charm hadrons

This article presents the first measurement of the interaction between charm hadrons and nucleons. The two-particle momentum correlations of $\mathrm{pD^-}$ and $\mathrm{\overline{p}D}^+$ pairs are measured by the ALICE Collaboration in high-multiplicity pp collisions at $\sqrt{s} = 13~\mathrm{TeV}$. The data are compatible with the Coulomb-only interaction hypothesis within (1.1-1.5)$\sigma$. The level of agreement slightly improves if an attractive nucleon(N)$\overline{\mathrm{D}}$ strong interaction is considered, in contrast to most model predictions which suggest an overall repulsive interaction. This measurement allows for the first time an estimation of the 68% confidence level interval for the isospin $\mathrm{I}=0$ inverse scattering length of the $\mathrm{N\overline{D}}$ state ${f_{0,~\mathrm{I}=0}^{-1} \in [-0.4,0.9]~\mathrm{fm^{-1}}}$, assuming negligible interaction for the isospin $\mathrm{I}=1$ channel.


Introduction
The study of the residual strong interaction among hadrons is a very active field within nuclear physics. This interaction can lead to the formation of bound states, such as nuclei, or molecular states as, for example, the Λ(1405), which is considered as being generated from the attractive forces in the nucleon(N)K-Σπ channels [1][2][3][4]. One of the most fervent discussions in this context is nowadays revolving around systems involving charm mesons (D, D * ). Studies of their interaction are motivated by the observation of several new states with hidden charm and/or beauty (so-called XYZ states) [5][6][7][8][9], as well as with open charm such as the T cc + [10,11], and also of pentaquark states like P c (4380) and P c (4450) [12,13]. These exotic hadrons can be described as compact multiquark states in the context of the constituent-quark model [14], but are also considered as natural candidates for loosely bound molecular states [5,6]. For example, the structure of the χ c1 (3872) (formerly X(3872)) has been interpreted as a DD * /DD * molecular state or as a tetraquark [15]. Currently, definite conclusions are difficult to draw because of the lack of any direct experimental information on the DD * strong interaction. Strong support for the molecular nature of the Λ(1405) came not least from low-energy NK scattering data and information on the pK scattering length from kaonic hydrogen atoms [16][17][18][19]. Hence, a determination of the scattering parameters of systems involving D and/or D * mesons are pivotal to advance in the interpretation of the many observed states. The first step in this direction is the investigation of the interaction between the p(uud)D − (cd) pair and its charge conjugate. This interaction does not couple to the lower energy meson-baryon channels since no qq annihilation can occur. A measurement of this interaction is also an essential reference for the study of the in-medium D-and D * -meson properties [20]. Similarly to kaons and antikaons, it is theoretically predicted that possible modifications of the charm-meson spectral function at large baryonic densities can be connected to a decrease of the chiral condensate, thus providing sensitivity to chiral-symmetry restoration [21].
So far, the topic of the strong interaction between hadrons containing charm quarks was addressed only from a theoretical point of view [22][23][24][25] by employing different effective models anchored to the successful description of other baryon-meson final states, such as the NK and NK systems, while data are missing. Scattering experiments [26] and systematic studies of stable and unstable nuclei [27], accompanied by sophisticated calculations achieved within effective field theories [28,29], allowed us to reach a solid comprehension of the interaction among nucleons. When extending these studies to interactions including strange hadrons, the average properties of the interactions of some strange nucleon-hadron combinations (pK ± [30][31][32], pΛ, and pΣ 0 [33][34][35]) could be gauged with the help of scattering data and measurements of kaonic atoms [36]. The study of Λ hypernuclei [37] led to the extraction of an average attractive potential. The situation has drastically changed in recent years, thanks to the novel employment of the femtoscopy technique [38] in pp and p-Pb collisions at the LHC applied to almost all combinations of protons and strange hadrons [39]. The ALICE Collaboration could precisely study the following interactions: pp, pK ± , pΛ, pΛ, pΣ 0 , ΛΛ, ΛΛ, pΞ − , pΩ − , and pφ [39][40][41][42][43][44][45][46][47]. Since conventional scattering experiments cannot be performed with D mesons and charm nuclei [48] have not been discovered yet (searches for charm nuclear states are included in the scientific program of the Japan Proton Accelerator Research Complex [49]), the femtoscopy technique can be employed to study the ND and ND interactions. In this article, the first measurement of the strong interaction between a D − meson and a proton is reported. This pioneering analysis employs D − instead of the more abundantly produced D 0 mesons because of the smaller contribution from decays of excited charm states and the possibility to separate particles and antiparticles without ambiguity.

Experimental apparatus and data samples
The analysis was performed using a sample of high-multiplicity pp collisions at √ s = 13 TeV collected by ALICE [50,51] during the LHC Run 2 (2016-2018). The main detectors used for this analysis to reconstruct and identify the protons and the D-meson decay products are the Inner Tracking System (ITS) [52], the Time Projection Chamber (TPC) [53] and the Time-Of-Flight (TOF) detector [54]. They are located inside a large solenoidal magnet providing a uniform magnetic field of 0.5 T parallel to the LHC beam direction and cover the pseudorapidity interval |η| < 0.9. The events were recorded with a high-multiplicity trigger relying on the measured signal amplitudes in the V0 detector, which consists of two scintillator arrays covering the pseudorapidity intervals −3.7 < η < −1.7 and 2.8 < η < 5.1 [55]. The collected data sample corresponds to the 0.17% highest-multiplicity events out of all inelastic collisions with at least one charged particle in the pseudorapidity range |η| < 1 (denoted as INEL > 0). Events were further selected offline in order to remove machine-induced backgrounds [51]. The events were required to have a reconstructed collision vertex located within ±10 cm from the center of the detector along the beam-line direction to maintain a uniform acceptance. Events with multiple primary vertices (pileup), reconstructed from track segments measured with the two innermost ITS layers, were rejected. The remaining undetected pileup is of the order of 1% and therefore negligible in the analysis. After these selections, the analyzed data sample consists of about 10 9 events. The Monte Carlo (MC) samples used in this analysis consist of pp collisions simulated using the PYTHIA 8.243 event generator [56,57] with the Monash-13 tune [58] and GEANT3 [59] for the propagation of the generated particles through the detector. The proton candidates are selected according to the methods described in [39]. Charged-particle tracks reconstructed with the TPC are required to have transverse momentum 0.5 < p T < 4.05 GeV/c and pseudorapidity |η| < 0.8. Particle identification (PID) is conducted by measuring the specific energy loss and the time of flight with the TPC and TOF detectors, respectively. The selection is based on the deviation n σ between the measured and expected values for protons, normalized by the detector resolution σ . For proton candidates with a momentum p < 0.75 GeV/c, only the TPC is used by requiring |n TPC σ | < 3, while for larger momenta the PID information of TPC and TOF are combined and tracks are accepted only if the condition (n TPC σ ) 2 + (n TOF σ ) 2 < 3 is fulfilled. With these selection criteria, the purity of the proton sample averaged over p T is P p = 98% [39]. The contribution of secondary protons originating from weak decays or interactions with the detector material is assessed by using MC template fits to the measured distribution of the distance of closest approach of the track to the primary vertex. The estimated average fraction of primary protons is 86% [39].
The D ± mesons are reconstructed via their hadronic decay channel D ± → K ∓ π ± π ± , having a branching ratio BR = (9.38 ± 0.15)% [60]. D-meson candidates are defined combining triplets of tracks reconstructed in the TPC and ITS detectors with the proper charge signs, |η| < 0.8, p T > 0.3 GeV/c, and a minimum of two (out of six) hits in the ITS, with at least one in either of the two innermost layers to ensure a good pointing resolution. To reduce the large combinatorial background and the contribution of D ± mesons originating from beauty-hadron decays (non-prompt), a machine-learning multi-class classification algorithm based on Boosted Decision Trees (BDT) provided by the XGBOOST library [61,62] is employed. The variables utilized for the candidate selection in the BDT are based on the displaced decay-vertex topology, exploiting the mean proper decay length of D ± mesons of cτ ≈ 312 µm [60], and on the PID of charged pions and kaons. Before that, a preselection of the D ± candidates based on the PID information of the decay products is applied by requiring a 3σ compatibility either with the TPC or the TOF expected signals of the daughter tracks. Signal samples of prompt (originating from charm-quark hadronization or decays of excited charm states) and non-prompt D ± mesons for the BDT training are obtained from MC simulations. The background samples are obtained from the sidebands of the candidate invariant mass distributions in data. The BDT outputs are related to the candidate probability to be a prompt or non-prompt D meson, or combinatorial background. D-meson candidates are selected in the p T interval between 1 and 10 GeV/c by requiring a high probability to be a prompt D ± meson and a low   probability to be a combinatorial-background candidate.
A selection on the candidate invariant mass (M(Kππ)) is applied to obtain a high-purity sample of D ± mesons. To this end, the M(Kππ) distribution of D ± candidates is fitted in intervals of p T of 1 GeV width in the range 1 < p T < 10 GeV/c with a Gaussian function for the signal and an exponential term for the background. The left panel of Fig. 1 shows the M(Kππ) distribution for D ± with 2 < p T < 3 GeV/c. The width of the Gaussian function used to describe the signal peak, σ D ± , increases from 6 to 10 MeV/c 2 with increasing p T as a consequence of the p T dependence of the momentum resolution. The D ± -meson candidates in the invariant mass window |M(Kππ)| < 2σ D ± are selected to be paired with proton candidates. This selection, displayed by the two vertical lines in Fig. 1, leads to a purity that is P D − = (61.7 ± 0.9(stat) ± 0.7(syst))% on average. The systematic uncertainty of P D − is evaluated by repeating the invariant mass fits, varying the background fit function and the invariant mass upper and lower limits.
The contributions of prompt and non-prompt D ± mesons are depicted in the left panel of Fig. 1 with the red and blue distributions, respectively, They are obtained with a data-driven method based on the sampling of the raw yield at different values of the BDT output score related to the probability of being a non-prompt D ± meson [63]. The yields of prompt and non-prompt D ± mesons can be extracted by solving a system of equations that relate the raw yield value Y i (obtained with the i-th threshold on the BDT output score) to the corrected yields of prompt (N prompt ) and non-prompt (N non-prompt ) D ± mesons via the corresponding acceptance-times-efficiency factors for prompt ((Acc × ε) prompt i ) and non-prompt The δ i factors represent the residuals that account for the equations not holding exactly due to the The system of equations can be solved via a χ 2 minimization, which leads to the determination of N prompt and N non-prompt . The right panel of Fig. 1 shows an example of a raw-yield distribution as a function of the BDT-based selection used in the minimization procedure for D ± mesons with 2 < p T < 3 GeV/c. The leftmost data point of the distribution represents the raw yield corresponding to the loosest selection on the BDT output related to the candidate probability of being a non-prompt D ± meson, while the rightmost one corresponds to the strictest selection, which is expected to preferentially select non-prompt D ± mesons. The prompt and non-prompt components, obtained for each BDT-based selection using the procedure described above, are represented by the red and blue filled histograms, respectively, while their sum is reported by the magenta histogram.
The obtained corrected yields N prompt and N non-prompt can then be used to compute the fraction of nonprompt D ± mesons f j non-prompt for a given selection j, The f non-prompt factor for the selections used to build the pD − and pD + pairs is estimated to be (7.7 ± 0.5(stat) ± 0.2(syst))%. The systematic uncertainty of f non-prompt is evaluated by repeating the procedure with different sets of selection criteria and varying the fitting parameters in the raw-yield extraction. In addition, since the efficiency depends on the charged-particle multiplicity, the multiplicity distribution in the MC sample used for the efficiency computation was weighted in order to reproduce the one in data.
Differently from the component originating from beauty-hadron weak decays, D ± mesons originating from excited charm-meson strong decays cannot be experimentally resolved from promptly produced D ± mesons due to their short lifetime. The two largest sources are the D * ± → D ± π 0 and D * ± → D ± γ decays, having BR = (30.7 ± 0.5)% and BR = (1.6 ± 0.4)% [60], respectively. Their contribution is estimated from the production cross sections of D + and D * + mesons measured in pp collisions at √ s = 5.02 TeV [63,64] and employing the PYTHIA 8 decayer for the description of the D * ± → D ± X decay kinematics. The fraction of D ± mesons in 1 < p T < 10 GeV/c originating from D * ± decays is estimated to be f D * − = (27.6 ± 1.3(stat) ± 2.4(syst))%, where statistical and systematic uncertainties are propagated from the measurements of the D + and D * + production cross sections.

The correlation function
The proton and D − candidates are then combined and their relative momentum k * is evaluated as where p p p * p,D are the momenta of the two particles in the pair rest frame. The k * distribution of pD − pairs, N same (k * ), is then divided by the one obtained combining proton and D − candidates from different events, N mixed (k * ), to compute the two-particle momentum correlation function, which is defined as C exp (k * ) = N × N same (k * )/N mixed (k * ) [65]. The latter provides a correction for the acceptance of the detector and the normalization for the phase space of the particle pairs. To ensure the same geometrical acceptance as for N same , the mixing procedure is conducted only between particle pairs produced in events with similar z position of the primary vertex and similar charged-particle multiplicity. Since the correlation functions for pD − and pD + are consistent with each other within statistical uncertainties, they are combined and in the following pD − will represent pD − ⊕ pD + . The measured two-particle momentum correlation function can be related to the source function and the two-particle wave function via the Koonin-Pratt equation C(k * ) = d 3 r * S(r * )|Ψ(k * , r * )| 2 [65], where S(r * ) is the source function, Ψ(k * , r * ) is the two-particle wave function, and r * refers to the relative distance between the two particles. The source function for the pD − pairs is estimated by employing the hypothesis of a common source for all hadrons in high-multiplicity pp collisions at the LHC corrected for strong decays of extremely short-lived resonances (cτ 5 fm) feeding into the particle pairs [66]. This is the case for resonances strongly decaying into protons. In contrast, both beauty-hadron and D * ± decays occur at larger distances than the typical range for the strong interaction [60]. This implies that the correlation function for D − mesons originating from these decays will only carry the imprint of the interaction of the parent particle with the proton without impacting the size of the emitting source. The core source determined in [66] features a dependence on the transverse mass m T of the particle pair, which can be attributed to a collective expansion of the system [65,[67][68][69]. The collective behavior has been studied in high-multiplicity pp collisions by the CMS Collaboration and found to be comparable for light-flavor and charm hadrons [70]. Hence, the core source of pD − pairs with k * < 200 MeV/c is estimated by parameterizing the measured m T dependence of the source radius extracted from pp correlations in [66] and evaluating it at the m T = 2.7 GeV/c 2 of the pD − pairs. Since the production mechanism of charm mesons might not be identical to that of light-flavor baryons, the emission of the pp and pD − pairs is studied by simulating pp collisions with PYTHIA 8.301 [57] and computing their relative distance in the pair rest frame, r * , considering only D − mesons originating directly from charmquark hadronization. These studies indicate that the core source of pD − at the pertinent m T is smaller by about 25% compared to that of pp pairs. This is included in the systematic uncertainty of the source radius. The resulting overall source is parametrized by a Gaussian profile characterized by an effective radius R eff = 0.89 +0.08 −0.22 fm, where the uncertainty includes both the one arising from the m T -dependent parametrization and the PYTHIA 8 study.
The correlation function due to the genuine pD − interaction can be extracted from the measured C exp (k * ) by estimating and subtracting the contributions of D − mesons originating from beauty-hadron and D * − decays, protons originating from strange-hadron decays, as well as misidentified protons and combinatorial-background D-meson candidates. The experimental correlation function is decomposed as The combinatorial (K + π − π − ) background below the D − peak and the final-state interaction among protons and D − from D * − decays play a significant role. All other contributions are assumed to be characterized by a C(k * ) compatible with unity and are therefore included in the C flat contribution. The relative weights, λ i , are evaluated considering the contributions to D − candidates described above and following the procedure explained in [39] for the protons. They are about 33.9% for C pD − (k * ) and 38.8%, 14.4%, and 13.4% for the p(K + π − π − ), pD * − , and flat contributions, respectively.
The correlation function C p(K + π − π − ) is extracted from the sidebands of the D − candidates, chosen as for the left and right sidebands, respectively. The contamination from D * − → D 0 π − → K + π − π − decays in the right sideband is suppressed by a 2.5σ D * − rejection around the mean value of the D * − invariant mass peak. The resulting correlation function is parametrized by a third-order polynomial in k * ∈ [0, 1.5] GeV/c and is displayed by the green curve reported in the right panel of Fig. 2. The observed behavior is determined by meson-meson and baryon-meson mini-jets and residual two-body interactions among the quadruplet, as obtained from previous studies [42,46].
The residual pD * − correlation function is computed employing the Koonin-Pratt formalism using the CATS framework [71] to obtain a two-particle wave function Ψ(k * , r * ) considering only the Coulomb interaction and assuming that the source radius is the same as for pD − pairs. The obtained pD * − correlation function is transformed to the momentum basis of the pD − relative momentum by considering the kinematics of the D * − → D − X decay [72]. The resulting correlation function is shown in the right panel of Fig. 2 as a red band. The purple band in the same figure represents the total background that includes all contributions with their corresponding weights. Finally, the genuine pD − correlation function is obtained by solving Eq. 3 for C pD − (k * ) and is shown in Fig. 3.
The systematic uncertainties of the genuine pD − correlation function, C pD − (k * ), include (i) the uncertainties of C exp (k * ), (ii) the uncertainties of the λ i weights, and (iii) the uncertainties related to the parametrization of the background sources, C p(K + π − π − ) (k * ) and C pD * − (k * ). In particular, as previously mentioned, the systematic uncertainty on C exp (k * ) is estimated by varying the proton and D − -candidate selection criteria and ranges between 0.5% and 3% as a function of k * . The uncertainties of the λ i weights are derived from the systematic uncertainties on the proton and D − purities (P p and P D − ), f D * − , and f non-prompt reported in Section 3.1. The systematic uncertainties of C p(K + π − π − ) (k * ) are estimated following the same procedure adopted for C exp (k * ) and, in addition, by varying the range of the fit of the correlation function parametrized from the sidebands regions of the invariant mass distribution. Additional checks are performed by varying the invariant mass interval used to define the sidebands region of up to 100 MeV/c 2 . The resulting systematic uncertainty ranges from 1% to 5%. The systematic uncertainty of C pD * − (k * ) is due to the uncertainty on the emitting source. Considering the small λ pD * − (k * ) this uncertainty results to be negligible compared to the other sources of uncertainty. The overall relative systematic uncertainty on C pD − (k * ) resulting from the different sources ranges between 3% and 10% and is maximum in the lowest k * interval.

Results
The resulting genuine C pD − (k * ) correlation function can be employed to study the pD − strong interaction that is characterized by two isospin configurations and is coupled to the nD 0 channel. First of all, in order to assess the effect of the strong interaction on the correlation function, a reference calculation including only the Coulomb interaction is considered. The corresponding correlation function is obtained using CATS [71]. Secondly, various theoretical approaches to describe the strong interaction are benchmarked, including meson exchange (J. Haidenbauer et al. [22]), meson exchange based on heavy quark symmetry (Y. Yamaguchi et al. [25]), an SU(4) contact interaction (J. Hoffmann and M. Lutz [23]), and a chiral quark model (C. Fontoura et al. [24]). The relative wave functions for the model of J. Haidenbauer et al. [22] are provided directly, while for the other models [23][24][25] they are evaluated by employing a Gaussian potential whose strength is adjusted to describe the corresponding published I = 0 and I = 1 scattering lengths listed in Table 1. The pD − correlation function is computed within the Koonin-Pratt formalism, taking into account explicitly the coupling between the pD − and nD 0 channels [73] and including the Coulomb interaction [74]. The finite experimental momentum resolution is considered in the modeling of the correlation functions [39].
The outcome of these models is compared in Fig. 3 with the measured genuine pD − correlation function. The degree of consistency between data and models is quantified by the p-value computed in the range k * < 200 MeV/c. It is expressed by the number of standard deviations n σ reported in Table 1, where the n σ range accounts, at one standard deviation level, for the total uncertainties of the data points and the models. The values of the scattering lengths f 0 for the different models are also reported in Table 1.
Here, the high-energy physics convention on the scattering-length sign is adopted: a negative value corresponds to either a repulsive interaction or to an attractive one with presence of a bound state, while a positive value corresponds to an attractive interaction. The data are compatible with the Coulombonly hypothesis within Finally, the scattering parameters can be constrained by comparing the data with the outcome of calculations carried out varying the strength of the potential and the source radius. In this case the interaction potential is parametrized by a Gaussian-type functional form with the range of ρ-meson exchange. In this estimation, it is assumed that the interaction in the I = 1 channel is negligible for simplicity. The correlation function C pD − (k * ) is computed including also the Coulomb interaction and the coupled channel.  This procedure is repeated for different values of the interaction potential for the I = 0 channel (V I=0 ). For all the correlation functions corresponding to the different interaction potentials, the agreement with the data is evaluated by computing the χ 2 using a bootstrap procedure. Both the statistical and systematic uncertainties of the data are considered in the bootstrap procedure, as well as the uncertainty on the emitting source radius (R eff ) in the computed C pD − (k * ), which is varied within 1σ of its uncertainty. The resulting overall χ 2 distributions are shown in Fig. 4 as a function of f −1 0, I=0 and V I=0 in the left and right panels, respectively. The data are found to be consistent with a potential strength of V I=0 ∈ [−1450, −1050] MeV within 1σ . This corresponds to an inverse scattering-length interval of f −1 0, I=0 ∈ [−0.4, 0.9] fm −1 . Since the determined potential strength is always attractive, the positive values of the scattering length imply an attractive interaction without bound states, while the negative values are consistent with the presence of a ND bound state. The same procedure was repeated for fixed values of R eff in order to obtain the 1σ confidence interval as a function of the emitting source radius. Figure 5 shows the confidence interval as a function of the source radius varied within 1σ of its uncertainty. The dashed interval corresponds to the radius uncertainty due to only the m T dependence while the full-shaded interval shows the total radius uncertainty. The most probable value reported in Fig 5 with the star symbol corresponds to an attractive interaction with the formation of a bound state. Given that most models predict a repulsive I = 1 interaction, in reality the I = 0 interaction might have to be even more attractive. The herewith presented limits provide valuable guidance for further theoretical studies advancing the understanding of  the strong interaction in the charm sector.

Summary
In conclusion, this article presents the first measurement of correlation functions involving charm hadrons, which allows one to access to the strong interaction between a proton and a charm meson. The genuine pD − correlation function reflects the pattern of an overall attractive interaction. The data are compatible within (1.1-1.5)σ with the correlation function obtained from the hypothesis of a Coulomb-only interaction. The degree of consistency improves when considering, in addition, state-of-the-art models that predict an attractive strong ND interaction with or without a bound state. Finally, assuming no interaction for the I = 1 channel, the scattering length of the ND system in the isospin I = 0 channel is estimated as f −1 0, I=0 ∈ [−0.4, 0.9] fm −1 . This exploratory study paves the way for precision studies of the strong interactions involving charm hadrons, facilitated by about one order of magnitude larger pp data samples expected to be collected in the next years during the LHC Runs 3 and 4 [75].