Bethe-Salpeter amplitudes of Upsilons

Based on lattice non-relativistic QCD (NRQCD) studies we present results for Bethe-Salpeter amplitudes for $\Upsilon(1S)$, $\Upsilon(2S)$ and $\Upsilon(3S)$ in vacuum as well as in quark-gluon plasma. Our study is based on 2+1 flavor $48^3 \times 12$ lattices generated using the Highly Improved Staggered Quark (HISQ) action and with a pion mass of $161$ MeV. At zero temperature the Bethe-Salpeter amplitudes follow the expectations based on non-relativistic potential models. At non-zero temperatures the interpretation of Bethe-Salpeter amplitudes turns out to be more nuanced, but consistent with our previous lattice QCD study of excited Upsilons in quark-gluon plasma.

Based on lattice non-relativistic QCD (NRQCD) studies we present results for Bethe-Salpeter amplitudes for Υ(1S), Υ(2S) and Υ(3S) in vacuum as well as in quark-gluon plasma. Our study is based on 2 + 1 flavor 48 3 × 12 lattices generated using the Highly Improved Staggered Quark (HISQ) action and with a pion mass of 161 MeV. At zero temperature the Bethe-Salpeter amplitudes follow the expectations based on non-relativistic potential models. At non-zero temperatures the interpretation of Bethe-Salpeter amplitudes turns out to be more nuanced, but consistent with our previous lattice QCD study of excited Upsilons in quark-gluon plasma.

I. INTRODUCTION
Potential models give a good description of the quarkonium spectrum below the open charm and bottom thresholds, see e.g. Refs. [1,2] for reviews. Even some of the states above the threshold are also reproduced well within this model. Potential models can be justified using an effective field theory approach [3,4]. This approach is based on the idea that for a heavy quark with mass m, there is a separation of energy scales related to the quark mass, inverse size of the bound state, and binding energy, m mv mv 2 , with v being the velocity of the heavy quark inside the quarkonium bound state. The effective field theory at scale mv is the nonrelativistic QCD (NRQCD), where the heavy quark and anti-quark are described by non-relativistic Pauli spinors and pair creation is not allowed in this theory [5]. The effective theory at scale mv 2 is potential NRQCD (pN-RQCD) and the quark anti-quark potential appears as a parameter of the pNRQCD Lagrangian. Potential model appears as the tree level approximation of pNRQCD [4]. Non-potential effects are manifest in the higher order corrections. For very large quark mass, v ∼ α s 1. Therefore, the large energy scales can be integrated out perturbatively [3,4]. However, for most of the quarkonium states realized in nature this condition is not fulfilled. If Λ QCD mv 2 , all the energy scales can be integrated out non-perturvatively and the potential is given in terms of Wilson loops calculated on the lattice [3,4]. So in this limit too the potential description is justified. However, for many quarkonia, Λ QCD mv 2 , and it is not clear how to justify the potential models.
In potential models one can also calculate the quarkonium wave function. On the other hand, in lattice QCD we can calculate the Bethe-Salpeter amplitude, which in the non-relativistic limit reduces to the wave function. Thus, one can use the Bethe-Salpeter amplitude for further tests of the potential models. In particular, one can also reconstruct the potential from the Bethe-Salpeter * rlarsen@bnl.gov amplitude [6][7][8][9][10]. Most of these studies focused on quark masses close to or below the charm quark mass, though in Ref. [9] quark masses around the bottom quark have also been considered. The resulting potential turned out to be similar to the static potential calculated on the lattice, but some differences have been found. The potential description is expected to work better for larger quark masses and therefore bottomonium is best suited for testing this approach. Studying bottomonium on the lattice using a fully relativistic action is more difficult because of the large cutoff effects and the rapid fall-off of the correlators. One of our aims is to test the potential model by calculating the bottomonium Bethe-Salpeter amplitude using lattice NRQCD [11,12], which is very well suited for studying bottomonium [13][14][15][16][17][18][19][20].
The existence and the properties of quarkonia in the hot medium attracted a lot of attention in the last 30 years. It has been proposed a long time ago that quarkonium production in heavy-ion collisions can be used to probe quark-gluon plasma (QGP) formation [21]. The study of in-medium properties of quarkonia and their production in heavy ion collisions is an extensive research program, see e.g. Refs. [22][23][24] for reviews. The in-medium properties of quarkonia as well as their dissolution (melting) are encoded in the finite temperature spectral functions. Quarkonium states show up as peaks in the spectral function that become broader as the temperature increases and eventually disappear above some temperature (T ). The temperature above which no peaks in the spectral function can be identified is often called the melting temperature. Reconstructing quarkonium spectral functions from lattice calculations at nonzero temperature appeared to be very challenging (see e.g. discussions in Refs. [25][26][27][28]). The study of Bethe-Salpeter amplitudes has been proposed as an alternative method to address this problem. The idea behind this approach is to compare the Bethe-Salpeter amplitude calculated on the lattice with the expectations of the free field theory that would indicate an unbound heavy quark anti-quark pair. Bethe-Salpeter amplitudes at non-zero temperature for charmonium have been calculated in previous lattice QCD studies [29][30][31][32][33][34], but presently our un-derstanding regarding the interpretations of quarkonia Bethe-Salpeter amplitudes at T > 0 remains murky. Although using a weak-coupling approach it is possible to generalize the potential description to non-zero temperature [35,36], it is unclear if such an approach and the interpretations of quarkonia Bethe-Salpeter amplitudes is applicable in the temperature regime of interest. In this paper we focus on lattice NRQCD based determinations of Bethe-Salpeter amplitudes of Υ(1S), Υ(2S) and Υ(3S) states at T > 0. By comparing with the corresponding T = 0 results, where the interpretations of Bethe-Salpeter amplitudes are more straightforward, we point out and discuss subtleties associated with interpretations of Bethe-Salpeter amplitudes at T > 0.

