QCD in a moving frame: an exploratory study

The framework of shifted boundary conditions has proven to be a very powerful tool for the non-perturbative investigation of thermal quantum field theories. For instance, it has been successfully considered for the determination of the equation of state of SU(3) Yang-Mills theory with high accuracy. The set-up can be generalized to QCD and it is expected to lead to a similar breakthrough. We present first results for QCD with three flavours of non-perturbatively O($a$)-improved Wilson fermions and shifted boundary conditions.


Introduction
The understanding of a large variety of phenomena involving the strong interactions, ranging from the dynamics inside a nucleon star to the evolution of the early Universe, crucially depends on the accurate knowledge of the equation of state (EoS) of QCD. In particular, those extreme conditions are now being reproduced and investigated at heavy-ion colliders, where the EoS is an essential input for the analysis of the data [1].
First principles determinations of the EoS of QCD are on the other hand very challenging, as one needs non-perturbative control of QCD over a wide range of temperatures. Even at relatively high temperatures, indeed, standard perturbative methods are characterized by very poor convergence and elaborated techniques are necessary to improve convergence (see e.g. [2][3][4]). Most importantly, due to the asymptotic nature of the perturbative expansion, it is difficult (if not impossible) to access the accuracy of the results within perturbation theory itself (even if apparent convergence is seen), unless one can compare with non-perturbative data over a wide range of temperatures. Therefore, lattice QCD is at present the only known framework that allows us to tackle this problem from first principles, in a fully systematic and predictive way.
Results for the EoS of QCD with N f = 2 + 1 quark flavours at zero chemical potential have been recently obtained using well-established lattice techniques [5][6][7]: the so-called integral method [8] and its variants. However, although many interesting results have been obtained with these methods, they become computationally very demanding as the temperature is increased, thus limiting in practice the accessible range of temperatures. Most calculations are in fact confined to T 500 MeV. Only very recently results at higher temperatures, i.e. up to 2 GeV, began to appear [9], but reliable continuum extrapolations are still difficult at the highest temperatures.
Speaker, e-mail: mattia.dalla.brida@desy.de For the very same reason of high computational cost, all state of the art determinations have been exclusively obtained in the framework of staggered fermions. It is, of course, of crucial importance at this point to provide independent results with also other discretizations, in order to increase our confidence in the continuum results obtained so far.
In the last few years, a big effort has been invested into devising new methods to address this problem, and overcome the limitations of the current state of the art techniques. In this contribution we consider, for the first time, the framework of shifted boundary conditions (SBC) [10][11][12] applied to QCD. This set-up has been successfully applied to the case of SU(3) Yang-Mills theory [13][14][15][16][17], leading above all to a determination of the EoS over two orders of magnitude in the temperature with half a per-cent accuracy [14,16,17]. A similar breakthrough is expected in the case of QCD.
In the following, the feasibility of a determination of the EoS in QCD with N f = 2 + 1 nonperturbatively O(a)-improved Wilson fermions and SBC is investigated, and it showed that precise results are possible over a wide range of temperatures. Two orders of magnitude in the temperature, 0.36 GeV T 45 GeV, could easily be covered.

