Extrapolation for meson screening masses from imaginary to real chemical potential

We first extend our formulation for the calculation of $\pi$- and $\sigma$-meson screening masses to the case of finite chemical potential $\mu$. We then consider the imaginary-$\mu$ approach, which is an extrapolation method from imaginary chemical potential ($\mu=i \mu_{\rm I}$) to real one ($\mu=\mu_{\rm R}$). The feasibility of the method is discussed based on the entanglement Polyakov-loop extended Nambu--Jona-Lasinio (EPNJL) model in 2-flavor system. As an example, we investigate how reliable the imaginary-$\mu$ approach is for $\pi$- and $\sigma$-meson screening masses, comparing"screening masses at $\mu_{\rm R}$ in the method"with"those calculated directly at $\mu_{\rm R}$". We finally propose the new extrapolation method and confirm its efficiency.


I. INTRODUCTION
T and µ dependence of hadron masses are closely related with those of the ground-state structure of hot QCD matter, where T is temperature and µ means quark-number chemical potential. In fact, medium modification of vector and η ′ mesons has been measured in heavy-ion collisions [1,2]. These results indicate the chiral and the effective U (1) Asymmetry restoration. It is, therefore, important to determine T and µ dependence of light hadron masses.
Lattice QCD (LQCD) is powerful tool to investigate the QCD matter at finite T and µ. In fact, many LQCD calculations have been done for low density (µ/T < 1). The calculation in high density region is still challenging because of well-known "sign problem". Several methods were proposed so far to circumvent the sign problem; the Taylor expansion method [3,4], the reweighting method [5], the imaginaryµ method [6][7][8][9], the canonical approach [10], the complex Langevin method [11][12][13][14], and the Lefschetz thimble theory [15,16]. These have made great progress, but all the results are consistent only in µ/T < ∼ 1 at the present stage. Among them, we pick up the the imaginary-µ method in the present paper. When one considers µ as complex variable, this method corresponds to the analytic continuation from the imaginary chemical potential (µ = iµ I ) to the real one (µ = µ R ).
Meson masses can be classified into "meson pole mass" and "meson screening mass". In LQCD simulations at finite T , the derivation of meson screening mass is easier than that of meson pole mass, since the spatial lattice size is larger than the temporal one; see Appendix of Ref. [42] for the further explanation. Meanwhile, in NJL-type effective models, timeconsuming calculations were needed for the meson screening mass compared with that of the meson pole mass. Recently, this problem was solved by our previous works [42][43][44] for the case of µ = 0.
In this paper, for simplicity, we concentrate on the π-meson and σ-meson screening masses, M scr π and M scr σ , in the framework of 2-flavor EPNJL model. We apply the method of Ref. [42][43][44] for the case of finite µ R and µ I , and then investigate how reliable the imaginary-µ method is for M scr π and M scr σ . For this purpose, we compare " the M scr ξ extrapolated from iµ I (extrapolating result)" with "the M scr ξ calculated directly at µ R (direct result)" for ξ = π, σ mesons.
In Sec. II, we explain a way of calculating the meson screening mass at finite µ. Numerical results are shown in Sec. III. Section IV is devoted to a summary.

