Revealing the nature of hidden charm pentaquarks with machine learning

In the last two decades, we have witnessed the emergence of tens of candidates for the so-called exotic hadrons which are beyond the conventional meson and baryon configurations, such as the hidden charm pentaquarks [1, 2, 3], the fully charmed tetraquark $X(6900)$ [4, 5, 6], or the doubly charmed tetraquark $T_{cc}^{+}$ [7, 8]. The first hidden charm pentaquarks were reported by the LHCb Collaboration in 2015 in Ref. [1] by observing the $J/\psi p$ invariant mass via the process $\Lambda_{b}\to J/\psi pK^{-}$. Four years later, the $P_{c}(4312)$, $P_{c}(4440)$ and $P_{c}(4457)$ were reported with an order of magnitude larger luminosity in Ref. [3]. Since then, many interpretations have been proposed for these states, e.g. compact multi-quark, hadronic molecule, hybrid, and triangle singularity [9, 10, 11, 12, 13, 14, 15, 16, 17, 16, 18, 19, 20]. Traditionally, judging whether a model works or not is to compare it with the experimental data. However, such a top-down approach makes the conclusion model-dependent. Physicists are thus trying to find bottom-up approaches [21] to obtain a definitive conclusion. One direction is to use machine learning to explore hadron properties, which benefits from the evolution of computing capabilities during the last decades. For instance, genetic algorithms [22, 23, 24] and neural networks [25, 26, 27, 28] can be utilized to explore the nature of hadrons or the properties of nuclei. Machine learning has been widely used in nuclear physics [29, 30, 31, 32, 33, 34, 35], high energy nuclear physics [36, 37, 38], experimental data analysis [39, 40], and theoretical physics [41, 42], even to discover the physical principles underneath the experimental data [21, 43, 44]. However, its application in hadron physics is in the early stage. A preliminary attempt was made in Refs. [45, 46, 21, 47] for the hadronic molecular picture [12], which is one of the natural explanations of near-threshold peaks. Refs. [48, 47] demonstrate that deep learning can be applied to classify poles and regress the model parameters in the one-channel coupling case.

Here, extending and improving the works [45, 46], we consider the multi-channel coupling case. We take the hidden charm pentaquarks, i.e. $P_{c}(4312)$, $P_{c}(4440)$ and $P_{c}(4457)$ [3] as examples, to achieve the following goals:

• (1)

  In the $\Sigma_{c}^{(*)}\bar{D}^{(*)}$ hadronic molecular picture [49, 50, 51], both the $P_{c}(4440)$ and the $P_{c}(4457)$ are related to the $\Sigma_{c}\bar{D}^{*}$ threshold. Their spin-parity ($J^{P}$) could be either $\frac{1}{2}^{-},~{}\frac{3}{2}^{-}$ (solution A) or $\frac{3}{2}^{-},~{}\frac{1}{2}^{-}$ (solution B) [49, 50, 51], from the viewpoint of pionless effective field theory (EFT). In solution A, the lower mass hidden charm pentaquark has lower spin, and vice versa in solution B. Although the $\chi^{2}$ fitting cannot distinguish the two solutions, we try to find out which solution is preferred in machine learning.

• (2)

  In the traditional fitting approach, the physics is embedded in the model parameters. One has to extract the model parameters from the experimental data and further extract the physical quantities of interest. We use a neural network (NN) based approach to extract the pertinent physical quantities directly and illustrate the role of each experimental data point (here the bins in the mass distributions) in the physics.

We have used one particular NN for the above goals and demonstrated that the NN-based approach consistently favors solution A over B compared to the fitting approach, and clearly gives the properties of the poles in the multi-channel case.