QCD in a moving frame
In a relativistic thermal field theory one can relate the entropy of the system in its local rest frame to the momentum measured by an observer in a moving reference frame (see e.g. [18]). Hence, if one is able to measure the momentum of the system in a moving frame, one can directly access its entropy in its rest frame. Any other thermodynamic potential can be obtained from standard thermodynamic relations.
Remarkably, the Euclidean path integral formulation of a thermal field theory in a moving frame is rather simple: it amounts to imposing shifted boundary conditions to the fields in the compact (time) direction [12]. In the case of lattice QCD this means requiring: where L 0 is the physical extent of the compact direction, and ξ is the shift vector corresponding to the imaginary velocity of the system; periodic boundary conditions are assumed in the three spatial directions of extent L.
The entropy density s(T ) at the temperature T can now be computed as [12]: where · · · ξ denotes the path-integral expectation value in presence of SBC, and T R 0k corresponds to the momentum k-component of the renormalized energy-momentum tensor (EMT).
From the results of the analysis of [19][20][21], the renormalized momentum components of the EMT can be written as, where T F 0k and T G 0k are the bare fermionic and gluonic components of the EMT; Z T F and Z T G are two finite (i.e. renormalization scale independent) renormalization constants, which depend on the exact discretization of the QCD action and of the bare fields T F 0k and T G 0k . For the gluonic component of the EMT, we choose (see [16,19] for a precise definition), where g 0 is the bare coupling and F a µν are the color components of the usual clover discretization of the field strength tensor.
For the fermionic part, we take, where ∇ µ , ∇ * µ and ← ∇ µ , ← ∇ * µ are the usual forward and backward covariant lattice derivatives, acting on either the quark or anti-quark fields, respectively. (Once again, we refer the reader to [19] for a precise definition.) We note that in this exploratory investigation we may neglect any O(a) operator counterterm to these fields.
Finally, as anticipated, for the lattice action we consider N f = 2 + 1 non-perturbatively O(a)improved Wilson fermions and, if not stated otherwise, the Wilson (plaquette) gauge action.

The entropy in free-field theory
The information from lattice perturbation theory can give one useful insight into the lattice artifacts of the theory, and allows educated guesses on which set of kinematic parameters might lead to small discretization errors in non-perturbative investigations. Hence, it is useful to first study Eq. (2) in the limit of free quarks and gluons i.e. at the lowest order in perturbation theory. Of course, the information inferred from perturbation theory is expected to be more accurate the higher is the temperature, and the closer one is to the continuum limit. In particular, as the results are expected to be more trustworthy at high temperatures where quark mass effects are strongly suppressed, below we shall focus on the case of massless quarks. In Figure 1 we show the continuum limit of the ratio of the entropy density computed at the lowest order in lattice perturbation theory through Eq. (2), and the expected continuum limit, namely the Stefan-Boltzmann (SB) result for QCD with N f quark-flavours: Specifically, we present data for the relevant case of N f = 3, and different choices of shift vector ξ. We considered lattices with L 0 /a = 4, . . . , 16, and a fixed ratio L/L 0 = 32; with this choice finite volume effects can be ignored in the present discussion. The results are plotted against the temperature T in lattice units (cf. Eq. (2)). It is interesting to note that, in the free case and for N f = 3, the fermionic contribution to the entropy density is twice as larger than the gluonic one. In particular, not only the physical component of the entropy is dominated by the quark sector, but also its lattice artifacts. Nonetheless, the cases with ξ = (1, 0, 0) and (1/2, 1/2, 0) appear to have in general very small lattice artifacts: below 5% for (aT ) 2 < 0.02. For these shift vectors, the latter condition translates into having lattices with L 0 /a ≥ 6, for which discretization errors are about 2% or below. Moreover, we note that the approach to the continuum limit is clearly O(a 2 ). This is expected since, in perturbation theory, the O(a) effects in this quantity are all of O(am q ), where m q is the subtracted bare quark mass. These effects are hence absent in the chiral limit.

