Critical behaviors of the $O(4)$ and $Z(2)$ symmetries in the QCD phase diagram

In this work we have studied the QCD phase structure and critical dynamics related to the 3-$d$ $O(4)$ and $Z(2)$ symmetry universality classes in the two-flavor quark-meson low energy effective theory within the functional renormalization group approach. We have employed the expansion of Chebyshev polynomials to solve the flow equation for the order-parameter potential. The chiral phase transition line of $O(4)$ symmetry in the chiral limit, and the $Z(2)$ line of critical end points related to the explicit chiral symmetry breaking are depicted in the phase diagram. Various critical exponents related to the order parameter, chiral susceptibilities and correlation lengths have been calculated for the 3-$d$ $O(4)$ and $Z(2)$ universality classes in the phase diagram, respectively. We find that the critical exponents obtained in the computation, where a field-dependent mesonic nontrivial dispersion relation is taken into account, are in quantitative agreement with results from other approaches, e.g., the conformal bootstrap, Monte Carlo simulations and $d=3$ perturbation expansion, etc. Moreover, the size of the critical regime in the QCD phase diagram is found to be very small.


I. INTRODUCTION
Significant progress has been made in studies of QCD phase structure over the last decade, both from the experimental and theoretical sides; see, e.g. [1][2][3][4][5][6][7][8][9][10]. One of the most prominent features of the QCD phase structure is the probable presence of a second order critical end point (CEP) in the phase diagram spanned by the temperature T and baryon chemical potential µ B or densities, which separates the first order phase transition at high µ B from the continuous crossover at low µ B [1]. The existence and location of CEP are, however, still open questions, whose answers would definitely help us to unravel the most mysterious veil related to the properties of strongly interacting matter under extreme conditions. The Beam Energy Scan (BES) Program at the Relativistic Heavy Ion Collider (RHIC) is aimed at searching for and locating the critical end point, where fluctuation observables sensitive to the critical dynamics, e.g., high-order cumulants of net-proton, net-charge, net-kaon multiplicity distributions, have been measured [11][12][13][14]. Notably, a non-monotonic dependence of the kurtosis of the netproton multiplicity distribution on the beam energy with 3.1σ significance in central collisions has been reported by the STAR collaboration recently [15].
On the other hand, lattice QCD simulations have provided us with a plethora of knowledge about the QCD phase structure, e.g., the crossover nature of the chiral phase transition at finite T and vanishing µ B with physical current quark mass [16], pseudo-critical temperature [17,18], curvature of the phase boundary [9,19], etc. Because of the notorious sign problem at finite chemical potential, the reliability regime of lattice calculations is restricted to be µ B /T 2 ∼ 3, where no CEP has been found. Free from the sign problem, the firstprinciple functional approaches, e.g, the functional renor- * wjfu@dlut.edu.cn malization group (fRG) and Dyson-Schwinger equations (DSE), could potentially extend the regime of reliability to µ B /T ∼ 4 [5,7]. With benchmark tests of observables at finite T and low µ B in comparison to lattice calculations, e.g., the quark condensate, curvature of the phase boundary, etc., functional approaches, both fRG and DSE, have predicted a CEP located in a region of 450 MeV µ B 650 MeV [5,7,[20][21][22] recently.
An alternative method used to circumvent the possible location of CEP, is to determine the critical temperature T c of the chiral phase transition in the chiral limit, more specifically, i.e., massless light up and down quarks and a physical strange quark mass. Since it is believed that the value of T c sets an upper bound for the temperature of CEP [23,24]. Very recently, the critical temperature T c in the chiral limit has been investigated and its value is extrapolated from both lattice simulations [25] and functional approach [26]. Moreover, further lattice calculations indicate that axial anomaly remains manifested at T ≈ 1.6 T c , which implies that the chiral phase transition of QCD in the chiral limit is of 3-d O(4) universality class [27]; see, e.g., [28] for more discussions about the relation between the axial anomaly and the symmetry universality classes.
In this work, we would like to study the QCD phase structure in the chiral limit and finite current quark mass, i.e., with a finite pion mass, in the two-flavor quarkmeson low energy effective theory (LEFT) within the fRG approach. For more discussions about the fRG approach, see, e.g., QCD related reviews [29][30][31][32][33][34][35][36]. In contrast with the lattice simulation and the first-principle fRG-QCD calculation [25,26], the chiral limit could be accessed strictly in the LEFT. Furthermore, we would also like to study the critical behaviors of the 3-d O (4) and Z(2) universality classes, including various critical exponents, which belong to the second-order chiral phase transitions in the chiral limit and at the critical end point with finite quark mass, respectively. To that end, we expand the effective potential of order parameter as a sum of Chebyshev polynomials in the computation of fRG arXiv:2101.08484v1 [hep-ph] 21 Jan 2021 flow equations; see [37] for more details. The Chebyshev expansion of solutions to a set of integrodifferential equations is, in fact, a specific formalism of more generic pseudo-spectral methods [38], and see also, e.g., [39][40][41] for applications of pseudo-spectral methods in the fRG.
In fact, another two numerical methods are more commonly used in solving the flow equation for the effective potential: one is the Taylor expansion of the effective potential around some value [42,43], and the other discretization of the effective potential on a grid [44]. The (dis)advantages of these two methods are distinct. The former is liable to implementation of the numerical calculations, but short of global properties of the effective potential, that is, however, indispensable to studies of chiral phase transition in the chiral limit or around CEP; the latter is encoded with global information on the potential, but it loses numerical accuracy near the phase transition point which is necessary especially for the computation of critical exponents. The Chebyshev expansion used in this work combines the merits from both approaches, i.e., the global potential and the numerical accuracy, and thus it is very suitable for the studies of critical behaviors in the QCD phase diagram. Remarkably, a discontinuous Galerkin scheme has been applied in the context of fRG recently [45], which is well-suited for studies of the first-order phase transition. This paper is organized as follows: In Sec. II we briefly introduce the flow equations in the quark-meson LEFT and the method of the Chebyshev expansion for the effective potential. The obtained phase diagram and QCD phase structure are presented and discussed in Sec. III. In Sec. IV scaling analyses for the the 3-d O(4) and Z(2) universality classes are performed, and various critical exponents are obtained. We also discuss the size of the critical regime there. In Sec. V we give a summary and conclusion. Some threshold functions and anomalous dimension in the flow equations, and some relations for the Chebyshev polynomials are collected in Appendix A and Appendix B, respectively.

II. FUNCTIONAL RENORMALIZATION GROUP AND THE LOW ENERGY EFFECTIVE THEORIES
Thanks to the Wilson's idea of the renormalization group (RG), see, e.g., [46], it has been well known that usually the active degrees of freedom are quite different, when the energy scale of a system evolves from a hierarchy into another. The relevant dynamics in different hierarchies are connected with each other through the evolution of RG equations. To be more specific, in QCD the partonic degrees of freedom, i.e., the quarks and gluons, in the high energy perturbative regime are transformed into the collective hadronic ones in the nonperturbative region of low energy, with the RG scale evolving from the ultraviolet (UV) to infrared (IR) limits [47], and see also, e.g., [7,30,[48][49][50][51][52][53][54] for recent development of the relevant ideas within the fRG approach. When the momentum or RG scale is below, say ∼ 1 GeV, which is related to a narrow transition region from the perturbative to nonperturbative QCD, calculated results of Yang-Mills theory and QCD in Landau gauge indicate that the gluons develop a finite mass gap and decouple from the system, and see, e.g. [7,52,55,56] for more details. As a consequence, contributions to the flow equations of effective action from the glue sector could be safely neglected, if the initial evolution scale is set at a UV scale Λ 1 GeV.
Hence, within the fRG approach, one is left with the flow equation for the low energy effective theory, which reads with the RG scale k and the RG time defined as t = ln(k/Λ). Apparently, Eq. (1) is an ordinary differential equation for the k-dependent effective action, Γ k [Φ], the arguments Φ = (q,q, φ) of which are the quark and mesonic fields in the LEFT. The equation in Eq. (1), which describes the evolution of the effective action with the RG scale, is also well known as the Wetterich equation [57], see also [58,59]. The flow receives contributions from both the quark and mesonic degrees of freedom, as shown on the r.h.s. of Eq. (1), where G qq,k and G φφ,k are the k-dependent full quark and meson propagators, respectively, and are related to the quadratic derivatives of Γ k [Φ] with respect to their respective fields, viz.
In this work, we adopt a truncation for the effective action in Eq. (1) as follows with the shorthand notation x = 1/T 0 dx 0 d 3 x, where the quark field q = (u , d) T and the meson field φ = (σ, π) are in the fundamental and adjoint representations of SU (N f ) in the flavor space with N f = 2, respectively. They interact with each other via a Yukawa coupling with a coupling strength h Y,k , where the subscript Y is used to distinguish it from the reduced external field h in Eq. (15). Here T i (i = 1 , 2 , 3) are the generators of SU (2) with Tr(T i T j ) = 1 2 δ ij and T 0 = Note that both the effective potential V k (ρ) and the mesonic wave function renormalization Z φ,k (ρ) in Eq. (3) depend on the meson field by means of ρ = φ 2 /2, which are O(4) invariant. Z q,k is the quark wave function renormalization. Notice that the term linear in the order parameter field, i.e., −cσ in Eq. (3), breaks the chiral symmetry explicitly, and thus here c is essentially an external "magnetic" field in the language of magnetization. Moreover,μ = diag(µ u , µ d ) is the matrix of quark chemical potentials in the flavor space, and µ = µ u = µ d is assumed throughout this work, which is related to the baryon chemical potential via µ = µ B /3. For more discussions about the quark-meson LEFT in Eq. (3) or its extensions, e.g., Polyakov-loop quark-meson LEFT, QCD assisted LEFT, etc., and their applications in calculations of QCD thermodynamics and phase structure, fluctuations and correlations of conserved charges, etc., see, e.g., [10,43,44,[64][65][66][67][68][69][70][71][72][73][74][75][76][77][78][79][80][81][82].

A. Flow equations
Substituting the effective action in Eq. (3) into the Wetterich equation in Eq. (1), one readily obtains the flow equation of the effective potential as follows with the threshold functions l which are RG invariant and dimensionless. The meson and quark anomalous dimensions in the threshold functions in Eq. (4) are defined as follows where the meson anomalous dimension is obtained by projecting the flow equation in Eq. pion propagator, to wit, the explicit expression of which is presented in Eq. (A9). Note that η φ,k is dependent on the meson field via ρ.
In comparison to the effects of the meson wave function renormalization on the chiral phase transition at finite temperature and density, it has been found that those of quark wave function renormalization and the running Yukawa coupling are relatively milder, see, e.g., [42,43,75]. Therefore, in this work we adopt the simplification as follows with the renormalized Yukawa coupling given in Eq. (A10), and use two different truncations: one is the usual local potential approximation (LPA), where the mesonic anomalous dimension is vanishing as well, and the k-dependent term in Eq. (3) is just the effective potential; the other is the truncation with the fielddependent mesonic anomalous dimension in Eq. (8) taken into account besides the potential, which is denoted as LPA in this work. Note that the notation LPA in literatures, e.g., [75,83], usually stands for the truncation with a field-independent mesonic anomalous dimension which is, strictly speaking, different from the case in this work.
As an illustrative example, we show the mesonic wave function renormalization Z φ ≡ Z φ,k=kIR as a function of the renormalized sigma fieldσ = Z 1/2 φ σ obtained in LPA in Fig. 1, where k IR is the RG scale in the IR limit, and one would has k IR → 0 in principle, which, however, is impossible to realize in numerical calculations. In our calculation the value of k IR is reduced as small as possible, and we find the convergence is obtained when k IR = 1 MeV. Note that the mesonic wave function renormalization at the scale of UV cutoff Λ, see Sec. III in the following, is assumed to be identical to unity, i.e., Z φ,k=Λ = 1. In Fig. 1, we choose several values of temperature T = ∆T + T c at and above the critical temperature that is T c = 143.6 MeV in the chiral limit and at vanishing µ B . One observes that with the increase of the temperature, the peak structure of Z φ as a function of the renormalized sigma fieldσ becomes smoother.

B. Chebyshev expansion of the effective potential
In this work we solve the flow equation in Eq. (4) by expanding the effective potential as a sum of Chebyshev polynomials up to an order N v , to wit, where quantities with a bar denote renormalized variables. The Chebyshev polynomial T n (ρ) is given in Eq. (B6), and the superscript [0,ρ max ] in Eq. (B6) denoting the interval ofρ is neglected for brevity here. Differentiating Eq. (10) with respect to the RG time t with ρ fixed, one is led to where we have used the Chebyshev expansion for the derivative of the effective potential as shown in Eq. (B10) and coefficients d n,k 's are the respective expanding coefficients. Employing the discrete orthogonality relation in Eq. (B4) by summing up the N + 1 zeros of T N +1 (ρ) in Eq. (B7), one arrives at which is the flow equation for the expansion coefficients in Eq. (10).

III. PHASE DIAGRAM
It is left to specify the parameters in the LEFT, prior to presenting our calculated results. The UV cutoff of flow equations in the LEFT is chosen to be Λ = 500 MeV, and the effective potential in Eq. (3) at k = Λ reads with λ Λ = 20 and ν Λ = 0 . The Yukawa coupling is k-independent as shown in Eq. (9) and is given byh y = 6.4. Concerning the Chebyshev expansion, we choose N = 81 for the number of zeros and N v = 21 for the maximal order of Chebyshev polynomials. We have also checked that there is no difference when the value of N v is increased. Moreover, the upper bound ofρ is chosen to bē ρ max = 9 × 10 3 MeV 2 , well above the value of minimum of the potential in the IR. In the LPA, these values of parameters lead to the pion decay constant f π =87 MeV and the constituent quark mass m q =278.4 MeV in the vacuum and in the chiral limit. While if the explicit breaking strength of the chiral symmetry in Eq. (3) is increased up to c = 1.85 × 10 −3 (GeV) 3 , one obtains the physical pion mass m π =138 MeV, as well as f π =93 MeV and m q =297.6 MeV in the vacuum. Note that in order to facilitate the comparison between the calculation with the truncation LPA and that with LPA , we use the same values of parameters above in the LPA computation as in LPA.
In Fig. 2 we show the phase diagrams of LEFT in the T −µ B plane, calculated within the fRG approach with the truncations LPA and LPA , in the left and right panels, respectively. The black dashed lines in both panels denote the second-order O(4) chiral phase transition of N f = 2 flavor in the chiral limit. The black circles indicate the location of the tricritical point, beyond which the second-order phase transition evolves into a discontinuous first-order one, which are shown by the solid lines. Note that the solid lines of different colors in the left panel denote the first-order phase transitions with different pion masses in the vacuum, i.e. different values of c in Eq. (3), and in the right panel, we only give the first-order phase transition line in the chiral limit, since numerical calculations become quite difficult in the region of high µ B and low T with the truncation LPA . The red dashed lines in both panels are the trajectories of the critical end points with the change of the strength of explicit chiral symmetry breaking c, which belong to the 3-d Z(2) Ising university class.
The critical temperature at vanishing baryon chemical potential is found to be T c = 144 MeV in LPA and 143 MeV in LPA in the chiral limit.

IV. CRITICAL BEHAVIOR AND CRITICAL EXPONENTS
A variety of scaling analysis has been performed for the O(4) universality class, e.g., in the O(N ) model [84][85][86][87][88][89] and two-flavor quark-meson model [90][91][92][93]. The dynamics of a system in the critical regime near a secondorder critical point is governed by long-wavelength fluctuations, and the correlation length tends to be divergent as the system moves towards the critical point. Critical exponents play a pivotal role in studies of the critical dynamics, which are independent of micro interactions, but rather universal for the same symmetry class, dimension of the system, etc., and see [88,93] for more details. In the following, we follow the standard procedure and give our notations for the relevant various critical exponents.
To begin with, from the effective action in Eq. (3) one readily obtains the thermodynamic potential density, which reads where the order parameter field σ ≡ σ or ρ = σ 2 /2 is on its equation of motion. We then introduce the reduced temperature and reduced external "magnetic" field as follows where T c is the critical temperature, and they are normalized by T 0 and c 0 , i.e., some appropriate values of T and c. In the language of magnetization under an external magnetic field, the order parameter σ here is just the corresponding magnetization density, i.e., M ≡ σ, and the explicit chiral symmetry breaking parameter c is equivalent to the magnetic field strength H ≡ c. We will not distinguish them in the following any more. In the critical regime the thermodynamic potential in Eq. (14) is dominated by its singular part f s , i.e., where the second term on the r.h.s. is the regular one, and the notation for the baryon chemical potential is suppressed. In what follows we adopt the notations in [94], and the scaling function f s (t, h) on the r.h.s. of Eq. (16) satisfies the scale relation to leading order, viz.
where is a dimensionless rescaling factor. The scaling function in Eq. (17) leads us to a variety of relations for various critical exponents [88,90,91,95], e.g., with the spacial dimension d. The critical exponents β and δ describe the critical behavior of the order parameter in the direction of t or h, respectively, i.e., The exponent γ is related to the susceptibility of order parameter χ, and ν to the correlation length ξ, which reads The scaling relation in Eq. (17) allows us to readily obtain the critical behavior for various observables. For instance, the order parameter and its susceptibilities read where χ σ ≡ χ l and χ π ≡ χ t are also called as the longitudinal and transverse susceptibilities, respectively.
Choosing an appropriate value of the rescaling factor such that h y h = 1 in Eq. (17), one is led to with the scaling variable z = t/h 1/(βδ) . Inserting Eq. (23) into the first equation in Eq. (22), one arrives at where we have introduced which is a scaling function dependent only on z. With appropriate values of H 0 and T 0 in Eq. (15), it can be shown that the scaling function in Eq. (25) has the properties f (0) = 1 and f (z) (−z) β with z → −∞ [94]. Consequently, it is straightforward to express the longitudinal and transverse susceptibilities in Eq. (22) in terms of the scaling function f (z), to wit, with and Alternative to the choice of h y h = 1 in Eq. (17), one can also employ t yt = 1, which is equivalent to the Widom-Griffiths parametrization [96,97] of the equation of state by means of the scaling variables, as follows which are obviously related to the other parametrization by the relations which read Hence the scaling function y(x) has the properties y(0) = 1 and y(−1) = 0. In the same way, one readily obtains the expressions of susceptibilities in this parametrization, which read A. Order parameter The flow equation of effective potential in Eq. (4) is solved by the use of the Chebyshev expansion as discussed in Sec. II B, i.e., evolving the flow equations of the expansion coefficients in Eq. (12) from the UV cutoff Λ to the infrared limit k → 0, and then the expectation value of the order parameter σ is determined by minimizing the thermodynamic potential in Eq. (14). Note that two different truncations, i.e., LPA and LPA as shown in Sec. II A, are employed in the calculations.
The critical exponents β and δ are given in Eqs. (19) and (20), which are related to the scaling behavior of the order parameter as the phase transition is approached towards in the temperature or external field direction, respectively. Note, however, that in the case of Z(2) phase transition as indicated by the blue cross in the phase diagram in Fig. 2, the order parameter should be modified slightly and we introduce the reduced order parameter which readsσ where f π is the pion decay constant in the vacuum and σ is the expectation value of sigma field at the phase transition point, which is nonvanishing on the red dashed lines of Z(2) in the phase diagrams in Fig. 2. Correspondingly, the reduced external field in Eq. (15) is modified into where c is the σ -related external field on the Z(2) phase transition line. Notice that both c and σ are vanishing on the O(4) phase transition line, viz., the black dashed lines in Fig. 2. In our calculations below, the normalized external field strength c 0 in Eq. (34) is chosen to be the value corresponding to the physical pion mass, and the normalized temperature in Eq. (15) is to be the critical one T 0 = T c . In Fig. 3 we show the log-log plots of the reduced order parameterσ versus the reduced temperature −t or external field h for the second-order O(4) and Z(2) phase transitions. The calculations are performed in the quarkmeson LEFT with the fRG in both LPA and LPA . The phase transition points are chosen to be the locations of the red and blue crosses in the phase diagrams in Fig. 2 for the O(4) and Z(2) universality classes, respectively. A linear relation is used to fit the calculated discrete data points in Fig. 3, and as shown in Eq. (19) and Eq.
In the same way, the values of δ are obtained as follows It is found that the critical exponents β and δ of the O(4) and Z(2) phase transitions in 3-d systems calculated in this work are consistent with previous results, e.g., Monte Carlo simulation of spin model [98] and d = 3 expansion for Z(2) [99]. Comparing the relevant results in LPA and LPA , one observes that both β and δ obtained in LPA are slightly smaller than those in LPA.

B. Preliminary assessment of the size of the critical region
It is well known that critical exponents are universal for the same universality classes. The size of the critical region is, however, non-universal and depends on the interactions and other details of system concerned. Furthermore, there has been a longstanding debate on the size of the critical region in QCD. Lattice QCD simulations show that the chiral condensate, i.e., the order parameter in Eq. (24), for physical quark masses are well described by Eq. (24) plus a small analytic regular term [25,100,101], which, in another word, implies that the size of the critical regime of QCD is large enough, such that QCD with physical quark mass is still in the chiral critical regime. On the contrary, it is found in [88,94,102] that the pion mass required to observe the scaling behavior is very small, at least one order of magnitude smaller than the physical pion mass. Moreover, it is also found that the critical region around the CEP in the QCD phase diagram is very small [103]. In Tab. I we present the values of the critical exponent β extracted from different ranges of temperature. One observes that when the temperature range is away from the critical temperature larger than 0.01 MeV, the value of β deviates from its universal value pronouncedly. This applies for both the O(4) and Z(2) universality classes. Given the systematic errors in the computation of this work, one could safely conclude that our calculation indicates that the critical region in the QCD phase diagram is probably very small, and it is smaller than 1 MeV in the direction of temperature.

C. Chiral susceptibility
According to Eq. (29), the reduced order parameter readsσ Moreover, it has been shown in [97] that given x > 0 and M > M 0 for some value M 0 in Eq. (29), the scaling function can be expanded as Inserting the leading term in Eq. (40) into Eq. (39) and utilizing the relation γ = β(δ − 1) as shown in Eq. (18), one is led to the reduced order parameter with t > 0 and h → 0, which readsσ Consequently readily obtained as follows which is in agreement with Eq. (21) in the limit h → 0 and in the symmetric phase, as it should be. Equation (41) also allows us to extract the value of the exponent γ, by directly investigating the scaling relation ofσ and t in the chiral symmetric phase with a fixed, small value of h. In Fig. 4 we show the logarithm of the longitudinal susceptibility χ σ versus that of the reduced temperature, where h = 3.5 × 10 −10 is chosen in the calculations. We have checked that this value of h is small enough to make sure that the value of γ obtained from the linear fit of ln(χ σ )-ln(t) is convergent. In the same way, the flow equations of fRG are resolved with two truncations LPA and LPA , and the phase transition points are chosen to be at the locations of the red and blue crosses in the phase diagrams in Fig. 2 Once more, one observes that these values, in particular those obtained in the LPA , are in good agreement with the values of γ for the O(4) and Z(2) symmetry universality classes, respectively; see, e.g., [98,99]. Calculations are done within the fRG approach with the truncations LPA and LPA . The phase transition is chosen to be near the location of the red cross in the phase diagrams in Fig. 2 for the O(4) symmetry universality class.
In Fig. 5 the longitudinal susceptibility of the order parameter χ σ , as shown in Eq. (22), is depicted versus the reduced temperature with several different values of the reduced external field. Here we only focus on the case of O(4) symmetry, and thus choose the phase transition to be near the location of the red cross in the phase diagrams in Fig. 2, i.e., the phase transition with vanishing baryon chemical potential. When the external field h that breaks the chiral symmetry explicitly is nonzero, the second-order phase transition becomes a continuous crossover, as shown in Fig. 5. One can define a pseudo-critical temperature T pc , which is the peak position of the curve χ σ versus T , and thus the reduced pseudo-critical temperature reads One observes from Fig. 5 that with the increasing h, the peak height of the susceptibility decreases and the pseudo-critical temperature t pc increases. The rescaling relation between t pc and h as well as that between the peak height of χ σ and t pc reads and see, e.g., [104] for more details.
In Fig. 6 we show the logarithm of the reduced pseudocritical temperature versus the logarithm of the reduced external field strength, and the logarithm of the peak height of the susceptibility versus the logarithm of the reduced pseudo-critical temperature in the left and right panels, respectively. The phase transition is also chosen to be near the location of the red cross in the phase diagrams in Fig. 2 for the O(4) symmetry universality class, where the baryon chemical potential is vanishing. Linear fitting to the calculated discrete data in Fig. 6 yields β = 0.403 (19) and γ = 1.543 (15) for the LPA, and β = 0.405 (22) and γ = 1.454 (17) for the LPA , which are in agreement with the relevant values in Eq. (35) and Eq. (43) within errors for the O(4) second-order phase transition in 3-d space. In turn, the agreement of critical exponents obtained from different scaling relations also provides us with the necessary check for the inner consistency of computations. Note, however, that the critical exponents β and γ determined from the scaling relations in Eq. (46) are significantly less accurate than those in Eq. (35) and Eq. (43).
As another check for the consistency, we consider the susceptibilities in the chiral broken phase near the coexistence line, i.e., x = −1, with t < 0 and h → 0. Inserting Eq. (39) into Eqs. (31) (32), one is led to when the system is near the coexistence line, one has x → −1 and y ∼ h/(−t) βδ . Hence, the transverse susceptibility is readily obtained as follows In order to obtain a similar expression for the longitudinal susceptibility, one needs further information on the equation of state y(x). As the system is located in the broken phase near the coexistence line, the dynamics is dominated by Goldstone modes, which are massless in the chiral limit. The relevant critical behavior in this regime is governed by a Gaussian fixed point, and thus the corresponding exponents are as same as values of mean fields [105,106], which leaves us with and see, e.g., [88,93] for more relevant discussions. Substituting equation above into Eq. (47), one arrives at As Eqs. (49) (51) show, the transverse and longitudinal susceptibilities are proportional to the external field with different powers in the broken phase, i.e., −1 and −1/2 for the former and latter, respectively. In Fig. 7 we show ln(χ π ) and ln(χ σ ) versus ln(−t) with a fixed value of the reduced external field h = 8.4 × 10 −9 in the chiral broken phase near the coexistence line. Similarly, here we only consider the phase transition of O(4) symmetry with µ B = 0 in the phase diagrams in Fig. 2. As shown in Eqs. (49) (51), the ratios of the linear fitting to ln(χ π )-ln(−t) and ln(χ σ )-ln(−t) are just the values of β and β − (βδ/2), respectively. Consequently, one arrives at β = 0.3979 (41) and δ = 4.984 (74) in LPA and β = 0.3832 (54) and δ = 4.86 (10) in LPA , which agree very well with the relevant values in Eq. (35) and Eq. (37).

D. Correlation length
It is well known that the correlation length ξ, plays a pivotal role in the critical dynamics, since fluctuations of wavelength ∼ ξ are inevitably involved in the dynamics. As a system is approaching towards a second-order phase transition, the most relevant degrees of freedom are the long-wavelength modes of low energy, and the correlation length is divergent at the phase transition [114].
The critical behavior of correlation length is described by the critical exponent ν, as shown in Eq. (21). In the symmetric phase t > 0, it reads which illustrates the scaling relation between the correlation length and the reduced temperature. Moreover, one can also define another critical exponent ν c related to the scaling relation between the correlation length and the reduced external field, to wit, In our setup in the quark-meson LEFT, cf. Sec. II, the correlation length is proportional to the inverse of the renormalized σ-meson mass, viz., where m σ is related to the dimensionless k-dependent sigma massm σ,k in Eq. (5) via the relation as follows where the scale k is chosen to be in the IR limit k → 0, and the mass is calculated on the equation of motion of the order parameter field. In Fig. 8 we show the scale relation between the correlation length and the reduced external field strength, and that between the correlation length and the reduced temperature, respectively. In the same way, we adopt the two different truncations: LPA and LPA . The phase transition points are also chosen to be at the locations of the red and blue crosses in the phase diagrams in Fig. 2 for the O(4) and Z(2) universality classes, respectively. By the use of the linear fitting to the calculated data, one obtains values of the critical exponent ν as follows   [110,111]. Moreover, results from other approaches, such as the conformal bootstrap for the 3-d conformal field theories (CFTs) [112,113], Monte Carlo simulation [98], and d = 3 perturbation expansion [99], as well as the mean-field values of exponents are also presented. Note that values with an asterisk are obtained with scaling laws in Eq. (18).
Finally, we close Sec. IV with a summary of various critical exponents calculated in this work in Tab. II. Respective results for the O(4) and Z(2) symmetry universality classes with truncation LPA or LPA are presented in the first several rows in Tab. II. As we have discussed in Sec. II, the effective potential is expanded as a sum of Chebyshev polynomials in our calculations, which captures global properties of the order-parameter potential very well. In Tab. II we also present values of critical exponents obtained from other computations, e.g., scalar theories calculated within the fRG with the effective potential expanded in a Taylor series [88,[107][108][109], or discretized on a grid [92], quark-meson LEFT within the fRG in LPA [93], derivative expansion of the fRG up to orders of O(∂ 4 ) and O(∂ 6 ) [110,111], the conformal bootstrap for the 3-d conformal field theories [112,113], Monte Carlo simulation [98], and the d = 3 perturbation expansion [99]. One observes that our calculated results are in good agreement with the relevant results from previous fRG calculations as well as those from the conformal bootstrap, Monte Carlo simulation, and the d = 3 perturbation expansion. Remarkably, the calculation with the truncation LPA is superior to that with LPA, and the former has already provided us with quantitative reliability for the prediction of the critical exponents in comparison to other approaches.

V. SUMMARY
QCD phase structure and related critical behaviors have been studied in the two-flavor quark-meson low energy effective theory within the fRG approach in this work. More specifically, we have expanded the effective potential as a sum of Chebyshev polynomials to solve its flow equation. Consequently, both the global properties of the effective potential and the numerical accuracy necessary for the computation of critical exponents are retained in our calculations. Moreover, we have employed two different truncations for the effective action: one is the usually used local potential approximation and the other is that beyond the local potential approximation, in which a field-dependent mesonic wave function renormalization is encoded.
With the numerical setup within the fRG approach described above, we have obtained the phase diagram in the plane of T and µ B for the two-flavor quark-meson LEFT in the chiral limit, including the second-order phase transition line of O(4), the tricritical point and the first-order phase transition line. Furthermore, we also show the Z(2) line in the phase diagram, which is the trajectory of the critical end point moving with the successive variance of the strength of explicit chiral symmetry breaking, or the varying pion mass.
In the phase diagram, we have performed detailed scaling analyses for the 3-d O(4) and Z(2) symmetry universality classes, and investigated the critical behaviors in the vicinity of phase transition both in the chiral symmetric and broken phases. Moreover, the transverse and longitudinal susceptibilities of the order parameter have been calculated in the chiral broken phase near the coexistence line.
A variety of critical exponents related to the order parameter, chiral susceptibilities and correlation lengths have been calculated for the 3-d O(4) and Z(2) symmetry universality classes in the phase diagram, respectively. The calculated results are also compared with those from previous fRG calculations, either employing the Taylor expansion for the order-parameter potential or discretizing it on a grid, derivative expansion of the effective action, the conformal bootstrap, Monte Carlo simulations, and the d = 3 perturbation expansion. We find that the critical exponents obtained in the quarkmeson LEFT within the fRG approach, where the orderparameter potential is expanded in terms of Chebyshev polynomials and a field-dependent mesonic wave function renormalization is taken into account, are in quantitative agreement with results from approaches aforementioned. Furthermore, we have also investigated the size of the critical regime, and it is found that the critical region in the QCD phase diagram is probably very small, and it is smaller than 1 MeV in the direction of temperature. withh y,k = h y,k Z q,k (Z φ,k ) 1/2 . (A10) Note that threshold functions BB (2,2) , F (2) and F (3) in Eq. (A9) can be found in e.g., [43,75].
Appendix B: Some relations for the Chebyshev polynomials In this appendix we collect some relations for the Chebyshev polynomials, which are used in solving the flow equation for the effective potential in Eq. (10). The Chebyshev polynomial of order n reads T n (x) = cos n arccos(x) , with nonnegative integers n's and x ∈ [−1, 1]. The explicit expressions for the Chebyshev polynomials could be obtained by the recursion relation as follows with T 0 (x) = 1 and T 1 (x) = x.
The N + 1 zeros of T N +1 (x) in the region −1 ≤ x ≤ 1 are given by A discrete orthogonality relation is fulfilled by the Chebyshev polynomials, to wit, Therefore, the zeros in y corresponding to Eq. (B3) read y k = y max − y min 2 cos π(k + 1 2 ) N + 1 + y max + y min 2 , with k = 0, 1, · · · N . Then, a function f (y) with y ∈ [y min , y max ] can be approximated as where the coefficients could be readily obtained by the use of the orthogonality relation in Eq. (B4), which yields with i = 0, 1, · · · N . With the Chebyshev approximation of the function f (y) in Eq. (B8), it is straightforward to obtain its derivative, viz.
where the coefficients d i 's can be deduced by the recursion relation, that reads