Spectral-statistics properties of the experimental and theoretical light meson spectra

We present a robust analysis of the spectral fluctuations exhibited by the light meson spectrum. This analysis provides information about the degree of chaos in light mesons and may be useful to get some insight on the underlying interactions. Our analysis unveils that the statistical properties of the light meson spectrum are close, but not exactly equal, to those of chaotic systems. In addition, we have analyzed several theoretical spectra including the latest lattice QCD calculation. With a single exception, their statistical properties are close to those of a generic integrable system, and thus incompatible with the experimental spectrum.


Introduction
Since Yukawa first proposed the pion to explain the internucleon force [1] and mesons started to populate the Review of Particle Physics (RPP) [2], meson spectroscopy has played a central role in the study of the strong interaction helping on its understanding and the development of Quantum Chromodynamics (QCD). Mesons constitute the simplest bound states for quarks and gluons, and their accurate description is one of the principal aims of QCD. However, so far a quantitative and predictive theory of confined states has not been achieved, hence, in order to study the properties of mesons we have to rely on models which have to be consistent with the underlying QCD. Constituent quark models [3] are examples of this kind of modeling.
In the last decade an enormous experimental effort has been made in meson spectroscopy with several facilities conducting research programs [4] whose main goal has been to find exotic mesons [5] which do not fit within the quark-antiquark picture of quark models. This search has been fruitless so far but has put meson physics at the forefront of scientific research, becoming a thriving research area with experimental collaborations in several facilities -i.e. BES (China) [6], CLAS at JLab (USA) [7], COMPASS at CERN (Switzerland) [8], J-PARC (Japan) [9] and Hall D under construction at JLab (USA) [10].
Theoretical research has not been oblivious to this experimental interest and several quark models of mesons have made their appearance in the literature [11,12,13] trying to match the low-lying experimental spectrum and complementing the classic calculation by Godfrey and Isgur [14]. Among the theoretical developments, it is noteworthy the lattice QCD calculation Email addresses: laura@nuc5.fis.ucm.es (L. Muñoz), cesar@nuc2.fis.ucm.es (C. Fernández-Ramírez), armando.relano@gmail.com (A. Relaño) of the meson spectrum by the Hadron Spectrum Collaboration (HSC) at JLab [15,16], although with the drawback of being computed at a high pion mass of 396 MeV.
Mesons can be considered as aggregates of quarks and gluons. Therefore, the mass spectrum of low-lying mesons can be understood as the energy spectrum of a quantum system, like an atomic nucleus, and it consists on all the possible states which stem from an interacting quantum system. Since Wigner discovered that the statistical properties of complex nuclear spectra are well described by the Gaussian Orthogonal Ensemble (GOE) of Random Matrix Theory [17], statistical methods have become a powerful tool to study the energy spectra of quantum systems [18,19]. The most striking result in this field is that the statistical properties of the energy-level fluctuations determine if a system is chaotic, integrable or intermediate. Moreover, the energy-level fluctuations of integrable and chaotic systems are universal. The former display a non-correlated sequence of levels, which follows the Poisson distribution [20], whereas the latter are characterized by a correlation structure described by GOE [21]. This kind of analysis has been already applied to the hadron mass spectrum in [22] obtaining a chaotic-like behavior. In [23] the spectral-statistic techniques have been used to compare the experimental baryon spectrum with theoretical ones, focusing on the problem of missing resonances. The main conclusion of this work is that quark models give rise to integrable spectra, whereas the experimental one is chaotic. This result is not compatible with the existence of missing resonances, predicted by the models and not observed in the experiments. On the contrary, the lack of chaos in the models shows that some ingredients are missing in them, as they are unable to reproduce the experimental spectrum. As a corollary of this work, it is stated that the theoretical models must reproduce the spectral fluctuations of the experimental spectra, since they determine the dynamical regime and the complexity of the real interac-tions.
In this Letter we employ an improved version of the approach in [23] to extend the work to mesons. We infer the dynamical regime of light mesons from the statistical properties of their spectrum. We also study a number of theoretical models based upon constituent quarks [11,12,13,14], and the lattice QCD calculation by the HSC at JLab [15,16].

