Topological susceptibility at $T>T_{\rm c}$ from master-field simulations of the SU(3) gauge theory

The topological susceptibility is computed in the SU(3) gauge theory at temperatures $T$ above the critical temperature $T_{\rm c}$ using master-field simulations of very large lattices, where the infamous topology-freezing issue is effectively bypassed. Up to $T=2.0\,T_{\rm c}$ no unusually large lattice effects are observed and the results obtained in the continuum limit confirm the expected rapid decay of the susceptibility with increasing temperature. As a byproduct, the reference gradient-flow time $t_0$ is determined in the range of lattice spacings from $0.023$ to $0.1\,{\rm fm}$ with a precision of 2 per mille.


Introduction
The temperature dependence of the topological susceptibility χ t in QCD is of interest in connection with the dark-matter candidacy of the axion, a hypothetical particle related to the so-called strong CP problem [1][2][3][4]. Computations of χ t in numerical lattice QCD are however not straightforward for various reasons. A direct sampling of the topological charge is often impractical, for example, because the simulation algorithms tend to get trapped in a fixed-charge sector of field space. Another source of difficulty is the fact that the susceptibility decreases rapidly at high temperatures and consequently becomes more and more sensitive to lattice effects.
Most computations of the topological susceptibility at temperatures T larger than the critical temperature T c performed to date [5][6][7][8][9][10][11][12][13][14] rely on some form of reweighting or the so-called integral method, where χ t is obtained by integrating its derivative with respect to T from low to high temperatures. The systematic uncertainties and the statistical errors are generally fairly large in these calculations, particularly so when the light quarks (which lead to an additional chiral suppression of χ t ) are included.
Master-field simulations [15] bypass the topology freezing issue by simulating lattices with four-dimensional volumes V satisfying (1.1) Fixed-topology effects are of order 1/V in this case [16,17] and are thus parametrically smaller than the statistical errors, which decrease like V −1/2 at large V . In the present paper, master-field simulations are used to calculate the topological susceptibility in the SU(3) gauge theory at temperatures approximately equal to 1.5 T c and 2.0 T c . The study also serves as a first test of the feasibility of such simulations at non-zero temperatures, where having a physically large three-dimensional volume may be of some general interest.
In the next section, the theoretical framework is described in more detail. Since the topological susceptibility is rapidly varying with temperature, its extrapolation to the continuum limit requires a highly accurate scale setting. A separate computation of the reference gradient-flow time t 0 [18] was therefore performed using master-field simulations at vanishing temperature. The computation of χ t is discussed in sect. 3 and conclusions are drawn in sect. 4.

Lattice theory
The SU(3) Yang-Mills theory studied in this paper is set up on hyper-cubic L 0 × L 3 lattices with spacing a and periodic boundary conditions in all directions. At high temperatures T = 1/L 0 , the time extent L 0 of the lattice is always taken to be much smaller than its spatial size L. For the gauge action the Wilson plaquette action [20] with bare coupling g 0 is chosen.

Definition of χ t
Since the correlation function of the topological density (where F µν denotes the field strength of the gauge potential) has a non-integrable short-distance singularity, the topological susceptibility is only formally given by A sensible definition of the susceptibility in the continuum theory must therefore be provided before it can be computed on the lattice. In the present context, the susceptibility is tied to the flavour-singlet U(1) chiral symmetry of QCD, which becomes a non-anomalous symmetry when the axion field is included in the theory. The soft breaking of the symmetry by the quark masses then leads to the well-known formula relating the axion mass to χ t , provided the latter is defined consistently with the chiral Ward identities. When this condition is met, χ t is unambiguously determined and can be shown to be given by a singularityfree expectation value of "density chains" [21][22][23].
Far easier to evaluate than the density chains is the topological charge at positive gradient-flow time [18]. The associated susceptibility does not require any subtraction or renormalization [19] and is known to coincide with the susceptibility defined through the density chains, at least in the pure gauge theory [24]. All this holds in the continuum limit of the lattice theory, provided the flow time is held fixed in physical units when the lattice spacing is taken to zero. In the present paper, the topological susceptibility is measured in this way, the implementation of the gradient flow and other technical details being the same as in ref. [18].