A. Model setting
The Lagrangian density of 2-flavor EPNJL model is defined by with u-and d-quark fields ψ = (u, d) T and the isospin matrix τ . We assume isospin symmetry, i.e., u and d quarks have the same mass m 0 . The gluon field A ν is introduced through the covariant derivative where the matrices λ a are the Gell-Mann matrices in color space and g is the gauge coupling. Here, we consider only the time component A 4 of A ν and assume that the A 4 is a homogeneous and static background field.
In the EPNJL model, the Polyakov loop Φ and its Hermitian conjugateΦ are defined by The trace tr c is taken in color space. The relation between A jj 4 and Φ orΦ is summarized in Appendix . The coupling constant G S of the four-quark interaction is assumed to depend on the Polyakov loop Φ andΦ; We set the parameters α 1 , α 2 to α 1 = α 2 = 0.2 to reproduce LQCD data on T dependence of chiral condensate [45] and Polyakov loop [46]; see Sec. III A for the further explanation. The Polyakov loop Φ and its Hermitian conjugateΦ are mainly governed by the Polyakov-loop potential U in Eq. (1). We use the logarithm-type Polyakov-loop potential U of Ref. [26]. The parameter set in U is determined from LQCD data on thermodynamic quantities in the pure gauge limit. The U has one dimensionful parameter T 0 and the value is often set to T 0 = 270 MeV since the deconfinement transition occurs at T = 270 MeV in the pure gauge limit. When one considers the dynamical quarks, the typical energy scale T 0 depends on the number of flavors (N f ). Hence we treat T 0 as an adjustable parameter and determine the value to reproduce the pseudocritical temperature T χ c = 173 ± 8 MeV for chiral transition in 2-flavor LQCD simulations at zero chemical potential [45][46][47]. The parameter thus obtained is T 0 = 200 MeV.
Applying the mean field approximation to Eq. (1) leads to the linearized Lagrangian density where the dressed quark propagator S is defined by with the effective quark mass M = m 0 − 2G S (Φ)σ and the chiral condensate σ = ψ ψ . One can make the path integral over the quark fields analytically, and the thermodynamic potential Ω per unit volume is obtained by with The meanfield variables σ, Φ,Φ are determined so as to minimize the potential Ω. For real µ, we take the approximation Φ =Φ for simplicity. This approximation is pretty good for µ R /T < 1 and not so bad even for µ R /T > 1 [38].
In the µ I region, this thermodynamic potential Ω has the RW periodicity [32,33,37]. The RW periodicity stems from the fact that Ω is invariant under the extended Z 3 transformation [37] defined by for integer k.
The three-dimensional momentum p integral in Eq. (6) has ultraviolet divergence and needs to be regularized. In this paper, we use the Pauli-Villars (PV) regularization [48,49].
where M 0 = M and the M α (α = 1, 2) mean masses of auxiliary particles. The parameters M α and C α are determined so as to satisfy the condition 2 α=0 C α = 2 α=0 C α M 2 α = 0 in order to remove the quartic, the quadratic and the logarithmic divergence in Ω F . We then set (C 0 , C 1 , C 2 ) = (1, −2, 1) and The parameter Λ should be finite even after the regularization (8), since the present model is non-renormalizable.
The EPNJL model has three parameters m 0 , G S (0), Λ in addition to T 0 , α 1 , α 2 . We set m 0 to m 0 = 6.3 MeV and determine G S (0), Λ to reproduce the experimental values of pion mass M π = 138 MeV and its decay constant f π = 93.3 MeV at vacuum. The EPNJL model parameters are summarized in Table I. Following the previous work [43], we first consider π and σ mesons at T = µ = 0. The current operator is expressed by with x = (t, x) for meson species ξ = π, σ, where Γ σ = 1 for σ meson and Γ π = iγ 5 τ 3 for π meson. The mesonic correlation function in coordinate space is defined by Here, the symbol T stands for the time-ordered product. The for an external momentum q = (q 0 , q) andq = ±|q|. When we take the random-phase approximation, we can get χ ξξ as for ξ = π, σ. The one-loop polarization function Π ξ is explicitly calculated by for σ meson and for π meson, where the trace tr c,d is taken in color and Dirac spaces. Three functions in Eqs. (13) and (14) are defined by . (17) These functions are regularized with the same procedure as shown in Eq. (8).
In the two cases of (a) finite T and µ = µ R and (b) finite T and µ = iµ I , one can get the final equations by taking the following replacement The meson screening mass M scr ξ for ξ meson is defined by where the correlation function ζ ξξ (0, x) in coordinate space is obtained by the Fourier transformation of the correlation function χ ξξ (0,q 2 ) in momentum space as see Fig. 1 to understand the meaning ofq integral. NJL-type effective models have two problems in the calculation of Eq. (20). The first problem stems from the regularization. The three-dimensional-momentum cutoff regularization commonly used explicitly breaks Lorentz invariance, and induces unphysical oscillations in ζ ξξ (0, x) [48]. This problem can be solved by taking the PV regularization [49]. We then use the PV regularization in this paper. The second problem is the fact that direct numerical calculations ofq integral is quite difficult because the integrand is highly oscillating at large r where M scr ξ is defined. In order to overcome this problem, one can rewrite theq integral to the complexq integral by using the Cauchy's integral theorem. However, it is shown in Ref. [48] that the complex function χ ξξ (0,q 2 ) has logarithmic cuts in the vicinity of the realq axis. The evaluation of the cuts still demands time-consuming numerical calculations. Our previous works [43,44] showed that the emergence of these logarithmic cuts is avoidable by making the p integration analytically before taking the Matsubara (n) summation in Eqs. (12)- (18).
Consequently, we obtain the regularized function I reg 3 as an infinite series of analytic functions: with a complex valued thermal mass where we take the principle value for logarithm in Eq. (21) and the square root in Eq. (22). Each term in last line of Eq. (21) has four cuts starting atq = ±2iM(M α , ω n , A jj 4 , µ) andq = ±2iM(M α , ω n , −A jj 4 , −µ), as shown in Fig. 1. For later convenience, we define the threshold mass M th and the decay width Γ th by the M located at the lowest branch point in the upper-half plane: Namely where M th (Γ th ) is the real (imaginary) part of 2M lowest . Meson screening mass M scr ξ is a pole of χ ξξ and is calculated by when the pole is located below the lowest branch point. This condition leads to [43] M scr ξ ≤ M th .

