Excited and exotic bottomonium spectroscopy from lattice QCD

We explore the spectrum of excited and exotic bottomonia using lattice QCD. Highly excited states are identified with masses up to 11,000 MeV, many of which can be grouped into supermultiplets matching those of the quark model while exotic spin-parity-charge-conjugation quantum numbers JPC = 0+−, 1−+, 2+− that cannot be formed from q¯q\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{q}q $$\end{document} alone are also identified. Single-meson operator constructions are used that have good JPC in the continuum, these are found to overlap well onto heavy quark states up to J ≤ 4. A continuum JPC is assigned to each level, based on the distribution amongst lattice irreps and dominant operator overlaps. States with a dominant gluonic component are identified and form a hybrid supermultiplet with JPC = (0, 1, 2)−+, 1−−, approximately 1500 MeV above the ground-state ηb, similar to previous computations with light, strange and charm quark systems.


Introduction
The spectrum of bottomonium mesons has been extensively studied using effective field theories, sum rules, models and ab initio theoretical methods. More recently, the discovery of unexpected, narrow charmonium-like states -the so-called XYZs -has led to a resurgence of interest in heavy quark spectroscopy fuelling theoretical discussions as to their nature and structure. Suggestions include hydrid mesons, molecular or tetraquark configurations or simple quark-antiquark mesons albeit with unexpected masses. In addition, in the plethora of newly discovered states no observation of a state with manifestly exotic quantum numbers, that cannot be produced by quark-antiquark constructions alone, has been made. To date an explanation of the charmonium-like XYZs as well as the absence of charmonium and bottomonium states with exotic J PC has not emerged.
A similar abundance of new states has yet to emerge in the bottomonium sector, where the excited and exotic meson spectrum is less thoroughly explored, although this will be addressed in the near future by LHCb and Belle II [1,2]. A recent review [3] summarises the current theoretical and experimental results in both charm and bottomquark systems. Quark models appear to account for many of the low-lying states, although many are yet to be found. Nevertheless, five hadronic states containing a bb component, which are inconsistent with quark models have already been discovered. In a similar pattern to the charmonium XYZs, three of these states have conventional quantum numbers JHEP02(2021)214 Υ(10580, 10860, 11020) while two states, the (charged) Z b (10610) and Z b (10650) appear to require a four-quark configuration, for example in a compact tetraquark [4,5] or a molecular arrangement [6][7][8]. Calculations of the spectrum of excited and exotic bottomonium mesons are therefore timely.
Relatively little is known about the presence of exotic hadrons in bottomonium from lattice QCD. However, exotic hadrons have been predicted in both integer and half-integer spins, in a series of studies considering light, strange [9][10][11][12][13][14], and charm quarks [15][16][17][18][19]. Evidence has also emerged from exploratory lattice QCD studies of compact tetraquark constructions containing b-quarks that some channels may admit bound-state solutions [20][21][22][23]. The large separation of scale between bottom and light quarks makes reliable firstprinciples lattice QCD computations challenging. Demanding that discretisation effects of O(am Q ) for both light and bottom quarks are small, while simultaneously ensuring that the volume is large enough (m π L > 4) requires very large and very fine lattices.
Effective field theories computed on the lattice, including NRQCD and the Fermilab approach, have played an important role in bridging the gap between maintaining a reasonable computational cost and accessing the physics of interest. In recent work heavy quarks have also been simulated with a relativistic highly-improved staggered quark action (heavy-HISQ) and a range of heavy quark masses are used to enable a reliable extrapolation to the bottom quark mass. An impressive array of precision b-physics quantities has been determined and the methods themselves are now proven to be robust and reliable. It is relatively commonplace to extract ground-state quantities precisely, and several studies have also reported success in computing excited state spectra [24][25][26][27][28][29][30][31][32].
In this study we present results from an exploratory relativistic dynamical anisotropic lattice calculation of the excited spectrum of bottomonium. On an anisotropic lattice the lattice spacing in the temporal direction, a t , is finer than in the spatial direction, a s providing enhanced resolution for spectroscopy. Earlier work using anisotropic lattices for heavy quark physics is described in refs. [33,34]. This work follows refs. [15,19] where charm quarks were simulated using a relativistic Wilson-clover quark action (the same action as used for the light quarks) on an anisotropic lattice with stout-smeared [35] spatial links, and distillation for quark smearing [36]. We apply this relativistic action to bottom quarks for the first time noting that while a s m b > 1, in the temporal direction a t m b remains less than unity. The effect of spatial discretisation is investigated through dispersion relations and is found to be mild, at least within the scope of this calculation. Nevertheless, the mass parameter is of the order of the lattice spacing and so no firm conclusions about the continuum limit can be drawn from this work. A calculation at a second, finer lattice spacing to or similarly comprehensive calculations using other heavy quark methods would be extremely useful to better understand the robustness of the results presented here.
In this work we demonstrate that the techniques developed and used in refs. [9][10][11][12][13][14][15][16][17][18][19] for spin identification from large correlation matrices formed of large bases of operators and analysed variationally, allow for the identification of extensive spectrum of bottomonium mesons, including quark-model-like and hybrid supermultiplets spanning conventional qq and exotic quantum numbers. Some preliminary results from these calculations were presented in ref. [37].