The $P_{c}(4312)$ and $P_{c}(4440)/P_{c}(4457)$ states are close to the $\Sigma_{c}\bar{D}$ and $\Sigma_{c}\bar{D}^{*}$ thresholds, respectively, making them prime candidates for hadronic molecules. In the heavy quark limit, these hidden charm pentaquarks are related to the scattering between the $(\Sigma_{c},\Sigma_{c}^{*})$ and the $(\bar{D},\bar{D}^{*})$ heavy quark spin doublets. Note that the heavy quark symmetry breaking effect is about $(m_{\rho}/m_{D})^{2}\sim 25\%$ in the charm sector, with $m_{\rho}$ and $m_{D}$ the masses of the $\rho$ and $D$ mesons. Since both the $P_{c}(4440)$ and $P_{c}(4457)$ are related to the $\Sigma_{c}\bar{D}^{*}$ channel, their mass difference comes from the spin-spin interaction between the $\Sigma_{c}$ baryon and the $\bar{D}^{*}$ meson for a given total spin, which should vanish in the exact heavy quark limit. Thus in the heavy quark limit the $P_{c}(4440)$ and $P_{c}(4457)$ should coincide with each other. Their hyperfine splitting stems from the heavy quark symmetry breaking effect. Consequently, we consider the heavy quark symmetry breaking effect from the mass difference in the same doublet. By constructing the interactions between the two doublets, i.e. the $(\Sigma_{c},\Sigma_{c}^{*})$ and the $(\bar{D},\bar{D}^{*})$ doublets, in the EFT respecting the heavy quark spin symmetry (HQSS), Refs. [51, 50] extracted the pole positions of the seven hidden charm pentaquarks by fitting the $J/\psi p$ invariant mass distribution for the process $\Lambda_{b}\to J/\psi pK^{-}$ ¹¹1Note that there is additional dynamic $\Lambda_{c}\bar{D}^{(*)}$ channel in Ref. [51] in comparison with Ref. [50]. However, the parameters of the $\Lambda_{c}\bar{D}^{(*)}$ channel are out-of-control due to the absence of the experimental data in the $\Lambda_{c}\bar{D}^{*}$ channel and the properties of the poles are driven by the $\Sigma_{c}^{(*)}\bar{D}^{(*)}$ channel. As the result, we only consider the elastic $\Sigma_{c}^{(*)}\bar{D}^{(*)}$ channel and the inelastic $\eta_{c}p$, $J/\psi p$ channels in our framework. . Following the same procedure, we consider the $\Sigma_{c}^{(*)}\bar{D}^{(*)}$ and $J/\psi p$, $\eta_{c}p$ channels as elastic and inelastic channels, respectively. The potential quantum numbers of hidden charm pentaquarks are $\frac{1}{2}^{-},~{}\frac{3}{2}^{-}$ and $\frac{5}{2}^{-}$ which correspond to the scattering of the $\Sigma_{c}\bar{D}-\Sigma_{c}\bar{D}^{*}-\Sigma_{c}^{*}\bar{D}^{*}$, $\Sigma_{c}^{*}\bar{D}-\Sigma_{c}\bar{D}^{*}-\Sigma_{c}^{*}\bar{D}^{*}$ and $\Sigma_{c}^{*}\bar{D}^{*}$, respectively. In the heavy quark limit, their interaction only depends on the light degrees of freedom [52]. To extract their interactions, one can expand the two-hadron basis in terms of the heavy-light basis $|H\otimes L\rangle$, with $H$ and $L$ the heavy and light degrees of freedom, respectively,

Here, the subscripts $\frac{1}{2}$, $\frac{3}{2}$ and $\frac{5}{2}$ give the total spin of the two-hadron system. In the HQSS, as the interaction only depends on the light degrees of freedom [52], the scattering amplitude is described by the two parameters [51]

with the subscripts corresponding to the spin of the light degrees of freedom. The coupling for the $S$-wave and $D$-wave inelastic channels, i.e. the $J/\psi p$ and $\eta_{c}p$ channels, are described by two parameters

respectively. The weak bare production amplitude can also be parameterized in terms of the $|H\otimes L\rangle$ basis

with $(H\otimes L)_{i}^{J}$ the $i$th state in the $|H\otimes L\rangle$ basis for spin $J$. With this definition, the bare decay amplitudes read

with the upper index $J=\frac{1}{2},\frac{3}{2},\frac{5}{2}$ for the total spin of the two-hadron system. Putting pieces together, the inelastic amplitude of the $J/\psi p$ for the decay process $\Lambda_{b}\to J/\psi pK^{-}$ can be expressed as [51, 50]

