Effects of quark chemical potential on analytic structure of gluon propagator

We perform complex analyses of the gluon propagator at non-zero quark chemical potential in the long-wavelength limit, using an effective model with a gluon mass term of the Landau-gauge Yang-Mills theory, which is a Landau-gauge limit of the Curci-Ferrari model with quantum corrections being included within the one-loop level. We mainly investigate complex poles of the gluon propagator, which could be relevant to confinement. Around typical values of the model parameters, we show that the gluon propagator has one or two pairs of complex conjugate poles depending on the value of the chemical potential. In addition to a pair similar to that in the case of zero chemical potential, a new pair appears near the real axis when the chemical potential is roughly between the effective quark mass and the effective gluon mass of the model. We discuss possible interpretations of these poles. Additionally, we prove the uniqueness of analytic continuation of the Matsubara propagator to a class of functions that vanish at infinity and are holomorphic except for a finite number of complex poles and singularities on the real axis.


I. INTRODUCTION
For a long time, it has been expected that quark degrees of freedom would dominate in a highly dense matter of quantum chromodynamics (QCD) rather than hadrons, although details of the phase structure are still unclear mainly due to the sign-problem [1]. Studying the analytic structure of the gluon propagator is of importance to this end since this structure provides information on the in-medium behavior, e.g., whether or not a quasi-particle description is appropriate. We thus explore the analytic structure of the gluon propagator in this article.
The main difficulty in a continuum approach for the quark matter is the breakdown of the perturbation theory in the infrared QCD. Indeed, perturbative calculations of dense QCD matter suggest that the quark matter is already strongly correlated at the quark chemical potential µ q 1 GeV [2]. Therefore, a method valid in infrared is required to study relatively low density region of the quark matter. About ten years ago, an effective model of the Landaugauge Yang-Mills theory has been proposed [3,4] to understand recent numerical lattice results that support the decoupling (massive-like) solution of the Dyson-Schwinger equation (DSE) [5,6]. This model consists of the Faddeev-Popov Lagrangian and the simple gluon mass term, i.e., the Landau gauge limit of the Curci-Ferrari model [7], which we call the massive Yang-Mills model. This mass deformation could be a consequence of generating the dimension-two gluon condensate [8][9][10][11][12] or avoiding the Gribov ambiguity [13,14]. The massive Yang-Mills model has the modified Becchi-Rouet-Stora-Tyutin (BRST) symmetry and is multiplicatively renormalizable to be proved through the modified Slavnov- * Electronic address: yhayashi@chiba-u.jp † Electronic address: kondok@faculty.chiba-u. jp Taylor identities (at all orders of the perturbation theory) [7]. Moreover, there exists the "infrared safe" renormalization scheme, which respects the non-renormalization theorems [15,16], in that the running gauge coupling constant g is finite at all scales on some renormalization group (RG) flows [4,14,17]. This model provides the gluon and ghost propagators that agree strikingly with the numerical lattice results just in the one-loop level [3,4]. The three-point functions [18] and two-point correlation functions at finite temperature [19] were compared to the numerical lattice results with good accordance. Furthermore, the two-loop corrections improve the agreement for the gluon and ghost propagators [20]. Therefore, the effective mass captures some nonperturbative aspects of the Yang-Mills theory.
For unquenched cases with the number of quark flavors N F = 2, 2 + 1 + 1, the massive Yang-Mills model with dynamical quarks gives the gluon and ghost propagators consistent with the numerical lattice results as well [21]. For the quark sector, higher loop corrections are important in this model [21][22][23]. Also, QCD phases with heavy quarks have been studied in a similar model in the Landau-DeWitt gauge [24]. Despite the shortcoming of the massive Yang-Mills model for describing the quark sector, this model will be useful for the analyses of the gluon propagator.
One might worry about the absence of the nilpotent BRST symmetry. The massive Yang-Mills model certainly suffers from the physical non-unitarity [7,25] as a consistent theory. However, as this model gives the wellapproximating propagator given by a mass-deformation that has some background refereed above, we can still consider the massive Yang-Mills model as a model for describing the gluon and ghost propagators. Therefore, it is interesting to investigate the analytic structure of the gluon propagator of this model.
In [26], the gluon propagator of the massive Yang-Mills model at finite chemical potential (and zero temperature) obtained in the vanishing momentum renormalization scheme, which is not infrared safe, has been com-pared with numerical lattice data [27] for the gauge group SU (2). While the singlet diquark gap can improve the consistency with the lattice results, the agreement with the lattice results is not quite satisfactory for parameters (g, M ) that are independent of chemical potential. If one enables the gluon mass parameter M to depend on the chemical potential, one can obtain a fair agreement between the massive Yang-Mills model and the numerical lattice results. Therefore, although this model may lack some important aspects, it is still worthwhile studying the analytic structure of the gluon propagator at a finite chemical potential by utilizing the massive Yang-Mills model with various model parameters.
In the vacuum case, i.e., vanishing chemical potential µ q = 0, we have investigated the analytic structures of the gluon, quark, and ghost propagators and revealed that the gluon and quark propagators have one pair of complex conjugate poles while the ghost propagator has no complex poles [28][29][30]. Other several models and a reconstruction method also predict such complex poles of the gluon propagator, e.g. [13,[31][32][33][34][35][36][37][38][39]. The DSE with the ray technique had provided the gluon propagator holomorphic except for timelike momenta [40], but the recent study [41] has updated this conclusion and strongly suggested a singularity on the complex momentum plane. Complex poles invalidate the Källén-Lehmann spectral representation [42] and might correspond to unphysical degrees of freedom in an indefinite metric state space [43]. Therefore, the complex poles represent deviations from observable particles and are expected to be related to the confinement mechanism. For example, the connection between complex poles of the fermion propagator and confining potential in three-dimensional quantum electrodynamics has been discussed in [44]. Incidentally, another generalization of the spectral representation taking unphysical degrees of freedom into account is proposed in [45].
In this article, we employ the massive Yang-Mills model with quantum corrections being included within the one-loop level and investigate the analytic structure of its in-medium gluon propagator at finite quark chemical potential µ q . Since we are interested in the longdistance behavior and analytic structure of the gluon propagator on the complex frequency plane, we perform complex analyses on the gluon propagator in the longwavelength limit k → 0. In addition, we consider the uniqueness of the analytic continuation in the presence of complex poles, since we use the Matsubara propagator in the zero-temperature limit and the analytic continuation is in principle not unique before taking the limit.
This article is organized as follows. In the next section, the definition of complex poles of an in-medium propagator and the method for counting complex poles are presented. A proof of the uniqueness in a class of functions having a finite number of complex poles is provided in Appendix A. The massive Yang-Mills model and its one-loop expressions are presented in Sec. III. We detail the vacuum part of the one-loop expressions in Appendix 1: Schematic picture of singularities on the complex z(= k0) plane. We analytically continue a Matsubara propagator D(z = iωn, k) defined at the Matsubara frequencies z = iωn (shown as the dots) to D(z, k) on the complex z plane. The two sets of complex conjugate poles in the z plane represent a pair of complex conjugate poles with respect to z 2 in (1).
B. In Sec. IV and Appendix C, we determine the number of complex poles and their locations in the space of the model parameters and the spectral function at a specific set of model parameters. It turns out that the gluon propagator has one pair of almost real poles in addition to the other pair of complex conjugate poles similar to the vacuum ones at intermediate quark chemical potential. In Sec. V, we discuss possible interpretations of these almost real poles and estimates for slightly large µ q . In Sec. VI, a summary of these findings and future prospects are given.

II. COMPLEX POLES OF IN-MEDIUM PROPAGATORS
In this section, we define complex poles of propagators in medium. Then, we introduce a method to count the number of complex poles, which is utilized in the subsequent sections.

A. Definition
In medium, we compute a Matsubara propagator D(iω n , k) within the imaginary-time formalism, where ω n is the Matsubara frequency and k is the spatial momentum. We consider the analytic continuation D(z, k) on the complex z plane for a fixed k from the Matsub-ara frequencies on the imaginary axis z = iω n . This provides information on the spectrum and is useful for studying linear response, in which the retarded propagator, namely the propagator analytically continued to the real axis from the upper-half plane, plays an important role [46].
For a field describing a physical observable particle, the usual spectral representation holds. The spectral condition forces analytically-continued Matsubara propagator D(z, k) to have singularities only on real axis z ∈ R.
However, the spectral condition may be violated for confined degrees of freedom, since not all states have to be physical. Thus, we can consider the possibility of complex spectra, which need not be excluded in an indefinite metric state space [43]. If a state with complex energy exists, this should correspond to a confined state. Further formal aspects will be discussed elsewhere.
Here, we assume the following generalized spectral representation allowing complex poles for the gluon propagator D(z, k), which is a propagator obtained by the analytic continuation from the Matsubara propagator D(iω n , k) defined at points on the pure imaginary axis of the complex z plane: where ρ(σ, k) is the spectral function, w ( k) is a position of a complex pole, and Z ( k) is its residue for arbitrary but fixed k. Notice that, in the vacuum case, there is a one-to-one correspondence between the propagator D(z, k) analytically continued to the upper-half plane in z and the analytic continuation in the complex k 2 planeD(k 2 ), which has been considered in the previous articles [28,30], in the sense that D(z, k) =D(z 2 − k 2 ).
Since the set of Matsubara frequencies {ω n } has no accumulation points, uniqueness of the analytic continuation is an important problem to be proved. Indeed, there is a well-known theorem saying that the uniqueness holds in a class of functions satisfying (i) D(z) → 0 as |z| → ∞ and (ii) D(z) is holomorphic except for the real axis, i.e., these two conditions are sufficient to determine the unique continuation [47]. Although this theorem cannot be applied to our case due to the existence of complex poles, we can generalize this theorem in a straightforward way. In Appendix A, we present a proof of the uniqueness under the weaker conditions allowing complex poles: Therefore, the uniqueness of the analytic continuation from the Matsubara propagator is valid in a similar sense even in the presence of complex poles. Note that complex poles defined here do not correspond to poles of quasi-particles. This is because the complex poles defined here yield poles in both of the upper-half and lower-half planes in z. While a quasiparticle pole is in the second Riemann sheet in z 2 , the complex pole is in the first Riemann sheet.

B. Counting complex poles
Let us introduce a method to count the number of complex poles based on the argument principle [28,30] to be used in the following sections. We can relate a propagator at real frequencies to complex poles and zeros.
In the vacuum case, we have applied the method to a propagator on the complex k 2 plane. For an in-medium propagator, we can take k 2 as the squared complex frequency z 2 . The statement is as follows.
Then the winding number, which is the difference between the number of complex zeros (N Z ) and poles (N P ) with respect to z 2 , reads Thus the number of complex poles N P is given by For details of the derivation, see [30]. When the three conditions (i), (ii), and (iii) hold, we can numerically compute the number of complex poles N P from the number of zeros (N Z ) and data at the real frequencies {D(x n + i )} Throughout this article, N P denotes the number of complex poles on the z 2 plane, i.e., the number of poles on the (upper-)half plane on the z plane, and "complex conjugate poles" denote those on the z 2 plane. The propagator has 2N P complex poles on the whole z complex plane.

III. MODEL
In this section, we introduce the massive Yang-Mills model, which is regarded as an effective model of the Landau gauge Yang-Mills theory, or the Landau gauge limit of Curci-Ferrari model, and review the one-loop expressions.

A. Massive Yang-Mills model
The Euclidean Lagrangian of the model at N colors with N F flavors is given by [3,4,21] where we have introduced the bare gluon, ghost, antighost, Nakanishi-Lautrup, and quark fields denoted by , and ψ i (i = 1, 2, · · · , N F ) respectively, the bare gauge coupling constant g b , the bare gluon mass M b , and the bare quark mass (m b ) q,i , while f ABC (A, B, C = 1, 2, · · · , N 2 − 1) stand for the structure constants with the generators t A of the fundamental representation of the group G = SU (N ).
The renormalization factors (Z A , Z C , ZC = Z C , Z ψi ), Z g , Z M 2 , Z mq,i for the gluon, ghost, anti-ghost, and quark fields (A µ , C ,C , ψ i ), the gauge coupling constant g, and the gluon and quark mass parameters M 2 , m q,i are introduced respectively as follows: In this article, we consider the two flavor case N F = 2 and employ this model with degenerate quark masses, m q := m q,i , and therefore Z ψ := Z ψi and Z mq := Z mq,i . Notice that the quark mass parameter m q of this model is chosen to fit the propagators obtained from other methods, e.g., numerical lattice results. In particular, the quark mass parameter m q is non-zero even for massless quarks due to the spontaneous breakdown of the chiral symmetry.
The general tensorial structure of the gluon propagator D µν (k E ) reads, from the spatial rotational symmetry and the transversality of the Landau gauge, where k E = (k 1 , k 2 , k 3 , k 4 ) = ( k, k 4 ) is the Euclidean momentum, P T µν and P L µν are the transverse and longitudinal projectors respectively, i.e., and, We define the vacuum part of the gluon and ghost twopoint vertex functions Γ (2) A ,vac , Γ (2) gh,vac as the zero temperature T = 0 and the zero chemical potential µ = 0 limit, where ∆ gh is the ghost propagator. As a renormalization scheme, we adopt the "infrared safe scheme" [4,21] respecting the nonrenormalization theorem Z A Z C Z M 2 = 1 [15]. For the gluon and ghost sector, we impose combined with the Taylor scheme [16] Z g Z

1/2
A Z C = 1 for the coupling. 1 In this renormalization scheme, it turns out that there exist RG flows on which the running coupling constant is always finite in a whole momentum region, which implies that the perturbation theory is valid to some extent.

B. One-loop expressions
Here we review the results of one-loop calculations of the in-medium gluon propagator.
Beforehand, we decompose the vacuum polarization Π µν (k E ) into the vacuum part Π vac µν (k E ) and the matter part Π mat µν (k E ), Π vac µν (k E ) had been calculated in [3,4,21]. For completeness, the vacuum part is presented in Appendix B.
The relation between Π µν (k E ) and D µν (k E ) is given by the further decomposition of Π mat µν (k E ) as follows: in general, the spatial rotational symmetry yields where the last term δΠ µν is spanned by the tensorial structures k E,µ k E,ν and (P µρ t ρ )k E,ν + (P νρ t ρ )k E,µ with t µ = ( 0, 1) and does not contribute to the propagator due to the transversality of the Landau gauge, while the vacuum part can be written as Π vac µν (k E ) = Π vac (k 2 E )P µν . The gluon propagator is thus of the form (9) with the components of the vacuum polarization: .
The matter part Π mat µν (k E ) at zero temperature T = 0 and non-zero quark chemical potential µ q > 0 is the quark-loop contribution; for µ q > m q , [46] Π mat ρρ = 2 s (k E ). Note that, with the RG functions determined by these renormalization conditions, the parameter dependence of the analytic structure of the RG-improved gluon propagator is qualitatively the same as that of the strict one-loop gluon propagator for N F ≤ 9 in the vacuum case. For N F = 3, 6, see Fig. 6 and Fig. 8 of [30].
and Re denotes the real part when k 4 is real, namely, for any function f (ik 4 ). Now, since we are interested in complex mass and longdistance behavior, let us take the long-wavelength limit k → 0 symmetrically. This limit reduces technical difficulties on the analytic continuation significantly.
In the long-wavelength limit k → 0, we have and, where θ(µ q − m q ) is the step function. Then, the gluon propagator D µν ( k → 0, k 4 ) can be written as where Π vac (s) is the vacuum part given in Appendix B (B8), and,Π with Notice that

IV. RESULTS
In this section, we study the analytic structure of the gluon propagator with the one-loop quantum corrections presented in the previous section.
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 and the up and down quark mass parameters in the case of N F = 2.
An important advantage of this model is the existence of the infrared safe scheme, in which there exist RG trajectories whose running coupling constant is finite for all scales in the one-loop level. Thus, we can implement a one-loop RG improvement, which will give a better fitting result.
However, in the vacuum case µ q = 0 [30], the parameter dependence of the analytic structure of the RGimproved gluon propagator is qualitatively the same as that of the strict one-loop gluon propagator. Therefore, while we adopt the infrared safe scheme, we employ only the strict one-loop gluon propagator to study its analytic structure to simplify analyses. to the whole z 2 plane. In terms of (1), Let us check the prerequisites for the claim of Sec. IIB. The gluon propagator takes the form (1), since it has no branch cut except for the real axis as can be confirmed from (22). Thus, it can have only complex poles in the complex plane excluding the real axis. Also, this gluon propagator satisfies the conditions (i) and (ii) in Sec. II B and has no zeros N Z = 0: M 2 ) plane at the set of parameters (28a), which gives the number of complex poles through the relation NP = −NW (C). In the NP = 2, 4 regions, the gluon propagator has one pair and two pairs of complex conjugate poles, respectively. We used N = 8 × 10 5 , xn = (n + 1) × 10 −5 M 2 (n = 0, · · · , N ), and xN+1 = 50M 2 for the discretization (5) and = 10 −9 M 2 for the infinitesimal imaginary part. For larger ζ or ξ, the gluon propagator has one pair of complex conjugate poles.
• The gluon propagator has no zeros N Z = 0, since the inverse of the propagator (22) does not diverge.
Therefore, the number of complex poles can be calculated according to (5) and N P = −N W (C). For the condition (iii) in Sec. II B, we numerically check convergence of the refinement of the discretization. plane with (25). At the vacuum case µ q = 0, the gluon propagator has one pair of complex conjugate poles, namely two complex poles (N P = 2), irrespective of the value m q . The novel N P = 4 region appears for light quarks (ξ 0.5, or m q 0.30 GeV). As the quark chemical potential µ q increases for such light quarks, the number of complex poles becomes four (N P = 4) at slightly above the quark mass m q and backs to two (N P = 2) at ζ ≈ 0.6, or µ q ≈ 0.33 GeV. In the intermediate quark chemical potential, the gluon propagator has four complex poles in complex z 2 plane. For large m q or µ q , the gluon propagator has two complex poles as in the vacuum case.
B. Analytic structure at a specific set of parameters Next, we take a further look at the analytic structure of the gluon propagator at a specific set of parameters. As the N P = 4 region with intermediate µ q will be interesting, let us choose (28a), (28b), µ q = 0.25 GeV, i.e.,  In what follows, we use k 0 to denote the complex variable z: To take a look at the analytic structure of the gluon propagator, let us see its modulus on the complex k 2 0 plane. The modulus of the gluon propagator D T (k 2 0 ) is plotted in Fig. 3. We can observe that the gluon propagator at the given parameters (32) has indeed two pairs of complex conjugate poles. One pair that is clearly visible in the top panel of Fig. 3 is located at k 2 0 /M 2 ≈ 1.4 ± 2.6i, or k 0 ≈ ±0.62 ± 0.36i GeV. The other pair of complex conjugate poles is at The latter pair has very small imaginary part, while the former one is similar to that in the vacuum case. This smallness of the imaginary part is a universal feature not only around the transition, but also on the whole N P = 4 region, as we will see in the next subsection.
The gluon propagator (22) with real k 2 0 and its spectral function ρ(ω) := ρ(ω, k → 0) := 1 π Im D T (ω 2 + i ) (34) are displayed in Fig. 4. The propagator shows a rapid oscillation at k 2 0 /M 2 = −k 2 4 /M 2 ≈ 1.4, or k 0 ≈ 0.5 GeV ≈ (2µ q ). The negative peak of the spectral function has a larger value than the positive one: max ρ ∼ 2.7 GeV −2 and min ρ ∼ −29 GeV −2 . The rapid change is consistent with existence of almost real complex poles. Apart from the sharp peak, the gluon propagator is similar to the vacuum one. The quark chemical potential affects the gluon propagator significantly only around k 0 ≈ 2µ q .

