Instanton liquid properties from lattice QCD

We examined the instanton contribution to the QCD configurations generated from lattice QCD for $N_F=0$, $N_F=2+1$ and $N_F=2+1+1$ dynamical quark flavors from two different and complementary approaches. First via the use of Gradient flow, we computed instanton liquid properties using an algorithm to localize instantons in the gauge field configurations and studied their evolution with flow time. Then, the analysis of the running at low momenta of gluon Green's functions serves as an independent confirmation of the instanton density which can also be derived without the use of the Gradient flow.


Introduction
Instantons and other semi-classical solutions of the QCD Lagrangian are believed to play an essential role in the low energy dynamics of QCD, where the crucial phenomena of confinement and chiral symmetry breaking take place [1]. Instantons are by definition topologically non-trivial solutions of the classical field equations in Euclidean space with finite action. These solutions play an interesting role from Quantum Mechanics, where they describe tunneling processes, in Minkowski spacetime, from one vacuum to another vacuum all the way to Yang-Mills theories where they describe tunnelling processes between different vacua which are labeled by a different value of the topological charge (winding number). Additionally, they play a crucial role in the explanation of the mechanism of spontaneous breaking of chiral symmetry [2] but their relation to confinement in four dimensional theories is much less clear despite the initial success in three-dimensional gauge theories [3]. Apart from tunnelling paths in Minkowski spacetime, instantons are intimately related to many interesting phenomena such as the U (1) problem, the strong CP problem, they have interesting counterparts in the electroweak sector (sphalerons) that can lead to violation of the baryon and lepton number conservation and could possibly describe rare processes of baryon decay. Also, in supersymmetric gauge theories the exact β-function can be computed by instanton methods. For more details we refer the interesting reader to [4]. QCD topology is a fully non-perturbative topic, and naturally the best way to study it systematically is in ab-initio lattice simulations. Strictly speaking, topology is not mathematically well defined on the lattice. One can employ different definitions of the topological charge either gluonic (which are based on some discretization of the topological charge density after some smoothing procedure) or fermionic (based on the spectrum of the Dirac operator). However, for small values of the lattice spacing were the gauge fields are smooth enough, it is still meaningful to study the topology of gauge fields and there are intensive explorations in lattice simulations [5][6][7][8]. While not demonstrated theoretically, it is expected and also confirmed by numerical simulations that all definitions of the topological charge agree in the continuum limit (of course can be affected by different sizes of cutoff effects) [8].
The main reason that a smoothing procedure is needed when one employs a gluonc definition of the topological charge is that both the topological charge and the instanton contribution to the lattice gauge fields are "hidden" by the presence of short-range (UV) fluctuations and most studies reveal the presence of instantons only after the application of such a filtering technique. Cooling [9], different smearing methods and, more recently, Gradient flow are efficient techniques that remove short-range fluctuations in the gauge fields that have been widely exploited for the study of QCD topology [6].
After UV fluctuations have been removed, different methods have been used to recognize instantons in the remaining gauge fields [10][11][12] and measure their density and size distributions. The application of a filtering technique may nevertheless introduce biases in the characteristics of the underlying semiclassical configuration. For example, instanton/antiinstanton pair annihilation would lead to a reduction of the instanton density that would therefore depend on the amount of UV filtering applied. The fact that the instanton size may be modified by the filtering is also known, at least with the use of cooling, in such a way that depending on the gauge action used, instantons shrink or grow with cooling. A third phenomenon, the disappearance of small instantons of size comparable to the lattice spacing may introduce uncontrolled and worse effects because it would affect the value of the topological charge.
A different approach to the determination of the instanton nature of QCD vacuum was presented in [13], where the IR running of some gluon Green's functions was asserted to be related to some properties of the instanton ensemble. Much effort has been devoted by the non-perturbative QCD community to the understanding of the deep IR running of the correlation functions among the fundamental fields of the theory. Quite remarkably, the combination of results from lattice, Dyson-Schwinger equations and some other continuum approaches has led to the firm conclusion that the gluon propagator acquires a dynamical mass in the deep IR while the ghost propagator remains massless [14][15][16][17][18][19][20][21][22][23][24][25][26][27][28]. An appealing possibility to describe the large distance (low momenta) correlations among gluon fields is the use of an instanton liquid model with the advantage that it can be applied both before and after the removal of UV fluctuations.
Our proposal in this paper is to apply both the direct recognition of instantons and the analysis of the running of gluon Green's functions, using the comparison between both methods in order to quantify the systematic uncertainties present in each method. We will furthermore analyze the instanton properties before applying the Gradient flow from i) an extrapolation of the results to zero flow time and ii) the analysis of the IR running without Gradient flow. Overall we expect to obtain a precise image of the instanton description of the QCD vacuum, and its evolution with Gradient flow. Using this combination of methods we will explore the dependence of the instanton liquid parameters on the number of dynamical fermions of the theory, by exploiting quenched ensembles (N F = 0) and unquenched ensembles (with N F = 2 + 1 and N F = 2 + 1 + 1 dynamical flavors).
The action density and topological charge density obtained after removing UV fluctuations show the presence of large lumps that are thought to be the realization of the IR structure of the initial gauge fields. One possible description of those structures are instantons. Section 2 will be devoted to a brief description of the Ansätze commonly used for the instantonic description of the QCD vacuum. A fit of the topological charge density obtained from the lattice to the instanton prediction has been used to quantify not only the topological charge but the instanton density as well (see section 3.3).
Although instantons are believed to be the basis of long-distance physics of the QCD vacuum, few studies have searched instanton traces without the application of any filtering technique. A few years ago, the IR running of the gluon Green's functions was proposed as a way to study the instantonic content of lattice gauge field configurations without filtering short-range fluctuations [29][30][31].
The aim of the paper is twofold. The fist purpose of the paper is therefore to combine the methods based in the analysis of the IR running of gluon Green's functions and the localization of instantons by fitting the lumps after application of the Gradient flow. The second motivation is to describe the parameters of the instanton liquid model before the application of any filtering technique. This article starts with the definition of the multiinstanton background which is introduced in Section 2, along with all the Landau gauge Green's functions formulae of the instanton liquid model. In Section 3 the local instanton recognition method is introduced along with all the basics of the Gradient flow. In Section 4 we present a detailed analysis of the IR running of the MOM coupling constant, with and without the effects of the Gradient flow. In the same Section the fate of instantons, under the Gradient flow, is studied in a simple toy model. Section 5 summarizes our conclusions.

