Università degli Studi di Parma Dipartimento di Fisica M. Melloni Dottorato di Ricerca in Fisica XXIII ciclo Numerical Stochastic Perturbation Theory for Lattice QCD: computation of renormalization constants and prospects for the study of the Dirac operator spectrum Supervisore: Dr. Francesco Di Renzo Candidato: Michele Brambilla CONTENTS Introduction iii Motivations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi 1 Numerical Stochastic Perturbation Theory 1 1.1 Stochastic perturbation theory . . . . . . . . . . . . . . . . . . 2 1.2 Stochastic quantization for gauge theories . . . . . . . . . . . 4 1.3 Lattice Langevin equation for Wilson action . . . . . . . . . . 6 1.4 NSPT for LGT . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2 An environment for NSPT 13 2.1 Hide perturbation theory . . . . . . . . . . . . . . . . . . . . . 14 2.2 Multilayer code . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.3 Fields and parallelization . . . . . . . . . . . . . . . . . . . . . 22 3 Renormalization coefficients 27 3.1 RI’-MOM scheme . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.1.1 Method . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.1.2 Deep into RC’s . . . . . . . . . . . . . . . . . . . . . . 32 3.2 Matching αTLS to αMS . . . . . . . . . . . . . . . . . . . . . . 35 3.2.1 Method . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.2.2 Computation . . . . . . . . . . . . . . . . . . . . . . . 38 3.2.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.3 Chiral limit . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 3.4 Hypercubic Taylor expansion and continuum limit . . . . . . . 44 3.5 Infinite volume approximation . . . . . . . . . . . . . . . . . . 48 i CONTENTS 3.6 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.6.1 Strategies . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.6.2 The simplest measurement: critical mass . . . . . . . . 52 3.6.3 Field renormalization coefficients . . . . . . . . . . . . 55 3.6.4 Current renormalization coefficients . . . . . . . . . . . 56 3.6.5 Conclusions, i.e. summing the series . . . . . . . . . . 65 4 The Dirac operator spectrum 69 4.1 The Dirac spectrum and Perturbation Theory . . . . . . . . . 71 4.2 The Dirac Spectrum: what to look at . . . . . . . . . . . . . . 71 4.3 The Dirac Spectrum in NSPT . . . . . . . . . . . . . . . . . . 73 4.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 4.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 Conclusions 81 Ringraziamenti 83 ii Introduction Quantum field theory (QFT) has proven to give a correct picture of elemen- tary particles and their interactions. While elementary "particles” (leptons and quarks) are the building blocks of matter, the interaction are described in the theory in terms of gauge fields. The most common approach to the theory is the path integral formalism, where the expectation value of a generic observable O is given by 〈Ω| O |Ω〉 = ∫ DADψ¯DψO[A, ψ¯, ψ]e−SE [A,ψ¯,ψ]∫ DADψ¯Dψe−SE [A,ψ¯,ψ] (1) where A is the gauge field, ψ are the fermionic fields and SE [A, ψ¯, ψ] the Euclidean action. The interpretation of such a formula is the following: the expectation value of the observable is given by the average over all the field configuration, where the weight is given by the exponential of the action. Since the weight is given by e−S most important contribution are given by the path which minimizes the action, that is where δS = 0. (2) This is indeed the principle of least action which leads to classical equation of motion, and quantization is given by the fluctuation around the classical vacua. The path integral formalism (1) reveals a deep resemblance with statisti- cal mechanic, where the partition function is given by Z = ∫ DADψ¯Dψe−SE [A,ψ¯,ψ]. (3) This allows to approach the study of a quantum field system with the tech- niques of statistical mechanic. iii Introduction As it is well known, QFT shows a lot of subtleties, in primis UV diver- gences: the naive formulation gives infinities for almost all the observables. In order to avoid such kind of problems a regulator is required. Nevertheless the regulator is not physical, and a procedure that allows to remove the reg- ulator keeping quantities finite should be implemented. Such a mechanism is the renormalization procedure to which we will come back later. Different regularizations have been used over the years. One of the possible regulator is the lattice. Lattice has proven not only to be a useful regulator, but a viable alterna- tive to perturbative computations on continuum enabling the study of QCD (and field theory in general) also at strong coupling, where perturbation the- ory does not hold. From a practical point of view, this is implemented via MonteCarlo simulations. This does not mean that one can not approach e.g. lattice QCD from the point of view of perturbation theory. One has to compute and all the possible Feynman graphs contributing at a given order to a process as in continuum perturbation theory. The regulator appears in terms of the lattice spacing a and is removed by means of renormalization. A drawback of this approach in the case of lattice QCD is the appearance of a huge number of Feynman graphs. This comes out of decompactifing the Wilson formulation of gauge theories. These (irrelevant) contributions are in the end lattice artifacts, but they result in making very difficult to compute the overall finite contribution: stated in a very naive way, a term proportional to a can nevertheless diverge, resulting in a finite contribution in the continuum limit. This makes perturbative computations on the lattice a very hard task: most of the result are known only up to one loop. MonteCarlo simulations are based on the statistical interpretation of path integral in the Euclidean formulation. One has to construct the correct form of an observable in terms of field content, generate a quite large number of field configurations according to a correct MonteCarlo procedure and finally measure the observable. In extracting all the statistical information one has to work out results and errors in terms of means, standard deviations, au- tocorrelations and so on. Unlike perturbation theory, this approach works also in the case of strong coupling, giving access to information to the low energy regime of QCD. While the general picture of the renormalization pro- cess which is at work in a lattice simulation is quite simply stated, there are a lot of technical subtleties which in the end have to do with the same subtleties which are inherent to every renormalization scheme for a field the- ory. In practice, one has to face the finiteness of the lattice spacing and of the lattice volume. This makes very difficult to perform simulations at the iv physical point. In all this framework, a key role is played by the generalized couplings of the theory: in the end, one is probing a complicated statistical system which has a non trivial phase structure. In particular, simulating QCD at the physical value of quark mass is very difficult. Numerical Stochastic Perturbation Theory (NSPT) is something which lies in between the two worlds. One could say that it allows perturbative computations of quantities by means of a MonteCarlo simulation. In this spirit, it requires statistical treatment of data as any MonteCarlo simulation. The range of applicability of NSPT for Lattice QCD is from a formal point of view broad, but in practice one has to expect it is limited to the high energy region, as any QCD perturbative computation. Motivations This thesis accounts for my research project in NSPT. From a technical point of view, it actually covers three quite different fields. As in any numerical investigation, my research project required an ex- tensive programming work. My PhD activity was the chance for what is technically known as a refactoring of Parma group LGT (and in particular NSPT) codes. This was in particular true in a phase in which the group was moving to the usage of new multi-cores architectures and (even more) of a new parallel platform (Aurora). As it is often the case for such a work, there is no obvious way to account for this in a thesis without letting this numerical work actually take over Physics. The solution which I chose is to give a brief methodological account of the package which is the result of my activities (PRlgt). I would regard the computation of three loop renormalization constants as the core business of this work. Perturbative versus non-perturbative renor- malization has been an issue for quite a long time. As a matter of fact, such a comparison has always been limited by the low order at which perturbative results were available. The ambitious goal of my work was to get to three loop results for quark bilinears (in a given regularization) with a fair account of all the systematics. This in turn enables us to better assess the level of confidence of non-perturbative results. The bottom line is that there are in general more indeterminations in non-perturbative determinations than it is usually stated: truncation errors are not the end of the story. As a last subject, the work pins down the prospects for an NSPT study of the Dirac operator spectrum. The spirit is to investigate the pattern of chiral symmetry breaking from Banks-Casher relation: can one inspect the v Introduction reshuﬄing of Dirac operator eigenvalues due to color interactions? This part of the work has by far a different status with respect to renormalization constants computation. While I put the technical basis for a computation that has not even been attempted before with traditional techniques, there is a long way to go to get quantitative results. I try to give a qualitative account of what emerges from one loop bare perturbation theory: eigenvalues do repel each other and the spectrum is well rearranged with respect to the free field. Next step will be the computation of the relevant renormalization constant (which is currently work in progress and is not reported at all). In the end, I think it emerges that (not without a lot of technical work. . . ) NSPT can both provide phenomenologically relevant high loop computations and enable to embark in quite non-standard computations. Outline Chapter 1 introduces the theoretical environment we are dealing with. We will introduce Stochastic Quantization, a quantization procedure which offers an alternative to canonical quantization by letting the fields stochastically evolve by means of a Langevin equation. The technique can used to carry out a perturbative approach which can be proved to converge to the per- turbation theory of canonical quantization. Using the technique to quantize gauge theories requires particular care, because of the appearance of zero modes which give rise to divergences. Moreover, in the case of a numeri- cal approach a gauge fixing is required. After the general introduction to Stochastic Quantization of gauge theories and Stochastic Perturbation The- ory we will be ready to introduce Numerical Stochastic Perturbation Theory, the framework of our computations. We will introduce the general algorithm and show how the inclusion of dynamical fermions is easier than in the case of nonperturbative simulations. Chapter 2 gives a brief account of the concrete computational side of my research. During these three years I implemented the package used to obtain all the result which follows. The spirit of the package is to have a general, flexible environment, easy to read and to use, with a modular implementation, ready for parallelization. In particular, we will see how the features of C++ can be exploited to hide the perturbative side of the program and a structure based on classes allows the writing of operations in a “human form”. The main classes will be presented, together with just a few examples of operators. In the end we will see how data can be rearranged in order to obtain “ready to parallelize” structures, without increasing the complexity of the (high level) programming task. vi Chapter 3 contains the results obtained by applying NSPT to the com- putation of three loop renormalization constants of quark bilinears in the regularization defined by tree level Symanzik gauge/Wilson fermions. Par- ticular care is devoted to the analysis of the errors: are truncation errors in perturbative determination the main source of indetermination? Does trun- cation errors really make non-perturbative determinations preferred with re- spect to perturbative ones? After describing the RI’MOM scheme we are working with we derive the explicit form of the Z of quark bilinears. In order to carry out our computation to three loop, a two loops matching be- tween the MS scheme and the lattice (for tree level Symanzik gauge/Wilson fermions) is required. We will perform this matching thanks to the interme- diate matching with the potential scheme. We then describe how we can take into account different sources of errors: how we can stay in the chiral regime and approach continuum and infinite volume limit. Finally, a comparison to non-perturbative results is discussed. In chapter 4 we discuss the prospects for the NSPT study of the spec- trum of Dirac operator. Via the Banks-Casher relation, this can be related to chiral symmetry breaking, a key feature of QCD. We will introduce a perturbative approach to the study of the spectrum of the Dirac operator, and qualitatively discuss the first, one loop results obtained in a naive ap- proach, which only accounts for bare PT. Next step will be the computation of the relevant renormalization factor (which is not discussed). The spirit of this chapter is not to show any definitive result, but to introduce the envi- ronment developed for a non standard computation in Numerical Stochastic Perturbation Theory. vii CHAPTER 1 Basics of Numerical Stochastic Perturbation Theory Stochastic Quantization, first introduced by Parisi and Wu [1], provides a viable alternative to canonical quantization. From a practical point of view, a quantum field theory is solved once one can compute any correlation function of the form 〈φ(x1) . . . φ(xn)〉 = 1 Z ∫ Dφ e−S[φ]φ(x1) . . . φ(xn) (1.1) To formulate Stochastic Quantization, one introduces a new degree of free- dom in the system, the “stochastic time” t, so that φ(x)→ φη(x, t) (1.2) This can be thought of as a time since the system “evolves in it” according to the Langevin equation ∂ ∂t φη(x, t) = − δS δφη(x, t) + η(x, t). (1.3) The equation is made out of a drift term (given by the equations of motion) plus a white noise η(x, t) 〈η(x, t)〉η = 0 〈η(x, t)η(x′, t′)〉η = 2δ(x− x′)δ(t− t′) (1.4) The Gaussian distribution of the noise is fully specified by 〈. . .〉η = ∫ dP [η] . . . dP [η] = [Dη] exp ( −1 2 ∫ dxdt η2(x, t) ) . (1.5) 1 Numerical Stochastic Perturbation Theory The main assertion of stochastic quantization is the equivalence with canon- ical quantization: lim t→∞ 〈φ(x1, t) . . . φ(xn, t)〉η = 〈φ(x1) . . . φ(xn)〉. (1.6) This can be stated in a very simple way: in the limit t → ∞ the stochastic system reaches an equilibrium such that equal time correlation functions approach Green functions of quantum theory. Put in another form, eq. (1.6) states that the equilibrium distribution of the process converges to the probability distribution which defines the functional integral. Stochastic quantization has been widely used as an al- ternative to MonteCarlo simulations of quantum field theory (see for exam- ple [2] [3]). A formal proof of such equivalence can go through the so called Fokker- Plank equation, which describes the time evolution of the probability density function for a stochastic process. It can be indeed shown that the Langevin equation (1.3) can be traded for an equation for the quantity P [φ, t] defined by P [φ, t] = 〈 ∏ x δ(φ(x)− φη(x, t))〉η (1.7) The relevant equation is the Fokker-Plank equation P˙ [φ, t] = ∫ dx δ δφ(x) ( δS δφ(x) + δ δφ(x) ) P [φ(x), t], (1.8) It is important to stress that since the Langevin equation (1.3) describes a diffusive process, the action involved must be a smooth function. 1.1 Stochastic perturbation theory One can obtain a perturbative formulation of stochastic quantization observ- ing that the action can be written S[φ] = S0[φ] + SI [φ, g], where S0[φ] is the free (non interacting) term and the interaction is contained in the SI [φ, g] term alone. Here g is the coupling of the theory. For the free fields the solu- tion of Langevin equation exists, and can be written as a Green function. In the case of the full theory the solution is obtained in form of a power series in the coupling. One starts expanding the stochastic fields φη(x, t) = φη(x, t) (0) + ∑ n>0 gnφη(x, t) (n) (1.9) 2 Stochastic perturbation theory and plugging this expansion into the Langevin equation (1.3) obtains a hi- erarchy of equations which can be truncated at a given order. It is straight- forward to notice that the noise only only appears in the free term, since it doesn’t depend on the coupling. The stochastic effect at higher orders enters due to the dependence on lower orders. We only sketch the proof of the equivalence between Stochastic Quanti- zation Theory and standard quantization in perturbation theory. A rigorous proof can be found in [4]. Relying on the Fokker-Plank formalism (1.8) the proof can be summarized as in following: • P [φ(x), t] is expanded in power series in the coupling P [φ(x), t] = P [φ(x), t](0) + ∑ n>0 gnP [φ(x), t](n) (1.10) so that the Fokker-Plank equation turns into a hierarchy of equations; • one can prove that at tree level an equilibrium distribution Peq[φ(x)](0) exists P [φ(x), t](0) t→∞−−−→ Peq[φ(x)](0) = Z−10 e−S0[φ] (1.11) which is the free path integral measure; • in a weak sense the same holds for all the perturbative orders: P [φ(x), t](n) t→∞−−−→ Peq[φ(x)](n); (1.12) • the various Peq[φ(x)](n) are related by a set of equations in which one can recognise the Swinger-Dyson equations. Since the solutions of the Schwinger-Dyson equations are unique in perturbation theory, one has recovered a standard perturbative solution for the quantization proce- dure. We have just sketched a perturbative proof for Stochastic Quantization to hold in terms of an asymptotic distribution. It does not come as a surprise that a practical implementation of perturbation theory comes out of Langevin equation. This is known as Stochastic Perturbation Theory [1]. It can be derived in the most direct way by going back to the plugging of (1.9) into (1.3). The resulting tower of equations can be solved at any finite order. The contributions which arise order by order can be represented by some stochas- tic diagrams, similar to Feynman graphs. Numerical Stochastic Perturbation Theory (NSPT) [5] is simply a numerical implementation of what we have 3 Numerical Stochastic Perturbation Theory just described: one simply put the tower of stochastic differential equations on a computer and integrate them, eventually computing observables order by order at asymptotic times, in the spirit of (1.6). We will later add a few details of NSPT, specifying to the case of the Wilson action. For the moment it is important to emphasize that Numerical Stochastic Perturbation Theory is not a traditional MonteCarlo simulation, which is most often a discontin- uous process: it is instead the integration of a tower of stochastic differential equations. 1.2 Stochastic quantization for gauge theories Since NSPT is a numerical implementation of Stochastic Perturbation The- ory, the existence of the limit distribution is not the only concern: the prop- erties of convergence also matter. Since we are interested in the case of gauge theories, we need to devote some attention to the stochastic quantization of gauge fields. In canonical quantization the redundant degrees of freedom connected to gauge invariance should be properly handled to avoid infinities. In stochastic quantization (and as consequence in NSPT as well), something similar happens. We start recalling the action in the case of pure gauge S = 1 4 ∫ d4xF aµν(x)F a µν(x), (1.13) where F aµν(x) = ∂µA a ν(x) − ∂νAaµ(x) − gfabcAbµ(x)Acν(x). If we consider the Langevin equation (1.3) in case of gauge fields, we have to substitute the derivative with respect to the field with a covariant derivative. Then ∂ ∂t Aaµ(x, t) = Dabν F bνµ(x, t) + ηaµ(x, t). (1.14) Once one has gone through the perturbative expansion in terms of Aa(n)µ , the solution in Fourier space reads, Aa(n)µ (k, t) = T ab µν ∫ t 0 ds e−k 2(t−s)f b(n)ν (k, s) + L ab µν ∫ t 0 ds f b(n)ν (k, s) (1.15) where T abµν and L ab µν are the transverse the longitudinal projectors T abµν = ( δµν − kµkν k2 ) δab L ab µν = kµkν k2 δab (1.16) 4 Stochastic quantization for gauge theories satisfying T abµνT bc νρ = T ac µρ L ab µνL bc νρ = L ac µρ T ab µνL bc νρ = 0 (1.17) and f b(n)ν represent the interaction terms. To understand what is going on it suffices to notice that in particular f b(0)ν (k, t) = ηbν(k, t). Projecting the solution (1.15) by the longitudinal projector Labµν the Langevin system presents no damping factor e−k 2t for any component A(n)aµ (k, t) LabµνA (n)a µ (k, t) = L ab µν ∫ t 0 ds f (n)bν (k, s). (1.18) Only the transverse modes are subjected to an attractive force which drives towards the correct solution; longitudinal modes evolves randomly. In particular the field A(0)aµ diverges like a random walk. A second source of divergence is due to the zero mode. Inspecting eq. (1.15) it is obvious that also the zero has no damping factor: e−k 2(t−s) = 1 for k2 = 0. In the case of a finite lattice the integral is replaced by a finite sum, and the relative weight of the zero mode increases with the decreasing of the lattice size. This divergence is then enhanced in small lattices. Both the divergences are associated with the presence of zero eigenvalues in the action, which correspond to a degree of freedom without an attractive force. It is important to stress that this is not a problem in the analytical computation of gauge invariant quantities in Stochastic Quantization, which was the spirit of the original proposal of Parisi and Wu [1]: every non-gauge invariant contribution actually cancels. This can be indeed a problem in a numerical integration of the stochastic process. To overcome the problem of the longitudinal modes we make use of Stochastic Gauge Fixing, first introduced by Zwanziger [6]. The idea of stochastic gauge fixing is in some way similar to gauge fixing in canonical quantization. One introduces a term into the Langevin equation such that the evolution of gauge invariant quantities is not affected. For instance we can write ∂ ∂t Aaµ(x, t) = Dabν F bνµ(x, t)−Dabµ V b[A, t] + ηaµ(x, t) (1.19) where V b[A, t] is a generic non invariant functional. Zwanziger’s choice −Dabµ V b[A, t] = α−1Dabµ ∂νAbν (1.20) introduces a force that keeps limited the norms of the gauge fields. In this 5 Numerical Stochastic Perturbation Theory case the formal solution (1.15) becomes Aa(n)µ (k, t) = T ab µν ∫ t 0 ds e−k 2(t−s)f b(n)ν (k, s) +Labµν ∫ t 0 ds e− k2 a (t−s)f b(n)ν (k, s). (1.21) In this case also the longitudinal component shows a dumping factor. As for the problem of zero mode, we notice that this occurs also in stan- dard finite volume lattice perturbation theory: the generic form of a given Feynman graph is a finite sum which in general entails a singularity for zero momentum. A common prescription to overcome the problem of the zero mode consists in simply dropping the k = 0 contribution [7]. A different strategy can be to impose particular boundary conditions, such as twisted boundary conditions [8]. A simple approach in NSPT is the subtraction of zero modes, that can be obtained imposing all over the trajectory the condition∫ dx Aaµ(x) = 0. (1.22) This means that there is no zero mode contribution to the field, and should be regarded as a prescription on a single mode, becoming irrelevant in the infinite volume limit for a quantity which does not have IR problems. We are the ready for next step: put everything on the lattice. 1.3 Lattice Langevin equation for Wilson ac- tion Wilson formulation of lattice QCD [9] has a gauge action that reads SG = 1 g20 ∑ n,ν 6=µ Tr (1− Uµν(n)) (1.23) where Uµν(n) is the elementary plaquette in the µ − ν plane at the site n, defined as Uµν(n) = Uµ(n)Uν(n+ µˆ)U † µ(n+ νˆ)U † ν(n) (1.24) 6 Lattice Langevin equation for Wilson action in which U(n) = eg0Aµ(n) where Aµ(n) belongs to su(3) algebra. The fermionic contribution to the action is given by SF = ∑ n,m ψ¯(n)Mm,n[U ]ψ(m) = ∑ n (m0 + 4)ψ¯(n)ψ(n)− 1 2 ∑ n,µ [ ψ¯(n + µˆ)(1 + γµ)U † µ(n)ψ(n) +ψ¯(n)(1− γµ)Uµ(n)ψ(n+ µˆ) ] , (1.25) where m0 is the bare quark mass. The action of a gauge transformation on the fields is the following: ψ(n) → G(n)ψ(n) (1.26) ψ¯(n) → ψ¯(n)G(n)−1 (1.27) Uµ(n) → G(n)Uµ(n)G(n+ 1)−1 (1.28) U †µ(n) → G(n+ µ)U †µ(n)G(n)−1, (1.29) where G(n) is an element of SU(N). The Langevin equation (1.3) for the U fields reads ∂ ∂t Uµ(n, t) = (i∇n,µS[U ]− iηn,µ)Uµ(n, t). (1.30) In order to find an explicit form of the derivative ∇n,µS[U ] (notice that this is a Lie derivative), we start integrating out the fermionic degrees of freedom in the action. The result of this integration appears as a determinant and one is then left with a path integral which is integrated only on the gauge degrees of freedom with a new weight given by e−SG[U ]detM [U ] = e−SG[U ]+Tr lnM [U ] = e−Seff [U ]. (1.31) The derivative of the action is then the derivative of the effective action Seff . One simply needs to replace ∇an,µS →∇an,µSeff = ∇an,µSG −∇an,µTr lnM = ∇an,µSG − Tr [ (∇an,µM)M−1 ] (1.32) To carry out the derivative on the gauge action we will need to make use of the definition of the Lie derivative ∇aV f(V ) = lim α→0 1 α [ f(e−iαT a V )− f(V )] (1.33) T a being the generators of the algebra associated with the Lie group. 7 Numerical Stochastic Perturbation Theory Discretizing the Langevin time with step ǫ (and assuming t = τǫ), one can implement a simple Euler scheme for (1.30) in the form U(n, τ + 1) = e−Fµ[U ]U(n, τ) (1.34) where Fµ[U ] = ǫ∇an,µSeff + √ ǫηµ(n). (1.35) Since effective action is a sum of the pure gauge action SG and a term which depends on the fermionic matrix M , we have to compute the sum of the two derivative. The gauge action depends on the link variables Uµ(n) only via the plaque- tte Uµν(n). The derivative then acts on this quantity. Making use of (1.33), we find that ∇n,µTr Uµν(n) = T a∇an,µTr Uµν(n) = iT aTr (T aUµν(n)) . (1.36) Since the plaquette can be decomposed as Uµν = IuI + T bub, (1.37) the element under trace can be computed to be Tr (T aUµν(n)) = Tr ( T auI + T aT bub ) = 1 2 ua (1.38) and one obtains ∇n,µTrUµν(n) = iT a1 2 ua = i 2 ( Uµν(n)− 1 N Tr Uµν(n) ) (1.39) Obtaining the derivative of the fermionic contribution requires a similar procedure. Referring to eq. (1.32) one needs to compute ∇an,µMmβb,lγc = i 2 [ δm,l+µˆδn,l(1 + γµ)βγ ( U †µ(n)T a ) bc −δm,l−µˆδn,l−µˆ(1− γµ)βγ (T aUµ(n))bc ] . (1.40) The difficulty in this case is not the derivative itself, but the inversion of the matrix. We still need to discuss the lattice version of stochastic gauge fixing. Redundant degrees of freedom are associated with large fluctuations. In order to avoid that one has to perform a gauge fixing transformation according to (1.28). We consider the form U ′µ(n) = e iw(n)Uµ(n)e −iw(n+µˆ) (1.41) 8 NSPT for LGT where a choice corresponding to (1.20) is w(n) = −α ∑ µ ∂LµAµ(n) 0 < α < 1 (1.42) where ∂L = is the left derivative defined on the lattice. The strategy is to alternate a Langevin evolution step and a gauge transformation in the form (1.42). This choice results in attracting the fields toward Landau con- dition ∑ µ ∂LµAµ(n) = 0. (1.43) 1.4 NSPT for LGT We now see in some practical detail how Numerical Stochastic Perturbation Theory is implemented. We have already said that we need to numerically implement the substitution of a power expansion in Langevin equation, which in the case of Lattice Gauge Theories we take in the form of (1.30). The starting point is the expansion of link variables in a power series in the coupling U = I + β− 1 2U (1) + β−1U (2) + . . . (1.44) This is not the usual continuum perturbation theory expansion in the algebra. Though, since A = logU the two expansions are related: A = logU = log ( I + β− 1 2U (1) + β−1U (2) + . . . ) = β− 1 2U (1) + β−1 ( U (2) − 1 2 U (1) 2 ) + . . . (1.45) U = eA = e(β − 1 2A(1)+β−1A(2)+...) = I + β− 1 2A(1) + β−1 ( A(2) + 1 2 A(1) 2 ) + . . . (1.46) Both expansions can be used. While the expansion in the algebra makes easier to keep unitarity, the expansion in the group (1.44) is more natural in the evolution step (1.30). We now need to implement the NSPT version of the evolution step (1.34). Since Fµ is the sum of two pieces (gluonic and fermionic), a first order con- venient numerical scheme is the one in which we write the exponential of the sum as the product of two exponentials. In practice, we split the evolution 9 Numerical Stochastic Perturbation Theory step in a pure gauge step plus a fermionic step. Let us now discuss the two contributions which are in place. The pure gauge contribution is given by (1.39) Fµ[U ] = ǫβ 4N ∑ n,ν [ Uµν(n)− 1 N Tr Uµν(n) ] + √ ǫηµ(n) (1.47) This can in turn be expanded in a power series Fµ[U ] = ∑ k β− k 2F (k)µ [U ] (1.48) once the time step is rescaled according to ǫβ = τ . The significance of this rescaling is the following: higher values of β requires smaller steps ǫ in order to keep discretization errors under control. It is important to notice that, since the plaquette contribution Uµν(n) is a product of link variables, which are in turn expanded in series, the product must be taken order-by-order, defined as (A ∗B)(n) = ∑ i≤n A(i)B(n−i). (1.49) As for fermions, we saw that their contribution in (1.32) is non trivial because of the computation of the inverse operator. It turns out that such a computation in NSPT is quite trivial. The fermionic operator M can be expanded in power series M =M (0) + ∑ k>0 β− k 2M (k), (1.50) and the inverse operator is obtained inverting the series: M−1 = (M−1)(0) + ∑ k>0 β− k 2 (M−1)(k). (1.51) The notation enlightens the fact that the zeroth-order of the inverse is the inverse of the zeroth-order. At higher orders a recursive equation holds (M−1)(1) = −M (0)−1M (1)M (0)−1 (M−1)(2) = −M (0)−1M (2)M (0)−1 −M (0)−1M (1)(M−1)(1) . . . (M−1)(n) = −(M−1)(0) n−1∑ j=0 M (n−j)(M−1)(j) (1.52) 10 NSPT for LGT This shows that there is no need of any explicit inversion but for zero order: what is left comes out of a recursion, in which M (n) are in place. Moreover, the inverse of the zeroth order has, in the case of Wilson fermions, an easy analytical form: (M (0))−1 = m0 + 2 ∑ µ sin 2 ( kµ 2 ) − i∑µ γµ sin2 (kµ)[ m0 + 2 ∑ µ sin 2 ( kµ 2 )]2 + ∑ µ sin 2 (kµ) (1.53) which is diagonal both in momentum and in color space. As we saw be- fore, the evaluation of the derivative is an easy task; one is only left with the evaluation of the trace. The idea [2] is to introduce another Gaussian (nonperturbative) source ξ, 〈ξiξj〉ξ = δij, (1.54) which allows to write the trace which appears in (1.32) as Tr [ (∇n,µM)M−1 ] = Re [ ξ†k(∇n,µM)kl(M−1)lnξn ] . (1.55) Introducing a vector ψ, solution of the equation Mklψl = ξk (1.56) which turns out to be perturbative, the perturbative fermionic contribution has a more compact form which reads at k-th order k∑ j=0 [ ξ†k(∇an,µM (j))klψ(k−j)l ] . (1.57) The recursive form of the inversion (1.52) translates into the recursive con- struction of the field ψ: ψ(0) = M (0) −1 ξ ψ(1) = −M (0)−1M (1)ψ(0) ψ(2) = −M (0)−1 [M (2)ψ(0) +M (1)ψ(1)] . . . ψ(n) = −M (0)−1 n−1∑ j=0 [ M (n−j)ψ(j) ] (1.58) This is the closed formula which is actually numerically implemented. In order to obtain the fermionic contribution one needs to construct the pertur- bative ψ field, apply the matrix (∇n,µM) on it and finally contract with ξ, 11 Numerical Stochastic Perturbation Theory provided that every operation is intended as order by order. Moreover, this inversion scheme is not limited to field evolution. Every times a fermionic observable is in place, such inversion in required. From eq. (1.58) we see that the inverse of zero order has to be applied several times in the recursion. In order to do that one goes back and forth from Fourier space: in one space the “tree level propagator” is trivial, in the other the direct application of M is easy. The only problem associated to the expression (1.53) is its computation at k = 0 in case of zero bare mass. Various infrared regularizations are available: leave a small mass, introduce anti-periodic boundary conditions for the fermionic fields in time direction or constrain the zero mode degree of freedom to zero. The last one is the choice we adopted. Let us now go back to the splitting of the Langevin dynamics. As we said, since we are using an Euler integration scheme, we can proceed in the following way: • evolution of the pure gauge contribution according to a dynamics dic- tated by (1.47); • iterative construction of ψ by means of the application of multiple FFT; • evolution by the fermionic contribution F (ferm)µ = ǫRe [ ξ†k(∇n,µM)klψl ] ; • perform a gauge fixing step according (1.41). • subtraction order-by-order of the zero momenta of the Aµ field. This is simply obtained by A(i)µ (k = 0) = ∑ n A(i)µ (n). (1.59) Since the equivalence between canonical and stochastic quantization is true in continuous time, the discretization errors due to the finite size of ǫ must be extrapolated to zero. Such extrapolation is performed after having simulated at different integration time steps. 12 CHAPTER 2 PRlgt: an environment to implement NSPT Numerical Stochastic Perturbation Theory is something which lies in between a traditional (analytical) perturbation theory computation and a MonteCarlo simulation. It is a perturbation theory computation, which means that we compute something that is equivalent to Feynman diagrams computation. But this comes out of the integration of a stochastic process. We are dealing with numerical computations like in nonperturbative MonteCarlo, like for example Hybrid MonteCarlo. Indeed, we need data analysis to take into ac- count statistics, we have to deal with thermalization and correlation between data just like in a typical MonteCarlo. Being NSPT “Numerical”, an appropriate package is required. We devel- oped a library called PRlgt, which is not only a perturbative library, but allows traditional computation too. The main ideas of the implementation are the following: • a “multilayer” C++ implementation which could account for the per- turbative structure of the technique; • a modular implementation which allows development without any mod- ification at the rest of the code; • data are organized in order to optimize parallelization. 13 An environment for NSPT 2.1 Hide perturbation theory The main problem we have to deal with are perturbative computations. This means that the fields are power series in the coupling, and any operation must be performed taking into account the correct power of the coupling itself. One possibility is of course to write down the expansion of the action in a “verbose” form; however this wold require a dramatic grows of the number of terms. Let’s consider for example two perturbative matrices A and B: A = A(0) + gA(1) + g2A(2) + . . . B = B(0) + gB(1) + g2B(2) + . . . (2.1) then the simplest (non trivial) operation one can think of is the product, which reads A ∗B = (A(0) ∗B(0))+ g (A(1) ∗B(0) + A(0) ∗B(1))+ g2 ( A(2) ∗B(0) + A(1) ∗B(1) + A(0) ∗B(2))+ . . . (2.2) where the product A(i) ∗B(j) (2.3) is the traditional matrix product. If we keep in mind that the basic element in the computation of the gauge action, the plaquette, is the product of four matrices, we easily understand that a straightforward, “verbose” expansion of action and observables would become impossible after just a few orders. Moreover, one should compute those expansions for every quantity. A useful reorganization of data will help us in this task. If we consider the expansions (2.1) as vectors of matrices, it is easy to define perturbative operations. For example the multiplication (2.2) becomes A ∗B → for i ∈ 0 . . . k (2.4) for j ∈ 0 . . . k− i (A ∗ B)[i] = (A ∗ B)[i] + A[i− j] ∗ B[j] end end The C++ implementation of this strategy is quite natural: one defines perturbative structures (typically, classes) on which for example multiplica- tions (or any other operation) are defined by “overloading” standard opera- tors. In our example (2.4), we will overload the definition of the operator ∗ between perturbative matrices. 14 Multilayer code This allows to write a code which is very similar to our “natural” way of writing. Any operation can be implemented by means of the same math- ematical symbol one would have on the paper. All the knowledge of the details of the computation is hidden inside the definition of the objects and operators. 2.2 Multilayer code In order to have a more flexible code which would allow extensions without effort, the library is “multi layer”. There are three main layers, corresponding to the physical (or mathematical) object one is dealing with: • basic background mathematics (complex numbers, matrices, vectors); • low level physical entities (gluons and fermions); • fields. This multilayer system enables one to forget about the underlying or over- lying structures when one needs to improve one particular layer. It is then possible to optimize the lowest level (in order to obtain a faster computation or to specialize it for a given architecture) without any further modification in the code. At the same time it is possible to implement new features; for example, we can have SU(2) or SU(3) fields. The innest building block of the whole data structure are the complex numbers, implemented by the class Cplx which is made of two members, real and imaginary part respectively: class Cplx { public: double re,im; Cplx (double r = 0., double i=0.) { re = r; im = i;} Cplx operator+(const Cplx& b) const { return Cplx(re+b.re, im+b.im); } . . . 15 An environment for NSPT }; By means of Cplx overloading of all the useful operators, living with the complex nature of all the computations for all (higher level) structures gives no problem. Since the mathematical objects required in any LGT computation are matrices and vectors, the next step in the library is to have them imple- mented: class SU3 { public: Cplx whr[NC*NC]; SU3 (){}; SU3 operator+(const SU3& A) const { SU3 res; for( i = 0; i < NC*NC; i++) { res.whr[i] = whr[i] + A.whr[i]; } return res; } . . . }; where NC is the number of colors. The reader will notice that as it often happens to computer programmers, names can be misleading: SU3 was the name chosen when we only had the number of colours NC = 3, but the same name remained after the generalization to SU(N) ... NC is a parameter that can be chosen at compilation time, enabling simulations for different gauge groups. As already pointed out, SU(3) and SU(2) simulations are currently being performed. Similarly to SU(3) matrices the low level library contains colour vectors, which are the main ingredient for spinors: 16 Multilayer code class CVector { public: Cplx whr[NC]; CVector (){}; CVector operator+(const CVector& V) const { CVector res; for( i = 0; i < NC; i++) { res.whr[i] = whr[i] + A.whr[i]; } return res; } . . . }; The objects implemented in the low level library are quite general: they are not strictly related to perturbative or nonperturbative computations. Only at this point one has to think about this problem. The solution, as pointed out before, is indeed quite easy. If we carry on this strategy, the most obvious thing is to define the perturbative objects equivalent to matrices and vectors: class ptSU3{ private: public: Cplx flag; SU3 ptU[allocORD]; ptSU3(Cplx i = 1.0) { flag = i; } ptSU3 operator*(const ptSU3& A) const { ptSU3 B; 17 An environment for NSPT for(int i = 0; i < (PTORD-1); i++) { for(int j = 0; j < (PTORD-1-i); j++) { B.ptU[i+j+1] += ptU[j]*A.ptU[i]; } } if( (flag == 0) && (A.flag == 0) ) { B.flag = 0; return B; } else { B.flag = flag * A.flag; for(int i = 0; i < PTORD; i++) { B.ptU[i] += A.flag*ptU[i] + flag*A.ptU[i]; } return B; } } . . . }; Notice that both in the attributes (look at the flag) and e.g. in complex multiplication, we can inspect some non trivial consequence of NSPT. This is nothing but the implementation of (2.4), in which a trivial vacuum back- ground has been chosen. In this implementation, the class ptSU3 is not only made of perturbative matrices: the first element (zero order) is a single com- plex number. This account for the possible expansion of an element of the SU(N) gauge group around the ZN center of the group itself. In the case of SU(3), for example, the possible flags are 1, e−ipi/3, e−2ipi/3. Moreover some- times expansion of an element of the algebra are required. This special case is taken into account to avoid unnecessary multiplications and consequent waste of computing power: this is the reason for the if condition in the multiplication showed above. Restricting to trivial vacuum is of course not compelling (think about the Schrödinger Functional scheme): having this specific implementation for the (most common) case of standard PT is a matter of efficiency. 18 Multilayer code class ptCVector { public: CVector ptCV[allocORD + 1]; ptCVector operator*(ptSU3& A) const { ptCVector cv; for(int i = 0; i < (PTORD-1); i++) { for(int j = 0; j < (PTORD-1-i); j++) { cv.ptCV[i+j+1] += ptCV[j]*A.ptU[i]; } } if(A.flag != 0) { for(int i = 0; i < PTORD; i++) { cv.ptCV[i] += A.flag*ptCV[i]; } } return cv; } . . . }; The case of perturbative color vector is similar. However in this case the first element in the expansion is not trivial, since there are no fermions at tree level, they enter only at one loop and the first element of the fermionic expansion must be nontrivial. The reported example of multiplication be- tween a perturbative matrix and perturbative color vector is the building block of any simulation involving fermions. All the mathematical entities involved have been created: it is possible to switch to physics. In QCD gluons and quarks are involved, and the next step is to create these objects both in the nonperturbative and perturbative case. 19 An environment for NSPT Gluons are made of SU(3) matrices, in a number equal to the number of space-time dimensions, typically four. This information is encoded in the parameter dim, which can be eventually chosen to be different from four as well. Once again, simply changing a parameter makes life easy, from the point of view of both code implementation and reading ease. class Gluon { public: SU3 U[dim]; Gluon(){}; Gluon dag(const Gluon &A) { Gluon res; for(int i = 0; i < dim; i++) { res.U[i] = dag(A.U[i]); return res; } . . . }; The only operator we report is the dagger operator. In the perturbative case the structure is similar class ptGluon { public: ptSU3 U[dim]; ptGluon() {}; ptGluon dag(const ptGluon &A) { ptGluon res; for(int i = 0; i < dim; i++) { res.U[i] = dag(A.U[i]); 20 Multilayer code } return res; } . . . }; The last “single” item we have to set up is the “quark”, or more generally the fermion. In agreement with what done so far we want perturbative and non perturbative to be similar. class SpinColor { public: CVector psi[dim]; SpinColor(){}; SpinColor& operator*=(const Cplx& z) { for (int mu = 0; mu < dim; mu++) { psi[mu] *= z; } return *this; } . . . }; class ptSpinColor { public: ptCVector psi[dim]; ptSpinColor(){}; ptSpinColor& operator*=(const Cplx& z) { for (int mu = 0; mu < dim; mu++) 21 An environment for NSPT { psi[mu] *= z; } return *this; } . . . }; Let us stress once again one of the big advantages in all this procedure: we can write operations between different entities in a “human” form, just like what one would write on a piece of paper. A multiplication between two object A and B will be written A ∗B, regardless of the nature of A and B. 2.3 Fields and parallelization The successive layer of the library is strictly related to data arrangement. There are three main classes responsible for the final physical entities one wants to deal with, which are fields. Fields are nothing but a collection of a certain number of gluons or fermions, the number depending on the size of the lattice. There are two main different kinds of perturbative and nonperturbative fields: ptGluon_fld, ptSpinColor_fld, Gluon_fld and SpinColor_fld. On top of these we need two extra fields to account for gauge transformations: SU3_fld and ptSU3_fld. In order to show again the symmetry between perturbative and nonperturbative we report as an example the two gluonic fields. class Gluon_fld{ public: latt *Z; Gluon *W; Gluon_fld(latt *z){ Z = z; W = new Gluon[z->Size]; } . . . 22 Fields and parallelization }; class ptGluon_fld{ public: latt *Z; ptGluon *W; ptGluon_fld(latt *z) { Z = z; W = new ptGluon[z->Size]; } . . . }; The object latt which appears in both definitions (as well as in all the other fields) is the data structure which is responsible for data arrangement and neighbourhood relations. Without entering into details, we will see it at work later on in a basic sample code. We stress that it holds all the infor- mations on “who lives where” in memory and provides the environment for the basic “translation” operators which we will encounter. The importance of a clever data organization has mostly to do with the parallelization issue, which will now discuss in brief. A typical parallel lattice simulation is accomplished by spreading the whole lattice configuration among computing units. In the simplest update step one needs informations about next neighbours, which in general (and in particular in the case of a boundary variable) can be in the memory of a different computing unit. This requires communications among nodes. As an example, we discuss the case of a two dimensional parallelization, i.e. the case in which the lattice is spread in two dimensions. Figure 2.1 shows the scheme for such a parallelization It is possible to recognise three main areas: • the bulk, whose data do not require to access anything outside the local memory; • inner boundaries where part of the information required has to come from one (or more) neighbour. On the other hand, the neighbours require these data in order to make their computations; 23 An environment for NSPT Forward edge Backward edge Backward edge Forward edge R ig h t ed g e L ef t ed g e L ef t ed g e R ig h t ed g e Bulk Forward edge Backward edge Backward edge Forward edge R ig h t ed g e L ef t ed g e L ef t ed g e R ig h t ed g e Bulk Forward edge Backward edge Backward edge Forward edge R ig h t ed g e L ef t ed g e L ef t ed g e R ig h t ed g e Bulk Forward edge Backward edge Backward edge Forward edge R ig h t ed g e L ef t ed g e L ef t ed g e R ig h t ed g e Bulk Figure 2.1: Example of parallelization in two dimensions. • outer boundaries which are data coming from the neighbours. These data are not updated by the node, but only accessed. The arrows in figure 2.1 represent internode transfers required before any upgrade step. Figure 2.2 shows how data are usually allocated in the computer, following a row (or column) major ordering. Overlapping this with figure 2.1 suggest that part of the data to be sent (and received) are not contiguous. This would require either many send/receive operations or a “collect” operation which would store data to be sent in a temporary location, which results in an overhead. A useful organization of memory can be the following: one can place different sublattice portions into contiguous tranches of memory (figure 2.3). Since each node has four neighbours, it has to divide its sublattice into many portions, namely the bulk, the inner edges and the outer edges. Particular care has to be taken to the overlapping parts on the corners (darker red in figure 2.1) where two edges meet, which have to be sent to more than one neighbour node. Since we want to perform only one transfer, the regions 24 Fields and parallelization 0 1 2 3 Nx Nx + 1 . . . 2Nx . . . . . . x y Figure 2.2: Non optimized data organization. which have to be sent to more neighbours must be duplicated. Bulk Inner boundaries Outer boundaries Figure 2.3: Data reorganization allows the reduction of the number of com- munications. In order to fulfil such a data arrangement, a “lookup table” contains all the informations concerning the positions. This is exactly the object latt whose reference is contained in the field classes. The core of the latt is the lookup table which contains both the translation between “logic” position in the lattice and “physical” position in memory and the the memory position of all the next neighbours, accessible without requiring any computation. It is important to underline that all the fields refers to one and the same lattice structure. int main(int argc, char** argv){ . . . latt LL(act_pars.sz); ptGluon_fld Umu(&LL); ptSpinColor_fld Pmu(&LL); SpinColor_fld Xi(&LL); 25 An environment for NSPT . . . }; Once the fields have been created, one can completely forget about the memory organization. The access to a specific lattice site or link, the access to a next neighbour or to a site/link along a specific path is performed via specific methods of the field classes: • access to a specific site: ψ(x) Pmu.get(x) • access to a specific link: Uµ(x) Umu.get(x,mu) • access to a link at the thuple: ψ(~n) Pmu.get(n) • access to a next neighbour site: ψ(x+ µˆ) Pmu.near(x,1,mu) • access to a next-next neighbour link: Uν(x− µˆ+ νˆ) Umu.get(x,-1,mu,1,nu,nu) • access to a site along a given (two dimensional) path: ψ(x+ 3µˆ− 5νˆ) Pmu.get(x,3,mu,-5,nu) We thus have a general library for both perturbative and nonperturba- tive QCD, which enables us to write a variety of specific programs with ease (different people have actually used it). As said, we can stay as close as possible to “human writing”. Since we had to deal also with data obtained on the APE machine in Rome, we also developed all the interfaces which are needed to load configurations generated on the other architecture. All the results in the following chapters have been obtained with this code, whose implementation was one of the objective of the doctoral research. 26 CHAPTER 3 Renormalization coefficients Renormalization of operators for Lattice QCD has two main purposes. Like in any regularization of a field theory, a renormalization condition is required to get physical results which are both finite and independent from the regulator. As for any regulator, infinities are still present if the regulator is removed in a naive way, but in the renormalization procedure they are reabsorbed in the renormalization coefficients (RC’s): OR = Z−1Obare (3.1) On top of this, one should keep in mind that finite renormalizations are needed to connect different renormalization schemes, and in the end lattice results are always to be matched to continuum QCD (whose standard schemes are actually perturbative) ORlat = Zcont Zlat ORcont (3.2) For a long time perturbation theory was the only available tool for com- puting RC’s on the lattice. As it is well known (see for example [10]), per- turbative computations on the lattice require a big effort, because of the appearance of many irrelevant contributions. This prevents computations at very high orders with standard (diagrammatic) techniques. The introduc- tion [11] [12] of non-perturbative renormalization drastically changed this scenario: by now non-perturbative computations of RC’s are in general pre- ferred. In view of this, it is worth to spend a few words on the relevance of a perturbative computation like the one we are going to report on. First of all, 27 Renormalization coefficients it is important to underline that there is no theoretical obstacle to the compu- tation of both finite and logarithmically divergent RC’s, like quark bilinears (or their ratios). Moreover, renormalization is proved to be multiplicative only in perturbation theory, even if there is no serious concern against the generally accepted expectation that this holds also in the non-perturbative case. Having said this, the real issue is the control of systematics. The choice of any tool which is at work in a lattice computation has always had to face a precision issue. But now the lattice community regards the current one as the era of high precision measurements: the goal is to achieve percent accuracy in the computation of many observables. It thus goes without saying that the assessment of systematic errors is a crucial issue. A typical example is the computation of quark masses: the values changed quite a lot over the years, and the discrepancies between different results have often been ascribed to the issue of perturbative vs non-perturbative techniques in the computation of RC’s. This is a snapshot for s-quark mass results as in 2007 [13]a: MILK(2004) mMSs (2GeV ) = 76(3)(7)MeV (PT 1 loop) MILK(2006) mMSs (2GeV ) = 90(5)(4)MeV (PT 2 loop) ETMC(2007) mMSs (2GeV ) = 105(3)(9)MeV (non-PT) Table 3.1: Different determinations of the s-quark mass [13]. One could argue that things changed a little bit since then, but the overall picture still holds: it looks like determinations coming from different methods for the computation of the relevant RC’s can not find an agreement within the errors they quote. The very point is that as long as a comparison is made taking into account one loop (and even two), it is virtually impossible to assess the systematics of both results. Non-perturbative computations are in fact supposed to get rid of trunca- tion errors, which are indeed an obvious matter of concern with perturbation theory. Still, the general landscape of systematic errors is wide for both per- turbative and non-perturbative techniques. The shopping list for the control of systematics includes at least the following: 1. PT definitely has to face truncation errors. 2. Most of the schemes are defined in the massless limit, and thus one most often has to deal with chiral extrapolations. aThe reader will notice that this is the start date of this research project. 28 3. In any case one needs to face extrapolations toward the continuum limit. 4. If the renormalization scheme is not a finite volume one, one has to keep finite size effects under control. It is important to stress that while the first item of this list only pertains to a perturbative computation, all the other stand both for perturbative and non-perturbative computations. Our goal is to compute three loops quark bilinears RC’s for the action defined by Symanzik tree level improved gauge action and nf = 2 Wil- son fermions. The former gauge action attracted much attention in recent times because it has been adopted by the European Twisted Mass Collabo- ration [14]. As a matter of fact, since RC’s are defined in the massless limit, our results for Wilson fermions can be directly compared to the Twisted- Mass fermionic action. In particular, results from non-perturbative ETMC determination of RC’s are available in the literature [14], with a choice for renormalization scheme which the same we adopt here, i.e. RI’-MOM. The claim is that all the systematics can be taken under fairly good control. 1. In subsection 3.1 we show how the knowledge of anomalous dimensions in RI’-MOM enables us to reach three loop order. With this respect, the only missing item is the two loop matching between lattice and con- tinuum couplings for tree level Symanzik improved gauge action: this is filled in subsection 3.2 We thus set a new standard for perturbative results for RC’s: what used to be one loop is now three loop order. It goes without saying that the usage of NSPT is the technical break- through for this achievement, as it has already been showed to be the case in other computations [15,16,17,18]. In particular, three (and even four) loops order has already been achieved in a similar computation of quark bilinears for the Wilson gauge action. 2. In subsection 3.3 we show how we are able to stay in the chiral limit. 3. In subsection 3.4 we show how we take the continuum limit. 4. Finally, in subsection 3.5 we discuss the control on finite size effects, which is an important issue, since RI’-MOM is defined in infinite vol- ume. We stress that this is a valuable add on with respect to previous 29 Renormalization coefficients computations of quark bilinears for Wilson gauge actionb. The relevant question we finally address is to what extent our results can contribute to better assess the overall systematic error which plagues the RC’s at hand: to this issue we devote subsection 3.6.5 3.1 RI’-MOM scheme The first step one has to go through in the computations of RC’s is to choose the scheme among different possibilities. This depends on the nature of the computation one is interested in, and each scheme presents advantages as well as disadvantages. One of the most common, in particular for perturbation theory on the continuum, is the min- imal subtraction (MS) scheme (the bar means that the actual scheme is a slightly modification, in particular in the choice of the scale, with respect to the traditional MS scheme), which has the advantage to keep coefficients of the perturbative expansion fairly small, and this leads of course to a better behaviour of the series. The drawback is that the scheme depends on the regularization procedure [19]. Another common scheme (in particular since the introduction of non- perturbative renormalization) is the RI’-MOMc, which presents expansion coefficients bigger than the previous, but, as the name RI suggests, is Regula- tor Independent, which means that it can be implemented with any regulator (for example, dimensional or lattice). We will adhere to RI’-MOM scheme for the following reason: in order to perform the computation of RC’s at three loop, anomalous dimensions are required at the same order. Now, anomalous dimensions are dependent on the scheme, but not on the regulator. By choosing a scheme in which they are known at the required order [20] one can get them from computations performed with a different regulator. In general, it would be very hard to compute logarithmic contributions in NSPT (this is a matter of numerical accuracy). One could say that within standard, diagrammatic techniques, it is easy to compute log’s, while it can be very hard to compute finite con- stants. In NSPT, it’s just the other way around. bThe same technology at work here will soon fill this gap in the case of Wilson gauge action cThe prime denotes a renormalization condition for the quark field which is slightly different from the original one. 30 RI’-MOM scheme In the RI’-MOM scheme the quark wave function renormalization condi- tion is given by Zq = −i 1 12 Tr ( /pS−1(p) ) p2 ∣∣∣∣∣ p2=µ2 , (3.3) where S−1(p) is the two point massless correlation function and p is the external momenta. For quark bilinears of the form OΓ = ψ¯Γψ, where Γ = 1, γµ, γ5, γµγ5, σµν stands for scalar, vector, pseudoscalar, axial and tensor currents respectively, the renormalization condition in the RI’-MOM scheme is given by ZOΓ(µa, g(a))Z −1 q (µa, g(a))OΓ(pa) ∣∣ p2=µ2 = 1. (3.4) More detail on the definition of the scheme can be found in [21]. Being these quantities gauge dependent, a choice of gauge condition has to be made. Using Landau gauge, easy to fix on the lattice [22], one gets an extra bonus: the anomalous dimension for quark field is zero at one loop. 3.1.1 Method The basic quantity we need to compute are the quark bilinears at fixed (off- shell) external momentum p:∫ dx〈p|ψ¯(x)Γψ(x)|p〉 = GΓ(pa). (3.5) Given the quark propagator S(pa), one can obtain the amputated functions GΓ(pa)→ ΓΓ(pa) = S−1(pa)GΓ(pa)S−1(pa), (3.6) to be eventually projected on the tree level structure by a suitable operator POΓ to obtain the observable OΓ(pa) = Tr ( PˆOΓΓΓ(pa) ) . (3.7) In order to develop a mass independent renormalization scheme, one im- poses Renormalization conditions (3.3) and (3.4) on a massless quark. This requires, in perturbation theory, to introduce the relevant counterterms for the so-called critical mass. One and two loops results for this quantity are known from the literature [23], while the third loop is a byproduct of the current computations. The situation in the non-perturbative framework is more cumbersome: it is not possible to have a massless quark simulation, and the critical mass is 31 Renormalization coefficients the prototype of a non-perturbative computation of an additive renormal- ization constant. Being a power divergent quantity, it cannot be computed in perturbation theory in the continuum limit, and the massless limit is ob- tained by means of a chiral extrapolation. This is one of the possible sources of errors in the determination on non-perturbative RC’s. So far we used the notation pa and µa, which accounts for the presence of the lattice regulator. Starting from here we introduce the definition pˆ = pa and µˆ = µa. Similarly, one has to remember that any observable on the lattice should be multiplied for the correct power of the spacing a in order to get the correct dimension. In order to account for that, we use the notation “hat” also for the other quantities, avoiding to put the correct number of a. In the continuum the two point function is given by Γ2(p 2) = S−1(p2), (3.8) while on the lattice it contains also irrelevant terms, and can be decomposed according to Γˆ2(pˆ 2, mˆcr, β −1) = Sˆ−1(pˆ, mˆcr, β −1) = i/ˆp+ mˆW (pˆ)− Σˆ(pˆ, mˆcr, β−1) (3.9) where the dependence on the coupling β and on the critical mass mˆcr is made explicit. Here mˆW (pˆ) = O(p2) is the irrelevant mass term generated by the Wilson prescription and Σˆ(pˆ, mˆcr, β−1) = O(β−1) is the quark self- energy which we will discuss in detail later in this section. The term mˆW (pˆ) is present also in the case of massless fermions, due to the chiral symmetry breaking in the Wilson regularization. This is the content of the Nielsen- Ninomiya no-go theorem: it is not possible to have a lattice fermionic action which is local, chirally invariant and free from unphysical degrees of freedom. In the case of Wilson fermions one relaxes the condition of chiral symmetry, and is then left with an additive mass counterterm (in other cases, one gives up with other conditions: for example with staggered fermions one accept to have unphysical degrees of freedom - the taste of the fermions - in order to keep chiral invariance). 3.1.2 Deep into RC’s Since the RC’s must reabsorb the divergences coming from the propagator or the currents, they must be divergent too. Moreover, since the redefinition of the fields from bare to renormalized must not change the form of the La- grangian, the multiplicative RC’s must be dimensionless. 32 RI’-MOM scheme The continuum limit form of Zq with a lattice cutoff reads Zq(µa) = 1 + ∞∑ n=1 αndn dn = n∑ i=0 d(i)n L i (3.10) where α is a renormalized coupling and L = log(µa). The dn coefficients contain the dependence on the anomalous dimension γq, defined by γq = −∂logZ 1/2 ∂L = ∞∑ n=1 γ(n)q α n (3.11) Being the RI’-MOM scheme regulator independent, one can get the coeffi- cients of the expansion from continuum computations [21] [20]. In order to plug them in our lattice computation at the order we are interested in, we will see that one needs to know the matching between the couplings in the two schemes up to two loop. This matching was known only up to one loop, and so we needed to compute the next loop (section 3.2). Plugging (3.10) in the definition of the anomalous dimension and keeping in mind that when the derivative with respect to L acts on α one gets the β-function β(α) = ∂α ∂L2 (3.12) = −b0α2 − b1α3 − b2α4 + . . . . (3.13) one obtains an expression which has to equal the known expansion of γq − ∂logZ 1/2 ∂L = −1 2 d (1) 1 α + [ b0d (0) 1 + 1 2 (d (0) 1 d (1) 1 − d(1)2 ) (3.14) + (b0d (1) 1 + 2d (1) 1 2 − d(2)2 )L ] α2 + . . . One can solve for the relevant coefficients by canceling all the log’s in the final expression for γq d (1) 1 = −2γ(1)q (3.15) d (1) 2 = 2(b0d (0) 1 − d(0)1 γ(1)q − γ(2)q ) d (2) 2 = −2(b0γ(1)q − γ(1)q 2 ) . . . 33 Renormalization coefficients Notice that we have till now worked with expansions in a continuum cou- pling, in which the computation of anomalous dimension has been performed. By the way, the practical choice made in the literature is that this coupling is αMS. Since we are interested in obtaining our result in terms of the bare cou- pling α0, we need to substitute the expansion of the renormalized α according to α = α0 + c1α 2 0 + c2α 3 0 + . . . (3.16) which we will discuss in detail in section 3.2. For the moment we are only interested in the log’s content of the ci factors, which have the form c1 = −2b0L+K1 (3.17) c2 = 4b0L 2 − (2b1 + 4b0K1)L+K2, (3.18) In the previous formulaK1 andK2 have known forms in terms of Λ-parameters and coefficients of β-function (even if the numerical value of K2 is to be de- termined). Since we are working in Landau gauge one loop anomalous dimension vanishes (γ(1)q = 0), so that Zq in terms of α0 reads Zq = 1 + Z (1) q α0 + ( Z(2)q − 2γ(2)q L ) α20 (3.19) + [ Z(3)q − 2 ( γ(3)q + 2γ (2) q K1 + Z (1) q γ (2) q ) L+ 4b0γ (2) q L 2 ] α30 +O(α40) The complete expressions for the finite parts read Z(1)q = d (0) 1 (3.20) Z(2)q = d (0) 2 +K1Z (1) q (3.21) Z(3)q = d (0) 3 − 2K21Z(1)q +K2Z(1)q + 2K1Z(2)q (3.22) In order to compute RC’s for composite operators (in our case bilinears), we start from the relation between bare and renormalized two point function: Γ0 = Z −1 O (µa)ΓR (3.23) with ZO(µa) = 1 + ∞∑ n=1 αnen en = n∑ i=0 e(i)n L i (3.24) 34 Matching αTLS to αMS Manipulating once again the definition of the relevant anomalous dimension, one obtains ZO = 1 + ( Z (1) O − 2γ(1)O L ) α0 (3.25) + [ Z (2) O − 2 ( γ (2) O + γ (1) O K1 + Z (1) O γ (1) O ) L+ 2γ (1) O ( γ (1) O + b0 ) L2 ] α20 + [ Z (3) O − 2 ( γ (3) O + 2γ (2) O K1 + γ (1) O K2 + γ (2) O Z (1) O − 2γ(1)O K1Z(1)O − γ(1)O (Z(2)O −K1Z(1)O ) ) L+ 2 ( b1γ (1) O + 2b0γ (2) O + 2γ (1) O γ (2) O + 2b0γ (1) O K1 + 2γ (1) O 2 K1 + b0γ (1) O Z (1) O + γ (1) O 2 Z (1) O ) L2 − 4 3 ( b20γ (1) O + 3b0γ (1) O 2 + γ (1) O 3 ) L3 ] α3 +O(α40) where Z (1) O = e (0) 1 (3.26) Z (2) O = e (0) 2 +K1Z (1) O (3.27) Z (2) O = e (0) 3 − 2K21Z(1)O +K2Z(1)O + 2K1Z(2)O (3.28) Again, the various γ(i)O are known from literature [21] as expansions in αMS. 3.2 Matching αTLS to αMS As we have already said, the only missing item in our previous formulas was K2 in (3.18). The general form of the matching between two schemes, in which the coupling constants are respectively α and α′, is well knownd α(sµ) = α′(µ) + c1(s)α ′(µ)2 + c2(s)α ′(µ)3 + . . . , (3.29) where s is a coefficient that accounts for the choice of a different scale, and the coefficients c1(s) and c2(s) are given by c1(s) = 2b0 log Λ Λ′ − 2b0 log s (3.30) c2(s) = c1(s) 2 − 2b1 log s+ 2b1 log Λ Λ′ + b2 − b′2 b0 . (3.31) dWe stop at the order we are interested in. 35 Renormalization coefficients b0, b1 and b2 are coefficients of the β function, and Λ is the parameter asso- ciated to the regularization. Comparing to (3.17) and (3.18) one can recognise that K1 = 2b0 log Λ Λ′ (3.32) K2 = 4b 2 0 ( log Λ Λ′ )2 + 2b1 log Λ Λ′ + b2 − b′2 b0 . (3.33) It is worth to remember that b0 and b1 are universal, while Λ and b2 depend on the scheme. In order to compute the finite coefficients of the RC’s we need the two loop matching of αTLS (starting from here we will refer to αTLS as the coupling defined by the tree level Symanzik improved gauge/Wilson nf = 2 fermions action) to αMS which entails the knowledge of b TLS 2 , since bMS2 is known. For one loop matching see for example [24]. 3.2.1 Method One could of course envisage a direct matching of MS to the tree level Symanzik scheme. Another strategy, the one we will go through, is to match- ing first to an intermediate scheme: with this approach we need no compu- tation of logarithm, which would require a terrific precision in NSPT. Following [25] [26] we exploit the knowledge of the static quark poten- tial, which defines a renormalized coupling called αV (R), which has been computed in the MS scheme [27] at two loops. The static quark potential is the quantity which describes the interaction energy of a infinitely heavy qq¯ pair. This is the non abelian equivalent of the Coulomb potential; it is a quantity which is expected to play a key role in the understanding of the theory, in particular the confinement mechanism. In its full (non-perturbative) form, it is in first approximation just the sum of a long distance term (string tension), which is responsible for confinement, and a R−1 contribution, whose interpretation is different in different IR/UV regimes V (R) = C R + σR. (3.34) In the non-perturbative region the first term is the leading order −π/12 correction to string energy (see for example [28]), while in the perturbative (UV) regime is just the Coulomb potential. Because of the large time extent (T ) behaviour of the Wilson loop for a heavy pair QQ¯ W (R, T ) ≡ 〈WC [A]〉 −−−→ T→∞ F (R)e−V (R)T (3.35) 36 Matching αTLS to αMS the static quark potential as a function of the distance R between the quarks can be defined via the vacuum expectation value of a Wilson loop W around a contour C: V (R) = − lim T→∞ 1 T ln〈WC [A]〉 (3.36) In the computation of such a quantity one has to deal with three sources of divergences. A first one is a logarithmic divergence which is absorbed in the definition of the coupling constant; and additional logarithmic divergence is due to the presence of the corners and a third is linear and related to the so-called residual mass of the quarks [29]. In order to remove the corner divergences one can consider a Creutz ratio in the form VT (R) = log ( W (R, T − 1) W (T,R) ) (3.37) which in turn can yield the static potential via V (R) = lim T→∞ VT (R). (3.38) Such a technique has already been applied in NSPT in order to compute the perturbative expansion of the residual quark mass [25] [26]. In fact, the ratio (3.37) in the limit of (3.38) leaves us with potential, including the linearly divergent term which gives the residual quark mass: V (R) = δm− CF αV (R) R . (3.39) Eq (3.39) defines the renormalized coupling in the so-called potential scheme αV (R) we will be concerned with. Till now we have made no use of the fact that we are going to compute all the quantities in the lattice scheme, apart from having explicitly written the linear divergence δm. The strategy is the following: we compute the potential via ratios of Wil- son loop in lattice perturbation theory. Then, by separating the linear from the logarithmic divergence, we actually compute the coupling in the potential scheme (αV ) as an expansion in the lattice coupling (αTLS), i.e. (3.29). The perturbative computation of the static quark potential in our lattice scheme reads V (R) = 2δm− CF αTLS R ( 1 + c1(R)αTLS + c2(R)α 2 TLS +O(α3TLS) ) , (3.40) 37 Renormalization coefficients where the coefficients (3.30) and (3.31) now read c1(R) = 2b0 log ΛV ΛTLS + 2b0 logR (3.41) c2(R) = c1(R) 2 + 2b1 logR + 2b1 log ΛV ΛTLS + b (V ) 2 − b(TLS)2 b0 . (3.42) As a byproduct of the computation we also obtain the residual mass δm as an expansion in αTLS ∑ n≥0 δm(n)αn+1TLS. (3.43) 3.2.2 Computation In order to compute the potential we measured the Wilson loopsW (R, T ) for all the values of R and T up to 16 on a 324 lattice and averaged results over ∼150 lattice configurations. The (sea) quark mass was set to zero by plugging the appropriate mass counterterm for the perturbative critical mass [23]. We extracted both the residual quark mass and the coefficient we are interested in by fitting (order-by-order) our data to the expected form defined by eq. (3.40). The extraction of the potential magnify the statistical errors in the com- putation of the Wilson loops, which must be taken into account together with the systematic errors due to the finite lattice size nature of our compu- tation: we can not actually take neither the T →∞ limit nor the continuum limit. In other words, we expect that results are distorted by lattice artifacts, and taking these effect under control is the major issue of the computation. While the overall volume was kept fixed, both the R/T and R scales were relevant and there is no simple finite size scaling that could be set up within these constraints. In order to minimize these effects we consider an interval of R such that • T > R (T/R ∼ 2.5); • R itself is not too small (R ≥ 3); • the fitting intervals themselves are from 3 up to 7 points long. The results so obtained were regarded as approximations of the infinite vol- ume, lattice artifacts free results. Results are of course affected by both systematic (lattice artifacts) and statistical errors. The relative weight of these effects is different for different 38 Matching αTLS to αMS orders. This actually opens the way to looking for a careful tradeoff between the errors. We adopt the following strategy: when systematic effects are clearly dis- tinguishable (i.e. statistical effects are relatively small), we only consider T = 16 data. This is the case of the tree level potential (figure 3.1): the different VT (R) are separated curves (though they intersect within errors). 6 7 8 9 10 113.44 3.46 3.48 3.5 3.52 3.54 3.56 3.58 R V( R) (0) T = 16 T = 15 T = 14 T = 13 Figure 3.1: Tree level potential. The different lines represent the VT (R) ex- tracted form different values of T , from T = 13 (higher) to T = 16 (lower). When statistical errors are significant on their own, we need a different approach. In figure 3.2 we show the two loop potential. VT (R) for different values of T widely intersect within errors; in other terms, the systematic (finite T ) effect is not that clear. In this case we decided to neglect this systematic effect and tame the statistical noise by averaging over different values of T (starting from T = 14 to T = 16). In other words, we thus obtain smoother curves. 3.2.3 Results In order to verify the reliability of this computation we first checked known results. The correctness of this is not trivial, since at any order residual mass is unknown, and we get it as a byproduct. We estimate 3 ≤ R ≤ 7 as the best fitting interval for tree level and one loop; the same interval was also taken for two loop computation. In figure 3.3 we show tree level data, to be fitted to the functional form V (R)(0) = 2δm(0) − CF R . (3.44) 39 Renormalization coefficients 4 6 8 10 12 14 16 40 45 50 55 60 65 70 R V( R) (2) T = 16 T = 15 T = 14 T = 13 Figure 3.2: Two loop potential. In this case statistical errors are higher than systematic one, making hard to distinguish them. 2.5 3 3.5 4 4.5 5 5.5 6 6.5 3.1 3.15 3.2 3.25 3.3 3.35 3.4 3.45 3.5 R V( R) (0) fit T = 16 Figure 3.3: Data and fit (continuous line) for tree level potential. We obtain δm(0) = 1.84±0.01, while CF is reconstructed to a few percent. This gives a rough idea of the impact of systematic effects. In figure 3.4 we show one loop data. At one loop we are able to extract the constant term log ΛV Λ(TLS) , which we can compare to the analytical result: V (R)(1) = δm(1) − CF R 2b0 ( logR + log ΛV ΛTLS ) . (3.45) 40 Matching αTLS to αMS 2.5 3 3.5 4 4.5 5 5.5 6 6.5 8.2 8.4 8.6 8.8 9 9.2 9.4 9.6 9.8 10 R V( R) (1) fit T = 16 Figure 3.4: One loop potential: data and fit. We obtain log ΛV Λ(TLS) = 2.8 ± 0.1, to be compared to the analytical result 2.8191 [24]. We also obtain δm(1) = 5.71± 0.01. At two loop we finally tackle the determination of the quantity we are interested in ( b (V ) 2 −b (TLS) 2 b0 ): V (R)(2) = δm(2) − CF R ( c1(R) 2 + 2b1 logR + 2b1 log ΛV ΛTLS + b (V ) 2 − b(TLS)2 b0 ) . (3.46) In figure 3.5 one can see that two loop fluctuations are larger than at lower orders, and as a consequence the fit will suffer of a larger indetermination. In this case we obtain δm(2) = 30± 1 and b (V ) 2 − b(TLS)2 b0 = 4± 1. (3.47) Though the relative error in this value is high, we must emphasize that we are interested in the whole coefficient of the α3 term, in particular in the final matching to αMS αMS = αTLS + 2.79866α 2 TLS + 11.5(±1.0)α3TLS +O(α4TLS), (3.48) where the relative error on the second coefficient is less than 10%. While we are trying to increase our accuracy, one should keep in mind that what we are in the very end interested in is the impact of this matching in the logarithmic contribution to the ZO. 41 Renormalization coefficients 2.5 3 3.5 4 4.5 5 5.5 6 6.5 39 40 41 42 43 44 45 46 47 48 49 R V( R) (2) fit mean Figure 3.5: Two loop potential: data and fit. 3.3 Chiral limit The renormalization procedure requires to fix a normalization point in the space of parameters. In order to keep the computation simpler one typically chooses the point at zero quark mass. This choice, in the case of lattice simulations, has a drawback. Construction of fermionic observables requires the inversion of Dirac ma- trix. In non-perturbative simulation this is a heavy computational task, since the computational time required in the inversion of a matrix is related to the distance between minimum and maximum eigenvalue. This distance grows with the decreasing of the quark masse. In practice one computes operators at different values of the mass and then extrapolate toward the chiral limit. In the perturbative case the situation is different: the mass enters as a parameter which is possible to set to zero at the end of the computation with- out introducing indeterminations. Is this the end of the problem related to the chiral limit? Not exactly. As previously enlighted (3.9) Wilson fermions break chiral invariance introducing an (irrelevant) mass term. In order to keep the chiral limit, a mass counterterm (the critical mass) is required at all orders. The critical mass can be computed from the inverse quark prop- agator (3.9). The self energy Σˆ(pˆ, mˆcr, β−1) can indeed be expanded in the eUntil a few years ago was supposed to exist an effective “limit”, the so-called Berlin wall, in decreasing the quark mass. Nowadays the situation is not that bad as we used to think. 42 Chiral limit form Σˆ(pˆ, mˆcr, β −1) = Σˆc(pˆ, mˆcr, β −1) + ΣˆV (pˆ, mˆcr, β −1) + Σˆo(pˆ, mˆcr, β −1) (3.49) where Σˆc(pˆ, mˆcr, β−1) is the contribution along the identity and ΣˆV (pˆ, mˆcr, β−1) the component along the γ matrices: Σˆc(pˆ, mˆcr, β −1) = 1 4 TrspinΣˆ(pˆ, mˆcr, β −1) (3.50) ΣˆV (pˆ, mˆcr, β −1) = 1 4 ∑ µ γµTrspin ( γµΣˆ(pˆ, mˆcr, β −1) ) . (3.51) Finally, Σˆo(pˆ, mˆcr, β−1) contains all the irrelevant contribution along all the other elements of the Dirac basis, which we are not interested in. The term Σc contains the contribution to the critical mass: Σˆ(0, mˆcr, β −1) = Σˆc(0, mˆcr, β −1) = amcr, (3.52) and by restoring the physical dimension one can recognise its a−1 divergence a−1Σˆc(0, mˆcr, β −1) = mcr. (3.53) 0 1 2 3 4 5 6 7 8 −2 −1.8 −1.6 −1.4 −1.2 −1 −0.8 −0.6 −0.4 −0.2 0 (pa)2 m cr(1) L = 32 L = 16 L = 12 (a) 0 1 2 3 4 5 6 7 8 −2 −1.8 −1.6 −1.4 −1.2 −1 −0.8 −0.6 −0.4 −0.2 0 (pa)2 m cr(2) L = 32 L = 16 L = 12 (b) Figure 3.6: 1 (a) and 2 (b) loop Σˆc(pˆ, mˆcr, β−1) vs (ap)2. Since 1 and 2 loop critical mass for the scheme defined by Symanzik gauge action and nf = 2 quark are analytically known we plug their values as counterterms in the computation. Figure (3.6) shows that we can keep the chiral limit at these orders. At 3 loop the mass counterterm for such regularization is not yet known, and in figure (3.7) new result for this value is shown. We will enter deeper into the results later on. Here we were only concerned with the conceptual point that there is no practical difficulty in staying in the chiral limit. 43 Renormalization coefficients 0 1 2 3 4 5 6 7 8 0 0.5 1 1.5 2 2.5 3 3.5 4 (pa)2 m cr(3) L = 32 L = 16 L = 12 Figure 3.7: 3 loop Σˆc(pˆ, mˆcr, β−1) vs (ap)2. 3.4 Hypercubic Taylor expansion and contin- uum limit In a lattice computation the finite size a of the lattice spacing gives rise to artifacts which must vanish in the continuum limit. Such artifacts appear in a form which is dictated by the breaking of the continuum Poincaré symme- try, in particular the O(4) group of rotation. Having less constraints to fulfill, the number of terms which satisfy the residual hypercubic symmetry is larger. Let us now discuss the expected form of a renormalization constant in a lattice computation. To keep notation simple, we consider a one loop computation of Zq f Zq = 1 + α0(c (1) 1 log(pa) + c (0) 1 (pa)) + . . . (3.54) c (0) 1 (pa) entails lattice artifacts. What we can expect is that these arti- facts comply to the lattice H4 symmetry. If we can constrain their expected form, we can fit them and extract the continuum limit, as it will be clear in the following. The form of these terms obviously depend on the observable we are in- terested in. For example, scalar invariants relevant for observables like the fWe write momentum p for the scale µ. 44 Hypercubic Taylor expansion and continuum limit critical mass are S (0) ≡ 1 (3.55) S (2) 1 ≡ pˆ2 = a2 ∑ µ p2µ S (4) 1 ≡ pˆ4 = a4 ∑ µ p4µ S (4) 2 ≡ (pˆ2)2 = a4 (∑ µ p2µ )2 S (6) 1 ≡ pˆ6 = a6 ∑ µ p6µ S (6) 2 ≡ pˆ2pˆ4 = a6 (∑ µ p2µ )(∑ ν p4ν ) S (6) 3 ≡ (pˆ2)3 = a6 (∑ µ p2µ )3 . . . where dots means that one can goes to higher powers in p. Actually, due to 45 Renormalization coefficients the form of Feynman integral (ratios of trigonometric functions) one expects S (0) ≡ 1 (3.56) S (2) 1 ≡ pˆ2 = a2 ∑ µ p2µ S (2) 2 ≡ pˆ 4 pˆ2 = a2 ∑ µ p 4 µ∑ µ p 2 µ S (4) 1 ≡ pˆ4 = a4 ∑ µ p4µ S (4) 2 ≡ (pˆ2)2 = a4 (∑ µ p2µ )2 S (4) 3 ≡ pˆ 6 pˆ2 = a4 ∑ µ p 6 µ∑ µ p 2 µ S (6) 1 ≡ pˆ6 = a6 ∑ µ p6µ S (6) 2 ≡ pˆ2pˆ4 = a6 (∑ µ p2µ )(∑ ν p4ν ) S (6) 3 ≡ (pˆ2)3 = a6 (∑ µ p2µ )3 . . . As a consequence, if one looks for a formal Taylor expansion for a typical scalar observable S(pˆ), one has to expect a form like S(pˆ) = c0S(0) + c1S(2)1 + c2S(2)2 + c3S(4)1 + c4S(4)2 + c5S(4)3 + . . . (3.57) Because of the presence of the powers of a, all but the first term (c0S(0)) vanish as a → 0. This defines a convenient strategy to get the continuum limit: by fitting the expansion of S(p) we get the term we are interested in. A different set of invariants must be considered when also vector repre- sentations of the hypercubic group are in place. For example, the lattice generalization of the continuum γµpµ structure reads V(pˆ) = i ∑ µ γµpˆµ (V(0)(pˆ) + pˆ2µV(1)(pˆ) + pˆ4µV(2)(pˆ) + . . .) . (3.58) At this point the V(i) coefficients can be expanded by means of Hypercubic 46 Hypercubic Taylor expansion and continuum limit Taylor expansion with the invariant defined before in a form similar to (3.57) V(i)(pˆ) = c(i)0 S(0)+c(i)1 S(2)1 +c(i)2 S(2)2 +c(i)3 S(4)1 +c(i)4 S(4)2 +c(i)5 S(4)3 +. . . (3.59) The main difference with respect to the scalar case is that we need to be more careful with power-counting (in a). For example O(1) = c(0)0 S(0) (3.60) O(a2) = c(0)1 S(2)1 + c(1)0 p2µS(0) (3.61) O(a4) = c(0)2 S(4)1 + c(0)3 S(4)2 + p2µc(1)1 S(2)1 (3.62) (3.58) is just the relevant form for Zq. This has to be computed according to (3.3), which in turns requires to get ΣˆV (pˆ, mˆcr, β−1) out of the self energy (3.49). By computing the inverse propagator at many different values of pa we can extract V(0)(pˆ) + pˆ2µV(1)(pˆ) + pˆ4µV(2)(pˆ) + . . . (3.63) This is shown in figure (3.8). One can recognise different “families” of points, corresponding to the different lengths of pµ 0 1 2 3 4 5 −0.7 −0.6 −0.5 −0.4 −0.3 −0.2 (pa)2 Z q( 1) L = 32 L = 16 L = 12 Figure 3.8: 1 loop ΣˆV (pˆ, mˆcr, β−1) vs (ap)2. As a general recipe, all the possible covariant polynomials can be found via a character’s projection of the polynomial representation of the hypercubic group onto the defining representation of the group. The only term which does not vanish in the continuum limit is c(0)0 , the first coefficient of the expansion of V(0). 47 Renormalization coefficients 3.5 Infinite volume approximation In the previous section we have shown how one can overcome the effects due to the discretization. The second effect we have to take into account is the finite volume of the lattice. Such technique has already been applied to the computation of the gluon and ghost propagators [30] [31]. Let us take ΣˆV (pˆ, mˆcr, β −1) as an example. In the spirit of [32] consider the ansatz Σˆ (0) V (pˆ, pL) = Σˆ (0) V (pˆ) + ( Σˆ (0) V (pˆ, pL)− Σˆ(0)V (pˆ) ) = Σˆ (0) V (pˆ) + ∆Σˆ (0) V (pˆ, pL), (3.64) where Σˆ(0)V (pa) = Σˆ (0) V (pa,∞) and we omit the dependence on mˆcr and β−1 in order to keep the formula simple. Dimensional arguments suggest the pL dependence. Note that this is not an irrelevant effect: these effects would be there also in a continuum computation on a finite volume. One can then apply this ansatz to the expansion (3.58) Σˆ (0) V (pˆ, pL) = i ∑ µ γµpˆµ (V(0)(pˆ) + pˆ2µV(1)(pˆ) + pˆ4µV(2)(pˆ) + . . .)+∆Σˆ(0)V (pˆ, pL) (3.65) which is, however, not yet useful for a fitting procedure. To make the ansatz useful for a fit we should in principle perform an Hypercubic Taylor expansion both on the V(i)(pˆ) and on the term responsible for the finite volume effect. As a first approximation one can only consider ∆Σˆ(0)V (pˆ = 0, pL) ≡ ∆Σˆ(0)V (pL). Neglecting pa corrections is motivated by regardime them as corrections on corrections. With such an approximation (3.65) becomes Σˆ (0) V (pˆ, pL) = i ∑ µ γµpˆµ (V(0)(pˆ) + pˆ2µV(1)(pˆ) + pˆ4µV(2)(pˆ) + . . .)+∆Σˆ(0)V (pL) (3.66) which is amenable to a fit. The key observation is the following pµL = 2πnµ L L = 2πnµ (3.67) where nµ is one element of the 4-tuple ~n = (n1, n2, n3, n4) with integer nµ. The interpretation is the following: neglecting pa effects in∆Σˆ(0)V (pˆ, pL) term, measurements at momenta defined by the same 4-tuple ~n at different lattice sizes N = L/a are affected by the same pL effect. The practical implementation is the following. First of all one has to select a collection of lattice sizes (e.g. N = 32, 20, 16, 12) and an inter- val [(pa)2min, (pa) 2 max]. Then one considers all the tuples ~n = (n1, n2, n3, n4) 48 Infinite volume approximation which fall in [(pa)2min, (pa) 2 max] for all value of N . Actually one chooses one representative for all the tuples connected by an H4 transformation. This choice so far still does not fix the overall normalization of the data. For that purpose we notice that finite volume effects decrease with both increasing momentum squared and increasing lattice size. Therefore we further add to the selected data the measurement taken at the highest value of (pa)2 which falls in the interval [(pa)2min, (pa) 2 max] on the lattice with the biggest value of N = L/a. Assuming that this tuple (which we will call n∗) is a good approx- imation to pL =∞, this measurement will be considered as a normalization point. In the end one fits data coming from different lattice sizes (N) and gets both the parameters in the relevant Hypercubic Taylor expansion and the finite volume corrections ∆Σˆ(0)V (pL), which come in the same number of the considered tuples. To be explicit, (3.66) can then be formulated in terms of the (~n,N) de- pendence: keep in mind that in the combined fit pa depends on n and N . At a low order Σˆ (0) V (~n 6= ~n∗, N) = c(0)0 S(0) + c(0)1 S(2)1 (~n,N) + c(0)2 pµ(nµ, N)2S(2)2 (~n,N) +c (1) 0 p 2 µ(nµ, N) 2 S (0)(~n,N) + ∆Σˆ (0) V (~n) +O(a4),(3.68) where ~n∗ is the tuple taken as normalization point. For the measure taken at the tuple ~n∗ the formula is slightly different: Σˆ (0) V (~n ∗, N) = c (0) 0 S (0) + c (0) 1 S (2) 1 (~n ∗, N) + c (0) 2 pµ(nµ, N) 2 S (2) 2 (~n,N) +c (1) 0 pµ(n ∗ µ, N) 2 S (0)(~n∗, N) +O(a4). (3.69) The technique can be applied without modification to higher loop cases. Before fitting the data in a chosen momentum range, all logarithmic pieces (supposed to be known) have to be subtracted. As a general recipe, in order to make this procedure effective, some re- quirement have to be fulfilled: • a number of measurements on different lattice sizes (N = L/a) should be available; • the interval [(pa)2min, (pa)2max] should allow an effective Hypercubic Tay- lor expansion; • the number of data points which enter in the procedure has to be large enough with respect to the number of fit parameters. Being an effective fitting procedure, based on the assumptions (3.64) and (even more) (3.66), its effectiveness has to be assessed a posteriori. In par- ticular one should check that: 49 Renormalization coefficients • the data point associated to ~n∗ is virtually free of finite size effects; • all the procedure results in a sufficiently stable fit. As a last remark, we can go back to eq. (3.65). We discussed the case in which one neglects pa effects in ∆Σ(0)V (pˆ, pL). This is not a necessary step: one can of course take into account a pa expansion also for this term. This would of course increase the number of parameters in the fit. 3.6 Results Aim of this work is to compute renormalizations constants for the field and scalar (S), vector (V), pseudoscalar (P) and axial (A) currents. As a byprod- uct we obtain the quark critical mass counterterm. Partial results for the tensor current are available, not yet analyzed at all. All the results are at three loop. Measurements were taken after having fixed Landau gauge. These renormalization constants are actually known in the literature from the ETMC non-perturbative determinations [14]. We want to understand to which extent these results and the errors that go on top of them can be fully trusted. In particular, the cited paper presents different determina- tions of the RC’s, which are supposed to differ for irrelevant effects. One of the method put forward in [14] is actually supposed to deliver a continuum limit result. Among the different current, the pseudoscalar (which is a very relevant one, being connected to the quark mass in the Twisted Mass regu- larization) appears to be both quite distant from the one loop perturbative result and quite unstable with respect to the choice of different methods to determine it. 3.6.1 Strategies From a practical point of view, we measure the observables Oq = ΣV and the OΓ which have been defined in (3.7). Recalling the renormalization conditions on themg Z−1q Oq = 1 Z −1 q ZΓOΓ = 1 (3.70) one can consider a variety of strategies. 1. The measurement of the field renormalization coefficient is straightfor- ward from the definition, Oq = Zq. (3.71) gWe move to a sloppy, but easy notation. 50 Results In this case the logarithmic divergences are related to the field anoma- lous dimension. 2. One can directly consider the observable OΓ: OΓ = Zq ZΓ . (3.72) One thus obtains a ratio of RC’s and the logarithmic divergences derive from both the field and the current anomalous dimensions. 3. One can aim at a direct measurement of the current renormalization constant ZΓ = Zq OΓ = Oq OΓ . (3.73) In this case logarithmic contributions, if present, are due to the current anomalous dimension. 4. Observables which are always logarithmic safe can be computed from ratios between scalar and pseudoscalar currents or vector and axial ZPS = ZP ZS = OS OP ZV A = ZV ZA = OA OV (3.74) These different approaches can be used as a test of consistency for the mea- surements. One can for example extract the renormalization coefficients from (3) and compare with the ratio of (1) and (2), or similarly compute coeffi- cients for two currents from (3) and compare to the ratio (4). A lot of combinations are possible. We will apply both Hypercubic Taylor expansions alone, and the com- bined fit Hypercubic Taylor expansions/pL effects in order to show the ef- fectiveness of the technique. As a general rule, we must keep in mind that a Hypercubic Taylor expansion is viable only once one has subtracted the logarithms coming from anomalous dimensions. In the case of the scalar current, for example, the observable one measures is OS = Zq ZS . (3.75) Inspecting the right hand side one can recognise the presence of logarithms coming both from the field (3.19) and the scalar (3.25) renormalization con- stant. One loop can be inspected in figure (3.9). What one measures is 51 Renormalization coefficients 0 0.2 0.4 0.6 0.8 1 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 pa2 Z q( 1) − Z S( 1) unsubtracted subtracted analytic result Figure 3.9: Comparison between unsubtracted O(1)s and subtracted O (1) s − γ (1) S log (pa) 2 measurement of scalar current. supposed to be O(1)s = Z (1) q − Z(1)S + γ(1)S L2 (3.76) where Z(1)q and Z (1) S have to be intended as in (3.19) and (3.25) and we traded γ(1)S L 2 for 2γ(1)S L. Before we attempt a Hypercubic Taylor expansion, the observable must be subtracted of the log: O(1)s − γ(1)S L2 = Z(1)q − Z(1)S (3.77) This is a quantity which is supposed to be amenable to a fit. We recall that vector (V) and axial (A) are finite, so (2) suffers only from logarithmic contributions from the field anomalous dimension and (3) is finite. 3.6.2 The simplest measurement: critical mass We first start computing the simplest scalar quantity, the critical mass. In this case the renormalization coefficient is additive, and does not involve all the intricate structure we discuss before: the measured quantity is directly the counterterm we need in order to stay in the chiral limit. We will show the effect of finite lattice size and volume, seeing the strategy explained in section 3.4 and 3.5 at work. 52 Results As a first step we will reproduce the 1 and 2 loop results: in this case the measured critical mass should be zero, since the relevant counterterms are already present, having been plugged into our simulation. 0 0.2 0.4 0.6 0.8 1 −0.5 −0.4 −0.3 −0.2 −0.1 0 (pa)2 m cr(1) L = 32 L = 20 L = 16 L = 12 fit pa fit pa + pL (a) 0 0.2 0.4 0.6 0.8 1 −0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 (pa)2 m cr(2) L = 32 L = 20 L = 16 L = 12 fit pa fit pa + pL (b) Figure 3.10: One (a) and two (b) loop Σˆc(pˆ, mˆcr, β−1) vs (ap)2. In figure 3.10 we show in detail the results for critical mass at one and two loop for the different lattice sizes we are considering, L/a = N = 324, 204, 164 and 124. Different shapes of blue points correspond to measurements on the different lattice sizes, the green line is the extrapolation toward (pa)2 = 0 for N = 324 and the red line is the combined fit over the different lattice sizes. Results are presented in tab. (3.2). The first column is relative to the lattice sizes of the data used for the fit. The first three results are obtained on a 324 lattice. With “pL” we refer to fit on 324, 204 and 164 lattice sizes, and in “pL1” we further added some data from the 124 lattice. The second column refers to the (pa)2 interval of data to be fitted. The choice depends on the observable and on the order in a of the fit. A higher order fit requires many points, but holds at higher values of (pa)2 and delivers more precise results. It is clear from both figures that the simple Hypercubic fit is not enough to reproduce the analytic result, even if the green line fits the blue circles very well. We can see the effect of the ∆Σm(pL) terms. The lowest value of (pa)2 corresponds to the (1, 1, 1, 1) tuple on 324 (it is actually a circle): it is actually below the red (fit) curve. The square point which is around (pa)2 = 0.4 corresponds to the same (1, 1, 1, 1) tuple on 204: the deviation from the red (fit) curve is to a good approximation the same. This validates our ∆Σm(pa, pL) ∼ ∆Σm(pL) assumption (and as a matter of fact the fit is quite good). Moving toward higher values of (pa)2 reduces the finite volume effects, and 53 Renormalization coefficients Fit (pa)2 interval Order m(1)cr m (2) cr 32 (0-0.6) a2 -0.0307 -0.0534 32 (0-0.9) a4 -0.0173 -0.0376 32 (0-1.9) a6 -0.0102 -0.0286 pL (0-1.55) a4 -0.0057 -0.0112 pL (0-2.1) a6 0.0018 0.0027 pL1 (0-1.9) a4 -0.0061 -0.0137 pL1 (0-2.2) a6 0.0018 0.0027 Table 3.2: 1 and 2 loop critical mass, vanishing because of the presence of the counterterms. Precision increases as the order of the pa fit is increased and as the combined pa, pL fit is taken into account. the two lines (green and red) virtually overlap. Deviations from zero of these one and two loop results are a measure of the indeterminations which are in place (quality of data, effectiveness of fits). Having shown the effectiveness of the technique we can estimate three loop critical mass (figure 3.11). 0 0.2 0.4 0.6 0.8 1 2.8 3 3.2 3.4 3.6 3.8 4 (pa)2 m cr(3) L = 32 L = 20 L = 16 L = 12 fit pa fit pa + pL Figure 3.11: Three loop critical mass vs (ap)2. We performed different determinations, using N = 324 lattice alone or applying “pL” and “pL1” strategies. Such different values can be used to estimate the accuracy. Results of the different fitting procedures are shown in tab. (3.3). 54 Results Fit (pa)2 interval Order m(3)cr 32 (0-0.6) a2 3.818 32 (0-0.9) a4 3.860 32 (0-1.1) a6 3.878 pL (0-1.55) a4 3.919 pL (0-2.1) a6 3.958 pL1 (0-1.9) a4 3.921 pL1 (0-2.2) a6 3.958 Table 3.3: Three loop critical mass. Taking into account pL effect, the value is slightly higher than the value obtained on the 324 alone. Notice that the 324 determinations are systematically lower than results coming from fitting many lattice sizes. In the pL case, there is a systematic effect coming from the order of the fit. χ2 values are always reasonable. We can assume that the indetermination coming from different fit orders is a fair estimate of errors, so that we can state m(3)cr = 3.94(2). Fitting the critical mass has been a laboratory to understand our fitting procedures, in particular with respect to single/multiple lattice sizes and pa and pL effects. As a matter of fact, we will see that some extra care is needed when anomalous dimensions are in place. The bottom line is that the log’s that we subtract are strictly speaking not supposed to be there on a finite lattice, and as a consequence pL effects are in a sense expected to be more severe. As a general, pragmatic rule, one has to face the clash of two requirements. On one side, one is interested in staying at the lowest possible value of (pa)2 (in order to have a more effective Taylor expansion), but at the same time one should avoid the points which mostly suffer from pL effects, which would drive us just the other way around. As said, these pL effects are mainly present in divergent observables (S and P), while finite ones (V and A) are less affected. As in the combined fit the pL effects are taken into account, one wants to consider data also at low (pa)2. 3.6.3 Field renormalization coefficients Measurement of Zq at two and three loop is affected by large pL effects, while at one loop one can safely fit the low momentum region as one can inspect by observing figure 3.12. As a matter of fact, one loop anomalous dimension is zero in Landau gauge. In the figure the red curve is V(0)(pˆ) in (3.63). 55 Renormalization coefficients One should keep this explanation in mind every time a plot with different “families” is in place. 0 0.2 0.4 0.6 0.8 1 −0.66 −0.64 −0.62 −0.6 −0.58 −0.56 −0.54 −0.52 −0.5 −0.48 pa2 Z q( 1) measured Oq (1) analytic low (pa)2 Figure 3.12: One loop field renormalization coefficient fitted on the 324 lattice. In black the analytical result. See the main text for explanation of the red curve. The comparison between the different fit is shown in tab. 3.4. At lowest order analytic result can be compared to our measurement. Once again, we trust the pL determinations, which are quite close to each other (we leave a2 determination out) and result in good χ2 values: I quote Z(2)q = −0.614(2) and Z(3)q = −1.184(4). Fit (pa)2 interval Order Z(1)q Z (2) q Z (3) q 32 (0.8-1.1) a2 -0.6234 -0.5709 -1.035 32 (0.6-1.1) a4 -0.6544 -0.5630 -1.204 pL (0-.7) a2 -0.6453 -0.6043 -1.166 pL (0-1.2) a4 -0.6509 -0.6148 -1.187 pL1 (0-1.55) a4 -0.6528 -0.6137 -1.181 pL1 (0-1.86) a4 -0.6545 -0.6134 -1.183 Analytic -0.6598 Table 3.4: Result of the fit for Zq. 3.6.4 Current renormalization coefficients We have just seen in the case of one loop field renormalization constant that, if the anomalous dimension is zero, one can safely fit data at low (pa)2 values. This turns out not to be true anymore when one considers a logarithmically 56 Results divergent quantity. As an example, we can consider once again the one loop scalar current (3.77) after the subtraction of the log. 0 0.2 0.4 0.6 0.8 1 −0.02 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 pa2 Z q( 1) − Z S( 1) measured Oq (1) analytic fit H4 (a) 0 0.2 0.4 0.6 0.8 1 −0.02 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 pa2 Z q( 1) − Z S( 1) measured Oq (1) analytic fit H4 (b) Figure 3.13: One loop scalar current fitted at low ( 0 < (pa)2 < 0.4 ) and high ( 0.6 < (pa)2 < 1.1 ) values of (pa)2 on the 324 lattice. In black the analytical result. Figure (3.13) shows this effect in a fit according to the (3.72) method. The low-momenta data are clearly affected by pL effects: a naive fit (figure (a)) wouldn’t return correct results. A “quick and dirty” procedure should be to consider only data at high momenta (figure (b)). On the other hand, we know that the Hypercubic Taylor Expansion fitting procedure has a worse behaviour in this region. A combined fit accounts for pL effects at low values of (pa)2, reducing the indeterminacy (figure 3.14) 0 0.2 0.4 0.6 0.8 1 −0.02 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 pa2 Z q( 1) − Z S( 1) measured OS (1) analytic fit H4 Figure 3.14: One loop scalar current fitted by means of combined fit on 324, 204, 164 and 124 lattice. In black the analytical result. 57 Renormalization coefficients In the following we will actually prefer the measurement ZΓ = Oq OΓ , leav- ing (3.72) only as a consistency cross-check. Scalar and pseudoscalar current suffer from logarithmic divergences from current anomalous dimension, while vector and axial current are finite and thus they do not need any subtraction. Once again, we try to combine different approaches. We first of all fit data measured on lattice 324. In this case, S and P currents can be fitted taking into account only the high (pa)2 region (this is the “quick and dirty” procedure we have already referred to), while V and A can rely on the safe procedure of staying on fairly low (pa)2 values. Actually the “quick and dirty” procedure can be implemented in two ways: either we only consider 324 data or we collect the same highest momenta (remember the discussion on pL effects: what actually matter are the “plain” tuples) on all the sizes. We notice that the resulting fits often make a reasonable sense. One example is in figure (3.15) : one and two loop for P current in the (3.72) approach. 0 1 2 3 4 5 6 7 8 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 Figure 3.15: One and two loop pseudoscalar current fitted by means of com- bined fit on 324, 204, 164 and 124 lattice.Only high tuples are taken into account to tame pL effects. Consistent results are recovered. On top of this, we perform combined fits to many lattice sizes according to the “pL” and “pL1” strategies we described when we discussed the critical mass fits. In the case of axial current results are available only on 324 and 204. Nev- ertheless, in tab. (3.8) is possible to notice that the results of the combined fit are consistent with result on 324. As already said, as a general recipe I tried to cross-check all the results by comparing numbers coming from different strategies. We already saw there 58 Results Fit (pa)2 interval Order Z(1)S Z (2) S Z (3) S 32 (0.8-1.1) a2 -0.6690 -0.7469 -1.706 32 (0.6-1.1) a4 -0.6866 -0.7462 -1.703 pL (0-0.7) a2 -0.6804 -0.7829 -1.855 pL (0-1.2) a4 -0.6889 -0.8121 -2.022 pL1 (0-1.55) a4 -0.6849 -0.7794 -1.866 pL1 (0-1.86) a4 -0.6844 -0.7769 -1.823 32,20 (0-0.7) a2 -0.6816 -0.7802 -1.793 32,20 (0-1.78) a4 -0.6846 -0.7855 -1.801 pL (0-1.9) a4 -0.6876 -0.7968 -1.914 pL1 (0-2.0) a4 -0.6835 -0.7730 -1.823 20,16,12 (0-2.01) a4 -0.6801 -0.7524 -1.737 Analytic -0.6893 Table 3.5: ZS Fit (pa)2 interval Order Z(1)P Z (2) P Z (3) P 32 (0.8-1.1) a2 -1.067 -1.192 -2.637 32 (0.6-1.1) a4 -1.094 -1.135 -2.531 pL (0-0.7) a2 -1.095 -1.184 -2.603 pL (0-1.2) a4 -1.110 -1.201 -2.739 pL1 (0-1.55) a4 -1.102 -1.184 -2.632 pL1 (0-1.86) a4 -1.103 -1.171 -2.561 32,20 (0-0.7) a2 -1.097 -1.165 -2.467 32,20 (0-1.78) a4 -1.095 -1.193 -2.593 pL (0-1.9) a4 -1.107 -1.203 -2.684 pL1 (0-2.0) a4 -1.105 -1.188 -2.661 20,16,12 (0-2.01) a4 -1.083 -1.182 -2.737 Analytic -1.101 Table 3.6: ZP 59 Renormalization coefficients 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 −0.8 −0.7 −0.6 (pa)2 Z V( 1) analitic L = 32 L = 20 L = 16 L = 12 pµ = 0 (a) 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 −0.65 −0.6 −0.55 −0.5 −0.45 (pa)2 Z A( 1) analitic L = 32 L = 20 pµ = 0 (b) 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 −0.76 −0.74 −0.72 −0.7 −0.68 −0.66 −0.64 −0.62 −0.6 −0.58 −0.56 (pa)2 Z S( 1) analitic L = 32 L = 20 L = 16 L = 12 pµ = 0 (c) 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 −1.1 −1.05 −1 −0.95 −0.9 −0.85 (pa)2 Z P( 1) analitic L = 32 L = 20 L = 16 L = 12 pµ = 0 (d) Figure 3.16: One loop data and fit results for RC’s measured by means of (3.73). 60 Results Fit (pa)2 interval Order Z(1)V Z (2) V Z (3) V 32 (0 - 0.6) a2 -0.8294 -0.8849 -1.864 32 (0 - 0.85) a4 -0.8299 -0.8968 -1.884 pL (0 - 1.55) a4 -0.8360 -0.9173 -1.897 pL1 (0 - 1.86) a4 -0.8383 -0.9107 -1.913 pL1 (0 - 2.1) a4 -0.8376 -0.8957 -1.874 Analytic -0.8411 Table 3.7: ZV Fit (pa)2 interval Order Z(1)A Z (2) A Z (3) A 32 (0 - 0.6) a2 -0.6349 -0.6131 -1.210 32 (0 - 0.85) a4 -0.6292 -0.6168 -1.221 pL (0 - 1.6) a2 -0.6460 -0.6097 -1.178 pL (0 - 1.8) a4 -0.6376 -0.6151 -1.202 Analytic -0.6352 Table 3.8: ZA are a number in place. All in all, χ2 were often very good. It was reassuring to notice that when χ2 values were comparable, results were comparable too. In particular, I could extract the values of the coefficients ZPS, ZSP , ZV A and ZAV either from direct measurement of ratio of observables (3.74) or from the coefficients obtained before. Notice that in the case of these ratios (which are finite) the approximation inherent in our pL fits were actually larger than finite size effects. In the end, measurements of ratios from the 324 lattice were the most reliable. Actually, we have some evidence that pL fits do improve when we also try to take into account (pa)2 corrections on top of pL corrections. Figure 3.17 shows a one loop example of (Zq/ZO)(1). After all the cross-checks, I quote as results • Z(2)S = −0.782(5) and Z(3)S = −1.82(7); • Z(2)P = −1.19(4) and Z(3)P = −2.59(9); • Z(2)V = −0.896(12) and Z(3)V = −1.88(3); • Z(2)A = −0.618(9) and Z(3)A = −1.21(3). 61 Renormalization coefficients 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 (pa)2 (Z q/Z S)( 1) analitic L = 32 L = 20 L = 16 L = 12 pµ = 0 (a) 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0.18 0.19 0.2 0.21 0.22 0.23 0.24 (pa)2 (Z q/Z V)( 1) analitic L = 32 L = 20 L = 16 L = 12 pµ = 0 (b) Figure 3.17: Data and fit result for RC’s measured by means of (3.72) Fit (pa)2 interval Order (Zq/ZS)(1) (Zq/ZS)(2) (Zq/ZS)(3) 32 (0.7-1.1) a2 0.0417 0.2123 0.863 32 (0.5-1.1) a4 0.0137 0.1383 0.619 pL (0-1.5) a4 0.0378 0.2007 1.004 pL1 (0-1.55) a4 0.0339 0.1947 0.889 pL1 (0-1.86) a4 0.0330 0.1825 0.888 Analytic 0.0295 Table 3.9: Zq/ZS Fit (pa)2 interval Order (Zq/ZP )(1) (Zq/ZP )(2) (Zq/ZP )(3) 32 (0.7-1.1) a2 0.4422 1.223 3.617 32 (0.5-1.1) a4 0.4227 1.203 3.399 pL (0-1.5) a4 0.4587 1.206 3.759 pL1 (0-1.55) a4 0.4526 1.182 3.804 pL1 (0-1.86) a4 0.4516 1.172 3.670 Analytic 0.4422 Table 3.10: Zq/ZP 62 Results Fit (pa)2 interval Order (Zq/ZV )(1) (Zq/ZV )(2) (Zq/ZV )(3) 32 (0 - 0.6) a2 0.1678 0.3907 1.126 32 (0 - 0.85) a4 0.1644 0.3744 1.268 pL (0 - 1.55) a4 0.1826 0.4317 1.237 pL1 (0 - 1.86) a4 0.1839 0.4271 1.231 pL1 (0 - 2.1) a4 0.1818 0.4296 1.238 Analytic 0.1813 Table 3.11: Zq/ZV Fit (pa)2 interval Order (Zq/ZA)(1) (Zq/ZA)(2) (Zq/ZA)(3) 32 (0 - 0.6) a2 -0.0263 -0.0365 -0.0494 32 (0 - 0.85) a4 -0.0372 -0.0591 -0.0715 pL (0 - 1.6) a2 -0.0009 0.0142 0.0937 pL (0 - 1.8) a4 -0.0169 -0.0142 0.0319 Analytic -0.0246 Table 3.12: Zq/ZA Fit (pa)2 interval Order Z(1)V A Z (2) V A Z (3) V A 32 (0 - 0.6) a2 -0.1944 -0.3957 -1.024 32 (0 - 0.85) a4 -0.1996 -0.4012 -1.042 pL (0 - 1.6) a2 -0.1842 -0.3846 -0.994 pL (0 - 1.8) a4 -0.1990 -0.4024 -1.046 Analytic -0.2059 Table 3.13: ZV A Fit (pa)2 interval Order Z(1)AV Z (2) AV Z (3) AV 32 (0 - 0.6) a2 0.1945 0.4328 1.184 32 (0 - 0.85) a4 0.1976 0.4452 1.214 pL (0 - 1.6) a2 0.1840 0.4162 1.133 pL (0 - 1.8) a4 0.1969 0.4411 1.211 Analytic 0.2059 Table 3.14: ZAV 63 Renormalization coefficients Fit (pa)2 interval Order Z(1)PS Z (2) PS Z (3) PS 32 (0-0.6) a2 -0.4088 -0.7790 -1.992 32 (0-0.9) a4 -0.4103 -0.7809 -1.998 32 (0-1.1) a6 -0.4018 -0.7870 -1.997 pL (0-1.55) a4 -0.4171 -0.8026 -2.073 pL (0-2.1) a6 -0.4165 -0.8357 -2.073 pL1 (0-1.9) a4 -0.4171 -0.8034 -2.114 pL1 (0-2.2) a6 -0.4221 -0.8267 -2.149 Analytic -0.4117 Table 3.15: ZPS Fit (pa)2 interval Order Z(1)SP Z (2) SP Z (3) SP 32 (0-0.6) a2 0.4088 0.9449 2.692 32 (0-0.9) a4 0.4102 0.9489 2.707 32 (0-1.1) a6 0.4110 0.9476 2.705 pL (0-1.55) a4 0.4176 0.9769 2.779 pL (0-2.1) a6 0.4204 0.9741 2.818 pL1 (0-1.9) a4 0.4183 0.9807 2.844 pL1 (0-2.2) a6 0.4207 0.9806 2.883 Analytic 0.4117 Table 3.16: ZSP 0 0.5 1 1.5 2 2.5 3 −0.24 −0.22 −0.2 −0.18 −0.16 −0.14 −0.12 −0.1 (pa)2 (Z V/ Z A )(1 ) analitic L = 32 L = 20 pµ = 0 (a) 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 −0.42 −0.4 −0.38 −0.36 −0.34 −0.32 −0.3 −0.28 (pa)2 (Z P/ Z S )(1 ) analitic L = 32 L = 20 L = 16 L = 12 continuum (b) Figure 3.18: Data and fit result for RC’s measured by means ratios of ob- servables. 64 Results 3.6.5 Conclusions, i.e. summing the series We have now access to three loop RC’s, and we claimed that this would help to keep under control truncation errors in summing the series. Is this true? We can compare the summations of our results with non-perturbative computations performed by the European Twisted Mass Collaboration [14]. Actually Twisted Mass fermion action differs from Wilson fermion action, but in the chiral limit, where the renormalization scheme is defined, they reduce to the same. In figure (3.19) we display the value of summation at each lattice mo- menta, i.e. we leave out all the extrapolation process which we had at work before in extracting the continuum limit. These figures should actually be compared to figure (7) of [14] (they are actually obtained at the same value for the coupling). In order to improve convergence properties, one can rephrase the series in terms of different couplings. In figure (3.20) the summation of RC’s in different coupling is presented. The different coupling constants we used are obtained in terms of the basic plaquette P. We define x0 = β−1, x1 = − 1 1.4649 log (P ), x2 = −β(−1)P . As a general recipe, one would like to see better and better convergence as order increases: differences of loop i from loop i−1 has to decrease as i increases. On the other hand, at some point determina- tions from different couplings should be the same. Notice with this respect that the x2 coupling tends to destabilize the series for certain observables: as a thumb rule, one sees this effect when “smaller” couplings already display good apparent convergence properties. All in all, there are cases in which our three loop results match very well (possibly by boosting the coupling) to non-perturbative results: this is the case for ZV , ZA and ZQ. Notice that in this case the extrapolation that the authors of [14] attempt in their figure (9) seems to be well under control. On the other side, we obtain quite different results in the case of ZS and ZP . Once again, we nevertheless notice that comparisons have ac- tually to be performed keeping figure (9) of [14] in mind. In these cases, the extrapolation does not seem in such a brilliant shape. Most probably there is some distance still to go in summing the series, i.e. one needs even higher loopsh. Nevertheless, if we have to gain insight from our experience, figure (9) of [14] should be better understood in terms of both finite volume hExploratory four loops results for finite quantities are actually on their way 65 Renormalization coefficients 0 0.5 1 1.5 2 2.5 3 0.7 0.72 0.74 0.76 0.78 0.8 0.82 0.84 0.86 pa2 Z V 0 0.5 1 1.5 2 2.5 3 0.78 0.8 0.82 0.84 0.86 0.88 0.9 pa2 Z A 0 0.5 1 1.5 2 2.5 3 0.7 0.72 0.74 0.76 0.78 0.8 0.82 0.84 0.86 0.88 0.9 pa2 Z S 0 0.5 1 1.5 2 2.5 3 0.55 0.6 0.65 0.7 0.75 0.8 pa2 Z P 0 0.5 1 1.5 2 2.5 3 0.8 0.82 0.84 0.86 0.88 0.9 0.92 0.94 pa2 Z P S 0 0.5 1 1.5 2 2.5 3 0.76 0.78 0.8 0.82 0.84 0.86 0.88 0.9 0.92 0.94 pa2 Z q Figure 3.19: Resummation in the bare coupling of observables at one (blue), two (red) and three (black) loop for lattices 324 (circles) and 204 (diamonds) 66 Results effects and higher orders (pa)i effects. Notice that if one tries to trust all their momentum window with an higher order extrapolation, their results get quite closer to ours. As a further caveat, notice that the authors of [14] do not try to keep the effects of what we called “families” under control (they simply perform some cuts to control hypercubic deviations from continuum; this is a common strategy, but also this effect should be assessed). As a last remark, one should assess also the effects in (pa)i of the chiral extrapolations. In the end, my understanding is that the errors quoted in the continuum limit of [14] are quite underestimated. All in all, truncation errors are once again not the end of the story. 67 Renormalization coefficients 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 0.42 0.66 0.68 0.7 0.72 0.74 0.76 0.78 0.8 couplings re su m m e d Z V 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 0.42 0.74 0.76 0.78 0.8 0.82 0.84 0.86 couplings re su m m e d Z A 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 0.42 0.68 0.7 0.72 0.74 0.76 0.78 0.8 0.82 0.84 0.86 couplings re su m m e d Z S 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 0.42 0.54 0.56 0.58 0.6 0.62 0.64 0.66 0.68 0.7 0.72 0.74 couplings re su m m e d Z P 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 0.42 0.76 0.78 0.8 0.82 0.84 0.86 0.88 0.9 couplings re su m m e d Z P S 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 0.42 0.88 0.89 0.9 0.91 0.92 0.93 0.94 0.95 couplings re su m m e d Z V A Figure 3.20: Resummation of the RC’s to one (red), two (green) and three (blue) loop. We show resummations for different couplings: on the x-axis, the (different) values of the different couplings x0, x1 and x2 68 CHAPTER 4 The Dirac operator spectrum: a perturbative approach Spontaneous symmetry breaking is a common phenomenon in many branches of Physics. Probably the most known example of symmetry breaking is spon- taneous magnetization. Let us consider a 3-d spin system with two (continuum) degrees of free- dom. At high temperature this system is rotationally invariant. Even if an external magnetic field is applied, its effect disappears as it is removed. De- creasing the temperature something happens: the spins promptly polarise in the direction of a magnetic field and rotational invariance is lost even after the magnetic field is removed. The mathematical interpretation of this phe- nomena is that the ground state is no longer invariant under the whole O(3) group of rotations in space, but only under the O(2) subgroup of rotations in the plane perpendicular to the applied field. A similar phenomena is a key feature of QCD, the so-called chiral sym- metry breaking. The presence of a mass m in the QCD Lagrangian L = ∑ f ψ¯f ( /D + ig /A ) ψf + ∑ f mf ( ψ¯fψf ) , (4.1) breaks the invariance under the independent left (L) and right (R) transfor- mations ψfL → ∑ f ′ Uff ′ L ψ f ′ L ψ f R → ∑ f ′ Uff ′ R ψ f ′ R (4.2) 69 The Dirac operator spectrum Which one would have in the massless limit. In the notation above f runs over the flavours and ψL = 1 2 (1− γ5)ψ ψR = 1 2 (1 + γ5)ψ. (4.3) There is a very intuitive way of capturing the physics of such a phe- nomenon: we can say that a small quark mass leads to a macroscopic re- alignment of the QCD vacuum (this is a strict quotation from [33]). Since the QCD partition function reads Z = 〈 ∏ f det(D +mf)〉 = 〈 ∏ f ∏ n (iλn +mf )〉 (4.4) one would naively conclude that this is hardly possible. In a more precise way: in order for this to be possible there must be an accumulation of Dirac operator eigenvalues near zero; otherwise the effect of a small quark mass would be overwhelmed by much larger eigenvalues. This is just the message encoded in the Banks-Casher relation [34] 〈ψ¯ψ〉 = πρ(0) V (4.5) This is the explicit relation of the chiral condensate (the order parameter of the transition associated to spontaneous symmetry breaking) to the density of eigenvalues of the Dirac operator spectrum, defined as ρ(λ) = 〈 ∑ n δ(λ− λn)〉. (4.6) Although not a natural observable in Field Theory, the Dirac operator spectrum has in force of (4.5) become a natural probe for the chiral transi- tion. Recent work [35] has investigated the field theoretic status of spectral observables, in particular with respect to their renormalization properties. From a numerically point of view, it should be pointed out that Lattice QCD can quite naturally compute (4.6), once a lattice regularization of the Dirac operator is given. A typical study that has been extensively performed is the comparison of Dirac lattice spectra to what is conjectured by Random Matrix models. It is quite amazing that there has been no attempt at com- puting the Dirac operator spectrum in Perturbation Theory. This is from one point of view easy to understand: this is kind of Physics that is assumed to be non-perturbative. On the other hand, Perturbation Theory can at least provide a natural mechanism for the (re)arrangement of the spectrum. 70 The Dirac spectrum and Perturbation Theory 4.1 The Dirac spectrum and Perturbation The- ory Since the free Dirac operator has a vanishing eigenvalues density near zero, one is lead to the conclusion that the small eigenvalues are due to gauge interactions. There is actually a natural candidate: any quantum interac- tion produces a repulsion among the eigenvalues. This is just what we were anticipating a few lines above. Basically, one starts with a highly degen- erate eigenvalue problem and look at how degeneracy is lifted: we expect that eigenvalues repel each other. One should say that with this respect Perturbation Theory is in a tantalizing situation: • on one side, PT sits (deep) in the chirally restored regime, while one looks for an effect which lives at its border; • on the other side, PT could give in principle a unique opportunity to follow the fate of eigenvalues in their mutual repulsion. Notice that at least a couple of other concerns should be well taken into ac- count. First of all, perturbative series always show an asymptotic behaviour. Second, we are working with Wilson fermion, whose chiral proprieties are not brilliant at all. We want to emphasize that our work is still at a very preliminary stage. In particular, we do not want to address here the subtleties which arise in properly defining a perturbative expansion of (4.6). We will discuss a quick and dirty procedure in which we first compute the perturbative corrections to the free spectrum eigenvalues λn = λ (0) n + β −1/2λ(1)n + β −1λ(2)n + . . . and then resum the expansion at given values of the coupling β. Given these summations, we can proceed to compute a density of eigenvalues. This is actually much the same work as in non-perturbative computations of the spectrum. 4.2 The Dirac Spectrum: what to look at Since the operator is antihermitean, what one usually do is to compute not actually its eigenvalues, but the eigenvalues of the hermitean operator γ5D or D†D. We choose the second form, even if nothing prevents us from computing the former. 71 The Dirac operator spectrum We recall the form of the (massive) Wilson Dirac operator in configuration space on the lattice: D = (m+ 4) ∑ x ψ¯xψx −1 2 ∑ µ [ ψ¯x (1− γm)Ux+µψx+µ + ψ¯x+µ (1 + γm)U †xψx ] , (4.7) where we choose the Wilson parameter to be one. It is straightforward to see the degeneracy of the free operator D†D we referred to before in momentum space. The operator reads M(p) ≡ D(p)†D(p) = ( m+ 1 2 p˜2 − iγp¯ )( m+ 1 2 p˜2 + iγp¯ ) = ( m+ 1 2 ∑ µ p˜2µ )2 + ∑ µ p¯2µ (4.8) where p˜µ = sin (pµ/2) and p¯µ = sin (pµ). To explicitly recognise the presence of degeneracies, we first point out that, since lattice simulations are performed in a finite volume, one should choose boundary conditions. The typical choice is to have periodic boundary conditions in spacial directions and antiperiodic boundary conditions in the temporal direction a. Since the components of momentum pµ in (4.8) appear squared, the oper- ator should be invariant under the transformation pµ → −pµ both in spatial and temporal direction. In force of the boundary condition it holds that −pµ = L − pµ, where L the extent of the lattice in the µˆ direction. Then it should be true that, in the free case, there is a degeneracy between all the eigenvectors containing the values pµ and L − pm. Such a degeneracy can be graphically computed by drawing a polygon in which vertices are lattice points, placed in slightly different positions depending on a direction being periodic or not. In (figure 4.1(a)) and (figure 4.1(b)) the couples with the same x value will be degenerate. A second source of degeneracy can be easily understood because of the absence of a µ dependence in (4.8). This means that all the momenta which differ for a permutation of the periodic components are degenerate. For example, the tuples (nt, nx, ny, nz) and (nt, ny, nx, nz) correspond to the same eigenvalue. aThis is indeed the case of fermionic field: gluonic field is periodic in the spatial direction too. 72 The Dirac Spectrum in NSPT 0 12 3 4 5 (a) 0 1 2 3 4 5 (b) Figure 4.1: Example of pattern of degeneracy for periodic (a) and antiperi- odic (b) boundary conditions in the case L=6. The dashed lines connect the degenerate couples. 4.3 The Dirac Spectrum in NSPT Making use of the intrinsic perturbative nature of NSPT, one has to solve a typical eigenvalue/eigenvector problem in PT. In practice, by plugging the NSPT expansions of gluon field U (n)µ into (4.7) one finds the order-by-order perturbation to be solved. The operator can then be decomposed in the free term M0 and the perturbation N , which is expanded in orders: M = M0 +N =M0 + ∑ i>0 N (i). (4.9) In order to solve for the spectrum we have to set up the eigenvalue problem M |α〉 = ǫ |α〉 (4.10) which has to be solved by ǫ = ǫ0 + gǫ1 + g 2ǫ2 + . . . |α〉 = ∣∣α(0)〉 + g ∣∣α(1)〉 + g2 ∣∣α(2)〉 + . . . (4.11) Due to the (huge) degeneracy of the free field solution, for every eigenvalue we need to explicitly separate components inside and outside the starting (degenerate) eigenspace |α〉 = P ∗ |α〉 + Pin |α〉 + Pout |α〉 (4.12) 73 The Dirac operator spectrum where P ∗ |α〉 = |α0〉 is the projection along the free (degenerate) eigenspace singled out as the zeroth order of the solution; Pin is the projector onto the component of the free eigenspace which is orthogonal to |α0〉 ; Pout projects instead outside the free eigenspace. In order to compute corrections to eigenvalues and eigenvectors, one should write the eigenvalue equation (4.10) in terms of the projectors P , i.e. (ǫ− ǫ0 −N) |α0〉 + (ǫ− ǫ0 −N)Pin |α〉 + (ǫ−M0 −N)Pout |α〉 = 0 (4.13) Applying once again the projectors (ǫ− ǫ0 − P ∗N) |α0〉 = 0 (4.14) (ǫ− ǫ0 − PinN)Pin |α〉 = (PinNPout |α〉 + PinN |α0〉) (4.15) (ǫ−M0 − PoutN)Pout |α〉 = (PoutNPin |α〉 + PoutN |α0〉) (4.16) one can find the formal solution for the eigenvectors correction. Once multiplied on the left for 〈α0| equation (4.10) gives the correction to the eigenvalues: ǫ− ǫ0 = 〈α0|N |α〉 ⇒ ǫn = n∑ k=0 〈α0|Nk−n |αk〉 . (4.17) Eigenvectors correction can be computed by means of (4.15) and (4.16): Pout |α〉 = (ǫ−M0 − PoutN)−1 (PoutNPin |α〉 + PoutN |α0〉) (4.18) Pin |α〉 = (ǫ− ǫ0 − PinN)−1 (PinNPout |α〉 + PinN |α0〉) . (4.19) This is the (closed) solution only provided degeneracy is lifted at first order. Should this not be the case, the formalism should be generalized by intro- ducing a new projector for each level of degeneracy still present (the solution is nevertheless closed also in such a situation, which actually occurs in our computations). In standard non-perturbative LGT computations of the Dirac operator spectrum one gets distributions of eigenvalues by generating configurations and computing the spectrum on each of them. The density of eigenvalues is then simply obtained by plain histograms of the results. We stress once again that we will adhere to the naive recipe of first computing the eigenvalues in PT, then summing the expansions at given values of the coupling and finally constructing histograms much the same way as in the non-perturbative case. 74 Results 4.4 Results In figure 4.2 we plot examples of our results: we collect all the measurements for first (trivial) and second (one loop) order corrections to free field results for the second lowest lying eigenspace on a 64 lattice. We stress that this eigenspace is degenerate (the dimension of this eigenspace is 144), but on top of this degeneracy the histograms entail the multiplicity which comes from the number of measurements. −1.5 −1 −0.5 0 0.5 1 1.5 0 20 40 60 80 100 120 140 160 180 O(β−1/2) correction to second eigenvalue on 64 −4 −3.8 −3.6 −3.4 −3.2 −3 0 50 100 150 200 250 300 350 O(β−1) correction to second eigenvalue on 64 Figure 4.2: First (trivial) and second (one loop) corrections to the second (lowest lying) free field eigenvalue on a 64 lattice (overall distributions of the measures). Figure 4.3 displays data once the average over all the measurements has been taken. In this case we plot first order corrections (as one expects, they average to zero) for lowest lying and second lowest lying eigenvalues. There are issues which are worth stressing. First of all, one can inspect degenera- cies which are not lifted. Second, the distributions of corrections in the two eigenspaces differ quite a lot. Figure 4.4 displays another interesting feature. In this case we plot a third order correction, which enlights how higher orders display long tails. One probably needs to carefully assess when the free field degeneracy is ac- tually lifted, as it is clear from the impact of denominators in (4.18). With this respect we point out that one can always check the accuracy of the computation by considering quantities like 75 The Dirac operator spectrum 0 5 10 15 20 25 −0.25 −0.2 −0.15 −0.1 −0.05 0 0.05 0.1 0.15 0.2 0.25 O(β−1/2) corrections to lowest lying free eigenvalue (64; free eigenspace degeneracy is 24) 0 50 100 150 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 O(β−1/2) corrections to second lowest lying free eigenvalue (64; free eigenspace degeneracy is 144) Figure 4.3: First (trivial) corrections to the first (lowest lying) and second free field eigenvalue on a 64 lattice (averages over Langevin histories). Free eigenspace degeneracies are 24 (left) and 144 (right). 〈Tr(D†D)k〉 = . . . 〈Tr(D†D)−k〉 = . . . They can be both computed directly and reconstructed from the eigen- values distribution, eventually validating the latter. We can now go back to our quick and dirty procedure to inspect the im- pact of the perturbative corrections. Basically, we can sum the contributions at any given value of the coupling and try to follow the resulting distribution of eigenvalues as the gauge interaction comes into play. We plot in figure 4.5 what we get at one loop. Figure 4.5 is something like a sequence of pictures taken while the interac- tion is switched on. One starts at zero coupling, where the key feature of the free field is on display: bins are centered where free field eigenvalues sit, and 76 Results −10 −8 −6 −4 −2 0 2 4 6 8 10 0 50 100 150 200 250 300 350 400 450 O(β−2/3) correction to second eigenvalue on 64 Figure 4.4: Third correction to the second free field eigenvalue on a 64 lattice: while it is centered in zero (as expected), it displays long tails. bins height simply entails the degeneracy of the various eigenspaces. Notice anyway that at this resolution some bins actually result from the contribu- tion of two free field eigenvalues sitting very close to each other. While the interaction is switched on (i.e. the value of the inverse coupling β decreases) the bins spread and overlap and eventually a non-zero density near zero is generated. A natural question arises: where do eigenvalues moving to zero come from? One should remember the point we made on repulsion among eigenvalues. Figure 4.6 displays an example of how this takes place: we plot the contribution to ρ coming from two eigenvalues starting very close to each other in free field. It is worth better assessing the impact of the repulsion among the couple of free eigenvalues we have just looked at. It actually turns out that they give a substantial contribution to the rearrangement of the eigenvalues density: one can recognize their splitting on the right of figure 4.7. Black line in figure 4.7 is nothing but another way of plotting the first 77 The Dirac operator spectrum 0 5 10 15 20 25 30 0 2000 4000 6000 8000 β =∞ 0 5 10 15 20 25 30 0 2000 4000 6000 8000 β = 18 0 5 10 15 20 25 30 0 2000 4000 6000 8000 β = 14 0 5 10 15 20 25 30 0 2000 4000 6000 8000 β = 10 0 5 10 15 20 25 30 0 2000 4000 6000 8000 β = 6 Figure 4.5: The evolution of the eigenvalues density: from free field limit (β =∞) to the interacting case (at different values of the coupling β). 0 2 4 6 8 10 12 14 16 18 0 50 100 150 β =∞ 0 2 4 6 8 10 12 14 16 18 0 50 100 150 β =60 0 2 4 6 8 10 12 14 16 18 0 50 100 150 β =30 0 2 4 6 8 10 12 14 16 18 0 50 100 150 β =10 Figure 4.6: Following the repulsion of two eigenvalues on 64. They start very close in free field limit and then strongly repel each other. row of figure 4.5: we plot the first 5122 free eigenvalues and the length of each plateaux is just the degeneracy of each free field eigenspace. The su- 78 Conclusions 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 −5 0 5 10 15 20 The values of the first 5122 free field eigenvalues (black) on 64 and how they are corrected at one loop (red) at β =7.5 Figure 4.7: The first 5122 free eigenvalues on a 64 lattice (black line): the lengths of each segment is the degeneracy in free field. Red curve displays how they move at one loop at β = 7.5. perimposed red line shows the summation (at first loop) of the perturbative series for these eigenvalues at β = 7.5. Some caveats are of course in order: • first of all: this bare Perturbation Theory and one should renormalize results. Next step is actually the computation of the relevant renor- malization constant. • Is this a finite-volume effect? At the moment we have actually got the same qualitative picture at any (still moderate) size we studied. • Is this a finite a effect? Testing this is more difficult. • One should carefully take care of the order of limits which is in place in the Banks Casher relation. 4.5 Conclusions There is actually no conclusion I can draw at the moment. Still, even though at a very preliminary stage, I showed some results which display that I set up an environment for a perturbative computation of the Dirac operator spectrum by means of NSPT. The only thing that I can state is that bare 79 The Dirac operator spectrum corrections to the spectra quantitatively support the picture of the repul- sion among eigenvalues being responsible for a rearrangement of eigenvalues. Whether this ultimately contributes in any way to Banks Casher is a ques- tion that will be addressed only after the relevant renormalization constant has been computed. 80 Conclusions My PhD research project was in NSPT. As it was clear, the levels of com- pleteness of different subjects were very different. Even if I resisted to the temptation of letting a report on the numerical work take over Physics, I think I was able to suggest that the package which was the result of my activities (PRlgt) is something quite complete and useful. The computation of three loop renormalization constants was a hard task. The issue of perturbative versus non-perturbative renormalization deserved nevertheless attention. My conclusion is that there are in general more in- determinations in non-perturbative determinations than it is usually stated: truncation errors are not the end of the story. In particular, our three (or even higher) computations are valuable to better assess the overall system- atics. I also pinned down the prospects for an NSPT study of the Dirac oper- ator spectrum. The relevant question (to which extent can one inspect the reshuﬄing of Dirac operator eigenvalues due to color interactions?) can be probably really addressed once we have computed the relevant renormaliza- tion constant. All in all, I think I can fairly state that NSPT can both provide fenomeno- logically relevant high loop computations and enable to enbark in quite non- standard computations. 81 Ringraziamenti Questa tesi segna la conclusione non solo di una esperienza scientifica, ma soprattutto una forte esperienza umana. Ho ancora bene in mente le sen- sazioni provate vedendo sfilare Milano al mio fianco il giorno in cui mi sono trasferito per iniziare questo Dottorato, così come la sorpresa nello scoprire questa fantastica città che è Parma. É quindi doveroso ringraziare le persone che, nei modi più diversi, hanno contribuito a farmi crescere. Prima di tutti, desidero ringraziare il Prof. Francesco Di Renzo, senza di lui niente di tutto questo sarebbe stato possibile. Fin dal primo giorno ha mostrato una grandissima cordialità, che si è rinnovata giorno dopo giorno. Gli insegnamenti che ho tratto da lui sono sicuramente molto più vasti di quanto esposto in questa tesi, ma spero soprattutto di essere riuscito ad assorbire almeno un po’ della sua carica umana. Grazie veramente di tutto! Ringrazio il Prof. Giuseppe Marchesini, che ho avuto modo di ammirare sotto tutti gli aspetti e conoscere durante i miei precedenti studi universitari a Milano. Credo di essere stato profondamente segnato dalla sua grande cu- riosità ed intuizione scientifica, e non smetterò mai di considerarlo un esempio da “inseguire”. Non posso certamente esimermi dal ringraziare il Dott. An- drea Banfi, che mi ha insegnato molto durante e dopo la Laurea. Nel poco tempo che ho avuto per poter collaborare con lui ho conosciuto un’ottima persona ed un bravissimo ricercatore. Spero un giorno di poter avere una nuova occasione per poter lavorare insieme, sarebbe certamente stimolante in quanto avrei ancora molto da imparare da lui. Un ringraziamento è poi dovuto a tutte le persone del gruppo di Parma con cui ho avuto modo di interagire, in primis il Prof. Roberto Alfieri ed il Dott. Roberto Covati con i quali ho avuto modo di interagire e collaborare, il Prof. Enrico Onofri e la signora Liliana Superchi, che si sono sempre mostrati disponibili e cordiali. 83 Ringraziamenti Tante altre persone sono entrate nella mia vita parmense in questi tre anni, ma nessuno prepotentemente come i compagni di ufficio: Riccardo, Daniele, Alessio, non è stato certo facile convivere per gran parte della gior- nata tutti insieme, ma ci siamo decisamente divertiti. Penso (e con una punta di orgoglio spero) che ne passerà di tempo prima che quest’ufficio veda altri anni come questi! Grazie a voi e a tutti gli altri amici, in particolare Anteo, Yuri e Scion, che hanno reso questa esperienza indimenticabile. Spero che la vita possa regalarci altre occasioni per passare del tempo insieme. Grazie a Stefano, mio coinquilino da tre anni, e a tutti gli amici in giro per il mondo che ho avuto modo di conoscere a scuole e conferenze (la maggior parte è il gruppo dell’indimenticabile Les Houches), anche se difficilmente avranno modo di leggere questi ringraziamenti. La forzata lontananza da Bellusco mi ha fatto sempre più apprezzare la fortuna che ho ad avere amici come Mirko (ragazzo, sono 25 anni che siamo amici, mica roba da poco!), Daniele e tutti gli altri che nominare sarebbe ingiusto nei confronti di chi non verrebbe nominato. Se vivendo distanti e vedendosi una volta ogni tanti mesi si riesce a rimanere così amici, finché sarete a Bellusco quella sarà la mia casa. Un grazie ai miei genitori, solo grazie ai sacrifici che hanno fatto e con- tinuano a fare sono potuto arrivare fin qui ed essere la persona che sono. Un grazie a mio nonno, mancato ormai da diversi anni ma sempre con me: fin da quand’ero piccolo dicevi che sarei andato sulla luna, magari non ci arriverò, ma quello che ho raggiunto mi ha dato (quasi) altrettanta soddisfazione. Una parte della dedica di questa tesi non può essere che per te. Un grazie speciale, infine, va a te, Greta: stare al tuo fianco è stata la forza che mi ha permesso di andare avanti nei momenti difficili, gustare veramente i momenti di gioia ed imparare cos’è la felicità. 84 BIBLIOGRAPHY [1] G. Parisi and Y. Wu. Perturbation theory without gauge fixing. Sci. Sin., 24:483 – 496, 1981. [2] G. G. Batrouni, G. R. Katz, A. S. Kronfeld, G. P. Lepage, B. Svetit- sky, and K. G. Wilson. Langevin simulations of lattice field theories. Phys. Rev. D, 32(10):2736 – 2747, 1985. [3] A. S. Kronfeld. Dynamics of Langevin simulations. Prog. Theor. Phys. Suppl., 111:293 – 312, 1993. Invited chapter to appear in the special supp. ‘Stochastic Quantization’ of Progress of Theoretical Physics. [4] J. Floratos et al. Equivalence of stochastic and canonical quantization in perturbation theory. Nucl. Phys. B, 214(3):392 – 404, 1983. [5] F. Di Renzo, G. Marchesini, P. Marenzoni, and E. Onofri. Lattice per- turbation theory by Langevin dynamics. 1993. [6] D. Zwanziger. Covariant quantization of gauge fields without Gribov ambiguity. Nucl. Phys. B, 192(1):259 – 269, 1981. [7] U. Heller and F. Karsch. One-loop perturbative calculation of Wilson loops on finite lattices. Nucl. Phys. B, 251(2):254 – 278, 1985. [8] M. Luscher and P. Weisz. Efficient numerical techniques for perturbative lattice gauge theory computations. Nucl. Phys. B, 266(2):309 – 356, 1986. 85 BIBLIOGRAPHY [9] K. G. Wilson. Confinement of quarks. Phys. Rev. D, 10(8):2445 – 2459, Oct 1974. [10] H. J. Rothe. Lattice gauge theories: an introduction. World Scientific Pub Co Inc, 2005. [11] G. Martinelli, C. Pittori, Christopher T. Sachrajda, M. Testa, and A. Vladikas. A General method for nonperturbative renormalization of lattice operators. Nucl. Phys. B, 445:81 – 108, 1995. [12] M. Luscher, R. Narayanan, P. Weisz, and U. Wolff. The Schrodinger functional: A Renormalizable probe for nonAbelian gauge theories. Nucl. Phys. B, 384:168 – 228, 1992. [13] B. Blossier et al. Light quark masses and pseudoscalar decay constants from Nf=2 Lattice QCD with twisted mass fermions. JHEP, 04:020, 2008. [14] M. Constantinou et al. Non-perturbative renormalization of quark bilin- ear operators with Nf=2 (tmQCD) Wilson fermions and the tree- level improved gauge action. JHEP, 08:068, 2010. [15] F. Di Renzo, A. Mantovi, V. Miccio, L. Scorzato, and C. Torrero. Two and three loops computations of renormalization constants for lattice QCD. Nucl. Phys. Proc. Suppl., 140:716 – 718, 2005. [16] F. Di Renzo, V. Miccio, L Scorzato, and C. Torrero. Renormalization constants for Lattice QCD: New results from Numerical Stochastic Per- turbation Theory. PoS, LAT2006:156, 2006. [17] F. Di Renzo, V. Miccio, L. Scorzato, and C. Torrero. High-loop pertur- bative renormalization constants for Lattice QCD. I. Finite constants for Wilson quark currents. Eur. Phys. J. C, 51:645 – 657, 2007. [18] F. Di Renzo, L. Scorzato, and C. Torrero. High loop renormalization constants by NSPT: A Status report. PoS, LAT2007:240, 2007. [19] W. Celmaster and R. J. Gonsalves. The Renormalization Prescription Dependence of the QCD Coupling Constant. Phys. Rev. D, 20(6):1420 – 1434, Sep 1979. [20] K.G. Chetyrkin and A. Retey. Renormalization and running of quark mass and field in the regularization invariant and MS-bar schemes at three loops and four loops. Nucl. Phys. B, 583:3 – 34, 2000. 86 BIBLIOGRAPHY [21] J.A. Gracey. Three loop anomalous dimension of nonsinglet quark cur- rents in the RI-prime scheme. Nucl. Phys. B, 662:247 – 278, 2003. [22] C. T. H. Davies, G. G. Batrouni, G. R. Katz, A. S. Kronfeld, G. P. Lepage, K. G. Wilson, P. Rossi, and B. Svetitsky. Fourier acceleration in lattice gauge theories. i. landau gauge fixing. Phys. Rev. D, 37(6):1581 – 1588, Nov 1988. [23] A. Skouroupathis, M. Constantinou, and H. Panagopoulos. Two-loop additive mass renormalization with clover fermions and Symanzik im- proved gluons. Phys. Rev. D, 77(1):014513, Jan 2008. [24] S. Aoki, K. Nagai, Y. Taniguchi, and A. Ukawa. Perturbative renormal- ization factors of bilinear quark operators for improved gluon and quark actions in lattice QCD. Phys. Rev. D, 58(7):074505, Sep 1998. [25] F. Di Renzo and L. Scorzato. The Residual mass in lattice heavy quark effective theory to alpha**3 order. JHEP, 0102:020, 2001. [26] F. Di Renzo and L. Scorzato. The N(f) = 2 residual mass in pertur- bative lattice-HQET for an improved determination of m(b)**anti-MS (m(b)**anti-MS). JHEP, 0411:036, 2004. [27] Y. Schroder. The Two loop static potential. Nucl. Phys. Proc. Suppl., 86:525 – 528, 2000. [28] M. Luscher and P. Weisz. String excitation energies in SU(N) gauge theories beyond the free-string approximation. JHEP, 0407:014, 2004. [29] G. Martinelli and Christopher T. Sachrajda. Computation of the b quark mass with perturbative matching at the next-to-next-to-leading order. Nucl. Phys. D, 559:429 – 452, 1999. [30] F. Di Renzo, E. M. Ilgenfritz, H. Perlt, A. Schiller, and C. Torrero. Two-point functions of quenched lattice QCD in Numerical Stochas- tic Perturbation Theory. (I) The ghost propagator in Landau gauge. Nucl. Phys. B, 831:262 – 284, 2010. [31] F. Di Renzo, E.-M. Ilgenfritz, H. Perlt, A. Schiller, and C. Torrero. Two-point functions of quenched lattice QCD in Numerical Stochas- tic Perturbation Theory. (II) The Gluon propagator in Landau gauge. Nucl. Phys. B, 842:122 – 139, 2011. 87 BIBLIOGRAPHY [32] H. Kawai, R. Nakayama, and K. Seo. Comparison of the lattice param- eter with the continuum parameter in massless QCD. Nucl. Phys. B, 189(1):40 – 62, 1981. [33] J. J. M. Verbaarschot and T. Wettig. Random matrix theory and chiral symmetry in QCD. Ann. Rev. Nucl. Part. Sci., 50:343 – 410, 2000. [34] T. Banks and A. Casher. Chiral symmetry breaking in confining theories. Nucl. Phys .B, 169(1-2):103 – 125, 1980. [35] L. Giusti and M. Luscher. Chiral symmetry breaking and the Banks- Casher relation in lattice QCD with Wilson quarks. JHEP, 0903:013, 2009. 88

Gell-Mann - Low Function in QED

Abstract

Extracted data

Gell-Mann - Low Function in QED

Abstract

Extracted data

Related items

Related items