JHEP02(2021)214
Volume a s /fm m π /MeV N cfgs. N tsrcs N vecs.  Table 1. Details of the lattice gauge field ensembles used in this study. The volume is given as (L/a s ) 3 × (T /a t ) where L and T are the spatial and temporal extents of the lattice. The number of gauge field configurations, the number of (perambulator) time-sources per configuration and the number of eigenvectors employed in the distillation method are given as N cfgs. , N tsrcs and N vecs. respectively.
The paper is organised as follows: sections 2 and 3 describe the details and parameters of this lattice calculation. The bottomonium spectrum is presented in section 4, where the determination of exotic states and a supermultiplet of hybrid mesons is also discussed. A summary and outlook is presented in section 5.

Calculation details
We use 2 + 1 flavours of dynamical quarks on an anisotropic lattice with a temporal lattice spacing, a t that is finer than the spatial lattice spacing a s , with anisotropy ξ = a s /a t ≈ 3.5. The light quarks correspond to a heavier-than-physical pion of approximately 391 MeV. A tree-level Symanzik-improved action is used in the gauge sector while the fermion sector is described by a tadpole-improved Sheikholeslami-Wohlert (clover) action, including stoutsmeared spatial gauge fields [35] and distillation quark smearing [36]. Full details of the action and parameters can be found in refs. [38,39]. Table 1 summarises the lattice ensembles used in this study. The final results are produced using two time-sources on an (L/a s ) 3 × T /a t = 20 3 × 128 volume at a single value of the lattice spacings (a s , a t ). To test the volume dependence, a single time-source on a 24 3 volume is used. Dispersion relations are extracted from energy levels at finite momentum computed with a single time-source on the 20 3 volume.

