Bottomonium spectrum with a Dirac potential model in the momentum space

We study the bottomonium spectrum using a relativistic potential model in the momentum space. This model is based on a complete one gluon exchange interaction with a momentum dependent screening factor to account for the effects due to virtual pair creation that appear close to the decay thresholds. The overall model does not make use of nonrelativistic approximations. We fit well established bottomonium states below the open charm threshold and predict the rest of the spectrum up to $\approx 11200$ MeV and $J^{PC}=3^{--}$. Uncertainties are treated rigorously and propagated in full to the parameters of the model using a Monte Carlo to identify if which deviations from experimental data can be absorbed into the statistical uncertainties of the models and which can be related to physics beyond the $b\bar{b}$ picture, guiding future research. We get a good description of the spectrum, in particular the Belle measurement of the $\eta_b(2S)$ state and the $\Upsilon(10860)$ and $\chi_b(3P)$ resonances.

In this paper we develop a relativistic quark model for bottomonia based on a complete one gluon exchange. The approach is completely relativistic and does not rely on nonrelativistic approximations. In this way the standard spin-orbit, spin-spin, and tensor interactions are automatically included. We also incorporate a relativistic scalar interaction and a momentum dependent screening factor to account for the effects due to virtual pair creation that appear close to the decay thresholds. All the calculations are performed in the momentum space. The same model was successfully applied to reproduce the charmonium spectrum in Ref. [55] which we refer the reader to for technicalities. We fit the model to all the known states of each J P C below the BB threshold except for the recently measured χ b1 (3P ) and χ b2 (3P ) which we prefer to predict in order to gain insight on their nature and the η b (2S) which we exclude of our fit owing to the disagreement between CLEO [56] and Belle [57] measurements. We perform a rigorous error estimation * djmolinab@unal.edu.co † cesar.fernandez@nucleares.unam.mx ‡ santopinto@ge.infn.it that allows us to assess if the inclusion of a new effect in the phenomenological model is necessary or not, and we compute the parameter correlations which provide insight on how independent are the different pieces of the model among them. A full error analysis is mandatory to identify which deviations from experimental data can be absorbed into the statistical uncertainties of the models and which can be related to physics beyond the bb picture, guiding future research.
The paper is organized as follows: In Sec. II we provide the relativistic quark model and the employed solution method; In Sec. III we describe the fitting procedure as well as the statistical method used to compute the uncertainties; In Sec. IV we report the computed bottomonium spectrum up to J P C = 3 −− and ≈ 11200 MeV as well as the comparison to the available experimental information. We obtain a very good description of both fitted and nonfitted bottomonia and also predict many unobserved states; Sec. V summarizes the conclusions.

A. Hamiltonian model
We apply to bottomonia the same model developed in [55] for charmonia. The total interaction Hamiltonian model in the momentum space is given by the sum of vector (H (v) ) and scalar (H (s) ) interactions where p a and p b represent the three-momenta of both quark and antiquark in the center of mass of the bottomonium system. The vector interaction is based on one gluon exchange, which in the Coulomb gauge reads where represents the standard four-current of the quarks, σ i stands for the Pauli matrices, and γ µ i are the gamma matrices, where i = 1, 2 is the particle label. We also introduce the quark energy difference and the squared (positive) four momentum transfer where q = p b − p a represents the three momentum transfer. The scalar interaction is defined as where I i is a scalar vertex. Finally, the vector and the scalar effective potentials have the following form: Equation (7a) represents a regularized Cornell potential, where α st is the coupling constant and β v corresponds to the vector confinement strength. Additionally, Eq. (7b) contains a phenomenological constant term A plus a β s term which corresponds to the scalar confinement strength. The constant parameter b has been introduced to avoid the divergence when | q| → 0. As in [55] for charmonia, we use two different prescriptions for the scalar interaction: potential I → model using Eqs. (7) with β s = 0, potential II → model using Eqs. (7) with β s = 0. (8) In this way we can check if the two forms of the scalar interaction (with or without the confinement term) have the same effect on the spectrum, as in the case of the charmonium system. Besides, in order to take into account the effects of the virtual [58,59] pair creation that appear close to the decays thresholds, we include a screening momentum dependent factor. Hence, the total Hamiltonian takes the final form where the factor F s (p) is defined as In this way, the model, with potential I and potential II, depends on seven and eight parameters, respectively.

