Mass sensitivity of the three-flavor chiral phase transition

The mass sensitivity of the chiral phase transition of QCD with and without axial $U_A(1)$-symmetry breaking at vanishing and finite quark chemical potential is investigated. To focus on the low-energy sector of QCD, a quark-meson model with three dynamical quark flavors is employed. Non-perturbative quantum fluctuations are taken into account with the functional renormalization group. The inherent ambiguities in fixing the low-energy model parameters away from the physical mass point and their consequences for spontaneous chiral symmetry breaking are discussed in detail and a heuristic parameter fixing scheme motivated by chiral perturbation theory is proposed. The influence of vacuum and thermal fluctuations of quarks and mesons on the order of the chiral phase transition is additionally assessed with a mean-field analysis.

The mass sensitivity of the chiral phase transition of QCD with and without axial UA(1)-symmetry breaking at vanishing and finite quark chemical potential is investigated. To focus on the lowenergy sector of QCD, a quark-meson model with three dynamical quark flavors is employed. Nonperturbative quantum fluctuations are taken into account with the functional renormalization group. The inherent ambiguities in fixing the low-energy model parameters away from the physical mass point and their consequences for spontaneous chiral symmetry breaking are discussed in detail and a heuristic parameter fixing scheme motivated by chiral perturbation theory is proposed. The influence of vacuum and thermal fluctuations of quarks and mesons on the order of the chiral phase transition is additionally assessed with a mean-field analysis.

