Keywords: Quantum Chromodynamics, Gribov copies, Non-perturbative QCD

Foremost, I would like to express my gratitude to my advisors Nicolás Wschebor and Matthieu Tissier for their continuous support, for their patience and their guidance. I can not imagine writing this thesis without the confidence they gave me.

I would like to thank my family Alicia, Fernando and Paula for their love, their support and attention and for always prioritizing my education. I am also grateful to my partner, Pablo Lessa, for all these beautiful years together, for the constant encouragement to continue learning physics

Fonctions de corrélation infrarouges de la

Chromodynamique Quantique

Résumé

Le but de cette thèse est l'étude des fonctions de corrélation des théries Yang-Mills dans le régime infrarouge. Il est connu que, à cause de l'invariance jauge, il est nécessaire de fixer la jauge pour calculer des valeurs moyennes analytiquement. La procedure de fixation gauge standard est la procedure de Faddeev-Popov (FP). Le Lagrangien de FP permet de faire des calculs perturbatifs pour la Chromodynamique Quantique dans le régime de hautes énergies dont les résultats sont comparés avec succès avec des expériences. Cependant, dans le régime de basses énergies, il se trouve que la constante de couplage, calculée avec la procedure antérieure, diverge. En conséquence, la théorie des perturbations standard n'est plus valide. D'autre part, les simulations du réseau trouvent que la constante de couplage est finie avec une valeur modérée même dans le régime infrarouge. Ceci suggère qu'il devrait exister une manière de faire des calculs perturbatifs également dans le régime infrarouge.

Cette différence dans la constante de couplage peut être due au fait que la procedure de FP n'est pas bien justifiée dans ce régime parce qu'elle ignore l'existence des copies de Gribov. Afin de prendre en compte les copies de Gribov, nous proposons de modifier le Lagrangien de FP avec un terme massif pour les gluons. Ce Lagrangien est un cas particulier du Lagrangien de Curci-Ferrari. Cette modification est également justifiée par le fait que le réseau trouve un propagateur du gluon qui paraît massive aux basses énergies.

Nous utilisons cette version massive pour calculer à une boucle les fonctions de corrélations à deux et trois points pour une configuration cinématique générale et en dimension quelconque dans la jauge de Landau. En particulier, nous calculons les propagateurs du gluon et du fantôme sans la présence des quarks et on compare nos résultats avec du réseau. Une fois fixés les paramètres de la théorie, on utilise les mêmes pour comparer les vertex fantôme-gluon et trois-gluons en obtenant un très bon accord avec le réseau. Finalement, on inclut les effets des quarks dynamiques pour étudier les propagateurs du gluon, du fantôme et du quark. L'analyse est aussi faite pour le vertex quark-gluon.

On trouve que les comparaisons de notre calcul à une boucle avec les résultat du réseau donnent, en géneral, un très bon accord.

Mots-clefs

Chromodynamique Quantique, Copies de Gribov, QCD non-perturbative.

Funciones de correlación en Cromodinámica Cuántica

Resumen

La presente tesis investiga el comportamiento de las funciones de correlación de las teorías Yang-Mills en el límite infrarrojo. Como es sabido, para poder calcular valores esperados en el marco de estas teorías invariantes gauge, debe realizarse algún procedimiento que permita fijar el gauge. El procedimiento estándar, ultilizado exitosamente para hacer cálculos perturbativos a grandes energías, es el procedimiento de Faddeev-Popov (FP). Los cálculos perturbativos usando el Lagrangiando de Faddeev-Popov para Cromodinámica Cuántica permiten reproducir con éxito los valores obtenidos en los experimentos que involucran la fuerza fuerte a grandes energías.

Sin embargo, los cálculos perturbativos usando el Lagrangiano de FP muestran una constante de acoplamiento divergente en el régimen de bajas energía. Por ende, la teoría de perturbaciones en este régimen no es válida. Por otra parte, las simulaciones Monte Carlo, que no implementan dicho procedimiento, encuentran una constante de acoplamiento que se mantiene finita tomando valores no muy grandes incluso a bajas energías. Este hecho sugiere que debe existir alguna forma de teoría de perturbaciones que pueda ser aplicada en el régimen infrarrojo. Esta diferencia encontrada en el comportamiento de la constante de acoplamiento puede deberse al hecho de que el procedimiento de FP no está completamente justificado a bajas energías. Este último ignora la presencia de copias de Gribov en este régimen.

Para tener en cuenta la existencia de copias de Gribov, proponemos usar una pequeña modificación del Lagrangiano de FP. Esta consiste en agregar un término de masa para los gluones, caso particular del Lagrangiano de Curci-Ferrari. Otra motivación para considerar el término de masa para los gluones es el hecho de que las simulaciones Monte Carlo muestran que el propagador del gluon se comporta a bajas energías como un propagador masivo.

Dicha extensión del Lagrangiano de FP es usada para calcular a un loop las funciones de correlación a dos y tres puntos en el gauge de Landau para una configuración cinemática genérica y en dimensión arbitraria. El primer cálculo está enfocado en el propagador del fantasma y del gluón despreciando el efecto de los quarks. Ambos propagadores son renormalizados usando dos esquemas diferentes. Teniendo en cuenta los efectos del grupo de renormalización comparamos nuestros resultados con los obtenidos por las simulaciones numéricas. Esto permite fijar los dos parámetros de nuestra teoría para luego usarlos en los cálculos del vértice gluón-fantasma y del vértice de tres gluones. Para finalizar, incluimos los efectos de los quarks dinámicos en el estudio a primer orden del propagador de gluón,

Introduction

The strong nuclear force is one of the fundamental forces in nature. It is responsible for holding neutrons and protons together within the atomic nucleus in spite of the electric charge of protons. The particles that feel the strong nuclear force are generically called hadrons, and include not only protons and neutrons but also many other particles such as pions, deltas, etc. The most accepted theory to describe microscopically the strong force is Quantum Chromodynamics (QCD). QCD is a quantized non-abelian gauge field theory invariant under the SU (3) symmetry group. This theory predicts the masses of hadrons and the interaction between them with great accuracy. However, the building blocks of Quantum Chromodynamics are not the hadrons themselves. The fundamental constituents of QCD are the quarks and gluons. Quarks are the fundamental fermions which bound together to form the hadrons.

|u ↑, u ↑, u ↑ However, there is a problem in this state since it is formed by three identical fermions in the same state contradicting Pauli's principle. Therefore, these quark should have another quantum number with at least three different values to be able to distinguish the mentioned up quarks.

Fortunately the study of non-abelian gauge theories was in progress since the Yang and Mills's work in 1954 [START_REF] Yang | Conservation of Isotopic Spin and Isotopic Gauge Invariance[END_REF]. In 1967 Faddeev and Popov quantized this kind of theories [START_REF] Faddeev | Feynman Diagrams for the Yang-Mills Field[END_REF], and in 1971 't Hooft proved they were renormalizable [tH71a]. In 1973, Politzer, Gross and Wilczek computed the renormalization group equation and noticed that the coupling constant falls down at high energies [START_REF] Politzer | Reliable Perturbative Results for Strong Interactions?[END_REF][START_REF] Gross | Ultraviolet Behavior of Nonabelian Gauge Theories[END_REF]. This property was called asymptotic freedom. Politzer, Gross and Wilczek were awarded with a Nobel prize in 2004 for this discovery.

The same year, Coleman and Gross probed that the only theories with asymptotic freedom in dimension four are the non-abelian gauge theories, see [START_REF] Coleman | Price of asymptotic freedom[END_REF]. Gell-Mann and Fritzsch proposed to consider SU (3) as the group to describe the strong interaction [START_REF] Bardeen | Light cone current algebra, pi0 decay, and e+ e-annihilation[END_REF][START_REF] Fritzsch | Advantages of the Color Octet Gluon Picture[END_REF]. Therefore, the corresponding non-abelian gauge theory has an extra non-abelian symmetry which will give to the quarks the extra quantum number needed, the color. For this reason the theory acquired the name Quantum Chromodynamics. Gluons are the bosons who play as mediators of the strong force. Gluons in QCD have the equivalent role as photons in the quantum electrodynamics. Evidence of their existence was found in 1978 in electron-positron annihilation.

Every quark has an extra quantum number called color with values red, green and blue. And gluons are essentially associated with a pair of colors. Neither the flavour or the color are related to the colloquial meaning of the word. Experimental evidence for the existence of color was found, for instance, in the electron-positron annihilation and in the decay rate of the pion into photons. The color is used to describe why quarks and gluon are not found isolated in nature. They always appear forming bound states (hadrons). This property is called confinement. Quantum Chromodynamics describes confinement imposing that only combinations of quark and gluons which together form singlet states (that means that the state has all the three colors in the same quantity, or a quark and an anti-quark with opposite color charge) are the possible hadrons. How this affirmation is deduced from first principles using the QCD Lagrangian is still unsolved. Nowadays, explaining confinement is one of the most ambitious tasks for physicists.

Another important feature in nature that is supposed to be explained by QCD is the spontaneous chiral symmetry breaking. The chiral symmetry allows to interchange the left handed and the right handed quarks when no masses for the quarks are considered. The chiral symmetry is not a real symmetry of the QCD Lagrangian since quarks are massive but it is a good approximation in what concerns the lightest quarks. This symmetry seems to be spontaneously broken and the consequence of this fact is that, in the broken phase, the mass of the hadrons is much larger than the mass of the quarks forming the hadron. For example, as the proton is made of two up quark and one down the sum of their masses gives approximately 9 MeV (m u ∼ 2 MeV and m d ∼ 5 MeV, see [O + 14]) while the mass of the proton is 938 MeV. Most of the mass of the proton arrives from the interactions. The main consequence of chiral symmetry breaking is the generation of 98% of the hadron masses.

We have mentioned three important features of nature that should be explained by QCD. Now, we want to summarize the techniques used to work with this theory. First of all, we have mentioned that the strength of the interaction between quarks and gluons (represented by the coupling constant) goes to zero when the energy scale grows. When the coupling constant is small the system behaves as a weakly interacting theory and therefore its study using perturbation theory is well justified. Perturbation theory is the most useful technique to compute cross sections when the momentum scale is much bigger that the mass of the proton. In this regime a variety of predictions of perturbation theory, such as the running of the coupling constant and the rates of the electron-proton annihilation, were successfully compared with the experimental result.

On the other hand, the low momentum regime is indeed difficult to access and phenomena with characteristic momentum scale less than 1 GeV are not well described using usual perturbation theory. In fact, the coupling constant computed perturbatively in this limit grows significantly. For this reason most of the processes can not be calculated directly with perturbative QCD. Generally the calculation of the cross section is separated in two parts. The first part is the perturbative process at high energies where the computation with quark and gluons are justified. The other one is how to relate the starting hadrons with the quarks and gluons that will be actors of the perturbative process. This is generally done through the universal long-distance functions such as the parton distribution functions, fragmentation functions, generalized parton distribution and other types of form factors. This techniques allowed to test QCD in several processes such as the Drell-Yan process, the rate of hadron production from the annihilation of a electron-positron pair, etc.

To study the infrared (low momenta) regime some alternatives to perturbation theory have been explored. The most reliable tool to study this regime are Monte Carlo simulations [START_REF] Montvay | Quantum fields on a lattice[END_REF]. It consists in simulating the system in a discrete lattice with finite volume. The accuracy of the simulations grows when the volume increases and the lattice spacing decreases. However, increasing the volume and reducing the lattice space increase the cost of the simulations. The challenge is to choose these parameters to make the simulations treatable while keeping the discretization errors small. The lattice simulations allowed to compute observables far away from the perturbative regime. It was possible, for example, to calculate the hadron masses (protons, neutrons, pions, etc.) and the results are in good agreement with experimental results.

Moreover, the results of lattice simulations agree with the principal features of the strong force. For instance, the simulations are consistent with confinement. When a static quark-antiquark pair is studied, the simulations show that when the quarks are far apart, the energy of the configuration is proportional to their distance and therefore grows to infinity. In this case, as the energy needed to keep them apart is huge, it is energetically favourable to create another pair quark-antiquark and form again two pairs that are color singlets as confinement imposes. There is another important characteristic of QCD well reproduced by lattice simulations, which is the spontaneous chiral symmetry breaking. Lattice simulations observe that even when the quark mass in really small in the ultraviolet the mass is enhanced strongly in the infrared.

Even if lattice simulations have had an important role in our understanding of QCD, they are difficult to produce and only the simplest quantities have been computed with them. Moreover, the lattice results can guide us to know how QCD works but in general it does not explain why or how its features depend on each other or which contributions are more important than others.

In consequence, several semi-analytical approaches have been developed to access the infrared properties. The most popular are Dyson-Schwinger equations (DS) [START_REF] Lorenz Von Smekal | The Infrared behavior of gluon and ghost propagators in Landau gauge QCD[END_REF][START_REF] Alkofer | The Infrared behavior of QCD Green's functions: Confinement dynamical symmetry breaking, and hadrons as relativistic bound states[END_REF][START_REF] Aguilar | A Dynamical gluon mass solution in a coupled system of the Schwinger-Dyson equations[END_REF] and the non perturbative renormalization group (NPRG) [START_REF] Berges | Nonperturbative renormalization flow in quantum field theory and statistical physics[END_REF]. Both approaches consist of a infinite set of coupled equations relating correlation functions with Introduction a different number of fields. These sets match with the perturbative expansion in the limit where perturbation theory is valid. However, in the infrared regime, some approximations have to be done in order to solve the system of equations. The idea is to use approximations schemes consistent with the information obtained from lattice simulations. The only problem is that the most accurate lattice simulations are done for observables (i.e. gauge invariant quantities) but the most treatable DS or NPRG equations relate correlation functions which are gauge dependent and therefore gauge fixing is needed.

For this reason, some simulations on the lattice have been done for the correlation functions [START_REF] Cucchieri | Constraints on the IR behavior of the gluon propagator in Yang-Mills theories[END_REF][START_REF] Bogolubsky | Lattice gluodynamics computation of Landau gauge Green's functions in the deep infrared[END_REF][START_REF] Cucchieri | Landau-gauge propagators in Yang-Mills theories at beta = 0: Massive solution versus conformal scaling[END_REF]. These simulations are restricted to the propagators and the three-points correlation functions. As the correlation functions depend on the gauge choice, most of the simulations have been done in the Landau gauge. This gauge is preferred because in the continuum it presents many symmetries and because there is a well established non perturbative definition of this gauge which is easy to implement on lattice simulations.

Fixing the gauge in QCD is necessary to define the correlation functions and in particular the gluon propagator. In the end, for calculations in the continuum we have to use a gauge-fixed Lagrangian instead of the original QCD one. The use of a gauge-fixed Lagrangian to compute gauge invariant quantities gives the same results as the one expected for the standard QCD Lagrangian. The standard procedure to fix the gauge in the continuum is due to Faddeev and Popov [START_REF] Faddeev | Feynman Diagrams for the Yang-Mills Field[END_REF]. The Faddeev-Popov Lagrangian is, therefore, the one used to do perturbation theory at high energies.

However, this procedure is ill-defined in the infrared regime. The issue lies on the existence of Gribov copies [START_REF] Gribov | Quantization of Nonabelian Gauge Theories[END_REF], which are not considered in the Faddeev-Popov procedure. Gribov copies refer to the different configurations of the fields, related through a gauge transformation, that fulfil the same gauge condition. For instance, in Landau gauge we impose ∂ µ A µ = 0 (where A µ represents the gluon field) and this condition presents a family of possible solutions. To fix the gauge we have to choose only one of these solutions. Lattice simulations do not need to fix the gauge but when they do, they can avoid the Gribov problem. Therefore, the fundamental difficulty on fixing the gauge comes when an analytical analysis is performed because it is not known in general how to build a gauge fixed Lagrangian without Gribov copies.

This last problem was faced first by Gribov and then by Zwanziger. They proposed to restrict the functional integrals appearing in the correlation function to what is known as the first Gribov zone. In this way the number of Gribov copies decreases considerably, however, it was shown that there still are Gribov copies inside this region. The restriction to this zone is done through the introduction of new auxiliary fields in an action called the Gribov-Zwanziger action, S GZ . This action was shown to be renormalizable even though it is not known if it is unitary. The results obtained with this action were not completely satisfactory with respect to lattice results when treated in some naive approximation. A different treatment was proposed, which leads to the so called the refined Gribov-Zwanziger action, which reproduces the behaviour of the propagators observed by lattice simulations. The extra fields included by this action make the calculations more complicated and therefore additional approximations have to be used.

Let us stress the importance of the lattice results to guide us in the study of the infrared regime. For example, to solve the system of DS equation the first solution proposed was to consider that every correlation function behaves in the infrared as a power law in the momentum scale. This assumption gives a solution of the DS equation called the scaling solution. However, with the improvement of lattice simulations it was seen that this solution was not coherent with lattice results which instead find a finite gluon propagator and a tree-level like ghost propagator in the infrared regime. The scaling solution cannot reproduce both result simultaneously. Therefore, another solution was looked for called the decoupling or massive solution. This solution gives the same qualitative behaviour for the propagators found by lattice simulations.

Another important observation of lattice simulations is that the coupling constant does not present a Landau pole, this means that it does not go to infinity for a finite value of momentum as it was believed from the perturbative analysis. Moreover, the coupling constant observed at low momenta point out that higher orders in the perturbative analysis should be of order of 20% of the leading contributions even in that regime. This result suggests that it should be possible to consider a perturbative analysis able to reproduce the infrared results. This fact opens the question of which is the most simple Lagrangian that allows to perform a perturbative expansion even in the low momenta regime, to reproduce the lattice results and match with the standard Lagrangian in the ultraviolet. An answer to this question was proposed in [START_REF] Tissier | Infrared propagators of Yang-Mills theory from perturbation theory[END_REF][START_REF] Tissier | An Infrared Safe perturbative approach to Yang-Mills correlators[END_REF]. The Lagrangian consists on a simple modification of the Faddeev-Popov Lagrangian by adding a gluon mass term. This Lagrangian corresponds to a particular case of the Curci-Ferrari Lagrangian proposed in the seventies [START_REF] Curci | A Theory of Massive Yang-Mills Field[END_REF]. This model was successfully proved to be renormalizable [START_REF] Curci | On a Class of Lagrangian Models for Massive and Massless Yang-Mills Fields[END_REF][START_REF] De Boer | On the renormalizability and unitarity of the Curci-Ferrari model for massive vector bosons[END_REF].

The Curci-Ferrari model was discarded at that moment because there were unsatisfactory aspects of theory. For instance, the standard definition of the physical space could not be applied to it because it was shown that it includes states with negative norms. States with negative norm should not belong to the physical space. It is an open question to determine a satisfactory physical space of this model. However, states with negative norm were observed by lattice simulations and also by Gribov-Zwanziger approach in the infrared regime. In particular the transverse gluon propagator present positivity violations, a fact that could be related to confinement [Gri78, Zwa93, DGS + 08a]. We insist on not discarding a model with positivity violation given that the lattice simulations (which is a non perturbative technique) observe the same behaviour. This only shows that the usual definition of physical space is not appropriate.

The aim of this thesis is to prove that, generally speaking, one loop calculation of this model reproduces the results obtained by lattice simulations. For that reason in the first chapter we introduce the basic features of QCD and the standard quantization of it. In the second chapter we treat in detail the problems found at low momenta and possible solutions. In chapter 3 we are going to present the analytical results for the quenched propagators in the model with a gluon mass which compare very well with lattice data [START_REF] Tissier | Infrared propagators of Yang-Mills theory from perturbation theory[END_REF][START_REF] Tissier | An Infrared Safe perturbative approach to Yang-Mills correlators[END_REF]. We also discuss the order of magnitude of higher loop corrections. In the next chapter we compute the quenched three-point correlation functions at one loop. The results fit with great accuracy the lattice data [START_REF] Pelaez | Three-point correlation functions in Yang-Mills theory[END_REF]. In chapter 4 the unquenched gluon propagator and the quark propagator are studied [START_REF] Peláez | Two-point correlation functions of QCD in the Landau gauge[END_REF]. The last chapter complements the unquenched analysis by presenting the quark-gluon vertex and its matching with the lattice results. At the end we present the conclusions of our work.

Chapter 1

Quantum Chromodynamics 1.1 Non-abelian gauge theories

Quantum Chromodynamic is a theory based on the non-abelian gauge group SU (3) which describes the interaction of colored particles. Every particle is described by a local field. Fermions interact between each other by exchanging a gauge boson called gluon. There are 8 such particles which arise from the local gauge symmetry of the theory. The Lagrangian governing this theory is SU (3) invariant. In order to find the Lagrangian of QCD we will work in a more general case considering the gauge group SU (N ) where N generalizes the number of colors.

Let us start discussing which are the terms invariants under the SU (N ) group (a detailed analysis can be found e.g. in [START_REF] Peskin | An Introduction to quantum field theory[END_REF]). We are allowed to include fields belonging to different representations. In each case we have to say how the SU (N ) group acts over those fields. An element of local SU (N ) can be written as U (x) = e -igθ a (x)t a where θ a (x) is an arbitrary function. The Latin letter a goes from 1 to N 2 -1 and t a are the generators of the corresponding Lie algebra which satisfy [t a , t b ] = if abc t c with f abc the structure constants. We choose the generators to be hermitian and normalized as Tr(t a t b ) = δ ab 2 when t a are in the fundamental representation. Contrarily to the deduction presented in [START_REF] Peskin | An Introduction to quantum field theory[END_REF] where the non-abelian gauge invariant Lagrangian is obtained in the Minkowski space, we are going to work in the Euclidean space. The equivalent in Minkowski space can be obtained through a Wick rotation, see [START_REF] Peskin | An Introduction to quantum field theory[END_REF]. We choose to work in the Euclidean space since our calculations will be compared with lattice simulations which are done in the Euclidean space. It is worth mentioning that in this case the position of the "Lorentz" indices is irrelevant.

The fundamental representation of the group is formed by the N -component fields, ψ. When ψ describes a fermion field, this vector is also a spinor under space-time transformations. An element of the SU (N )-group acts over ψ(x) in the fundamental representation with the standard matrix product. The ψ(x) field transforms as

ψ(x) → ψ (x) = U (x)ψ(x)
under the local gauge symmetry SU (N ). Consequently we have that ψ(x) † → ψ † (x) = ψ(x) † U † (x). Therefore, we can build a gauge invariant terms as ψ † ψ. However, this is not Lorentz invariant. Instead, we will use ψ = ψ † γ 0 , where γ 0 is one of the Euclidean Dirac matrices γ µ which satisfy {γ µ , γ ν } = 2δ µ,ν . The first term gauge and Lorentz invariant that we can build is ψψ.

To include terms with derivatives in the Lagrangian we have to introduce the parallel transport, T (y, x), to compare two fields at different points of the space-time. The parallel transport in general is related to a curve C xy from x to y, as

T (y, x) = Pe ig Cxy dxµAµ(x)
where P denotes the path ordering and A µ which is a N × N matrix is usually called the connection and in physics the gauge field. The gauge field belongs to the adjoint representation of the Lie algebra of SU (N ) and can be written as A µ (x) = A a µ (x)t a . The gauge field will describe the gluon field. The parallel transport of ψ(x) gives a new field ψ(y) which is required to transform like ψ(y). That means that

ψ(y) = T (y, x)ψ(x) → U (y) ψ(y). Therefore, T (y, x) → U (y)T (y, x)U † (x)
and the gauge bosons under a gauge transformation transform as

A µ (x) = U A µ U † - i g (∂ µ U )U † (1.1)
which in a infinitesimal gauge transformation can be written as

A µ (x) ∼ A µ (x) -∂ µ θ a + gf abc A b µ θ c t a + O(θ 2 ).
Let us define the covariant derivative, D µ , in the direction of the unit vector n using the parallel transport as

n µ D µ ψ = lim →0 1 (ψ(x + n) -T (x + n, x)ψ(x)) .
A straightforward calculation shows that the covariant derivative takes the form

D µ = ∂ µ -igA a µ t a .
Over the spinor ψ it takes the form

(D µ ψ) b = (D µ ) bc ψ c = ∂ µ ψ b -igA a µ (t a r ) bc ψ c
where we have made explicit the color index b and c for the fermions. The generators t a in this case are in the fundamental representation denoted by the subindex r. If, on the other hand, we consider a field X in the adjoint representation (subindex G) the generators have to be in the same adjoint representation and therefore

(D µ X) b = (D µ ) bc X c = ∂ µ X b + gf bac A a µ X c
where we have used that (t a G ) bc = if bac . These derivatives are called covariant since they transform as the field itself, meaning that

D µ ψ(x) → U (x)D µ ψ(x).
Then, a good gauge and Lorentz invariant candidate is ψ(γ µ D µ )ψ which is a similar term of the one in the Dirac Lagrangian by changing the standard derivative to the covariant one. This term includes a kinetic term for the fermions and also an interaction term between fermions and gluons.

A gauge invariant term governing the dynamics of the gauge bosons can be obtained by observing that [D µ , D ν ]ψ(x) is also covariant and therefore

[D µ , D ν ]ψ(x) → U (x)[D µ , D ν ]ψ(x).
Using some algebra we can compute

[D µ , D ν ]ψ(x) = -ig ∂ µ A a ν -∂ ν A a µ + gf abc A b µ A c ν t a ψ(x) = -igF a µν t a ψ(x)
where we have defined

F a µν = ∂ µ A a ν -∂ ν A a µ + gf abc A b µ A c ν .
With this definition Tr F a µν t a F b µν t b = 1 2 F a µν F a µν is also gauge invariant. We could also include a term of the form Θ ab µνρσ F a µν F b ρσ , where Θ ab is a constant matrix and µνρσ is the Levi-Civita tensor, but this term is a total derivative and therefore it does not affect the equation of motion or the Feynman rules. For this reason it is not included in our Lagrangian.

To conclude, we can write the SU (N ) invariant Euclidean Lagrangian as

L inv = 1 4 F a µν F a µν + N f i=1 ψi (-γ µ D µ + M i )ψ i , (1.2)
where

F a µν = ∂ µ A a ν -∂ ν A a µ + gf abc A b µ A c ν , D µ ψ = ∂ µ ψ -igA a µ t a ψ.
The index i runs over the N f quark flavours. The coupling constant g represents the strength of the interaction between quarks and gluons as well as the gluon self-interaction.

The gauge invariance imposes us to describe both interactions with the same bare coupling constant. It is important to mention that the sign for each term on the Lagrangian was chosen so this Lagrangian matches through a Wick rotation with the Minkowski Lagrangian used by, e.g. [START_REF] Gracey | Three loop MS-bar renormalization of the Curci-Ferrari model and the dimension two BRST invariant composite operator in QCD[END_REF].

It is possible to consider other terms invariant under this gauge group, however, these terms would include a larger number of fields or derivatives. These other terms are known as non-renormalizable and they are associated with irrelevant operators, that means that at low energies with respect to a microscopical natural scale, the contribution of this terms can be ignored.

Quantization

In the previous section we have introduced the classical Lagrangian density for the SU (N )invariant renormalizable theory. In this section we would like to quantize the theory.

In order to quantize the theory there are a few standard methods commonly used in physics in the Minkowski space, see e.g. [START_REF] Muta | Foundations of quantum chromodynamics: An introduction to perturbative methods in gauge theories[END_REF]. We are going to mention only two of them here. The first one is the canonical quantization based on considering the fields as operators with some commutation rules. In this case, the Green's functions are computed as the expectation value of the operators acting in the vacuum as

0| T [ φ(x 1 )... φ(x n )] |0
where φ represents any possible field and φ the corresponding operator in Heisenberg picture, T is the time ordering, and |0 is the vacuum of the theory. The second equivalent method is the functional formalism based on Feynman integrals. It relies on the idea that a physical system going from state a to state b can take any possible path between those configurations. Every path contributes with an amplitude proportional to the exponential of iS M with S M the classical action in the Minkowski space. In this formalism the Green's functions can be written, using the functional integral, as

Dφ φ(x 1 )...φ(x n )e iS M Dφ e iS M
computing S M considering a slightly imaginary temporal axis. Both formalisms give the same Green's functions, therefore,

0| T [ φ(x 1 )... φ(x n )] |0 = Dφ φ(x 1 )...φ(x n )e iS M Dφ e iS M .
The right hand side integral can be computed in the Euclidean space through a Wick's rotation which consists in rotating the temporal coordinate to the imaginary axis by changing x 0 M to ix 0 E . This deformation of the contour transform the Minkowski path integral into an Euclidean one:

Dφ φ(x 1 )...φ(x n )e iS M Dφ e iS M = Dφ φ(x 1 )...φ(x n )e -S inv Dφ e -S inv
where S inv = d 4 xL inv is the Euclidean classical action. The last term can be thought as the expectation value of the product φ(x 1 )...φ(x n ) with a Boltzman-like weight e -S inv . The correlation functions that we are interested in calculating are

0| T [ φ(x 1 )... φ(x n )] |0 ≡ φ(x 1 )...φ(x n ) = Dφ φ(x 1 )...φ(x n )e -S inv Dφ e -S inv .
The functional or Feynman formalism [START_REF] Feynman | Space-time approach to nonrelativistic quantum mechanics[END_REF] shows the symmetries involved in the theory explicitly. That is an important reason to choose it when we want to study features related to the symmetries. Moreover, any correlation function can be easily obtained by derivation of the generating functional, Z[J] , which is analogue to the partition function defined as

Z[J] = Dφ e -S+ d 4 xJ i (x)φ i (x)
where φ i denotes all the field involved. The correlation function

φ in (x n )...φ i 1 (x 1 ) = 1 Z[J] Dφ φ in (x n )...φ i 1 (x 1 )e -S+ d 4 xJ i (x)φ i (x)
is obtained by successive derivatives of the generating function

φ in (x n )...φ i 1 (x 1 ) = 1 Z[J] δ n Z[J] δJ in (x n )...δJ i 1 (x 1 )
where the first derivative was taken with respect to J i 1 and so on. Generally it is convenient to define W [J] as Z[J] = e W [J] , W [J] is the generating function of the connected Green's function. That means that

