Mesonic dynamics and QCD phase transition

We study the nontrivial dispersion relation of mesons resulting from the splitting of the transversal and longitudinal mesonic wave function renormalizations, and its influences on the QCD phase transition, equation of state, and the fluctuations of baryon number. The calculations are performed in a two-flavor low energy effective model within the functional renormalization group approach. We find that influences of the splitting mesonic wave function renormalizations on the equilibrium thermodynamical bulk properties are mild. Furthermore, we have compare the fixed and running expansion for the effective potential in detail, through which the role of the field-dependent meson and quark wave function renormalizations could be inferred.


I. INTRODUCTION
Recent years have seen remarkable progresses in the studies on the QCD phase structure from both the experimental and theoretical sides. Non-monotonic behavior of non-Gaussian cumulants of the net-proton distribution have been observed in Phase I of the Beam Energy Scan (BES) Program at the Relativistic Heavy Ion Collider (RHIC) [1][2][3]. This non-monotonic behavior is potentially related to the critical end point (CEP) in the QCD phase diagram in terms of the temperature and the baryon chemical potential or the baryon density. The CEP is usually believed to be a key feature of the QCD phase structure, which separates the first order phase transition at high density from the continuous crossover at low density.
On the theoretical side, lattice QCD simulation is the first-principle approach to explore the properties of QCD at finite temperature and densities. There have been significant progresses in the lattice studies in recent years, such as the phase boundary [4,5], QCD equation of state at finite baryon chemical potential and the fluctuations and correlations of the conserved charges [6][7][8], the chiral phase transition temperature [9], etc. Because of the sign problem at finite chemical potential, computation of lattice QCD is still restricted to the regime of small chemical potential in the phase diagram, e.g., µ B /T ≤ 2 ∼ 3. Functional continuum field approaches, such as the Dyson-Schwinger equation [10], the functional renormalization group (FRG) [11], etc., are also very suited for the studies of QCD at finite temperature and densities. Functional approaches are not hampered by the sign problem, thus can be employed to investigate the whole phase diagram [10,12,13], in particular to predict the location of CEP, e.g., see [14][15][16][17][18].
The low-energy effective models, e.g. the quark-meson (QM) model [12], Nambu-Jona-Lasinio (NJL) model [47], and their Polyakov-loop improved variants: PQM and PNJL, are suitable to be employed to study the QCD phase transitions. They have been investigated quite a lot in literatures, see, e.g., [48][49][50][51][52][53] for more details. In this work, we adopt the scale-dependent effective action for the two-flavor PQM model, as follows with µ = (0, 1, ..., 3) and i = (1, 2, 3). In Eq. (1) we have used notation x = 1/T 0 dx 0 d 3 x, where the imaginary time formalism for the field theory at finite temperature is used, and the temporal length reads β = 1/T . Apparently, when the temperature is nonzero, the O(4)symmetry in the Euclidean space is broken into that of Z 2 ⊗ O(3), which leads to the split of the magnetic and electric components of correlations functions. They correspond to the components transversal and longitudinal to the heat bath, respectively. In this work, we take this split into account in the two-point correlation function for the mesons, as shown in the second line on the r.h.s. of Eq. (1), where Z φ,k and Z ⊥ φ,k indicate the longitudinal and transversal wave function renormalizations for the temporal and spatial components, respectively.
The reason why we concentrate on the splitting of the wave function renormalization, especially for the mesons, is due to the facts as follows. Firstly, in comparison to the quark wave function renormalization Z q,k and the scale dependent Yukawa coupling h k in Eq. (1), it is found that the meson wave function renormalization Z φ,k plays the most significant role beyond the local potential approximation (LPA) [35,54]. In the LPA, the propagators are classical, i.e., Z q,k = Z φ,k = 1 and the Yukawa coupling h is a constant and independent of the scale k. Secondly, In Ref. [55] calculations based on the full momentum-dependent two-point correlation functions of mesons are compared with those from LPA and LPA , where a momentum-independent Z φ,k is included in LPA in comparison to LPA, and it is found that there is a good agreement between the full momentum calculation and the LPA , rather not LPA, which indicates that the dispersion relation for the meson, resulting from a scale dependent Z φ,k , have already captured most momentum dependence of the two-point correlation function.
Considering the importance of the wave function renormalization for the mesons and the success of LPA , in this work we would like to investigate the effects of the splitting of Z φ,k in the LPA as shown in Eq. (1), which is a natural choice at finite temperature as discussed above. Furthermore, we will also study the interplay between the splitting of Z φ,k and other truncation approaches, e.g., the field dependent Yukawa coupling h k (ρ) which encodes higher order quark-meson scattering processes [54], fixed point expansion for the effective potential V k (ρ) in Eq. (1) versus the running point expansion, etc. Their influences on the QCD phase transition and observables, e.g. fluctuations of the baryon number, will be investigated in detail.
This paper is organized as follows. In Sec. II we give the flow equations for the effective potential, anomalous dimensions, and the Yukawa coupling. The equation of state and the baryon number fluctuations are presented in Sec. III. In Sec. IV we show our numerical results, and finally a summary and outlook is presented in Sec. V. Explicit expressions for the anomalous dimensions and the flow of the Yukawa are given in Appendix A, and some useful threshold functions are collected in Appendix B.

