K-matrix Analysis of $e^+e^-$ Annihilation in the Bottomonium Region

We perform the first global and unitary analysis of $e^+e^-\to b\bar{b}$ cross sections. We analyze exclusive cross sections in the $B\bar{B}$, $B^* \bar{B}(+c.c.)$, $B^*\bar{B}^*$, $B_s^*\bar{B}_s^*$, $\Upsilon(nS)\pi^+\pi^-$ and $h_b(nP)\pi^+\pi^-$ channels as well as the total inclusive cross section for $b\bar{b}$ production. Pole positions and residues are determined for four vector states, which we associate with the $\Upsilon(4S)$, $\Upsilon(10750)$, $\Upsilon(5S)$ (or $\Upsilon(10860)$), and $\Upsilon(6S)$ (or $\Upsilon(11020)$). We find strong evidence for the new $\Upsilon(10750)$ recently claimed by Belle, although with parameters not well constrained by the data. Results presented here cast doubt on the validity of branching ratios reported earlier using Breit-Wigner parameterizations or ratios of cross sections. We also compare our results with a selection of theoretical calculations for the vector bottomonium spectrum.


I. INTRODUCTION
The spectrum of vector (J P C = 1 −− ) bottomonium states above BB threshold has been the source of a series of surprises and unresolved issues. The initial exploration of this region using inclusive e + e − annihilation to hadrons [1,2] showed evidence for the production of two states with masses heavier than the Υ(4S), consistent with potential model expectations for the Υ(5S) and Υ(6S). More recent measurements of the same process [3][4][5] have revealed more complex structure. While the putative Υ(5S) and Υ(6S) states (also called the Υ(10860) and Υ(11020), respectively) still appear as prominent peaks in the inclusive cross section, the effects due to coupled-channel scattering and the opening of a variety of open bottom thresholds (e.g., B * B , B * B * , B sBs , B * sBs , B * sB * s ) are now more apparent and complicate the observed spectrum. Extracting vector bottomonium masses, total widths, and partial e + e − widths from these spectra has posed serious challenges. While fits to the inclusive e + e − spectrum using a coherent sum of Breit-Wigner amplitudes are possible [5], the fits violate unitarity and the results are expected to be unreliable.
Recent measurements of the energy dependence of exclusive e + e − → B ( * )B( * ) cross sections [6] confirm the importance of coupled-channel scattering. Rather than showing distinct peaks for the Υ(5S) and Υ(6S), the cross sections are marked by dramatic peaks and valleys at various open-bottom thresholds. These non-trivial features in the open-bottom cross sections undermine older measurements of the Υ(5S) and Υ(6S) branching fractions, such as those currently listed in the Particle Data Group's Review of Particle Properties (RPP) [7]. Previous branching fractions of the Υ(5S) to open-bottom final states, for example, were estimated by first measuring the cross section of e + e − to a given open-bottom final state at an energy near the presumed mass of the Υ(5S) and then dividing by the inclusive bb cross section at the same energy [8][9][10]. This ratio of cross sections would approximate an Υ(5S) branching fraction only if the Υ(5S) were produced in isolation, an assumption we now know to be false.
Besides strong coupled-channel effects in the openbottom final states, anomalously large cross sections for e + e − to closed-bottom channels, such as ππΥ(nS) (n = 1, 2, 3) and ππh b (nP ) (n = 1, 2), have been observed [11][12][13][14][15]. Their production rates were later found to be enhanced by the presence of the exotic isovector bottomonium-like states, the Z b (10610) and Z b (10650) (also called the Z b and Z b ), which decay to πΥ(nS) and πh b (nP ) [16,17]. In contrast to the complications in the open-bottom channels, the Υ(5S) and Υ(6S) appear to be well-isolated in the closed-bottom channels, which allows for a more reliable extraction of their mass and width [15]. These relatively well-behaved cross sections also provide evidence for an additional state, the Υ(10750) [15], which may be the Υ(3D) bottomonium state.
With the recent publication of inclusive e + e − cross sections [5] and exclusive e + e − cross sections to openbottom [6] and closed-bottom [15] final states, we are now in a position to perform the first global and unitary analysis of the vector bottomonium system above BB threshold. We use the K-matrix formalism for this analysis. In the literature, the K-matrix is regularly used in the spectroscopy of hadrons containing u, d and s quarks, for example in studies of scalar mesons in Refs. [18][19][20]. A first application of a K-matrix in the analysis of heavy quarkonia was made by Uglov et al. [21], performing a unitary coupled channel analysis of e + e − annihilation to the DD, D * D and DDπ channels for energies up to 4.7 GeV. Masses and partial decay widths for the charmonium resonances ψ(3770), ψ(4040), ψ(4160) and ψ(4415) were extracted. Comparing to Ref. [21], our study of the bottomonium system includes three-body final states, allows for non-resonant scattering, is further constrained by the total inclusive cross section and allows for an analytic continuation to the complex energy plane.
In the following, we describe our model (Sec. II), the datasets used in the analysis (Sec. III), and the fit procedure (Sec. IV). In Sec. V, we present our results for the pole positions of the Υ(4S), Υ(10750), Υ(5S), and Υ(6S), as well as estimates for their partial widths to all considered channels. We discuss the results in Sec. VI, comparing them to a variety of theoretical calculations, and we share our conclusions in Sec. VII. Finally, a more detailed description of the K-matrix formalism is given in the Appendix (App. A) together with a discussion of an analysis that omits three-body channels (App. B).