The corresponding diagrams are shown in Fig. 1a. Here, $U^{J}_{i}$ is the $i$th $J/\psi p$ inelastic channel amplitude with the quantum number $J$. $\mathcal{V}^{J}_{i\beta}$ is the transition vertex between the $\beta$th elastic and the $i$th inelastic channel, and described by the coupling constants $g_{S}$ for the $S$-wave and $g_{D}$ for the $D$-wave.

is the two-body propagator of the $\beta$th elastic channel [51, 50], with $\mu_{\beta}$ and $m_{\mathrm{th}}^{\beta}$ the reduced mass and threshold of the $\beta$th channel. The physical production amplitude of the $\beta$th elastic channel $U^{J}_{\beta}$ for total spin $J$ is obtained from the Lippmann-Schwinger (LSE)

with $P^{J}_{\alpha}$, constructed by seven parameters $\mathcal{F}^{J}_{n}$ defined as in Eq. (8), the production amplitude for the $\alpha$th elastic channel. The corresponding diagrams are shown in Fig. 1b.

Note that the above integral equations can be reduced to algebric equations in the framework of a pionless EFT. For a general introduction to this type of EFT, see [53]. $V^{J}$ is the effective potential

which contains both the scattering between the elastic channels $V_{C}$ described by two parameters $C_{1/2}$, $C_{3/2}$ and that $V^{J}_{\mathrm{in}}(E)$ between the elastic and inelastic channels described by the two parameters $g_{S}$, $g_{D}$. In total, eleven parameters are combined in a vector

To describe the $J/\psi p$ invariant mass distribution, a composite model is constructed via a probability distribution function (PDF) as:

with $U^{J}$ the $\Lambda_{b}\to J/\psi pK^{-}$ amplitude (14) for total spin $J=\frac{1}{2},~{}\frac{3}{2},~{}\frac{5}{2}$ and $\mathrm{p.s.}(E)$ is the phase space of the corresponding process. A Gaussian function $G(E^{\prime}-E)$ which represents the experimental detector resolution is convoluted with the physical invariant mass distribution. Here, we take the energy resolution as a constant of $\sigma\sim 2.3$ MeV [3] without considering its energy dependence. Further, ${\rm Chebyshev_{6}}$ is the sixth-order Chebyshev polynomial which represents the background distribution, with $1-\alpha$ its fraction, $\alpha\in[0,1]$. The background fraction in the data is determined to be (96.0$\pm$0.8)% as discussed in the Supplementary materials. The coefficients ($c_{0},...,c_{6}$) for the background component are obtained by fitting the Chebyshev polynomials to data [54], as listed in the Supplementary materials. Note that the area enclosed by the signal (background) distribution, i.e. the remnant of the first term after removing the factor $\alpha$ (the second term after removing the factor $1-\alpha$ ) in Eq. (17) is normalized to one.