Multi-instanton background
Our main purpose is to make the semiclassical multi-instanton background resulting from the non-trivial QCD vacuum manifest, basically by computing its contribution to the gluon gauge field, assuming it to be dominant in the deep infrared domain, and identifying then this contribution from results obtained from lattice simulations. In order to do so, we can proceed in two different manners: (i) a direct recognition of local structures in configuration space either for the gauge action or for the topological charge density, due to the instanton brackground, and (ii) by the scrutiny of the low-momenta gluon correlations expressed by the appropriate Green's functions. Let us then start by introducing a practical description of the gauge field within the classical multi-instanton background.

The gauge field within a multi-instanton ensemble
We base our description of the QCD vacuum in a reliable approximation of a multi-instanton solution for the SU (3) gauge action, built on the ground of the appropriate Ansatz for the minimization of the resulting action per particle. The following trial function, has been proposed in ref. [32] for the gauge-field classical solution from an ensemble of instantons, B a µ ; where y i = (x − z i ) and η α µν is the 't Hooft tensor, that should be replaced by η α µν when summing over anti-instantons as i = A. R aα (i) represents the color rotations embedding the canonical SU (2) instanton solution in the SU (3) gauge group (i.e., α = 1, 2, 3 and a = 1, 2, . . . 8). f (x) is a shape function, to be obtained by the minimization of the action, that should obey f (0) = 1 in order not to spoil the field topology at the instanton centers and which additionally provides sufficient cut-off at large distances guaranteeing convergence of the sum. The length parameter ρ is known as the instanton radius or size and, of course, remains transparent for a classical solution of the field equations. It can be only fixed, phenomenologically, after a successful comparison of the instanton-based prediction of a given observable with its empirical result, or related to the lattice spacing when describing lattice results. In terms of this length scale, regimes of large and small distances can be asymptotically defined from Eq. (2.1) and, as discussed in [29], in both cases, the gauge field can be effectively described by the following independent-pseudoparticle Ansatz [33] 2) which appears to be a linear superposition of pseudo-instantons at different positions, and where the profile function φ behaves as where f (x) is the shape function defined in Eq. (2.1). As will be seen below, a linear superposition makes possible a simple and closed result for the gluon two-and three-point correlation functions, reliable both in the large as in the small-distance domains. The profile function φ appears to match the behaviors from both domains, expressed by the shape function, and breaks explicitly the scale independence required by the (small-distance) limit of an isolated instanton, provided by the BPST solution [34] φ BPST (z) = 1 1 + z 2 .
(2.4) Thus, assuming a BPST profile near any instanton center, one can therein approximate the topological charge density by (2.5) with "+" sign for instantons and "−" sign for anti-instantons.
On the other hand, while the one-instanton contributions dominate at small distances (near any instanton center) over the non-linear effects from other instantons in the background, the latter are dominant for large distances and their average define the drop of the shape function which, according to [33], can be approximated as being also independent of the instanton size.
It is worth emphasizing that, even when using a profile function φ(z) that takes into account instanton interactions, the Ansätze given by Eqs. (2.1) and (2.2) are only approximations to the solutions for the classical equations of motion. 1 Later below, aiming at a careful examination of the classical background underlying the gauge fields produced by lattice QCD simulations, we will apply a prescription based on the local recognition of the peaks induced in the topological charge density by the pseudo-instantons of Eq. (2.2). Such a goal requires that both Ansätze provide a good description of the gauge field locally around each pseudo-instanton center, such that the topological charge density results fairly approximated by Eq. (2.5). Furthermore, a non-local prescription for the scrutiny of the classical background in the gauge field, based on the analysis of two-and three-point gluon Green's functions, has been recently investigated in [29] and will be extensively used in the following. For this latter method to work, Eq. (2.1) needs to be a good approximative representation for the classical gauge field within the multi-instanton background, and the deviations in Eq. (2.2) with respect to Eq. (2.1) are to be neglected from the gluon correlations in momentum space, at least in the low-momentum regime. Its main advantage is however that, as will be later explained, it can be directly applied to the lattice gauge fields without any filtering to deprive the fields from quantum fluctuations; being thus a cheaper and less distorting prescription for the examination of the quasi-classical background.