II. K-MATRIX FORMALISM
Our fit to the data will employ a K-matrix that includes resonant and non-resonant scattering terms and is modeled as follows: The index R refers to resonances while Greek indices denote continuum channels. Except where noted below, the resonant couplings (g) and non-resonant terms (f ) incorporate energy-dependent form factors (see Eqs. A31 and A32) that are meant to capture the hadronic nature of the relevant interactions. Form factors are parameterized in terms of a universal scale β, while the resonant and non-resonant coupling strengths are denotedĝ and f , respectively. Derivations and further details can be found in Appendix A.
The scattering amplitude is written as where C is the Chew-Mandelstam function, with the property and ρ is the diagonal phase space matrix with elements Here, S ν is a symmetry factor and k ν is the center of mass momentum in channel ν. A detailed description of the model choices and our reasons for them is given in Appendix A. Some of the fits described below employ Aitchison's P -vector formalism for production [22]. In our case this is implemented as where P ν = K ν,ee , andK andĈ are defined in the restricted channel space that does not include the initial channel. No form factors were used in the e + e − channels of the production vector, as these are not relevant for leptons. Additional modeling is required since we are fitting three-body channels in a two-body formalism. We have chosen two approaches: (i) Create mock two-body channels for (bb)ππ consisting of the bottomonium state and a "ππ" state of mass 2m π with relative angular momentum set to zero. This is a common method used in hadronic scattering analyses.
(ii) We note that all the three body channels have cross sections of order 10 pb, while the open bottom two-body channels are order 100 pb. Thus it is reasonable to treat the three-body channels perturbatively, neglecting terms of order g 2 R:3body . This can be realized by placing the production (e + e − ) channel in the K-matrix and the threebody channels into a "final state" matrix, which we call F . Thus (see Appendix A for details) where ∆ denotes a three-body channel,K is defined in the restricted channel space that includes two-body channels and the production channel, and We are at liberty to choose an appropriate model for the F vector. Two different model choices will be tested: (iia) Resonances R feed the three-body channels ∆, yielding (Notice that the units of g R:µ are dimension one while g R:∆ are dimension zero.) (iib) The two-body channel µ feeds the three-body channel ∆ directly, As a further extension, which we did not pursue, intermediate states saturating ππ (say an f 0 ) in a three-body channel like ππΥ could be included if efficiency corrected Dalitz plots were available. Form factors were not incorporated in the modeling of three body channels since this was judged to be an unnecessary complication.
Once the scattering amplitude is computed, cross sections are obtained by integrating over the relevant Dalitz plot region. This method comprises a rigorous, if perturbative, approach to solving the three-body problem in the K-matrix formalism. The general problem requires solving integral equations and has a long history of investigation [23]. Results for both of our methods (iia) and (iib) will be reported below.
An additional novel aspect of our approach is the use of the inclusive scattering process. This is achieved by fitting the sum of all channels to σ bb = R bb · σ(e + e − → µ + µ − ). An immediate problem that this procedure raises is that the sum over the measured channels falls short of the inclusive cross section for high center-of-mass energies √ s 10.75 GeV, implying that missing channels can be important. Possible missing open bottom channels that are relevant to the energy region of the fit are listed in Table I. Because this represents a substantial source of modeling ignorance, we choose to represent this dynamics with a single dummy channel chosen to correspond to the lightest channel. Thus the channel masses are set to m(B s ) and the relative angular momentum is = 1. We have tested sensitivity to these choices by adjusting the dummy threshold to m(B) + m(B * ) + m(π) and 2 · m(B * ) + m(π) and find negligible difference. Additional error can arise in the model due to neglected resonances with masses much greater than √ s in the fit region. If considered, these would contribute to the K-matrix. This is equivalent to a non-resonant interaction, and is therefore accounted for in the fit to a large extent. Lastly, we consider neglected high mass channels. In this case the relative channel momentum will be zero and there will be no effect in the scattering amplitude for P-waves and beyond. For S-waves the imaginary part of the amplitude denominator will not change. The real part will shift byĝ 2 [1/(8π 2 )−s/(96π 2 M 2 )+. . .] and thus will be largely absorbed in the bare parameters, and will not affect the physical properties of the amplitude.