The hidden charm pentaquarks are described as generated from the scattering of the $(\Sigma_{c},\Sigma_{c}^{*})$ and $(\bar{D},\bar{D}^{*})$ doublets. As the thresholds of the inelastic channels $J/\psi p$ and $\eta_{c}p$ are far away from those of elastic channels, their effects on the classification of the poles relevant to the physical observables are marginal. Thus, we consider the classification of poles by considering the elastic channels. In this case, the $\frac{1}{2}^{-}$ and $\frac{3}{2}^{-}$ states are given by a three-channel case with $2^{3}=8$ Riemann surfaces denoted as $\mathrm{R}_{\pm\pm\pm}$, see [20] for a pictorial. Here, the “+” and “$-$” signs in the $i$th position mean the physical and unphysical sheet of the $i$th channel, respectively. The physical observables, i.e. the $J/\psi p$ invariant mass distribution in our case, are real, reflecting the amplitude in the real axis. Those are related to the poles on the complex plane. We consider the poles most relevant to the physical observables, i.e. the poles on the physical sheet $\mathrm{R}_{+++}$ and those close-by ones. As we aim at extracting the most important information of those poles, it is sufficient to work with leading order contact interactions. In this case, the poles accounting for the near threshold structures can be classified as in Fig. 2. For the classification ²²2For the classification of poles, strictly, bound states do not decay to any final particles, which means them stable. However, in the standard model, there are only 11 stable particles, i.e. photon, the three generations of neutrino and antineutrino, electron, positron, proton, and antiproton. Other particles can decay into final states with lower masses. In our case, bound state, virtual state, and resonance are defined based on the most coupled elastic channels. of the states, we start from the one-channel case, i.e. $J^{P}=\frac{5}{2}^{-}$ state related to the $\Sigma_{c}^{*}\bar{D}^{*}$ channel. The “bound state”, “resonance”, and “virtual state” are defined as a pole on the physical sheet below threshold, a pole on the unphysical sheet above the threshold, and a pole on the unphysical sheet below the threshold, respectively (Fig. 2a). Those poles are labeled as $0$, $1$, $2$ in order. For the three-channel case (Fig. 2b), i.e. $J^{P}=\frac{1}{2}^{-}$ and $J^{P}=\frac{3}{2}^{-}$ channels, the poles are defined according to the most strongly coupled channel. The case with three poles on the $\mathrm{R}_{+++}$, $\mathrm{R}_{-++}$ and $\mathrm{R}_{--+}$ sheets slightly below the first, second, and third thresholds, respectively, is defined as a “bound state” which is labeled as $0$. The case with three poles on the $\mathrm{R}_{-++}$, $\mathrm{R}_{--+}$ and $\mathrm{R}_{---}$ sheets slightly above the first, second, and third thresholds, respectively, is defined as a “resonance” which is labeled as $1$. The case with three poles on the $\mathrm{R}_{-++}$, $\mathrm{R}_{--+}$ and $\mathrm{R}_{---}$ sheets slightly below the first, second, and third thresholds, respectively, is defined as a “virtual state” which is labeled as $2$. Note that these three states are defined according to their most strongly coupled channel, which the word “slightly” means. For instance, on the sheet $\mathrm{R}_{-++}$, both the “resonance” (blue circle point) and “bound state” (green crossed point) locate between the first and the second thresholds.

However, the former one strongly couples to the first channel and is located close to the first threshold, defined as a resonance of the first channel. The latter one strongly couples to the second channel and is located close to the second threshold, defined as a bound state of the second channel. The situation for the other states is similar.

We use a four-digit label $``olll"$, with $l=0,1,2$ and $o=0,1$, to denote the situation of the poles. The last three labels are for the nature of the poles in the $J^{P}=\frac{1}{2}^{-},~{}\frac{3}{2}^{-}$ and $\frac{5}{2}^{-}$ channels in order. The first label “$o$”, which we put in front, is used for the mass order. That is due to the indeterminacy of the quantum numbers of the $P_{c}(4440)$ and $P_{c}(4457)$. As they are close to the $\Sigma_{c}\bar{D}^{*}$ threshold, the $S$-wave interaction can give both $J^{P}=\frac{1}{2}^{-}$ and $\frac{3}{2}^{-}$ states. In the limit of HQSS, there are two solutions, i.e. Solution A and Solution B defined in Refs. [49, 50, 51], as discussed before. In Solution A, the higher spin state has larger mass, i.e. $J_{P_{c}(4440)}=\frac{1}{2}$ and $J_{P_{c}(4457)}=\frac{3}{2}$. In Solution B, the situation is reversed, i.e. $J_{P_{c}(4457)}=\frac{1}{2}$ and $J_{P_{c}(4440)}=\frac{3}{2}$. To distinguish these two scenarios, another label $o=1$ and $o=0$ for Solution A and Solution B, respectively, is required. In total, we have a four-digit label “$olll$” to denote the situation of the poles. Taking the “1012” case for example, it means that the poles for $J^{P}=\frac{1}{2}^{-},\frac{3}{2}^{-}$ and $\frac{5}{2}^{-}$ channels are “bound state”, “resonance” and “virtual state”, respectively. The first label $1$ means Solution A, i.e. $J_{P_{c}(4440)}=\frac{1}{2},~{}J_{P_{c}(4457)}=\frac{3}{2}$. Note that we can also have cases with the poles different from the three situations discussed above, e.g. the poles are far away from the thresholds. However, their probabilities are almost zero and are thus neglected.