II. FUNCTIONAL RENORMALIZATION GROUP AND FLOW EQUATIONS
To proceed, we describe other notations in the effective action in Eq. (1). φ = (σ, π) is a meson field with four components. The effective potential V k (ρ) with ρ = φ 2 /2 is O(4) invariant and the c-term cσ breaks the chiral symmetry explicitly. The mesons interact with quarks through the scalar and pseudo-scalar channels with a mesonic field dependent Yukawa coupling h k (ρ), and T 0 and T i are the generators in the flavor space with the convention as follows: Tr(T i T j ) = 1 2 δ ij and Besides the wave function renormalization for the meson, we also introduce one for the quark, i.e., Z q,k . Since it plays a minor role in the chiral dynamics in comparison to Z φ,k , the splitting of Z q,k into the transversal and longitudinal components are neglected for simplicity in calculations. Finally,μ = diag(µ u , µ d ) in the first line on the r.h.s. of Eq. (1) denotes the matrix of the quark chemical potential, and µ = µ u = µ d is assumed in the following unless stated specifically. A 0 is the temporal gluon background field, through which the Polyakov dynamics is taken into account.
The RG scale k in Eq. (1) is an infrared cutoff, below which quantum fluctuations are suppressed in the effective action. The evolution of the effective action with k is described by the Wetterich equation [56], which reads with the RG time t = ln(k/Λ) and the initial ultraviolet (UV) cutoff Λ, where R φ k and R q k are the regulators for the meson and quark fields, respectively, and they are given in Eqs. (B1) and (B2). The scale dependent meson and quark propagators are given by with Φ = (q,q, φ) denoting all species of fields. Inserting the effective action in Eq. (1) into the flow equation (2), one arrives at the flow equation for the effective potential, which reads with the RG invariant dimensionless meson and quark masses as follows The threshold functions l Note that z φ ≡ Z φ,k /Z ⊥ φ,k enters into the flow of V k (ρ) through the mesonic fluctuations as shown in Eq. (4).
The transversal anomalous dimension for the π-meson is obtained by employing the projection as follows and the longitudinal one reads Their explicit expressions are given in Eqs. (A1) and (A2), respectively. The difference of the anomalous dimension between the π and σ mesons is neglected here. Note that even they are different, choosing the π-meson anomalous dimension for η φ , as done in this work, could minimize the errors of calculation, at least when the baryon chemical potential is not very high, since the mass of the pion is less than that of σ-meson, because of its nature of Goldstone particle, and the number of the degrees of freedom for the pions is also larger. Therefore, the mesonic degrees of freedom are dominated by the pions. But we should mention that, in the region of high baryon chemical potential in the phase diagram, especially near the CEP, the σ-mode is the most relevant collective mode and the mass of the σ-meson is vanishing, it is necessary to distinguish the anomalous dimensions of π-and σ-mesons, which will be investigated in the future.
The anomalous dimension for the quark is obtained by projecting the inverse quark propagator onto the vector channel, as follows where the spatial component of the external momentum is chosen to be vanishing, as same as the mesonic anomalous dimensions in Eqs. (8) and (9). Note that because of the fermionic property of the quark, its lowest Matsubara frequency is nonvanishing, and we denote it here as p 0,ex , which is described in Appendix A. Projection of the inverse quark propagator onto the scalar channel leads us to the flow of the Yukawa coupling [54], which reads The analytic expressions of Eq. (10) and Eq. (11) are given in Eq. (A3) and Eq. (A4), respectively.

