Localized rainbows in the QCD phase diagram

We study the phase diagram of strongly interacting matter with light quarks using a recently proposed, small parameter approach to infrared QCD in the Landau gauge. This is based on an expansion with respect to both the inverse number of colors and the pure Yang-Mills coupling in the presence of a Curci-Ferrari mass term. At leading order, this leads to the well-known rainbow equation for the quark propagator with a massive gluon propagator. We solve the latter at nonzero temperature and chemical potential using a simple semi-analytic approximation known to capture the essence of chiral symmetry breaking in the vacuum. In the chiral limit, we find a tricritical point which becomes a critical endpoint in the presence of a nonzero bare quark mass, in agreement with the results of nonperturbative functional methods and model calculations. This supports the view that the present approach allows for a systematic study of the QCD phase diagram in a controlled expansion scheme.


I. INTRODUCTION
Hadronic matter is expected to present a rich phase structure when submitted to sufficiently high energy and baryonic densities, large magnetic fields, etc., as encountered in various environments such as the early universe, ultradense astrophysical objects, or relativistic heavy ion collisions in the laboratory [1][2][3]. Unravelling the phase diagram of quantum chromodynamics (QCD) at nonzero temperature T and baryonic chemical potential µ B is a major challenge both experimentally and theoretically. At vanishing chemical potential, first principle lattice simulations unambiguously demonstrate a smooth crossover from a mostly confined to a mostly deconfined phase, accompanied by a restoration of chiral symmetry [4,5]. The crossover region sharpens for increasing quark masses and turns in a second order phase transition for critical values of the quarks masses, above which the transition is first order. The same is expected to happen with the chiral transition in the opposite limit of decreasing quark masses: the transition turns first order below a critical value of the quark masses. Although not firmly established by lattice calculations yet [6], this is the expected behavior of a theory with at least three light quark flavors. The situation with two light quarks is more subtle due to the possible role of the axial anomaly [7].
The situation is even less clear at µ B = 0 in the low quark mass region (including the physical point), where standard Monte Carlo algorithms are plagued by the infamous sign problem [8]. The typical expectation is that of a line of first order chiral transition at low temperatures ending at a critical point [9]. Firmly establishing the existence of the latter and studying its possible experimental signatures has been the topic of intense theoretical work [10][11][12][13][14][15][16] and is among the major physics goals of various present and upcoming experiments [17][18][19][20]. Methods to circumvent the sign problem on the lattice have been devised but remain, so far, limited to µ B /T 1, and no critical endpoint has been firmly established [21]. To go beyond, a fruitful proposal has been to supplement lattice results for the quenched gluon dynamics with explicit quark contributions by means of nonperturbative functional methods [22][23][24]. Existing works neglect the mesonic degrees of freedom-although baryons have been included [24]-and find a critical endpoint at a relatively large µ B /T 3. A complementary approach uses phenomenological, Nambu-Jona-Lasinio (NJL) or quarkmeson models with various degrees of sophistication [25][26][27][28][29]. These typically predict a critical endpoint at relatively large µ B /T , whose precise location, however, varies significantly from one study to another. One common weakness of such approaches is that the employed approximations lack a systematic ordering principle. One typically explores the whole phase diagram with truncations adjusted against lattice data at µ B = 0, far from the region where a critical point is found.
In the present article, we study this question using a semi-perturbative approach to the infrared dynamics of QCD based on a simple massive extension of the Faddeev-Popov (FP) Lagrangian in the Landau gauge, known as the Curci-Ferrari (CF) model [30,31]. In this context, the gluon mass term is motivated both by the results of lattice simulations [32] and by the necessity to modify the FP Lagrangian in the infrared due to Gribov ambiguities [33]. The CF model is the simplest renormalizable deformation of the FP Lagrangian and remains under perturbative control down to the deep infrared: The gluon mass screens the standard perturbative Landau pole and the (running) gauge coupling remains moderate at all scales [31,34], as observed in lattice simulations. A series of recent studies has shown that the perturbative Curci-Ferrari model gives an accurate description of the phase structure of pure Yang-Mills the-ories and of QCD with heavy quarks [35][36][37]. The case of light quarks is more delicate because, unlike the couplings in the pure gauge sector, the quark-gluon coupling becomes significant in the infrared [38]. A systematic approximation scheme, nonperturbative in the quark-gluon vertex, has been proposed in Ref. [39], based on a double expansion in powers of the pure gauge coupling and of the inverse number of colors 1/N c . At leading order, this leads to the well-known rainbow equation and, therefore, successfully describes the dynamics of chiral symmetry breaking in the vacuum. There are however too notable differences with respect to the usual implementations of the rainbow equation: first, the gluon exchange is described in terms of a tree-level Curci-Ferrari propagator, and second, higher order corrections are controlled by small parameters, which allows in particular to include renormalization group effects in a consistent way [39].
The purpose of the present work is to extend this approach to nonzero temperature and nonzero chemical potential, paving the way for a systematic study of the predictions of the CF model for the QCD phase diagram.