Based on the PDF given in Eq. (17), 240184 samples corresponding to physical states, i.e. with poles not too far away from the real axis in the considered energy range, are selected among 1.85 million samples uniformly generated in the space of $\mathcal{P}$, distributed in the mass window from 4.25 $\mathrm{GeV}$ to 4.55 $\mathrm{GeV}$. These samples are represented as a set: $\{\mathcal{H}_{j},\mathcal{L}_{j}\}$ where the $j$ index refers to a given sample. A histogram $\mathcal{H}$ with 150 bins (see Fig. 3) denotes the invariant mass spectrum of the $P_{c}$ states (background included) and the label $\mathcal{L}$ indicates the state label defined in Sec. III.

The NN is implemented with an infrastructure of ResNet [56] as discussed in the Supplementary materials, and solved with the Adam [57] optimizer parametrized by the cross-entropy loss function. In order to balance the impact of unequal sample sizes, the cross-entropy loss function is weighted by the reciprocal of the number of samples. To avoid the network falling into a local minimum solution during the training process, the weights of the neurons are randomly initialized with a normal distribution centered at 0 and with a width of 0.01, and the biases of the neurons are set to zero. A reasonable solution can be obtained after training 500 epochs, using an initial learning rate value of 0.001.

Here the learning rate of the next interval is reduced to half of the previous one, controlled by the stepLR method. The interval is set to 100 epochs in our case. Note that 70% of MC simulation samples are used for training and the rest 30% for a benchmark. To avoid overfitting in the training, dropout layers are introduced in the network during the training. The dropout probability is chosen to 0.3 to optimise performance. In case of overfitting, the obtained NN works well on training datasets while probably does not on testing datasets. Thus the loss function calculated with training datasets should obviously differ from the one calculated with testing datasets. However, Fig. 9 in the Supplementary materials indicates that there is no obvious difference. The training details are summarized in the Supplementary materials. In short, the NNs successfully retrieve the state label with an accuracy (standard deviation) of 75.91(1.18)%, 73.14(1.05)%, 65.25(1.80)% and 54.35(2.32)% for MC simulation samples $\{\mathcal{S}^{90}\}$, $\{\mathcal{S}^{92}\}$, $\{\mathcal{S}^{94}\}$ and $\{\mathcal{S}^{96}\}$. The prediction accuracy decreases as the background increases, as expected.

Turning to the mass relation, we verified 100 $\{\mathcal{S}^{90}\}$ samples corresponding to the label “1xxx”, and another 100 samples corresponding to the label “0xxx”. Fig. 4a/b shows the predicted probability distribution for the label “1xxx/0xxx” samples, in which all labels are successfully retrieved. In other words, the NN can accurately predict the mass relationship. For comparison, we also perform an analysis with the $\chi^{2}$/ndf-fitting method. Fig. 4c/d shows the normalized $\chi^{2}$/ndf distributions corresponding to the label “1xxx/0xxx”. A 3% misidentification is observed for the label “1xxx” samples.

We now use the well-trained NNs to analyse the experimental data. The results are listed in Tab. 2. To reduce the systematic uncertainties arising from the NNs, we trained five and ten NN models, which have an identical structure under different initialization, for each group of samples with different background fractions. The probabilities of the sum of all the other labels are smaller than 1%. The labels with the top three probabilities are “1000”, “1001” and “1002”, with the sum of them larger than $99\%$. The standard deviation decreases with the increasing number of NN models, as expected. It means that our NNs can well determine the first three labels and favors Solution A, i.e. $J_{P_{c}(4440)}=\frac{1}{2}$ and $J_{P_{c}(4457)}=\frac{3}{2}$, as the first label is $1$.