III. EQUATION OF STATE AND BARYON NUMBER FLUCTUATIONS
In the PQM model (1) the thermodynamical potential density reads where all the fields are on their respective equations of motion, see, e.g., [35] for more details. L is the traced Polyakov loop andL is its complex conjugate. They are related to the temporal gluonic background field A 0 in Eq. (1) through the equations as follows with the Polyakov loop P(x) which reads where the path-ordering operator P has been employed.
In this work we employ the parametrization of the glue potential in [57], which reads with the Haar measure whereV glue = V glue /T 4 is dimensionless. The temperature dependence of the glue potential is encoded in the coefficients in Eq. (15), which are given by for x ∈ {a, c, d}, and with the reduced temperature t r = (T − T c )/T c . The parameters in the glue potential are determined by fitting the thermodynamics, the Polyakov loop and its fluctuations in the Yang-Mills (YM) theory, and see [57] for their values. It has been found that the unquenched effect on the glue potential in QCD is well captured from that in the YM theory [58], by employing the rescale for the reduced temperature as follows 0 100 200 300 400 500 600 700 = 250 MeV is adopted in this work. The pressure and the energy density are given by with the entropy density and quark number density reading respectively. The interaction measure or the trace anomaly is given as follows In this work we will also investigate the fluctuations of the baryon numbers, which are obtained through highorder derivatives of the pressure w.r.t. the baryon chemical potential µ B = 3µ. The n-th order baryon number fluctuation reads  which is also called as the n-th order generalized susceptibility of the baryon number. χ B n 's in Eq. (25) are related to the cumulants of the baryon number distributions, e.g., up to the fourth order. Here ... denotes the ensemble average and δN B = N B − N B . In terms of observables in experiments, one has the mean value

