QCD phase transitions in the light quark chiral limit

We investigate the order of the QCD chiral transition in the limit of vanishing bare up/down quark masses and variations of the bare strange-quark mass $0 \le m_{\mathrm{s}} \le \infty$. In this limit and due to universality long range correlations with the quantum numbers of pseudoscalar and scalar mesons may dominate the physics. In order to study the interplay between the microscopic quark and gluon degrees of freedom and the long range correlations we extend a combination of lattice Yang--Mills theory and a (truncated) version of Dyson--Schwinger equations by also taking back-reactions of mesonic degrees of freedom into account. Both this system and the meson backcoupling approach have been studied extensively in the past but this is the first work in a full $(2 + 1)$-flavor setup. Starting from the physical point, we determine the chiral susceptibilities for decreasing up/down quark masses and find good agreement with both lattice and functional renormalization group results. We then proceed to determine the order of the chiral transition along the left-hand side of the Columbia plot, for chemical potentials in the range $-(30 \,\textrm{MeV})^2 \le \mu_q^2 \le (30 \,\textrm{MeV})^2$. We find a second-order phase transition throughout and no trace of a first-order region in the $N_{f} = 3$ corner of the Columbia plot. This result remains unchanged when an additional Goldstone boson due to a restored axial $\mathrm{U_A}(1)$ is taken into account.