I. INTRODUCTION
Quantum Chromodynamics (QCD) at finite temperature and density predicts a phase transition at low energies from confined hadronic matter to a deconfined quark-gluon plasma. The nature of the QCD transition is subject of intense studies and relevant to ongoing and planned heavy-ion experiments [1][2][3][4][5][6]. The global QCD symmetries associated to this transition are the chiral and center symmetries, which are realized only in two antipodal extreme limits of the quark masses. In the limit of infinitely heavy quarks, i.e. the pure gauge theory, the QCD vacuum obeys center symmetry. In this limit a first-order transition to a phase with spontaneously broken center symmetry occurs at a critical temperature of T c ∼ 270 MeV [7]. In the opposite limit of vanishing quark masses the QCD action is invariant under global chiral symmetry and, at least for three massless flavors, also a first-order transition is expected [8].
For finite quark masses both chiral and center symmetry are explicitly broken. As a consequence, the associated phase transitions weaken away from both quark mass limits and finally terminate at critical points of second-order transitions that belong to the threedimensional Z 2 Ising universality class. For physical quark masses the chiral and deconfinement transitions are smooth analytic crossovers with approximately coinciding pseudocritical temperatures [9][10][11][12].
For the case of two degenerate light (up and down) quarks and one strange quark the mass dependence of the order of the phase transitions is summarized in the Columbia plot [13]. In the heavy quark regime the Columbia plot, and even its extension to finite baryon density, is by now well understood both from continuum and lattice studies, e.g. [14][15][16][17][18]. In the latter case, the sign problem becomes treatable due to the large quark masses. The region of small masses is less well understood, see e.g. [19][20][21][22][23][24][25][26] for lattice studies of the location of the chiral critical line for three degenerate quark flavor. However, the majority of lattice results confirm the existence of a first-order region close to the three-flavor chiral limit. Although improved calculations indicate a small first-order region, the results still show huge discrepancies due to the notoriously difficult implementation of chiral fermions on the lattice.
Another longstanding problem concerns the order of the chiral transition in the chiral limit of two-flavor QCD. In the Columbia plot, this is the point of massless light quarks and infinitely heavy strange quarks. There is also a history of conflicting lattice results in this case, e.g. [27][28][29][30]. The two-flavor case is particularly intriguing since it is expected that the order of the phase transition is most sensitive to the fate of the global U A (1) symmetry [8]. While this symmetry is anomalously broken in the QCD vacuum [31], its restoration is expected at large temperatures [32]. If the restoration occurs close to the two-flavor chiral phase transition in the chiral limit, the transition can be of first order. Otherwise, a secondorder transition is expected and depending on the axial anomaly strength at the chiral transition, the associated critical behavior can be SU (2) L × SU (2) R ∼ = O(4) or U (2) L ×U (2) R /U (2) V whose corresponding critical exponents are very similar [8,[33][34][35][36]. In case of a second-order transition in the two-flavor and a first-order transition in the three-flavor chiral limit, the O(4) critical line of the light quark chiral limit and the Z 2 critical line at finite quark masses meet at a tricritical strange quark mass m * s . The existence and its relative location of this tricritical point with respect to the physical mass point is not established yet [37] and even the possibility m * s → ∞ cannot be excluded [38].
An additional open problem is the interplay between the three-and two-flavor chiral theories at finite chemical potential µ. It is not known whether a tricritical point m * s at µ = 0, if it exists, is analytically connected by a tricritical line to the tricritical point at finite µ in the two-flavor theory. For a positive curvature of the chiral critical surface, i.e. when the first-order region around the chiral limit increases with µ, both theories are connected analytically most probably and hence a critical endpoint exists for physical quark masses. However, if the size of the first-order region shrinks for increasing µ the situation can completely change. Even if the two-flavor chiral theory exhibits a tricritical point at some finite µ it need not be (analytically) connected to the chiral critical surface which thus hampers any conclusions about the existence of a critical endpoint in the phase diagram [39].
In summary, the nature of the chiral transition in the small mass region of the Columbia plot is still controversial concerning several lattice QCD studies for two and three quark flavors. However, functional continuum approaches such as the functional renormalization group (FRG) or Dyson-Schwinger equations (DSE) do not suffer a sign problem at finite density and the implementation of chiral fermions is straightforward. Renormalization-group methods such as the FRG are also best suited to investigate universal features of phase transitions. Much progress towards an understanding of the QCD phase structure from first principles has been made in recent years with these methods, e.g. [40][41][42][43][44][45] and references therein. For example, FRG results predict a second-order transition in the two-flavor chiral limit [46].
Valuable qualitative and quantitative insights into the chiral phase structure can be gained from effective lowenergy QCD models such as the linear sigma (LSM), quark-meson (QM) or Nambu-Jona-Lasinio (NJL) models. These models share the same global symmetries with QCD by construction and exhibit similar or even the same symmetry breaking pattern as the chiral transition in QCD. Their criticality belongs to the same universality class as QCD. Initiated by the first renormalization group analysis of the LSM in [8], the investigation of the Columbia plot with low-energy models has a long history. For example, an investigation of the influence of the axial U A (1) anomaly on the chiral phase transition within the Cornwall-Jackiw-Tomboulis (CJT) formalism in Hartree approximation of the LSM yields a large first-order chiral transition region in absence of the anomaly. Including the anomaly the results depend drastically on the sigma meson mass and no first-order transition is found at all in the most realistic case [47]. In [48,49] a zero temperature parameterization for the mass dependence of the couplings based on chiral perturbation theory (χPT) in an optimized perturbation theory scheme for the LSM results in a first-order region around the chiral limit which persists up to large strange quark masses. In contrast, NJL models on a mean-field level predict a very small first-order region [50]. In optimized perturbation theory and the mean-field treatment of the QM model a large first-order region in the chiral and in the light chiral limit is obtained [51,52]. Interestingly, in such studies, the Columbia plot exhibits features which are qualitatively independent of the axial anomaly [53].
However, one drawback of these model investigations is the parameter dependency which in parts can be reduced by combining effective low-energy QCD models with the FRG such that an improved non-perturbative truncation can be obtained in a systematical manner, see e.g. [54][55][56][57].
In this work, we go beyond these previous studies and perform a non-perturbative FRG study of the Columbia plot in a three-flavor quark-meson truncation at vanishing and finite density as well as with and without the axial U A (1) anomaly for the first time. The FRG allows us to incorporate arbitrarily high loop orders as well as genuine non-perturbative effects [58,59]. A recent FRG study for a quark-meson truncation with N f = 2 + 1 dynamical quark flavors has been done in [60] in a leadingorder derivative expansion, the local potential approximation (LPA), where the effective potential is treated as energy scale dependent, while all other correlation functions are kept constant. It was found that the chiral transition in the light chiral limit at physical strange quark masses is of second-order belonging to the O(4)universality class in the presence of the axial anomaly. Without the anomaly, the transition was of first order in this limit. The influence of scale dependent correlation functions beyond the leading LPA has been worked out in [61]. However, the full Columbia plot has not been studied with the FRG yet. Our approach allows us to systematically study the impact of thermal and vacuum fluctuations of quarks and mesons on the order of the chiral phase transition for arbitrary quark masses.
Any investigation of the Columbia plot with a lowenergy QCD description faces the problem to fix the initial parameters unambiguously. At the physical mass point the parameters can be adjusted to reproduce certain experimental observables. Away from the physical point no such phenomenological guidance can be found. We discuss the physical consequences of different parameter fixing procedures and propose a heuristic fixing scheme motivated by chiral perturbation theory (χPT). This paper is organized as follows: after the introduction of our effective low-energy description of QCD in terms of a quark-meson model in Sec. II, the incorporation of non-perturbative quantum fluctuations by the FRG is summarized in the following Sec. III. In Sec. III B the vacuum fluctuations of quarks, which are relevant for QCD thermodynamics, are included into the conventional mean-field approximation. Owing to the ambiguity in the model parameters away from the physical point, two different parameter fixing schemes are introduced in Sec. III C. Our numerical findings on the chiral critical lines in the Columbia plot at vanishing and finite chemical potential with and without the axial anomaly are presented in Sec. IV. A summary with an outlook can be found in Sec. V. Details of our numerical implementation are given in the App. A.

II. THREE-FLAVOR QUARK-MESON MODEL
In this section we briefly describe our setup of the three flavor quark-meson model as a low-energy effective model for QCD. It is based on a linear sigma model where the (pseudo)scalar meson nonets are represented by a colorsinglet N f ×N f matrix field Σ. It shares the same global chiral symmetry properties as QCD. The inclusion of constituent quarks then facilitates the study of finite chemical potential and significantly improves the chiral dynamics such that, e.g., the chiral critical temperature at vanishing density is much closer to lattice QCD, cf. e.g. [53,62]. In general, quark-meson models arise from QCD by successively integrating out quark and gluon fluctuations towards low energies, see e.g., [44,45,56,63,64].
Under the assumption that gluon fluctuations decouple below a certain UV scale Λ, the U L (N f ) ⊗ U R (N f ) chiral invariant part of the Euclidean effective action is given by with the three flavor quark field q = (u, d, s) T [65,66]. The scalar and pseudo-scalar meson fields σ a and π a are comprised in the complex matrix fields and Σ 5 = T a (σ a +iγ 5 π a ) where T a label the U (N f ) group generators. They are given explicitly by T 0 = 1 3 / √ 6 and T i = λ i /2, with the eight Gell-Mann matrices λ i , i = 1, . . . , 8. The action includes standard kinetic terms for the quark and meson fields as well as a Yukawa type interaction vertex ∼ h. A quark chemical potential matrixμ = diag(µ u , µ d , µ s ) is included. In the following we restrict the discussion to symmetric quark matter µ ≡ µ u = µ d = µ s . We consider chirally invariant meson-meson interactions of arbitrary order within an effective chiral potential U χ (ρ 1 , ρ 2 ) which is a function of the chiral invariants ρ n = Tr (Σ † Σ) n , n ≤ N f . In the effectve potential a third invariant ρ 3 has been dropped for the sake of simplicity. If the chiral potential U χ would only depend on one invariant ρ 1 , it would exhibit an enhanced SO(18) flavor symmetry. Adding a second invariant ρ 2 breaks this symmetry down to the desired chiral [65,67]. Hence, the general dependency on ρ 1 and ρ 2 represents a minimal requirement to capture the correct flavor symmetry.
To effectively account for finite current quark masses as well as the axial U A (1) anomaly, we add the following symmetry breaking terms to the action: where the first term models the explicit breaking of the chiral SU L (N f ) ⊗ SU R (N f ) due to finite quark masses.
The anomalous breaking of the U A (1)-symmetry is described by the second term which represents a bosonized version of the 't Hooft determinant, ξ = det(Σ † ) + det(Σ) [31,68,69]. The determinant gives a N f -meson interaction vertex which leaves SU L (N f ) ⊗ SU R (N f ) ⊗ U V (1) intact but explicitly breaks the U A (1)-symmetry down to Z A (N f ). Later, we will effectively explore the effect of the anomalous breaking and consider the two constant cases: c A = 0 and c A = 0. However, in QCD this coupling has also an explicit temperature and/or density as well as a scale dependency which we ignore in the present work. For exploratory studies see e.g. [54,[70][71][72][73]. The vacuum expectation values (VEV) for the fields are found at the stationary point of the action. Due to parity and (approximate) isospin symmetry as well as flavor neutrality of the vacuum only the diagonal components of the scalar meson sector related to the generators T 0 and T 8 can have finite expectation values. Hence, we assume light isospin symmetry and write q = (l, l, s) T . At the stationary point we can therefore express the chiral invariants and the explicit symmetry breaking action, Eq. (3), in terms of the light and strange scalar meson VEVsσ l = σ l andσ s = σ s as follows: σ l andσ l are directly related to the light and strange chiral condensates l l and ss . The fields denoted by the indices l, s are expressed in the light-strange basis. They are related to the ones in the singlet-octet basis, Eq. (2), via the rotation The pion and kaon decay constants can be expressed in terms of the light and strange meson VEVs as [62]: which yield the quark masses The current quark mass sensitivity of the chiral phase transition can be controlled by varying the explicit symmetry breaking parameters j l and j s . Since we have no direct access to the perturbative QCD regime within the used low-energy description, it is more sensible to express the Columbia plot in terms of the purely light pion-and the open strange kaon masses, m π and m K , instead of the corresponding current quark masses. It is therefore useful to derive relations between the explicit symmetry breaking sources and the pion and kaon masses, The light-quark chiral limit j l = 0 always implies m π = 0, while the strange quark chiral limit j s = 0 leads to the relation m 2 K = f π m 2 π /(2f K ). As a consequence, kaons are always massive als long as the pions are. Eq. (8) can be derived by using the equations of motion ∂(U χ + L SB )/∂σ l,s = 0 together with the explicit expressions for m 2 π and m 2 K in App. D of [61].

III. NON-PERTURBATIVE QUANTUM FLUCTUATIONS
For an accurate description of the chiral phase transition all quantum, thermal and density fluctuations, in particular in the vicinity of the transition, should be taken into account in a non-perturbative manner. As mentioned in the introduction this can be achieved with the FRG, which is a semi-analytical method providing a non-perturbative regularization and renormalization scheme for the resummation of an infinite class of Feynman diagrams. For QCD related FRG reviews, see e.g. [58,59,[74][75][76][77].

A. Functional Renormalization Group
The central object of the FRG is a scale dependent effective action Γ k that is identified with the classical action S = Γ Λ at some UV-cutoff scale Λ. By successively integrating out fluctuations momentum shell by momentum shell, the full quantum effective action Γ = Γ 0 , i.e. the generating functional of all one-particle irreducible correlation functions, is recovered in the infrared (IR) at k → 0. The renormalization group flow of Γ k is governed by the Wetterich equation [78], which for the quark-meson model reads Here, Γ k denotes the second functional derivative w.r.t. the corresponding meson or quark fields. Fluctuations with momenta p 2 < k 2 are suppressed by appropriate regulator functions R Σ,q k . The super trace STr indicates the trace over discrete (color, flavor, spinor) and continuous (loop momenta) indices as well as various fields. The Wetterich equation (9) represents an exact functional equation whose solution involves an infinite tower of partial differential equations. Hence, in practice Γ k has to be truncated. In this work we choose the local potential approximation (LPA) which amounts to a scale dependence of the meson potential U χ,k (ρ 1 , ρ 2 ). All other parameters are k-independent. Thus, the scale-dependent effective action reads (10) where S χ,k is given by Eq.
(1) with a scale-dependent effective potential U χ → U χ,k . The symmetry breaking part L SB , Eq. (3), is unchanged. In this approximation, and by employing optimized regulators for the spatial momenta [79][80][81], the Wetterich equation turns into the following partial differential equation for the effective potential: The quark and meson quasi particle energies are abbreviated as E i = k 2 + m 2 i with the corresponding light quark and strange quark masses, m l = h σ l /2 and m s = h σ s / √ 2. The meson (curvature) masses m b are obtained from diagonalizing the Hessian of the full effective potential U χ,k +L SB , for explicit expressions we refer to the appendix of Ref. [61]. Details on the truncation and regularization of the effective action can be found in Ref. [60]. We remark that improved LPA truncations as well as a different regulator choice are likely to have some quantitative effects on the presented results. In the present context the former has been studied in [61,82] and the latter has been discussed in [83]. However, the qualitative statements in this work are expected to be unaffected by truncation and regularization effects.

B. Extended Mean-Field Approximation
In order to understand the effect of the mesonic fluctuations we also solve the quark-meson model in a meanfield approximation, which can conveniently be achieved by switching off the meson contributions to the flow equation (11). The mesonic quantum fields are replaced by their vacuum expectation values on the level of the action. The remaining contribution is quadratic in the quark fields and thus yields the conventional Gaussian path integral for the quark contribution to the effective potential.
The quark contribution contains an UV-divergent vacuum part which requires regularization. In standard (or no-sea) mean-field approximation (MFA) this contribution is simply ignored. However, as shown in [84] the  vacuum contribution has a significant impact on the chiral dynamics of the QM model. In the present context, the order of the chiral transition in the light chiral limit depends crucially on fermionic vacuum fluctuations [85]. For these reasons, we also explore the Columbia plot in the extended mean-field approximation (eMFA), which takes these vacuum fluctuations into account.
We want to emphasize that the vacuum contribution is automatically included, and properly regularized, in the FRG approach. This also includes the mesonic vacuum fluctuations. Hence, by switching off the meson loops in Eq. (11) we directly work in the eMFA. Alternatively, the vacuum contribution of the QM model in eMFA can also be regularized with standard techniques like e.g. dimensional regularization yielding a renormalized QM model potential, e.g. [85]. However, performing a MFA within the FRG it is guaranteed that the same non-perturbative regularization and renormalization schemes are employed.

C. Parameter Fixing
Here, we discuss the initial effective action we use to initiate the RG flow of our model. The general idea is to fix the couplings of the effective action Γ k=Λ , Eq. (10), at a UV scale Λ such that the resulting effective action Γ k=0 reproduces phenomenologically known quantities in the infrared, like meson masses and decay constants. The RG evolution is determined by the flow equation Eq. (9). For initial scales Λ larger than the scale of chiral symmetry breaking (in our case around k χSB ≈ 500 MeV) the VEVs of the mesons are small. The meson masses are large and consequently their fluctuations are negligible. Therefore it is reasonable to expand the initial effective meson potential about small VEVs of the mesons and keep only the leading terms. It turns out for Λ 700 MeV finite initial values for marginal and relevant terms are sufficient. Furthermore, choosing Λ 700 MeV also guarantees that UV cutoff artifacts at the relevant temperature and chemical potentials are negligible. For a recent discussion regarding UV cutoff effects in low-energy effective models we refer to [86]. Irrelevant couplings, i.e. those with a negative mass dimension, are suppressed due to small initial meson fluctuations. Their initial values can be set to zero. This yields the initial potential where we have introduced a modified chiral invariant ρ 2 = ρ 2 − ρ 2 1 /3. For this ansatz seven parameters {a 10,Λ , a 20,Λ , a 01,Λ , c A , j l , j s , h} have to be fixed in the UV. At the physical point, i.e. for realistic masses, we can fix the initial parameters such that we reproduce wellknown masses and decay constants in the IR. Note, that without anomalous breaking, c A = 0, one has m η = m π and m η is completely fixed by the kaon and pion masses and their decay constants, cf. [62]. For the mean-field approximation an identical parameter fixing procedure with the same UV cutoff scale Λ is employed, see [53] for further details. The used parameters are summarized in Tab. I.
We want to emphasize that irrelevant couplings, even though they are zero initially, are generated during the RG flow towards the IR. This has been demonstrated explicitly within a QM model in [82]. In QCD, this follows from the fact that the QCD effective action just above the chiral symmetry breaking scale is dominated by effective four-quark interactions whose resonance signals the spontaneous breaking of chiral symmetry [44,45,54,56]. Based on this, we assume that higher order mesonic operators are generated only after the decoupling of the gluons. These operators are quantitatively and qualitatively relevant for the non-universal physics considered here. We therefore solve the flow equation for the full effective potential ∂ k U χ,k (ρ 1 , ρ 2 ), Eq. (11), on a two-dimensional grid instead of resorting to an expansion in terms of only relevant and marginal operators, see App. A.
In order to investigate the mass sensitivity and the influence of the axial anomaly on the chiral phase structure the corresponding symmetry breaking parameters j l,s and c A are varied and thus the system is tuned away from the physical point. We cannot rely on IR observables to fix the initial conditions in this case. Since the initial UV parameters are fixed in general by the underlying fundamental theory of the strong interactions, it is a priori not clear how the remaining parameters of the low-energy model change when the current quark masses and/or the 't Hooft coupling are tuned away from the physical point.
For the axial symmetry parameter c A , two scenarios are considered (see Tab. I): one is where a constant c A is chosen such that the η -meson mass is fixed to its experimental value at the physical point. In the other case we assume that the axial anomaly is absent, i.e. c A = 0. The corresponding initial UV parameters are fixed in such a way that the same experimental IR values (as listed in Tab. I) are obtained. We want to emphasize that the latter case is not directly connected to QCD, at least for small temperatures, since the QCD vacuum always breaks the U A (1)-symmetry through quantum fluctuations. Only for high temperatures [32] or in the large-N c limit [87] this symmetry is restored. Thus, we implicitly assume the limit of infinite colors for the c A = 0 computations. Nonetheless, it is interesting to investigate the influence of a U A (1)-symmetry on the phase structure and study the consequences for the mass sensitivity from a more general point of view. For exploratory FRG studies about U A (1)-symmetry breaking, see [54,[70][71][72][73].
While varying the explicit chiral symmetry breaking j l,s is less delicate, it gives rises to ambiguities for the parameter fixing of low-energy models. In principle, the initial action of the low-energy model at a given scale Λ could be determined for varying current quark masses by solving the corresponding QCD flow: Starting with a microscopic QCD action at a perturbatively large energy scale k 1 GeV the RG evolution towards Λ would yield the desired initial conditions in an unique way. This is beyond the scope of the present work, see, e.g., [88] for an explicit example of a model parameter fixing procedure within an NJL model framework.
Owing to the scarce information about QCD away from the physical point, we explore two strategies to fix the parameters in this case: • fixed-UV scheme: Away from the physical point only the explicit symmetry breaking sources j l,s are assumed to change while all other parameters of the initial effective action are identical to the ones at the physical point. So the symmetric part of the initial effective action does not change in the UV.
• fixed-f π scheme: For each j l,s the initial effective action is adjusted such that the IR pion decay constant is always fixed to its physical value, i.e. f π = 92.4 MeV at every mass point in the Columbia plot.
The fixed-UV scheme is practically the simplest choice and it has been used frequently in the literature, e.g., [47,50,53,89,90]. In this scheme the initial action is given by the parameters listed in Tab. I at the cutoffscale Λ = 700 MeV and only j l,s are varied to explore the quark-mass sensitivity of the phase transition. The underlying assumption is that a change in the current quark masses of QCD can be mapped directly onto a change of the symmetry breaking sources j l,s in the effective low-energy model. The fixed-f π scheme is motivated by the findings of chiral perturbation theory [91]. In the light chiral limit, j l = 0, the pion decay constant only slightly decreases. In the chiral limit when also j s = 0, it does not change at leading order. Since χPT suggests a change of f π only on the order of 10% towards the chiral limit, we always adjust the initial conditions such to yield f π = 92.4 MeV in the IR, independent of j l,s . A more accurate adaption to a more precise quark-mass dependence of f π (and also other observables) will only lead to minor quantitative changes which are irrelevant in the present case. So in the fixed-f π scheme we demand (approximate) compliance with χPT. A similar parameter fixing procedure for a linear sigma-model based on χPT has been proposed in [48]. It is a priori possible that both schemes give the same results. However, we will demonstrate that there are drastic differences and discuss the physical reasons in the next section. The Yukawa coupling h and the anomaly coupling c A are not running in our current approximation. In addition, we assume that they do not change away from the physical point. The parameters that can potentially be adjusted in the fixed-f π scheme are a 10,Λ , a 20,Λ , a 01,Λ , j l and j s . Clearly, there is no unique choice for the initial parameters in this case since we assume that m π , m K and f π are the only IR quantities we know away from the physical point. An exception is the chiral limit, where also the relation between f π and f K as well as the masses of the pseudoscalar meson octet are fixed by chiral symmetry. Our model consistently reproduces this case. It turns out that for decreasing j l,s it is always possible to find a larger initial scale Λ > Λ with a 10,Λ = a 10,Λ , a 20,Λ = a 20,Λ and a 01,Λ = a 01,Λ (cf. Tab. I) such that the fixed-f π scheme is realized. This procedure is equivalent to a change in the initial parameters for a fixed Λ. It is numerically much more convenient since everything is controlled by only one parameter at given j l,s . In Tab. II we exemplify the initial conditions within the two employed schemes for three different choices of explicit symmetry breaking parameters j l,s , parametrized as

IV. NUMERICAL RESULTS
Lattice QCD calculations establish that for vanishing quark chemical potential the chiral three-flavor phase transition is a smooth crossover at the physical mass point [9]. This is also a generic outcome of linear sigmaand quark-meson models. In this work we extend the investigation of the chiral phase transition for arbitrary symmetry breaking sources j l and j s as well as finite chemical potential. We begin by reviewing various limiting cases: • j s → ∞: Increasing j s above its physical value, the strange quark and all mesons with strangeness become heavier and eventually decouple from the flow towards the infrared. This limit is equivalent to the case of only two degenerate light flavors and, depending on the fate of the U A (1) anomaly, will either resemble a O(4)-symmetric or an U (2) L ⊗ U (2) R -symmetric system.
• j l = 0 and j s > 0: We call this limit the light chiral limit. It corresponds to two massless light quark flavors but finite strange quark masses and coincides with the (m π = 0)-axis in the Columbia plot.
• j l > 0 and j s = 0: This limit corresponds to the strange chiral limit. The kaons then obey m 2 K = f π m 2 π /(2f K ). This defines the line of smallest possible kaon masses in the Columbia plot.
• j l = j s = 0: All quark masses vanish and chiral SU (3) L ⊗ SU (3) R -symmetry is not explicitly broken. Since this is the N f = 3 chiral limit, we refer to it as the chiral limit.

A. Spontaneous Symmetry Breaking
As discussed in Sec. III C, the model parameter fixing away from the physical point is not unique. In principle the effective action Γ Λ at a given initial UV scale Λ is determined by an RG flow from the perturbative QCD regime down to Λ. At least for the physical point this ambiguity in the initial action of the low-energy model can be circumvented by adjusting the initial action to known IR observables. As soon as the physical point is left this procedure is not constructive anymore.
In order to proceed we apply the in Sec. III C already introduced two parameter fixing methods when the sources j l,s are varied. The fixed-UV scheme, where all other parameters remain fixed (cf. Tab. I), has been used in a FRG study of the light chiral limit [60] and in a similar MFA study for the whole Columbia plot [53]. In the latter case a rather heavy sigma mass of the order of 800 MeV was required in order to observe spontaneous chiral symmetry breaking in the chiral limit, which is in conflict with recent experimental data [92].
While the fixed-UV scheme yields a reasonable value for f π ≈ 87 MeV in the light chiral limit which is in agreement with chiral perturbation theory, we find that both light and strange condensates vanish towards the chiral limit and spontaneous chiral symmetry breaking is lost. This is displayed in Fig. 1 where the light (left panel) and the strange (right panel) condensates are shown as a function of the temperature for different points on a path connectiong the physical point to the chiral limit in the Columbia plot. This path is parametrized by a factor α ≥ 0 as defined in Eq. (13), such that α = 1 corresponds to the physical mass point and α = 0 to the chiral limit. Both condensates melt with decreasing α and at the chiral limit they vanish. Furthermore, the (pseudo)critical temperature T c drops, in particular T c → 0 towards the chiral limit. The same behavior is seen for the RG scale of spontaneous chiral symmetry breaking, k χSB . Hence, the system remains in the symmetric phase in the chiral limit. This observation can be understood already on the mean-field level. In this case the initial meson potential Eq. (12) matches the full meson potential. For the light sector of the potential in the chiral limit, a finite expectation valueσ l can only be generated for a 10 < 0 if the potential is bounded from below, which in turn requires a positive quartic coupling (a 20 + a 01 /3) > 0. With the fixed-UV scheme we find for sigma masses in the range m σ ∈ [400, 600] MeV a positive a 10 and thus spontaneous symmetry breaking cannot occur. In the MFA this can be circumvented by choosing unphysically large sigma masses (m σ 800 MeV) [53]. The above mean-field argument cannot directly be transferred to the FRG case, since the full effective potential as a function of arbitrary powers of the chiral invariants is evaluated beyond mean-field. For instance, the FRG flow only allows for sigma masses in the range m σ ∈ [400, 600] MeV. Basically, m σ is controlled by the quartic meson coupling. Vacuum stability of the initial potential requires a positive quartic coupling which fixes the lower bound for m σ . On the other hand, the upper bound is given by the requirement that mesons have to decouple at large energy scales. Since the initial meson masses decrease with increasing quartic coupling, this coupling cannot be arbitrarily large. In order to understand why spontaneous chiral symmetry breaking is lost in the fixed-UV scheme, we need to understand the fluctuations that eventually drive the system to criticality.
We observed that the values of the condensates, as well as k χSB and T c , decrease with decreasing α, i.e. with less explicit symmetry breaking. At first glance, this seems to be counter-intuitive since smaller j l,s imply lighter current quark masses that yield a potential enhancement of symmetry breaking fermionic fluctuations. Furthermore, the quark masses are always lighter than the RG scale above k χSB , i.e. m l/s,k k for k ∈ [k χSB , Λ], so that any variations in the current quark masses should have less impact on k χSB . However, in the fixed-UV scheme the meson masses increase with decreasing sources j l,s at the initial scale Λ. This also includes the critical sigma mode which is related to the correlation length via ξ = 1/m σ . Conversely, the initial correlation length decreases with smaller j l,s such that (pseudo) criticality is reached later in the RG flow. Hence, k χSB , the condensates and, as a consequence, T c decrease with decreasing j l,s . The reason for the sensitivity of the initial meson masses on the explicit symmetry breaking parameters is that they change the evaluation point of the effective potential that defines the physical parameters and not the global form of potential itself, see, e.g., [82]. Thus, different j l,s can change the meson masses substantially.
This observation contradicts our current understanding of chiral dynamics in the chiral limit. For example, χPT predicts in the chiral limit only minor modifications of the decay constants as discussed in Sec. III C. In addition, this observation might propose that the lower bound of the conformal window, N c f , where spontaneous chiral symmetry breaking is lost for N f ≥ N c f chiral quarks in the fundamental representation, would be N c f = 3. However, functional continuum methods, supersymmetryinspired all-orders beta function approaches and lattice QCD studies suggest N c f ≈ 8 − 12, e.g., [55,[93][94][95][96] and references therein. It is therefore reasonable to conclude that the fixed-UV scheme is not applicable within the present low-energy theory. We note that this statement is unlikely to be affected by truncation or regularization scheme effects within the QM model. The main effect discussed above is the drastic increase of the initial meson masses with decreasing j l,s . It follows directly from the form of the initial effective action. Only if the RG flow is initiated at larger scales where gluon dynamics are relevant, our results would be independent of the form of the initial meson potential (as long as mesons are decoupled from the dynamics). This clearly goes beyond the lowenergy effective model employed here and is discussed in detail in the context of dynamical hadronization in [44,54,56].
In the following we assume that chiral symmetry of the QCD vacuum is spontaneously broken in the chiral limit. By demanding compliance with χPT we introduced the fixed-f π scheme in Sec. III C. It follows from Eq. (6) and the symmetric form of the initial effective potential (12) that fixing f π to a finite value for any j l,s always implies spontaneously broken chiral symmetry in the infrared. As discussed above, if we fix a 10,Λ , a 20,Λ and a 01,Λ at a given UV cutoff Λ, a decrease of the explicit symmetry breaking parameters j l,s not only decreases the pion and kaon masses in the IR but also leads to an unphysically large reduction of the initial correlation length. Since the quark fluctuations that drive the system to (pseudo) criticality do not change much for scales k > k χSB and decreasing j l,s , we expect that the same behavior holds for the correlation length. Due to the intimate relation between the decay constants and the chiral condensates, a larger UV correlation length also entails larger IR decay constants. So the fixed-f π scheme resolves the problems of the fixed-UV scheme by construction. At the physical mass point our initial conditions are fixed such that f K /f π ≈ 1.2, in agreement with lattice QCD [97]. As we will demonstrate below, our procedure naturally leads to f K = f π in the chiral limit, in accordance with chiral symmetry. Thus, the fixed f π -scheme, backed up in particular by χPT, seems to be a physically reasonable procedure to fix the initial conditions away from the physical mass point as it provides a natural interpolation between the physical point and the chiral limit. A more accurate procedure to fix the initial conditions that also takes the scaling, e.g., of f π , f K , m π and m K near the chiral limit as predicted by χPT into account is beyond the scope of the present work. But since the physical picture presented here is expected to hold in general, we anticipate only minor quantitative changes.
Our findings within the fixed f π -scheme are shown in Fig. 2 where the light and strange condensates (left and middle panel) as well as the pion and sigma meson masses (right panel) as a function of temperature for different UV cutoffs towards the chiral limit from Λ = 700 MeV at the physical point to Λ ∼ 1143 MeV in the chiral limit are shown. Again, a straight path in the Columbia plot from the physical point to the chiral limit has been chosen, see Eq. (13). The light condensateσ l is always fixed to f π = 92, 4 MeV at T = 0 by construction for every α. With decreasing α the chiral phase transition becomes steeper and turns into a first-order transition below a critical value α c 0.04. At α c (not shown in the figure) the sigma mass drops to zero at the critical temperature signaling a second-order phase transition. Since the sigma meson is the only critical mode at α c , the transition lies in the three-dimensional Ising universality class.
The constituent quark mass consists of the current quark mass and a contribution from spontaneous chiral symmetry breaking proportional to the chiral condensates. Thus, towards the chiral limit no significant changes for the light quarks are expected, while the strange quark constituent mass decreases substantially. Sinceσ s is connected to the strange constituent quark mass via Eq. (7), it also decreases significantly with decreasing α. In fact, in the chiral limit chiral symmetry implies the relationσ l = √ 2σ s , which results in f K = f π , cf. Eq. (6). The pseudoscalar meson octet with the pions, kaons and the η-meson becomes massless. The pseudoscalar singlet η -meson remains massive in the presence of the U A (1) anomaly. The masses of the scalar meson nonet with the sigma, f 0 , a 0 and the kappas, are generated solely by spontaneous chiral symmetry breaking. This is exactly what we find.

B. Columbia Plot in Mean-Field Approximation
We begin with an investigation of the Columbia plot in the extended mean-field approximation as outlined in Sec. III B. We express the symmetry breaking in terms of the pion and kaon masses, m π and m K , by using the identities in Eq. (8). We will concentrate on lower and moderate mass values in the Columbia plot since the opposite quenched mass limit is unreachable without direct access to gauge degrees of freedom. As parameter fixing we also apply the fixed f π -scheme by adjusting the UV cutoff Λ such that f π is fixed in the vacuum for any values of m π and m K . This extends the analysis in [53], where the Columbia plot has been computed in standard MFA where vacuum quark fluctuations are omitted. Together with the following FRG investigation, we are now able to distinguish between effects stemming from fermionic (vacuum as well as thermal) and/or mesonic fluctuations.
In standard MFA the chiral transition is of first order in the light chiral limit, apparently independent of m K and of the presence of the axial anomaly. Particularly, this would imply a first-order chiral transition in the twoflavor limit.
In Fig. 3 the mass sensitivity of the chiral phase transition in the (m π , m K )-plane as well as for finite chemical potentials in eMFA are shown. The results differ significantly from the MFA results in [53]. Including the anomalous U A (1)-symmetry breaking (left panel) we find a distinctive first-order region around the chiral limit for vanishing chemical potentials. The transition is of second-order in the light chiral limit for kaon masses larger than the physical ones and in particular in the two flavor limit. These findings agree with the predictions in [8] for the three and two flavor chiral limits. Since standard MFA predicts a first-order transition [84], our results imply that the second-order transitions in the light and two-flavor chiral limits are induced by fermionic vacuum fluctuations. Similarly, and in contrast to a standard MFA study [53], a tricritical point m * K = 169 MeV on the m K -axis for µ = 0 is found when quark vacuum fluctuations are taken into account. At this point the second-order chiral transition line in the light chiral limit terminates and the transition turns from second-to firstorder with decreasing m K . In the opposite direction for finite pion masses the chiral critical line terminates at the j s = 0 axis (dashed orange line in the figure) corresponding to m π ≈ 110 MeV. Everywhere else the transition is a smooth crossover for µ = 0. For µ > 0 the size of the first-order region around the chiral limit increases such that the chiral critical surface has a positive curvature, cf. [98]. Consequently, a chiral critical endpoint exists in the (T, µ)-phase diagram for physical masses which is marked as a red point at the vertical line in Fig. 3 and represents the intersection point with the chiral critical surface [22]. In this scenario the found tricritical point m * K at µ = 0 is expected to be analytically connected by a tricritical line to the tricritical point at some finite µ in the two flavor chiral limit by increasing m K [99]. This behavior is also confirmed in Fig. 3: for increasing µ and larger m K the chiral critical line has a positive slope and saturates for m K 500 MeV at µ c = 239 MeV and is continuously connected to the critical endpoint in the two-flavor chiral limit. For chemical potential below the chiral critical line the chiral transition is always of second order for m π = 0 and first-order above.
Without the U A (1)-symmetry breaking (right panel in Fig. 3) no first-order phase transition for µ = 0 is found at all. This is in disagreement to [8] as well as standard MFA results [53]. A second-order transition is found in the light chiral limit independent of m K when the quark vacuum fluctuations are included and a first-order one if not. Hence, fermionic vacuum fluctuations induce a second-order transition in the (light) chiral limit. A first order region is only seen at larger chemical potentials. In the chiral limit at µ c = 232 MeV which rises to µ c = 256 MeV for a kaon mass around m K = 243 MeV and then slowly declines again. At the physical point, we find a critical endpoint at µ c ∼ 298 MeV.
Note that in the standard MFA analysis [53] an extended first-order region along the m K -axis is seen independent of the U A (1)-symmetry breaking. However, in order to obtain a spontaneous symmetry breaking in the chiral limit a large sigma meson mass was chosen. We have also verified that with vacuum fluctuations and a comparatively large sigma meson mass similar results as in Fig. 3 are obtained. This confirms that the differences we found can indeed be attributed to fermionic vacuum fluctuations.

C. Columbia Plot with the FRG
In Fig. 4 the Columbia plot for µ = 0 with vacuum and thermal quark and meson fluctuations within the FRG framework is presented. On the m K -axis a second-order chiral phase transition line is found that terminates at the tricritical point m * K ≈ 23 MeV, which is significantly smaller than the physical mass and the eMFA result. For smaller kaon masses surrounding the chiral limit a small first-order region is formed that is limited by a secondorder line in agreement with the predictions of [8]. The physical mass point is far outside the plot range and lies deep within the crossover region. In comparison to the eMFA findings, the first-order region around the chiral limit shrinks drastically upon the additional inclusion of  It is important to stop the RG flow as deep in the IR as possible towards the chiral limit. The reason is that the lowest IR mass scale in the system determines when the RG flows freeze out and the mesons in the pseudoscalar octet become very light close to the chiral limit. Integrating out fluctuations at very small scales, k m π /2, is numerically expensive due to the convexity of the effective potential where the flow in the flat potential region approaches a singularity [82,100]. We used k IR = 70 MeV throughout this work. In particular, the location of the chiral critical line close to the tricritical point is sensitive to IR-cutoff effects. The gray shaded area in Fig. 4 gives a rough estimate for the sensitivity of the critical line on k IR . The upper boundary of this area was extracted at k IR = 100 MeV while the lower boundary (red line) denotes the critical line at k IR = 70 MeV. Between k = 100 MeV and k = 70 MeV the tricritical point moves towards lower kaon mass by ∆m K ∼ 6 MeV while the location of the critical line is already frozen for k 100 MeV and m π 10 MeV. We therefore expect only very minor changes for even smaller k IR so that our main finding, the surprisingly small first-order region around the chiral limit, will not be altered. Scaling properties at the critical lines/points and a precise determination of m * K are not in the scope of the present work.
The Columbia plots obtained with the FRG for vanishing as well as finite µ with and without the axial anomaly are shown in Fig. 5. With the axial anomaly (left panel in the figure) the first-order region at µ = 0 around the chiral limit is significantly smaller than the corresponding region in eMFA. To be more precise: to-wards the SU (3)-symmetric chiral limit, i.e. where the light and the strange quark masses coincide and consequently m π = m K , the boundary of the first-order region at µ = 0 is located at a critical pion mass m c π ∼ 17 MeV compared to m c π ∼ 86 MeV in the eMFA and m c π ∼ 150 MeV in the standard MFA. Thus, the conclusion is that vacuum quark as well as meson fluctuations significantly reduce the size of the first-order region in the Columbia plot around the chiral limit at µ = 0.
The second-order chiral critical line extends from the tricritical point at m * K ≈ 23 MeV along the m K -axis to the two-flavor chiral limit. For µ > 0 the critical line steeply rises to µ ≈ 200 MeV around m K ≈ 50 MeV and then slowly grows to µ ≈ 295 MeV at physical kaon masses. Similar to the eMFA case the chiral critical surface has a positive slope such that the size of the firstorder region around the chiral limit increases with increasing µ. Also shown in the left panel of Fig. 5 is the extrapolation to the two-flavor chiral limit. An explicit FRG calculation with a chiral SU (2) ⊗ SU (2)-symmetric quark-meson truncation with the same sigma meson mass and f π = 92.4 MeV yields a critical µ c = 281 MeV which is also shown in the left panel of the figure. However, this comparison should be viewed with some caution since in a chiral SU (2)⊗SU (2)-symmetric setup the axial U A (1)symmetry is maximally broken while in our N f = 2 + 1 setup the U A (1)-symmetry breaking, albeit large, is still finite.
Below the chiral critical line the transition for m π = 0 is always of second-order as in the eMFA case but the area is clearly larger in the FRG computation which can be attributed to the mesonic fluctuations. In both calculations the tricritical point and the two-flavor critical point are continuously connected by the chiral critical line.
Our result without explicit U A (1)-symmetry breaking is shown in the right figure of Fig. 5. The transition is always first-order in the light chiral limit at µ = 0, independent of m K . This is in particular true for the chiral limits in the three-and two-flavor cases, in agreement with [8]. Hence, while we concluded in the last section that fermionic vacuum fluctuations induce a second-order transition in the light chiral limit in absence of the axial anomaly, we observe here that additional meson fluctuations induce a first-order transition again. In the light chiral limit the transition is always of first order, also at finite µ. We find a chiral critical line at vanishing density at m π ∼ 20 MeV for all m K except very small ones, where we see a slight bending towards larger m π .
In Fig. 6 the m K -dependency of the CEP in the (T, µ)plane is shown for c A = 0 in the light chiral limit. Dashed green open symbols correspond to eMFA and solid purple symbols are the FRG results. Starting from the tricritical point at µ = 0 which has the highest critical temperature, the CEP moves towards larger µ and smaller T with increasing kaon mass. The m K -dependency of the CEP is similar to the movement of the CEP when the sigma meson mass is increased as seen in [53]. We summarize our results in the chiral limits and compare with the available literature in Tabs. III and IV. The obtained order of the chiral transition in different chiral limits with and without the axial anomaly is confronted to other works in Tab. III. The important influence of quantum fluctuations and the axial anomaly is evident.
In Tab. IV the critical pion mass m c π at µ = 0 and c A = 0 for three degenerate quark flavors is compared to available lattice results. Even recent lattice results show enormous discrepancies depending on the implementation of fermions. The general trend is that the closer one approaches the continuum limit, by decreasing the lattice spacing or by improving the action, the lighter m c π becomes. It was argued in [101] that converged lattice results for the critical pion mass are expected to be very small. Our continuum results are certainly in favor of a very small critical m c π .

V. SUMMARY
We have studied the mass sensitivity of the chiral phase transition, the left side in the Columbia plot, within a low-energy effective description of QCD. To fully account for non-perturbative corrections to the effective potential caused by quark and meson fluctuations, we applied the functional renormalization group. In addition, we also performed an extended mean-field computation within our framework in order to discern between effects from different fluctuation sources. This allowed us to systematically assess the impact of vacuum and thermal fluctuations of quarks and mesons.
In the presence of the axial U A (1) anomaly in the light chiral limit, fermionic vacuum fluctuations induce a second-order phase transition for larger kaon masses at vanishing density. A rather extensive first-order region around the three-flavor chiral limit emerges. It is separated from a crossover region by a chiral critical line. Meson fluctuations significantly reduce the size of the first-order region. In both the eMFA and the FRG computations the tricritical point at µ = 0 and the two-flavor CEP at finite µ are continuously connected by a chiral critical line. The critical surface bends away from the origin, i.e., the first-order region increases for increasing µ. The physical point always lies deep within the crossover region at vanishing density. Consistent with the curvature of the chiral critical surface a CEP is found at rather large µ. light chiral limit, N f = 3 strange chiral limit, N f = 3 chiral limit, N f = 3 chiral limit, N f = 2 cA > 0 cA = 0 cA > 0 cA = 0 cA > 0 cA = 0 cA > 0 cA = 0 -exp. LSM [8] X X X X 1 st 1 st 2 nd 1 st CJT LSM [47] 2 nd ∀ mK 1 st ∀ mK crossover 1 st ∀ m l 12 2 nd 1 st 2 nd 1 st MFA QM [53] 1  TABLE III. Chiral phase transition order in different chiral limits with and without the axial anomaly from various methods (masses in MeV units). No predictions are labeled as X. The transition is always of second-order above a given upper bound on mK in the light chiral limit and similarly a crossover in the strange chiral limit when an upper bound on mπ is given. In the CJK LSM [47] the Columbia plot was explored in the (m l , ms)-plane and the results with the lightest available sigma mass (600 MeV) has been chosen.
For the U A (1)-symmetric system, we have shown that fermionic vacuum fluctuations, on top of thermal fluctuations, induce a second-order phase transition in the light chiral limit. As a consequence, there is no first-order region in the Columbia plot. Meson fluctuations then turn this second-order transition in the light chiral limit into first order and move the chiral critical line to finite m π values. There is no CEP for any µ for very small pion masses since the transition is always first-order. Thus, the effects of quantum fluctuations on the order of the phase transition drastically depend on the fate of the axial anomaly.
We have also demonstrated the crucial importance of a proper definition of the initial action of the employed low-energy effective model. Simply varying the current quark masses at the initial scale of the effective theory results in the loss of spontaneous chiral symmetry breaking in the chiral limit. This contradicts our current knowledge about gauge theories with light fermion flavors. We argued that the variation of the explicit chiral symmetry breaking parameters leads to sizable modifications of the initial correlation length and therefore concluded that the loss of spontaneous symmetry breaking can indeed be atmethod m c π [MeV] Year standard staggered, Nt = 4 [20] ∼ 290 2001 p4 staggered, Nt = 4 [21] ∼ 67 2004 standard staggered, Nt = 6 [22] ∼ 150 2007 stout staggered, Nt = 6 [23] could be zero 2014 Wilson-clover, Nt = 6, 8 [24] ∼ 300 2014 HISQ staggered, Nt = 6 [25] 50 2017 Wilson-clover, Nt = 8, 10 [26] 170 2017 MFA QM [53] ∼ tributed to an improper choice for the initial action. We proposed the fixed f π -scheme for fixing the initial action for arbitrary values of the current quark masses. Its motivation is drawn from χPT results on the quark mass and flavor dependence of the decay constants. It implies that the initial action has to be modified when the explicit chiral symmetry breaking parameters are varied.
Our results outline the crucial role that fluctuations play for the order of the chiral phase transition. We were able to identify how different sources of quantum fluctuations affect the phase transition. A major drawback of the present analysis is the absence of gluon fluctuations. This can conveniently be taken into account with the FRG by using dynamical hadronization along the lines of, e.g., [56]. The advantage of this approach is that the RG flow can be started in the perturbative regime and the low-energy physics is uniquely fixed by the parameters of microscopic QCD. In addition to the access to the influence of gauge degrees of freedom on the phase transition, it also cures the problem of a heuristic determination of the initial effective action. Another drawback is the implementation of the axial anomaly via a constant coupling to the 't Hooft determinant. In a more realistic scenario, the anomaly coupling c A depends on the medium parameters such as T and µ as well as on the RG scale k.