Higher bottomonium zoo

In this work, we study higher bottomonia up to the nL=8S\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$nL=8S$$\end{document}, 6P, 5D, 4F, 3G multiplets using the modified Godfrey–Isgur (GI) model, which takes account of color screening effects. The calculated mass spectra of bottomonium states are in reasonable agreement with the present experimental data. Based on spectroscopy, partial widths of all allowed radiative transitions, annihilation decays, hadronic transitions, and open-bottom strong decays of each state are also evaluated by applying our numerical wave functions. Comparing our results with the former results, we point out difference among various models and derive new conclusions obtained in this paper. Notably, we find a significant difference between our model and the GI model when we study D, F, and G and n≥4\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n\ge 4$$\end{document} states. Our theoretical results are valuable to search for more bottomonia in experiments, such as LHCb, and forthcoming Belle II.


Introduction
Since J/ψ and ϒ(1S) were observed in 1974 and 1977 [1][2][3][4], respectively, a heavy quarkonium has become an influential and attractive research field because its physical processes cover the whole energy range of QCD. This energy range provides us an excellent place to study the properties of perturbative and non-perturbative QCD [5]. An exhaustive and deeper study of a heavy quarkonium helps people better understand the QCD characteristics. Since the bottomonium family is an important member of heavy quarkonia, people have made great effort to investigate a bottomonium experimentally and theoretically in the past years. Thirty years ago, a couple of higher excited bottomonia have been expera e-mail: wangjzh2012@lzu.edu.cn b e-mail: sunzhif09@lzu.edu.cn c e-mail: xiangliu@lzu.edu.cn d e-mail: matsuki@tokyo-kasei.ac.jp imentally observed in succession, e.g., ϒ(nS) with the radial quantum numbers n from 2 to 6 and P-wave spin-triplet states χ bJ (1P) and χ bJ (2P) with J = 1, 2, 3 [3,4,[6][7][8][9][10][11]. With the upgrade and improvement of two B-Factories BaBar and Belle together with LHC, new breakthroughs in experiments have been achieved in recent years. First, the pseudoscalar partners η b (1S) and η b (2S) were identified by the BaBar Collaboration in 2008 and 2011, respectively [12,13]. Subsequently, the follow-up studies by the CLEO [14,15] and Belle Collaborations [16] confirmed the existence of these two bottomonia, where Belle did the most accurate measurements of the mass to date with values of 9402.4 ± 1.5 ± 1.8 MeV and 9999.0 ± 3.5 +2.8 −1.9 MeV, respectively [16]. The χ b1 (3P) state is a new state discovered by LHC in the chain decay of χ b1 (3P) → γ ϒ(1S)/γ ϒ(2S) → γ μ + μ − [17]. This state has been confirmed by the D∅ Collaboration [18]. BaBar discovered a possible signal of the P-wave spin-singlet h b (1P) state in 2011 [19], and subsequently, Belle firstly and successfully confirmed the observation of the h b (1P) in the process of ϒ(5S) → π + π − h b (1P), where the first radial excited state h b (2P) has been also observed with mass value 10259.76 ± 0.64 +1. 43 −1.03 MeV [20]. Thus, both the ground and first radial excited states of the P-wave bottomonium have been fully established experimentally. However, this is not enough to make people satisfied because many of experimental information is still incomplete, including the total width and branching ratios of some significant decay channels, which still require both of experimental and theoretical efforts. With regard to the search for a D-wave bottomonium in the experiments, there also has some progresses by the CLEO and BaBar Collaborations [21,22], where spintriplet ϒ(1 3 D 2 ) was observed with a 5.8 σ significance while confirmation of ϒ(1 3 D 1 ) and ϒ(1 3 D 3 ) states is dubious on account of lower significance [22].
Such plentiful experimental achievements in the bottomonium family not only increase our motivation to search for more higher excited bottomonia in future experiments but also provide an excellent opportunity to test various theories and phenomenological models. To better understand the properties of a bottomonium, both the progress of experimental measurements and the calculation of theoretical and phenomenological methods are necessary. The lattice QCD is usually considered to be the most promising solution to the non-perturbative difficulty in low energy regions. In principle, the spectrum of the heavy quarkonium is directly obtained from the lattice QCD, but actually, it is quite troublesome because an additional large heavy quark mass m Q needs a tremendous computational effort than that of the light quarkonium [23]. In spite of this, the lattice QCD still has an advantage in the calculations of bottomonia [24,25]. Before the lattice method becomes more comprehensive than the past, the mass spectra of hadrons are usually treated by the phenomenological model, namely, potential model, which takes into account the non-perturbative effects of QCD in the interaction potentials and can give satisfactory results consistent with the experiments.
In the past decades, a variety of potential models have been used to study the bottomonium system [26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42][43][44][45], among which the most well-known model is the Godfrey-Isgur relativistic quark model. This model has been widely used in the study from light mesons to heavy mesons [26,46]. Recently, Godfrey and Moats used the Godfrey-Isgur relativistic quark model for systematically investigating the properties of the bb system. This system contains a large number of higher excited states, and they have given the results of the production in e + e − and pp collisions so that experimentalists can look for and observe the most promising new bottomonia [26] in future. However, this is similar to the case of other meson families with abundant experimental information. That is, high-lying states are related to the higher mass values predicted by the GI model. In comparison with the measured results of bb states, one of the most salient particles of the model is the ϒ(11020) or ϒ(6S) whose theoretical mass value is about 100 MeV larger than the experimental data. In Ref. [47], the present authors proposed a modified GI model with color screening effects to investigate higher excited charm-strange mesons. Namely, they have taken into account the effect that the confinement potential is softened by the influence of light quarks induced from the vacuum in the long-range region [48], and have found that the addition of this effect can well describe the properties of higher charmstrange mesons. Hence, the recent study of bottomonia by the GI model [26] motivates us to explore what different conclusions can be drawn after considering screening effects in the GI model. With this motivation, we will carry out a most comprehensive study on the properties of the bottomonium family by the modified GI model with screening effects, primarily focusing on the properties of higher excited states.
In this work, we calculate the mass spectra and decay behaviors of bottomonia using the modified GI model, including computations of radiative transition, hadronic transition, annihilation process, and OZI-allowed two-body strong decay, where the meson wave functions also come from the modified GI model. This model will be thoroughly introduced in the next section. After the χ 2 fit with abundant experimental information, we can give a fairly good description for the mass spectra of bottomonia, where the masses of the higher excited states have been dramatically improved compared with the estimates by the GI model [26]. At the same time, we also predict the mass values of higher bottomonia, which are valuable for the experiments such as BaBar, Belle, and LHC to search for these missing particles. From the mass spectra of bb states, we discuss bottomonia by dividing them into those below and above open-bottom thresholds.
For the states below the threshold, since strong decay channels are not open, radiative transitions and annihilation decays are usually dominant. Dipion hadronic transitions are, of course, also major decay modes. Comparing calculated partial widths of these modes with those of the GI model, we find that the results between the two models are almost the same for the low-lying states but are different from them for the higher excited states such as 1F or 1G states. It is shown that the screening effect has little influence on the low-lying states but it is crucial for the higher states. For the states above the threshold, we calculate the strong decays using the 3 P 0 model, where we adopt a meson wave function numerically obtained rather than the SHO wave function with a corresponding effective β value. Because the states above the threshold have relatively large masses, an influence of the screening effect becomes more obvious for these states and our results that are largely different from the GI model should give a better prediction for the higher excited states.
In addition, it is emphasized that the discussion and prediction on the higher bottomonia are the main points of our work derived from our calculations. We hope that the present investigation can not only reveal the inherent properties of the observed bottomonia but also provide valuable clues to the experimental search for more missing bb states in the future.
This paper is organized as follows. In Sect. 2, we will give a brief introduction to the modified GI model with the screening effect and compare the results of the GI model and modified GI model. With our mass spectrum, the bottomonium spectroscopy will be analyzed in this section. Combining information of mass spectra and decay behaviors, we will study properties of the states bellow the BB threshold and those above the threshold in Sects. 3 and 4, respectively, where plentiful predictions will be given. In Sect. 5, we compare the results between the modified GI model and a coupled-channel quark model. We make a summary in Sect. 6. Finally, all the theoretical tools of various decay processes and physical quantities, like partial decay width, branching ratio, annihilation decay, radiative transition, hadronic transition, and total width, will be presented in Appendix A. A complete list of interaction potentials of the modified MG model will be given in Appendix B.

Modified GI model with a screening effect
In this work, we use the modified Godfrey-Isgur model with a screening potential [47] to calculate bottomonium mass spectrum and wave function, which will be employed in the calculation of decay processes. The Godfrey-Isgur relativized quark model [46] is one of the most successful model in predicting mass spectra of mesons. Even though the GI model has achieved great success, for the higher orbital and radial excitations, the predicted masses are larger on the whole than the observed masses of newly discovered particles in recent decades. Coupled-channel effects usually play an important role for the higher excitations and a screening effect could partly substitute the coupled-channel effect [49,50]. In 1985, Godfrey and Isgur proposed a relativized quark model motivated by QCD [46]. Compared with other models, relativistic corrections and universal one-gluon-exchange interaction as well as linear confinement potential are the most important features of this model. In the following, before presenting our modified GI model, we will introduce the GI model one by one. The Hamiltonian of a meson system is composed of kinetic energy and interaction between quark and anti-quark, and kinetic energy adopts a relativistic form as where m 1 and m 2 are constituent quark masses corresponding to quark and anti-quark, respectively. The interaction between quark and anti-quark includes short-range onegluon-exchange potential G(r ) and long-range confinement S(r ), spin-orbit interaction, color hyperfine interaction (contact and tensor term). G(r ) and S(r ) have the following forms as and S(r ) = br + c, where α s (r ) is a running coupling constant. Relativistic contributions are divided into two parts. Firstly, the model makes a smearing transformation of G(r ) and S(r ). To express it shortly, we use a general function symbol V (r ) instead of G(r ) and S(r ). A special operation follows asṼ (r ) = V (r )ρ(r − r )dr 3 (4) with ρ(r − r ) = σ 3 π 3 2 e −σ 2 (r − r ) 2 , where σ is a smearing parameter. Second, an important reflection of relativistic effects lies in the momentum dependence of interactions between quark and anti-quark. Therefore, a momentum-dependent factor is brought into the interactions. In this situations, a one-gluon-exchange potential G(r ) is modified as Tensor, contact, scalar spin-orbit, and vector spin-scalar terms should be changed as where E 1 and E 2 are the energies of quark and anti-quark and ε i corresponds to i-th type of interactions (ε c , ε t , ε sov , and ε sos ). We readily notice that in the non-relativistic limit the factors become unity, and particular values of ε i can be obtained from Table 1.
After a brief review of the GI model, we will introduce a modified GI model with a screening potential. Our previous work [47] presented the modified GI model with a color screening effect and revisited the properties of charm-strange mesons, subsequently, charm mesons have been also studied in the framework of the modified GI model [51], where higher excitations have been greatly improved. A screening effect can be achieved by modifying a confinement potential as where μ is a parameter which expresses the strength of a screening effect. As described in the preceding paragraphs, the modified confinement potential Eq. (8) also needs relativistic corrections. Details of techniques can be found in Ref. [47], and a complete list of the interaction potentials in the modified GI model are listed in Appendix B. It is worth mentioning that the modified GI model also employs simple harmonic oscillator (SHO) wave functions which can be considered as a complete basis set to expand exact meson wave functions. In the momentum space, an SHO wave function has the form   Table 1 with where Y L M L ( p ) is a spherical harmonic function, L L+1/2 n−1 p 2 β 2 is an associated Laguerre polynomial, and β is a parameter of an oscillator radial wave function. In the next subsection, we calculate bottomonium mass spectra via the modified GI model, which help us well understand the bb family and are also used for the subsequent decay processes. Fig. 1 shows a comparison of confinement potentials of modified GI model with screening effects and the GI model.

Mass spectrum
Due to the introduction of a new parameter μ in our modified GI model, we need to combine experimental information to scale all available model parameters. At the same time, the richness of bottomonia in the experiment also provides us with an excellent opportunity to determine these model parameters, where the χ 2 fitting method is adopted. The χ 2 fitting method is to find the minimum χ 2 value, thereby, to obtain a set of corresponding fitting parameters in which the theoretical predictions of phenomenological model and experimental results are most consistent on the whole. χ 2 is defined as where the V E x p (i), V T h (i) and V Er (i) are experimental value, theoretical value and error of i-th data, respectively. We select eighteen bottomonia to fit our model parameters as shown in Table 2 that have been established in the experiments. In this Table, all experimental masses are taken from PDG [52] and a uniform value of V Er (i) = 5.0 MeV is chosen as a fitting error for all the states, which is larger than their respective experimental uncertainty. The reason is that the experimental errors of these particles are relatively small and unevenly distributed. In other words, if the error of a particle is much smaller than the others, its mass will be too precise so that it is unfavorable to the overall fitting. The final fitting χ 2 value is given as 11.3 in conjunction with a uniform error and all corresponding fitting parameters of the modified GI model are listed in Table 1, where parameters of the GI model are also given for comparison.
In Table 2, the theoretical mass values of the GI model are also listed and the corresponding χ 2 value is calculated to be 31.4. Comparing χ 2 values, one can easily see that the fitted modified GI model well improves the whole description of the bottomonium spectrum. Although the GI model was successful in investigating the bottomonium spectrum, comparing with the observed masses, there are two main shortcomings in the theoretical estimation of the GI model. The first one is that the theoretical mass values of low-lying states are universally 10-20 MeV smaller than the expeimental values. The second is that the theoretical predictions for higher excited states become larger than the experiments. This is because the screening effect is not taken into account, especially for the ϒ(6S) state whose mass difference is close to 100 MeV. It is worth mentioning that the experimental mass value of ϒ(10860) or ϒ(5S) state is overestimated. That is, the recent BaBar experiment tends to give a 30 MeV higher than the initial experimental value [52]. This means that the theoretical prediction of the modified GI model for this state will be inevitably small and this can also explain why the fitted χ 2 value can not be much smaller. In fact, other studies of non-relativistic potential models that consider a screening effect [28,29,43] are also not very satisfactory for the description of this state, ϒ(5S). Except for ϒ(10860), our results for the bottomonium spectrum are pretty good that can be seen from Table 2. Our model not only solves the two above shortcomings of the GI model but also gives a fairly precise theoretical estimate, where the calculated mass values of many particles are almost the same as experimental results within the fitting errors. Based on our excellent description for the bottomonium spectrum with experimental information, we also predict the mass values for higher bottomonia from S-wave to G-wave that have not yet been observed in Table 3. In order to facilitate readers to understand the bottomonia more intuitively, the mass spectra are also shown in Fig. 2, which could give an overall grasp of the spectroscopy and is convenient for the comparison among the results of two models and experiments. As seen from Fig. 2, after introducing the screening effect, the impact of a screening potential on the ground states and low-lying states is not so obvious. In contrast, the screening effect begins to be important in the region of higher excited states, especially for those with larger radial quantum numbers n. The greater the n value the more prominent this effect is because of the more and more obvious mass differences between the GI model and modified GI model. This feature can also be reflected from the wave function of bottomonia and we take the S-wave ϒ family as an example, whose radial wave functions are shown in Fig. 3. Here the wave functions of the ϒ(mS) with m = 1, . . . , 4 are almost undistinguishable but the higher radial excited states appear to be visibly different in the positions of nodes and radial distributions. Of course, to thor-oughly understand the impact of the screening effect and the nature of the bottomonium family, we need to analyze their decay behaviors, which is also the main task of the later sections.

The η b states
The η b family with a quantum number 1 S 0 is the partner of spin-triplet ϒ states. The ground state η b (1S) and radial excited state η b (2S) have been established in the experiment, and their measured average masses are 9399.0 ± 2.3 MeV and 9999 ± 3.5 +2.8 −1.9 MeV [52], respectively, which are in pretty good agreement with our theoretical values in Table 2. There is a more interesting physical quantity, i.e., the hyperfine mass splitting between the spin-singlet and spin-triplet , which reflects the spin-dependent interaction and can be used to test various theoretical models. For the 1S state, our theoretical result of hyperfine mass splitting is 65 MeV within the upper limit of error, which is well consistent with the experimental result of 62.3 ± 3.2 MeV [52] and the lattice calculation of 60.3 ± 7.7 MeV [53]. The m(2S) from our modified GI model is estimated as 28 MeV which is also equally well consistent with the measured value of 24.3 +4.0 −4.5 MeV [16] and the lattice computation of 23.5 − 28.0 MeV [53]. It is worth mentioning that the predicted mass splittings from the GI model and our modified GI model are exactly the same, although there are some differences in the respective masses. Together with the successful description for the mass of η b system, we also explore their decay behaviors. We list partial widths and branching ratios of electromagnetic decays, annihilation decay, and hadronic transitions for η b states below the open-bottom thresholds from η b (1S) up to η b (4S) in Table 4. Here, the OZI-allowed two-body strong decay channels do not yet open so that the annihilation into the two gluons is dominant, whose corresponding branching ratio is almost 100%. From Table 4, we see that the decay widths from our modified GI model are not so different from those of the GI results [26]. This is because screening effects are weak for the wave function of low-lying states as discussed in Sect. 2. The decay predictions of η b (1S) and η b (2S) are also consistent with those of nonrelativistic constituent quark model [27].
For the two unobserved η b (3S) and η b (4S) state, we predict the masses of 10336 and 10597 MeV and total widths of 5.52 and 4.07 MeV, respectively. The corresponding hyperfine mass splitting m(3S) is 20 MeV while the result of nonrelativistic constituent quark model [27] is 19 MeV, and our estimate for m(4S) is 15 MeV. These predictions require further validation by the future experiments. In addition to the dominant two gluon annihilation decay, some other possible main decay channels of η b (3S) are η b (1S)π π and h b (2P)γ , both of which have almost the same branching ratios. The decay mode h b (1P)γ is also estimated to be important in our calculation, but the corresponding prediction of nonrelativistic constituent quark model [27] is two orders of magnitude smaller than our result. Similarly, the hadronic transition to η b (1S)π π and E1 radiative transition to h b (3P)γ are predicted to be important decay modes of η b (4S) with the branching ratios of (9.1 × 10 −3 ) and (3.7 × 10 −4 ), respectively.

The ϒ states
From Fig. 2, we can see clearly that ϒ(4S) state is slightly above the open-bottom threshold. Hence, the states below the thresholds are only ϒ(1S), ϒ(2S), and ϒ(3S), which were first discovered bottomonia together in the E288 Collaboration at Fermilab by studying produced muon pairs in a regime of invariant masses larger than 5 GeV [3,4]. At present, these three particles do not bring much controversy due to the high accuracy of the experimental measurements for their masses, which can be usually matched by the calculation of most of the potential models and the lattice QCD. The mass differences between our theoretical estimates and experimental center values for these three states are 3, 5 and 1 MeV, respectively, which also indicates the reliability of our modified GI quark model in the mass spectrum. In addition, the BaBar collaboration [54]  Partial widths and branching ratios of electromagnetic decays, annihilation decay, and hadronic transitions and the total widths for ϒ(1S), ϒ(2S), and ϒ(3S) are listed in Table 5. Compared with the S-wave spin-singlet η b family, the experimental information of the spin-triplet ϒ states is clearly much more abundant, which includes total widths and partial rates of most of the decay processes. Next, we   Table 5, the annihilation decay to three gluons is dominant since each branching ratio is 89%, 63.8% and 58.7% from our model, respectively, and the contributions to the total width from other three annihilation decay modes + − , γ gg, and γ γ γ are much smaller. Especially, the γ γ γ decay is difficult to search for in the experiment due to the very small branching ratio, the 10 −5 ∼ 10 −6 order of magnitude, and leptonic annihilation decay and γ gg decay modes have almost the same predicted partial widths. All the experimental widths and branching ratios basically agree with our theoretical calculations although the errors in ϒ(3S) state are large. Additionally, it needs to be emphasized that the experimental partial widths marked as b shown in Table 5 are obtained by combining measured total widths and branching ratios of the PDG [52].
There are some difficulties in the theoretical description of ϒ(3S) as a whole. Our theoretical total width of ϒ(3S) is 35.8 keV, which is larger than the PDG result of 20.32 ± 1.85 keV. The excess mainly comes from the annihilation mode of ggg and hadronic mode of ϒ(1S)π π . Considering the uncertainty of the phenomenological models, we make a comparison in radiative transitions by using experimental partial widths rather than directly measured branching ratios.  [52]. The theoretical results of the GI model [26] and the nonrelativistic constituent quark model [27]  a From the summation of partial widths of Ref. [27] b From the calculation of combining experimental total widths and branching ratios of the PDG [52] Our predictions for χ bJ (2P)γ are satisfactory and in spite of some deviations, the χ bJ (1P)γ and η b (1S)γ modes can be ensured in the order of magnitude.
In the same year, 2012, the Large Hadron Collider brought good news about the confirmation of the first observed 3P bottomonium state χ b1 (3P) in the chain decay of the radiative transition to ϒ(1S)γ and ϒ(2S)γ and to ϒ(1S, 2S) → μ + μ − . This experiment reported the measured mass of 10530 ± 5 ± 9 MeV [17], whose central value is consis-tent with our theoretical prediction of 10527 MeV in Table  2 although statistical and systematic uncertainties are relatively large. Several subsequent measurements from LHCb give a mass value about 15 MeV less than that in Ref. [17], which are incompatible with the predictions of most of the potential models. Therefore, further study on χ b1 (3P) in more experiments is necessary. Furthermore, the predicted masses of other 3P states are listed in Table 3. Table 2 shows that our theoretical masses for all of 1P and 2P states are in remarkable agreement with the experiments with the error of less than 5 MeV, which again proves the validity of the modified GI model with a screening effect. Partial widths and branching ratios of radiative transition, annihilation decay, and hadronic transition and the total widths for 1P to 3P bottomonium states are successively given in Tables 6,7,8,9. Because there are currently no measured total widths or partial widths available for the observed 1P and 2P states, we will analyze the decay properties of Pwave bottomonium states directly from the branching ratios. For 1P states, there are actually not many open decay channels and the annihilation processes of multi-gluon or hybrid qqg final state have dominant contributions to the total width. As for radiative transitions, the CLEO collaboration released the most accurate measurement on electromagnetic process of χ bJ (1P) → ϒ(1S)γ with J = 1, 2, 3, whose branching ratios are 1.76 ± 0.30 ± 0.78%, 33.1 ± 1.8 ± 1.7%, and 18.6 ± 1.1 ± 0.9% [55], respectively.
In the past few years, the Belle collaboration has studied the decay mode of h b (1P) → η b (1S)γ , whose branching ratio is 49.2 ± 5.7 +5.6 −3.3 % in 2012 [16] and 56 ± 8 ± 4% in 2015 [56]. Combining these results with the average ratios of PDG [52], we can see that the theoretical estimates shown in Table 6 are also generally supported by experimental data.
Different from the case of 1P states, the measurements of E1 radiative decays χ bJ (2P) → ϒ(1S, 2S)γ with J = 1, 2, 3 and h b (2P) → η b (1S, 2S)γ are much larger than theoretical results as shown in Table 7. It may be due to an overestimation of annihilation decay calculation or the experimental measurement error. After all, the uncertainties of some experimental data are quite large.
In addition, it is worth noting that our results of critical hadronic decay channel of χ bJ (2P) → χ bJ (1P)π π with non-flip J and h b (2P) → h b (1P)π π calculated by the QCD multi-pole expansion approach are roughly consistent with those of the GI model [26]. However, our results for the process of χ bJ (2P) → χ bJ (1P)π π with J = J are incompatible with the GI model, where our estimates are quite suppressed similar to the calculations of nonrelativistic constituent quark model [27]. The analogous situation also exists in the 3P and D-wave bottomonium states, which needs to be identified by the future experimental measurements.
Although the first 3P state has been experimentally observed, the experimental information of its decay behaviors is still lacking. Our predictions for various decay channels of 3 1 P 1 , 3 3 P 0 and 3 3 P 1 , 3 3 P 2 states are listed in Tables  8 and 9, respectively. By comparing our results with those of the GI model [26], we find that although most of decay modes have not changed much, the screening effects have demonstrated the power in the M1 radiation transitions. Taking the mode h b (3P) → χ bJ (1P)γ as an example, there is an order of magnitude difference between the two models. Such a large difference mainly comes from the change of bottomonium wave functions rather than the phase space. To illustrate it, we compare the square of the M1 radiation matrix element ψ f j 0 where ω is almost the same for both. For the final states χ bJ (1P)γ with J =0, 1, 2, their numerical results calculated by the GI model are 2.56×10 −4 , 9.61×10 −4 and 1.02×10 −3 [26], and those of our modified GI model are 1.38 × 10 −3 , 1.08 × 10 −4 and 1.13 × 10 −4 , respectively. This comparison once again proves the importance of the screening effects in describing the inner structure of meson states. In addition, the total widths of the 3P states predicted by us are higher on the whole.
In addition to the annihilation decay ggg, the decay is confirmed in the future, we suggest experiments to search for the missing η b (3S) state by studying this radiative process of h b (3P). Likewise, compared to other radiative transition processes, the decay process to the S-wave ϒ state is very significant for χ bJ (3P) states. It should be noted that the χ b1 (3P) state is just found in the chain decay of the radiative transition to ϒ(1S, 2S)γ and to ϒ → μ + μ − . It should be also emphasized that our prediction in partial width of hadronic transition χ b0 (3P) → χ b0 (1P)π π gives an opportunity to observe χ b0 (3P) in this channel.
Finally, we discuss the predicted total widths of P-wave bottomonia below the open-bottom threshold, where all of total [h b (n P)] with n=1,2,3 are about 80-90 keV. The total widths of χ b1 (n P) and χ b2 (n P) states with n = 1, 2, 3 are around 110-145 keV and 196-264 keV, respectively. However, the total widths of χ b0 (n P) states with n=1,2,3 are around 3.13-3.54 MeV, and these differences could possibly distinguish the P-wave spin-triplet χ bJ states from each other in experiments.

The 1D and 2D states
The ground states and first radial excitations of the D-wave bottomonia with orbital angular momentum L = 2 were estimated as 10150-10170 MeV and 10440-10460 MeV, respectively, which are below the BB thresholds. For the D-wave Table 6 Partial widths and branching ratios of annihilation decay and radiative transition and total widths for the 1P bottomonium states. Experimental results are taken from the PDG [52]. The theoretical results of the GI model [26] and the nonrelativistic constituent quark model [27] are summarized in the rightmost columns. The width results are in units of keV State Channels This work Expt. [52] G I [ 26] R e f . [ 27] Width Total 196 100 · · · · · · 180 100 122.84 100 Table 7 Partial widths and branching ratios of annihilation decay, radiative transition, and hadronic transition and total widths for the 2P bottomonium states. Experimental results are taken from the PDG [52]. The theoretical results of the GI model [26] and the nonrelativistic constituent quark model [27] are summarized in the rightmost columns. The width results are in units of keV State Channels This work Expt. [52] G I [ 26] R e f . [ 27] Width     Table 7 continued State Channels This work Expt. [52] G I [ 26] R e f . [ 27] Width   Table 8 Partial widths and branching ratios of annihilation decay, radiative transition, and hadronic transition and total widths for 3 1 P 1 and 3 3 P 0 bottomonium states. The width results are in units of keV states, the annihilation decays to three-gluon or two-gluon are generally suppressed than the S-wave. Therefore, all of the 1D and 2D states of bottomonium are expected to be narrow states. In 2010, the BaBar Collaboration observed the D-wave spin-triplet ϒ(1 3 D J ) states through decays into ϒ(1S)π + π − [22] and the member of J = 2 was confirmed with a significance of 5.8 σ , although other two states ϒ 1 (1 3 D 1 ) and ϒ 3 (1 3 D 3 ) have lower significances of standard deviations 1.8 and 1.6, respectively. In general, the experimental studies on the D-wave bottomonia are presently not enough, just as the total width and the branching ratio of typical decay channels are still unknown. For the ϒ 2 (1 3 D 2 ) state, the predicted masses by the GI model [26] and the nonrelativistic constituent quark model [27] are 10147 and 10122 MeV, respectively, which are much Table 9 Partial widths and branching ratios of annihilation decay, radiative transition, and hadronic transition and total widths for 3 3 P 1 and 3 3 P 2 bottomonium states. The width results are in units of keV Total 145 100 Total 264 100 lower than the measured value 10163.7 ± 1.4 MeV [52]. Our theoretical mass is 10162 MeV and such a consensus again shows an excellent description of bottomonium mass spectrum in this work. Furthermore, we also predict the mass splittings m[ MeV. Ignoring delicate differences of hyperfine spin splittings, the mass of the 2D bottomonium is estimated to be about 10450 MeV, which is in good agreement with that calculated in Ref. [26]. Experimental study of 2D states is also an interesting issues. Partial widths and branching ratios of annihilation decay, radiative transition, and hadronic transition and total widths for the 1D and 2D bottomonia are shown in Tables 10, 11 and 12, respectively. We want to emphasize that our results for most of decay channels are comparable with those given by the GI model [26]. However, there is still a significant difference for the M1 radiation transition widths calculated from the present work and the GI model, which is mainly due to the difference of the wave functions obtained under the quenched and unquenched quark models. This situation can be happen when comparing the result from the screened potential model and the GI model for the decays of the 1F and 1G states, which will be discussed later. For the 1D bottomonium states, we notice that the process ϒ 2 (1 3 D 2 ) → ϒ(1S)π + π − can determine the constant C 2 of Eq. (A11) and the predicted branching ratios of the identical hadronic processes for other two members with J = 1, 3 are also consistent with different theoretical models [26,27]. The radiative transition η b2 (1D) → h b (1P)γ may be an optimal channel to detect the η b2 (1D) state on account of the estimate of the branching ratio 96% by us. Accordingly, the radiative process of ϒ(1 3 D J ) → χ bJ (1 3 P J )γ with J = J + 1 is the dominant decay mode for spin-triplet states of 1D bottomonium, which has the estimated branching ratios of 44.1 %, 74 % and 90.5 % for J = 1, 2, 3, respectively. Surely, the nonforbidden radiative decay to P-wave χ bJ states is also significant for ϒ(1 3 D J ) state. There are similar decay behaviors in 2D bottomonium states, where the radiative transitions to P-wave bottomonium states are still dominant except for the ϒ 1 (2 3 D 1 ) state, whose half of the total width is contributed by the annihilation mode ggg. The analysis of other decay channels of 2D bottomonium states from Tables 11 and 12 can be summarized as follows: 1. The branching ratios of M1 radiative transitions of 2D → 1D and 2D → 2D with spin-flip are about 10 −5 ∼ 10 −7 , which indicate that it is difficult to observe these decay modes in experiments. 2. The E1 spin non-flip radiative transitions of 2D → 1F with J i = J f − 1 with the total angular momentum of initial (final) state J i( f ) can be used to study first Fwave bottomonium, due to the predicted partial width of 0.16 ∼ 1.4 keV and branching ratio of 0.7% ∼ 6%. 3. Compared with radiative transitions, the contributions of the dipion hadronic decays are much smaller, where the decays to S-wave η b or ϒ states and D-wave η b2 or ϒ J states occupy about 10 −3 ∼ 10 −1 % and 10 −10 ∼ 10 −2 %, respectively. 4. It is easy to see that the width of leptonic annihilation decay of D-wave ϒ state is much smaller than that of Swave ϒ states with three orders of magnitude by comparing Table 5 with Tables 10, 11, which can be used to distinguish two kinds of ϒ states with the same J PC = 1 −− in experiments.

Table 10
Partial widths and branching ratios of annihilation decay, and radiative transition and total widths for the 1D bottomonium states. Experimental results are taken from the PDG [52]. The theoretical results of the GI model [26] and the nonrelativistic constituent quark model [27] are summarized in the rightmost columns. The width results are in units of keV State Channels This work Expt. [52] G I [ 26] R e f . [ 27] Width    Table 11 Partial widths and branching ratios of annihilation decay, radiative transition, and hadronic transition and total widths for 2 1 D 2 and 2 3 D 1 bottomonium states. The width results are in units of keV Total 43.5 100

The 1F and 1G states
As the high spin states, the F-wave and G-wave bottomonia have no experimental signals at present. If these states can be experimentally observed, it could be a good confirmation for the theoretical calculation of the potential model. The predicted decay properties of 1F and 1G bottomonia for partial widths and branching ratios of annihilation decay, electromagnetic transition, and hadronic transition and the total widths are listed in Tables 13 and 14, respectively, where one can find a very interesting common decay feature. That is, the dominant decay modes of eight particles are all sorts of the radiative transition of 1F (J ) → 1D (J −1) or 1G (J ) → 1F (J −1) with almost more than 90% branching ratios, where subindices J / J − 1 represent the total angular momentum of a particle. This also indicates that experiments are likely to observe the first F-wave bottomonia via γ once all the ground states of D-wave bottomonium are experimentally established. In the same way, search for the first G-wave bottomonia also requires the confirmation of the F-wave states in experiment. This forms a chain-like relationship, which means that search for these high spin states is probably achieved step by step.
Our predicted mass values of 1F and 1G states are shown in Table 3, where the average mass values of spin triplet states for 1F and 1G are estimated as 10366 MeV and 10535 MeV, respectively. These are almost the same as those of spin singlet states. We hope that these predicted results will be helpful for the future experimental studies on high spin F and G-wave states. In addition, it can be seen from Tables 13, 14 that the 1F (1G) states could also transit to the specific P (D)wave states by emitting two light mesons ππ. However, par- Total 23.6 100 tial widths of these hadronic decay processes are calculated to be relatively small. The predicted total widths of the ground states of F and G-wave bottomonia are all about 20 keV, which are consistent with the estimates of the GI model [26].

Higher η b radial excitations
As the pseudoscalar partner of the S-wave ϒ states, the established members of the η b family are much less than those of the ϒ family. Based on this fact, we need to predict some theoretical results and simultaneously make a systematic study on the higher excited states of the η b family that will not only provide meaningful clues to experimentally search for them above thresholds but also further reveal their inner nature (Tables 15 and 16).
For the bb states above thresholds, they usually decay mainly into a pair of positive and negative bottom mesons or bottomed-strange mesons. The OZI-allowed strong decay behaviors could be studied by the QPC model, whose details Total 21.8 100  Table 17. For the input parameters of bottomonium states, we generally use the experimental masses and spatial wave functions of our modified GI model. However, for the predicted bottomonia, the principle of the input masses is consistent with Ref. [26]. The input masses of bottom and bottomed-strange mesons are summarized in Table 18, where two 1 + bottom mesons B(1P 1 ) and B(1P 1 ) are the mixture of 3 P 1 and 1 P 1 states, and the mixing angle θ 1P = − arcsin( √ 2/3) = −54.7 • is adopted as in the heavy quark limit [57][58][59]. Finally, the wave functions of bottom and bottomed-strange mesons are taken from the calculations of Ref. [60] and the quark mass parameters are given in Appendix A.
In this subsection, higher radial excitations from η b (5S) to η b (8S) are systematically studied, whose mass values are predicted in Table 3. Their corresponding decay properties including total widths and partial widths and branching ratios of OZI-allowed strong decay, annihilation decay and radiative transition are listed in Tables 19, 20, 21. Here, we make a comparison of the decay results of η b (5S) and η b (6S) pre-

Table 15
Partial widths and branching ratios of the OZI-allowed strong decay, annihilation decay, and radiative transition and the total width for ϒ(4S). Experimental results are taken from the PDG [52]. The theoretical results of the GI model [26] and the nonrelativistic constituent quark model [27] [52] sented in this work and given by the GI model [26], which generally reflect that these two quark models can achieve the same goal. The mass of η b (5S) state is calculated as 10810 MeV by the modified GI model. We give another estimate of the mass value of η b (5S) of 10870 MeV utilizing the results of theoretical mass splitting and the measured mass of ϒ(10860). At the same time, this mass estimate is also used as the model parameter of strong decay calculations. In Table 19, one can see that the total width of η b (5S) is predicted to be 28.    In comparison with the ϒ(nS) states with n=7, 8, which are likely to become a new research area for bottomonium physics in the future as their lower radial excited states have been observed experimentally, we also study the η b (7S) and η b (8S) state here. We predict their mass values as 11149 and 11289 MeV, and hyperfine mass splittings m(7S)\ m(8S) as 8\7 MeV, respectively. There are obviously large differences between the predicted mass by the GI model and modified GI model, which can be seen in Table 3. From Tables  20 and 21

Higher ϒ radial excitations
The S-wave ϒ states with J PC = 1 −− have always been the significant research field of bottomonium physics. Up to now, it is also the unique place to experimentally study some interesting effects above open-bottom thresholds, which includes ϒ(4S), ϒ(10860), and ϒ(11020). Therefore, further theoretical exploration for higher ϒ bottomonia is necessary. Additionally, the strong decay into two bottom mesons or bottom-strange mesons will be dominant for the three ϒ's mentioned above. Hence, their experimental information will determine the model parameter γ of the strong decay calculation in the QPC model.
There is no controversy in experimental measurements of the ϒ(4S) state, whose total width was first measured as 25 ± 2.5 MeV and 20 ± 2 ± 4 MeV by the CUSB [7] and CLEO [6] collaborations, respectively. On the other hand, BaBar's measurements in 2005 gave the similar resonance parameters [61] and the current average total width of PDG is 20.5 ± 2.5 MeV [52].
For the ϒ(5S) and ϒ(6S), however, the situation is not so optimistic since the early measured data and the recent measurements of resonance parameters are incompatible with each other. In short, after 2010, the Belle collaboration released consistent mass values of ϒ(10860) by the different production cross sections of e + e − → ϒ(nS)π + π − with n=1,2,3 and e + e − → bb, as well as e + e − → h b (n P)π + π − with n=1,2 [62][63][64]. Their values are larger than the theoretical calculation and early experimental values. In addition, the measurement of the total width by Belle also disagreed with the previous experimental conclusion of a broad state of ϒ(10860) [6,7]. Although there are a few experimental results of ϒ(11020), the center mass values from the BaBar [65] and the recent Belle collaboration [63,64] are consistently about 20 MeV smaller than the earlier experimental results [6,7]. It can be seen clearly from the PDG [52] that the total width of ϒ(11020) is measured to be about 30-60 MeV by integrating the present experimental information. At the same time, because the amount of data is not sufficient, the measurements of these two resonances by the analysis of the R value in the previous experiments are probably unconvincing to a great extent, where R = σ (e + e − → hadrons)/σ (e + e − → μ + μ − ). Therefore, we hope to see the more accurate experimental results for the ϒ(10860) and ϒ(11020) in the forthcoming Belle II experiment. It is worth   mentioning that in the measurements of the Belle [63], the resonance parameters from R ϒ(nS)π π and R b are completely consistent within the error range [63] although the application of a flat continuum in the R b fit brings some inconsistencies between the fitted amplitudes for R ϒ(nS)π π and R b . Furthermore, the measured mass values and total widths from R ϒ(nS)π π have too large statistical and systematic errors. , respectively. Based on the above consideration, in the following, the systematical study will be performed on the ϒ(nS) states above the thresholds. Higher bottomonia ϒ(7S) and ϒ(8S) states are also discussed in this subsection.
The total widths of the ϒ(10860) and ϒ(11020) are estimated to be 45 [26,27,31], where there is even a difference of 50 MeV on the prediction of ϒ(11020). Experiments have released the branching ratios of leptonic annihilation decays in the 10 −6 orders of magnitude both for ϒ(5S) and ϒ(6S) [52]. This is in accord with our theoretical predictions. Branching ratios of all six kinematically-allowed open bottom decay modes B B, B B * , B * B * , B s B s , B s B * s , and B * s B * s of ϒ(5S) have been experimentally measured as 5.5±1.0%, 13.7±1.6%, 38.1±3.4%, 0.5 ± 0.5%, 1.35 ± 0.32%, and 17.6 ± 2.7% [52], respectively. One can see that B * B * is the leading decay channel and the following B * s B * s and B B * are next-to-leading decay modes. This is not only contradictory with our calculations but most of other theoretical estimates of potential models [26,27,31]. The possible reasons are the following: firstly, as mentioned before, experimental measurement is not enough and one cannot rule out the latent error in experiment; secondly, the calculated partial width of strong decay is dependent on particle's real mass and a more accurate mass value may reduce the possible deviation from the phenomenological decay model. Finally, there may be some other physical effects and mechanisms in the strong decay process of ϒ(10860). In Ref. [66], authors just discussed the application of the Franck-Condon principle, which is common in molecular physics on the anomalous high branching ratios of B * s B * s versus B s B * s . For the ϒ(11020) state, there are no experimental data in the open-bottom decay channel at present. We predict that the dominant modes of ϒ (6S) are B B  *  , B B(1P 1 ), B B, B * B * with the corresponding branching ratios as 43 %, 21.6 % ,20.4 % and 11.5 %, respectively. These partial widths are quite different from the predictions by the GI model [26], though both of the predictions of total widths are contiguous, which are expected to be tested by the forthcoming Belle II.

ϒ(7S) and ϒ(8S)
In Ref. [26], authors did not study the properties of ϒ(7S) and ϒ(8S) in the framework of the GI model. Considering the importance of a screening effect for higher bottomonia, we employ the modified GI model with a screening effect to study the nature of the two particles in this work. In Table 3, we predict the masses of ϒ(7S) and ϒ(8S) as 11157 MeV and 11296 MeV, which are raised by 154 and 293 MeV, respectively, compared with the mass value of ϒ(11020) from Belle [63]. We suggest future experiments to search for these two bottomonia in the vicinity of their corresponding energies mentioned above. The decay information of partial widths and branching ratios of strong decay, annihilation decay, and radiative transition and the total widths for the  Table 20, we find that the combination of a vector meson B * and a P-wave bottom meson accounts for a large portion of the total width of ϒ(7S). Contributions from the B + B(1P) and bottomstrange meson modes can be almost ignored. Some typical ratios of partial widths are given by which can be tested by future experiments. In Table 21, the total width of ϒ(8S) is predicted to be 66.  Table 3, the mass values of higher radial P-wave bottomonia 5P and 6P states are estimated to be about 10940 and 11100 MeV, respectively, which are reduced by about 70 and 110 MeV compared to the predictions of the GI model, respectively. In this subsection, we will give a theoretical analysis of their decay properties including higher radial Pwave bottomonia 5P and 6P states in the framework of the modified GI model. We hope to provide useful information for the search for these bottomonia in future experiments.
The decay behaviors of 4P, 5P, and 6P bottomonia are given in Tables 22, 23, 24, 25, 26 in succession. For the 4P bottomonium states and even higher bottomonia with higher radial quantum number or higher orbital angular momentum, the screening effect begins to manifestly show its power in the description of corresponding decay properties because the theoretical calculations of various decay processes are directly dependent on the mass value and the wave function of the related hadrons. Therefore, in the following discussion, including the subsequent higher D, F, and G-wave states, we will perform a detailed comparison of the predicted decay properties with and without a screening effect. In Ref. [26], authors studied the P-wave bottomonia only up to 5P states by the GI model. Thus, we focus only on the comparison of 4P and 5P bottomonia, and our predicted decay properties of 6P bottomonia will be illustrated at the end.
According to the numerical results in Tables Table 23 Partial widths and branching ratios of the OZI-allowed strong decay, annihilation decay, and radiative transition and the total widths for 5 1 P 1 and 5 3 P 0 bottomonium states. The width results are in units of keV

D-wave states
In this subsection, we will discuss features of D-wave bottomonia. Among them, the D-wave vector bottomonium with J PC = 1 −− is quite interesting because unlike charmonium system, there is no clear signal to show the existence of Dwave vector states although many S-wave vector bottomonia have been discovered by experiments. Hence it is helpful for understanding this puzzle and further understanding the behaviors of bottomonia that we study the intrinsic properties of these missing ϒ 1 (n 3 D 1 ) states. There are only 1D and 2D bottomonia below the BB threshold. The mass of 3D states is located at around 10680 MeV according to our estimates, which is 120 MeV larger than the threshold. We also predict that the masses of 4D and 5D bottomonia are about Table 27 Partial widths and branching ratios of OZI-allowed strong decay, annihilation decay, and radiative transition and total widths for 3 1 D 2 and 3 3 D 1 bottomonium states. The width results are in units of keV  10880 and 11050 MeV, respectively. In Table 3, we notice that the predicted mass of ϒ 1 (4 3 D 1 ) is 10871 MeV, which is very close to the measured value of the ϒ(10860) state [63]. However, the possibility of a such candidate for ϒ(10860) can be basically excluded by the following analysis. Next, we discuss the decay properties of D-wave bottomonia up to 5D states in detail.
In Tables 27, 28, 29, 30, we list the numerical results of 3D and 4D bottomonium decays. Comparing ours with the calculation of Ref. [26], we can conclude 1. Compared with the GI model [26], almost all of our estimated partial widths of radiative transition for 3D and 4D bottomonia become smaller except for the electro-       (3D). Additionally, its branching ratio of annihilating to the leptonic pair + − is three orders of magnitude smaller than ϒ(4S) state. Hence, it is hard to detect this D-wave vector bottomonium in the final state events of   Although the mass of ϒ(10860) is compatible with our prediction of ϒ 1 (4 3 D 1 ) state in Table 3, we can safely rule out this allocation since the two orders of magnitude difference on the branching ratios of the leptonic pair decay + − between the theoretical and measured data. We predict that the main decay channels of ϒ 1 (4D)  The detailed decay information of 5D bottomonium states are presented in Tables 31 and 32. Our results indicate that the spin-singlet η b2 (5D) and triplet ϒ J (5D) with J = 1, 2, 3 are all broad states, whose predicted total widths are 110.7, 121.7, 101.6, and 86.0 MeV, respectively. We also notice that the dominant decay channels of η b2 (5 1 D 1 ), ϒ 2 (5 3 D 2 ) and ϒ 3 (5 3 D 3 ) are all B * B * and B B * and other relatively important decay modes are provided by the strong decay channels containing P-wave bottom mesons. It is very interesting to study the properties of ϒ 1 (5 3 D 1 ) because the mass of ϒ 1 (5D) is only 40 MeV larger than that of ϒ(6S) according to our prediction in Table 3, which is nearly 100 MeV smaller than that of the GI model. Our numerical results show that the ϒ 1 (5 3 D 1 ) is a broad state, for which there are ten openbottom decay modes. Furthermore, its main decay channels are B * B * , B B, B B * , B B(1P 1 ) and B B(1 3 P 2 ), and the corresponding branching ratios are 38.7 %, 16.4 %, 15.8 %, 14.8 %, and 7.59 %, respectively. The contributions of bottomstrange mesons are too low to exceed 1 %. Finally, we hope that these results can provide valuable clues for future experiments to search for more D-wave bottomonia.

F-wave states
In the following, we will focus on the higher F-wave bottomonia, i.e., 2F, 3F, and 4F bottomonia. The theoretical mass values of F-wave bottomonia are presented in Table 3. Their average mass values are the same as those of the spinsinglet states of 10609, 10812, and 10988 MeV with radial quantum number n = 2, 3, 4, respectively. From Fig. 2, it is also easily to find that the mass values of these particles Table 31 Partial widths and branching ratios of OZI-allowed strong decay, annihilation decay, and radiative transition and total widths for 5 1 D 2 and 5 3 D 1 bottomonium states. The width results are in units of keV    Table 32 Partial widths and branching ratios of the OZI-allowed strong decay, annihilation decay, and radiative transition and total widths for 5 3 D 2 and 5 3 D 3 bottomonium states. The width results are in units of keV  are close to the S-wave bottomonium states corresponding to radial quantum number n = 4, 5, 6, respectively.
In Tables 33, 34, 35, we present the numerical decay results of 2F and 3F bottomonia, which are also studied by Ref. [26] in the framework of the GI model. By analyzing and comparing our calculation with the GI model [26], we conclude 1. The radiative transitions of F-wave states also have a rule very similar to P-wave and D-wave bottomonia, which can be seen by comparing the predictions of the modified GI model and GI model. That is, our partial widths of radiative decays are smaller than those of the GI model on the whole except for the radiative processes 2F → 1D and 3F → 2D. 2. 2F bottomonia are states near the threshold, and their corresponding open-bottom channels are either B B or B B * , whose partial widths from the GI model are much larger than our predictions. We predict that the total widths of h b3 (2F) and χ b3 (2F) states are 0.413 and 0.543 MeV, respectively, which are about one-30th times the results of the GI model. The total widths of the other two particles χ b2 (2F) and χ b4 (2F) are calculated as 41.4 and 1.03 MeV, respectively. This result is about 2 to 3 times smaller than those of the GI model. From the above results, it can be easily noticed that the χ b2 (2F) is a relatively broad state but the other 2F states are all narrow states, which means that it is not difficult to distinguish the spintriplet with J = 2 of 2F bottomonia from the experiment. Some important radiative decay modes of these narrow states may also have large contributions to the total width. Hence, the processes h b3 (2F) → η b2 (2D)γ ,            which have an extra calculation of 3G bottomonium states compared to Ref. [26]. The hyperfine mass splittings among the spin-singlet and spin-triplet states of 2G and 3G are quite small and the average mass values of 2G and 3G bottomonia in Table  3 are 10747 and 10929 MeV, respectively, which are close to those of 4P and 5P bottomonia, respectively. The decay information on partial widths of several decay types and total

Comparison between the screening potential and coupled-channel quark models
The screening potential model adopted in this work and the coupled-channel quark model are typical unquenched quark models. In Ref. [50], the authors compared the results from the screening potential model and the coupled-channel model by taking a charmonium spectrum as an example, which, to some extent, reflects the equivalence between the screening potential model and the coupled-channel quark model. We notice that Ferretti and Santopinto [31] studied the bottomonium spectrum below 10.6 GeV by the coupledchannel quark model, which makes us have a chance to test the equivalence of two unquenched quark models further. In Fig. 4, we compare the results, which are from the present work obtained by the screening potential model and from a coupled-channel quark model [31]. In general, these two unquenched quark models present the comparable results, which again supports the conclusions of Ref. [50]. Surely, we also find the differences in the results from two unquenched quark models for the 3P and 1D states. This fact shows some differences between the screening potential and the coupled-channel quark models. After all, they are two approaches to phenomenologically describe the unquenched effects. Besides, the effects of nearby thresholds are important for some higher bottomonia, which cannot be normally reflected in the screening potential model. To explore the strength of this effect in the bottomonium system, we will take ϒ(4S) and ϒ(5S) states as an example to calculate their Fig. 4 A comparison of the results of our screening potential model and those obtained by a coupled-channel quark model [31] mass shifts in the coupled-channel model [67]. In Ref. [67], the inverse meson propagator can be presented as where m 0 is the mass of a bare state, and Re n (s) and Im n (s) are real and imaginary parts of the self-energy function for the n-th decay channel, respectively. The bare mass m 0 can be obtained from the GI model [46]. To obtain the mass of a physical state, the P −1 (s) can be expressed in the Breit-Wigner representation where m(s) 2 = m 2 0 + n Re n (s) and tot (s) = n Im n (s) m BW . The physical mass m BW can be determined by solving an equation m(s) 2 − s = 0. Furthermore, the imaginary of the self-energy function Im n (s) can be related to the amplitude of the 3 P 0 model by applying the Cutkosky rule and the corresponding Re n (s) can also be obtained by the dispersion relation. The parameter γ 0 in the 3 P 0 model can be fixed by fitting the experimental width of ϒ(4S) and ϒ(5S), which gives us γ 0 = 0.337 for the bottomonium system. By calculation of the coupled-channel model, we obtain the physical mass m BW (4S) = 10.592 GeV and m BW (5S) = 10.838 GeV. Comparing the above results and those from our screening potential model, we find that mass differences between two models are 20 MeV and 16 MeV for ϒ(4S) and ϒ(5S) states, respectively. These results show that the effects of nearby thresholds are essential for the bottomonium system, but on the other hand, they are also within our expectations. Generally, the contributions of effects of nearby thresholds to higher bottomonia should be systematically calculated in the coupled-channel model, which can be considered for further research in the future.
In addition, when comparing our results of average mass of these 3P states with those from other models [26,27,31], we notice that our results are in agreement with those given in Refs. [26,27,31]. In Ref. [31], the authors specified the coupled-channel effect on the mass splittings for the χ bJ (3P) states. If further focusing on the mass splittings of χ bJ (3P) states, there exist differences among different model calculations [26,27,31]. It is obvious that the theoretical and experimental study on this mass splittings of χ bJ (3P) states will be an intriguing research topic in future.

Summary
As we all know, the screening effect usually plays a very important role for highly excited mesons. It affects the mass values, wave functions of mesons and hence, estimates of corresponding decay behaviors. Motivated by the recent studies of bottomonium properties in the framework of the Godfrey-Isgur relativized quark model [26], we have performed the most comprehensive study on the properties of bottomonium family by using the modified Godfrey-Isgur model with a screening effect. We have studied radiative transition, annihilation decay, hadronic transition, and OZI-allowed two-body strong decay of bb states. Our calculated numerical results indicate that our predictions for the properties of higher bottomonia are quite different from the conclusion of GI model. Hence, we also expect that this work could provide some valuable results to the future research of bottomonium in experiment.
Our work in this paper can be divided into two parts, study of mass spectrum and study of decay behaviors of bottomonium states. Furthermore, we have focused our main attention on the prediction and analysis of higher bottomonia due to significant reflection of the screening effect. Firstly, we have taken advantage of the measured mass of 18 observed bottomonium states in Table 2 to fit eight undetermined parameters of the modified GI model in Table 1. It can be found that our theoretical mass values have been greatly improved compared to those of the GI model. At the same time, our results have been well matched with experimental results. Based on the above preparation, the predicted mass spectrum of bottomonium states has been given in Table 3. It is interesting to note that the mass values of ϒ(7S) and ϒ(8S) are predicted as 11157 and 11296 MeV, respectively, which are higher than the measured mass of ϒ(11020) only by 154 and 293 MeV, respectively.
Classifed in L, the decay properties of bottomonium states have been separately discussed in accordance with mass values above and below the open-bottom threshold. We have found that a screening effect is weak for the decay behav-iors of the most bottomonium states below the threshold, whose estimates are similar to those of the GI model. For the higher bottomonium states above the threshold, the screening effect has become important. We have obtained fairly inconsistent conclusion on characteristic decay behaviors of bottomonium mesons between the GI model with and without screening effects. Here, we have provided results to check the validity of our model in future experiments.
In the following years, exploration for higher bottomonium states will become a major topic in the future LHCb and forthcoming Belle II experiments. Until then, the highly excited states that are still missing are likely to be found. Moreover, some hidden experimental information of observed bb states can be further perfected. In this work, we have provided abundant theoretical information for higher bottomonia, which is helpful for piloting experiments to search for these missing bottomonium states.

Radiative transitions
The radiative transitions of a heavy quarkonium are important in the sense that radiative decays not only are main decay channels of some particles below the open-flavor threshold, but help us better understand the inner structure of a quarkonium, i.e., wave functions and interactions of QQ. For the E1 transition process, n 2S+1 L J → n 2S+1 L J + γ of charmonium, the partial width is given by [68], with where the e b is a bottom quark charge in units of |e| and α is a fine-structure constant, ω is an emitted photon energy and ψ f |r |ψ i is the transition matrix element which has the integral form ∞ 0 R n L (r )r R nL (r )r 2 dr. Here the radial wave function R nL (r ) is obtained from the modified GI model using parameters listed in Table 1, and they are the same in the calculation of M1 radiative transitions.
The partial width of the M1 radiative transion with the spin flip from the initial state n 2S+1 L J to the final state n 2S +1 L J can be written as [69] where m b is the mass of a bottom quark, j 0 ( ωr 2 ) is the spherical Bessel function and other parameters have been defined above. The results of radiative transition will be discussed in Sects. 3 and 4.

Annihilation decays
The annihilation decay of a heavy quarkonium is important especially for those low excited states below the open-flavor threshold, where the annihilation decay to gluons is dominant. Furthermore, the decay mode is generally available in experiments, in which a vector meson 3 S 1 or 3 D 1 generates lepton pairs, i.e., e + e − , μ + μ − , and τ + τ − . The measured branching ratios of lepton-pair decays can also be used to judge whether experimental XY Z exotic states are treated as conventional mesons or not because this ratio is usually very small for the multiquark states [70]. It is important to study the annihilation decays of bottomonium states into gluons, light quarks, leptons, and photons in this paper, and the general formulas for annihilation decays of a heavy quarkonium have been extensively studied by using perturbative QCD methods [71][72][73][74][75][76][77][78][79][80][81][82][83]. The most important feature of annihilation decays is that the probability of annihilation is related to the zero point of the meson wave function or its n-th order derivative where n corresponds to the orbital angular momentum of a meson, i.e., n = L. For the lepton pair annihilation decay of 3 S 1 or 3 D 1 bottomonium state, this process occurs via a virtual photon in the tree level, and their width expression with the first-order QCD radiative corrections can be found in Refs. [77,83]. The formulas for annihilation processes n 3 S 1 → ggg, n 3 S 1 → γ γ γ and n 3 S 1 → ggγ and the annihilation modes of P-wave states including QCD radiative corrections are given in Ref. [77]. The general expressions for the decays into two gluons and two photon of spin singlet state with an arbitrary total angular momentum are given in Ref. [78]. Finally, for the bottomonium states with the higher angular momentum, the authors of Refs. [79,80] have given complete expressions for the annihilation decay to three gluons for the D-wave spin-triplet states. In Refs. [81,82], the annihilation decay of the F-wave spin-triplet states into two gluons was also studied. Fortunately, the authors of Ref. [26] have summarized all these lowest-order annihilation decay formulas with the first order QCD corrections of a heavy quarkonium from S-wave to G-wave states. Hence, we do not specifically list these formulas here. It is worth noticing that according to the formula of Ref. [77] and our model input, we recalculate each coefficient of first order radiative corrections for the processes χ b0 (n P) → gg and χ b2 (n P) → gg. The radiative correction terms for χ b0 (5P, 6P) and χ b2 (5P, 6P) states are also added in this work. We find that the coefficient C(n P) in the radiative correction term C(n P)α s /π of the above processes are modified only for the radial excited states, χ b0 (3P) and χ b2 (2, 3, 4P). The corresponding constants C(n P)'s are 10.4, 0.89, 1.59 and 2.10, respectively. The correction constants for χ b0 (5, 6P) and χ b2 (5, 6P) are also estimated to be 10.6, 10.7, 2.50 and 2.85, respectively. Finally, we need to emphasize that some of parameters in the annihilation decay calculations are given by m b =5.027 GeV, α s (bb)=0. 18 and α = 1/137.

Hadronic transitions
The hadronic transition of a heavy quarkonium usually refers to the release of a light hadron when the QQ state moves to a lower energy level, which is very critical in the search for some bottomonium particles below open-bottom thresholds. In this work, the hadronic transitions of bottomonia will be studied in the framework of the QCD multipole expansion method, in which the hadronic transition is described as first emitting one gluon from a heavy quark to form the intermediate hybrid state with a color octet QQ pair and then recombine themselves into a light hadron with another emitted gluon via the hadronization process. Since the mass difference of a heavy quarkonium between before and after the transition is usually small, the wavelengths of emitted gluons are generally far larger than the radius of the heavy quarkonium. Similar to electromagnetic radiation, gluon field can be treated in the multipole expansion form in this situation, which was first proposed in Ref. [84]. Next, we briefly introduce the QCD multipole expansion method and more details can be found in the review article Ref. [85].
The quark and gluon fields are assumed to be ψ(x) and A a μ (x), respectively and are transformed as by the introduction of an operator U (x, t) defined as where P is the path-ordering operator and x is the mass center of a quarkonium. In fact, this transformation indicates that (x, t) dressed by gluons plays a role of a constituent quark in the effective Lagrangian of system which can be obtained in Ref. [86]. As mentioned above, the process of launching gluons of heavy quarkonium can be treated by the multipole expansion method where in the zero position, the emitted transformed gluon field A μ (x, t) can be expanded as A (x, t) = − 1 2 x × B(0, t) + · · · (A7) on the basis of effective Lagrangian, the corresponding Hamiltonian can be derived as the follow form [86] H e f f where the H QC D is part of kinetic and potential energy of heavy quarkonium field which is not the simple Hamiltonian of free field but already contains relatively strong interaction, and the H (1) QC D and H (2) QC D are usually seen as a perturbation which consist of the interactions of color charge, colorelectric dipole moment and color-magnetic dipole moment as well as higher order multipole momentum of quarkonium field, respectively. The general formula of S-matrix element in the QCD multipole expansion is given in Ref. [87]. Here, we only focus on the spin-nonflip ππ transitions which were dominated by double electric-dipole (E1-E1) transitions since other transitions including spin-flip ππ processes where E1-M1 transition is main and spin-nonflip η decays which are contributed by E1-M2 and M1-M1 transitions are usually suppressed compared with the E1-E1 transitions. Starting from the general formula of S-matrix, the amplitude of spin-nonflip ππ transitions can be written as [86][87][88] where | I and F h| are initial quarkonium and final quarkonium and light hadron, respectivrly.x is the separation of heavy quark and anti-quark and (D 0 ) bc = δ bc ∂ 0 − g s f abc A a 0 . After inserting a complete set of intermediate states, this transition amplitude can be divided into two parts which are a heavy quark MGE (multipole gluon emission) factor and an H (hadronization) factor, respectively and the concrete form is given by [88] M E1E1 = i g 2 As for the MGE factor which has two electric dipole factors, the initial state | I first transforms to the intermediate vibrational state K L| formed by the color-octet quarkonium and gulon called as a hybrid state. Since this three-body bound state cannot be solved by the QCD, we use the quark confining string (QCS) model [89][90][91] as a viable approach to calculate the intermediate hybrid, which will be mentioned later. The part of an MGE factor can be calculated by appllying the eigenvalue and the wave function of an intermediate hybrid, initial and final quarkonium states. The H factor ππ|E a k E a l |0 clearly reflects the process of two emitted gluons transforming to the light hadrons after hadronization. This H factor is highly nonperturbative due to the low scale of energy. Hence, this matrix element can not be also directly obtained by the QCD, however, a phenomenological approximation can be given by using the partially conserved axialvector current and soft pion theorem [88,92]. Based on the above treatments, the final transition rate is given by [88] (φ I → φ F + ππ) =δ l I l F δ J I J F G|C 1 | 2 − 2 3 H |C 2 | 2 |A 1 | 2 + (2l I + 1)(2l F + 1)(2J F + 1) where C 1 and C 2 are parameters which are determined by the processes of ϒ(2S) → ϒ(1S)π π and ϒ(1D) → ϒ(1S)π π , respectively. The symbols l I (l F ), J I (J F ) are the orbital and total angular momentum of initial (final) state, respectively and the spin s does not change after the reaction. The f L P I P F I F has the following structure f L P I P F where M K L and R K L (r ) are the mass and radial wave function of the intermediate state, respectively. The phase-space factors G and H are written as with K given by The intermediate hybrid states can be described by the quark confining string (QCS) model [89][90][91], in which we consider that the quark and anti-quark are connected by an appropriate color electric flux tube or string. If the string is in the ground state, the system of a quark-antiquark pair is a meson where the string corresponds to the strong confinement interaction. The vibration of the string means a new state with gluon excitation effects, which is composed of the excited gluon field and quark-antiquark pair, i.e., the socalled hybrid state. For this vibrational mode, assuming both ends of a string are fixed because of too heavy quark masses, then the effective vibrational potential can be given by [90] V n (r ) = σ (r )r 1 − where d is the correction of finite heavy quark mass and n indicates the excitation level. The α n related to the shape of the vibrational string [90] is taken as √ 1.5, which is consistent with Ref. [27] and is insensitive to our mass spectrum of hybrid states.
The potential of a hybrid meson can be expressed as [91] V hyb (r ) = V G (r ) where V G (r ) is one-gluon exchange potential and V S (r ) is a color confining potential. It is easy to see that the above potential becomes a general QQ interaction when n = 0 for the vibrational potential V n . For theoretical self-consistency, forms of V G (r ) and V S (r ) are taken from our modified GI model and due to a screening effect, the effective string tension σ (r ) is not a constant but rather a function of a distance r of Q andQ. The specific expressions of potentials V G (r ) and V S (r ) can be written as Solving the Schrödinger equation for a hybrid meson, one obtains the mass spectrum and corresponding wave function of a hybrid state, which are used to calculate the width of hadronic transition by Eq. (A11). Nevertheless, we have to emphasize that the QCD multipole expansion is dependent on the inputs and has its own error due to theoretical uncertainties. Hence, we should regard the calculated width of hadronic transition as rough estimates rather than precise results. In addition, considering that there may be a more complex mechanism of hadronic transition for highly excited states, we focus on the hadronic transition of the bottomonium states only below open-bottom threshold in this work. The numerical results of hadronic decay will be discussed in Sect. 3 and it should be noted that the GI's results of hadronic transition Ref. [26] are derived from the reduced matrix elements, which are obtained by measured transition rates rather than direct calculations. Here, we adopt the direct QCD multipole-expansion method to calculate a partial width of hadronic transition of bottomonium states and it is also useful to make a comparison with the results of Ref. [26]. 4. Two-body OZI-allowed strong decays Quark-Pair Creation (QPC) model is applicable to the calculation of OZI allowed hadron strong decays. This model is proposed by the Micu [93] at the earliest in 1968 and it has been further developed by the Orsay Group [94][95][96][97] which is one of the most popular phenomenological method to deal with the OZI allowed strong decays and has been greatly used in the calculation of strong decay. The model assumes that a created quark-antiquark pair qq from the vacuum is a 3 P 0 state which has spin-parity J PC = 0 ++ , so the model also known as 3 P 0 model. In the following, we will briefly introduce this model. For the OZI allowed strong decay process A → B + C, the transition operators T can be expressed as T = −3γ m 1m; 1 − m|00 dp 3 dp 4 δ 3 (p 3 + p 4 ) where Y lm (p) = p l Y lm (θ p , φ p ) is a solid spherical harmonic function and b † 3 (d † 4 ) is quark (antiquark) creation operator, φ 34 0 = (uū + dd + ss)/ √ 3 and ω 34 0 are SU(3) flavor and color wave function of vacuum quark pair respectively, and the dimensionless parameter γ describes the strength of creating a quark-antiquark pair from the vacuum. γ value for ss pair creation is generally more than a factor of 1/ √ 3 compared to that of uū/dd pair creation. The reason of the existence of factor 1/ √ 3 is in order to show the SU(3) symmetry breaking [51,[94][95][96][97][98]. The transition matrix of decay process can be written as ×φ D ω D dp 1 dp 2 δ 3 (P D − p 1 − p 2 ) where the front factor E is particle energy and where J and L are the total and orbital angular momenta between final state B and C respectively and P = P B . Finally, the partial width of the A → BC can be written as where m A is the mass of the initial state A. In addition, in the calculations of strong decay, the constituent quark mass of bottom, up/down and strange quark are taken as 5.027, 0.22 and 0.419 GeV, respectively.