III. DATA SETS
With Belle, BaBar, CLEO, and CUSB, four different experiments made contributions to the study of vectorbottomonium states in electron-positron annihilation in the past, with the latter three focused on the inclusive cross section, while the former also studied exclusive cross sections at various center-of-mass energies above the open-bottom threshold.
For this work, we use measurements of the ratio of inclusive bottom anti-bottom quark production over muon production by BaBar and Belle [3,4]. This ratio R bb is obtained from a measurement of the total hadronic cross section, subtracting u, d, s, c-quark contributions extrapolating a model fitted to experimental data of σ(e + e − → hadrons) below the open-bottom region (see Ref. [24] for details). A detailed comparison between the BaBar and Belle measurements and a determination of radiative corrections has been performed in Ref. [5]. We will use the combined data reported in Ref. [5] for the dressed R bb ratio and multiply by σ(e + e − → µ + µ − ) = 4πα 2 3s to obtain the total dressed cross section for inclusive bottom-quark pair production σ bb .
In addition, the Belle collaboration recently published cross section measurements of many exclusive open-(e + e − → (bq)(bq)) and hidden-bottom (e + e − → (bb)ππ) processes. The open-bottom processes include BB, B * B and B * B * [6] and B * sB * s [25]. Here and in the following, B s , only "visible" (directly measured) cross sections were made available. Initial state radiation effects are in general reaction-dependent and can have an influence on the reconstruction efficiency, so that there is no straightforward way to convert the visible cross section for e + e − → B * sB * s into a dressed cross section. Instead, we use the visible cross section in our fit for this process alone and thus implicitly assume that the ISR correction is small compared to the sizable uncertainties of the measurement. The hidden-bottom channels measured by Belle were Υ(nS)π + π − (n = 1, 2, 3) [15] and h b (nP )π + π − (n = 1, 2) [14]. We also consider the measurements of Υ(1S)π + π − and Υ(2S)π + π − at the Υ(4S) resonance [15,26]. For these hidden-bottom channels, we abbreviate π + π − as ππ in the following. In case of all hidden-bottom processes, Born cross sections-observed cross sections corrected for initial-state radiation (ISR) and vacuum polarization (VP) effects-were published. As VP effects are independent of the final state, we multiply the hidden-bottom cross sections by the VP correction factors given in Ref. [5] to turn the Born cross sections into dressed cross sections. If no value of the VP correction factor is available at a given center-of-mass energy, we interpolate linearly between two neighbouring points. Given the large density of data for R bb and thus the VP correction factor in Ref. [5], we only have to interpolate over small center-of-mass energy regions. Ref. [15] shows additional data points in the Υ(10753) region obtained using the ISR method (e + e − → γ ISR Υ(nS)ππ). These are not publicly available because the effective luminosity changes rapidly with center-of-mass energy for these measurements [27]. These data points are thus ignored here.
In principle, as was mentioned in Sec. II (see Eq. 10), our model could accommodate intermediate resonances in three-body processes. However, no information is available on the center-of-mass energy dependence of the Υ(nS)ππ and h b (nP )ππ Dalitz plots. These would be of particular interest for future high statistics measurements by Belle II.
We define a total of nine different models, using three different form factor scales β = 0.8, 1.0, and 1.2 GeV and three different treatments of the three-body hiddenbottomonium channels. In the first set of models, threebody channels are described in a quasi two-body approach including a (ππ) quasi-particle with mass m ππ = 2m π . In the other two sets of models, three-body channels are treated perturbatively, using either a resonant or non-resonant coupling as defined in Eqs. 8 and 9. Given the lack of data on the center-of-mass energy dependence of the Υ(nS)ππ and h b (nP )ππ Dalitz plots, we do not use F (∆) µ given in Eq. 10. In addition, we also consider three two-body models applied to the open-bottom channels only as a test of the importance of the three-body data and the robustness of the conclusions (see Appendix B).
The difficulty in estimating initial parameters, the large number of parameters and, by construction, strong correlation between them, forced the adoption of a stepby-step fit procedure. We start using a single pole and a single channel (BB) to describe the total cross section σ bb in the Υ(4S) region below the B * B threshold, where σ bb ≈ σ(e + e − → BB) and Br(Υ(4S) → BB) ≈ 100%. Then, we add the remaining three poles together with the BB data. In an iterative procedure, we then add a new channel to absorb differences between the exclusive and the inclusive cross section, re-fit, and add the exclusive data for that channel. These steps are repeated until all open-bottom channels have been included, with the "B sBs " channel, used to absorb missing intensity, added last. For the "B sBs " channel, we fix the non-diagonal background parametersf µ,ν (µ = ν) to zero. For all other open-bottom channels, all background termsf µ,ν are free parameters (see Eq. A32 for the definition of the couplingf ).
In addition, the bare masses m R and couplingsĝ R:µ are free parameters in the fit, apart from the couplings of the lowest mass pole for which we fixĝ R:µ = 0 for µ = BB owing to the experimental observation that the Υ(4S) is fully saturated by its decay to BB.
Based on these fits, we add the three-body hiddenbottom processes one-by-one and re-fit the data for each step. Here, we use the three different approaches of adding either quasi two-body channels or treating the three-body channels perturbatively in one of the final state matrices F detailed above. In the former case, we fix the background terms in the K-matrixf µ,ν to zero if µ or ν is a three-body channel, but allow for a free background term in the production P -vector.
To determine statistical uncertainties of the fit, we use a bootstrapping method, generating 1000 pseudo datasets by randomly varying all data points following Gaussian distributions G(σ, δσ) where σ and δσ are the measured cross section and its uncertainty at any given s. Given that Belle did not use a prior requiring σ to be positive, we allow random variations to negative cross sections as well. For each pseudo dataset, we repeat the fit starting from the nominal solution for a given model choice. To cross check our result, each fit result to a pseudo dataset is tested as a new set of starting parameters in the fit to the actual data, and if that fit yields a better result, we replace the previous solution and repeat the determination of uncertainties with the new fit result. In each channel, we obtain confidence levels for each center-of-mass energy containing the central 68% (90%) of all fit variations for a given model. We remark that the nominal fit result does not necessarily lie in the center of the interval, in large part because the fit result is positive by definition, whereas the (pseudo-)data are allowed to have negative cross section values.
To extract pole locations, in those cases where the e + e − channel is contained in the production vector, we extend K and M by one dimension using K µ,e + e − = P µ and K e + e − ,e + e − = R g 2 R:ee m 2 R −s and systematically scan M e + e − ,e + e − (s) in the complex plane. We look for poles on the nearest sheets defined such that the sign of the imaginary part of the breakup momentum k µ is negative for channels whose thresholds are below (m thr < (m pole )) and positive for channels whose thresholds are above the pole mass. If a candidate for a pole in the complex plane is found, we repeat the scan around the pole with a finer grid. This procedure is repeated multiple times and only those candidates that can still be identified as poles on a grid with eV binning are examined further.
Once pole candidates are identified, we determine the residues Res µ using a discrete version of Cauchy's residue theorem (this aids in producing stable results): where the sum over n is chosen such that the discrete values of s n run along a closed circle around the position of the pole candidate m 2 pole in the complex plane. These residues are related to the partial width of the relevant resonance for channel µ via Γ R,µ = |Res µ | R .
It is possible for "ghost poles" to appear in the formalism. This is particularly true for the Υ(4S) region where the fit tends to arrange the background terms such that [1 + KC] is not invertible for s = m 2 ghost pole . We identify ghost poles by comparing their total widths as determined by the residues to twice the imaginary part of the pole location. Lack of agreement is evidence for a ghost pole. We also adjust the strengths of the couplingŝ g andf by a factor λ < 1. As this factor is reduced to zero, poles should approach the bare mass on the real axis, while ghost poles can exhibit other behavior. We have found that both methods yield valuable information; for example, most ghost poles have total widths that are O(10 5 ) times too small with respect to the imaginary part of the pole location.

