One-loop unquenched three-gluon and ghost-gluon vertex in the Curci-Ferrari model

In this article we study the unquenched three-gluon and ghost-gluon vertex in all momentum range going from the ultraviolet to the infrared regime using the Curci-Ferrari model at one-loop in Landau gauge as an extension of the results presented in \cite{Pelaez:2013cpa}. Results are compared with recent lattice data for $SU(3)$ in the unquenched case. This calculation is a pure prediction of the model because it does not require fixing any parameter once two-point functions are fitted. An analysis of the influence of dynamical quarks in the position of the zero crossing of the three-gluon vertex is presented. Due to the recent improvements in infrared lattice data for the quenched three-gluon correlation function \cite{Aguilar:2021lke} we also redo the comparison between our one-loop results in this limit and the lattice obtaining a very good match.


I. INTRODUCTION
The infrared sector of QCD is usually called the Nonperturbative regime due to the fact that standard perturbation theory based on Faddeev-Popov Lagrangian presents a Landau pole in the infrared. This implies that perturbation theory cannot be applied together with this particular gauge-fixed Lagrangian to study the low energy region. For these reasons different semianalytical alternatives have been developed in order to approach this regime. For instance, approaches using non-perturbative functional techniques as treatments based on Schwinger-Dyson equations (SD) (see e.g. ), the functional renormalization group [25][26][27][28][29][30][31] or the variational Hamiltonian approach [32]. Other approaches focus on the fact that Faddeev-Popov procedure, generally used to fixed the gauge, is not completely justified in the infrared due to the existence of Gribov copies [33]. However, until now, it is not known how to build a gaugefixing procedure from first principles taking into account properly the problem of Gribov copies in the infrared. There are some interesting attempts to reach this gaugefixed Lagrangian based on Gribov-Zwanziger action and the refined Gribov-Zwanziger approach [34][35][36].
On top of these semi-analytical studies there are lattice simulations. Lattice simulations can deal with the problem of Gribov copies so they are a good tool to obtain information about the infrared behavior of Yang-Mills theory. Two important observations of lattice simulations are, first, that the gluon propagator reaches a finite nonzero value in the infrared, thus behaving as a massivelike propagator in this region. [37][38][39][40][41][42]. Second, that the relevant expansion parameter obtained through the ghost-gluon or the three-gluon vertex does not present a Landau-pole and in fact it does not become too large [15,38,43,44]. These points have motivated us to study the infrared regime using a gauge-fixed Lagrangian with a gluon mass term [45,46]. This Lagrangian is a particular case of Curci-Ferrari Lagrangians in Landau gauge (CF) [47]. On a different approach, the mentioned lattice results also motivates a screened massive perturbation theory where the mass term is added and subtracted changing the starting theory for expansion [48][49][50][51].
Even though we do not know how to justify the CF Lagrangian from first principles it is important to observe that it can reproduce a great variety of correlation functions using the first order in perturbation theory. It is important to mention that we do not attempt to reproduce all infrared quantities of QCD perturbatively. In particular the perturbative expansion for correlation functions involving quarks within CF near the chiral limit fails. Other approach using CF model was proposed in [52,53] in order to explore the chiral limit. See [54] for a detail summary of the already studied properties of the model. In particular, one-loop corrections within the CF model were computed for propagators, ghost-gluon vertex and the quenched three-gluon vertex [1,46,[55][56][57]. In addition to this, two-loop corrections were studied for propagators [58,59] and the ghost-gluon vertex with a vanishing gluon momentum [60] and compared with lattice data with great accuracy. It is important to mention that vertices are obtained as a pure prediction of the model, in the sense that the free parameters are fixed by minimizing the error between propagators and the corresponding lattice data and therefore there are no free parameters left when studying vertices.
In this article, we study the one-loop effects of dynamical quarks in the three-gluon vertex using the CF model. The unquenched results are compared with lattice data from [81]. Moreover, recent simulations of the quenched three-gluon vertex show a better handling of the infrared regime, yielding more precise data in this limit [2]. For this reason it is worth to extend the results presented in [1] for SU (2) to SU (3) gauge-group and compare it with the newest lattice data. For both cases, quenched and unquenched, the parameters used in the plots were chosen to minimize the error (understood as discrepancy with the lattice) of the propagators previously computed in [55,58]. In this sense, the results shown in this article are a pure prediction of the model that reproduces with great accuracy the lattice data. Due to the presence of massless ghosts, CF model also features a zero crossing as it is observed in [1,2,20,30,[65][66][67][68][69][70][71][72][73][74][75][76]. We also find that the inclusion of dynamical quarks shifts the zero crossing towards the infrared in a way consistent with what is observed in [68,69].
The article is organized as follows. In Sec. II we describe in more detail the Curci-Ferrari model in Landau gauge. We give some details on the one-loop calculations of the three-gluon vertex in Sec. III in terms of the Ball-Chiu components. In Sec. IV we present the renormalization conditions and the renormalization group equations. We present our results in Sec. V. and compare them with lattice data. At the end of the article we present the conclusions of the results.

II. CURCI-FERRARI MODEL WITH QUARKS
We start by introducing the Curci-Ferrari Lagrangian [47] in the presence of dynamical quarks in Euclidean space: where g is the coupling constant, and the flavor index i spans the N f quark flavors. The covariant derivative D µ acting on a ghost field in the adjoint representation of SU(N) reads (D µ c) a = ∂ µ c a + gf abc A b µ c c , while when applied to a quark in the fundamental representation it reads D µ ψ = ∂ µ ψ − gt a A a µ ψ. The latin indices correspond to the SU (N ) gauge-group, the t a are the generators in the fundamental representation and the f abc are the structure constants of the group. Finally, the field strength is given by F a µν = ∂ µ A a ν − ∂ ν A a µ + gf abc A b µ A c ν . The Feynman rules associated to this Lagrangian are the standard Feynman rules for Euclidean-QCD in Landau gauge except for the gluon's free propagator, which reads where we have introduced the transverse projector: It is important to mention that the gluon mass term added to the Faddeev-Popov Lagrangian breaks the BRST symmetry. However, it can be shown that (1) satisfies a modified (non-nilpotent) BRST symmetry that can be used to prove its renormalizability [82].
Probably the most interesting aspect of this model is the fact that, as it has been shown in various previous articles (see [45,46,57,58,88] for instance), the addition of a gluon mass term regularizes the theory in the infrared, allowing for a perturbative treatment of the theory in this region. More specifically, it is possible to find a renormalization scheme without an infrared Landau pole for particular choices of the initial condition of the renormalization-group flow. These features have made it possible to use this model to compute various two and three-point functions to 1-loop and 2-loop order, obtaining a very good match with lattice simulations [1,46,55,56,[58][59][60]. It is important to mention that this model has also been studied at finite temperature and chemical potential in [83][84][85]. In this work we extend the one-loop computation of the gluon three-point function obtained in [1] for SU (N ) Yang-Mills theory to unquenched QCD. In order to calculate the gluon's three-point function at one-loop order we need to compute the Feynman diagrams shown in Fig  5. As shown in [86], the color structure of the three-gluon vertex is simply the structure constant f abc of the SU (N ) group, so it can be factored out. Furthermore, we follow the usual convention [87] of factorizing the coupling constant, and thus we define: .
The tensor structure of Γ µνρ (p, k, r) can be easily deduced: it must depend on three Lorentz indices (one for each gluon) and on two independent momenta due to momentum conservation. As a consequence, we can have only two types of tensor structures: the ones made up of three momenta (p µ p ν p ρ , p µ p ν k ρ ,...) and the ones made up of one momentum and the Euclidean metric tensor (p µ δ νρ , k ν δ µρ ,...). It's not hard to convince oneself that there are eight possible terms of the first kind and six of the second, adding up to a total of 14 possible terms in the vertex's tensor structure. However, the vertex is symmetric under the exchange of any pair of external legs, and this ends up reducing the total number of possible independent coefficients in the vertex's tensor structure to six.
A cleaner way of exploiting these symmetries is following the decomposition proposed by Ball and Chiu in [89], where they parametrize the vertex using six scalar functions: The scalar functions have the following symmetry properties: A, C and F are symmetric under permutation of the first two arguments; B is antisymmetric under permutation of the first two arguments; H is completely symmetric and S is completely antisymmetric.
It is important to note that only some of these scalar functions are accessible through lattice simulations, since they have access to the vertex function only through the correlation function, i.e. the vertex contracted with the full external propagators. Since the propagators in Landau gauge are transverse, the longitudinal part of the vertex function is lost in the process when requiring the conservation of momentum at the vertex. In particular this means that the B and S functions are not accessible through lattice computations.
We decomposed every diagram contributing to the three gluon vertex into the Ball-Chiu tensorial structure. In this way we obtained the contribution of each diagram to each of the scalar functions A, B, C, F, H and S. To perform our computations we expressed the integrals in Feynman diagrams of Fig. 5 in terms of only three Master Integrals defined following the convention on [90] as: whereC = 16π 2μ 2 (2π) d , and the regularization scaleμ is related to the renormalization scale µ by µ 2 = 4πe −γμ2 . The A and B-Master Integrals can be solved analytically in d = 4 − 2 in terms of the external momentum and the masses, but the C-Master Integral must be treated numerically except for particular kinemetics. We chose the FIRE5 algorithm [91] to perform the Master Integral re-duction, thus obtaining analytic expressions for each of the scalar functions in terms of the three Master Integrals for arbitrary momentum configurations. The expressions are complicated and not very enlightening, however, the explicit expressions appear in the suplemental material of [1] for the quenched case while the quark contribution is written in [87]. In the case of one vanishing external momentum the computation becomes much simpler and the result for the quenched vertex function in this configuration is given in [1], while the unquenched case is presented in the Appendix A.

B. Checks
Various checks for the Yang-Mills part of the result for the three-gluon vertex function had already been made in [1]. We only need to check the quark triangle diagram to test our unquenched results. To do this we compared our results to those of [87], verifying that they yield the same expressions when properly transformed to the Euclidean space. This was expected, as the quark triangle diagram is independent from the mass of the gluons and therefore its contribution in the Curci-Ferrari model is the same as in standard QCD.

IV. RENORMALIZATION AND RENORMALIZATION GROUP
In this section we introduce the renormalization scheme that we used in this work and we explain how we implemented renormalization-group ideas to improve our perturbative calculation.

A. Renormalization
To take care of the divergences appearing in the 1loop quantities we took the usual approach of redefining the fields, masses and coupling constants through renormalization factors that absorb the infinities. The renormalized quantities are defined in terms of the bare ones (denoted with a "0" subscript) as follows: From now on, all expressions will refer to renormalized quantities unless explicitly stated otherwise. The renormalized propagators and the three-gluon 1PI correlation function are thus defined as: B. Infrared Safe renormalization scheme To fix the renormalization factors we used the Infrared Safe (IS) renormalization scheme defined in [46]. It is based on a non-renormalization theorem for the gluon mass [92][93][94], and is defined by where Γ ⊥ (p) is the transversal part of Γ (p) and J(p) is the ghost dressing function. To fix the renormalization of the coupling constant we used the Taylor scheme, in which the coupling constant is defined as the ghost-gluon vertex with a vanishing ghost momentum. Requiring that the renormalized vertex is finite leads to a relation among the renormalization factors Z A , Z c and Z g : The divergent part of the renormalization factors for the quenched case were presented in [46]. Here we show the extension to the unquenched Curci-Ferrari model already computed in [55,95]. In d = 4 − 2 they read: Finally, the quantity we are interested in is actually Γ µνρ as defined earlier. Since in it's definition we factorized a factor of g, the relation between the renormalized and bare quantities is the following: where in the last equality we used equation 9.

C. Renormalization Group
After the renormalization procedure we obtain a finite expression for the three-gluon vertex, but it comes with the usual loop corrections of the form log( p µ ). To handle this situation we implemented the renormalization-group flow to take care of the large logarithms coming from loop corrections. First we define the β functions and anomalous dimensions of the fields: where χ can take the role of the coupling constant, the gluon mass or the quark mass and φ represents the different fields A, c, ψ.
The renormalization group equation for the vertex function with n A gluon legs and n c ghost legs reads: This equation allows us to obtain a relation for the vertex function renormalized at scale µ 0 and the same vertex function renormalized at a different scale µ: where g(µ), m 2 (µ) and M (µ) are obtained by integration of the β functions with initial conditions given at some scale µ 0 and z A and z c are obtained repectively from: We can then use the non-renormalization theorems of Eq.(8) and Eq.(9) to relate the anomalous dimensions and the β functions. It is simple to check that one obtains the following relations: Finally we use these relations to integrate Eq.(15), obtaining analytical expressions for z A and z c in terms of the running gluon mass and coupling constant, being this feature another of the advantages of the Infrared Safe scheme: We are able now to express the three-gluon vertex renormalized at scale µ 0 in terms of the same quantity using a running scale µ = p, thus avoiding the large logarithms problem. Taking into account the factor of g on the definition of Γ µνρ (p, r) this reads:

V. RESULTS
We now present our results for the different scalar functions associated to the three-gluon vertex introduced in the previous section. All our results correspond to four dimensions and the SU (3) gauge group, and we evaluate the scalar functions in different momentum configurations in order to compare them with the available lattice data.

A. Fixing Parameters
The only fitting parameters we need to adjust to compare our results with the lattice are the overall normalization constant of the gluon three-point function and the initial conditions of the renormalization-group flow, i.e. the values of the mass of the gluon, the mass of the quark and the coupling constant at some renormalization scale µ 0 .
The initial conditions for the renormalization-group are best obtained by looking for the set of parameters (m 0 , M 0 , g 0 ) that produce the best fit between the gluon and ghost propagators computed using the Curci-Ferrari approach and the lattice data, since the lattice results are much more precise for propagators than for threepoint functions. This task was done in [1,58] for different gauge groups and renormalization schemes in the quenched case, and in [55] including dynamical quarks. For the SU (3) group and the IS scheme the initial conditions for the R-G flow at µ 0 = 1 GeV obtained are the ones listed in table I .
In this work we use these values to compute the one loop three-gluon vertex, which means that up to the overall normalization constant our results are a pure prediction of the model.
For the quenched case, we compared our results with the lattice data from [2]. Following their definitions, in the symmetric configuration we work with the scalar functionsΓ sym 1 andΓ sym 2 : (20) where µ ν ρ (p 1 , p 2 , p 3 ) defined as the perturbative treelevel tensor of the three-gluon vertex, and λ µνρ On the other hand the asymmetric configuration of the vertex is parametrized by a single scalar functionΓ asym 3 defined by with λ µνρ 3 (p, −p, 0) = 2p ρ P ⊥µν (p).
We compared our unquenched results with the lattice data from [81]. They work in the orthogonal configuration, and define the usual scalar function G 1 , which consists on contracting the external legs of the vertex with transverse propagators and the tree-level momentum structure of the 3-gluon vertex, normalized to the same expression at tree-level. This reads: The results of the model are shown below using the different scalar functions defined in this section including the renormalization group effects.

C. SU (3) Yang Mills results
We first present our results for SU (3) Yang Mills theory and we compare them with lattice results from [2]. As stated before, we integrate the beta functions with initial conditions at µ 0 = 1 GeV using the initial conditions listed in table I.
In Fig. 2 we show the results for the scalar func-tionΓ asym 3 in the asymmetric configuration (one vanishing momentum), and in Fig. 3 we do the same for the functionsΓ sym 1 andΓ sym 2 the symmetric configuration (all momenta equal).
In all cases the agreement is very good, specially considering that the initial conditions for the renormalization-group flow were not fitted for the three-point function but for the propagators. It is also noticeable that in all cases the different scalar functions become negative at low energies, a qualitative feature that was observed in many lattice simulations as well as in different analysis [65,71,[73][74][75][76]. While the scalar functions associated to the tree-level tensor diverges logarithmically, theΓ sym 2 goes to a constant value in the infrared as stated in [2]. The simplicity of one-loop CF model allows to write the infrared behaviour ofΓ sym 2 analytically: which is indeed finite in the infrared, and where Cl 2 is the Clausen function satisfying Cl 2 π 3 = 1.0149417. It is worth mentioning that this behaviour is not modified by the effects of the renormalization group.
These results also show the divergent behavior ofΓ 1 , which can be easily understood due to the presence of (bottom) scalar functions as a function of momentum for all momenta equal (symmetric configuration). The points are lattice data from [2]. The plain red line corresponds to our 1-loop computation. massless ghosts as stated in [1].

D. Unquenched QCD results
If we want to include the influence of dynamical quarks to the previous computation we must add the quark triangle diagram to the vertex and use the running of the coupling obtained in the unquenched analysis [55]. The contribution of that diagram can be computed with no difficulty in arbitrary dimension and for arbitrary number of quarks. Our explicit expressions match the ones presented in [87] when continuing them to Euclidean space. In order to be more specific we show as an example the explicit bare contribution of the quark-loop diagram to the factor G 1 in d = 4 − 2 dimensions in Appendix A.
The total result for G 1 is shown in Fig. 4 where it is compared with lattice data from [81]. The data available corresponds to the G 1 scalar function in the case of two mass-degenerate quark flavors (N f = 2) in the orthogonal configuration (two momenta orthogonal to each other and of equal magnitude). In this case, the agreement is still very good in the infrared but worsens in the UV. More precisely, the model and the lattice results start separating at a scale of about 2.5 GeV. This scale is of the order of magnitude of the inverse of the lattice spacing used in most lattice simulations, and therefore lattice results beyond this scale are subject to hypercubic lattice spacing artefacts. Taking this fact into account and also considering that perturbation theory must work at one-loop in the UV, we suspected that the decrease in the values of G 1 after the inverse lattice spacing scale must be caused by finite lattice artefacts such lattice effects in the UV.
To confirm this statement, we computed analytically the high-energy limit of G 1 , finding that it behaves in the UV as Ln which is compatible with our results. The idea behind this computation is that since the UV limit of G 1 (p, p) is equal to 1, the high-energy behavior of G 1 (µ 0 , p) must be given by z to the three-gluon vertex arising from the different oneloop Feynman diagrams, which are shown in the asymmetric configuration in Fig. 5. Our results match the behaviour observed in [18] for the Yang-Mills contributions, and indeed show that the divergence of the vertex in the infrared is caused by the ghost triangle's logarithmic divergence.

Zero crossing and the number of flavours.
In this section we study the influence of dynamical quarks in the position of the zero crossing of the threegluon vertex. In the one-loop CF-model the influence of quarks in the renormalized vertex can be isolated as the term proportional to N f . At one loop, the quarks contribution to the renormalized vertex is proportional to: represent the finite part of the coefficient proportional to g 2 N f T f in the Z A -and Z g -renormalization factor and in G 1 respectively.
We studied the sign of the factor F in order to analyze in which direction the zero crossing is shifted. As F depends on the finite part of the renormalization factor, it is expected that its sign depends on the chosen renormalization scheme. For the IS-scheme we observe that even though the contribution of the quark triangle diagram is positive (see Fig.5) the renormalization factors together with the renormalization group flow makes F negative. This means that the whole contribution arising from considering dynamical quarks is negative when compared to the quenched quantity using the same flow. However, as flows should be different in each situation, it is better to study the influence of dynamical quarks in the infrared region using the same set of initial conditions of the renormalization group flow at an ultraviolet scale (for instance µ = 3 GeV). In this context, we can see in Fig. 6 that dynamical quarks shift the zero crossing to the infrared as it is also observed by [68,69].

Unquenching the ghost-gluon vertex.
It is also interesting to observe the influence of dynamical quarks in the ghost-gluon vertex. Eventhough quarks do not contribute directly in its one-loop diagrams the inclusion of dynamical quarks affects its renormalization  Table I in the orthogonal configuration (top), symmetric configuration (middle) and gluon vanishing momentum (bottom) group flow. The function G ccA (p, k, r) defined through the vertex as: where is shown in Fig. 7 for different kinematical configurations using the parameters of Table I which are the parameters that give a better fit to one-loop propagators. It is important to mention that G ccA (p, k, r) is renormalized by the combination Z g √ Z A Z c which is set to one accordingly to Taylor scheme (9). Therefore the value of ghost-gluon function only depends on the values of g and m on each momentum scale. In particular, initial conditions of the renormalization flow are initialized at 1GeV, therefore the vertex function at that scale is purely determined by the values from Table I. As the value of g 0 is larger in the unquenched case, the unquenched vertex function will be above the quenched one as is observed in Fig. 7.
To study the influence of dynamical quarks on the ghost-gluon vertex is therefore convenient to raise the flavor number but using the same initial conditions of the flow. In Fig. 8 we compare the ghost-gluon vertex when raising N f using the same set of parameters at 3GeV. It can be seen that the unquenching effects reduces the vertex contribution in the infrared. Another observation is that the unquenched ghost-gluon vertex is almost insensitive to the value of the quark mass as it is seen in Fig. 9.

VI. CONCLUSIONS
With the aim of studying the infrared properties of the gluon and ghost-gluon three-point correlation functions we presented a one-loop calculation using Curci-Ferrari model in Landau gauge for arbitrary kinematical configurations. The results are an extension of a previous work [1] to the unquenched case. In particular, we compared the results for the vertex with the available lattice data including dynamical quarks corresponding to the kinematical configuration with a vanishing gluon momentum and two degenerate flavors. A study of the position of the zero crossing of the vertex was also done, observing that the position of the zero crossing is shifted to the infrared due to the presence of dynamical quarks when compared to the quenched case.
We also studied the quenched case because some infrared properties observed by the model in the previous work were not clear in lattice simulations at that time. However, the infrared lattice study of the three-gluon vertex has improved in the last years and now error bars are good enough to understand its infrared behavior as discussed in [2]. In particular, lattice simulations show a change of sign in the deep infrared that is easily understood by the Curci-Ferrari model. As it has been discussed in [1,65] it can be explained as a consequence of the diagram with a loop of massless ghosts. In the quenched case we compare the results with lattice simulations for the completely symmetric and the antisymmetric configuration. The results show an excellent match to lattice results, specially considering that the free parameters of the model were already fixed by fitting the propagators.