II. BETHE-SALPETER AMPLITUDES AT T = 0
To define the Bethe-Salpeter amplitude for bottomonium we consider the correlation functioñ whereÕ α is the meson operator that has a good overlap with a given quarkonium state α and O r qq is a pointsplit meson operator with the quark and antiquark fields separated by distance r, Here, Γ fixes the quantum number of the meson. Furthermore, in the present work we use Coulomb gauge fixed ensembles to define the expectation value. Inserting a complete set of states we obtain the following spectral decomposition of the correlator: Assuming that only one state |α contributes at large τ , which is correct for an appropriately chosenÕ α , at large Euclidean time we havẽ where is called the Bethe-Salpeter (BS) amplitude and describes the overlap of the quarkonium state |α with the state that is obtained by letting the two field operators at distance r act on the vacuum. In the non-relativistic limit it reduces to the wave function of the given quarkonium state. Thus, up to normalization factor the Bethe-Salpeter amplitude is given by the large τ behavior of exp(E α τ )C α (τ ), with E α being the energy of quarkonium state |α , which is also calculated on the lattice. In the following we will use the terms BS amplitude and wave function interchangeably. As mentioned in the introduction, we aim to calculate the bottomonium BS amplitudes using NRQCD. We performed calculations using 2+1 flavor gauge configurations generated by HotQCD with the highly improved staggered quark (HISQ) action [37,38]. The strange quark mass was fixed to its physical value, while the light quark masses correspond to the pion mass of 161 MeV in the continuum limit [37,38]. We use the same NRQCD formulation as in our previous study [39,40]. For the calculations at zero temperature we use 48 4 lattices and β = 10/g 2 0 = 6.74 corresponding to lattice spacing a = 0.1088 fm. We use 192 gauge configurations in our analysis with eight sources per configuration.
To construct the meson operators that have the optimal projection we start with the source [40] Here, ψ i (r) is the trial wave function of the ith bottomonium state obtained by solving the Schrödinger equation with the Cornell potential modified by discretization effects [15].
to obtain the optimized operator for bottomonium state αÕ Thus, to obtain the BS amplitude we consider the large τ behavior of the following combination: In practice, the value of τ does not have to be very large. We find that τ > 0.3 fm works for all states, i.e. the resulting BS amplitudes are time independent. For τ = 0 the BS amplitude will be equal to the trial wave function ψ i (r). To obtain the proper normalization of the BS amplitude we require that ∞ 0 drr 2 |φ α (r)| 2 = 1. In Fig. 1 we show the BS amplitude φ α (r) for Υ(1S), Υ(2S) and Υ(3S) states compared to the corresponding trial wave functions ψ α (r) used to construct the optimized meson operators. We see that the r-dependence of the BS amplitudes is in qualitative agreement with the expectations of non-relativistic potential model. However, the details of the r dependence are different from the input trial wave function. We also note that the orthogonalization procedure is important for getting the correct r dependence of the BS amplitudes.
If the potential picture is valid the BS amplitude should satisfy the Schrödinger equation with m b being the b-quark mass of the potential model. Note that the reduced mass in the bb system is m b /2, hence the absence of factor two in the above equation.
Using the BS amplitude and the energy of at least two bottomonia states determined in NRQCD from the above equation we can obtain m b and the potential V (r). We determine the b-quark mass using Υ(1S) and Υ(2S) states as follows To evaluate ∇φ α we use the simplest difference scheme. The value of m b determined from the above equation for different values of quark antiquark separation r is shown in Fig. 2. The r-range was chosen such that it does not include the node of Υ(2S) and large distances, where the statistical errors are large. We see some modulation of the extracted m b in r, which may indicate that the BS amplitude cannot be completely captured by the Schrödiner equation, but there is no clear tendency of m b as function of r. Therefore we fitted the values of m b obtained for different r to a constant. This resulted in This value of the effective bottom quark mass obtained by us is not very different from the one used by the original Cornell model, m b = 5.17 GeV [41] but is significantly larger than the b-quark mass used in most of the potential models (see e.g. Ref. [42]). We also determined the value of m b using the BS amplitudes and the energy levels of Υ(1S) and Υ(3S) and obtained m b = 5.82(0.51) GeV. This agrees with the above result within the errors.
Having determined m b , we can also calculate the potential, V (r), using the BS amplitudes and the bottomo- nium energy levels as The results are shown in Fig. 3. Given our findings for m b it is not surprising that the values of the potential obtained using Υ(1S), Υ(2S) and Υ(3S) states agree within errors. In the figure we also compare the value of V (r) determined from the different states to the phenomenological potential of the original Cornell model [41] and the energy of static quark antiquark pair obtained from Wilson loops at lattice spacing a = 0.06 fm [38]. It is quite non-trivial that all these potentials agree with each other. A similar conclusion has been reached in Refs. [7][8][9] when the limit of quark mass going to infinity was taken. We note that the relativistic corrections to the spin-dependent part of the potential are quite small for the b quark mass [43] and, thus, are not visible given our statistical errors. The discussion above ignored spin-dependent effects. To address the spin-dependent part of the potential we also calculated the BS amplitude for η b (nS) states, n = 1, 2, 3. We have found that the corresponding BS amplitudes agree with the ones of the Υ(nS) states within errors. Therefore, with the present statistics we cannot resolve the spin-dependent part of the potential.
As discussed above, the r-dependence of the BS amplitudes qualitatively follow the r-dependence of the trial wave function ψ i (r) obtained from the potential model. But at qualitative level significant differences can be seen, c.f. Fig. 1. This potential model used m b = 4.676 GeV [15], which is smaller than the effective quark masses determined above. Therefore, we calculated the wave functions of (nS) bottomonium states using the static quark anti-quark energy [38] as a potential and m b = 6 GeV. The results are shown in Fig. 4 and we see that the agreement between the BS amplitude and the wave  [41] shown as solid line as well as to the the energy of the static quark antiquark pair obtained from Wilson loops using a = 0.06 fm lattice [38]. All the lattice results were normalized to coincide with the Cornell potential at r = 0.4 fm. functions is significantly improved. We also note that the dependence of the energy levels on m b is rather mild, e.g. changing m b from 4 GeV to 6 GeV only reduces the spin-averaged 2S-1S splitting by 3.5%. Thus, using large values of m b in the potential model is a viable option.