Statistical analysis and results
Prior to any statistical analysis of the spectral fluctuations one has to accomplish some preliminary tasks. First of all, it is necessary to take into account all the symmetries that properly characterize the system. Mixing different symmetries, i.e. energy levels with different values of the good quantum numbers, spoils the statistical properties deflecting them towards the Poisson statistics (see [17,19] for generic reviews and [24] for a recent work where the effects of both mixing symmetries and missing levels in the same sequence are surveyed). Hence, it is necessary to separate the whole spectrum into sequences of energy levels involving the same symmetries (good quantum numbers). The usual symmetries associated to mesons are spin (J), isospin (I), parity (P ), C−parity (C) and strangeness. Strangeness can be dropped due to the assumption of flavor SU (3) invariance and because strange mesons correspond to I = 1/2 while the rest of the mesons considered in this Letter are I = 0, 1 (we do not include in our analysis mesons with c, b or t quark content). Therefore, the meson spectrum is split into sequences with fixed values of J, I, P and C. In this work we take d = 126 energy levels from the RPP experimental spectrum up to 2.5 GeV of energy, which are distributed in n = 23 sequences of lengths l i ≥ 3, so that The energy spectrum of a quantum system is fully characterized by its level density ρ(E). It can be split into a smooth part ρ(E), giving the secular behavior with the energy, and a fluctuating part ρ(E), which is responsible for the statistical properties of the spectrum. The standard procedure to remove the smooth component of the level density is the unfolding. Once this procedure has been performed the fluctuations from different systems or from different parts of a single spectrum can be compared. Since the experimental meson spectrum has been divided in very short sequences of levels, we have to use the so-called local unfolding procedure [25] which considers that the variation of ρ(E) along a short sequence of energy levels must be negligible, a reasonable assumption in the present case. The procedure is as follows. Let {E i , i = 1, 2, . . . , l x } X be an energy-level sequence characterized by the set X of good quantum numbers. Then, the distances between consecutive levels, S i = E i+1 −E i , are rescaled using their average value S = (l x − 1) S i to obtain the quantities s i = S i / S , called generically nearest neighbor spacings (NNS). For the rescaled spectrum the mean level density ρ(E) = 1, and therefore s = (ρ(E)) −1 = 1. In this Letter, the statistical properties of the NNS are studied by means of the nearest neighbor spacing distribution (NNSD) [26], denoted P (s), which gives the probability that the spacing between two consecutive unfolded levels lies between s and s + ds. The NNSD follows the Poisson distribution P p (s) = exp(−s) for generic integrable systems [20] while chaotic systems with time reversal and rotational invariance are well described by the GOE of random matrices, whose NNSD follows the Wigner surmise P W (s) = πs 2 exp − πs 2 4 [21]. Figure 1 compares the P (s) distribution for the experimental spectrum with the Poisson distribution P P (s) and the Wigner surmise P W (s). The experimental distribution is intermediate between the Poisson and Wigner limits, albeit it seems closer to the latter.
Before we proceed with a quantitative analysis, it is convenient to study the possible spurious effects of the unfolding because the combination of the local unfolding with very short sequences of levels may distort the actual P (s) distribution [27]. Since s = 1 for every spacing sequence, no spacing can be greater than l − 1, where l is the sequence length, and therefore the P (s) distribution must exhibit a sharp cutoff at s = l − 1. When l is large enough this cutoff is irrelevant due to the exponential and Gaussian decays of the Poisson and Wigner distributions. Obviously, this is not the case for smaller values of l; in order to take into account this problem we will, as in [23], "distort" Wigner and Poisson predictions including the effects of the unfolding procedure in the same way as the experimental distribution. To generate these distributions we divide two generic GOE and Poisson spectra of dimension d = 126 in 23 level sequences with the same lengths of the experimental spectrum, and calculate their spacings. Then, we take a step further with respect to [23] and instead of using just one spectrum in each case, we obtain the smooth behavior of the distorted distributions by averaging over 1000 realizations of the spectra. The corresponding distorted Poisson and Wigner NNSDs are denoted by P DP (s) and P DW (s), and it is important to stress here that they will play from now on the role of theoretical distributions with which the data have to be compared.  displays the shape of these curves and a comparison with figure  1 shows that the unfolding effect is quite noticeable for P P (s). Figure 2 confirms that, after taking care of the unfolding distorsions, the statistical properties of the experimental spectrum are intermediate between the Poisson and Wigner predictions. In order to assess how close the fluctuations are to one limit or the other, some distributions can be used to describe intermediate NNSD. One of the most frequently used is the Berry-Robnik distribution P BR (f, s) [28], where f represents the fraction of phase-space volume dominated by chaotic orbits in the semiclassical limit. Of course, for a proper comparison, we need the corresponding distorted Berry-Robnik distribution P DBR (f, s), which we compute following the same procedure employed to obtain P DP (s) and P DW (s). Fitting the experimental P (s) to P DBR (f, s) we obtain for the fractional density f = 0.78 ± 0.13. Figure 2 also displays P DBR (0.78, s).
Additional information can be gained by performing the K-S distribution test [29]. The reference distribution is either P DP (s), P DW (s) or P DBR (0.78, s), and the null hypothesis is that the experimental P (s) distribution coincides with the one selected as reference, against the hypothesis that both distributions are different. The results for the p-value obtained in each case are p DP = 0.13, p DW = 0.47 and p DBR = 0.68, respectively. This result is fully compatible with our previous analysis of the NNSD. The statistical properties of the experimental meson spectrum are intermediate between the Poisson and Wigner limits. However, they are closer to the latter since a Berry-Robnik distribution with f = 0.78 fits well the experimental NNSD and p DBR = 0.68. It is also worth to note that p DP = 0.13 is close to the usual limit for the null hypothesis to be rejected (p 0.10).
As a complementary test, we calculate the moments of the distributions. Gathering together all the spacings s i , the k-th moment, M (k) is calculated as M (k) = (d − n) −1 d−n i=1 s k i , where d stands for the spectrum dimension and n is the number of spacing sequences. Figure 3 displays the moments M (k) , 1 ≤ Considering the experimental meson energies as Gaussian random variables with mean equal to the RPP estimation and variance equal to the error bar, we have generated 1000 "realizations" of the experimental spectrum and performed for each of them the K-S test. If the energy levels are allowed to fluctuate independently (in our case these fluctuations are induced by the error bars), the correlations are usually weakened, displacing the statistics towards Poisson. Therefore, it is appropriate to use as reference distributions P DP (s) and P DBR (0.78, s). Figure  4 shows that the histograms of the resulting p-values are separated with almost no overlap. The distribution of p DBR -values is concentrated in the upper half with p DBR = 0.77 ± 0.07, and the histogram of the p DP -values lies in the lower half with centroid p DP = 0.28 ± 0.07. It is important to notice that for almost every "realization" of the experimental spectrum p DBR > p DP , sustaining the good agreement of the experiment with the Berry-Robnik distribution for f = 0.78. For the sake of completeness we have also used as reference distribution the Wigner surmise, obtaining p DW = 0.46 ± 0.09.
To summarize, our analysis is fairly robust against experimental errors and allows us to conclude that the statistical properties of the experimental spectrum are intermediate between the Wigner and Poisson limits, closer to the former and safely incompatible with the latter. That is, mesons are much closer to chaotic systems than to integrable ones. Moreover, a Berry-  Robnik distribution with 78% of chaos provides the best description of the experimental NNSD.