That shows one of the advantages of machine learning comparing to the normal fitting approach. However, this conclusion is only based on the pionless EFT. One cannot obtain a solid conclusion about the hadronic molecules without the one-pion-exchange (OPE) potential. Thus it makes little sense to compare our results with that of Refs. [51, 50]. However, the hidden charm pentaquarks in pionless EFT are sufficient to illustrate the difference between machine learning and the normal fitting approach, as well as the potential merit of machine learning. The second label $0$ and the third label $0$ mean that the poles for the $J^{P}=\frac{1}{2}^{-}$ and $J^{P}=\frac{3}{2}^{-}$ channels behave as “bound states”, which is different from the virtual state conclusion for the $P_{c}(4312)$ in Ref. [21] also using machine learning. The pole situation for the $J^{P}=\frac{5}{2}^{-}$ is undetermined, i.e. all of the three labels $0$, $1$, $2$ appearing for the fourth label, as the structure around the $\Sigma_{c}^{*}\bar{D}^{*}$ is not significant. To resolve this issue, precise data in this energy region would be needed. In the normal fitting approach, the $J$=1/2, 3/2 and 5/2 channels are described by the same set of parameters in the heavy quark limit. As the result, once the set of parameters is determined, all the poles of the three channels are determined as well. However, in our NN approach, we extract the types of poles without extracting the intermediate parameters directly from the experimental data. This makes the poles of the $J$=1/2, 3/2 well determined, but the $J$=5/2 channel undetermined. Of importance is to understand how the NN learns information from the input mass spectrum $\mathcal{H}$. The SHapley Additive exPlanation (SHAP) value is widely used as a feature importance metric for well-established models in machine learning. One of the prevailing methods for estimating NN model features is the DeepLIFT algorithm based on DeepExplainer. A positive (negative) SHAP value indicates that a given data point is pushing the NN classification in favor of (against) a given class. A large absolute SHAP value implies a large impact of a given mass bin on the classification. This method requires a certain amount of dataset as a reference for evaluating network features. We choose 3000 {$\mathcal{S}^{90}$} samples in the dataset as a reference to test the SHAP values of both 1000 {$\mathcal{S}^{90}$} samples and the experimental data on our network. The upper panel of Fig. 18 shows the SHAP values of experimental data corresponding to 5 labels, i.e. “0000”, “1001”, “0100”, “1000” and “1002”. The lower panel of Fig. 18 shows the SHAP values of MC samples corresponding to the “1200” label, with others presented in the Supplementary materials. Both panels of Fig. 18 indicate that data points around the peaks in the mass spectrum have a greater impact.

A comparison is performed by checking the traditional $\chi^{2}$/ndf-fitting approach, which minimizes the

between the theoretical prediction $N(k;\mathcal{P})$ and experimental data $N^{d}(k)$ by summing over all bins of $\mathcal{H}$, where $\mathcal{P}$ represents the model parameters to be determined and $k$ is the bin index. We check how the $j$th bin of $\mathcal{H}$ impacts on the parameters $\mathcal{P}$ and the pole positions.

For this purpose, we define the uncertainty $\Delta U_{i}\equiv|\frac{\mathcal{P}_{i}(\mathrm{on})-\mathcal{P}_{i}(\mathrm{off})}{\mathcal{P}_{i}(\mathrm{on})}|$ for the $i$th parameter of $\mathcal{P}$, which is determined with the central values obtained by switching on/off the $j$th bin of $\mathcal{H}$ in the fitting. The bigger uncertainty observed by excluding a bin means a larger impact power of the bin on a parameter. Fig. 6 illustrates the $\Delta U_{i}$ distributions for two of the eleven parameters. The distributions for other parameters can be found in the Supplementary materials. The bins near the peaks in the mass spectrum do not show stronger constraints on the uncertainties of the parameters $C_{1/2}$ and $\mathcal{F}^{1/2}_{1}$ than the others. That is because that the model parameters are correlated to each other and do not reflect the underlying physics directly. In addition, to check the impact of the $j$th experimental point on the pole positions of the $J^{P}$ channel, the uncertainty