φ in (x n )...φ i 1 (x 1 ) c = δ n W [J] δJ in (x n )...δJ i 1 (x 1 )
(perturbatively only connected diagrams contribute to the connected Green's functions). If we perform the Legendre transformation of W [J] we construct Γ[ φ ] called the effective action or the generating functional of the proper vertex, explicitly defined as

Γ[ φ ] = -W [J] + d 4 xJ i (x) φ i (x) .
From this definition it can be shown that

δΓ[φ] δφ i (z) = J i (z)
where

φ i = φ i = δW [J] δJ i .
The functional Γ is also called the one particle irreducible (1-PI) generating functional, because the n-derivative of Γ, Γ (n) can be represented by the 1-PI Feynman diagrams with n external legs. We say that a connected diagram is 1-PI if, when cutting any internal line, the diagram remains connected. The second and third derivative of the vertex functional are related to the two and three-point correlation functions as it is shown in the Appendix A.1. This will be very useful for our calculations.

Functional integral quantization

In this formalism it is easy to observe the problems arising from the gauge symmetry. For instance, the bare propagator of the gluon is ill defined. Remember that the propagator is the two-point connected correlation function using only the quadratic part of the Lagrangian, L 0 , and it can be calculated as the inverse of the second derivative of the action as

φ(x 1 )φ(x 2 ) 0 = Dφ φ(x 1 )φ(x 2 )e -d 4 xL 0 Dφ e -d 4 xL 0 = δ 2 S δφ(x 1 )δφ(x 2 ) -1 .
The term of the quadratic Lagrangian involved in the gluon propagator is

1 4 ∂ µ A a ν -∂ ν A a µ 2
. The second derivative in the Fourier space takes the form p 2 δ µν -p µ p ν (see Appendix A.2.2) which is degenerate and is thus not invertible. To define the gluon propagator we have to fix the gauge, by imposing a gauge condition.

With this purpose we are going to define what is called a gauge orbit. A gauge orbit consists of iterations of the gauge transformation, e.g. Eq.(1.1), of a chosen configuration of the field. Let us write φ = φ θ , where φ represents an arbitrary field, in order to make explicit the dependence on the gauge parameter θ. The gauge orbit of φ(x) is the set {φ θ } for all θ. Along the orbit the physics of the systems remains the same due to the gauge invariance of the action and of the physical observables. We want to choose only one representative of the gauge orbit, φ θ, for example the one satisfying a certain gauge condition f [φ θ, x] = 0.

Faddeev and Popov in [START_REF] Faddeev | Feynman Diagrams for the Yang-Mills Field[END_REF] proposed a mechanism for fixing the gauge in gauge theories which can be expressed in a Lagrangian formalism. The Faddeev-Popov procedure is really useful to determine the gauge-fixed Lagrangian which allow to compute with success several quantities in QCD. However, it relies on the hypothesis of the existence of a unique θ which satisfies the gauge condition. Constraints with that property are called ideal gauge conditions [START_REF] Williams | Nonperturbative QCD, gauge fixing, Gribov copies, and the lattice[END_REF]. However, there are some gauges for which this hypothesis is not true. That means that there exists different φ θ in the same orbit, usually called Gribov copies, that satisfy the same gauge condition. In this case, the Faddeev-Popov procedure is not well justified. We will come back to this issue in chapter 2.

Faddeev and Popov procedure

As explained in the above section, in order to be able to compute the Green's function we would like to choose in the integral one representative of the gauge orbit satisfying the gauge condition, f [φ θ , x] = 0.

For that purpose (see e.g. [START_REF] Weinberg | The quantum theory of fields[END_REF]) Faddeev and Popov proposed to consider the integral

J = Dφ Y[φ]B[f [φ]]DetF[φ]
where

Dφ will represent DψD ψDA µ , Y[φ] is a gauge invariant functional, B[f [φ]
] is a weight that impose the gauge condition and

F ax,by = δf a [φ θ , x] δθ b (y) θ=0 .
If we consider the same integral J but we use instead as integration variable the transform of φ with parameter Θ, φ Θ we have

J = Dφ Θ Y[φ Θ ]B[f [φ Θ ]]DetF[φ Θ ]
and assuming that the measure and Y are gauge invariant we obtain

J = Dφ Y[φ]B[f [φ Θ ]]DetF[φ Θ ].
We can integrate now with respect to Θ with an auxiliary weight ρ[Θ] that will be characterised latter. We can observe that

J x,a dΘ a (x)ρ[Θ] = Dφ Y[φ]C[φ]
where

C[φ] = x,a dΘ a (x)ρ[Θ]B[f [φ Θ ]]DetF[φ Θ ].
We would want to show that C[φ] actually does not depend on φ. Let us define Θ to be the parameter associated with the product of the gauge transformation with parameter θ and Θ. If we apply the chain rule to F[φ Θ ], we obtain

F ax,by [φ Θ ] = δf a [(φ Θ ) θ , x] δθ b (y) θ=0 = d 4 z δf a [φ Θ, x] δ Θc (z) Θ=Θ δ Θc (z) δθ b (y) θ=0 (1.3)
where we are going to denote

δf a [φ Θ, x] δ Θc (z) Θ=Θ = δf a [φ Θ , x] δΘ c (z) = J ax,cz [φ, Θ] (1.4) and δ Θc (z, Θ, θ) δθ b (y) θ=0 = R cz,by .
Therefore, we can write

DetF[φ Θ ] = DetJ[φ, Θ] DetR[Θ].
Now, it is convenient to choose ρ[Θ] = 1/DetR[Θ] and therefore when these expression are introduced in the definition of C[φ] we obtain

C[φ] = x,a dΘ a (x)B[f [φ Θ ]]DetJ[φ, Θ]. (1.5)
If we change the variable integral from Θ to f , DetJ[φ, Θ] acts as the Jacobian determinant and the integral can be simplified to

C[φ] = a,x df a (x)B[f ] = C.
We can conclude that

Dφ Y[φ] = 1 C J x,a dΘ a (x)ρ[Θ] (1.6)
If we want to compute the expectation value of a gauge invariant operator O inv , we can use the preceding result taking Y = O inv e -S inv and Y 0 = e -S inv to write

O inv = Dφ O inv e -S inv Dφ e -S inv = Dφ Y[φ] Dφ Y 0 [φ] = Dφ Y[φ]B[f ]DetF[φ] Dφ Y 0 [φ]B[f ]DetF[φ] (1.7)
For simplicity, we can choose

B[f ] = e -1 2ξ d 4 xfa(x)fa(x)
where ξ is an arbitrary real parameter. Moreover, we can consider the simplest covariant gauge condition:

f a [A µ , x] = ∂ µ A a µ (x)
. Therefore, the only effect of B[f ] is to modify the action by a gauge-fixing term S GF =

d 4 x 1 2ξ ∂ µ A a µ 2 .
In consequence, including this result in Eq.(1.7), we obtain

O inv = Dφ Y[φ]e -S GF DetF[φ] Dφ Y 0 [φ]e -S GF DetF[φ] = Dφ O inv e -(S inv +S GF ) DetF[φ] Dφ e -(S inv +S GF ) DetF[φ] .
This means that gauge invariant operator can be calculated with or without gauge fixing obtaining the same result. This can be verified in lattice simulations where both calculations are well defined. However, as we mentioned earlier the integrals without fixing the gauge in the continuum carry some problems due to the gauge invariant property and they are not properly defined. Therefore, in the continuum it is necessary to compute the expected values fixing the gauge.

The ghost invasion

With the help of Grassman variables (see e.g. [START_REF] Peskin | An Introduction to quantum field theory[END_REF]) we can express the determinant of an arbitrary matrix M as detF = DcDc exp -d 4 xd 4 y c(x)F(x, y)c(y)

(1.8)

We will use this trick to write DetF[φ] explicitly . The auxiliary Grassman field, c and c are usually called ghosts because they are not physical fields.

For the covariant gauge condition ∂ µ A θ,a µ (x) the matrix is

F ax,by = δ(∂ x µ A θ,a µ (x)) δθ b (y) θ=0 = ∂ x µ δA θ,a µ (x) δθ b (y) θ=0 = -∂ x µ (∂ µ δ ab + gf acb A c µ )δ 4 (x -y). (1.9)
which is called the Faddeev-Popov operator. In summary the factor DetF[φ] can be expressed as

DetF[φ] = DcDc exp -d 4 x ∂ µ ca (∂ µ δ ab + gf acb A c µ )c b .
(1.10)

The final effective Lagrangian density in a linear covariant gauge is

L = L inv + L GF + L FP = 1 4 F a µν F a µν + 1 2ξ (∂ µ A a µ ) 2 + ∂ µ ca (∂ µ δ ab + gf acb A c µ )c b (1.11)
where the first term is the gauge invariant Lagrangian, the second term arises from fixing the gauge and the third one is the Faddeev-Popov ghost Lagrangian. The ghost lives in the adjoint representation. Accordingly the covariant derivative over c takes the form, (D µ c) a = ∂ µ c a + gf abc A b µ c c . The Faddeev-Popov procedure remains untouched for linear covariant gauges with the addition of the quarks. To conclude, the effective Lagrangian might be written as

L = 1 4 F a µν F a µν + N f i=1 ψi (-γ µ D µ + M i )ψ i + ∂ µ ca (D µ c) a + 1 2ξ (∂ µ A µ ) 2 .
(1.12)

Landau gauge

There is still a little modification of the specific form of the Lagrangian that we are going to use. In this thesis the gauge fixing condition chosen is the Landau gauge. Landau gauge condition consists in constraining the gauge field to the ones fulfilling exactly ∂ µ A a µ = 0 which is equivalent to imposing ξ = 0.

The term of the functional integral corresponding to L

GF = 1 2ξ ∂ µ A a µ 2
can be expressed in a more convenient form by using an auxiliary bosonic field h a . For that purpose let us note that the Gaussian integral

Dh a e -d 4 x ξ 2 h a + i √ 2ξ ∂µA a µ 2 = const
gives a constant value which can be omitted in the generating functional. Thus, we can do the following replacement e

-d 4 x 1 2ξ (∂µA a µ ) 2 ∝ e -d 4 x 1 2ξ (∂µA a µ ) 2 Dh a e -d 4 x ξ 2 h a + i √ 2ξ ∂µA a µ 2 = Dh a e -d 4 x( ξ 2 (h a ) 2 +ih a ∂µA a µ ) . (1.13)
That is, one can add the auxiliary field h a and replace the term of the action

d 4 x 1 2ξ ∂ µ A a µ 2 by d 4 x ξ 2 (h a ) 2 + ih a ∂ µ A a µ .
For now on we work in Landau gauge (ξ → 0). In consequence, the Lagrangian we are going to work with is

L = 1 4 F a µν F a µν + N f i=1 ψi (-γ µ D µ + M i )ψ i + ∂ µ ca (D µ c) a + ih a ∂ µ A a µ .
(1.14)

The Feynman rules associated with this Lagrangian are presented in section 1.3.

BRST symmetry

After fixing the gauge the action is no longer gauge invariant. However, it still has a symmetry called BRST symmetry. The BRST symmetry was discovered by Becchi, Rouet and Stora [START_REF] Becchi | Renormalization of the Abelian Higgs-Kibble Model[END_REF][START_REF] Becchi | Renormalization of Gauge Theories[END_REF] and independently by Tyutin [Tyu75]. This symmetry acts on the fields introducing variations of the form:

δ s Φ = η sΦ
where η is an x-independent Grassmann number and s is the BRST operator that acts over the fields as:

sA a µ = (D µ c) a , sψ = -igt a c a ψ, sc a = - g 2 f abc c b c c , sc a = ih a and sh a = 0.
Moreover for any product of two set of fields Φ 1 and Φ 2 the BRST operator acts as

s (Φ 1 Φ 2 ) = (sΦ 1 )Φ 2 ± Φ 1 sΦ 2
where the ± sign is plus if Φ 1 is bosonic and minus if Φ 1 is fermionic. Knowing this, one of the most useful properties of BRST symmetry can be demonstrated: BRST symmetry is nilpotent meaning that s 2 = 0, see e.g. [START_REF] Weinberg | The quantum theory of fields[END_REF].

In order to verify that the above transformation is a symmetry of the action it is convenient to note that the BRST transformation for the matter and gauge field is like an infinitesimal gauge transformation with gauge parameter ηc. Therefore, the gauge invariant terms of L inv are also BRST invariant. On the other hand, it can be shown that the remaining terms of the action corresponding to d 4 x (L GF + L F P ) are of the form sS(A, ψ, ψ, c, c, h). In the case of Landau gauge, we can write

d 4 x (L GF + L F P ) = d 4 x ih a ∂ µ A a µ + ∂ µ ca (D µ c) a = s d 4 xc a ∂ µ A a µ .
More generally, see [START_REF] Weinberg | The quantum theory of fields[END_REF], the BRST-invariant action takes the form

S = S inv + sS. (1.15)
The first term is BRST invariant because it is gauge invariant and the second term is also BRST invariant because s is nilpotent. We can say that S inv ∈ Ker(s) and sS ∈ Im(s).

Fixing the gauge with the Faddeev-Popov procedure brings new fields which are not expected to have physical meaning. Hence, we need to clarify which is the Hilbert space containing only physical states. To define physical states we are going to consider the asymptotic fields, long-time before or after the interaction. Physical states should have positive norm and do not depend on the gauge. That means that any change in S, ∆S, should not change the matrix element α|β for any pair of physical states, |α and |β .

To define the physical space let us introduce the BRST charge. As the BRST is a global symmetry it has a conserved current. This current allows us to define an associated charge Q s which is the generator of the symmetry. The variation of an arbitrary set of fields Φ can be expressed as:

i sΦ = [Q s , Φ] ∓
where the sign beingor + depending on whether Φ is bosonic or fermionic. For physical states we expect that

α| s∆S |β = 0 This implies, α| [Q s , ∆S] |β = 0.
As this happens for any possible change ∆S we can conclude that one must require α| Q s = Q s |β = 0, both states belong to the kernel of Q s . Two states differing by a state in the image of Q s give the same S-matrix element (because Q s is also nilpotent), so those states are equivalents. That is the reason why the physical space, V phys , is taken as the cohomology of the BRST symmetry,

V phys = Ker(Q s ) Im(Q s ) .
It can be shown, see [START_REF] Weinberg | The quantum theory of fields[END_REF], that only gluons with transverse polarization belong to the cohomology of the BRST symmetry. The antighost does not belong to the KerQ s and the ghost belongs to Im(Q s ). Therefore none of them belong to V phys .

To conclude that Ker(Qs) /Im(Qs) is a good definition for the physical space it is necessary to prove that the S-matrix, S, restricted to V phys is unitary. The S-matrix is pseudo-unitary for an action of the form (1.15) that means that S † S = SS † = I but we must assure that the space is a Hilbert space with positive norm. In fact, we build V phys in order to ensure that any matrix element between two physical states are gauge independent. Therefore, we can study the S-matrix in a more physical gauge like Coulomb or axial gauge where there is no problem of positivity of the states, see [START_REF] Weinberg | The quantum theory of fields[END_REF]. Therefore in this gauges the restriction of the S-matrix to the physical space V phys is unitary and, as in this space the elements of the S-matrix do not depend on the gauge, we can conclude that the S-matrix is unitary in V phys .

Another important feature of BRST symmetry is that it allows us to prove that the theory is renormalizable (concept introduced in Sect.1.4) to all orders of perturbation theory [START_REF] Hooft | Renormalizable Lagrangians for Massive Yang-Mills Fields[END_REF][START_REF] Hooft | Renormalization of Massless Yang-Mills Fields[END_REF][START_REF] Gerard 't Hooft | Regularization and Renormalization of Gauge Fields[END_REF].

Feynman rules

Let us return to the computation of Green's functions. We are going to use the gauge-fixed Lagrangian presented in Eq.(1.14). The Green's function of n-arbitrary fields, φ, in the functional formalism can be expressed as

G(x 1 , ..., x n ) = φ(x 1 )...φ(x n ) = DADhDcDcD ψDψ φ(x 1 )...φ(x n )e -d 4 xL
DADhDcDcD ψDψ e -d 4 xL

Obviously the problem is not completely solved: those integrals can not be easily computed analytically. However, perturbation theory and Feynman diagrams allows us to manipulate them in a clever way.

To use perturbation theory it is convenient to separate the Lagrangian into the free or quadratic Langrangian, L 0 , and the interacting Lagrangian, L int .

L = L 0 + L int where L 0 = 1 4 (∂ µ A a ν -∂ ν A a µ ) 2 + N f i=1 ψi (-γ µ ∂ µ + M i )ψ i + ∂ µ ca ∂ µ c a + ih a ∂ µ A a µ (1.16)
and

L int = 1 2 gf abc (∂ µ A a ν -∂ ν A a µ )A b µ A c ν + 1 4 g 2 (f abc A b µ A c ν ) 2 + ig N f i=1 ψi γ µ A a µ t a ψ i + gf abc ∂ µ ca A b µ c c .
(1.17) If the coupling constant, g, is small enough we can expand the exponential

e -d 4 xL ∼ e -d 4 xL 0 1 -d 4 xL int + 1 2 d 4 xL int 2 - 1 3! d 4 xL int 3 + . . . .
The advantage is that every term in the expansion can be seen as an expectation value with a Gaussian weight. Using Wick theorem, see e.g. [START_REF] Peskin | An Introduction to quantum field theory[END_REF], the Gaussian expectation value can be computed easily as the sum of all possible pairwise contractions. That means,

φ(x 1 )...φ(x n ) 0 = φ(x 1 )φ(x 2 ) 0 ... φ(x n-1 )φ(x n ) 0 + φ(x 1 )φ(x 3 ) 0 ... φ(x n-1 )φ(x n ) 0 + all other product of pairwise Gaussian expectation values (1.18)
where the subindex 0 denotes the expectation value with L 0 . This property allows us to describe the correlation functions with the Feynman diagrams. The Wick theorem is a very useful tool since the Gaussian expectation value of two fields is really easy to compute, it takes the value

D 0 (x 1 -x 2 ) = φ(x 1 )φ(x 2 ) 0 = Dφ φ(x 1 )φ(x 2 )e -d 4 xL 0 Dφ e -d 4 xL 0 = δ 2 S δφ(x 1 )δφ(x 2 ) -1
. The Green's function φ(x 1 )φ(x 2 ) 0 is called the free propagator and it will be represented in Feynman diagrams with a line (full, dashed, curled, etc depending on the field) from x 2 to x 1 .

Chapter 1. Quantum Chromodynamics

The interacting Lagrangian which will define the interaction vertex of the theory can be written in a symmetric way as follows

d d xL int = d d x 1 2 gf abc (∂ µ A a ν -∂ ν A a µ )A b µ A c ν + 1 4 g 2 (f abc A b µ A c ν ) 2 + gf abc ∂ µ ca c b A c µ = d d xd d yd d z 1 3! A a µ (x)A b ν (y)A c ρ (z)S abc µνρ (x, y, z) + d d xd d yd d zd d w 1 4! A a µ (x)A b ν (y)A c ρ (z)A d σ (w)S abcd µνρσ (x, y, z, w) + d d xd d yd d zA c µ (z)c b (y)c a (x)S abc µ (x, y, z) + d d xd d yd d zA c µ (z) ψb (y)S abc µ,ψ (x, y, z)ψ a (x) (1.19)
where

S abc µνρ (x, y, z) = δ 3 S δA a µ (x)δA b ν (y)δA c ρ (z) A,c,c=0 , S abcd µνρσ (x, y, z, w) = δ 4 S δA a µ (x)δA b ν (y)δA c ρ (z)δA d σ (w) A,c,c=0 , S abc µ (x, y, z) = δ 3 S δc a (x)δc b (y)δA c µ (z) A,c,c=0
and

S abc µ,ψ (x, y, z) = δ 3 S δψ a (x)δ ψb (y)δA c µ (z) A,ψ, ψ=0
.

All these derivative are presented in the Appendix A.2. The Fourier transform of, for instance, S abc µ (x, y, z), is given by

d d xd d yd z e -ipx e -iqy e -ikz S abc µ (x, y, z) = (2π) d δ(p + q + k)S abc µ (p, q, k)
and in a similar way for the others (all momenta are incoming to the vertex). The inverse Fourier transform is defined as

S abc µ (x, y, z) = d d p (2π) d d d q (2π) d d d k (2π) d (2π) d δ(p + q + k)S abc µ (p, q, k) e ipx e iqy e ikz .
At the lowest order in the perturbation theory (tree level) the Legendre transform of the effective action, Γ, has the same form as the action, that means that

Γ tree [ φ ] = S[ φ ]
Therefore, Γ (3),tree-level and Γ (4),tree-level will represent the third and fourth derivative of the action, usually called the three-gluon and four-gluon vertex. The number of Lorentz indices reflects the number of derivatives with respect to gauge fields.

To summarise the Feynman diagrams are a diagrammatic representation of the expression to be computed. Each diagram represents an integral in the perturbative expansion. They consist in putting together: external lines (associated with φ(x 1 ), ..., φ(x n )), propagators (internal lines) and vertices; obeying the following rules (in momenta space, with all momenta incoming) of table 1.1.

Each diagram comes with a symmetry factor, which can be calculated as follows:

• There is a minus for each fermion loop. (3),tree-level

ψ ψA a µ = igt a γ µ r p k µ, a three-gluon vertex Γ (3),tree-level abc µνρ (p, q, r) = igf abc [(q -r) µ δ νρ +(r -p) ν δ µρ + (p -q) ρ δ µν ] q r p µ, a ν, b ρ, c
four-gluon vertex Γ (4),tree-level abcd µνρσ = g 2 f eab f ecd (δ µρ δ νσ -δ µσ δ νρ )

+ f eac f ebd (δ µν δ ρσ -δ µσ δ νρ ) +f ead f ecb (δ µρ δ νσ -δ µν δ σρ ) µ, a ν, b ρ, c σ, d
• There is a (-1) v /v! (where v is the number of vertex) that come from the expansion of the exponential.

• There is a 1/3! for each three gluon vertex.

• There is a 1/4! for each four gluon vertex.

• There is an extra factor for the diagrams that have different kind of vertex that come from expanding (L int ) v . We will have a factor v!/(v 1 !v 2 !v 3 !) where v 1 is the number of the vertex of one type and v 2 for the other type, and v 3 for the last one.

• Finally, we have to add the usual combinatorial factor due to the possible contractions coming from the Wick theorem.

The extraction of Feynman rules was done using the following convention of the Fourier transformation:

f (p) = e -ipx f (x)d 4 x.
We want to recall that the second derivative of Γ (2) is the inverse of the propagator, see Appendix A.1. While, the three-point vertex, Γ (3) , is related to the expectation value of the mentioned involved fields, G(p 1 , p 2 , p 3 ; g, M ), by

G(p 1 , p 2 , p 3 ; g, M ) = -D(p 1 )D(p 2 )D(p 3 )Γ (3) (p 1 , p 2 , p 3 ; g, M )
where D(p i ) are the propagators of the external fields.

For a detailed deduction of the Feynman rules from the action see Appendix A.2.

Regularization and Renormalization

When we want to explicitly calculate the Feynman diagrams some of the integrals present ultra-violet divergences. To avoid these infinities, a regularization method is commonly used. The most popular ones are: the inclusion of an ultra-violet cut-off, Λ, as the upper bound of the integrals; Pauli-Villars regularization; lattice regularization and the dimensional regularization (see [START_REF] Bollini | Lowest order divergent graphs in nu-dimensional space[END_REF][START_REF] Bollini | Dimensional Renormalization: The Number of Dimensions as a Regularizing Parameter[END_REF][START_REF] Gerard 't Hooft | Regularization and Renormalization of Gauge Fields[END_REF]). The last two have the advantage of respecting gauge invariance. In the analytical calculation of Feynman diagrams the dimensional regularization is the most convenient regularization to use and is the one that we will use in chapters 3,4,5 and 6. The dimensional regularization consists in computing the integrals appearing in the Feynman diagrams analytically as a function of the space dimension, d. Usually the integral in four dimension momentum space

d 4 q (2π) 4 is replaced by d d q (2π) d .
This behaves as an ultraviolet regulator because it eliminates some ultraviolet divergent integrals. For instance, the integral

d d q (2π) d 1
(q 2 +m 2 ) n can be computed as

d d q (2π) d 1 (q 2 + m 2 ) n = Γ(n -d/2) (4π) d/2 Γ(n) (m 2 ) d/2-n
and therefore the integral

d d q (2π) d 1 q 2 ,
which is a priori UV-divergent in four dimensions, is set to zero in the dimensional regularization.

Even though the cut-off regularization is not gauge invariant and therefore it will not be used in our calculations, for the introduction of the renormalization concept it is better to work with the cut-off regularization, that means that the integrals are replaced by

∞ dq → Λ dq
We would like to redefine the field, masses, and coupling constant in order to absorb the dependence on the cut-off. The new field and couplings are said to be renormalized and the ones appearing at the beginning in the Lagrangian are called bare quantities.

A a µ B = Z A A a µ , c a B = Z c c a , ca B = Z c ca , ψ a B = Z ψ ψ a , g B = Z g g and M B = Z M M (1.20)
The index B represents the bare quantities. In this section, all quantities are assumed to be renormalized unless the subindex B is present.

With the new normalization the bare Green's function are related to their corresponding renormalized one, i.e. the expected value of the renormalized fields, through the product of Z-factors of the intervening fields, Z G ,

G R (p 1 , .., p n ; g, M ) = Z G G B (p 1 , .., p n ; g B , M B )
The renormalized quantities are defined by choosing a renormalization scheme. That means, in a given momentum scale µ we choose the value of some renormalized Green's functions in order to define the renormalization factors. For a renormalizable theory it is enough to fix the Green's functions appearing at tree level, i.e. the ones that explicitly appear in the Lagrangian.

However, ignoring the gauge-fixing term, there are seven different Green's functions with tree-level contribution; the inverse of gluon, ghost and quark propagators and the three-gluon, four-gluon, ghost-gluon and quark-gluon vertex, and only five Z-factors to define. Fortunately, thanks to BRST invariance it is possible to find some identities, called Slavnov-Taylor identities, which relate the Green's functions and allow to show that the renormalization can be successfully done.

There are different renormalization schemes, all of which absorb the divergence in the Z-factors but they differ in the modification of the finite part. The most common schemes are the minimal subtraction schemes MS and MS. The scheme that we used in our calculations is presented in section 3.2.1.

Let us continue with an example of the procedure to obtain the Z-factor of a given renormalization scheme. The bare ghost propagator computed from the Feynman digramas has the form

D ab B (p) = δ ab J B (p) p 2
where J B (p) is the bare dressing function of the ghost propagator. The renormalized ghost propagator is obtained by including the Z c and expressing all the masses and coupling constant in term of the renormalized ones, thus

D ab (p) = Z c δ ab J B (p) p 2 = δ ab J(p) p 2
where J(p) is called the ghost dressing function. In the scheme that we are going to use, a particular case of the MOM scheme, Z c is defined by requiring that the propagator acquires its tree-level form at some momentum scale µ. That means,

D ab (p = µ) = δ ab 1 µ 2 .
This allow us to obtain the Z c . Obviously Z c depends on the scale µ where the renormalization scheme was defined, and so does the propagator.

Renormalization group

In this section we are going to introduce the renormalization group equation in order to study the dependence on the renormalization scale of the correlation function. When we renormalize the theory at the scale µ, generally in four dimensions, the calculation of the correlation function brings logarithms of the form log( µ p ) where p is the typical momentum of the function studied. In order to avoid large logarithms the renormalization scale µ should be chosen near p. The problem is that we are usually interested in studying a large spectra of momenta. The idea is to think µ as a sliding scale and compare how the evolution of µ affects the correlation functions.

For the derivation of the renormalization group equation it is convenient to note that the bare correlation functions does not depend on the scale µ, that is

0 = µ d dµ G B (p, g B , {M i } B ).
On the other hand, when we express the bare function in terms of the renormalized correlation function the last one depends explicitly and implicitly on the scale µ. It depends implicitly on the renormalization scale through the dependence of the renormalized coupling constant, masses and field. The dependence of these quantities is determine by the β-functions and anomalous dimensions defined as:

β g (g, {M i }) = µ dg dµ g B ,{M i } B , γ A (g, {M i }) = µ d log Z A dµ g B ,{M i } B , γ c (g, {M i }) = µ d log Z c dµ g B ,{M i } B , β M i (g, {M i }) = µ dM i dµ g B ,{M i } B , γ ψ (g, {M i }) = µ d log Z ψ dµ g B ,{M i } B . (1.21)
The evolution of g(µ) and M i (µ) are obtained by integration of the beta functions with initial conditions given at some scale µ 0 . That gives origin to the renormalization group equation:

  µ ∂ ∂µ + 1 2 nγ + β g ∂ ∂g + N f i β M i ∂ ∂M i   G(p, µ, g(µ), {M i (µ)}) = 0. (1.22)
where n is the number of fields appearing in the correlation function G. For convenience we are going to write down the renormalization group equation for the vertex function with n A gluon legs, n c ghost legs and n ψ quark legs explicitly:

  µ∂ µ - 1 2 (n A γ A + n c γ c + n ψ γ ψ ) + β g ∂ ∂g + N f i β Mi ∂ ∂M i   Γ (n A ,nc,n ψ ) (p, µ, g(µ){M i (µ)}) = 0, (1.23)
The solution of the renormalization group equation (1.23)has the form

Γ (n A ,nc,n ψ ) (p, µ, g(µ), {M i (µ)}) = z A (µ) n A /2 z c (µ) nc/2 z ψ (µ) n ψ /2 Γ (n A ,nc,n ψ ) (p, µ 0 , g(µ 0 ), {M i (µ 0 )}). (1.24)
where

log z A (µ) = µ µ 0 dµ µ γ A g(µ ), {M i (µ )} , log z c (µ) = µ µ 0 dµ µ γ c g(µ ), {M i (µ )} , log z ψ (µ) = µ µ 0 dµ µ γ ψ g(µ ), {M i (µ )} .
(1.25)

This last equation relates the vertex function renormalized at the scale µ 0 , Γ (n A ,nc,n ψ ) (p, µ 0 , g(µ 0 ), {M i (µ 0 )}), to the vertex function Γ (n A ,nc,n ψ ) (p, µ, g(µ), {M i (µ)}) renormalized at the scale µ. In this thesis we are going to compute some particular cases of Γ (n A ,nc,n ψ ) (p, µ 0 , g(µ 0 ), {M i (µ 0 )}) when we choose µ 0 = 1GeV. In order to avoid the large logarithms we use the relation (1.24) to relate it to Γ (n A ,nc,n ψ ) (p, µ, g(µ), {M i (µ)}) which ensure us avoiding large logarithms if µ ∼ p.

Asymptotic freedom

The evolution of the coupling constant is of great importance when studying the nonabelian gauge theories. For instance, we know that when the coupling constant is small enough we are allowed to do perturbation theory. In that regime, even though the calculations could be tedious there is a well known strategy to follow.

The calculations of the β g -function in QCD was first done by Politzer, Gross and Wilczek in 1973, see [START_REF] Politzer | Reliable Perturbative Results for Strong Interactions?[END_REF][START_REF] Gross | Ultraviolet Behavior of Nonabelian Gauge Theories[END_REF]. It can be computed from the ghost-gluon, quarkgluon, three-gluon or four-gluon vertex. The first mentioned is the easiest to compute, and shows no divergence in four-dimensions in Landau gauge. In that gauge Z g can be easily obtained from Z A and Z c of the gluon and ghost propagators. The relation in Landau gauge is given by Taylor's theorem [START_REF] Taylor | Ward Identities and Charge Renormalization of the Yang-Mills Field[END_REF]:

Z g Z A Z c = 1
Taylor's theorem (see Appendix B.5.1) is based on the fact that, except for the treelevel expression, the contribution for the ghost-gluon vertex is proportional to the ghost momentum to all orders of the perturbation theory. Therefore, the ghost-gluon vertex at vanishing ghost momentum, Γ abc µ,B (p = 0, r), is equal to the tree-level vertex form, g B f abc r µ (with r µ the anti-ghost momentum). In the MS scheme, the renormalized quantities satisfy the relation

Z A Z c Γ abc µ (p = 0) = g Z g f abc r µ .
Taylor's scheme generalize this property also to the finite parts even without taking the M S scheme. Once, the renormalized coupling constant is defined we can compute the perturbative β-function at µ M , for M the mass of the quarks, obtaining [START_REF] Caswell | Asymptotic Behavior of Nonabelian Gauge Theories to Two Loop Order[END_REF][START_REF] Jones | Two Loop Diagrams in Yang-Mills Theory[END_REF]:

β g (g(µ)) = -β 0 g(µ) 3 16π 2 -β 1 g(µ) 5 256π 4 + O(g 7
) were the first two coefficients in the expansion, β 0 and β 1 , are renormalization-scheme independent and have the expressions.

β 0 = 11 3 C a - 4 3 T f N f , and β 1 = 34 3 C 2 a -4C f T f N f - 20 3 C a T f N f where C a = N ,C f = N 2 -1

2N

and T f = 1/2. For SU (3), as the number of quark flavours does not exceed six, the β-function is negative. This sign is really important because it means that the coupling constant falls down at high momentum. On the contrary, it seems to blow up at low momentum. This fact forbids the use of standard perturbation theory at low momenta, where others techniques need to be implemented.

Chapter 2

Infrared Regime of QCD At very high energies, QCD is asymptotically free, meaning that quarks and gluons behave like weakly interacting particles. On the other hand, at low energies quarks and gluons strongly interact. In fact, when we try to separate a pair quark-antiquark the gluon field forms a narrow tube with energy proportional to the distance. As the quarks are taken further apart the energy stocked in this tube increases until nature prefers to use the energy to create a new quark-antiquark pair instead of isolate the quarks. In fact, the way to study the strong force is through the bound states of quarks and gluons, the hadrons, as no free quark or gluon have ever been detected. This phenomenon is called confinement. Why confinement takes place and how it can be explained from the QCD Lagrangian is one of the most interesting questions of modern physics since there is no fully-understood analytical demonstration of confinement from first principles.

Although some progress in the subject have been made we are not going to treat the problem of confinement in this thesis. However, we are going to focus on a previous analysis trying to understand the problems and properties of the infrared regime. In this chapter, we would like to mention the problems which appear in the infrared regime, and the relevant features of the different methods used to attack them.

Gribov region

In chapter 1 we introduced the classical Lagrangian density invariant under the SU (N )gauge group. In order to work with a quantized theory we presented the functional Feynman formalism. The quantization of the theory implies to define the propagators of the involved fields. In order to properly define the propagators of a gauge-invariant theory it is necessary to fix the gauge. The standard procedure to fix the gauge is the Faddeev-Popov procedure which was presented in chapter 1. However, for instance in covariant gauges (Lorentz invariant gauges), the Faddeev-Popov procedure is not fully justified. This quantization assumes that every gauge orbit intersect the gauge condition only once. However, this assumption is not always true. For simplicity we are going to focus on a particular covariant gauge, the Landau gauge, which satisfies ∂ µ A µ = 0.

If we go back to the FP procedure there are some observations that we have ignored. Faddeev-Popov quantization is based on performing the change of variable Θ → f [φ Θ ] in the Eq. (1.5) with the goal of proving that C[φ] is constant. Nevertheless, this can be done if the DetJ defined in Eq. (1.4) is positive and f [φ] invertible. The procedure is not justified if there are more than one field configuration satisfying the same gauge condition. The field configurations on the same gauge orbit which satisfy the same gauge condition are known as Gribov copies [START_REF] Gribov | Quantization of Nonabelian Gauge Theories[END_REF][START_REF] Singer | Some Remarks on the Gribov Ambiguity[END_REF].

We are going to study the condition for a gauge orbit to intersect the Landau gauge condition twice, for simplicity we are going to focus on the gauge field, φ = A µ . Let us consider A µ and A θ µ fulfilling ∂ µ A µ = 0 and ∂ µ A θ µ = 0 both in the same gauge orbit. That means,

A θ µ = U A µ U † - i g (∂ µ U )U †
with U = e -igθ(x) an element of the SU (N )-group where θ = θ(x) represents θ a (x)t a . Assuming that the gauge fields are close enough we can approximate U ∼ 1 -igθ(x), then the relation becomes:

A θ µ = A µ -D µ θ where D µ θ = (D µ θ) a t a and D µ is the covariant derivative in the adjoint representation (D µ θ) a = ∂ µ θ a + gf abc A b µ θ c .
Both gauge fields can satisfy Landau gauge condition if and only if there is a solution to the equation

-∂ µ D µ θ = 0. (2.1)
The term in the left hand side of this equation is the FP operator already introduced in section 1.2.3. It is worth mentioning that in the perturbative regime (A µ ∼ 0) as well as in the abelian theory (f abc = 0) the gauge field does not present Gribov copies in Landau gauge. The reason is that in both cases the previous equation simplifies to -∂ 2 θ = 0. The Laplace equation does not have any solution except for the trivial one, with the border condition that θ(x) vanishes at infinity. However, Gribov showed that eq. (2.1) has non trivial solutions for SU (N ) and that means that for some values of the field A the Faddeev-Popov operator has zero modes.

To find a well defined gauge fixing procedure one option is to restrict the integration domain of the functional integral to a region containing only one representative of each gauge orbit. This aim can be divided in two difficult tasks. The first one is to identify the mentioned region and the second is to find a clever way to constrain the functional integral to it.

Eq. (2.1) gives a hint of how to determine the region without Gribov copies. It is natural to think that this region should avoid the zone in which FP operator has zero modes.

In fact, Gribov in [START_REF] Gribov | Quantization of Nonabelian Gauge Theories[END_REF] suggested to consider the so-called first Gribov region, Ω.

Ω = A a µ , ∂ µ A a µ = 0 and F positive definite
where the matrix F was introduced in Eq. (1.9) and reads

F ax,by = -∂ µ D ac µ δ(x -y) = -∂ 2,x µ δ ac + ∂ x µ f abc A b µ δ(x -y).
For the case when ∂ µ A µ = 0, F is hermitian. Therefore its eigenvalues are real, and the condition of being in the Gribov region implies that all the eigenvalues of F are positive.

At first sight the Gribov region seems to be the good choice to restrict the domain of the functional integrals. However, even though the Gribov region does not have Gribov copies which are infinitesimally close, it has been proved to still have Gribov copies [STSF82, DZ89, DZ91, vB92]. To see this, it is better to study an equivalent definition of Ω.

Let us consider the set of relative minima of:

A θ 2 = 1 2 d d xA θ,a µ (x)A θ,a µ (x) (2.2)
as a function of the gauge parameter θ at fixed A µ . We are going to show that the region of A θ µ defined by this set is the first Gribov region. First of all, we are going to demonstrate that the extrema of A θ 2 are the configurations in Landau gauge.

The extrema are obtained by taking δ A 2 = 0 where δ represents an infinitesimal gauge variation around a configuration A µ .

δ A 2 = d d x (δA a µ (x))A a µ (x) = -d d x (D µ θ) a A a µ (x) = -d d x ∂ µ θ a (x)A a µ (x) = d d x θ a (x)∂ µ A a µ (x) = 0 (2.3)
where in the second line we have used that the second term of the covariant derivative vanishes since f abc A b µ A a µ is zero by interchanging a and b. The relative extrema fulfils the above condition for all infinitesimal gauge variations θ and thus it implies that ∂ µ A a µ = 0, the Landau gauge condition required by the first definition.

Secondly, for the minima the matrix of second derivatives should be positive. This condition forces the FP operator to be positive definite because:

δ(δ A 2 ) = dxθ a (x)(-∂ µ (D µ θ) a ) > 0
This new definition of the Gribov region i.e. the relative minima of { A θ 2 } makes us doubt on the non-existence of Gribov copies inside Ω. It seems natural to think that there may be many local minima of { A θ 2 }. In fact, an explicit example is given in [vB92] and, moreover, in lattice simulations a lot of Gribov copies have been found.

In order to consider only one representative configuration for every gauge orbit one should take, for example, the absolute minimum of { A θ 2 }. This set is called the fundamental modular region, FMR [START_REF] Zwanziger | Fundamental modular region, Boltzmann factor and area law in lattice gauge theory[END_REF]. The FMR is included in the Gribov region and is the ideal restriction of the functional integral domain. However, it is more difficult to work with the FMR than with the first Gribov region. The difficulty lies in how to include this constraint in the functional integral. Some interesting ideas were developed by Gribov using the semiclassical approach [START_REF] Gribov | Quantization of Nonabelian Gauge Theories[END_REF] and by Zwanziger in [START_REF] Zwanziger | Nonperturbative Faddeev-Popov formula and infrared limit of QCD[END_REF] including a new local term in the action.

Properties of the Gribov region

It is important to mention that every gauge orbit intersect the Gribov region at least once, see [START_REF] Dell'antonio | Every gauge orbit passes inside the Gribov horizon[END_REF]. Therefore, the Gribov region contains informations about all configurations.

Another important property is that A µ = 0 belongs to the Gribov region. We have already seen that, for A µ = 0, the FP operator is positive definite, see [START_REF] Zwanziger | Nonperturbative Faddeev-Popov formula and infrared limit of QCD[END_REF]. This means that the perturbative region is included in Ω. The Gribov region is also convex and bounded in every direction, see [START_REF] Zwanziger | Nonperturbative Faddeev-Popov formula and infrared limit of QCD[END_REF].

In relation with Gribov copies, it is known that every configuration in the first Gribov region near enough to its boundary, δΩ, presents Gribov copies outside Ω but near the δΩ, [START_REF] Gribov | Quantization of Nonabelian Gauge Theories[END_REF]. Moreover, it was proposed by Zwanziger in [START_REF] Zwanziger | Nonperturbative Faddeev-Popov formula and infrared limit of QCD[END_REF] that the important configurations lie on the intersection of the Gribov boundary, δΩ, with the boundary of the fundamental modular region, δΛ. Therefore, if it were confirmed, the copies in the interior of the Gribov region should not play any role but there is no proof for the moment for that.

Restriction to the Gribov region

The aim of this section is to briefly present some ideas used to restrict the functional integral in the generating function to the Gribov region,

Z[J] = Ω Dφ e -S inv -S FP -S GF +Ssource
where Dφ represents the product of DA µ DψD ψDcDc and S source is the part of the action that includes the external sources as in Sec. 1.2.1.

To implement the restriction it is better to modify the integrand using a factor V(Ω), to be determined later. So the generating function can be expressed as

Z[J] = Ω Dφ e -S inv -S FP -S GF +Ssource = Dφ V(Ω)e -S inv -S GF +Ssource-d 4 x{-c a ∂µ(Dµc) a })
(2.4)

The restriction of the integral Z to Ω implies the restriction to fields in which the FP operator is positive definite thus it has no zero modes. This implies that the ghost propagator in presence of a background gauge field in Landau gauge, that we denote D ab (p 2 , A µ ), does not have a non-trivial pole inside the first Gribov region (remember that the ghost propagator is related to the inverse of the FP operator which does not have zero modes inside the Gribov zone). This is known as the Gribov no-pole condition because in [START_REF] Gribov | Quantization of Nonabelian Gauge Theories[END_REF] Gribov proposed using this condition of the ghost propagator in order to define the constraint for the Gribov region.

In order to study the no-pole condition, Gribov proposed to compute the ghost propagator in presence of a background gauge field using a semi-classical approximation. The ghost propagator in presence of a background gauge field can be obtained as the second derivative with respect to the external sources J and J of the log I [A, J, J], where I [A, J, J] is defined as

I [A, J, J] = DcDc exp -d 4 xd 4 y ca (x)(-∂ µ D ab µ )c b (y) + d 4 x Ja c (x)c a (x) + ca (x)J a c (x) = Cdet (-∂ µ D µ ) exp -d 4 xd 4 y Ja c (x)(∂ µ D ab µ ) -1 J b c (y) (2.5)
where the character C represents a constant that will be ignore because no multiplicative constant survives in the correlation functions.

Gribov in [START_REF] Gribov | Quantization of Nonabelian Gauge Theories[END_REF] computed D ab (p 2 , A µ ) to the second order and found that

D ab (p 2 , A µ ) ∼ 1 p 2 1 1 -σ(p, A µ )
where σ is given by

σ(p, A µ ) = 1 V 1 p 2 N g 2 N 2 -1 d d q (2π) d A a µ (-q)A a ν (q) (p -q) µ p ν (p -q) 2 .
The no-pole condition for the ghosts forces σ to be less than one. Therefore, one obtains

V(Ω) = Θ(1 -σ(0, A µ ))
where Θ represents the step function.

In the way to rewrite the restriction to the Gribov region in a local action Zwanziger, in [START_REF] Zwanziger | Renormalizability of the critical limit of lattice gauge theory by BRS invariance[END_REF][START_REF] Zwanziger | Fundamental modular region, Boltzmann factor and area law in lattice gauge theory[END_REF], described this change using the horizon function, S H ,

S H = γ 4 d 4 xH(x)
where H is proportional to

F -1 H(x) = g 2 f abc A b µ F -1 ad f dec A c µ
and γ has dimension of mass and it can be determined as function of g and Λ QCD through the horizon condition [Zwa89, Zwa93]

H(x) = 4(N 2 -1).
The restriction to the Gribov region was done under the assumption that the relevant Gribov copies are located near the border of the first Gribov region and therefore we can approximate the step function into a δ-function. Using this conjecture the functional integral restricted to the first Gribov region can be modified as:

Ω → e -S H .
This introduces a non-local modification in the action that suppresses the probability of configurations near the boundary δΩ. (A pedagogical introduction to this subject can be found in [START_REF] Sobreiro | A Study of the Gribov copies in linear covariant gauges in Euclidean Yang-Mills theories[END_REF][START_REF] Vandersickel | A Study of the Gribov-Zwanziger action: from propagators to glueballs[END_REF].) It is worth mentioning that the horizon function, which is non-local, can be localized through the help of additional fields [Zwa89, Zwa93],

e -S H = DϕD φDωD ωe -S loc .

The term S loc has the form

S loc = d 4 x -φac µ F ab ϕ bc µ -ωac µ F ab ω bc µ -γ 2 g d 4 x f abc A a µ ϕ bc µ + f abc A a µ φbc µ ,
where the fields φµ and ϕ µ are a pair of complex conjugate bosonic fields, and ωµ and ω µ are fermionic fields. Each of them with 4(N 2 -1) 2 components.

To conclude, under some assumptions, it is possible to express the restriction to the Gribov region with a local action, usually called Gribov-Zwanziger action, S GZ as follows:

Ω DA µ δ(∂ µ A µ )det(F)e -S inv = DA µ DψD ψDϕD φDωD ωDbDcDce -S GZ
where

S GZ =S inv + d 4 x b a ∂ µ A a µ -ca F ab c b + d 4 x -φac µ F ab ϕ bc µ -ωac µ F ab ω bc µ -γ 2 g d 4 x f abc A a µ ϕ bc µ + f abc A a µ φbc µ .
(2.6)

The behaviour of the propagator using the Gribov-Zwanziger action at tree level shows an enhanced ghost propagator and a gluon propagator which vanishes at the origin. Unfortunately, this is not the behaviour found by lattice simulation, discussed in Sec.2.3.

The study of the restriction to the Gribov region has gone beyond Gribov-Zwanziger action with the introduction of the refined Gribov-Zwanziger action, S RGZ , [DGS + 08b].

It consists of taking into account the two-dimension non-trivial condensate 0| φab

µ ϕ ab µ -ωab µ ω ab µ |0 by the introduction of the operator φab µ ϕ ab µ -ωab µ ω ab µ directly in the action. The action includes the term,

S φϕ = -d 4 xµ 2 φab µ ϕ ab µ -ωab µ ω ab µ
where the µ parameter can be fixed by a variational principle. The refined Gribov-Zwanziger action is therefore,

S RGZ = S GZ + S φϕ .
This new action shows a non-vanishing gluon propagator in the infrared and a ghost propagator that goes as 1 p 2 when p → 0 [START_REF] Sorella | Gribov horizon and BRST symmetry: A Few remarks[END_REF] as it is expected by lattice simulations, see [START_REF] Cucchieri | Constraints on the IR behavior of the ghost propagator in Yang-Mills theories[END_REF].

Both the GZ and RGZ models are renormalizable. However, there are a few remarks that need to be made. The GZ and RGZ actions are no longer BRST invariant. In fact, this symmetry is softly broken. This is enough to prove renormalizability but the definition of the physical Hilbert space using the cohomology of the BRST symmetry is no longer valid. The problem of the unitarity of the theory is unsolved. Another inconvenient with this approach is that the calculations in this model are cumbersome since there are several more fields in the action, however some corrections for the propagators were computed in [START_REF] Gracey | Alternative refined Gribov-Zwanziger Lagrangian[END_REF][START_REF] Gracey | Loop calculations in the three dimensional Gribov-Zwanziger Lagrangian[END_REF].

Dyson-Schwinger approach

There exists a non-perturbative analysis of quantum field theories introduced by Dyson and Schwinger [START_REF] Dyson | The S matrix in quantum electrodynamics[END_REF][START_REF] Schwinger | On the Green's functions of quantized fields. 1[END_REF] that can be used to study the infrared region of QCD. The Dyson-Schwinger equations is an infinite set of coupled non-linear integral equations relating the different Green's functions.

The Dyson-Schwinger equations can be found through the observation that the integral of a total derivative is zero. If we apply this to the generating function Z[J], we obtain

0 = Dφe -S[φ]+ d 4 xφJ -∂S[φ] ∂φ + J = - ∂S ∂φ -J J
as usual, φ represents any arbitrary field and J the source associated with it. The Dyson-Schwinger equation is therefore,

∂S ∂φ -J J = 0. (2.7)
It is worth mentioning that the source J can be related to the effective action Γ[ φ ], by

J(x) = δΓ δ φ(x)
.

The set of infinite equations can be obtained from successive derivations of the Dyson-Schwinger equation (2.7), and by setting the sources to zero.

In order to clarify how to obtain these equations let us study an example. We are going to focus on the ghost-gluon sector, that means that we are going to ignore the quarks. This restriction is usually called quenched approximation. The action has the form

S = d 4 x 1 4 F a µν F a µν -ca ∂ µ (∂ µ + gf abc A b µ c c ) + ih a ∂ µ A a µ .
Note that this action corresponds to the Lagrangian of Eq. (1.14) without including quarks.

In this case, we have to consider sources J a µ , χa , χ a and R a for the gluon, ghost, antighost and h a respectively, so the action associated with the sources is

S sources = d 4 x J a µ A a µ + χa c a + ca χ a + R a h a .
Then, as an example, the Dyson-Schwinger equation taking φ = ca is

- ∂S ∂c a (x) + χ a (x) J = 0.
The next step is to take the derivative respect to χ b (y). We have to be careful with the signs because the derivative is made with respect to a fermionic field. At this point one obtains

- ∂S ∂c a (x) cb (y) = -δ ab δ(x -y).
The derivative ∂S ∂c a (x) can be calculated explicitly

∂S ∂c a (x) = -∂ µ D µ c a .
Thus,

∂S ∂c a (x) cb (y) = (-∂ 2 ) c a (x)c b (y) -gf adc ∂ x µ d 4 z d 4 z A d µ (z)c c (z )δ(z -x)δ(z -z )c b (y) = (-∂ 2 )D ab (x -y) + d 4 z d 4 z(-gf adc ∂ x µ δ(z -x)δ(z -z ) c c (z )c b (y)A d µ (z) . (2.8)
In the last line we have identify the ghost propagator D ab (x -y) = c a (x)c b (y) .

Furthermore, the connected correlation function c c (z )c b (y)A d µ (z) can be expressed in term of the three-point vertex as (see the Appendix for a perturbative explanation)

c c (z )c b (y)A d µ (z) = -d 4 ud 4 vd 4 wD de µν (z -u)D cf (z -v)Γ ef g ν (u, v, w)D gb (w -y).
Here, D de µν (z -u) = A d µ (z)A e ν (u) is the gluon propagator and Γ ef g ν (u, v, w) represents the ghost-gluon vertex function defined as:

Γ ef g ν (u, v, w) = δ 3 Γ δc g (w)δc f (v)δA e ν (w)
.

To conclude with the example one can include these elements in the last line of Eq. (2.8) to obtain

(-∂ 2 )D ab (x -y) -d 4 ud 4 vd 4 wd 4 zd 4 z (-gf adc (∂ z µ δ(z -x))δ(z -z )D de µν (z -u) × D cf (z -v)G ef g ν (u, v, w)D gb (w -y) = d 4 wδ(x -w)δ ag (-∂ 2 ) -d 4 ud 4 vd 4 wd 4 zd 4 z (-gf adc (∂ z µ δ(z -x))δ(z -z ) ×D de µν (z -u)D cf (z -v)G ef g ν (u, v, w) D gb (w -y) = δ ab δ(x -y) (2.9)
Therefore the expression within the braces is the inverse of the full ghost propagator. Moreover, we can identify in the first term the tree-level of the ghost propagator, (-∂ 2 ). In the second term the -gf

adc (∂ z µ δ(z -x))δ(z -z ) = gf acd ∂ z µ δ(z -x)δ(z -z )
is the tree level of the ghost-gluon vertex (see the Appendix A.2). Taking into account all these observations, Eq. (2.9) can be represented diagrammatically as shown in figure 2 We can see that nth-derivative of Γ is related to higher order derivatives. This gives a infinite system with coupled equations. It is worth noting that until now the procedure to obtain the DS equation is exact, no approximation was done. However, to go further and obtain more information some approximations are needed. In general, one implements a truncation scheme and give an ansatz for the unknown higher order vertex functions.

= -1 6 -1 -1 2 + -1 -1 2 -1 2 -1 2
The first approximation made for the gluon DS equation was proposed by Mandelstam in [START_REF] Mandelstam | Approximation Scheme for QCD[END_REF]. It consists in neglecting all ghost contributions, neglecting the two loop diagrams and taking the vertices as the bare ones. However, for the Landau gauge the ghosts are of great importance in the infrared [START_REF] Lorenz Von Smekal | A solution to coupled Dyson-Schwinger equations for gluons and ghosts in Landau gauge[END_REF][START_REF] Lorenz Von Smekal | The Infrared behavior of gluon and ghost propagators in Landau gauge QCD[END_REF].

In fact, there exists a solution of the DS equation where the ghost propagator is enhanced in the infrared. This solution is the so called scaling solution [START_REF] Fischer | Uniqueness of infrared asymptotics in Landau gauge Yang-Mills theory[END_REF]. It is obtained by assuming that all the one-particle irreducible Green's function with 2n c external ghost legs and n A external gluon legs follow an infrared scaling law [START_REF] Alkofer | Vertex functions and infrared fixed point in Landau gauge SU(N) Yang-Mills theory[END_REF][START_REF] Huber | The Infrared behavior of Landau gauge Yang-Mills theory in d=2, d=3 and d=4 dimensions[END_REF]

of the form Γ (nc,n A ) (p 2 ) ∼ (p 2 ) (nc-n A )κ+(1-nc)(d/2-2)
where d is the space-time dimension and κ the anomalous dimension. For the ghost dressing function and the gluon propagator, J(p) and D(p) respectively, defined by

D ab µν (p) = δ ab δ µν - p µ p ν p 2 D(p), D ab (p) = δ ab J(p) p 2
the scaling solution implies that in the infrared for d=4

D(p) ∼ (p 2 ) κ A -1 , J(p 2 ) ∼ (p 2 ) -κc .
In four dimensions the relation κ A = 2κ C holds where κ C = κ. For κ A > 1 the gluon propagator goes to zero in the infrared while the ghost propagator goes to infinity quicker than the tree-level propagator. This solution satisfies the Gribov-Zwanziger confinement condition [START_REF] Zwanziger | Nonperturbative Faddeev-Popov formula and infrared limit of QCD[END_REF] and the Kugo-Ojima confinement scenario [START_REF] Kugo | Local Covariant Operator Formalism of Nonabelian Gauge Theories and Quark Confinement Problem[END_REF]. Among others, the scaling solution has been studied by [START_REF] Lorenz Von Smekal | The Infrared behavior of gluon and ghost propagators in Landau gauge QCD[END_REF][START_REF] Alkofer | The Infrared behavior of QCD Green's functions: Confinement dynamical symmetry breaking, and hadrons as relativistic bound states[END_REF][START_REF] Fischer | Nonperturbative propagators, running coupling and dynamical quark mass of Landau gauge QCD[END_REF][START_REF] Jacques | Two loop improved truncation of the ghost gluon Dyson-Schwinger equations: Multiplicatively renormalizable propagators and nonperturbative running coupling[END_REF][START_REF] Fischer | On the infrared behavior of Landau gauge Yang-Mills theory[END_REF].

On the other hand, lattice simulations [BIMPS07, CM07, SvSLW07, CM08a, CM08b] show an infrared finite gluon propagator and more important an infrared finite ghost dressing function. A infrared finite gluon can also be obtained by the scaling solution considering κ A = 1, but even in this case the ghost propagator exponent is enhanced with respect to the tree level one.

A new solution called the decoupling solution was studied in [Blo03, AN04, BBL + 06, AP08, ABP08, BLLY + 08, RQ11a, HvS13]. The decoupling solution succeeds in recovering the lattice prediction of the propagators. Its name is due to the fact that the gluons decouple below a scale corresponding to their mass. The generation of an effective dynamic gluon mass was first elaborated by [START_REF] Cornwall | Quark Confinement and Vortices in Massive Gauge Invariant QCD[END_REF][START_REF] Cornwall | Dynamical Mass Generation in Continuum QCD[END_REF]. The infrared behaviour imposed by the decoupling solution for the dressing functions are

D(p) ∼ const J(p) ∼ const.
It is worth mentioning that the scaling solution is unique while there is a family of decoupling solutions [START_REF] Fischer | On the infrared behavior of Landau gauge Yang-Mills theory[END_REF]. One of the decoupling solutions seem to reproduce all lattice results even for the simulation of the vertex functions, see e.g. [START_REF] Aguilar | Ghost propagator and ghost-gluon vertex from Schwinger-Dyson equations[END_REF][START_REF] Aguilar | A new method for determining the quark-gluon vertex[END_REF].

Lattice simulation

All the techniques presented before in this chapter have the aim of studying non-perturbative QCD analytically. However, we have seen that in the infrared regime some approximations have to be done in order to go further in the investigation. The best way to choose the approximation is to use the results obtained by lattice simulations as a benchmark.

Lattice simulations are the most reliable tool to investigate the infrared regime. The aim of lattice studies is to perform the functional integral explicitly with the help of computers. The functional formalism contains an integral over all possible field configurations with a weight analogous to a statistical mechanics Boltzman one. To compute the path integral on the lattice a finite number of configuration are taken into account. We are going to return to the path integral in the subsection 2.3.3.

In order to compute the action it is necessary to discretize the space considering only a finite set of points {x i } on it. Lattice simulations are realized on an hypercubic fourdimensional lattice where the field configurations are simulated (the form of the lattice could also be changed).

Let us define the four dimensional lattice of length L with site-spacing a, as it is shown is figure 2.4. The space-time integral becomes a sum over all lattice sites, d 4 x → a 4

x and we use an equivalent discretized form of the Lagrangian. That means that Langrangian depends on the value of the fields evaluated at the lattice sites, φ(x i ). In the lattice approximation we only know the fields at the lattice sites. Thus all derivatives in the Lagrangian must be converted to finite differences.

The size of the lattice spacing has to be chosen small enough so that the discretization remains a good approximation. But for smaller lattice spacing the cost of the simulation is larger. The difficulty is to choose a lattice spacing as large as possible to make the simulations treatable while keeping the discretization errors small. The lattice spacing also acts as a ultraviolet regulator. The gluons and quarks with momenta larger than π/a are not reachable by the lattice simulations. This cut-off introduces an error because in a quantum field theory like QCD ultraviolet modes strongly affect infrared ones. This effect is usually reduced by the renormalization of the parameters or by adding new local interactions.

Gluon sector

In this subsection we are going to present the basic ideas of lattice QCD in the gluon sector, a pedagogical review can be found in [START_REF] Lepage | Lattice QCD for novices[END_REF][START_REF] Davies | Lattice QCD[END_REF][START_REF] Sternbeck | The Infrared behavior of lattice QCD Green's functions[END_REF]. A detailed explanation and more modern strategies can be found e.g. in [START_REF] Montvay | Quantum fields on a lattice[END_REF].

In this case we are going to focus on the dynamical term of the gluons

S YM = d 4 x 1 4 F a µν F a µν .
Let us remark that it is extremely difficult to formulate a lattice version of QCD directly in terms of the gluon field A µ , because its gauge transformation is complicated. For the gluon sector it is better to describe the theory using the link variable U µ (x),

U µ (x) = Pe ig 1 0 Aµ(x+atμ)dt
, where U µ (x) is an element of SU (3) and P denotes the path ordering along the integration from x to x + μ. The link variable U µ (x) is represented pictorially by a directed line from x to x + aμ. The gauge transformation of U is simply,

U µ (x) → U (x)U µ (x)U † (x + aμ).
Gauge invariant objects can easily be made out of links variables forming a closed loop.

The smallest nontrivial gauge invariant object built from the link operators is the plaquette operator P x,µν which involves the product of link variables around the smallest square at site x.

P x,µν = 1 N ReTr U µ (x)U ν (x + aμ)U † µ (x + aμ + aν)U † ν (x + aν)
The plaquette allows us to describe the action in terms of the link variables since

P x,µν ∼ 1 -c 1 a 4 T r(gF µν (x)) 2 + O(a 6 )
where c 1 is a constant. The simpler gauge invariant action that can be built, is the so called Wilson action

[Wil74] S Wil [U] = β x 1≤µ<ν≤4
(1 -P x,µν )

where β = 2N g 2 and N = 3 for QCD. The continuum expression for the Wilson action takes the form of S YM up to order a 2 , in fact,

S Wil = d 4 x 1 4 T r(F µν ) 2 + O(a 2 )
The error in the discretization of the action can be improved by adding a rectangle operator which cancels the order a 2 of the Wilson action [START_REF] Lepage | Lattice QCD for novices[END_REF]. In this way the discretization of the classical action is achieved to order a 4 . There are different ways to define the lattice action of QCD. All of them are required to be gauge invariant and to be equal to the continuum action when the lattice spacing vanishes.

It is important to remark that lattice simulations respect gauge invariance but the space-time discretization softly breaks Lorentz symmetry, rotational symmetry, etc. Of course those symmetries are recovered once the continuum limit is taken. However, it is necessary to impose gauge invariance in order to have the same coupling constant for all the interactions. In this way we reduce the number of parameters to work with.

Quark sector

Working with fermions on the lattice is much more difficult. The dynamical term of the continuum action of the fermions is

S f = d 4 x ψ (-γ µ ∂ µ + M ) ψ.
A naïve discretisation leads us to consider

S latt,naïve = a 4 x    ψx 4 µ=1 -γ µ ψ x+μ -ψ x-μ 2a + M ψx ψ x   
This action gives problems in the computation of the propagator. In order to see them let us consider the propagator in Fourier space:

-iγ µ Sin(p µ a) a + M which in the continuum takes the form -iγ µ p µ +M . The problem arises in the Sin function because it goes to zero when p ∼ 0 but also when p ∼ π a , this is called the doubling problem. There are various possible definitions for the fermion action, e.g. the Wilson fermion action, the Staggered fermion action, the clover-improved Wilson action, etc. The important clue is to impose the fermionic behaviour in the functional integral, which means that if the quark action has the form S f = a 4

x ψx Q xy ψ y , where Q is the fermion matrix that depends on the gauge field A µ , then

DψD ψe -S f = det(Q).
Having integrated the fermions the QCD lattice action can be written as

S lattice = β x 1≤µ<ν≤4 (1 -P x,µν ) -log (det(Q)) .
Even in the unquenched case, to compute the path integral we have to integrate over all configurations of U (the gluon field) and the role of quarks is to introduce the extra term log (det(Q)). Let point out, however, that the det(Q) is more difficult to compute than the Wilson action.

Path integral

A mean value in the functional formalism can be computed as the path integral over all possible configuration of the field weighted by the exponential of the action properly normalised,

O = DU O e -S DU e -S
To compute this integral in the lattice we can generate the configurations U α with a density probability e -S lattice and substitute the correlation function by

O = 1 N conf N conf α O[U α ] (2.10)
where N conf is the number of configurations. The observable is calculated as the average of the operator evaluated in the field configurations which were generated with a probability e -S .

Landau gauge

We want to emphasize that in lattice simulations it is not necessary to fix the gauge in order to calculate gauge invariant physical observables. However, as the lattice simulations is the most controlled tool to investigate the infrared zone it is useful to use these simulations to compare with analytical computations. Generally the analytical calculations are done in a given gauge. In the lattice approach it is possible to fix the gauge and the procedure avoid the Gribov problem.

For fixing the gauge in the lattice, the basic idea is to generate the set of configurations {U α } using some Monte Carlo method with gauge-invariant action. This set does not respect any gauge condition in particular. Then, one applies a gauge transformation to take the gauge configuration through the gauge orbit until it satisfies the gauge condition. Let us call U θ µ (x) the gauge transformation of U µ (x) that satisfies a given gauge condition,

U θ µ (x) = U (x)U µ (x)U † (x + μ).
With these configurations we compute the expectation value of a given quantity by

O = 1 N conf N conf α O[U θ α ].
It is important to observe that this procedure does not affect gauge invariant quantities because

O[U α ] = O[U θ α ].
In order to find the gauge configuration U θ µ (x) satisfying the Landau gauge condition we are going to consider the functional

F U [θ] = 1 4 x 4 µ=1 ReTrU θ µ (x).
This functional is the lattice analogous to Eq. (2.2). In this case all the extrema of the functional F U [θ] fulfil the Landau gauge condition. There exists many extrema, each of them represent a Gribov copy. If we want to fix the gauge it is enough to consider only one extremum. Contrarily to the analytical analysis, the choice of one and only one maximum on the lattice can be achieved without much problem using standard algorithms and this avoids the Gribov problem. There exist different algorithms for gauge-fixing in the lattice, for instance the over relaxation [START_REF] Mandula | Efficient gauge fixing via overrelaxation[END_REF] and Fourier-accelerated gauge-fixing [DBK + 88], a comparison of both methods can be found in [START_REF] Cucchieri | Critical slowing down in SU(2) Landau gauge fixing algorithms[END_REF][START_REF] Sternbeck | The Infrared behavior of lattice QCD Green's functions[END_REF]. The problem is that the non-perturbative corresponding continuum limit of the action is unknown.

The Landau gauge condition on the lattice can be expressed as ∇ µ A θ µ (x) = 0 where ∇ denotes the finite difference,

∇ µ A µ (x) = 4 µ=1 {A µ (x + μ/2) -A µ (x -μ/2)}
The A µ can be defined, for instance, at the midpoint of a link as

A µ (x + μ/2) = 1 2i U µ (x) -U † µ (x) - 1 6i Tr U µ (x) -U † µ (x)
having a precision of O(a 2 ).

In this thesis we are going to compute correlation functions, that are gauge-dependent quantities. Our calculations will be done in Landau gauge because this is the gauge condition in which more simulations have been done and accordingly, the results can be better tested. To conclude with the basic ideas of lattice simulation we are going to present how to compute the gluon and ghost propagator on the lattice.

On the lattice, the gluon propagator is the Monte Carlo average of the corresponding two-point function

D ab µν (x, y) = A a µ (x)A b ν (y) here A a µ (x) = A µ (x + μ/2) using U µ (x) satisfying the Landau-gauge condition.
The discretised Fourier transform of the propagator is

D ab µν (q) = 1 V x,y A a µ (x)A b ν (y)e iq.(x+μ/2) e -iq.(y+ν/2) .
In the continuum the tensor structure of the gluon propagator takes the form

D ab µν (q) = δ ab δ µν - q µ q ν q 2 D(q 2 ).
Assuming that lattice propagator has the same tensor structure as in the continuum limit, D(q 2 ) can be computed as

D(q 2 ) = 1 N 2 -1 1 d -1 a,µ D aa µµ (q).
The gluon propagator has been simulated in several opportunities, see e.g. [BCLM04, FN04a, CM06, CM07, CMM08, SMP12] for quenched (without quarks) SU (2), [LSWP98, BBLW00, FN04a, BLLY + 05, BIMPS07, SvSLW07, CMM08] for quenched SU (3), [BHL + 04, FN06, IMPS + 07, SMMPvS12] for full QCD simulations. The figure 2.5 shows the typical behaviour of the gluon propagator which is qualitatively the same for both quenched and unquenched cases and for SU (2) and SU (3). The gluon propagator is found to be finite at low and zero momentum [OS09, BIMPS07, CM07, SvSLW07, CM08a, CM08a, CM08b, CM10, BIMPS09, DOV10]. Another interesting feature is the violation of the positivity of the spectral function. In the Euclidean space, the spectral function, ρ(σ), is related to the gluon propagator, by

D(p) = 1 π ∞ 0 dσ ρ(σ) p 2 + σ .
The spectral function of physical fields should be positive. In Landau gauge lattice simulations found that the gluon spectral function has a negative part [BPS94, LSWP98, FN04b, CMT05, SIMP + 06, BHL + 07, DSVV08, Cor13]. This may be related to gluon confinement as the gluon propagator does not correspond (even in the transverse sector) to physical particles.

In Landau gauge, the algorithms to fix the gauge can be thought as a standard algorithm for maximizing a function. Therefore there is no need for introducing explicitly the ghosts. However we can obtain information about the ghost propagator through the FP operator F. The FP operator is defined only as a function of gluons. On the lattice the FP operator is computed by

F[U] = -∇.D[U]
where D[U] denotes the covariant derivative. It can be shown that in the lattice representation

(D µ [U]) ab xy = 2ReTr[t a t b U x µ ]δ x+μ,y -2ReTr[t a t b U x µ ]δ
x,y Thus, the ghost propagator in momentum space is obtained by

D ab (q) = 1 V x,y F -1 ab xy e iq.(x-y) .
If we extract the color factor D ab (q) = δ ab J(q) q 2 then J(q)

q 2 = 1 N 2 -1 a D aa (q) = 1 N 2 -1 TrF -1 .
A typical behaviour of the ghost dressing propagator in d = 4 is shown in figure 2.6, [START_REF] Cucchieri | Three-point vertices in Landaugauge Yang-Mills theory[END_REF]. The infrared divergence of the ghost propagator has the same order as the one appearing at tree level in three and four dimensions [START_REF] Cucchieri | Constraints on the IR behavior of the ghost propagator in Yang-Mills theories[END_REF]. Contrarily to the d = 4 and d = 3 case where in the infrared the ghost propagator diverges as the tree-level, in two dimension the ghost propagator diverges faster as p -2-2κ with anomalous dimension κ between 0.1 and 0.2. As we mentioned earlier there is no Gribov problem on the lattice. One could want to try to implement in analytical calculations the same gauge fixing as the one used on the lattice. The problem is that no representation of this fixing in terms of an action is known.

Another option could be to try to implement directly a BRST invariant action on the lattice. This was proposed by Neuberger in [START_REF] Neuberger | NONPERTURBATIVE BRS INVARIANCE[END_REF] where a general BRST invariant Lagrangian is presented. As in the continuum it is the solution to Q s L = 0 where Q s is the charge generator of BRST symmetry. The solution found has the form L 0 + Q s L where L 0 is gauge invariant and does not depend on the antighost and the ghost fields. This fact makes the functional integral to vanish [START_REF] Neuberger | Nonperturbative BRS Invariance and the Gribov Problem[END_REF] due to the antighost integral. In order to see this, let us define in the functional integral a BRST invariant measure, which has the form

dµ = x,a dc a (x)dc a (x)dh a (x) × x dψ(x)d ψ(x) µ dU µ (x) ρ(h)
where U µ are the link variables and ρ(h) is a function that ensures the convergence of the h-integral. This measure has the property that any integral of the form dµQ s R vanishes for any operator R, see [START_REF] Neuberger | Nonperturbative BRS Invariance and the Gribov Problem[END_REF]. If we consider a BRST invariant operator O = O 0 + Q s Ô where O 0 is gauge invariant we can define the function

f (t) = dµe -L 0 +tQs LO.
For t = 1 the function matches with the numerator of the expectation value of O. We can observe that the derivative df dt is zero since

df dt = dµQ s e -L 0 +tQs L LO = 0.
Therefore,

f (1) = f (0) = dµe -L 0 O = dµe -L 0 O 0 + dµe -L 0 Q s Ô
as Q s e -L 0 = 0 the second integral is zero. Moreover, the first integral is also zero, dµe -L 0 O 0 since the integrand does not depend on the ghost field. This implies that

dµe -L 0 +Qs LO = 0.
This zero is independent of the physics and should be factorized from the integral and would have to cancel with the corresponding zero in the denominator of the expectation values. This problem is similar to FP problem where we want to factorize the gauge volume, in this case we wonder how to factorise these zeros both in the numerator and in the denominator of a Green's function. Perturbatively these zeros are associated with disconnected vacuum diagrams and in the average of gauge invariant quantities this spurious 0/0 cancels. Therefore, to any order of perturbation theory the Green's functions of gauge invariant quantities are the same as the ones computed with a gauge invariant action. On the other hand, if we want to compute the expectation value non perturbatively we do not know a way to eliminate these zeros. In consequence the problem of finding a continuum action in the non-perturbative (infrared) regime remains unsolved.

Curci-Ferrari model in Landau gauge

The Dyson-Schwinger equations allowed to study the infrared regime reproducing with reasonable accuracy the lattice data. However, it could be nice to have a more manageable method to study the infrared regime, for example, pertubation theory. We wonder which is the simplest Lagrangian that allows one to make perturbative expansion at low momenta, reproduces the lattice data, and matches QCD at high momenta. The aim of this thesis is to prove that the simplest modification of the Lagrangian consisting of adding a gluon mass term reproduces the lattice results with good accuracy and at the same time, gives the standard QCD results in the ultraviolet. That means that a standard one-loop calculation of the propagators and vertices gives the expected results obtained by lattice simulations for both the quenched and the unquenched case in the infrared regime. The following chapters will show a detailed computations of the propagators and vertices and the corresponding matching with lattice data.

The Lagrangian we are going to work with is a particular case of Curci-Ferrari model [START_REF] Curci | A Theory of Massive Yang-Mills Field[END_REF][START_REF] Curci | On a Class of Lagrangian Models for Massive and Massless Yang-Mills Fields[END_REF] but setting Landau gauge condition. In the Euclidean space the Lagrangian density has the form

L = 1 4 F a µν F a µν + ∂ µ ca (D µ c) a + ih a ∂ µ A a µ + m 2 2 A a µ A a µ + N f i=1 ψi (-γ µ D µ + M i )ψ i , (2.11)
where, we recall that

F a µν = ∂ µ A a ν -∂ ν A a µ + gf abc A b µ A c ν , (D µ c) a = ∂ µ c a + gf abc A b µ c c , D µ ψ = ∂ µ ψ -igA a µ t a ψ.
The Curci-Ferrari model was introduced in the mid-seventies as an alternative to the Higgs mechanisms for describing massive vector mesons. However, even though the model did not succeed at that time, it was useful to study the infrared problem of massless Yang-Mills theories.

The problem of the change in the Lagrangian is that the addition of a mass term breaks BRST symmetry. BRST symmetry and the fact that it is nilpotent was used to prove renormalizability at all orders using the Slavnov-Taylor identities and to define the physical space by the cohomology of the symmetry. However, the Curci-Ferrari Lagrangian breaks the BRST symmetry softly and it is possible to prove that the theory is renormalizable [START_REF] Curci | On a Class of Lagrangian Models for Massive and Massless Yang-Mills Fields[END_REF]. The multiplicative renormalization factors were computed at one loop explicitly in [START_REF] De Boer | On the renormalizability and unitarity of the Curci-Ferrari model for massive vector bosons[END_REF][START_REF] Andriyash | One loop renormalization of the Yang-Mills theory with BRST invariant mass term[END_REF] and more recently at two and three loops in M S in [START_REF] Gracey | Two loop MS-bar renormalization of the Curci-Ferrari model[END_REF][START_REF] Gracey | Three loop renormalization of QCD in various nonlinear gauges[END_REF].

It is worth mentioning that the procedure used in chapter 1 to define the Hilbert space of physical states cannot be applied to the Curci-Ferrari Lagrangian. Furthermore [START_REF] Ojima | Comments on Massive and Massless Yang-Mills Lagrangians With a Quartic Coupling of Faddeev-popov Ghosts[END_REF][START_REF] De Boer | On the renormalizability and unitarity of the Curci-Ferrari model for massive vector bosons[END_REF] investigated the BRST algebra and explicitly constructed a negative norm state, showing that the theory is not unitary with the standard definition of the physical space. The only hope to prove the unitarity is to define another physical space with the properties mentioned in Sec.1.2.4. It is important to mention that the same difficulty takes place each time that BRST symmetry is broken as happens in lattice simulations, in Gribov-Zwanziger action and in Dyson-Shwinger approach, that is in all non-perturbative gauge-fixed analysis of QCD.

Symmetries and identities

In this subsection we want to discuss the remaining symmetries of the Lagrangian and the Slavnov-Taylor identities produced by them. These identities will be useful later on when we have computed the vertices to test the validity of the calculations.

The Slavnov-Taylor identities [Tay71a, Sla72] are a generalization of Ward identities [START_REF] Clive | An Identity in Quantum Electrodynamics[END_REF][START_REF] Takahashi | On the generalized Ward identity[END_REF]. They are obtained by using the symmetries of the action in the functional integral. In order to derive them let us consider that the action and the measure are invariant under the infinitesimal transformation φ → φ = φ + δφ. Thus, the generating functional, Z[J], can be written as

Z[J] = e W [J] = Dφe -S[φ]+Ssources[φ] = Dφ e -S[φ ]+Ssources[φ ] = Dφe -S[φ]+Ssources[φ]+δSsources ∼ Dφe -S[φ]+Ssources[φ] (1 + δS sources ) .
Therefore the Slavnov-Taylor identity is

δS sources = d d x δΓ δ φ δφ = 0.
(2.12)

In the particular case where δφ is linear in the fields we recover the Ward identities.

For simplicity we are going to focus on the quenched case. As we mentioned above the Lagrangian of Eq. (2.11) does not respect the nilpotent BRST symmetry but it is invariant under the soft-BRST symmetry, which for an arbitrary field φ is φ → φ + ηsφ where η is a Grassmann variable. The soft-BRST symmetry is defined as

sA a µ = (D µ c) a , sc a = - g 2 f abc c b c c , sc a = ih a and si(h a ) = m 2 c a .
(2.13)

For the soft-BRST symmetry is enough to consider the term of sources

S sources = d d x J a µ A a µ + χa c a + ca χ a + R a h a + Ka sA a µ + La sc a .
Remembering the relation between Γ and the action

Γ( A a µ , c a , ca , h a , Ka , La ) = -Dφe -S+Ssources + d d x J a µ A a µ + χa c a + ca χ a + R a h a (2.14)
we obtain that

δΓ δ A a µ = J a µ , δΓ δ c a = -χa , δΓ δ ca = χ a , δΓ δ h a = R a , δΓ δ Ka µ = -sA a µ and δΓ δ La = -sc a .
If we include these relations in the Slavnov-Taylor equation (2.12) we can write

d d xJ a µ sA a µ -χa sc a + sc a χ a + R a sh a = 0 = d d x J a µ sA a µ -χa sc a + sc a χ a + R a sh a = d d x - δΓ δA a µ δΓ δ Ka µ - δΓ δc a δΓ δ La + ih a δΓ δc a -im 2 δΓ δh a c a (2.15)
We can conclude that the Slavnov-Taylor identity associated with the soft-BRST is

d d x δΓ δA a µ (x) δΓ δ Ka µ (x) + δΓ δc a (x) δΓ δ La (x) -ih a (x) δΓ δc a (x) + im 2 δΓ δh a (x) c a (x) = 0 (2.16)
The second symmetry we are going to briefly mention is the conservation of the ghost number. One must associate a ghost number 1 to the ghost field, a -1 to the anti-ghost field and to Ka µ and a ghost number -2 to L while other fields are consider to have zero ghost number so the action has zero ghost number. The same property is fulfilled by Γ.

Another important relation between the derivatives of the effective action is the equation of motion for ca . It consists in translating the anti-ghost ca (x) → ca (x) + ¯ a (x), where is an infinitesimal x-dependent Grassman variable. This transformation is not a symmetry of the action and therefore we have to modify a little the condition (2.12) to δS -δS sources = 0.

Which in terms of reads

d d x(-∂ µ ¯ a (D µ c) a + ¯ a χ a ) = 0,
we can integrate by parts and express the source as the derivative of Γ to obtain

d d x¯ a (-∂ µ δΓ δ Ka µ + δΓ δ ca ) = 0.
As this relation holds for an arbitrary ¯ a , we can deduce the equation of motion related to ca as

∂ µ δΓ δ Ka µ (x) = δΓ δ ca (x)
.

(2.17)

More relations between the different derivatives of Γ can be obtained by taking the derivatives with respect to other fields. Moreover, the symmetry group of the action is bigger that the soft-BRST [START_REF] Delduc | A Note on Some Nonlinear Covariant Gauges in Yang-Mills Theory[END_REF][START_REF] Tissier | Gauged supersymmetries in Yang-Mills theory[END_REF], this group and its Slavnov-Taylor equations will be presented and studied in the Appendix together with the equation of motions of the theory.

Chapter 3

One loop corrections: Ghost and Gluon propagator

In the last chapter we have motivated the use of Curci-Ferrari model based on phenomenological arguments. The gluon propagator behaves as if it was massive at low momentum. In the Landau gauge the Curci-Ferrari model consists in adding a gluon mass term to the QCD Lagrangian,

L = 1 4 F a µν F a µν + ∂ µ ca (D µ c) a + ih a ∂ µ A a µ + m 2 2 A a µ A a µ + N f i=1 ψi (-γ µ D µ + M i )ψ i , (3.1)
where,

F a µν = ∂ µ A a ν -∂ ν A a µ + gf abc A b µ A c ν , (D µ c) a = ∂ µ c a + gf abc A b µ c c , D µ ψ = ∂ µ ψ -igA a µ t a ψ.
The aim of this thesis is to analyse to what extent this little modification can reproduce the lattice simulations of Landau gauge correlation functions by doing perturbation theory. The theory was already shown to be renormalizable [CF76, dBSvNW96, DS89, BM96] but the unitarity of the theory remains unproven. More than that lattice simulations show that at low momenta the standard definition of the physical space also includes negative norm states. On the opposite regime, this model matches with Yang-Mills theory in the limit of large momentum since, for momentum scales much larger than the mass of the gluons, this mass can be neglected. Some ideas to justify the presence of this new term from first principles can be found in [START_REF] Serreau | Lifting the Gribov ambiguity in Yang-Mills theories[END_REF].

In this chapter we are going to work in the quenched approximation meaning that quarks will be ignored. We are going to compute one loop corrections for the gluon and ghost propagators and for the ghost-gluon and three-gluon vertices. These corrections were computed several times in standard Yang-Mills theory, see e.g. [START_REF] Davydychev | Three gluon vertex in arbitrary gauge and dimension[END_REF][START_REF] Davydychev | One loop results for three gluon vertex in arbitrary gauge and dimension[END_REF]. The presence of a mass term yields a different propagator for the gluons and this fact makes the calculation of the vertices considerably harder.

In Sect 3.1 we will present the techniques used to compute the Feynman diagrams using the ghost propagator as an example. The calculation can be done by hand for the propagators and for the ghost-gluon vertex setting the gluon momentum to zero. More general cases are difficult to manage manually and the Mathematica software [START_REF][END_REF] is used. The expressions for the one-loop gluon and ghost propagator are presented in Sect 3.2. All the calculations are done in the Euclidean case for arbitrary space-time dimension, d.

In section 3.3 it is shown that our calculations compare very well with lattice simulations in d = 3, d = 4 both for SU (2) and SU (3). Later the renormalization group improvement is considered and we show the comparison with lattice data in this situation. Furthermore, an estimation of the error depending on the initial values of the parameters in the flow of the renormalization group equation is studied.

Feynman diagram techniques: Ghost propagator

In this section we are going to compute the one loop correction to the ghost propagator. The calculation will be presented in detail in order to introduce the techniques used to compute correlation functions in general. In Sect. 1.3 we presented the perturbative approach where it was established that when the coupling constant g is small it is convenient to separate the Lagrangian in two, the quadratic part and the interacting Lagrangian L int . The quadratic part gives Gaussian expectation values which are easier to compute due to Wick's theorem.

To compute a propagator in Fourier space, D(p), we are going to consider only the 1-PI diagrams since the propagator is related with Σ(p), the sum of all amputated 1-PI contributions, as

D -1 (p) = D -1 0 (p 2 ) -Σ(p 2
) where D 0 (p 2 ) is the free propagator. This relation holds if Σ commutes with the free propagators. This is clearly the case for ghosts and gluons but it is not so clear for quarks. However, for the quarks propagator the corresponding 1-PI diagrams depend only on / p, that means Σ = Σ( / p) and the propagator, 1 -i / p+M , also commutes with Σ. Therefore in the computation of the ghost propagator, we only need terms associated with one-particle irreducible diagrams which correspond at one loop to the diagram of figure 3.1. In Fourier space it takes the form The factors of the structure constants are equal to f aAc f Bac = -N δ AB . Furthermore, the last line can be improved observing that q ν P ⊥ µν (q -p) = (q ν -p ν )P ⊥ µν (q -p)+p ν P ⊥ µν (qp) = p ν P ⊥ µν (q -p) which has the advantage of having one internal momentum component less to compute in the integral. This trick will be used wherever possible. Therefore, the diagram of figure 3.1 represents the integral

d d q (2π) d ip µ gf aAc 1 q 2 iq ν gf Bac P ⊥ µν (q -p) (q -p) 2 + m 2 = -g 2 f aAc f Bac p µ d d q (2π) d q ν 1 q 2 P ⊥ µν (q -p) (q -p) 2 + m 2 . p p q q -p A B a c µ ν
N g 2 δ AB p µ p ν d d q (2π) d 1 q 2 P ⊥ µν (q -p) (q -p) 2 + m 2

Feynman trick

It is important to remark that, as some of the integral are divergent, we have to regularize the theory. In the previous example we have already regularized the theory using the dimensional regularization and we will continue using it so the integrals that we are going to present are well-defined.

The gluon propagator always comes with a factor

P ⊥ µν (q)
q 2 +m 2 which has a term like 1 q 2 1 q 2 +m 2 . The more fractions of the form 1 q 2 +m 2 we have, the harder the corresponding integrals are to compute. The first trick we are going to present allows us to reduce the number of propagators in each integral at the cost of increasing the number of integrals to be done. It consists in the substitution

1 q 2 1 q 2 + m 2 = 1 m 2 1 q 2 - 1 q 2 + m 2 . (3.2)
This trick will be implemented as many times as needed to reduce the number of propagator in the integrals. As a consequence the number of Feynman parameters (presented below) needed will be reduced and the integrals will be easier to compute.

The integral

d d q (2π) d 1 (q + p) 2 P ⊥ µν (q) q 2 + m 2
which is equivalent to the integral of Eq. (3.1) will be separated in

d d q (2π) d 1 (q + p) 2 δ µν q 2 + m 2 - 1 m 2 q µ q ν q 2 + 1 m 2 q µ q ν q 2 + m 2 .
It will be conveniently expressed in terms of

I 1 = d d q (2π) d 1 (q + p) 2 δ µν q 2 + m 2 , I 2 = - 1 m 2 d d q (2π) d 1 (q + p) 2 q µ q ν q 2 and I 3 = 1 m 2 d d q (2π) d 1 (q + p) 2 q µ q ν q 2 + m 2
The denominators can be combined using the Feynman trick, see e.g [START_REF] Peskin | An Introduction to quantum field theory[END_REF], which will simplify the computation of the angular integral. Feynman trick consists in using x 1 ,...,x n , Feynman parameters, one for each different propagator and making the substitution

1 D ν 1 1 D ν 2 2 ...D νn n = 1 0 dx 1 ...dx n δ( x i -1) x ν i -1 i (x 1 D 1 + ...x n D n ) ν i Γ(ν 1 + ... + ν n ) Γ(ν 1 )...Γ(ν n ) (3.3)
where D i denotes the ith-propagator and Γ in this case is the Euler Gamma function.

To compute the three integrals mentioned earlier we only need one Feynman parameter, x. If we apply the Feynman trick to the denominator of I 1 we have

1 (q + p) 2 (q 2 + m 2 ) = 1 0 dx 1 (x(q + p) 2 + (1 -x)(q 2 + m 2 )) 2 = 1 0 dx 1 (q 2 + 2xq.p + xp 2 + (1 -x)m 2 ) 2 = 1 0 dx 1 (l 2 + ∆) 2
where in the last line we change of variable l = q + xp and ∆ = x(1 -x)p 2 + (1 -x)m 2 does not depend on the integration momentum variable.

The first integral I 1 in dimensional regularization is given by

I 1 = 1 0 dx d d l (2π) d δ µν (l 2 + ∆) 2 = 1 0 dx δ µν Γ(2 -d/2) (4π) d/2 ∆ d/2-2 .
Applying the same trick to the third integral I 3 we obtain

I 3 = 1 m 2 1 0 dx d d l (2π) d (l -xp) µ (l -xp) ν (l 2 + ∆) 2 = 1 m 2 1 0 dx d d l (2π) d l µ l ν + x 2 p µ pν (l 2 + ∆) 2
where the odd power of l were eliminated because they integrate to zero. Moreover, integrals with two power of l can be substituted by

d d l (2π) d l µ l ν f (l 2 ) = δ µν d d d l (2π) d l 2 f (l 2 ).
Including this into the calculation of I 3 we have

I 3 = 1 m 2 1 0 dx δ µν d d d l (2π) d l 2 (l 2 + ∆) 2 + x 2 p µ p ν d d l (2π) d 1 (l 2 + ∆) 2 .
The integrals over l are known and they can be found in e.g. [START_REF] Bellac | Quantum and statistical field theory[END_REF]. We are interested in the value of

d d l (2π) d l 2 (l 2 + ∆) 2 = d 2 Γ(1 -d/2) (4π) d/2 ∆ d/2-1 .
Putting this result in I 3 we obtain

I 3 = 1 m 2 1 (4π) d/2 1 0 dx δ µν 2 ∆ d/2-1 Γ(1 -d/2) + x 2 p µ p ν ∆ d/2-2 Γ(2 -d/2) = 1 m 2 Γ(2 -d/2) (4π) d/2 1 0 dx δ µν (2 -d) ∆ + x 2 p µ p ν ∆ d/2-2 .
The second integral I 2 has a similar form

I 2 = - 1 m 2 Γ(2 -d/2) (4π) d/2 1 0 dx δ µν (2 -d) ∆(m = 0) + x 2 p µ p ν ∆(m = 0) d/2-2 .
To conclude, the steps we are going to use for computing a Feynman diagram are:

• Apply the substitution of Eq. (3.2) in order to decrease the number of propagators in each integral.

• Expand the sum and perform the Feynman trick.

• Ignore the integrals with an odd number of l in the numerator by using parity symmetry.

• Perform the integrals of the form

d d l (2π) d l n (l 2 +∆) 2 .
All these tasks are done by the Mathematica algorithm. The result is an integral over the Feynman parameters that in some cases can be computed analytically in integer dimensions. In other cases, this integral can not be calculated analytically so the integration is done numerically.

To conclude Eq. (3.1) is given by

N g 2 δ AB p µ p ν (I 1 + I 2 + I 3 ) = N g 2 δ AB p µ p ν Γ(2 -d/2) (4π) d/2 1 0 ∆ d/2-2 δ µν + 1 m 2 ( δ µν (2 -d) ∆ + x 2 p µ p ν ) -∆ d/2-2 1 m 2 δ µν 2 -d ∆ + x 2 p µ p ν = N g 2 δ AB Γ(2 -d/2) (4π) d/2 p 2 m 2 1 0 (1 -x)m 2 + x(1 -x)p 2 d/2-2 × m 2 + p 2 x 2 + (1 -x)m 2 + x(1 -x)p 2 2 -d -x(1 -x)p 2 d/2-2 p 2 x 2 + x(1 -x)p 2 2 -d .
The bare one-loop inverse of the ghost propagator is given by Γ

(2)

c a cb (p) = δ ab p 2 -N g 2 δ ab Γ(2 -d/2) (4π) d/2 p 2 m 2 1 0 dx (1 -x)m 2 + x(1 -x)p 2 d/2-2 × m 2 + p 2 x 2 + (1 -x)m 2 + x(1 -x)p 2 2 -d -x(1 -x)p 2 d/2-2 p 2 x 2 + x(1 -x)p 2 2 -d .
This integral can be done analytically for integer dimension (see for some examples Sect.3.3).

Gluon propagator

Let us consider now the self energy for gluons. The diagrams which contribute at one loop are the three shown in figure 3.2. In order to calculate those diagrams we are going to follow the same procedure explained for the ghost propagator. The computation of these diagrams by hand using the Feynman parameters are cumbersome but they were done in order to compare the results with the ones given by the Mathematica algorithm. The calculation is not particularly enlightening and we do not give the details here.

D AA 1 = - 1 2 g 2 d d q (2π) d 1 q 2 + m 2 P ⊥ ρσ (q) f eab f ecc (δ µρ δ νσ -δ µσ δ νρ ) + f eac f ebc (δ µν δ ρσ -δ µσ δ νρ ) +f eac f ecb (δ µρ δ νσ -δ µν δ ρσ ) = - 1 2 g 2 f eac f ebc d d q (2π) d 1 q 2 + m 2 P ⊥ ρσ (q) {2δ µν δ ρσ -δ µσ δ νρ -δ µρ δ νσ } = - 1 2 g 2 N δ ab d d q (2π) d 1 q 2 + m 2 2δ µν (d -1) -2P ⊥ µν (q) = -N g 2 δ ab d d q (2π) d 1 q 2 + m 2 δ µν (d -1) -P ⊥ µν (q) = -g 2 N δ ab (4π) d/2 (d -1) 2 d m d-2 Γ(1 -d/2)δ µν
where the last integral was done by standard textbook formulas.

The next diagram that we are going to compute is more interesting and is represented in figure 3.4. There is a factor 1/2 coming from the expansion of the exponential and a factor 2 from Wick's theorem. On the other hand, the ghost loop brings a minus sign. Feynman rules lead to the expression D AA 2 given by

D AA 2 = - d d q (2π) d igf acd q µ 1 (q -p) 2 igf bdc (q ν -p ν ) 1 q 2 = -g 2 N δ ab d d q (2π) d q µ (q -p) ν (q -p) 2 q 2 = g 2 N p 2 d/2-2 (d -2)(4π) d/2 Γ 2 -d 2 Γ d 2 2 Γ(d) (d -2)p µ p ν + p 2 δ µν . (3.4)
The last diagram is the most complex due to the presence of two three-gluon vertices. It is represented in figure 3

D AA 3 = 1 2 d d q (2π) d if acd [(2q + p) µ δ σρ + (p -q) σ δ µρ -(2p + q) ρ δ µσ ] P ⊥ σ (p + q) (q + p) 2 + m 2 × if bdc [(2q + p) ν δ η + (p -q) δ νη -(2p + q) η δ ν ] P ⊥ ηρ (q) q 2 + m 2 = 1 2 f acd f bdc d d q (2π) d i [(2q + p) µ δ σρ + (2p -(p + q)) σ δ µρ -((p + q) + p) ρ δ µσ ] P ⊥ σ (p + q) (q + p) 2 + m 2 × i [(2q + p) ν δ η + (2p -(p + q)) δ νη -(2p + q) η δ ν ] P ⊥ ηρ (q) q 2 + m 2 .
Eliminating the trivial projections of the form q ν P ⊥ µν (q) we have

D AA 3 = N 2 d d q (2π) d [(2q + p) µ δ ρσ + 2p σ δ µρ -2p ρ δ µσ ] [(2q + p) ν δ η + 2p δ νη -2p η δ ν ] × P ⊥ σ (p + q) (q + p) 2 + m 2 P ⊥ ηρ (q) q 2 + m 2 .
After doing the Feynman trick and performing the integral over momentum space the third diagram reads

D AA 3 = g 2 N Γ 2 -d 2 2(4π) d/2 1 0 dx   4p 2 ∆ d/2 0 + ∆ d/2 1 -2∆ d/2 2 (d -2)dm 4 + 4p 2 ∆ d/2-2 1 2 + p 2 /m 2 2x 2 -2x + 1 -8∆ 2 1 p 4 /m 2 ∆ d/2 2 + (2 -d) (d -1)∆ d 2 -1 1 - p 4 m 4 2x 2 -2x + 1 ∆ d/2-1 2 -(x -1) 2 ∆ d/2-1 1 + p 4 m 4 (x -1) 2 ∆ d/2-1 0 + 3p 2 m 2 ∆ d 2 -1 1 -∆ d 2 -1 2 δ µν + g 2 N Γ 2 -d 2 8(4π) d/2 m 4 1 0 dx   32 ∆ d/2 0 + ∆ d/2 1 -2∆ d/2 2 (d -2)d - (d -2) x -1 -(x -1) ∆ d 2 -1 1 8m 2 + p 2 -20x 2 + 4x + 3 -2∆ d 2 -1 2 4m 2 -5p 2 (1 -2x) 2 -3p 2 (x -1)∆ d 2 -1 0 -4(5x -1)∆ d/2 0 + 4 1 ∆ 2 1 ∆ 2 2 ∆ 2 3 x 2 ∆ 2 2 ∆ 2 3 x 2 ∆ d/2 1 m 4 d(1 -2x) 2 -4x 2 + 4x -9 +m 2 p 2 -4x 2 + 4x -3 + p 4 (1 -2x) 2 x 2 -∆ 2 1 p 2 ∆ d/2 2 × ∆ 2 2 m 2 1 -4x 2 + p 2 8x 4 -16x 3 + 14x 2 + 1 -2∆ 2 3 x 3 2m 2 + 3p 2 + (1 -2x) 2 ∆ d/2 0 (x -1) 2 p µ p ν (3.5)
where

∆ 0 = -p 2 (-1 + x)x, ∆ 1 = m 2 -p 2 (-1 + x)x, ∆ 2 = m 2 x -p 2 (-1 + x)x, ∆ 3 = m 2 -p 2 (-1 + x).
It is possible to perform the integral over the Feynman parameter x analytically in integer dimension but since the result is really long we will not write it here (for special dimensions see Sect.3.3). The results of all diagrams agree with those presented in [START_REF] Tissier | An Infrared Safe perturbative approach to Yang-Mills correlators[END_REF].

Both the ghost and gluon two-points vertex functions are divergent when d = 4. The divergent part of the ghost two-point vertex function is

Γ (2),1-loop c a cb (p) d→4 ∼ - 3 2 g 2 N δ ab p 2 16π 2 (4 -d)
and for the gluon

Γ (2),1-loop A a µ A b ν (p) d→4 ∼ g 2 N 16π 2 (4 -d) 3 2 m 2 - 13 3 p 2 δ ab P ⊥ µν (p). (3.6)
The divergences found are in agreement with previous calculations [START_REF] De Boer | On the renormalizability and unitarity of the Curci-Ferrari model for massive vector bosons[END_REF][START_REF] Gracey | Three loop MS-bar renormalization of the Curci-Ferrari model and the dimension two BRST invariant composite operator in QCD[END_REF]. The renormalization factors in MS scheme in dimensional regularization, that means the ones that do not modify the finite parts while absorbing the divergences, are

Z m 2 = 1 - 35 6 g 2 N 16π 2 1 4 -d , Z c = 1 + 3 2 g 2 N 16π 2 1 4 -d , Z A = 1 + 13 3 g 2 N 16π 2 1 4 -d .
(3.7)

Renormalization scheme and renormalization group

We have already computed the renormalization factors in the MS scheme. However, we prefer to introduce another scheme with more physical meaning which leads to different finite parts for the renormalization factors. Moreover, to estimate the dependence on the scheme we are going to work with two different schemes, the vanishing momentum (VM) and the infrared safe(IS), already presented in [START_REF] Tissier | An Infrared Safe perturbative approach to Yang-Mills correlators[END_REF]. Physical quantities are supposed to be scheme independent. However, as our calculations are only the leading order in perturbation theory a dependence is expected (hopefully small).

The relations between bare and renormalized vertices are the following:

Γ (2) A a µ A b ν (p) = Z A Γ (2),bare A a µ A b ν (p) Γ (2) c a cb (p) = Z c Γ (2),bare c a cb (p) Γ (3) c a cb A c µ (p, r) = Z c Z A Γ (3),bare c a cb A c µ (p, r) Γ (3) A a µ A b ν A c ρ (p, r) = Z 3/2 A Γ (3),bare A a µ A b ν A c ρ (p, r) (3.8)
For convenience in the definition of the renormalization scheme we are going to parametrize the two-point vertex function as

Γ (2) A a µ A b ν (p) = δ ab Γ ⊥ (p)P ⊥ µν + Γ (p)P µν , (3.9) Γ (2) c a cb (p) = δ ab p 2 J(p) (3.10)
It is important to note that only the transverse part of the gluon vertex function is related with the gluon propagator in Landau gauge,

A a µ A b ν (p) = δ ab P ⊥ µν (p) Γ ⊥ (p)
, and therefore the longitudinal part is not directly accessible by lattice simulations. However, some information about the longitudinal part can be obtained from Ward identities.

The relation for the renormalization of the coupling constant is based on the Taylor scheme [START_REF] Taylor | Ward Identities and Charge Renormalization of the Yang-Mills Field[END_REF]. The relation Z g √ Z A Z c = 1 is automatically fulfilled by the divergent parts of the renormalization factors to all orders of perturbation theory since no loop corrections appear when the ghost momentum is zero. The Taylor scheme extends this property also to the finite parts. The Taylor scheme used in both schemes makes the β-function of the coupling constant easy to compute in terms of the gluon and ghost anomalous dimensions, as

β g = µ ∂g ∂µ g B = µ ∂ g B Zg ∂µ = -µ g B Z g 1 Z g ∂Z g ∂µ = -g µ ∂ log Z g ∂µ = -g µ ∂ log( √ Z A Z c ) -1 ∂µ = g µ ∂(1/2 log(Z A ) + log(Z c )) ∂µ = g γ A 2 + γ c . (3.11)
The β-functions and the anomalous dimensions are important in order to know the evolution of the renormalized parameter with the scale. All the parameters will be expressed in term of the renormalized ones, which evolve following Eqs. (1.21). Once the initial condition are imposed at a scale µ 0 and taking the running scale µ ∼ p, the two-point vertex functions have the form

Γ (2) AA (p, µ 0 , g, m 2 ) = p 2 + m 2 (p) z A (p) ,
and

Γ (2) cc (p, µ 0 , g, m 2 ) = p 2 z c (p)
where z was defined in Eq. (1.25).

These last relations can be simplified using the precise choice of the renormalization scheme. We are going to separate the rest of the section in two, one for each scheme in order to give more details of the behaviour of each one.

The vanishing momentum scheme

The vanishing momentum (VM) scheme is defined as

Γ ⊥ (p = µ) = m 2 + µ 2 , J(p = µ) = 1, Γ ⊥ (p = 0) = m 2 and Z g Z A Z c = 1.
(3.12)

Unfortunately it presents a Landau pole, that means that the coupling constant goes to infinity at finite renormalization scale µ. The infrared Landau pole that appears in this scheme has no physical meaning. It is just an artificial pole which is consequence of the renormalization scheme chosen, which is not compatible with the fact that Γ ⊥ is not monotonic.

The evolution of the gluon mass is given by

β m 2 = µ ∂m 2 ∂µ g B ,m B = µ ∂ m 2 B Z m 2 ∂µ = -m 2 µ ∂ log Z m 2 ∂µ = -m 2 µ ∂ log Z -1 A ∂µ = m 2 γ A . (3.13)
where in the last line we have used that

∂ log Z m 2 ∂µ = ∂ log Z -1 A ∂µ
. This last relation is a consequence of the non-renormalization theorem for the mass presented in Appendix B.4. It can be deduced observing the equation (B.23):

Γ ⊥ B (0) J B (0) = m 2 B
and imposing the vanishing momentum condition as follows

Γ ⊥ B (0) J B (0) = Z -1 A Γ ⊥ R (0) J B (0) = Z -1 A m 2 R J B (0) = Z m 2 m 2 R therefore, Z -1 A J B (0) = Z m 2 .
If we take the logarithm we obtain

log Z -1 A -log J B (0) = log Z m 2 .
As the bare quantities do not depend on the renormalization scale µ, we can conclude that

∂ log Z m 2 ∂µ = ∂ log Z -1 A ∂µ
The relation (3.13) in addition to the Taylor scheme simplifies the factors z considerably, its expressions are given by

log z A (µ) = µ µ 0 γ A (µ ) µ dµ = µ µ 0 β m 2 (µ ) µ m 2 (µ ) dµ = µ µ 0 1 m 2 (µ ) ∂m 2 ∂µ dµ = m 2 (µ) m 2 (µ 0 ) dm 2 m 2 = log m 2 (µ) m 2 (µ 0 ) and log z c (µ) = µ µ 0 γ c (µ ) µ dµ = µ µ 0 β g (µ ) µ g(µ ) - 1 2
β m 2 (µ ) µ m 2 (µ ) dµ = log g(µ) g(µ 0 ) - 1 2 log m 2 (µ) m 2 (µ 0 )
where in the first line we have used that γ c = βg g -

β m 2
2m 2 deduced from combining Eq. (3.11) with Eq. (3.13). Therefore

z A (µ) = m 2 (µ) m 2 (µ 0 )
, and (3.14)

z c (µ) = g(µ) g(µ 0 ) m 2 (µ 0 ) m 2 (µ) . (3.15)
The infrared-safe scheme

The other scheme that we are going to consider has an advantage. For some initial conditions, it does not present a Landau pole so we can always chose the scale equal to the relevant momentum scale p. The infrared safe (IS) scheme is characterize by

Γ ⊥ (p = µ) = m 2 + µ 2 , J(p = µ) = 1, Z m 2 Z A Z c = 1. (3.16)
in addition to the Taylor scheme. The IS renormalization scheme differs from the VM scheme by replacing the constraint at zero momentum of the gluon two-point vertex function by

Z A Z c Z m 2 = 1. (3.17)
This relation, for the divergent parts, was first observed through perturbative calculations in [START_REF] Doria | Landau gauge and anti-BRS symmetry[END_REF][START_REF] Gracey | Three loop MS-bar renormalization of the Curci-Ferrari model and the dimension two BRST invariant composite operator in QCD[END_REF] and it can be seen as a consequence of the nonrenormalization theorem presented in [DGL + 05, Wsc08, TW11] and explained in Appendix B.4. The equation (B.24) shows that, as the renormalized quantities are finite, the divergent parts of the renormalization factors satisfy Z A Z c Z m 2 = 1. This relation is also automatically satisfied by the leading order finite parts. Therefore it is possible to choose to satisfy this identity, also for the finite parts, for all order of perturbation theory. This is the case of the IS scheme. In this case, the renormalized mass m 2 R does not correspond to the value of the renormalized propagator at zero momentum any more.

From the relation Z A Z c Z m 2 = 1, it is easy to relate the β m 2 -function with the anomalous dimension of the ghost and gluon propagators, having

β m 2 = -m 2 µ ∂ log Z m 2 ∂µ = -m 2 µ ∂ log(Z A Z c ) -1 ∂µ = m 2 µ ∂(log(Z A ) + log(Z c )) ∂µ = m 2 (γ A + γ c ) . (3.18)
In this scheme, the gluon anomalous dimension change with respect to the VM scheme while the ghost anomalous dimension remains identical at one loop. The expressions relating z A (µ) and z c (µ) to the coupling constants can be deduced using the relations

γ A = 2 - βg g + β m 2 m 2 and γ c = 2 βg g - β m 2
m 2 which were obtained by Eq. (3.11) and Eq. (3.18) as follows

log z A (µ) = µ µ 0 γ A (µ ) µ dµ = µ µ 0 {-2 β g (µ ) µ g(µ ) + 2 β m 2 (µ ) µ m 2 (µ ) }dµ = -2 log g(µ) g(µ 0 ) + 2 log m 2 (µ) m 2 (µ 0 ) and log z c (µ) = µ µ 0 γ c (µ ) µ dµ = µ µ 0 {2 β g (µ ) µ g(µ ) - 1 2
β m 2 (µ ) µ m 2 (µ ) }dµ = 2 log g(µ) g(µ 0 ) - 1 2 log m 2 (µ) m 2 (µ 0 ) therefore z A (µ) = m 4 (µ) m 4 (µ 0 ) g 2 (µ 0 ) g 2 (µ) and (3.19) z c (µ) = g 2 (µ) g 2 (µ 0 ) m 2 (µ 0 ) m 2 (µ) .
(3.20)

Results

Results in d = 4 imposing the VM-renormalization scheme

In general, two-point vertex functions containing only two types of propagators can be computed analytically in an integer dimensional space. Here, we present the results in four space dimension which is the dimension with the most physical interest. Imposing the vanishing momentum condition within the renormalization scheme we find that the renormalized Γ ⊥ and J(p) are finite and take the form

Γ ⊥ (p)/m 2 = s + 1 + g 2 N s 384π 2 111s -1 -2s -2 + (2 -s 2 ) log s + 2(s -1 + 1) 3 s 2 -10s + 1 log(1 + s) + (4s -1 + 1) 3/2 s 2 -20s + 12 log √ 4 + s - √ s √ 4 + s + √ s -(s → µ 2 /m 2 ) , (3.21) J -1 (p) = 1 + g 2 N 64π 2 -s log s + (s + 1) 3 s -2 log(s + 1) -s -1 -(s → µ 2 /m 2 ) , (3.22)
where s = p 2 /m 2 . These expressions can be compared with lattice simulations. In order to do so, we have to fix the value of the coupling constant and gluon mass at the renormaliztion scale µ 0 where the renormalization scheme was established. The fitting is done by using a multiplicative constant to relate our renormalization scheme with the unrenormalized lattice data. The fitting parameters were chosen to minimize simultaneously the squared errors defined as

χ 2 AA = 1 4N i (Γ ⊥ lt. (µ 0 ) 2 + Γ ⊥ lt. (p i ) 2 ) 1 Γ ⊥ lt. (p i ) - 1 Γ ⊥ th. (p i ) 2 χ 2 cc = 1 4N i (J -2 lt. (µ 0 ) + J -2 lt. (p i )) (J lt. (p i ) -J th. (p i )) 2 . (3.23)
This definition of the error corresponds to a sort of average between the absolute error and the relative error.

We expect that if we implement the RG flow the accuracy of the fit in the ultraviolet will increase allowing us to fit simultaneously the UV and IR regimes. In four dimensions, imposing the vanishing momentum scheme, the anomalous dimensions read as

γ A = - g 2 N 192π 2 t 3 t(34t 2 -175t + 6) -2t 5 log t + 2(t + 1) 2 (2t 3 -11t 2 + 20t -3) log(t + 1) + 2t 3/2 √ t + 4 t 3 -9t 2 + 20t -36 × log √ t + 4 - √ t √ t + 4 + √ t , (3.24)
and

γ c = - g 2 N 32π 2 t 2 2(t + 1)t -t 3 log t + (t + 1) 2 (t -2) log(t + 1) , ( 3.25)
where t = µ 2 /m 2 . From these expressions we can deduce the β-functions for the coupling constant and mass. The function β m 2 is obtained by multiplying the gluon anomalous dimension by the squared mass. The evolution of the coupling constant is given by Eq. (3.11) which at high momentum with respect to the gluon mass matches with the standard β-function. The ultraviolet and infrared limits are

β g ∼ -g 3 N 16π 2 11 3 if µ m, -g 3 N 16π 2 1 12 if µ m.
(3.26)

In the VM scheme β g remains negative for renormalization scale smaller than m. This property results in the existence of a Landau pole meaning that the coupling constant diverges at a finite renormalization scale.

However, as we established earlier there is no need for the renormalization group at small momentum in d = 4. At low momentum the large logarithm could only arrive from the ghost massless propagator since the gluon propagator is regularized in the infrared by the mass. In four dimensions this effect is not as important as in lower dimensions. Consequently, one could reproduce the results considering the evolution of the flow equation only for p ≥ m. This can be implemented by using a running scale of the form

µ = p 2 + αm 2 , (3.27)
where α is a positive parameter to be chosen. In this way, µ ∼ p for scales larger than the mass (where large logarithm are important) and µ freezes out for small momenta.

Results for the SU (2) gauge group

Let us start comparing our results of equations (3.21) and (3.22) without considering the effect of the renormalization group. Therefore we fix µ = µ 0 = 1 GeV and use the value of g and m that simultaneously minimize the errors defined in eq. (3.23). The best fits are obtained using the parameters g = 7.5, m = 0.68 GeV. The corresponding fits are shown in figures 3.6 for the gluon propagator and in 3.7 for the ghost dressing function. The results match the lattice data with surprisingly good accuracy. The lattice data does not reach large values of the momenta therefore we can not test the necessity of the renormalization group in the ultraviolet. In spite of this, we are going to consider the renormalization group effect and study the impact of the value α in the definition of the running scale (3.27). We show in figure 3.8 the results obtained with RG improvement for different values of α. The initial conditions of the flow equation, g 0 = g(µ 0 ) and m 0 = m(µ 0 ), were imposed at µ 0 = 1GeV and are presented in the table 3.1 for different values of α in SU (2). As the dependence on α is not so big we are going to work with α = 1 every time the VM scheme is implemented.

We want to stress that the parameter chosen as initial condition of the RG flow are not unique in the sense that there is a region in the space of parameters which produce fits with the same level of precision, i.e. less than 10%. This region for the case of α = 1 is shown in figure 3.9. It can be seen that the corresponding region for χ AA is narrower than the region for χ cc (which expands far beyond what is shown in the figure). We can conclude that the gluon propagator is more difficult to reproduce than the ghost dressing function.

The evolution of the running coupling constant and the running gluon mass are shown Some results for the SU (3) gauge group are presented in [START_REF] Tissier | An Infrared Safe perturbative approach to Yang-Mills correlators[END_REF].

Results in d = 4 imposing the IS-renormalization scheme

When we change the vanishing momentum condition to the one based on the non-renormalization of the mass, the ghost dressing function and its anomalous dimension remain equal to (3.22) and (3.25) respectively. The modification acts only in the gluon sector where the gluon propagator is

Γ ⊥ (p) m 2 = s + 1 + g 2 N 384π 2 2s -2 s 2 -10s + 1 (s + 1) 3 log(s + 1) -s s 2 -2 log(s)
-s -1/2 (s + 4)

√ t + √ t + 4 √ t + 4 - √ t -2t -3 (t + 1) 3 s t 2 -10t + 4 -3t log(t + 1) - 2s t 3 (t + 4) 7t 4 + 30t 3 -54t 2 -200t + 192 log 1 - t t + 4 -2s -1 -6t -1 + 96 (3.28)
where s = p 2 /m 2 and t = µ 2 /m 2 . For the renormalization group, the ghost anomalous dimension remains the same as in the previous scheme and the gluon anomalous dimension takes the form

γ A = g 2 N 96π 2 t 3 -(t -2) 2 (2t -3)(t + 1) 2 log(t + 1) + (-17t 2 + 74t -12)t + t 5 log(t) -t 3/2 √ t + 4 t 3 -9t 2 + 20t -36 log √ t + 4 - √ t √ t + 4 + √ t , (3.29)
Here, the ultraviolet behaviour (µ m) of the β g function is the same as in standard Yang-Mills, see (3.26), but for low renormalization scale it changes its sign:

β g ∼ g 3 N 16π 2 1 6 if µ m. (3.30)
Depending on the initial condition, this can prevent the flow from having a Landau pole as a function of the renormalization scale. Let us recall that if the initial condition of the gluon mass is set to zero we recover the standard Yang-Mills beta function and the infrared Landau pole is recovered. That is why depending on the initial conditions the flow can be infrared safe or not. When m is large enough the solution of Eq. (3.30) gives a coupling constant that goes to zero at low momenta. This property justifies the use of perturbation theory in the infrared range. It is worth mentioning that this is not in disagreement with the universality of the β-function since the scheme independence property of the two first orders of β-functions is valid only for mass-independent schemes or in the ultraviolet limit when the masses can be neglected.

Results for the SU (2) gauge group

We are going to start studying the results considering the effects of the renormalization group for α = 0. The fit of the analytical results with the lattice data is shown in figure 3.11. The initial conditions of the RG flow imposed at 1GeV were g = 5.2 and m = 0.44 GeV. These parameters were chosen to minimize the errors defined in Eq.

(3.23) demanding a compromise between improving the gluon fits with the ghost one.

The contour plot of these errors are presented in figure 3.12 where it can be appreciated that there is a region in the parameter space that lead to an error small than 10%. The region associated with the gluon propagator is narrower than the one associated with the ghost dressing function which continues beyond what is shown in the figure. This make us conclude that the gluon propagator is more exigent to fit than the ghost dressing function.

The evolution of the running coupling constant and the running gluon mass are shown in figure 3.13. The coupling constant does not reach a Landau pole but instead goes to zero again in the infrared. Similarly the running mass goes to zero in the infrared. This is not incompatible with the fact that the gluon propagator is finite at low momenta because, in the IS scheme, the mass is not defined through the value of the propagator at zero momenta as in the VM scheme.

Although this scheme is infrared safe and it is not necessary to consider a positive value of α in (3.27), we study the effect of changing α in figure 3 Henceforth when working with the IS scheme we are going to take α = 0. This choice is based on the curves that we will show for d = 3 where to reproduce the gluon propagator the RG flow has to reach the deep infrared. Figure 3.15 -Contour plots for the error functions χ AA for the gluon propagator (elliptic region, for 10%, 15% and 20%) and χ cc for the ghost propagator (diagonal region, for 4%, 7% and 10%).

Results for the SU (3) gauge group

We complete our study in four dimensions presenting the results of the gluon propagator and the ghost dressing function when the IS-renormalization scheme is done and when the group effects are taking into account. In this case we take α = 0 and therefore the running scale µ = p. With this purpose we analyse the corresponding errors defined in (3.23). The contour lines of them are shown in figure 3.15.

We choose the parameters in the intersection of the darkest region of 3.15 which correspond to consider g 0 = 4.1 and m 0 = 0.43 GeV at µ 0 = 1 GeV. The fits for the gluon propagator and the ghost dressing function with the lattice data of [START_REF] Bogolubsky | Lattice gluodynamics computation of Landau gauge Green's functions in the deep infrared[END_REF] are shown in figure 3.16. It can be observe that both set of lattice data are well reproduce with our one-loop results.

Results in d = 3

In this part we are going to present the results for the three dimensional case. The twopoint vertex functions can be written as

Γ (2 c a cb (p) = δ ab p 2 + g 2 mN 32π √ s πs 2 + 2 √ s(1 -s) -2(s + 1) 2 arctan √ s , ( 3
Γ ⊥ (p) = p 2 + m 2 + g 2 mN 128πs 3/2 -4 5s 2 + 7s -1 √ s + π s 2 -2 s 2 -4(s + 1) 2 s 2 -6s + 1 arctan √ s + 2s(s + 4) s 2 -12s + 8 arctan √ s 2 , (3.32)
where we used the same notation introduced in the previous sections. It is worth noting that in three dimensions the coupling constant has dimensions of √ GeV. We are also going to renormalize the theory with the two schemes discussed above. As these diagrams do not diverge the only consequence of the renormalization scheme is to modify the finite part of them.

In figure 3.17 and 3.18 we compare our next to leading order expression with the d = 3 lattice simulations for SU (2) (without including the RG effects). The quality is much lower than in four dimensions. The fitting that match better seems to be the VM scheme although the peak in the lattice data is not reproduced. Both schemes were done using g = 3.7 √ GeV, m = 0.89 GeV and µ = 1 GeV. .17 -Three dimensional renormalize ghost dressing function for SU (2). In red the results of Eq. (3.31) without including the effects of the RG using g = 3.7 √ GeV, m = 0.89 GeV and µ = 1 GeV. The gray dots are the lattice data from [START_REF] Cucchieri | Constraints on the IR behavior of the gluon propagator in Yang-Mills theories[END_REF] We observed that no set of parameters reproduces the gluon's pronounced growth in the IR range even though all of them give a finite contribution at low momenta. The renormalization group approach is important in the infrared for d < 4 since the radiative corrections become larger in the infrared due to the presence of the massless ghosts at lower dimensions.

In figure 3.19 we show how the RG can considerably improve the results at low momenta in d = 3. Moreover the enhancement of the gluon propagator in the infrared is reproduced using the IS scheme and the evolution of the corresponding RG equation. The initial condition of the IS flow g 0 = 2.4 √ GeV and m 0 = 0.55 GeV were fixed at µ 0 = 1 GeV. On the other hand for the vanishing momentum flow the curve was done using g 0 = 4 √ GeV and m 0 = 1 GeV at 1 GeV. As this scheme has an infrared Landau pole, we used α = 1 in the running scale µ = p 2 + αm 2 in order to freeze the flow at scales lower than the mass.

We study the dependence of the VM scheme on α. For α = 0.5, 1.0, 2.0, 3.0 we determine the best fitting parameters, however, none of them help to reproduce the pronounced gluon peak, see figure 3.20. The same analysis was done for IS scheme even though the inclusion of a positive α is not needed since this scheme does not present a Landau pole. We observe in figure 3.21 that a RG evolution for the IR is important in order to reproduce the enhancement of the gluon propagator, therefore the flow is required to run into the deep infrared. That is probably why when using the RG for the IS scheme the parameter α is set to zero. The fitting parameters used in both cases are summarized in table 3.3. For now on, we will focus on α = 1 for the study of the VM scheme and α = 0 for the IS scheme. We want to stress that the choice of parameters is not unique meaning that there is a region in the parameter space which gives the same accuracy in the fits. This effect can be seen by studying the contour level of the errors defined in Eq. (3.23). The left hand side figure of 3.22 shows that we are able to chose any pair of parameters belonging to the intersection of the contours in order to have an error less that 10%. However, for the (VM) case the right hand side of figure 3.22 has no intersection so we choose a set of parameters near both regions. This is in agreement with the fact that the IS scheme work much better in d = 3. We do not have an explanation why in d = 4 case the preferred renormalization is VM while when d = 3 it is IS.

Estimation of higher loop corrections

The values obtained for the coupling constant could possibly make us doubt on the validity of perturbation theory. Let us analyse the value of the coupling constant in figure 3.23 in the IS scheme. The value of g goes to zero at vanishing momentum but it takes large values around µ ∼ 1GeV. This situation is better than having a Landau pole but it may invalidate perturbation theory. The following discussion was obtained from [START_REF] Tissier | An Infrared Safe perturbative approach to Yang-Mills correlators[END_REF] and allows us to justify the use of perturbation theory even at momenta of order of one GeV.

First of all it is worth noting that the real expansion parameter is not g 2 but it is

u(p) = g 2 (p)N p d-4 (4π) d/2 Γ(d/2)
which is dimensionless for any value of the space dimension and includes N and the angular factor appearing in the loop integral. The full line in Figure 3.24 shows how the angular factors reduce the expansion parameter with respect to the value of the coupling constant. Moreover, there is actually another effect that suppresses the corrections in the infrared. As the gluons are massive, for momenta lower than their mass all the internal gluon propagators are suppressed by the inverse of the squared mass and therefore the internal gluon freezes out at low momenta. So we can think of an effective theory with only a dynamical massless ghost. But as the ghost interacts through the exchange of gluons all the vertices of the form shown in figure 3.25 behaves as p 1 p 2 g 2 /m 2 , where p 1 and p 2 are the corresponding ghost momenta. Consequently, the effective expansion parameter is suppressed by p 2 /m 2 where p is a typical external momentum. In order to estimate this effect for arbitrary momenta we can consider an interpolation of the ultraviolet expansion parameter u(p) with the effective infrared parameter u(p)p 2 /m 2 as performing the expansion with respect to u(p) p 2 p 2 + m 2 which is shown in figure 3.24 with dashed lines. This new parameter is small enough to do perturbation theory even in the region where the coupling constant is not so small. It is important to mention that this expansion parameter at low momentum is not as small as its value in the ultraviolet or in comparison with the one used in QED where perturbation theory has been completely tested. In our case, the two loop contribution seems to be of order 0.4 2 ∼ (15%) which is approximately the error that we will have with the three point correlations functions computed in next chapter.

Quenched three-point vertex functions

In the preceding chapter we made the first steps towards investigating the theory based on the Lagrangian (2.11). The one-loop calculations done for the gluon propagator and the ghost dressing function compare very well with lattice simulations both in d = 4 and d = 3 [TW11, PTW13] at least in some schemes. It is surprising that such accuracy can be obtained with a modest one-loop calculation. In the previous chapter we reproduced the analysis of the estimation of two-loops corrections done in [START_REF] Tissier | An Infrared Safe perturbative approach to Yang-Mills correlators[END_REF]. There, it was found that higher loop corrections to the propagators are suppressed by powers of the external momenta in the infrared regime. This fact supports the idea that perturbation theory in this massive version of QCD in Landau gauge is under control.

In this chapter we want to go further in the investigation of one loop corrections introduced by the Lagrangian (2.11). This time we are going to study more complex quantities, such as the three-point ghost-gluon correlation function and the three-gluon correlation function. The analysis of these vertices is extremely interesting since their tensor structure is not as simple as for the propagators. They depend on two independent momenta or, equivalently on three independent scalar variables. Therefore, the computation of the three-point vertices is a non-trivial test to validate our model. Moreover, the value of the coupling constant and gluon mass at 1 GeV is fixed by the analysis of the previous chapter, and, therefore the three-point vertices are completely determined and we do not have any free parameter to fix.

In the past few years, several techniques were used with the aim of describing these three-point correlation functions. For instance, semi-analytical methods such as truncated Schwinger-Dyson equation or Non-perturbative RG [HvS13, RQ11b, DORQ12, AIP13] were employed. Some studies [ABIP14b, EWAV14, VAEW14] specially interested in the zero-crossing of the three-gluon vertex appeared after our results were published. Fortunately several simulations on the lattice were also done, see [Par94, AHP + 97, CMM04, IMPSS06, CMM06, CMM08]. However, as the vertex functions are more complex quantities, lattice data are not as precise as for the propagators.

In this chapter, we present the Feynman diagrams and their evaluation for the ghostantighost-gluon vertex and the three-gluon vertex. The explicit calculation was performed by implementing a Mathematica algorithm. We performed the calculation in two different ways. The results for the one-loop expressions using the method of Feynman parameters are too large to be presented here. Instead we used Passarino-Veltman [START_REF] Passarino | One Loop Corrections for e+ e-Annihilation Into mu+ mu-in the Weinberg Model[END_REF] integrals (presented in Appendix C) to present them in a more compact way within the complementary material of [START_REF] Pelaez | Three-point correlation functions in Yang-Mills theory[END_REF]. Some checks were done on the calculations by com-paring for example the zero mass limit with the expressions on standard Yang-Mills of [START_REF] Davydychev | Three gluon vertex in arbitrary gauge and dimension[END_REF][START_REF] Davydychev | One loop results for three gluon vertex in arbitrary gauge and dimension[END_REF]. Moreover, we deduced some Slavnov-Taylor identities which allow us to do an independent check of the expressions. Both algorithms, the one performing the Feynman Trick and the one which reduces the expressions to Passarino-Veltman integral fulfil the mentioned checks, and, furthermore, they give the same results.

Towards the end of the chapter we implement the renormalization group approach in order to reproduce the lattice data given by [START_REF] Cucchieri | Three-point vertices in Landaugauge Yang-Mills theory[END_REF]. As we already mentioned, the comparison was done using the same initial condition obtained in the previous chapter for the two-point correlation functions. In this sense, the calculation of the three-point correlation functions becomes a pure prediction without any parameter to adjust. The results match the three different kinematic configurations provided by lattice simulations very well.

Ghost-gluon vertex

In this section we are going to introduce the one loop contribution to the ghost-gluon vertex. All the calculations presented in this section are in terms of bare quantities even though it is not specified. The renormalization procedure is deferred to section 4.4. The one loop contribution to the ghost-gluon vertex function is diagrammatically represented by the two diagrams of figure 4.1. These diagrams can be computed by using the methods described in the previous chapter. The corresponding expressions in Landau gauge, in arbitrary dimension d and number of colors N is described in the next paragraphs.

Our convention of momenta and color indices are elucidated in fig 4 .2. The expression of this diagram, D1 ccA , can be obtained following the Feynman rules described in Sec.1.3, and it reads

D1 ccA = -(ig) 3 d d q (2π) d f lmc p ν 1 (q -p) 2 f anm (q -p) µ 1 (q + r) 2 f lbn (q + r) ρ P ⊥ νρ (q) q 2 = i N 2 g 3 f abc p ν r ρ d d q (2π) d (q -p) µ 1 (q -p) 2 1 (q + r) 2 P ⊥ νρ (q) q 2
where in the last line we used that f bac = -i(t a G ) bc where t a G is the generator in the adjoint representation and therefore

f anm f lbn f lmc = f man f nbl f lcm = iTr(t a G t b G t c G ) = -N 2 f
abc and where we also cancel trivial contractions of the type q ρ P ⊥ νρ (q). The second digram is a bit longer and it is represented in figure 4.3. The Feynman rules leads us to the expression

D2 ccA = -(ig) 3 d d q (2π) d f amn [(q + p -(r -q)) µ δ νρ + (r -q + p + r) ν δ µρ + (-(p + r) -(q + p)) ρ δ µν ] × P ⊥ νσ (q + p) (q + p) 2 + m 2 P ⊥ ρη (q -r) (q -r) 2 + m 2 f lcm (-p σ ) 1 q 2 f bln q η = -ig 3 f amn f bln f lcm p σ d d q (2π) d [(2q + p -r) µ δ νρ + (-q + p + 2r) ν δ µρ + (-q -r -2p) ρ δ µν ] q η × 1 q 2 P ⊥ νσ (q + p) (q + p) 2 + m 2 P ⊥ ρη (q -r) (q -r) 2 + m 2 = -i N 2 g 3 f abc p σ r η d d q (2π) d [(2q + p -r) µ δ νρ + 2(p + r) ν δ µρ -2(p + r) ρ δ µρ ] × 1 q 2 P ⊥ νσ (q + p) (q + p) 2 + m 2 P ⊥ ρη (q -r) (q -r) 2 + m 2
Therefore we can compute the ghost-gluon vertex or even better a richer tensor Γ νµ .

Our analysis of the three-point ghost-gluon correlation function was split in two. First, we calculated directly the one-loop contributions to the three-point vertex. Independently we performed the computation of Γ (3) Aµc Kν is -gf abc δ µν ). This allowed us to compute all the scalar functions that we present in the suplemental material of [START_REF] Pelaez | Three-point correlation functions in Yang-Mills theory[END_REF]. The relation between both independent calculations was checked using the identity (4.4). The ghost-gluon vertex is found to have no ultraviolet divergence for d = 4 in Landau gauge in agreement with the non-renormalization theorem of the coupling constant. It is important to stress that this does not mean that g is unrenormalized but that its renormalization is fixed in terms of the renormalization field factors Z A and Z c . In particular the one loop contribution is zero when the ghost or the anti-ghost momentum vanishes. The limit of vanishing gluon momentum is easy to compute analytically and can be expressed in d = 4 as

q p ν, c b r q + r q -p k = -(p + r) µ, a m n ν ρ l q p ν, c b r r -q q + p k = -(p + r) µ,
Γ µν (-p, p,0) = δ µν 1 + g 2 0 N 128π 2 9/2 + s + 5s -1 -(7s -1 + 5s -2 ) log(s + 1) -(s -1)s log s -1 + 1 (4.6)
were we recall that s = p 2 /m 2 .

In three dimensions this tensor reads

Γ µν (-p, p,0) = δ µν 1 + g 2 0 N 384πms 2(6s 2 -5s -21) -3π √ s(2s 2 -s + 288) + 6s -1/2 (2s 3 -s 2 -68s + 7) arctan( √ s) (4.7)
These expressions are regular even in the infrared for p m. In fact, for arbitrary momentum configuration if we multiply all the momenta by a coefficient κ we have that

Γ µν ({κp i }) -δ µν ∼ κ d-2
which is regular in the infrared for d > 2.

Three-gluon vertex

In this section we present the calculations of the three-gluon vertex at one loop. There are three diagrams that contribute to this vertex shown in figure 4.5. The first diagram, D1 AAA , described in detail in figure 4.6 has a loop of ghosts which contributes with a minus sign. Taking into account the minus coming from the expansion of the exponential and the combinatorial factor the global symmetry factor is one, so no factor will precede the integral. We have to take into account that there are two possible directions for the ghost loop. For convenience, we are going to compute separately the diagram with the anti-clockwise loop, D1

[1] AAA , and the clockwise loop,D1

[2] AAA . Therefore, the expression of those diagrams can be computed as

D1 [1] AAA = d d q (2π) d (igf alm q µ ) 1 (q -p) 2 (igf cnl (q -p) ρ ) 1 q + r (igf bmn (q + r) ν ) 1 q 2 = -ig 3 N 2 f abc d d q (2π) d q µ (q -p) ρ (q + r) ν 1 q 2 1 (q -p) 2 1 (q + r) 2 (4.8)
and

D1 [2] AAA = ig 3 N 2 f abc d d q (2π) d (q + p) µ (q -r) ρ q ν 1 q 2 1 (q + p) 2 1 (q -r) 2 = -ig 3 N 2 f abc d d q (2π) d (q -p) µ (q + r) ρ q ν 1 q 2 1 (q -p) 2 1 (q + r) 2 (4.9)
It is worth noting that this is the only infrared singular diagram that we have since all internal lines correspond to ghost propagators (which are massless). The other diagrams only have internal gluon propagators which are regular in the infrared.

The second diagram of figure 4.5 is denoted as D2 AAA . We have to consider the possible contribution from the three different external momentum configurations. That is made explicit in figure 4.7, where the left diagram is called D2

[1]
AAA , the middle one is D2 AAA , in the middle D2

[2]

AAA , and in the right D2

[3] AAA .

D2

[1]

AAA = 1 2 d d q (2π) d igf mlb (q -r) β δ αν + (2r + q) α δ βν + (-2q -r) ν δ α β P ⊥ γα (q) q 2 + m 2 P ⊥ δβ (q + r) (q + r) 2 + m 2 × g 2 f eac f eml (δ µδ δ ργ -δ µγ δ ρδ ) + f eac f ecl (δ µρ δ γδ -δ µγ δ µγ δ ρδ ) + f eal f emc (δ µδ δ ργ -δ µρ δ γδ ) = 3 2 ig 3 f abc (δ µν δ ργ -δ µγ δ ρδ ) d d q (2π) d (r β δ αν -r α + q ν δ αβ ) P ⊥ γα (q) q 2 + m 2 P ⊥ δβ (q + r) (q + r) 2 + m 2
where in the last line we have simplified the momentum contracted with the perpendicular projector. The other diagrams are obtained by permutation of the external legs and take the form

D2 [2] AAA = 3 2 ig 3 f abc (δ µδ δ νγ -δ µγ δ νδ ) d d q (2π) d ((p + r) β δ αρ -(p + r) α δ βρ -q ρ δ αβ ) × P ⊥ γα (q) q 2 + m 2 P ⊥ δβ (q -p -r) (q -p -r) 2 + m 2 D2 [3] AAA = 3 2 ig 3 f abc (δ νδ δ ργ -δ νγ δ ρδ ) d d q (2π) d (-p β δ αµ + p α δ βµ -q µ δ αβ ) P ⊥ γα (q) q 2 + m 2 P ⊥ δβ (q + p) (q + p) 2 + m 2 .
The last diagram that we have to consider is formed by three three-gluon vertices. This is by far the most difficult diagram to compute due to the complex structure of the three-gluon vertex in addition to the three gluon propagators. It is presented in figure 4

D3 AAA = - d d q (2π) d f aml [(-2q + p) µ δ αβ + (q -2p) β δ µα + (p + q) α δ µβ ] P ⊥ γβ (q) q 2 + m 2 × igf nmb [(q -r) δ δ νγ + (2r + q) γ δ νδ + (-2q -r) ν δ γδ ] P ⊥ δη (q + r) (q + r) 2 + m 2 P ⊥ λα (q -p) (q -p) 2 + m 2 × igf lnc [(q + p + 2r) λ δ ρη + (q -r -2p) η δ λρ + (p -2q -r) ρ δ λη ] = -i4N g 3 f abc d d q (2π) d [p α δ µβ -q µ δ αβ -p β δ µα ] [r γ δ δν -r δ δ γν -q ν δ γδ ] × [(p + r) λ δ ρη -(p + r) η δ λρ -(q + r) ρ δ λη ] P ⊥ γβ (q) q 2 + m 2 P ⊥ δη (q + r) (q + r) 2 + m 2 P ⊥ λα (q -p) (q -p) 2 + m 2
The unrenormalized three-gluon vertex function is obtained, at one loop, by summing all these diagrams plus the contribution at tree level. That is means that

Γ (3) A a µ A b ν A c ρ (p, r, k) = -igf abc {(r -k) µ δ νρ + (k -p) ν δ µρ + (p -r) ρ δ µν } + D1 [1] AAA + D1 [2] AAA + D2 [1] AAA + D2 [2] AAA + D2 [3] AAA + D3 AAA (4.10)
As for the ghost-gluon vertex, we ignore that for a general gauge group, contributions at higher loop have more complicate color structure involving, for instance, the completely symmetric color structure d abc . However, for SU (2) or N = ∞ as well as for the one-and two-loop corrections the color structure is f abc and can be factorized as

Γ (3) A a µ A b ν A c ρ (p, r, k) = -igf abc Γ µνρ (p, r, k) (4.11)
(for the two-loop color structure see [START_REF] Davydychev | Two loop three gluon vertex in zero momentum limit[END_REF]). This vertex has a richer momentum structure due to the three different Lorentz indices. In fact, it involves six independent scalar components. As suggested by Ball and Chiu in [START_REF] Ball | Analytic Properties of the Vertex Function in Gauge Theories. 2[END_REF], it can be expressed in terms of scalar functions as

Γ µνρ (p, r, k) = A(p 2 , r 2 , k 2 )δ µν (p -r) ρ + B(p 2 , r 2 , k 2 )δ µν (p + r) ρ -C(p 2 , r 2 , k 2 )(δ µν p.r -p ν r µ )(p -r) ρ + 1 3 S(p 2 , r 2 , k 2 )(p ρ r µ k ν + p ν r ρ k µ ) + F (p 2 , r 2 , k 2 )(δ µν p.r -p ν r µ )(p ρ r.k -r ρ p.k) + H(p 2 , r 2 , k 2 ) -δ µν (p ρ r.k -r ρ p.k) + 1 3 (p ρ r µ k ν -p ν r ρ k µ ) + cyclic permutations.
(4.12)

The three-gluon vertex is symmetric with respect to any interchange of the external legs. That is the reason why the scalar functions A, B, C, S, F, H have extra symmetries. A, C and F are symmetric under a permutation of the first two arguments; B is antisymmetric under a permutation of the first two arguments while H is completely symmetric and S completely antisymmetric. In the three-gluon correlation function in Landau gauge neither B nor S will appear since the vertex is contracted with the transverse gluon propagator which makes the longitudinal part of the vertex disappear.

The full expression for the scalar functions is extremely large and will be not presented here. It can be found together with the previous vertex in the supplemental material of [START_REF] Pelaez | Three-point correlation functions in Yang-Mills theory[END_REF]. For general momenta it requires to perform the integrals numerically.

The expression considerably simplifies when we set one gluon momentum to zero. For

d = 4 -it reads Γ µνρ (p, 0, -p) = 1 - N g 2 0 768π 2 - 136 (1 -log m) + 1 3 (36s -2 -594s -1 + 319 + 6s) + (3s 2 -2) log s -4s -3 (1 + s) 3 (s 2 -9s + 3) log(1 + s) + (4 + s) 3/2 s 3/2 24 -30s + s 2 log √ 4 + s + √ s √ 4 + s - √ s (p µ δ νρ + p ρ δ µν ) -2 + N g 2 0 384π 2 - 136 (1 -log m) + 1 3 (18s -2 -321s -1 -97 + 24s) + (s -1)(s 2 -2s -2) log s -2s -3 (1 + s) 2 (s -1)(s 3 -7s 2 + 7s -3) log(1 + s) + √ 4 + s s 3/2 48 + 16s + 22s 2 -11s 3 + s 4 log √ 4 + s + √ s √ 4 + s - √ s p ν δ µρ - N g 2 0 384π 2 m 2 (-36s -3 + 278s -2 -74s -1 -10) -s 2 log s + s -3 (1 + s) 2 (36s -1 -44 -4s -12s 2 + 2s 3 ) log(1 + s) + s(4 + s)s -3 (-144s -1 + 80 + 4s + 10s 2 -s 3 ) log √ 4 + s + √ s √ 4 + s - √ s p µ p ν p ρ + O( ) (4.13)
where m2 = m 2 e γ /(4π) with γ the Euler constant. The divergent term ∝ 1/ will be absorbed within the renormalization factors. On the other hand, this quantity is not ultraviolet divergent in d = 3 where it reads

Γ µνρ (p, 0, -p) = 1 - N g 2 0 128πm √ s π 2 (2 -3s 2 ) + 2 3s 3/2 (3s 3 + 23s 2 + 56s -15) -s -1 (4 + s)(16 -18s + s 2 ) arctan( √ s/2) + 2s -2 (s -5)(s + 1) 2 arctan( √ s) (p µ δ νρ + p ρ δ µν ) -2 + N g 2 0 128πm √ s π(-s 3 + 3s 2 -1) + 2 3s 3/2 (15s 3 -51s 2 + 53s -15) + 2(-s 3 + 6s 2 + 2s -16 -32s -1 ) arctan( √ s/2) + 2(2s 3 -9s 2 -5)(1 -s -2 ) arctan( √ s) p ν δ µρ - N g 2 0 128πm 3 √ s π(s -1 + s 2 ) + 2 3s 5/2 (-21s 3 + 5s 2 -139s + 75) + 2s -2 (s 4 -5s 3 -16s 2 -40s + 96) arctan( √ s/2) -2s -3 (s + 1)(2s 4 -7s 3 -9s 2 -15s + 25) arctan( √ s) p µ p ν p ρ (4.14)
The infrared limit when all momenta are smaller that the mass scale shows that

Γ µνρ ({κp i }) ∼ κ d-4 (4.15)
where all momentum have been multiplied by a common coefficient κ. As a consequence of this behaviour, G AAA diverges as log κ in d = 4 and diverges as 1/κ in d = 3 when κ → 0. This behaviour is entirely due to the ghost loop diagram D 1 because it is the only diagram with massless internal propagators allowing a singular behaviour in the IR.

Checks

In order to avoid errors in the computation of the diagrams we run different tests to check the correctness of the big expressions for the vertices. First of all, we have compared our expression in the limit of vanishing gluon mass with the one already obtained in the literature for standard Yang-Mills theory [START_REF] Davydychev | Three gluon vertex in arbitrary gauge and dimension[END_REF]. This check allows to compare independently all the scalar functions. The comparison was done numerically since the relation (3.2) introduces factors of the form 1/m that compensate in a very subtle way.

Even though the diagrams are regular in the limit of vanishing gluon mass, the analytic comparison is cumbersome.

There are other tests related to Slavnov-Taylor identities that our results have to fulfil. The first identity is the relation between both vertices and the propagators,

[Γ ⊥ (k)P ⊥ µρ (k) + Γ (k)P µρ (k)]Γ µν (p, k, r) -[Γ ⊥ (r)P ⊥ µν (r) + Γ (r)P µν (r)]Γ µρ (p, r, k) = p µ J -1 (p)Γ ρνµ (k, r, p) (4.16)
already known in the literature (see Appendix B.5.2). It is interesting that this is the standard identity and is not modified by the presence of the gluon mass. This relation was verified analytically with our one loop results. Another identity useful to check the results of the ghost-gluon vertex is deduced in Appendix B.5.3, and reads as

Γ µ (p, k, r) + Γ µ (k, p, r) - r µ r 2 p ν p 2 Γ ν (k, r, p) + k ν k 2 Γ ν (p, r, k) = 0 (4.17) where Γ µ (p, k, r) = k ν Γ νµ (p, k, r)r 2 J -1 (r). (4.18)
To our knowledge, this identity was not known before our study.

It is worth noting that the color structure for the three-point vertices was supposed to be the same as in the the one-loop one. The identity (4.17) was checked numerically for 50 arbitrary momentum configurations. This new identity gives a new non-trivial constraint for the ghost-gluon vertex.

Results

It is important to mention that some of the scalar functions cannot be obtained from lattice simulations. Lattice simulations compute correlations functions, and in consequence, the three-point vertex functions have their legs contracted with the full propagators, see relation (A.5). Since the gluon propagator is transverse, the longitudinal part of the vertex is not accessible. This is the case for the function c(r 2 , p 2 , k 2 ) defined in (4.5). However, we are going to compute the whole diagram because we can obtain information on this function through Slavnov-Taylor identities. The quantity simulated in the lattice for the ghost-gluon vertex by [START_REF] Cucchieri | Three-point vertices in Landaugauge Yang-Mills theory[END_REF] is a scalar function called G ccA , defined as

G ccA (p, r, k) = p ν P ⊥ µν (k)p ρ Γ ρµ (r, p, k) p ν P ⊥ µν (k)p µ (4.19)
where G ccA is the bare ghost-gluon vertex contracted with the corresponding propagator of the external legs normalized by the same tree level expression. This function depends on a unique linear combination of the scalar functions a, b, c, d, e, given by a(r 2 , p 2 , k 2 ) + p • k b(r 2 , p 2 , k 2 ) + d(r 2 , p 2 , k 2 ) + p 2 e(r 2 , p 2 , k 2 ). (4.20)

For the three-gluon vertex, lattice simulations on [START_REF] Cucchieri | Three-point vertices in Landaugauge Yang-Mills theory[END_REF] compute the scalar expression G AAA (p, r, k) defined as

G AAA (p, r, k) = [(k -r) γ δ αβ + cyclic perm]P ⊥ αµ (p)P ⊥ βν (r)P ⊥ γρ (k)Γ µνρ (p, r, k) [(k -r) γ δ αβ + cyclic perm]P ⊥ αµ (p)P ⊥ βν (r)P ⊥ γρ (k)[(k -r) ρ δ µν + cyclic perm] (4.21)
where the external legs of the vertex function are contracted with the propagators and the tree-level of the vertex, normalized to the same expression at the bare level.

Once we have computed the bare expressions (4.19) and (4.21) which will be denoted respectively G ccA B and G AAA B we have to renormalize the fields and coupling constants with the renormzalization factors defined in (1.20). In this way the lattice quantities are related to renormalized ones as follows

G ccA B (p, k, r, g B , m B ) = G ccA (p, k, r, g R , m R ) G AAA B (p, k, r, g B , m B ) = Z c Z A G AAA (p, k, r, g R , m R ). (4.22)
The lattice value of G AAA will differ with ours by a multiplicative factor due to the fact that our is renormalized and the lattice one is not (but their renormalization factors are finite because of the lattice regularization). The quantity G AAA depends on the lattice parameters so the multiplicative factor change for each set of lattice data. When we Figure 4.10 shows the results for G AAA for the same three momentum configurations. We want to stress that the comparison was done in the momentum configuration where lattice data is available although our expressions are valid for arbitrary kinematic configuration. The statistical errors for G AAA are larger than for the ghost-gluon correlation function and even much larger than for the propagators. Our expression for G AAA becomes negative and diverges like log(p) for very low momenta, this is not seen by lattice data because the simulations does not reach sufficiently small momenta.

It can be observed that the preferred scheme to reproduce the propagators in d = 4, which was the VM, is also the scheme that better reproduce these vertices.

d=3

The same analysis was done in three dimensions. In figure 3.22 the error contours show the region of parameters for which the gluon propagator and the ghost dressing function match the lattice data with less than 10% of error. In particular, the IS scheme with α = 0 reproduces with more accuracy the increase of the gluon propagator at low momentum. However, this effect is not reproduced by the VM scheme. The parameters chosen for the IS scheme are g 0 = 2.4 √ GeV and m 0 = 0.55 GeV which belong to the intersection of the contour errors for the gluon propagator and the ghost dressing function. As mentioned above this scheme is preferred to reproduce three-dimension lattice data. For the VM scheme the error contours do not intersect but we have considered a set of parameters in the middle, g 0 = 4.0 √ GeV and m 0 = 1.0 GeV. The results for the ghost-gluon factor, G ccA , is shown in figure 4.11, for both schemes. While in figure 4.12, the comparison for three-gluon factor, G AAA , is shown.

In all cases our results reproduce the qualitative behaviour of the lattice data. We find that G AAA becomes negative at small momenta and diverges like 1/p for vanishing momenta as found in [START_REF] Cucchieri | Three-point vertices in Landaugauge Yang-Mills theory[END_REF] for the lattice data.

In d = 3 the scheme that better reproduces the vertices is the IS scheme which is consistent with the study of the propagators. It is important to remark that in the following studies we will restrict ourself to the IS scheme. The reason is that it has been shown that this scheme considerably improves the propagators at d = 3 while the difference between both scheme in d = 4 is not very important. Moreover, the IS scheme do not have a Landau pole so we can choose the renormalization scale equals the momentum which correspond to α = 0. Chapter 5

Unquenched two-point correlation functions

In the previous chapters we have shown that a modest one loop calculation using the quenched approximation of the Lagrangian introduced in (2.11) remarkably reproduces the lattice results obtained for the quenched ghost dressing function, gluon propagator, ghost-gluon vertex and three gluon vertex. The comparison of all these quantities has been done by fixing just two parameters, the initial condition in the renormalization group flow of the coupling constant and the gluon mass.

In this chapter we extend the study to the unquenched case, i.e. taking into account the dynamics of the quarks. We will start by including the contribution of the quarks to the gluon propagator. The ghost propagator do not have contributions with internal quarks at one loop. In Sect. 5.1 we present the analytical results in arbitrary dimension in terms of Passarino-Veltman integrals [START_REF] Passarino | One Loop Corrections for e+ e-Annihilation Into mu+ mu-in the Weinberg Model[END_REF] and the simplified expressions in d = 4. Sect. 5.2 is devoted to the quark propagator. In the following section we implement the renormalization group approach in the IS scheme and we present the comparison with lattice data for the ghost and gluon sector. The calculation is done for arbitrary SU (N ) gauge group and for arbitrary number of fermions. For comparison with the lattice data, the fit is performed for QCD (N = 3 and d = 4) with two light degenerate quarks (N f = 2) and two light degenerate quarks in addition to two heavier quarks (N f = 2 + 1 + 1). In Sect. 5.5 the quark propagator, parametrized in the usual way in terms of two scalar functions is compared with the lattice data in the case with N f = 2 + 1. In the last section we present an estimation of the two-loop contributions in order to estimate error bars. All the results presented in this chapter were published in [START_REF] Peláez | Two-point correlation functions of QCD in the Landau gauge[END_REF].

Unquenched gluon propagator

The one-loop contribution to the gluon propagator is represented by the diagrams in Fig. 5.1.

A a µ A b ν = -(-g 2 ) d d q (2π) d Tr γ ν t b i / q + M q 2 + M 2 γ µ t a i( / q + / p) + M (q + p) 2 + M 2 = g 2 Tr t b t a d d q (2π) d -Tr (γ ν γ ρ γ µ γ σ ) q ρ (q + p) σ (q 2 + M 2 ) ((q + p) 2 + M 2 ) +Tr (γ ν γ µ ) M 2 1 (q 2 + M 2 ) ((q + p) 2 + M 2 ) = g 2 T f δ ab -4 (δ νσ δ ρµ -δ νµ δ ρσ + δ νρ δ µσ ) d d q (2π) d q ρ (q + p) σ (q 2 + M 2 ) ((q + p) 2 + M 2 ) +4M 2 δ µν d d q (2π) d 1 (q 2 + M 2 ) ((q + p) 2 + M 2 )
where T f is defined by Tr(t a t b ) = T f δ ab ( e.g. in the fundamental representation T f = 1/2). The mass of the internal quark is denoted by M . A diagram of this kind should be summed over all quark flavour, changing M for the corresponding mass.

This result can be expressed in terms of Passarino-Veltman integrals [PV79] A(M ) and B 0 (p, M, M ) defined as Following the procedure explained in the Appendix C, the above diagram can be expressed in terms of Passarino-Veltman integrals [START_REF] Passarino | One Loop Corrections for e+ e-Annihilation Into mu+ mu-in the Weinberg Model[END_REF] as

Γ 1loop,4 A a µ A b ν (p) = 2g 2 T f (d -1) [4M 2 -(d -2)p 2 ]B 0 (p, M, M ) + 2(d -2)A(M ) P ⊥ µν (p). (5.1)
where the integrals A(M ) and B 0 (p, M, M ) are defined as (see Appendix C)

A(m) = d d q (2π) d 1 q 2 + m 2 B 0 (p, m 1 , m 2 ) = d d q (2π) d 1 (q 2 + m 2 1 ) (q + p) 2 + m 2 2 .
(5.2)

The presence of quarks modifies the perpendicular part of the gluon self energy and therefore the propagator, however, the parallel part is not affected. The above expression was successfully checked with [START_REF] Davydychev | Three gluon vertex in arbitrary gauge and dimension[END_REF].

In d = 4 -the perpendicular part defined in 3.9 reads

Γ ⊥ 1loop,4 (p) = g 2 T f p 2 6π 2 - 1 + log M e γ/2 √ 4π - 5 6 + 2t + (1 -2t) √ 4t + 1 coth -1 √ 4t + 1 + O( )
where t = M 2 p 2 . It is important to note that for the ghost self-energy the only diagram that contributes at one loop, even in the unquenched case, is the one presented in figure 3.1. Therefore, the one-loop calculation of the ghost propagator is not directly modified with the dynamic quarks. Note, however, that the ghost dressing function is indirectly influenced by them through the renormalization group flow.

Quark propagator

The quark propagator is another fundamental quantity in QCD. The two-point quark vertex function, i.e. the inverse of quark propagator, can be parametrized in terms of two scalar components, A(p) and B(p), as

Γ (2) ψ ψ(p) = A(p)(-i / p) + M B(p).
We would like to compute the one-loop contribution to the quark propagator using the Lagrangian (2.11). It is given by the diagram of figure 5.3. Following the Feynman rules we can compute the diagram as µ, α ν, β q p + q p p (2)

ψ ψ,1loop = -g 2 d d q (2π) d γ ν t β i( / p + / q) + M (q + p) 2 + M 2 γ µ t α P ⊥ µν (q)δ αβ q 2 + m 2 = -g 2 C f i (γ ν γ ρ γ µ ) d d q (2π) d (p + q) ρ P ⊥ µν (q) ((q + p) 2 + M 2 ) (q 2 + m 2 ) +M (γ ν γ µ ) d d q (2π) d P ⊥ µν (q) ((q + p) 2 + M 2 ) (q 2 + m 2 )
where C f is defined by t a t a = C f Id. In the fundamental representation,

C f = N 2 -1
2N . The M here denotes the mass of the quark line. In this case, no sum has to be performed. The one loop contribution to the quark propagator only depends on the mass of the external quark flavour and the other quarks are not involved. This expression for the quark propagator can be inserted in our Mathematica algorithm in order to compute the integral by doing the Feynman's trick. Once the integrals are done we can identify the tensor structure and contract the external momentum with the Dirac matrices of each term.

For the first integral,

d d q (2π) d (p + q) ρ P ⊥ µν (q) ((q + p) 2 + M 2 ) (q 2 + m 2 )
we obtain terms proportional to p µ p ν p ρ , p µ δ νρ , p ν δ µρ and p ρ δ µν , which, contracted with γ ν γ ρ γ µ become p 2 / p, d / p, d / p and (2 -d) / p respectively. All of them contribute to the scalar function A(p). For the second integral,

d d q (2π) d P ⊥ µν (q) ((q + p) 2 + M 2 ) (q 2 + m 2 )
the tensor structure is simpler and there are only two terms. One, proportional to δ µν and the other to p µ p ν . The first term contracted with γ ν γ µ gives a d and the other one a p 2 . Since no Dirac gamma remains, this integral contributes to B(p).

The expression of the diagram can be written in terms of Passarino-Veltman integrals as

Γ (2) ψ ψ,1loop (p 2 ) = -g 2 C f (d -1)M B 0 (p, m, M ) + i / p g 2 C f 2m 2 p 2 (2 -d)m 4 + (d -3)m 2 (M 2 -p 2 ) + (M 2 + p 2 ) 2 B 0 (p, m, M ) +(M 2 + p 2 ) 2 B 0 (p, 0, M ) + A(m) (2 -d)m 2 -M 2 -p 2 + (d -2)m 2 A(M ) .
The previous expression coincides with that of [START_REF] Davydychev | Quark gluon vertex in arbitrary gauge and dimension[END_REF] when the gluon mass m is set to zero. It is worth to mention another test: in each case the divergent part when d = 4 is compared with the results of [START_REF] Gracey | Three loop MS-bar renormalization of the Curci-Ferrari model and the dimension two BRST invariant composite operator in QCD[END_REF].

In d = 4 -, the expression can be written analytically as

Γ (2) ψ ψ,1loop (p 2 ) = i / pg 2 C f 64π 2 m 2 p 4 k 2 2m 4 + m 2 (p 2 -M 2 ) -(M 2 + p 2 ) 2 Q -2m 2 p 2 (-2m 2 + M 2 + p 2 ) -2[2m 6 + 3m 4 (p 2 -M 2 ) + (M 2 + p 2 ) 3 ] log M m -2(M 2 + p 2 ) 3 log M 2 + p 2 M 2 + 3g 2 C f M 8π 2 - 1 + log me γ/2 √ 4π - 2 3 - k 2 4p 2 Q + 1 2p 2 (m 2 -M 2 + p 2 ) log M m + O( )
where

k 2 = m 4 + 2m 2 (p 2 -M 2 ) + (M 2 + p 2 ) 2 and Q = log (k 2 -p 2 ) 2 -(M 2 -m 2 ) 2 (k 2 + p 2 ) 2 -(M 2 -m 2 ) 2 .
(5.

3)

It can be noted that the part proportional to / p is finite as it is well known in Landau gauge at one loop. Therefore, the renormalization factor of the fermion field is also finite. Moreover, even the finite part proportional to / p vanishes in the limit of vanishing gluon mass. In order to complete the test of this expression we computed the divergent part in an arbitrary linear gauge and checked it with [START_REF] Gracey | Three loop MS-bar renormalization of the Curci-Ferrari model and the dimension two BRST invariant composite operator in QCD[END_REF].

Renormalization and Renormalization group

As discussed in the previous chapters, in order to absorb the divergences in the vertex functions we redefine the coupling constant, masses and fields. The renormalization factors were already introduced in (1.20). Until now we have used the factors related to the ghost -gluon sector but here, we include the renormalization factor for the quark field and quark masses, defined as

ψ a B = Z ψ ψ a and M B = Z M M
where the subindex B denotes the bare quantities. The renormalized Γ (2) is related with the bare one as

Γ (2),R ψ ψ (p) = Z ψ Γ (2),B ψ ψ (p) = Z ψ A(p)(-i / p) + Z M M B(p) .
In this chapter we focus on the infrared safe IS scheme since although the VM scheme seems give better results in d = 4 the difference between the schemes is small while in d = 3 the IS scheme really improves the results. The unquenched IS scheme is defined by Eqs (3.16) together with

Γ (2),R ψ ψ (µ) = -i / µ + M. (5.4) Therefore, Z ψ = A -1 (µ) and Z M = A(µ)B -1 (µ). (5.5)
In order to avoid large logarithms the running scale µ has to be chosen as µ p for µ m for a correlation function with typical momentum p. The β-function for the quark mass and quark anomalous dimension are given in (1.21).

In this chapter we will focus on the physical dimension d = 4 and all the results will take into account the renormalization group treatment. In four dimension the contribution of the quarks to the gluon anomalous dimension is

γ quarks A = N f i=1 g 2 T f π 2    t i -6 6t i - 2 t 3 i (t i + 4) log √ t i + 4 - √ t i √ t i + √ t i + 4   
where t i = µ 2 /M 2 i . Even though the ghost anomalous dimension remains untouched, this change in the gluon anomalous dimension also affects the β-functions in the IS scheme,

β m 2 = m 2 (γ A + γ c ) and β g = g(γ A /2 + γ c ).
This time we have to consider the fermion anomalous dimension, given by

γ ψ = g 2 C f 16π 2 m 2 µ 4 m 2 µ 2 4m 2 -2M 2 + µ 2 + -2 2m 6 -3m 4 M 2 + M 6 -3µ 2 M 4 + m 4 + µ 6 log M m + -2M 6 -3M 4 µ 2 + µ 6 log M 2 + µ 2 M 2 + 1 2k 2 2(m 2 -M 2 ) 3 (2m 2 + M 2 ) + µ 2 (m 2 -M 2 )(7m 4 + 6m 2 M 2 + 5M 4 ) +µ 4 3m 4 -2m 2 M 2 -3M 4 + µ 6 m 2 + M 2 + µ 8 Q and β M = M γ ψ + 3g 2 C f M 8π 2 -µ 2 + (m 2 -M 2 ) log M m - 1 2k 2 m 4 + m 2 µ 2 -2M 2 + M 2 M 2 + µ 2 Q
where k and Q are the expressions appearing in Eq. (5.3) with p replaced by the renormalization-group scale µ.

Results for the gluon and ghost sectors

We would like to compare the unquenched gluon propagator and the ghost dressing function with the lattice data available. Lattice data shows that the qualitative behaviour remains the same with the inclusion of quarks. Our analytical expression has been computed for arbitrary number of colours N , number of flavours N f and dimension. In this section we compare the results in four dimensions for SU (3) with two degenerate light quarks N f = 2 (both quarks have the same mass), and, with two degenerate light quarks and two heavier quarks N f = 2 + 1 + 1 (there are three different quark masses).

Our results depend on the choice of the initial conditions for the renormalization group flow. We use the coupling constant and masses at the scale µ 0 = 1 GeV as fitting parameters. The best fit between our results and lattice data take place when the value of the fitting parameters minimize simultaneously the relative errors defined as

χ 2 AA,Unq = 1 N i Γ ⊥ lt. (p i ) 2 1 Γ ⊥ lt. (p i ) - 1 Γ ⊥ th. (p i ) 2 χ 2 cc,Unq = 1 N i J -2 lt. (p i ) (J lt. (p i ) -J th. (p i )) 2 .
(5.6)

We can observe that there is a small difference between the errors used in the quenched analysis (3.23) and the last ones (5.6). In the quenched analysis we choose to balance the relative error with the absolute error. However, we have seen that the information given by the relative error is practically the same as the defined in Eq.(3.23), considering that there is a expanded region of possible values, so we decided to use the relative error instead.

We first study the case with two light fermions with the same mass. Therefore, we have three fitting parameters to choose: the coupling constant, the gluon mass and the quark mass. First, we fix the quark mass with the approximate lattice value, M u,d 0 = 0.13 GeV at 1 GeV. Second, we plot the error contours in terms of the coupling constant and gluon mass, see figure 5.4. Third, we choose the values of the parameters belonging to the intersection of the contours. In this case g 0 = 4.5 and m = 0.42 GeV. The corresponding gluon and ghost propagators are shown in figure 5.5 where we can see an excellent match with the lattice data of [START_REF] Sternbeck | Determination of LambdaMS from the gluon and ghost propagators in Landau gauge[END_REF]. Finally, we fix the value of the coupling constant of the third step and make the error contour lines as function of the gluon and quark masses. The last plot displayed in figure 5.6 shows us that the quality of the fits are almost independent of the quark mass as long as M 1 GeV. In the case N f = 2 + 1 + 1, we have a priori five parameter to fit. They correspond to the initial conditions of the renormalization group flow at 1 GeV of the coupling constant, gluon mass and three quark masses. However, we impose at the scale µ 0 = 1 GeV the Figure 5.4 -Contour levels for the quantities χ AA,Unq (forward, oval form) and χ cc,Unq (backward) for d = 4, for N f = 2 and M u,d 0 = 0.13 GeV. The contour lines correspond to 4%, 7% and 10 % for the gluon and 2%, 4% and 6 % for the ghost. same ratio for the quark masses as those used in lattice simulations. That means that the initial condition of the middle quark mass is twice the lightest mass and the heavier quark mass is twenty times the value of the light mass at 1 GeV. These ratios give at 2 GeV the ratios used by lattice simulations in [ABB + 12] which are the data we are going to compare with.

The error contours in four dimensions with quark masses respecting the ratios mentioned above and M u,d 0 = 0.13 GeV fixed at 1 GeV are shown in figure 5.7. The fitting parameter that we choose are g 0 = 5.3 and m 0 = 0.48 GeV which give the fits depicted in figure 5.8. Again in this case, lattice result are reproduced with great accuracy. The comparison with lattice data shows that the new diagram with an internal loop of quarks is not strong enough to modify the qualitative behaviour of the propagators.

As in the previous case, we are not paying attention of the quark mass since the results seem to be independent of small changes of this value (as long as the lighter quarks have a mass much smaller than 1 GeV). This independence with respect to the quark mass value can be seen in figure 5.9.

It is important to note that in this case we have three parameters to fit but, as discussed for the quenched case, the coupling constant should match with the known value after the appropriate running. For instance, in [BBB + 12] the coupling constant using Taylor scheme is computed at the Z 0 -boson mass scale where g(m Z 0 ∼ 96) = 1.2274. If we perform the renormalization group flow of the coupling constant obtained by us at 1 GeV, g 0 = 5.3, it should match to their value at µ = m Z 0 . At the momentum scale m Z 0 ∼ 96 GeV five quarks must be considered, since the mass of the up, down, charmed, strange and bottom quarks is smaller than m Z 0 . Therefore, we use the renormalization group flow for the for N f = 2 and g 0 = 4.5. The contour lines correspond to 6%, 10% and 12 % (left), 3%, 5% and 8 % (right). The quality of the fits is almost insensitive to the value of the quark mass.

N f = 2 + 1 + 1 case up to the scale of the bottom mass and above that scale we include the bottom quark with a frozen mass in the renormalization group equations. The coupling constant obtained at m Z 0 is g = 1.4. This gives a 17% error on the coupling constant that is the typical error of the one loop approximation. Even though the coupling constant is not a free parameter it is convenient to fit it at 1 GeV in order to study the infrared regime better. We always choose the initial condition of the coupling constant within this error bar.

It is interesting to analyse the effect of increasing the number of flavours. Figure 5.10 shows the gluon and ghost dressing function for both N f = 2 and N f = 2 + 1 + 1. The comparison is done by normalizing the curves such that they coincide at 4 GeV so their infrared behaviour can be easily compared. It can be seen that the gluon dressing function decreases with the addition of more quarks as was observed in lattice simulations [BHL + 04, ABB + 12]. It is important to remark that both curves have been done using the corresponding fitting parameters in each case. We also observed an enhancement of the ghost dressing function.

Results for the quark sector

In this section we compare our one loop calculations for the quark propagator shown in section 5.2 with the lattice data available for this quantity. The data published in [BHL + 04, BHL + 05] for the gluon and the light quark propagator respectively correspond to the SU (3) group in four dimensions and N f = 2 + 1, i.e. two degenerate light quarks and one heavy quark. Therefore we restrict our calculations to that case. It is important to note that in the previous analysis we restricted ourself to the ghost-gluon sector because no quark data was available.

In this section we have three quantities available from the lattice data to compare with. They correspond to the gluon propagator and the two scalar functions related to the light To compare with the lattice data we include the effects of the renormalization group and therefore we can relate

Γ (2) ψ ψ(p, µ 0 ) = 1 z ψ (µ) Γ (2) ψ ψ(p, µ) = Z ψ z ψ (µ) A(p)(-i / p) + Z M M (µ)B(p) = A -1 (µ) z ψ (µ) A(p)(-i / p) + A(µ)B -1 (µ)M (µ)B(p) µ=p = 1 z ψ (µ = p) -i / p + M (µ = p) (5.7)
where we have used the relations (5.5). We note that M u,d (p) corresponds to the M u,d (µ) (the light quark running mass) changing µ by p. Therefore, at 1 GeV, M u,d (1GeV) = M u,d 0 . We observe also that Z u,d (p) = z -1 ψ (µ = p), however, due to the different renormalizations when comparing with the lattice data we have to include a multiplicative factor.

As in the previous case, we first fix the parameters g 0 and m 0 to obtain a good fit for the gluon propagator with M u,d 0 = 0.13 GeV. For the heavy quark mass, the initial condition at 1 GeV is fixed to be five times the value of M u,d 0 which is equal to the ratio used by lattice simulations. Figure 5.11 shows the contour lines for the error quantities χ AA,Unq and χ CC,Unq . We mention that the analysis for the ghost dressing function in this section is done by using the lattice data for the case with N f = 2 + 1 + 1 as no lattice data is available for N f = 2 + 1. This is supported by the idea that this function is rather insensitive to the inclusion of a heavy quarks. Under this hypothesis, it is interesting to compare our findings with the ghost dressing function for N f = 2 + 1 + 1.

Figure 5.11 -Contour levels for the quantities χ AA Unq (forward) and χ cc Unq (backward, oval form) for d = 4 with N f = 2 + 1. The contour lines correspond to 5.5%, 6.5% and 7.5 % for the gluon and 2%, 3% and 4 % for the quark mass. The best fit for both dressing functions was obtained choosing at 1 GeV, g 0 = 4.8 and m 0 = 0.42 GeV and it is depicted in figure 5.12. As in the previous section both functions are almost independent of the value of M 0 , as can be seen in figure 5.13. However, the quark propagator itself strongly depends on the choice of M 0 . In figure 5.14 we show that those values for the coupling constant and gluon mass are not good enough for fitting the scalar function M u,d (p) using M u,d 0 = 0.13 GeV. In fact, figure 5.15 shows that in order to reproduce M (p) maintaining M u,d 0 = 0.13 GeV at 1 GeV the values for g 0 and m 0 needed are far away of the expected rank (being g 0 = 9.6 and m 0 = 1.4 GeV). Clearly, these values do not reproduce the gluon propagator, see figure 5.16.

It is important to note that the lattice results for the mass of the quark show the consequence of the chiral symmetry breaking. We can see that even if the starting value with g 0 = 4.8. The contour lines correspond to 6%, 7% and 8 % for the gluon and 2.5%, 3% and 4 % for the ghost.

of the quark mass at large momentum is small, the quark mass increases considerably in the infrared. In principle, this should happens if we set an arbitrary small value of the quark mass at high energy. However, our one loop calculation does not reproduce this feature. Instead, if the ultraviolet quark mass is set to zero, the mass remain zero at every momentum scale. However, when we considered a small non zero value for the quark mass the situation is similar enough to the lattice data.

In order to find the correct value of M u,d 0 to reproduce the lattice data we fixed g 0 = 4.8 (extracted from figure 5.11) and study the error dependence on the masses. The left figure of 5.17 shows the contour errors for χ M defined as

χ 2 M = 1 N i M -2 lt. (p i ) (M lt. (p i ) -M th. (p i )) 2 .
The right figure of 5.17 shows the superposition of the error contours χ M with the contours of χ AA,Unq . It can be observed that the fitting for the M (p) does not depend on the gluon mass for masses of the gluon bigger than 0.4 GeV. The parameters found to reproduce the gluon propagator in figure 5.12 belong to the region of smallest χ M . Therefore, we use the same coupling constant and gluon mass used to reproduced the gluon dressing propagator and we change the initial condition for the quark mass to M u,d 0 = 0.08 GeV. In order to complement our analysis we also compare the ghost dressing function. This comparison is just qualitative due to the fact that we used the ghost data for N f = 2 + 1 + 1 as no lattice data was available in the current case. The best fit is shown in figure 5.18 where g 0 = 4.8, m 0 = 0.42 GeV and M u,d 0 = 0.08 GeV. We also show in figure 5.19 the masses for the degenerate light quark and for the heavy quark.

We can conclude that a simple one loop calculation for the gluon and ghost dressing functions and the quark mass reproduces the data obtained in lattice simulations with great accuracy for gluon and ghost dressing functions and reasonable accuracy for the quark mass.

There remains to compare Z u,d (p). Let us recall that using the IS scheme Z u,d (p) matches with z ψ (µ), up to a multiplicative constant, replacing µ by p. The result for Z u,d (p) is shown in figure 5.20 where we can see that the results are not good at all. In fact, the behaviour is not even qualitatively reproduced. We attribute this mismatch to the fact that the one loop contribution to this function is really small. The reason for that is that in standard QCD the one loop contribution to this function vanishes. The important correction to the tree-level expression arrives from the two-loop calculations. Therefore, it is natural to expect that for the massive case also the two-loop contributions are not negligible for describing Z(p).

In order to test this explanation, in the next section, we estimate the contribution of the two-loop diagrams.

Estimate of the two loop contributions

The results of the previous section show that the function Z u,d (p) is not well reproduced by the one-loop calculations. Although, one-loop contributions are enough to describe with great accuracy the gluon and ghost dressing functions and with reasonably accuracy the light quark mass function. The Z u,d (p) function has no one-loop contribution when the gluon mass is set to zero. In that case, two loop contributions dominate the one loop corrections. In this section, we argue that two-loop contributions are important in our case in order to reproduce Z u,d (p), even if they do not influence in a significant manner the dressing functions and M u,d (p).

The completely calculation of the two-loops diagrams of momentum dependent quantities with massive propagators is extremely demanding so we are not going to compute them here. However, we make a first estimation of the two-loop contribution. For that purpose, we use the ultraviolet behaviour of the β functions and anomalous dimensions computed in [START_REF] Gracey | Three loop MS-bar renormalization of the Curci-Ferrari model and the dimension two BRST invariant composite operator in QCD[END_REF]. Since all the two-loop diagrams of the quark propagator have at least one gluon propagator, a two-loop calculation in Curci-Ferrari model would be suppressed in the infrared. We mimic this effect by modifying the anomalous dimension in a minimal way by multiplying them by where s is of order 1 GeV. Let us call the results hybrid expressions. In this way, the hybrid anomalous dimension matches with the one obtained at two loops by [START_REF] Gracey | Three loop MS-bar renormalization of the Curci-Ferrari model and the dimension two BRST invariant composite operator in QCD[END_REF] in the ultraviolet limit and goes to zero as µ 2 in the infrared, as expected.

µ 2 µ 2 + s 2
We repeat the same procedure as before in order to find the best fitting parameters. The best fits are obtained for g 0 = 3.4, m 0 = 0.4 GeV and M = 0.08 GeV at µ 0 = 1 GeV and lead to the curves presented in the figure 5.21.

We can estimate the error of our one-loop calculations to be of the order of the difference between our one loop results and the hybrid results. The shaded areas of figure 5.22 diagrammatically represent a region centred at the hybrid expression with a shaded region and limited by the one loop results and by the symmetric curve of one loop results with respect to the hybrid model. This figure shows that higher correction for the gluon dressing function are small and of order comparable with the difference between our one-loop calculations and the lattice data. On the other hand, the higher order contributions to the quark mass are small but they are not enough to reproduce the lattice data. However, we still see a quark mass which has small values in the ultraviolet but is abruptly enhanced in the infrared as is observed in lattice data, as a consequence of the chiral symmetry Chapter 6

One loop corrections to the quark-gluon vertex

The quark-gluon vertex is another fundamental ingredient of QCD which controls the interaction between quarks and gluons. It is believed that, at least in Landau gauge, it plays an important role in the confinement mechanism and in the dynamical chiral symmetry breaking [MR03, RW94, FA03, AP11, CR14]. It is also essential for the formation of physical bound states [MT99, BDRT02, BHK + 04, HKR05, CR09, Wil14]. For instance, the meson spectrum can be studied by using Bethe-Salpeter equation [START_REF] Hrayr | Quark-gluon vertex dressing and meson masses beyond ladder-rainbow truncation[END_REF], taking as an input this vertex and the propagators of the theory.

The strength of the interaction between quarks and gluons is measured through the running coupling constant α s = g 2 4π which can be determined from this vertex as explained in [START_REF] Bethke | Determination of the QCD coupling α s[END_REF]. As we have mentioned before, the coupling constant obtained through this vertex or any other of the vertices in the gluon sector give the same value in the ultraviolet at two-loop order. However, at low momenta the α s extracted from this vertex is larger than the one computed from the ghost-gluon vertex or from the three-gluon vertex as it is shown in [START_REF] Skullerud | Quark gluon vertex from lattice QCD[END_REF].

For these reasons, there is a huge interest in determining the behaviour of the quarkgluon vertex. For instance, it has been computed at one loop in the limit of vanishing gluon mass by [START_REF] Davydychev | One loop results for the quark gluon vertex in arbitrary dimension[END_REF] in arbitrary linear covariant gauge and dimension. Some partial computations in particular gauges and momentum configurations have been done at two and three loops in [START_REF] Chetyrkin | Two loop QCD vertices and three loop MOM beta functions[END_REF][START_REF] Chetyrkin | Three loop three linear vertices and four loop similar to MOM beta functions in massless QCD[END_REF] in the same zero gluon mass limit. It was also studied by Dyson-Schwinger methods in [START_REF] Aguilar | A new method for determining the quark-gluon vertex[END_REF][START_REF] Williams | The quark-gluon vertex in Landau gauge bound-state studies[END_REF] and by lattice simulations in [SK02, SBK02, SKB + 04, SBK + 05]. However, the understanding of this vertex is not complete since its kinematics has a complex tensor structure with 12 scalar components which depend on three independent kinematic variables. The lattice simulations for the general kinematic has been done only for a particular combination of the scalar functions but are not as precise as in the zero momentum gluon case.

In this chapter present the computation of the quark-gluon vertex at one loop using the Lagrangian (1.14) which is the standard QCD Lagrangian with the addition of a gluon mass term. The preceding chapters showed that this little change in the Lagrangian reproduces with good accuracy the results for the quenched and unquenched propagators and the three point correlation functions in the gluon sector. All these computations are done in arbitrary dimension, number of colors and flavours and in general kinematics. The comparisons between our results and the lattice simulations to be presented are done for N = 3 in d = 4 and published in [START_REF] Peláez | Quark-gluon vertex in the Landau gauge[END_REF].

Tensor structure

We begin by presenting the tensor structure of the quark-gluon vertex which is represented in figure 6.1 where p, r and k denote the quark, anti-quark and gluon momentum respectively. The momentum conservation tells us that r = -k -p. The Lorentz structure of the vertex consists of 12 independent tensors. They can be grouped in the transverse terms (8 components) and the non-transverse ones as it was done in [START_REF] Skullerud | Quark gluon vertex from lattice QCD[END_REF]. The vertex function can be written as follows

Γ ψ ψA a µ (p, r, k) = t a Γ µ (p, r, k) = -igt a 4 i=1 λ i (p 2 , r 2 , k 2 )L i,µ + 8 i=1 τ i (p 2 , r 2 , k 2 )T i,µ
where

L 1,µ = γ µ L 2,µ = -( / p -/ r)(p -r) µ L 3,µ = -i(p -r) µ L 4,µ = -iσ µν (p -r) ν T 1,µ = i k µ (r.k) -r µ k 2 T 2,µ = k µ (r.k) -r µ k 2 ( / p -/ r) T 3,µ = / kk µ -k 2 γ µ T 4,µ = -i k 2 σ µν (p -r) ν -2k µ σ νλ r ν k λ T 5,µ = iσ µν k ν T 6,µ = / k(p -r) µ -k.(p -r)γ µ T 7,µ = - i 2 k.(p -r) ( / p -/ r)γ µ -(p -r) µ -i(p -r) µ σ νλ r ν k λ
T 8,µ = -γ µ σ νλ r ν k λ -/ rk µ + / kr µ (6.1)

Symmetries

The scalar functions respect some symmetries due, for example, to the effect of the charge conjugation applied to the quark-gluon vertex, which implies that

CΓ µ (p, r, k)C -1 = -Γ T µ (r, p, k).
Remembering that the charge conjugation acts over the Dirac matrices as

Cγ µ C -1 = -γ T µ ,
we obtain that all the scalar functions are symmetric under the interchange of the quark and anti-quark momenta except for λ 4 , τ 4 and τ 6 which are odd. These properties were verified by our calculations.

Let us also study the constraint of the chiral symmetry over the scalar functions. With this purpose we consider the part of the Lagrangian involving quarks but taking the quark mass to zero L = -ψγ µ D µ ψ.

We can define the left-handed and the right-handed fermions as

ψ L = 1 -γ 5 2 ψ and ψ R = 1 + γ 5 2 ψ
and therefore ψL = ψ 1 + γ 5 2 and ψR = ψ 1 -γ 5 2 .

This gives the relations ψ = ψ L + ψ R and ψ = ψL + ψR . The Lagrangian can be expressed as

L = -ψL γ µ D µ ψ L -ψR γ µ D µ ψ R .
This Lagrangian is invariant under the chiral symmetry defined as

ψ L → U L ψ L = e iθ a L t a ψ L = ψ L + iθ a L t a ψ L + O(θ 2 L ) and ψ R → U R ψ R = e iθ a R t a ψ R = ψ R + iθ a R t a ψ R + O(θ 2 R )
which can be performed independently for each quark flavour. This symmetry leads to the following Ward identity Γ ψ ψA a µ (p, r, k)γ 5 + γ 5 Γ ψ ψA a µ (p, r, k) = 0 which implies that, in the case with vanishing quark mass, the quark-gluon vertex anticommutes with γ 5 .

In particular this implies that λ 3 , λ 4 , τ 1 , τ 4 , τ 5 and τ 7 vanish in the chiral limit. We checked that in our one-loop results, where these functions are proportional to M in the limit of small M . We will see that the running of M is really important for these quantities and they have a behaviour that is similar to the one of the chiral symmetry breaking when the renormalization group is considered. The other functions, instead, are of order one in this limit.

One loop contribution

We consider the one loop contribution for this vertex given by our model. The diagrams involved are shown in figure 6.2. The first diagram has a similar structure as the one loop contribution to QED except for a multiplicative factor. It is shown in detail in figure 6.3 and its expression, denoted Γ a,1loop µ,1 (p, r, k), is the following

(p, r, k) CB = -(ig) 3 d d q (2π) d γ ρ (t c ) CD i( / q -/ k -/ p) + M (q -k -p) 2 + M 2 γ µ (t a ) DE × i( / q -/ p) + M (q -p) 2 + M 2 γ ν (t b ) EB P ⊥ νρ (q) q 2 + m 2 δ bc = ig 3 t b t a t b CB d d q (2π) d γ ρ i( / q -/ k -/ p) + M (q -k -p) 2 + M 2 γ µ i( / q -/ p) + M (q -p) 2 + M 2 γ ν P ⊥ νρ (q) q 2 + m 2 = - i 2N g 3 (t a ) CB d d q (2π) d γ ρ i( / q -/ k -/ p) + M (q -k -p) 2 + M 2 γ µ i( / q -/ p) + M (q -p) 2 + M 2 γ ν P ⊥ νρ (q) q 2 + m 2
where the Latin capital letters going from 1 . . . N denotes color of the quark, the Latin small letters from 1 to N 2 -1 represent the color indices of the gluons and, as usual, the greek letters represent the Lorentz indices. In the last line we have used that t b t a t b = (C f -C a /2) t a = -1/(2N )t a , see [START_REF] Peskin | An Introduction to quantum field theory[END_REF]. Where we have used that, for SU (N ),

C a = N and C f = N 2 -1 2N .
We can split the integral according to the structure in terms of the Dirac gamma matrices. This leads to

Γ a,1loop µ,1 (p, r, k) = -i 2N g 3 t a (I 1 + I 2 + I 3 + I 4 )
where

I 1 = γ ρ γ γ µ γ η γ ν - d d q (2π) d (q -k -p) (q -k -p) 2 + M 2 (q -p) η (q -p) 2 + M 2 P ⊥ νρ (q) q 2 + m 2 I 2 = γ ρ γ γ µ γ ν iM d d q (2π) d (q -k -p) (q -k -p) 2 + M 2 1 (q -p) 2 + M 2 P ⊥ νρ (q) q 2 + m 2 I 3 = γ ρ γ µ γ η γ ν iM d d q (2π) d 1 (q -k -p) 2 + M 2 (q -p) η (q -p) 2 + M 2 P ⊥ νρ (q) q 2 + m 2 I 4 = γ ρ γ µ γ ν M 2 d d q (2π) d 1 (q -k -p) 2 + M 2 1 (q -p) 2 + M 2 P ⊥ νρ (q) q 2 + m 2
The expressions between squared brackets are computed with the Mathematica algorithm which is in charge of doing the Feynman trick and integrate over the internal momentum. We can then identify the structure in momentum variables and perform the contraction between them and the corresponding gamma matrices. We can finally identify the contribution to each scalar function defined in Eq.(6.1). The calculation reproduces the results of [DOS00] in the limit of vanishing gluon mass.

It can be seen that I 1 and I 4 contribute to λ 1 , λ 2 , τ 2 , τ 3 , τ 6 and τ 8 scalar functions while I 2 and I 3 contribute to the remaining ones. An interesting feature is that no contribution appears for λ 4 at one loop.

Moreover, it is known that the λ i scalar functions of the QED-like vertex are related to the quark propagator through an identity identical to the abelian Ward identity of QED up to a multiplicative factor. The identity, valid only for the QED-like diagram is

k µ Γ a,1loop µ,1 (p, r, k) = (C f -1/2C a ) C -1 f Γ (2) ψ ψ -1 (-r) -Γ (2) ψ ψ -1 (p) 1loop .
This identity was successfully checked by our result.

The calculation was also done using the Mathematica's package, FeynCalc [START_REF] Mertig | FEYN CALC: Computer algebraic calculation of Feynman amplitudes[END_REF], which allows us to rewrite the expression in terms of Passarino-Veltman integrals. To allow FeynCalc to reduce our equation to a manageable size it was necessary to simplify the input by using, for instance, the substitution (3.2). The result of the scalar functions using this package is presented in the complement material of [START_REF] Peláez | Quark-gluon vertex in the Landau gauge[END_REF].

The second diagram is shown in detail in figure 6.4. Its contribution to the quark-gluon vertex is given by This diagram does not satisfy a simple abelian Ward identity but the full vertex satisfy a Slavnov-Taylor identity (see for instance the presentation in [START_REF] Davydychev | One loop results for the quark gluon vertex in arbitrary dimension[END_REF]). The corresponding check is beyond the scope of the present calculation because we would need to compute other vertices appearing in that Slavnov-Taylor identity.

Γ a,1loop µ,2 (p, r, k) CB = -(-g) 2 d d q (2π) d γ η (t e ) CD i( / q -/ p) + M (q -p) 2 + M 2 γ σ (t d ) DB δ db P ⊥ σν (q) q 2 + m 2 δ ce × P ⊥ ηρ (k + q) (k + q) 2 + m 2 igf abc [(2q + k) µ δ νρ -(2k + q) ν δ µρ + (k -q) ρ δ µν ] = ig 3 t c t b CB f abc d d q (2π) d γ η i( / q -/ p) + M (q -p) 2 + M 2 γ σ P ⊥ σν (q) q 2 + m 2 P ⊥ ηρ (k + q) (k + q) 2 + m 2 × [(2q + k) µ δ νρ -2k ν δ µρ + 2k ρ δ µν ] = N 2 g 3 (t a ) CB d d q (2π) d γ η i( / q -/ p) + M (q -p) 2 + M 2 γ σ P ⊥ σν (q) q 2 + m 2 P ⊥ ηρ (k + q) (k + q) 2 + m 2 × [(2q + k) µ δ νρ -2k ν δ µρ + 2k ρ δ µν ] k µ, a B C p k + p q q -p -(k + q)
Again we can separate the integral in two parts according to the number of Dirac matrices.

Γ a,1loop µ,2 (p, r, k) CB = N 2 g 3 (t a ) CB (I 5 + I 6 )
where

I 5 = γ η γ γ σ i d d q (2π) d (q -p) (q -p) 2 + M 2 P ⊥ σν (q) q 2 + m 2 P ⊥ ηρ (k + q) (k + q) 2 + m 2 [(2q + k) µ δ νρ -2k ν δ µρ + 2k ρ δ µν ]
and

I 6 = γ η γ σ M d d q (2π) d 1 (q -p) 2 + M 2 P ⊥ σν (q) q 2 + m 2 P ⊥ ηρ (k + q) (k + q) 2 + m 2 [(2q + k) µ δ νρ -2k ν δ µρ + 2k ρ δ µν ] .
After performing the integral using the Mathematica algorithm we can observe that the integral I 5 contributes to λ 1 , λ 2 , τ 2 , τ 3 , τ 6 and τ 8 while I 6 contributes to λ 3 , λ 4 ,τ 1 ,τ 4 , τ 5 and τ 7 .

In order to have an idea of the order of magnitude of the scalar functions obtained in the calculation we present in figure 6.5 all the scalar functions for a fixed value of the gluon momentum k = 0.555 GeV and the anti-quark momentum r = 0.844 GeV as function of p (except λ 1 that requires to be renormalized and will be discussed below). The general expression of the scalar functions in the case of vanishing gluon mass was successfully compared numerically with the results obtained for the standard QCD by [START_REF] Davydychev | Quark gluon vertex in arbitrary gauge and dimension[END_REF].

Renormalization scheme

In order to compare our calculations with the lattice simulations we implement the renormalization scheme IS defined in (3.16) and (5.4). It is important to remark that no extra renormalization condition needs to be imposed to the quark-gluon vertex. Since the renormalization factors already defined are enough to renormalize this quantity. Therefore, we have that

Γ R ψ ψA a µ (p, r, k) = Z ψ Z A Γ B ψ ψA a µ (p, r, k) = t a Z ψ Z A Z g Γ B µ (p, r, k)
where the combination Z ψ √ Z A Z g absorbs the divergences of this vertex. The relation between the renormalized vertex in two different scales is given by

Γ R ψ ψA a µ (p, r, k, µ 0 ) = 1 z ψ (µ) z A (µ) Γ R ψ ψA a µ (p, r, k, µ)
but, as we want to represent the scalar quantities after factorizing (-it a g), the scalar functions that we are going to compare can be written, up to a multiplicative factor, as

λ R i (g(µ 0 ), m(µ 0 ), M (µ 0 )) = Z ψ √ Z A Z g z ψ (µ) z A (µ) g(µ) g(µ 0 ) λ B i (g(µ), m(µ), M (µ)) and τ R i (g(µ 0 ), m(µ 0 ), M (µ 0 )) = Z ψ √ Z A Z g z ψ (µ) z A (µ) g(µ) g(µ 0 ) τ B i (g(µ), m(µ), M (µ))
where we have made explicit that the effects of the renormalization group flow are considered.

Results

The one loop vertex was computed for generic number of colors, flavours, arbitrary dimension and general kinematics. Our one-loop calculations are renormalized using the IS-scheme at µ 0 = 1 GeV and we take the running scale µ = r 2 +p 2 +k 2

2

. The choice of this scale is motivated so it matches with µ ∼ p for vanishing gluon momentum and when p 2 = r 2 k 2 , µ ∼ p and in the opposite case the running scale is frozen out by the gluon momentum.

In this section we compare the scalar functions with the lattice data available in the literature. The lattice simulations for the quark-gluon vertex have been done for QCD in the quenched approximation in d = 4. It is important to remark that the gluon and ghost propagators were already studied with our methods for quenched QCD in section 3.3. In that study the coupling constant and the gluon mass were fixed to reproduce the lattice data of the propagators. We use here the same values of the parameters, g = 4.1 and m = 0.43 GeV at 1 GeV, for comparing the quark-gluon vertex. Therefore, we only have to fix the initial condition of the quark mass in the RG flow. In this section we show how the parameters used to fit the quenched propagators allow us to reproduce with good accuracy the lattice simulations for the quark-gluon vertex.

Simulations with small error have been done for the case corresponding to the vanishing gluon momentum by [SBK + 03]. We present the comparison of our results, also considering the quenched approximation, with the lattice data of the three scalar functions involved in this configuration. Another configuration treated in [SBK + 03] is the quark reflection configuration, which corresponds to r = -p and therefore k = -2p [SBK + 03]. In this case there are only two scalar functions directly accessible to the lattice. They correspond to the combination λ 1 = λ 1 -k 2 τ 3 and τ 5 . The comparison of these functions with our results is also presented in this section. Some more general kinematics have been studied on the lattice [SBK + 05] for one special combination of the scalar functions but as they present larger statistical errors we only compare some special cases. For instance, the completely symmetric configuration with r 2 = p 2 = k 2 and the case in which the quark and antiquark has the same momentum norm but the gluon momentum is arbitrary.

Vanishing gluon momentum

Let us present the result for the case of vanishing gluon momentum. In this particular configuration the tensor structure simplifies to

Γ µ (p, -p, 0) = -ig λ 1 (p 2 )γ µ -4λ 2 (p 2 ) / pp µ -2iλ 3 p µ . (6.2)
Moreover the expression for the scalar functions in this configurations simplify considerably and, in d = 4 -, they can be done analytically:

λ 1 (p 2 ) = 1 + 3N g 2 32π 2 1 -log m + g 2 192π 2 m 4 N p 4 log M m -6m 8 N 2 -1 + m 6 M 2 10N 2 -9 -14N 2 + 3 p 2 -3m 4 M 4 N 2 -14M 2 p 2 + 3 N 2 + 2 p 4 + 3m 2 M 6 -7M 4 p 2 + 15M 2 p 4 -p 6 -N 2 M 2 + p 2 4 - g 2 192π 2 m 4 N p 4 log M 2 + p 2 M 2 M 2 + p 2 3 N 2 M 2 + p 2 -3m 2 + g 2 384π 2 m 4 N p 4 Y 1 m 2 p 2 Y 1 12m 4 N 2 -1 + m 2 M 2 6 -8N 2 + m 2 p 2 19N 2 + 6 -6L 2 2m 2 -5M 2 + p 2 m 2 -M 2 + p 2 -2N 2 M 2 + p 2 2 + L 1 m 4 + 2m 2 p 2 -M 2 + M 2 + p 2 2 6m 6 N 2 -1 + m 4 M 2 3 -4N 2 + m 4 8N 2 -3 p 2 -M 4 m 2 N 2 -3 + 6m 2 M 2 p 2 + m 2 N 2 + 3 p 4 -N 2 M 2 + p 2 3 , (6.3)
where m2 = m 2 e γ /(4π) with γ the Euler constant as in (4.13),

λ 2 (p 2 ) = g 2 384π 2 m 2 N p 4 12m 4 N 2 -1 -m 2 4N 2 -3 2M 2 -p 2 -2N 2 M 2 + p 2 2 - g 2 384π 2 m 4 N p 6 log M m 12m 8 N 2 -1 -m 6 10N 2 -9 2M 2 -p 2 + 6m 4 M 4 N 2 + p 4 + m 2 -6M 6 -9M 4 p 2 + 6M 2 5N 2 -1 p 4 + 9p 6 +2N 2 M 8 + 4M 6 p 2 -9M 4 p 4 + 19M 2 p 6 + p 8 - g 2 384π 2 m 4 N p 6 log p 2 M 2 + 1 M 2 + p 2 2 m 2 3p 2 -6M 2 + 2N 2 M 2 + p 2 2 + g 2 768π 2 m 4 N p 8 L 1 Y 1 12m 10 N 2 -1 p 2 + m 8 p 2 M 2 30 -32N 2 + 22N 2 -21 p 2 + m 6 p 2 2M 4 13N 2 -9 + 3M 2 1 -2N 2 p 2 + 10N 2 -9 p 4 -m 4 150M 8 N 2 + 6M 6 N 2 + 1 p 2 + M 4 10N 2 -3 p 4 + 2M 2 2N 2 -3 p 6 + 3p 8 + m 2 300M 10 N 2 + 2M 8 3 -149N 2 p 2 + M 6 4N 2 + 15 p 4 + 9M 4 p 6 -M 2 4N 2 + 3 p 8 -2N 2 + 3 p 10 + 150M 8 N 2 Y 2 1 -2N 2 M 2 + p 2 2 75M 8 + M 6 p 2 + 3M 4 p 4 + 3M 2 p 6 + p 8 - g 2 128π 2 m 4 N p 2 L 2 m 2 + 5M 2 N 2 m 2 -M 2 + p 2 , (6.4) and λ 3 (p 2 ) = g 2 M 64π 2 m 4 N p 4 2 log M m 3m 4 N 2 -1 (m 2 -M 2 ) -N 2 p 4 m 2 -4M 2 -p 4 m 2 -p 2 (m 2 -M 2 ) m 2 -3M 2 N 2 -N 2 p 6 -6m 4 N 2 -1 p 2 - 3L 1 Y 1 m 4 N 2 -1 p 2 m 2 + M 2 + m 2 -M 2 2 -p 2 L 2 m 2 -M 2 + p 2 m 2 + N 2 p 2 -3M 2 , } (6.5)
where

L 1 = Log Y 1 -p 2 2 -M 2 -m 2 2 (p 2 + Y 1 ) 2 -(m 2 -M 2 ) 2 , L 2 = Log -m 2 + M 2 + p 2 2 -Y 2 1 (m 2 -M 2 + p 2 ) 2 -Y 2 1 and Y 1 = (m 2 -M 2 + p 2 ) 2 + 4M 2 p 2 .
In lattice simulations, these functions are obtained from the simulated vertex Γ µ (p, -p, 0) by projecting in the various tensor structures in the following way [SK02, SBK + 03]

λ 1 (p 2 ) = -1 4g B Im Tr γ µ Γ µ (p, -p, 0) pµ=0, pν =0 for µ =ν , λ 2 (p 2 ) = 1 4p 2 µ 1 4g B Im [Tr γ µ Γ µ (p, -p, 0)] + λ 1 (p 2 ) (6.6) and λ 3 (p 2 ) = 1 2p 2 µ p µ 1 4g B Re [Tr Γ µ (p, -p, 0)]
where no implicit sum over repeated indices is meant.

To compare our results with the lattice simulation of [SBK + 03] we renormalize the vertex using the IS scheme and implement the renormalization group as was presented in the section 6.4. As mentioned before, the SU (3) quenched results in d = 4 for the gluon and ghost propagators renormalized using the IS scheme were studied in section 3.3 where it was found that the best fitting parameter for the renormalization group flow at µ 0 = 1 GeV were g 0 = 4.1 and m 0 = 0.43 GeV. These values are fixed in our study of the quark-gluon vertex in the quenched approximation.

The only remaining free parameter is the initial condition of the quark mass (except for a multiplicative factor needed to match the different renormalizations). In order to fit the quark mass we start studying the scalar function λ 1 . With this purpose we define the error associated with it as

χ 2 λ 1 = 1 N i (λ 1 ) -2 lt. (p 2 i ) (λ 1 ) lt. (p 2 i ) -(λ 1 ) th. (p 2 i ) 2
and represent the error contour lines in figure 6.6. The analysis of the error shows that the fit for λ 1 improves when a smaller quark mass is considered. However the dependence on the quark mass is not so important. This fact can be seen in figure 6.8 where the comparison with the lattice data and our analytical results for M 0 = 0.02 GeV, M 0 = 0.04

GeV and M 0 = 0.08 GeV is done. The same procedure was done for the scalar function λ 3 using the error defined as

χ 2 λ 3 = 1 N i (λ 3 ) -2 lt. (p 2 i ) (λ 3 ) lt. (p 2 i ) -(λ 3 ) th. (p 2 i ) 2 .
The contour lines of χ λ 3 are represented in figure 6.7 where it can be noted that the fit improves for masses near M = 0.08 GeV. This is shown explicit in figure 6.9 where the lattice data is compared with the analytical result using this mass and also M 0 = 0.02 GeV and M = 0.04 GeV. We can also see that this scalar function depends strongly on the value of the quark mass as all functions that tend to zero in the chiral limit. Moreover, the rapid growth of the scalar function λ 3 looks like a consequence of the chiral symmetry breaking. The contourlines correspond to 4.3%, 4.4% and 4.6%.

We can conclude that the comparison of λ 1 is very good with an error less than 5% and λ 3 is also good when we use M 0 = 0.08 GeV. It is important to stress that once the fitting parameters are fixed we use the same values to reproduce all the other kinematic configurations. We consider, however, three values of the quark mass in order to see the dependence on it.

The situation is drastically different for λ 2 . Instead our analytical expression seems to be really small in comparison with the lattice data as seen in Fig. 6.10. The same discrepancy was obtained by Dyson-Swinger methods [START_REF] Aguilar | A new method for determining the quark-gluon vertex[END_REF][START_REF] Williams | The quark-gluon vertex in Landau gauge bound-state studies[END_REF]. In fact, λ 2 is not constrained by the chiral symmetry and therefore there is no reason to expect the growth found in the lattice data. We believe that this can be an artefact introduced by the way λ 2 is obtained from the lattice. We can see in (6.6) that λ 2 is computed through the addition of dλ 1 with λ2 (p 2 ) = 1 4g B Im (Trγ µ Γ µ (p, -p, 0)) = 4p 2 λ 2 -dλ 1 (which is directly measured on the lattice) divided by 4p 2 . The errors of each term sum to give the error of λ 2 which can be larger than expected. In fact, our calculations see that λ 2 tends to a constant in the infrared. This can be interpreted as a compensation of λ2 and dλ 1 in the limit p → 0 with the difference being of order p 2 . Therefore, a small error in one of these two terms would lead to a spurious divergence of λ 2 .

The divergence of λ 2 (p 2 ) observed by lattice simulations implies a non-analytic behaviour of the vertex function at small momenta. In our calculations, the only non-analytic behaviour comes from the ghost loops which do not contribute at one loop to this vertex. Therefore, our vertex at one-loop is analytic.

For all these reasons, we decided to compare instead the expression λ2 (p 2 ) which is directly obtained by the simulations and is not expressed as the subtraction of large numbers. The results are shown in figure 6.11 where the one loop results are a good approximation of the lattice data. That seems to confirm that the lattice data of λ 2 may be altered by the way it has been extracted from the lattice data.

The quark reflection kinematics

In this section we present the results for the kinematic configuration given by r = p and k = -2p which was studied in [SBK + 03] lattice simulations. In this kinematics the quark-gluon vertex also simplifies considerably to

Γ µ (p, p, -2p) = -ig λ 1 (p 2 )γ µ + τ 3 (p 2 )(/ kk µ -k 2 γ µ ) + iτ 5 (p 2 )σ µν k ν (6.7)
and the transverse projected vertex, Γ P µ (p, p, -2p) = P ⊥ µν Γ µ (p, p, -2p), has the form

Γ P µ (p, p, -2p) = -ig λ 1 (p 2 )(γ µ -/ kk µ /k 2 ) + iτ 5 (p 2 )σ µν k ν (6.8)
with λ 1 = λ 1 -k 2 τ 3 . In spite of this simplification, some of the integrals in terms of Feynman parameters can not be done analytically in this configurations.

In the lattice the scalar functions λ 1 and τ 5 in this configuration are obtained as

λ 1 = - 1 3 µ 1 4g B Im Trγ µ Γ P µ (p, p, -2p) and τ 5 = 1 3k 2 µ,ν k µ 1 4g B Re Trσ µν Γ P µ (p, p, -2p) .
In this case, we use a running scale µ = √ 3p (coherent with µ = p 2 +r 2 +k 2

2

). In figure 6.12 we compare the results with the lattice data using the same initial condition of the coupling constant and masses as in the previous case, g 0 = 4.1 and m 0 = 0.43 GeV for different values of the quark mass. We observe that the function λ 1 is well reproduced by our results. However, we do not succeed in reproducing quantitatively the behaviour of the scalar function τ 5 . In particular, our one-loop calculation does not reach such large values in the infrared. This observation was already done in the Dyson-Schwinger results of [START_REF] Aguilar | A new method for determining the quark-gluon vertex[END_REF][START_REF] Williams | The quark-gluon vertex in Landau gauge bound-state studies[END_REF], and in [START_REF] Aguilar | A new method for determining the quark-gluon vertex[END_REF] was found that they can reproduce the lattice data multiplying their results by a factor different for the factor needed to reproduce λ 1 . We also observe that this function is very sensitive to the choice of the quark mass, as it was expected from the chiral symmetry.

Completely symmetric configuration

In this section we study the particular case of the completely symmetric configuration, meaning that the norm of the three momentum are equals, p 2 = r 2 = k 2 . To do the comparison we had to interpolate the lattice data in order to obtain the results for this configuration as it was done in [SBK + 05]. The matching is done by considering the renormalization group with running scale µ = 3/2p. The fitting parameters that we use correspond to the ones that better fit the quenched propagators and they correspond to g 0 = 4.1 and m 0 = 0.43 GeV at 1 GeV while the value of the quark mass varies between 0.02 GeV and 0.08 GeV. The fit is shown in figure 6.13 where we can observe that most of the lattice data is well reproduced except for the only value which is considered in the extreme infrared where an error of the order for 25% is offered.

Symmetric configuration for arbitrary gluon momentum

The last two cases we investigate correspond to a gluon momenta fixed to k = 0.277 GeV and to k = 0.838 GeV and to equals modulus for the quark and anti-quark momentum.

As in the preceding cases the renormalization scheme used is the IS scheme and the initial conditions of the renormalization group equations correspond to g 0 = 4.1 and m 0 = 0.43 GeV at µ 0 = 1 GeV. The fits are shown in figure 6.14. The results are better for the lower gluon momentum while for k = 0.838 GeV there is an infrared discrepancy of order 25% between our results and the lattice data. One possible explanation to this fact is that we worked with the coupling constant obtained though the ghost-gluon vertex instead of using the quark-gluon vertex. As the coupling constant obtained through this vertex is larger in the infrared than the one we used it can be expected that perturbation theory requires to go to the next to leading order. . In red the one-loop results with M = 0.02 GeV (plain line), M = 0.04 GeV (large dashes) and M = 0.08 GeV (dotted line). The dots the lattice data of [SBK + 05].

The influence of dynamic quarks

In this section we analyse how the effects of dynamic quarks influence the results we have obtained. The inclusion of dynamic quarks on the lattice is cumbersome for Monte Carlo simulations and therefore no lattice data is still available for the unquenched vertex. However, in our model taking into account the dynamic quarks is really simple.

First of all it is important to note that the one-loop contribution of the quark-gluon vertex does not directly depend on the presence of dynamic quark. The first contribution of dynamic quarks come from two loops diagrams. However, the dynamic quarks modify the results through the β functions and anomalous dimensions. In fact, dynamic quarks affect the renormalization factor for the gluon Z A and for the gluon mass Z m 2 , which both have a contribution from a fermion closed loop. The unquenched effects of the RG was already studied in chapter 5 where we performed the calculations of the unquenched propagators.

In this section we compare our one-loop result of the quark-gluon vertex renormalized with the IS scheme taking into account, on one hand, the quenched analysis that we have presented in this chapter and on the other hand the influence of the unquenched effects. In particular, we compare the scalar functions in the vanishing gluon momentum case considering in the unquenched case two light fermions. The fitting parameters for the unquenched propagator with N f = 2 are g 0 = 4.5 and m 0 = 0.42 GeV (see figure 5.5). We choose to compare both cases using M 0 = 0.08 GeV which has proven to be a reasonable value for the quark mass.

In order to see how is the influence of dynamic quarks in the infrared we normalized both cases so the λ 1 functions coincide in the ultraviolet at 6 GeV. The comparison is shown in figure 6.15.

The inclusion of dynamic quarks seems to enhance the absolute value of λ 1 , λ2 and λ 3 in the infrared. We conclude that the presence of dynamic quarks enhance the quark-gluon vertex.

Conclusions and perspectives

In this thesis we studied a simple Lagrangian which allows us to do perturbative calculations in the infrared and reproduces the behaviour of the gluon propagator given by lattice simulations. This Lagrangian is a particular case of the Curci-Ferrari model [START_REF] Curci | A Theory of Massive Yang-Mills Field[END_REF] in Landau gauge. It consists on the standard Faddeev-Popov Lagrangian with an extra term that gives mass to the gluons. This model had been studied several times in [START_REF] Curci | A Theory of Massive Yang-Mills Field[END_REF][START_REF] Curci | On a Class of Lagrangian Models for Massive and Massless Yang-Mills Fields[END_REF][START_REF] De Boer | On the renormalizability and unitarity of the Curci-Ferrari model for massive vector bosons[END_REF][START_REF] Gracey | Two loop MS-bar renormalization of the Curci-Ferrari model[END_REF][START_REF] Gracey | Three loop MS-bar renormalization of the Curci-Ferrari model and the dimension two BRST invariant composite operator in QCD[END_REF][START_REF] Gracey | Three loop renormalization of QCD in various nonlinear gauges[END_REF], where it was demonstrated that it has the property of being renormalizable. Also it matches with the standard Faddeev-Popov Lagrangian at high energy which assures us that all the properties already verified by the standard perturbative QCD are reproduced also by this Lagrangian.

There are three questions to answer about this Lagrangian. The first one is if the theory is unitary. This question comes together with the definition of the physical space in which the theory has to be unitary. The standard definition of the physical space relies on the cohomology of the BRST symmetry [KO78, BRS75, BRS76, Tyu75]. Usually, this definition is well-defined because BRST symmetry is nilpotent. However, the addition of a gluon mass term breaks the BRST symmetry. BRST in no longer a symmetry of the theory based on the Lagrangian (2.11). Fortunately this model presents other symmetries studied along this thesis. The most important is called soft-BRST. Unfortunately this new symmetry is not nilpotent so the standard definition of the physical Hilbert space can not be used. The problem with the definition of the physical space is shared with the GZ approach and for the moment there is no answer to this question (see however [START_REF] Schaden | BRST Cohomology and Physical Space of the GZ Model[END_REF] for recent development).

The second question is what is the origin of this mass. We believe that this term could encode the problem of the Gribov copies. That means that if a proper definition of the gauge fixing procedure is considered, a mass term for the gluon would appear as is observed on the lattice. A possible explanation from first principles is presented in [START_REF] Serreau | Lifting the Gribov ambiguity in Yang-Mills theories[END_REF].

The third question and aim of this thesis is to verify if this model reproduces the lattice data. The first two questions become relevant once this third question is answered. To address it we used the standard perturbative analysis with Feynman diagrams. In chapter 3 we presented the computation of the quenched gluon propagator and ghost dressing function at one loop for arbitrary number of color and dimension. We renormalized the theory using two different renormalization schemes. The first scheme, called the vanishing momentum scheme presents a Landau pole as a function of the renormalization point. However, the symmetries of the theory allow us to prove some non-renormalization theorems. The extension of these theorems to the finite parts of the renormalization factors gives a renormalization scheme without a Landau pole. This scheme is called the infrared-safe scheme. In order to reproduce the lattice data we considered the effects of the renormalization group. The comparison between our analytical results and the lattice simulations has only two free parameters to fix. These parameters are the initial condition of the renormalization group for the coupling constant and the gluon mass. We choose these parameters in order to reproduce the gluon propagator and the ghost dressing function simultaneously maintaining good precision in both fits with a maximum error of ∼ 10%. The matching between our results and the lattice data in four and three dimension are very good [START_REF] Tissier | Infrared propagators of Yang-Mills theory from perturbation theory[END_REF][START_REF] Tissier | An Infrared Safe perturbative approach to Yang-Mills correlators[END_REF][START_REF] Pelaez | Three-point correlation functions in Yang-Mills theory[END_REF].

It is important to mention that the coupling constant is not entirely a free parameter. It should match with the lattice coupling constant at the ultraviolet scale of the order of the inverse of the lattice spacing. In the same chapter we verified that the coupling constant we obtained at the renormalization scale 1 GeV evolves to the expected value for higher scales up to an error of the order of 20%.

In chapter 4 we computed the one-loop quenched three-point correlation functions (corresponding to the ghost-gluon and the three gluon vertex) in Landau gauge for arbitrary dimension and number of colors in all momenta regime [START_REF] Pelaez | Three-point correlation functions in Yang-Mills theory[END_REF]. These calculations are much heavier than the previous ones because they depend on three kinematic independent variables and the tensor structure is more complex. Few analytical results were known for the infrared behaviour of the ghost-gluon vertex ([RQ11b, DORQ12, ABCP12, AIP13, HvS13] and later [START_REF] Williams | The quark-gluon vertex in Landau gauge bound-state studies[END_REF]). For the three-gluon functions some ansatz have been proposed before this calculation [START_REF] Huber | On the influence of three-point functions on the propagators of Landau gauge Yang-Mills theory[END_REF]. Our one-loop results were renormalized with the renormalization schemes already presented and also the renormalization group effects were taking into account. We used the same value of the fitting parameters as the ones obtained by fitting the two-point correlation functions where the lattice error are small and the quality of the fits are remarkably good. Therefore the comparison of the three-point functions with the lattice simulations of [START_REF] Cucchieri | Three-point vertices in Landaugauge Yang-Mills theory[END_REF] were done without any extra free-parameter (with an exception of a global multiplicative factor for the three-gluon function). It was surprising that without fitting any new parameter the fits for these quantities in d = 4 and d = 3 and for the three kinematic configurations available on the lattice were very good. We observed that the qualitative properties of these functions observed in the lattice were reproduced by our results. For instance, lattice simulations observe that the three-gluon function becomes negative in d = 3 at low momenta and it seems to diverge when all the momenta goes to zero. Our calculations showed the same behaviour which can be explained as a consequence of the infrared singular behaviour of the diagram with a ghost loop.

In the second half of the thesis we include dynamic quarks in our calculations. We started by computing the unquenched corrections for the gluon propagator and the one loop contributions for the quark propagator. These results were presented in chapter 5 and published in [START_REF] Peláez | Two-point correlation functions of QCD in the Landau gauge[END_REF]. As in the previous chapters, the calculation were done for Landau gauge, arbitrary number of colors, number of fermions and dimension. In the unquenched case we found that the gluon and ghost propagator are reproduced with high precision for d = 4 and N f = 2 and N f = 2 + 1 + 1. We also see that these quantities do not depend on small changes on the quark mass. On the other hand, the quark sector is not as well reproduced as the ghost-gluon one. The quark mass seems to have the behaviour found in lattice simulations where the chiral symmetry breaking makes the quark mass to increase quickly in the infrared even when the ultraviolet mass is small. However, the one loop contribution of the scalar function, Z u,d , associated with the tensor / p is not well reproduced, as it happens in standard QCD where the two loop contribution is dominant. For this reason, we estimate the two loop contributions using a hybrid model. This model consists on adding the ultraviolet two loop results obtained by Gracey in [START_REF] Gracey | Three loop MS-bar renormalization of the Curci-Ferrari model and the dimension two BRST invariant composite operator in QCD[END_REF] (adapted to the infrared regime) to our one loop calculations. This analysis allows us to estimate how important are the two loop corrections. We observed that for Z u,d this correction can not be neglected but it is reasonable small for the other quantities.

In the last chapter, we extended the analysis of the quark sector studying the quarkgluon vertex. We computed the one loop contribution in Landau gauge for this quantity in arbitrary dimension. For this calculation we used the package of the Mathematica software called FeynCalc which helped us to manipulate the Dirac matrices. For arbitrary kinematics the quark-gluon vertex depends on twelve scalar functions. We computed all of them but we compared only some, limited by the lattice data available. The comparison was satisfactory for most of the scalar functions, in particular we obtained very good results for λ 1 in all the configurations studied. This scalar function is directly obtained from lattice simulations and it is the one containing the tree-level structure. Our analysis was done prioritizing the fits for this quantities. However there are a pair of scalar functions that we could not quantitatively reproduce. They are λ 2 in the vanishing gluon momentum case and τ 5 in the quark reflection kinematic. We believe that the disagreement with λ 2 , also found by DS study [START_REF] Aguilar | A new method for determining the quark-gluon vertex[END_REF], may be due to the way that this function is extracted from the lattice simulations. If, instead, we compare the function which is directly obtained by lattice simulations we obtain a very good agreement. For the scalar function τ 5 , we observed that its absolute value increases in the infrared as is shown by the lattice data. However, we can not reproduce the quantitative behaviour of this quantity which strongly depends on the quark mass as is expected by the chiral symmetry. The same quantitative problem is found by [START_REF] Aguilar | A new method for determining the quark-gluon vertex[END_REF][START_REF] Williams | The quark-gluon vertex in Landau gauge bound-state studies[END_REF] where they note that for reproducing this scalar function an extra ad-hoc multiplicative factor is needed. The problem is that this multiplicative factor is not coherent with the matching of λ 1 .

Independently of this thesis, the investigation of the infrared properties including a gluon mass term was extended to finite temperature. In [START_REF] Reinosa | Yang-Mills correlators at finite temperature: A perturbative perspective[END_REF] it is shown that the properties of the gluon and ghost propagators qualitative reproduced at finite temperature. The same framework was used to study the deconfinement transition [START_REF] Reinosa | Deconfinement transition in SU(N) theories from perturbation theory[END_REF][START_REF] Reinosa | Deconfinement transition in SU(2) Yang-Mills theory: A two-loop study[END_REF]. It was observed that second order transition for SU (2) and the first order transition for SU (3) are found using a one-loop calculation. Moreover if a nextto leading order calculation is performed in SU (2) the quantitative value of the critical temperature is obtained. This work can be extended for the SU (3) and also including the effects of the renormalization group. The study adding quarks at non zero chemical potential was also done [START_REF] Reinosa | Perturbative study of the QCD phase diagram for heavy quarks at nonzero chemical potential[END_REF].

All these results support the idea that at least an important part of the infrared effects present in QCD can be reproduced with a gluon mass term. This makes us think that this mass encode the effects of Gribov copies [Gri78, Zwa89, Zwa93, DGS + 08a]. Contrarily to the previous analysis of the infrared effects of Gribov copies, our approach does not require extra fields. It has the advantage that the Feynman rules remains practically unchanged except for a mass term in the gluon propagator and therefore the calculation of many quantities become treatable. Therefore our analysis can be extended for different situations, for instance we have in mind to compute more physical quantities like observables (gauge invariant quantities) such as the glueball masses [START_REF] Cornwall | Glueballs as Bound States of Massive Gluons[END_REF]. Moreover, in order to check completely the validity of the perturbative expansion we would like in the future to compute also some two-loop corrections. The first task in this direction is to compute the quenched propagators. It would be also interesting to study deeper the chiral symmetry breaking and the possibility of explaining it through our model. We have already mentioned some of the fundamental open questions produced with this modification of the Lagrangian. A big task in our future is to investigate the possible answer for those questions including the definition of a the physical space and the origin of the gluon mass. From this definition it can be shown that

δΓ[φ] δφ i (z) = J i (z)
where φ i = φ i = δW [J] δJ i . If another derivative δ δJ j (x) is taken we find δ δJ j (x)

δΓ[φ] δφ i (z) = δ ij δ(z -x) = dy δφ k (y) δJ j (x)
δ 2 Γ(φ) δφ k (y)δφ i (z) .

We can insert the definition of φ dy δ 2 W [J] δJ j (x)δJ k (y)

δ 2 Γ(φ) δφ k (y)δφ i (z) = δ ij δ(z -x)
where we can identify the propagator, then dy φ j (x)φ k (y) c δ 2 Γ(φ) δφ k (y)δφ i (z) = δ ij δ(z -x).

The last equation allows us to affirm that φ i (x)φ j (y) c = δ 2 Γ(φ) δφ i (x)δφ j (y) -1 .

(A.2)

To obtain the tree-level propagator for the Feynman rules one has to differentiate the action twice and take the inverse. In order to find the relation between full connected Green's functions and the vertex we will take the derivative δ δJ k (z) of the above equation. There is another symmetry, similar to the gauge symmetry, presented in [TW09a] which we will not present here. The same identities were obtained in [START_REF] Binosi | Anti-BRST symmetry and background field method[END_REF] and related to the background gauge formalism.

B.2 Equations of motion

Associated with h a

We consider the transformation of h a given by ih a → ih a + a (x), with an arbitrary function. The expected value of the variation of the action has to be zero, so The effective action has the same dependence on h a as the bare action. In Landau gauge we do not need an independent renormalization factor for h a . It is worth noting that there are only two two-point vertex involving h a non zero a when the limit of vanishing field and sources is taken. They are

Γ βc Γ β K Γ Kc Γ k β Γ cc Γ L L (B.16)
The derivatives of Slavnov-Taylor identities and the equations of motion allow us to relate all the two-point vertex functions between them. The final expressions will be evaluated in zero field and sources so they the ghost number should be zero. Moreover they should respect Lorentz invariance.

The ghost number of each field is field

A µ c c h J µ χ χ R K µ Kµ L L M α µ β β Ghost number 0 1 -1 0 0 -1 1 0 1 -1 2 -2 0 0 1 -1
We present a list of all possible derivatives of Slavnov-Taylor identities and the equations of motion that relate the two-point vertex function. Some of them can be reach from more than one way.

Figure 2

 2 Figure 2.1 -Dyson-Schwinger equation for ghost propagator, original figure from [AHS10]. The bold line refers to the full propagator while the biggest point is the full vertex.

Figure 2 Figure 2

 22 Figure 2.2 -Dyson-Schwinger equation for gluon propagator, original figure from [AHS10]. The bold line refers to the full propagator while the biggest point is the full vertex.

  Figure 2.4 -Lattice

Figure 2

 2 Figure 2.5 -Gluon propagator, lattice data from [CMM08].

  Figure 2.6 -Ghost dressing function for d = 4, lattice data from [CMM08].

Figure 3 . 1 -

 31 Figure 3.1 -One-loop correction for the ghost propagator.

Figure 3 . 2 -

 32 Figure 3.2 -One-loop correction for the two-point gluon vertex function.

Figure 3

 3 Figure 3.3 -Snail

Figure 3 . 4 -Figure 3

 343 Figure 3.4 -One loop diagram with a ghost loop.

Figure 3 . 6 -Figure 3 . 7 -

 3637 Figure 3.6 -Four dimensional gluon propagator for SU (2). In red the results of Eq.(3.21) renormalized in the VM scheme without including the renormalization group flow. The value of the parameters used are g = 7.5 and m = 0.68 GeV and µ = 1 GeV. The gray dots are the lattice data from [CM08a].

Figure 3

 3 Figure3.8 -Four dimensional SU (2) ghost dressing propagator (top) and gluon propagator (bottom) as a function of momentum renormalized in the VM scheme by including the renormalization group effects. The points are lattice data of[START_REF] Cucchieri | Constraints on the IR behavior of the gluon propagator in Yang-Mills theories[END_REF]. The large dashed line to α = 1.0, while the small and tiny dashed curves correspond to α = 2.0 and α = 3.0 respectively. The initial conditions of the RG flow are the presented in table 3.1.

  Scheme α g 0 m 0 (GeV) VM 1.0 7.5 0.77 VM 2.0 9.0 0.78 VM 3.0 9.1 0.75 Table 3.1 -Fitting parameters retained for computing correlation functions in d = 4 for different schemes.

Figure 3 . 9 -Figure 3

 393 Figure 3.9 -Contour levels for the quantities χ AA and χ cc for the VM scheme in d = 4 for α = 1. The large diagonal region corresponds to χ cc and the small elliptic one to χ AA . From dark to light: 4%, 7% and 10%.

Figure 3

 3 Figure 3.11 -Gluon propagator (top) and ghost dressing function (bottom) as a function of momentum in d = 4. The plain blue line corresponds to the infrared safe scheme with α = 0 considering the effects of the renormalization group. The initial condition of the parameters corresponds to g 0 = 5.2 and m 0 = 0.44 GeV at µ 0 = 1 GeV. The points are lattice data of [CM08a].

Figure 3 . 12 -Figure 3 Figure 3

 31233 Figure3.12 -Contour levels for the quantities χ AA and χ cc for the IS scheme for d = 4. The large diagonal region corresponds to χ cc and the small elliptic one to χ AA . From dark to light: 4%, 7% and 10%.

Figure 3

 3 Figure 3.16 -Four dimensional gluon propagator and ghost dressing function for SU (3). In blue the results of Eq.(3.21) renormalized in the IS scheme including the renormalization group flow using µ = p. The value of the parameters used are g 0 = 4.1 and m 0 = 0.43 GeV and µ 0 = 1 GeV. The gray dots are the lattice data from [BIMPS09].

  Figure3.17 -Three dimensional renormalize ghost dressing function for SU (2). In red the results of Eq. (3.31) without including the effects of the RG using g = 3.7 √ GeV, m = 0.89 GeV and µ = 1 GeV. The gray dots are the lattice data from[START_REF] Cucchieri | Constraints on the IR behavior of the gluon propagator in Yang-Mills theories[END_REF]

Figure 3 Figure 3 Figure 3

 333 Figure 3.19 -Gluon propagator (top) and ghost dressing function (bottom) as a function of momentum in d = 3 considering the effects of the RG. The points are lattice data of [CM08a]. The plain blue line corresponds to the infrared safe scheme with α = 0,g 0 = 2.4 √ GeV and m 0 = 0.55 GeV were fixed at µ 0 = 1 GeV. The dashed red line corresponds to the vanishing momentum scheme with g 0 = 4 √ GeV and m 0 = 1 GeV at µ 0 = 1 GeV.

Figure 3 .

 3 Figure 3.22 -Contour levels for the quantities χ AA and χ cc for the IS scheme (left) and VM scheme (right) both for d = 3 . The large diagonal region corresponds to χ cc and the elliptic one to χ AA . From dark to light: 4%, 7% and 10%.

Figure 3

 3 Figure 3.23 -Running coupling constant in the IS scheme with initial conditions at µ 0 = 1 GeV corresponding to g 0 = 5.2 and m 0 = 0.44 GeV in d = 4 (left) and g 0 = 2.4 √ GeV and m 0 = 0.55 GeV in d = 3 (right) for SU (2)

Figure 3

 3 Figure 3.24 -Expansion parameter u(p) (full line) and the effective expansion parameter u(p) p 2 p 2 +m 2 , in d = 4 (left) and d = 3 (right)

∼ ∼ p1p2g 2 m 2 p 1 p 2 Figure

 2 Figure 3.25 -Effective vertex

Figure 4

 4 Figure 4.1 -One-loop Feynman diagrams for the ghost-gluon vertex.

Figure 4

 4 Figure 4.2 -Description of the diagram, D1 ccA , with three ghost-gluon vertices.

Figure 4

 4 Figure 4.3 -Description of the diagram, D2 ccA , with one three-gluon vertex and two ghost-gluon vertices.

  Kc ν using the diagrams associated with the one loop contributions (see diagrams of figure 4.4, where the tree-level vertex for Γ

Figure 4 . 4 -

 44 Figure 4.4 -Description of the one loop diagrams of the three ghost-K-gluon vertex. Here, the left bottom wavy external leg correspond to Kc ν .

Figure 4

 4 Figure 4.5 -One-loop Feynman diagrams for the 3-gluon vertex.

[ 2 ]Figure 4 . 6 -Figure 4

 2464 Figure 4.6 -Description of the diagram with three ghost-gluon vertices. In the left D1 [1] AAA

Figure 4

 4 Figure 4.8 -Description of the diagram with three three-gluon vertices, D3 AAA .

Figure 4 Figure 4 Figure 4 Figure 4 . 12 -

 444412 Figure 4.9 -Ghost-antighost-gluon correlation function G ccA for one vanishing momentum (top figure), all momenta equal (middle figure), two momenta orthogonal, of equal norm (bottom) as a function of momentum, in d = 4. The lattice data of [CMM08] are compared with our calculations. The plain blue line corresponds to the IS-scheme with α = 0 and the dashed red line corresponds to the VM-scheme with α = 1.

Figure 5 . 1 -Figure 5 . 2 -

 5152 Figure 5.1 -Diagrams contributing to the gluon self-energy.

Figure 5 . 3 -

 53 Figure 5.3 -One loop contribution to the quark two-point vertex

Figure 5 . 5 -Figure 5 . 6 -

 5556 Figure 5.5 -Gluon dressing function (top) and ghost dressing function (bottom) as a function of momentum in d = 4 for N f = 2. Using g 0 = 4.5, m 0 = 0.42 GeV and M u,d 0 = 0.13 GeV at 1 GeV. The points are lattice data of [SMMPvS12].

Figure 5

 5 Figure 5.7 -Contour levels for the quantities χ AA,Unq (forward) and χ cc,Unq (backward) for d = 4, for N f = 2 + 1 + 1 and M u,d 0 = 0.13 GeV. The contour lines correspond to 3%, 4% and 5% for the gluon and 1%, 1.5% and 2% for the ghost.

Figure 5 . 8 -Figure 5

 585 Figure 5.8 -Gluon dressing function (top) and ghost dressing function (bottom) as a function of momentum in d = 4 for N f = 2 + 1 + 1. Using g 0 = 5.3, m 0 = 0.48 GeV and M u,d 0 = 0.13 GeV at 1 GeV. The points are lattice data of [ABB + 12].

Figure 5

 5 Figure 5.10 -Gluon dressing function (left) and ghost dressing function (right) for different number of flavours in d = 4.

Figure 5 . 12 -

 512 Figure 5.12 -Gluon dressing function (left) and ghost dressing function (right) as a function of momentum in d = 4 for N f = 2+1. Using g 0 = 4.8, m 0 = 0.42 GeV and M u,d 0 = 0.13 GeV at 1 GeV. The points are lattice data of [BHL + 04, ABB + 12] respectively.

Figure 5

 5 Figure5.13 -Contour levels for the quantities χ AA Unq (left) and χ CC Unq (right) for d = 4 with g 0 = 4.8. The contour lines correspond to 6%, 7% and 8 % for the gluon and 2.5%, 3% and 4 % for the ghost.

Figure 5

 5 Figure 5.14 -The scalar function M (p) for the light quark as a function of momentum in d = 4 for N f = 2 + 1. Using g 0 = 4.8, m 0 = 0.42 GeV and M u,d 0 = 0.13 GeV at 1 GeV. The points are lattice data of [BHL + 05].

Figure 5

 5 Figure 5.15 -Contour levels for the quantity χ M . The contour lines correspond to 10%, 20% and 30 %.

Figure 5

 5 Figure 5.16 -Gluon dressing function (left) and the scalar function M (p) for the light quark (right) as a function of momentum in d = 4 for N f = 2 + 1. Using g 0 = 9.5, m 0 = 1.4 GeV and M u,d 0 = 0.13 GeV at 1 GeV. The points are lattice data of [BHL + 04, BHL + 05].

Figure 5 Figure 5 Figure 5 Figure 5 Figure 5 Figure 5

 555555 Figure5.17 -Contour levels for the quantities χ M (left) and contour levels for the quantities χ M and χ AA (right) for d = 4, for N f = 2 + 1, maintaining g 0 = 4.8 at 1 GeV. The contour lines correspond to 7%, 8% and 10 % for the gluon and 50%, 70% and 90 % for the quark mass.

Figure 6

 6 Figure 6.1 -Quark gluon vertex.

Figure 6 . 2 -Figure 6 . 3 -

 6263 Figure 6.2 -One loop contribution to the quark-gluon vertex.

Figure 6

 6 Figure 6.4 -One loop contribution to the quark-gluon vertex; diagram with a three-gluon vertex.

Figure 6 . 5 -

 65 Figure 6.5 -Unrenormalized scalar functions as function of p for k = 0.555 GeV and r = 0.844 GeV.

Figure 6

 6 Figure6.6 -Contour levels for the quantity χ λ 1 for d = 4, maintaining g 0 = 3.7 at 1 GeV. The contourlines correspond to 4.3%, 4.4% and 4.6%.

Figure 6 . 7 -Figure 6 . 8 -Figure 6 Figure 6

 676866 Figure6.7 -Contour levels for the quantity χ λ 3 for d = 4, maintaining g 0 = 3.7 at 1 GeV. The contourlines correspond to 8.8%, 9% and 9.2%.

Figure 6

 6 Figure 6.11 -The figure shows the function λ2 (p 2 ) for d = 4. The initial conditions of the RG equation correspond to g 0 = 4.1, m 0 = 0.43 GeV at µ 0 = 1 GeV. The red line are the one-loop results with M = 0.02 GeV (plain line), M = 0.04 GeV (large dashes) and M = 0.08 GeV (dotted line). The dots the lattice data of [SBK + 03].

Figure 6 . 12 -

 612 Figure 6.12 -Scalar function of the quark-gluon vertex: λ 1 (top) and τ 5 (bottom) for the quark reflection configuration. The red line are the one-loop results with M = 0.02 GeV (plain line), M = 0.04 GeV (large dashes) and M = 0.08 GeV (dotted line). The dots the lattice data of [SBK + 03].

Figure 6

 6 Figure 6.13 -Scalar function of the quark-gluon vertex λ1 = λ 1 -k 2 τ 3 with k 2 = r 2 = p 2 .The initial conditions of the RG equation correspond to g 0 = 4.1 and m 0 = 0.43 GeV at µ 0 = 1 GeV while the running scale was taken as µ = 3/2p. In red the one-loop results with M = 0.02 GeV (plain line), M = 0.04 GeV (large dashes) and M = 0.08 GeV (dotted line). The dots the lattice data of [SBK + 05].

Figure 6

 6 Figure 6.14 -Scalar function of the quark-gluon vertex λ 1 = λ 1 -k 2 τ 3 for r 2 = p 2 and k = 0.277 GeV (top) and k = 0.838 GeV (bottom). The initial conditions of the RG equation correspond to g 0 = 4.1 and m 0 = 0.43 GeV at µ 0 = 1 GeV while the running scale was taken as µ ∼ p 2 +r 2 +k 2 2

Figure 6

 6 Figure 6.15 -Scalar functions of the quark-gluon vertex: λ 1 (top), λ2 (middle) and λ 3 (bottom) for a vanishing gluon momentum, as a function of the quark momentum. The full line represent the quenched results using initial conditions of the RG equation g 0 = 4.1, m 0 = 0.43 GeV and M = 0.08 GeV at µ 0 = 1 GeV. The dashed line represent the unquenched results for N f = 2 using g 0 = 4.5, m 0 = 0.42 GeV and M = 0.08 GeV at µ 0 = 1 GeV.

  If we perform the Legendre transformation of W [J] we construct Γ, the 1-PI generating function Γ[ φ ] = -W [J] + d 4 xJ i (x) φ i (x) .

1 =a m 2 c a - δΓ δ βa -m 2 α a µ δΓ δ Ka µ +

 1µ For example,δ 3 W [J] δJ k (z)δJ i (x)δJ j (y) = δ δJ k (z) δ 2 Γ(φ) δφ i (x)δφ j (y) -1 = d d w δφ l (w) δJ k (z) δ φ l (w) δ 2 Γ(φ) δφ i (x)δφ j (y) Using the derivative of the inverse of a matrix, δ ∂x M -1 (x) = -M -1 ∂M ∂x M -1 , we get δ 3 W [J] δJ k (z)δJ i (x)δJ j (y) = -d d w δφ l (w) δJ k (z) d d ud d v δ 2 Γ(φ) φ i (x)φ m (u) -1 δ 3 Γ(φ) δφ l (w)δφ m (u)δφ n (v) δ 2 Γ(φ) φ n (v)φ j (y) --d d wd d ud d v φ l (w)φ k (z) c φ i (x)φ m (u) c φ j (y)φ n (v) c δ 3 Γ(φ) δφ l (w)δφ m (u)δφ n (v) (A.4)To sumarize,φ k (z)φ i (x)φ j (y) c = -dwdudv φ l (w)φ k (z) c φ i (x)φ m (u) c φ j (y)φ n (v) c φ l (w)φ m (u)φ n (v) 1PI(A.5) Similar identities can be obtained for general correlation functions. which in terms of Γ takes the formd d xβa m 2 δΓ δ La -2m 2 δΓ δM a β a = 0. (B.10)For the anti-soft BRST the corresponding identity isδ βa -m 2 c a -2M a δΓ δβ a + m 2 α a µ δΓ δK a µ -βa m 2 ih a -2m 2 δΓ δL a β a = 0 (B.11)It is worth mentioning that the identity (B.10) reduces to the standard Slavnov-Taylor identity for the BRST symmetry (2.16) when K µ = L = M = α = β = β = 0.

d

  d x -a ∂ µ A a µ -iR a a + α a µ D µ ( a ) + βa (-gf abc b c c ) -gf abc b cc β a = 0As this is valid for all we have-∂ µ A a µ -i δΓ δh a -D µ α a µ -gf abc βc c b -gf abc cb β c = 0.Therefore, the explicit dependence of the vertex function Γ with respect to the field h a isδΓ δh a = i ∂ µ A a µ + D µ α a µ + gf abc ( βc c b + cb β c )

  δΓ δA b µ (y)δh a (x) = iδ ab ∂ µ (δ(x -y)) (B.12) and δΓ δα b µ (y)δh a (x) = iδ ab ∂ µ (δ(x -y)). (B.13)

Table 1

 1

				.1 -Feynman rules.
	ghost propagator	D ab 0 (p) = δ ab p 2			p
	gluon propagator	D ab 0,µν = δ ab 1 p 2 P ⊥ µν (p)	p
	quark propagator	D ab 0,ψ	p 2 +M 2 ψ(p) = δ ab i / p+M	p
						µ, a
	ghost-gluon vertex	Γ (3),tree-level ccA	bca µ	(r, p, k) = -igf abc p µ	k
					p	r
	quark-gluon vertex	Γ

  .14.

	Scheme	α	g m (GeV)
	IS	0.0 5.2	0.44
	IS	1.0 5.2	0.43
	IS	2.0 5.8	0.48
	IS	3.0 6.3	0.53
	Table 3.2 -Fitting parameters retained for computing correlation functions in d = 4 for
	different schemes.

δΓ δM c + ih c (B.15)B.3 Two-point vertex functionsAll the two-point vertex functions areΓ AA Γ Ah Γ αh Γ Aα Γ AM Γ αM Γ M M Γ αα Γ β β Γ K K Γ c β Γ c K

Acknowledgements

where, as in the previous diagram, we have cancelled the trivial vanishing projections and used f amn f bln f lcm = N 2 f abc . Thanks to momentum conservation we can express the diagrams as functions of only two momenta, r and p.

These expressions are the input for the Mathematica algorithm which performs the simplification mentioned in the Eq. (3.2), does the Feynman trick, integrates over the momentum variable and performs the corresponding integrals (numerically or analytically when possible) over the Feynman parameters. One difference with the calculation of the propagators is that the Feynman trick brings two different Feynman parameters.

In Fourier space, the bare ghost-gluon vertex including tree-level and one loop contributions then reads

The tensor structure of the ghost-gluon vertex is rather simple. As we have seen in the one loop calculation the color dependence is proportional to the structure constant f abc , where a, b, c are respectively the color indices of the gluon, ghost and antighost external legs. This is the exact color structure in the case of SU (2) and for large N . For general SU (N ) this is known to be true at one-and two-loop level [START_REF] Davydychev | Two loop three gluon vertex in zero momentum limit[END_REF]. However, in a general situation other color structures could appear for instance the completely symmetric color structure d abc .

We would like to describe the vertex tensor structure explicitly in momentum space. The ghost-gluon vertex has two independent tensor structure. However, it is useful to relate this vertex to Γ The same Ward identity gives a similar relation for Γ

(2) A a µ c b cc that matches with the derivative of Eq. (4.1) with respect to A a µ . This expression takes the form

Assuming the color dependence discussed above (that is valid at one loop), the threepoint vertex function Γ This justify the decomposition introduced by Ball and Chiu in [START_REF] Ball | Analytic Properties of the Vertex Function in Gauge Theories. 2[END_REF] where the ghostgluon vertex is expressed as a function of a rank two tensor Γ νµ ,

where the Lorentz structure can be explicitly written using scalar functions a, b, c, d, e as Γ νµ (r, p, k) = δ µν a(r 2 , p 2 , k 2 ) -k ν r µ b(r 2 , p 2 , k 2 ) + p ν k µ c(r 2 , p 2 , k 2 ) + k ν p µ d(r 2 , p 2 , k 2 ) + p ν p µ e(r 2 , p 2 , k 2 ). (4.5)

Chapter 4. Quenched three-point vertex functions consider the renormalization group, the vertex functions renormalized a different scales are related as Γ

(3)

and Γ

(3)

Therefore the quantities, G ccA (p, k, r) and G AAA (p, k, r), that we want to compare are related with our calculations as follows

where we have used the abbreviation tl to note that only the tree-level contribution is considered and therefore we can conclude that no z(µ) factor has to be consider for the quantity G ccA . On the contrary, G AAA can be expressed as

.

To compare our results with the lattice data we need to fix the multiplicative factor and two other parameters. These last two correspond to the value of the coupling constant and mass at the renormalization scale µ 0 , let us call them g 0 and m 0 . We are going to use the same renormalization scheme introduced in the previous chapter, the VM scheme (3.12) and the IS scheme (3.16). The coupling constant and the gluon mass will then run obeying the renormalization group equations.

The best fitting parameters were obtained in the study of the quenched propagators in chapter 3 and we will use the same values in our study of the vertices.

d=4

In this subsection we compare the renormalized scalar functions G ccA and G AAA with the lattice data simulated for SU (2) in d = 4 and published in [START_REF] Cucchieri | Three-point vertices in Landaugauge Yang-Mills theory[END_REF]. The value of the coupling constant and mass at 1 GeV are g 0 = 7.5 and m 0 = 0.77 GeV with α = 1 for the vanishing momentum scheme and g 0 = 5.2 and m 0 = 0.44 GeV with α = 0 for the infrared safe scheme. In figure 4.9 we show the comparison for the ghost-gluon factor, G ccA , for different momentum configurations. It is a non-trivial result that the set of parameters chosen for describing the two-point functions simultaneously give a very good agreement for the three-point functions.

Appendix A

Feynman Rules

A.1 Relation between the vertex and the connected correlation function

The generating functional, Z[J], is analogue to the partition function defined as

where φ i denotes all the field involved. The correlation function

is obtained by successive derivatives of the generating functional

where the first derivative was taken with respect to the source J i 1 and so on. Generally it is convenient to define W [J] as Z[J] = e W [J] . Moreover, W [J] is the generating function of the connected Green function because

.

The relation between the Green function and the connected ones is

The propagator between φ i (x) and φ j (y) is

.

A.2 Feynman Rules

In this section we present the derivation of Feynman rules for the Curci-Ferrari model in the Landau gauge. The action considered in a d-dimensional Euclidean space is

From now on, the Fourier transform convention chosen is

A.2.1 The ghost propagator

The tree-level ghost propagator is:

where the first derivative is with respect to cb (y) and the second one with respect to c a (x). It is important to remark that these are Grassman variables so the order is important. In order to derive the Feynman rules we consider the second derivative of the action with respect to the ghost and the anti-ghost, which is

Making the Fourier transformation, we have:

A.2.2 The gluon propagator

The second ingredient for the Feynman rules is the gluon propagagor. To compute it let us focus on the gluon quadratic term of the Lagrangian,

To obtain the gluon propagator we have to compute the second derivative of this Lagrangian and take the inverse. The derivative is given by

In the momentum space we get

where P T µν (p) is the transverse projector P T µν (p) = δ µν -pµpν p 2 . We can not ignore that there is another quadratic term in h a and A: L hA = ih a ∂ µ A a µ . The derivative of the action with respect to those fields is

which in Fourier space takes the form

The matrix for the second derivatives of the action in terms of A and h is:

We would like to take the inverse of this matrix, in order to doing so we consider an ansatz for the inverse given by

Imposing the product of the ansatz with the second-derivative matrix to be the delta function we obtain

We conclude that the solution is

Therefore the gluon propagator is

It is worth noting that the gluon propagator is transverse even though Γ (2) is not. We can note that h a A b µ and h a h b are not zero, however, as no vertex involving the field h a exists they do not play any role in the perturbative calculations.

A.2.3 The quark propagator

The quark propagator is ψ(y)

. Therefore we are going to compute the second derivative of quadratic part of S. It is easy to see that

Taking the Fourier transformation we obtained

Therefore the propagator in the Fourier space is

where the momentum goes in the opposite direction of the charge.

A.2.4 The quark-gluon vertex

The third derivative of the action with respect to anti-quark, a quark and a gluon is

When we take the Fourier transformation we obtain

Therefore, the quark-gluon vertex is igt a γ µ .

A.2.5 The ghost-gluon vertex

We want to calculate

where the first derivative is taken with respect to A a µ (z), then the second one respect to cc (y) and the last one with respect to c c (x).

Performing the Fourier transformation e -ipy e -irx e -ikz

We obtain that the ghost-gluon vertex with all momenta incoming, shown in figure A.1, is given by

Therefore, the tree-level ghost-gluon vertex with all incoming momenta is given by

A.2.6 The three-gluon vertex

Let us consider the terms in the Lagrangian that involves only the interaction of three gluons.

The three-gluon vertex is diagramaticaly represented as it is shown in figure A.2. For calculate the vertex we take the third derivative:

Once all the interested derivatives are done we make the Fourier transform considering all momenta incoming to the diagram as is shown in the picture.

With this convention of signs we arrive at:

A.2.7 The four-gluon vertex

For derive the four-gluon vertex let us focus in the only term of the action containing four gluons

The first step in order to deduce the four-gluon vertex is to take the derivative with respect to A a µ (x). The symmetry of S AAAA allows us to write this derivative as

The second derivative is given by

To perform the last two derivatives it is better to take them together so

As the above expression does not depend on the momenta the Fourier transform is easier to calculate, then, grouping conveniently the terms we obtain

Symmetries and identities B.1 Slavnov-Taylor identities

The aim of this section is to present the group of symmetries of the theory and some of the consequently Slavnov-Taylor identities that we are going to use as checks of our calculations. For simplicity we are going to focus on the quenched case. As we mentioned above the Lagrangian of Eq. (2.11) does not respect the nilpotent BRST symmetry but it is invariant under the soft-BRST symmetry, which for an arbitrary field φ is φ → φ + ηsφ where η is a Grassmann variable, and under the soft-anti-BRST, φ → φ + ηsφ, defined as

In order to study the equation of motion of this theory and the implications of soft-BRST symmetry it is better to first introduce other simpler symmetries. The ghost-charge conjugation given by

while the other fields are invariant. Another symmetry is the symplectic symmetry [START_REF] Delduc | A Note on Some Nonlinear Covariant Gauges in Yang-Mills Theory[END_REF] with generators t, t and N defined as

3)

The algebra of these symmetries is given by the commutation relations

The generator N is nothing but the counting of the number of ghost. The Slavnov-Taylor identities associated with these symmetries were already studied for the explicitly symmetric form of the Lagrangian in [START_REF] Tissier | Gauged supersymmetries in Yang-Mills theory[END_REF]. To recover the analysis in the asymmetric Lagrangian of the form of (2.11) we have to change ih a sym → ih a Asym + g 2 f abc cb c c . Then the Slavnov-Taylor identity can be obtained by considering the sources J a µ , χa , χ a , R a , Ka µ , K a µ , La , L a , M a , α a µ , βa and β a coupled to A a µ , c a , ca , h a , sA a µ , sc a , sc a , sc a + sc a , ssA a µ , ssc a and ssc a respectively. The number of sources is obviously larger than for studying standard BRST.

The action of the sources is

It is important to note that those are all the sources needed. First of all, the variations associated with t and t are already considered, due to the fact that t(ih a ) = sc a and t(ih a ) = -sc a . Secondly, it is not necessary to consider sources for s 2 ,s 2 or ss + ss. This relies on the equalities s 2 = m 2 t, s2 = m 2 t and {s, s} = -m 2 N . Third, in relation to the sources for sc a and sc a we can consider a source for sc a -sc a which is already included in R a and another for sc a + sc a . And last, ss(ih a ) = sc a , therefore it is not necessary to consider an extra source for this term.

To write Slavnov-Taylor identities it is better to consider the effective action Γ, with

which can be related to the expected value of the fields as Doing some straightforward algebra we can computed the Slavnov-Taylor identity for the t-symmetry, which yields

B.2.1 Associated with ca

To deduce the equation of motion for the anti-ghost we are going to consider the transformation ca → ca + ¯ a (x).

In consequence,

The above relation holds for all a (x) so

B.2.2 Associated with c a

A similar procedure can be done considering the transformation for the ghost field c a → c a + a (x) in addition to ih a → ih a -gf abc b cc The variation of the action is given by

Therefore the equation of motion is

From the equation of motion related to ca

From the equation of motion of c a

From the Slavnov-Taylor equation from gauge transformation

From the Slavnov-Taylor equation of t-symmetry

From the Slavnov-Taylor equation of s-symmetry

B.4 The non-renomalization theorem for the mass

One of the important consequence of the relations presented in the preceding section is the non-renormalization theorem for the mass. This implies, in addition with the Taylor scheme, that to renormalize the theory, in the quenched approximation, we only need to define two renormalization factors.

To deduce the non-renormalization theorem let us start observing the first equation of (B.17) which can be expressed as

where Γ c a Kb µ (p) can be written as Γ c a Kb µ (p) = δ ab p µ f (p) because it is the only tensor structure allowed. Therefore, we can invert the relation using that

On the other hand the second equation of (B.21), deduced from the Slavnov-Taylor identity associated with the soft-BRST symmetry, is

Using the assumption that the gluon two-point vertex is sufficiently regular at low momentum one must have that Γ ⊥ (0) = Γ (0). Therefore, Γ A a µ A b ν (0) = δ ab δ µν Γ ⊥ (0). Putting all the information together we obtain

where we have introduced the subindex B to make explicit that we are dealing with bare quantities. If we change to the renormalized fields it holds that

Therefore, the renormalization of the mass is not independent of the renormalization of the ghost and gluon fields. To absorb the divergences it is not necessary to include an independent renormalization factor of the gluon mass, that is why we called it a nonrenormalization theorem.

B.5 Three-point vertex functions

We extended our study made for the two-point vertex functions to the three-point vertex functions by finding automatically all the relations given by the symmetries that involves the two and three-point vertex functions. The number of relations is really large and therefore they are not presented here. However, some of them have a particular interest and will be described in detail.

B.5.1 Taylor non-renormalization theorem

In this thesis we have used the Taylor scheme based on the relation √ Z A Z c Z g = 1 to renormalized the coupling constant (also known as minimal MOM scheme [START_REF] Lorenz Von Smekal | The Strong coupling and its running to four loops in a minimal MOM scheme[END_REF]). This relation can be hold due to the non-renormalization theorem for the coupling constant shown by Taylor in [START_REF] Taylor | Ward Identities and Charge Renormalization of the Yang-Mills Field[END_REF] where it was proven that the coupling constant can be renormalized using the ghost and gluon renormalization factors. The theorem is a consequence of the equation of motion for c. Let us start from the equation (B.15). If we take the derivative of this equation with respect to cb (y) follow by the derivative with respect to A c ν (z) we obtain the relation:

We can use equation (B.12) to replace

In the Fourier space this last relation become

If we set the ghost momentum p to zero the relation reads

where we have include the subindex B to remember that we are working with the bare quantities.

If we renormalized this expression, we obtain

As the renormalized quantities are finite the divergent part of the renormalization factor have to fulfil the relation

As the tree level contribution trivially verifies this relation, it is possible to extend this property to all orders of perturbation theory. Therefore, Taylor scheme propose to extend this relation also for the finite parts of the renormalization factors.

B.5.2 Relation between the ghost-gluon vertex and the three-gluon vertex

The ghost-gluon and the three-gluon vertex are related to each other as a consequence of the Slavnov-Taylor equation (B.10). In order to prove this we take the derivative of equation (B.10) first with respect to c b (y) and then with respect to A c ν (z) and A d ρ (w). This procedure takes as to the following relation

which in the Fourier space can be written as

We can use the color structure defined in (3.9), (4.3) and (4.11) and the relation (B.22) to simplify the above equation, this gives

This identity is the same identity obtain in standard QCD already known in the literature.

B.5.3 Constraint for three-gluon vertex

There is another interesting relation arriving from the Slavnov-Taylor identity which makes a constraint for the three-gluon vertex. This relation is obtained by taking the derivative of Eq. (B.10) with respect to two ghosts and one antighost and taking the Fourier transformation. This leads to:

The second step consist in taking the derivative of the Slavnov-Taylor identity associated with the t symmetry Eq. (B.7) with respect to a gluon and twice with respect to the ghost field. If we take the Fourier transform we obtain

Now we can obtain Γ

(3)

26) and substitute it in Eq. (B.27). These steps lead us to the identity:

This new identity is used to check our calculation for the one-loop three-gluon vertex.

Appendix C

Reduction to Passarino-Veltman integrals.

The analysis of the diagrams led us to integrals of the form

where p 1 and p 2 represent two independent momenta and the m i takes the value 0 or m in our case. At this moment, we can choose to follow the techniques mentioned in the previous chapter. That means, to apply the Feynman trick, shift the momentum q to l, eliminate the odd powers of l and perform the integral over l. However, when the integral over the momentum space can be done the remaining integrals over the Feynman parameters are cumbersome. The difficulty increases with the number of momenta at the numerator and with the number of propagators. That is why we would like to relate the integral (C.1) to simpler integrals such as,

These are the analogue of Passarino-Veltman integrals [START_REF] Passarino | One Loop Corrections for e+ e-Annihilation Into mu+ mu-in the Weinberg Model[END_REF] but in Euclidean space. To succeed we have to face two different issues. The first one is to make explicit the tensor dependence of the numerator so that the integrals become scalars. For doing so we are going to follow the idea presented in [START_REF] Passarino | One Loop Corrections for e+ e-Annihilation Into mu+ mu-in the Weinberg Model[END_REF]. The second task is to reduce the number of propagators in the denominator following the steps mentioned in [START_REF] Davydychev | Recursive algorithm of evaluating vertex type Feynman integrals[END_REF].

To explain how things work, let us start with the integrals B µ and B µν defined as

The integral B µ has only one external momentum p and one Lorentz index µ, therefore it can be written as

where B 1 is scalar integral that can be determined for example if we contract B µ with p µ ,

The next step is to relate the numerator q.p to the propagators, using, for instance, that

We can then cancel the inverse of the propagators which appear in the numerator with the propagators itself. This procedure simplify the integrals at the cost of increasing the number of integrals

For the other integral B µν , the tensor structure is a bit more complicated, we can describe the tensor structure as

We follow the same strategy to obtain B 21 and B 22 . For example, if we contract B µν with δ µν we have

where the numerator q 2 was written as (q 2 + m 2 1 ) -m 2 1 . The contraction with p µ p ν gives The study for

q µ ; q µ q ν ; q µ q ν q ρ ; q µ q ν q ρ q σ (q 2 + m 2 1 ) (q + p 1 ) 2 + m 2 2 (q + p 1 + p 2 ) 2 + m 2 3 (q + p 1 + p 2 + p 3 ) 2 + m 2 4 (C.6) can be found in [START_REF] Passarino | One Loop Corrections for e+ e-Annihilation Into mu+ mu-in the Weinberg Model[END_REF].

We also implement the algorithm for

The reduction to Passarino-Veltman integrals in simple cases was already implemented in the Mathematica package FeynCalc. However, the treatment done by FeynCalc was not enough to simplify our expressions so we implemented a Mathematica algorithm to do it. We implemented a reduction with a bigger number of propagators and Lorentz indices but just for the case where the value of the masses are 0 or m.

The final expression simplifies considerably because of the symmetries of Passarino-Veltman integral. To mention some of them, A(0) = 0 because of the use of dimensional regularization,

Once the numerator was simplified and the calculation has been reduced to the scalar integrals, we want to reduce the number of propagators to be able to express them in terms of A, B 0 and C 0 . We follow the strategy presented in [START_REF] Davydychev | A Simple formula for reducing Feynman diagrams to scalar integrals[END_REF]. Here, we change the notation a little in order to be clearer. We define the scalar function J(ν 1 , ν 2 , ν 3 ) as

for simplicity we omit the indices m 1 , m 2 , m 3 and p 1 , p 2 .

We would like to reduce the number of propagators. Let us start by observing that the integral of a divergence is zero.

= dJ(ν 1 , ν 2 , ν 3 ) -2ν 1 (q + p).(q + p)J(ν 1 + 1, ν 2 , ν 3 ) -2ν 2 (q + p 2 ).(q + p 1 )J(ν 1 , ν 2 + 1, ν 3 ) -2ν 3 q.(q + p 1 )J(ν 1 , ν 2 , ν 3 + 1).

(C.9)

We can now use the scalar products to reduce the number of propagators. To this end, we express the scalar products as follows:

(q + p 1 ).(q + p 1 ) = (q + p 1 ) 2 + m 2 1 -m 2 1 (q + p 2 ).(q + p 2 ) = (q + p 2 ) 2 + m 2 2 -m 2 2 q 2 = (q 2 + m 2 3 -m 2 3 (q + p 1 ).(q + p 2 ) = q 2 + q.p 1 + q.p 2 + p 1 .p 2 = 1 2 (q + p 1 ) 2 + m 2 1 + 1 2 (q + p 2 ) 2 + m 2 2 -1 2 (p 1 -p 2 ) 2 + m 2 1 + m 2 2 q.(q + p 1 ) = 1 2 (q + p 1 ) 2 + m 2 1 + 1 2 q 2 + m 2 3 -

Using the above relations, Eq. (C.9) now reads 0 = dJ(ν 1 , ν 2 , ν 3 ) -2ν 1 J(ν 1 , ν 2 , ν 3 ) -m 2 1 J(ν 1 + 1, ν 2 , ν 3 ) -ν 2 J(ν 1 , ν 2 , ν 3 ) + J(ν 1 -1, ν 2 + 1, ν 3 ) -(p 1 -p 2 ) 2 + m 2 1 + m 2 2 J(ν 1 , ν 2 + 1, ν 3 ) -ν 3 J(ν 1 , ν 2 , ν 3 ) + J(ν 1 -1, ν 2 , ν 3 + 1) -(p 2 1 + m 2 1 + m 2 3 )J(ν 1 , ν 2 , ν 3 + 1) .

Therefore, we can relate on function J with another function J with a smaller number of propagators, as 2ν 1 m 2 1 J(ν 1 + 1, ν 2 , ν 3 ) + ν 2 (p 1 -p 2 ) 2 + m 2 1 + m 2 2 J(ν 1 , ν 2 + 1, ν 3 ) +ν 3 p 2 1 + m 2 1 + m 2 3 J(ν 1 , ν 2 , ν 3 + 1) = (2ν 1 + ν 2 + ν 3 -d) J(ν 1 , ν 2 , ν 3 ) + ν 2 J(ν 1 -1, ν 2 + 1, ν 3 ) + ν 3 J(ν 1 -1, ν 2 , ν 3 + 1).

(C.10)

A similar procedure can be done starting with (q + p 2 ) µ in the numerator of Eq. (C.9). In this case we get the following relation ν 1 (p 1 -p 2 ) 2 + m 2 1 + m 2 2 J(ν 1 + 1, ν 2 , ν 3 ) + 2ν 2 m 2 2 J(ν 1 , ν 2 + 1, ν 3 ) +ν 3 p 2 2 + m 2 2 + m 2 3 J(ν 1 , ν 2 , ν 3 + 1) = (ν 1 + 2ν 2 + ν 3 -d) J(ν 1 , ν 2 , ν 3 ) + ν 1 J(ν 1 + 1, ν 2 -1, ν 3 ) + ν 3 J(ν 1 , ν 2 -1, ν 3 + 1).

(C.11)

Finally, we can do the same thing with q µ in the numerator of (C.9) ν 1 p 2 1 + m 2 1 + m 2 3 J(ν 1 + 1, ν 2 , ν 3 ) + ν 2 p 2 2 + m 2 2 + m 2 3 J(ν 1 , ν 2 + 1, ν 3 ) +2ν 3 m 2 3 J(ν 1 , ν 2 , ν 3 + 1) = (ν 1 + ν 2 + 2ν 3 -d) J(ν 1 , ν 2 , ν 3 ) + ν 1 J(ν 1 + 1, ν 2 , ν 3 -1) + ν 2 J(ν 1 , ν 2 + 1, ν 3 -1).

(C.12)

The three equations (C.10),(C.11) and (C.12) form the follow coupled system

(C.13) Accordingly, J(ν 1 + 1, ν 2 , ν 3 ), J(ν 1 , ν 2 + 1, ν 3 ) and J(ν 1 , ν 2 , ν 3 + 1) can be obtained as function of J with one less propagator.

With this iterative procedure we can take all the integrals to the form J(1, 1, 1), J(0, 1, 1), J(1, 0, 1), J(1, 1, 0), J(0, 0, 1) which can be easily related to Passarino-Veltman integrals.