Physical regimes at high temperatures
The topological susceptibility is a potentially complicated function of the temperature T and the spatial volume L 3 , particularly so when L is less than 1 fm, where the effective gauge coupling is small and the semi-classical approximation becomes asymptotically exact †. If L is much larger than the correlation lengths in the pseudoscalar sector, χ t is independent of L up to exponentially small terms. This regime † In the case of a four-dimensional spherical space-time, χ t can be worked out analytically in this limit and is found to be a steep function of V [25]. At non-zero temperatures, the situation is far more complicated already at the classical level [26,27]. 6.45 11 12.196(21) sets in at values of L of a few fermi, for all temperatures, but at high temperatures the bound (1.1) only holds at much larger spatial sizes. At these temperatures there is then an interesting intermediate regime, in which L is large while the variance of the distribution of the topological charge Q is much smaller than 1. It is plausible that χ t is dominated by the sectors with charge Q = ±1 in this case. Moreover, if their contribution is assumed to be suppressed by the factor exp{−S min }, S min being the minimum of the gauge action in these sectors (the instanton action), the renormalization group implies that with a logarithmically varying proportionality constant. It goes without saying that this argumentation is quite crude and that eq. (2.4) should not be taken as a solid theoretical result.

Computation of the reference flow time t 0
The extrapolation to the continuum limit of lattice results for the topological susceptibility requires a precise scale-setting. When the limit is taken, the temperature must be held fixed in units of some physical scale such as the Sommer radius [28]. Moreover, since χ t has mass dimension 4, its value must also be expressed in such units. In view of the steep temperature dependence of χ t , a relative numerical error (8t 0 ) 1/2 / r 0 Fig. 1. Plot of the simulation results for ln(t 0 /a 2 ) (diamonds) and the interpolation (2.5),(2.6). As shown by the plot on the right, setting the scale with t 0 or the available data for the Sommer radius r 0 [29,30] comes to the same within a margin of about 1% (grey band; r 0 was computed using different methods above and below β = 6.5). The sinusoidal curve is obtained from the fit function (2.5) and the one published by Necco and Sommer for r 0 /a [30].
in the reference scale thus results in an approximately 11 times larger error of the converted values of χ t .
The target statistical precision of χ t in the present paper is a few percent and the reference scale must therefore be known with errors less than a few per mille to permit unbiased continuum-limit extrapolations. This level of precision is generally difficult to reach in practice, but can be attained with a limited computational effort if the reference gradient-flow time t 0 [18] is used to set the scale.
The values of t 0 /a 2 quoted in table 1 were obtained from master-field simulations of physically large lattices. In the range of β = 6/g 2 0 considered, the lattice spacing decreases from about 0.10 to 0.023 fm. The lattice sizes L are at least 6 fm and reach values above 9 fm in some cases. On all these lattices, χ t V is in the thousands and frozen-topology effects are therefore expected to be neglible. The numbers N mf of master fields included in the measurement of t 0 /a 2 were adjusted so as to have approximately constant statistical errors of about 2 per mille. Further details of the simulations are reported in appendix A. As shown in fig. 1, the data for ln(t 0 /a 2 ) rise roughly linearly with β and can be well represented by a polynomial for the coefficients. The fit approximates t 0 /a 2 in the range 5.96 ≤ β ≤ 7.01 with an estimated error of 2 per mille. A comparison with more precise results previously obtained on small lattices [24] confirms this up to β = 6.42 and the fit also reproduces the values at β = 6.3, 6.4, . . . , 7.0 quoted in ref. [31] within errors varying from 0.2 to 1.1 percent.

