Complex poles and spectral functions of Landau gauge QCD and QCD-like theories

In view of the expectation that the existence of complex poles is a signal of confinement, we investigate the analytic structure of gluon, quark, and ghost propagators in the Landau gauge QCD and QCD-like theories by employing an effective model of Yang-Mills theory with a gluon mass term, which we call the massive Yang-Mills model. In this model, we particularly investigate the number of complex poles in the parameter space of the model consisting of gauge coupling constant, gluon mass, and quark mass for the gauge group $SU(3)$ and various numbers of quark flavors $N_F$ within the asymptotic free region. This investigation extends the previous result obtained for the pure Yang-Mills theory with no flavor of quarks $N_F=0$ that the gluon propagator has a pair of complex conjugate poles and the negative spectral function while the ghost propagator has no complex pole. The gluon and quark propagators at the best-fit parameters for $N_F=2$ QCD have one pair of complex conjugate poles as in the zero flavor case. By increasing quark flavors, we find a new region in which the gluon propagator has two pairs of complex conjugate poles for light quarks with the intermediate number of flavors $4 \lesssim N_f<10$. However, the gluon propagator has no complex poles if very light quarks have many flavors $N_f \geq 10$ or both of the gauge coupling and quark mass are small. In the other regions, the gluon propagator has one pair of complex conjugate poles. Moreover, as a general feature, we argue that the gluon spectral function of this model with nonzero quark mass is negative in the infrared limit. In sharp contrast to gluons, the quark and ghost propagators are insensitive to the number of quark flavors within the current approximations adopted in this paper. These results suggest that details of the confinement mechanism may depend on the number of quark flavors and quark mass.

In view of the expectation that the existence of complex poles is a signal of confinement, we investigate the analytic structure of gluon, quark, and ghost propagators in the Landau gauge QCD and QCD-like theories by employing an effective model of Yang-Mills theory with a gluon mass term, which we call the massive Yang-Mills model. In this model, we particularly investigate the number of complex poles in the parameter space of the model consisting of gauge coupling constant, gluon mass, and quark mass for the gauge group SU (3) and various numbers of quark flavors NF within the asymptotic free region. This investigation extends the previous result obtained for the pure Yang-Mills theory with no flavor of quarks NF = 0 that the gluon propagator has a pair of complex conjugate poles and the negative spectral function while the ghost propagator has no complex pole. The gluon and quark propagators at the best-fit parameters for NF = 2 QCD have one pair of complex conjugate poles as in the zero flavor case. By increasing quark flavors, we find a new region in which the gluon propagator has two pairs of complex conjugate poles for light quarks with the intermediate number of flavors 4 N f < 10. However, the gluon propagator has no complex poles if very light quarks have many flavors N f ≥ 10 or both of the gauge coupling and quark mass are small. In the other regions, the gluon propagator has one pair of complex conjugate poles. Moreover, as a general feature, we argue that the gluon spectral function of this model with nonzero quark mass is negative in the infrared limit. In sharp contrast to gluons, the quark and ghost propagators are insensitive to the number of quark flavors within the current approximations adopted in this paper. These results suggest that details of the confinement mechanism may depend on the number of quark flavors and quark mass.

