The topological susceptibility from grand canonical simulations in the interacting instanton liquid model: zero temperature calibrations and numerical framework
aa r X i v : . [ h e p - ph ] J a n The topological susceptibility from grand canonicalsimulations in the interacting instanton liquid model:zero temperature calibrations and numerical framework.
Olivier Wantz a a Department of Applied Mathematics and Theoretical Physics, Centre for MathematicalSciences,University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, United Kingdom
Abstract
This is the first in a series of papers which ultimately aims on improving onthe present estimates on the axion mass by modelling the topological non-perturbative QCD dynamics. Axions couple to instantons and their mass isset by the topological susceptibility whose temperature dependence we estimatewith the interacting instanton liquid model (IILM). Since accurate finite tem-perature instanton calculations have problems and do not consider fluctuationsin the topological charge, we develop an improved grand canonical version ofthe IILM to study topological fluctuations in the quark gluon plasma. In thisfirst paper we will calibrate the model against the topological susceptibility atzero temperature, in the chiral regime of physical quark masses.
1. Introduction
The strong interactions at finite temperature are believed to display a num-ber of interesting, non-perturbative phenomena, among which are the confine-ment/deconfinement transition, spontaneous P and CP violation and chiralsymmetry restoration. The latter is believed to have its origin in topolog-ical fluctuations. Lattice simulations, e.g. [8, 52], and phenomenology, e.g.[14, 40, 41, 42], have shown that the chiral dynamics of QCD is well describedby instanton models.Another interesting question related to topological fluctuations is the de-termination of the axion mass. Axions couple to instantons and their mass isdirectly proportional to the topological susceptibility. The latter turns out tobe a chiral property of QCD and can thus also be expected to be well describedby the IILM. The main physical question that we want to address is the compu-tation of the topological susceptibility and the systematic effects that pertain to Email address:
[email protected] (Olivier Wantz)
Preprint submitted to Elsevier November 23, 2018 ts determination with instanton methods in the chiral regime of light, physicalquark masses.Instanton models are based on a combination of semi-classical expansion[47] and variational approach [13, 14]. Taking this variational path integral asa starting point, Shuryak investigated what has become known as the interact-ing instanton liquid model (IILM) [43, 44]. In [39] many bulk properties werecomputed and seen to be consistent with the available lattice data and phe-nomenology. Some recent studies [20, 9] corroborate the earlier results that theIILM rather accurately describes the chiral properties of QCD, i.e. that instan-tons are the dominant degrees of freedom as far as the chiral regime of QCD isconcerned. However, the IILM fails to reproduce confinement.The topological susceptibility is a key parameter of QCD and has been in-vestigated in many lattice studies. Comparatively few studies have addressedthis quantity within the IILM [46]. One reason is that, so far, the IILM is basedon a canonical ensemble. Although one can extract the topological susceptibil-ity from the canonical ensemble through the decay of correlators, e.g. see [46]for IILM and [5, 6] for lattice simulations, it is most natural to use the grandcanonical ensemble to study the topological susceptibility. Recent investigationsof the IILM in the grand canonical ensemble [18, 17] were based on canonicalsimulations and a fugacity expansion, while we will set up a grand canonicalIILM that uses grand canonical Monte Carlo simulations; see also [10] for a‘mean-field’ study of the grand-canonical ensemble and [35] for an exploratoryinvestigation of grand canonical simulations in a simpler framework.While developing the grand canonical IILM, we have found that the existingfinite temperature IILM, [45], displays an unphysical behaviour in that it doesnot allow for a thermodynamic limit. Specifically the instanton–instanton inter-action, Eq. (3.11) in [45], contains a term that decays very slowly with instantonseparation R , ln (cid:18) βR (cid:19) R →∞ −−−−→ βR , (1)and is not integrable. In their original paper the authors do discuss this long-range interaction and point out that they found the O (1 /R ) dyon–dyon be-haviour for a wide range of intermediate separations. It might well be thatthe interactions are still well described by this ansatz for the simulation boxesused in subsequent numerical investigations, e.g. [39], but for studying the largevolume behaviour it is not appropriate.We remedy this deficiency by re-deriving the interactions and setting upa numerical framework that avoids using parametric fitting formulas, such as(1); instead we will integrate the Lagrangian density for a pair of instantonsexactly, i.e. numerically. To the best of our knowledge, the exact action densityfor a pair has not been published in the literature before; we will provide itfor Harrington–Shepard calorons [25]. The explicit form allows us to performthe numerical integration in an efficient way by exploiting the symmetries ofthe integrand, and an exact analytic computation for widely separated pairsis possible because of the localised nature of the integrands. Hence, the large2eparation interactions are under very good control. In particular, the largeseparation instanton–instanton interaction at finite temperature is not givenby (1), as we will see more explicitly in the final paper of this series [54]; inthis paper, however, we will restrict ourselves to zero temperature. We hopethat this framework can be build upon to include the non-trivial holonomycalorons [28, 30, 29], which may play part in the confinement/deconfinementphase transition, and for which good fitting formulas will be even harder tocome by because of their more complicated structure. As mentioned before, inthis series we will restrict ourselves to the Harrington–Shepard calorons.In section 2 we will review the standard strategy used to derive the parti-tion function for an ensemble of background gauge fields, i.e. the semi-classicalapproximation. We will then re-derive the interactions for the so-called ratioansatz, used to construct multi-instanton backgrounds from individual instan-tons, in section 3 and compare it with other available ans¨atze. In section 4 wepresent the numerical framework we have set up to deal with the simulations.Given that different ans¨atze are available, we will study their effect on somebulk properties in section 5 and we endeavour to get a handle on systematic un-certainties inherent in this approach. Finally we fix the free parameters of themodel and summarise our results in section 6. Finite temperature simulationswill be dealt with in [53] and [54].
2. Saturating the path integral
The IILM path integral is an approximation to the fundamental path integralby saturating the latter with a given ansatz for the multi-instanton background.The functional measure consists of small fluctuations around that classical con-figuration. To make analytical progress, the action is expanded to quadraticorder to define the ‘free’ part that is used in perturbation theory. In generalthe background induces non-Gaussian fluctuations that need to be treated ex-actly. The directions of these zero modes can (sometimes) be integrated up,and correspond to the tangent space of a generically non-trivial manifold. Thecoordinates on this so-called moduli-space can be interpreted as those degreesof freedom whose quantum mechanics approximates the low-energy dynamicsof the fundamental theory.In order to discuss the approximations that are eventually used, we willnow sketch the construction of the variational path integral, paying particularattention to the low lying modes. Details pertaining to the variational approach,gauge fixing and renormalisation can be found in the original papers [13, 14].We denote by φ the collection of bosonic fields. The classical action we write as S c and the classical interaction is defined as S int = S c − N S , (2)where N is the number of instanton constituents of the background and S theaction of an individual instanton. Assume that the background has N γ (quasi)3ero modes, which we denote collectively by γ . We can then write , using theeigenfunctions of the free part of the action δ S/δφ (cid:12)(cid:12) φ = φ c at zero and finite γ , φ c ( x,
0) + φ ( x ) = φ c ( x,
0) + ∞ X n =1 ζ n η n ( x, , = φ c ( x, γ ) + ∞ X n = N γ +1 ¯ ζ n η n ( x, γ ) + O ( γ ) . (3)This can be rearranged to (omitting the x -dependence for notational clarity) φ ( { γ, ¯ ζ } ) = φ c ( γ ) − φ c (0) + ∞ X n = N γ +1 ¯ ζ n η n ( γ ) , = ∞ X m =1 "Z d n x φ c ( γ ) − φ c (0) + ∞ X n =2 ¯ ζ n η n ( γ ) ! η m (0) η m (0) , (4)where we used the fact that η ( x,
0) forms a complete basis. Clearly φ ( γ = 0 , ¯ ζ ) ⊥ η i with i = 1 , . . . , N γ . The Jacobian for the variable change { ζ n } → { γ, ¯ ζ m } ,follows from the following partial derivatives ∂ζ n ∂γ i (cid:12)(cid:12)(cid:12)(cid:12) γ =0 = Z (cid:0) ∂ γ i φ c (0) η n (0) − φ (¯ ζ ) ∂ γ i η n (0) (cid:1) , (5) ∂ζ n ∂ ¯ ζ m (cid:12)(cid:12)(cid:12)(cid:12) γ =0 = δ mn . (6)Note the occurrence of the φ part. This will lead to new interactions which haveno classical counterpart but are purely quantum mechanical; to 1-loop order,we are allowed to discard them. From this matrix structure it follows that theJacobian is given by det( R ∂ γ n φ η m ).Now, we do not know the set { η } of exact low lying eigenfunctions. However,we can approximate it by constructing an orthonormal set of the known singleparticle zero modes that descend from the exact solutions used to build up thebackground field. With a slight abuse of notation, we substitute η → ¯ η = O B η ; O B is the matrix that generates an orthonormal basis from the original set { η } of single particle zero modes. The Jacobian is then given bydet (cid:18)Z ∂ γ n φ ¯ η m (cid:19) = det (cid:18)Z ∂ γ n φ η m (cid:19) det O B . (7)It corresponds to the quantum mechanical gluonic interactions. The high-frequency eigenvalues, encoded in the determinant of the fluctuation operator δ S/δφ are assumed to be N -fold degenerate, and so the fluctuation determi-nant factorises. This follows [27].
4n QCD we also need to introduce quarks and treat their interactions withthe background field. In the case where the Dirac operator admits quasi-zeromodes we can approximate the low frequency part in the same way as for thegluonic case; the high-frequency fluctuations will again be assumed to factorise.As for the Jacobian, we do not know the exact set of low-lying eigenfunctionsfor the superposition, but we approximate it by constructing an orthonormalset of the exact single particle zero modes ξ n , i.e. ¯ ξ = O F ξ . The Dirac operator,truncated to that subspace, is then given by( D/ + m ) low = ( D/ + m ) ij | ¯ ξ i i ⊗ h ¯ ξ j | = (cid:16) ¯ D/ ij + mδ ij (cid:17) | ¯ ξ i i ⊗ h ¯ ξ j | , (8)with ¯ D/ ij = h ¯ ξ i | D/ | ¯ ξ j i . The matrix of overlaps is related to the single particlezero mode overlaps by D/ low = ¯ D/ = O † F D/ O F , (9)with D/ ij = h ξ i | D/ | ξ j i . Note that this is not a similarity transformation because O F is not unitary.
3. Interactions in the IILM
We will now turn to instantons in QCD. In this paper, we will only discussBPST instantons [2]. In terms of the ’t Hooft potential Π( x, { y, ρ } ) = ρ r , (10)with r = ( x − y ) , the BPST instanton in singular gauge is given by A aµ = − O abi ζ bµν ∂ ν Π( x, { y, ρ } )1 + Π( x, { y, ρ } ) , (11)with ζ bµν = ¯ η bµν for instantons, ζ bµν = η bµν for anti-instantons and η the ’t Hooftsymbols. The collective coordinates are: y the centre, ρ the size and O thecolour orientation in the adjoint representation.The simplest background configuration is the sum ansatz, as used for in-stance in [13]. It was shown in [41], that the sum ansatz produces an unphysicalamount of repulsion; this is due to the fact that the field strength actuallydiverges at the individual centres, and is in sharp contrast to the individual sin-gular gauge instanton whose field strength is finite at the centre . In this case, Actually 1 + Π is the ’t Hooft potential. Note, however, that the field strength of the individual singular gauge instanton is notcontinuous at the centre and is only defined on the punctured Euclidean space. Incidentally,the winding about this singular point corresponds to the winding at infinity of the regularinstanton. R E Interactions for ratio ansatz as derived in this paper. R H Gluonic interactions are derived from the ratio ansatz whereasthe quark overlaps use a sum ansatz, [39]. S Interactions have been derived from the so-called streamlineansatz. These are only available at zero temperature, [39] [48].
Table 1: Several ans¨atze for the classical background field have been proposed. The followingtable summarises what they will be referred to throughout the rest of the paper. the author therefore proposed a different ansatz, inspired by ’t Hooft’s multi-instanton form, that stays finite at the centre of the instantons, and dubbed itthe ratio ansatz. It is given by A aµ = − P i O abi ζ bµν ∂ ν Π i ( x, { y i , ρ i } )1 + P i Π i ( x, { y i , ρ i } ) . (12)In what follows we will refer to R E as the interactions or the ensemble generatedby the ratio ansatz. We will compare the predictions from R E with those of thestreamline ansatz S [48] and another ‘hybrid’ ratio-sum ansatz R H [39]. Thisis summarised in Table 1.Phenomenological considerations have lead to the conclusion that the QCDvacuum consists of a dilute ensemble of instantons; a fact corroborated by latticestudies and self-consistency checks within the IILM. Diluteness and the localisednature of instantons render negligible contributions other than two-body inter-actions , given here for an instanton–anti-instanton pair, A aµ = − ¯ η aµν ∂ ν Π ( x, { x , ρ } ) + O ab η bµν ∂ ν Π ( x, { x , ρ } )1 + Π ( x, { x , ρ } ) + Π ( x, { x , ρ } ) , (13)with O = O t O . The formulas for an like-charged pairs follow trivially from theabove. The complete classical gluonic interaction is given by the sum over all thepossible pairings. It is clear from the structure of (13) that the colour degreesof freedom can be completely factorised out. After some lengthy algebra the A note on terminology: whenever we use the word interaction, we mean a quantity ‘nor-malised’ to the dilute gas, i.e. we subtract the dilute gas counterpart if the term naturallyoccurs in the exponential, as in the classical gluonic interactions, or we divide by the dilutegas counterpart if the interaction is a pre-exponential factor, as in the gluonic Jacobian or theDirac determinant. R / 0 0.5 1 1.5 2 2.5 3 3.5 g S -5051015 : instanton-instanton E R : instanton-instanton H R : instanton-anti-instanton E R : instanton-anti-instanton H R Figure 1: For instantons with equal sizes the interaction of R E agrees very well with R H foroppositely charged instantons. There is slight discrepancy for like-charged instantons, in thatthe repulsion is a bit steeper in the R E case. (We have set ¯ ρ = q ρ + ρ .) result for the squared field strength can be written in the form F aµν F aµν = I + (Tr O t O + (¯ ηOη ) µνµν ) J + (¯ ηOη ) ρµρν I µν + (¯ ηOη ) µρνσ I µρνσ + ( ηO t Oη ) µρνσ J µρνσ + (¯ ηOη ) αµαρ (¯ ηOη ) βνβσ K µρνσ . (14)The different contributions are given in appendix Appendix A.Factorising out the single instanton contributions and the coupling constant,the classical interaction between instantons is given by S g /S ≡ V ≡ ( S [ A ] /S − , (15) S [ A ] = 14 g Z F aµν F aµν , (16)where S [ A ] is the action of the background gauge fields and S = 8 π/g that of asingle instanton. For equal sizes the agreement with R H is very good, see Fig. 1.However, for unequal sizes there are noticeable differences, see Fig. 2. Thediscrepancy follows from the functional dependence on the size parameters beingof the form √ ρ ρ in R H . As can be seen from the asymptotic behaviours, seeappendix Appendix A.2, the sizes enter rather in the combination p ρ + ρ , atleast in the parameter regions of large and small separations; this is in agreementwith [13] . The authors of [13] have considered the sum ansatz; for large separations, however, everyansatz is equivalent to the sum ansatz. R / 0 0.5 1 1.5 2 2.5 3 g S -4-20246810121416 : instanton-instanton E R : instanton-instanton H R : instanton-anti-instanton E R : instanton-anti-instanton H R Figure 2: For unequal size parameters, e.g. ρ /ρ = 3 in this case, large differences start tobecome apparent. The reason is that the dependence on the sizes is more complicated thanthe functional form √ ρ ρ used in R H . Note that the attractive well is much deeper in the R H case which will eventually lead to a denser ensemble. (We have set ¯ ρ = q ρ + ρ .) We will a adopt a couple of simplifications in the practical implementationthat have been introduced in previous work. These are, on the one hand, theapproximation of the high-frequency quantum interaction by an inverse runningcoupling constant evaluated at the scale of the mean instanton size and, on theother hand, the neglect of the Jacobian that introduces the collective coordinatesand represents the low-frequency quantum interaction.In the single instanton case the high frequency quantum fluctuations leadto charge renormalisation and the coupling constant is replaced by the runningcoupling at the scale given by the instanton size [47], S ( ρ ) = 8 π/g ( ρ ). Thesame calculation has never been performed for a pair. In the original paper [13],the interaction part of high-frequency quantum fluctuations have been estimatedto be subdominant to the classical interactions . In that paper the authors arguethen that the quantum interaction can be estimated by modulating the (total)classical interaction with the inverse running coupling constant, a slowly varyingfunction of the background field, at the scale provided by the mean instantonsize ¯ ρ . We will adopt the parametrisation put forth in [41, 39] that estimatesthe scale of the running coupling constant on a pair-by-pair basis and uses thegeometrical mean of the sizes to set this scale. The full gluonic interaction isthen given by S g = S ( √ ρ ρ ) V . (17)The Jacobian interaction is positive by definition and can therefore be inter-preted as a repulsive (low-frequency) quantum interaction. A rough estimate of8he large distance behaviour suggests that the asymptotic power-law decay is O (1 /R ), with R the separation between the pair. This is a faster decay thanthe well-known dipole–dipole interaction that follows from the classical action.For strong overlaps the Jacobian matrix will become approximately degenerate,and its determinant small, essentially because the matrix elements of the pairwith the other instantons will be roughly equal. For complete degeneracy thesingularity should be logarithmic because one singular value will tend to zero asthe rank of the matrix decreases by one. The repulsion will thus be of the formln J sing12 ∼ ln (cid:18) R ρ + ρ (cid:19) , (18)with proportionality factor of order O (1) − O (10) because, as we argued, thedegeneracy is due to one overlapping pair and should not get contributionsfrom other instantons. In (4.1) we will discuss the small separation asymptoticbehaviour for the ratio ansatz; the analytical expressions are given in appendixAppendix A.2.2, and we note that the singular behaviour is also repulsive andlogarithmic, I sing IA ∼ ln (cid:18) ρ I + ρ A R (cid:19) , (19)with a proportionality factor that is again of order O (1) − O (10). In the in-termediate region it is harder to estimate the Jacobian interaction, but thelogarithm should make its contribution subdominant. Also, the classical inter-action is boosted by the quantum contribution through charge renormalisation.We conclude that the Jacobian interaction is probably negligible compared tothe classical interactions.Thus, the gluonic interactions we will use in this work will be given solelyby the classical interaction. The quark interaction arises from (8), as is clear from our discussion insection 2, and is purely quantum mechanical. As for the gluonic interaction,some further approximations have been used in the literature; we will adoptthese, albeit rephrased sightly differently.We will assume that the single instanton zero modes { ξ } form a functionalorthonormal basis, i.e. we neglect contributions arising from non-vanishing over-laps among the ξ i . With this in mind, the finite dimensional low-frequency Diracoperator is then given by( D/ + m ) ij = h ξ i | D/ + m | ξ j i = D/ ij + mδ ij . (20)To reiterate, we attribute the diagonal mass term to the requirement of or-thonormality rather than the degree of dilution of the instanton ensemble, e.g. Writing H = H ij | ψ i i ⊗ h ψ j | makes only sense if { ψ i } forms an orthonormal system, giventhe scalar product h·|·i . R / 0 1 2 3 4 5 6 | I A | T ρ / ρ : E R =1 ρ / ρ : H R =3 ρ / ρ : E R =3 ρ / ρ : H R Figure 3: The relatively large discrepancy is due to the fact that R E uses the full ratio ansatzin the Dirac operator whereas R H uses the sum ansatz. (We have set ¯ ρ = q ρ + ρ .) [37]. On the practical level this is irrelevant in as far as we recover the samedeterminantal interaction as used in previous works.The quark zero mode, in singular gauge and in the chiral representation, isgiven by [24] ξ I = 12 πρ I p I ∂/ Π I I (cid:18) U I ϕ (cid:19) , (21) ξ A = 12 πρ A p A ∂/ Π A A (cid:18) U A ϕ (cid:19) , (22)with ϕ αa = ǫ αa , normalised according to ǫ = 1. Finally, U i is the 3 × O ab = 1 / U τ a U † τ b ), with U = U † I U A .The Dirac operator, as defined above, is anti-hermitian. Eventually we needto diagonalise it, but, since readily available routines work with hermitian ma-trices, we display here the matrix elements of iD/ . Within the ratio ansatz, thematrix elements T IA = R ξ † I iγ µ D/ µ ξ A are as follows T IA = Z d x π ρ I ρ A
12 Tr(
U τ + β ) I β . (23)The concrete realisation of I β is given in appendix Appendix B.The rather large difference between R E and R H , see top of Fig. 3, is dueto the fact that the latter use a sum ansatz. The ratio ansatz was introduced10 R / 0 2 4 6 8 I A q S -25-20-15-10-50 =0.4 s =0.02 m d =0.01 m u : m E R =0.4 s =0.02 m d =0.01 m u : m H R =0.4 s =0.2 m d =0.1 m u : m E R =0.4 s =0.2 m d =0.1 m u : m H R Figure 4: On the level of the effective interaction, the difference between R E and R H isnot as pronounced, i.e. the relative difference has decreased substantially. We can clearlysee that light quark masses lead to a stronger attractive interaction between instantons andanti-instantons. Note that the relative difference between the ans¨atze R E and R H does notseem to depend strongly on the quark masses. The instantons have been set up with equalsizes. (We have set ¯ ρ = q ρ + ρ .) to remove the unphysical divergence in the field strength; no such problemafflicts the overlap matrix elements. On top of that the quark determinant isa pre-exponential factor and as an effective interaction the extra logarithmicfactor should make it rather insensitive to its exact form, see [45]. Withinour numerical framework, the full ratio ansatz does not produce any additionaloverhead and has the merit to be more consistent with the gluonic interactions.We have checked that upon neglect of the contributions special to the ratioansatz, i.e. simplifying the overlaps so as to recover the sum ansatz, our resultsagree very well with those of R H , apart form the aforementioned discrepancy inthe instanton size parametrisation. As for the gluonic interactions, the colourmatrices could again be completely factorised out.Note that the Dirac operator only connects instantons to anti-instantonsdue to the extra γ -matrix factor as compared to the mass operator, which van-ishes between instantons and anti-instantons. Therefore the quark fluctuationoperator has the following form m I − i (cid:18) TT † (cid:19) , (24)with T the N I × N A matrix of overlaps T IA , and N I ( N A ) is the number ofinstantons (anti-instantons); the 0-matrices are N I × N I and N A × N A dimen-sional, respectively; finally, I is the identity operator on the quasi-zero mode11pace of dimension ( N I + N A ) × ( N I + N A ). To diagonalise iD/ , it suffices toknow the singular-value-decomposition of T . The left and right singular vectors, ψ L and ψ R , are defined by T ψ Rn = λ n ψ Ln , (25) T † ψ Ln = λ n ψ Rn . (26)The singular eigenvalues λ n are always positive. The kernel of the Dirac operatoris spanned by the λ = 0 singular eigenvectors, ψ K , of either T or T † , dependingon whether N I < N A or N I > N A . We can then construct the eigenvaluedecomposition of the Dirac operator. The non-zero eigenvalue part has thefollowing eigensystem (cid:26) (cid:20) λ n , (cid:18) ψ Ln ψ Rn (cid:19)(cid:21) , (cid:20) − λ n , (cid:18) − ψ Ln ψ Rn (cid:19)(cid:21)(cid:12)(cid:12)(cid:12)(cid:12) n ∈ { , . . . , min( N I , N A ) } (cid:27) . (27)Finally, the kernel is spanned by the eigensystem (cid:26) (cid:20) , (cid:18) ψ Kn (cid:19)(cid:21)(cid:12)(cid:12)(cid:12)(cid:12) n ∈ { , . . . , N A − N I } (cid:27) , N I < N A , (cid:26) (cid:20) , (cid:18) ψ Kn (cid:19)(cid:21)(cid:12)(cid:12)(cid:12)(cid:12) n ∈ { , . . . , N I − N A } (cid:27) , N I > N A . (28)Note that the non-zero eigenvalues come in pairs. Together with the zeroeigenvalues, the determinant of the Dirac operator can be written asdet( iD/ ) = m | Q | min( N I ,N A ) Y n ( m + λ n ) , (29)with Q = N I − N A the topological charge. If we are only interested in the deter-minant, and not so much in the eigensystem, this can be put in the equivalentform m | Q | (cid:26) det( T T † + m ) , Q < T † T + m ) , Q > . (30)Upon placing this term into the exponential, the normalised determinant ofquark zero mode overlaps leads to an effective interaction. The normalisationconsists of dividing (30) by m N I + N A . The quark interaction is thus given by S qN f = − N f X n =1 (cid:26) ln det( T T † + m n ) − N I ln m n , Q < T † T + m n ) − N A ln m n , Q > , (31)with N f the number of active quark flavours. Note that the quark interactionis always attractive. This follows from the fact that we can write the overlapmatrix for each flavour as I + T m n , and this form makes it explicit that thedeterminant is bounded from below by unity because the smallest eigenvalue iseasily seen to satisfy λ min ≥
1. 12his exhausts the interactions in the IILM because the fluctuation operatorof the ghost part is positive definite and its lack of zero modes prevents theconstruction of the low frequency part of the spectrum within the moduli-spaceapproximation. We are thus left with the high-frequency part which, as in theother cases, is assumed to factorise and cannot lead to interactions.
4. Numerical Implementation
The decoupling of the colour degrees of freedom is a computational benefit:by using global SO (4) transformations, without loss of generality, we place thefirst instanton at the origin and the second along the z -direction. The initialorientational dependence is then factored out of the integrand and combineswith the colour matrices as in [13]. These integrations are too time consumingto perform during actual simulations; instead, they are computed beforehandto fill interpolation tables that are, in turn, used during the simulations. Theinterpolation grid is three-dimensional, and depends on ρ , ρ and R = | x − x | .For numerical stability we choose to use simple linear interpolation.A uniform grid can, of course, only extend over a finite region and we mustdecide which portion of the parameter space to cover. We took the single in-stanton moduli-space measure as a guide for the size grid because, suitablynormalised, it can be interpreted as a probability density. We choose the lowerlimit, ρ min ≈ Λ, to be a fairly small quantile . Here, Λ is the scale at whichQCD starts to become strongly coupled. The upper limit is set to ρ max = Λ.Larger instantons cannot be treated consistently in the IILM because it usesperturbation theory, which breaks down below Λ.We believe that these choices cover the relevant parameter space, and wesample the sizes from the interval [ ρ min , ρ max ]. As a consistency check we moni-tored the actual size distribution and did not find any evidence for a significantweight at the edges of the sample interval. We therefore conclude that thisprocedure is well-defined.The classical gluonic interaction in the ratio ansatz suffers from gauge sin-gularities that prevent us from extending the grid down to vanishingly smallinstanton separations, R →
0. The opposite limit, R → ∞ , cannot be coveredeither unless we use a non-uniform measure on R + . In principle this would seemlike the most elegant approach, however, it is not feasible practically because thenumerical integration becomes inaccurate at larger separations; the only rem-edy would be to set very small error tolerances for the numerical integrations,but that is computationally prohibitive. Therefore, we decided to use matchingformulas for both the large and small separation regimes.The rationale is not to derive accurate formulas in absolute terms but to getthe absolute value from the interpolation results at a matching point R m . The It corresponds to less than the millionth quantile. (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:1) (cid:0)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) I I R Figure 5: The instantons I and I are so far apart that, within the shaded region that givethe dominant contribution to the field strength of each, the partner’s field strength is roughlyconstant and fixed at x µ − R µ ≈ − R µ . We can then safely extend the integration regionto be all of R , with a negligible error due to the rather strong localisation of the individualinstantons. matching formulas are thus to be understood as accurate in a relative sense, i.e.the asymptotic interactions f asy should behave, asymptotically, like the exactnumerical interactions f ex . This ensures that we reproduce the correct fall-offor singularity behaviour. Thus, we compute the interactions according to f ( R ) = f asy ( R ) f ex ( R m ) f asy ( R m ) , (32)whenever they fall out of the grid. Since the localisation of the instantons is setby the sizes, it is natural for the matching point to be proportional to the for-mer. Eventually, the exact proportionality factor follows from an ‘optimisation’procedure, given that we aim for the interpolated interactions to be correct atthe one percent level.The full gluonic interaction consists of different pieces that are added to-gether, (14). We could use (32) for these subinteractions term by term butit turns out that such a matching is numerically rather unstable. Thus, eventhough we are only interested in asymptotic relations, we need a systematicprocedure that insures that the different asymptotic subinteractions are addedup with the correct magnitude relative to each other.For the large separation case we want the instantons to be so far apart fromeach other that within the region in which the field strength for I is strongthe field strength of I hardly changes: we can approximate x µ − R µ ≈ − R µ .Since the field strength is negligible at and beyond R µ , we can safely extendthe integration region to cover all of R . The field strength of I behaves as aconstant, and we can use the rather simple rational expression for the interactionin terms of the ’t Hooft potential to find exact results. We add to this theanalogous contribution from I ↔ I . The configuration is illustrated in Fig. 5.We shall call this the zeroth order approximation, and it is clear that, to thisorder, terms odd in derivatives of I i will vanish due to O (4) symmetry. However, We use a translation to place I at the origin. (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) R Large separation asymptotic integration up to R (II) I I Proper small separation asymptotic integration down to R (I)neglected Figure 6: The instantons I and I are strongly overlapping. We approximate the integralby, first, integrating over I keeping I fixed at R µ , as in the large separation case but withupper limit R/
2; to this we add the analogous contribution from I . Secondly, the, possibly,singular behaviour is picked up by integrating from infinity down to R , and approximatingthe arguments to be x µ − R µ / ≈ x µ and x µ + R µ / ≈ x µ respectively. it turns out numerically that, upon combining all the different terms from (14),the zeroth order terms are sufficient. In particular, no non-integrable termsare present; as we will demonstrate in [54] , the finite temperature interactionsdo not include terms such as (1) that prohibit the thermodynamic limit. Ourformulas are given in appendix Appendix A.2.1.In principle, we can compute the neglected terms by going to first order, i.e. g ( x − R ) ≈ g ( − R ) + x µ ∂ µ g ( − R ), or beyond. Such higher order contributionswill typically no longer converge on R . It seems natural to cut them off at R , and this will generally lead to logarithms, ln(1 + R / ( ρ + ρ )), togetherwith rational functions. However, in contrast to the fitting formulas of [39],the Taylor expansion of our asymptotic formulas produce only power-law likedecays for large separations; in addition they fall off more strongly than thezeroth order terms and thus will not produce non-integrable interactions either.We have thus achieved our goal of deriving interactions that allow us to studythe thermodynamic limit of the IILM.We now turn to the case of asymptotically small separations. A typicalsituation is depicted in Fig. 6. The rationale is to split the integration into 2regions. At finite temperature, the present framework remains virtually unchanged. The onlymodifications are that the ’t Hooft potential, Π, changes and that the integration regionbecomes S × R . -3 -2 -1
10 1 10 I n t e r ac t i on -12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1
10 Gluonic interaction (approximation)Gluonic interaction (exact)Quark interaction (approximation)Quark interaction (exact)
Figure 7: The gluonic interaction is well approximated by the combination of interpolationand asymptotic matching. The quark overlap is very poorly approximated by the zeroth ordersmall separation asymptotic formula; it tends to zero with too high a power as compared to theexact result. The correct behaviour can in principle be obtained from higher orders, and we’veestimated that the second order contribution will suffice. In practice, the quark interaction inthis region is completely irrelevant as compared to the gluonic interaction.
I The far-field region beyond both centres, placed symmetrically around theorigin; we approximate the arguments by x µ ± R µ / ≈ ± x µ .II The region around each instanton up to R/
2, with R the pair separation.We integrate around ± R µ / x µ ∓ R µ / ≈ ∓ R µ /
2. This is similar to the large separationcase, but here we only integrate up to R/ -2 -1
10 1 10 I n t e r ac t i on -100102030405060 ApproximationExactR -3 -2 -1
10 1 10 R e l a t i ve d i ff e r e n ce -3 -2 -1 zero crossingMatching Large SeparationMatchingSmall Separation Interpolation Figure 8: The total pair interaction is accurate on the one percent level. Note that the errorin region I from the quark interaction is completely negligible. The spike in the bottom plotis due to the sensitivity to zero crossings.
Previous studies, lattice results and phenomenology indicate that the instan-ton ensemble is fairly dilute. Therefore, we organise the partition function intoa dilute gas measure times the exponential of interactions, Z = ∞ X N I ,N A N I ! 1 N A ! N I Y i d ( ρ i ) N A Y j d ( ρ j ) exp ( − ( S g + S q )) , (33) ≡ ∞ X N I ,N A N I ! 1 N A ! Z N I ,N A , (34) S g = X ij S gij , (35)where S gij is given in (17) and S q is given in (31). We follow [39] and use one-loop accuracy for the charge renormalisation factor that modulates the classicalgluonic interaction, i.e. S ( √ ρ ρ ) → β ( √ ρ ρ ), with β given below in (41).Although not really consistent, the single instanton density is given at two-loopin order to replace the pre-exponential bare by running coupling constants [13];the former were induced by the transformation to collective coordinates. Thetwo-loop single instanton measure is then given by d ( ρ ) = d g ( ρ ) d q ( ρ ) N f , (36) d g ( ρ ) = C N c ρ − β ( ρ ) N c exp (cid:20) − β ( ρ ) + (cid:18) N c − b ′ b (cid:19) b ′ b ln β ( ρ ) β ( ρ ) (cid:21) , (37) d q ( ρ ) = mρ exp (cid:18) −
13 ln mρ (38)+ ln mρ + 2 α − (6 α + 2 β )( mρ ) + 2 A ( mρ ) − A ( mρ ) − mρ ) + B ( mρ ) + B ( mρ ) + B ( mρ ) (cid:19) . (39)18or the quark term, d q , we use the generalisation of ’t Hooft’s [47] result validfor arbitrary mass [11]. The different terms in d g are given by C N c = 0 . e − . N c ( N c − N c − , (40) β ( ρ ) = − b ln( ρ Λ) , b = 113 N c − N f , (41) β ( ρ ) = β ( ρ ) + b ′ b ln (cid:18) b β ( ρ ) (cid:19) , b ′ = 343 N c − N c N f + N f N c . (42)Note that the above has been derived in Pauli-Villars regularisation.Being an interacting many-body system, the partition function cannot beevaluated analytically, and we choose Monte Carlo methods to cope with itnumerically. More precisely, we will use the Metropolis algorithm to samplethe important integration regions of the partition function. This is, of course,all well known, but it seems appropriate to introduce the, possibly less known,Monte Carlo moves corresponding to insertion and deletion of instantons neededfor grand canonical simulations.Following the usual strategy of imposing detailed balance, the simplest in-sertion/deletion algorithm consists of randomly placing an instanton in the boxand randomly selecting an instanton to be removed. Imposing detailed balanceand considering the case of an instanton, we arrive at1 V p eq N I ,N A A N I ,N I +1 = 1 N I + 1 p eq N I +1 ,N A A N I +1 ,N I . (43)As usual, p eq N I ,N A = Z N I ,N A /Z is the probability to be in the state { N I , N A } .The acceptance probability A ij is implicitly defined through (43), and theMetropolis algorithm defines it to have the following form, [19], A N I ,N I +1 = min(1 , A ) , (44) A N I +1 ,N I = min(1 , A − ) . (45)Plugging this into (43) we finally arrive at A = VN I + 1 p eq N I +1 ,N A p eq N I ,N A . (46)The difference to ordinary Monte Carlo moves, as used in the canonicalensemble , is that the proposal probabilities do not cancel and the transition Note that we neglect the factorial terms in the definition of the equilibrium probabilitydensity p eq because they are an artifact as far as the measure is concerned. They have beenintroduced to render the integration volume simple, i.e. the product of the single instantonmoduli-spaces M N I . During the integration process all the permutations of a given set ofcoordinates are generated, but, since the instantons are indistinguishable, they really corre-spond to only one configuration. To correct for this overcounting, we then have to divide bya factor of N I !. The important point is that for the transition probabilities these factorialfactors are irrelevant. That is, updates for the positions, sizes and colour orientations C δ GC h N i ξ N h Q i ξ Q h S int i ξ int . . . . − . . . . . − . . . . . − . . . . . − . Table 2: The sample size is roughly equivalent for each set, with 200 independent configura-tions generated according to the autocorrelation time ξ N . Considering some bulk properties,we see that the sampling does not really depend on the a-priori-probabilities δ i . Even thoughthe autocorrelation times are only rough estimates, we will take the data at face value andchoose δ C = 0 . δ GC = 0 . matrix is not symmetric. In this specific case, the proposal probability for aninsertion is P ins P = 1 /V , corresponding to the probability to place the instantonrandomly within the box, whereas the proposal probability for a deletion is P del P = 1 / ( N I + 1), corresponding to the probability to select an instantonamong the N I + 1 available.When we perform the standard updates, it is easy to monitor the acceptancerates and tune the the proposal probabilities to achieve good rates, i.e. 50% say.For the move described by (46) we do not have a parameter to tune though. At T = 0, this is not really a big issue because it turns out that the acceptancerate is ≈ .
4, even for the rather small quark masses that we will use. This isstill acceptable and does not really justify the overhead of more sophisticatedupdate algorithms.Also, note that such an insertion/deletion step is a fairly large change ascompared to the normal coordinate updates, and so these grand canonical movesactually help to sweep through phase space more quickly.Finally, we need to decide how many grand canonical moves we perform percoordinate update, that is we need to fix the a-priori-probabilities δ C and δ GC .We found that, for T = 0, the ensemble is not sensitive at all to this parameter,see Table 2. Since we will ultimately be interested in computing the topologicalsusceptibility, we will aim to achieve low autocorrelation times for the instantonnumber, N = N I + N A , and topological charge, Q = N I − N A , i.e. we willperform rather more insertion/deletion moves than less. In practice we performcanonical moves only 60% of the time. As mentioned in the introduction, we want to study the IILM for ’physical’quark masses. In that case, we must make sure that the simulation box is largeto be insensitive to finite size effects. In the lattice community it is commonpractice to use a box length that corresponds to 4 − V → ∞ . 20s compared to fitting formulas, our combination of interpolation and asymp-totic matching results in a rather substantial computational overhead. This isparticularly so in the quenched case. For unquenched simulations the situationis less drastic as the computationally most demanding part is the evaluation ofthe determinant and/or the determination of the eigensystem of the Dirac op-erator. Increasing the simulation box, i.e. increasing the number of instantons,this becomes the bottleneck to large volume simulations.The Monte Carlo changes are, however, of a rather simple form, changingonly one column of the overlap matrix T at a time. We can therefore use de-composition update techniques to reduce the complexity from O ( N ) to O ( N ).For the update step we only need to evaluate the determinant (30). Giventhe fact that m + T T † , respectively m + T † T , is a positive definite hermitianmatrix, the fastest evaluation will be achieve by using the Cholesky decompo-sition. An added bonus is that the Cholesky decomposition and its algorithmare known to be very stable.Focusing on M = m + T T † = LDL † , an update T ′ = T + ∆ T can bewritten as two rank 1 updates for M , of the form M ′ = M + ΦΦ † − ΨΨ † , (47)with Φ, Ψ vectors. Details are given in appendix Appendix C, where we alsodiscuss more efficient ways to deal with adding and removing instantons, and thecorresponding updates. The Cholesky decomposition can be updated efficientlywhen it only changes by rank 1 matrices, that is transformations of the form M ′ = L ′ D ′ L ′ † = M + αzz † = L ( D + αww † ) L † , (48)where Lw = z . The algorithms then compute the decomposition of D + αww † =˜ L ˜ D ˜ L † , which can be achieved in O ( N ) because D is diagonal. Furthermore,the matrix ˜ L has a special form which allows an efficient matrix multiplication, L ′ = L ˜ L , in O ( N ). Details can be found in [22]. The algorithm we use inpractice is known to be unstable for downgrading, α <
0, unless the resultingmatrix, M ′ , is known to be positive definite. Since upgrading, α >
0, is alwaysstable, it is important to perform the two consecutive updates in the order givenby (47).In general we will be performing grand canonical simulations, and need tokeep track of two decompositions, one for m + T T † and one for m + T † T .Furthermore, we deal with 3 active quarks so that each Monte Carlo updateentails 2 · ·
5. Different Ensembles
To be predictive, the IILM should not depend too sensitively on the chosenansatz. Given the fact that, for instance, the streamline and the ratio ansatz21ave quite different functional forms for overlapping pairs, the insensitivity ofthe model to specific background ans¨atze can only be determined a posteriori.On a heuristic, level we expect insensitivity to emerge if the ensemble stabilisesin a rather dilute form so that the precise functional form of the repulsionis irrelevant. The large separation limit is a priori unproblematic because allans¨atze are constructed such that, asymptotically, they approach the simplesum ansatz, i.e. A = A + A , with A the gauge field.First, we will frame our discussion on the pair interactions. The total effec-tive interaction of a pair of oppositely charged partners is given by S = S ( √ ρ ρ ) V − N f X n =1 ln (cid:18) | T | + m n m n (cid:19) . (49)Identically charged pairs only feel the gluonic interaction as T II = T AA =0. As expected, the ratio and streamline ansatz are markedly different onlyfor strongly overlapping pairs, see Fig. 9, where strongly overlapping pairs arecharacterised by R ≤ p ρ + ρ .In the quenched case, we notice that the R E ansatz has a higher absoluteminimum as compared to the R H ansatz, occurring roughly for the same sepa-ration. So we expect the ensemble to become slightly more dilute because it willnot be as favourable, energetically, for instantons to come close. For unequalsizes, however, the repulsion is weaker in the R E case which would favour adenser ensemble, as less volume is excluded. The streamline ansatz will leadto a substantially more dilute system because the core repulsion is broader,excluding more volume for the instantons to move through.In the unquenched case, the difference in the absolute interaction strength ismuch more pronounced between R E and R H . We therefore expect that the R E ansatz should be quite a bit more dilute as compared to the R H ansatz. Con-sidering that the streamline ansatz has a deeper minimum than the R E ansatz,the former will favour instantons to come closer. However, it has more excludedvolume. Both trends work in opposite directions, and there is a possibility thatthey lead to roughly identical ensembles, at least on the level of the instantondensity.We will address these issues in more detail by performing canonical simula-tions and minimising the free energy. This follows closely [39], and also serves tovalidate our code against their results. Note that the simulations are performedin the topologically trivial sector, for which N I = N A , i.e. the topological charge Q = 0. In Fig. 10 we plot the free energy F = − ln Z/V against the instantondensity n = ( N I + N A ) /V . As expected from our considerations of the pairinteractions, in the quenched case the R E and R H ans¨atze are only slightlydifferent, with R E leading to a slightly denser ensemble. Also, the interactionsstored in that ensemble are a bit lower, again as could be anticipated from the The difference in the free energies is directly related to the difference of the interactionper instanton R / 0 1 2 3 4 S -15-10-5051015 : quenched E R : quenched H RS: quenched : unquenched E R : unquenched H RS: unquenched ρ R / 0 1 2 3 4 5 S E R : quenched H RS: quenched : unquenched E R : unquenched H RS: unquenched
Figure 9:
Top:
Most attractive colour orientation.
Bottom:
Random colour orientation. Forboth graphs the instantons have different size parameters. We expect the R H ensemble tobe denser than the R E ensemble because the attraction well is deeper, whereas the excludedvolume due to repulsion is not that different. Along the same lines, the S ansatz should leadto a rather more dilute system in the quenched case. For unquenched simulations, the deeperattraction well, and the steeper and broader repulsion, of the S interactions might lead to anensemble roughly equivalent to the R E one. u m d m s M . . . M .
05 0 .
05 0 . M .
012 0 .
022 0 . Table 3: We use three different sets of quark masses and investigate how the instanton liquiddepends on them. pair interactions. The S ansatz leads to a much more dilute instanton ensemble,and our data reproduces well that of [39].Ultimately, we will be interested in smaller quark masses. It is clear from(49) and Fig. 4 that smaller masses increase the quark interaction strength ascompared to the gluonic counterpart, which stays constant and is responsible forthe core repulsion. From a purely energetic point of view, smaller quark massesshould then lead to denser ensembles. However, we clearly see in Fig. 10 thatthe ensembles become more dilute. The reason is that the small quark massesenter the instanton size distribution; in turn, the density, in the dilute gas limit,is entirely set by the size distribution, i.e. n = 2 R dρd ( ρ ). Bringing it into theaction, we can interpret the size distribution as the energy cost needed to insertan instanton into the box. This is a well-known fact, namely that small quarkmasses suppress instanton contributions to the QCD vacuum because the differ-ent topological vacua become equivalent in the limit of vanishing quark masses;phrased differently, the energy barrier has disappeared, and only field configu-rations with topological charge Q = 0 survive. In the dilute gas approximationthis leads to the disappearance of instantons altogether. As Fig. 10 shows, thisis not true for an interacting instanton ensemble, where the instanton densityconverges to a finite limit as the quark mass is lowered . The results fromFig. 10 also show that the R E ansatz generates an ensemble that differs moreand more from the R H ansatz, as was anticipated from our considerations ofthe pair interactions. We can also clearly see that the R E ensemble does notconverge to the S ensemble.So far we have framed the discussion essentially in terms of the instantondensity. To investigate the similarities and differences in more detail, we willlook at a few bulk properties and their dependence on the different ans¨atze inthe thermodynamic limit and the grand canonical ensemble. We clearly seehow the density decreases with the quark masses, but approaches a finite limit,Fig. 11. As we have discussed before, the quark masses will suppress fluctuationsto inequivalent topological sectors and in the limit of vanishing masses only thetrivial Q = 0 sector will survive. The fluctuations between topological sectorsare encoded in the topological susceptibility χ = h Q i /V , which vanishes withthe quark masses, see Fig. 12. Both the instanton number and the topological Remember that the simulations take place in the topologically trivial sector, i.e. N I = N A or Q = 0. Λ n ( -2 -1
10 1 ) Λ F ( -2-1.5-1-0.500.5 E R H RSQuenchedM1M2M3
Figure 10: The different quark masses are summarised in Table 3. The simulations wereperformed in the topologically trivial sector, i.e. N I = N A . For the quenched and the M case we fixed N = 64 as in [39]. The other two unquenched simulations have N = 200. Weclearly see that small quark masses suppress instanton contributions to the QCD vacuum,but also that there exists a finite limit for the instanton density as the quark masses vanish;this is in contrast to the dilute gas approximation which suppresses instanton contributionscompletely for zero quark masses. For unquenched simulations the free energy for the R E ensemble roughly agrees with that of the S ensemble, although the equilibrium densities arerather different. Still, the approximate equality between the free energies might be interpretedas evidence of an approximate equivalence between both ensembles for bulk properties, e.g.equivalent pressure, since it is directly related to the free energy. However, the R E liquid doesnot seem to converge towards the S ensemble as we lower the quark masses. -4 Λ V (0 500 1000 1500 2000 2500 3000 3500 4000 < N > E R H RSQuenchedM1M2M3
Figure 11: As anticipated from the canonical study, the instanton number for both ratioans¨atze is very similar in the quenched case. There seems to exist a finite limit for theinstanton density as the quark masses vanish, and instantons will be present in the QCDvacuum even in the chiral limit; this is in sharp contrast to dilute gas approximations. charge fluctuations exhibit a nice scaling with the volume.We will now turn to an intensive quantities, the mean instanton size ¯ ρ Fig. 13.For the quenched and the first unquenched, M , simulations, which have beencalibrated to achieve N ≈
64 for the smallest volume as in the canonical simu-lations, we see that the simulation boxes are not large enough, even though thedensity and the charge fluctuations seem to suggest otherwise, i.e. display goodscaling with V . For the other two unquenched simulations, tuned to N ≈ R E ansatz as compared to the R H ansatz dominates over the deeper attractivewell of the latter. Therefore, the total interaction in the R E ensemble is slightlylower, leading to a denser system. The stronger repulsion for the S interactionsleads to more excluded volume; this, in turn, leads to lower interactions anda more dilute ensemble. These conclusions are in agreement with the directmeasurement of the instanton density Fig. 11.We are mostly interested in unquenched results, and the following commentsrelate this sector. From the data of Table 4, we can infer that results for χ / -4 Λ V (0 500 1000 1500 2000 2500 3000 3500 4000 > < Q E R H RSQuenchedM1M2M3
100 200 300 400 50050100150200250
Figure 12: The topological susceptibility, the slope of the graphs, is very sensitive to thequark masses. It is screened by small quark masses and will vanish in the chiral limit. This isexpected as QCD with massless quarks does not have topologically inequivalent vacua; in thiscase the so-called θ parameter is not physical and can be rotated away by a chiral rotation ofthe quark fields. See also [46] for another work on the topological susceptibility in the IILM. ) -4 Λ V (0 500 1000 1500 2000 2500 3000 3500 4000 ) - Λ > ( ρ < E R H RSQuenchedM1M2M3
Figure 13: For the quenched and the M simulations, which have been fixed to N ≈
64 for thesmallest volume, the simulation boxes are still a too small, as can be seen by the systematicdrift. For the other two unquenched simulations ( N ≈
200 for the smallest volume) we aremuch closer to the thermodynamic limit, although there are still systematic deviations forthe R H ansatz. In any case, the different ans¨atze give rather similar results. Also the meaninstanton size approaches a unique limit for small quark masses. -4 Λ V (0 500 1000 1500 2000 2500 3000 3500 4000 > q + S g < S -6-5-4-3-2-101 E R H RSQuenchedM1M2M3
Figure 14: As anticipated from the pair interaction considerations, the interactions are verysimilar in the quenched sector for R E and R H . For full simulations, the differences betweenthe ans¨atze stay constant as the quark masses vary, with R H leading to the strongest attractiveinteractions, as was expected. Note that the pair interactions are less sensitive to finite sizeeffects compared to the mean instanton size. n χ h ρ i h S g + S q i /V Quenched 0 . R E ]0 . R H ]0 . S ] 0 . . . . . . . . . M . . . . . . . . . − . − . − . M . . . . . . . . . − . − . − . M . . . . . . . . . − . − . − . Table 4: Thermodynamic extrapolations for the instanton density, the topological suscepti-bility, the mean instanton size and the mean interaction. The data has been obtained fromFigs. 11, 12, 13 and 14 respectively. R E and R H ansatz in the unquenched sector. Themean instanton size is indeed a rather robust quantity and only affected on the3% level by these systematics. Also, note that the mass dependence on h ρ i israther small, with differences not larger than 5%. The instanton interactionsand the n / agree within 20%, and the latter converges to a fixed limit as thequark masses vanish.
6. Fixing parameters
In the quenched case, the IILM has only one freely adjustable constant, thelambda parameter Λ. We need one observable, from the lattice say, to fix it.Different approaches can be chosen. In the early works, Λ was determined byfixing the instanton density to 1 fm − at T = 0. To compare this with the latticeis not straightforward as the classical instanton content is convoluted with thequantum mechanical fluctuations. With the discovery of the KvBLL calorons,there is a renewed interest in studying the topological structures on the lattice,see for instance [3]. Since the topological susceptibility is well measured on thelattice and is easily accessible within the IILM, it is a natural candidate. Thelambda parameter is then given byΛ = r χ lat χ IILM . (50)We will use χ / = 193 MeV, [12]. The topological susceptibility in the IILM isextracted from Fig. 15 by using the definition χ top = lim V →∞ h Q i V . (51)This yields Λ = 234(1) MeV. The error is purely statistical. The instantondensity turns out to be n = 0 .
543 Λ = 1 . − , fairly close to the usuallyquoted phenomenological value of n = 1 fm − . We find that even for theselarger volumes the mean instanton size is still evolving towards lower values,as in Fig. 13. The largest volume then leads to the upper bound ¯ ρ < .
57 fm.Using a simple fit to ¯ ρ = ¯ ρ ∞ + αV − . to extrapolate to the asymptotic value,we find ¯ ρ ∞ ≈ .
53 fm; this is rather large compared to the phenomenologicalvalue of ¯ ρ ≈ .
33 fm.To estimate the systematic error due to the dependence on the ansatz, wewill use the data from Table 4. The fact that we take a fourth root reduces therather large differences in χ IILM to about 15% for Λ, i.e. Λ = 234(35) MeV. Ourvalue has been obtain through simulations in PV regularisation. To compareit with lattice data, we will convert it to the MS scheme, [26], Λ MS / Λ PV =exp( − / MS = 224(33) and compares well with the latticeresult Λ MS = 259(20) [23]. 29 -4 Λ V (100 200 300 400 500 600 700 > < Q
5; chiral perturbation theory can then beused to extrapolate to physical masses. Ultimately, the lattice wants to testthe predictions of chiral perturbation theory as well, and in recent years, thecomputing power and, most importantly, the algorithms have improved to suchan extent that physical quark mass simulations are becoming feasible; however,these are still immensely costly simulations, and 2 + 1 flavour simulations wererare until recently.We follow a rather more modest rationale by simply demanding that thequark mass be at least so small as to be comparable to the lowest eigenvalueof the Dirac operator, h λ min i , see Fig. 17. This sets the smallest box we usein our simulations. We then use ever larger volumes and extrapolate to thethermodynamic limit.In [9] the lambda parameter is fixed by computing the meson and nucleonmasses, through current correlators of the interpolating fields and their asymp- Actually, the authors fixed the mean instanton size. But it is trivial to relate the latterto the lambda parameter. h ¯ qq i . The chiral condensate has been studied within chiral perturbation the-ory and, more recently, it has been precisely determined on the lattice [7, 6, 5].We will take it to be h ¯ qq i MS0 ( µ = 2 GeV) = 250 MeV.To extract the chiral condensate from the IILM, we will use the procedureadopted in [7, 6, 5]: we compute the topological susceptibility for different setsof quark masses and extrapolate to the chiral limit, m i →
0. The condensatecan then be determined by chiral perturbation theory [33], χ = m eff h ¯ qq i + O ( m ) , (52) m eff = N f X n m n − . The chiral condensate has an anomalous dimension and, therefore, dependson the scale. Furthermore, the IILM is set up with a PV regulator, whereasthe quoted result is computed in dimensional regularisation. It is well knownthat within an unphysical renormalisation scheme such as MS the results de-pend on the regulator (for unphysical quantities like masses, coupling constantsand amplitudes). We therefore need to compute the finite counterterms thatrelate the PV to the MS regularised results. Deferring the details to appendixAppendix D, we find that h ¯ qq i PV0 ( µ = 2 GeV) ≈
244 MeV. This is a one-loopresult. The two-loop correction can be estimated very roughly to be on the 10%level as is typical for computations around the scale of µ = 2 GeV .We will define the scale of the IILM by µ Λ = Λ / ¯ ρ , as suggested in [9], anddetermine Λ from the self-consistency equation h ¯ qq i PV0 ( µ Λ ) = Λ h ¯ qq i IILM0 , (53)where we run the chiral condensate h ¯ qq i PV0 ( µ ) at one loop. To that order, thereis no difference between schemes and we can use the MS results, e.g. [49]. To reiterate, we implicitly rely on the fact that the IILM is describing well the chiralproperties of QCD, as has been checked in numerous studies, the most convincing being [9]. ’t Hooft’s computation of the one-loop instanton measure, using Pauli-Villars regular-isation, is also unphysical because, instead of poles, logarithms of the regulator mass aresubtracted. Strictly speaking, we should use the two-loop result because the simulations in the IILMhave been obtained using the two-loop improved instanton measure. However, Pauli-Villarsregularisation is not straightforward for non-Abelian gauge theories beyond the one-loop level,and we do not have the expertise to embark on this endeavour. In any case, the differenceshould still be on the 10% level.
31o get an estimate of the quark mass ratio dependence, we have used twodifferent sets of quark masses, one inspired by the chiral perturbation theoryand the other by the quark masses extracted from the lattice [34]. The two setshave the following ratios m i m j = (cid:26) .
83 : 36 . M )1 : 2 .
32 : 45 . M ) . (54)For each set we perform 5 simulations with ever smaller absolute masses, seeFig. 16. This data is fitted to (52) to extract the chiral condensate. The resultsfor the two sets agree on the 1 σ level, and we can argue that the chiral condensatedepends only weakly on the quark mass ratios, given that the latter vary byroughly 25%, see (54). This is as it should be since the exact chiral condensatedoes not depend on the quark masses at all. From an operational point of view,the robustness against quark mass ratios makes the chiral condensate a goodquantity to set Λ.Solving (53) we find that the lambda parameter is given byΛ i = (cid:26) , (55)where the errors follow from the fit, h ¯ qq i PV0 and the systematic on χ . This leadsto an overall error of 44 MeV, or roughly 11%, and is strongly dominated bythe one-loop result for h ¯ qq i PV0 . Since we run at one loop in the self-consistencyequation (53), such a large error is certainly realistic, if not underestimated.We found that running at two-loop gives results consistent with (55). Usingthe prescription of [9], the scale and the mean instanton size for the IILM is µ Λ i = (cid:26) , ρ Λ i = (cid:26) . . . (56)This is in very good agreement with the precision study [9]. Given that bothworks use chiral properties for the calibrations, the nice overlap is probably nottotally unexpected.Note that (55) is a prediction for the lambda parameter with 3 active quarkflavours. To compare our result with experimental data we run down the cou-pling constant α MS s = 0 . M Z to µ Λ and convert it to a lambdaparameter. This is a rather big difference in scales and it is appropriate to usetwo-loop running, although not entirely consistent when we compare it to the i.e. taking the limit from different directions in quark mass space. We use β -functions and anomalous dimensions from the MS scheme, [49], since we donot know them for PV regularisation. However, both regularisations are thought to giveroughly similar results, for instance h ¯ qq i PV0 and h ¯ qq i MS0 agree on the 3% level at one-loop and µ = 2 GeV, and for the purpose of estimating errors in the one-loop running this procedureshould be fine. Remember that in the unquenched case the instanton size is fairly independent of thequark masses. ff m0.002 0.004 0.006 χ eff2 + 5.87e+01 m eff =1.96e-01 m χ eff m0.002 0.004 0.006 χ eff2 + 5.62e+01 m eff =2.30e-01 m χ Figure 16: Computing the topological susceptibility allows for the extraction of the chiralcondensate by using chiral perturbation theory, (52). The rationale is the same as used inrecent lattice studies to extract the chiral condensate [7, 6, 5]. In order to get a rough estimateon the systematic error introduced by the chiral limit, two sets of masses have been used. Theupper plot corresponds to the set M and the lower plot to M , as given in (54). The chiralcondensate for both mass ratios agrees on the 1 σ level, and we conclude that it depends onlyweakly on the quark mass ratios. Mathematica package
RunDec , [4]. The conversion between the MS and PV lambda parame-ters is given by [26], [1]Λ PV = Λ MS exp (cid:18) − N f / (cid:19) . (57)This leads to Λ (3)PV = 325(40) and the IILM result agrees on the 1 σ level. Trustingthe perturbative running down to the rather low scale µ Λ is a leap of faith.However, earlier studies have seen good agreement between IILM and latticepredictions for physical quantities, such as meson masses, and so the agreementbetween the lambda parameters might not just be a fluke.To determine the physical quark masses, we will use (52) rewritten in termsof the pion mass χ = m π f π m u m d ( m u + m d ) + O (cid:18) m s (cid:19) = (cid:26) (77 . (75 . . (58)We used m π = 135 MeV and f π = 93 MeV. Together with the fits, Fig. 16,we can compute the corresponding quark masses. We convert them into MSmasses at 2 GeV, run at one-loop, in order to compare them more easily withother sources. Our results are m PV i ( µ = 0 . (cid:26) . . . . , (59) m MS i ( µ = 2 GeV) = (cid:26) . . . . . (60)The errors include an estimate from the 2-loop running. These masses comparewell with the particle data group masses [15], i.e. m u = 1 . − . m d =3 . − . m s = 70 −
130 MeV, and to the lattice masses [34], i.e. m u = 1 . m d = 4 . m s = 87(6) MeV.Very large volume simulations are expensive even in the IILM. We have seenin section 5 that the instanton density becomes independent of quark masses inthe chiral limit. Therefore, in the physical region of parameter space that weare considering the volume is directly proportional to the number of instantonsin the box. The computation of the fermionic determinant is the most costlypart of the simulations. In a naive implementation it would scale as O ( V ).However, in section 4.3 we were able to drastically reduce this cost to O ( V )by implementing fast update algorithms for the Cholesky decomposition. Theabsolute scale for the volumes is set by considerations regarding finite size effects.To have these under control the lightest propagating particle, in our case thepion, should fit into the box. For masses beyond the chiral limit the instantondensity increases, see again section 5. However, the pion mass increases too andit turns out that for the same level of control over finite size effects the systemsize can be smaller, i.e. larger quark masses are computationally cheaper.We have studied the thermodynamic limit on four volumes, in the range 2 . Lm π .
3. Even though the data has displayed a nice scaling with the volume, it34 -4 Λ V (2000 4000 6000 8000 10000 12000 14000 16000 m i n λ -5 -4 m m Figure 17: The simulations are performed for different simulation boxes. The average of thesmallest Dirac eigenvalue, h λ min i , is smaller than the quark masses for all but the smallestsimulation box. is important to check whether the thermodynamic limit was consistent. To thisend we will run large volume simulations, Lm π ∈ [2 . , . m u = 0 . m d = 0 . m s = 0 . ρ ∈ [0 . , . ρ = ¯ ρ ∞ + αV − . gives¯ ρ ∞ = 0 . − = 0 . n = 1 . − , and like thetopological susceptibility displays a nice thermodynamic limit.
7. Conclusion
With the discovery of the new non-trivial holonomy calorons [28, 30], [31],there is renewed interest in studying the role of non-trivial field configurations inQCD, especially their role in the confinement/deconfinement phase transition.A lattice based approach [21] is well suited for the pure gauge sector becauseit easily incorporates many-body instanton interactions in the classical action.However, the introduction of fermions will be plagued by the same problems that35 -4 Λ V (2000 4000 6000 8000 10000 12000 14000 16000 > < Q 200 to h N i ≈ ) -4 Λ V (5000 10000 15000 > ρ < ρ < Figure 19: Even for the largest volumes the mean instanton size ¯ ρ is still decreasing and doesnot seem to converge to a constant. We are rather lucky that, although the effect is clearlysystematic, the variation of ¯ ρ is small in absolute terms. O ( N ) to O ( N ); this is achieved by rewriting theupdates in a form suitable for fast matrix modifications. Secondly, we study thethermodynamic limit and monitor some bulk quantities to guarantee a consis-tent large volume extrapolation.The topological susceptibility is easy to compute in the IILM and has beenstudied extensively on the lattice. It represents a natural candidate to fix unitsin the IILM. For quenched simulations we found rather good agreement be-tween the IILM and the lattice in this way. We are, however, mainly interestedin the unquenched case. Instead of using the topological susceptibility directly,we decided to use the chiral condensate, extracted through the topological sus-ceptibility and chiral perturbation theory, to set units. The reason we decidedagainst a direct use of the topological susceptibility is that it depends stronglyon the quark masses. We found that the chiral condensate has a very weak de-pendence on the chiral limit, and use it to set the units in the unquenched sector.We achieve good agreement with previous work and also with experimental dataon the strong coupling constant α s . Using chiral perturbation theory, we areable to determine physical quark masses. These turn out to compare well withexperimental bounds and lattice simulations.Finally, we investigated the uncertainties introduced by the large volumeextrapolations, and found that our procedure, of bounding the volume in sucha way that the quark masses are smaller than the smallest Dirac eigenvalues,allows for a systematic thermodynamic limit.In a further publication we will use these input parameters to study theIILM at finite temperature. Acknowledgements We are very grateful for many informative discussions with P. Faccioli andE.P.S. Shellard. Simulations were performed on the COSMOS supercomputer(an Altix 4700) which is funded by STFC, HEFCE and SGI. This work wassupported by STFC grant PPA/S/S2004/03793 and an Isaac Newton TrustEuropean Research Studentship. Appendix A. Gluonic Interactions The ratio ansatz for an instanton–anti-instanton pair is defined by A aµ = − ¯ η aµν ∂ ν Π ( x, { x , ρ } ) + O ab η bµν ∂ ν Π ( x, { x , ρ } )1 + Π ( x, { x , ρ } ) + Π ( x, { x , ρ } ) , (A.1)where Π( x, { y, ρ } ) = ρ r and O = O t O , with O i the respective colour embed-dings. A global colour rotation has been performed to bring the gauge potential38nto this form, which is irrelevant since the action is gauge invariant. Instanton–instanton and anti-instanton–anti-instanton pairs differ by having either only ¯ η or η in the above formula. A brute force computation then gives F aµν F aµν = I + (Tr O t O + (¯ ηOη ) µνµν ) J + (¯ ηOη ) ρµρν I µν + (¯ ηOη ) µρνσ I µρνσ + ( ηO t Oη ) µρνσ J µρνσ + (¯ ηOη ) αµαρ (¯ ηOη ) βνβσ K µρνσ . (A.2)The different terms have the following form I = 4(1 + Π + Π ) [( ∂ µ ∂ ν Π )( ∂ µ ∂ ν Π ) + ( ∂ µ ∂ ν Π )( ∂ µ ∂ ν Π )] − + Π ) [( ∂ µ ∂ ν Π )( ∂ µ Π )( ∂ ν Π ) + ( ∂ µ ∂ ν Π )( ∂ µ Π )( ∂ ν Π )+2( ∂ µ ∂ ν Π )( ∂ µ Π )( ∂ ν Π ) + 2( ∂ µ ∂ ν Π )( ∂ µ Π )( ∂ ν Π )]+ 4(1 + Π + Π ) [3( ∂ µ Π ∂ µ Π )( ∂ µ Π ∂ µ Π ) + 3( ∂ µ Π ∂ µ Π )( ∂ µ Π ∂ µ Π )+3( ∂ µ Π ∂ µ Π ) + 3( ∂ µ Π ∂ µ Π ) + 2( ∂ µ Π ∂ µ Π )( ∂ µ Π ∂ µ Π )+( ∂ µ Π ∂ µ Π ) (cid:3) . (A.3) J = 2(1 + Π + Π ) ( ∂ µ Π ∂ µ Π )( ∂ µ Π ∂ µ Π ) . (A.4) I µν = 4(1 + Π + Π ) ( ∂ µ ∂ σ Π )( ∂ µ ∂ σ Π )+ 4(1 + Π + Π ) [( ∂ µ ∂ ν Π )( ∂ σ Π ∂ σ Π ) + ( ∂µ∂ ν Π )( ∂ σ Π ∂ σ Π ) − ∂ µ ∂ σ Π )( ∂ ν Π )( ∂ σ Π ) − ∂ µ Π )( ∂ σ Π )( ∂ ν ∂ σ Π ) − ∂ µ ∂ σ Π )( ∂ σ Π )( ∂ ν Π ) − ∂ µ Π )( ∂ ν ∂ σ Π )( ∂ σ Π )]+ 4(1 + Π + Π ) [ − ( ∂ µ Π )( ∂ ν Π )( ∂ σ Π ∂ σ Π ) − ( ∂ µ Π )( ∂ ν Π )( ∂ σ Π ∂ σ Π )+3( ∂ µ Π )( ∂ ν Π )( ∂ σ Π ∂ σ Π ) + 3( ∂ µ Π )( ∂ ν Π )( ∂ σ Π ∂ σ Π )+3( ∂ µ Π )( ∂ ν Π )( ∂ σ Π ∂ σ Π )] . (A.5) I µρνσ = 4(1 + Π + Π ) ( ∂ µ ∂ ν Π )( ∂ ρ ∂ σ Π )+ 8(1 + Π + Π ) [( ∂ µ Π )( ∂ ρ ∂ ν Π )( ∂ σ Π ) + ( ∂ µ Π )( ∂ ρ ∂ ν Π )( ∂ σ Π )]+ 8(1 + Π + Π ) ( ∂ µ Π )( ∂ ρ Π )( ∂ ν Π )( ∂ σ Π ) . (A.6)39 µρνσ = 2(1 + Π + Π ) ( ∂ µ Π )( ∂ ρ Π )( ∂ ν Π )( ∂ σ Π ) . (A.7) K µρνσ = 2(1 + Π + Π ) ( ∂ µ Π )( ∂ ρ Π )( ∂ ν Π )( ∂ σ Π ) . (A.8) Appendix A.1. Exact Interactions When computing the look-up tables, we use global translations and rota-tions in R to place one instanton at the origin and the partner at y ′ = R = p R µ R µ = | y I − y I | , where y i are the instanton centres. The rotation willreemerge in contractions of R µ with the colour structure, as we will now see.The relation between the position vector R µ and R ′ µ ≡ (0 , , , R ) is given bythe following rotation matrix R ′ µ = O tµν R ν , O µ = R µ R , (A.9)and the other components of the rotation matrix are irrelevant.Note that, with the choice of R ′ µ , the integrands are O (3) symmetric in thesubspace orthogonal to the 4-direction. Denoting the arguments of the ’t Hooftpotentials Π( x, { y, ρ } ) by x µ and ˜ x µ ≡ x µ − R µ , we can extract extract the R µ from the integrands with help of the following formulas, which we orderaccording to the tensor structure of the x µ -dependence on the integrand. Z x µ = O µ Z x ′ . (A.10) Z x µ x ν = δ µν Z x ′ + O µ O ν Z ( x ′ − x ′ ) . (A.11) Z x µ x ν x κ = ( δ µν O κ + δ κµ O ν + δ νκ O µ ) Z x ′ x ′ (A.12)+ O µ O ν O κ Z ( x ′ − x ′ x ′ ) . (A.13) Z x µ x ν x κ x δ = ( δ µν δ κδ + δ µκ δ νδ + δ µδ δ κν ) Z x ′ x ′ (A.14)+ ( δ µν O κ O δ + perm . ) Z ( x ′ x ′ − x ′ x ′ ) (A.15)+ O µ O ν O κ O δ Z ( x ′ − x ′ x ′ + 3 x ′ x ′ ) . (A.16)40erms with ˜ x can be constructed from these. Incidentally, splitting thedifferent integrands according to the above formulas is the most stable procedurenumerically. Taking into account the antisymmetry of the ’t Hooft symbols, weend up with the following integrands. I = 4(1 + Π + Π ) (cid:2) (Π ′′ ) + 3(Π ′ /r ) + (Π ′′ ) + 3(Π ′ / ˜ r ) (cid:3) − + Π ) (cid:20) ′′ (Π ′ ) + 2Π ′′ (Π ′ ) + x ˜ xr ˜ r (Π ′′ Π ′ Π ′ + Π ′ Π ′′ Π ′ ) (cid:21) + 4(1 + Π + Π ) (cid:20) ′ ) + 12(Π ′ ) + 8(Π ′ ) (Π ′ ) + 4( x ˜ xr ˜ r Π ′ Π ′ ) +12(Π ′ ) ( x ˜ xr ˜ r Π ′ Π ′ ) + 12( x ˜ xr ˜ r Π ′ Π ′ )(Π ′ ) (cid:21) . (A.17)Note that to achieve good numerical precision, we need to subtract theone-instanton integrands from the above before performing the numerical inte-gration. J = 2(1 + Π + Π ) (Π ′ ) (Π ′ ) . (A.18) I µν = δ µν ˜ I µµ + R µ R ν R ˜ I µν . (A.19)˜ I µµ = 4(1 + Π + Π ) (cid:20) x ′ r (Π ′′ − (Π ′ /r ))(Π ′ / ˜ r ) + x ′ ˜ r (Π ′ /r )(Π ′′ − (Π ′ / ˜ r ))+(Π ′ /r )(Π ′ / ˜ r ) + x ′ r ˜ r x ˜ xr ˜ r (Π ′′ Π ′′ − Π ′′ (Π ′ / ˜ r ) − (Π ′ /r )Π ′′ + (Π ′ /r )(Π ′ / ˜ r )) (cid:21) + 1(1 + Π + Π ) (cid:2) ′ /r )(Π ′ ) + (Π ′ ) (Π ′ / ˜ r ))+ x ′ r (4(Π ′′ − (Π ′ /r ))(Π ′ ) − ′ ) (Π ′ / ˜ r ))+ x ′ ˜ r (4(Π ′ ) (Π ′′ − (Π ′ / ˜ r )) − ′ /r )(Π ′ ) )+ x ′ r ˜ r ( − x ˜ xr ˜ r ((Π ′′ − (Π ′ /r ))(Π ′ ) + (Π ′ ) (Π ′′ − (Π ′ / ˜ r ))) − ′′ Π ′ Π ′ − ′ Π ′′ Π ′ ]+ 1(1 + Π + Π ) (cid:20) − x ′ r (Π ′ ) (Π ′ ) − x ′ ˜ r (Π ′ ) (Π ′ ) +12 x ′ r ˜ r Π ′ Π ′ ((Π ′ ) + (Π ′ ) + x ˜ xr ˜ r Π ′ Π ′ ) (cid:21) . (A.20)41 I µν = 4(1 + Π + Π ) (cid:20) x ′ − x ′ r (Π ′′ − (Π ′ /r ))(Π ′ / ˜ r )+ ( x ′ − R ) − x ′ ˜ r (Π ′ /r )(Π ′′ − (Π ′ / ˜ r ))+ x ′ ( x ′ − R ) − x ′ r ˜ r x ˜ xr ˜ r (Π ′′ Π ′′ − Π ′′ (Π ′ / ˜ r ) − (Π ′ /r )Π ′′ + (Π ′ /r )(Π ′ / ˜ r )) (cid:21) + 1(1 + Π + Π ) (cid:20) x ′ − x ′ r (4(Π ′′ − (Π ′ /r ))(Π ′ ) − ′ ) (Π ′ / ˜ r ))+ ( x ′ − R ) − x ′ ˜ r (4(Π ′ ) (Π ′′ − (Π ′ / ˜ r )) − ′ /r )(Π ′ ) )+ x ′ ( x ′ − R ) − x ′ r ˜ r ( − x ˜ xr ˜ r ((Π ′′ − (Π ′ /r ))(Π ′ ) + (Π ′ ) (Π ′′ − (Π ′ / ˜ r ))) − ′′ Π ′ Π ′ − ′ Π ′′ Π ′ ]+ 1(1 + Π + Π ) (cid:20) − x ′ − x ′ r (Π ′ ) (Π ′ ) − x ′ − R ) − x ′ ˜ r (Π ′ ) (Π ′ ) +12 x ′ ( x ′ − R ) − x ′ r ˜ r Π ′ Π ′ ((Π ′ ) + (Π ′ ) + x ˜ xr ˜ r Π ′ Π ′ ) (cid:21) . (A.21) I µρνσ = δ µν δ ρσ ˜ I µνµν + δ µν R ρ R σ R ˜ I µρµσ . (A.22)˜ I µνµν = 0 (analytically) . (A.23)˜ I µρµσ = 4(1 + Π + Π ) (cid:20) x ′ − x ′ r (Π ′′ − (Π ′ /r ))(Π ′ / ˜ r )+ ( x ′ − R ) − x ′ ˜ r (Π ′ /r )(Π ′′ − (Π ′ / ˜ r )) x ′ R ( r ˜ r ) (Π ′′ − (Π ′ /r ))(Π ′′ − (Π ′ / ˜ r )) (cid:21) − + Π ) (cid:20) x ′ − x ′ r (Π ′ ) (Π ′ / ˜ r ) + ( x ′ − R ) − x ′ ˜ r (Π ′ /r )(Π ′ ) + x ′ R ( r ˜ r ) ((Π ′′ − (Π ′ /r ))(Π ′ ) + (Π ′ ) (Π ′′ − (Π ′ / ˜ r ))) (cid:21) + 8(1 + Π + Π ) (cid:20) x ′ R ( r ˜ r ) (Π ′ ) (Π ′ ) (cid:21) . (A.24)˜ J µρνσ = δ µν R ρ R σ R + Π ) x ′ R ( r ˜ r ) (Π ′ ) (Π ′ ) . (A.25)42 K µρνσ = (cid:2) ( δ µν δ ρσ + δ µρ δ νσ + δ µσ δ ρν ) x ′ x ′ + δ µν R ρ R σ R (( x ′ − R ) x ′ − x ′ x ′ ) + δ ρσ R µ R ν R ( x ′ x ′ − x ′ x ′ )+ ( δ µρ R ν R σ R + δ µσ R ν R ρ R + δ νρ R µ R σ R + δ νσ R µ R ρ R )( x ′ ( x ′ − R ) x ′ − x ′ x ′ )+ R µ R ν R ρ R σ R ( x ′ ( x ′ − R ) + 3 x ′ x ′ − x ′ x ′ − ( x ′ − R ) x ′ − x ′ ( x ′ − R ) x ′ ) (cid:3) + Π ) r ˜ r ) (Π ′ ) (Π ′ ) . (A.26) Appendix A.2. Asymptotic Interactions As explained in the main text, the small separation asymptotic formulasget contributions which have the same functional form as those for the largeseparation asymptotics; the difference lies in the integration limit. We willtherefore start with the large separation formulas and leave the integrals explicit. Appendix A.2.1. Large Separation The upper integration limit z follows from variable substitution and has thethe following for integration over I , with I held fixed, z = 1 + Π ρ r . (A.27)Apart from the dependence of z on Π, the rational form of the ’t Hooftpotential allows for a complete factoring out of Π under the above mentionedvariable substitution. For the large separation formulas it is understood that z → ∞ because the initial integration variable r extends to infinity.The integral over I contains terms that do not mix the ’t Hooft potential Π and Π except for the denominators. At zeroth order in our expansion, theseterms can be transformed to exactly match the single instanton contributionsby exploiting scale invariance. Remembering that we actually subtract theone-instanton contributions to get the interactions, we can neglect these termsaltogether. We then end up with the following formulas. Z I = 72 π ρ ∂ µ Π ∂ µ Π(1 + Π) Z z s ds ( s + 1) + sym . (A.28) Z J = 16 π ρ ∂ µ Π ∂ µ Π(1 + Π) Z z s ds ( s + 1) + sym . (A.29)43 I µν = 16 π ρ ∂ µ ∂ ν Π(1 + Π) Z z s ds ( s + 1) − (cid:18) π ρ δ µν ∂ σ Π ∂ σ Π(1 + Π) + 8 π ρ ( ∂ µ Π)( ∂ ν Π)(1 + Π) (cid:19) Z z s ds ( s + 1) + sym . (A.30)At zeroth order, partial integration and the antisymmetry of the ’t Hooftsymbols can be used to simplify I µρνσ → + Π ) ( ∂ µ Π )( ∂ ρ Π )( ∂ ν Π )( ∂ σ Π ) , (A.31)with asymptotic behaviour Z I µρνσ = 16 π ρ δ µν ( ∂ ρ Π)( ∂ σ Π)(1 + Π) Z z s ds ( s + 1) + sym . (A.32) Z J µρνσ = 4 π ρ δ µν ( ∂ ρ Π)( ∂ σ Π)(1 + Π) Z z s ds ( s + 1) + sym . (A.33)For K µρνσ no ’t Hooft symbols can be used to exchange the index pairs( µ, ν ) ↔ ( ρ, σ ), and so we cannot simplify with a symmetry argument anymore. Z K µρνσ = 4 π ρ δ µν ( ∂ ρ Π )( ∂ σ Π )(1 + Π ) Z z s ds ( s + 1) + 4 π ρ δ ρσ ( ∂ µ Π )( ∂ ν Π )(1 + Π ) Z z s ds ( s + 1) . (A.34) Appendix A.2.2. Small Separation As explained in the main text, the small separation asymptotic formulas getcontributions from the large asymptotics. Also, in this case we have performeda global translation so that the instantons sit at ± R µ / 2. Therefore, in the largeseparation formulas we need to put z = ρ ( R/ .We now turn to the proper small separation asymptotic formulas that encodethe repulsion through the gauge singularity. We will again introduce an explicitupper limit for the integrals; abusing notation we will use the same letter asbefore, but here the meaning becomes z = R ρ + ρ , z i = R ρ i . (A.35)To derive these formulas, we approximate the arguments x µ ± R µ / → x µ .We have, therefore, explicitly restored O (4) symmetry, which can be exploitedto set several integrals to zero. Eventually, we arrive at44 I = 384 π (cid:20) ρ + ρ ( ρ + ρ ) Z z dss ( s + 1) − (cid:18) ρ ρ ( ρ + ρ ) + 2 ρ + ρ ( ρ + ρ ) (cid:19) Z z dss ( s + 1) + ρ + ρ + ρ ρ + ρ ρ + ρ ρ ( ρ + ρ ) Z z dss ( s + 1) − Z z s dss ( s + 1) − Z z s dss ( s + 1) (cid:21) . (A.36) Z J = 64 π ρ ρ ( ρ + ρ ) Z z dss ( s + 1) . (A.37) Z I µν = δ µν (cid:20) π ρ ρ ( ρ + ρ ) Z z dss ( s + 1) − π ρ ρ ( ρ + ρ ) Z z dss ( s + 1) +32 π ρ ρ + 3 ρ ρ + 3 ρ ρ ( ρ + ρ ) Z z dss ( s + 1) (cid:21) . (A.38) Z I µρνσ = δ µν δ ρσ (cid:20) − π ρ ρ ( ρ + ρ ) Z z dss ( s + 1) +32 π ρ ρ ( ρ + ρ ) Z z dss ( s + 1) (cid:21) . (A.39) Z J µρνσ = 0 . (A.40) Z K µρνσ = 83 π ( δ µν δ ρσ + δ µρ δ νσ + δ µσ δ νρ ) ρ ρ ( ρ + ρ ) Z z dss ( s + 1) . (A.41) Appendix B. Fermionic Interactions The Dirac overlap matrix elements are given by T IA = Z d x π ρ I ρ A 12 Tr( U τ + β ) I β . (B.1)Note that Tr( U τ + β ) ≡ iu β is the colour four-vector used for instance in [39].After some straightforward algebra, we find that I β has the following form I β = − I + Π A )(1 + Π I ) / (1 + Π A ) / (cid:18) Π A I ( ∂ µ Π I ∂ µ Π I ) ∂ β Π A + ( ∂ µ Π A ∂ µ Π A ) ∂ β Π I (cid:19) . (B.2)45 ppendix B.1. Exact Interactions Using the same rotations (A.9) as for the gluonic interactions to marry thespace-time with the colour indices, we get I β = R β R − I + Π A )(1 + Π I ) / (1 + Π A ) / (cid:26) x ′ r Π ′ I (Π ′ A ) + x ′ − R ˜ r (Π ′ I ) Π ′ A Π A I (cid:27) . (B.3) Appendix B.2. Asymptotic InteractionsAppendix B.2.1. Large Separation In order to get rather simple formulas, we make the following additionalsimplification1 + Π I + Π A → (cid:26) I : Integration over Π I A : Integration over Π A . (B.4)Given these further assumption, we can proceed as for the gluonic interac-tions. Finally, caution needs to be taken in the case of an anti-instanton becauseit sits at − R µ so that ∂ β Π A generates an extra minus sign. Z I β = 8 π ρ I Π A ∂ β Π A (1 + Π A ) / Z z I s ds ( s + 1) / − π ρ A ∂ β Π I (1 + Π I ) / Z z A s ds ( s + 1) / . (B.5) Appendix B.2.2. Small Separation At zeroth order, i.e. x µ ± R µ / → x µ , the contribution to I β vanishes becauseof O (4) symmetry. It turns out that the large separation asymptotics falls offtoo quickly as R → 0. However, this is not important because in this regimethe gluonic interaction is dominant. Appendix C. Cholesky decomposition update In this appendix we look in detail at how the structure suitable for theCholesky decomposition update comes about. We will also see that insertioncan be performed faster whereas deletions will be the most costly.46 ppendix C.1. Canonical Moves Upon updating instanton I , we have that T → T +∆ T , with (∆ T ) ij = δ iI ξ ∗ j .This induces the following changes( T † T ) ij → ( T † T ) ij + T † iI ξ ∗ j + ξ i T Ij + ξ i ξ ∗ j , (C.1) ψ i ≡ T ∗ Ii , (C.2) φ i ≡ ξ i + ψ i , (C.3) T † T → T † T + φφ † − ψψ † . (C.4)( T T † ) ij → ( T T † ) ij + δ iI ξ ∗ k T † kj + T ik ξ k δ Ij + δ iI δ jI , (C.5) ψ i ≡ | ξ | ( T ξ ) i , (C.6) φ i ≡ δ Ii | ξ | + ψ i , (C.7) T † T → T † T + φφ † − ψψ † . (C.8)Changes in an anti-instanton will have analogous formulas. Appendix C.2. Insertion We focus on inserting an instanton. Insertion of an anti-instanton is thensimilar. Since in the code we always add an instanton at the end of the arrays,an insertion corresponds to adding a row to T . T → (cid:18) Tξ † (cid:19) , (C.9) T T † → (cid:18) Tξ † (cid:19) (cid:0) T † ξ (cid:1) = (cid:18) T T † T ξξ † T † ξ † ξ (cid:19) . (C.10)On the level of the Cholesky decomposition L → = (cid:18) L χ † (cid:19) , (C.11) D → = (cid:18) D d (cid:19) , (C.12) LDL † → = (cid:18) LDL † LDχχ † DL † χ † Dχ + d (cid:19) . (C.13)Remembering that the insertion also adds a mass term in the diagonal, wehave to solve the following system (cid:26) LDχ = T ξd = ξ † ξ + m − χ † Dχ , (C.14)47hich can be solved in O ( N ) by using backsubstitution. The case for T † T issimply given by T † T → (cid:0) T † ξ (cid:1) (cid:18) Tξ † (cid:19) = T † T + ξξ † , (C.15)which is a rank-1 update. Appendix C.3. Deletion We focus again on an instanton. Deletion will be a two step process. Wefirst delete the last instanton and then swap it with that instanton that has beenselected for deletion. The swapping is similar to a canonical move, where now ξ is given by the difference between the last instanton and the selected instanton.The proper deletion part is given by T T † → = (cid:18) T T † 00 0 (cid:19) , (C.16) L → = (cid:18) L 00 1 (cid:19) , (C.17) D → = (cid:18) D 00 0 (cid:19) . (C.18)The T † T part is again simply related to a rank-1 update because, uponrearranging the result from the insertion part, we get T † T → T † T − ξξ † . (C.19) Appendix D. MS to PV Operators with anomalous dimensions run, and for mass independent renor-malisation prescriptions they depend on the scheme. The IILM makes predic-tions within a subtraction scheme that uses Pauli-Villars regularisation. How-ever, the lattice results have been quoted in MS, and so we need to work outthe relation between the two.It is not hard to convince ourselves that the quark masses run inverselyto the chiral condensate: note that chiral perturbation theory, (58), relatesthe topological susceptibility to the pion mass and decay constant, which arephysical quantities; it also relates the chiral condensate and the quark massesto the topological susceptibility through (52), and therefore the renormalisationscheme dependence must exactly cancel among the two.Eventually, we will also relate the quark masses of the two schemes, and sohere we will focus on mass renormalisation. We will only work at one-loop ,i.e. we will need to evaluate (Fig. D.20) in both schemes. Maintaining manifest gauge-invariance in Yang-Mills theories using Pauli-Villars regular-isation is not straightforward beyond one-loop. Figure D.20: Feynman diagram needed to compute the difference between the MS and PVscheme at one-loop. This is a textbook computation, [38]. After subtracting off the divergences,we end up withΣ PV = α s π C (3) (cid:26) − m + 14 p/ + Z dx (2 m − (1 − x ) p/ ) ln µ xm − x (1 − x ) p (cid:27) , (D.1)andΣ MS = α s π C (3) (cid:26) − m + 12 p/ + Z dx (2 m − (1 − x ) p/ ) ln µ xm − x (1 − x ) p (cid:27) . (D.2)Relating both through the pole mass and using that, in our case C (3) = 4 / m MS = m PV (1 − α s π 53 ) . (D.3) References [1] A. Armoni, M. Shifman, and G. Veneziano. QCD Quark Condensate fromSUSY and the Orientifold Large-N Expansion. Phys. Lett. B , 579:384–390,2004.[2] A. A. Belavin, A. M. Polyakov, A. S. Schwartz, and Yu. S. Tyupkin. Pseu-doparticle solutions of the Yang-Mills equations. Phys. Lett. B , 59:85–87,1975.[3] F. Bruckmann, C. Gattringer, E.M. Ilgenfritz, M. M¨uller-Preussker,A. Schafer, and S. Solbrig. Quantitative comparison of filtering methodsin lattice QCD. Eur. Phys. J. A , 33:333–338, 2007.[4] K. Chetyrkin, J. K¨uhn, and M. Steinhauser. RunDec: a Mathematica pack-age for running and decoupling of the strong coupling and quark masses. Comput. Phys. Commun. , 133:43–65, 2000.495] T. Chiu, S. Aoki, H. Fukaya, S. Hashimoto, T. Hsieh, T. Kaneko, H. Mat-sufuru, J. Noaki, K. Ogawa, T. Onogi, and N. Yamada. Topological sus-ceptibility in 2-flavor lattice QCD with fixed topology. PoS Lattice 2007 ,068, 2007.[6] T. Chiu, S. Aoki, S. Hashimoto, T. Hsieh, T. Kaneko, H. Matsufuru,J. Noaki, T. Onogi, and N. Yamada. Topological susceptibility in (2+1)-flavor lattice QCD with overlap fermion. PoS Lattice 2008 , 072, 2008.[7] T. Chiu, T. Hsieh, and P. Tseng. Topological susceptibility in (2+1) flavorslattice QCD with domain-wall fermions. Phys. Lett. B , 671:135–138, 2008.[8] M.C. Chu, J. Grandy, S. Huang, and J. W. Negele. Evidence for the Role ofInstantons in Hadron Structure from Lattice QCD. Phys.Rev. D , 49:6039–6050, 1994.[9] M. Cristoforetti, P. Faccioli, M.C. Traini, and J.W. Negele. Exploringthe Chiral Regime of QCD in the Interacting Instanton Liquid Model. Phys.Rev. D , 75:034008, 2007.[10] D. Diakonov, M. Polyakov, and C. Weiss. Hadronic matrix elements ofgluon operators in the instanton vacuum. Nucl. Phys. , B461:539–580, 1996.[11] G.V. Dunne, J. Hur, C. Lee, and H. Min. Calculation of QCD InstantonDeterminant with Arbitrary Mass. Phys.Rev. D , 71:085019, 2005.[12] S. Durr, Z. Fodor, C. Hoelbling, and T. Kurth. Precision study of theSU(3) topological susceptibility in the continuum. JHEP , 04:055, 2007.[13] D.I. Dyakonov and V.Yu. Petrov. Instanton-based vacuum from the Feyn-man variational principle. Nucl. Phys. B , 245:259–292, 1984.[14] D.I. Dyakonov and V.Yu. Petrov. A theory of light quarks in the instantonvacuum. Nucl. Phys. B , 272:457–489, 1986.[15] C. Amsler et al. Phys. Lett. B , 667:1, 2008.[16] K. Hagiwara et al. Quantum Chromodynamics. Phys. Rev. D , 66:010001–1,2002.[17] P. Faccioli. Strong CP breaking and quark-antiquark repulsion in QCD, atfinite θ . Phys. Rev. D , 71:091502, 2005.[18] P. Faccioli, D. Guadagnoli, and S. Simula. The Neutron Electric DipoleMoment in the Instanton Vacuum: Quenched Versus Unquenched Simula-tions. Phys.Rev.D , 70:074017, 2004.[19] D. Frankel and B Smit. Molecular Simulation: From Algorithm to Appli-cations . Academic Press, second edition edition, 2002.5020] A. M. Garc´ıa-Garc´ıa and J. C. Osborn. Chiral phase transition and An-derson localization in the Instanton Liquid Model for QCD. Nucl. Phys.A , 770:141–161, 2006.[21] P. Gerhold, E.M. Ilgenfritz, and M. M¨uller-Preussker. An SU (2) KvBLLcaloron gas model and confinement. Nucl. Phys. B , 760:1–37, 2007.[22] P. Gill, G. Golub, W. Murray, and M. Saunders. Methods for ModifyingMatrix Factorizations. Mathematics of Computations , 28:505–535, 1974.[23] M. G¨ockeler, R. Horsley, A. Irving, D. Pleiter, P. Rakow, G. Schierholz, andH. St¨uben. A Determination of the Lambda Parameter from Full LatticeQCD. Phys. Rev. D , 73:014513, 2006.[24] B. Grossman. Zero energy solutions if the Dirac equation in an N-pseudoparticle field. Phys. Lett. A , 61:86–88, 1977.[25] B.J. Harrington and H.K. Shepard. Periodic Euclidean solutions and thefinite-temperature Yang-Mills gas. Phys. Rev. D , 17:2122–2125, 1978.[26] A. Hasenfratz and P. Hasenfratz. The scales of euclidean and hamiltonianlattice QCD. Nucl. Phys. B , 193:210–220, 1981.[27] H. Kleinert. Path Integral in Quantum Mechanics, Statistics, PolymerPhysics,and Financial Markets . World Scientific, 2006.[28] T.C. Kraan and P. van Baal. Exact T -duality between Calorons and Taub-NUT spaces. Phys.Lett. B , 428:268–276, 1998.[29] T.C. Kraan and P. van Baal. Monopole Constituents inside SU ( N )Calorons. Phys.Lett. B , 435:389–395, 1998.[30] T.C. Kraan and P. van Baal. Periodic Instantons with non-trivial Holon-omy. Nucl.Phys. B , 533:627–659, 1998.[31] K. Lee and C. Lu. SU (2) Calorons and Magnetic Monopoles. Phys.Rev.D , 58:025011, 1998.[32] F. Lenz, J.W. Negele, and M. Thies. Confinement from merons. Phys. Rev.D , 69:074009, 2004.[33] H. Leutwyler and A. Smilga. Spectrum of Dirac operator and role of wind-ing number in QCD. Phys. Rev. D , 46:5607–5632, 1992.[34] Q. Mason, H. Trottier, R. Horgan, C. Davies, and G. Lepage. High-precisiondetermination of the light-quark masses from realistic lattice QCD. Phys.Rev. D , 73:114501, 2006.[35] G. Munster and C. Kamp. Distribution of instanton sizes in a simplifiedinstanton gas model. Eur. Phys. J. , C17:447–454, 2000.5136] J.W. Negele, F. Lenz, and M. Thies. Confinement from Instantons orMerons. Nucl. Phys. Proc. Suppl. , 140:629, 2005.[37] M.A. Nowak, M. Rho, and I. Zahed. Chiral Nuclear Dynamics . WorldScientific, 1996.[38] M. Peskin and D. Schroeder. An Introduction to Quantum Field Theory .Perseus Books, 1995.[39] T. Sch¨afer and E.V. Shuryak. Interacting instanton liquid in QCD at zeroand finite temperatures. Phys.Rev. D , 53:65226542, 1996.[40] E.V. Shuryak. Toward the quantitative theory of the instanton liquid (I).Phenomenology and the method of collective coordinates. Nucl. Phys. B ,302:559–579, 1988.[41] E.V. Shuryak. Toward the quantitative theory of the instanton liquid (II).The SU (2) gluodynamics. Nucl. Phys. B , 302:574–598, 1988.[42] E.V. Shuryak. Toward the quantitative theory of the instanton liquid (III).Instantons and light fermions. Nucl. Phys. B , 302:599–620, 1988.[43] E.V. Shuryak. Instantons in QCD (I). Properties of the ”instanton liquid”. Nucl. Phys. B , 319:521–540, 1989.[44] E.V. Shuryak. Instantons in QCD (II). Correlators of pseudoscalar andscalar currents. Nucl. Phys. B , 319:541–569, 1989.[45] E.V. Shuryak and J.J.M. Verbaarschot. QCD Instantons at finite temper-ature. Nucl. Phys. B , 364:255–282, 1991.[46] E.V. Shuryak and J.J.M. Verbaarschot. Screening of the topological chargein a correlated instanton vacuum. Phys. Rev. D , 52:295–306, 1995.[47] G. ’t Hooft. Computation of the quantum effects due to a four-dimensionalpseudoparticle. Phys.Rev. D , 14:3432–3448, 1976.[48] J.J.M. Verbaarschot. Streamlines and conformal invariance in Yang-Millstheories. Nucl. Phys. B , 362:33–53, 1991.[49] J. Vermaseren, S. Larin, and T. Ritbergen. The 4-loop quark mass anoma-lous dimension and the invariant quark mass. Phys. Lett. B , 405:327–333,1997.[50] M. Wagner. Classes of confining gauge field configurations. Phys.Rev. D ,75:016004, 2007.[51] M. Wagner. Fermions in the pseudoparticle approach. Phys. Rev. ,D76:076002, 2007. 5252] Z.Q. Wang, X.F. Lu, and F. Wang. Dilute liquid of instanton and itstopological charge dominate the QCD vacuum. AIP Conf.Proc , 865:242–247, 2006.[53] Olivier Wantz. The topological susceptibility from grand canonical simula-tions in the interacting instanton liquid model: strongly associating fluidsand biased Monte Carlo. Nucl. Phys. B , 829:91–109, 2010.[54] Olivier Wantz and E. P. S. Shellard. The topological susceptibility fromgrand canonical simulations in the interacting instanton liquid model: chi-ral phase transition and axion mass.