III. NUMERICAL RESULTS
A. Deconfinement and chiral transition lines in θ-T plane Figure 2 shows T dependence of σ and |Φ| for the case of θ = 0. The EPNJL-model results with the parameter set of Table I well simulate LQCD data [46,47]. This means that the present EPNJL model is reliable at least for θ = 0. Figure 3 shows the deconfinement and chiral transition lines in the imaginary-µ region, where the transition temperatures are determined from peak positions of chiral and Polyakov-loop susceptibilities. θ dependence of the transition lines are well fitted in 0 ≤ θ ≤ π/3 by using where the superscript "X = d" means the deconfinement transition and "X = χ" corresponds to the chiral transition. The results of the fitting are summarized in Table. II. . The σ is normalized by the value (σ0) at T = 0. LQCD data are taken from Refs. [46,47]. Note that the 10 % errors come from those of Tc. B. θ dependence of π and σ meson screening masses First we have confirmed that πand σ-meson screening masses have the RW periodicity and charge symmetry:  In the next subsection, we will extrapolate the meson screening masses from µ = iµ I to µ = µ R . For this purpose, we first fit our model results with the polynomial function, in 0 ≤ θ ≤ π/3. We take n max = 1, 2, 3, 4 in order to confirm convergence of the expansion. θ dependence of M scr π and M scr σ is well fitted with n max = 4. In this procedure, θ is varied in the range 0 ≤ θ ≤ π/3, although T is fixed.
We consider the following two cases: Fig. 3: The system is in both the deconfinement and the chiral-symmetry restored phase for any θ, since T ≥ T χ c (π/3). (B) T = 180 MeV in Fig. 3: This case satisfies T χ c (0) ≤ T ≤ T RW . The system is in the deconfinement phase for 0 ≤ θ ≤ 0.697 but in the confinement phase in 0.697 ≤ θ ≤ π/3. The system is in the chiralsymmetry restored phase for 0 ≤ θ ≤ 0.403 but in the chiral-symmetry broken phase for 0.403 ≤ θ ≤ π/3.  Figure 5 shows θ dependence of σ-meson screening masses for two cases (A) and (B). The M scr σ have non-monotonic θ dependence for case (B). As for case (A), the πand σ-meson screening masses agree with each other due to the chiral symmetry restoration.

C. Extrapolation from µI to µR region
We compare the extrapolating result with the direct one for finite µ R in order to confirm applicability of the analytic continuation. One can easily make the analytic continuation by replacing θ with −iµ R /T : Figure 6 explains µ R dependence of π-meson screening masses for two cases (A) and (B). In µ R /T < 0.4, the M scr π converge to the direct results as n max increases for both the two cases.  restored in case (A), and θ dependence of M scr σ is almost same as that of M scr π . The extrapolating results tends to the direct ones for 0 ≤ µ R /T < 0.4, and the deviation in 0.4 ≤ µ R /T can not be improved by taking the higher order terms.
The origin of the deviation can be understood when one considers the relation between σ-meson screening mass and chiral susceptibility. Equation (19) indicates that the inverse of M scr σ corresponds to the correlation length in the fluctuation of ψ (x)ψ(x) ; see Ref. [50] for the further explanation, and note that screening mass is referred to be the frequency of "sound mode" there. Hence M scr σ is related to the chiral susceptibility χ σ as Particularly for the chiral limit, µ R and µ I dependence of M scr σ is non-analytic on the chiral phase transition line in µ R - T and µ I -T plane, since χ σ is non-analytic on the chiral phase transition line. As for finite quark mass, a remnant of the nonanalycity makes the accuracy of the analytic continuation less accurate.

