Excited bottomonia in quark-gluon plasma from lattice QCD

We present the first lattice QCD study of up to $3S$ and $2P$ bottomonia at non-zero temperatures. Correlation functions of bottomonia were computed using novel bottomonium operators and a variational technique, within the lattice non-relativistic QCD framework. We analyzed the bottomonium correlation functions based on simple physically-motivated spectral functions. We found evidence of sequential in-medium modifications, in accordance with the sizes of the bottomonium states.


Introduction
Quarkonium suppression has been proposed as a signature of quark-gluon plasma (QGP) formation in heavyion collisions [1]. The main idea behind this proposal was the observation that color screening within a deconfined medium can make the interaction between the heavy quark and anti-quark short ranged, leading to the dissolution of quarkonia in QGP. At a given temperature, different quarkonium states are expected to be affected differently by QGP-a more tightly bound quarkonia having a smaller size is less influenced by the medium than a relatively loosely bound, larger one. Therefore, following the hierarchy of their binding energy and sizes, a sequential pattern of in-medium modification is expected [2,3]. Some evidence of sequential in-medium modification of quarkonia comes from lattice QCD studies of S-wave and P-wave quarkonium correlators along the temporal [4][5][6][7][8] and spatial [9,10] directions. Recent studies have revealed that inclusion of dissipative effects lead to a more complex theoretical picture of in-medium heavy quark and antiquark interactions [11,12]. However, the main conclusion, i.e., quarkonium dissolve in QGP when the temperature is high enough compared to its inverse size and binding energy, have remain unchanged [13][14][15][16][17]. Hints for sequential in-medium modification of bottomonia have also been observed in heavy-ion collision experiments [18][19][20][21][22]. While the connection between the observed hierarchy of the Υ(nS) yields in heavy-ion collisions and the expected sequential melting of these states in QGP is complicated by dynamical effects, such a link is expected to exist [23][24][25]. For this reason, the study of sequential in-medium quarkonium modifications in heavy-ion collisions is a subject of extensive experimental and theoretical efforts; for recent reviews see Refs. [26,27].
There have been many attempts to study in-medium properties of charmonium [28][29][30][31][32][33][34][35][36] and bottomonium [4-8, 33, 37, 38] in lattice QCD, almost entirely focused on in-medium modifications of ground states of S-and P-wave quarkonium. Previous lattice QCD studies of in-medium quarkonium used point meson operators, i.e., operators with quark and anti-quark fields located in the same spatial point, which are known to have non-optimal overlap with the quarkonium wave-functions, especially, with the excited states. As a result, these correlators are largely dominated by the vacuum continuum parts of the spectral function, and isolating the contributions of in-medium bottomonium become quite difficult [16,17,39]. Recently, we explored the possibility of studying in-medium bottomonium properties using correlators of extended meson operators [40]. We found that such operators have very good overlap with the lowest S-and P-state bottomonia, thereby, allowing us to cleanly isolate the vacuum continuum contributions to the bottomonium correlators. We showed that these correlators are more sensitive to the inmedium bottomonium properties than with point sources. Analyzing these correlators we found evidence for thermal broadening of the 1S-and 1P -state bottomonia, however, excited bottomonium states remained elusive. In this letter we introduce novel extended meson operators within the lattice non-relativistic QCD (NRQCD) formalism, which, for the first time, allow us to probe in-medium modifications of up to 3S and 2P bottomonium.