I. INTRODUCTION
It is one of the main goals of the Beam Energy Scan program at RHIC/BNL [1] and the ongoing and future HADES/CBM experiment at GSI/FAIR [2][3][4] to unravel the possible existence and location of a critical endpoint (CEP) in the chiral phase transition line of the QCD phase diagram.The quest of extracting signals for such a CEP from the experimental data is quite delicate and much work is currently being invested to improve the rigorousness of theory-experiment connections (see, e.g., [5,6] for reviews).
On the theoretical side, there is wide-spread consensus on the crossover nature of the chiral transition of QCD at zero chemical potential.The corresponding pseudocritical temperature has been localized around T c ≈ 155 MeV [7,8] with a couple of MeV difference between different definitions of the chiral order parameter.Furthermore, thermodynamic properties of the hot matter in a broad temperature range around T c have been determined with great accuracy [9][10][11][12][13].
Different theoretical approaches to QCD agree with each other that no CEP is found in the region of the temperature-baryon-chemical-potential plane (T, µ B ) with µ B /T < 2.5.This region is excluded by recent studies on the lattice, see, e.g., Refs.[14,15] and references therein, as well as studies using functional methods [16][17][18][19][20][21].Beyond this region, errors in lattice extrapolations accumulate rapidly and no definite statements can be made.On the other hand, functional approaches, i.e., approaches via Dyson-Schwinger equations (DSE) and/or the functional renormalization group (FRG), do in principle allow for a mapping of the whole QCD phase diagram but inherently depend on approximations and truncations necessary to make the equations tractable.
One way to learn more on the behavior of QCD for physical quark masses is to map out the behavior at unphysical up-, down-and strange-quark masses and track structures like critical lines and surfaces at zero, imaginary and real chemical potential.The variation of these reveal areas of different type of transitions, sketched in Fig. 1, the 'Columbia plot' [22].Each of these transitions is related to an underlying symmetry of QCD: chiral symmetry and center symmetry.Their explicit breaking due to nonvanishing (chiral) or non-infinite (center) quark masses generates possible patterns for the order of the transition at finite temperature and vanishing chemical potential as a function of the quark masses displayed in Fig. 1.In the upper right corner of each of the three plots, we find the first order deconfinement transition in the pure gauge limit of infinite quark masses, separated by a second-order critical line from the crossover region.The second-order separation line in the upper right corner of the Columbia plot is in the Z(2) universality class and its location in the u/d-s-quark-mass plane has been mapped out by lattice gauge theory [23][24][25][26][27][28], effective models [29,30], the Dyson-Schwinger approach [31] and background-field techniques [32,33].Thus, although the precise location of the secondorder critical line may differ between the approaches, the qualitative picture is undisputed.This is different for the chiral upper left and lower left corners of the Columbia plot, and the left-hand side of varying strange-quark masses in the light chiral limit.This region is governed by the chiral transition and the corresponding axial symmetries U A (1) × SU A (N f ).Whereas the latter one is broken dynamically at low temperatures (and always explicitly by nonzero quark masses), the for- mer one is broken anomalously.Both the dynamical and anomalous breaking can be restored at large temperatures, albeit the corresponding transition temperatures may very well differ from each other.Whether U A (1) remains broken at the chiral SU A (N f ) transition is an open question with conflicting indications in both directions [34][35][36][37][38][39][40][41][42] The fate of the U A (1) symmetry is expected to affect the order of the chiral SU A (N f ) transition.With an anomalously broken U A (1) at all temperatures, it has been conjectured that the chiral transition for the twoflavor theory (upper left corner) is second order and in the universality class of the O(4) theory, whereas the chiral three-flavor theory (lower left corner) is expected to be first order [43], since no three-dimensional SU(N f ≤ 3) second-order universality class is known [44,45].Consequently, these regions are connected and the left-hand side of the Columbia plot features a tricritical strangequark mass m tri s where the first-order region around the chiral three-flavor point merges into the second-order line connected to the chiral two-flavor point.This is the 'standard' plot seen on the left of Fig. 1.The middle diagram of Fig. 1 shows a possible scenario with restored U A (1).Then the upper left corner may remain first order [43] and the two first-order corners are expected to be connected along the left-hand side of the plot.
It is currently an open question which of these scenarios is realized in QCD.The situation in the upper left corner and, related, in the light-quark chiral limit of the N f = 2+1-theory with strange-quark mass fixed is not clear and indications from lattice simulations vary between favoring either of the two left scenarios of Fig. 1 [40,[46][47][48][49][50][51][52][53].Both scenarios of Fig. 1 can be also realized in effective low energy QCD models such as the PQM or PNJL model, see, e.g., [54][55][56][57][58][59][60][61] and Refs.therein and FRG approaches to QCD [62,63].In Ref. [61], it has been demonstrated that results on the Columbia plot from mean-field approaches are substantially modified once fluctuations have been included.
For the theory with three degenerate flavors, lattice studies support the existence of a first-order transition for light-quark masses on coarse lattices [64][65][66][67][68][69][70].However, the size of the first-order region depends strongly on the formulation of the lattice action and the temporal extend of the lattice and has not yet been determined unambiguously.Thus, it has been conjectured [45] that the third option for the Columbia plot show in the right diagram of Fig. 1 is a realistic possibility.Indeed, recent results on the lattice by Cuteri, Philipsen and Sciarra clearly point in this direction [71] and have been followed up in [72] with similar results.In Ref. [73], it has been suggested that a second-order N f = 3 transition may not be at odds with previous FRG results, see [61] and references therein.
It is the purpose of this work to re-examine the situation in a functional continuum framework that takes both, microscopic quark and gluon degrees of freedom and effective, long-range degrees of freedom with the quantum numbers of pseudoscalar and scalar mesons, into account.In the framework of Dyson-Schwinger equations (DSE) that we employ, these appear naturally as part of fermion four-point functions in the DSE for the quark-gluon vertex [74,75].At nonzero temperature and chemical potential, the corresponding framework has been already explored for physical quark masses [21,76,77] and has led to a prediction of the location of the critical endpoint in agreement with recent FRG studies [18][19][20].Here, quarks have been taken into account on the N f = 2 + 1 level but the meson sector remained N f = 2 [21,76,77].In this work, we extend the framework to a consistent N f = 2 + 1 level and therefore make it suitable for a study of the left-hand side of the Columbia plot.
The paper is organized as follows.In the next section II, we detail the framework of Dyson-Schwinger equations including the quark and gluon DSEs as well as the abovementioned fluctuation effects on the quark-gluon vertex.
2. General form of the DSE for the quark propagator (top) and truncated gluon DSE (bottom).Large grey and white dots indicate dressed quantities; solid and curly lines represent quark and gluon propagators, respectively.There is a separate quark DSE for the up, down and for the strange quarks.The large shaded dot denotes the quenched gluon propagator that is taken from the lattice while the quark loop is evaluated explicitly.The latter contains an implicit flavor sum over up, down and strange.
We discuss our treatment of the corresponding meson masses and decay constants for varying strange-quark mass, taking particularly care of the limits m s → ∞ and m s → 0. In Section III, we then present our results for the order of the phase transition along the left hand side (i.e., chiral up/down quarks but varying strange-quark mass) of the Columbia plot for zero and small real and imaginary chemical potential.We discuss critical temperatures, the resulting Columbia plot and the dependence of our result on the restoration temperature of the U A (1) symmetry.We conclude in Section IV.

A. Dyson-Schwinger equations
All necessary quantities for our investigation of the Columbia plot can be obtained directly from dressed (i.e., full) quark propagator S f .For a given quark flavor f ∈ {u, d, s}, its inverse at nonzero temperature T and quark chemical potential µ f is given by (1) Here, p = (p, ωn ) represents the four-momentum, while ωf n = ω n + iµ f denotes a combination of the fermionic Matsubara frequencies ω n = (2n + 1)πT , n ∈ Z, with the chemical potential.All non-perturbative information such as the non-trivial momentum dependence is carried by the quark dressing functions A f , B f and C f .
To calculate the quark propagator, we solve its associated Dyson-Schwinger equation (DSE) which reads (2) Above, m f denotes the current-quark mass while Z 2 and Z m labeling the wave function and mass renormalization constants, respectively, which are calculated in vacuum using a momentum-subtraction scheme.
The quark DSE is displayed pictorially in the top row of Fig. 2. The quark self-energy Σ f comprises both the gluon propagator and the quark-gluon vertex.In our framework, we calculate the gluon propagator explicitly by solving its DSE albeit in a truncated form (for details see below).For the quark-gluon vertex we rely on a welltested vertex model constructed to (approximately) satisfy Slavnov-Taylor identities and preserving all perturbative and renormalization constraints (see below) [16,17,[78][79][80].
As the correlation length diverges in the vicinity of a second-order phase transition, long-range correlations in the quark-gluon vertex become important.These arise from a specific diagram in the DSE for the quark-gluon vertex that involves a four-quark kernel.In pole approximation, this diagram is shown in the upper equation of Fig. 3.This diagram provides contributions to all tensor components of the quark-gluon vertex [74].In the quark DSE, the resulting two-loop diagram can be simplified to a one-loop diagram using a homogenous BSE, see Refs.[21,74] for details.The effect of this specific contribution to the quark-gluon interaction has been studied in a number of works at zero temperature/chemical potential including a discussion of the analytic structure of the quark propagator [75], a discussion of its effect onto the meson spectrum [81], and an exploratory study of meson-cloud effects in baryons [82].In all these studies, it has been noted that meson-backcoupling effects typically provide contributions of the order of 10-20 % as compared with other components of the quark-gluon interaction.The effect of this contribution on the location of the CEP has been studied in Ref. [21].
The quark-gluon vertex, split into a non-hadronic (NH) and a mesonic (M) part and inserted into the quark-DSE leads to the following expression for the analogous splitting of the quark self-energy: The non-hadronic part of the quark self-energy corresponds to the usual quark self-energy from previous works: Here, k = p − q indicates the gluon momentum, g labels the strong coupling constant, Z3 represents the ghost renormalization constant and D νρ is the dressed gluon propagator.The prefactor of 4/3 originates in the trace over color space which has already been carried out for The new element of the truncation used in this work as compared to Ref. [21] is the inclusion of the strangequark contributions to the mesonic part of the vertex.The necessary technical steps are discussed in detail in Subsection II B below.
In Eq. ( 4), Γ f ρ labels the non-hadronic quark-gluon vertex for which we employ the ansatz: This is a combination of the leading Dirac tensor structure of the Ball-Chiu vertex [83] -which leads to backcoupling effects of the quarks onto the vertex -with a phenomenological vertex dressing function: The first term is an IR enhancement inspired by Slavnov-Taylor identities while the second term ensures the correct perturbative behavior of the propagators in the UV.The running coupling is given by α = 0.3, the scales d 2 = 0.5 GeV 2 and Λ = 1.4 GeV are fixed to match the ones in the gluon lattice data.The anomalous dimension reads For the argument, we have x = k 2 in the quark DSE while x = p 2 +q 2 in the gluon DSE in Eq. ( 8) to ensure multiplicative renormalizability.More details on the vertex and the choice of the parameters can be found in Refs.[16,17] and references therein.Since the meson-backcoupling diagrams originate in a modification of the vertex, we need to adjust the vertex-strength parameter d 1 .We tune d 1 such that the pseudocritical temperature obtained from the chiral susceptibility at the physical point, i.e., the maximum of m π = 139 MeV in Fig. 6, corresponds to the one from the lattice T p c = 156.5 MeV.This yields d 1 = 8.98 GeV 2 as opposed to d 1 = 8.49 GeV 2 without meson contributions [17].
The gluonic part of our truncation is unchanged as compared to previous works, i.e., in the Yang-Mills sector we do not take the meson effects explicitly into account 1 .As stated earlier, the gluon propagator is calculated using 1 The main reason is feasibility: in the quark-loop diagram of a simplified version of the full gluon DSE (illustrated in the bottom row of Fig. 2): Here, D YM νρ denotes the quenched gluon propagator given by a combination of all pure Yang-Mills for which we use temperature-dependent fits to results of quenched lattice calculations [79,85,86] as input.The quark loop Π νρ , however, is calculated explicitly within our framework.This leads to an unquenching of the gluon and consequently a non-trivial backcoupling of the chiral dynamics of the quarks onto the gluon.It is given by The prefactor of 1/2 stems from the color trace, the flavor sum f ∈ {u/d, s} runs over the investigated 2 + 1 quark flavors.We work in the isospin-symmetric limit of degenerate up and down quarks (m u = m d , µ u = µ d ).Therefore, we will from now on only refer to a 'light' quark ℓ = u = d for the sake of simplicity.At the physical point, the quark masses are fixed using results for the pion and kaon masses in vacuum obtained from the Bethe-Salpeter formalism developed in Ref. [87].This leads to values of m ℓ = 0.8 MeV and m s = 20.56MeV at a renormalization point of 80 GeV.

B. Meson-backcoupling self-energy
In this subsection, we discuss the mesonic part of the quark self-energy in some detail.To this end, we begin again with the two-loop expression in the quark DSE in the lower equation of Fig. 3.This diagram originates from the gluon-DSE the meson-exchange diagram remains two-loop and is therefore too expensive in terms of CPU time.However, these contributions are also irrelevant when it comes to critical exponents [84].the meson pole approximation of a fermion four-point function in the DSE for the quark-gluon vertex [74].The corresponding meson propagator in the diagram is therefore bare and accompanied by two Bethe-Salpeter amplitudes that connect the quark lines with the exchanged meson in question.In Ref. [74], it has been realized that the left half of the two-loop diagram displayed in Fig. 3 can be interpreted as the interaction diagram in a homogeneous Bethe-Salpeter equation (BSE) and therefore can be replaced with a Bethe-Salpeter amplitude (BSA).This way, the mesonic part of the quark self-energy reduces to a one-loop diagram illustrated in Fig. 4. Of course, this simplifies calculations tremendously.In principle, this diagram contains mesons with all quantum numbers that can be build from a quark-antiquark pair.In practice, we are only interested in those mesons that have the potential to become massless at phase transitions, i.e., the pseudoscalar meson octet, its critical chiral partner modes and the pseudoscalar singlet in case the axial U A (1) is restored at the transition temperature.All these are potentially long ranged and are expected to become the dominant degrees of freedom at second-order phase transitions.All other meson contributions are subleading due to their large masses in the meson propagator and are therefore omitted in our approach.
We thus end up with the lightest pseudoscalar octet, i.e., pions, kaons and the η 8 , as well as the η 0 in a crosscheck calculation (see Section III).Additionally, we consider the scalar σ meson (i.e., the f 0 (500)) as it is vital for the correct O(4)-scaling behavior in the upper left corner of the Columbia plot and the ss-partner of the σ (which we identify with the f 0 (980)) that may be important in the N f = 3 chiral limit.As detailed below, Eq. ( 11), we assume the f 0 to be massless in the chiral N f = 3 limit, since it has the quantum numbers of the strange-quark condensate.In order to obtain a consistent N f = 3 limit, we alter its flavor factor by hand to match the one of the σ meson thus obtaining three identical DSEs for the up, the down and the strange quark in this limit (cf.Tab.I). 2estricting to the pions and σ meson, this type of meson backcoupling was discussed in detail, e.g., in Refs.[21,84].
Building on the explanations therein, we generalize this to the N f = 2 + 1 case to arrive at the following mesonic contribution to the quark self-energy: Here, P = p − q denotes the meson's total momentum while l i represent the relative momenta of the Bethe-Salpeter amplitudes.Using an appropriate momentum routing, we can identify l 1 = p and l 2 = q.
Most quantities in Eq. ( 9) are specific to the flavor of the external quark f ∈ {ℓ, s} and the exchanged meson X ∈ {π, K, η 8 , σ, f 0 , (η 0 )}.First, there are the multiplicities F f X of the respective meson-backcoupling diagram obtained by performing the trace over flavor space.Second, we have the quark propagator of the internal quark S f X which differs from the external one for the kaon, i.e., S ℓ K = S s and S s K = S ℓ .The (inverse) meson propagator at nonzero temperature is given by [88]: with u X = f s X /f t X being the meson velocity which is given by the ratio of the spatial and temporal meson decay constants, f s X and f t X , respectively.Again, we restrict ourselves to potentially critical modes and consider only the zeroth Matsubara frequency of the meson propagator, i.e., P 4 = 0.All other Matsubara frequencies act as an effective meson mass that leads to suppression of the respective contribution.This restriction implies ω q = ω p in Eq. ( 9) and consequently cancels the Matsubara sum [89].
For the meson masses m X , we choose the vacuum values that have been obtained in the following way: We have solved the coupled system of N f = 2 + 1 DSEs without meson-backreaction effects in the vacuum and for complex quark momenta and then extracted the corresponding meson masses from solving their corresponding Bethe-Salpeter equations for different up/down-and strangequark masses.The resulting mass curves for the pion and the kaon have been fitted with the expressions above, which correspond to the expected behavior from Gell-Mann-Oakes-Renner relations.The remaining masses dynamically.In this case, adjusting flavor factors by hand would not be necessary.
are expressed in terms of these for the sake of simplicity in such a fashion that the correct massless modes appear in the N f = 2 and N f = 3 chiral flavor limits.Note that this treatment overestimates the effects of the critical modes in the low-temperature, chirally broken phase since the critical modes are already massless by construction instead of becoming massless at the critical temperature.We have checked that this simplification does not affect the order of the transition, but it may affect the location of the transition, i.e., the critical temperature.We discuss this further in Section III, when we present our results.In principle, one could solve the temperature-dependent Bethe-Salpeter equations also at nonzero temperature along the lines of Ref. [89].There, it has been shown explicitly that the pion and sigma modes follow the correct pattern of symmetry breaking and restoration in the N f = 2 chiral limit.In practice, this would add an extra layer of complications and an order of magnitude more in computing time to an already demanding endeavor and we there resort to the simplifications expressed in Eq. (11).
The central unknown quantities are the meson Bethe-Salpeter amplitudes ΓX .In the chiral limit, it is an exact property of QCD [90,91] that the leading BSA of the Goldstone boson can be expressed through the scalar dressing function B of the quark propagator and the Goldstone-boson decay constant via Γ X (l) = γ 5 B(l)/f X , with relative momentum l between the quark and the antiquark, see [92] for a review and a detailed explanation of this property.This behavior persists approximately also away from the chiral limit with the caveat that the quark dressing function then develops a logarithmic tail at large momenta, whereas the Bethe-Salpeter amplitude always falls like a power in momentum.We therefore adopt the following prescription for our meson amplitudes: where f f,t X labels the respective temporal meson decay constant and γ X = γ 5 for the pseudoscalar mesons and γ X = I for the scalar mesons.Additionally, we also apply these relations to mesons comprising non-chiral strange quarks.To account for the correct power-law behavior in the large momentum limit, we supplement the quark dressing function with a Pauli-Villars-like term with a scale that matches our renormalization point [93].As a consequence, this also renders the meson-backcoupling diagrams ultraviolet finite so that no further regularization is necessary.Note that for mesons with mixed flavor content we always use the B function of the quark external to the loop in which the BSA appears.This turned out to be numerically advantageous for a consistent N f = 3 limit.
Additionally, we introduce the closely related quantity ΓX which originates in the two-loop diagram of the vertex expansion.Apart from a non-trivial sign which is a result 5. Pictorial representation of the generalized Pagels-Stokar formula in Eq. ( 14) we use to calculate the pseudoscalar and scalar meson decay constants.
of its two-loop origin, we identify it with the BSA ΓX : where (−1) X = −1 for the pseudoscalar mesons and (−1) X = +1 for the scalar mesons.Last, we also require the meson decay constants.These are calculated using a generalized Pagels-Stokar formula [94], (14) which merits some explanations.First, we use the abbreviation Pµ = (ω P , u X P ), while x and y label the flavor indices of the quarks contributing to the meson in question.Consequently, the temporal and spatial decay constants are obtained in the limit P µ → 0 from the temporal and spatial momentum component, respectively.The index z labels the quark flavor of the external quark of the backcoupling diagram where the meson appears and Y ∈ {ps, sc} represents the parity of the meson.This relation is pictorially displayed in Fig. 5.This notation is necessary for the following reason: the decay constants in our diagrams do not depend on the mesons directly but rather on the contributing quark propagators of the backcoupling diagram.Furthermore, they have to match the type of scalar dressing function used in the BSA.For hidden-flavor mesons this is unproblematic and Eq. ( 14) matches the usual Pagels-Stokar approximations.For open-flavor mesons, however, the four kaons in our approach, this has the consequence that the decay constant appearing in the light-quark DSE is different from the one in the strange-quark DSE since the external quark is different.Eq. ( 14) accounts for this.Furthermore, we symmetrize Eq. ( 14) in a mathematically well-defined manner 3 with the arithmetic mean of the exchanged quark flavors.This way, we arrive at the following definitions: Information of multiplicities, internal quark propagators and decay constants for all considered meson backcoupling diagrams.
from which we can obtain the required (temporal) meson decay constants used in our calculations: Of course, this procedure is only relevant for finite, nonzero strange-quark masses along the left-hand side of the Columbia plot.In the chiral N f = 2 and N f = 3 limits (lower and upper left corners), it becomes immaterial.In total, we summarize all necessary information for our meson-backcoupling procedure compactly in Tab.I.

III. RESULTS AND DISCUSSION
In this section, we present the numerical results obtained in our framework.Our investigation of the Columbia plot is based on monitoring the behavior of the light-quark condensate as the order parameter for chiral symmetry breaking.It can be obtained from the quark propagator via the relation where the factor three stems from the color trace and, as always, f ∈ {ℓ, s}.The quark condensate is quadratically divergent for all flavors with a nonzero bare-quark mass due to a contribution proportional to m f Λ 2 where Λ denotes the ultraviolet (UV) momentum cutoff.As a consequence, the divergent part of the light-quark condensate can be cancelled with the corresponding one of the strangequark condensate when the latter is multiplied with the light-to-strange mass ratio.In our (2 + 1)-flavor setup, a regularized expression for the light-quark condensate can therefore be obtained by considering the difference Note that we are working with renormalized quantities, hence the appearance of the renormalization constants Z f m in order to preserve multiplicative renormalizability.In the case of massless light quarks, the subtracted condensate reduces to the unsubtracted one which is then well-defined, i.e., UV finite.
The chiral susceptibility is then defined as the derivative of the regularized condensate with respect to the lightquark mass: Up to normalization factors and for m ℓ = m u = m d , this definition is equivalent to the ones used in Refs.[63,95].The remainder of this section is structured as follows.First, we study the line between the physical point and the left-hand side of the Columbia plot by analyzing the dependence of the chiral susceptibility on the pion (and thus the light-quark) mass.Second, we investigate the type of the chiral phase transition across the whole lefthand side of the Columbia plot, i.e., for chiral light quarks and strange-quark masses between m s ∈ [0, ∞).We also quantify the dependence of the critical temperature on the strange-quark mass.Third, we analyze the scaling behavior of the light-quark condensate.Last, we extend our analysis to small but nonzero chemical potential, both real and imaginary.

A. Towards the chiral limit
We start our investigation of the Columbia plot with the line between the physical point and the left-hand side.That is, we keep the strange-quark mass physical m s = m p s and decrease the light-quark mass from its physical value down to zero.This path has been explored already by the HotQCD collaboration [95] with lattice QCD methods, the fQCD collaboration using the FRG [63] as well as with a combined FRG-DSE approach in Ref. [96].
In the left diagram of Fig. 6, we show the chiral susceptibility as a function of temperature for four different pion masses compared to the lattice results of Ref. [95].Analogously to Ref. [63], we normalize the susceptibilities to the maximal value at a physical pion mass: Qualitatively, we find similar results as both the lattice, the FRG and the FRG-DSE approach in Refs.[63,95,96].That is, for decreasing pion masses, the peak of the susceptibilities increases in height and moves towards lower temperatures monotonically.Quantitatively, the results are also very similar.Namely, the peak heights are comparable with the FRG and FRG-DSE findings for all investigated pion masses and the (pseudo)critical temperatures away from the chiral limit do not deviate more than 2 MeV and 3.5 MeV, respectively, see Tab.II.It is, however, interesting to note that the linear extrapolation of our results towards the chiral limit, shown in right diagram of Fig. 6  temperature calculated explicitly in the chiral light-quark limit by about two MeV.The linear extrapolation results in 145. 4 MeV, whereas we find the value which is about 5 MeV-5.5 MeV larger than the FRG and FRG-DSE results, respectively, and more than ten MeV larger than the extrapolated lattice value T HotQCD c = 132 +3 −6 MeV.We attribute this to our vertex construction, which contains a strength parameter d 1 which is fixed at the physical point and not changed with quark mass.We therefore slightly overestimate the interaction strength in the chiral limit leading to slightly too large transition temperatures.This will be discussed again also in the next section.

B. Left edge of the Columbia plot
Next, we turn our analysis to the left edge of the Columbia plot.To this end, we display the temperature dependence of the quark condensate for chiral light quarks and a set of six selected strange-quark masses between m s ∈ [0, ∞) in the left diagram of Fig. 7.One can immediately notice that for all investigated strangequark masses we observe a second-order phase transition.
That is, the quark condensate continuously changes from a nonzero value to zero with increasing temperature with no (apparent) jumps.As we will see in Sec.III C, it is indeed a genuine second-order transition since the condensate exhibits a scaling behavior in the vicinity of the respective critical temperatures.We emphasize that this also holds true for the N f = 3 corner where we consequently find no first-order region at all.In general, the condensate is smaller for smaller strange-quark masses at all temperatures.The only exception occurs close to the three-flavor limit at around m s ∼ 10 −9 MeV where we do find a sudden and small increase in the condensate for all temperatures which then remains constant until m s → 0. In the left diagram of Fig. 7, this increase is visible when comparing the m s = 10 −3 MeV result with the the one at m s = 0 MeV.Since we neither found a technical nor a physical reason for this glitch we attribute it (for the moment) to a pure numerical artefact of the three-flavor limit.
The dependence of the critical temperature on the strange-quark mass is illustrated in the right diagram of Fig. 7. Qualitatively, we find that T c varies very little for very small and very high strange-quark masses but increases monotonously in the range m s = 1 MeV-10 4 MeV.
In Tab.III, we compare our findings of the critical temperature quantitatively to the most-recent lattice re- sults available for zero, physical and infinite strange-quark masses.We observe that our values for T c are consistently larger than the ones found on the lattice, with the smallest difference at the physical strange-quark mass.Our explanation for this discrepancy is based on the discussion above: we fix the interaction strength d 1 for the non-hadronic part of our quark-gluon interaction, Eq. ( 6), at the physical point and do not take into account any changes of the vertex strength with variation of the quark masses.Presumably, this leads to the small discrepancy in transition temperature in the light chiral limit with physical strange quarks already discussed above and larger discrepancies in the chiral corners of the Columbia plot.We have explicitly checked what happens in the N f = 3 limit, when we adapt the vertex strength.With d 1 = 7.13 GeV 2 we reproduce the transition temperature of Ref. [72] while the transition is still second order.Thus, the value of d 1 (at least within the range studied here) had no material influence on the order of the transition.
Finally, we investigated whether the fate of the U A (1) symmetry has any influence on the order of the transition.On the complexity level of the present truncation, an anomalously broken U A (1) at the critical temperature results in a massive η ′ meson, whereas a restored U A (1) renders the η ′ meson massless.So far, we assumed the first case and neglected the η ′ meson in the backreaction diagrams completely together with all other mesons that remain massive in the chiral N f = 3 limit.In order to gauge the influence of the η ′ on the order of the transition in this limit, we repeated our calculation with the massive η ′ and, even more importantly, with a massless η ′ explicitly present in the loops.As a result, we find only a very small reduction of the transition temperature of about 0.3 MeV when including a massive η ′ with no changes in the second-order nature of the transition.This result confirms our notion that additional massive mesons barely have any influence on our results and their omission, therefore, is a good approximation.Including a Adapting our vertex strength parameter to match the lattice critical temperature in the ms = 0-limit does not change the second-order nature of the chiral phase transition.
massless η ′ reduces the transition further by about 1.5 MeV but again does not change the second-order nature of the phase transition.We therefore conclude that within the framework presented here the fate of the anomalously broken axial U A (1) symmetry with temperature has no material effect on the order of the chiral phase transition.

C. Scaling behavior
In the vicinity of a second-order phase transition the order parameter, obeys a power law with respect to some (universal) reduced quantity.For the regularized condensate, we expect the following behavior: labels the reduced temperature, β indicates the critical exponent depending on the underlying universality class and c denotes a non-universal constant.The scaling properties of the quark DSE and, related, that of the condensate have been studied in Ref. [84] in the chiral N f = 2 limit.
Here we expect a second-order phase transition in the O(4) universality class of the Heisenberg antiferromagnet, see, e.g., Ref. [43].Indeed, it has been shown in Ref. [ scaling of the temporal meson decay constants is taken into account.Therefore, to obtain self-consistent scaling from the quark DSE, one would need to explicitly include the BSE for the relevant long-range degrees of freedom, the pion and the sigma, as well as the corresponding equations for their decay constants in a self-consistent manner.This is beyond the present truncation and would require extensive additional work.A shortcut, also used in Ref. [84], is to assume the critical scaling of the decay constants and only check for consistency in the quark DSE using an appropriate ansatz for this scaling.Indeed, one then finds that the diagrams with long-range massless meson degrees of freedom dominate over the gluon ones, i.e., universality sets in and scaling is observed.On the other hand, any setup without scaling decay constants delivers the usual mean-field scaling observed in rainbowladder type truncations already very early and reviewed in Ref. [98].
For completeness, we have checked both, mean field scaling without and O(4) scaling with proper scaling decay constants.For the former we use the Pagels-Stokar approximation for the decay constants discussed in Eq. ( 14), for the latter we use the following ansatz4 Here, we use the decay constants from the Pagels-Stokar setup at some temperature T 0 = 100 MeV as input, while the critical temperatures T c are the ones from Fig. 7.The scaling law for f Y ℓℓ is taken from Ref. [84], whereas the scaling law for f ps ℓs , valid for m s ̸ = 0 is obtained from an analogous scaling analysis as in Ref. [84] for the kaon diagrams.All other decay constants do not exhibit any critical scaling.In the limit of m s → 0, we just assume the same critical scaling for f Y ℓℓ , f ps ℓs and f ps ss .This is, however, for simplicity only since the universality class in this limit is not known.
In Fig. 8, we display our results.In the left diagram, we show the (logarithm of the) regularized condensates in the chiral light-quark limit and for the same strange-quark masses as in Fig. 7 as functions of the (logarithm of the) reduced temperature t.For the sake of comparability, we fit each dataset to the relation in Eq. ( 23) and divide by their respective constant c.As can be seen, all curves collapse nicely onto each other and align along the line t 0.5 for log(t) ≤ −1.5.The spread of data points below log(t) ≤ −2 is entirely due to our numerical uncertainty of ∆T c = 0.1 MeV in the determination of T c .We observe the expected mean-field scaling behavior for all investigated strange-quark masses.
In the right diagram of Fig. 8, we show corresponding results under the assumption that the decay constants scale according to Eq. ( 24) with the correct O(4) critical exponents of QCD, where β = 0.73/2.As expected, the scaling behavior of the decay constants induces a consistent corresponding scaling of the order parameter.Furthermore, with induced scaling, the collapse of all datasets for log(t) ≤ −1.5 is almost perfect and the spread for low log(t) vanishes completely.This is due to the fact that the additional appearance of T c in the scaling ansatz stabilizes the numerics considerably.
Of course, since β = 0.73/2 is an input, this setup reveals nothing about the universality class in the chiral three-flavor limit.In order to study this issue from DSEs, as discussed above, one needs to solve the corresponding Bethe-Salpeter equations and the defining equations for the decay constants without further approximations.This is possible, in principle, and should be attempted in future work.In practice, it may however be more straightforward to perform this calculation in the framework of the functional renormalization group, where scaling properties are more directly approachable [57,61,63,[99][100][101][102].

D. Three-dimensional Columbia plot
Finally, we would like to explore the fate of the secondorder phase transition along the left-hand side of the Columbia plot when we switch on chemical potential.Is there a second-order critical sheet?And if yes, does it bend at some point towards nonzero quark masses and is it connected to the CEP that we find at physical quark masses [17,21,78]?
Unfortunately, these questions are difficult to study.Our current approximation for the meson Bethe-Salpeter amplitudes, Eq. ( 12), is known to be accurate at vanishing chemical potential.From the explicit calculation in Ref. [76,77], however, it is known that the amplitudes are modified substantially at large chemical potential.We can therefore only trust the approximation Eq. ( 12) at small chemical potentials.
We therefore restrict ourselves to real and imaginary baryon chemical potentials of |µ B | = 30 MeV.The corresponding results are shown in Fig. 9 together with results for µ B = 0 as a comparison.The calculations have been performed for m s = 0, m p s , ∞.In the top panel, we display the condensate as a function of temperature in the N f = 3 chiral limit, i.e., m ℓ = m s = 0. We find no significant changes within this range of chemical potential.Similar results are obtained for all investigated strangequark masses in m s 0, m p s , ∞.In total, we find little quantitative and no qualitative difference between the results for vanishing and small chemical potentials.That is, one can still observe a second-order phase transition with identical scaling behavior and an almost unchanged critical temperature.We therefore arrive at the slice of the three-dimensional Columbia plot shown in the bottom panel of Fig. 9.For zero chemical potential, this ties in with the lattice results of Ref. [71,72] and for imaginary chemical potential with Ref. [103].It also agrees with one of the scenarios displayed in the FRG approach in Ref. [61] (their right diagram of Fig. 3), but disagrees with the other scenarios they give.This needs to be reexamined in some detail.In any case, the analyticity of the second-order transition plane from small imaginary to small real chemical potential visible in Fig. 9 is, to our knowledge, shown for the first time.
Last, we perform a first exploration of the 'bottom' line of the Columbia plot, i.e., the zero chemical-potential line with a massless strange quark and up/down-quark masses in the range m ℓ = 0, ..., m p ℓ , ..., ∞.Then, the strangequark condensate becomes the corresponding order parameter for the chiral transition at finite temperature, so we include one massless Goldstone boson due to the dynamical breaking of a U A (1) subgroup of the flavor SU A (3) and we expect the isoscalar scalar f 0 with sscontent to be the only additional massless mode at the critical temperature.Furthermore, we assume a restored anomaly at the critical temperature.The corresponding meson masses, flavor coefficients and decay constants are detailed in Appendix A. For this setup, we indeed find again a second-order phase transition, also indicated in our three-dimensional Columbia plot in the bottom panel of Fig. 9.This second-order transition persists to large but not infinite up/down-quark masses.At some large up/down-quark mass, m ℓ > 100 GeV, our calculation breaks down indicating that we are approaching the one-flavor limit of QCD with different symmetries.The detailed study of this limit is non-trivial and deferred to future work.

IV. SUMMARY AND CONCLUSIONS
In this work, we studied the order of the phase transition in the light chiral limit of massless up/down quarks as a function of the mass of the strange quark and at zero and small values for the baryon chemical potential.Using a truncation of Dyson-Schwinger equations that takes into account microscopic degrees of freedom as well as potential long-range correlations with the quantum numbers of pseudoscalar and scalar mesons, we obtain a chiral crossover as long as the light-quark masses remain nonzero, but a second-order phase transition in the light chiral limit.This behavior persists along the left-hand side of the Columbia plot, i.e., for all strange-quark masses in 0 ≤ m s ≤ ∞ and also for (small) imaginary and real chemical potential.It persists regardless whether we include a massive η ′ meson (in case the axial U A (1) remains anomalously broken at T c ) or a massless η ′ meson (in case the axial U A (1) is restored at T c ).Our findings do not support the long-standing notion of a chiral first-order N f = 3 corner in the Columbia plot [43], but agree with recent findings from lattice QCD [71,72] and notions from effective models [73].
N f = 1 + 2 → 3 is consistent with the corresponding limit N f = 2 + 1 → 3 in Tab.IV under inclusion of the η 0 .In the scalar meson sector, a massless isoscalar ss (which we  call f 0 ) at the chiral transition temperature reflects the fact that the strange-quark condensate is the appropriate order parameter.The kaons are massive away from the N f = 3 limit, similar to the N f = 2 + 1 case, but this time due to the non-chiral up/down quarks.In total, we reuse the parametrizations of Eq. ( 11)) but adjust the quark masses in the kaon argument and set the f 0 and η s masses to zero:

FIG. 1 .
FIG.1.Three different versions of the Columbia plot[22] of phase-transition orders at nonzero temperature and vanishing chemical potential as functions of quark masses.The 'standard plot' with anomalously broken UA(1)-symmetry is shown on the left, in the middle we display a possible version with restored UA(1) and on the right we show an alternative version without chiral first order regions.Here, we assume mass-degenerate up and down quarks, m ℓ = mu = m d .

FIG. 3 .
FIG. 3. Top: Vertex ansatz for including long-range correlations originating the in the skeleton expansion of the quark-gluonvertex DSE.Bottom: Resulting DSE for the quark propagator.

1 l 2 FIG. 4 .
FIG. 4. One-loop meson-backcoupling diagram in the quark self-energy as an approximation of the two-loop diagram of Fig. 3.

15 −∆ ℓs [GeV 3 ] 1 FIG. 7 .
FIG.7.Left: Regularized quark condensate as a function of temperature for different strange-quark masses in the up/down quark chiral limit.Right: Dependence of the critical temperature on the strange-quark mass in the same limit.The dashed vertical line indicates the physical strange-quark mass.

FIG. 9 .
FIG.9.Regularized quark condensate for (small) chemical potentials in the up-and strange-quark chiral limits.Top: Regularized quark condensate as a function of temperature.Bottom: Illustration of the three-dimensional Columbia plot we find.
m π = 156.525MeV 1/2 • √ m ℓ , m σ = 2m π , m K = 74.2MeV 1/2 • √ m ℓ + 1.54 • m ℓ , m η ℓ = 2m K , m ηs = m f0 = 0 .(A1) [95]6.Left: Chiral susceptibility as a function of temperature for fixed, physical strange-quark mass but different up/down-quark masses corresponding to different pion masses.The susceptibilities are normalized to the maximal value of mπ = 140 MeV.The lattice data have been taken from Ref.[95].Right: Critical temperatures for the same pion masses.Shown in a linear extrapolation of our data at finite pion masses to zero mass compared with the result of an explicit calculation at zero pion mass.

TABLE II
. Comparison of critical temperatures for different up/down-quark masses corresponding to different pion masses and fixed physical strange-quark masses between our DSE findings, the FRG, FRG-DSE and the lattice results, respectively.

TABLE III .
Comparison of critical temperatures for different strange-quark masses between our DSE findings and lattice results.
Scaling behavior of the regularized quark condensate as a function of the reduced temperature for different strange-quark masses in the up quark chiral limit without scaling decay constants.Right: Same scaling behavior with scaling decay constants (see main text for discussion).

TABLE IV .
Information of multiplicities, internal quark propagators and decay constants for all considered meson backcoupling diagrams in the N f = 1 + 2 setup.