Topological susceptibility of pure gauge theory using Density of States
aa r X i v : . [ h e p - l a t ] F e b Topological susceptibility of pure gauge theory using Density of States
Szabolcs Bors´anyi and D´enes Sexty Department of Physics, Wuppertal University, Gaussstr. 20, D-42119 Wuppertal, Germany J¨ulich Supercomputing Centre, Forschungszentrum J¨ulich, D-52425 J¨ulich, Germany University of Graz, Institute for Physics, A-8010 Graz, Austria
February 26, 2021
Abstract
The topological susceptibility of the SU(3) pure gauge theory is calculated in the deconfined phase at temperaturesup to 10 T c . At such large temperatures the susceptibility is suppressed, topologically non-trivial configurations areextremely rare. Thus, direct lattice simulations are not feasible. The density of states (DoS) method is designed tosimulate rare events, we present an application of the DoS method to the problem of high temperature topologicalsusceptibility. We reconstruct the histogram of the charge sectors that one could have obtained in a naive importancesampling. Our findings are perfectly consistent with a free instanton gas. The past decade has witnessed an immense progress in the theoretical description of the thermodynamics of stronglyinteracting 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 ofQCD 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 experimentaldesigns, e.g. by shining light through a wall [8], helioscopes [9, 10] and haloscopes [11, 12]. The search can be narroweddown by constraints on the axion mass, e.g by the requirement, that axions have no more contribution to the dark matterthan the total observed abundance [13–15]. For the latter cosmological input to be effective we have to obtain informationfor the axion potential at the temperatures of the Early Universe, where these were produced. This strategy was pursuedin 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 ) = 132 π ǫ µνρσ Tr( F µν ( x ) F ρσ ( x )) . (1)This quantity is called topological charge, as its integral Q evaluates to integer numbers and this charge is topologicallystable. Evaluating this term at temperature T in a space-time volume of Ω we can express the topological susceptibilityas χ ( T ) = Z d x h q ( x ) q (0) i T, Θ=0 = lim Ω →∞ h Q i T, Θ=0 Ω . (2) χ ( T ) determines the quadratic term of the axion potential, while higher order fluctuations of the charge control the detailsof the shape of the potential.The determination of the topological susceptibility using lattice methods has a long history [17–25]. Slightly below theQCD 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 DiluteInstanton Gas Approximation (DIGA) [26]. The rapid drop of the susceptibility with the temperature manifests in thefinite volume lattice simulations such that calorons (the finite temperature localized objects that carry a topologicalcharge) are extremely rare. Lattice QCD simulations would need to sample such rare events with sufficient statistics todetermine at least the variance of Q . The problem of freezing topological sectors in the lattice update algorithms posesan 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 fromimaginary Θ parameters was used to map out χ ( T ) and other parameters of the free energy, offering a way to calculate1igher moments of the topological charge. Refs. [16, 29] addressed the rarity of topological configurations in the Markovchain of simulation updates. It was observed that at sufficiently high temperatures configurations with | Q | ≥ g 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 ) athigh T . While this method was applied successfully, one had to rely on the cancellation of a quartic divergence in the traceanomaly. A very different, reweighting based approach was advocated by Refs. [30,31] where a modified update algorithmwas introduced to enhance the production of dislocations that may grow into calorons. Ref. [32] makes a further stepand includes the enhancement force into the microcanonical update rather than deferring it into a Metropolis step. Acommon feature of these reweighting based methods is the use of a proxy charge, which is an easily accessible non-integerfunction of the gauge fields that strongly correlates with the integer charge Q . Most of these methods can or have beengeneralized for the case of dynamical fermions [16, 23].In many lattice studies configurations with Q = ± Q = 2 configurations have arguably small weight at high temperature, one may not want toexclude 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 ata temperature as low as T ≈ .
05 MeV.In our work we corroborate the DIGA picture in the high temperature Yang-Mills theory. We calculate the topologicalsusceptibility using lattice simulations, not ignoring the very rare Q = ± Q = ± First we state the idea of the Density of States method generally, before we specify to the problem of interest. Given anaction S [Φ] of space-time dependent fields Φ( x, t ), we are interested in the partition function: Z = Z D Φ e − S [Φ] . (3)Now using the fact that the value of the Gaussian integral Z ∞−∞ dc e − P ( c − a ) = r πP (4)is independent of the constant a , we can rewrite the integral expression for the partition function such that up to anirrelevant normalization factor we have Z = Z D Φ Z ∞−∞ dc e − P ( c − F [Φ]) e − S [Φ] , (5)where F [Φ] is an arbitrary functional of the fields. Swapping the order of the integrations we can write Z = Z ∞−∞ dc ρ ( c ) (6)where we defined the ’density of states’: ρ ( c ) = Z D Φ e − S [Φ] − P ( c − F [Φ]) . (7)2f we choose F to be the energy functional and consider the limit P → ∞ than ρ ( c ) indeed describes the density of theenergy states of the system, and we get the partition function as an integral over energy. The formulas remain correctfor 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 h A i = 1 Z Z ∞−∞ dc Z D Φ A [Φ] e − S [Φ] − P ( c − F [Φ]) = R ∞−∞ dc ρ ( c ) h A i c R ∞−∞ dc ρ ( c ) (8)where we have defined the notation h· · · i c , which is an average with the action S c [Φ] = S [Φ] + P ( c − F [Φ]) : h A i c = R D Φ e − S c [Φ] A [Φ] R D Φ e − S c [Φ] = 1 ρ ( c ) Z D Φ e − S c [Φ] A [Φ] (9)The density of states is reconstructed by measuring the derivative of its logarithm: ∂ ln ρ ( c ) ∂c = 1 ρ ( c ) Z D Φ e − S [Φ] − P ( c − F [Φ]) ( − P ( c − F [Φ])) = h− P ( c − F [Φ]) i c , (10)we can thus measure ∂ c ln ρ ( c ) on a predetermined set of points and reconstruct ln ρ ( c ) (and thus ρ ( c )) using numericalintegration (with e.g. the trapezoid rule). Using this prescription we obtain ln ρ ( c ) with an error magnitude approximatelyindependent of c , thus we get ρ with approximately constant relative errors in the whole c range. An alternativedetermination of ρ ( c ) in the P → ∞ limit is possible by simply measuring the histogram of the observable F [Φ] ina simulation with P = 0. In this case, however, the statistical errors are proportional to p ρ ( c ), which can be prohibitive.Using the former method, thus, allows the determination of the probability of certain rare events in the configurationspace of the theory, which one could not hope to reach in a naive importance sampling simulation. A similar setup wasused 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 usedin a U (1) gauge theory with a Θ-term. There the authors used open boundary conditions to avoid quantization of thetopological charge and make the system amenable to the DoS treatment. Here we take a different route, see below.We will calculate the topological susceptibility χ = h Q i / Ω, where Q is the topological charge and Ω is the four-volume. At large temperatures in the deconfined phase, the topological susceptibility is known to be very small [26]. Inan importance sampling simulation the theory is almost always in the zero charge sector, thus the value of the susceptibilityis given by the probability of rare visits to the ± P → ∞ limit), and the procedure described above is not applicable. To work around thisproblem, 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], similarlyto the cooling techniques introduced in [42]. (For other equivalent definitions, see [43].) One evolves the gauge fieldconfigurations using the flow equations given by the Wilson plaquette action, and measures a discretised version of thefield strength tensor appearing in Eq. (1). One observes that this discretised definition tends to integer numbers at largeflow times, and one can carry out the continuum limit by fixing the flow time (at which the measurement of the chargeis 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 fromthe integer topological sector of the configuration. The idea is that if we define the proxy charge Q P to be the chargeat small flow time, then it will not be restricted to integer numbers. The integer values are approached after a longerflow time fixed e.g. to the temperature scale. In order to be able to use a Hybrid Monte Carlo algorithm, we need to beable to calculate the derivative of Q P with respect to the gauge fields. This suggests to use the analytic stout smearingprocedure [44], so we define Q P [ U ] = Q clov [ U ′ n,ρ ] , (11)where Q clov is the clover discretisation of the topological charge (1) and U ′ n,ρ are the stout smeared link variables using n smearing steps with stepsize ρ . A similar use of the proxy charge was described in Ref. [32]. We thus use F [ U ] = Q P [ U ]in the following to constrain the action. 3 -14 -12 -10 -8 -6 -4 -2
0 0.5 1 1.5 2 2.5 3 3.5 ρ ( c ) c T/T c =1.5T/T c =2 T/T c =3 T/T c =4 T/T c =6 T/T c =10 10 -16 -14 -12 -10 -8 -6 -4 -2
0 0.2 0.4 0.6 0.8 1 1.2 1.4N T =6, T=6T c P=1000, n=4, ρ =0.1 ρ ( c ) c disc. 1disc. 2 Figure 1: The density of states for various temperatures (left), and at T = 6 T c with two different discretisations in c (right). The topological susceptibility is often normalized to the transition temperature’s fourth power. On an N S × N T latticeis given by χ ( T ) T c = h Q i Ω T c = h Q i ( N S /N T ) (cid:18) TT c (cid:19) . (12)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 appropriaterange of the proxy charge is covered. Typically we have used 30-60 c values to measure ∂ c ln ρ ( c ) and expectation valuesas a function of c . As the theory is symmetric in c at Θ = 0, we only used non-negative c values, except for a test at T = 3 T c where we observed good agreement of the results of a simulation using c ∈ [0 , .
5] and an other, independent oneusing values c ∈ [ − . , .
2] (10 χ/T c = 1 . . ∼ Q , we use an improved discretisation for the topological charge density including the 1 × t = 1 / (8 T ) = N T a /
8. We round the obtained charge values tointeger values, subsequently.Unless stated otherwise, results for N T = 6 , N S = 24 are presented. We use Hybrid Monte Carlo for updating theconfigurations, such that the force of the fixing term of the action is calculated in every 3-4th step. For the calculationof Q P we typically use n = 4 stout smearing steps with ρ = 0 .
1. The algorithmic parameter P required hand-tuningto P = 1000. Too small values don’t constrain the dynamics enough to allow extrapolation of Q = 0 sectors, too largevalues 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 keepingthe integral R ∞−∞ ρ ( c ) normalized to 1. In practice we restrict the integral over positive values of c , and the integral has acut-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 thetemperature is increased. At larger temperatures one sees a second local maximum around c ∼ . − .
9, correspondingto configurations which have topological charge Q = 1. At larger c values one can find similar peaks corresponding to the Q = 2 , , . . . sectors.There is a further ingredient in the used algorithm with the goal to reduce auto-correlation times. We use paralleltempering [48, 49] across the several simultaneously run ensembles, each working on a different c parameter. Paralleltempering adds a further update step between the HMC trajectories. The update consists of the swap of the gaugeconfigurations between ensembles at neighbouring points of the c -grid. The change in the total action is taken intoaccount by a Metropolis step. The tempering update allows dislocations produced at some c to travel in the c space andenhance the variety of the configurations at any given c parameter. This comes at the price of having correlated errorson the ρ ( c ) curve. These correlations are correctly kept when the c -integrals are calculated.4 〈 Q 〉 c cT/T c =1.5T/T c =2 T/T c =3 T/T c =4 T/T c =6 T/T c =10 Figure 2: The average of the topological charge Q as a function of c . -14 -12 -10 -8 -6 -4 -2
0 0.5 1 1.5 2 2.5 3 3.5 ρ ( c ) 〈 Q 〉 c c T/T c =1.5T/T c =2 T/T c =3 T/T c =4 T/T c =6 T/T c =10 -5x10 -7 -7 -6 -6 -6 -6
0 0.5 1 1.5 2 2.5 3N T =6, T=1.5T c 〈 Q 〉 c ρ ( c ) / Ω c LT=6LT=4 Figure 3: The average of the topological charge squared times the density of states ρ ( c ) as a function of c . The integralof this function gives the topological susceptibility. On the left panel we show the logarithm of the integrand for varioustemperatures with LT = 4, on the right panel we fix T = 1 . T c and vary the volume in a linear plot.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 theintegral under control. Thus, we can keep the systematic error coming from the finite c grid below the magnitude of thestatistical errors. In the right panel of Fig. 1 two reconstructed ρ ( c ) functions are compared: “disc. 1” has twice as manygrid points as “disc. 2”. They differ only in the range where ρ ( c ) drops below 10 − as observed in the Figure. Thecorresponding χ ( T ) /T c values are 1 . · − and 1 . · − , respectively for “disc. 1” and “disc. 2”. We usedthe coarser grid for the result plots.For each c ensemble we determine the h Q i c average, this we show in Fig. 2 for several temperatures. As expected itroughly follows the Q = c line, and it has plateaus at integer values: if c is close to an integer number, then the systemstays in the topological sector picked out by c . In the region c ≈ . − . ρ ( c ) h Q i c is shown. The integral of this quantity gives the topological susceptibility (wenormalize the integral of ρ ( c ) to 1). We see that this quantity gets more and more sharply peaked with increasingtemperature, and there the | Q | > Q = 1 sector carries by farthe largest contribution to h Q i . On the right panel we clearly see that the Q = 1 peak is located at a value c <
1, whichthen translates to a Q P value significantly less than one. This was expected since there is a multiplicative renormalizationbetween Q P and Q . After performing the c integral this renormalization factor does not enter our susceptibility results.5 e-071e-061e-050.00010.0010.010.1 1 2 4 8 10 χ / T c T/T c brute forceDoS Figure 4: The topological susceptibility measured with the brute force method and with the DoS approach. We show apower law fit with a slope of -6.3(1) with dashed line. N T c vals. c range c parameter as well as the statistics used at each c value to calculate the continuumextrapolation in Fig. 5. We start with the calculation of the topological susceptibilities by performing the c integrals at each temperature. In Fig. 4we show the resulting χ ( T ) function for 24 × P = 0) are also included for comparison. Direct simulations get increasingly difficult as the temperatureis increased, as in that case the system is almost always in the Q = 0 sector, with rare visits to the Q = ± Q = 0 makes direct measurements of the topological susceptibilityabove T ∼ T c practically impossible [22].The density of states approach does not depend on rare tunnelings and can be used at high temperatures, though athigher temperatures one has to deal with increasing thermalization and autocorrelation times. We observe a power lawdependence of the susceptibility with the exponent − . χ ( T ) as well as the exponent are affected bydiscretization 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 . T c , using simulations on N T = 6 , , , N S = 4 N T lattices, as visible in Fig. 5. The temperature is chosen to facilitate comparison with the result from [31] which reads: χ ( T = 4 . T c ) /T c = 4 . e ± . − , where they performed the continuum extrapolation for the quantity ln χ . (Note thatRefs. [30, 31] used the plaquette gauge action.) Our result is χ/T c = (5 . ± . − and it’s close to the result we getfrom continuum extrapolating the logarithm: χ/T c = 6 . e ± . − . The difference of the two extrapolations one cantake as the systematical error of the continuum extrapolation. One can wonder whether rounding Q to integer values ata certain flow time has an influence on our results. In Fig. 5 we show also the values which one gets without the finalrounding step. We observe that for decreasing lattice spacing the effect of this rounding steadily decreases. In Tab. 1 thenumber of configurations we used to calculate the topological charge Q for this continuum extrapolation is listed. Notethat we measured the charge after every ∼
50 HMC trajectories, but the quantity ∂ c ln ρ ( c ) is measured more often as itis 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 thatthe probability that two calorons appear close to each other is small, therefore the DIGA is expected to be a goodapproximation. In this case it is expected that the value of h Q i is proportional to the volume and thus the susceptibilityin Eq. (12) is independent of the volume.To investigate this behavior we have performed simulations at T = 1 . T c at two different spatial box sizes L = 4 /T and L = 6 /T . In Fig. 3 (right) the quantity h Q i c ρ ( c ) is shown for the two volumes. Since the ρ ( c ) function is normalizedto one, the integral of this function gives h Q i . One sees that at LT = 4 mainly the Q = 1 peak contributes, as in the6 c , N s =4N t , Symanzik action * χ / T c T2 roundednot roundedcontinuum extrapolation Figure 5: Continuum extrapolation of the topological susceptibility at T = 4 . T c using N T = 6 , ,
10. We show resultsfrom the prescription when at flow time t = 1 / (8 T ) the values of Q are rounded to integer values or when this roundingstep is not performed. A linear fit for both sets of points is indicated. -20 -15 -10 -5
0 1 2 3 4 p r obab ili t y QT=1.5T c , LT=6T=1.5T c , LT=4T=2T c , LT=4T=3T c , LT=4T=4T c , LT=4T=6T c , LT=4T=10T c , LT=4 10 -10 -8 -6 -4 -2
0 1 2 3 4 5N T =6, T=1.5T c p r obab ili t y QLT=6Skellam( λ )LT=4Skellam( λ *(4/6) )bruteforce LT=4 Figure 6: The reconstructed histogram of the topological charge in a simulation without the forcing term in the action.The lines represent fitted Skellam distributions: p k = e − λ I k ( λ ), where the λ parameter is consistent with h Q i .smaller volume the appearance of two calorons is relatively rare. In contrast, in the larger volume the Q = 2 sector gives anon-negligible contribution to the susceptibility. The results for the susceptibility are consistent: χ/T c = 0 . χ/T c = 0 . χ values remain volume independent, showing thatneglecting calorons’ interaction is indeed a good approximation and we could get fairly accurate results in the smallervolume by restricting our simulations to the 0 ≤ c ≤ . b , defined by b = − h Q i − h Q i h Q i (13)which characterises the anharmonicity of the axion potential. It is expected that at large temperatures, where the DIGAapproximation holds, b assumes the value -1/12. In small volumes where only the Q = 0 , ± h Q i = h Q i and h Q i is small. In larger volumes, where calorons and anticaloronsappear independently, their probability distribution follows the Skellam distribution: p k = e − λ I k ( λ ) (see also below),which leads to b = − /
12 as well. Earlier results show that starting from
T > . T c , b is very close to -1/12 [19]. Infact in all our simulations b is consistent with -1/12 within errors. As argued above, this is nontrivial only in the casewhere h Q i has a sizeable contribution from | Q | ≥ LT = 6 , T = 1 . T c allows testing thiscase and we get a result compatible with the DIGA predicition: b = − . h ( n ) = h δ Q,n i , (14)where the overall normalization is chosen such that P n h ( n ) = 1. Using the DoS simulations we can reconstruct this as h ( n ) = Z dcρ ( c ) h δ Q,n i c . (15)In Fig. 6 we show reconstructed histograms. Note that some of the probabilities are so low that it would be practicallyimpossible to measure them using naive simulations. In the right panel of Fig. 6 we show the reconstructed histogramsat T = 1 . T c for two different spatial volumes T L = 4 and
T L = 6. We can model these histograms under theassumption that instantons are independent of each other: this means that the number of calorons and anticaloronsfollow a Poisson distribution. Thus the total topological charge is distributed as the difference of two Poisson distributedrandom variables, known as the Skellam distribution: p k = e − λ I k ( λ ), with I k ( λ ), the modified Bessel function of the firstkind. (This distribution is equivalent to takeing the total number of calorons and anticalorons as a Poisson distributedvariable and then assigning a random sign to each defect [50].) We observe that the fitted Skellam distribution describesthe histograms well. This verifies that instantons appear independently of each other in our simulations. As expectedthe same λ parameter scaled with the volume describes the different spatial volumes, as one observes on the right panelof Fig. 6. The results of a “bruteforce” calculation with the unconstrained gauge action are also shown, using ≈ In this study the Density of States method is applied to the SU(3) pure gauge theory in order to calculate its topologicalsusceptibility at high temperatures. In practice, we introduce a force term on a proxy charge, which may assume non-integer values, but correlates with the integer topological charge. This allows the mapping out of the proxy charge densityusing DoS.The topological susceptibility is calculated in a temperature range
T /T c = 1 . . . .
10. The DoS method also allowsreconstructing the histograms of the topological sectors one would get in a naive importance sampling simulation. Weobserve that these histograms follow the Skellam distribution, which implies that instanton interactions are negligible atthe spatial volumes used in this study. This is additionally verified by performing simulations at two different spatialvolumes at T = 1 . T c , such that on the larger volume the | Q | = 2 charge sector has a non-negligible contribution to thesusceptibility and to the b parameter.In this exploratory study we mostly use N T = 6 ensembles, except for one test at T = 4 . T c , where the continuumextrapolation 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 opposedto the integral method in Refs. [16, 29]. Also the absence of large cancellations between the subtracted free energycontributions may open the way towards simulations at finer lattices.For eventual applicability to axion phenomenology, fermionic degrees of freedom are also required. Performing acontinuum 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 conceptualchallenge to the DoS method described here. Acknowledgements
This project was partially funded by the DFG grant SFB/TR55. We acknowledge funding by the Gauss Centre forSupercomputing (GCS) for providing computer time on the JUWELS supercomputer at the J¨ulich SupercomputingCentre (JSC) under the GCS/NIC project ID HWU16. Some parts of the numerical calculations were done on the GPUcluster at the University of Wuppertal, and on the cluster at the University of Graz.
References [1] J. Ghiglieri, A. Kurkela, M. Strickland, and A. Vuorinen, “Perturbative Thermal QCD: Formalism and Applications,”Phys. Rept. (2020) 1–73, arXiv:2002.10188 [hep-ph] .[2] C. S. Fischer, “QCD at finite temperature and chemical potential from DysonSchwinger equations,”Prog. Part. Nucl. Phys. (2019) 1–60, arXiv:1810.12938 [hep-ph] .[3] F. Gao and J. M. Pawlowski, “QCD phase structure from functional methods,” Phys. Rev. D no. 3, (2020) 034027, arXiv:2002.07500 [hep-ph] .
4] J. N. Guenther, “Overview of the QCD phase diagram - Recent progress from the lattice,” arXiv:2010.15503 [hep-lat] .[5] R. Peccei and H. R. Quinn, “CP Conservation in the Presence of Instantons,” Phys. Rev. Lett. (1977) 1440–1443.[6] S. Weinberg, “A New Light Boson?,” Phys. Rev. Lett. (1978) 223–226.[7] F. Wilczek, “Problem of Strong P and T Invariance in the Presence of Instantons,” Phys. Rev. Lett. (1978) 279–282.[8] R. B¨ahre et al., “Any light particle search II —Technical Design Report,” JINST (2013) T09001, arXiv:1302.5647 [physics.ins-det] .[9] CAST
Collaboration, S. Aune et al., “CAST search for sub-eV mass solar axions with 3He buffer gas,”Phys. Rev. Lett. (2011) 261302, arXiv:1106.3919 [hep-ex] .[10] E. Armengaud et al., “Conceptual Design of the International Axion Observatory (IAXO),” JINST (2014) T05002, arXiv:1401.3233 [physics.ins-det] .[11] ADMX
Collaboration, N. Du et al., “A Search for Invisible Axion Dark Matter with the Axion Dark Matter Experiment,”Phys. Rev. Lett. no. 15, (2018) 151301, arXiv:1804.05750 [hep-ex] .[12]
MADMAX Working Group
Collaboration, A. Caldwell, G. Dvali, B. Majorovits, A. Millar, G. Raffelt, J. Redondo, O. Reimann,F. Simon, and F. Steffen, “Dielectric Haloscopes: A New Way to Detect Axion Dark Matter,”Phys. Rev. Lett. no. 9, (2017) 091801, arXiv:1611.05865 [physics.ins-det] .[13] J. Preskill, M. B. Wise, and F. Wilczek, “Cosmology of the Invisible Axion,” Phys. Lett. B (1983) 127–132.[14] L. Abbott and P. Sikivie, “A Cosmological Bound on the Invisible Axion,” Phys. Lett. B (1983) 133–136.[15] M. Dine and W. Fischler, “The Not So Harmless Axion,” Phys. Lett. B (1983) 137–141.[16] S. Borsanyi et al., “Calculation of the axion mass based on high-temperature lattice quantum chromodynamics,”Nature no. 7627, (2016) 69–71, arXiv:1606.07494 [hep-lat] .[17] B. Alles, M. D’Elia, and A. Di Giacomo, “Topological susceptibility in full QCD at zero and finite temperature,”Phys. Lett. B (2000) 139–143, arXiv:hep-lat/0004020 .[18] C. Gattringer, R. Hoffmann, and S. Schaefer, “The Topological susceptibility of SU(3) gauge theory near T(c),”Phys. Lett. B (2002) 358–362, arXiv:hep-lat/0203013 .[19] C. Bonati, M. D’Elia, H. Panagopoulos, and E. Vicari, “Change of θ Dependence in 4D SU(N) Gauge Theories Across theDeconfinement Transition,” Phys. Rev. Lett. no. 25, (2013) 252003, arXiv:1301.7640 [hep-lat] .[20] E. Berkowitz, M. I. Buchoff, and E. Rinaldi, “Lattice QCD input for axion cosmology,” Phys. Rev. D no. 3, (2015) 034507, arXiv:1505.07455 [hep-ph] .[21] R. Kitano and N. Yamada, “Topology in QCD and the axion abundance,” JHEP (2015) 136, arXiv:1506.00370 [hep-ph] .[22] S. Borsanyi, M. Dierigl, Z. Fodor, S. Katz, S. Mages, D. Nogradi, J. Redondo, A. Ringwald, and K. Szabo, “Axion cosmology, latticeQCD and the dilute instanton gas,” Phys. Lett. B (2016) 175–181, arXiv:1508.06917 [hep-lat] .[23] C. Bonati, M. D’Elia, M. Mariti, G. Martinelli, M. Mesiti, F. Negro, F. Sanfilippo, and G. Villadoro, “Axion phenomenology and θ -dependence from N f = 2 + 1 lattice QCD,” JHEP (2016) 155, arXiv:1512.06746 [hep-lat] .[24] Y. Taniguchi, K. Kanaya, H. Suzuki, and T. Umeda, “Topological susceptibility in finite temperature ( 2+1 )-flavor QCD using gradientflow,” Phys. Rev. D no. 5, (2017) 054502, arXiv:1611.02411 [hep-lat] .[25] P. Petreczky, H.-P. Schadler, and S. Sharma, “The topological susceptibility in finite temperature QCD and axion cosmology,”Phys. Lett. B (2016) 498–505, arXiv:1606.03145 [hep-lat] .[26] D. J. Gross, R. D. Pisarski, and L. G. Yaffe, “QCD and Instantons at Finite Temperature,” Rev. Mod. Phys. (1981) 43.[27] M. D’Elia and F. Negro, “Phase diagram of Yang-Mills theories in the presence of a θ term,” Phys. Rev. D no. 3, (2013) 034503, arXiv:1306.2919 [hep-lat] .[28] C. Bonati, M. D’Elia, and A. Scapellato, “ θ dependence in SU (3) Yang-Mills theory from analytic continuation,”Phys. Rev. D no. 2, (2016) 025028, arXiv:1512.01544 [hep-lat] .[29] J. Frison, R. Kitano, H. Matsufuru, S. Mori, and N. Yamada, “Topological susceptibility at high temperature on the lattice,”JHEP (2016) 021, arXiv:1606.07175 [hep-lat] .[30] P. T. Jahn, G. D. Moore, and D. Robaina, “ χ top ( T ≫ T c ) in pure-glue QCD through reweighting,”Phys. Rev. D no. 5, (2018) 054512, arXiv:1806.01162 [hep-lat] .[31] P. T. Jahn, G. D. Moore, and D. Robaina, “Improved Reweighting for QCD Topology at High Temperature,” arXiv:2002.01153 [hep-lat] .[32] C. Bonati, M. D’Elia, G. Martinelli, F. Negro, F. Sanfilippo, and A. Todaro, “Topology in full QCD at high temperature: amulticanonical approach,” JHEP (2018) 170, arXiv:1807.07954 [hep-lat] .[33] R. A. Vig and T. G. Kovacs, “Ideal topological gas in the high temperature phase of SU(3) gauge theory,” arXiv:2101.01498 [hep-lat] .[34] F. Wang and D. Landau, “Efficient, Multiple-Range Random Walk Algorithm to Calculate the Density of States,”Phys. Rev. Lett. no. 10, (2001) 2050, arXiv:cond-mat/0011174 .[35] N. Garron and K. Langfeld, “Anatomy of the sign-problem in heavy-dense QCD,” Eur. Phys. J. C no. 10, (2016) 569, arXiv:1605.02709 [hep-lat] .[36] K. Langfeld, “Density-of-states,” PoS LATTICE2016 (2017) 010, arXiv:1610.09856 [hep-lat] .[37] Z. Fodor, S. D. Katz, and C. Schmidt, “The Density of states method at non-zero chemical potential,” JHEP (2007) 121, arXiv:hep-lat/0701022 .[38] G. Endrodi, Z. Fodor, S. D. Katz, D. Sexty, K. K. Szabo, and C. Torok, “Applying constrained simulations for low temperature latticeQCD at finite baryon chemical potential,” Phys. Rev. D98 no. 7, (2018) 074508, arXiv:1807.08326 [hep-lat] .
39] C. Gattringer and O. Orasch, “Density of states approach for lattice gauge theory with a θ -term,” Nucl. Phys. B (2020) 115097, arXiv:2004.03837 [hep-lat] .[40] M. L¨uscher, “Properties and uses of the Wilson flow in lattice QCD,” JHEP (2010) 071, arXiv:1006.4518 [hep-lat] . [Erratum:JHEP 03, 092 (2014)].[41] M. Luscher and P. Weisz, “Perturbative analysis of the gradient flow in non-abelian gauge theories,” JHEP (2011) 051, arXiv:1101.0963 [hep-th] .[42] P. de Forcrand, M. Garcia Perez, and I.-O. Stamatescu, “Improved cooling algorithm for gauge theories,”Nucl. Phys. B Proc. Suppl. (1996) 777–780, arXiv:hep-lat/9509064 .[43] C. Alexandrou, A. Athenodorou, K. Cichy, A. Dromard, E. Garcia-Ramos, K. Jansen, U. Wenger, and F. Zimmermann, “Comparison oftopological charge definitions in Lattice QCD,” Eur. Phys. J. C no. 5, (2020) 424, arXiv:1708.00696 [hep-lat] .[44] C. Morningstar and M. J. Peardon, “Analytic smearing of SU(3) link variables in lattice QCD,” Phys. Rev. D69 (2004) 054501, arXiv:hep-lat/0311018 [hep-lat] .[45] G. D. Moore, “Improved Hamiltonian for Minkowski Yang-Mills theory,” Nucl. Phys. B (1996) 689–728, arXiv:hep-lat/9605001 .[46] P. de Forcrand, M. Garcia Perez, and I.-O. Stamatescu, “Topology of the SU(2) vacuum: A Lattice study using improved cooling,”Nucl. Phys. B (1997) 409–449, arXiv:hep-lat/9701012 .[47] S. O. Bilson-Thompson, D. B. Leinweber, and A. G. Williams, “Highly improved lattice field strength tensor,”Annals Phys. (2003) 1–21, arXiv:hep-lat/0203008 .[48] R. H. Swendsen and J.-S. Wang, “Replica monte carlo simulation of spin-glasses,” Physical review letters no. 21, (1986) 2607.[49] D. J. Earl and M. W. Deem, “Parallel tempering: Theory, applications, and new perspectives,” Physical Chemistry Chemical Physics no. 23, (2005) 3910–3916.[50] M. van der Meulen, D. Sexty, J. Smit, and A. Tranberg, “Chern-Simons and winding number in a tachyonic electroweak transition,”JHEP (2006) 029, arXiv:hep-ph/0511080 ..