Topological susceptibility of pure gauge theory using Density of States

The topological susceptibility of the SU(3) pure gauge theory is calculated in the deconfined phase at temperatures up to $10T_c$. At such large temperatures the susceptibility is suppressed, topologically non-trivial configurations are extremely rare. Thus, direct lattice simulations are not feasible. The density of states (DoS) method is designed to simulate rare events, we present an application of the DoS method to the problem of high temperature topological susceptibility. We reconstruct the histogram of the charge sectors that one could have obtained in a naive importance sampling. Our findings are perfectly consistent with a free instanton gas.


Introduction
The past decade has witnessed an immense progress in the theoretical description of the thermodynamics of strongly interacting matter through the advances in the solution strategies of the underlying theory, Quantum Chromodynamics (QCD). New insights came from a wide range of first principle approaches ranging from resummed perturbation theory [1] through functional methods [2,3] to direct simulations on the lattice [4]. One of the remaining less understood aspects of QCD is related to the role of instantons. One example is the strong CP problem for which the Peccei-Quinn mechanism [5] offers a solution by introducing the axion particle [6,7]. The hypothetical axion is searched for in various experimental designs, e.g. by shining light through a wall [8], helioscopes [9,10] and haloscopes [11,12]. The search can be narrowed down by constraints on the axion mass, e.g by the requirement, that axions have no more contribution to the dark matter than the total observed abundance [13][14][15]. For the latter cosmological input to be effective we have to obtain information for the axion potential at the temperatures of the Early Universe, where these were produced. This strategy was pursued in the framework of lattice QCD in Ref. [16].
The QCD axion effectively couples to the gauge invariant but CP breaking combination of the strong fields q(x) = 1 32π 2 ǫ µνρσ Tr(F µν (x)F ρσ (x)) . ( This quantity is called topological charge, as its integral Q evaluates to integer numbers and this charge is topologically stable. Evaluating this term at temperature T in a space-time volume of Ω we can express the topological susceptibility as χ(T ) determines the quadratic term of the axion potential, while higher order fluctuations of the charge control the details of the shape of the potential. The determination of the topological susceptibility using lattice methods has a long history [17][18][19][20][21][22][23][24][25]. Slightly below the QCD transition temperature susceptibilities close to the zero temperature value were observed. At high temperatures, on the other hand, χ(T ) drops with an approximate power law. The power law was actually expected from the Dilute Instanton Gas Approximation (DIGA) [26]. The rapid drop of the susceptibility with the temperature manifests in the finite volume lattice simulations such that calorons (the finite temperature localized objects that carry a topological charge) are extremely rare. Lattice QCD simulations would need to sample such rare events with sufficient statistics to determine at least the variance of Q. The problem of freezing topological sectors in the lattice update algorithms poses an additional challenge. Thus, brute force approaches e.g. in Ref. [22] are naturally limited to a short temperature range.
Several ideas have been proposed to circumvent aforementioned problems. In Refs. [27,28] analytic continuation from imaginary Θ parameters was used to map out χ(T ) and other parameters of the free energy, offering a way to calculate higher moments of the topological charge. Refs. [16,29] addressed the rarity of topological configurations in the Markov chain of simulation updates. It was observed that at sufficiently high temperatures configurations with |Q| ≥ 2 practically never occur even if the volume is kept large enough to contain the g 2 T or even T c scale. Simulating in the Q = 1 and Q = 0 sectors separately and determining their relative free energy provides for an indirect method to calculate χ(T ) at high T . While this method was applied successfully, one had to rely on the cancellation of a quartic divergence in the trace anomaly. A very different, reweighting based approach was advocated by Refs. [30,31] where a modified update algorithm was introduced to enhance the production of dislocations that may grow into calorons. Ref. [32] makes a further step and includes the enhancement force into the microcanonical update rather than deferring it into a Metropolis step. A common feature of these reweighting based methods is the use of a proxy charge, which is an easily accessible non-integer function of the gauge fields that strongly correlates with the integer charge Q. Most of these methods can or have been generalized for the case of dynamical fermions [16,23].
In many lattice studies configurations with Q = ±2 are mostly missing. Thus, possible interactions between calorons are not described. While Q = 2 configurations have arguably small weight at high temperature, one may not want to exclude such interactions from the begining. In Ref. [33] the statistics of both calorons and anti-calorons were considered (as opposed to the net charge). The distribution of the topological objects was perfectly consistent with an ideal gas at a temperature as low as T ≈ 1.05 MeV.
In our work we corroborate the DIGA picture in the high temperature Yang-Mills theory. We calculate the topological susceptibility using lattice simulations, not ignoring the very rare Q = ±2 sectors. We find that the weight of the latter supports the ideal gas description.
We approach the problem of rare calorons with the Density of States (DoS) method. Originally, the DoS method studies several small energy ranges separately to determine the energy dependence of the density of states [34]. The DoS is designed to simulate the physics of rare events. A prominent example for its successful use is the solution of the sign problem in heavy dense QCD [35]. For a detaled review of the recent progress see [36] and references therein. This study aims to show for the first time that the DoS is applicable to the problem of measuring the topological susceptibility.
In a nutshell, our strategy uses a micro-canonical force on a proxy charge in the well defined framework of the DoS method. The rare events with a large proxy charge are sampled and provide for configurations with Q = ±1 and beyond. As an additional measure to reduce auto-correlation times we also deploy the parallel tempering method. We consider here pure gauge theory, expanding the method to QCD with quark degrees of freedom involves no further conceptual problems.
In Section 2 after a brief overview of the DoS method, we explain the details of the application of DoS to the pure gauge system. We introduce our lattice setup in Section 3 and present the numerical results in Section 4. Finally, conclusions are offered in Section 5.

Density of States
First we state the idea of the Density of States method generally, before we specify to the problem of interest. Given an action S[Φ] of space-time dependent fields Φ(x, t), we are interested in the partition function: Now using the fact that the value of the Gaussian integral is independent of the constant a, we can rewrite the integral expression for the partition function such that up to an irrelevant normalization factor we have where F [Φ] is an arbitrary functional of the fields. Swapping the order of the integrations we can write where we defined the 'density of states': If we choose F to be the energy functional and consider the limit P → ∞ than ρ(c) indeed describes the density of the energy states of the system, and we get the partition function as an integral over energy. The formulas remain correct for an arbitrary functional F [Φ] and also for finite P . In this case we call ρ(c) the generalized density of states. We can measure observables using the formula where we have defined the notation · · · c , which is an average with the action S c The density of states is reconstructed by measuring the derivative of its logarithm: we can thus measure ∂ c ln ρ(c) on a predetermined set of points and reconstruct ln ρ(c) (and thus ρ(c)) using numerical integration (with e.g. the trapezoid rule). Using this prescription we obtain ln ρ(c) with an error magnitude approximately independent of c, thus we get ρ with approximately constant relative errors in the whole c range. An alternative determination of ρ(c) in the P → ∞ limit is possible by simply measuring the histogram of the observable F [Φ] in a simulation with P = 0. In this case, however, the statistical errors are proportional to ρ(c), which can be prohibitive.
Using the former method, thus, allows the determination of the probability of certain rare events in the configuration space of the theory, which one could not hope to reach in a naive importance sampling simulation. A similar setup was used in [37,38] at nonzero chemical potential µ, where an additional sign problem is also present. For the present study, we stay at Θ = 0, so the theory has no sign problem. In a recent paper [39] the density of states method was used in a U (1) gauge theory with a Θ-term. There the authors used open boundary conditions to avoid quantization of the topological charge and make the system amenable to the DoS treatment. Here we take a different route, see below.
We will calculate the topological susceptibility χ = Q 2 /Ω, where Q is the topological charge and Ω is the fourvolume. At large temperatures in the deconfined phase, the topological susceptibility is known to be very small [26]. In an importance sampling simulation the theory is almost always in the zero charge sector, thus the value of the susceptibility is given by the probability of rare visits to the ±1 charge sectors. The idea of this paper is to use the density of states method to measure the probability of these unlikely visits to nonzero charge sectors. However, the topological charge of the configurations is given by integer numbers such that the density of states in this case would be a sum of delta functions on the integer values (in the P → ∞ limit), and the procedure described above is not applicable. To work around this problem, we need a proxy charge Q P [U ] (a function of the link variables U ) which is a continuous value such that Q P [U ] is close to the integer topological charge Q[U ] [32]. Recall the topological charge of a Yang-Mills theory defined in Eq. (1). On the lattice, one can e.g. choose a field theoretical definition of the charge based on the Wilson flow [40,41], similarly to the cooling techniques introduced in [42]. (For other equivalent definitions, see [43].) One evolves the gauge field configurations using the flow equations given by the Wilson plaquette action, and measures a discretised version of the field strength tensor appearing in Eq. (1). One observes that this discretised definition tends to integer numbers at large flow times, and one can carry out the continuum limit by fixing the flow time (at which the measurement of the charge is to be carried out) in physical units.
At zero flow time the gluonic definition of the charge is not close to integer numbers, and typically it can be far from the integer topological sector of the configuration. The idea is that if we define the proxy charge Q P to be the charge at small flow time, then it will not be restricted to integer numbers. The integer values are approached after a longer flow time fixed e.g. to the temperature scale. In order to be able to use a Hybrid Monte Carlo algorithm, we need to be able to calculate the derivative of Q P with respect to the gauge fields. This suggests to use the analytic stout smearing procedure [44], so we define

Simulation setup
The topological susceptibility is often normalized to the transition temperature's fourth power. On an N 3 S × N T lattice is given by The gauge action with tree-level Symanzik improvement, and the clover discretised topological charge density in Q P is used in the simulations.
As discussed in Section 2 a separate simulation must be performed for several c values, such that the appropriate range of the proxy charge is covered. Typically we have used 30-60 c values to measure ∂ c ln ρ(c) and expectation values as a function of c. As the theory is symmetric in c at Θ = 0, we only used non-negative c values. We measured the topological charge of every ∼ 50-th configuration using the gluonic definition after the Wilson flow. For the measurement of the exact topological charge Q, we use an improved discretisation for the topological charge density including the 1 × 2 plaquettes in the clover formula [45][46][47], which we evaluate at the flow time t = 1/(8T 2 ) = N 2 T a 2 /8. We round the obtained charge values to integer values, subsequently.
Unless stated otherwise, results for N T = 6, N S = 24 are presented. We use Hybrid Monte Carlo for updating the configurations, such that the force of the fixing term of the action is calculated in every 3-4th step. For the calculation of Q P we typically use n = 4 stout smearing steps with ρ = 0.1. The algorithmic parameter P required hand-tuning to P = 1000. Too small values don't constrain the dynamics enough to allow extrapolation of Q = 0 sectors, too large values lead to large force terms that require small HMC step sizes and thus slow down the simulation.
Finally ρ(c) is reconstructed by integrating ∂ c ln ρ(c). Expectation values of observables are calculated from Eq. (8) using the trapezoid rule. The undetermined overall factor of ρ(c) (which drops out of observables) can be fixed by keeping the integral ∞ −∞ ρ(c) normalized to 1. In practice we restrict the integral over positive values of c, and the integral has a cut-off at the largest c value simulated. Since ρ(c) = ρ(−c) and ρ(c) typically has an overall exponential decay for large c this is well justified. Statistical errors are calculated using the Jackknife procedure. In Fig. 1 we show the reconstructed ρ(c) function for several temperatures. One observes a roughly exponential decay for large c which gets faster as the temperature is increased. At larger temperatures one sees a second local maximum around c ∼ 0.8 − 0.9, corresponding to configurations which have topological charge Q = 1. At larger c values one can find similar peaks corresponding to the Q = 2, 3, . . . sectors.
There is a further ingredient in the used algorithm with the goal to reduce auto-correlation times. We use parallel tempering [48,49] across the several simultaneously run ensembles, each working on a different c parameter. Parallel tempering adds a further update step between the HMC trajectories. The update consists of the swap of the gauge configurations between ensembles at neighbouring points of the c-grid. The change in the total action is taken into account by a Metropolis step. The tempering update allows dislocations produced at some c to travel in the c space and enhance the variety of the configurations at any given c parameter. This comes at the price of having correlated errors on the ρ(c) curve. These correlations are correctly kept when the c-integrals are calculated.
The grid of c values for a simulation is chosen such that there is sufficient overlap between neighbouring simulations. A dense grid helps maintaining a high acceptance rate for the tempering updates and keeps the systematic error of the  integral under control. Thus, we can keep the systematic error coming from the finite c grid below the magnitude of the statistical errors. In the right panel of Fig. 1 two reconstructed ρ(c) functions are compared: "disc. 1" has twice as many grid points as "disc. 2". They differ only in the range where ρ(c) drops below 10 −12 as observed in the Figure. The corresponding χ(T )/T 4 c values are 1.00(18) · 10 −6 and 1.03(17) · 10 −6 , respectively for "disc. 1" and "disc. 2". We used the coarser grid for the result plots.
For each c ensemble we determine the Q c average, this we show in Fig. 2 for several temperatures. As expected it roughly follows the Q = c line, and it has plateaus at integer values: if c is close to an integer number, then the system stays in the topological sector picked out by c. In the region c ≈ 0.5 − 0.7 we see that the average charge smoothly goes from 0 to 1. By checking the Monte Carlo history of the topological charge in this region one can see that here the system goes through many tunnelings and the autocorrelation time of the topological charge remains small, though with the decrease of the lattice spacing the situation gradually worsens as expected.
In Fig. 3 the quantity ρ(c) Q 2 c is shown. The integral of this quantity gives the topological susceptibility (we normalize the integral of ρ(c) to 1). We see that this quantity gets more and more sharply peaked with increasing temperature, and there the |Q| > 1 sectors have negligible contribution to the topological susceptibility. Note that a logarithmic scale is used which means that even at the smallest temperature the peak of the Q = 1 sector carries by far the largest contribution to Q 2 . On the right panel we clearly see that the Q = 1 peak is located at a value c < 1, which then translates to a Q P value significantly less than one. This was expected since there is a multiplicative renormalization between Q P and Q. After performing the c integral this renormalization factor does not enter our susceptibility results.

Results
We start with the calculation of the topological susceptibilities by performing the c integrals at each temperature. In Fig. 4 we show the resulting χ(T ) function for 24 3 × 6 lattices. Measurements of the susceptibility using direct simulations (i.e.  simulations with P = 0) are also included for comparison. Direct simulations get increasingly difficult as the temperature is increased, as in that case the system is almost always in the Q = 0 sector, with rare visits to the Q = ±1 sector. The quick decay of the probability of configurations with Q = 0 makes direct measurements of the topological susceptibility above T ∼ 4T c practically impossible [22].
The density of states approach does not depend on rare tunnelings and can be used at high temperatures, though at higher temperatures one has to deal with increasing thermalization and autocorrelation times. We observe a power law dependence of the susceptibility with the exponent −6.3 (1). The presented χ(T ) as well as the exponent are affected by discretization effects, an agreement with the perturbative result [26] is expected to hold in the continuum limit only.
The continuum extrapolation is performed at one point: at T = 4.1T c , using simulations on N T = 6, 8, 10, N S = 4N T lattices, as visible in Fig. 5. The temperature is chosen to facilitate comparison with the result from [31] which reads: 24 10 −6 , where they performed the continuum extrapolation for the quantity ln χ. (Note that Refs. [30,31] used the plaquette gauge action.) Our result is χ/T 4 c = (5.5 ± 2.8)10 −6 and it's close to the result we get from continuum extrapolating the logarithm: χ/T 4 c = 6.61e ±0.30 10 −6 . The difference of the two extrapolations one can take as the systematical error of the continuum extrapolation. One can wonder whether rounding Q to integer values at a certain flow time has an influence on our results. In Fig. 5 we show also the values which one gets without the final rounding step. We observe that for decreasing lattice spacing the effect of this rounding steadily decreases. In Tab. 1 the number of configurations we used to calculate the topological charge Q for this continuum extrapolation is listed. Note that we measured the charge after every ∼50 HMC trajectories, but the quantity ∂ c ln ρ(c) is measured more often as it is a much cheaper observable.  Next we turn to the question of instanton interactions and the applicability of the dilute instanton gas approximation (DIGA), which assumes interactions are negligible. At large temperatures the appearance of a caloron is so rare that the probability that two calorons appear close to each other is small, therefore the DIGA is expected to be a good approximation. In this case it is expected that the value of Q 2 is proportional to the volume and thus the susceptibility in Eq. (12) is independent of the volume.
To investigate this behavior we have performed simulations at T = 1.5T c at two different spatial box sizes L = 4/T and L = 6/T . In Fig. 3 (right) the quantity Q 2 c ρ(c) is shown for the two volumes. Since the ρ(c) function is normalized to one, the integral of this function gives Q 2 . One sees that at LT = 4 mainly the Q = 1 peak contributes, as in the smaller volume the appearance of two calorons is relatively rare. In contrast, in the larger volume the Q = 2 sector gives a non-negligible contribution to the susceptibility. The results for the susceptibility are consistent: χ/T 4 c = 0.00716(56) on the smaller lattice and χ/T 4 c = 0.00755(48) on the larger lattice. The χ values remain volume independent, showing that neglecting calorons' interaction is indeed a good approximation and we could get fairly accurate results in the smaller volume by restricting our simulations to the 0 ≤ c ≤ 1.5 range. Because of the strong suppression of the calorons at larger temperatures the contribution from the higher sectors is even smaller.
We also study the observable b 2 , defined by which characterises the anharmonicity of the axion potential. It is expected that at large temperatures, where the DIGA approximation holds, b 2 assumes the value -1/12. In small volumes where only the Q = 0, ±1 sector contributes, the value -1/12 follows from the fact that Q 4 = Q 2 and Q 2 is small. In larger volumes, where calorons and anticalorons appear independently, their probability distribution follows the Skellam distribution: p k = e −λ I k (λ) (see also below), which leads to b 2 = −1/12 as well. Earlier results show that starting from T > 1.15T c , b 2 is very close to -1/12 [19]. In fact in all our simulations b 2 is consistent with -1/12 within errors. As argued above, this is nontrivial only in the case where Q 4 has a sizeable contribution from |Q| ≥ 2 sectors. The simulation using LT = 6, T = 1.5T c allows testing this case and we get a result compatible with the DIGA predicition: b 2 = −0.0804(58).
The histogram of the topological charge is defined with the help of the Kronecker-delta function as where the overall normalization is chosen such that n h(n) = 1. Using the DoS simulations we can reconstruct this as In Fig. 6 we show reconstructed histograms. Note that some of the probabilities are so low that it would be practically impossible to measure them using naive simulations. In the right panel of Fig. 6 we show the reconstructed histograms at T = 1.5T c for two different spatial volumes T L = 4 and T L = 6. We can model these histograms under the assumption that instantons are independent of each other: this means that the number of calorons and anticalorons follow a Poisson distribution. Thus the total topological charge is distributed as the difference of two Poisson distributed random variables, known as the Skellam distribution: p k = e −λ I k (λ), with I k (λ), the modified Bessel function of the first kind. (This distribution is equivalent to takeing the total number of calorons and anticalorons as a Poisson distributed variable and then assigning a random sign to each defect [50].) We observe that the fitted Skellam distribution describes the histograms well. This verifies that instantons appear independently of each other in our simulations. As expected the same λ parameter scaled with the volume describes the different spatial volumes, as one observes on the right panel of Fig. 6. The results of a "bruteforce" calculation with the unconstrained gauge action are also shown, using ≈ 3000 independent configurations.

Conclusions
In this study the Density of States method is applied to the SU(3) pure gauge theory in order to calculate its topological susceptibility at high temperatures. In practice, we introduce a force term on a proxy charge, which may assume noninteger values, but correlates with the integer topological charge. This allows the mapping out of the proxy charge density using DoS. The topological susceptibility is calculated in a temperature range T /T c = 1.2 . . . 10. The DoS method also allows reconstructing the histograms of the topological sectors one would get in a naive importance sampling simulation. We observe that these histograms follow the Skellam distribution, which implies that instanton interactions are negligible at the spatial volumes used in this study. This is additionally verified by performing simulations at two different spatial volumes at T = 1.5T c , such that on the larger volume the |Q| = 2 charge sector has a non-negligible contribution to the susceptibility and to the b 2 parameter.
In this exploratory study we mostly use N T = 6 ensembles, except for one test at T = 4.1T c , where the continuum extrapolation is carried out. The continuum extrapolation over all temperatures is to be carried out in a follow-up study, allowing a quantitative comparison with perturbative results.
Using the method presented here the topological susceptibility at individual temperatures can be addressed as opposed to the integral method in Refs. [16,29]. Also the absence of large cancellations between the subtracted free energy contributions may open the way towards simulations at finer lattices.
For eventual applicability to axion phenomenology, fermionic degrees of freedom are also required. Performing a continuum extrapolation with dynamical fermions can be highly non-trivial and often very fine lattices are required. Nevertheless, since the modified dynamics affects the gauge sector only, the inclusion of fermions presents no conceptual challenge to the DoS method described here.