II. THE RAINBOW EQUATION IN THE CURCI-FERRARI MODEL
For simplicity, we study a theory with N f degenerate quark flavors. At nonzero temperature T and quark chemical potential µ = µ B /3, the rainbow equation for the (Euclidean) quark propagator S reads where M 0 and g 0 denote the bare quark mass and quarkgluon coupling, respectively. We have introduced Euclidean momenta P = (ω p , p), Q = (ω q , q) and K ≡ P − Q = (ω k , k), with ω n = 2πnT andω n = 2π(n + 1/2)T (n ∈ Z) the bosonic and fermionic Matsubara frequencies respectively. Correspondingly, the bosonic and fermionic Matsubara sums will be denoted where it is implicitly understood that integrals over the norm of three-dimensional momenta are cut off at a scale Λ. The matrices γ µ stand for the Euclidean Dirac matrices, with {γ µ , γ ν } = 2δ µν , which we choose in the Weyl basis, such that γ * 0,2 = γ t 0,2 = γ 0,2 and γ * 1,3 = γ t 1,3 = −γ 1,3 . Finally, the tree-level gluon propagator is with P ⊥ µν (K) = δ µν − K µ K µ /K 2 the transverse projector and m the Curci-Ferrari mass.
In Appendix A, see also [13], we recall that the quark propagator decomposes as (we keep the µ-dependence explicit) with any of the componentsX =Ã 0 ,Ã v ,B orC depending on p only through its norm p ≡ | p| and such that X(−ω p , p ; µ * ) * =X(ω p , p ; µ) .
These considerations apply also to the inverse propagator S −1 (ω p , p ; µ) which we parametrize as We have X = ∆X, with Projecting Eq. (1) onto the various tensor components, one arrives at a non-linear system of integral equations: whereÂ 0 (P ) ≡ (ω p + iµ)A 0 (P ),Â v (P ) ≡ pA v (P ) and C(P ) ≡ p C(P ), as well asp ≡ p/p. We have also introduced the quadratic Casimir in the fundamental representation C F = 4/3. In the chiral limit, corresponding to M 0 → 0, an unbroken chiral symmetry implies B = C = 0, which obviously solve the (homogeneous) equations (10) and (13). This also means that a solution with either B = 0 or C = 0 signals the spontanous breaking of chiral symmetry. In what follows, we use B as our order parameter for chiral symmetry breaking since C = 0 remains an allowed solution (which we stick to) also away from the chiral limit. In order to keep the discussion as simple as possible, we also set the other functions to their tree-level values, With this ansatz, the rainbow equation for the quark mass function B, Eq. (10), reads where we have defined Q iµ ≡ (ω q + iµ, q) and we recall that, in the case of a real chemical potential, as follows from Eqs. (6) and (7). In particular, B is real for µ = 0 but becomes a priori complex when µ is nonzero. 1

III. LOCALIZATION
In this article, we analyse the solutions to Eq. (14) by using an approximation scheme called localization [40,41], which we now recall and extend to the problem at hand. In a given model, the value of the mass function at a given momentum depends on the value of the mass function at any other momentum. However, there can be cases where, in some range of parameters and to a reasonable level of accuracy, the value of the mass function at a particular scale decouples from the rest and obeys, therefore, a simpler, "localized" equation. For instance, in Ref. [40], the behavior of the mass function at zero momentum was essentially controlled by the zeromomentum mass itself. The self-consistent equation for this zero-momentum mode could be obtained by expanding the mass function about this zero-momentum value in the corresponding integrals.
Even in cases where there is not a clear argument of why a certain scale could decouple, the localized equa-tions often provide a good qualitative guide. In the present case, it correctly captures the phenomenology of chiral symmetry breaking in the vacuum.
Of course, rainbow equations similar to Eq. (14) have been solved before without the need to resort to localization, even at finite temperature. Localization is however a convenient approach that could be used as a first investigation of more intricate settings such as the rainbow equation in the presence of non-trivial gluonic background (accounting for the interplay between chiral and center symmetry), or the corrections beyond the rainbow approximation, within the systematic expansion scheme alluded to in the Introduction.
One of the goals of this work is to investigate this strategy in a simpler setting before applying it to these more complicated cases. In order to test the robustness of the approach, we shall consider two types of localizations, referred to below as Euclidean and physical respectively, and which we now define more precisely.

