Thermal Broadening of Bottomonia: Lattice Non-Relativistic QCD with Extended Operators

We present lattice non-relativistic QCD calculations of bottomonia correlation functions at temperatures $T \simeq 150-350$ MeV. The correlation functions were computed using extended bottomonia operators, and on background gauge-field configurations for 2+1-flavor QCD having physical kaon and nearly-physical pion masses. We analyzed these correlation functions based on simple theoretically-motivated parameterizations of the corresponding spectral functions. The results of our analyses are compatible with significant in-medium thermal broadening of the ground state S- and P-wave bottomonia.


I. INTRODUCTION
At temperatures above the chiral cross-over temperature T c = 156.5 ± 1.5 MeV [1] hadrons start to melt and chiral symmetry gets restored. However, quarkonia, bound states of charm or bottom quark and anti-quark can exist at higher temperatures. It has been conjectured by Matsui and Satz that color screening will eventually lead to melting of quarkonia, and suppression of quarkonium production in heavy ion collisions can serve as signals of formation of deconfined medium in such collisions [2]. The study of in-medium properties of quarkonia and their production in heavy ion collisions is an extensive research program, see Refs. [3,4] for recent reviews. For interpretation of experimental results on quarkonium production in heavy ion collisions it is necessary to know quarkonium properties at non-zero temperature.
Reconstructing the spectral function from discrete set of lattice data points for C(τ, T ) is a difficult task. Early attempts in this direction have been reported in Refs. [5][6][7][8][9][10][11][12][13]. It has been pointed out that Euclidean correlation functions have limited sensitivity to the in-medium quarkonia properties and/or their melting [14,15], because the maximal time extent is limited to τ max = 1/(2T ) and becomes small at high temperature (c.f., Eq. 1). Studying bottomonium on the lattice is also challenging because of the large bottom quark mass, M b ∼ 5 GeV, which implies large discretization effects ∼ aM b , with a being the lattice spacing. To circumvent this problem an effective field theory approach can be used by integrating out the energy scale related to the bottom quark mass. The corresponding theory is called the lattice nonrelativistic QCD (NRQCD) [16,17]. The heavy quark fields in this theory are represented by Pauli spinors coupled to gauge fields, and quark-anti-quark pair creation of the heavy quark is not allowed. Lattice NRQCD has been used to study bottomonia properties at zero temperature [18][19][20][21][22][23][24][25]. Attempts to calculate bottomonia spectral functions at non-zero temperatures using lattice NRQCD have also been reported [26][27][28][29][30][31][32]. Since pair creation is not allowed in NRQCD, the quark fields do not satisfy periodic boundary condition [33]. Therefore, the relation of the Euclidean time correlator and the spectral function has the form [26] C(τ, T ) = +∞ −∞ dωρ(ω, T )e −τ ω . ( As a result, the maximal time extent at non-zero temperature is τ max = 1/T , i.e., twice larger than in the relativistic case. For this reason, the meson correlators in NRQCD are more sensitive to in-medium properties of quarkonium [26,[30][31][32]. Till now, all the lattice NRQCD studies at non-zero temperature have used point-like meson operators for computations of the quarkonia correlation functions. The use of point-sources for the non-relativistic quarkonia is appealing, since the corresponding spectral function has the physical interpretation in terms of di-lepton rate, and the melting of quarkonia is easily understood as disappearance of peaks in the corresponding spectral functions. However, it is well-known that point-like hadron operators, usually, do not have good overlap with a particular meson state in that quantum number channel, and precise results C(τ, T ) for large τ are needed to isolate the contribution of a particular meson to the correlation function. As a result, point meson operators are not very sensitive to in-medium modifications of quarkonia.
Extended operators are widely used to study quarkonia properties at zero temperature [18][19][20][21][22][23][24][25]. Correlators of extended operators have a better projection onto quarkonium states and are less sensitive to the continuum, i.e., the bound-state contribution to these correlators is sig-nificantly larger. For this reason, we expect that correlators of extended meson operators will be more sensitive to in-medium quarkonia modifications and/or their dissolution. In this paper we explore the temperature dependence of bottomonium correlators corresponding to extended meson sources within lattice NRQCD. As we will see later, compared to the correlators of point meson operators, the correlators of extended sources have significantly larger temperature dependence. We try to understand the observed temperature dependence of the correlation functions of extended meson operator using simple theoretically-motivated parameterization for the in-medium spectral functions, and find evidence for thermal broadening of bottomonium states.
The rest of the paper is organized as follows. In section II we describe our setup of lattice calcualtions. We present our results at zero temperature in section III. Non-zero temperature case and bottomonia properties in the deconfined medium are discussed in section IV. Section V contains our conclusions.

II. DETAILS OF LATTICE QCD CALCUALTIONS
In NRQCD the heavy quarks and anti-quarks are described in terms of two component Pauli spinors, ψ and χ. In our study, we use tree-level tadpole-improved NRQCD action. The tadpole improvement means that the gaugelink variables that enter the NRQCD Lagrangian on the lattice are divided by the average value of the link, u 0 . The NRQCD action used in this study includes all terms of order v 4 as well as v 6 spin-dependent terms, with v being the heavy-quark velocity inside quarkonium.
The calculations in NRQCD are set-up as an initial value problem. The heavy-quark propagator, is calculated as follows Here, t = τ /a and U 4 (t) is the temporal link. The is the leading order non-relativistic Hamiltonian, and is the covariant lattice Laplacian. Here, δH is the part of the NRQCD Lagrangian that contains corrections of order v 4 and spin dependent v 6 correction. The explicit form of δH is given in Ref. [20]. The parameter n plays the role of the discretization time step and is called the stability parameter or the Lepage parameter. For the tadpole parameter u 0 we use the forth root of the plaquette. The anti-quark propagator is also determined by Eq. (4), since it can be related to quark propagator as G χ (x, t) = −G † ψ (x, t). In the above equations we use the same notation as in Ref. [20].
In this paper we are mostly interested in the correlators of extended meson operators of the form O =χ † Γψ, obtained from smeared quark and anti-quark fields, ψ = W ψ , and χ = W χ , with For sufficiently large number of the smearing steps, N , the meson operator has a Gaussian shape having a width σ. The root mean square size of the meson source is r RM S = √ 3σ/2. To avoid instabilities, one has to choose N > 3σ 2 /2. Spin and color structures of the meson are determined by the matrix Γ. Along with the corresponding lattice irreducible representations (irrep), the choices for Γ used in the present study of bottomonia are given in Table I, in terms of the Pauli matrices σ i and the covariant derivative operator We are interested in the meson two-point function, which can be expressed in terms of the quark propagator as In the above equations the gauge links entering the smearing operator W and its conjugate are defined on time-slice indicated in parenthesis, i.e. the gauge link entering W are defined on time slice t, while the gauge link entering W † are defined on time-slice zero. We choose N and σ such that r RM S is about 0.21 fm for all lattice spacings. This choice provides a good overlap with the ground state for both S-state and P-state bottomonia correlators. In our calculations at non-zero temperature we also considered other choices for r RM S , varied between 0.1 fm and 0.4 fm. This is discussed in Appendix A. To make contact with earlier studies we also calculated the correlators of point meson operators that are obtained from Eq. (8) by setting W to unit matrix.
In our study we use 2+1 flavor gauge configurations generated by the HotQCD collaboration using highly improved staggered quark (HISQ) action at the physical strange quark mass and light quark masses corresponding to the pion mass of 161 MeV in the continuum limit [34,35]. Thus, our calculations are performed almost at the physical point. We perform calculations on 48 3 × 12 lattices in the temperature range T = 151 − 334 MeV. These correspond to lattice spacings a = 0.05 − 0.11 fm. For each temperature we perform the corresponding calculations at T = 0. We use multiple sources when evaluating meson correlation functions. The parameters of the lattice calculations, including the lattice sizes, lattice spacings number of gauge configurations, number of sources per gauge configurations, tadpole parameters and the values of the Lepage (stability) parameter are given in Tab. II.
For a given value of the heavy quark mass parameter, M b , we determine the masses of quarkonium states by fitting the large τ behavior of the meson correlation functions with a single exponential. As already mentioned above, correlators of extended meson operators have a better projection onto the lowest state in a given channel. In order to judge to what extent correlators are dominated by the lowest energy state we consider the effective mass, defined as In Fig. 1 we show the effective masses from Υ correlators at T = 0 for extended (smeared) and point sources. For point sources the effective mass reaches a plateau only for τ > 1.2 fm, while for the extended sources the plateau is reached already for τ 0.4 fm, implying that the correlation function is dominated by the ground state contribution at relatively small τ . This will be important for the analysis of the correlation functions at T > 0, where the time extend is limited by the inverse temperature. The highest temperature used in our study is T = 333.5 MeV, which corresponds to τ max 0.6 fm. Thus, we expect to probe quarkonia properties even at the highest temperature using correlators with extended sources. To obtain the T = 0 bottomonium energy levels we fit correlators with smeared sources to a single exponential form in the interval [t min : t max ]. The value of t min is chosen such that it corresponds to τ min > 0.5 fm to eliminate the possibility of excited state contamina- tion. The value of t max is determined by the signal-tonoise ratio of the correlators and the desire to have well controlled statistical errors. In Tab. II we give the values of t min and t max use in the fits.

III. NUMERICAL RESULTS AT ZERO TEMPERATURE
In this section we discuss our results on bottomonia masses and correlators at T = 0 [36], which provide an essential baseline for the study of bottomonia properties in the medium. Before we can study the bottomonia properties we have to fix the mass parameter, M b , in the NRQCD Lagrangian. This can be done by studying the dependence of the energy of the bottomonium state as function of the spatial momentum p, defining the so-called kinetic mass, M kin , of the meson [20]. In this study we use Eq. (10) for the kinetic mass, however, the non-relativistic definition gave very similar numerical results. We determine the kinetic mass of η b meson for a given input mass parameter M in b . Then we interpolate the kinetic mass of the η b in M in b and determine the physical value of the mass parameter by requiring that M kin = M P DG being the experimental value of the η b mass from Particle Data Group (PDG) [37]. This procedure is demonstrated in Fig. 2. In practice, we find that linear interpolation in M in b works well. The values of the tuned M b for different β are given in Tab. II. NRQCD is expected to break down for aM b < 1, because radiative corrections become very large when the inverse lattice spacing is larger than the mass parameter, see e.g. Ref. [22]. In our setup this happens for lattice spacing of about 0.04 fm. This breakdown can be seen in the behavior of the kinetic mass as function of the mass  parameter aM b . The kinetic mass decreases monotonically with decreasing aM b . This decrease is usually well described by a linear dependence, c.f. Fig. 2. For sufficiently small aM b we start seeing deviations from the linear behavior, and eventually the kinetic mass does not decrease with decreasing aM b , but starts to saturate at some value. We see that for β = 7.825 corresponding to a = 0.04 fm we cannot reach the physical η b mass by lowering aM b . For this reason, β = 7.596 was the largest gauge coupling used in our study. For charmonia the largest β value that can be used is β = 6.74. Therefore, in-medium properties of charmonia cannot be studied within the present NRQCD framework and using HotQCD lattices. Attempts to study charmonia at nonzero temperature using NRQCD have been presented in Ref. [32]. In that work the NRQCD parameter was fixed to be M c = 1.275 GeV for all values of β. Our analysis shows that this choice of the mass parameter results in η c mass above 4 GeV for β > 6.74.
Having determined M b we can study the spectrum of bottomonium states at T = 0. In NRQCD absolute values of the meson masses cannot be calculated, only the corresponding energy levels that are related to the masses by a lattice spacing dependent constant can be determined. Therefore, we will consider the differences between the energy levels of various quarkonium states and the energy levels of η b . These are equal to the corresponding differences in meson masses. In what follows, we will use the energy of η b states to calibrate the energy scale at different lattice spacings, i.e., we will set the energy of η b state to be zero. With this in mind we will use the term mass and energy level interchangeably. In Table III we show the difference of the masses of various bottomonium states and the η b mass. These differences are compared to the experimental values from PDG [37]. For P-wave bottomonia, namely for χ b0 , χ b1 , χ b1 and h b the mass difference agrees well with the experimental results. The difference between the Υ and the η b mass, i.e., the hyper-fine splitting is smaller than the experimental value. This is similar to the findings of Ref. [20], where order v 6 NRQCD action was used. One needs the radiative corrections in the NRQCD Lagrangian in order to reproduce the hyper-fine splitting [22].
Before studying bottomonium properties at non-zero temperature we would like to understand the spectrum of energy levels encoded in the extended meson operators used in this study. On general grounds, we expect that the spectrum of energy levels consist of a ground state, one or two excited states below the open bottom threshold, and many higher lying states that in the infinite volume limit form a continuum. In the case of finite volume, the spectral function is always given by the sum of δ-functions. However, even in this case the density of states at large ω is very large [32] and, in practice, it is possible to approximate the spectral function by a continuum. Thus, the spectral function can be approximated as Here, s 0 is the open beauty threshold. For point sources with n = 1/2 for S-wave quarkonia and n = 3/2 for P-wave quarkonia. For large values of ω ∼ 1/a the spectral function is distorted by lattice artifacts and vanishes above some ω max [32]. For extended sources the form of ρ cont is not known, but the above general form of the spectral function is still valid. Since only a few lattice data points in the correlation function are sensitive to the high energy part of the spectral function, and excited state contribution is suppressed when extended sources are used, we can employ a simplified parameterization of the spectral function i.e., we can consider only the ground state contribution and some continuum contribution, which has support only for ω > M . The explicit form of ρ high is not important for our analysis. If the above equation holds, the correlation function can be written as As discussed in the previous section, the correlators of extended operators are dominated by the ground state for τ > 0.5 fm and, therefore, it is possible to determine the parameters A and M from the single exponential fit in this τ region. From this we can obtain C high (τ ). In Appendix B we show the determination of C high (τ ) for different lattice spacing, which will be used in the following sections to analyze the correlators at non-zero temperature. It turns out that the results for C high (τ ) can be well fitted using a simple step function for ρ high (ω) at some threshold ω = ω 0 .

IV. FINITE TEMPERATURE RESULTS
The temperature range which we explore in this paper goes from 151 MeV up to 334 MeV. As discussed in the previous section, the temperature range is limited by the fact that the NRQCD approximation starts to break down for aM b < 1. The last temperature we look at is right on the boundary of this limit. We have explored all temperatures with both point to point correlators, and smeared to smeared correlators. The analyses of the T > 0 correlators are performed in terms of the effective mass, M ef f , calibrated with respect to the T = 0 mass of η b meson.

A. Correlators with point sources at T > 0
We start our discussion with the results obtained with point sources. In Fig. 3 and Fig. 4 we show the effective masses corresponding to η b and χ b0 correlators, respectively. For the η b correlator we see very little temperature dependence. In fact, no medium effect can be seen at the lowest temperature. Somewhat larger temperature dependence is observed for the χ b0 correlators. We see that the effective masses become smaller compared to the T = 0 effective masses with increasing temperature and τ . This implies enhancement of the finite temperature correlators compared to the corresponding one at zero temperature with increasing temperature and Euclidean time, which is consistent with the previous studies of correlators with point sources [30][31][32].
At sufficiently high temperatures we expect that bottomonium states will be dissolved and the correlator, as well as the spectral function will be reasonably well described by the free theory. This will happen at smaller temperature for the P-wave bottomonia because of their larger size. In absence of lattice artifacts the free correlator of the P-states will be proportional to τ 5/2 , c.f. Eq. (12). Therefore, the effective masses will have the form M ef f = b + 5/(2τ ). In Fig. 4 we compare this expectations with the lattice results for the χ b0 effective masses at the highest temperature, T = 333.5 MeV. The constant was fixed to match with the lattice results at the largest available τ . The free-theory form seems to capture the behavior of the lattice results except at  small τ . This is in agreement with the conclusion reached in Ref. [26], and suggest that the lattice results are not incompatible with the dissolution of P-state quarkonia. The free-theory form of M ef f does not describe the lattice results for τ < 0.4 at quantitative level. This is not completely surprising, since the spectral function at large ω that dominates the small τ behavior of the correlator is strongly influenced by cutoff effects.
From the results shown in Figs. 4 and 3 we see that the effective masses corresponding to correlators with point sources are dominated by the large ω part of the spectral function and, thus, show little temperature dependence. Similar conclusion has been reached in Ref. [32]. Therefore, in the following we will study in-medium bottomonium properties using smeared Gaussian sources and sinks.

B. T = 0 with smeared sources
We study the effective masses for the smeared correlators to get some insight about in-medium modification of the bottomonium states. Our results for the effective masses of Υ and χ b0 at two representative temper- atures are shown in Fig. 5. In the figure we also show the corresponding zero temperature results for the reference. At small τ we see little to no temperature dependence in the effective masses. For the lower temperature, T = 199 MeV, we see an approximate plateau in the Υ effective masses for τ < 0.8 fm and a rapid drop at larger τ . For the highest temperature the Υ effective masses show a mild approximately linear decrease in the τ region, where the corresponding T = 0 effective mass has an approximate plateau. At large Euclidean time, τ > 0.45 fm we see again a rapid drop in the effective mass of the Υ correlator. The effective mass of χ b0 correlator at T = 199 MeV shows a behavior that is very similar to that of the Υ effective mass at the highest temperature. This is expected, since χ b0 state being larger is more affected by the medium. The effective mass of χ b0 at the highest temperature does not show any remnant of a plateau, but a significant decrease with a slope that is increasing with increasing τ .
A similar behavior in the effective masses to the one observed above has also been found in the calculations of static meson correlation functions at finite temperature, for sufficiently large separation between the static  quark Q and static anti-quarkQ [38,39]. In the case of static mesons, the temperature and τ dependence of the effective masses can be more easily understood. For sufficiently high temperatures, the spectral function of a static meson can be calculated in Hard Thermal Loop (HTL) re-summed perturbation theory [40]. Perturbative calculations show that the energy of static QQ pair is complex [33,41]. This energy is also often referred to as the complex QQ potential. Therefore, the spectral function of a static meson has a peak at small ω. The position of the peak is given by the real part of the potential, while the imaginary part of the potential determines the width of the peak. We also expect that there is a continuum part of the spectral function at large ω [42]. Around the peak position, the shape of the spectral function is well described by a Lorentzian (or skewed Lorentzian) [40]. However, the peak has also long tail at small ω [40], which is not related to the imaginary part of the potential [43]. This tail determines the large τ behavior of the correlation function of static QQ meson. The qualitative features of the spectral function of static QQ meson obtained in HTL re-summed perturbation theory can also help to explain the behavior of the effective masses of static meson correlators observed in lattice calculations [38,39]. In particular, it can be shown that the large τ behavior of the effective masses is determined by tails of the spectral function at small ω. Given the above discussion it is reasonable to assume that the spectral function at T > 0 has a broadened peak around ω that corresponds to a quarkonium state, and a high frequency part that is identical to the zero temperature one, i.e., we can assume that the spectral function has the form ρ(ω, T ) = ρ med (ω, T ) + ρ high (ω) .
Here, ρ med (ω, T ) describes the spectral functions for ω M and/or ω < M . A natural parameterization of ρ med (ω, T ) would be a Breit-Wigner (Lorentzian) form with M and Γ being the temperature dependent mass and width of the bottomonium state, respectively. The Lorentzian form is expected to capture well the main features of the spectral function for ω ∼ M , but also has a long tail for ω M , where it is not adequate. As discussed above in the case of static mesons, the Lorentzian only works in the vicinity of the peak and the functional dependence on ω is very different from the Lorentzian for ω well below the peak position. Since we do not know  the functional form of the spectral function at very small ω, we assume that The parameter ω cut is an effective way to parameterize the tail of the spectral function at small ω. It turns out that the Euclidean correlation function and the effective masses are sensitive to the choice of ω cut at large τ . Thus, we need at least three parameters to describe the medium dependent part of the spectral function, ρ med : the peak position, M , the thermal width Γ, and the ω cut that parameterize the tail of the spectral function at small ω.
The above discussion holds in the infinite volume limit. In lattice calculations the volume is finite and, therefore, the number of available energy levels is also finite. Therefore, ρ med (ω, T ) should be given by a sum of delta functions. Furthermore, for typical volumes used in present day lattice calculations the number of energy levels is not very large for ω ∼ M , see discussions in Ref. [32]. In particular, we should not expect that there are many delta functions in the lattice volume that will effectively parameterize the low ω tail of the spectral function.
At small τ the correlators and the effective masses are mostly sensitive to the high energy part of the spectral function. Since the high energy part of the spectral function is temperature independent, the effective masses at small Euclidean time are not sensitive to the effects of the medium on quarkonium states. One should consider only large τ values to gain sensitivity to in-medium bottomonium properties. On the other hand, as discussed above, at large τ the behavior of the correlators is sensitive to the small ω tail of the spectral function, which is also unrelated to bottomonium properties at T > 0. Therefore, we consider the subtracted correlator If Eq. (15) is valid, the subtracted correlator at small τ should be sensitive to ρ med (ω, T ) at ω M .
We can also define an effective mass corresponding to the subtracted correlator M ef f,sub (τ, T ), which is shown in Fig. 6 for η b and χ b0 . There are three key features of M ef f,sub (τ, T ): at small τ it is close to T = 0 bottomonium mass, at intermediate τ we see a mild, approximately linear decrease, and finally, there is a sharp drop at large τ values. The fact that M ef f,sub (τ, T ) is close to the T = 0 mass at small τ may suggest that the quarkonium mass is not significantly shifted relative to its vacuum value. The mild linear decrease in the difference of the effective masses is closely related to the width of the bottomonium state. To extract in-medium bottomonium properties we fit the data on the subtracted effective masses using the following Ansatz for ρ med ρ med (ω, T ) = A cut δ(ω − ω cut ) The parameters A cut and ω cut effectively describe the tail of the spectral function at low ω. If A cut is small, the subtracted effective mass corresponding to the above equation has the form M ef f,sub = M (T ) − Γ 2 τ , i.e., the Gaussian form naturally explains the approximate linear dependence of the effective masses. At large τ the tail of the spectral function at small ω also becomes important, and we see deviations from the linear behavior. As can be seen from Fig. 6, fits using Eq. (19) welldescribe the lattice results on M ef f,sub (τ, T ) . The χ 2 /dof lies in the range of 0.2 to 0.7 for η b and Υ, and 0.1 to 0.6 for χ b . From the fits we obtain the in-medium bottomonium mass M (T ) and the width parameter Γ(T ). To test the robustness of the fit procedure, we also performed fits omitting two data corresponding to the two largest τ values, and by setting A cut = 0. Such fits work well for the effective masses of η b and Υ at all temperatures as well as for χ b0 and χ b1 except at the highest temperature. We began the fits from τ = a or 2a, and used the average of the two fit results; the difference of these two fit results were chosen as the systematic errors. This  systematic error was added in quadrature to the statistical error calculated using jackknife sampling. Our results for the medium masses are shown in Fig. 7 in terms of the mass differences (shifts) ∆M (T ) = M (T ) − M (T = 0) for 1S and 1P bottomonium states. The mass shift is compatible with zero when statistical and systematic uncertainties are taken into account. The in-medium width parameter for different bottomonium states is shown in Fig. 8. The fits generally work well, but some systematic errors does appear, as seen in the right plot for χ b at T = 199 MeV. At this temperature the linear behavior in the effective mass generated by the Gaussian form can equally be well produced by the negative peak, which was introduced to explain the strong drop off at τ = 1/T . For the same reason in the fits for η b and Υ at T ≤ 251 MeV and for χ b at T ≤ 173 MeV, the negative peak has not been included in the fit, since the data points appear to follow a linear behavior to a very accuracy. The fit can not see the difference between a fit with a Gaussian or a fit with two delta functions. We see from the figure that at the lowest temperature the width is compatible with zero for both 1S and 1P bottomonia. The width of η b and Υ is small for T = 173'MeV, while it becomes significant for χ b0 and χ b1 . In general, the width is larger for 1P bottomonia than for 1S bottomonia. This is expected since 1P bottomonia being larger are more affected by the deconfined medium. The extracted values of the width parameter are, in most cases, not sensitive to the simplified description of the low ω tail, since the results from the fits using all data points agree mostly with the results from the fits without the last two data points and A cut = 0.
It is important to ask the question to what extent the above results for the in-medium mass shift and width depend on the interpolating operator used in the analysis. The continuum part of the spectral function clearly depends on the choice of the interpolating operator. However, this part is, to a large extent, temperature independent and does not show up in the subtracted effective masses. If the in-medium bottomonium mass and width are physical then they should not depend on the choice of the interpolating operator. On the other hand, the low energy tail of the spectral function may depend on the choice of the interpolating operator. We discuss this issue in Appendix A, where we study the bottomonium correlators corresponding to interpolating operators of different sizes. As shown in this Appendix, the in-medium bottomonium mass and width only mildly depends on the size of the interpolating operators, and this dependence is likely due to the imperfection of the fit Ansatz of the spectral function.
The absence of a significant bottomonium mass shift at finite temperature should not come as a complete surprise. The energy levels that enter the spectral decomposition are the same as in the vacuum, i.e., are not effected by the temperature, and include the vacuum bottomonium state. An effective in-medium shift of the peak position of the spectral function can come from additional energy levels corresponding to multi-particle states containing quarkonium. The density and distribution of these energy levels could be such that they correspond to a peak of some width, with a maximal density at an ω value that is different from the vacuum quarkonium mass. We do not see any indication of such shift in our analysis. One may wonder if this is due to the small volume used in our calculations, since the density of the additional states could be too small and the location of the corresponding energy levels could be highly distorted. Therefore, for the highest temperature T = 334 MeV we performed calculations on a smaller 36 3 × 12 lattice. The results are discussed in Appendix C. As shown in this Appendix, we do not see any significant volume effects in the bottomonium correlators, thus, it is not clear if the absence of in-medium mass shift is due to finite volume effects. Instead, it is possible that the absence of an effective mass shift is due to the use of a simple Ansatz for ρ med .
For finite volume the most natural representation of the in-medium spectral function is the sum of δ-functions. If we want to represent ρ med (ω, T ) by a sum of N δfunctions, we have to fit 2N parameters and, thus, N cannot be too large. We have seen from the previous analysis that two to three parameters are sufficient to describe the τ -dependence of the subtracted effective masses, which show an approximate linear behavior at small and moderate τ and a fast drop-off at large τ that corresponds to the low ω part of ρ med (ω, T ). We have also seen that bottomonium masses are not modified with respect to their vacuum values. The simplest representation of ρ med (ω, T ) in terms of δ-functions consistent with the above features is where M 0 is the vacuum bottomonium mass, and the parameters A cut and ω cut describe the low ω tail of the spectral function. It is easy to see that for small δM the above Ansatz gives an effective mass that decreases linearly in τ , with the slope equal to δM 2/3 if the contribution proportional to A cut can be neglected. Therefore, δM 2/3 can be interpreted as a width parameter. We performed fits of the lattice results using Eq. (20), and obtained the values δM , A cut and ω cut . Our results for δM 2/3 are shown in Fig. 9. We see from the figure that the width parameters obtained from the fits with Eq. (20) agree well with the estimate of the bottomonium states obtained using the Gaussian form of ρ med . This suggest that our estimates of the width are robust.

V. CONCLUSION
We studied S-and P-wave bottomonia correlators at non-zero temperature using extended Gaussian, as well as point meson sources within the framework lattice NRQCD including v 6 spin-dependent terms. Correlators of point sources show little temperature dependence because of their limited sensitivity to bottomonia properties at small Euclidean time. Therefore, we focused on analyses of the correlators with extended Gaussian meson sources. We identified the contributions of the high energy part of the spectral functions to the correlators of extended Gaussian sources at T = 0. This contribution was then subtracted from the correlators at finite temperature, significantly simplifying the analyses. The τ -dependence of the subtracted correlators with extended Gaussian sources could be understood in terms of a simple theoretically-motivated spectral functions, consisting of a single broadened peak. We identified two prominent structures in the corresponding effective masses: an approximately linear decrease at small and intermediate values of τ , which is related to the width of the peak, and a rapid drop in the spectral function around τ 1/T , corresponding to the small ω tail of the peak. At present statistical accuracy these features are very well described just by three parameters.
We estimated the thermal width of bottomonium properties using two different forms of the spectral functions that capture the above behavior of the effective masses. We found that the thermal broadening of P-states is larger than the thermal broadening of the S-states, and sets in at lower temperatures. This is expected based on the difference in size of S-state and P-states. The inmedium bottomonium masses, defined as peak positions, did not change relative to their vacuum values. This may appear somewhat unexpected, especially for the P-states. We did not find any indications that the absence of mass shift being due to finite volume of the lattice. At current level of statistical accuracy the lattice results of the correlators cannot resolve further details on the shape of the spectral function beyond its overall width. Therefore, to resolve the detailed shape of the spectral function, including its asymmetric nature, and possible shift in the peak position, more precise lattice results are needed. It is possible that with increase statistical precision one will see more sensitivity to finite volume effects too.

Appendix A: Smearing Dependence
In the main analysis we use Gaussian sources of size 0.21 fm, which we will refer here to as the standard sources size. We studied to what extent our conclusions about in-medium bottomonium properties depend on the choice of the source size. Therefore, we calculated η b correlators at T = 334 MeV using source sizes 0.15 fm, 0.21 fm, 0.25 fm and 0.41 fm, which are 71%, 122%, 141% and 200% of the standard source size used in the main analysis, respectively. The corresponding numerical results for the η b effective masses are shown in Fig. 10. The qualitative behavior of the effective masses obtained for different smearing sizes is similar, except for the largest sources size (200% of the standard source size), where no plateau like structure can been seen. At small τ the effective mass is the smallest for the standard source size implying that the standard sources size is close to the optimal one.
We also calculated the P-wave correlators at T = 199 MeV for different source sizes, namely for sources sizes that are 50%, 141% and 200% of the standards source size. The results for the effective masses are shown in Fig. 11. We see that the source size of 0.1 fm is too small, the corresponding effective mass is very large and does not approach a plateau. We also see no clear plateau for the largest source size corresponding to 200% of the standard source size. The effective masses obtained using source size of 0.21 fm and 0.29 fm look similar, but the effective masses are smaller at small τ for the source size of 0.29 fm. Therefore, for P-waves the optimal source size is around 0.29 fm.
To check the dependence of in-medium bottomonia properties on the smearing size, we performed fits on the effective masses using the following Ansatz for the spectral function ρ(ω, T ) = A cut δ(ω − ω cut ) The last term in this equation parameterize the continuum part of the spectral function. We performed fits on the effective masses obtained with source sizes of 71% to 141% of the standard source size and extracted the in-medium masses and width of the bottomonium states. The in-medium parameters obtained from the effective masses with different source sizes are similar but not exactly the same. This is due to the fact that here we fit the continuum part of the spectral function, which is different for different source sizes. The continuum parameters are correlated with the mass and width parameter and thus the latter are also affected. It is possible that the simple fit form given by Eq. (A1) is less appropriate for source sizes that are not around the optimal value. In order to have a more robust comparison of the results obtained with different source sizes, for the source size of 0.29 fm we also performed calculations of the Pwave correlator at zero temperature. This enabled us to calculate the high energy part of the correlator, C high (τ ), and subtract it from the corresponding T > 0 result. The resulting subtracted effective mass is shown in Fig. 12, and compared to the result obtained with source size of 0.21 fm discussed in the main text. The subtracted effective masses for these two source sizes appear to be quite similar, though the effective mass corresponding to source size of 0.29 fm shows a faster fall-off at large τ . This means that the low ω tail is more important for the larger source size. We performed fits on the subtracted effective masses for the above two source sizes using Eq. (19) for the spectral function, and obtained the in-medium mass and width. We find that the inmedium mass agrees with the vacuum mass also for the The width parameter is 36 ± 59 MeV larger for the larger source size. This is likely due to the correlation between the Gaussian part of the spectral function and the low ω tail.  similar to the ones obtained on 48 3 × 12 lattice. In fact, the largest deviation we have seen between the effective masses obtained on two different lattice volumes was 1.5σ. As an example in Fig. 14 we show the difference of the effective masses calculated on 36 3 × 12 lattice and 48 3 ×12 lattice as function of τ for Υ and χ b0 . At present statistical accuracy we do not see volume dependence of quarkonium correlators at the highest temperature.