C. Locations of complex poles
Let us investigate locations of complex poles of the gluon propagator for various parameters (ζ = µ 2 q M 2 , ξ = m 2 q M 2 ) with fixed (g, M ) of (28a). We present the ratio ω I /ω R of the real and imaginary parts of a complex pole on the (ζ, ξ) plane and a trajectory of poles for varying µ q and at fixed m q . First, we compute the ratio ω I /ω R to obtain an overview on positions of complex poles on the parameter space (ζ, ξ). We can restrict ourselves to ω R > 0, ω I > 0 without loss of generality from the Schwarz reflection principle and the symmetry k 0 → −k 0 . As the gluon propagator has at most two pairs of complex conjugate poles with respect to k 2 0 , it is sufficient to find max ω I /ω R and min ω I /ω R .
Contour plots of the ratios (max ω I /ω R and min ω I /ω R ) are shown in Fig. 5.
This result is consistent with Fig. 2 as max ω I /ω R = min ω I /ω R only on N P = 4 region, where the gluon propagator D T (k 2 0 ) has two pairs of complex conjugate poles with respect to k 2 0 . These figures indicate that the gluon propagator has a pair of almost real complex poles in the N P = 4 region shown, while the pair with max ω I /ω R is always of the same order of magnitude.
Moreover, in general, the ratio ω I /ω R tends to increase as the quark chemical potential µ q increases, except for the almost real poles. In other words, the gluon propagator becomes "less particlelike" for large µ q .
In the previous subsection, we observed that both the sharp spectral peak and almost real poles appear at Re k 0 ≈ 2µ q (≈ 0.5 GeV) at µ q = 0.25 GeV. This feature is not limited to the specific parameter but universal. Let us examine locations of complex poles at the parameter (28a) and (28b) with varying µ q . Figure 6 plots a trajectory of complex poles on the complex k 0 plane and µ q -dependence of the real parts of the complex poles. As µ q increases, a new pole appears from the branch cut (at µ q ≈ 0.16 GeV), then moves along the real axis, and is finally absorbed into the branch cut (at µ q ≈ 0.33 GeV). On the other hand, the other pole increases its imaginary part gradually. This feature is consistent with the number of complex poles of Sec. IV A.
The bottom panel of Fig. 6 clearly indicates that the real part of the new almost real pole can be approximated by 2µ q : ω R ≈ 2µ q . We have also checked that the almost real poles are at Re k 0 ≈ 2µ q for different values of m q . Before concluding this section, let us consider (g, M ) dependence of the above results, especially, the number of complex poles. For details of these analyses, see Appendix C. We have found that the N P contour plot is not sensitive to a detailed choice of the parameters (g, M ).