Green's functions in Landau gauge
Let us assume the independent pseudo-particle Ansatz in (2.2) for the gauge field within the multi-instanton ensemble, i.e. an instanton liquid model (ILM) for the gauge field. Then, if we assume that instanton positions and color orientations are uncorrelated, the gluon propagator dressing function reads [31,36,37] in Landau gauge (always in Euclidean space), with √ V B a µ (k) = d 4 ke ik·x B a µ (x), V being formally the volume of spacetime and where the usual average over the gauge group space, · · · , is recast by the ILM as the average over all the possible configurations of pseudoinstantons within the statistical ensemble defining the classical background of the QCD vacuum. All the dependence on the profile function φ(z) is captured by where J 2 is a Bessel function of the first kind. Analogously, albeit requiring more algebra, the scalar form factor for the tree-level component of the symmetric Landau-gauge threegluon Green's function, the kinematical configuration of momenta defined by p 2 = q 2 = r 2 = k 2 .
A quite remarkable feature of the strong coupling defined in the so-called symmetric MOM scheme [38], computed from the form factors (2.6) and (2.8), happens to be its independence on the instanton profile in both the limits of small and large momenta [13,29], as a result of whereρ = ρ 2 expresses the mean instanton size and δρ 2 = (ρ −ρ) 2 stands for the mean square width of the instanton size distribution. Then, at large and small momenta, the contribution from the classical instanton background to the running of the coupling defined by Eq. (2.11) remains insensitive to both the instanton profile function and to the statistical distribution of the multi-instanton ensemble represented by the independent pseudo-particle Ansatz in (2.2). Indeed, it will only keep track of the instanton density, as immediately results from plugging (2.12) into Eq. (2.11), one thus obtaining that the coupling approximately behaves according to a k 4 -power law wherein the classical background appears to be dominant. This was clearly shown in Ref. [13] to happen at low momenta, where the ultraviolet (UV) quantum fluctuations have little effect on the gluon correlations, by a careful examination of the coupling computed from lattice QCD. The same also happens, after filtering the short-range UV fluctuations out from the gauge fields, at large momenta [29,31,39]. In particular in [29], we applied the so-called Gradient flow as filtering technique and thus harvested a solid confirmation of the picture emerging from Eqs. (2.11,2.12). This is especially so, because the effect of the size-distribution width yields a different factor in front of the power for low and large momenta (see Eq. (2.12) above), and this is precisely what the analysis made in [29] appeared to show. We will supplement here this non-local technique, in some representative cases, with a local one based on the shape recognition of peaks in the topological charge density which, albeit more expensive from the point of view of computing time, will be of great help for pinpointing the instanton sizes and obtaining a complete picture of the instantonic description.

Instanton local recognition from lattice QCD
As explained in the previous subsection, one of the main objectives of this paper is to examine the lattice gauge fields in configuration space, in order to reveal their underlying multi-instanton content by applying the local recognition recipe that will be introduced and discussed below. A consistent comparison of the properties so obtained for the quasiclassical multi-instanton ensemble with those extracted with the Green's function method will be finally presented.

The topological charge density
As already advanced, the examination of the classical background underlying the lattice gauge fields in lattice configurations is based on the recognition of the shape of the peaks, locally around their centers, induced in the topological charge density by the pseudoparticles employed to build the multi-instanton Ansantz for the gluon gauge field, as given by Eq. (2.2). The first step needs thus to be the computation of the topological charge density, which is obtained, as well as the gauge action, as follows.
Let us call Π µν (x) the Hermitean traceless part of the plaquette starting at site x in the µ − ν plane, with µν (x) the average of the four plaquettes in the µ − ν plane. Then, the action density can be obtained as where the 8π 2 normalization factor comes from the action of a single instanton, such that d 4 xS(x) = 1 either for the instanton or anti-instanton cases. While, for the topological charge density one is left with which, again, keeps the appropriate normalization factor to result in d 4 xQ(x) = +1(−1) for instantons (anti-instantons). Eq. (3.3) is the so-called clover definition of the topological charge density, which is expressed by a combination of the plaquettes given in Eq. (3.1). More elaborated definitions, which include rectangular Wilson loops of size 1 × 2 in order to reduce discretization errors, can be found in literature (see, e.g. [8]). Notwithstanding this, our use of Eq. (3.3), supplemented with the prescription and acceptance criteria that will be described below, as well as with the comparison with the Green's function results, will be enough to achieve our purpose of revealing the quasi-classical multi-instanton background underlying the gauge fields and making its properties manifest. However, the quantum UV fluctuations are obviously present in any ensemble of lattice gauge fields obtained by Monte Carlo techniques and implementing the QCD action, therefore hiding the quasi-classical background. In computing the gluon correlations, the UV fluctuations have been shown to be negligible in the low-momentum regime, basically dominated by the long-distance correlations [13,29,31]. On the contrary, gauge fields are locally dominated by short-range contributions, which, if they are not previously suppressed, would spoil the topological charge density defined by (3.3). In the following subsection, we will briefly explain the use of the Gradient flow to remove such short-range contributions, as it has been also done in [29], in that case to unveil the quasi-classical k 4 -power law also at large-momenta.