Statistical significance of the analysis
The statistical analysis performed in the previous section accounts for the shortness of the sequences of levels through comparing to the theoretical ad hoc distorted NNS distributions. Also the size of the whole spectrum is not too large and, though spectra of this kind can be found in the literature (i.e. [27] and references therein), a careful analysis of the statistical significance is in order.
In the previous section we have found that the experimental spectrum is closer to a chaotic system than to an integrable one, actually, obtaining the best description employing a Berry-Robnik distribution with a 78% of chaos. To test the validity of the analysis one should start from a spectrum with the same size and structure (organized with the same number and length of the sequences) of the experimental one, explicitly built with 78% of chaos, and perform the analysis in order to see how safely can we arrive to the same conclusion. Then we build a random ensemble (RE) of 1000 spectra with a 78% of chaos, construct the NNSD in each case and run K-S tests against the NNS distributions P DW (s), P DP (s) and P DBR (0.78, s). We obtain the mean p-values and error bars (standard deviations) shown in table 1 together with the ones obtained for the actual experimental spectrum from the previous section. In Fig.  5 we show the NNSD for one of the spectra of the RE, where it can be seen that is similar to the experimental one, that is, the histogram is not very smooth due to the small size of the sample and the p-values are: p DW = 0.49, p DP = 0.20 and p DBR = 0.74.
From the mean p-values one can conclude that a typical spectrum of the ensemble is closer to chaotic systems than to integrable ones and is best fitted with a mixed distribution with 78% of chaos. From the error bars we can see that the range of p-values for Wigner is the broadest but the values inside the interval are such that the Wigner null hypothesis cannot be rejected. The Poisson hypothesis, although the range of p-values is lower (the lower edge is close to rejection), cannot be dismissed either. Thus, most probably the actual distribution lies in between both, and hence the case for trying a Berry-Robnik distribution, as we did for the experimental one. When the K-S test is run for a Berry-Robnik with f = 0.78, the highest pvalues of all three are obtained within a narrow range. Then, we can deduce that in most of the cases we would arrive to the same conclusion as for the experimental spectrum. But in order to better quantify the statistical significance one should look at each individual realization.
It is important to compare the three p-values in each individual case to see what would be the final decision in each one. If we do this we find that in 80% of the cases p DW > p DP and in 90% p DBR is the highest of the three. The Poisson hypothesis can be rejected (p < 0.1) in 10% of the cases, while the Wigner in only a 0.5% and the Berry-Robnik in none of them. Hence, at a confidence level of approximately 90% we are recovering the correct result: the spectrum lies in between the Wigner and Poisson extremes and the Berry-Robnik distribution with f = 0.78 is the best description of the spectrum. Additionally, we have run K-S tests for other values of f be- low and above f = 0.78 checking that in all the cases lower p-values are obtained. In summary, the methodology employed to analyze the experimental spectrum is reliable and robust.