A. Euclidean localization
At finite temperature, it is important to stress that the Euclidean mass function B is defined a priori on the fermionic Matsubara frequencies, which never vanish. We shall therefore localize the mass function at the smallest momentum available, that is, B(ω 1 ) ≡ B(ω 1 , 0; µ). Since B(ω 1 ) = B(−ω 1 ) * , it is natural -and in fact crucial, as we illustrate below -to localize B(ω 1 ) together with B(−ω 1 ). This means that we should consider two regions in the integrals, the one where it makes sense to expand the mass function aboutω 1 and the one where it makes sense to expand it about −ω 1 . Therefore, we approximate (14), with P = (ω 1 , 0) or P = (−ω 1 , 0) by and which are easily checked to be compatible with B(ω 1 ) = B(−ω 1 ) * . It is convenient to work with the real quanti- In terms of B r and B i , the rainbow equations for B(±ω 1 ) read and After performing the Matsubara sums and the angular integrals, Eqs. (20)-(21) rewrite in the simple form where, for notational convenience, we have defined with with ε x y ≡ x 2 + y 2 and where n (±) x = (e x ± 1) −1 denote the Bose-Einstein and Fermi-Dirac distributions. We note that F (B * ) is not the complex conjugate of F (B), due to the dependence on iω 1 which we leave, however, implicit in what follows.

B. Physical localization
One inconvenient aspect of the Euclidean localization is that one has to deal with two variables, B r and B i . This prevents the definition of a potential associated to the localized equations, due to the fact that the latter do not comply with the Cauchy-Riemann conditions. Although many features of the phase diagram do not require the existence of an underlying potential, it is convenient to find a setting where one deals with only one variable instead of two. 2 Another drawback of the Euclidean localization is that it involves the first Matsubara frequency πT . Therefore, we expect its quality to decrease as the temperature is increased.
One way to try to cope with these limitations is to consider a localization based on the retarded mass function evaluated for p 0 = 0 and p = 0. The reason for the presence of µ in the prescription to obtain the physical retarded Green's function is recalled in Appendix B, see also Ref. [42]. We find that, for B real and smaller than m, the corresponding analytic continuation of F (B) in Eq. (24) is real in the limit p 0 → 0. 3 Therefore, the equation for B i becomes compatible with the solution B i = 0, which we assume from now on, and the equation for B r ≡ B reduces down to where F (B) is to be evaluated withω 1 → 0 and µ → 0 in the energy denominators. While F vac (B) remains the same as above, we now have 2 We will see below that, after renormalization, there is one particular way to achieve this within the Euclidean localization. There are also certain limits, such as T → 0 or µ → 0, where this becomes possible due to the fact that B i → 0. 3 One should pay attention to the fact that what is really continued are not Br and B i , but rather B(ω).
One price to pay for this choice of localization is the singularity at B = m, which is regulated in the retarded self-energy through B 2 − m 2 → B 2 − m 2 + i0 + . We ignore this issue and restrict our analysis to cases where B < m. This is justified both because the physically relevant values in the vacuum fall in this range (see below) and because most of our discussions below concern the vicinity of the symmetric point B = 0. 4 We note finally that, interestingly (although not surprisingly), the localized equation (28) reduces to the corresponding meanfield gap equation of a NJL-type model with an effective nonlocal four-fermion vertex corresponding to a massive gluon exchange.

