The bottomonium spectrum at finite temperature from $N_f=2+1$ lattice QCD

We present results on the bottomonium spectrum at temperatures above and below the deconfinement crossover temperature, $T_c$, from dynamical lattice QCD simulations. The heavy quark is treated with a non-relativistic effective field theory on the lattice and serves as a probe of the hot medium. Ensembles with a finer spatial lattice spacing and a greater range of temperatures below $T_c$ than those previously employed by this collaboration are used. In addition, there are $N_f=2+1$ flavours of Wilson clover quark in the sea with $M_\pi\approx400$ MeV and we perform a more careful tuning of the bottom quark mass in this work. We calculate the spectral functions of S and P wave bottomonium states using the maximum entropy method and confirm earlier findings on the survival of the ground state S wave states up to at least $2T_c$ and the immediate dissociation of the P wave states above $T_c$.


Introduction
The dissociation of heavy quarkonia in a deconfined medium may provide a valuable thermometer for relativistic heavy-ion collisions [1]. Dissociation, or melting, contributes to the suppression of the yield of these states in nuclear collisions compared with hadronic ones [2]. Suppression patterns [3,4] are complicated by the statistical recombination of heavy quarks in the quark-gluon plasma (QGP). However, competing effects are expected to be less pronounced for bottomonium than for charmonium [5]. Indeed sequential suppression has been observed recently in the Υ system by CMS at the LHC [6]. It is therefore desirable to understand the dissociation of mesons in the bottomonium system from ab initio QCD. Analytic results from effective field theories [7][8][9][10][11][12] can aid the identification and interpretation of the relevant physical phenomena, such as the role of the finite width of such states in the plasma as well as the familiar colour-Debye screening mechanism [9]. However, in order to achieve a separation of scales, they generally rely on a choice of the parametric size of the temperature, for instance T ∼ gm b or T ∼ g 2 m b , where m b is the heavy quark mass and g is the coupling which is assumed to be sufficiently weak for a hierarchy of scales to emerge. In contrast, by directly simulating the non-relativistic field theory for the heavy quark (NRQCD) non-perturbatively, all that is required is that m b T , which is certainly satisfied for the bottom quark in the relevant regime up to 4 or 5T c . Numerical lattice simulations are well suited to investigate the properties of the strongly-coupled QGP and to capture the non-perturbative dynamics of the hot medium formed around the crossover temperature, T c . However, it remains a challenge to control systematic errors in dynamical lattice simulations. This paper continues the work of the fastsum collaboration's programme to investigate phenomenologically relevant observables in the QGP with improved control over uncertainties [13][14][15][16][17]. Here, we calculate bottomonium spectral functions from N f = 2 + 1 ensembles using the maximum entropy method. The determination of the spectral functions from the lattice aims to complement other approaches such as weakcoupling effective field theory and potential models to give a complete understanding of the behaviour of heavy quarkonium through the deconfinement crossover. Previous studies of quarkonium spectral functions from lattice QCD, mostly for charmonium, include [18][19][20][21][22][23]. An earlier analysis of the spectral functions by this collaboration from N f = 2 flavour ensembles suggested the ground state S waves (η b and Υ channels) survive up to at least 2T c while the first excited state may be suppressed in the deconfined phase, both at zero [14] and non-zero momentum [15]. Those conclusions are consistent with predictions from effective field theory [9] and experimental data [6], although a detailed comparison with the latter is beyond the scope of this paper. An examination of the correlation functions indicated [13] that the P waves (h b , χ b1,2,3 channels) are greatly modified directly above T c . By comparison with the expectations from the free continuum effective theory the observed thermal modification of the correlators provided evidence in favour of the hypothesis that the P waves dissociate above the crossover temperature. This interpretation is supported by the examination of the spectral functions in those channels [16]. Here, we extend that work by using new N f = 2 + 1 ensembles with improved lattice parameters. The results for the spectral functions are compatible with earlier findings but in the analysis of the correlators we note that the finite threshold plays a more prominent role, discussed in section 3.1. The following section outlines the simulation details and zero temperature calibration. We present the bottomonium correlators and spectral functions at finite temperature in section 3, discuss systematic effects in the spectral function reconstruction in section 4 and conclude in section 5.