V. RESULTS
In this section, we first present fit results for our nine model variations (Sec. V A). We then analyze the structure of the solutions in the complex plane. This allows us to determine resonance masses and total widths by extracting pole locations (Sec. V B) and to determine branching fractions and partial widths by calculating the residues of each pole in each channel (Sec. V C). As a byproduct, we also plot two-body to two-body hadronic cross sections (Sec. V D).
A. Fits Fig. 1 presents experimental data and fit results for the quasi two-body model with β = 1.0 GeV for the 10 channels under consideration. The inner (red) and outer (green) shaded regions show the central 68% and 90% confidence levels as determined with the bootstrap method described in Section IV. As can be seen, the fit is quite good, with a χ 2 per degree of freedom of 1.17. Note, however, that the paucity of data above √ s = 11 GeV leads to ambiguity in the large energy behavior of the exclusive channels, subject to the constraint provided by the inclusive reaction. Furthermore, the dummy channel (panel (d)) takes up substantial strength, especially near the resonances at (approximately) 10.9 and 11.0 GeV. The question of what this channel represents will be considered below. We remark that the spike seen in the left of panel (a) is the trailing edge of the Υ(4S).
The fit significantly underestimates the data near 10.65 GeV in the B * B channel (panel (b)). This is caused by the low value of σ bb in this energy region (see panel (l)), where evidently the inclusive data does not agree with exclusive measurements [28].
Fit results for all nine models are displayed in Fig. 2. One sees that the fits are quite consistent through the regions in which they are constrained. Many of the channels show the ambiguity in the high energy region men-tioned above, which is to be expected due to the absence of exclusive data. Chi-squared values for all the fits are given in Table II, showing similar fit quality for all the models (with perhaps a slight preference for the threebody non-resonant class of models).

B. Poles
A summary of the extracted pole positions for each model is presented in Table III. The poles are labelled Υ(4S), Υ(10750), Υ(5S), and Υ(6S) because they clearly reproduce the states reported in the RPP [7]. (The RPP names for these states are Υ(4S), Υ(10753), Υ(10860), and Υ(11020) respectively.) We also report 68% confidence intervals. In a few instances the sheet structure of the corresponding pole is ambiguous. In those cases, multiple pole locations are reported with the sheet being indicated by the signs on the imaginary part of the breakup momentum of the five open-bottom channels (ordered by threshold). Intervals reported in square brackets correspond to alternative solutions found during the bootstrap procedure.
Those T -matrix poles that correspond to poles of the K-matrix are displayed in color in Fig. 3. Each point represents a pole from a fit to the bootstrap pseudo-data for different model choices. Poles that we identify as spurious, but which nevertheless have sizable residues, are shown in gray. The RPP estimates of mass and width (corresponding to points M − iΓ/2) are also shown as stars with error bars.
All models and bootstrap variations agree quite well for the Υ(5S) and Υ(6S) pole positions, with the Υ(5S) agreeing with the RPP estimate, while our fits obtain a total width for the Υ(6S) that is approximately twice that of the RPP. In contrast, the situation for Υ(4S) and Υ(10750) is less clear.  Turning attention to the Υ(4S), we see that most pole positions cluster near the nominal RPP value, although 10-20 MeV higher in mass. There is, however, a region of the β = 1.0 GeV, three-body non-resonant model poles that lies near ( √ s) = −0.03 GeV. The number of data points is sufficiently sparse (and the model is sufficiently general) that perhaps a variety of nearly degenerate minima of the objective function exist. In this case, the secondary group of poles appear to be associated with BB threshold. We note that the secondary group of poles has substantial overlap with a group of ghost poles (indicated in gray). These points tend to move towards the BB threshold upon rescaling the couplings -indicative of their spurious nature, and hint that the three-body non-resonant secondary group of bootstrap poles should not be considered as viable bottomonium resonance candidates.
For the Υ(10750), the main accumulation of poles seems to agree with the RPP value. However, other solutions cannot be ruled out. The high model dependence found in these fits, yielding poles with a large range of masses and widths, reflects the lack of data in the energy region around 10.7 GeV. Some models contained ghost poles in this region that move towards the thresholds of either the dummy or the B * s B * s channel as the couplings are decreased. New data from an upcoming Belle II measurement in this energy region will be key in determining the properties of the Υ(10750) with higher precision.
Comparing fit results for different model choices: Red lines correspond to the quasi two-body models, yellow lines to the three-body models with resonant and yellow lines to three-body models with non-resonant couplings. Solid lines correspond to β = 1.0, whereas dashed and dotted lines correspond to β = 0.8 and β = 1.2, respectively. The order of plots (a) through (l) is the same as in Fig. 1.