Conversion to physical units
The SU (3) Yang-Mills theory is unphysical and any assignment of physical units is therefore a bit arbitrary. Often the Sommer radius r 0 is taken as the reference scale and its physical value is set to 0.5 fm. In the range 5.96 ≤ β ≤ 6.92 of validity of the fit curves of both r 0 /a [30] and t 0 /a 2 , the ratio of scales plotted in fig. 1 averages to 0.950. The traditional choice r 0 = 0.5 fm thus amounts to setting Throughout this paper the conversion to physical units is performed using eq. (2.7) and the values of t 0 /a 2 given by the interpolation (2.5).

Computation of the topological susceptibility
The computations reported in this section follow the lines of refs. [15,18] except for the fact that lattices at high temperatures are simulated.

Master-field simulations
In total six lattices were simulated, at two temperatures and three lattice spacings at each temperature, so as to allow for an extrapolation of the results to the continuum limit (see table 2). The critical temperature T c in the SU(3) gauge theory is 294 MeV [32] and the chosen temperatures T are thus about 1.5 T c and 2.0 T c . As will become clear below, the bound (1.1) is well satisfied on all lattices. Moreover, the relevant correlation lengths are much smaller than the spatial sizes L, so that the master-field simulation strategy is expected to work out.
At high temperatures, the Polyakov loop is the Polyakov loop at positive flow time, since its distribution does not require renormalization [19] and unambiguously shows the increasingly strong polarization of the loop with increasing temperature (see fig. 2). Like the freezing of the topological charge, the spontaneous breaking of the center symmetry is associated with very long autocorrelation times if the standard simulation algorithms are used.
Master fields representative of the theory in a pure phase can be built up in several steps from approximately thermalized configurations on smaller lattices. If L is not very much larger than L 0 , the simulation algorithm rapidly evolves the gauge field to a field with definite polarization of the Polyakov loop. Periodic extensions of the field in space to larger lattices preserve the polarization and long equilibration times caused by large domains with different polarization are avoided. Reflections in space preserve the distribution of the Polyakov loop too and additionally ensure that the topological charge of the field and thus its effects on the correlation functions [16,17] remain small.

Simulation results
In the continuum limit, the topological susceptibility is independent of the flow time t at which the charge density q(x) is computed, provided t is held fixed in physical units when the limit is taken. The choice of the flow time however has an influence on the size of the lattice effects. In the calculations reported here, two values of t given  (32) in units of t 0 were chosen corresponding to smoothing ranges √ 8t [18] approximately equal to 0.28 fm and 0.47 fm.
As explained in ref. [15], χ t can be obtained in master-field simulations by integrating the two-point correlation function of the charge density, up to some sufficiently large radius R, where the integral reaches its asymptotic value within statistical errors (see fig. 3 for illustration). Reflection positivity implies that the asymptotic value is approached from above with an exponential rate given by the screening lengths in the pseudo-scalar channel. The bumps in the data shown in fig. 3 and the plateaus at R ≥ 1.2 fm are characteristic features of χ t (R) on all lattices listed in table 2. At large T , small R and small flow times t, χ t (R) probes the two-point function of the topological density at short distances, where perturbation theory applies. The bumps in the data are in fact roughly matched by leading-order perturbation theory (appendix B). This computation also shows that χ t (R) is suppressed already at small R by the gradient-flow smoothing of the charge density and then gets further suppressed at larger radii by the negative (non-perturbative) long-distance contributions.
The results for the topological susceptibility quoted in table 3 coincide with the calculated values of χ t (R) at R ≃ 1.4 fm, where the asymptotic plateaus are, in all cases, safely reached within errors.

Continuum limit
The calculated values of t 2 0 χ t must be expected to depend on the lattice spacing, the leading effects near the continuum limit being of order a 2 . Statistically significant lattice effects are, however, only observed at the larger temperature considered (see   fig. 4). As further elucidated in subsect. 3.4, it is in fact no suprise that the relative size of the effects increases with temperature, since the lattice expression for the topological charge density includes non-topological contributions of order a 2 .
Linear extrapolation in a 2 /t 0 of the data listed in table 3 to the continuum limit yield results for t 2 0 χ t with errors ranging from 5.3 to 14 percent. The values obtained at the two flow times considered agree within errors, as should be the case, the ones at the larger flow time, being a bit more precise. These figures are orders of magnitude smaller than the susceptibility t 2 0 χ t = 6.67(7) × 10 −4 [24] at zero temperature and the observed rapid decrease from T = 1.5 T c to T = 2.0 T c is in rough agreement with the power law (2.4). The agreement might however be somewhat fortuitous in view of the fact that the derivation of eq. (2.4) assumes the effective gauge coupling to be small, which is not the case at these temperatures.