Methodology
The lattice NRQCD Lagrangian employed in this study is exactly the same as in Refs. [40,41]-tree-level tadpole improved, accurate up to order v 4 , but also includes v 6 spin dependent terms. For the background gauge fields we used 2 + 1-flavor gauge configurations on 48 3 × 12 lattices with bare gauge couplings β = 6.74, 6.88, 7.03, 7.  pling we also carried out the corresponding vacuum T = 0 calculations. All gauge configurations were generated by the HotQCD collaboration [42,43], with the physical value of the strange quark mass, and up/down quark masses corresponding to the pion mass of 160 MeV in the continuum limit. The lattice spacings, a, were determined using the r 1 scale from the static quark anti-quark potential and the value r 1 = 0.3106 (18) fm [44]. The mass parameter in NRQCD Lagrangian was fixed through the kinetic mass of the η b meson, described in detail in Ref. [40].
To calculate bottomonium correlators we used novel extended meson operators in Coulomb gauge of the form The different choices of Γ used in this work can be found in Table 1 of Ref. [40]. The index i refers to the different states in a given channel, e.g., 1S, 2S, 3S etc . The shape-functions Ψ i were obtained by solving the discretized Schrodinger equation with a Cornell potential on a 3-dimensional lattice having a lattice spacing and a volume exactly same as that of the corresponding QCD background. Spin interactions were neglected. For the lattice Schrodinger equation we used an O(a 2 )-improved Laplace operator. For the Cornell potential the string tension was chosen to be (468 MeV) 2 , and the Coulomb part was computed at tree-level in lattice perturbation theory for the Symanzik-improved lattice gluon action with a fixed strong coupling constant α S = 0.24. The bottom-quark mass was set to m b = 4.676 GeV. More details can be found in the Appendix D of Ref. [41]. The Ψ's used for calculations of η b and Υ correlators for lattices with a = 0.0655 fm are shown in Fig. 1.
Since Ψ i is a good approximation for the wave-function of the i th vacuum bottomonium, as expected, the corresponding operator O i was found to have a good overlap with the i th state. However, the off-diagonal correlators, were found to be nonzero, though small. Thus, we resorted to the variational analysis by considering linear combinationsÕ α = Ω αj O j such that Õ α (t)Õ † β (0) ∝ δ α,β . The matrices Ω αj were obtained using the generalized eigenvalue problem [45][46][47][48][49] G ij (t)Ω αj = λ α (t, t 0 )G ij (t 0 )Ω αj . We calculated the cor-relators of optimized operators C α (t) = Õ α (t)Õ † α (0) for 1S, 1P, 2S, 2P and 3S. When calculating Ω αj , the value of t was chosen such that it corresponds to physical extent τ = t · a 0.5 fm and t 0 was chosen to be τ = 0 fm. Choosing t 0 /a to be 1 or 2 does not change the results significantly.

Vacuum case
In Fig. 2 we show some examples of the effective masses, aM eff (t) = ln[C α (t)/C α (t + 1)], at T = 0. Those were found to reach a plateau for τ 0.4 fm for all states. For the excited states we see stable plateau up to τ of around 1.2 fm, beyond which signal deteriorated. Performing single-exponential fits in the plateau region we extracted the energy levels of different bottomonium states. Since the NRQCD energy levels contain an additional lattice spacing dependent constant, the differences of these energy levels have the physical interpretation of mass differences of different bottomonium states. Thus, in this work we set the spin-average energy of 1S bottomonium, E 1S = (E η b + 3E Υ )/4, to be the zero of the mass/energy, and quote masses/energies of the rest of the bottomonium states with respect to these baselines for each lattice spacing. (Note that, in Ref. [40] this baseline was set by the η b energy level.) The reason beingĒ 1S remains unaffected by the spin-spin interaction, which is difficult to reproduce accurately using tree-level NRQCD Lagrangian used in this study. The energy differences shows a very mild dependence on the lattice spacing, and in most cases can be fitted with a constant to obtain our final estimate for the energy differences. The only exceptions are the energy of χ b0 and the 2S hyperfine splitting, where the a-dependence cannot be neglected. We also fitted the resulting energy differences with constant plus a term proportional to a 2 to remove the remaining discretization effects, this resulted in small differences in the central value but larger statistical errors for the final estimate. We used this procedure to compare with the experimental results, except for 3S state, where it gives too large statistical errors. The comparison of the zero temperature results on the mass differences with the experimental result from Particle Data Group  (PDG) [50] is shown in Table 1. We see very good agreement between our lattice NRQCD calculations and the experimental results within the estimated errors in most cases. The only exception is the χ b0 (1P ) state, where a small tension between our result and the PDG value is observed. Using the result for the 3S hyperfine splitting we can also predict the mass of the η b (3S) state to be 10341.8 ± 6.2 MeV, that has not yet been observed experimentally.