III. BETHE-SALPETER AMPLITUDES AT T > 0
We can also consider the mixed correlatorC r α (τ, T ) defined in Eq. (1) for T > 0 by evaluating the expectation value over a thermal ensemble at a temperature T = 1/β with the thermal partition function Z(β) = Tr e −βH .
Using energy eigenstates to evaluate the trace and inserting a complete set of states we obtain the following expression for the correlatorC r α (τ, T ) Since we perform calculations in NRQCD the sum over m should be restricted to states that do no contain the heavy quark anti-quark pair; heavy quark pair creation is not allowed in NRQCD. We denote those states as |m . If we write the states |n as |n γ , where index n labels the light degrees of freedom and γ labels the quarkonium states, the above expression forC r α (τ, T ) can be re-written as If we write E m γ = E γ + E m + ∆E m γ and assume that the operatorÕ α mostly projects on to quarkonium state |α we can obtain a simplified form with A * α = A * α /Z(β). In the above equation we separated out the the m = 0 vacuum contribution in the sum corresponding to the thermal trace. At small temperature the first term in the above equation is the dominant one and the correlator is approximately given by the T = 0 BS amplitude. In general, however, there is no simple interpretation of the correlatorC r α (τ, T ) in terms of some finite temperature quarkonium wave function. The temperature dependence of this correlator crucially depends on the value of the matrix elements m |O r qq |m α and m α|Õ α |m . The size of m |O r qq |m α depends on the separation r and, therefore, also the size of the thermal effect will be r dependent. For values of r that are about the size of the bottomonium state of interest the matrix elements m |O r qq |m α and m α|Õ α |m should be of similar size and, thus, the temperature dependence of C r α (τ, T ) is expected to be comparable to the correlator ofÕ α explored in Ref. [40].
We performed calculations ofC r α at six different temperatures using 48 3 × 12 lattices from HotQCD collaboration. The parameters of the calculations including the  gauge coupling β = 10/g 2 0 and number of configurations are summarized in Table I. As at zero temperature, we used 8 sources per gauge configuration.
We could use the same approach as in Ref. [40] to explore the temperature dependence of the correlator C r α (τ, T ) and define the effective mass for a fixed r aM r eff (τ, T ) = ln C r α (τ, T ) C r α (τ + a, T )) .
Now the effective mass also depends on the distance r between the quark and antiquark in the point-split current. In Fig. 5 we show the effective mass of Υ(1S) correlator as function of r and τ at the lowest and the highest temperature. The errors of the effective masses are not shown to improve the visibility. Since the energy levels in NRQCD are only defined up to a lattice spacing dependent constant, as in Ref. [39] we calibrate the effective masses with respect to the energy level of η b (1S) state at zero temperature. At large τ and r the errors are quite large and within these errors we do no see any medium effects in the effective mass at the lowest temperature. For small r the effective mass quickly reaches a plateau with increasing τ . For large r the effective mass at 151 MeV reaches the plateau from below. At the highest temperature, T = 335 MeV the r and τ dependence of the effective masses looks similar for not too large values of r. However, the behavior of the effective mass is qualitatively different for large r. In particular the effective mass does not reach a plateau with increasing τ .  [40].
Since the correlatorC r α does not correspond to a positive definite spectral function it is difficult to infer inmedium properties of bottomonia from M r eff . The large statistical errors make this even more complicated. Another way to analyze the temperature dependence ofC r α is to consider the integral At zero temperature this quantity should be proportional to exp(−2E α τ ) for sufficiently large τ . This is also expected to be true below the crossover temperature. The combination should be independent of τ and can be interpreted as the normalization of the BS amplitude. In Fig. 7 we show N norm (τ, T ) as function of τ for different temperatures normalized to one at t = τ /a = 3. For the lowest temperature as well as for T = 0 we see that N norm (τ, T ) is approximately constant as expected. Here we note that the τ range in Fig. 7 is different for Υ(1S), Υ(2S) and Υ(3S) states. This is due to the fact that the correlators C r Υ(2S) and C r Υ(3S) will be contaminated by the lowest Υ(1S) state at large τ as the projection is not perfect due to the small operator basis of only three operators used in this study. As the temperature increases we see that N norm (τ, T ) no longer approaches a constant but increases at large τ . This implies that the correlatorC r α is no longer dominated by the first term in Eq. (17). The τ -dependence of N norm (τ, T ) is larger for high temperatures and is also more pronounced for excited states, as expected.
We could also analyze the τ -dependence of N α (τ, T ) in terms of the corresponding effective masses aM Nα eff (τ, T ) = ln At large τ these effective masses should reach a plateau equal to 2E α . Our results for M Nα eff for the different Υ(nS) states are shown in Fig. 8. As before the effective masses have been calibrated with respect to the energy of η b state at T = 0. We see that at T = 0 as well as at the lowest temperature the effective masses reach a plateau corresponding to the physical mass (energy), but at higher temperatures this is not the case, in general. For the ground state the errors are large enough so that no clear medium shift can be seen, except at the highest temperature, T = 334 MeV. For the Υ(2S) the corresponding effective masses decrease with increasing τ for T ≥ 251 MeV. For the Υ(3S) we see a significant shift in M Nα eff (τ, T ) already for T > 191 MeV. The behavior of M Nα eff (τ, T ) is qualitatively similar to the behavior of the effective masses of the correlator of optimized operators studied in Ref. [40]. This corroborates the findings of Ref. [40] on the in-medium modifications of the bottomonium spectral functions. For the Υ(1S) state our findings are also consistent with other studies of bottomonium at non-zero temperature using NRQCD [44][45][46][47][48][49].
Before concluding this section we mention that so far we only discussed Υ(nS) states but very similar results have been obtained for η b (nS) states as well.