for the real part of the $i$th pole position is defined by switching on/off the $j$th bin of $\mathcal{H}$ in the fitting. A similar behavior can be found for the imaginary part of the poles, see Supplementary materials. Fig. 7 illustrates the $\Delta U_{i}^{J^{P}}$ distributions for the $J^{P}=\frac{1}{2}^{-},\frac{3}{2}^{-}$ and $\frac{5}{2}^{-}$ channels, respectively, without considering the correlation among the parameters. Although those values are of the order $10^{-5}$, one can still see that the experimental data around the $\Sigma_{c}\bar{D}$, $\Sigma_{c}\bar{D}^{*}$, and $\Sigma_{c}^{*}\bar{D}$ thresholds are more important than the others. The reason why the data around the $\Sigma_{c}^{*}\bar{D}^{*}$ threshold is not important is that the small production amplitude of the $\Sigma_{c}^{*}\bar{D}^{*}$ channel and the experimental data have little constraint about the corresponding poles [50, 51]. One can also obtain the same conclusion from Fig. 7g, where the experimental data around the $J^{P}=\frac{5}{2}^{-}$ dynamic channel $\Sigma_{c}^{*}\bar{D}^{*}$ do not show any significance. Due to the coupled-channel effect, the experimental data around the coupled channels still have strong constraints on the physics, e.g. the seven poles, around the other coupled channels. Taking the first pole of the $J^{P}=\frac{1}{2}^{-}$ channel as an example, e.g. Fig. 7a, the data around the $\Sigma_{c}\bar{D}^{*}$ threshold are as significant as those around the $\Sigma_{c}\bar{D}$ threshold. As the sample data of Fig. 7 corresponds to Solution A, the data around $P_{c}(4440)$ ($P_{c}(4457)$) are more important than those around the $P_{c}(4457)$ ($P_{c}(4440)$) in Fig. 7b (Fig. 7e) as expected.

We have investigated the nature of the famous hidden charm pentaquarks with a NN-based approach in a pionless EFT, which strongly favors $J_{P_{c}(4440)}=\frac{1}{2}$ and $J_{P_{c}(4457)}=\frac{3}{2}$, i.e. solution A in Refs [49, 50, 51, 58, 59]. From this NN, we find that solution A is systematically preferred over solution B. This conclusion is based on the pionless EFT, which is used to illustrate the difference between the machine learning and the normal fitting approach. Furthermore, we also performed checks on both the NN-based approach and the $\chi^{2}$/ndf-fitting approach. Our conclusion is that both approaches work well on the MC simulation samples. In the NN-based approach, the role of each data bin on the underlying physics is well reflected by the SHAP value. For the $\chi^{2}$/ndf-fitting approach, such a direct relation does not exist. This further explains why the two solutions can be better distinguished in the NN-based approach than in the $\chi^{2}$/ndf-fitting approach. At the same time, the effect of the background fraction on the accuracy of the network is accurately obtained, which provides a paradigm for NNs to study exotic hadrons by first extracting the background fraction from the data, and then analyzing the physics of the possible states from the data with the network trained from the corresponding background fraction data. This study provides more insights about how the NN-based approach predicts the nature of exotic states from the mass spectrum.

Acknowledgements:  We are grateful to Meng-Lin Du, Feng Kun Guo and Christoph Hanhart for the helpful discussion. This work is partly supported by the National Natural Science Foundation of China with Grant No. 12035007, Guangdong Provincial funding with Grant No. 2019QN01X172, Guangdong Major Project of Basic and Applied Basic Research No. 2020B0301030008. Q.W. and U.G.M. are also supported by the NSFC and the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through the funds provided to the Sino-German Collaborative Research Center TRR110 “Symmetries and the Emergence of Structure in QCD” (NSFC Grant No. 12070131001, DFG Project-ID 196253076-TRR 110). The work of U.G.M. is further supported by the Chinese Academy of Sciences (CAS) President’s International Fellowship Initiative (PIFI) (Grant No. 2018DM0034) and Volkswagen Stiftung (Grant No. 93562).

Author contributions:  Zhenyu Zhang and Jiahao Liu did the calculations. Jifeng Hu, Qian Wang and Ulf-G. Meißner focused on how to use machine learning to solve the problems in hadron physics. All the authors made substantial contributions to the physical and technical discussions as well as the editing of the manuscript. All the authors have read and approved the final version of the manuscript.