Lattice heavy quarks
In an anisotropic lattice computation with dynamical quarks, the bare gauge and fermion anisotropies, which enter in the Monte Carlo ensemble generation, must be simultaneously tuned to produce a target value when measured non-perturbatively [38]. For heavy valence quarks, an additional tuning of the fermion anisotropy is useful to ensure a consistent measured anisotropy [15]. The tuning condition is the same as for light quarks namely that the speed of light, c = 1, or in other words that a relativistic dispersion relation is measured. For charm quarks, it was demonstrated that a single valence anisotropy parameter, tuned with the pseudoscalar (η c ) dispersion relation, also yields a relativistic dispersion relation for the vector (J/ψ) and for the heavy-light (D-meson) sector, consistent within uncertainties.
The dispersion relation also plays an important role in the Fermilab treatment of lattice heavy quarks [40]. In this approach a mass-dependent tuning of parameters guarantees JHEP02(2021)214 a smooth interpolation between the light and heavy sectors. The bottom quark mass in the simulation is tuned so that the kinetic mass measured from the quarkonium dispersion relation takes its physical value and need not be the same as the rest mass. In a relativistic simulation the rest and kinetic masses are equal but this is not required in the Fermilab approach, while in non-relativistic formulations of QCD (NRQCD) the rest mass is discarded.
In this work the bottom quark is treated relativistically, computed with the same tadpole-improved Wilson-clover action used for the light quarks. The bottom quark mass parameter in the simulation is determined by imposing that the pseudoscalar meson mass takes the physical η b value [41]. We investigate the occurrence and significance of discretisation errors using dispersion relations, and mass splittings of the pseudoscalar and vector ground states. The relativistic form of the dispersion relation for a meson A can be written and since the quark fields in the meson satisfy periodic boundary conditions for the finite extent of the lattice, the momenta are quantised as a s p = 2π L (n x , n y , n z ) with n i ∈ Z and for simplicity we label momenta by [n x n y n z ] in the text. The renormalised anisotropy, ξ A can be determined from the slope of the dispersion relation. The input anisotropy parameter in the valence heavy-quark action was tuned so that the pseudoscalar (η b ) dispersion relation produces ξ ≈ 3.5. Keeping the heavy-quark mass and the input anisotropy fixed, the anisotropy determined from the heavy-heavy vector (Υ) and in heavy-light (B-meson) systems was then compared with the η b result.
The dispersion relations for the η b and Υ mesons are shown in the left panel of figure 1, for a range of momenta up to p = [211]. 1 At each momenta a large basis of interpolating operators has been used in a variational analysis to determine the energy level in the appropriate lattice irreducible representation, as described in section 3. The rest mass and anisotropy determined from fits to eq. 2.1 are summarised in table 2. The goodness-of-fit (χ 2 /N dof ) values obtained show the data are well-described by eq. 2.1. Little evidence of deviation from the relativistic form of the dispersion relation is apparent below momentum [210], while [210] and [211] appear to show a small differences. We note that for these latter high momenta points the correlation functions only plateau at large euclidean times and the energy levels are less robustly determined. The measured anisotropies are ξ ηc = 3.590 (15) and ξ Υ = 3.574 (26), within 3% and 2% of the target value respectively.
To probe the reliability of the simulations further we determine the pseudoscalar and vector dispersion relations in the B-meson sector using the same input parameters as for the η b and Υ. The measured anisotropy is determined from the slope of the dispersion relation and comparing this parameter in the heavy-heavy and heavy-light sectors effectively isolates momentum-dependent discretisation effects providing a clean test of the robustness of heavy-quark simulations. This was highlighted in an analysis of differences in binding energies between heavy-heavy and heavy-light systems in refs. [43,44] Table 2. Parameters determined from dispersion relation fits, using eq. 2.1, to rest and movingframe energies in J P = 0 − , 1 − for both heavy-heavy and heavy-light systems, with n as defined in the text below eq. 2.1. The energies used and the fit for n 2 ≤ 4 are shown in figure 1.
anisotropies of the pseudoscalar and vector heavy-light mesons are ξ B = 3.360 (29) and ξ B * = 3.293(31), within 4% and 6% of the target value respectively. As noted, there is some discrepancy between the anisotropy measured in the bottomonium and B meson systems but in each case the measured values are at most 6% from the target anisotropy. The results presented here are not sensitive to the discrepancy between heavy-light anisotropies and the target value, which becomes relevant only in scattering calculations. A similar, although smaller, effect was observed between charmonium and D mesons. Improved determinations of the energy levels at finite momentum and tuning on a larger volume may JHEP02(2021)214  Table 3. A comparison of the 1S hyperfine splitting determined with the tree-level clover coefficient and with increased values, c s = 2.0 and 2.2 to investigate the possible effects of nonperturbative improvement and to provide an approximate scale for the discretisation error reduce this discrepancy, which will be investigated in future work. In conclusion, the dispersion relations for both heavy-heavy and heavy-light vector and pseudoscalar mesons yield precise determinations of the measured anisotropy demonstrating relativistic behaviour including for large values of spatial momentum. Further calculations at multiple lattice spacings would allow for a fuller accounting of the sources of discretisation uncertainties.

