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

TL;DR

This study calculates the topological susceptibility in SU(3) gauge theory at temperatures above the critical point using large-scale master-field simulations, effectively overcoming topology-freezing issues and confirming expected decay patterns.

Contribution

It introduces a novel application of master-field simulations to compute topological susceptibility at high temperatures, bypassing topology-freezing problems and providing precise continuum results.

Findings

Susceptibility decreases rapidly with temperature above T_c

No significant lattice effects up to 2 T_c

Precise determination of gradient-flow scale t_0

Abstract

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.

Equations

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad\chi_{t}V\gg 1.\hfill\hbox to0.0pt{\hss{\tenrm(1.1)}}\crcr}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad q(x)=-{1\over 32\pi^{2}}\,\epsilon_{\mu\nu\rho\sigma}{\tenrm tr}\{F_{\mu\nu}(x)F_{\rho\sigma}(x)\}\hfill\hbox to0.0pt{\hss{\tenrm(2.1)}}\crcr}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad\chi_{t}=\int{\tenrm d}^{4}x\,\langle q(x)q(0)\rangle.\hfill\hbox to0.0pt{\hss{\tenrm(2.2)}}\crcr}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad\langle Q^{2}\rangle=\chi_{t}V,\qquad V=L^{3}/T,\hfill\hbox to0.0pt{\hss{\tenrm(2.3)}}\crcr}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad\chi_{t}\mathrel{\mathop{\propto}_{T\to\infty}}T^{-7}\hfill\hbox to0.0pt{\hss{\tenrm(2.4)}}\crcr}}$$

$$\beta$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad\ln(t_{0}/a^{2})=\sum_{k=0}^{4}c_{k}(\beta-6)^{k}\hfill\hbox to0.0pt{\hss{\tenrm(2.5)}}\crcr}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad(c_{0},\ldots,c_{4})=(1.16390,3.37888,-1.36231,1.20666,-0.45672)\hfill\hbox to0.0pt{\hss{\tenrm(2.6)}}\crcr}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad(8t_{0})^{1/2}=0.475\,{\tenrm fm}.\hfill\hbox to0.0pt{\hss{\tenrm(2.7)}}\crcr}}$$

$$\beta$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad P({\tenmib x})=\hbox{$1\over 3$}{\tenrm tr}\{W(x)\}\hfill\hbox to0.0pt{\hss{\tenrm(3.1)}}\crcr}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad\chi_{t}(R)=a^{4}\sum_{x_{0}}\sum_{|{\tenmib x}|\leq R}\langle q(x)q(0)\rangle,\hfill\hbox to0.0pt{\hss{\tenrm(3.2)}}\crcr}}$$

$$R/a$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad t_{0}^{2}\chi_{t}=2.25(12)\times 10^{-5}\quad\hbox{at}\quad T\sqrt{8t_{0}}=1.081,\hfill\hbox to0.0pt{\hss{\tenrm(3.3)}}\cr\vskip 10.76385pt\cr 0.0pt{\hfil$\displaystyle\qquad t_{0}^{2}\chi_{t}=3.43(27)\times 10^{-6}\quad\hbox{at}\quad T\sqrt{8t_{0}}=1.434,\hfill\hbox to0.0pt{\hss{\tenrm(3.4)}}\crcr}}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad q(x)=-{1\over 32\pi^{2}}\,\epsilon_{\mu\nu\rho\sigma}{\tenrm tr}\bigl{\{}F_{\mu\nu}(x)F_{\rho\sigma}(x)-\hbox{$2\over 3$}a^{2}F_{\mu\nu}(x)[F_{\mu\rho}(x),F_{\mu\sigma}(x)]\bigr{\}}\hfill\hbox to0.0pt{\hss{\tenrm(3.5)}}\crcr}}$$