Theoretical spectra
Next we analyze six theoretical calculations of the light meson spectrum and compare them to the results from the previous section. These are: (i) The classic model by Godfrey and Isgur (set GI) [14], which is a relativized quark model where the interaction is built employing a one gluon exchange potential and confinement is achieved through a spin-independent linear potential; (ii) and (iii) are the fully relativistic quark models by Koll et al. (sets K1 and K2 which correspond, respectively, to models A and B in [11]) based on the Bethe-Salpeter equation in its instantaneous approximation, a flavor-dependent two-body interaction and a spin-dependent confinement force, being the latter the difference between the two models; (iv) the relativistic quark model by Ebert et al. (set E) [13] based on a quasipotential (this calculation has the disadvantage that isoscalar and isovector mesons composed by u and d quarks are degenerate); (v) the effective quark model by Vijande et al.
(set V) [12], based upon the effective exchange of π, σ, η and K mesons between constituent quarks; and (vi) the lattice QCD calculation by the HSC (set LQCD) [16]. Lattice QCD calculation does not include strange mesons as the previous models, but it includes exotics such as the isoscalar J P C = 2 +− states. Table 2 displays relevant information on the six theoretical spectra, like their dimension d, the number n of pure sequences included in the analysis and the total number of spacings, which is equal to d − n. It also provides the p-values obtained by applying the K-S test to their NNSDs, taking as null hypotheses that the NNSD coincides either with P DW (s) or with P DP (s). The first relevant outcome is that, according to the K-S test, the NNSDs of sets K1, K2, E and LQCD are incompatible with the Wigner correlations and closer to the Poisson statistics. Thus, the dynamics predicted by these models is essentially regular, while the statistical properties of the experimental light meson spectrum show that the dynamical regime should be chaotic. This fact resembles the results obtained for baryons: while the fluctuations of the experimental baryon spectrum are well reproduced by Wigner predictions, the theoretical calculations give rise to spectra with Poisson statistics [23]. Figure 6 shows the NNSD of the six theoretical spectra. Sets K1 and E provide flat NNSDs with a cut at s = 2. The cut is expected as was explained in section 2.1. When the Poisson distribution is distorted it flattens due to the small amount of levels, so actually the NNSDs that we find for sets K1 and E are the ones we expect from a Poisson distribution, confirming that these sets have less correlations than the experimental data as the K-S test suggested. The comparison between models K1 and K2 by Koll et al. is particularly interesting because they only differ on the confinement interaction and show how important that interaction can be for the spectral statistics, hinting that it should be revised to obtain a better agreement with the experiment.
The result for set LQCD is particularly interesting because lattice QCD is currently the only tool available to compute lowenergy observables employing QCD directly. We find that the current state-of-the-art calculation in [16] does not describe properly the statistical properties of the meson spectrum. Lattice QCD NNSD is relatively close to the P DP (s) as it is shown in figure 6. This is evident at zero spacing where P LQCD (s = 0) ≈ 0.6. The value at zero spacing is critical to distinguish the Wigner and Poisson cases, being characteristic of a correlated spectrum and has an important impact in the K-S test. Our results remain unaltered if the statistical errors of the lattice QCD calculation are taken into account. Thus, the LQCD calculation should be considered a step forward in lattice calculations but still far away from being a description of the data or their structure. It is not something unexpected given that the LQCD set has been obtained at a pion mass of 396 MeV, far away from its physical mass, and it is reasonable to expect a drastic change in the statistical properties when calculations get the pion mass closer to its actual value. However, the fact that the lattice QCD calculation has a lot less correlations than the experimental data, being practically uncorrelated, demands, besides the need of bringing the calculation closer to the physical pion mass, a careful examination of the approximations employed.
The results of the K-S test for sets GI and V are inconclusive because they suggest that the models are compatible with both chaotic and integrable dynamics. It is thus mandatory to take a close look to the NNSDs (figure 6) before obtaining any conclusion. Set GI NNSD has a strange shape with peaks and dips, completely different from any of the usual distributions (Poisson, Wigner or Berry-Robnik), and therefore not close at all to experiment (see figure 2). It only has some similarity with some very particular integrable systems whose P (s) is equal to a sum of δ functions, constituting an exception to the rule of Poisson distribution [30]. On the contrary, set V displays a smooth NNSD, which can be very well fitted to a distorted Berry-Robnik distribution with f = 0.63±0.19 (also displayed in figure 6). Then we can conclude that the model by Vijande et al. gives a better account of the dynamical regime of the light meson spectrum.  Figure 6: (color online.) NNSDs for the theoretical spectra: i) Top left, set GI by Godfrey and Isgur [14]; ii) middle left, set K1 by Koll et al. [11]; iii) bottom left, set K2 by Koll et al. [11]; iv) top right, set E by Ebert et al. [13]; v) middle right, set V by Vijande et al. [12] and ad hoc distorted Berry-Robnik with f = 0.63; vi) bottom right, set LQCD by Dudek et al. [16]. Distorted Wigner (D-Wigner) is represented with diamonds, distorted Poisson (D-Poisson) with crosses and distorted Berry-Robnik (D-Berry-Robnik) with dots.