Lattice set-up
In this work we employ ensembles with anisotropic lattice spacings, using a tadpole-improved Wilson clover quark action for the light and strange quarks and a tadpole-and Symanzikimproved gauge action. Tree-level improvement coefficients are used for both the fermion and gauge actions. The tuning of the lattice parameters was performed by the Hadron Spectrum Collaboration and further discussion can be found in ref. [24]. A range of temperatures above and below the deconfinement crossover is accessible from 0.76T c to 1.90T c . The fixed-scale approach is used whereby the temperature is varied by changing the number of temporal lattice sites while the lattice spacings are kept constant. This reduces the overhead of zero temperature simulations required to tune the lattice parameters. The renormalized anisotropy, ξ ≡ a s /a τ , is tuned to ξ = 3.5, which allows us to maintain an adequate resolution of the correlation functions to obtain the spectral functions at high temperatures. These "second generation" ensembles [17] represent multiple improvements over the "first generation" ones used in the previous work by this collaboration [13][14][15][16]. In particular, the spatial lattice spacing, a s = 0.1227 (8) Table 1. Summary of the ensembles used in this work. The crossover temperature is determined from the renormalized Polyakov loop [25]. The zero temperature tuning of the lattice parameters was completed by the Hadron Spectrum Collaboration [24].  is reduced. Moreover, the strange quark is now included in the sea. Details are given in tables 1 and 2. On these ensembles the pion is relatively heavy, M π ≈ 400 MeV, while the kaon is roughly physical, M K ≈ 500 MeV [26]. The highest accessible temperature is slightly reduced in this study with respect to the earlier work.

Lattice NRQCD
NRQCD is an effective field theory with power counting in the heavy quark velocity in the bottomonium rest frame, v ∼ |p|/m b [29]. The finite lattice spacing cuts off the relativistic modes of the heavy quark in the discretized theory. The heavy quark and anti-quark fields decouple and their numbers are separately conserved. Their propagators, S(x), solve an initial-value problem whose discretization leads to the following choice for the evolution equation with U τ (x) being the temporal gauge link at site x and e τ the temporal unit vector. The leading-order Hamiltonian is defined by The higher order covariant finite differences are written in terms of the components of the usual forward (∇ + i ) and backward (∇ − i ) first order ones. The correction terms are given by parameter is sufficient for these anisotropic lattices in order to satisfy the stability criterion max |1 − a τ H 0 /2k| < 1. Other choices of k = 2, 3 were investigated but their effects are negligible. The clover definition of the field-strength tensor is used to define the unimproved chromoelectric and magnetic fields and ∇ ± is the symmetric covariant finite difference operator. Tadpole improvement is implicit and implemented by dividing all links by the mean link determined from the fourth root of the average plaquette. Only energy differences are physically significant in NRQCD because the rest-mass energy can be removed from the heavy quark dispersion relation by performing a field transformation. The determination of the lattice heavy quark mass, a s m b , is achieved by tuning a hadronic kinetic mass, M 2 , at zero temperature which appears in the hadronic lattice dispersion relation where n i = −N s /2 + 1, . . . , N s /2. In this work we tune the heavy quark mass by requiring the spin-averaged 1S kinetic mass, M 2 (1S) = (M 2 (η b ) + 3M 2 (Υ))/4, to be equal to its experimental value. Using the spin-averaged kinetic mass mitigates the systematic error due to the determination of the hyperfine splitting in the kinetic mass [30]. The fits to the dispersion relations, shown in figure 1, at the tuned value of the heavy quark mass a s m b = 2.92, are given by with a statistical error determined from a bootstrap analysis. Higher order terms in the dispersion relation could not be resolved within the statistical precision. The tuned value of the heavy quark mass corresponds to M 2 (1S) = 9560(110) MeV which is consistent with the experimental value, M expt (1S) = 9444.7(8) MeV. The error includes a contribution from the statistical uncertainty in the scale set from the Ω baryon mass [31]. This tuning is an improvement over the ad hoc choice of the heavy quark mass in the previous study. Figure 2 shows correlation functions in the Υ (red crosses) and χ b1 (blue circles) channels and the energies determined from single exponential fits. The experimental Υ(1S) mass is used to fix the absolute energy shift, E 0 , of spectral quantities