$$n_{\text{tenrm}step}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad\chi_{t}(R)=\int_{0}^{1/T}{\tenrm d}x_{0}\int_{|{\tenmib x}|\leq R}{\tenrm d}^{3}{\tenmib x}\,\left\langle q(x)q(0)\right\rangle\hfill\hbox to0.0pt{\hss{\tenrm(\hbox{\tenrm B}.1)}}\crcr}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad 8tT^{2}\pi^{2}\gg 1,\hfill\hbox to0.0pt{\hss{\tenrm(\hbox{\tenrm B}.2)}}\crcr}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad\left\langle q(x)q(0)\right\rangle=\alpha_{s}^{2}{T^{2}\over\pi^{5}(8t)^{3}r^{3}}\gamma(\hbox{$3\over 2$},r)\bigl{\{}3\gamma(\hbox{$3\over 2$},r)-4\gamma(\hbox{$5\over 2$},r)\bigr{\}},\hfill\hbox to0.0pt{\hss{\tenrm(\hbox{\tenrm B}.3)}}\crcr}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad\gamma(a,r)=\int_{0}^{r}{\tenrm d}s\,s^{a-1}{\tenrm e}^{-s}\hfill\hbox to0.0pt{\hss{\tenrm(\hbox{\tenrm B}.4)}}\crcr}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad\left\langle q(0)^{2}\right\rangle=\alpha_{s}^{2}{4T^{2}\over 3\pi^{5}(8t)^{3}},\hfill\hbox to0.0pt{\hss{\tenrm(\hbox{\tenrm B}.5)}}\crcr}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad\chi_{t}(R)=\alpha_{s}^{2}{4T\over\pi^{4}(8t)^{3/2}}f(\rho),\qquad\rho={R\over(8t)^{1/2}},\hfill\hbox to0.0pt{\hss{\tenrm(\hbox{\tenrm B}.6)}}\cr\vskip 10.76385pt\cr 0.0pt{\hfil$\displaystyle\qquad f(\rho)=\rho^{-3}\gamma(\hbox{$3\over 2$},\rho^{2})^{2},\hfill\hbox to0.0pt{\hss{\tenrm(\hbox{\tenrm B}.7)}}\crcr}}}$$

Full text

CERN-TH-2018-259

Topological susceptibility at $T>T_{\fourteenrm c}$ from master-field

simulations of the SU(3) gauge theory

Leonardo Giustia and Martin Lüscherb,c

aDipartimento di Fisica, Università di Milano-Bicocca and

INFN, Sezione di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy

bCERN, Theoretical Physics Department, 1211 Geneva 23, Switzerland

cAlbert Einstein Center for Fundamental Physics

Institute for Theoretical Physics, Sidlerstrasse 5, 3012 Bern, Switzerland

Abstract