IV. COMPARISONS BETWEEN T > 0 AND T = 0 BETHE-SALPETER AMPLITUDES
If we insist on the interpretation of the correlator C r α (τ, T ) in terms of the wave function we could simply divide it by N α (τ, T ) and study the r-dependence of the corresponding ratio for sufficiently large τ . At small temperatures this ratio will have an r-dependence that closely follows the r-dependence of the BS amplitude at T = 0. In Fig. 9 we compare φ α (τ, T ) = C r α (τ, T )/N α (τ, T ) for the lowest temperature, T = 151 MeV and τ = 0.65 fm with the corresponding zero temperature BS amplitudes. For the Υ(1S) and the Υ(2S) we do not see any difference between the zero temperature BS amplitude and φ α (τ, T ). For the Υ(3S) some difference between the zero temperature and finite temperature result for φ α (τ, T ) can be seen at large r, though it is not statistically significant. In any case the r-dependence of φ α at T = 0 and T = 151 MeV is quite similar even for the Υ(3S). The lack of medium effects in the BS amplitude for T = 151 MeV is not surprising since at this temperature all bottomonia should exist as well defined states. Next, we compare φ α (τ, T ) at the lowest and the highest temperature at τ = 0.4 fm. This comparison is shown in Fig. 10 and no temperature effect can be observed. This is presumably due to the fact that for this τ value the contribution of the second term in Eq. (17) is too small. Therefore, in Fig. 11 we show our results for φ α (τ, T ) at T = 251 MeV and several values of τ . As one can see from the figure for the Υ(2S) and Υ(3S) there is a significant τ -dependence of φ α for large r. This suggests that the normalized BS amplitude cannot be interpreted simply as the wave function of inmedium Υ in potential model picture. Yet, the r dependence of φ α (τ, T ) does not change much from one τ value to another. In summary, the correlationC r α (τ, T ) shows significant temperature dependence as one would expect based on the previous studies. However, the r depen-  τ and therefore, in Fig. 11 we only show the numerical results for τ = 0.4 fm. This lack of τ -dependence indicates that Υ(1S) can exist in the deconfined medium at T = 251 MeV as a well defined state with little medium modification, in agreement with the previous studies of bottomonium at non-zero temperature based on NRQCD [44][45][46][47][48][49].
The lack of temperature dependence of the normalized BS amplitude at T > 0 at τ = 0.4 fm demonstrated in Fig. 10 has an interesting consequence. It means that φ α (r, T ) can be used as a proxy for the T = 0 BS am-  plitude at zero temperature. Since the two temperatures shown in Fig. 10 correspond to two different lattice spacings this result also implies that the lattice spacing dependence of the BS amplitude is small. Therefore, the comparison of the wave function obtained from potential model and BS amplitude obtained on the lattice with a = 0.1088 fm in Section II seems justified.