Zero temperature results
The resulting zero temperature spectrum is given in table 3. We note that in ref. [32] a heavy quark action with O(v 6 ) corrections was used in conjunction with the same zero temperature ensemble used in this work. Slight discrepancies between the spectra could arise from a different heavy quark action including a different choice for the matching coefficients. Precision studies of the zero temperature bottomonium spectrum using NRQCD have been performed by the HPQCD collaboration [30,33]. The maximum entropy method (MEM) with Bryan's algorithm [35] was used to obtain the spectral functions, ρ(ω), from the hadronic correlation functions, which are related through Although this is a typical ill-posed problem given the discrete and noisy estimator for the correlation function, G(τ ), MEM gives a unique solution after specification of the default Correlation functions in the Υ and χ b1 channels with the corresponding ground state energies extracted from exponential fits. Filled symbols denote negative data excluded from the fit.  Table 3. Bottomonium spectrum from standard exponential fits where the Υ mass has been used to set the absolute energy shift, E 0 . The error quoted in the third column is statistical, while the error in the fourth column includes a contribution from the statistical uncertainty in the scale [26,31]. The experimental masses are taken from the Particle Data Group booklet [34].
model. Further details on the implementation used can be found in ref. [14], where also the default model dependence is studied in detail. Figure 3 shows the spectral functions in the S wave Υ (left) and P wave χ b1 (right) channels along with energies extracted from multi-exponential fits directly to the correlators using the CorrFitter package [36] shown with vertical black dotted lines. Good agreement is observed between the energies extracted directly from the correlators and the peak positions in the reconstructed spectral functions. We note that while as many as six well-defined peaks can be discerned in the S wave channel presented here, the third and higher peaks are not compatible with the experimental spectrum [34]. A more thorough investigation of the zero temperature spectrum could be undertaken following the HPQCD approach [30,33].

Bottomonium at finite temperature
The heavy quarks are valence quarks which propagate through the thermal medium according to the solution of the initial value problem in NRQCD. Their propagators do not satisfy anti-periodic thermal boundary conditions so the heavy quarks are not in thermal equilibrium with the medium. This is illustrated in the representation of the correlation function in eq. (2.7) which is manifestly not symmetric about τ = 1/2T , while the kernel, K(τ, ω), is independent of the temperature. The thermal modification of the correlators can therefore be directly attributed to the modification of the associated spectral function. The asymmetry of the hadronic correlation functions can be seen explicitly in figure 2. These simplifications result from replacing ω → 2m b + ω and taking the m b /T → ∞ limit in the standard kernel at finite temperature where n B (ω) = 1/(exp(ω/T )−1) is the Bose-Einstein distribution. The inversion of eq. (2.7) to obtain the spectral function is simpler with this asymmetric kernel, since the correlation functions are accessible to much greater temporal separations than for relativistic quarks. Moreover, unlike a formulation in which the quarks are in equilibrium with the medium there is no constant contribution to the hadronic correlation functions. This reflects the fact that NRQCD is an effective field theory around the two-quark threshold.

Thermal modification of correlation functions
The ratios of the correlation functions at finite temperature to those at zero temperature are shown in figure 4. With increasing temperature the correlators are enhanced at earlier temporal separations. We observe temperature dependence already below T c and stronger dependence above T c . An enhancement of approximately 4% is seen in the S wave Υ  channel (left) above T c but there is a greater effect of almost 20% in the P wave χ b1 channel (right). The enhancement relative to the low temperature correlators in the N f = 2 studies [14,16] was approximately 2% and 25% in the S wave and P wave channels, respectively. The correlators in the other S wave (η b ) and P wave (h b , χ b0 , χ b2 ) channels show analogous modifications to the Υ and χ b1 channels. The S wave effective mass displays little temperature dependence (figure 5, left) but a clear effect is seen in the P wave channel effective mass (right). In ref. [13] it was also observed that the S wave effective mass showed little variation with temperature while the temperature dependence in the P wave channel effective mass was even more pronounced than visible here. It is useful to compare the effective mass at high temperatures with that in the non-interacting infinite temperature limit. In continuum NRQCD the spectral functions are known for free heavy quarks [9], and are given by

