Spectral properties of bottomonium at high temperature: a systematic investigation

An important task in QCD phenomenology is to determine the spectral functions of hadronic and current–current correlators. These encode physical information about inter alia the existence and properties of bound states and resonances, transport properties, and multiparticle thresholds. The spectral functions are related to the euclidean correlators computable on the lattice by an integral transform, but computing the spectral function $\rho(\omega)$ given a finite set of noisy correlator data $G(\tau)$ is known to be an ill-posed problem. Numerous methods have been developed to attempt to handle this problem (see [1, 2] for recent reviews), which all have their inherent limitations. The purpose of the study that is reported here is to obtain a better understanding of a subset of these methods by applying them to a specific problem, namely determining the temperature dependence of the mass and width of the 1S and 1P bottomonium states computed in lattice non-relativistic QCD (NRQCD).

The bottomonium system has been studied in lattice NRQCD for a long time [3, 4, 5, 6, 7, 8]. The NRQCD approach has the benefit of a simple spectral representation with a temperature-independent kernel,

Since the heavy quark is not in thermal equilibrium, the full range of points $\tau/a_{\tau}=0,\dots,N_{\tau}-1$ are available for the spectral reconstruction. Finally, since at least the $\Upsilon(1S)$ state is expected to survive as a bound state up to quite high temperatures, the peak position and thermal width remain well-defined quantities which can be used to benchmark the methods.

Earlier results from this project were presented in [9]. Some of the results shown here have previously been presented in [10, 11, 12].

Our simulations are carried out using anisotropic lattices with an ${\cal O}(a^{2})$ improved gauge action and an ${\cal O}(a)$ improved Wilson fermion action with stout smearing, following the parameter tuning and ensembles generated by the Hadron Spectrum Collaboration [13, 14]. The results presented here were produced using the “Gen2L” ensembles, which have $N_{f}=2+1$ active quark flavours with $m_{\pi}\approx 240\,$MeV and an approximately physical strange quark. The spatial lattice spacing is $a_{s}=0.112\,$fm and the anisotropy $\xi=a_{s}/a_{\tau}=3.45$. The temperature is given by $T=(a_{\tau}N_{\tau})^{-1}$ and is varied by changing the number of sites $N_{\tau}$ in the temporal direction. For more details about the ensembles, see Ref. [15, 16] and references therein.

The $b$ quarks were simulated using an NRQCD action including ${\cal O}(v^{4})$ and leading spin-dependent corrections, with mean-field improved tree-level coefficients [4].

The idea behind the moments method is to consider the correlator in (1) as the integral of a weighted spectral function $\rho^{\tau}(\omega)=e^{-\omega\tau}\rho(\omega)$. The temporal derivatives of the correlator are then the moments of $\rho^{\tau}(\omega)$, and the salient features of $\rho^{\tau}(\omega)$, such as the locations and widths of its peaks, and hence those of $\rho(\omega)$, can be inferred from these derivatives.

If $\rho(\omega)$ is approximated by a sum of Gaussians,

where $B_{i}(\tau)$ can be approximated by a constant if $\Delta\Gamma_{i}^{2}\tau\ll\Delta m_{i}$. Specifically, from (5) we see that the ground state width can be extracted in the large-$\tau$ limit assuming the same condition holds, while the first derivative (4) can be used to determine the ground state mass.

The results for the vector ($\Upsilon$) channel are shown in fig. 1. We see that both the first and second log-derivative approach a plateau at large $\tau$. To determine the mass and width, the log-derivatives are fitted to a simplified version of (4,5),

The results for $m_{0}$ and $\Gamma_{0}$ for $\Upsilon$(1S) are shown in fig. 2 as a function of temperature (red points). At first sight, both $m$ and $\Gamma$ appear to increase with temperature. However, this increase may in part be an effect of the shorter temporal range available at higher temperatures, where $M(\tau)$ and $\Gamma^{2}(\tau)$ have not yet reached a plateau at the maximum value of $\tau$. To disentangle this from real, thermal effects, the same analysis has been carried out on the $N_{\tau}=128$ ensemble (which may be considered to be at $T=0$), with the temporal range restricted to $\tau_{\max}=1/T$. The results of this are shown by the blue points in fig. 2. We see that there are no substantial thermal effects for $T\lesssim T_{pc}\approx 170\,$MeV. Above this temperature, the mass from the thermal ensemble lies below that from the corresponding truncated $T=0$ analysis, while the thermal width appears to be increasing from zero at $T\gtrsim 200\,$MeV (albeit still consisten with zero within errors). For further details about the method and results, see [11].