E. Summary of results
In summary, we have observed the following points in this section.
• There is a N P = 4 region, where the gluon propagator has two pairs of complex conjugate poles with respect to k 2 0 . See Fig. 2.
• In N P = 4 region, the gluon propagator has an almost real pair of complex conjugate poles at Re k 0 ≈ 2µ q . See Fig. 6 • With almost real poles, the real part and imaginary part (to be identified with the spectral function) have narrow peaks at k 0 ≈ 2µ q . See Fig. 4 • The ratio ω I /ω R of a complex pole k 0 = ω R + iω I , (ω R > 0, ω I > 0) tends to increase as µ q increases, except for the almost real poles. See Fig. 5.

V. DISCUSSION
In this section, we discuss implications of the results shown in the previous sections, especially the appearance of the almost real pole in the N P = 4 region, and comment on estimates of the analytic structure of the gluon propagator for relatively large µ q .

A. Almost real complex poles and spectral function
For the gluon propagator, we found a new pair of complex conjugate poles at Re k 0 ≈ 2µ q with quite small imaginary parts (Im k 0 ≈ 0). Together with the narrow peaks shown in Fig. 4, the quark chemical potential affects the gluon propagator significantly around k 0 ≈ 2µ q .
The importance of the scale 2µ q can be understood by the fact that 2µ q is the lowest energy for the quark pair creation to occur, which contributes to the spectrum of the gluon, due to the Fermi degeneracy. Moreover, in the massive model, the gluon "decouples" at low energies. Thus, quark loop dominates the low-energy region of the gluon spectral function. On the other hand, in the high energy region, the gluon and ghost loops win against the quark loop for the gluon spectral function to yield ρ < 0 in the large frequency limit for N F < 10 [48]. Therefore, 2µ q will be quite an important scale for relatively small µ q (but larger than m q ), while less important in the highenergy region. This might explain the appearance and disappearance of the almost real complex poles as varying µ q .
Since complex poles never appear in the physical spectrum, they should correspond to confined degrees of freedom. The transition between N P = 2 and N P = 4 regions indicates that timelike spectra transform to confined complex degrees of freedom, or vice versa. Therefore, the transition between N P = 2 and N P = 4 regions might have a physical significance on the dynamics of the strong interaction.
Note that, however, the appearance of the almost real pole may be an artifact of the approximation: where the vacuum polarization Π is replaced by the oneloop expression Π 1−loop . For example, in this approximation, even the propagator of the Higgs field in U (1) Higgs model with the small gauge-fixing parameter has complex poles with tiny imaginary parts [49]. The new . This plots ρ(ω) = 1 π Im DT (k 2 0 = ω 2 + i ) with /M 2 = 10 −3 , which is larger than the imaginary part of the almost real pole. This shows that the spectral function has a long-lived quasi-particle peak, if the complex pole is an error.
pole reported in the previous section may be similar to this one. In this case, the almost real pole should be interpreted as a long-lived collective mode with frequency 2µ q .
If the almost real pole is an artifact, an estimate of the spectral function will be given by ρ(ω) = 1 π Im D T (k 2 0 = ω 2 + i ), where is small but larger than the imaginary part of the almost real pole ω I . This estimate is displayed in Fig. 7 at (32) and N F = 2. We take /M 2 = 10 −3 because the complex poles are at k 2 0 /M 2 ≈ 1.4 ± (1.8 × 10 −4 )i. This plot implies that the new "complex pole" may correspond to actually an long-lived quasi-particle.
Finally, note that N P = 4 region is located in the region less than µ q ≈ 0.33 GeV, which is approximately the matter threshold, if exists. Therefore, in any case, the new complex pole or the quasi-particle pole will be in the confined dynamics. Although it might be accidental, it could be interesting that the right side of the boundary between N P = 2 and N P = 4 regions locates approximately at the liquid-gas threshold µ q ≈ 0.33 GeV for all m q 0.33 GeV.
In summary, we again emphasize the following points, • The chemical potential influences the gluon propagator significantly around k 0 ≈ 2µ q . This can be explained by the facts, (1) it is the least energy for the quark pair production without momentum transfer k = 0 and (2) the quark loop is important in the energy scale less than the effective gluon mass in this model.
• If the new pair of complex conjugate poles indeed emerges as µ q increases, there may be a transition on the boundary between N P = 2 and N P = 4 phase.
• On the other hand, the almost real pole may be an artifact of the approximation (36). Then, the gluon propagator would have a quasi-particle spec- tral peak instead of the complex poles, which correspond to confined states.