3/2, P wave. (3.2)
To facilitate the comparison with the interacting effective theory we have included a threshold, ω 0 , to account for the additive shift in the quarkonium energies 1 . This threshold may depend on the lattice parameters and the temperature. In the infinite temperature limit the correlation functions then have the following behaviour and the effective mass becomes Pure power law decay in the correlator or the absence of plateaus in the effective mass may be less evident in the presence of a larger threshold. In the earlier N f = 2 studies comparisons between the correlators and effective mass and their non-interacting counterparts  were possible without including a finite threshold. However, the threshold appears to play a more significant role here which is reflected in the fact that there is a change in the energy shift, E 0 , of about 300 MeV between the two studies, see table 2.
The thermal modification of the spectrum is expected to depend on the heavy quark mass. In simple potential models the binding radius is related to the typical inverse momentum transfer, r H ∼ (m b v) −1 , and effective colour Debye screening occurs when the screening length, r D , is comparable with or smaller than the binding radius, r D r H . At a given temperature, we would therefore expect such a mechanism to become less effective for heavier, more tightly bound states. Figure 6 shows the modification of correlators for various lattice heavy quark masses and illustrates the mass dependence in the Υ (left) and χ b1 (right) channels. Correlators in the χ b1 channel exhibit greater thermal modification than in the Υ channel at each of the lattice heavy quark masses investigated. At smaller values of the heavy quark mass, approaching the charm quark mass, a large enhancement is seen even in the Υ channel correlation function, while for large values of the heavy quark mass some enhancements are still seen in the χ b1 channel. The mass dependence has also

Spectral functions
Figures 7 and 8 depict the spectral functions in the Υ and χ b1 channels respectively at temperatures from 0.76T c up to 1.90T c . For clarity each panel displays just two neighbouring temperatures. In the Υ channel the ground state peak is clearly visible and coincides with the energy extracted from the exponential fit to the correlation function at zero temperature, see figure 3 (left). The ground state peak persists at all accessible temperatures demonstrating the survival of the ground state to at least T = 1.90T c . We observe a broadening in the peak and a decrease in its height above T c . Below T c the second peak may be identified with the first excited state. Its interpretation above T c is less clear which may be due to melting as well as the possible presence of lattice artefacts in the high frequency part of the spectral function, which are discussed further in appendix A. In the χ b1 channel the ground state peak can be discerned at temperatures below the T c and agrees with the energy from the exponential fit at zero temperature. This peak is observed to disappear immediately in the deconfined phase which we suggest indicates the dissociation of this state almost as soon as the deconfined phase is reached. We note that the ground state peak in the P wave channels is harder to distinguish than in the S wave channels, even below T c .