All the linear methods considered here (Tikhonov [17], Backus–Gilbert [18], and HLT [19]) regularise the inverse problem by avoiding a pointwise estimate of $\rho(\omega)$ and instead constructing a “smeared” spectral function,

where the smearing kernel $\bar{\Delta}(\omega,\omega_{0})$ plays the role of a finite-width approximation to $\delta(\omega-\omega_{0})$. It may in turn be expressed as a linear combination of the NRQCD kernel functions,

where $b_{\tau}(\omega_{0})$ are coefficients which encode the location of the sampling point $\omega_{0}$ and are determined by the method. Inserting (8) into (7) and using the spectral definition of the correlator (1) then yields

The coefficients $b_{\tau}$ are determined by minimising a functional $W[g_{\tau}]=A[g_{\tau}]+\lambda B[g_{\tau}]$, where the three methods differ by their choices for these functionals:

and $\Delta_{\sigma}$ is a gaussian with width $\sigma$ centred on $\omega_{0}$. An important distinction between HLT and the other two is that in HLT, the smearing width $\sigma$ is an input parameter, while in the other two methods the smearing kernel $\bar{\Delta}(\omega,\omega_{0})$ and hence its width is entirely determined by the data and the hyperparameter $\lambda$. The stability of the results with respect to variations in $\lambda$ is an important consideration with these methods.

Results for the $\Upsilon$ spectral functions from the three methods are shown in the left panel of fig. 3. We see that the results from all three methods are consistent, but that HLT yields a considerably more well-determined spectral function than the other two, as can be seen from the width of the respective error bands.

In the right panel of fig. 3 we show the spectral functions from the HLT method at all temperatures. We see a progressive broadening and apparent positive mass shift as the temperature is increased. However, in light of the discussion in the previous section, in the absence of an equivalent analysis on the zero-temperature ensemble, we are not yet in a position to draw any firm conclusion on this. For more details about the methods and results, see [10].

Bayesian methods, including the maximum entropy method (MEM) [20] and the BR method [21], regulate the inverse problem by introducing prior information, such that the probability of the spectral function $\rho$ given the data $D$ and the prior information $I$ is given by

where $L$ is the standard likelihood ($\chi^{2}$) and $S[\rho]$ parametrises the prior probability $P(\rho|I)$. $S$ is written in terms of a default model $m(\omega)$ which is the most likely spectral function in the absence of data. The main difference between MEM and BR lies in the form for $S[\rho]$,

The two methods also differ in that in the usual (Bryan) implementation of MEM, $\rho(\omega)$ is restricted to the singular subspace of the kernel ($e^{-\omega\tau}$ in the case of NRQCD), while BR constructs the solution from the full $N_{\omega}$-dimensional space.

In fig. 4 we compare the results from MEM and BR for the $\Upsilon$ spectral function in the vicinity of the $\Upsilon$(1S) peak. With both methods we see that the peak shifts to the left and broadens as the temperature increases. However, the quantitative differences are significant: MEM consistently gives a larger width than BR, and the shift in the peak position also tends to be larger in the MEM results.

Fig. 5 shows spectral functions from the BR method for the $\Upsilon$ and $\chi_{b1}$ at all temperatures. The left hand plots show the results from the thermal ensembles, while the right hand plots show results from the equivalent analysis on the zero-temperature ensemble. We see that the zero-temperature control always gives a positive shift in the peak position and a broadening. It is only when the thermal results differ from the control that a real thermal effect can be concluded. In this case we can infer thermal effects in both channels; in particular, for the $\chi_{b1}$, comparing the thermal and control results we can infer a negative mass shift as well as a thermal broadening, while the thermal results on their own would suggest a positive mass shift.

Fig. 6 shows the results for the thermal mass shift $\Delta M$ and thermal width $\Delta\Gamma$ of the $\Upsilon$ and $\chi_{b1}$ for all the methods discussed here. In addition, we show results from the use of a generalised eigenvalue problem (GEVP) [12], where masses are determined from standard exponential fits to the resulting optimised correlators, while the widths are determined by applying the time-derivative moments methods to the same correlators.

The mass shifts and thermal widths in fig. 6 are obtained by subtracting the zero-temperature mass $M(T_{0})$ and width $\Gamma(T_{0})$,