IV. RESULTS
We now investigate our predictions for the phase diagram using the two types of localizations, starting with the physical localization.
where we see that symmetry-breaking solutions exist only if g 2 0 > π 2 /F vac (0). We use Eq. (30) to trade the bare quark-gluon coupling g 0 for the (ultraviolet finite) quark mass B 0 . The rainbow equation rewrites as where we note that the cut-off can now be sent to infinity: It is useful to interpret Eq. (31) as deriving from a chirally symmetric potential . The absolute minima of W (B 2 ) then determine the state of the system. At µ = 0, one easily checks that the nontrivial minimum, equal to B 0 in the vacuum, decreases with increasing The first-order line is then determined from where B min is the nontrivial minimum at the transition. The associated lower and upper spinodals are respectively defined by with B sp the location of the nontrivial metastable state at the upper spinodal. The two spinodals flank the first order line and merge at the tricritical point, beyond which the lower spinodal becomes the critical line. The equation governing the critical and lower spinodal lines is easily obtained as the approximate expression This is similar to the result obtained in the quark-meson model with a large-N f approximation [25].
We note that there is an ambiguity in the definition of the potential W (B 2 ) since neither the solutions of Eq. (31) nor their convexity are altered by the replacement W (B 2 ) → f + (B 2 )R(B 2 ) with f + (B 2 ) a differentiable and strictly positive function. Interestingly, because the conditions (33) and (35) only involve W and its derivatives, they are, in fact, independent of the function f + and so are, thus, the spinodal lines, the line of second-order transition, and, of course, the tricritical point where these lines meet. Only the line of first-order transition explicitly depends on f + through the second condition in Eq. (34). However, it always lies in between the two spinodals.
In Fig. 1, we show our results for the phase diagram in the chiral limit. We use the typical values m = 500 MeV and B 0 = 300 MeV, motivated by the study of dynamical chiral symmetry breaking in the Curci-Ferrari model in the vacuum [39]. We find a tricritical point located at (µ, T ) ≈ (237 MeV, 69 MeV). The transition at zero chemical potential occurs at T c (µ = 0) ≈ 141 MeV and the two spinodals meet the T = 0 axis for µ ≈ 268 MeV, and µ ≈ 305 MeV, respectively. This gives an estimate of the first-order transition line at most at the 10% level, independently of the function f + (the estimate improves as one approaches the tricritical point). For the choice f + = 1, the T = 0 transition point is at µ ≈ 287 MeV.
if µ ≥ B and F th (B) = 0 otherwise. In particular, we see that the value of the order parameter below the transition point is independent of µ, B(T = 0, µ) = B 0 , until it jumps to B = 0 in the symmetric phase. This is known as the Silver Blaze property [43,44], which we illustrate in Fig. 2. 5 Also, using Eqs. (32) and (38) Let us now move away from the chiral limit. Equation (31) now reads 2BR(B 2 ) = H, where H ≡ π 2 M 0 /g 2 0 needs to be seen as a finite parameter controlling the departure from the chiral limit. For H = 0, the second- 5 We mention that the Silver Blaze property should in princple extend only up to the the first singularity on the T = 0 axis, namely the nuclear liquid-gas transition. Here, however, our level of description does not capture the corresponding dynamics, and the Silver Blaze property extends further. We also mention that the µ-independence below the first singularity applies in principle only to 0-point functions. For higher n-point functions, it takes a more general form as shown in [44]. However, due to the presence of µ in the retarded prescription (27), it can be argued that the Silver-Blaze property applies to retarded Green functions as it does for 0-point functions.
from which one extracts B c , T c , and µ c for each H. To this aim, it is convenient to vary B, determine T c (B) and µ c (B) from the last two conditions in Eq. (39)that do not involve H-and then deduce H(B) from the first condition. Inverting this relation, one then gains access to T c (H) and µ c (H). In particular, we find that the approach to tricriticality is governed by mean field exponents; see Fig. 3. This is expected because the potential is regular around B = 0. The trajectory of the critical endpoint in the phase diagram, shown in Fig. 1, exhibits a nonmonotonous behavior of T c as a function of µ, similar to that observed in Ref. [12] using an approach based on the Cornwall-Jackiw-Tomboulis effective potential. Finally, Fig. 4 shows the interval [x − (H), x + (H)] compatible with a critical point for each value of H. Interestingly, in the physical localization considered here, the Curci-Ferrari mass should be neither too large nor too small for a critical end point to exist.