In-medium case
In NRQCD the spectral function, ρ(ω, T ), is related to the Euclidean time correlation function: Here, α labels the bottomonium operator of interest. Bottomonium states correspond to peaks in the spectral function having some in-medium width. At large ω many states contribute to the spectral function forming a continuum. Therefore, we can write the spectral function as with the second term parameterizing the continuum part of the spectral function. In the zero-temperature limit ρ med α (ω, T ) = A α δ(ω − M α ), M α being the mass of the corresponding bottomonium state. Here we note that the use of extended operators reduces the relative contribution of ρ high α [51, 52]. For this reason the effective masses in Fig.  2 approach a plateau at relatively small τ . The continuum part, ρ high α (ω), is expected to be temperature independent. This was seen to fit with our calculation, as the temperature dependence of C α (τ, T ) for τ 0.3f m fm was very small, with the small deference being in agreement with changes to the medium. Thus, following Ref. [40], for each lattice spacing we can identify the contribution of ρ high α (ω) to the correlator, C high α (τ ), as Here, A α and M α are the amplitude and mass of the corresponding bottomonium state, and C high α (τ ) is the Laplace transform of ρ high α . Using the single-exponential fits to the vacuum correlators for τ 0.6 fm, and subtracting off this contribution from C α (τ, T = 0) we isolated C high α (τ ) for each value of β. Further, following Ref. [40], for each temperature we then defined the continuum-subtracted correlator as Thus, the continuum-subtracted correlator, C sub α (τ, T ), is mostly sensitive to ρ med α (ω, T ), encoding the in-medium bottomonium properties. We then studied the in-medium bottomonium properties using the continuum-subtracted effective masses, aM sub eff (τ, T ) = ln C sub α (τ, T ) /C sub α (τ + a, T ) .
In Fig. 3 we show typical examples of M sub eff as a function of τ -for the Υ states at T = 251 MeV and for the   χ b0 states at T = 199 MeV. At small τ the M sub eff are equal to the vacuum masses. As τ increases we see an approximately linear decrease of M sub eff . Finally, for τ 1/T we see a rapid drop-off. Similar behaviors of M sub eff for the ground states was also observed in the previous study using Gaussian smeared meson operators [40]. As discussed in Ref. [40], the slope of the linear decrease of M sub ef f can be understood in terms of a thermal width. We see that the slope is larger for higher excited bottomonium states, i.e., the thermal width of different bottomonium states follows the expected hierarchy of their sizes. Higher excited states have larger size and therefore are more affected by the medium, leading to larger width. The behavior of the effective masses at τ 1/T is related to the tail of the spectral function at small ω, and may depend on the choice of the meson operator [40]. Therefore, it is important to compare the results on the subtracted effective masses obtained with different meson operators. In Fig. 4 we compare M sub eff of subtracted Υ(1S) at T = 250 MeV with the corresponding results obtained with Gaussian smeared sources of Ref. [40]. Good agreement was found between the present results and those obtained of Ref. [40], especially for τ 1/T . Therefore, our conclusion regarding the in-medium modification of the spectral functions is not af-fected by the choices of meson operators. For τ 1/T we see a smaller drop-off in M sub eff compared to that observed in Ref. [40]. Thus, the small-ω tail of the spectral functions plays a less prominent role here. We found that the behaviors of M sub eff for η b (nS) are very similar to that of the Υ(nS), and that of the χ b1 (nP ), χ b2 (nP ) and h b (nP ) are very similar to that of the χ b0 (nP ).
As introduced in Ref. [40], the simplest theoretically motivated parameterization of the in-medium spectral function that can describe the generic behavior of M sub eff observed here is as follows The first term in the above equation provides a simple parameterization of the low-ω tail of the spectral function.
As explained in Ref. [40], this tail is important for understanding the behavior of the effective masses around τ 1/T . The second term gives rise to the linear behavior in τ of M sub eff , with the slope given by Γ 2 α . For each temperature, we fitted M sub eff (τ ) with the Ansatz given by Eq. (7), and using Eqs. (2), (6), to determine the inmedium masses, M α (T ), and width, Γ α (T ), of different bottomonium states. Since the the tail of the spectral function play a less prominent role in the present study, for T ≤ 173 MeV and all temperatures for Υ 1S, we fitted by setting A cut α = 0, and omitting 1-3 data points for the largest values of τ . Only for higher temperatures was the term proportional to A cut α included. We generally find good fits with χ 2 divided by degrees of freedom being around 0.5. Still we do see some cases where the fluctuations in the data were sometimes larger than the error bars. This was mostly the case for the 1S state of Υ, but was also observed in χ b0 at 199 MeV as can be seen on the right in Fig. 3. In these cases we found that χ 2 divided by degrees of freedom was around 2. The fit still seem to work nicely, so it is most likely the errors that were a bit too small.
The change of the in-medium mass parameter com-   Figs. 5 and 6, respectively. The in-medium masses of different states obtained from the fits turned out to be very similar to the vacuum masses. In fact, we do not see any statistically significant deviations from the T = 0 results. On the other hand, Γ α (T ) shows a clear increase with increasing temperature. For large enough temperatures, Γ α (T ) appears to approximately rise linearly with T . Nearly for the entire T -range, the in-medium width was found to follow the sequential hierarchical pattern according to the increasing sizes of the bottomonium states: As a result, at these temperatures 2S and 3S, as well as the 1P and 2P states will together appear as broad structures in their respective spectral functions. These observations lead us to conclude that, similar to what have been observed in the experiments [18,19], for T 200 MeV it will become difficult to individually identify the 2S, 3S, 1P and 2P states within the experimentally measured line shapes of the invariant-mass distributions.
Lastly, we address the question to what extent the estimated thermal widths of bottomonium states depend on the model for the spectral function used to interpret our lattice QCD results presented here. Since M sub eff (τ ) show a linear in τ behavior and we do not observe any significant thermal mass shift, following Ref. [40], we also used the following model for the spectral function: Here, M 0 α is the vacuum bottomonium mass, and the parameters A cut α and ω cut α describe the low-ω tail of the spectral function. For T ≤ 173 MeV again we used A cut α = 0. The equivalent thermal width in this case is Γ α (T ) = 2/3∆ α (T ). Carrying out fits with the above Ansatz we obtained thermal widths that, within errors, agreed with the ones obtained by using the Gaussian Ansatz. Therefore, our estimates of thermal width do not depend very much on the precise functional form of the fit Ansatz.

Conclusion
For the very first time, we studied in-medium properties up to 3S and 2P excited bottomonium states using lattice QCD at temperatures T 150 − 350 MeV. This lattice QCD study was made possible through the introduction of novel bottomonium operators within the lattice NRQCD framework, and implementation of a variational analysis based on these novel operators. We found that the effective masses constructed out of the continuumsubtracted bottomonium correlation functions drop off linearly in Euclidean time. We argued that the behaviors of the continuum-subtracted effective masses can be understood in terms of a couple of theoretically-motivated, simple models of the bottomonium spectral functions. For all of the models considered, we found indications of thermal broadening of bottomonium states in QGP. For the entire temperature range, the magnitudes of the thermal broadening were found to follow the expected sequential hierarchical pattern according to the increasing sizes of the bottomonium states. Further, we found that for T 200 MeV the thermal broadening of the 2S, 3S, 1P and 2P states becomes large enough that it would be difficult to identify these states separately within the corresponding spectral functions.