C. Residues
To determine partial widths and branching fractions, we use the method discussed in Sec. IV to calculate the residues for each channel at each pole in the complex plane. Considering the significant model dependence and large statistical uncertainties, we do not quote central values, but only estimate ranges for each measurement that take into account both the statistical and model spread in our solutions.
We first consider electronic widths, which are reported in the top panels of Figs. 4, 5, 6, and 7 for the Υ(4S), Υ(10750), Υ(5S), and Υ(6S), respectively. The results are summarized and compared to theoretical expectations in Tab. V of the next section.
Extracted electronic widths for the Υ(4S) state (Fig. 4) range from 0.003 to 0.62 keV, with most values clustering near 0.15 keV, somewhat lower than the RPP value of 0.272 ± 0.029 keV. A first measurement of the e + e − partial width for the Υ(10750) is reported in Fig. 5, with values ranging from 0.004 to 0.10 keV. These values are somewhat smaller than those for the Υ(4S), an issue to which we return in Section VI C. Figures 6 and 7 show that the situation is somewhat cleaner in the cases of the Υ(5S) and Υ(6S), with all extracted partial widths for both in a range between roughly 0.04 and 0.07 keV.
Our extracted values for both the Υ(5S) and Υ(6S) electronic widths are substantially smaller than the values reported in the RPP, which are 0.31 ± 0.07 and 0.13 ± 0.03 keV, respectively. The original measurements of the inclusive e + e − cross sections and their subsequent parameterization [1,2], which are the basis for the RPP values, were based on very little data and unconstrained models. In Ref. [1], for example, the inclusive e + e − cross section was modeled using a sum of Gaussian distributions, with a single Gaussian distribution covering the entire Υ(5S) region. Our model is more fine grained and better constrained by the addition of more experimental data. These circumstances, and the proximity of the recently discovered Υ(10750), drive the large deviations from the RPP values. The implications of this deviation will be explored in Section VI C.  Branching ratios for hadronic two-and three-body channels are also displayed in Figs. 4 to 7 and are summarized and compared to expectations in Tabs. VI to IX. Once again, the model variants permit moderately large variation in the extracted branching ratios. This can be due to a number of effects. For example, we report a nonzero branching fraction of 0.6 (0.4-0.9)% for Υ(4S) → B * B for the β = 1 GeV quasi two-body model. This is a result of the pole being found very close to the B * B threshold with a sizable width. In other models, the Υ(4S) → B * B branching fraction can grow to around 30%, depending on the location of the pole with respect to the threshold and on the width of the state.
A more subtle example of model variation is visible in Fig. 6, where three branching fractions (circles) for Υ(1S)ππ lie above the other six values. These three points are all models in which β = 1 GeV, which happen to have a narrower Υ(10750) than the other models, hence less overlap with the Υ(5S), and hence larger Υ(1S)ππ branching fractions.
Branching fractions to the "missing" channel show some dependence on the state, with values of (70-90)% (Υ(6S)), (31-77)% (Υ(5S)), and anywhere between zero and 79% (Υ(10750)). It is perhaps expected that these branching fractions diminish as the mass of the resonance decreases. Nevertheless it is disconcerting that so much of the Υ(5S) and Υ(6S) states disappear into unseen channels. Of course, this may be due to increasing phase space and channels with higher √ s. It may also be due to the lack of data above 11 GeV. For example, the putative Υ(6S) signal is truncated in the B * B channel by missing data (see panel b of Fig. 1), and therefore that 6S branching fraction may be missing some intensity. Lastly, we note that both the Υ(5S) and Υ(6S) are reported to decay to Z b π and Z b π, with the two Z b states decaying to B * B and B * B * , respectively. This implies that B * B π and B * B * π could be important components of the "missing" channel.

D. Hadronic Scattering Cross Sections
With the K-matrix in hand it is a simple matter to obtain cross sections for all channels. Here we examine the set of 25 two-body channels to confirm that the model is behaving in a reasonable way. The result for the twobody variant with β = 1.0 GeV is displayed in Fig. 8 along with the central 68% and 90% intervals that are obtained from fits to the pseudo-data samples. The cross sections obtained from all the fit models are shown in Fig. 9.
We see that the cross sections have the scale of typical hadronic cross sections, which is reassuring. All cross sections also display a sharp threshold rise, which is not unexpected. This is typically followed by a rapid drop over a range of √ s ≈ 100 MeV. Again, this is expected [29] as this scale is set by the relevant hadronic wavefunctions. Beyond the threshold region there is scant evidence of the bottomonium resonances, with small peaks visible at Υ(10750), Υ(5S) and Υ(6S) in some channels. This is  likely due to the large background scattering present and the large couplings to the dummy channel. It is interesting that the range of fits in Fig. 8 tend to follow those of Fig. 9. Evidently the different fit models tend to span a similar region in parameter space as the minima obtained in the variant model space. This is a strong indication that the modeling is consistent and that the fluctuations present are due to data quality.

VI. DISCUSSION AND INTERPRETATION
The large bottom quark mass makes bottomonium an ideal system for the application of constituent quark models. It is therefore informative to compare the results we find to model predictions and lattice field theory computations where they exist.  sive channels for large energy. As a result the strength in the inclusive measurement can be shared in different ways throughout the available channels. This observation is reinforced by the results of Fig. 2, which show that the large energy ambiguity is reflected in the model variant optimal fits. Evidently, model variation recapitulates data variance, which is a strong indication of the consistency of our approach. As we have noted, the paucity of data above 11 GeV makes it difficult to pin down branching fractions for the Υ(6S); hence measurements in this region would be most welcome. The other region of greatest model variation occurs near √ s = 10.75 GeV. As might be guessed, this is correlated with ambiguity in the pole location of the Υ(10750). We find very strong evidence for this state with a significance well above 10 σ compared to models using only three poles, where the significance is estimated using the difference in the minimum χ 2 of the respective fits. However, its location varies with model and under bootstrapping. The imaginary part of the pole runs between 10 and 70 MeV (with the largest cluster around 20 MeV). More importantly, the real part of the pole lies above or below the nominal missing channel threshold, which leads to a significant spread of the measured branching