B. Relativistic equation and solution method
The relativistic equation we use is obtained performing a three dimensional reduction of the Bethe-Salpeter equation and keeping only the contributions of the positive energy Dirac spinors [55]. In the center of mass of the bb system, the relativistic integral equation takes the form where we have introduced the energy and M 0 represents the phenomenological zero point energy of the spectrum, M is the resonance mass (i.e. the eigenvalue of the integral equation) and Ψ( p) is the resonance wave function. The wave function Ψ n,{ν} ( p) ({ν} = L, S, J) can be written as where R n,L (p;p) corresponds to the radial function in the momentum space with n the principal quantum number, p the variational parameter (with dimensions of momentum), Y L,M L (p) are the spherical harmonics, and χ S,M S is the spin function. To solve Eq. (11) we use the variational method. As trial functions we use a combination of a finite subset of three dimensional harmonic oscillators. Hence, we can write the Hamiltonian matrix as The eigenvalues and the eigenstates are found through the variational method, diagonalizing and minimizing the M {ν},n b ,na matrix in Eq. (14) [55,60]. The angular part is solved analytically and the radial part numerically. The details can be found in the Appendix of Ref. [55].

III. PARAMETER DETERMINATION
To determine the values of the parameters, the uncertainties, and the theoretical bottomonium spectrum we fit the experimental masses given in Table I, i.e. all the known states of each J P C below the BB threshold except for the recently measured χ b1 (3P ) and χ b2 (3P ) which we try to predict in order to gain insight on their nature and the η b (2S) which we prefer to exclude of our fit owing to the disagreement between CLEO and Belle masses [49]. From CLEO data a mass of 9974.6 ± 2.3 ± 2.1 MeV [56] was obtained while Belle measures 9999.0 ± 3.5 +2.8 −1.9 MeV [57]. BABAR reports a range value between 9974 and 10015 Mev [61]. The PDG favors the Belle measurement [12], therefore we show this experimental value in Table I and Figs. 1 and 2. To perform the fits and the error analysis we use the bootstrap technique [62][63][64][65] and proceed as follows: 1. We randomly choose values for the masses of the resonances by sampling a Gaussian distribution according to their uncertainties (systematic and statistical summed in quadrature), obtaining a resampled bottomonium spectrum; 2. We use the least-squares method to minimize the squared distance where M i are the resampled experimental bottomonia, i.e. the states 0 . The E th i represents the theoretical states calculated by solving the eigenvalue Eq. (14) with potentials I and II. The fit is performed using MINUIT [66].
This procedure is repeated 1000 times in order to obtain enough statistics to compute the expected values of the parameters as well as their uncertainties at a 1σ (68%) confidence level (CL). The expected value of the parameters (Table II) are computed as the mean value of the 1000 samples. The uncertainties are obtained as the the differences between the mean value and the highest and lowest masses of the best 68% of the fits. Hence, our uncertainties can be asymmetric. Once the parameters have been determined, we can compute the bottomonium spectrum and the associated uncertainties (Table I). We find an excellent agreement between theory and fitted states within uncertainties. We note that the values of the common parameters of the two potentials are very similar. These results show that, unlike for charmonia, the scalar confinement term of the interaction does not seem to be relevant in the bottomonia description. To gain further insight on this issue we compute the correlation matrices, Tables III and IV, for potentials I and II,  respectively. For potential I (Table III) we find a strong correlation between the parameters of the vector interaction (α s and β v ) and the scalar interaction parameter A, which indicates that vector and scalar interactions are physically correlated in this model. The screening parameter p s is weakly correlated with the vector interaction parameters but strongly correlated with the scalar interaction ones. For potential II, we have the additional parameterβ s . In this potential, the parameters are less correlated as shown in Table IV with one exception, the additional scalar interaction parameter β s is noticeble correlated with the vector interaction parameter β v . Consequently, we find a significant correlation between the confinement terms of the vector and the scalar interactions. The parameter p s of the screening factors is weakly correlated with the other parameters of the interactions except with the phenomenological parameter A in the scalar interaction. This sizeable correlation highlights how the screening factor impacts more on the scalar interaction.
Using the values obtained in the fitting procedure we plot the screening function F s (p) in Fig. 3 for the two potentials. As mentioned above we introduce the screening momenta p j 1/2 (j = I, II labels potentials I and II) which are given by F s (p j 1/2 ) = 1/2 (we recall that F s (0) = 1). Through the fitting values, we find p I 1/2 = 3.38 GeV and p II 1/2 = 3.34 GeV. These values correspond to the screening kinetic energȳ which amount toĒ I = 11.281 GeV for potential I and E II = 11.260 GeV for potential II. This result show that the screening effect is active above the open bottom threshold as in charmonia. Nevertheless, due to the high values ofĒ I,II , we find that the screening effect is less relevant for the low-lying part of the bottomonium spectrum than for charmonia [55].