Spatial discretisation effects and the hyperfine splitting
The hyperfine splitting inferred from the masses in table 2 2 is significantly smaller than the experimental value. A similar effect was noted using this approach in charmonia, as reported in ref. [15]. This underestimate is a common feature of calculations at finite lattice spacing. Nevertheless, several other approaches have determined the hyperfine splitting using lattice QCD [24-26, 28, 30, 32], in many cases finding good agreement when all systematic uncertainties are accounted for. The (spatial) clover term c s appears in eq. 5 of ref. [38], and its tree-level tadpoleimproved value is used consistently for the dynamical quarks and the valence bottom quark computed here. However, this is expected to underestimate the non-perturbative value [45] required for full O(a) improvement. The one-loop perturbative value of c s is also known to significantly underestimate the non-perturbative value, in particular for coarser lattice spacings, such as the spatial lattice spacing used in this study [46]. The effect of increasing c s was explored for charmonia in ref. [15] where it was found that by replacing the tree level value c s = 1.35 with c s = 2 the 1S hyperfine splitting increased from ∆ cs=1. 35 = 80(2) MeV to ∆ cs=2.0 = 114(2) MeV, significantly closer to the physical value of 117 MeV.
A similar analysis is repeated here for bottomonium and summarised in table 3. A clover coefficient c s = 2.2 was chosen to approximate the nonperturbative value at the bottom quark mass, and the η b -Υ mass difference was determined. Empirically, the hyperfine splitting increases with increasing c s , ∆ cs=1.2 = 19.5±0.6 MeV to ∆ cs=2.2 = 55.3±0.4 MeV, much closer to the experimental value of ∼ 61 MeV [41] Disconnected contributions are not included in this analysis, as was also the case in charmonium, however these are expected to be small due to OZI suppression and are unlikely to play a significant role in the discrepancies with experimental values discussed above.
We conclude that for the anisotropic action with stout smearing used here, the O(a s p) discretisation errors are under reasonable control for bottomonium and heavy-light simulations at zero and finite momenta up to [200]. The latter is important as in subsequent work we will extend our analysis to include meson-meson operators for states above de-JHEP02(2021)214 cay threshold. Finally, the analysis of the effect of enlarging the clover term suggests a scale for the discretisation errors on the masses determined in this study of approximately 35-40 MeV.

Operator construction, fitting and spin identification
Following the notation and methods developed and described in refs. [9,10,15], meson energies are determined from two-point correlation functions where the operators O † j (0) create the state of interest at t = 0 which is later annhilated by O i (t) at euclidean time t. The correlation functions have a spectral decomposition where Z n i = n|O † i |0 are the operator overlaps. Distillation is used to create correlation functions, facilitating the efficient construction of large operator bases. A derivative-based construction is used, combining gauge-covariant forward-backward, spatial derivatives and gamma matrices in a fermion bilinear to form operators of definite J PC . The discrete lattice breaks rotational symmetry, and the relevant symmetry group for a lattice calculation is O h whose five, ten including parity, irreducible representations (irreps) are labelled The distribution of a continuum spin state across the lattice irreps is shown in table 4, for J ≤ 4. The assignment of continuum spin values to lattice energy levels, characterised by lattice irreps, is not straightforward, the procedure we follow is discussed in refs. [9,10]. Finally, the lattice operators used in this work are constructed to transform in a definite lattice irrep and row, Λ and λ respectively, and are derived from continuum operators O J,M by subduction.
where M is the J z component of spin. In each irrep we include operators in the basis that are proportional to the commutator of two covariant derivatives, the field-strength tensor, to investigate the pattern of hybrid mesons. The number of operators in each lattice irrep, Λ PC , used in this study is given in table 4. The energy levels and operator overlaps are extracted from two-point correlation functions using the variational method [47,48]. In each lattice irrep energies and overlaps are extracted from a matrix of two-point correlators, C ij (t) by solving a generalised eigenvalue problem where i and j label the operators and t 0 is the reference timeslice. The generalised eigenvalues, or principal correlators, λ n yield the energies from fits to  with free parameters m n , m n , A n , where m n is the energy of the state of interest. An example of the principal correlators determined in this study is shown in figure 2. In general we find this method continues to work well with bottom quarks, as it also did in studies of the excited charmonium and open charm mesons and for light (isovector and isoscalar) mesons and baryons.