fractions.
The Υ(10750) was discovered in three-body decays [15] to Υ(nS)ππ. Because of this we thought it prudent to determine if evidence for this resonance exists in the twobody open bottom channels. This investigation is reported in Appendix B. We again see strong evidence for the Υ(10750) (similar to Ref. [5]).
Panel (d) of Fig. 2 shows the rate for e + e − to the missing dummy channel. One sees that the reaction is dominated by the Υ(5S) and Υ(6S), with a more ambiguous contribution from the Υ(10750). This result is driven by the lack of agreement between σ bb (panels k and l) and the sum over the measured exclusive channels. As we have discussed above, we interpret the discrepancy as being due to missing multipion final states such as B * B π and B * B * π, including Z b π and Z b π contributions. Panel k of Fig. 2 focuses on the Υ(4S) region in the inclusive cross section. One sees an asymmetric peak and a node near 10.6 GeV that is apparently associated with the B * B threshold (indicated as a dashed line in the panel). We have already noted that the RPP reports an Υ(4S) mass that is approximately 15 MeV below ours. In fact, the real part of our pole position does not lie at the peak of the cross section. The older model, from which the RPP mass is obtained, assumed a Breit-Wigner resonance amplitude (so that the mass necessarily lies near the peak position of approximately 10.58 GeV). Of course, the assumption of a Breit-Wigner form can be improved, as we have done with the K-matrix formalism. The node appears to be due to destructive interfer- ence between the resonance portion of the amplitude and the background scattering, which turns on just at B * B threshold. In the end, unitarity of the model and this effect combine to yield a pole position that is somewhat removed from the Breit-Wigner mass. Lastly, Table II gives the global chi-squared for the nine models considered here. All values lie between 1.08 and 1.19, indicating that the fits are both good, given the disagreement between the inclusive and the sum of exclusive cross sections around 10.65 GeV, and comparable with each other. Again, this shows that model variation is consistent with data variation as revealed in the bootstrap fits, which demonstrates both the reliability of the approach and its limits.

B. Masses
Pole positions are the most basic quantities extracted from the fit, it is thus of interest to compare our results to quark models and lattice field theory computations.
There is a large model literature, which we cannot possibly summarize here. We therefore focus on five models, two recent examples from the literature and three further models that are due to us. These are 1. (GM) a version of the Godfrey-Isgur "relativized" quark model, which assumes relativistic kinetic energies, smeared potentials based on the one gluon exchange interaction, and a Lorentz scalar linear confining potential [30].
2. (SOEF) a nonrelativistic constituent model with a frozen coupling, regulated one gluon exchange short range interactions, and a mixture of Lorentz scalar and vector screened confinement potentials [31]. Use of a screened confinement potential limits the model to low-lying states.
3. (NR) a nonrelativistic model with perturbative spin-dependent interactions reported in Ref. [32]. This model was originally used to describe the charmonium spectrum and was re-fit to 17 well-known bottomonium states here.

(ARM) a relativized model reported in
Ref. [33]. The model uses relativistic kinetic energy and spindependent interactions with a Lorentz structure informed by lattice field computations. The spin-dependent interactions are given by (this corrects typographic errors in the original) The model reported here uses = 0.25 as obtained by lattice computations and was fit to 59 meson masses across all flavor sectors.
5. (bbg) a simple spin-independent model of mesons and constituent gluon hybrid states, first reported in Ref. [34], where it was applied to charmonia and charmonium hybrids. We have refit the model to 17 bottomonium states here.
The final comparison will be made to the recent lattice field calculation of Ryan and Wilson [35]. This work uses 2+1 dynamical quark flavors on anisotropic lattices of size 20 3 × 128 and 24 3 × 128. The pion mass is relatively heavy at 391 MeV. The authors are able to identify the likely composition of states by evaluating overlaps with a variety of operators. Their preferred identifications are shown in the last column of  The table reveals that model predictions have typical deviations of tens of MeV from each other, and from experimental masses, with the situation getting worse as one moves up the spectrum. Certainly, it appears that the Υ(6S) is anomalously light (by 60-100 MeV) with respect to theory expectations. The lightest four masses extracted from lattice field theory line up fairly well (within several tens of MeV) with quark models. The authors of Ref. [35] identify the next two vector states as the (2D) and a hybrid meson, as indicated in the table.
On the experimental side, the authors of Ref. [15] speculate that the Υ(10750) is either a (3D) bottomonium state or a bottomonium hybrid. The approximately 50 MeV deviation from the (3D) quark model is perhaps not unusual; on the other hand, it is rather far removed from the (2D) state as predicted by the quark models. In contrast, the lattice computation did not resolve a D-wave near 10420 MeV and the authors assigned the 10718 state to (2D), which is a possible attribution for the Υ(10750) if the lattice results are confirmed at lighter quark masses. Finally, the lattice and model calculations yield a vector hybrid near 11000 MeV, which is too heavy for a reasonable interpretation of the Υ(10750) signal.