Miscellaneous remarks
Scaling behaviour. If both T and L are held fixed in physical units, the computational effort required for the generation of a single master field is expected to increase like a −6 when the continuum limit is approached. With respect to the integral method, which scales approximately like a −10 , this behaviour is rather mild. However, if T is increased at fixed a, L must grow too for the inequality (1.1) to remain true. While the computational effort then scales like T 7 or so, the higher cost of the simulations should be balanced against the fact that the effective statistics provided by a single master field increases proportionally to T 8 .
Improved topological charge. In all computations reported here, the standard symmetric expression was used for the topological charge density on the lattice, in which the field tensor F µν (x) is given by the so-called clover formula. A classically O(a 2 )improved expression is then up to derivative terms that do not contribute to the total charge Q. Contrary to what may be expected, the a 2 -correction in eq. (3.5) tends to increase the lattice-spacing dependence of the topological susceptibility. A complete O(a 2 )-improvement of the theory [33] and the gradient flow [34] is thus presumably required if the convergence to the continuum limit is to be accelerated.
Finite-volume effects in traditional simulations. At high temperatures T , the basic screening lengths are expected to decrease proportionally to 1/T . The approximate susceptibility χ t (R) therefore approaches its asymptotic value at large R more and more rapidly, but as suggested by fig. 3, a significant R-dependence may persist in a core range of R extending up to R = 1.2 fm or so. In traditional high-temperature simulations, where the topology freezing is overcome in ways other than through a large volume, spatial sizes L ≥ 2.4 fm are thus required to be safe of finite-volume effects.

Conclusions
Dimensional analysis suggests that the topological susceptibility grows proportionally to T 4 at high temperatures T , but instead it decreases rapidly as a result of a nearly perfect cancellation of short-and long-distance contributions. This behaviour is commonly attributed to the topological nature of the charge density q(x), i.e. to the fact that variations of q(x) with respect to the gauge field are total derivatives. None of the non-perturbatively well-defined expressions for the susceptibility known to date however embodies this property of the charge density to the extent that the smallness of the susceptibility at high temperatures would be explained. Master-field simulations provide new opportunities for non-perturbative studies of QCD. At non-zero temperatures below T c , for example, the physically large volumes that become accessible in this way allow the theory to be studied in kinematic regimes close to the thermodynamic limit, where multi-hadron states make important contributions to the partition function. Another motivation for the use of this new type of simulations is the fact that the topology-freezing issue (which tends to become severe at lattice spacings a ≤ 0.05 fm) can be bypassed in a conceptually transparent manner.
The computations of the topological susceptibility reported in the present paper could proceed straightforwardly for this reason and led to results with unprecedented precision. At temperatures higher than the ones considered here, master-field simulations however require larger and larger lattices to be simulated and thus become impractical at some point. Moreover, the topological susceptibility must be expected to be increasingly sensitive to lattice effects. To be able to control these effects, the lattice spacing must then be decreased. This second problem is, however, not specific to master-field simulations and will persist until an expression for the susceptibility is found which is naturally small at high temperatures.
All simulations were performed on a HPC cluster at CERN and on the Marconi machine at CINECA through agreements of INFN and the University of Milano-Bicocca with CINECA. We gratefully acknowledge the computer resources and the technical support provided by these institutions.

Appendix A. Simulation algorithm and other implementation details
Apart from some specific technical details related to the very large sizes of the simulated lattices, the master-field simulations reported in this paper followed established lattice-QCD strategies.