Conclusions
In this letter we have studied the spectral fluctuations of the light meson energy spectrum. These properties give information about the dynamical regime of the physical system, i.e., whether the system is chaotic, integrable or intermediate, and thus provide insight on the underlying interactions. Our analysis seems to be fairly robust against to the experimental errors, and unveils that the dynamical regime of the light mesons is near to the fully chaotic limit -our results show that the spectral fluctuations of light mesons resemble those of chaotic systems and indeed their NNSD can be well fitted to a Berry-Robnik distribution with a 78% of chaos.
On the other hand, the analysis of the spectra of the LQCD calculation and all the quark models but one shows that they are incompatible with chaos. This result is quite similar to that obtained for the baryons in Ref. [23]. Only the model by Vijande et al. predicts an intermediate spectrum with a 63% of chaos. The other quark models and the LQCD calculation [11,13,14,16] predict a regular or nearly regular dynamics in clear contradiction with the experiment.
We would also like to point out that the experimental spectrum can be totally chaotic, but the possible existence of missing levels, responsible for losing correlations, would deflect the actual NNSD towards a distribution between Wigner and Poisson. Ongoing experimental research will help to establish how large is the effect of the missing states in the spectral-statistics properties of the meson spectrum.
Further work is needed to study the origin of the failure of theoretical models. In the case of the LQCD spectrum it is probably due to the fact that the calculation is made at an unrealistic pion mass. For the quark models, the ratio between the quark-antiquark interaction and the confinement potential should be surveyed as well as the existence of states that go beyond quark-antiquark content.