Gradient flow
The Gradient flow can be conceived as a smoothing procedure which diminishes the shortdistance fluctuations which, within the context of a quantum field theory, correspond to UV quantum fluctuations. Therefore, depriving the gauge fields from them, potentially, implies to isolate the underlying non-trivial classical solutions which minimize the gauge action. Despite of sharing the capacity of eliminating the short-distance fluctuations with other filtering techniques such as cooling or smearing, the Gradient flow has a theoretical ground with attractive features such as the simple renormalization of the flown fields [40].
In continuum language, the Gradient flow field B µ (τ, x) of SU (N ) gauge field results from the solution of the following first order differential equation that introduces the flow time τ and the field strength tensor for the flown fields and covariant derivative, respectively. As initial condition for the differential equation, one chooses B µ (0, x) = A µ (x), so that the fundamental gauge field A µ (x) evolves from τ = 0 and, flown after any τ , reads as a formal expansion that, at tree-level, makes apparent an effective suppression of shortdistance fluctuations up to a distance of √ 8τ . Its lattice formulation, the SU (N ) matrices are flown by obeying a discretized counterpart of the first-order differential equation (3.5), with the initial condition V µ (0, x) = U µ (x), and where S is a discretization of the gauge action (herein, we will mainly use the Wilson gauge action for this, and will only make use of the Iwasaki gauge action to perform some tests, as discussed below), g 0 is the bare coupling and ∂ x,µ are the link derivatives that can be found defined in [41]. Despite its important advantages as it is possessing a well defined continuum limit or the existence, uniqueness and smoothness of its solutions, the effect of Gradient flow over the gauge fields is in the practical level equivalent to that of cooling, at least for the purpose of dealing with the topological charge [5,6,8]. A remarkable property of the Gradient flow is however the steadiness of any exact solution of D µ G µν = 0 under the flow, thereupon implying that quantum UV fluctuations can be filtered out without distorting exact solutions such as an isolated instanton. In practice, this property of the Gradient flow does not prevent it from modifying the parameters of the multi-instanton ensemble (2.2), representing the quasi-classical background underlying the lattice gauge fields. And the latter is true for two different reasons: first because the whole argument about steadiness corresponds to the continuum formulation of both instantons and flow, and second because the multi-instanton representations are only approximated solutions to the classical equations of motion, the isolated instantons never being an adequate description of the QCD vacumm, where a rather dense distribution is expected [12]. For comparative purposes, the flow time should be expressed in units of t 0 , defined by √ 8t 0 = 0.3fm following [41]. This allows to settle physical units for the Gradient flow times, given that t 0 = a 2 τ 0 = 0.01125fm 2 and t = τ t 0 /τ 0 .