IV. BOTTOMONIUM SPECTRUM
Using the relativistic model interaction, with either potential I or II, we obtain the bottomonium spectrum. Through the bootstrap method, the errors in the fitted states are carried in full to the computed uncertainties in the parameters and to the spectrum. We provide the computed spectrum in Tables I (fitted states) and V (predicted states). The computed and the experimental spectra are compared in Figs. 1 (potential I) and 2 (potential II). In general, the spectrum is reproduced by the model within the experimental uncertainties. We note that the parameters obtained with both potentials are very similar, leading to closely akin spectra. This result shows that the confining part of the scalar potential does not impact the bottomonium spectrum. However, the presence of the scalar interaction is necessary for an optimum fit, i.e. the parameter A contribution in Eq. (7b). In what follows we look into the states that were not included in the fit as well as the predicted higher-lying spectrum. These resonances belong to the family with quantum numbers 1 −− . They were discovered by means of e + e − collisions in the mid-eighties [68,69] and were more recently measured by the Belle collaboration [70]. The Υ(4S) is regarded as a 4 3 S 1 state; its experimental mass is M Υ(4S) = 10579.4 ± 1.2 MeV and is not well reproduced by either potential I or II. This resonance is generally considered as a bb state, but its mass is overestimated by models that make use of different approaches: e.g., the nonrelativistic model in Ref. [41] pro-vides M Υ (4S) 10630 MeV, the semirelativistic model of Ref. [49] finds M Υ(4S) = 10607 MeV, and the nonrelativistic coupled channels model in Ref. [43] reports M Υ(4S) = 10603 MeV. Our computations provides approximately 10642 ± 40 MeV, with both potentials. This result is compatible with the other models, but far away from the experimental value, even when the uncertainties are taken into account. Consequently, our result combined with non-relativistic calculations suggest that there must be beyond the qq picture effects that need to be included to properly describe the state.
The Υ(10860) resonance is generally interpreted as a Υ(5 3 S 1 ), e.g. in [41,43,47,49,50]. However, the theoretical calculations for the pion emission decay widths, to Υ(1S), Υ(2S) and Υ(3S) are two orders of magnitude [46] greater than the measurement [71] leading to different possible interpretations, such as that Υ(10860) is a mixing of a standard Υ(5S) with a P hybrid state [72], Finally, in Ref. [42] this state is interpreted as a Υ(6S), and, hence, the Υ(5S) becomes a missing resonance of the experimental spectrum In our model, this mass state can be reproduced as a Υ(5S) (5 3 S 1 ) (see Table V Table II. We highlight the values of the screening momenta p I 1/2 (potential I) and p II 1/2 (potential II) introduced in [55].
with both potentials. We do not find support the Υ(6S) interpretation. Actually, our predicted mean value mass, with Υ(5S), is only 5 (1) MeV away from the experimental value with potential I (II). Consequently, we identify this state as a bb with Υ(5S) quantum numbers. However, any final conclusion requires the explanation of the before mentioned pion emission decay widths which we leave for a future work.
Finally, the Υ(11020) state is mostly described as a bb meson in a 6 3 S 1 state except for [42] which interprets it as a 7 3 S 1 state. We do not find a satisfactory description of the mass of this state, neither as 6 3 S 1 nor 7 3 S 1 with either potential. In fact, our results are similar to that of a nonrelativistic model in [41]. Hence, we favor the existence of additional physics to explain the mass of this state, such as coupled channel effects as shown in [43] where a mass of 11023 MeV is obtained, very close to the experimental value.
The χ b (3P ) states have been the focus of several experimental collaborations during the last years. An estimation of the χ b (3P ) system barycenter (i.e. spinweighted mass average of the χ b0 (3P ), χ b1 (3P ), and χ b2 (3P ) states) was reported by ATLAS [1] and D0 [2] collaborations, yielding 10530 ± 5(stat) ± 9(syst) MeV and 10551 ± 14(stat) ± 17(syst) MeV, respectively. More recently, two out of the three state masses were measured; χ b1 (3P ) by the LHCb collaboration obtaining 10515.7 +2.2 −3.9 (stat) +1.5 −2.1 (syst) MeV, and χ b1 (3P ) and χ b2 (3P ) by the CMS collaboration [7] yielding 10513.42± 0.41(stat) ± 0.18(syst) MeV and 10524.02 ± 0.57(stat) ±    [44]. All of the results overestimate the mass of χ b1 (3P ). In our calculation (which porpously does not fit this state) we obtain 10540±30 MeV with both potentials whose central value also overestimates the mass of the state. When the uncertainties are taken into account, the experimental mass falls within our error bars and no indication of the need for additional physics is called for. This shows how important it is to perform a rigorous error estimation when performing a level-by-level comparison between theory and experiment, as differences that can be accounted by the error analysis can be mistaken by physics beyond the bb picture. Regarding χ b2 (3P ), 10532.4 Mev is obtained in Ref. [43] using the coupled channels formalism and 10578 Mev under the unquenched quark model [44]. We obtain 10554 +25 −28 and 10557 +22 −42 with potentials I and II, respectively. The CMS value falls well within our uncertainites for potential I and slightly out of them for potential II, although certainly within 2σ uncer-   10992.9 +10.0 −3.1 [12] tainties. Hence, the individually measured χ b (3P ) states are well reproduced by our model. Finally, we obtain the barycenter mass 10545 +24 −27 MeV for potential I and 10549 +23 −41 MeV for potential II, both compatible with the previously quoted ATLAS and D0 estimations. Recalling that not all the individual states of the χ bJ (3P ) system have been measured, we provide in Tables VI (poten-tial I) and VII (potential II) the n = 1, 2, 3 barycenter masses, given by [73,74] along with the available experimental measurements and estimates from PDG values Given that both potentials Theoretical results obtained, using Potential I, for the states χ bJ (nP ) compared with the available experimental data; n = 1, 2, 3 is the principal quantum number;Mn stands for the barycenter of the system for each n. The experimental states for n = 1, 2 are taken from Ref. [7]. The statistical and systematic errors of the experimental states have been summed in quadrature in order to obtain the errors of the experimental barycenter masses. The theoretical uncertainties of the barycenters were propagated from the parameters through the bootstrap technique. produce similar spectra, the χ b (nP ) barycenters are very similar. In summary, we find a good agreement, within errors, between the models and the experimental barycenters.
Finally, we would like to mention that it has been theorized that some of the states in the χ b (3P ) system could be the bottomonia counterparts of the X(3872) charmonium [45,75], i.e. states closely related to the opening if the BB, BB * , and B s B s thresholds. Our results do not support such hypothesis, as the model reproduces the χ b (3P ) system within (large) uncertainties, contrary to the X(3872) case which was overestimated using the same model [55], and whose description (both mass and width) calls for additional dynamics beyond the cc picture. Along the same ideas, according to Ref. [76], the χ b1 (4P ) state could significantly couple to the BB * and B * B * channels. The measurement of this particular state combined with the comparison to quark model calculations, like the one presented in this work, can provide insight on the impact in the masses of the dynamical effects due to the open bottom thresholds.

C. Missing resonances
Besides reproducing the experimentally established states, in Table V we provide predictions of states both above and below the open bottom thresholds (≈ 10.6 GeV). In total, we predict 38 states up to 11.3 GeV for 0, 1, 2 (with either ± combinations for P and C) and 3 −− quantum numbers. These predicted states are of interest for future analysis at LHCb [15,77,78] and Belle II [13,14,[79][80][81]. In particular, pinning down the Υ(6S) would provide further insight on bottomoniumlike states [81].
The missing η b (nS) sector (n 1 S 0 states) can be studied through their relation to their angular momentum partners Υ(nS) (n 3 S 1 ) -known from experiment-, by computing the ∆S n = n 3 S − n 1 S mass splitting. This difference should decrease as n increases in the potential model context [82]. The experimental data for ∆S 1 and ∆S 2 shown in Table VIII support this theoretical results. Thereby, we consider our mass estimations for both η b (nS) and Υ(nS) reasonable.
We also provide predictions for states of the n 1,3 D 1,2,3 family, which remain undetected except for the 1 3 D 2 resonance. The predicted missing states (with uncertainties) provide useful information to guide the forthcoming spectroscopy programs in Belle II [13,14] and LHC [15,16]. However, the production rate of these states should be low, hence, difficult to detect [77].

V. CONCLUSIONS
We have developed a relativistic quark model in momentum space to study the bottomonium spectrum. The model closely follows the one used in Ref. [55] to study charmonium. It combines vector and scalar interactions with a momentum dependent screening factor to account for the effects due to virtual pair creation that appear close to the decay thresholds. We fitted our model to all the known states of each J P C below the BB threshold except for the recently measured χ b1 (3P ) and χ b2 (3P ) which we prefer to predict in order to gain insight on their nature and the η b (2S) which we exclude of our fit owing to the disagreement between CLEO and Belle measurements. Our prediction for η b (2S) mass agrees with the Belle result.
We have performed a full statistical error analysis using the bootstrap technique, that provides a rigorous treatment of the statistical uncertainties. In this way we obtain the uncertainties of the parameters and their correlations and we can propagate both to the predicted spectrum. Previous error analysis within phenomenological models have been very limited and incomplete. The rigorous error estimations allow us to assess if the inclusion of a new effect in the phenomenological model is necessary or not, and the correlations provide insight on how independent are the different pieces of the model among them. A full error analysis is mandatory to identify which deviations from experimental data can be absorbed into TABLE VIII. Differences ∆Sn = n 3 S − n 1 S. We observe that these differences decrease when n is increased. All the differences reported in this the statistical uncertainties of the models and which can be related to physics beyond the bb picture, guiding future research. We find that the model reproduces very well the fitted states as well as the nonfitted ones within uncertainties.
To asses the importance of a confining term in the scalar interaction, i.e. β s = 0 in Eq. (7b), we fitted the data with and without such contribution. The results obtained with the two potentials are very similar for the fitted and the predicted states, both in the low and the high parts of the spectrum. Therefore, such confining contribution to the scalar interaction can be disregarded in a bottomonium relativistic model. Even so, the correlations found among the parameters belonging to the scalar interaction and the rest of the model parameters, show that the scalar interaction A in Eq. (7b) is strictly necessary to reproduce the spectrum. The screening factor F s (p) included in the interaction Hamiltonian begins to impact the predictions in significant way at ≈ 11200 MeV, i.e. further away from the open bottom decay thresholds. Hence, the screening effect is not particularly intense and has a slight impact on the bottomonium spectrum, contrary to what it was found for the charmonium one [55].
We have also studied the χ b (3P ) resonances. In particular we have calculated the mass of each state of this system and its barycenter. The experimental mass value of the χ b1 (3P ) falls into the theoretical uncertainty calculated with both potentials. Whereby, we conclude that the model is able to properly predict this state. Also, the model, with both potentials, reproduces the χ b1 (3P ) state. Our result indicates that the χ b1,2 (3P ) states are more likely to be bb mesons than the hypothetical X b states.
Our model overestimates the Υ(4S) mass and is consistent with results obtained by semirelativistic quark models, within errors. This is an indication of physics beyond the bb picture for this state. We identify the Υ(10860) as a 5 3 S 1 state and the model fails to reproduce the Υ(11020), although it is well reproduced by other potential models that take into account coupled channel effects [43]. Hence, the first can be considered (mostly) a bb state while the latter is up for discussion.
Finally, we report some states that, up to now, have not been observed experimentally but the confirmation of their existence is part of the experimental plans at LHC B factories and Belle II.