B. Slightly larger µq
To obtain a fair agreement with lattice results in twocolor QCD, the effective gluon mass parameter M is chosen of order µ q for µ q ∼ 0.6 -1 GeV [26]. As an attempt to obtain an estimate of the analytic structure of the gluon propagator for the slightly large µ q , we investigate it at µ q = M .
Beforehand, let us see how the in-medium modification of the effective gluon mass affects the analytic structure. The real and imaginary parts of the gluon propagator at µ q = 0.8 GeV for M = 0.42 GeV and M = 0.8 GeV are plotted in Fig. 8. The change of M does not largely modify the location of the spectral peak, k 2 0 ≈ (2µ q ) 2 , while the direction of the peak is inverted. The real and imaginary parts of the gluon propagator at M = µ q = 0.6, 0.8, 1.0 GeV are plotted in Fig. 9. The spectral function has a negative peak at k 2 0 ≈ (2µ q ) 2 . The magnitude of this peak decreases as µ q increases. The gluon propagator has one pair of complex conjugate poles as the vacuum one. The effect of the chemical potential around k 0 ≈ 2µ q is less significant for large µ q in this model.
For complex poles, we have numerically confirmed N P = 2 in this set up for µ q ∼ 0.6 -1 GeV as inferred from Fig. 2. Their positions k 0 = ω R + iω I , (ω R > 0, ω I > 0) are plotted in Fig. 10. The apparent linearity of ω R and ω I with respect to µ q (= M ) suggests that M and µ q are the dominating scales in the propagator. A comparison with Fig. 6 indicates that the in-medium modification of the gluon mass makes ω R and ω I larger.