B. Euclidean localization
Let us now investigate the Euclidean localization, that allows us to test the robustness of the previous features.
To renormalize the equations in this case, we proceed as follows. First, we note that the function F vac (B) possesses a logarithmic UV divergence which however does not depend on B. It is then readily checked that the equation for B r has a divergence proportional to B r . This divergence can be absorbed into a redefinition of the bare coupling as follows. We divide the corresponding equation by g 2 0 and set 1 where the renormalized coupling g should be interpreted as being defined at the (real) renormalization scale B . Introducing H ≡ π 2 M 0 /g 2 0 as before, the renormalized equation for B r takes the form As far as the equation for B i is concerned, it is easily checked that it is finite, for any fixed g 0 . Therefore, using the redefinition (40) or the bare coupling is problematic here since it leads to a spurious cut-off dependence. This problem is rooted in the localization procedure which does not commute with the renormalization procedure, see [41] for more details. However, we can always define a renormalized localized scheme by replacing g 0 by g in the equation for B i : We note that, in this equation, one can interchangeably use F orF .
With this choice of renormalization, we can eventually express the equations in terms of the vacuum mass in the chiral limit, B 0 , see Eq. (30). It is related to g and B by Replacing g 2 in that form, one checks that the first term in the RHS of (41) disappears while the scale B inF (B) is replaced by B 0 . Thus, Eq. (41) does not depend on B . However a scale dependence remains in equation (42), as expected at a given order of approximation. Below, we shall test the dependence of our results on the renormalization scale B .
Let us also mention that, because of the localization procedure, the coupled gap equations (41) and (42), which we denote formally as H = R r (B r , B i ) and 0 = R i (B r , B i ) in what follows, cannot be seen as deriving from a potential, because, in general the Cauchy-Riemann condition ∂R r /∂B i = ∂R i /∂B r is not satisfied. However, certain features of the phase diagram can be defined without the need of a potential because they correspond to the merging of different solutions of the gap equations. For instance, suppose that we want to investigate whether B r becomes critical. To this purpose, we solve for B i as a function of B r from its gap equation: and construct a potential for B r by integrating In the chiral limit, the criticality condition reads 0 = V (0). Simple algebra using (44) and (45) leads to the condition and therefore ∂R i /∂B i | Bi=0 = π 2 /g 2 . From these remarks, it follows that the condition for a critical point in the chiral limit simplifies to The same equation defines the lower spinodal in the case of a first order phase transition. The upper spinodal is also determined from (46) but without evaluating it for B r = B i = 0 and coupling it to the gap equation for B r . Finally, the tricritical point is determined from where the indices u, v, and w take the values r or i, and ∆B = (∂R i /∂B i , −∂R i /∂B r ) Br=Bi=0 . We have again made use of ∂R r /∂B i | Bi=0 = 0. This formula simplifies further because ∂ 3 R r /∂B 2 r ∂B i | Bi=0 = ∂ 3 R r /∂B 3 i | Bi=0 = 0. We note however that we are not able to fully eliminate B i , contrary to what happened for the critical point [see below for a particular limit where this becomes possible].
Our results in the chiral limit (with B = 1 GeV) are summarized in Table I and compared to the results in the physical localization as well as to the results in other approaches. The last column shows the values of µ at which the lower and upper spinodals (µ ls and µ us ) and the first order transition (µ first ) are reached for T = 0. As already discussed in the previous section, the value of µ first has the largest uncertainty since it depends on the potential which is not uniquely defined in our approach. Also, it may look surprising that we could obtain a value µ first in the Euclidean localization case since there is no potential in this case compatible with the gap equations. However, in the T → 0 limit, it is readily checked using (15) that B i vanishes. We note also that, along the T = 0 axis, B ≡ B r is not constant below the transition. This is of course not in contradiction with the Silver-Blaze property since only 0-point functions should be constant. The case of the physical localization is a bit peculiar since the retarded prescription (27) with µ included makes the retarded function behave like a 0-point function as far as the Silver-Blaze property is concerned.
We mentioned above that it was crucial to localize Chiral limit (H = 0) µ tric T tric T c µ ls µ first µ us simultanously in B(ω 1 ) and B(−ω 1 ). Let us illustrate this point here. In Fig. 5, we show the curves R r (B r , B i ) = 0 and R i (B r , B i ) = 0 for decreasing temperatures and a large enough chemical potential. The crossings correspond to the various possible solutions in the chiral limit and, because the chemical potential is chosen large enough, we should observe a first order transition pattern. Let us now see how this comes about. The first plot is at a temperature right above the upper spinodal, that is the appearance of a new crossing between the curves R r (B r , B i ) = 0 and R i (B r , B i ) = 0, at which two new extrema are about to appear. 6 Below this temperature, the various branches making the curve R r (B r , B i ) = 0 fuse and reorganize, in such a way that, at an even lower temperature, a second spinodal occurs at B = 0 where two extrema merge. We observe that the proper realization of this scenario requires not only the various branches to fuse at some temperatures, but also that the number of intersections of the curves R r (B r , B i ) = 0 and R i (B r , B i ) = 0 changes from one to five, and then to three, as the temperature is decreased. Had we performed the localization only with respect to B(ω 1 ), that is by writing only the first term in Eq. (16) and taking the real and imaginary part of the corresponding equation, this second requirement would not be fulfilled, as we have checked explicitly.
Regarding the renormalization scale dependence of our results, we observe numerically that the critical/lower spinodal line does not depend on the scale B . This is no surprise since the corresponding equation (46) depends only on R r (B r , 0) which, has we have already argued above is B independent. In contrast, the position of the tricritical point along (46), or even the other spinodal emerging from this point do depend on B . We note, however, that the inverse coupling 1/g 2 diverges positively as the renormalization scale is taken to infinity, see Eq. (43). Since there is no other dependence with respect to B in Eq. (42), it follows that B i should approach 0 in this limit and all relevant features (boundary lines, tri/critical points, . . . ) should converge to a certain limit, obtained by considering a single gap equation (41) in which one sets B i = 0 from the start. 7 Take for instance the tricritical point. Because ∆B r = −π 2 /g 2 ∆B i in the limit B → ∞, Eq. (48) becomes which is indeed the condition for a tricritical point if one restricts from the beginning to Eq. (41) with B i = 0. The relative difference between the tricritical values for B = 1 GeV and B → ∞ are found to be about a few percent. This indicates a controlled renormalization 7 This equation does not become trivial in the limit B → ∞ because, in this case, the B -dependence of 1/g 2 is cancelled by the corresponding B -dependence hidden inF . We also mention that the equation obtained in the B → ∞ limit is nothing but the one we would have obtained by applying the naïve renormalization and sending the cut-off to infinity. Indeed the remaining cut-off dependence in the equation for B i , only present in the term containing 1/g 2 0 , would enforce B i → 0 as Λ → ∞.
scale dependence and allows us from now on to work in a simplified picture in the B → ∞ limit. In particular, in this limit, we can associate a potential to the Euclidean localization, such that V (B) = R r (B, 0). We shall now employ this simpler setting to move away from the chiral limit. 8 For a non-zero bare mass, chiral symmetry is explicitly broken. As already mentioned, the second order transition line turns into a crossover and the tricritical point into a critical point. Unlike the critical endpoint, the crossover line has no unique definition. Moreover, there exist as many crossover lines as there are order parameters. In what follows, we define the crossover temperature by the inflection of B as a function of the temperature. With this choice, we can determine for which value of H the crossover temperature of the quark mass function becomes T χ = 170 MeV in the limit of vanishing chemical potential, which is the value found by lattice simulations at the physical point for two flavors [59]. Within our approach, we treat this particular value of H to correspond to the "physical point", H phys . We can determine it conveniently as follows. Denoting by H = R(B) the gap equation in both localization schemes and taking two T -derivatives, we have Imposing the inflection condition d 2 B/dT 2 = 0 and upon plugging (50) into (51), we arrive at the following condition  certainly fall within the group of lower temperatures and larger chemical potentials.
Finally, one can study how our findings for the phase diagram depend on the gluon mass of the Curci-Ferrari model. While m = 500 MeV is the value that globally works best in both the pure Yang-Mills as well as the unquenched sector, it is nonetheless insightful to vary it as a free parameter. Thereby, for each value of m, we always insist on fixing the coupling such that we keep the T = µ = 0 solution B 0 fixed at 300 MeV, in the chiral limit. Away from the chiral limit, we only vary H, without further changing g.
In Fig. 6, we display the position of the tricritical point in the chiral limit as the Curci-Ferrari mass parameter is varied. As can be seen, the obtained trajectories are qualitatively quite different depending on the considered localization scheme, although for m = 500 MeV, the tricritical points are not so far apart, in particular in temperature values. Interestingly, while in all localization schemes considered, the gluon mass can never exceed an upper limit m + , in the Euclidean localization, one might take m → 0 without losing the tricritical point, so m − = 0 in this case.
As before the definition of m ± is trivially extended to the case of non-zero H , where the tricritical point is replaced by a critical end-point. We then show our results for these values in dependence of the symmetry breaking parameter H in Fig.7 and in comparison with the findings in the physical localization. Here, we iterate our observation that the existence of a CEP puts an upper bound on the allowed for values for the gluon mass in both localization schemes considered, whereas a lower bound only exists for the physical one.

