Higher bottomonium zoo

In this work, we study higher bottomonia up to the $nL=8S$, $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\ge 4$ states. Our theoretical results are valuable to search for more bottomonia in experiments, such as LHCb, and forthcoming Belle II.


I. 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 experimentally 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 Pwave 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 spin-triplet Υ(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 nonperturbative 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 Sec. II, 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 Secs. III and IV, respectively, where plentiful predictions will be given. In Sec. V, we compare the results between the modified GI model and a coupled-channel quark model. We make a summary in Sec. VI. 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 GI model will be given in Appendix B.

II. SPECTRUM
A. 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 one-gluonexchange 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 G(r) = − 4α s (r) 3r (2) 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 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 asG 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 I. 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  Table I. 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 with 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.

B. Mass spectrum
Due to the introduction of a new parameter µ in our modified GI model, we need to combine experimental information

Mass MeV
BB 9500 10 000 10 500 11 000 11 500 Mass spectrum of bottomonium. Here, bottomonia are classified by the quantum number 2S +1 L J , from left to right successively in each category , black, pink, and purple lines stand for our results from the modified GI model, the predictions of the GI model [26,46], and experimental values taken from PDG [52], respectively. The position of the open-bottom threshold is identified by the dashed line. In addition, the notation J means all the total angular momentum numbers of triplet states such as J = 0, 1, 2 for the P-wave states and J = 1, 2, 3 for the D-wave states and so on. 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

R r
where the V Exp (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  [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 I, where parameters of the GI model are also given for comparison.
In Table II, 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 nonrelativistic 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 II. 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 III. 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 thoroughly 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 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 II. There is a more interesting physical quantity, i.e., the hyperfine mass splitting between the spin-singlet and spintriplet states ∆m(nS ) = m[Υ(nS )] − m[η b (nS )], 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  VI. 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 VI, 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 Section.  II. 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.
B. 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 de-cays, annihilation decay, and hadronic transitions and the total widths for Υ(1S ), Υ(2S ), and Υ(3S ) are listed in Table VII. 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 firstly start from the analysis of the common decay properties of Υ(1S ), Υ(2S ), and Υ(3S ). From Table VII, 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 VII are obtained by combining measured total widths and branching ratios of the PDG [52]. The M1 radiative transition Υ(1S ) → η b (1S )γ is the unique electromagnetic decay of Υ(1S ) state, which has no experimental information until now. The calculations of three models, our work and Refs. [26,27] give a consistent estimate, 0.01 keV. Numerous radiative decay modes are opened for Υ(2S ) including the transition to χ bJ (1P), all of which are measured by experiment. By comparison, we can see that the experimental data of radiative transition of Υ(2S ) can be well reproduced by our model and that of Refs. [26,27] except for the decay channel η b (1S )γ. It should be mentioned that the calculations of branching ratios by the nonrelativistic constituent quark model [27] adopt the measured total widths, which is the reason why our theoretical results are smaller but the estimates of partial widths are close to the vlaues of Ref. [27]. At last, the hadronic transition Υ(2S ) → Υ(1 3 S 1 )ππ has been used to fix the unknown model parameter C 1 in the QCD multipole expansion approach, Eq. (A11) in Appendix A.
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. 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. The P-wave states χ bJ (1P) and χ bJ (2P) with J = 1, 2, 3 were firstly discovered in the search for the radiative processes Υ(2S ) → χ bJ (1P)γ [10] and Υ(3S ) → χ bJ (2P)γ [11] in 1982. The corresponding spin-singlet partners h b (1P) and h b (2P) VI: Partial widths and branching ratios of annihilation decay, radiative transition, and hadronic transition and total widths for the S -wave η b below the open-bottom thresholds. Experimental results are taken from the PDG [52]. The theoretical results of the Godfrey-Isgur relativized quark model [26] and the nonrelativistic constituent quark model [27] are summarized in the rightmost columns. The width results are in units of keV.

This work
Expt. [52] GI [26] Ref. [27] State Channels Width   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].
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], VIII: 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.
In addition, it is worth noting that our results of critical hadronic decay channel of χ bJ (2P) → χ bJ (1P)ππ with nonflip J and h b (2P) → h b (1P)ππ calculated by the QCD multipole 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 Dwave 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 X and XI, 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 demon-strated 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 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 modes η b (1S , 2S , 3S )γ and h b (1P)ππ are also important for h b (3P) states. Especially, the process of h b (3P) → η b (3S )γ has a predicted branching ratio of 10.4%. If h b (3P) 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 empha- IX: 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.

D. 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 states, the annihilation decays to three-gluon or two-gluon XII: 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.

This work
Expt. [52] GI [26] Ref. [27] State 96 · · · · · · 24.9 91.5 17.23 96.58 Total 25.3 100 · · · · · · 27.2 100 17.84 100 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 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 Table XII and  Tables XIII-XIV, 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 op- timal channel to detect the η b2 (1D) state on account of the estimate of the branching ratio 96% by us. Accordingly, the radiative process of Υ( 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 non-forbidden 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 chan- Total 25.9 100 nels of 2D bottomonium states from Tables XIII and XIV 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 VII with Tables XII-XIII, which can be used to distinguish two kinds of Υ states with the same E. 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 pre-  dicted 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 XV and XVI, 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 these radiative processes of h b3 (1 1 F 3 ) → η b2 (1 1 D 2 )γ and χ bJ (1 3 F J ) → Υ J−1 (1 3 D J−1 )γ 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 Fwave 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 III, where the average mass values of spin triplet states for 1F and 1G are estimated as 10366 MeV and 10535 MeV, respectively. These are al-most 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 XV-XVI that the 1F (1G) states could also transit to the specific P (D)-wave states by emitting two light mesons ππ. However, partial widths of these hadronic decay processes are calculated to be relatively small. The predicted total widths of the ground states of F and Gwave bottomonia are all about 20 keV, which are consistent with the estimates of the GI model [26].

IV. ANALYSIS AND PREDICTIONS OF HIGHER BOTTOMONIA
A. 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.
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 are given in Appendix A. In this work, the parameter γ for the bottomonium system in the QPC model can be determined by fitting to experimental widths of Υ(4S ), Υ(5S ), and Υ(6S ) states and the fitted results are shown in Table V. For the input masses of bottomonium states, the experimental masses will be considered. If the members of a bottomonium multiplet are partially discovered by experiments, the input masses of the missing states can be obtained by combining the measured mass and predicted splitting. For the predicted bottomonia, we calculate the behaviours of strong decays by utilizing our theoretical masses in Table III as an input. Besides the input masses, the spatial wave functions of bottomonium states can be directly derived from our modified GI model. The input masses of bottom and bottomed-strange mesons are summarized in Table IV, 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 III. 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 XIX-XXI. Here, we make a comparison of the decay results of η b (5S ) and η b (6S ) presented 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 XIX, one can see that the total width of η b (5S ) is predicted to be 28.4 MeV, which can mainly decay into the BB * , gg, B * s B * s , B s B * s , and B * B * , and corresponding branching ratios are 71.5 %, 11.5 %, 10.2 %, 3.91 %, and 2.71 %, respectively. The dominant radiative decay mode is h b (4 1 P 1 )γ with a partial width of 33.5 keV. We should notice that the notation BB * refers to BB * +B * B , B * B * to B * B * , and so on. The predicted mass of the η b (6S ) is 10991 MeV, which should also be improved compared to the estimates of the GI model, just like Υ(6S ). The total width of η b (6S ) is estimated to be 20.4 MeV, which is about the same as that of η b (5S ). The dominant decay channel of η b (6S ) is still BB * with a branching ratio of 73 % and other decay modes gg, BB(1 3 P 0 ), B s B * s , and B * B * also have large contributions to the total width. In general, our decay results indicate that the η b (5S ) and η b (6S ) are most likely to be detected in the BB * mode.
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 III. From Tables  XX and XXI 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. Addi- XVII: 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] are summarized in the rightmost columns. The width results are in units of keV.
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 (nP)π + π − 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. So in the calibration of parameter γ, we select the latest Belle's measurements of the process e + e − → bb, which give the resonance parameters for 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. XVIII: Partial widths and branching ratios of the OZI-allowed strong decay, annihilation decay, and radiative transition and the total widths for Υ(10860) and Υ(11020). 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  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 BB * , BB(1P 1 ), BB, 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  III, 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 Υ(7S ) and Υ(8S ) are shown in Tables XX and XXI, respectively. The total width of Υ(7S ) is estimated to be 101.4 MeV, which means a broad state. There are thirteen open-bottom modes, among which the important decay channels are B * B(1 3 P 2 ), BB * , B * B * , B * B(1P ′ 1 ), B * B(1P 1 ), and BB with the corresponding partial widths, 28.1, 22.0, 20.4, 9.26, 9.03, and 6.79 MeV, respectively. From Table XX, 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 bottom-strange meson modes can be almost ignored. Some typical ratios of partial widths are given by which can be tested by future experiments. In Table XXI, the total width of Υ(8S ) is predicted to be 66.6 MeV and the dominant decay channels are B * B * and BB * with branching ratios 40.9 % and 31.0 %, respectively. The B * B * and BB * modes are excellent decay channels to detect the Υ(8S ) bottomonium state in the future experiments. All the decay modes BB, BB(1P ′ 1 ), and BB(1 3 P 2 ) have a consistent partial width of about 4−5 MeV, which almost occupies the remaining contributions to the total width. It is difficult to experimentally observe the configuration of bottom-strange mesons in the Υ(8S ) decay.
C. P-wave states From Fig. 2, the first above-threshold P-wave bottomonia are 4P states, which include spin-singlet h b (4P) and spintriplet χ bJ (4P) with J = 0, 1, 2 and exceed the BB threshold by about 200 MeV according to our estimates. These abovethreshold particles have not yet been found by experiments. In Table III, 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 P-wave 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 XXII-XXVI 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 XXII-XXIV and Ref. [26], we conclude the followings. Although the difference in total width is not conspicuous for most of the 4P and 5P states, but the predicted dominant decay modes of 4P and 5P states from two different models are almost all opposite, which also illustrates that the influence of the screening effect on the higher bottomonia is quite significant. Thus the present work should be of great value for revealing nature of bottomonium. The dominant decay channels of the h b (5P) state are BB * and B * B * , and corresponding branching ratios are 75.7 % and 21.1 % respectively. We predict that the total width of χ bJ (5P) is about 40 ∼ 50 MeV, and decay mode BB and B * B * are critical for the χ b0 (5P) state and the dominant decay channels of χ b1 (5P) are BB * and B * B * , furthermore, for the χ b2 (5P), there are three important decay modes BB * , B * B * and BB which all have a great contribution for its total width. The above conclusions could be examined by experiment in the future.
The predicted average mass of 6P bottomonium states by modified GI model is 11099 MeV, which is about 80 ∼ 100 MeV above the experimental mass of the observed Υ(11020) state. Hence, more strong decay channels of 6P states are opened. From Table XXV and XXVI, Our results indicate that the h b (6P) and χ bJ (6P) all are broad bottomonium mesons because of the predicted total width of about 107 ∼ 140 MeV. Additionally, BB * , B * B * , and B * B(1 3 P 2 ) just simultaneously are the main decay channels for spin-singlet h b (6P) and triplet χ bJ (6P) with J = 1, 2, and the sum of their contributions all are more than 70 % to the total decay width. The decay modes B * B(1 3 P 2 ), B * B * and BB are dominant for the χ b0 (6P) state with branching ratios of 28.8 %, 25.7 % and 24.5 %, respectively. A common decay feature of 6P bottomonium states can clearly be seen, that is, the role of mode B * B * and B * B(1 3 P 2 ) are quite considerable.

D. 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 10880 and 11050 MeV, respectively. In Table III, 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 XXVII-XXX, 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 electromagnetic processes of 3D → 1P and 4D → 2P, which become larger from 20 % to 300 % range.
2. Similar to the situation of χ b0 (4P), our prediction of the total width for the Υ 1 (3D) is 54.1 MeV, which is largely different from an estimate of 103.6 MeV by the GI model. Overall, the 3D bottomonia are broad resonances, where the predicted total widths of other three particles η b2 (3D), Υ 2 (3D) and Υ 3 (3D) are 143.0, 96.3, and 223.8 MeV, respectively. Our results of 3D and 4D states on the partial widths and branching ratios of strong decay channels are still quite different from those of the GI model. The strong decay modes of η b2 (3D) are only B * B * and BB * with close predicted ratios, which accounts for almost all the contributions to the total width. For the Υ 1 (3D), the dominant decay mode is  (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 µ + µ − . This situation is also applied to the Υ 1 (4D) and Υ 1 (5D) states. Finally, one can easily find that the modes B * B * and BB * are dominant both for the Υ 2 (3D) and Υ 3 (3D) at the same time. 3. The total widths of 4D bottomonia including the spin singlet η b2 (4D) and three triplets Υ J (4D) are estimated to be about 70 ∼ 90 MeV. Our results are generally 10 ∼ 20 MeV larger than those of the GI model. The η b2 (4D) mainly decays into BB * and B * B * with the branching ratios 63.2 % and 33.4 %, respectively. The corresponding partner Υ 2 (4D) has the similar decay behavior. Although the mass of Υ(10860) is compatible with our prediction of Υ 1 (4 3 D 1 ) state in Table III, 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) are B * B * , BB and BB * with corresponding     annihilation widths are too small as mentioned before. Hence, this can also explain why the S -wave Υ's are experimentally found in succession but the D-wave Υ 1 's always have no movement in experiments.
The detailed decay information of 5D bottomonium states are presented in Table XXXI and XXXII. 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 BB * 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 III, 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 * , BB, BB * , BB(1P ′ 1 ) and BB(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.

E. 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 III. 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 are close to the S -wave bottomonium states corresponding to radial quantum number n = 4, 5, 6, respectively.
In Tables XXXIII-XXXV, 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 BB or BB * , 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 states are all narrow states, which means that it is not difficult to distinguish the spin-triplet 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)γ, χ b3 (2F) → Υ 2 (2D)γ,    Tables XXXVI and XXXVI and the predicted total widths  of 4F states are located at near 70 MeV. In their openbottom decay channels, the contributions from bottom-strange mesons are all small, among which the largest is not more than 3 %. The dominant decay modes of h b3 (4F) and χ bJ (4F) with J = 2, 3 are B * B * and BB * , whose the sum of branching ratios is more than 80 %. Only the mode B * B * with branching ratio 92.1 % governs the χ b4 (4 3 F 4 ) state. The BB channel is also important for the χ b2 (4 3 F 2 ) state but is negligible for the χ b4 (4 3 F 4 ) state on account of the branching ratios 17.5 % and 1.27 × 10 −4 , respectively.
F. G-wave states G-wave bottomonia are typical high spin states, which are difficult to be experimentally found in recent years. However, we still give our predictions of decay behaviors for the higher G-wave bottomonia with radial quantum number up to 3, 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 III 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 widths of 2G and 3G bottomonia are shown in Tables XXXVIII-XL. For the radiative transition of 2G states, we obtain the behaviors similar to P, D, F-wave states. Hence, we no longer explain it here. Furthermore, the total widths of 2G bottomonia are estimated to be about 110 ∼ 160 MeV, which means that they are all broad states although these values are significantly small relative to those of the GI model. Since the mass values of 2G bottomonia do not reach the threshold line of including a P-wave bottom meson, the main decay modes of 2G states are BB, BB * and B * B * . The mode BB * is dominant for the Υ 3 (2 3 G 3 ) and Υ 4 (2 3 G 4 ) with branching ratios 50.9 % and 62.5 %, respectively. The mode B * B * with a branching ratio 76.8 % is critical for the Υ 5 (2 3 G 5 ). Finally, the modes BB * and B * B * have the same importance for the η b4 (2 1 G 4 ) state.
Similarly, the 3G bottomonium states cannot decay into Pwave bottom mesons on account of the mass values. In Tables XXXIX and XL, the predicted total widths of η b4 (3G) and Υ J (3G) with J = 3, 4, 5 states are 53.3, 39.8, 50.4, and 67.5 MeV, respectively. All the contributions of bottom-strange mesons to the total widths of 3G states occupy about 4 %. In addition, the strong decay behaviors of channels BB, BB * , and B * B * of 3G bottomonia are almost exactly the same as the case of 2G bottomonia, which also indicates that the dominant decay modes of each 2G bottomonium state and corresponding radially excited 3G state are similar to each other.  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 coupled-channel 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 mass shifts in the coupled-channel model [67]. In Ref. [67], the inverse meson propagator can be presented as P −1 (s) = m 2 0 − s + n (ReΠ n (s) + iImΠ n (s)) , 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 selfenergy 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 coupledchannel 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.

VI. 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 II to fit eight undetermined parameters of the modified GI model in Table I. 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 III. 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 behaviors 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.
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) 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 so-called 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 − 2nπ 2nπ + σ(r)[(r − 2d) 2 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) + V S (r) + [V n (r) − σ(r)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 Secs. III 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].

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 ) ×Ψ D nLM L (p 1 , p 2 )|q 1 (p 1 )q 2 (p 2 ) , where the front factor E is particle energy and Ψ D nLM L (p 1 , p 2 ), χ D , φ D and ω D donate spatial, spin, flavor and color wave function of meson D, respectively, and LM L ; S M S |JM J is Clebsch-Gordan coefficients. For the spatial wave function of initial and final states, we use the exact eigenfunctions by solving Schrödinger equation in potential models rather than simple harmonic oscillator (SHO) wave function. Combined Eqs. (A23)-(A24) and Eq. (A25), the decay amplitude M M J A M J B M Jc can be derived.
For the convenience of experimental measurement the decay amplitudes could be related to the helicity partial wave amplitudes by Jacob-Wick formula [100] M 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. momentum-independent correction to the one-gluon exchange potential, which can be obtained from Eq. (7) and their specific expressions arẽ respectively, with constant parameters ǫ t and ǫ c . The last term V so is the spin-orbit coupling and it includes vector and scalar spin-orbit potentials, namely, whereṼ ii are also the momentum-dependent corrections for the vector and scalar spin-orbit interactions, respectively, which are given bỹ where the subscripts i( j)=1, 2 denote bottom and anti-bottom quark, respectively.