A.1 Simulation algorithm
All simulations were performed using the HMC algorithm [35] with trajectory length τ = 2. The molecular-dynamics equations were integrated by applying the forthorder integrator given by eqs. (63) and (71) in ref. [36]. This scheme proves to be highly efficient and an only mild adjustment of the step number n step was required on the larger lattices in order to preserve a good acceptance rate P acc (see table 4).
Using standard MPI communication functions, the computational work was distributed over up to 32768 processing units. Most demanding from the point of view of the memory requirements was the measurement program for the topological susceptibility, which occupied a total memory of about 16 TB in the case of the largest lattice.

A.2 Thermalization
As already indicated in sect. 3, the master fields were generated in several steps from smaller lattices, where thermalizations of the gauge field alternate with extensions to the next larger lattice through reflections at the lattice planes. The plaquette action per point is unchanged after a reflection and the topological charge vanishes, but the gauge-field tensor changes abruptly across the reflection planes, which can give rise to a low acceptance rate in the early phase of the subsequent thermalization. A few update cycles with a more accurate integration of the molecular-dynamics equations may be required in this case to get the thermalization started.
The lengths τ th of the final thermalization runs listed in table 4 are much longer than the relevant autocorrelation times. A drift in the single-field expectation values [15] has in fact never been seen after these long thermalization phases (see fig. 5 for an example). It may be worth noting in passing that outliers, such as the measurement number 25 in fig. 5, must occur with some non-zero probability, as in traditional simulations, where whole ensemble averages may be similarly outlying.
The separation ∆τ mf in simulation time of the master fields included in the computations of expectation values need not be particularly large, since any statistical correlations among the fields are automatically taken into account [15]. Autocorre- lations however lead to larger statistical errors relative to what they would be for uncorrelated fields. On the B 3 lattice, for example, the separation was duplicated with respect to the other runs for this reason.

A.3 Use of quadruple-precision arithmetic
On the simulated lattices, significance losses of up to 11 decimal places occur when the energy deficit ∆H is computed at the end of the molecular-dynamics evolution of the fields. Standard IEEE 754 double-precision data and arithmetic may be barely good enough under these conditions and it is, therefore, advisable to use quadrupleprecision artithmetic in the summation of the action densities over all lattice points. ∆H is then obtained with absolute precision given by the now practically exactly accumulated numerical errors of the densities. Assuming these are randomly distributed, their sum scales like (V /a 4 ) 1/2 and the accumulated inaccuracies are then far below any statistically relevant level.
A convenient portable implementation of quadruple-precision numbers is through pairs of double-precision numbers. Algorithms for the associated arithmetic operations were published by Dekker [37] many years ago. The subject is also discussed in a book of Knuth [38] and more extensively in an article by Shewchuk [39].

A.4 Parallel I/O
In master-field simulations, the computer time spent for field configuration I/O may not be negligible. Current HPC systems however permit the storage facilities to be accessed concurrently and thus offer a high aggregate I/O bandwidth.
In the I/O programs used in the present study, the lattice is logically divided into fairly large rectangular blocks. The part of the gauge field residing on a given block is then written out in a portable format by one of the processing units. A single field is thus stored in several files and advantage of the parallel capabilities of the storage facility is taken by having many processing units write their blocks concurrently.

Appendix B. Calculation of χ t (R) in perturbation theory
In the continuum theory and at flow time t > 0, the integrated correlation function for the approximate susceptibility (B.1). To this order of perturbation theory, χ t (R) thus depends on the summation radius R roughly like the data plotted in fig. 3 (see fig. 6). In particular, at the flow times chosen in the simulations, the maxima of the bumps in fig. 3 are at R = 0.30 fm and 0.45 fm, while the leading-order expression (B.6) has its maximum at R = 0.35 fm and 0.58 fm in these cases. The plateaus in fig. 3, on other hand, occur at distances, where perturbation theory is not expected to apply and instead goes to zero consistently with the vanishing of χ t to all orders.