Functional renormalization group analysis of the soft mode at the QCD critical point

We make an intensive investigation of the soft mode at the quantum chromodynamics (QCD) critical point on the basis of the functional renormalization group (FRG) method in the local potential approximation. We calculate the spectral functions $\rho_{\sigma, \pi}(\omega,\, p)$ in the scalar ($\sigma$) and pseudoscalar ($\pi$) channels beyond the random phase approximation in the quark--meson model. At finite baryon chemical potential $\mu$ with a finite quark mass, the baryon-number fluctuation is coupled to the scalar channel and the spectral function in the $\sigma$ channel has a support not only in the time-like ($\omega\,>\,p$) but also in the space-like ($\omega\,<\, p$) regions, which correspond to the mesonic and the particle--hole phonon excitations, respectively. We find that the energy of the peak position of the latter becomes vanishingly small with the height being enhanced as the system approaches the QCD critical point, which is a manifestation of the fact that the phonon mode is the {\em soft mode} associated with the second-order transition at the QCD critical point, as has been suggested by some authors. Moreover, our extensive calculation of the spectral function in the $(\omega, p)$ plane enables us to see that the mesonic and phonon modes have the respective definite dispersion relations $\omega_{\sigma.{\rm ph}}(p)$, and it turns out that $\omega_{\sigma}(p)$ crosses the light-cone line into the space-like region, and then eventually merges into the phonon mode as the system approaches the critical point more closely. This implies that the sigma-mesonic mode also becomes soft at the critical point. We also provide numerical stability conditions that are necessary for obtaining the accurate effective potential from the flow equation.