C. Chiral condensate
As mentioned previously, the mostly used order parameter for the chiral transition is not the constituent quark mass but the chiral condensate. Within the localized schemes considered here, it is natural to approximate the chiral condensate as with At this level of description, we can adopt an adhoc renormalization of the condensate by removing the divergence in J B , up to the scale in the logarithm of the vacuum contribution of J B which makes the renormalized condensateσ a scale dependent quantity: To get some intuition on the behavior ofσ, consider the physically localized gap equation which we rewrite identically as In the limit T → ∞, the tadpole sum-integralJ m grows like ∼ T 2 /12, which is only counter-balanced provided B vanishes as 1/T 2 , in which case we also haveJ B ∼ −T 2 /24. Putting all the pieces together, we find and thusσ In the Euclidean localization the condensate is found to grow with T 2 at large T , which is clearly an artefact. 9 Despite this feature, in both cases, we can define a crossover temperature associated to the inflexion point of the chiral condensate as a function of the temperature. It is found to be lower than the corresponding crossover temperature for B. We can also try to fine-tune the value of H to bring the crossover temperature for the condensate to the lattice value of 170 MeV. To this purpose, consider the equation d 2σ /dT 2 = 0. Seeingσ as a function of B and T , this equation rewrites which we use again to determine the crossover value for B, assuming a crossover temperature T χ = 170 MeV. The corresponding value of H is then obtained from the gap equation. In the physical localization, we find 9 The coefficient is proportional to H however in such a way that the condensate becomes zero in the chiral limit, in the high temperature region. H phys = 117 MeV, and a critical end-point located at (504 MeV, 11 MeV), but the value of B gets suspiciously close to the bound B = m. In the Euclidean localization, it seems not to be possible to reach these transition temperatures, at least at this level of approximation.