The topological susceptibility is computed in the ${\ninerm SU(3)}$ gauge theory at temperatures $T$ above the critical temperature $T_{\ninerm c}$ using master-field simulations of very large lattices, where the infamous topology-freezing issue is effectively bypassed. Up to $T=2.0\,T_{\ninerm 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\,{\ninerm fm}$ with a precision of 2 per mille.

1. Introduction

The temperature dependence of the topological susceptibility $\chi_{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–4]. Computations of $\chi_{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_{\tenrm c}$ performed to date [5–14] rely on some form of reweighting or the so-called integral method, where $\chi_{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 $\chi_{t}$) are included.

Master-field simulations [15] bypass the topology freezing issue by simulating lattices with four-dimensional volumes $V$ satisfying

[TABLE]

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_{\tenrm c}$ and $2.0\,T_{\tenrm 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 $\chi_{t}$ is discussed in sect. 3 and conclusions are drawn in sect. 4.

2. Theoretical framework

2.1 Lattice theory

The SU(3) Yang–Mills theory studied in this paper is set up on hyper-cubic $L_{0}\times 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.

2.2 Definition of $\chi_{t}$

Since the correlation function of the topological density

[TABLE]

(where $F_{\mu\nu}$ denotes the field strength of the gauge potential) has a non-integrable short-distance singularity, the topological susceptibility is only formally given by

[TABLE]

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 ${\tenrm 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 $\chi_{t}$, provided the latter is defined consistently with the chiral Ward identities. When this condition is met, $\chi_{t}$ is unambiguously determined and can be shown to be given by a singularity-free expectation value of “density chains” [21–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].

2.3 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\,{\tenrm fm}$, where the effective gauge coupling is small and the semi-classical approximation becomes asymptotically exact††$\dagger$ In the case of a four-dimensional spherical space-time, $\chi_{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].. If $L$ is much larger than the correlation lengths in the pseudo-scalar sector, $\chi_{t}$ is independent of $L$ up to exponentially small terms. This regime 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

[TABLE]

of the distribution of the topological charge $Q$ is much smaller than $1$. It is plausible that $\chi_{t}$ is dominated by the sectors with charge $Q=\pm 1$ in this case. Moreover, if their contribution is assumed to be suppressed by the factor $\exp\{-S_{\tenrm min}\}$, $S_{\tenrm min}$ being the minimum of the gauge action in these sectors (the instanton action), the renormalization group implies that

[TABLE]

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.

Table 1. Lattice parameters and simulation results for $t_{0}/a^{2}$

[TABLE]

2.4 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 $\chi_{t}$ has mass dimension $4$, its value must also be expressed in such units. In view of the steep temperature dependence of $\chi_{t}$, a relative numerical error in the reference scale thus results in an approximately $11$ times larger error of the converted values of $\chi_{t}$.

The target statistical precision of $\chi_{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.

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 $\beta=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].

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 $\beta=6/g_{0}^{2}$ considered, the lattice spacing decreases from about $0.10$ to $0.023\,{\tenrm fm}$. The lattice sizes $L$ are at least $6\,{\tenrm fm}$ and reach values above $9\,{\tenrm fm}$ in some cases. On all these lattices, $\chi_{t}V$ is in the thousands and frozen-topology effects are therefore expected to be neglible. The numbers $N_{\tenrm 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 $\beta$ and can be well represented by a polynomial

[TABLE]

of degree $4$. A least-squares fit yields the values

[TABLE]

for the coefficients. The fit approximates $t_{0}/a^{2}$ in the range $5.96\leq\beta\leq 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 $\beta=6.42$ and the fit also reproduces the values at $\beta=6.3,6.4,\ldots,7.0$ quoted in ref. [31] within errors varying from $0.2$ to $1.1$ percent.

2.5 Conversion to physical units

The ${\tenrm 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\,{\tenrm fm}$. In the range $5.96\leq\beta\leq 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\,{\tenrm fm}$ thus amounts to setting

[TABLE]

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).

Table 2. Parameters of the high-temperature lattices

[TABLE]

3. 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.

3.1 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_{\tenrm c}$ in the SU(3) gauge theory is $294\,{\tenrm MeV}$ [32] and the chosen temperatures $T$ are thus about $1.5\,T_{\tenrm c}$ and $2.0\,T_{\tenrm 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.

Fig. 2. Normalized histograms of ${\ninerm Re}\,\{zP({\ninemib x})\}$ at flow time $t=0.35\,t_{0}$ measured on the $A_{1}$ (left) and $B_{1}$ (right) lattices. In both cases, the bin size is $1/60$ and the phase factor $z\in\{1,\exp(\pm i2\pi/3)\}$ is chosen so as to cancel the phase of the average value of the Polyakov loop.

At high temperatures, the Polyakov loop

[TABLE]

(where $W(x)$ denotes the Wilson line that passes through $x$ and wraps around space-time in the time direction) assumes a non-zero expectation value. The expectation value breaks the ${{Z}}_{3}$ center symmetry of the theory and its phase is spontaneously chosen to be [math], $2\pi/3$ or $-2\pi/3$. A technically attractive choice of order parameter 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.

Fig. 3. Values of $\chi_{t}(R)$ obtained on the $A_{3}$ lattice at two flow times corresponding to smoothing ranges equal to $0.28\,{\ninerm fm}$ (squares) and $0.47\,{\ninerm fm}$ (diamonds).

3.2 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 in units of $t_{0}$ were chosen corresponding to smoothing ranges $\sqrt{8t}$ [18] approximately equal to $0.28\,{\tenrm fm}$ and $0.47\,{\tenrm fm}$.

As explained in ref. [15], $\chi_{t}$ can be obtained in master-field simulations by integrating the two-point correlation function of the charge density,

[TABLE]

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\geq 1.2\,{\tenrm fm}$ are characteristic features of $\chi_{t}(R)$ on all lattices listed in table 2. At large $T$, small $R$ and small flow times $t$, $\chi_{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 $\chi_{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 $\chi_{t}(R)$ at $R\simeq 1.4\,{\tenrm fm}$, where the asymptotic plateaus are, in all cases, safely reached within errors.

Table 3. Simulation results for $\chi_{t}$

[TABLE]

3.3 Continuum limit

The calculated values of $t_{0}^{2}\chi_{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 table 3 and 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_{0}^{2}\chi_{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,

[TABLE]

being a bit more precise. These figures are orders of magnitude smaller than the susceptibility $t_{0}^{2}\chi_{t}=6.67(7)\times 10^{-4}$ [24] at zero temperature and the observed rapid decrease from $T=1.5\,T_{\tenrm c}$ to $T=2.0\,T_{\tenrm 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.

Fig. 4. Extrapolation of the values of $t_{0}^{2}\chi_{t}\times 10^{5}$ listed in table 3 to the continuum limit (left: $A$-lattices, right: $B$-lattices). The data at flow time $0.35\,t_{0}$ (squares) and $0.97\,t_{0}$ (diamonds) are extrapolated linearly in $a^{2}$, the grey points at $a=0$ being the extrapolated values.

3.4 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_{\mu\nu}(x)$ is given by the so-called clover formula. A classically ${\tenrm O}(a^{2})$-improved expression is then

[TABLE]

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 ${\tenrm 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 $\chi_{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\,{\tenrm 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\geq 2.4\,{\tenrm fm}$ are thus required to be safe of finite-volume effects.

4. 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_{\tenrm 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\leq 0.05\,{\tenrm 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 $\tau=2$. The molecular-dynamics equations were integrated by applying the forth-order 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_{\tenrm step}$ was required on the larger lattices in order to preserve a good acceptance rate $P_{\tenrm 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.

Table 4. Simulation parameters$\,{}^{\ast}$

[TABLE]

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 $\tau_{\tenrm 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 $\Delta\tau_{\tenrm 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]. Autocorrelations 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.

Fig. 5. Relative deviation from the ensemble average of the values of the topological susceptibility computed using single master fields (run $A_{3}$, flow time $0.97\,t_{0}$).

A.3 Use of quadruple-precision arithmetic

On the simulated lattices, significance losses of up to 11 decimal places occur when the energy deficit $\Delta 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 quadruple-precision artithmetic in the summation of the action densities over all lattice points. $\Delta 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 $\chi_{t}(R)$ in perturbation theory

In the continuum theory and at flow time $t>0$, the integrated correlation function

[TABLE]

of the topological charge density can be straightforwardly expanded in powers of the gauge coupling. The computation proceeds along the lines of ref. [18] except for the fact that the time components $p_{0}$ of the momenta are quantized in units of $2\pi T$.

At high temperatures, where

[TABLE]

the contributions of the $p_{0}\neq 0$ modes of the gauge field to the leading-order expression for the two-point function of the charge density are exponentially suppressed. Up to these terms the latter is then given by

[TABLE]

$\alpha_{s}$ being the strong coupling,

[TABLE]

the incomplete $\Gamma$-function and $r={\tenmib x}^{2}/(8t)$. The correlation function thus decreases monotonically from

[TABLE]

becomes negative at some point and eventually goes to zero with a rate proportional to $|{\tenmib x}|^{-6}$ at large distances $|{\tenmib x}|$.

Fig. 6. The function (B.7) assumes its maximal value $0.1583(1)$ at $\rho=1.236(1)$ and decays like $\rho^{-3}$ at large $\rho$.

Equation (B.3) leads to the expression

[TABLE]

for the approximate susceptibility (B.1). To this order of perturbation theory, $\chi_{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\,{\tenrm fm}$ and $0.45\,{\tenrm fm}$, while the leading-order expression (B.6) has its maximum at $R=0.35\,{\tenrm fm}$ and $0.58\,{\tenrm 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 $\chi_{t}$ to all orders.

References

[1] R. D. Peccei, H. R. Quinn, CP conservation in the presence of instantons, Phys. Rev. Lett. 38 (1977) 1440

[2] R. D. Peccei, H. R. Quinn, Constraints imposed by CP conservation in the presence of instantons, Phys. Rev. D16 (1977) 1791

[3] S. Weinberg, A new light boson?, Phys. Rev. Lett. 40 (1978) 223

[4] F. Wilczek, Problem of strong $P$ and $T$ invariance in the presence of instantons, Phys. Rev. Lett. 40 (1978) 279

[5] E. Berkowitz, M. I. Buchoff, E. Rinaldi, Lattice QCD input for axion cosmology, Phys. Rev. D 92 (2015) 034507

[6] S. Borsanyi et al., Axion cosmology, lattice QCD and the dilute instanton gas, Phys. Lett. B752 (2016) 175

[7] C. Bonati et al., Axion phenomenology and $\theta$-dependence from $N_{f}=2+1$ lattice QCD, JHEP 1603 (2016) 155

[8] P. Petreczky, H. P. Schadler, S. Sharma, The topological susceptibility in finite temperature QCD and axion cosmology, Phys. Lett. B762 (2016) 498

[9] S. Borsanyi et al., Calculation of the axion mass based on high-temperature lattice quantum chromodynamics, Nature 539 (2016) 69

[10] J. Frison, R. Kitano, H. Matsufuru, S. Mori, N. Yamada, Topological susceptibility at high temperature on the lattice, JHEP 1609 (2016) 021

[11] Y. Taniguchi, K. Kanaya, H. Suzuki, T. Umeda, Topological susceptibility in finite temperature (2+1)-flavor QCD using gradient flow, Phys. Rev. D95 (2017) 054502

[12] P. T. Jahn, G. D. Moore, D. Robaina, $\chi_{\tenrm top}(T\gg T_{\tenrm c})$ in pure-glue QCD through reweighting, Phys. Rev. D98 (2018) 054512

[13] F. Burger, E. M. Ilgenfritz, M. P. Lombardo, A. Trunin, Chiral observables and topology in hot QCD with two families of quarks, Phys. Rev. D98 (2018) 094501

[14] C. Bonati et al., Topology in full QCD at high temperature: a multicanonical approach, JHEP 1811 (2018) 170

[15] M. Lüscher, Stochastic locality and master-field simulations of very large lattices, EPJ Web Conf. 175 (2018) 01002

[16] R. Brower, S. Chandrasekharan, J.W. Negele, U.-J. Wiese, QCD at fixed topology, Phys. Lett. B560 (2003) 64

[17] S. Aoki, H. Fukaya, S. Hashimoto, T. Onogi, Finite volume QCD at fixed topological charge, Phys. Rev. D76 (2007) 054508

[18] M. Lüscher, Properties and uses of the Wilson flow in lattice QCD, JHEP 1008 (2010) 071 [Erratum: ibid. 1403 (2014) 092]

[19] M. Lüscher, P. Weisz, Perturbative analysis of the gradient flow in non-Abelian gauge theories, JHEP 1102 (2011) 051

[20] K. G. Wilson, Confinement of quarks, Phys. Rev. D10 (1974) 2445

[21] L. Giusti, G. C. Rossi, M. Testa, G. Veneziano, The $U_{\tenrm A}(1)$ problem on the lattice with Ginsparg--Wilson fermions, Nucl. Phys. B628 (2002) 234

[22] L. Giusti, G. C. Rossi, M. Testa, Topological susceptibility in full QCD with Ginsparg--Wilson fermions, Phys. Lett. B587 (2004) 157

[23] M. Lüscher, Topological effects in QCD and the problem of short-distance singularities, Phys. Lett. B593 (2004) 296

[24] M. Cè, C. Consonni, G. P. Engel, L. Giusti, Non-Gaussianities in the topological charge distribution of the SU(3) Yang--Mills theory, Phys. Rev. D92 (2015) 074502

[25] M. Lüscher, A semiclassical formula for the topological susceptibility in a finite space-time volume, Nucl. Phys. B205 [FS5] (1982) 483

[26] T. C. Kraan, P. van Baal, Periodic instantons with nontrivial holonomy, Nucl. Phys. B533 (1998) 627

[27] T. C. Kraan, P. van Baal, Monopole constituents inside SU(n) calorons, Phys. Lett. B435 (1998) 389

[28] R. Sommer, A new way to set the energy scale in lattice gauge theories and its applications to the static force and $\alpha_{s}$ in SU(2) Yang--Mills theory, Nucl. Phys. B411 (1994) 839

[29] M. Guagnelli, R. Sommer, H. Wittig (ALPHA collab.), Precision computation of a low-energy reference scale in quenched lattice QCD, Nucl. Phys. B535 (1998) 389

[30] S. Necco, R. Sommer, The $N_{f}=0$ heavy quark potential from short to intermediate distances, Nucl.Phys. B622 (2002) 328

[31] M. Asakawa et al., Determination of reference scales for Wilson gauge action from Yang--Mills gradient flow, arXiv:1503.06516

[32] G. Boyd et al., Thermodynamics of SU(3) lattice gauge theory, Nucl. Phys. B469 (1996) 419

[33] M. Lüscher, P. Weisz On-shell improved lattice gauge theories, Commun. Math. Phys. 97 (1985) 59 [E: ibid. 98 (1985) 433]

[34] A. Ramos, S. Sint Symanzik improvement of the gradient flow in lattice gauge theories, Eur. Phys. J. C76 (2016) 15

[35] S. Duane, A. D. Kennedy, B. J. Pendleton, D. Roweth, Hybrid Monte Carlo, Phys. Lett. B195 (1987) 216

[36] I. P. Omelyan, I. M. Mryglod, R. Folk, Symplectic analytically integrable decomposition algorithms: classification, derivation, and application to molecular dynamics, quantum and celestial mechanics simulations, Comp. Phys. Commun. 151 (2003) 272

[37] T. J. Dekker, A floating-point technique for extending the available precision, Numer. Math. 18 (1971) 224

[38] D. E. Knuth, Semi-Numerical Algorithms, in: The Art of Computer Programming, vol. 2, 2nd ed. (Addison-Wesley, Reading MA, 1981)

[39] J. R. Shewchuk, Adaptive precision floating-point arithmetic and fast robust geometric predicates, Discrete & Computational Geometry 18 (1997) 305