IV. NUMERICAL CALCULATIONS AND RESULTS
We employ the approach of Taylor expansion to solve the flow of the effective potential in Eq. (4) numerically. Expanding around the expansion point κ, one is led to where a subscript k is affixed to the expansion point κ, indicating that the expansion point can be RG scale dependent, and N v is the highest order which is included in the calculations. The convergency is obtained when the calculated results are almost no longer influenced by the increment of N v , after it is above some value. It is more convenient to work with the renormalized variables, to wit,V Two different expansion approaches are commonly used in the literature: one is the fixed bare expansion point with ∂ t κ k = 0, i.e., κ k is independent of the scale, which yields It follows immediately that the second term on the r.h.s. of Eq. (32) is vanishing for the fixed point expansion.
Consequently,λ n,k 's of different orders are decoupled and there is no feedback from the high order coupling to low order flows. Therefore, the convergency property is well controlled in the fixed point expansion. For more discussions about the fixed point expansion, see, e.g. [54]. Another approach is the running physical expansion, which demands that the expansion point is the solution of equation of motion for every value of k, to wit, withc k = c/(Z ⊥ φ,k ) 1/2 being the renormalized strength of the explicit chiral symmetry breaking in the effective action (1). Note that the bare c is scale independent, thus one has ∂ tck = (1/2)η ⊥ φ,kc k . Combining Eq. (32) and Eq. (34), one arrives at the flow of the physical expansion point, which reads In this work, we will investigate the influence of these two expansion approaches on the fluctuations of the baryon number. Furthermore, we also use the Taylor expansion to include the field dependence of the Yukawa coupling, which readsh  with where N h is the order, up to which the Taylor expansion ofh k (ρ) is done. Inserting Eq. (37) into Eq. (11), one is led to through which the Yukawa couplings of high orders can be solved.  Taylor expansion for the effective potential in Eq. (31) and the Yukawa coupling in Eq. (36) are chosen to be N v = 5 and N h = 5, respectively. In Sec. IV B one will see that these values are large enough to obtain convergent results in the fixed expansion, see also [54].
It is left to specify the parameters in the low energy effective models: the UV cutoff is chosen to be Λ = 700 MeV; the effective potential at k = Λ is parameterized as follows with λ Λ =10 and ν Λ = (0.556 GeV) 2 . Furthermore, the strength of the explicit chiral symmetry breaking in (1) is c = 1.97 × 10 −3 (GeV) 3 , and the Yukawa coupling is h Λ = 7.33. Employing the fixed point expansion and all the truncations discussed in this work, including the splitting of the meson wave function renormalization and the field dependent Yukawa coupling, etc., we obtain the physical observables in the vacuum, which read f π =92 MeV, m q =300 MeV, m π =136 MeV, and m σ =454 MeV.
In Fig. 1 we show the running of the transversal and longitudinal wave function renormalizations for the mesons, Z ⊥ φ,k , Z φ,k , with the RG scale k. The results obtained from the calculation without splitting are also presented for comparison. Note that values of the wave function renormalizations at the UV cutoff are irrelevant for the computation, since only the anomalous dimensions, rather the wave function renormalizations themselves, enter into the flow equations. One can see that the splitting of Z ⊥ φ,k and Z φ,k takes place even in the vacuum, which is attributed to 3d regulators used in this work, as shown in Eq. (B1) and Eq. (B2), which lead to partial loss of the O(4) symmetry in the vacuum. We also find that when the splitting is not taken into account,  which is usually adopted in the literature, the difference between Z φ,k and Z ⊥ φ,k is very small. Note that Z φ,k here is also derived from the spatial component, see Eq. (8). Even though, the small difference between Z φ,k and Z ⊥ φ,k has already indicated that the approximation neglecting the splitting is reasonable. In Fig. 2 we investigate the influence of the splitting of the meson wave function renormalization on the meson, quark masses and the π-meson decay constant f π . One observes that there is almost no difference between the two cases with and without splitting, except for a small effect for the σ-meson mass in the low temperature regime. In Fig. 3 we plot the trace anomaly as a function of the temperature with µ B = 0. In the same way, we compare the results with and without splitting of the meson wave function renormalization, and the results show that the effect of the splitting on the trace anomaly is small. In the left panel of Fig. 4, we investigate the splitting effect of the meson wave function renormalization on the pressure at finite temperature and vanishing baryon chemical potential. Only a mild difference between the two cases with and without splitting is observed in the low temperature regime, which is similar with the result of the trace anomaly. However, we should remind that, in comparison to mesons, the pressure is usually dominated by the degrees of freedom of quarks and gluons, here the gluons realized through the glue potential in the low energy effective theory. Therefore, it is necessary to just single out the contribution from the meson sector to the full pressure, where we denote the meson contribution by p φ , and investigate the splitting effect on p φ . The relevant result is presented in the right panel of Fig. 4. As expected, the magnitude of p φ is quite smaller than that of p in the left panel, and the splitting effect becomes more obvious for p φ .
In order to investigate whether the partial loss of the O(4) symmetry in the vacuum, arising from the 3d regulators employed in our work, results in unphysical results at finite temperature. We have also performed all the calculations with a larger splitting of the wave function renormalization, by increasing the magnitude of longitudinal anomalous dimension η φ,k artificially. We find that with the increase of the splitting, its influence on the thermodynamical observables grows as well. In another word, this indicates that although the O(4) symmetry is partially lost in the vacuum in our setup due to the use of the 3d regulators, this defect is not too harmful for extracting the thermal effects resulting from the nontrivial dispersion relations of mesons. Certainly, an improved regulator will be helpful to clarify this issue eventually.
In Fig. 5 we investigate the influence of the splitting meson wave function on the kurtosis of the baryon number distribution, i.e., κσ 2 = χ B 4 /χ B 2 . One observes that when the baryon chemical potential is vanishing, the two calculations give almost identical results. With the increase of the chemical potential, e.g., µ B = 400 MeV, as shown in the right panel of Fig. 5, the splitting meson wave function has begun to play a role, and there is difference between the two calculations, especially in the region near the phase transition. In order to find where the difference comes from. We present the fluctuations of the baryon number of different orders at µ B = 400 MeV in Fig. 6. One can see that the fluctuations of the baryon number of all the orders shown here, has been suppressed a bit by the splitting of the meson wave function renormalization during the chiral phase transition. But the magnitude of the difference is mild. The three columns correspond to three different truncations: "Full" includes all the ingredients of the truncation described in this paper; "η φ = 0" denotes the truncation with both the transversal and longitudinal meson anomalous dimensions be vanishing, i.e., η ⊥ φ,k = η φ,k = 0, and others as same as the "Full" one; "LPA" denotes the standard local potential approximation, where all the anomalous dimensions, including those for the meson and quark, are vanishing, and the Yukawa coupling is just a constant.
In a short summary, in this section we have investigated the influences of the splitting meson wave function renormalization on the phase transition, QCD thermodynamics and EoS, and the fluctuations of the baryon number. We find that the influences are small, except that for the baryon number fluctuations at large baryon chemical potential. It seems that the splitting of the meson wave function renormalization is negligible. However, we should remind that the situation may be different when one calculates, e.g., the spectral function of the meson correlator. Taking the π-meson for instance, the one particle irreducible two point function reads It is apparent that m π in the equation above as well as that in Fig. 2 is just the curvature mass, and the more physical relevant masses are the pole mass, and the screening mass, which are obtained by analytically continuing Eq. (40) from the Euclidean to Minkowski space time, and using the conditions as follow Based on this naïve argument, and inserting Eq. (41) and Eq. (42) into Eq. (40), one is led to m 2 π,pole = m 2 π /Z φ and m 2 π,screening = m 2 π /Z ⊥ φ . One can see that the pole and screening masses are related to the wave function renormalization, and accordingly the relevant spectral function is also affected by the wave function renormalization. Since studies on the influence of the splitting wave function renormalization on the spectral functions are beyond the scope of this paper, we will delay them to the work in the future.