VI. SUMMARY AND FUTURE PROSPECTS
Let us summarize our findings. We have performed complex analyses of the gluon propagator at non-zero quark chemical potential µ q in the long-wavelength limit k → 0, by using the massive Yang-Mills model. We have verified that the two conditions, (i) D(z) → 0 as |z| → ∞ and (ii) D(z) is holomorphic except for the real axis and a finite number of complex poles, are sufficient to sin-gle out the correct analytic continuation of a Matsubara propagator. Therefore, the uniqueness of the analytic continuation is concluded even if we allow the existence of complex poles. For the proof, see Appendix A.
We have found that there is N P = 4 region, where the gluon propagator has two pairs of complex conjugate poles with respect to the complex variable z 2 = k 2 0 . In this region, a new pair appears near the real axis in addition to the other pair similar to that in the vacuum case. At the typical parameters (Fig. 2), the N P = 4 region appears for light quarks (m q 0.30 GeV). As the quark chemical potential µ q increases, the number of complex poles becomes four (N P = 4) at slightly above the quark mass m q and backs to two (N P = 2) at µ q ≈ 0.8M ≈ 0.33 GeV. This structure is not sensitive to details of choice of the parameters (g, M ) as shown in Appendix C. Moreover, in this N P = 4 region, the new pair of complex conjugate poles has quite small imaginary part, and its location is approximately Re k 2 0 ≈ (2µ q ) 2 . On the other hand, in the N P = 2 region, the gluon propagator behaves less "particlelike" with larger ratio ω I /ω R of the complex pole at k 0 = ω R + iω I , as µ q increases.
The chemical potential influences the gluon propagator significantly around k 0 ≈ 2µ q , where the new pole appears and the spectral peak is observed. We can attribute this to the facts (i) it is the least energy for the quark pair production to occur at k = 0 and (ii) the quark loop dominates in the energy scale less than the gluon mass M .
Finally, we can interpret the new almost real poles in two ways. First, the results may imply that the gluon propagator indeed has a new pair of complex poles. This suggests a transition in confined degrees of freedom involving the gluon. Second, the almost real pole may be an artifact of "the one-loop approximation" (36). Then, the gluon propagator would have a long-lived quasi-particle spectral peak instead of the confined complex pole, which suggests a quasi-particle picture of the in-medium gluon.
To sum up, although the gluon propagator presents only mild changes on the Euclidean side [27], it might have a rich and interesting structure in the complex frequency plane.
As future prospects, there is plenty of room for improvement in the present work in many aspects. First, this work does not take into account the quark condensation, which is expected to be essential in the highly dense quark matter. The effect on the analytic structure of the quark gap would be interesting. Second, as remarked in the introduction, the one-loop level is not enough in the quark sector of the massive Yang-Mills model. A possible improvement is the double expansion that improves the quark mass function significantly [23]. Third, while a fair agreement with lattice results can be obtained by making the gluon mass M depend on µ q [26], the medium modification of the effective gluon mass should be determined in a more systematic way. Fourth, since the massive Yang-Mills model has the infrared safe renormalization scheme, it would be important to com-pare the RG improved Euclidean gluon propagator with the lattice one. This could improve the current unsatisfactory agreement. Lastly, when using lattice results, we have to keep in mind that the lattice gluon propagator has non-negligible systematic errors, e.g., finite lattice-spacing effect [50], at low momenta and how Gribov copies affect results because there is no reason of the coincidence between the minimal Landau gauge and the Euclidean version of Landau gauge of the well-known covariant operator formalism due to the Gribov ambiguity.
For other directions, it would be interesting to introduce temperature and to consider the physical sector and its transport properties in the massive Yang-Mills model and compare them with other approaches, e.g., [51]. Although it is very difficult, it is important to discuss implications of complex poles in the physical sector. The corresponding state should be confined and not itself have any physical impact, but its composite state might have physical significance [52]. Formal aspects of complex poles will be discussed in a future work.
In the absence of complex singularities, a theorem for the uniqueness of analytic continuation of the Matsubara propagator is well-known and proved in [47]. In this appendix, we shall extend the theorem to propagators with complex poles.
In practice, we have not faced the problem of the uniqueness of the analytic continuation as we have employed the gluon propagator at zero temperature T = 0. However, it is the low temperature limit T → 0 of the Matsubara propagator at finite temperature; it is conceptually essential to establish the uniqueness of the analytic continuation of the Matsubara propagator.

