The sphaleron rate from Euclidean lattice correlators: an exploration

We show how the sphaleron rate (the Minkowski rate for topological charge diffusion) can be determined by analytical continuation of the Euclidean topological-charge-density two-point function, which we investigate on the lattice, using gradient flow to reduce noise and provide improved operators which more accurately measure topology. We measure the correlators on large, fine lattices in the quenched approximation at $1.5\,T_c$ with high precision. Based on these data we first perform a continuum extrapolation at fixed physical flow time and then extrapolate the continuum estimates to zero flow time. The extrapolated correlators are then used to study the sphaleron rate by spectral reconstruction based on perturbatively motivated models.


I. INTRODUCTION
The classical vacuum of QCD is not unique; there are an infinite number of topologically distinct classical vacua defined by different integer values of the Chern-Simons number. Transitions between different classical vacua are important because such transitions change the axial quark number, which would otherwise be conserved for massless quarks and slowly changing for light quarks. The sphaleron rate is defined as the mean-squared change in Chern-Simons number per spacetime volume. It was first considered within QCD because of the role it may play in electroweak baryogenesis [1,2]. More recently it has raised interest because magnetic phenomena in heavy ion collisions might give axial charge density an important role to play in, for example, the chiral magnetic effect [3,4]. At high temperatures, where the QCD coupling is weak, there are semiclassical approaches to determine the sphaleron rate, but they falter long before one considers physically relevant temperatures [5]. Therefore we will make a first attempt to study the sphaleron rate nonperturbatively using the methods of lattice QCD and spectral reconstruction.
We begin by reviewing the definition and properties of the sphaleron rate in QCD. Topology in QCD is determined through the topological charge density q(x) and its integral, the topological charge Q, where µνρσ is the totally antisymmetric tensor and F µν is the field strength tensor. It can be proven that Q is an integer for smooth vacuum-to-vacuum configurations [6] and is equal to the number of left-handed zero modes minus the right-handed zero modes of the massless Dirac operator [7]. In Minkowski time, the topological charge density determines the violation of axial current conservation: neglecting quark masses, the axial current where N f is the number of light quark flavors. The sphaleron rate Γ sphal is defined as the rate of the mean squared change in the topological charge per Minkowski 4-volume, or equivalently the integration of the Wightman correlator of the topological charge density. It can be reexpressed in terms of the spectral function ρ(ω) of the topological-charge-density two-point correlator at zero spatial momentum, which in turn is related to the Euclidean topologicalcharge-density correlation function G(τ ) through an integral relation, Here T is the temperature and τ is the temporal distance between the two charge-density operators, and the unusual minus sign arises because q is a time reversal odd operator.
The sphaleron rate is of interest both in electroweak baryogenesis [2,5,8] and heavy ion collisions [3,4,9]. For electroweak interacting matter the sphaleron rate has been well understood and determined using Bödeker's effective theory [10][11][12][13] in the weak-coupling regime. These innovations allowed for a nonperturbative semiclassical real-time evaluation on a Minkowski lattice [14,15]. A similar study for the weak-coupling behavior of the SU(3) sphaleron rate is given in [5], while for the physically interesting coupling regime a direct evaluation in Minkowski space is impossible. The results of nonperturbative lattice studies are also rather limited. A recent lattice QCD study on the SU(3) sphaleron rate can be found in [16].
In this work we attempt to give a constraint on the sphaleron rate by fitting Euclidean correlators to perturbatively motivated models in different scenarios. Due to high-frequency fluctuations of the gauge fields on the lattice, the correlators considered in this work are noisy, and some noise reduction method is necessary. Gradient flow [17][18][19] is the preferred method because it has a valid definition in terms of continuum field theory and we have a good analytical understanding of how it affects the gauge fields. This makes it preferable to traditional approaches such as cooling [20,21], or APE [22], stout [23], or HYP [24] smearing. Related studies can be found in [25][26][27]. Since its introduction, gradient flow has proven to be useful for a variety of issues in lattice QCD [28][29][30][31][32]. It also provides access to investigate how instantons emerge in the continuum limit of lattice QCD [16,[33][34][35]. In this work the gradient flow is mainly employed as a noise reduction technique. The field smearing nature of the gradient flow enables the creation of smooth gauge configurations on the lattice, from which a well-defined topological charge can be obtained [18,36].
This study is closely related to our previous work [37], which shares the same implementation of the gradient flow and lattice setup. The strategy for the continuum and flow-time-to-zero extrapolation is also similar. The subtle differences will be addressed in later sections.
In the following we briefly introduce the lattice setup, including our implementation of the topological charge (density) and gradient flow on the lattice (Sec. II). In Sec. III we examine under what conditions our definition of topological density correctly returns the topology, and we investigate whether topological freezing will affect our results. Section IV is devoted to the continuum extrapolation at fixed flow time and the subsequent extrapolation to zero flow time. In Sec. V we review the perturbative spectral function and combine it with various models for low-frequency contributions to extract the sphaleron rate from the lattice data. We then discuss the results and conclude in Sec. VI.

II. LATTICE SETUP
The setup in this work is identical to that of our previous work [37]. We perform numerical calculations of SU(3) Yang-Mills theory in 4-dimensional spacetime with periodic boundary conditions for all directions. The configurations are generated on large, fine isotropic lattices in the quenched approximation using the standard Wilson gauge action. To make sure that the configurations are sampled from the thermal equilibrium we first perform 5000 heat-bath sweeps. Afterwards we save one configuration after every 500 combined sweeps, where a combined sweep consists of one heat-bath and four overrelaxation sweeps. We have checked that this procedure eliminates autocorrelations between saved configurations in all quantities we consider except for total topology, which we discuss later.
With the aim of providing results in the continuum limit, we measure the topological charge density correlators on five different lattices. For each lattice the β value is tuned according to r 0 T c = 0.7457 (45) [38] so that all are at almost the same temperature of 1.5 T c . The scale setting is done via the Sommer parameter r 0 [39] with parametrization from [38] and updated coefficients from [40]. We summarize the key information in Tab We use a gluonic definition of the topological charge density q(x) based on an a 2 -improved implementation of the field strength tensor. Rather than the standard "clover" sum of four square plaquettes, we use a mixture of square and 1 × 2 rectangular plaquettes [41] F Imp µν (n) = where C (m,n) µν (n) are field strength tensors constructed from the variant clover terms made up of plaquette rectangles W m×n µν , see [41] for details. This is, however, insufficient to provide a q(x) definition which captures topological information, as it suffers both from large renormalization issues and nontopological noise. In order to obtain a well-defined definition of q(x) that accurately captures topological information, we need to apply gradient flow to smooth out the highest-frequency fluctuations in the gauge fields.
The gradient flow introduces a flow-time variable τ F and a gauge field B µ (τ F , x) which is set equal to the physical gauge field A µ (x) at τ F = 0 and subsequently evolves with increasing flow time τ F by gradient descent under the Yang-Mills action: The superscript B indicates that D B ν G B νµ is evaluated using B µ (τ F ) rather than A µ . We adopt a Symanzik improved discretization of the gradient flow called Zeuthen flow [42], which introduces no new O(a 2 ) discretization errors. The discretized flow equation is integrated in small steps using a third-order Runge-Kutta algorithm with an adaptive step-size. The same numerical implementation of the gradient flow is also employed in [37]. The observables are measured on flow times √ 8τ F T ∈ {0, 0.001, . . . , 0.2, 0.205, . . . , 0.3}. Gradient flow can be understood as a modification of the operators used in the measurement, replacing the elementary links with unitarized averages over many paths, an extreme form of the use of "fat links." At lowest perturbative order it is equivalent to replacing the gauge fields with their averages over a Gaussian envelope with width √ 8τ F . But gradient flow is a nonperturbatively well-defined and gauge invariant procedure.
Using gradient flow to modify our operators is a twoedged sword. As we increase the amount of flow, the elimination of UV fluctuations makes the correlation functions less noisy. In addition, the lattice definition of q(x) is always contaminated by nontopological high-dimension operators suppressed by powers of the lattice spacing a, and if fluctuations with wave number k ∼ 1/a are present, these operators contaminate our measurement with nontopological effects. By eliminating such shortdistance fluctuations, gradient flow improves the topological behavior of q(x). As we will soon see, there is a lattice spacing dependent flow depth above which q(x) correctly captures topological information. Smaller values cannot be trusted and must be avoided. But the value of the correlation function is τ F -dependent and only the τ F → 0 limit corresponds to the desired correlation function. This limit must be extracted by extrapolation; but as we increase τ F , the difference eventually becomes enough that the finite τ F correlator is not useful in trying to extrapolate to τ F → 0, as very coarse lattices do not help us to take the continuum limit. On dimensional grounds, one expects that, for a correlation function measured over a separation τ , these corrections will be dependent on τ F /τ 2 . Only sufficiently small τ F /τ 2 values will be in a useful scaling window where they help determine the τ F → 0 extrapolation. For larger values physical information is lost; indeed, above a certain ratio of τ F /τ 2 the correlation function even has the wrong sign [43]. Therefore there is also a τ -dependent upper limit on available τ F values; Ref. [43] conservatively advises √ 8τ F < 0.33τ .

III. ROLE OF TOPOLOGY
Because the observable q is a topological density, it makes sense to first check what role topology plays in our calculations. This means, first, to see what τ F value is necessary before q returns accurate topological information, and second, to check that topological freezing does not significantly affect our results.
If we define q using the field strength from Eq. (5) but without applying any gradient flow, the desired topological density operator mixes with nontopological higherdimension operators which contribute large noisy contributions to the determined topology. This effect is ameliorated as we increase the amount of flow. To examine this in detail, we study the topological susceptibility χ ≡ Q 2 /Ω where Ω is the spacetime volume of the lattice. We compute this on each lattice at a series of flow times τ F . Our results, shown in Fig. 1  tioned artifacts significantly contaminate the determined topological susceptibility. We will therefore only rely on data with τ F ≥ 0.5a 2 in what follows, so that we can be confident that our operator faithfully represents the topological density.
In analyzing the susceptibility, we found that Q evolves slowly on our finer lattices, a phenomenon called topolog-ical freezing. This is visible in Figure 1 in the much larger error bars for the susceptibility for the finest lattices. If the correlation functions we are interested in are sensitive to the topological sector, then this can cause large autocorrelation problems and poor statistical power in our results. However, it is not clear that correlation functions of the topological density should be especially sensitive to the topological sector. After all, the sphaleron rate which we are seeking to compute is not related to the topological susceptibility; at weak coupling the sphaleron rate is completely unrelated to Euclidean topology and instantons [44]. Therefore, we will examine how much the topological charge density correlators are affected by topological sector, and therefore how much of a problem topological freezing may be. We do this by measuring the qq correlators separately in different topological sectors, and then analyzing how much a misweighting or exclusion of the Q = 0 sectors might affect our results. We use the N τ = 16 lattice for this study because the topological susceptibility is relatively well determined there. To improve our determination of the properties of the Q = 0 sector, we perform a Markov chain in which we prevent Q = 0 via an extra accept-reject step, to get a pure sample of Q = 0 configurations. We can then compare this with both the Q = 0 sector, and the proper mixture of Q = 0 and Q = 0 as determined from the topological susceptibility. The results are shown in Fig. 2. The Q = 1 sector disagrees with the Q = 0 sector by more than their error bars, but nevertheless the difference is so small that, given the small fraction of Q = 1 configurations in the full sample, the difference is smaller than our statistical errors even for the largest flow time where the errors are the smallest. This implies that a misweighting of the Q = 1 sector, even by a factor of 2 relative to its correct weighting, will change the correlation functions of interest by less than their statistical errors and can therefore be ignored. Since the finer lattices have a still smaller topological susceptibility and therefore even less Q = 1 configurations, we expect that this result holds for those lattices as well. Therefore we will proceed without considering topological freezing further.

IV. DOUBLE EXTRAPOLATION
Physics resides at zero lattice spacing, and Eq. (4) only holds for the correlation function at τ F = 0. Therefore it is necessary to perform both a continuum and a flowtime-to-zero extrapolation. We first perform the continuum extrapolation at each flow time, and then use these continuum-limit values of G τF (τ ) to extrapolate to τ F → 0.

A. Interpolation
In order to perform a continuum extrapolation, we first have to interpolate the discretized data of the coarser lattices to the distances of the finest one. For that we use cubic splines in which the first derivative is constrained to vanish at τ T = 0.5 which reflects the symmetry of the correlator on the lattice.
The data is only interpolated in the range τ T ∈ [0.166, 0.5]. The correlator at shorter distances is not very helpful in extracting the low-frequency spectral function, and our scaling window of useful τ F values closes up at small τ T , so we make no attempt to study smaller τ T values. A set of interpolations at one selected flow time is shown on the left side of Fig. 3.

B. Continuum extrapolation
The continuum extrapolation is performed by linear fits in 1/N 2 τ of the interpolated data at separations of the finest lattice (N τ = 36) using the Ansatz where b is the continuum correlator normalized by T 5 .
To estimate the statistical error, we do bootstrap resampling with 10000 samples, where each sample consists of 10000 configurations. The interpolations and extrapolations are done on each sample. In the next subsection, we will define minimum and maximum ranges for τ F to ensure that we are in a scaling window for the τ F extrapolation. In some cases that window will be in conflict with the condition τ F > 0.5a 2 obtained in Sec. III. For each τ value, we will then exclude those lattices where this conflict occurs from the continuum extrapolation; otherwise we extrapolate using all available lattices.
The left-hand side of Fig. 3 shows the topological charge density correlators at flow time √ 8τ F T = 0.15 for each lattice spacing. The error bands belong to the interpolations, and we see that the relative cutoff effects of the lattices are small. This is also visible on the righthand side of Fig. 3, where we show the correlator as a function of 1/N 2 τ for some selected separations at flow times chosen from a flow range which we will define in the next section (see Eq. (8)). The continuum extrapolated values are within the error bars of the lattice data in every case. This is not surprising since the gradient flow produces renormalized operators that are insensitive to lattice-scale fluctuations.

C. Flow-time-to-zero extrapolation
In order to perform an extrapolation to zero flow time, we first need to decide on the flow time range for the extrapolation. In Ref. [43] the authors propose a criterion for the upper limit of the flow time at given separation; it is determined by allowing the leading-order term of  this constraint appears to be overly strict: the correlator's τ F dependence is still moderate even for somewhat larger τ F values, which can therefore be used in the small flow-time extrapolation. We loosen the criterion to a 20% deviation, which leads to a maximum flow time within which a linear extrapolation is applicable, namely 8τ max F T = 0.5220τ T . Because very small τ F values are noisy, we then fit in the range A general analysis of continuum operators on gradientflowed field configurations shows that, in terms of unflowed operators, they can be expanded as an operator product expansion, and correspond to the desired operator, possibly with a renormalization factor, plus a series of high-dimension operators with compensating positive powers of τ F as determined by operator dimension; see for instance Refs [29,45]. In our case we know that the topological charge does not renormalize (as we easily verify by seeing that it integrates to an integer) and that the high-dimension operators must vanish on space integration, implying that they are of form, e.g., τ F D 2 q. Such a contaminating high-dimension operator does not affect the determined total topology qd 4 x, but it does affect the correlation functions, leading to corrections of order τ F /τ 2 based on dimensional reasoning. Therefore, an appropriate Ansatz for our correlation function at small τ F , incorporating the lowest-order corrections, is where d is the correlator at zero flow time, again normalized by T 5 . The extrapolation procedure is illustrated in Fig. 4. Each curve represents the continuumextrapolated G τF (τ ) at a fixed separation τ , as a function We see that the extrapolated value and the value at the largest-used τ F differ by at most 20%, indicating that a small-τ F expansion should still be valid and that higher-order τ 2 F terms should only give a few percent corrections and are not (yet) needed.
The topological charge density q(x) is odd under time reflections. Due to arguments based on reflection positivity [46,47], the correlation function has to be negative for all nonzero separations τ T = 0 in the continuum [46,48]. Fig. 5 shows our final continuum and flow time extrapolated correlator. We indeed observe that the correlator is negative in the range we are analyzing, that is τ T > 0.

V. SPHALERON RATE
The sphaleron rate is determined by the smallfrequency limiting behavior of the spectral function as shown in Eq. (3), but the Euclidean function receives large contributions from the high-frequency tail of the spectral function, particularly at smaller temporal separations, as shown in Eq. (4). Therefore we need to incorporate as much information as we can about this large-frequency region. Fortunately, the behavior at large frequencies should be more perturbative than the smallfrequency part, and the spectral function has been computed both at leading order (LO) and at next-to-leading order (NLO) in the coupling in Ref. [49]: Note that our definition of the spectral function differs from that in Ref. [49] by a relative minus sign. Here d A = N 2 c −1 = 8 is the dimension of the adjoint representation, which counts the number of gluon states. At leading order one does not obtain a prescription on determining the value of the running coupling. We therefore make an educated guess of the renormalization point choosing the value from the 1-loop order "EQCD" setup yielding (Eq.(5.26) in [49]) Using this relation the coupling is fixed to the value g 2 μ opt(T ) = 2.2346 at T = 1.5T c , where we use an updated relation T c = 1.24Λ MS [38]. We introduce an overall scaling coefficient B in the LO models to compensate for the possibly bad choice of value of the coupling constant. As determined in [49], at NLO, see Eq. (11), the optimization of the scaleμ and the running of the strong coupling constant with ω is possible in the regime ω πT . In this regime the function φ T (ω), defined in Eq.(4.4) of [49], is small which allows to define the optimized renormalization point as a function of ω as (Eq.(5.25) in [49]) For values of ω outside this regime one again falls back to the renormalization point given by Eq. (12). Following the prescription given in [49], one uses the larger value of Eq. (12) and Eq. (13) for given ω. The switch between the two formulations happens at ω/T = 19.456π. However, we also introduce the scaling parameter B in the NLO models in order to compensate for higher order corrections to the value of the renormalization point as well as other uncertainties in renormalization. Because the perturbative series is not rapidly converging at this temperature, we fit the lattice data using both the LO and the NLO spectral function, considering the difference as an estimate of the uncertainties which arise due to our incomplete knowledge of the spectral function's high frequency functional form.
First, we consider a fit to just the leading-order spectral function. The fit is poor, with χ 2 /d.o.f = 68.6. Replacing the leading-order spectral function with the NLO expression from Eq. (11), even allowing for an additional multiplicative overall rescaling B, also gives a poor fit, with χ 2 /d.o.f. = 33.2. Therefore it is necessary (and theoretically justified) to incorporate an additional structure, representing a low frequency contribution to the spectral function. Due to a lack of theoretical guidance, we will consider three possibilities. The first is a simple δ-peak in ρ/ω whose overall coefficient A/T 4 is treated as another fit parameter. This model is theoretically motivated by the appearance of a sharp feature in perturbative calculations [50]. We call this Model M1 when combined with the LO given in Eq. (10) times the scaling parameter B and M4 when combined with B times the NLO spectral function from Eq. (11). Because the actual coupling is rather large, the assumption of a sharp peak may not be reliable, so we also consider a model in which the peak is broadened into a Breit-Wigner distribution, e.g., ρ peak /(ωT 3 ) = (A/T 4 )CT 2 /(C 2 T 2 + ω 2 ). We treat A/T 4 as a fitting parameter and consider a few distinct C values, varying from a rather narrow to a quite wide structure. We call this model M2 when combined with Eq. (10) and M5 when combined with Eq. (11). Finally we consider the large-width limit of this form, ρ peak /(ωT 3 ) = A/T 4 , which is model M3 or M6. In summary, our models are: The sphaleron rate is 2T 4 times the value of the righthand side at ω = 0. For M2 and M5 the sphaleron rate can be calculated as Γ sphal /T 4 = 2A/CT 4 while for M3 and M6 it is Γ sphal /T 4 = 2A/T 4 . For M1 and M4, the δ-peak leads to an infinite sphaleron rate.
The fit results are summarized in Tab. II. In each case, our data rather tightly constraints the fit parameters, and the quality of the fit (χ 2 /d.o.f.) is fairly good. The fits based on the NLO high-frequency spectral function are slightly better than those using the LO spectral function, but no model for the low-frequency behavior is clearly preferred. We also show the ratio of fit correlators to the lattice data, and the resulting spectral functions in Fig. 6 and Fig. 7.
We use the Levenberg-Marquardt algorithm for χ 2fitting and stop the iteration when a relative tolerance of 10 −7 is reached for the fit parameters or the χ 2 /d.o.f. Only data points at or beyond τ T = 0.25 are used in the fit, since the transport information is mainly encoded in the correlator at large separations and the τ F -window over which we can perform the flow-time-to-zero extrapolation closes up as we go to shorter and shorter separations. Now we calculate the sphaleron rate in different models according to Eq. (3) and summarize them in Fig. 8. From our analysis we can see that a δ-like transport peak and a linear-in-frequency transport peak are just special cases of a Breit-Wigner transport peak with zero and infinite width. Within the set of fit models we have considered, we find that the sphaleron rate varies within the range Γ spha /T 4 ≥ 0.030(2), based on M1-M3, Γ spha /T 4 ≥ 0.024(2), based on M4-M6.
The first line is based on the LO high-frequency spectral function and the second line is based on the NLO spectral function. Note that our results indicate that a spectral function with only a high-frequency part does not provide a good fit; there must be a low-frequency structure in addition. The functional forms we have considered contain two extremes for such a form; an infinitely sharp peak (M1/M4) and an infinitely broad peak (M3/M6). Therefore we consider them to span the range of likely functional forms for a "peak," and we propose that the lower value found can be viewed as a lower limit on the sphaleron rate. We are aware that this claim depends somewhat on our choice of fitting functions considered, but again we have considered a broad range with peak structures ranging from sharp to perfectly wide, so we consider our bound to be reasonable.  The numbers in the parentheses are statistical uncertainties from a bootstrap analysis. The sphaleron rate Γ sphal /T 4 is calculated from the fitted parameters A/T 4 and C. The meaning of each fit parameter and how the sphaleron rate is determined in each ansatz can be found in Sec. V. In ansatz M1-M3 the LO spectral function has been used at large frequencies while for M4-M6 the NLO spectral function has been used.

VI. CONCLUSION
We have calculated the topological charge correlation functions on 5 different large and fine isotropic lattices in pure-glue QCD at a temperature T = 1.5T c . To improve the signal we used the gradient flow method. Within the framework of the gradient flow we developed a methodology to perform reliable double-extrapolations for the topological charge correlators. The correlators extrapolated to the a → 0 and τ F → 0 limit are then used in the spectral reconstruction, where we use perturbatively inspired models considering uncertainties from dif-ferent sources. Using either the LO or the NLO spectral function alone leads to a very poor fit; the data demand the addition of a low-frequency structure. By considering a range of such structures, the lowest value we can obtain for the sphaleron rate, for the case that the added structure is completely flat, is Γ spha /T 4 ≥ 0.024 at 1.5T c . Low-frequency structures containing an actual peak give higher values for Γ spha . If only the LO spectral function is used in the fit, then a stronger constraint, Γ spha /T 4 ≥ 0.030 is obtained.
In our opinion there are a few take-home messages from our work. First, the use of gradient flow is essential to en- sure that the topological charge correctly captures topological information. However, only a limited flow time range is useful for the flow-time-to-zero extrapolation. At small separations and coarse lattices these constraints cannot both be maintained, and we are driven to very fine lattice spacings and can only use the larger available separations. Second, although we were able to extract continuum and flow-time-to-zero extrapolated data with a few percent statistical error bars, the lack of clear theoretical guidance in fitting the spectral function leads to rather severe theory errors in the resulting sphaleron rate. We can place a lower bound on the sphaleron rate, but not a useful upper bound. It would be useful to improve our theoretical understanding of the expected form for the topological charge-density spectral function, particularly in the low-frequency region; such improvements might allow a reasonably good determination of the sphaleron rate. All data from our calculations, presented in the figures of this paper, can be found in [51].