Locating instantons
Let us start by describing the lattice simulations set-ups for the gauge field configurations both quenched (N F = 0) and unquenched (with N F = 2 + 1 and N F = 2 + 1 + 1 dynamical flavors) that will be either exploited in this section by applying the instanton local recognition, for some representative cases, or subsequently via the scrutiny of the low-momenta running of the coupling defined by Eq. (2.11).
For the quenched data, we consider two simulations exploited in [29] and [42] with lattice spacings and volumes appropriate for having both enough lattice sites in the vicinity of an instanton peak and enough data for momenta within the IR window, after Fourier transform. In the case of unquenched simulations, we will use two sets of lattice data, the first includes N F = 2 + 1 field configurations produced by the RBC/UKQCD collaboration using domain wall fermions very close to the physical pion mass (139 MeV), more details on the lattice set-up can be found in [43]; the second set corresponds to a N F = 2 + 1 + 1 simulation from the ETM collaboration with a pion mass of around 300 MeV. This last simulation has been carried out using Wilson twisted mass fermions at maximal twist (details of the lattice setup can be found in [44]). Details such as the gauge action, the values of the lattice spacing, the spacetime volumes or the total number of exploited configurations for each set-up appear gathered in Tab. 1.
As will be seen below, we will take the configurations for the tree-level Symanzik gauge action at β=4.2 (see Tab. 1) as a representative case and apply the Gradient flow for τ = 2, 3, 4, 5, 8, 10 and 15. Only for these fields, we will apply the instanton local recognition, while the Green's function method will be applied for all of the lattice ensembles displayed in Tab. 1. It should be particularly noticed that, by the combination of different gauge actions in the quenched case, we intend to guarantee that the discretization effects are under control or, at least, that the observed behavior appears irrespectively of the action considered. Finally, having access to N F = 0, N F = 2 + 1 and N F = 2 + 1 + 1 dynamical flavor configurations we will be in a position to check the dependence on the instanton liquid parameters with the number of fermionic flavors in the sea. Now, after applying the Gradient flow, the short distance fluctuations that are present in both the action and topological charge densities are suppressed and smooth lumps, presumably induced by the pseudo-particles in Eq. (2.2), appear unveiled. We will then explicitly check for those lumps in the topological charge density, locally around their centers, their shape similarity to the classical instanton solution in Eq. (2.5). In the rest of this section we will discuss the algorithm we have used for identifying instantons in the flown field configurations. The first step is to localize the extrema of the topological charge density that will be candidates for topological charges. Then, an estimate for their sizes has to be provided. The aim of the last step will be to eliminate candidates that are not topological charge structures but remain UV fluctuations.
We will define local maxima (minima) of the topological charge density as those sites x 0 where Q(x 0 ) is larger (smaller) than the closest 8 neighbors. We are aware that defining the extrema with respect to the 8 sites that are at ±a away in one direction may produce double counting when a single instanton is so much distorted by quantum fluctuations that another extremum takes place close the real one in the surrounding hypercube. Had we used a more strict definition of an extremum, such as Q(x) being larger than the 80 sites of the surrounding hypercube (see discussion in [45]), we would have avoided the double counting but would have been also left to disregard the distorted real instanton. Therefore, we have preferred to use a more modest definition of an extremum and apply a further filter to eliminate false candidates at the last step.
After localizing these instanton candidates, the next step will be to obtain their instanton size ρ. To this aim, we consider the ratio of the topological charge at the closest 8 neighbors to the value at the center that, according to Eq.(2.5) should be and fitted the lattice data to this expression using ρ as a free parameter. According to Eq. (2.5), the value of the topological charge at the center of a single instanton is related to the instanton size by Q(x 0 ) = 6 π 2 ρ 4 . To check this dependence, the peak value of the topological charge density has been plotted versus the fitted size in Fig. 1 for three typical quenched configurations after τ =4, 8 and 15 Gradient flow times. In this plot there is a point for each local extremum of the topological charge density, while the continuum line corresponds to Eq. (2.5). As can be seen, most of the extrema lie near the continuum line although a non-negligible number of outliers have to be filtered out as described below. Furthermore, the larger is the flow time, the more closely concentrated around the continuum line tend the extrema to be.
Then, as a first filtering prescription, we proceed first by checking that the topological charge and action densities for all the candidates behave as Eq. (2.2) for short distances, and reject those not satisfying and where R and Q are parameters that we fix in our algorithm. We have checked that the total number of instantons found is rather stable when varying them, and thus fixed R = 0.5 and Q =0.3 for the results shown below. Once we ensure that the functional form of the topological charge lump behaves locally as an isolated BPST instanton, and that it is selfdual, we are left with the above raised issue of preventing from the double counting of two extrema corresponding to a single but distorted instanton. To do so, we have finally included a filter for eliminating close pairs that works as follows, when an extremum is identified as an instanton of size ρ i at site x i , no other candidate is accepted if the site of the extremum lies within the hypersphere defined by x i + ρ i u, with u a unitary vector in physical units (and, in all the cases, the hypersphere radius is fixed to 2 in lattice units when the instanton size is smaller). The parameter is varied between 0.7 and 1, the difference between the densities obtained for those two values being used to estimate the systematic uncertainty associated to our localization method. This filter serves to exclude the possibility of finding one instanton in the core of another one 2 . The filtering of close pairs eliminates a number of candidates much larger than the functional form of Q(x), but many of them correspond to the outlier extrema of Fig. 1 that anyhow passed the criterion (3.9), namely those candidates for which instanton size ρ i and Q(x i ) are not well consistent with (2.5). Both parameters of the instanton ensemble, the average instanton radius and density, which are describing the semiclassical background underlying the gauge field, do depend on the Gradient flow time and, in general, evolve with any parameter controlling the suppression of the quantum fluctuations via a given smoothing prescription minimizing the action. Although this is well known [5,12,[46][47][48], few authors have studied this evolution in detail, quoting just the values of the instanton density or size distribution after a fixed flow time, cooling or smearing. In our case, after applying the Gradient flow and our instanton locating algorithm, the emerging picture is that of a very dense ensemble of instantons with a size of around 1/3 fm, which is more diluted as Gradient flow eliminates instantons and anti-instantons, while the remaining ones become larger with the flow. This is very apparent from the results displayed in Fig. 2 and also in Tab. 2. Furthermore, we have examined the correlations among instantons, both of like (II or AA) and unlike charges (IA), through the evaluation of the radial distribution function of pairs, g(r). This function is defined as the probability density, evaluated at a distance r = |x − x i | away from any given instanton located at the position x i , of finding one like-or unlike-charge instanton divided by the density of like-or unlike-charge instantons. Isotropy and translational invariance guarantee that, after averaging over all the possible configurations of the instanton ensemble, the distribution is indeed a radial function irrespective of the position of the instanton x i . At large distances and for the homogeneous distribution of an instanton gas or liquid, correlations among instantons tend to vanish and one would be thereupon left with g(r) → 1. This radial distribution g(r) can also be understood as the ratio of the number of instantons found in a spherical shell of infinitesimal width dr and a radius r away from each instanton divided by the shell volume and rescaled by the instanton density. Therefore, for a lattice evaluation, as the lattice volume is proportional to the number of lattice sites, one can in practice estimate the radial distribution of pairs by counting the number of like-or unlike-charge instantons at a given distance away from each instanton, perform then the average over all the instantons of the ensemble, divide next by the product of the total number of avalaible lattice sites at the same distance and the instanton density; and, finally, perform the average over the ensemble of lattice configurations. So proceeding, we have been left with the results displayed in Fig. 3 for several Gradient flow times. A first feature to be noticed from the plots is that there is an excluded volume for small distances which grows moderately with the flow time. This might be readily thought to be a consequence of applying the double-counting-preventing prescription that filters out the lumps lying within the core of an accepted instanton. How-ever, the prescription is only applied for like-charge instantons, while the same feature is also clearly manifest for those of unlike charge. Thus, we can conclude that the picture of "hard-core" instantons is fairly consistent with the results we obtain for the correlations of pairs. Furthermore, the increasing of the excluding volume reflects consistently the growing of the instanton radius which is manifest from Fig. 2. A second remarkable feature is that the distribution of instantons is not completely random: after Gradient flow, it is more likely to find a charge of the same sign than a particle of opposite sign at short distances. The latter is consistently accounted for the fact that the annihilation of opposite charge pairs during the flow increases the probability of finding instantons of the same charge close to each other.
A short comment about the effect of the filtering for small distances introduced in the localization algorithm is furthermore in order. Had we taken a smaller cut on the instanton minimum inter-distance, a large peak at short distances would have emerged for alike charges in Fig. 2, signaling the presence of very close pairs. The question of whether those instanton clusters are physical or an artifact of the localization algorithm was discussed long ago [10,47], but up to now there is no definite answer. Thus, our estimates for the instanton density with this algorithm might underestimate systematically the actual density of topological objects because of this. Later on, we will return to this point and justify, on the basis of the comparison of the densities obtained by instanton localization and by Green's functions, that this peak for close pairs ought to be indeed neglected as an artifact of the instanton localization algorithm. However, this implies that one needs, to admit a systematic deviation of around 20-30 %.
The instanton densities we have measured are of the order of n ∼ 1fm −4 for rather large flow times, a value that agrees with the one that can be set from the gluon condensate G 2 , and that has served as reference for decades [1]. Nevertheless, if we try to infer the density at zero Gradient flow time from the results in Fig.2, the density results much larger. Although it is difficult to extrapolate back to τ = 0 from our results, it points towards an order of magnitude larger (see next section).
Concerning the instanton size, Garcia-Perez et al. [46] found that when cooling with the Wilson action, individual instantons should shrink under a cooling procedure while with overimproved actions they stabilize or grow with the application of cooling. A similar discussion appears in the nice review by Creutz [49] who also analyzes the effects of these types of action modifications to reflection positivity [50,51]. Although it is not fully known how the use of different filtering techniques may modify the conclusions reached by the aforementioned authors, the fact of being in a dense liquid instead of being isolated is definitely a dramatic change compared to the setup of the previous studies. Indeed we found that most of the instantons we localized grow with the application of the Gradient flow, while only small instantons shrink (and eventually disappear). The latter appears to indicate that, at a first stage in the evolution with the flow time, a sort of effective interaction between the independent pseudo-instantons (used to describe the classical solution) dominates the process. This interaction takes into account the non-linear effects between the pseudoparticles of the Ansatz which minimize the total action and, at least when the density is large, tends to favor both the annihilation of opposite-charge instantons and their growing  when the density drops. We have checked that, if the Gradient flow is driven by the Iwasaki gauge action, the observed evolution does not differ from the one described above.

Momentum running, flow evolution and flavor effects
As we have discussed in Sec. 2, the running of two-and three-gluon Green's functions can be combined to yield a strong coupling definition, the estimate of which, computed from the gauge fields for a semiclassical multi-instanton ensemble leaves us with a very striking signature for the presence of instantons. In the following, we will use this coupling obtained from the lattice gauge fields to detect instantons by analyzing its running both at largeand small-momenta, as explained in Ref. [29], for all the set-up's of Tab. 1. We will then compare the obtained results with the ones stemming from instanton localization for the quenched case at β = 4.2 derived in the previous section, and will also study their evolution with the Gradient flow time. Furthermore, we will take advantage of the other lattice set-ups to analyze the evolution of the instanton ensemble parameters with the number of dynamical flavors.

Momentum running and the instanton density
Without the application of any filtering technique, the running at large momenta of the strong coupling defined by Eq. (2.11) is dominated by the perturbative prescription, giving rise to the appearance of asymptotic freedom, but the deep non-perturbative region (typically below ∼ 1GeV) exhibits a behavior that fits well to Eqs. (2.11,2.12). This finding serves as a confirmation that the low-momenta (large-distance) correlations are dominated by instanton-like objects, and allows to extract the instanton density without the need of any filtering technique. Moreover, after applying the Gradient flow, the UV fluctuations are filtered out and the gluon correlations become again dominated by the semiclassical multi-instanton ensemble and the large-momenta running of the strong coupling is also well accounted by Eq. (2.11). This can be clearly seen in Fig. 4 for the quenched case at β = 4.2, where the lattice estimates appear to be manifestly consistent with Eq. (2.11), after applying the Gradient flow, in both the large-and small-momenta domains; thus supporting a multi-instanton ensemble picture for the QCD vacuum at least when describing gluon correlations. It is indeed worthwhile to realize how the suppression of UV fluctuations from the gauge fields turns the well-known large-momenta logarithmic behavior of gluon correlations (brown solid circles in Fig. 4, depicting the results from non-flown fields) into the k 4 -law expressed by Eq. (2.11); law which becomes more and more apparent when the Gradient flow time evolves from τ =4 (red circles) to τ =15 (black circles). The same figure had been also presented in [29], where its main features were very well understood, in particular how the width for the instanton size distribution shifts the intercept of the logarithmic line up at low momenta, as dicated by (2.12).   Table 3. Comparisons of instanton densities obtained by applying the localization technique (third column) and the Green's function running method (fourth column), after the Gradient flow for three representative times. The errors here quoted incorporate, and are basically dominated by, some systematic uncertainties in the IR case (see text), but only statistical for the GF method.
Now, as the slope of α MOM (k) in k 4 , according to Eqs.(2.11,2.12), allows for a direct estimate of the instanton density from the lattice data at large momenta, we consider three representative flow times (τ =4,8 and 15) for the quenched simulation at β=4.2 and make a systematic comparison of the so obtained densities to the ones resulting from the direct instanton counting using the localization algorithm described in the previous section. The results from both methods lie in the same ballpark, as can be seen in Tab. 3 and Fig. 5. It appears however, that there are systematic deviations which, taking into account the statistical uncertainties, are of the order of 2-3 σ's. The deviations can be seen as an underestimation in the counting of instantons by localization of around a 30 % at τ =4 and of 24 % at τ =8 and 15. The latter can be consistently understood if one admits that some of the close pairs filtered out by the short-distance criterion correspond to real classical solutions. We can estimate that, had we accepted a 20-30 % of the rejected close pairs, the densities obtained by instanton location would have fairly agreed with the ones resulting from the analysis of the Green's function method. As a fully general classical solution describing such pairs of close instantons is not available, a systematic criterion to disentangle them from single distorted instantons that need to be rejected is not at hand. We have thus preferred to keep the short-distance criterion and assumed a systematic deviation which, anyhow, will not prevent the estimates from being in the same ballpark as those from the Green's function method. Moreover, either through the analysis of the direct counting of instantons or by the study of the running of α MOM (k), we reach a coherent image of an instanton density that drops with flow time. It is usually accepted that the results at a small fixed flow time represent the physics of the original gauge field configuration. Nevertheless, the fact that there is such a strong dependence on the flow time requires a deeper understanding. The same phenomenon is found with other filtering techniques such as cooling and we refer the reader to [52] for a simple approach for instanton density evolution through IA annihilation. Let us try here to build a toy model for this phenomenon and use it to extrapolate the instanton density down to zero flow time.

Gradient flow evolution of the instanton properties
If pair annihilation is the only phenomenon responsible for the observed density dropping, the rates at which both n I and n A decrease have to be equal. Under the hypothesis that this is a first order process, the time variation of instanton and anti-instanton densities will be given by Furthermore, if n I ≈ n A ≈ n/2 (i.e., assuming that the topological charge of a particular gauge configuration is much smaller than the instanton density), we arrive to The parameter λ may depend on the flow time τ via the instanton size, interparticle distance, etc. thus being not constant. In a first approximation, however, we assume that it is well described by a constant, and integrating out Eq. (4.2) results into This equation has been used to fit the measured instanton densities at different flow times with the instanton location algorithm for the quenched configurations at β = 4.2, and extrapolate to zero flow time. The result of the fit has been plotted in Fig. 5, thus obtaining an extrapolated density of 12.3(4)fm −4 . It is remarkable, despite of the simplicity of the toy model, how well the measured densities at different flow times agree with Eq. (4.3). Furthermore, albeit the systematic uncertainty associated to this extrapolation may be rather large and has not been quantified at this level, the so obtained density compares fairly well with the estimate from the Green's function method applied to non-flown gauge fields. The density, in this case, can be obtained once the instanton size distribution width is known and plugged into Eqs. (2.11,2.12). This width can be always estimated, after applying the Gradient flow, by the direct scrutiny of the size distribution for the localized instantons (see Fig. 2). However, as discussed in [29], it can be evaluated too, though indirectly, from the shift of the low-momentum line with respect to the large-momentum one in the logarithmic plot of Fig. 4, according to Eq. (2.12), thus providing with a width estimation fully consistently made within the Green's function approach. Owing to its consistency, we have preferred the latter (the estimates from both methods being anyhow in the same ballpark; namely, δρ 2 /ρ 2 around 0.02). Furthermore, whilst an extrapolation to zero flow time would be required, in both cases, the values extracted at low flow time remain very stable and we made the simple choice of keeping that at τ =2. Thus, a ratio δρ 2 /ρ 2 ≈ 0.014 is estimated and, applied to Eqs. (2.11,2.12), results in n ≈ 12(2)fm −4 . As can be also seen in Fig. 5, the four values for the density from Green's functions can be also well described by Eq. (4.3) but with a parameter λ slightly diminished.

Evolution with the number of flavors
We finally focus on the extensive analysis for the small-momenta behavior of α MOM (k) for the different number of flavors in Tab. 1. The running of α MOM (k) obtained from the lattice data in Tab. 1, has been plotted in Fig. 6 for N F = 0 (leftmost upper panel), N F = 2 + 1 (rightmost upper) and N F = 2 + 1 + 1 (lower) without Gradient flow. The same instanton picture that, so far, we have firmly grounded on the basis of our deep analysis of quenched lattice simulations ought to work also for the unquenched simulations and, indeed, the running of α MOM (k) between 0.3 and 0.9 GeV fits well to the power law given by Eq. (2.12) (straight line in the plots). When the lattice volume is large enough and there are lattice points at lower momenta (in particular in the quenched β=5.8 case), this instanton prediction breaks down, presumably due to the neighborhood of a zero crossing for the three gluon vertex [42,[53][54][55][56][57][58][59][60]. From the point of view of the ILM, the failure for small momenta is expected because for very large distances (typically larger than the instanton size) the model of uncorrelated instantons is not expected to work; in any case the purely quantum nature of the zero-crossing discussed in Refs. [42,59] seems to be in contradiction with an instantonic explanation. Now, as Eq. (2.11) tells, the IR running of α MOM (k) is basically determined by the instanton density and, thus, a direct comparison of the results of the low-momenta fits shown in Fig. 6 makes possible the scrutiny of the evolution of this density with the number of flavors. Indeed, it suggests that the instanton density increases monotonously with the number of dynamical flavors. If we furthermore accept that the ratio δρ 2 /ρ 2 is the same for N f = 0, N f = 2 + 1 and N f = 2 + 1 + 1, we obtain that the instanton density increases by a factor 1.3(1) from N f = 0 to N f = 2+1 and by a factor 1.5(1) from N f = 0 to N f = 2+1+1. These two ratios are nearly compatible within the errors and this suggests a very mild effect of the charm quark on the instanton density, certainly owing to its sizeable mass threshold. The role of the charm quark and the dependence on the quark masses has not been studied in detail yet but the fact that the instanton density grows with the presence of dynamical flavors has been previously reported, for instance, in Ref. [11]. In the previous reference it was interpreted as follows: it is well known that an isolated instanton produces a zero mode of the fermion determinant [11,48] and, as a consequence of this and the action of a given gauge field being weighted by [det(/ D + m)] N F , the existence of zero modes of the Dirac operator would be highly improbable in the chiral limit. Becoming more improbable the larger is the number of light fermion flavors. Thereupon, a large density of instantons can be understood as favoring the instanton superposition and thus suppressing the effect of isolated instantons. Or, in other words, it can be thought to enhance non-linear effects and pseudo-particle interactions destroying the zero-modes associated to the single instanton solution. However, one needs to be careful with the above argument (that for the massless theory, one can ignore any instanton effects since those configurations don't contribute to the partition function) since while it is correct that instantons do drop out of the partition function itself they do survive in physical correlation functions. We refer the reader to [61] for a detailed analysis of the topic and its relation to the 't Hooft vertex [62]. On the other hand, even admitting that the above qualitative argument might work, there is another source of explanation, at least partially, for this increasing with the number of dynamical flavors: any direct comparison of the instantons densities obtained from simulations with different number of flavors gets over the systematic deviations from the physical scale setting on phenomelogical estimates or in gauge quantities as the gluon correlation functions [63,64]. In setting the scale by imposing that a lattice estimate of a given quantity takes its physical value, the latter always incorporates effects from those dynamical quarks being active in the real world for the empirical determination of this quantity. In particular, the impact of this effect on a quenched theory is certainly non-negligible, but still for N f =2 simulations, there has been a recent work [65] finding a reduction of around a 5 % in the lattice spacing as a systematic deviation from the physical scale setting. A reduction of the lattice spacing ought also to be expected for the quenched lattice theory, as the lattice estimates for Λ MS appear to be lower at N f =0 [66] than at N f =2 or N f =2+1+1 [67,68], while the matching at the quark-mass thresholds 3 of theories with different number of dynamical flavors suggest the opposite trend [69]. Though it is very hard to quantify how this reduction amounts, one can readily conclude that a diminishing of a 10 % for the lattice spacing in the quenched theory would enhance by around a 40 % the instanton density and would thus make it lie on the same ballpark as N f =2+1 and N f =2+1+1. There might appear to exist also an impact from the charm quark in the lattice scale setting for the comparison of N f =2+1 and N f =2+1+1 results, however the scale setting involving quantities in the pion sector, such an impact should be nearly negligible [65].
Finally, one must also keep in mind that the estimates for the densities we compare have been obtained by invoking Eqs. (2.11) and (2.12), the last of which also needs the instanton size width to be plugged into. We have assumed, for the comparison to be made, that δρ 2 /ρ 2 is the same for any number of flavors. Whilst the latter seems to result, at low Gradient flow time, from the analysis of Ref. [29] and despite the fact that the width is only playing at the level of a subleading correction, any small deviation would become magnified by the factor 48 in (2.12) and so produce a sizeable impact on the comparison. However, as above stated, the picture we obtain for the impact of the dynamical flavors on the instanton density is consistent with the results of Ref. [11], obtained with a different technique for the instanton detection.

Conclusions
We have presented an analysis of the instanton contributions to the QCD vacuum using two complementary approaches. The first one is based on the fact that an independent pseudoparticle approach (i.e. an instanton liquid model without correlations among instantons) predicts a power law k 4 for the combination of gluon Green's functions used to define the coupling in the symmetric MOM scheme, α MOM (k). Moreover, the coefficient of k 4 can be used to infer the instanton density. The second approach uses the Gradient flow technique to remove short range fluctuations in the gauge field configurations and instantonlike structures are revealed in the topological charge density. We have proposed a method for identifying instantons among the extrema of the topological charge density and applied it to some quenched configurations.
Although GF is thought to modify minimally the topological structure of the gauge field configuration, it does not prevent the disappearance of instantons, neither of small instantons whose size is comparable to the lattice spacing, nor through instanton/anti-instanton annihilation. Indeed, we have observed that both methods (IR running of α MOM (k) and instanton localization) show a strong decrease in the instanton density with the flow time. Whilst the densities obtained from both methods differ by a ∼ 25 − 30 %, they can be seen to lie in the same ballpark and, on top of this, the fact that both evolve in the same manner with flow time strongly supports the instanton dominance after removing the short range fluctuations. The existence of noticeable discrepancies between the instanton densities obtained from both methods can be well interpreted as due to the systematic uncertainties related to the difficulties for identifying close instanton pairs (or clusters) with our localization algorithm owing to strong deviations locally around the instanton centers from the BPST profile used.
A simple model for instanton annihilation reproduces, qualitatively and quantitatively, the evolution of the instanton with GF, and allows to extrapolate the instanton density to zero flow time. The density obtained, ∼ 12fm −4 , is much larger than the traditionally quoted value of 1fm −4 , although some modern estimates [12] point towards larger densities. For zero flow time, i.e. without the application of the Gradient flow, we cannot localize instantons using our direct method, but can still fit the instanton density using the k 4 power law. We obtained that the slope of α MOM (k) supports the picture of a dense instanton liquid, with a density of n ≈ 12.3(2)fm −4 . This density is fully compatible with the extrapolation of the densities obtained from the localization method using the annihilation model.
After performing the comparison of both methods in the quenched case, we obtained the instanton density for unquenched (N F = 2 + 1 and N F = 2 + 1 + 1) lattice gauge field configurations without GF. From the fit to the power-law, we obtained larger instanton densities for the unquenched case, which although already reported [11] is a somehow controversial finding. The naive argument associates, via the Atiyah-Singer index theorem, instantons to zero modes of the Dirac operator and, thus, the presence of light dynamical quarks in the unquenched simulation should suppress instantons. Indeed our numerical findings with both methods employed seem to indicate the opposite behavior, namely the density of instantons grows so that they get distorted and can not any longer be associated to zero modes of Dirac operator. Two considerations can be made concerning this enhancement of the instanton density with the inclusion of dynamical quarks: first only the light quarks should have an effect and, therefore, the inclusion of the charm quark when passing from N F = 2 + 1 to N F = 2 + 1 + 1 is thought to have a negligible impact over the instanton density. We found ratios to the quenched case of 1.3 (1) for N F = 2 + 1 and 1.5(1) for N F = 2 + 1 + 1, nearly compatible within errors. The second consideration discussed in the text is the impact of the light quarks in the quenched lattice spacing setting. When using any observable to fix the lattice spacing, it incorporates the effect of any quark being active at the scale used for the determination of the physical observable. Therefore, in the comparison of quenched to unquenched results, it is conceivable that a systematic deviation of the lattice spacing is present. For our estimates of the instanton density, a ∼ 10% systematic deviation for the quenched lattice spacing would account for the differences between N F = 0 and N F = 2 + 1 and 2 + 1 + 1.
In summary, we feel that the results presented here draw a clear image of instanton dominance after application of gradient flow, when there is no trace of Λ QCD , and the running of gluon Green's functions is dominated by the semiclasical content of the gauge fields. Perhaps more interestingly, the IR behavior of gluon Green's functions without gradient flow seems to be well described by an instanton liquid model with a rather large density. Why instantons may serve for describing some properties of low energy QCD and not others, or to what extent instantons are relevant for the QCD phenomenology seem to be still open problems.