V. CONCLUSION
We have computed the phase diagram of QCD with light quarks in the context of a first-principle inspired approach to infrared QCD, the Curci-Ferrari model, using the double expansion in the pure gauge coupling and in 1/N c proposed in Ref. [39]. To allow for a semi-analytic grasp of the corresponding equations and since our first aim is a qualitative survey of what to expect when the equations are solved in full glory, we have used some simplifying approximations, in particular the localization scheme discussed in [40,41], which we have extended to the present context. For the parameters used here, the leading-order results agree well with those of effective quark-meson models when the chiral anomaly is neglected [25,29]. Although subleading in 1/N c , the latter and meson fluctuations are important to correctly deter-mine the phase structure in the Columbia plot [7,29,60]. In principle, they can be systematically included at nextto-leading order in the present expansion scheme.
We also mention that similar rainbow equations have already been considered in great detail to study the phase diagram in the context of nonperturbative functional approaches [12,22]. Not too surprisingly, our results agree qualitatively with those (the main differences may be attributed to the different treatments of the gluon sector). However, the key point of the present work is to systematically justify the employed approximation on the basis of identified small parameters in QCD.
Although we used here simplified versions of the complete rainbow equation (1), we expect, inspired by the vacuum case, that our main results are robust. In particular, we have tested that our results do not depend much on the type of localization we use. One notable exception is the fate of the (tri)critical point as the gluon mass is taken to zero. In some scenarios, the existence of a critical end-point seems to require a nonzero gluon mass. It would be interesting to investigate to what extent this is an artefact of the localization procedure by solving the original equation (14). Another, obvious extension of the present work will be to include the other scalar functions A 0 , A v , and C in Eq. (1), see Refs. [13,14].
Yet another interesting direction of investigation is to include the order parameter of the deconfinement transition, the Polyakov loop, in the spirit of Refs. [35][36][37], which would allow one to study the interplay between the chiral and deconfinement phase transition across the Columbia plot [26,61,62].