For the moments and BR methods, $M(T_{0})$ and $\Gamma(T_{0})$ are determined from the truncated zero-temperature correlators (control data) as described in section 3. This prescription attempts to correct for the main finite-$N_{\tau}$ effects, but it should be noted that it may be complicated by the presence of excited states at low $T$ which are not resolved in an analysis of the truncated data and may not be present in the thermal correlators. Further analysis will be necessary to disentange any such effects.

For the other methods (MEM, linear methods, GEVP) the ordinary analysis at the lowest temperature has been used to determine $M(T_{0})$. We expect the masses obtained from the GEVP method to be relatively insensitive to finite-$N_{\tau}$ artefacts, while the GEVP widths come from the moments method applied to the GEVP projected correlators and so may have similar artefacts. We expect the MEM and the linear methods to have a reduced $\Delta M$ and $\Delta\Gamma$ when the full control analysis is carried out. In all cases, uncertainties include both statistical and systematic uncertainties.

The first thing to note is that the linear methods have much larger uncertainties than the others. This is not surprising, as the width of the smearing kernel in these methods is in most cases larger than the physical width of the states under consideration. The Tikhonov and BG results for the $\Upsilon$ mass are too noisy to be shown. In HLT the width is an input parameter and can therefore not be used as an estimator (or upper bound) for the physical width.

With these provisos, we see an encouraging qualitative and, in some cases, even quantitative agreement between the different methods. Notably, there is agreement between BR, MEM, moments and GEVP about a small negative mass shift of up to 40 MeV for the $\Upsilon$(1S), for 150 MeV$\lesssim T\lesssim$250 MeV. A similar mass shift may be seen in the $\chi_{b1}$(1P), using BR, moments and GEVP. The current MEM analysis does not show this shift, but the effect may be obscured by the finite $N_{\tau}$, as observed with the BR method in fig. 5.

Concerning the width, there is rough agreement (within large uncertainties) for the $\chi_{b1}$, but for the $\Upsilon$ there is a discrepancy, with BR and moments giving a much smaller width ($\Delta\Gamma<100\,$MeV at $T=380\,$MeV) than the other methods. Although the full zero-temperature control analysis may change this, the GEVP method should be less sensitive to this, so this discrepancy needs to be studied further.

We have studied the mass and width of ground state S- and P-wave bottomonium ($\Upsilon$ and $\chi_{b1}$) with a range of different methods. We find agreement between correlator and Bayesian methods on a negative $\Upsilon$ mass shift of up to 20–40 MeV at $T\sim 250$ MeV, as well as qualitative agreement on the mass and width of $\chi_{b1}$. There is a discrepancy for the width of $\Upsilon$, which requires further investigation. Linear methods (Tikhonov, Backus–Gilbert, HLT) have intrinsically larger uncertainties when applied to this particular problem, but their results are consistent with those of the other methods.

Our results suggest that it is feasible to determine spectral features at high temperature by comparing results using different methods. A crucial part of understanding the systematics is the zero-temperature control analysis (see also [22, 5]). Applying this to all methods is currently in progress. The procedure outlined here may also be extended to study excited states [12], or to transport properties, where the balance of strengths and limitations of the different methods may be different.

The gauge ensembles were produced using openqcd-fastsum [23]; information about the ensembles can be found at [24]. The NRQCD correlators were produced using the FASTNRQCD package [25]. We anticipate making the full dataset and scripts for the results presented here available after a future publication.

We acknowledge EuroHPC Joint Undertaking for awarding the project EHPC-EXT-2023E01-010 access to LUMI-C, Finland. This work used the DiRAC Data Intensive service (DIaL2 and DIaL) at the University of Leicester, managed by the University of Leicester Research Computing Service on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The DiRAC service at Leicester was funded by BEIS, UKRI and STFC capital funding and STFC operations grants. This work used the DiRAC Extreme Scaling service (Tesseract) at the University of Edinburgh, managed by the Edinburgh Parallel Computing Centre on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The DiRAC service at Edinburgh was funded by BEIS, UKRI and STFC capital funding and STFC operations grants. DiRAC is part of the UKRI Digital Research Infrastructure. This work was performed using the PRACE Marconi-KNL resources hosted by CINECA, Italy. We acknowledge the support of the Supercomputing Wales project, which is part-funded by the European Regional Development Fund (ERDF) via Welsh Government, and the University of Southern Denmark for use of computing facilities. This work is supported by STFC grant ST/X000648/1 and The Royal Society Newton International Fellowship. RB acknowledges support from a Science Foundation Ireland Frontiers for the Future Project award with grant number SFI-21/FFP-P/10186. RHD acknowledges support from Taighde Éireann – Research Ireland under Grant number GOIPG/2024/3507.