B. Comparison between the fixed and running expansion approaches
In this subsection, we would like to investigate the properties of convergency for the fixed and running expansion approaches. We will vary the value of N v in Eq. (31) to investigate whether the convergency is arrived at for each expansion approach. Moreover, a more important and interesting question is that, when the two expansion approaches adopt the same initial conditions, do they produce the same convergent observables? especially the high order ones, such as non-gaussian cumulants of conserved charges, since any small difference for the two approaches could be amplified in the high order observables. The question proposed here is quite nontrivial, because it is not only related to the self-consistency of the theory, but also one can employ the comparison between the two approaches to investigate the importance of the missing effect in our calculation, such as the field dependence of the wave function renormalization, which is beneficial to the improvement of truncation.
In Fig. 7 we have performed the convergency test for both the fixed and running expansions for the effective potential. We find that the convergency for the fixed expansion is very well, and calculations with N v ≥ 3 have already produces reliable results. This behavior is also observed in [54]. This excellent convergency property for the fixed point expansion is easily understood, because there is no feedback from the high order expansion coefficients to the lower ones due to Eq. (33). The situation for the running expansion is more complicated, and we find that the numerical instability arising from the the feedback is harmful to the convergency, as shown in the left panel of Fig. 7, even the maximal expansion order N v grows up to eighth, oscillations are still observed there. However, the convergency for the running expansion could be improved significantly, if the system is away from the regime of the numerical instability, as we show in the right panel of Fig. 7, where the smaller initial Yukawa coupling suppresses the numerical instability remarkably, and a good convergency is observed for the running expansion approach.
In Fig. 8 we show the π-meson decay constant f π , constituent quark mass m q , and the kurtosis of the baryon number distribution χ B 4 /χ B 2 calculated in the fixed and running expansion approaches. Results obtained with N v = 5 for the fixed expansion, and those with N v =5 through 8 for the running expansion are presented for comparison. Note that N v = 5 is large enough to obtain the convergency for the fixed expansion, see Fig. 7. In the first column of Fig. 8 we present the results obtained with the "Full" truncation in which all the ingredients in Eq. (1) are encoded. One observes that the convergency for the running expansion has already been arrived at with N v ≥ 5, but the discrepancy between the fixed and running expansion is obvious, in particular for f π and χ B 4 /χ B 2 . We have included the field dependence for the effective potential in Eq. (31) and the Yukawa coupling in Eq. (36) in our calculations, but not for the wave function renormalizations or the anomalous dimensions. Therefore, it is reasonable to attribute the discrepancy to the missing field-dependence of the mesonic and quark wave function renormalizations in our computation. In order to make this point more manifest, we turn off the mesonic anomalous dimension, and the relevant results are given in the second column of Fig. 8. As we expect, the difference between the fixed and running expansion becomes smaller, especially for f π and χ B 4 /χ B 2 . The difference for the quark mass increases a bit, which indicates that the field dependence of the quark anomalous dimension also play a relatively small role there. In the third column of Fig. 8, we turn off all the anomalous dimensions and the flow of Yukawa coupling, and employ the standard local potential approximation (LPA). One can see that results from the fixed and running expansion approaches are in good agreement with each other. In Fig. 9 we separate the influence of the gluon dynamics, set the Polyakov loop L =L = 1, and compare χ B 4 /χ B 2 obtained from the fixed and running expansions. In the same way we adopt the LPA, and one can see the two expansions give consistent results.