I. INTRODUCTION
Color confinement, absence of color degrees of freedom from the physical spectrum, is one of the most fundamental and significant features of strong interactions. It is a long-standing and challenging problem in particle and nuclear physics to explain color confinement in the framework of quantum field theory (QFT). Understanding analytic structures of the correlation functions will be of crucial importance to this end because a QFT describing physical particles can be reformulated in terms of correlation functions [1], and there are some proposals of confinement mechanisms whose criteria are expressed by them, e.g., [2]. In particular, the analytic structures of propagators encode kinematic information as the Källén-Lehmann spectral representation [3], which will be useful towards understanding confinement.
In the past decades, numerous studies of both the lattice and continuum approaches have focused on the gluon, quark, and ghost propagators in the Landau gauge of the Yang-Mills theory and quantum chromodynamics (QCD). In the Yang-Mills theory, or the quenched limit of QCD, the so-called decoupling and scaling solutions of the gluon and ghost propagators are observed based on the continuum approaches [4]. The recent numerical lattice results support the decoupling solution [5].
The decoupling solution has an impressing feature that the running gauge coupling stays finite and nonzero for all non-vanishing momenta and eventually goes to zero in the limit of vanishing momentum, which cannot be predicted from the standard perturbation theory that is plagued by the Landau pole of the diverging running coupling.
A low-energy effective model of the Yang-Mills theory is proposed following the decoupling behavior of the gluon propagator, which provides the gluon and ghost propagators that show a striking agreement with the numerical lattice results by including quantum corrections just in the one-loop level [6,7]. This effective model is given by the mass-deformed Faddeev-Popov Langrangian of the Yang-Mills theory in the Landau gauge, or the Landau limit of the Curci-Ferarri model [8], which can be shown to be renormalizable due to the modified BRST symmetry, and we call it the massive Yang-Mills model for short. This effective mass term could stem from the dimension-two gluon condensate [9][10][11][12][13] or could be taken as a (minimal) consequence of avoiding the Gribov ambiguity [14,15]. Moreover, it has been shown that the massive Yang-Mills model has "infrared safe" renormalization group (RG) flows on which the running gauge coupling remains finite for all scales [7,15,16]. The three-point functions [17] and two-point correlation functions at finite temperature [18] in this model were compared to the numerical lattice results, showing good agreements. Moreover, the two-loop corrections improve the accordance for the gluon and ghost propagator [19]. These works indicate the validity of the massive Yang-arXiv:2001.05987v1 [hep-th] 16 Jan 2020 Mills model as an effective model of the Yang-Mills theory.
The gluon and ghost propagators for unquenched lattice QCD with the number of quark flavors N F = 2, 2 + 1, 2 + 1 + 1 have been studied, for instance [20], and exhibit the decoupling feature as well. The massive Yang-Mills model with dynamical quarks can reproduce the numerical lattice gluon and ghost propagators for QCD as well [21]. However, it is argued that higher-loop corrections are important for the quark sector in this model [22,23]. Despite this shortcoming, it appears that the effective gluon mass term captures some non-perturbative aspects of QCD. What is more, QCD phases have been extensively studied in a similar model with the effective gluon mass term of the Landau-DeWitt gauge [24].
Apart from the realistic QCD, it is also interesting to study gauge theories with many flavors of quarks. For many quark flavors, the infrared conformality is predicted [25,26] and well-studied in line with the walking technicolor of the beyond standard model [27], for example, [28][29][30]. Some argue that chiral symmetry restores while color degrees of freedom are "unconfined" in the conformal window. For a better understanding of the confinement mechanism, observing N F dependence will be thus extremely valuable.
All works on the correlation functions described above were implemented in the Euclidean space. Considerable efforts have been devoted to reconstructing the spectral functions from the Euclidean data based on the Källén-Lehmann spectral representation, e.g., [31,32]. On the other hand, several models of Yang-Mills theory [14,[33][34][35][36][37][38][39], including the (pure) massive Yang-Mills model [40,41], and a way of the reconstruction from the Euclidean data [42] predict complex poles in the gluon propagator that invalidate the Källén-Lehmann spectral representation. The existence of complex poles of the propagators of the confined particles is a controversial issue, e.g., [43]. The complex singularities invalidate the standard reconstruction from a Euclidean field theory to a relativistic QFT [1] and might correspond to unphysical degrees of freedom in an indefinite metric state space [44]. Therefore, complex poles are expected to be closely connected to the confinement mechanism.
In this paper, we investigate the analytic structure of the QCD propagators for various N F based on the massive Yang-Mills model, mainly focusing on complex poles by utilizing the general relationship between the number of complex poles and the propagator on timelike momenta from the argument principle [40].
This paper is organized as follows. In Sec. II, we present machinery to count the number of complex poles as an application of [40]. In Sec. III, we review the calculation of the massive Yang-Mills model with quarks of [21], consider the infrared safe trajectories of this model, and argue the infrared negativity of the gluon spectral function. Then, in Sec. IV, we analyze the analytic structures of the gluon, quark, and ghost propagators at the best-fit parameters of [21], and investigate the number of complex poles for various parameters and the number of quarks N F . Finally, Sec. IV is devoted to conclusions and discussions. In Appendix A, we provide a generalization of the proposition of Sec. II for various infrared behaviors.

II. COMPLEX POLES IN PROPAGATORS
To elucidate the starting point, we review a generalization of the spectral representation for a propagator so as to allow complex poles. We then develop a method for counting the number of complex poles from the data on timelike momenta, e.g. the spectral function, as a straightforward application of the general relation [40].

A. A generalization of the spectral representation
We introduce some definitions and underlying assumptions on propagators adopted in this article. Given a propagator defined in the Euclidean space, we analytically continue the propagator to the whole complex k 2 plane from the Euclidean momenta. In the complex k 2 plane, we call points on the negative real axis Euclidean momenta and points on the positive real axis timelike momenta. We will study the gluon, quark, and ghost propagators in the Landau gauge. We assume each propagator of the gluon, quark, or ghost has the following generalized spectral representation allowing the presence of complex poles: 1 where ρ(σ 2 ) is its spectral function, z stand for the positions of complex poles, and Z are the residues associated with the complex poles. This representation can be derived under the following conditions using the Cauchy integral formula for the closed contourC presented in Fig. 1 [40,45]: (i) D(z) is holomorphic except for singularities on the positive real axis and a finite number of simple poles.
(iii) D(z) is real on the negative real axis.
If we replace the condition (i) with the more strict condition FIG. 1: ContourC on the complex k 2 plane avoiding the singularities on the positive real axis and the complex poles. The contourC consists of the large circle C1, the path wrapping around timelike singularities C2, and the small contours gamma that clockwise surround the complex poles at z . The propagator D(k 2 ) is holomorphic in the region bounded by the contourC = C ∪ {γ } n =1 , where we denote the closed contour C1 ∪ C2 by C.
(i') D(z) is holomorphic except for singularities on the positive real axis, namely has no complex poles, the three conditions (i'), (ii), and (iii) lead to the Källén-Lehmann form [3], which propagators for unconfined particles are supposed to obey. In this sense, eq. (1) gives a generalization of the Källén-Lehmann spectral representation.

B. Counting complex poles
We present a procedure to count complex poles from the propagator on the timelike momenta based on [40].
We apply the argument principle to a propagator on the contour C presented in Fig. 1. Then the winding number N W (C) of the phase of the propagator D(k 2 ) along the contour C is equal to the difference between the number of zeros N Z and the number of poles N P in the region bounded by C, The winding number N W (C) can be calculated from the propagator on timelike momenta D(−k 2 + i ) and infrared (IR) and ultraviolet (UV) asymptotic forms.
In what follows, we assume the following asymptotic form for the propagator.
Let us comment on these assumptions. The first assumption is satisfied by the quark and gluon propagators in the Landau gauge, which follows from the RG analysis for asymptotic free theories. The RG argument of Oehme and Zimmermann [46] provides the following UV asymptotic form for the propagators, as |k 2 | → ∞, where Z U V is a positive constant and γ = γ 0 /β 0 is the ratio between the first coefficients γ 0 , β 0 of the anomalous dimension and the beta function, respectively. For the gluon propagator of Yang-Mills theories with N F quarks in the Landau gauge, γ is computed as follows.
where C 2 (G) and C(r) = N F /2 are the Casimir invariants of the adjoint and fundamental representations of the gauge group G, respectively. We restrict ourselves to asymptotically free theories: β 0 < 0, or N F < 11 2 C 2 (G), which is essential to derive the above UV asymptotic expression for the propagator (3).
Additionally, the one-loop corrections give no contribution to the quark anomalous dimension γ ψ in the Landau gauge because γ ψ = O(g 4 ). Therefore, the quark propagator behaves in the UV limit as the free one because of the asymptotic freedom and the vanishing of the first coefficient γ 0 ψ of the order g 2 of the quark anomalous dimension.
The second assumption (ii) indicates the massive-like behavior for the propagator, which corresponds to the decoupling solution for the Euclidean gluon propagator. For the general cases D(k 2 ) → Z IR (−k 2 ) α as |k 2 | → 0 with a real exponent α, e.g., the scaling solution for the pure Yang-Mills gluon and a massless propagator, there is an additional contribution to the expression of N W (C) described below (9). See Appendix A for the details of the additional term. From here on, we simply assume (ii) and will verify this assumption when we compute N W (C) for each propagator employed.
Let us add notes on the relationships between the number N W (C) and the spectral function [40]. With the assumptions (i) and (ii), the positive spectral function implies N W (C) = 0 and the negative one implies N W (C) = −2. Since the winding number is a topological invariant, N W (C) is invariant under continuous deformations. For example, if the spectral function is "quasi-negative", i.e. the spectral function is negative ρ(k 2 0 ) < 0 at all real and positive zeros k 2 0 of Re D(k 2 ) i.e., Re D(k 2 0 ) = 0 (k 2 0 > 0), then the propagator has N W (C) = −2. Actually, this is the case of the massive Yang-Mills model with N F = 2 quarks at the realistic parameters analyzed below.
In order to calculate the winding number N W (C) in a numerical way, we divide the interval [δ 2 , Λ 2 ] on the positive real axis into (N + 1) segments x 0 , x 1 , · · · , x N +1 such that the following condition on is sufficiently dense so that D(k 2 = x + i ) changes its phase at most halfwinding (±π) between x n + i and x n+1 + i , i.e., for n = 0, 1, · · · , N , where we denote sufficiently small x 0 = δ 2 > 0 and sufficiently large x N +1 = +Λ 2 , on which we will take the limits δ 2 → +0 and Λ 2 → +∞.
Let us now calculate N W (C) from the data {D(x n + i )} N n=1 under the above assumptions, by evaluating N W (C 1 ) and N W (C 2 ) separately, where C 1 stands for the large circle and C 2 for the path around the positive real axis depicted in Fig. 1.
The first assumption (i) yields For N W (C 2 ), from the Schwarz reflection principle, [D(z)] * = D(z * ), we have, Notice that we have used the second assumption (ii) to eliminate the contribution from the small circle around the origin. The third assumption (iii) transforms the integral (7) into a discrete sum as where Arg is the principal value of the argument (−π < Arg z < π).
To sum up, we have the expression for N W (C), Though we have to know the number of complex zeros for counting the exact number of complex poles N P from the winding number N W (C), note that it suffices to verify N W (C) < 0 in order to show the existence of complex poles because N P = N Z − N W (C) ≥ −N W (C). Moreover, since N W (C) is invariant under continuous deformations, we can expect that N W (C) should be robust under some approximations. We will consider complex poles of the QCD propagators using an effective model and the relation (9) based on this expectation.

III. MASSIVE YANG-MILLS MODEL AS AN EFFECTIVE MODEL
The massive Yang-Mills model to be defined shortly is an effective model of the Yang-Mills theory in the Landau gauge, which captures some non-perturbative aspects by introducing a phenomenological mass term for the gluon. The mass term may be generated by the effect of the dimension-two gluon condensate [9][10][11][12][13] or by avoiding the Gribov ambiguity [15]. For the latter effect, intuitively, some effect suppressing d D x A A µ A A µ should be taken into account due to the Gribov copies, and the most IR relevant term of this effect will be the mass term. The massive Yang-Mills model reproduces the numerical lattice results of the gluon and ghost propagators well and has a renormalization condition such that there are RG trajectories whose running gauge coupling remains finite in all scales. Moreover, the massive Yang-Mills model gives a correct UV asymptotic behavior [46] by RG improvement, as we will see in Sec. III C. Therefore, although we expect the massive Yang-Mills model is a low-energy effective model, this model can describe Yang-Mills theory in all scales to some extent. One might worry about the absence of the nilpotent BRST symmetry. Although the massive Yang-Mills model suffers from the problem of physical unitarity [47] as a consistent QFT, we insist that this model will suffice for our purpose of investigating analytic structures, as it gives the well-approximating propagators with a sensible and straightforward mass-deformation.
Moreover, by adding effective quark mass term, this model has good accordance with the numerical lattice results for the unquenched gluon and ghost propagators and can also reproduce the quark mass function qualitatively [21]. In this section, we review the results on the one-loop computations of the massive Yang-Mills model with quarks, prepare expressions that we will use in the investigation of the analytic structures of the propagators in the next section, and study asymptotic behaviors of the propagators of this model using RG.
In the Euclidean space, the Lagrangian of the model is given by We introduce the renormalization factors ( for the gluon, ghost, anti-ghost, and the quark fields (A µ , C ,C , ψ i ), the gauge coupling constant g, the gluon and quark mass parameters M 2 , m q,i respectively: Throughout this article, for simplicity, we employ this model with degenerate quark masses, m q := m q,i , and therefore Z ψ := Z (i) ψ .

A. Strict one-loop calculations
We review the strict one-loop results for quark and gluon propagators here and the RG functions in the next subsection [21].
For the gluon, ghost, and quark, we introduce the twopoint vertex functions Γ (2) A , Γ (2) gh , and Γ (2) ψ , the transverse gluon propagator D T , the ghost propagator ∆ gh , the quark propagator S, dimensionless gluon and ghost vacuum polarizationsΠ andΠ gh , and the scalar and vector part of the quark two-point vertex function Γ where k E is the Euclidean momentum, and The bare vacuum polarizations computed by the dimensional regularization read [7,21], for gluons, for ghosts, where ε := 2 − D/2, γ is the Euler-Mascheroni constant, and, witht For the quark sector [21], where and for G = SU (N ) and fundamental quarks.
Henceforth, we adopt the "infrared safe renormalization scheme" by Tissier and Wschebor [7], which respects the non-renormalization theorem Z A Z C Z M 2 = 1 [49], in the one-loop level, combined with the Taylor scheme [48] In this renormalization scheme, the running coupling of some RG flows turns out to be always finite, which implies that the perturbation theory will be valid to some extent. By imposing the above renormalization condition, we have the renormalized two-point vertex functions, Note that the gluon propagator exhibits the decoupling feature and satisfies the condition (ii) of the previous section, Indeed, where we have used h q (t → ∞) = O(1), h(s) = − 111 2s + O(ln s), f (0) = 5/2, and the fact that f (s) is a monotonically increasing function.
Finally, note that the one-loop expression of gluon propagator will have disagreement at momenta far from the renormalization scale µ. Indeed, in the UV limit, the strict one-loop expression has the wrong asymptotic form: while the RG analysis yields where we have analytically continued the gluon propagator from Euclidean momentum k 2 = −k 2 E to complex k 2 , γ 0 and γ = γ 0 /β 0 are given in (4), and Z U V > 0.
However, for γ = γ 0 /β 0 > 0, the phase of the gluon propagator on the ultraviolet momentum is qualitatively correct despite the wrong exponent of the logarithm. Furthermore, both of the propagators have negative spectral functions ρ(σ 2 ) := Im D T (σ 2 + i ) < 0 in the UV limit. Therefore, we expect the one-loop expression will provide a good approximation of the phase of the gluon propagator for γ > 0 based on the robustness of the winding number. In Sec. IV, we indeed confirm that the RG improved and strict one-loop gluon propagator yield the same N W (C) for qualitatively the same parameter region.

B. Renormalization group functions
Here we present the renormalization group functions for (g, M 2 , m q , Z A , Z C , Z ψ ) [7,21]. For later convenience, we put The beta functions β α and anomalous dimensions γ Φ are defined for the masses and coupling α = g, M 2 , m q , λ and for the fields Φ = A , C , ψ i renormalized at a scale µ as where the bare masses and the coupling are fixed in taking the derivative. From the nonrenormalization theorems, Z A Z C Z M 2 = 1 and Z g √ Z A Z C = 1, β λ , or equivalently β g , and β M 2 are expressed by γ A and γ C : The other RG functions, γ C , γ A , γ ψ and β mq are computed as follows within the one-loop approximation. By differentiating the counterterms, γ Φ = µ d dµ δ Φ , we have for ghosts, for gluons, where and, for quarks, Since Z mq = 1 + δ mq − δ ψ in the one-loop level, we obtain the beta function for the quark mass m q , where Q and K 2 are given in (26) Fig. 2. As the pure Yang-Mills case, the diagram has the "infrared safe" trajectories, on which the running gauge coupling λ is always finite in the all scales. One can confirm that the qualitative feature is the same for N F ≤ 16.
Beforehand, there is a caveat on this discussion. In the RG analysis of Oehme and Zimmermann [46], the parameter of the theory is only the gauge coupling g, which guarantees that the higher-loop effects are suppressed in the ultraviolet region precisely. Then, the one-loop RG gives a strong argument on asymptotic behaviors and enables us to establish the UV negativity of the gluon spectral function for N F < 10. However, The arrows indicate infrared directions µ → 0. The blue trajectories are infrared safe.
the massive Yang-Mills model has the mass parameter, which can potentially invalidate the perturbation theory in λ even though λ → 0 asymptotically. Here, we assume that the perturbation theory "works well". In particular, we assume that the spectral function is dominated by the one-loop contribution in both UV and IR limits.
In the pure Yang-Mills case, a similar discussion on the RG functions and Euclidean propagators can be found in e.g., [16]. In what follows, we study asymptotic behaviors of the RG functions, Euclidean propagators, and spectral functions in the massive Yang-Mills model described in the previous subsections.

UV limit
We consider the UV limit and confirm that the asymptotic behaviors of the massive Yang-Mills model is consistent with those of the Faddeev-Popov Lagrangian. Taking the limit u, t → 0, we have which reproduce the standard one-loop beta function: recovering the UV asymptotic behavior of the gauge coupling λ ∼ −β 0,λ / ln(µ/Λ) with some constant scale Λ. On the other hand, the beta functions for the masses read, for gluons, and for quarks, Notice that at G = SU (3), the gluon mass is suppressed logarithmically for β 0,M 2 < 0, or N F < 35 8 C 2 (G) = 105/8 ≈ 13.1 and enhanced logarithmically for 14 ≤ N F ≤ 16, while the correction to the quark mass always suppresses the quark mass in the logarithmic way.
One can see the UV negativity for the gluon spectral function and UV positivity for the ghost spectral function from (55) combined with the analyticity or RG improvement with respect to |k 2 | in the complex k 2 plane [46].
Thus, it turns out that the massive Yang-Mills model can describe the correct UV behavior by RG improvement from the above observations, although we have regarded the massive Yang-Mills model as a low-energy effective model.

IR limit
Next, we consider the IR limit of the infrared safe trajectories: λ → 0, u → ∞, t → ∞. Taking the limit u, t → ∞, we have as those of the pure Yang-Mills case shown in [16]. This indicates that the infrared gluon and ghost sector is insensitive to the quark flavor effect except the case of "massless quark" m q = 0. As shown in [16], we find β λ /λ = β M 2 /M 2 = λ/3, which leads as µ → 0 to, The quark mass m q runs according to which indicates absence of logarithmic correction to the quark mass in the infrared. Since u = M 2 /µ 2 → ∞ and t = m 2 q /µ 2 → ∞ exponentially increase more rapidly than λ ∼ 1 ln µ → 0, the above evaluation is consistent. From (57) and (54), the infrared safe trajectories describe the decoupling solution, i.e., the massive gluon and the massless free ghost in the IR limit as those in the pure massive Yang-Mills model of [16].
Finally, we examine the spectral functions. In the gluon vacuum polarization, the quark loop affects the gluon spectral function in the region of (2m q ) 2 < k 2 < ∞ as h q (t) has the branch cut only on 0 <t < 1/4. Therefore, the infrared spectra of the gluon and ghost are the same as those in the pure Yang-Mills case within the one-loop level because the quark mass will be finite in the infrared as can be seen from (58).
As s → 0, the vacuum polarizations arê From the assumption that the higher loop terms are suppressed as µ → 0, i.e. the running u and t do not invalidate the perturbation theory in λ → 0, the following approximation holds in the IR limit |k 2 | → 0, where Z A (|k 2 |, µ 0 ) and Z C (|k 2 |, µ 0 ) are the renormalization factors and Z A (µ 2 , µ 2 0 ) = λ0 M 2 (µ 2 ) from the non-renormalization theorems. As µ 2 → 0, Therefore, the spectral functions originating from the branch cut on timelike momenta have the IR asymptotic forms σ 2 → +0: for gluons, and for ghosts, These results demonstrate the infrared negativity of the gluon spectral function and the infrared positivity of the ghost spectral function for the infrared safe trajectories with m q > 0. The ghost spectral function has a delta function at σ 2 = 0 with a negative coefficient associated to the negative norm massless ghost state and shows the finite positive spectrum in the limit σ 2 → +0. Although we assume that the running u and t do not invalidate the asymptotic perturbative expansion in λ → 0, the infrared negativity of the gluon spectral function will be a remarkable consequence from the infrared safety.
As an alternative evidence for the infrared negativity of the gluon spectral function, we can utilize the infrared behavior of the Euclidean propagator. For example, although this relation gives a trivial result in this model, it is worthwhile to note that they are related . Let us consider the infrared asymptotic behavior of the gluon propagator. The Euclidean propagator (54) can be rewritten as [16], from which we have where we have used γ C → − λ 2 u −1 ln u + O(λu) and λ(µ) − 6 ln µ 2 as µ → 0. Therefore, as k E → 0, we have This logarithmic divergence is precisely consistent with the infrared spectral negativity of ρ(σ 2 ) ∼ −σ 2 < 0. Indeed, if ρ(σ 2 ) ∼ −σ 2 < 0, then, from the (generalized) spectral representation, as k E → +0, which shows the coherency between the approximation of (60) in the complex k 2 momentum and the infrared asymptotic behavior of the Euclidean propagator of (54). Both of the two arguments imply the infrared negativity of the gluon spectral function. As negativity of a spectral function in a weak sense leads to complex poles [40], the IR and UV negativity of the gluon spectral function supports the existence of complex poles in the gluon propagator.
For the quark propagator, the one-loop expression gives no non-trivial infrared spectrum from (25). This is rather clear from the facts that the Feynman diagram indicates Im Σ has non-vanishing value only for k 2 > m 2 q and that m q (µ) is constant in the infrared limit.
Incidentally, let us comment on the non-trivial infrared fixed point to give the Gribov-type scaling solution [16]. In the massive Yang-Mills model of the pure Yang-Mills theory, a non-trivial fixed point in (λ, u) plane appears, where the gluon and ghost exhibit the Gribov-type behavior. One can find a similar fixed point in the massive Yang-Mills model with "massless quarks" m q = 0. For instance, the fixed point lies at (λ, u, t) = (4.0, 1.6, 0) in N F = 3. However, this "scaling" fixed point has an infrared unstable direction. Therefore, with quarks m q > 0, the unstable "scaling" fixed point disappears due to the running of the quark mass.

IV. NUMERICAL RESULTS ON ANALYTIC STRUCTURES
In this section, we present results on the number of complex poles of the gluon and quark propagators in the massive Yang-Mills model with quarks at the oneloop approximation and its RG improvements. Without quarks, we found the gluon propagator has N P = −N W (C) = 2 for all parameters (g 2 , M 2 ) due to the negativity of the spectral function [40]. Surprisingly, it turns out that the model with many light quarks has N P = 0, 2, 4 regions, depending on the number of quarks N F and the parameters (g 2 , M 2 , m 2 q ). In particular, for 4 N F ≤ 9 light quarks, the N P = 4 region, where the gluon propagator has two pairs of complex poles, covers a typical value of coupling g, e.g., g ∼ 4 in the setting to be described shortly, and the models with 10 ≤ N F very light quarks have the region with no complex poles around the typical value of coupling g.
From here on, we set G = SU (3) and the renormalization scale µ 0 = 1 GeV. With the RG improvements, the best fit parameters reported in [21] are g = 4.5, M = 0.42 GeV, and the up and down quark masses m u = 0.13 GeV for the case of N F = 2 and g = 5.3, M = 0.56 GeV, and m u = 0.13 GeV for N F = 2 + 1 + 1 (with the assumption on the strange and charm quark masses of m s = 2m u and m c = 20m u ). Notice that the "quark mass" m q of this model should not be confused with the current quark mass. The "quark mass" parameter m q will be chosen to reproduce the propagators. Rather, m q will be of the same order of the constituent quark mass, since we renormalized this model at µ 0 = 1 GeV. In particular, the massless quarks will differ from the "massless quarks" m q = 0 due to the spontaneous breakdown of the chiral symmetry.
First, we investigate N W (C) of the strict one-loop gluon propagator, mainly at fixed "typical values" of the parameters g = 4 and M 2 = 0.2 GeV 2 to see a qualitative overview of the analytic structures of this model. Next, we consider the RG improvements of these results and evaluate the gluon spectral function at the best-fit parameters for N F = 2. Comparing the strict one-loop and RG-improved results could support the robustness of N W (C). Furthermore, we discuss the complex poles of the quark propagator and comment on the ghost propagator.

A. Gluon propagator: Strict-one-loop analysis
Based on the strict one-loop gluon propagator (13) combined with (29), (30), and (31) of Sec. III, we compute the winding number N W (C) for N F ≤ 9 according to the procedure (9) of Sec. II. Notice that the strict one-loop gluon propagator satisfies the conditions (i) and (ii) in Sec. II. Since the two-point vertex function Γ As a first step, we investigate the (N F , ξ = m 2 q /M 2 ) dependence of N P at the fixed "typical values" of the parameters to obtain an overview. istence of a Euclidean pole together with the fact that D T (0) > 0 and D T (k 2 ) has no complex zeros.
This result illustrates that the gluon propagator has four complex poles in the blue region (4 N F ≤ 9, 0.2 ξ 0.6). It predicts that the transition occurs from the N P = 2 region to N P = 4 region by adding light quarks since N W (C) is a topological invariant. In the N P = 4 region, the gluon propagator has two pairs of complex conjugate poles if we exclude the possibility of Euclidean poles on the negative real axis.
To see details of the N P = 4 region and its boundary, let us look into positions of complex poles.

Location of complex poles
Here we take a further look into positions of complex poles and the transition of N W (C). We will report the ratio w/v for a position of a complex pole k 2 = v + iw and trajectories of complex poles for varying ξ.
First, we focus on the ratio w/v at the complex pole at k 2 = v + iw, on which we set w ≥ 0 without loss of generality from the Schwarz reflection principle, in order to see whether or not the gluon presents a particlelike resonance. Since the gluon propagator has at most two pairs of complex conjugate poles, it is sufficient to find max w/v and min w/v.
The results are demonstrated in Fig. 4. In general, the ratio of complex poles w/v decreases as N F increases. This implies that the gluon exhibits more particlelike be-havior for larger N F . Note that the ratio min w/v rapidly decreases near the boundary between N W (C) = −2 and N W (C) = −4.
FIG. 4: Contour plots of min w/v (top) and max w/v (bottom) on the (NF , ξ) plane. We color the region of w/v > 1 with lightyellow, 1 > w/v > 0.1 with brown, 0.1 > w/v > 0.01 with blue, and 0.01 > w/v with darker blue. As NF increases, the ratio w/v tends to be small. Second, we examine trajectories of complex poles for N F = 3 and N F = 6. Figure 5 plots the pole location with varying ξ. In the case of N F = 3, no transition changing N P occurs; the pole moves gradually and lowers its real part v as increasing ξ. On the other hand, the trajectory of complex poles is completely different at N F = 6 where the transition occurs. The pole from the position of ξ = 0 moves continuously, but is absorbed into the branch cut at ξ ≈ 0.6. Beforehand, the new pole arises at ξ ≈ 0.2 from the branch cut.
These observations demonstrate that the gluon propagator has a pole at timelike momentum, which could correspond to a physical particle, on the boundary and that the pole bifurcates and becomes a new pair of complex conjugate poles in the N P = 4 region. This appearance of the new pair of complex conjugate poles from the branch cut on the real positive axis is compatible with the rapid decrease of min w/v on the boundary between the N P = 2 and N P = 4 regions in Fig. 4.

(g,M) dependence
We have seen that the gluon propagator changes its number of complex poles when many light quarks are incorporated at the fixed typical values of the parameters of g = 4 and M 2 = 0.2 GeV 2 . Therefore, we have to check whether or not the existence of the transition is insensitive to choices of the parameters (g, M ). We shall compute N W (C) in the three dimensional parameter space (λ = C2(G)g 2 16π 2 , u = M 2 /µ 2 0 , ξ). Figure 6 plots the boundary surfaces of equi-N W (C) volume at N F = 3 and N F = 6 in the three-dimensional parameter space (λ, u, ξ). The gluon propagator has four complex poles inside the green surface (0.2 ξ 0.6) and no complex poles inside the blue surface (with small λ and ξ), and two complex poles except these regions. In addition, since the gluon propagator acquires a pole at timelike momentum on the boundary of the N P = 4 region and its poles will move like Fig. 5, w/v for the pole of the gluon propagator decreases as N F increases. The gluon will become more particlelike as light quarks are added. Incidentally, to compare the results with that of the RG improved gluon propagator, we show in Fig. 7 the two-dimensional slices of the three-dimensional figure at u = M 2 /µ 2 0 = 0.2. Thus our qualitative conclusion will be valid: this model has the transition, and the gluon propagator has four complex poles in the presence of 4 N F ≤ 9 light quarks with m q satisfying 0.

B. Gluon propagator: RG improved analysis
So far, we have studied the gluon propagator in the strict one-loop level, relying on the robustness of the winding number. To make this robustness more reliable and to study the structure of the gluon propagator for N F ≥ 10, it is important to survey the winding number for the RG improved gluon propagator. Moreover, this model reproduces the numerical lattice results with the RG improved propagators at the "realistic" parameters g = 4.5, M = 0.42 GeV, and m q = 0.13 GeV at µ 0 = 1 GeV and N F = 2 [21]. In this subsection, we investigate the analytic structure of the one-loop RG improved gluon propagator for the realistic parameters and its parameter dependence for each number of quark flavors N F . on the positive real axis. This shows its spectral function is quasi-negative, from which NW (C) = −2. Notice that the spectral function exhibits the linear decrease with respect to k 2 in agreement with (62) in the infrared limit.

The RG equation for the gluon propagator is
where α denotes the set of coupling and masses α = (λ, u = M 2 /µ 2 , t = m 2 q /µ 2 ) and Z A (µ 2 , µ 2 0 ) is the renormalization factor computed by the anomalous dimension. We then approximate the two-point vertex function to avoid the large logarithms as Although the RG improvements only for the modulus |k 2 | on the complex k 2 plane may break the analyticity, this will give a better approximation than the strict oneloop approximation for the propagator on the timelike momenta.
First, let us see the gluon propagator for the realistic parameters. Figure 8  of N W (C) under continuous deformation [40]. Since the gluon propagator has no zeros N Z = 0 from the fact that the RG improvement only affects the renormalization factor and the running couplings, we deduce the gluon propagator of this model has one pair of complex conjugate poles as in the pure Yang-Mills case. Incidentally, notice that the spectral function ρ(σ 2 ) decreases linearly with respect to σ 2 in the IR limit σ 2 → +0, which is a general feature for the infrared safe trajectories with m q > 0 in this model as shown in (62).
Next, we investigate the number of complex poles N P in the whole parameter space for each N F . Since N Z = 0 within this approximation as before, N P can be computed as N P = −N W (C). Note that different points on a same renormalization group trajectory in the three di- mensional parameter space, e.g., Fig. 2, provide the same analytic structure. Indeed, they are exactly connected by the scale transformation from the dimensional analysis and the anomalous dimension: under a scaling κ, which shows that D T (k 2 , α(κµ 0 ), µ 0 ) and D T (k 2 , α(µ 0 ), µ 0 ) have the same number of complex poles. Therefore, it suffices to compute in the two-dimensional slice of the renormalization group flow, e.g., Fig. 2. From here on, we employ the two dimensional slice at u = M 2 /µ 2 = 0.2. Fig. 9 and Fig. 10 reveal the results for N F = 3, 6, 9, 10. Note that, since the analytic structure of the gluon prop-agator approaches to that of the pure massive Yang-Mills model as m q → ∞, the gluon propagator has one pair of complex conjugate poles N P = 2 for large ξ. At N F = 3, the N P = 2 region dominates the region of the parameters. At N F = 6, in the presence of light quarks, the N P = 4 region, on which the gluon propagator has two pairs of complex conjugate poles, occupies larger region, especially a typical coupling λ ∼ 0.3. On the other hand, at N F = 9, the N P = 0 region, on which the gluon propagator has no complex poles, expands and the N P = 4 region shrinks. At, N F = 10, notably, the N P = 4 region disappears, and the N P = 0 region appears in the larger region of λ for the very light quark masses m q /M 1. Note that the strict one-loop gluon propagator and RG improved one provide almost the same result for N F = 3, 6 by the comparison between Fig. 7 and Fig. 9. This supports the robustness of the winding number N W (C).
It is remarkable that the parameter dependence of the analytic structure drastically changes between N F = 9 and N F = 10. The extension of N P = 0 region for very light quark masses can be seen from the positivity of the UV spectral function for N F ≥ 10 [46], since the positivity in a weak sense indicates N W (C) = 0 [40]. In a similar sense, the disappearance of N P = 4 region could be related to the UV positivity of the spectral function, which might lower N W (C).

C. Quark propagator
As pointed out in [21], the one-loop RG computation is not enough to reproduce the quark wave function renormalization, although the one-loop quark mass function exhibits a qualitative agreement with lattice results. The higher loop corrections or other ignored effects are thus highly significant to describe the quark sector well. Here, we try to examine the quark propagator within the framework of the one-loop RG to catch the qualitative feature as a first attempt.
Let us consider the vector and scalar parts of the quark propagator: We define the spectral functions associated to these propagators, Note that the positivity of the state space would imply  We use the same improvement scheme as (69). The scalar and vector parts of the quark propagator at the realistic values of the parameters are plotted in Fig. 11. Both of the propagators have negative spectral functions and therefore have one pair of complex conjugate poles like the gluon propagator. Notice that both the conditions for the positivity (73) are violated. The former one can be seen from the vector part of the quark propagator in Fig. 11. Figure 12 shows the violation of the latter condition.
Next, we investigate the number of complex poles with various parameters (N F , g, M, m q ). The results on the winding number N W (C) at N F = 3 and N F = 12 are shown in Fig. 13. We numerically checked that both of the scalar and vector parts of the quark propagator (71) yield the same N W (C) on the region shown in Fig. 13, from which N Z = 0 for the both parts.
The quark sector does not exhibit a new structure by adding N F quarks. The qualitative insensitivity of the quark sector to N F is evident in the sense that the strict one-loop expression is, of course, independent of N F and that the quark propagator will only be influenced indirectly by N F via the running parameters. However, it is worthwhile to note that the region with N P = 0 extends slightly as N F increases, in particular for very light quarks ξ 1, which is a common feature to the gluon case.

D. Ghost propagator
Finally, let us add some comments on the ghost propagator. The strict one-loop propagator is not affected by the dynamical quarks. Therefore, from the proposition (Case III) of [40], the value N W (C) = 0 for ghosts will hold after the RG improvements unless the RG trajectory has a Landau pole. Therefore, we conclude the analytic structure of the ghost propagator within this approximation is insensitive to N F and have no complex poles.

V. CONCLUSION AND DISCUSSION
Let us summarize our findings. The argument principle for a propagator relates the propagator on the timelike momenta and the number of complex poles. The winding number N W (C) = N Z −N P can be computed numerically from a propagator on time like momenta according to (9), where N Z and N P stands for the number of complex zeros and poles respectively.
To study the analytic structures of the QCD propagators, we have employed an effective model of QCD, the massive Yang-Mills model. We have found that the infrared-safe trajectories, on which the running coupling is finite in all scales, reproduce the UV asymptotic behavior (55) originally obtained by Oehme and Zimmermann [46] and that the gluon spectral function is negative in the infrared limit ρ(σ 2 ) ∼ −σ 2 for quarks with the "quark mass" of this model m q > 0 (62). The UV negativity of Oehme-Zimmermann and IR negativity of the gluon propagator argued in this paper supports the existence of complex poles in the gluon propagator from the general relationship claiming that a negative spectral function in a weak sense leads to complex poles [40].
In the "realistic" parameters used in [21] to fit the numerical lattice results, the gluon and quark propagators have quasi-negative spectral functions and one pair of complex conjugate poles, while the ghost propagator has no complex poles.
Finally, we have investigated the number of complex poles in this model with many quarks by computing N W (C) for various parameters. For the gluon, there are several interesting features. First, the N P = 2 region dominates the parameter region if ξ = m 2 q /M 2 is not so small or N F 4. In this region, the gluon propagator has one pair of complex conjugate poles as in the pure Yang-Mills case [40]. Second, the N P = 4 region, where the gluon propagator has two pairs of complex conjugate poles, expands by adding light quarks. A typical set of the parameters (g ≈ 4, M 2 ≈ 0.2 GeV 2 ) is covered by the N P = 4 region when 4 N F ≤ 9, 0.2 ξ = m 2 q M 2 0.6. This feature holds in the computations both for the strict one-loop and the RG-improved one-loop results. Third, the strict one-loop result on the locations of complex poles indicates that the gluon tends to be "particlelike" as N F increases. Fourth, for N F ≥ 10, the one-loop RG-improved gluon propagator has no N P = 4 region. Moreover, the gluon propagator with N F ≥ 10 very light quarks ξ = m 2 q M 2 1 has no complex poles even for a "typical value" of the gauge coupling.
The existence of complex poles invalidates the Källén-Lehmann spectral representation, which is a fundamental consequence of QFT describing physical particles. Therefore, the existence of complex poles of the gluon and quark propagators could be a signal of confinement of the elementary degrees of freedom in QCD.
In conclusion, the above results imply that the confinement mechanism may depend on the number of quarks and their mass since complex poles represent a deviation from physical particles and will be related to confinement. Furthermore, the drastic change between N F = 9 and N F = 10 quarks could be possibly related to the "deconfinement" in line with the conformal window.
Several comments are in order. First, let us mention a comparison with a similar approach [38,39], where the analytic structures of gluon and quark propagators are investigated with a massivetype model in light of the variational principle and optimization. In there, the gluon propagator has two pairs of complex conjugate poles, while the quark propagator has a physical pole and no complex poles in N F = 2 QCD. These structures are different from ours, in which the gluon and quark propagators have one pair at the "realistic" parameter. Although the quark sector of both models lacks accuracy, the difference will be relevant in light of the confinement of quark degrees of freedom 3 ; a timelike pole might correspond to a physical one particle state even after some confinement mechanism works. In this sense, the absence of a timelike pole will be favored.
Second, we comments on the N F dependence of the condensation. We have studied the mass-deformed model as an effective model for QCD or QCD-like theories based on the facts that the gluon mass can minimally improve the Gribov ambiguity and that the effective potential for the operator A µ A µ by the local composite operator technique indicates the condensation of this operator [12,13]. The former one can give a mass-like effect independently upon N F . However, the latter argument will be substantially affected by the presence of quarks. The effective potential of [13] appears to be "unbounded" for γ 0 > 0, or N F ≥ 10 for G = SU (3). The "unboundedness" follows from the fact that the coefficient of the Hubbard-Stratonovich transformation ζ runs into the negative infinity ζ → −∞ in the UV limit µ → ∞ due to additive counterterms for N F ≥ 10 even if it is set to be a positive value at some scale. This problem might indicate the limitation of the perturbative treatment for the local composite operator. Therefore, we have no cogent argument supporting the dimension two gluon condensate for N F ≥ 10.
Third, related to the second remark, the validity of our results will be questionable for large N F . For instance, two-loop corrections will be important if the first coefficient of the beta function is small as in the original argument of the infrared conformality [25,26]. Moreover, the results from the truncated Schwinger-Dyson equation [29] shows that the gluon and ghost propagators seem to obey a scaling-type power law in a wide range of momentum in the conformal window; the description by the massive Yang-Mills model will be inappropriate above the critical value of N F . However, we can expect that the massive Yang-Mills model will be valid in the QCD-like phase, or below the critical value of N F . Therefore, the massive Yang-Mills model may capture some information on the transition from the QCD-like phase.
Fourth, N F = 10, where the first coefficient of the gluon anomalous dimension changes its sign, is the value at which the analytic structure of the gluon propagator changes drastically in our analysis. This value appears in various perspective. For example, there has been some proposal that the negativity of the gluon anomalous dimension is crucial for confinement [51]. N F ≈ 10 can be the critical value of the conformal phase transition [30].
Finally, the investigation of the analytic structures by model calculations is speculative and should be taken as an attempt towards capturing some aspects of the intricate dynamics of QCD. Many works of literature have different claims. For example, reference [50] claims that the quark propagator seems to have a pole at timelike momentum while the gluon propagator complex poles by using some parameterizations for the propagators and the numerical solution of truncated Dyson-Schwinger equation. The gluon propagator obtained by solving the truncated Dyson-Schwinger equation on the complex momentum plane in pure Yang-Mills theory is shown to have no complex poles [43]. A recent reconstruction technique indicates complex poles in the gluon propagator [42]. Pursuing the rich analytic structure of the QCD propagators deserves further investigations because QFT describing confined particles is not yet well understood.
On a formal side, local QFT cannot yield complex poles in the standard perspective, see e.g., [52]. One might assert that complex poles correspond to shortlived excitations and break the locality and unitarity in the level of propagators [34,35]. However, if we analytically continue 4 a propagator with complex poles not in the complex momentum but in the complex time from Euclidean space to Minkowski one, the resulted propa-gator can be interpreted as a propagator of a QFT with an indefinite metric state space having complex spectra that can satisfy local commutativity, e.g., [53]. Notice that such theories with complex spectra and an indefinite metric are out of the scope of the axiomatic quantum field theory because their Wightman functions are not tempered distribution due to complex energies, and the theorems derived by assuming the temperedness are not applicable to the "propagators with complex poles". Then, complex poles would not lead to the non-locality and just represent unphysical degrees of freedom. Further discussion on this issue is reserved for future works. As discussed in [36,54], it would also be interesting to study how the complex poles are "canceled" in the physical propagator.

2π
Arg Let us derive this expression. First, we decompose the path around the positive real axis C 2 into three pieces C 2 = C 2,+ ∪ C δ ∪ C 2,− , where C 2,± stands for the path along the positive real axis of C 2,± = {x ± i ; δ 2 < x < Λ 2 } and C δ for the small circle whose center is the origin k 2 = 0. Accordingly, the winding number can be decomposed into the integrals The integral N W (C 1 ) + N W (C 2,− ) + N W (C 2,+ ) can be evaluated as before, N W (C 1 ) + N W (C 2,− ) + N W (C 2,+ ) For the contribution from the small circle, note that D(k 2 ±i ) |D(k 2 ±i )| = e ∓iπα as k 2 → +0. Therefore the phase factor D/|D| varies from e +iπα to e −iπα , from which To sum up, we obtain (A1). Note that the infrared suppression contributes negatively to N W (C) = N Z − N P .