Systematic tests of MEM
A close examination of the reconstruction of the spectral functions is essential to have confidence in the interpretation of temperature effects. Here we address some pertinent issues due to the selection of the temporal range of the correlator and the frequency domain of the spectral function used in the MEM. Other effects such as the dependence on the default model and the statistical uncertainty have been investigated for similar data from the N f = 2 ensembles [14,16] where they were noted to be mild. Typical systematic effects in lattice studies such as the unphysical pion mass and finite-volume effects are not discussed although the latter are expected to be small for such heavy quarkonium states. The stability of the spectral function with the variation of the temporal range of the correlation functions used in the reconstruction is shown in figure 9. The spectral functions are stable as long as data at temporal separations close to τ /a τ = 0 or N τ are excluded on account of lattice artefacts. Effects due to the inclusion of temporal separations near N τ have also been discussed in ref. [16]. Although there are no temporal boundary conditions for the heavy quark fields, we recall that the gauge fields are periodic. Since the spatial lattice spacing is coarser than the temporal one, a s = 3.5a τ , effects at separations close to N τ may be expected at this scale. The effect is stronger in the P wave channels [16].
In the Υ channel (figure 9, top panels) the spectral function is stable when varying the temporal window as long as the correlator datum closest to τ /a τ = 0 is omitted. In the χ b1 channel (figure 9, bottom panels) the spectral function is stable as long as the largest temporal separation τ /a τ = N τ − 1 is also excluded. Therefore the reconstructed spectral function converges in all cases when the range of correlator data, [τ 1 /a τ , τ 2 /a τ ], is chosen Υ, N τ = 28 [2,25]   such that τ 1 /a τ 1 and τ 2 /a τ N τ − 1. We have also investigated using a subset of the available correlator data to ensure any inferences about the temperature dependence are due to physical effects and not artefacts of the analysis using fewer data. In figure 10, the reconstructed spectral function in the Υ channel at T = 0.76T c using only half the available correlator data is compared with that from using the entire correlation function. Only a small variation in the ground state peak height is observed. The frequency domain chosen for each reconstruction of spectral function is given in table 4. This interval must be chosen judiciously and may extend to negative frequencies as the energies can be shifted by an a priori unknown constant in the effective theory. Furthermore, this range must be sufficiently large to exclude unphysical spectral weight which has been observed to appear at the edges of the interval.

Conclusions
Calculating heavy quarkonium correlation functions using lattice QCD provides valuable input towards understanding the modification of the spectrum in the hadronic and plasma phases of QCD. We observe that with increasing temperature the correlation functions are enhanced relative to low temperature. In NRQCD this temperature dependence arises Υ χ b1 N τ a τ ω min , a τ ω max a τ ω min , a τ ω max  Table 4. Frequency ranges used in reconstruction of the spectral functions. The frequency interval is discretized into N ω = 1000 points for each N τ .
solely from the temperature dependence of the associated spectral function since the integral kernel relating the spectral function and the correlator is temperature independent. These enhancements are greater in the P wave than the S wave channels for each temperature below and above T c . There is significant mass dependence in these modifications, with lighter states showing increased temperature dependence. Further interpretation is aided by calculating the spectral functions using MEM. At zero temperature the reconstructed spectral functions have localised peaks coincident with bound state energies extracted directly from the correlation functions. The analysis of the spectral functions at finite temperature suggests the survival of the ground state Υ up to at least 1.90T c with some modification above T c . The Υ(2S) state appears to dissolve close to T c , although the proximity of lattice artefacts complicates the interpretation of this structure. The ground state χ b1 peak is suppressed immediately above the crossover temperature indicating significant alterations, compatible with the dissociation of this state in the QGP. These results are qualitatively consistent with the conclusions of the earlier N f = 2 studies and a systematic comparison between the temperature dependence of the peak positions and widths is underway. In the future we also plan to investigate the momentum dependence of the bottomonium spectral functions and examine new Bayesian approaches to the reconstruction of the spectral functions [38,39].

Acknowledgments
We thank D. K. Sinclair for code used in the calculation of the heavy quark propagators and A. Rothkopf for stimulating discussions about the Bayesian reconstruction of the spectral functions. We are grateful to the Hadron Spectrum Collaboration for the use of their zero temperature ensemble. This work is undertaken as part of the UKQCD collaboration and the STFC funded DiRAC Facility. We acknowledge the PRACE Grants 2011040469 and Pra05 1129, the Initial Training Network STRONGnet, Science Foundation Ireland grants 11-RFP.

A Free lattice spectral functions
The free lattice spectral functions are calculated by summing over the first Brillouin zone [14] according to A comparison between the free lattice spectral functions for the first and second generation parameters (see table 2) is shown in figure 11 for the S wave (left) and P wave (right) channels. The effect of the finer spatial lattice spacing is apparent as the cusp artefacts -which correspond to momenta reaching the corners of the Brillouin zone -are pushed to higher frequencies, and the support of the free lattice spectral function is correspondingly enlarged. Note that at high temperatures the support of the reconstructed spectral functions (figures 7 and 8) is comparable with that of the free lattice spectral function. The difference between the dotted lines representing the free spectral functions in the continuum is due to the slightly different choice of heavy quark mass between the two ensembles.