D. Phase-transition-line extrapolation
We propose the new extrapolation method by modifying a trajectory of (T, θ) in fitting. In standard extrapolation, θ is varied with fixed T . In new method, we also vary T so that the trajectory runs along the phase transition line. We then assume θ dependence of T as Up to (µ R /T) 4 Up to (µ R /T) 6 Up to (µ R /T) 8 Direct calc. with any constant R that is introduced to cover the θ-T plane; see Fig. 8 for the understanding. The symbol X means the chiral transition (X = χ) or deconfinement transition (X = d). In this paper, we refer to the modified extrapolation as "phase-transition-line (PTL) extrapolation". From now on, we consider the chiral transition (X = χ). We fit θ dependence of σ-meson screening masses with a polynomial series: In Eq. (31), the extrapolation line does not pass through the chiral transition line, we can use all range of θ for fitting, i.e., 0 ≤ θ ≤ π/3.  Up to (µ R /T PTL χ ) 4 Up to (µ R /T PTL χ ) 6 Up to (µ R /T PTL χ ) 8 Direct calc. We then extrapolate M scr σ (θ) and T χ PTL (θ) from finite µ I region to µ R region. Figure 9 shows the comparison between the direct result and the extrapolating ones for M scr σ , where we set T χ PTL (0) = 180 MeV. The extrapolating results rapidly converge to direct-calculated one in µ R /T χ PTL < 0.8. The PTL extrapolation yields better agreement than the standard extrapolation.
We also check the reliability of extrapolation by estimating the radius of convergence in Eq. (32) based on the d'Alembert ratio test. The coefficients b (n) σ in Eq. (32) are summarized in Table. III. The radius of convergence r σ is calculated by ≃ 0.84, whose value is consistent with the upper bound of the agreement region.  Parallel discussion is possible for M scr π , as shown in Fig.  10. We can obtain good agreement between direct results and extrapolating ones for µ R /T χ PTL ≤ π/3.

IV. SUMMARY
We first showed a method of calculating screening masses for finite µ R and µ I in the framework of the 2-flavor EPNJL model.
Next, we investigated how reliable the imaginary-µ approach is for M scr π and M scr σ by comparing "the results extrapolated from imaginary µ" with "those calculated directly in real µ". In the standard extrapolation, the agreement between Up to (µ R /T PTL χ ) 4 Up to (µ R /T PTL χ ) 6 Up to (µ R /T PTL χ ) 8 Direct calc. the direct and the extrapolating results is seen in µ R /T < 0.4 for M scr π and M scr σ for T = 180 and 250 MeV. Especially for σ meson, the disagreement in 0.4 < µ R /T can not be improved by taking higher order terms.
We can understand the difficulty of extrapolation when one remembers that M scr σ is nothing but the inverse of correlation length in fluctuation of local chiral condensate. The M scr σ is thus related with the chiral susceptibility χ σ as M scr σ ∝ χ −1/2 σ . When one set quark mass to zero, χ σ becomes nonanalytic on the chiral transition line T = T χ c (θ), and so does M scr σ . Even for finite quark mass, a remnant of this nonanalicity makes the accuracy of extrapolation less accurate, since quark mass is much smaller to temperature and negligible around the chiral phase transition. This indicates that the simple extrapolation is not useful for M scr σ (T, µ R ). In order to circumvent this problem, we propose the PTL extrapolation. In the method, the agreement between the direct and the extrapolating results is seen in µ R /T χ PTL < 0.8 for M scr σ and in µ R /T χ PTL < π/3 for M scr π with T χ PTL (0) = 180 MeV. The extrapolating results tend to the direct results as higher order terms are taken into account. The PTL extrapolation thus makes better extrapolating results than the standard one.
The difficulty of the simple extrapolation may be in common with other scalar, vector and pseudovector mesons composed of u and d quarks, since these meson masses are sensitive to the chiral transition. The application of PTL extrapolation to such mesons is thus interesting as a future perspective.