V. SUMMARY AND OUTLOOK
In this work we have studied the nontrivial dispersion relation of mesons resulting from the splitting of the transversal and longitudinal mesonic wave function renormalizations, and its influences on the QCD phase transition, equation of state, the fluctuations of baryon number, etc. The calculations are performed in a twoflavor low energy effective model within the functional renormalization group approach.
It is found that the influences of the splitting mesonic wave function renormalizations on the QCD phase transition, QCD EoS, and the fluctuations of baryon number are small, except that the fluctuations of the baryon number at high baryon chemical potential are suppressed a bit during the chiral phase transition by the splitting. Note that although the splitting effect on the equilibrium thermodynamical bulk properties is mild, it would potentially play a remarkable role for other observables related to the nonequilibrium dynamics, such as the spectral function of correlators, as discussed at the end of Sec. IV A, which will be investigated in our future work.
Furthermore, we have investigated the property of convergency for the fixed and running expansion for the effective potential. Since there is no feedback from the high order expansion coefficients to the lower ones, the convergency for the fixed point expansion is very well. On the contrary, the convergency could be spoiled, if there is numerical instability in the calculation by employing the running expansion for the effective potential, but the convergency can be improved significantly, once the calculation is performed far away from the regime of numerical instability.
Through the comparison between the results obtained from the fixed expansion and those from the running expansion, we find that the field dependence of the meson and quark wave function renormalizations, in particular the former, are important for the quantitative computations, which will be done in our future work.

Appendix B: Threshold functions
In this work we use the 3d-flat or Litim regulators [59,60], which are very suited for the computations at finite temperature and densities, since the summation for the Matsubara frequencies is not affected by the 3d regulators, and can be performed analytically. The regulators in Eq. (2) read with where Θ(x) is the Heaviside step function. Note that since Z ⊥ φ,k = Z φ,k , it is the transversal wave function renormalization for the meson appearing in Eq. (B1).
The definition of the threshold functions B (n) and F (n) is given by F (n) (m 2 q,k ; T, µ) = T k nq G q (q,m 2 q,k ) n .