Introduction
The phase diagram of quantum chromodynamics (QCD) is expected to have a rich structure and its clarification is one of the main topics in high-energy and nuclear physics [1]. One of the remarkable features in the expected phase structure is the existence of the first-order phase boundary between the hadronic phase and the quark-gluon plasma (QGP) phase at large baryon chemical potential μ. In particular, the end point of the first-order phase boundary is known as the QCD critical point, where the phase transition becomes second order. Although lattice QCD, which is a powerful nonperturbative first-principle method for QCD, has a limited predictive power in the case of large μ because of the sign problem [2][3][4], other QCD-motivated approaches, such as chiral effective models implementing relevant symmetries and functional methods with inputs from lattice QCD, support the existence of the QCD critical point [5][6][7][8][9][10][11][12][13][14][15][16].
While the location of the QCD critical point in these calculations strongly depends on details of the models and employed approximations [17][18][19][20], one may expect anomalous fluctuations of PTEP 2016, 073D01 T. Yokota et al. experimental observables in relativistic heavy ion collisions if the produced matter passes through critical region where the thermodynamic quantities are strongly influenced by the existence of the critical point [21]. In particular, the net-baryon number susceptibility and its higher-order ones, approximated by the net-proton ones in measurements [22], are expected to be sensitive to the critical behavior of the system [18,[23][24][25][26][27][28][29][30][31][32]. A recent summary of measurements in the beam-energy scan program at the Relativistic Heavy Ion Collider (RHIC) including intriguing behaviors of the net-proton fluctuations can be found in Ref. [33].
Although the mean-field theory seems to give a reasonable picture of the phase structure, a system near a critical point exhibits strong correlations among its constituents. Therefore methods beyond the mean-field theory are desirable for understanding the physics near the critical point. The functional renormalization group (FRG) [34][35][36][37] is one of the frameworks that enable us to incorporate fluctuation effects beyond the mean-field theory; see also Refs. [38][39][40]. Indeed, it has been pointed out that such fluctuations play an important role in the aforementioned observables [27,[29][30][31][32][41][42][43]. The FRG has been applied to a wide range of fields including hot and dense QCD and found to be useful in the description of chiral phase transition in QCD using effective chiral models [12,20,[44][45][46][47][48].
For a second-order transition, there exist specific modes that are coupled to the fluctuations of the order parameter and become gapless and long-lived at the critical point. Such a mode is called the soft mode of the phase transition. As for the QCD critical point, the nature of the soft modes changes depending on whether the quarks have a mass or not [49,50]: In the chiral limit where the quarks are massless, the theory has exact chiral symmetry and an O(4) critical line appears in the (T , μ) plane in the two-flavor case; the critical line terminates at a point (the tricritical point) in the (T , μ) plane and is connected to a first-order phase transition line [51]; the σ and π mesonic modes in the time-like region become massless on the critical line, implying that the quartet mesons are the soft modes of the chiral transition in the chiral limit [52][53][54]. At finite baryon chemical potential, however, the picture changes because the charge conjugation symmetry is lost. When the chiral symmetry is broken, the mixed correlator ψ ψψγ 0 ψ does not vanish. Then, when the system approaches the tricritical point in the broken phase, the density fluctuation (ψγ 0 ψ) 2 also shows a singular behavior. Furthermore, when the theory does not have chiral symmetry due to the current quark masses from the outset, the natures of the phase transition and the soft mode change drastically. In this case, the tricritical point becomes the critical point where the universality class belongs to that of the Z 2 critical point and the soft mode is considered to be the particle-hole mode corresponding to the density (and energy) fluctuations. It is noteworthy that not only the chiral susceptibility but also the susceptibilities of hydrodynamical modes such as the density fluctuation or the quark-number susceptibility diverge at the critical point, owing to the scalar-vector coupling caused by the finite quark mass at nonvanishing μ [55]. Such a picture is suggested in the time-dependent Ginzburg-Landau theory and random phase approximation (RPA) analysis of the Nambu-Jona-Lasinio (NJL) model in Ref. [49] and the Langevin equation in Ref. [50].
Then it would be intriguing to apply the FRG for investigating the nature and dynamical properties of the soft mode at the critical point, which is the purpose of the present work. Here it should be mentioned that the role of density fluctuations in the static critical properties of the QCD critical point has been investigated in FRG within a chiral quark-meson model [56], where static properties such as the quark-number susceptibility and the curvature (screening) masses in the scalar and vector channels are analyzed. Dynamical properties such as the dispersion relations of excitation modes 2 and those of the soft modes at the QCD critical point have not been touched upon. We note that these quantities can be read off from the spectral function in the relevant channel.
Needless to say, a real-time analysis is needed for the investigation of the spectral functions for excitation modes. Since analytic continuation of two-point functions from imaginary Matsubara frequencies to real frequencies is necessary to get real-time two-point Green's functions at finite temperature, it is in general difficult to obtain the spectral functions numerically with high accuracy [57][58][59][60]. Recently, a useful method for calculating the spectral functions within FRG has been developed, which adopts an unambiguous way of analytic continuation in the imaginary-time formalism and leads to reasonable results of meson spectral functions in the O(4) model in vacuum [61]. Furthermore, this method has also been successfully applied to the quark-meson model at finite temperature and chemical potential [62,63].
The spectral function in a specific channel with specific quantum numbers may have more than one peak and bump, the number of which can change according to that of the parameters characterizing the system such as temperature and baryon density. A peak or even bump in the spectral function in the channel may be identified with a particle excitation in the channel, and the width of the peak/bump shows the decay rate of the particle in a specific decay channel. This feature enables one not only to explore the appearance but also to analyze the nature of modes in the system using spectral functions in various channels. In particular, the spectral functions around a critical point give information on the soft modes. Actually, the spectral function in the scalar-isoscalar channel, i.e., the sigma channel, has been analyzed in the RPA using the Nambu-Jona-Lasinio (NJL) model in Ref. [49], where it is shown that the spectral function has a prominent bump in the space-like region corresponding to the phonon mode composed of particle-hole excitations, the peak position of which moves to vanishingly small frequency as the system approaches the critical point, whereas the mesonic sigma mode in the time-like region does not show such a behavior, retaining a finite mass. This result clearly shows that the phonon mode (or hydrodynamical mode in general) is the soft mode of the QCD critical point; see also Refs. [50,[64][65][66].
The purpose of the present paper is to investigate the nature of low-energy modes at the QCD critical point beyond the RPA using a two-flavor quark-meson model: We calculate the spectral functions in the σ and pion channels with FRG. Our results confirm the softening of the particlehole mode in the σ channel near the QCD critical point, but not in the pion channel. In addition, we find that the low-momentum dispersion relation of the sigma-mesonic mode penetrates into the space-like region and the mode merges into the bump of the particle-hole mode.
This paper is organized as follows. In Sect. 2, we recapitulate the method [62,63] and describe how to calculate the spectral functions in the mesonic channels numerically with FRG. The numerical results are shown in Sect. 3. The phase diagram, the critical region, and the precise location of the critical point are presented in Sect. 3.2. In Sect. 3.4 the results of the spectral functions are shown, and the soft mode at the QCD critical point is discussed. Section 4 is devoted to the summary and outlook. In Appendix A, the detailed forms of the FRG flow equations after Matsubara summation are shown. In Appendix B, we derive the conditions for numerically stable calculation of the flow equation as a nonlinear partial differential equation.

Method
In this section, we recapitulate the method developed in Refs. [62,63] for calculating spectral functions in the quark-meson model with FRG, and present some details of our numerical procedure.

Procedure to derive spectral functions in meson channels
The FRG is based on the philosophy of the Wilsonian renormalization group and realizes the coarse graining by introducing a regulator function R k , which has the role of suppressing lower-momentum modes than the scale k for the respective field. In this method, the effective average action (EAA) k is introduced such that it becomes bare action S at a large UV scale k = and becomes the effective action at k → 0 with an appropriate choice of regulators. The flow equation for EAA, the Wetterich equation, can be derived as a functional differential equation [34]: where (n) k is the nth functional derivative of k with respect to fields. This equation has a one-loop structure and can be represented diagrammatically ( Fig. 1(a)).
In principle, one can get the effective action k=0 by solving Eq. (1) with the initial condition = S . However, it is prohibitively difficult to solve Eq. (1) in an exact way and some simplifications are introduced for practical use. One of the simplifications is to truncate the form of EAA. Since our purpose is to reveal the behavior of the low-momentum modes around the QCD critical point, we adopt a truncation in which the low-momentum fluctuations are properly taken into account.
We should now specify the low-energy effective model of QCD; we employ the two-flavor quarkmeson model as such a model. Then we take the local potential approximation (LPA) for the meson flow part as our truncation. This truncation corresponds to considering only the lowest order of derivative expansion for the meson flow part.
In the imaginary-time formalism, our truncated EAA is as follows [46]: where φ = (σ , π). The quark field ψ has the indices of a four-component spinor, color N c = 3, and flavor N f = 2. The cσ represents the effect of the current quark mass, which explicitly breaks N f = 2 chiral symmetry. In this truncation, we also neglect the flow of g s and the wave function renormalization. Therefore, only the meson effective potential U k has a k dependence. Although the truncated EAA itself gives only the simplest information about the momentum dependence of two-point Green's functions, the nonperturbative effects are to be incorporated through Eq. (1) with the truncated EAA used as the initial condition. More specifically, we first calculate the effective potential U k (φ 2 ) using Eq. (1). Then the chiral condensate σ 0 is obtained as σ satisfying the quantum equation of motion (EOM) δ k=0 /δσ = 0. In our case, this condition corresponds to obtaining σ that minimizes U k (σ 2 )−cσ , under the assumptions that the condensate is homogeneous The diagrammatic expression of this equation is shown in Fig. 1 (2) k,σ (P) and (2) k,π (P) in the sigma and pion channels, respectively, which are defined as where σ (P) and π a (P) are the Fourier transforms of σ (x) and π a (x), respectively, and P = (iω n , p) with ω n being the bosonic Matsubara frequency. The RHS of Eq. (3) contains k . In general, the flow equation comprises an infinite hierarchy of differential equations such that the flow equation for . We can simplify the flow equation for the n by replacing these derivatives with those derived from the truncated EAA. Under the above procedure, the integration of the flow equation for (2) k,σ (π) (P) leads to the two-point Green's function in the sigma (pion) channel with the nonperturbative effects incorporated.
A real-time two-point Green's function at finite temperature is obtained by analytic continuation for the temperature Green's function to real time, i.e., from imaginary Matsubara frequencies to real frequencies in the case of momentum representation. In our case, such an analytic continuation is successfully carried out at the level of the flow equation after the Matsubara summations in the RHS of Eq. (3), as follows: When one derives a retarded two-point Green's function via analytic continuation in the frequency ω plane, the analyticity of the Green's function in the upper half-plane of ω must be retained [67]. In the present case, the flow equation itself should be analytic in the upper half-plane after the analytic continuation. One can retain the analyticity in the upper half-plane easily by taking into account the following points. First, by choosing ω n -independent regulators, one can avoid possible extra poles in the ω plane in the flow equation otherwise arising from ω n dependence of the regulators. The second point is about the analytic continuation of thermal distribution functions n B,F (E + iω n ) obtained for a discrete (multiple of 2πT ) frequency ω n , where the subscript B, F stands for a boson or fermion, respectively, and E is ω n independent. Such factors appear in the flow equation after the Matsubara summation. Because of the periodicity of the exponential function, n B,F (E + iω n ) is equal to n B,F (E). However if n B,F (E + ω) is substituted for n B,F (E + iω n ), such a factor breaks the analyticity of the flow equation in the upper half-plane. Therefore n B,F (E + iω n ) should be replaced by n B,F (E) before the analytic continuation. By taking into account these points, the substitution ω + i for iω n with being a positive infinitesimal gives the flow equation for the inverse of the retarded Green's function (2),R k,σ ,π (ω, p). Finally, the spectral functions in the meson 5 channels are given in terms of the thus-obtained retarded Green's functions as follows:

Flow equations
In the present work, we adopt the 3D Litim's optimized regulators for bosons and fermions [68] as ω n -independent regulators: Then the insertion of Eq. (2) into Eq. (1) leads to the following flow equation for U k : where According to the procedure presented in the previous subsection, the flow equations for (2) k,σ (P) and (2) k,π (P) become: respectively, where π,π ∈ {π 1 , π 2 , π 3 } and π =π. The loop-functions J k,αβ (P), I k,α , and J where Q = (iq n , q) and The three-and four-point vertices some of which are expressed in terms of U k : Analytic continuation in Eq. (12) and Eq. (13) is carried out after the Matsubara summation in Eqs. (14)- (16). The detailed forms of Eqs. (14)- (16) after Matsubara summation are presented in Appendix A.
To solve these flow equations, we employ the following initial conditions at the UV scale k = :

Numerical procedure
We employ the grid method to solve Eq. (10) numerically. This method reveals the global structure of U k (σ 2 ) on discretized σ . We employ the fourth-order Runge-Kutta method to solve Eq. (10). It is known that certain conditions between the intervals of discretization of variables must be fulfilled to solve partial differential equations numerically in a stable way [69]. We have derived such conditions for Eq. (10), which are presented in Appendix B. We fix the intervals of discretization of σ and t in Eq. (10) according to these conditions. As stated before, the flow equation (10)   to k = 0 from k = to get the effective action k=0 in principle. However, due to the conditions for stable calculation mentioned above, solving the flow equation to small k is quite time-consuming for some regions of the (T , μ) plane, such as the low-temperature region of the hadronic phase. In such a region, the curvature of U k , i.e., m 2 σ , can take a negative value, which leads to small E σ for some σ . This gives large F and |G| (Eqs. (B13) and (B14)) as k decreases and the condition (B12) becomes difficult to satisfy at small k. Thus, some infrared scale k = k IR is introduced in practice, at which the numerical procedure is stopped. Of course, k IR should be as small as possible so that sufficiently low-momentum fluctuations are taken into account to describe the system around the critical point where vanishingly low-momentum excitations exist. Thus we choose a much lower value of k IR than the 40 MeV adopted in Ref. [63], and set k IR = 1 MeV as being small enough to incorporate the low-momentum fluctuations. Therefore our calculation will be reliable in the vicinity of the critical point, except for the small surrounding region where excitation modes with momentum scales lower than 1 MeV are strongly developed. Although , which appears after the analytic continuation, is defined as a positive infinitesimal, we set it to 1 MeV in the present calculation, which should be small enough for present purposes.
3D momentum integrals remain after the Matsubara summation in Eq. (12) and Eq. (13). However, the integrals can be fully calculated analytically for zero external momentum. Even for a finite external momentum, they can be nicely reduced to 1D integrals, which are evaluated numerically. The numerical integrations involve a tricky point, and one has to take care of the poles of each term in the integrands. For example, the first and third terms of the integrand of the second integral in Eq. (A2) have the same poleẼ α = E α + ip 0 . If such terms are integrated separately, a large cancellation can occur, which then leads to big numerical errors. Therefore, we first combine such terms analytically in the integrand before numerical integrations.

Parameter setting
The truncated EAA Eq. (2) and the initial condition Eq. (25) have some parameters that are fixed so as to reproduce the observables in vacuum: We use the same values for the parameters as those in Ref. [63] and list them in Table 1.

Phase diagram and quark-number susceptibility
We show the phase diagram in Fig. 2, where a contour map of the chiral condensate is also given. One sees that chiral restoration occurs as the temperature is raised, and the phase transition is not a genuine  one but a crossover, except for the low-temperature and large chemical potential region, where the phase transition is of first order. This feature is qualitatively in accordance with the results given in the literature, although the location of the critical point here is in a somewhat smaller temperature region than that given in Ref. [63]. The detailed procedure for locating the critical point is described below.
At the QCD critical point, the chiral susceptibility diverges. Therefore, we locate the critical point by searching for the point where the sigma screening mass M σ , the square of which is the inverse of the chiral susceptibility, becomes the smallest: We seek the minimum position of M σ using the data points where M σ is greater than 1 MeV, because our choice of k IR = 1 MeV enables us to take into account fluctuations whose momentum scales are greater than k IR so as to make the result of M σ reliable when M σ is larger than 1 MeV. We also identify the first-order phase transition by a discontinuity of the chiral condensate. The results at T = 5. We find that the sigma screening mass becomes smallest between T = 5.0 MeV and T = 5.2 MeV and between μ = μ t (5.0 MeV) and μ = μ t (5.2 MeV). Therefore, the critical temperature T c and the critical chemical potential μ c are estimated as T c = 5.1 ± 0.1 MeV and μ c = 286.6 ± 0.2 MeV. The position of the critical point is quite different from the (T , μ) = (10 MeV, 292.97 MeV) given in Ref. [63]. Such a difference may be attributed to the different choice of k IR . In the following discussion, we regard T c and μ c as 5.1 MeV and 286.686 MeV, respectively. As seen in the behavior of the chiral condensate shown in the right panel of Fig. 3, the phase transition along the chemical potential is of first order when T = 5.0 MeV and a crossover when T = 5.2 MeV.
It is known [55] that the quark-number susceptibility χ q ≡ ∂ρ q /∂μ is coupled to the scalar susceptibility at finite μ and can be used to reveal the nature of the phase transition, where ρ q denotes the quark-number density. Indeed χ q shows a singular behavior at the critical point, and hence an enhancement of χ q is a useful measure of the critical region. We thus calculate χ q to map out the critical region, as was done in Refs. [23,56].   assumption that the chiral condensate is homogeneous and π = 0, which in our case reads The quark-number susceptibility χ q is given by differentiating Eq. (29) twice with respect to the chemical potential: We carry out the derivatives numerically. Figure 4 shows the contour map of the quark-number susceptibility, normalized by the value for the massless free quark gas:

Spectral function in the σ channel away from the critical point
Before entering into discussions on the spectral properties in the meson channel near the critical point, we first show the numerical result of the spectral function ρ σ (ω, p) in the sigma channel away from the critical point in the hadronic and QGP phases so that the peculiar behavior of the spectral functions near the critical point shown in the next subsection shall be prominent. In the former case, there is a sharp peak at ω = 290 MeV and a relatively small bump in the space-like region ω < p; these correspond to the sigma meson with a modified mass at finite temperature and the phonon mode composed of particle-hole excitations, respectively, which is in accord with the result in the RPA in Ref. [49].
The spectral function also tells us the decay and absorption processes of the particle excitations from the width of the corresponding peaks or bumps. In our energy scale, the following processes contribute to the spectral function ρ σ (ω, p): where σ * denotes a virtual state in the sigma channel with energy-momentum (ω, p). The energymomentum conservation gives constraints on the possible (ω, p) region for the first three processes as follows:  which are all in the time-like region. On the other hand, the second three processes are all collisional ones and possible only in the space-like region, 0 ≤ ω < p. In particular, the last process σ * ψ → ψ corresponds to the absorption process of the σ * mode into a thermally excited quark. In short, the width of the large bump at ω = 290 MeV in Fig. 5(a) comes from the 2π decay process, while the small bump arises from the space-like processes.
In the latter case, at T = 50 MeV and μ = 400 MeV, the peak position corresponding to the sigma meson is shifted to ω = 210 MeV, while the bump of the particle-hole excitations still persists in the space-like region.

Spectral functions near the QCD critical point
We calculate the spectral function in the σ channel near the QCD critical point, by increasing the chemical potential toward μ c along a constant temperature line T = T c . The results at μ = 286.00 MeV, μ = 286.50 MeV, and μ = 286.57 MeV are shown in Fig. 6(a). One can see the sigma-mesonic peak as well as bumps corresponding to 2σ and 2π decay in the time-like region.  The peak position of the sigma-mesonic mode shifts to lower energy as the system approaches the critical point. The position of the 2σ threshold also shifts to a lower energy while those of the 2π and ψψ thresholds hardly change. The spectral function in the space-like region is drastically enhanced as the system is close to the critical point. This behavior can be interpreted as the softening of the particle-hole mode, which is in accordance with the result in Ref. [49]. In Fig. 6(b), we show the results at chemical potentials much closer to the critical point. Because of numerical instability in 286.60 MeV ≤ μ 360 MeV, we choose μ = 286.58 MeV and μ = 286.59 MeV. For comparison, the result at μ = 286.57 MeV is also shown. These results are drastically different from those in μ ≤ 286.57 MeV. In μ > 286. 57 MeV, the peak of the sigma-mesonic mode penetrates into the space-like region and then merges into the particle-hole mode.
We can see the dispersion relations of the modes by making contour maps of the spectral functions as functions of ω and p. Figure 7 shows the dispersion relations of the sigma-meson and particlehole modes near the critical point. At μ = 286.3 MeV, the sigma-mesonic peaks can be seen in the time-like region as well as the particle-hole bump in the space-like region. As the chemical potential increases, the dispersion relation of the sigma-mesonic mode shifts downward and it touches the light cone near μ = 286.575 MeV. At μ = 286.59 MeV, in the low-momentum region the sigma-mesonic mode clearly penetrates into the space-like region and merges with the particle-hole bump, which has a flat dispersion relation in the low-momentum region. Our results indicate that the sigma-mesonic mode as well as the particle-hole mode can become soft near the critical point.
One of the possible triggers of this phenomenon is the level repulsion between the sigma-mesonic mode and other modes. In particular, the two-sigma (σ σ ) mode is considered to play an important role in the level repulsion since the threshold of the two-sigma mode shifts downward as the system approaches the critical point. Let us suppose that the particle-hole mode, the sigma-mesonic mode, and the two-sigma mode can each be described by a state having a single energy level. Then the system can be regarded as a three-level system, as depicted in Fig. 8(a): The interaction within the three states leads to a level repulsion: If the interaction between the σ -mesonic mode and the σ σ state becomes sufficiently strong as the system approaches the critical point, the energy level of the sigma meson will be so strongly pushed down that it penetrates into the space-like region. To show that this scenario can be the case, we change the strength of the three-point vertex (0, 3) k,σ σ σ by hand PTEP 2016, 073D01 T. Yokota et al.  to investigate the behavior of the sigma-meson peak. The results in the cases of multiplying (0,3) k,σ σ σ by factors 0.8 and 1.02 are shown in Fig. 8(b). The position of the sigma meson goes up when the three-point vertex is weakened, whereas it exhibits a downward shift to a lower energy when the three-point vertex is slightly enhanced. This result suggests that the above interpretation in terms of a level repulsion can be correct.
Here it should be noted that our results exhibit a superluminal group velocity of the sigma-mesonic mode near the critical point, as seen in Fig. 7 for p = 100 MeV at μ = 286.59 MeV. Such an unphysical extreme behavior may be an artifact of our truncation scheme, in which some of the higher-order terms in the derivative expansion, such as the wave-function renormalization and so on, are neglected, although a drastic softening of the sigma-mesonic mode may be true. Conversely speaking, such a drawback could disappear if one uses improved methods with higher-derivative terms being incorporated. One of the most important improvements is the inclusion of wave-function renormalization [70,71], since this may become important when additional modes emerge. However, this task is beyond the scope of the present work and will be left as a future project. 14 So far, we have concentrated on the spectral function in the sigma channel and seen interesting behaviors of it near the critical point. It would be intriguing to examine whether the spectral function ρ π (ω, p) in the pion channel shows any peculiar behavior near the critical point. The numerical result of ρ π (ω, p) near the critical point is shown in Fig. 9, from which one can clearly see the dispersion relation of the pion mode in the time-like region but not in the space-like region. In contrast to ρ σ , ρ π hardly changes near the critical point, indicating that there is no critical behavior in the isovector pseudoscalar modes in either the space-like or time-like region. This different critical behavior in the sigma and pion channels may be understood as follows: First of all, the finite current quark mass makes the would-be chiral transition cease to be of second order, and hence prevents the quartet mesons composed of the sigma meson and pion from becoming soft modes inherent in the second-order transition, as mentioned in Sect. 1. However, the finite chemical potential μ = 0 makes scalar-vector coupling possible and the critical point can get to exist with the Z 2 universality class of the second-order transition. Thus the soft modes inherent for the second-order transition appear due to the very scalar-vector coupling in the space-like region, which is primarily composed of particlehole excitations. In principle, the pion is coupled to fluctuations in the isovector pseudoscalar or axial vector channels, which are reduced to spin-isospin modes in the nonrelativistic limit [72]. Our result simply shows that such pionic modes do not develop in the space-like region, at least around the critical point. In the scalar channel, the coupling of the single and double sigma modes is so strong that a level repulsion between them drastically lowers the energy of the sigma-mesonic mode in the time-like region. In contrast to the scalar channel, the coupling of the single pion mode to other modes, such as the mode that consists of one sigma and one pion, would not be strong enough to cause the softening behavior of the isovector pseudoscalar mode in the time-like region.

Summary
We have calculated the spectral functions in the meson channels in the quark-meson model with the functional renormalization group method based on the local potential approximation (LPA). A particular emphasis is put on the behavior of the spectral function in the σ channel near the QCD critical point. Our results show that the particle-hole mode (phonon) is enhanced near the critical point, and thus imply that the density fluctuations are soft modes at the critical point, as was suggested in the RPA using the NJL model [49]. In addition, we have found that the low-energy dispersion curve of the sigma meson penetrates into the space-like region and the mode merges into the particle-hole mode near the critical point. This result may imply that the sigma meson also acts as a soft mode at the QCD critical point. We have also suggested that a possible level repulsion between the sigma meson and the two-sigma state leads to the anomalous softening of the sigma meson: An artificial variation of the strength of the σ three-point vertex (0, 3) k,σ σ σ strongly affects the position of the sigma meson. We have also investigated the spectral function in the pion channel near the critical point, which shows no softening in either the time-like region or the space-like region, in contrast to the isoscalar modes.
Since our result might provide a new picture in which the critical dynamics at the QCD critical point can be described by an effective theory composed of not only the hydrodynamical modes including the density fluctuation but also the sigma-meson mode, there might be some implications for the dynamical class of the QCD critical point [21,50,73].
In recent years, the possible existence of inhomogeneous chiral phases in dense quark matter has been intensively examined [74][75][76][77]. It would be interesting to investigate in the FRG how the existence 15 of the inhomogeneous phases affects the behavior of low-energy modes, including the Nambu-Goldstone modes. It is also worth emphasizing that our analysis of the spectral functions in the space-like region or particle-hole modes can be extended to that of precursors of such inhomogeneous phases: our results showing the softening and nonsoftening of the particle-hole modes in the sigma and pion channels might imply that the inhomogeneous phase with pion condensate does not come into existence as a result of a second-order transition.
Our calculation is based on the local potential approximation, in which some of the higher-order terms in the derivative expansion are neglected. Such a simple truncation scheme might be an origin of a superluminal group velocity in the close vicinity of the critical point encountered in Sect. 3. Thus, it is imperative to confirm the results by employing improved methods incorporating higherderivative terms, including the wave-function renormalization, since it is important to identify when composite collective modes emerge. This intriguing task will, however, be left for future work.
for which the external Matsubara frequency ip 0 is to be replaced to make analytic continuation.

B. Numerical stability conditions for solving the flow equation of the effective potential
In general, when one solves a partial differential equation numerically, the discretization of derivatives may cause numerical instability. Thus, one needs to impose numerical stability conditions to avoid the enhancement of the error due to accumulation. The derivation of such conditions is concretely demonstrated in the case of linear partial differential equations and briefly mentioned in the case of nonlinear partial differential equations in Ref. [69]. In this appendix we consider the numerical stability conditions for a nonlinear partial differential equation that is a generalized equation of Eq.
(10) to derive the numerical stability conditions for solving Eq. (10) in the grid method. Consider a nonlinear partial differential equation for some function u(t, σ ) of the following form: where f is an arbitrary real function. This equation is a generalized equation of Eq. (10). We derive the numerical stability conditions for this equation in the case of the forward difference for the t-derivative and the central three-point difference for the σ -derivative. We expect that the derived 18 where X ≡ cos(k x) and h(X ) is defined as Considering that h(X ) satisfies h(1) = 1, one can easily understand that the condition Eq. (B8) is equivalent to the following conditions: These conditions are rewritten as follows: We show later that G < 0 and t < 0 in the case of Eq. (10). Supposing that G < 0 and t < 0, Eq. (B11) is finally rewritten as In the case of Eq. (10), t and u(t, σ ) are identified with ln(k/ ) and U k (σ 2 ), respectively, and F and G are derived to be: (B14) G is negative definite, and t is also negative because the direction of flow is from k = to k = 0. Although the above discussion contains rough approximations, numerical instability can be avoided by adjusting t and σ so that Eq. (B12) is satisfied. Because the above condition of Eq. (B12) is too strict when σ is close to zero, we neglect the condition around σ = 0.