Simulation results
The numerical simulations have been performed using a customized version of the openQCD-1.6 package [22,23]. This allowed us to employ several efficient algorithms to speed up the simulations. More precisely, the simulation of the doublet of up and down quarks employed an optimized twisted-mass Hasenbusch preconditioning of the quark determinant [24], while the strange quark was simulated through a RHMC algorithm [25,26] with an optimized frequency splitting of the rational approximation. Even-odd preconditioning was used for both the light and strange quarks. The integration of the molecular dynamics equations was based on a three-level integration scheme. The gauge force was integrated on the finest level using a 4th-order Omelyan-Mryglod-Folk (OMF4) integrator [27], while the fermionic forces were integrated on the two coarser levels. On the finest of these we used a OMF4 integrator, while on the coarsest a 2nd-order OMF integrator [27]. The solution of the Dirac equation along the molecular dynamics evolution was obtained using a standard conjugate gradient with chronological inversion. More details on the exact implementation of these algorithms can be found in the references provided. In order to access the feasibility of a determination of the EoS in 3-flavour QCD with Wilson fermions and SBC, the first crucial issue to be addressed is the computational effort necessary to obtain a given statistical precision for the basic bare quantities T F,G 0k ξ ; particularly important is how this depends on the temperature. We considered three different ensembles of well separated temperatures; a detailed description is given below. For all simulations we chose ξ = (1, 0, 0), L 0 /a = 6, and L/a = 96, therefore having T L ≈ 11. We note that the leading finite volume effects in the entropy density are exponentially small in the product ML, where M is the lightest screening mass at the given temperature (see [12] for explicit formulas in presence of SBC). Having this noticed, based on the results of finite volume studies conducted in SU(3) Yang-Mills theory [14,17], and on perturbative estimates of the lightest screening masses in QCD [28], we expect finite volume effects to be well below the precision reached in our exploratory runs (see below). Of course, for a precise determination of the entropy density a systematic study of finite volume effects needs to be carried out. In this respect, we emphasize that considering larger spatial volumes is not a real issue in our set-up. Our computation of the entropy density is based on the measurement of a simple one-point function of the EMT (cf. Eq. (2)). The increase in computational effort due to a larger spatial volume is hence largely compensated by the statistical gain deriving from averaging the one-point function over a larger number of space-time points. Also, simulating large spatial volumes would help in controlling systematic effects arising from the freezing of topology at small lattice spacings [29,30].  Table 1. Bare action parameters for the simulated ensembles at temperatures: T 1 , T 2 , T 3 . As usual, β = 6/g 2 0 , and κ −1 = 2m 0 + 8, with m 0 the bare quark mass of either the up and down quarks (ud) or the strange quark (s).
c sw denotes the improvement coefficient of the Sheikholeslami-Wohler term.
We shall now describe in some detail the three ensembles we considered.
1. T 1 : The lattice bulk-action and action bare parameters of this ensemble match those of ensemble N203 of CLS [31,32]. Specifically, the simulations employ N f = 2+1 non-perturbatively O(a)improved Wilson fermions and the Symanzik tree-level O(a 2 )-improved (Lüscher-Weisz) gauge action (see the given references for more details). The bare action parameters are reported in Table 1 for the reader's convenience. These correspond to having a lattice spacing a ≈ 0.064 fm, and quark masses such that m π ≈ 340 MeV and m K ≈ 440 MeV at T = 0. The resulting temperature with our choice L 0 /a = 6 and ξ = (1, 0, 0) is T 1 ≈ 0.36 GeV.
2. T 2 : For this ensemble the bare coupling was fixed by requiring the Schödinger functional (SF) coupling at the energy scale µ 0 = T 2 × 3 √ 2/4 to have the valueḡ 2 SF (µ 0 ) = 2.012 [33,34]. From the results of [35] one infers that T 2 ≈ 4 GeV. The bare quark masses were all set to approximatively their critical value, hence simulating the theory at very small quark masses. This was feasible since at high temperature the massless Dirac operator develops a spectral gap λ ∝ T . 1 The critical mass was estimated by extrapolating to infinite volume the results obtained in finite volume SF simulations [34]. 2 The precise set of parameters considered is given in Table 1. In this respect, we must note that using the results for the quark-mass renormalization determined along the lines of constant physics defined through the SF couplings [36], it will be possible in future runs to tune the quark masses to the their physical values over a wide range of temperatures, i.e. up to T ≈ 70 GeV. 3 3. T 3 : Analogously to the case of the T 2 ensemble, also for this ensemble we fixed the bare coupling through the SF coupling, havingḡ 2 SF (µ) ≈ 1.27 with µ = T 3 ×3 √ 2/4. This corresponds to T 3 ≈ 45 GeV [33,35]. The quark masses were all set (approximatively) to their critical value, estimated from small volume SF simulations (cf. Table 1).  Table 2. Results for the bare expectation values T F,G 0k ξ at the three chosen temperatures: T 1 , T 2 , T 3 . A total of N ms measures was considered, each separated by 10 MDU. The estimated CPU time required for the generation of the ensembles is also given.
The results for the bare expectation values T F,G 0k ξ corresponding to the ensembles T 1 , T 2 , T 3 are given in Table 2. The table provides the total number of measurements gathered. These were collected every 10 MDU along the Monte Carlo histories, detecting no autocorrelations. In the table we also report the CPU time invested in the calculations. As one can see, the figures are modest, while the statistical precision reached on the bare quantities is remarkable: 0.5 − 1%, depending on the temperature, for the fermionic contribution T F 0k ξ , and slightly below 2% for the gluonic contribution T G 0k ξ . The results are particularly encouraging since, as expected for N f = 3 QCD, the fermionic contribution to the entropy is significantly larger than the gluonic one, but it can be determined more accurately for a given number of independent measurements, especially at high temperature.
A few miscellaneous observations are in order at this point. Firstly, if we had measured the observables more frequently, we would have obtained a significantly better precision at essentially the same computational cost. At T 2 and T 3 for example, we observed that the autocorrelations of T F,G 0k are in fact below 2 MDU: much shorter than we originally expected. As the generation of a gauge configuration is for all ensembles substantially more expensive than the measurement of the observables (about a factor 9 larger), measuring every 2 MDU instead of 10 MDU would reduce the computational effort necessary to obtain a given statistical precision by an interesting factor. Secondly, for these exploratory runs we did not consider any optimization of our code to run on KNL processors. Another interesting speed-up factor might thus be obtained even with a basic optimization strategy. 4 When moving from intermediate to high temperatures we observe that a precise determination of the bare expectation values T F,G 0k ξ becomes easier. The main reasons is because the simulations do become easier. This is so because the theory to be simulated approaches the simple free theory, and the spectral gap of the Dirac operator gets larger as the temperature increases; this is a similar situation to step-scaling simulations with the SF [33,34]. In addition, we note that the relative error of the fermionic contribution T F 0k ξ decreases for a given number of independent measurements. Finally, it is also interesting to note that the simulation algorithm appears to sample topology well at the lowest temperature T 1 , where a ≈ 0.064 fm. There we find: Q 2 ξ = 5.2(3). 5 At the higher temperatures T 2 and T 3 , however, quantitative conclusions cannot be drawn. Certainly, one would expect topology to be very much suppressed at these high temperatures, but concurrently the simulation algorithm must suffer from a severe freezing of the topology at these very small lattice spacings [38,39]. All we can say is that our runs seem to exclusively sample the Q = 0 sector.