C. Couplings
Constituent quark model computations of electronic widths and two-and three-body decay modes exist. (In principle these can also be obtained in lattice calculations, but this remains to be done.) Table V shows the electronic partial widths for the resonances considered in this study as reported in the RPP, determined by us, and the quark model computations due to Godfrey and Moats [30] and Segovia et al. [31]. We also give results from a quark model calculation that attempts to account for persistent over-estimation of decay constants in nonrelativistic models by softening the short range interaction with the aid of the running coupling [36] (labelled "LS" in the table). The table suggests that the quark models tend to agree on the magnitude of the coupling and that, perhaps not surprisingly, these roughly agree with the measured values as reported in the RPP. However, as discussed above, typical older model assumptions lead to widths that are higher than they should be. Indeed, our results are smaller by approximately a factor of two for the Υ(6S) and a factor of six for the Υ(5S). The value for Υ(10750) is new. The trend is that these widths (and hence, decay constants) fall more rapidly with radial quantum number than previously thought, which implies that short range dynamics in typical constituent quark models is not as well understood as hoped. Future lattice field theory computations of decay constants should help clarify the situation.
We extract an e + e − partial width of (4 -100) eV for the Υ(10750). For most models and bootstrap solutions, we find values substantially larger than the roughly 2 eV that quark models typically obtain for D-wave widths. Possible explanations are that (i) S-D wave mixing is larger than typical quark models predict, (ii) the D-wave identification is incorrect, or (iii) quark model descriptions of short range dynamics become worse as the an-gular momentum increases. In view of the large overall spread in our fit results, it is perhaps too early to speculate on this discrepancy; although we suspect that option (iii) should be considered further by modellers.
Tables VI -IX report branching fractions for the four resonances as determined here, RPP values, and theoretical expectations from the GM and SOEF quark models.
As expected, in most of the models and bootstraps, nearly the entire width of the Υ(4S) is in the BB mode (Tab. VI). As discussed above, we do however report a nonzero value for the B * B mode because a number of bootstrap poles lie above this threshold. Our threebody branching fractions are two to four times larger than those given in the RPP; however our ranges are quite large. These higher branching fractions are a consequence of the pole position being shifted towards higher masses. The three-body branching fractions are constrained by a single Belle data-point on the Υ(4S) resonance. If the pole is indeed shifted, the peak in the cross sections of the three-body processes is shifted as well, leading to a larger branching ratio when compared with the assumption, that the energy point is in the maximum of the resonance. The quark model three-body branching fractions are in the ranges of both the RPP and our measurements, and it is perhaps too early to make judgements on the efficacy of the models.    times larger than that reported in the RPP. We have attributed this substantial missing strength to B * B( * ) π, potentially through the Z ( ) b intermediate state; although we remind the reader that experiments reported ratios of cross sections rather than branching fractions in this case.
Other branching fractions also show substantial disagreement with the RPP values, with the largest discrepancy in the B * B channel. Comparison to the quark models is somewhat confused because they do not account for all possible channels in determining the total width. Of course, relative branching fractions remain meaningful. Both models predict that B * B is the dominant decay mode with B * B * an order of magnitude smaller. We find some evidence for the reverse situation, which agrees with the trend in the RPP.
Lastly, we consider the Υ(6S) branching fractions shown in Tab. IX. Once again, much of the full width appears to be unmeasured, with approximately 80% going to the missing channel. We attribute this to missing B * B( * ) π, as discussed above. In this case, the quark models agree on the ordering Bf (B * B ) Bf (B * B * ) > Bf (BB). We, on the other hand, find Bf (B * B ) Bf (B * B * ), although we stress that new data on exclusive cross sections above 11 GeV could lead to significant changes in these branching fractions.
Although branching fractions are not reported for the Υ(10750) and Υ(6S) in the RPP, they list ratios of widths, as shown in Tables X and XI. Our results for the Υ(6S) align very well with those reported in the RPP, which supports our estimated branching fractions in Tab. IX. Alternatively, our ratio for the Υ(3S)ππ decay mode of the Υ(10750) is smaller than that of the RPP over the full model range. The three-body branching fractions of Tables VII, VIII, and IX are all comparable in size, which we find reasonable.
Finally, we draw attention to the large branching fractions for closed bottom three-body channels for the Υ(5S) (Tab. VIII) and Υ(6S) (Tab. IX), which are three or more orders of magnitude larger than those seen for the Υ(4S). The formerly perplexing dipion decays of heavy quarkonia are now widely believed to be understood in terms of long-lasting open-bottom fluctuations and mixing with the novel Z b states (see, for example, Ref. [37]). The challenge here is that the Z b mesons have masses that are rather far removed from the Υ(5S) and Υ(6S) masses. Relevant box diagrams containB ( * ) B ( * ) mesons and are therefore also far removed in mass. Of course the possibility of higher mass resonances in the box diagram is open, and should be investigated.
In summary, it appears that the model predictions for open-bottom decay widths for Υ(5S) and Υ(6S) do not agree well with our extracted residues, although the uncertainties are large. Both groups (GM [30] and SOEF [31]) used the well-established "3p0" model for strong decays to obtain their results. Although the model is very simple, the ratios of partial widths that are related by spin symmetry should be robust. The lack of agreement could then imply that the assumed spin structure of the 3p0 model is incorrect. An alternative is that these partial widths also depend on wavefunction overlaps, and these can be particularly difficult to model high in the spectrum, where details such as the location of wavefunction nodes can become very important.