Spin identification
Lattice energy levels are assigned continuum spin values following the approach described in refs. [9,10]. The lattice operators used in this work are subduced from continuum operators with definite spin, as detailed in eq. 3.3, and empirically it has been observed and E −− . The masses and overlaps are found to be quite consistent across lattice irreps. We thus identify these four states as components of the same continuum 4 −− level. that they do not depend significantly on the specific lattice irrep into which they are subduced. Continuum spin assignments are then made by comparing the overlap values determined independently in different lattice irreps which have been subduced from the same continuum operator. This is a powerful tool to discriminate the spin of higher-lying energy levels especially in a dense spectrum with many states that are close to degenerate in mass. This method was shown to work well in the charmonium spectrum [15] and the same methodology is employed here.
An illustrative example is given in figure 3, which shows the overlap values Z n i determined in the variational analysis for the (a 1 × D and E −− when subduced and we find the operator produces a consistent value across irreps. In general we find good agreement between Z values that originate from the same continuum operator, and we use this approach to assign continuum spin values up to J = 4. An illustrative example of the relative operator overlaps for each state n is given in figure 4, here the overlap for operator i is plotted usingZ n i = Z n i /max ∀m (|Z m i |). Even in this densely packed T −− 1 irrep, the spin assignment is abundantly clear.

Results
In this section we present results for the spectrum of bottomonium mesons including candidate hybrid states. The results, organised by the relevant lattice irrep are shown first, followed by the final spin-identified spectrum labelled by the continuum quantum numbers J PC .

The bottomonium spectrum by lattice irrep and volume comparison
The spectrum of excited and exotic bottomonium states determined on the 20 3 volume and arranged by lattice irreps, labelled Λ PC , is shown in figure 5  relative to the energy of the lowest state in the system, the η b meson, to mitigate the uncertainty due to tuning the heavy quark mass. The vertical height of the boxes shows the one sigma statistical uncertainty about the mean from a determination of energy levels using the variational analysis, as described earlier. The energy levels are colour-coded according to the spin, identified as described above (black for J = 0, red for J = 1, green for J = 2, blue for J = 3 and gold for J = 4), as in ref. [15].
In figure 6 the volume dependence of the extracted spectrum is investigated, by comparing results from two spatial volumes, 20 3  observed although with some dependence on volume visible particularly in the higher-lying states, albeit with larger statistically uncertainties on the 24 3 volume.
In this study single-meson operators constructed from fermion bilinears are utilised. Since no multi-hadron operators are included a strong overlap onto multi-hadron states is not expected. Investigations of the consequences of this choice have found that when using only approximately-local operators a level usually arises within the width of the state [49,50]. We find that the calculated spectrum is well-described by single mesons as indicated by the only mild volume dependence and very consistent agreement for a state of definite continuum spin subduced into its relevant lattice irreps. However, we cannot rule out the presence of additional spectroscopic states that may be found by extracting and analysing the scattering amplitude [51][52][53]. In a finite volume calculation we expect there to be eigenstates associated with each channel allowed by the quantum numbers of each irrep. Some relevant combinations of two-and three-hadron strong decay thresholds include η b η, η b ππ and BB, the latter shown in figures 5, 6, 7, 9 and 10. Decay channels involving B-hadrons contain the leading quark-line connected diagrams and are expected to produce the most significant effects.

The spin-identified bottomonium spectrum
The spin-identified spectrum of bottomonium mesons labelled by continuum quantum numbers, J PC , is shown in figure 7. The numerical values relative to the η b , in MeV units, are provided in table 5. Using the operator overlaps to identify the different components of the given J = 2, 3, 4 states split across lattice irreps, we do not find statistically significant discretisation effects due to rotational symmetry breaking, and therefore continuum spins of states up to J = 4 are assigned. Clusters of states appear to follow quark model supermultiplets. For example, in the negative parity sector (leftmost on the plot) (1, 2, 3)S and (1, 2)D are visible with possibly a few states of the 1G supermultiplet. We find no 0 −− state in the energy region we consider. In the positive parity sector (middle of plot) (1, 2)P and 1F are clearly identified. The left panel of figure 8 shows an example of the 1F operator overlaps for (π, ρ) × D [3] J 13=2 , J=3 J , which are seen to be consistent across supermultiplet members. This enables us to identify the third state in 2 ++ as 1F , while the nearby second state forms part of the 2P supermultiplet.
The exotic quantum numbers, 1 −+ , 0 +− and 2 +− each contain a level, and there are additional levels not accommodated in the quark model supermultiplet counting which are determined in several non-exotic irreps.

Hybrid supermultiplets
Candidate hybrid mesons are identified by examining the operator-state overlaps for each energy level in a given lattice irrep. States with hybrid content are highlighted in red and blue in figure 7, while all others appear in green.

JHEP02(2021)214
A state is proposed to be a hybrid if it is characterised by a relatively large overlap onto an operator proportional to the field-strength tensor. In practice, this means that for a given state the dominant overlap factor (Z) is from operators of the form ((π, ρ) × D [2] J=1 ). A typical example can be seen in figure 4 where the 8th state in the T −− 1 channel is clearly dominated by these hybrid-like operators.
The same procedure was followed in previous analyses of mesons in the light, strange, open-charm and charmonium meson sectors. In each case a supermultiplet of hybrid states in [(0, 1, 2) −+ , 1 −− ] was identified at an energy scale ≈ 1.3 GeV above the lightest state in the spectrum. It was previously commented that such a pattern of states is predicted by bag models and models with constituent gluonic degrees of freedom [54][55][56], but appears to rule out some flux-tube models [57].
A similar pattern of hybrids is identified in this bottomonium system. The T −+ 1 irrep contains J PC = 1 −+ , 3 −+ and higher J, some of which are not exotic based on qq considerations. One difference to studies with lighter quarks is that a state dominated by a J PC = 4 −+ operator arises at a lower energy than the first level with J PC = 1 −+ , as can be seen in figures 5 and 6. The second state in T −+ 1 is dominated by operators containing a field-strength tensor in their constructions, singling it out as a hybrid candidate. Furthermore, all of the irreps with J PC = (0, 1, 2) −+ , 1 −− have a state dominated by such operators around a t m = 1.93. In T −− 1 this is state n = 8, whose relative overlap was shown earlier in figure 4.
In the right pane of figure 8, we show the overlap of the simplest hybrid operator construction, 4 (π, ρ) × D [2] J=1 J and find that, as for the 1F supermultiplet, we can identify a hybrid supermultiplet across lattice irreps that has a common origin, with masses in the range a t m = 1.93 − 1.94. These states are highlighted in red in figure 7. Other states identified as hybrids are highlighted in blue and may form parts of a higher-lying exotic supermultiplet. In figures 5 and 6, the hybrid candidates have been highlighted with a dotted outline. This qualitative picture is very similar to that found with charm quarks in refs. [15] and [19]. Figure 9 shows the similarity in the pattern and the mass scale of the lightest hybrid supermultiplet identified in charmonium [15] and in bottomonium as described here. Finally, we note that in refs. [58,59] a similar pattern of bottomonium hybrids is inferred using the charmonium lattice calculations from ref. [15] as input, although many more states are predicted there than are found in this calculation.

Comparison to experiment
The aim of this study has been to investigate the qualitative features accessible in bottomonium using a relativistic action and a basis of quark-antiquark operators in a proven methodology for extracting highly excited states. Although statistically precise, we do not necessarily expect these results to be quantitatively accurate for several reasons, including discretisation effects due to the heavy quark mass, the presence of additional levels due to hadronic decays, the larger-than-physical light quark mass and the lack of sea quark effects JHEP02(2021)214  for bottom and charm. While a comprehensive error budget is beyond the scope of this work it is still informative to make a comparison with the existing experimental spectrum.
In figure 10, we show the experimental energies listed in ref. [41] compared with those determined in this work. The hyperfine splitting is smaller in the lattice spectra, compared with experiment, as discussed above. The 2S states lie slightly below their experimental counterparts. The 1P states are in reasonable agreement with experiment while the 2P multiplet is consistently higher in energy. Comparing spin-averaged splittings, we find that ∆ 2S−1S = 564(3) MeV is in good agreement with the corresponding experimental value of 563.3 MeV. The spin-averaged (1P-1S) splitting is 495(3) MeV in this study compared to 455.12 MeV determined from experimental values [41].

JHEP02(2021)214
Experimentally, the 1P and 2P multiplets are observed to satisfy a simple nonrelativistic quark-model spin-splitting relation, 1 9 m 3 P 0 + 3m 3 P 1 + 5m 3 P 2 − m 1 P 1 ≈ 0, and these even appear to be relatively insensitive to the effects of nearby decay channels [60]. This splitting in the spectrum determined here is −6±5 MeV, which approximately satisfies this relation.
An extra level belonging to the 1D supermultiplet appears in 1 −− , similar to those seen in quark potential models [61]. This state has not been unambiguously observed experimentally although the 2 −− member is known [62,63]. The comparable state in charmonium, the ψ(3770), is clearly seen in both hadronic and e + e − processes, although the lighter quark mass and the proximity of DD are key differences that may mix the S and D-wave eigenstates and enlarge production rates.
Quark potential models typically predict a small electromagnetic coupling to the Υ(1D) [64] and it can in principle be computed on the lattice from a local vector current of the form 0|ψγ i ψ|Υ(1D) , for example as in ref. [65]. While such a calculation is beyond the scope of the present work, we already have the equivalent smeared vector operator, 5 visible as the fourth |Z| plotted from the eigenvectors of each state in figure 4. The Υ(1D) is identified as state n = 2 where this vector operator makes a tiny contribution, unlike the majority of states in the spectrum identified as having J PC = 1 −− . The pattern of overlaps suggests that state n = 5 is the Υ(2D), this receives a slightly larger contribution from this smeared vector operator. We note that the state n = 8, identified as a hybrid, also has a relatively small overlap with this smeared vector operator.
The mixing of the S-D eigenstates in 1 −− can be qualitatively assessed, once again from the operator overlaps. The final two non-hybrid vector operators shown in red in figure 4 ((ρ × D [2] J=2 ) J=1 and (ρ 2 × D [2] J=2 ) J=1 ) have previously been found to have dominant overlaps onto 2S+1 L J = 3 D 1 states [13]. The same is true here with a tiny overlap in states n = 0, 1 and a dominant overlap in state n = 2. This, and the previous comments regarding the vector operator, suggests very weak S-D mixing in the lowest three 1 −− states. 6 The same cannot be said of states n = 4, 5, 7 where significant contributions from both operators are present. These states are found on either side of the BB threshold, not dissimilar to the ψ(3770) and ψ(2S) in charmonium which is found above and below DD.
The operator basis used in this T −− 1 irrep contains 26 operators, many more than the 10 states we quote. Further levels are found at higher energies, but none lower in mass. We stress that this does not necessarily mean that further levels are not present, we know there must be hadron-hadron levels which are absent, and there is no guarantee the operators used will overlap well with all the physical eigenstates. Of the 26 operators used, 5 pairs are identical in the non-relativistic limit, which is a good approximation in bottomonium. The absence of any states belonging to the 4S and 5S supermultiplets is striking, but this should not be taken as evidence that they are not present, it is possible such states could be found in a more sophisticated calculation, even below the hybrids in 0 −+ and 1 −− . Figure 10. A selection of energy levels determined in this study for which there are corresponding experimental values. The energy levels obtained in this study, shown in green, compared with experiment [41], shown in black. The green (grey) dashed lines show relevant low-lying strong decay thresholds determined from this lattice data (experiment). Similar differences can be found in 2 ++ . Qualitatively, three states are found experimentally and three are present in this calculation, however, as discussed above and shown in the left panel of figure 8, the third level appears to be part of the 1F supermultiplet while the experimental candidate is assigned to 3P . No clear 3P states arise in our calculation, but again, this absence should not be taken as clear evidence that they are not present in reality.

Summary
A first study of bottomonium including exotic and hybrid mesons from lattice QCD has been presented. An extensive pattern of states comparable to that predicted by quark models is determined. In addition, several hybrid meson candidates are identified at an energy approximately 1500 MeV above the ground state η b , some of which can be grouped into a hybrid supermultiplet with J PC = (0, 1, 2) −− , 1 −+ . Similar qualitative features have previously been observed in work by our collaboration in studies spanning light, strange and charm quarks and now for the first time in bottom-quark physics.
This study was carried out with 2 + 1 dynamical quark flavours using distillation for quark propagation. A single lattice spacing and a b-quark mass parameter slightly less than unity are used so that discretisation uncertainties remain unquantified. Within the scope of the calculation we have assessed their impact by considering dispersion relations of heavy-

JHEP02(2021)214
heavy and heavy-light pseudoscalar and vector mesons and by varying the coefficient of the spatial clover term c s in the heavy-quark action. Rotation breaking effects are also found to be small and continuum spins are reliably assigned for states up to J = 4. A large basis of carefully-constructed single-meson operators is used however hadron-hadron or multihadron operators are not included. These are expected to be relevant above the indicated hadron-hadron thresholds where additional eigenstates should be present that would mix with states extracted here. Several studies have demonstrated that the methodology used here produces an eigenstate within the hadronic width of the resonance in the vast majority of cases. Such resonance effects are highly volume dependent and we have tested that the same qualitative picture holds on a larger volume where in the majority of cases good quantitative agreement is found.
A similar determination of the B meson spectrum and the B c spectrum is underway, and some preliminary results were presented in ref. [37]. These methods appear to be sufficiently reliable that investigations of near-threshold states can now begin, presenting many interesting possibilities.

Acknowledgments
We thank our colleagues in the Hadron Spectrum Collaboration, in particular Jozef Dudek and Christopher Thomas for useful comments. We also thank Tim Burns for useful comments. DJW acknowledges support from a Royal Society University Research Fellowship and from the U.K. Science and Technology Facilities Council (STFC) [grant number ST/P000681/1. SR acknowledges the hospitality of the HEP group at DAMTP where this project was started and the support and hospitality of the Technical University of Munich -Institute for Advanced Study through the award of an August Wilhelm Scheer Visiting Professorship. SR is supported by the European Union's Horizon 2020 research and innovation programme under grant agreement No. 824093. This work was performed using resources provided by the Cambridge Service for Data Driven Discovery (CSD3) operated by the University of Cambridge Research Computing Service (www.csd3.cam.ac.uk), provided by Dell EMC and Intel using Tier-2 funding from the Engineering and Physical Sciences Research Council (capital grant EP/P020259/1), and DiRAC funding from the Science and Technology Facilities Council (www.dirac.ac.uk). The DiRAC component of CSD3 was funded by BEIS capital funding via STFC capital grants ST/P002307/1 and ST/R002452/1 and STFC operations grant ST/R00689X/1. DiRAC is part of the National e-Infrastructure. Computations were also performed at Jefferson Laboratory under the USQCD Initiative and the LQCD ARRA project and on the Seagull cluster maintained by the Trinity Centre for High Performance Computing (TCHPC).
The software codes Chroma [66], QUDA [67,68], QPhiX [69], and QOPQDP [70,71] were used to compute the propagators required for this project. This research was supported in part under an ALCC award, and used resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725. This research is also part of the Blue Waters sustained-petascale computing project, which is supported by the National Science Foundation (awards OCI-0725070 and ACI-1238993) and the state of   (4) 1099 (17) 1042(21) 3 +− 1137 (12) 3 ++ 1133 (12) 4 ++ 1141 (15) 1 −+ 1536 (16)  Gauge configurations were generated using resources awarded from the U.S. Department of Energy INCITE program at the Oak Ridge Leadership Computing Facility, the NERSC, the NSF Teragrid at the TACC and the Pittsburgh Supercomputer Center, as well as at Jefferson Lab. In table 5 we summarise the masses found in physical units relative to the η b as measured in this study, as presented above in figure 7.

A Masses summary table
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.