V. CONCLUSIONS
Using lattice NRQCD in this paper we studied the correlation functions,C r α , between operators optimized to have good overlaps with the of Υ(1S), Υ(2S) and Υ(3S) vacuum wave functions and simple spatially nonlocal bottomonium operators, where the bottom quark and anti-quark are separated by distance r. This correlator has been calculated at zero as well as at nonzero temperature. At zero temperatureC r α can be interpreted in terms of the Bethe-Salpeter amplitude. We have found that the r-dependence of the Bethe-Salpeter amplitude closely resembles the corresponding potential model based bottomonium wave function. Moreover, by choosing the bottom quark mass used in the Schrödinger equation to be ∼ 5.5 GeV we estimated the heavy quark antiquark potential from Bethe-Salpeter amplitudes, and found agreement with the static quark potential calculated on the lattice. These findings support the potential model for bottomonium in vacuum.
We studied the temperature and Euclidean time dependence ofC r α in terms of effective masses. For Υ(1S) we see only very small temperature and Euclidean time dependence of the corresponding effective masses, except at the highest temperature of 334 MeV. For Υ(2S) and especially for Υ(3S) significant dependence on the Euclidean time were observed, making it difficult to draw parallels between Bethe-Salpeter amplitudes and potential model based in-medium wave functions. Since the r-dependence changes very little with varying Euclidean time and temperature, focusing solely on the rdependence ofC r α at a fixed τ might lead to misleading conclusions regarding existence of well-defined Υ(2S) and Υ(3S) in medium. On the other hand, we found that the behavior of the effective masses is similar to the one previously studied by us using correlators of optimized bottomonium operators [40], supporting the picture of thermal broadening of bottomonium states.