Theorem
Let D(z) be a complex function whose values at Matsubara frequencies z = iω n := i 2πn β are given. Then, its analytic continuation D(z) to the whole complex z plane is unique provided that an analytic continuation satisfies the following conditions, Let D 1 (z) and D 2 (z) be two analytic continuations satisfying the above two conditions that coincide at all the Matsubara frequencies: • ϕ(iω n ) = 0 for all Matsubara frequencies ω n , • ϕ(z) may have a finite number of poles, • ϕ(z) → 0 as |z| → ∞.
We shall show that ϕ(z) is identically zero, i.e., an assumption that ϕ(z) had only isolated zeros leads to a contradiction. The proof is a straightforward generalization of a proof of the Carleman theorem given in Titchmarsh's book [53].
Consider the integral where the contour C = (ρ, R) ∪ C R ∪ (−R, −ρ) ∪ C ρ is depicted in Fig. 11 and C R = {z; Im z > 0, |z| = R} and C ρ = {z; Im z > 0, |z| = ρ} are the semicircles with counterclockwise and clockwise directions respectively. In this integral, we are going to keep ρ finite and take a limit R → ∞. From here on, we omit +i for notational simplicity.
We take a sufficiently small ρ (or appropriate choice of branch cuts of ln ϕ(z)) so that C ρ does not intersect with any branch cut of the logarithm.
We evaluate this integral Im I(R) in two ways to obtain the contradiction.
First, we decompose the integral I(R) into four pieces following C = (ρ, R) ∪ C R ∪ (−R, −ρ) ∪ C ρ , Then, we have and, Thus, we obtain Im I(R) = Im I Cρ Note that Im I Cρ is O(1) as R → ∞. The other two integrals could diverge as R → ∞; however, then, Im I(R) would be negative infinity, since ϕ(z) → 0 as |z| → ∞ and the other parts of the integrands are positive, 1 On the other hand, the integral I(R) is closely related to zeros and poles inside C .
The first integral sums up 'discontinuities' from the branch cuts of the logarithm. Since we have assumed that the branch cuts of the logarithm do not intersect with C ρ , the first term contributes only from (ρ, R) ∪ C R ∪ (−R, −ρ), on which z R 2 + 1 z is real. Therefore, Moreover, as R → ∞, n iωn∈D ω n R 2 = 0<ωn<R ω n R 2 = O(1). (A11) These results indicate which contradicts the first evaluation: Im I(R) is bounded above. The assumption that ϕ(z) had only isolated zeros is false. Therefore, ϕ(z) = D 1 (z) − D 2 (z) is identically zero at least for the upper-half plane. In the same way, ϕ = 0 in the lowest-half plane follows by taking −z as z. This completes the proof. Incidentally, let us comment on the possibility of branch cuts. The uniqueness holds even if we allow the propagator to have a finite number of (non-closed) branch cuts that have finite length and represent finite discontinuities of D(z). Then, ϕ(z) = D 1 (z) − D 2 (z) could have branch cuts in addition to poles. We can still prove ϕ = 0 by (i) deforming C to avoid the branch cuts and (ii) taking the branch cuts of ln ϕ(z) so that they intersect with neither C ρ nor the path wrapping around the new branch cuts of ϕ.
Indeed, the first evaluation becomes I(R) = I ρ→R + I C R + I −R→−ρ + I Cρ + γ:cuts I γ , where γ is a path that surrounds a cut γ in |z| < R. This new contribution is finite for any R due to the finiteness of the branch cuts.
On the other hand, the second evaluation by the partial integration is the same as before, which leads to a contradiction again. Therefore, the conclusion is not changed in the presence of discontinuities on curves of finite length. while the asymptotic freedom and RG analysis yields where we have analytically continued the gluon propagator from the Euclidean momentum k 2 = −k 2 E to complex k 2 , Z U V > 0 is a positive constant, and γ 0 and β 0 are respectively the first coefficients of the gluon anomalous dimension and the beta function: Both the strict one-loop gluon propagator and RG improved one satisfy the condition (i) of Sec. II B. In spite of the wrong logarithmic exponent, the one-loop gluon propagator has qualitatively the same phase as the RG improved one (for N F < 10). Thus, the wrong logarithmic exponent will not change the value of N W (C) = N Z − N P , and hence the strict one-loop expression may be enough for our purpose. In the main text, we have investigated the analytic structure of the gluon propagator with the fixed parameters g = 4.5 and M = 0.42 GeV, as they give best-fit parameters to the lattice results [21]. In this appendix, we check that the qualitative features of the analytic structure are not sensitive to the model parameters (g, M ).
We have confirmed that the contour plots of N P on the (ζ =