Conclusions
The results we presented show that within the framework of SBC, a precision of 0.5 − 1% on the bare quantities entering the determination of the entropy density is reachable with modest computational effort and over a wide range of temperatures, easily covering two orders of magnitude in the temperature. This is very encouraging in view of a precise determination of the EoS of N f = 2 + 1 QCD using non-perturbatively O(a)-improved Wilson fermions. In particular, differently from the situation that occurs with more standard methods, the higher the temperature is, the easier the simulations become, and thus the easier is to obtain precise results. The fundamental reason for this favorable situation is that the problem of computing the bare expectation values and their renormalization are completely disentangled in our framework, and can therefore be tackled separately. The important question at this point is hence whether the necessary renormalization of the EMT can be obtained with a competitive level of precision. This issue is under current investigation and the details of our renormalization strategy are in preparation [40].
Regarding the renormalization of the EMT we must note that several interesting ideas, based on the Yang-Mills gradient flow [37], have been proposed [41][42][43] (see [44] for a review). Following some of these strategies promising results for the EoS have already been obtained, both in SU(3) Yang-Mills theory and QCD [45,46]. From a different perspective, other promising ideas for determining the EoS have been devised and tested [47,48].
Finally, another important issue we shall address in order to obtain an accurate determination of the EoS using Wilson fermions is a systematic study of the O(a)-improvement of the EMT. In this respect, it is interesting to note that chiral symmetry "restoration" at high temperature may have interesting consequences on the lattice artifacts of the theory. Following standard arguments [49], one expects indeed some degree of "automatic" O(a)-improvement at high temperature [40].