VII. CONCLUSIONS
We have provided the first comprehensive analysis of vector bottomonia using data on both the production of various exclusive open-and hidden-bottom channels and the inclusive cross section for bottom anti-bottom quark pair production in electron positron annihilation. In a unitary approach using the K-matrix formalism, we find a total of four poles, the Υ(4S), Υ(10750), Υ(5S) and Υ(6S) in accordance with the RPP. Allowing for a nonresonant contribution and a proper treatment of thresholds, we extract a pole position of the Υ(4S) that is about 10-20 MeV higher than previous values. While we find a significant contribution of the additional Υ(10750) state previously observed by Belle in e + e − → Υ(nS)π + π − , the paucity of the data in that energy region allows for a broad range of masses and widths, all describing the data similarly well. The large number of two-body thresholds in the Υ(4S) and Υ(10750) regions further complicate the situation. The pole positions of the Υ(5S) and Υ(6S) are well constrained by the (inclusive) data, with the Υ(5S) position agreeing well with the RPP average and the Υ(6S) being about twice as wide as previous measurements.
The K-matrix description allows us to determine absolute branching ratios of all four states for the first time, although relatively few data points in the exclusive measurements and (potentially as a consequence) a significant model dependence lead to large uncertainties. However, our method of extracting these branching fractions is more robust than the ratios of cross sections currently being reported in the RPP, especially given that none of the exclusive cross sections of open-bottom production exhibit clean peaks in the Υ(10750), Υ(5S) and Υ(6S) regions. For these three states, we find a large fraction of the intensity is still missing, with B ( * )B( * ) π channels, potentially populated by Z ( ) b contributions, identified as promising candidates to explain the discrepancy between the exclusive and the inclusive data. A measurement of these cross sections using Belle(II) data could resolve this issue in the future, thereby significantly reducing model uncertainties.
Additional data in the Υ(10750) region from Belle II is highly anticipated in order to more cleanly study this state, with a more precise measurement of its branching ratios aiding in identifying its nature. Similarly, additional exclusive measurements in the Υ(6S) region above 11 GeV would be beneficial to reduce current ambiguities caused by the total absence of such data. Finally, one might speculate about the nature of the Υ(5S) and Υ(6S) states, given that their decay rate to three-body hidden-bottom channels is significantly enhanced compared to the conventional Υ(4S) state. Whether this can be explained by a strong coupling to the exotic Z To facilitate using the conventional scattering formulas of quantum field theory (as exemplified, for example in the RPP), we choose to work with the invariant amplitude (M) defined by where the index f (i) refers to the final (initial) particles in the reaction. In this case Eq. A2 becomes where we now make the channel indices explicit and introduce a symmetrization factor, S γ . The right hand side comes from γ T µγ T * νγ where the notation is resolved as a sum over channels (denoted γ) and an integral over three-momentum. The result is written in terms of the invariant phase space dΦ (this is (2π) 4 larger than the RPP convention) defined by (A5) For two bodies A → C + D in the center of mass frame We have introduced the center of mass momentum k * (M 2 A ) which we generalize to We stress that this quantity will not be considered below threshold in the following (except when exploring the Tmatrix in the complex plane) for reasons to be discussed shortly. Thus it is implemented with a theta function forcing it to zero below threshold. In keeping with particle physics convention, the branch cut in the square root will be taken along the positive real axis (however, the usual branch cut along the negative real axis will be used for the square root in the denominator, which exhibits a non-physical singularity at s = 0). Now restrict attention to two-body scattering in the center of mass frame with two-body intermediate states γ. We get Finally, assume that the amplitudes are not functions of angle so that this equation reduces to an algebraic equation. Going to matrix form then gives: is the phase space matrix. Different conventions for the phase space exist -these can be absorbed into the definition of K and will not affect pole positions but will yield unconventional expressions for cross sections and other physical quantities. Eq. A9 is equivalent to which implies that the quantity in brackets is real and symmetric. It is therefore useful to identify a matrix with the same properties: where R is a real function along the diagonal of K. The function R can in turn be incorporated into the phase space by defining a new quantity, C, as Thus we have Many choices for the function C are possible (and many are made in the literature). However it is preferable to use one that maintains analyticity across thresholds (hence R = 0) since this smooths discontinuities that fitting routines find difficult to handle. In this paper we will therefore set C equal to the Chew-Mandelstam function [18,39]. The real part of the Chew-Mandelstam function is related to its imaginary part by a oncesubtracted dispersion integral, that ensures a smooth transition across threshold. Here s γ = (m 1 +m 2 ) 2 is the threshold for channel γ. The auxiliary function is chosen to guarantee that the imaginary part of C γ is ρ γ . Thus r is the continuation of ρ to the complex plane: The integral can be done yielding (we set C(s γ ) = 0) [39], where the additional function is defined by Fig. 10 displays the real and imaginary parts of the Chew-Mandelstam function. We remark that many authors choose to continue k * (and thus ρ(s)) below threshold in the name of "analyticity". This amounts to making a model choice of the real part of C that is particularly poor, as illustrated in the figure, and we recommend against this practice. Some authors choose to include angular momentum factors in the Chew-Mandelstam function [40]. The derivation we have presented leaves no room for modifications of this sort. And in fact, it is clear that these factors are associated with the scattering amplitude rather than the phase space, and are better modeled in the couplings, as described below.
The remaining task is to specify the elements of K. Gribov has argued that scattering amplitudes factorize near poles [41], hence it is sensible to write the elements of K in terms of the product of channel couplings and a "bare" pole. This is common practice in the field, which we follow here. Thus we parameterize K with resonant and non-resonant ("background") terms: The index R refers to a resonance and Greek indices refer to continuum channels. The couplings must be real but can be s-dependent (angular dependence is excluded by the assumptions leading to Eq. A15). Extracted couplings are spin-averaged for g R:ee and are spin-summed in other cases. Many choices for the s-dependence of the couplings exist. To facilitate the discussion and to gain insight into these choices it is helpful to compare to the propagator for a scalar field. To this effect we consider a heavy field φ with mass M coupled to a light scalar with mass m. The interaction lagrangian is taken to be 1 2 g d 4 x φϕ 2 . Performing the standard Dyson sum yields a scattering amplitude of where −iΠ is the φ self-energy given by (we use the M S scheme to renormalize): Π(s) = g 2 16π 2 S g 2 1 − s γ /s · θ(s − s thr,γ ) = −g 2 ρ γ (s).
(A23) At lowest order, the φ particle width is given by which can be immediately generalized to the case of twoto-two scattering in the center of mass frame as Γ(s) = g 2 k * (s) 8πSs .
Notice that, because the imaginary part of the selfenergy is (up to −g 2 ) the phase space, ρ, and because the loop integral is analytic, the Chew-Mandelstam function For two-to-three scattering, 12 → 345, the differential cross section is given by: then the perturbative relationship to the coupling is given by .
Similarly, the perturbative three-body partial width is given by where A is the area of the Dalitz plot for this decay.
In this work we do not rely on perturbative expressions, rather partial widths are extracted from residues, as we have described in Sec. IV. In this regard, two useful relationships are We present fit results for the case in which only twobody channels are considered. This is of interest because the evidence for the Υ(10750) was in three-body decay modes [15] and we wish to examine the possibility that this state can be discerned in the two-body data. We have therefore examined the two-body system with three models defined by β = 0.8, 1.0, and 1.2 GeV. Fit results are shown in Fig. 11. We note that the general features of the fits are similar to those of the full model. Paucity of data, however, is evidenced by nearly degenerate minima that are found for the β = 0.8 and 1.0 GeV cases. Global χ 2 values are given in Tab. XII, where we see very similar results (if slightly worse) to the full case.
Pole positions for the restricted model are given in Fig.  12. Once again, results are very close to those of the full model.
These results yield strong evidence for the Υ(10750), with a similar range of resonance parameters as when three-body channels are included in the dataset. The  minimum χ 2 in fits with three poles is at least 150 units worse with 6 fewer parameters than in fits with four poles, again leading to a statistical significance of larger than 10 σ. We conclude that the evidence for the Υ(10750) is robust and convincing.