Functional renormalization group study of the critical region of the quark-meson model with vector interactions

The critical region of the two flavour quark-meson model with vector interactions is explored using the Functional Renormalization Group, a non-perturbative method that takes into account quantum and thermal fluctuations. Special attention is given to the low temperature and high density region of the phase diagram, which is very important to construct the equation of state of compact stars. As in previous studies, without repulsive vector interaction, an unphysical region of negative entropy density is found near the first order chiral phase transition. We explore the connection between this unphysical region and the chiral critical region, especially the first order line and spinodal lines, using also different values for vector interactions. We find that the unphysical negative entropy density region appears because the s=0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s=0$$\end{document} isentropic line, near the critical region, is displaced from its T=0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T=0$$\end{document} location. For certain values of vector interactions this region is pushed to lower temperatures and high chemical potentials in such way that the negative entropy density region on the phase diagram can even disappear. In the case of finite vector interactions, the location of the critical end point has a non-trivial behaviour in the T-μB\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T-\mu _B$$\end{document} plane, which differs from that in mean field calculations.


I. INTRODUCTION
The phase diagram of Quantum Chromodynamics (QCD) is a widely studied topic by both experimental and theoretical physics and much has been learned about its properties since its first conjecture by N. Cabibbo and G. Parisi [1].However, the phase structure at low temperatures and high baryonic density remains a mystery, e.g., the existence of a first order phase transition and the critical end point (CEP).
Heavy ion collision (HIC) experiments conducted by the STAR Collaboration in the Relativistic Heavy Ion Collider (RHIC) at Brookhaven National Laboratory [2][3][4] and by NA61/SHINE Collaboration in the Super Proton Synchrotron (SPS) at CERN [5,6], are currently, not only studying the properties of the quark-gluon plasma (QGP), but also trying to map the phase boundary of QCD.In the future, other facilities like the Nuclotron based Ion Collider fAcility (NICA) at Joint Institute for Nuclear Research [7], Facility for Antiproton and Ion Research (FAIR) at GSI Helmholtzzentrum für Schwerionenforschung [8] and J-PARC Heavy Ion Project at Japan Proton Accelerator Research Complex (J-PARC) [9], will also join the collective effort to better understand the properties of nuclear and quark matter under extreme conditions of temperature, density and in the presence of magnetic fields.
The low temperature and high density region of the phase diagram, where the CEP might exist, is not only interesting for nuclear and particle physics studies, but also extremely important for astrophysical applications, namely to study the evolution and properties of neutron stars.Since the equation of state for nuclear matter derived from first principles is still unknown, the core composition of these objects is still an open question and several options have been proposed.Some model calculations propose different neutron star core compositions such as hyperon matter, pion or kaon condensates and quark matter [10,11].
From the theoretical side, lattice QCD (LQCD), a first principles method, is not able to shed light on these questions since it is not yet possible to do LQCD calculations at finite chemical potential due to the famous sign problem which renders the importance sampling needed in Monte Carlo simulations ineffective [12].Recently, different approaches have been tried to circumvent this problem like reweighing, Taylor series expansions, imaginary chemical potential and Complex Langevin dynamics [12][13][14].Due to these shortcomings, to seek for qualitative and increasingly quantitative understanding of QCD matter, other theoretical tools have been applied to study the phase diagram, such as Dyson-Schwinger equations and effective model calculations.Some of these calculations predict a first-order phase transition and a CEP in the low temperature and high density region of the phase diagram [15][16][17][18][19][20][21].
Model calculations, using the Nambu−Jona-Lasinio (NJL) model or the quark-meson (QM) model, can be improved by going beyond the common mean field approximation (MF) [22][23][24][25][26][27][28][29][30].Usually, when dealing with these chiral effective models, the quark contribution to the path integral can be reduced to a quadratic interaction and be exactly integrated out.The remaining path integrals, usually related to meson degrees of freedom, are not quadratic and some approximation has to be performed in order to obtain an effective action.In the MF approximation the remaining path integrals are calculated by the saddle-point approximation which effectively means that the only field configuration taken into account is the classical one, all quantum fluctuations to the remaining fields are left aside.
One way to go beyond the MF approximation is using an application of the renormalization group to continuous field theories, the Functional Renormalization Group (FRG), a powerful non-perturbative method which allows to incorporate quantum and thermal fluctuations in a field theory.The renormalization group is an important tool in theoretical physics since it allows the study of physical phenomena in different scales of distance and/or energy with enormous range of applications such as: studying the strong interaction, the electroweak phase transition, effective models of nuclear physics, condensed matter physics systems and quantum gravity [31][32][33].Some of its most important applications in the history of physics are the elimination of ultraviolet divergences in renormalizable quantum field theories and its application to explain the universality properties of continuous phase transitions.The FRG has been extensively used to study the QCD phase diagram using chiral effective models beyond the MF, like the NJL model [34][35][36][37] and the QM model [38][39][40][41][42][43][44].For detailed reviews on the FRG method see [45][46][47].
The application of the FRG method to the 2-flavour QM model leads to the presence of an unphysical negative entropy density region in the low temperature and high density region of the phase diagram, near the critical region where a first order chiral phase transition and CEP are predicted by the model.This behaviour was first discussed in detail by R. Tripolt et al. in [48], although previous FRG studies have reported decreasing pressures with increasing temperatures [39,40].The authors have put forward some explanations for this unphysical region: the truncation used to derive the QM flow equation is not enough to define a thermodinamically consistent model beyond the mean field approximation or the specific choice of regulator, used to account for fluctuations in the model, is not appropriate.They also discuss the possibility that the source for such behaviour is physical like a pairing transition to a color superconducting phase or to the existence of inhomogeneous phases.For more details see [48].
In this work, we will consider the 2-flavour QM model with vector interactions to explore the connection between these vector degrees of freedom with the critical region predicted by the model and the unphysical negative entropy density region.Vector interactions are very important to describe in-medium properties and are widely used to describe neutron stars [49][50][51], study the curvature of the critical line [52] and vector meson masses.The general effect of these interactions on the phase diagram, in MF calculations, is to drive the first order phase transition and CEP towards lower temperatures and higher chemical potentials.For high enough vector couplings, the critical region disappears leaving a smooth crossover for the chiral transition for all values of temperature and chemical potential.The ω 0 and ρ 3 0 vector mesons will be considered.The ω 0 vector is known to stiffen the equation of state of quark matter while the ρ 3 0 can be very important in isospin asymmetric systems, acting as an isospin restoring interaction.Hence the inclusion of these degrees of freedom can be essential to describe certain physical systems at high densities like neutron stars.
This paper is organized as follows.In Section II the 2-flavour QM model, including vector interaction and the FRG formalism are presented.The vector degrees of freedom are frozen and the flow equations for the effective potential and entropy density are laid out.In Section III the results are presented and the effect of the vector interactions on the critical region and on the unphysical negative entropy density are discussed.Finally, in Section IV conclusions are drawn and further work is proposed.

II. MODEL AND FORMALISM
The 2-flavour quark-meson model is invariant under chiral symmetry i.e., SU (2) L × SU (2) R .It can be built by considering a quark field ψ, interacting with dynamical meson fields via symmetry conserving terms at the Lagrangian level.Considering the scalar and pseudoscalar fields, σ, π and the isoscalar-vector and isovectorvector fields, ω µ and ρ µ , the following symmetry conserving Lagrangian density, in Minkowski spacetime, can be written: Here, the quark field ψ is a N c -component vector in flavour space, where each component is a Dirac spinor and τ are the three Pauli matrices.To study the system at finite density, a diagonal quark chemical potential matrix, μ = diag (µ u , µ d ), was also included.The field strength tensors F µν and R µν are used to define the kinetic terms for the ω µ and ρ µ fields, respectively, and are given by: The potential U (σ, π, ω µ , ρ µ ), must be invariant under chiral symmetry except for an explicit chiral symmetry breaking term, that tilts the potential to give a finite mass to the Goldstone mode, the pion.At the mean field level, for the the σ and π fields, this potential can include arbitrary powers of the chiral invariant φ 2 = σ 2 + π 2 .
Regarding the vector field contributions to the potential, several types of terms can be included, as long as the symmetries are respected.Due to the nature of the FRG calculation, one has only to specify the potential at the ultraviolet scale.The effect of current quark masses is to explicitly break chiral symmetry at the Lagrangian level, giving rise to a (slightly) massive Goldstone mode.This can be accomplished in the QM model by adding to the potential a non-vanishing expectation value for the σ field, This field will behave as an order parameter for the chiral transition.Finite temperature can be included using the Matsubara formalism in which a Wick rotation to Euclidean spacetime is applied to the action.To simplify the notation, we introduce the fields, φ i = {σ, π} and V i µ = {ω µ , ρ µ }.Using the Euclidean action S E = 1 /T 0 dτ d 3 x L, the generating functional of the fully connected Green's functions, for a given temperature (T ) and chemical potential (µ), is defined as: where we have included sources for the scalar fields J i and for the vector fields j i µ , omitting the sources for the fermion fields which can be integrated out.First we are only interested in dealing with the path integral over the vector fields hence, we write: We have defined the effective action for vector degrees of freedom as: Writing explicitly only the functional dependence on j i µ , the effective action can be computed by Legendre transforming W j i µ as follows: where Ṽ i µ is the expectation value of the vector fields V i µ , in the presence of an external source j i µ and it is defined as: The effective action can be written as [53]: In a chiral effective model, the most important dynamics comes from chiral symmetry breaking.This means that the dynamics of the more massive fields, will play a secondary role.Hence, following previous approaches [54], the vector fields will be used to model unknown degrees of freedom at short distances.This allows the use of the saddle-point approximation to solve the path integral in Eq. ( 10): the classical trajectories will be the most important for such fields, effectively freezing these heavier modes.
Using the saddle-point approximation, the main contribution to the integral will come from the minimum of the action S V [V i µ ] .Taylor expanding the action S V [V i µ + Ṽ i µ ] around Ṽ i µ one can get the effective action in the mean field approximation, Where the mean field configuration is calculated from: Due to rotational invariance, the spatial components of the mean fields Ṽ i j , vanish [54].Since we are not interested in studying the condensation of mean fields that change the properties of the vacuum, the non-diagonal elements ρ1 = ρ2 = 0 will also be zero.Therefore, only the fields ω0 and ρ3 0 can have non-zero values.These fields can be absorbed in the definition of the effective quark chemical potential matrix, μ, as: As expected, the mean field ρ3 0 , introduces an isospin asymmetry [54].Using Eqs.(7) and (11), we can write the effective action as: The same approximation could be performed in the remaining meson path integrals and the quarks can be integrated out exactly, yielding the quark-meson model in the mean field approximation.However, in the present work, we will go beyond mean field by taking into account quantum fluctuations of the σ and π fields using the FRG method.Modifying Eq. ( 14) with a regulator term, the effective average action can be defined through a modified Legendre transformation [45].
A. The FRG method In the formalism of the FRG, the central object is the average effective action, Γ k .This object depends explicitly on a momentum scale k and has well defined limits: at the momentum scale k = Λ, we have the classical action to be quantized, S, at the momentum scale k = 0, all quantum fluctuations have been included and we obtain the full quantum effective action, Γ i.e., Γ k→Λ = S, Γ k→0 = Γ.
The average effective action interpolates these regimes in the space of all possible actions.The behaviour of this quantity during the renormalization group flow is governed by the so-called Wetterich equation [55].For boson fields this equation is given by: while, for fermions, it can be written as: Here, t = ln k Λ is the adimensional renormalization time with respect to some cutoff momentum Λ.The derivatives of the average effective action, Γ k and Γ (1,1) k , follow the usual notations for boson and fermion fields derivatives.The function R k is the so-called regulator function and it can be interpreted as a scale dependent mass term.As long as the interpolation between the ultraviolet and the infrared is correct, the regulator can take any functional form since it will only interfere in the arbitrary path taken, between these points in the theory space.Of course, since from the numerical point of view it is impossible to reach k = 0 [45], a finite infrared cut-off has to be applied meaning that different regulators might lead to different infrared effective actions.
These equations are exact functional differential equations for the effective average action which, in principle, can be solved given a set of initial conditions.Solving exactly the Wetterich is an impossible task due to the infinitely high coupled behaviour of the equation and some approximation scheme is needed.There are two widely used approximations schemes to solve this equation: the vertex expansion and the operator expansion.In the present work we will use the latter approach in the so-called local potential approximation (LPA), by building an effective average action based on the operator expansion with increasing mass dimension.
To use Wetterich's equation, a regulator function, which respects the interpolating limits of the effective average action, has to be chosen.We employ the so-called optimized or Litim regulator function [56], for bosons and fermions, respectively given by: After solving the flow equation one can relate the effective action in the minimum, with the grand canonical potential, Γ k=IR (T, µ) min = βV Ω(T, µ), to calculate several thermodynamic quantities of interest such as the pressure (P ), particle (ρ i ), entropy (s) and energy densities ( ), using the following relations [57]: The constant P 0 is the vacuum pressure i.e., P 0 = P (0, 0).

B. The flow equations
In the lowest order of LPA, only the potential is scale dependent and using Eq. ( 14), the imaginary-time average effective action can be written as: with the scale dependent grand potential, U k , written in terms the chiral invariant φ 2 = σ 2 + π 2 and of the mean vector fields: The contribution U χ k is a function of the chiral invariant only and the term U V k represents the contribution from the vector degrees of freedom.While the functional dependence of the chiral part of the potential is calculated during the flow, a mean field approximation is performed in the vector channels, and a functional dependence for U V k must be chosen.In this work we use: Applying the stationary condition of Eq. ( 12) to the effective average action in Eq. ( 23) is equivalent to requiring that the potential U k σ, π, ω 0 , ρ 3 0 is minimal with respect to the vector fields, at each momentum scale k [54,58], i.e., Hence, the vector fields acquire an implicit dependence on the RG scale k.This requirement ensures that the flow equation follows a path, in theory space, where the effective potential is always in the minimum with respect to the vector fields.Since no pion condensation will be considered, only the radial direction of the field, φ = {σ, 0}, will contribute and we can switch variables to σ.
Calculating the scale derivative of Eq. ( 23) and using the stationary conditions for the vector fields given in Eq.( 26), yields: Hence, ensuring at each momentum shell that the vector fields stationary conditions hold, one can simply solve the flow equation for U χ k with effective quark chemical potentials modified by the vector fields.Putting everything together leads to the dimensionful LPA flow equation for the effective potential U χ k : Here, the effective chemical potential, μk,i , is defined as: With, i = 0 for up quarks and i = 1 for down quarks.The vector contribution, v k,i , is defined as: The functions, n B (E) and n F (E) are the Bose-Einstein and Fermi-Dirac distribution functions respectively and, After solving the above flow equation, one has access to U χ k=IR .The full potential in the infrared, U k=IR , containing the contribution coming from vector fields, can easily be calculated with k = IR: The contribution coming from the vector fields can be calculated using Eq. ( 25) in the infrared.Using Eq. ( 26) and the flow equation ( 28), applying the substitution ∂ µ → −(ηE/p)∂ p and performing an integration by parts [59], the following self-consistent equations for the vector fields can be written: Where we have defined: Since the equations depend only on the product g ω ω0,k = ωk and g ρ ρ3 0,k = ρk , we take this combined quantity as variables.Likewise, the equations depend only on the combination gω mω = G ω and gρ mρ = G ρ , we take these ratios as parameters.
We are also interested in studying the entropy of the system including quantum fluctuations in order to understand what happens to the low temperature behaviour of this quantity.Hence, a flow equation for the entropy must be derived.Using Eqs. ( 21) and ( 28) and considering that the temperature derivative commutes with the scale derivative, the following dimensionful flow equation for the chiral contribution to the average entropy density, s χ k , can be derived: If considering finite vector interactions, there is an extra contribution coming from the temperature dependence of the vector fields, at each momentum shell, ∂ T v k,i .For the detailed calculation of this quantity, see Appendix VI B.
The entropy density in the infrared, s k=IR , containing the contribution coming from vector fields, can easily be calculated after solving the system of flow equations through: The contribution coming from the vector fields can be calculated using the stationary conditions for the vector fields given by Eqs.(35) and (36).
The system of coupled, partial differential equations, for the effective average action and average entropy, given in Eqs. ( 28) and ( 38) alongside the self-consistent equations for the vector fields, (35) and (36), must be solved numerically.One way to do so, is to use a Taylor expansion around the scale-dependent minimum of the effective potential U k .This method however, is not well suited to study the low temperature and high density regime of the phase diagram, where for certain parametrizations, a first order chiral phase transition is expected and two minima co-exist.In the present work we use the grid method, a much more powerful technique that provides full access to the effective potential, in a given range of the σ field.This allows the study of the phase diagram around a first order phase transition.In this numerical approach, the field variable σ is discretized in an one-dimensional grid, and the first and second derivatives of the effective potential with respect to σ are calculated using finite differences.For more information about the used numerical approach, see the Appendix (VI A).
The value of the vector fields, ω0,k and ρ3 0,k , are calculated using Eqs.( 35) and ( 36) at each momentum shell k.In practice, by following this approach, the flow of the effective potential and the entropy density are automatically in the minimum with respect to to the vector fields.
In the MF calculation, the self-consistent equation for the ω0 field is directly related to the sum of the quark densities while the one for the ρ3 0 field is related to the difference of the quark densities.This means that the ρ3 0,k field is zero for symmetric matter (ρ 3 0,k = 0), i.e. if considering µ u = µ d .The FRG calculation leads to a similar scenario.Indeed, in [54] it was shown that neglecting the boson quantum fluctuations, the MF results can be recovered.However, the choice of non-zero ultraviolet value of the ρ3 0 vector field, ρ3 0,Λ , would lead to an explicit isospin breaking interaction and to a non-zero ρ3 0 field, even for symmetric matter.Indeed, in [58], non-zero values for ω0,Λ were considered and their effect on the phase diagram was studied.However there is no reason to consider a ultraviolet potential with explicit isospin breaking by the ρ3 0 field.Hence, in the present work, we will consider ρ3 0,Λ = 0.In order to investigate the effect of the ρ3 0 field on the phase diagram and the unphysical negative entropy density region, an asymmetry between the quark flavours has to be considered.Following [60], we allow for different chemical potentials for each quark flavour, In principle, upon considering a finite δµ, pion condensation could happen.This means that the effective potential would dependent on two chiral invariants, φ 2 = σ 2 + π 2 3 and ξ 2 = π 2 1 + π 2 2 [61].In such a scenario, not only the flow equations would be much more complicated but a two dimensional grid would have to be considered since there are two distinct chiral invariants.Following previous works [51,54], to simplify the calculations, we will neglect the possibility of pion condensation and work only with one chiral invariant.To make this approximation valid, a very small difference between quark chemical potentials of δµ = −30 MeV will be considered [60].For such a value of δµ, we will be describing matter with more down quarks then up quarks, a very important scenario to study neutron stars, for example.

III. RESULTS
In this section we present the phase diagram of the 2-flavour quark-meson model, calculated by solving the flow equation (28), for different values of temperature and chemical potential.Different vector couplings are considered in order to study their effect on the phase diagram.We also present, for the same scenarios, the results of solving the flow equation for the entropy, given in Eq. (38).From this calculation we are able to study the behaviour of the entropy density near the critical region with and without vector interactions where an unphysical region, of negative entropy density, is expected from previous calculations [48].We also test the thermodynamical consistency by checking if the numerical temperature derivative of the effective potential agree with the result coming from solving the flow equation for the entropy density.
Regarding the numerical calculation, solving the system of coupled flow equations in a grid is very computationally demanding.In fact, the computational time is not only dictated by the grid size and the infrared cutoff, k IR but, in the case of finite vector couplings, of consistently solving Eqs.(35) and (36) for each grid point at every momentum shell (see the Appendix (VI A)).
In order to make the numerical computations more efficient within the scope of the present work, we decided to use a higher infrared cutoff of k IR = 80 MeV then the one used in [48] of k IR = 40 MeV.We verified, by solving the flow equations for different values of k IR , that this change does not influence the results qualitatively, allowing the study of the qualitative effect of different vector interactions in the phase diagram and in the unphysical negative density entropy region, using less computing resources.
Different grid sizes were also studied and after some analysis we decided to use a 80-point grid size in σ ∈ [2,122] MeV.As for the infrared cut-off, using a thinner grid does not change the qualitative behaviour of the results and since the same grid is used for every scenario, considering a given grid size represents a systematic uncertainty.
The system of flow equations can then be solved from the UV scale, k = Λ, down to the infrared scale, k =IR, to yield U k=IR and s k=IR , the effective potential and entropy density in the infrared.In order to solve this system of coupled partial differential equations, a set of initial conditions has to be provided.In the case of Eqs. ( 28) and (38), these correspond to the effective action and entropy density in the k = Λ momentum shell.The effective potential in the UV, U k=Λ , is chosen in such a way that it respects the symmetries of the system and to yield, in the infrared, the experimental values for pion mass and decay constant.In this work we use the usual potential: with the parameters given in Table I.The vector fields in the UV, ω0,Λ and ρ3 0,Λ , in this work, were chosen to be zero, In [58] also the effect of a σ dependence for the ω0,Λ vector field was studied, which ended up not changing the phase structure significantly.We will consider the vector coupling constants, G ω = g ω /m ω and G ρ = g ρ /m ρ as free parameters and study the influence of different values on the structure of the phase diagram.In [54], these parameters were considered as bounded by G ω = G ρ = 0.001 − 0.01 MeV −1 .These bounds were obtained using vacuum properties, by considering the vector fields as massive, m ω , m ρ ∼ 1 GeV and g ω = g ω = 1 − 10 [54].However these parameters might be density dependent and in-medium modifications could change their magnitudes.
Due to the fact that there is no temperature dependence in the UV potential, the UV entropy density, s k=Λ , is simply given by: Since the UV scale is fixed at a finite value, there is no reason why the effective potential in the UV, should be temperature and chemical potential independent [41].Indeed, in [62], only the purely thermal flow equation was solved, effectively generating a temperature and chemical potential UV potential.
In the presence of a first order chiral phase transition, the effective potential has two minima.The phase transition in this case will be defined through the Maxwell construction: when the effective potential has several minima, the one with lowest energy represents the stable phase.In Fig. 1, we present such a construction at T = 20 MeV for the QM model using the FRG method.The dot is the chiral transition chemical potential while the squares are the chemical potentials of spinodal points.
We start our study by considering the QM model without vector interactions.In Fig. 2 we present the first order phase transition of the model, the s k=IR = 0 line and the CEP.One can see that the region in-between spinodal lines is very narrow, differently from MF calculations.Indeed, for T = 20 MeV one can analyse Fig. 1 and observe that the overall size of the region in-between spinodal points is less then 1.5 MeV.We also present the line s k=IR = 0 which, trivially, separates the region of positive and negative entropy densities.This result is very similar to the one presented in [48].However, in [48], a region of negative entropy density was only found on the right side of the first order phase transition.Here, we find such an unphysical region on both sides of the phase transition line.More, we see that the s k=IR = 0 line behaves like an isentropic line that crosses the first order phase transition: in [63,64] when an isentropic line crosses the first order phase transition it enters the critical region, touches each spinodal line once and exits the critical region.The branches entering from outside the spinodal region until touching the phase-transition line correspond to the entropy in the stable minimum of the potential.The two parts of the line between the phase transition line and touching the spinodal lines correspond to the entropy in the respective local minimum.Finally, the branch between touching both spinodal lines follows the solution in the maximum of the potential.This leads us to argue that the s k=IR = 0 isentropic line is displaced from its T = 0 location in this model within the FRG approach.
A possible origin for the displacement of this isentropic line and consenquently the existence of the negative entropy density region is that finite chemical potential effects are not correctly accounted in the model beyond mean field.Upon considering an UV potential which is independent of the temperature and chemical potential, one is considering that the initial conditions to solve the differential equations are the same for every point in the phase diagram.Such case may not be true and considering temperature and chemical potential dependences in the UV potential are known to modify the thermodynamics and the phase structure [62].Hence, building a temperature and/or chemical potential UV effective potential could provide more insights on the origin of the negative entropy density region.Such a study is beyond the scope of the present work and is left as future work.The next step in our study is to consider the effect of finite vector interactions.First we just consider the effect of the ω0 field, by setting G ρ = 0, and increasing G ω .The critical region with increasing G ω can be seen in Fig. 3.For increasing vector coupling in the range G ω = [0.001,0.004] MeV −1 (see Fig. 2, panels (a), (b), (c), and (d)), there are two main effects regarding the critical region: the CEP is moved to much higher temperatures and smaller chemical potentials and the extension of the region in-between spinodal lines increases.The low temperature first order phase transition line is slightly shifted to higher chemical potentials while for higher temperatures the first order line is dragged along with the CEP to smaller chemical potentials.Further increasing the vector coupling, G ω = [0.006,0.015] MeV −1 (see Fig. 2, panels (e), (f), (g), (h), and (i)), a very different behaviour is observed: the CEP is moved towards smaller temperatures and higher chemical potentials while the region in-between spinodal lines gets imperceptibly smaller.The behaviour of the CEP for these values of G ω is very similar to the one found in [58] even tough in that study, the chiral limit is used.
The behaviour of the negative entropy density region and s k=IR = 0 line is very interesting: increasing the vector coupling G ω pushes this region to lower values of temperature.In fact, there is a critical value of G ω to which there is no more negative entropy density region on the phase diagram of the model.
From MF studies one expects that the inclusion of repulsive vector interactions would push the CEP towards lower values of temperature, making it disappear for a high enough vector coupling.However we observe a rather different and complex behaviour when including quantum fluctuations with the FRG.Indeed the CEP and first order phase transition do not disappear for the range of considered vector couplings and the previous unphysical negative entropy density region disappears for increasing G ω .
As already stated, in order to study the effect of the ρ3 0 vector field on the first order phase transition, the two flavour quark system must be on an asymmetric state.As already discussed, we will consider a finite isospin chemical potential of δµ = −30 MeV.
In Fig. 4, we show the results of comparing the critical region of the model with δµ = 0 and δµ = −30 MeV, without vector interactions i.e., G ω = G ρ = 0.The effect of considering a finite isospin is the following: the first order line is shifted to higher chemical potentials (at lower temperatures) and the CEP is marginally moved to lower quark chemical potentials but its temperature remains the same (within our level of numerical accuracy).Since the s k=IR = 0 isentropic line is connected to the spinodal region, moving the first order line to higher chemical potentials also moves the unphysical negative entropy density region.The inclusion of a finite δµ also enlarges the region in-between the spinodal lines.MeV.The red, black and blue lines are the spinodals, first order phase transitions and the s k=IR = 0 lines, respectively.The CEPs are represented by the black dots.Entropy density is negative below the s k=IR = 0 line.
In order to study the isolated effect of the ρ3 0 vector field with δµ = −30 MeV, we set G ω = 0 and calculated the phase diagram for increasing values of G ρ .The results can be seen in Fig. 5. Increasing the coupling G ρ , has the opposite behaviour of considering a finite δµ: it shifts the first order line to smaller chemical potentials while the CEP is slightly moved to higher chemical potentials and low temperatures.The region in between spinodal lines is also larger with finite G ρ when compared to the case without vector interactions, even tough the effect is much less noticeable than when considering finite G ω .The first order phase transition line at low temperatures is very close to its original location with δµ = 0 for G ρ = 0.008 MeV −1 .Thus, increasing this coupling is effectively restoring the isospin symmetry, broken by the finite δµ.Indeed, in nuclear relativistic mean field models, the ρ3 0 vector field can be added to the theory as an isospin restoring interaction, mirroring the Bethe-Weizsäcker mass formula and the valley of beta stability in nuclear physics [65].
Finally in Fig. 6 we consider G ω = G ρ = 0.008 MeV −1 , with δµ = −30 MeV.In this scenario we are taking into account the combined effect of the ω0 and ρ3 0 vector fields.The obtained phase diagram is extremely similar to the one obtained in the Fig. 3 panel (f), with G ω = 0.008 MeV −1 and δµ = 0.The only difference is on the location of the first order line which is negligibly dislocated to smaller chemical potentials.Taking the previous results into account, this behaviour is expected: the ρ3 0 field is restoring the isospin symmetry while the influence of the ω0 field is identical to the one observed in the isospin symmetric case.

IV. CONCLUSIONS
We have calculated the critical region near the first order phase transition of the two flavour QM model with vector interactions, within the FRG approach to include quantum fluctuations.Besides the first order chiral transition and the CEP, the spinodal lines were presented.The unphysical region of negative entropy density reported by [48] was also found and its behaviour due to the presence of vector interactions was studied.
The behaviour of the critical region under finite vector interactions is different from mean field calculations: increasing the repulsive vector interaction pushes the CEP towards higher values of temperature and lower values of chemical potential.Further increasing the vector interaction, drives the CEP to smaller temperatures and higher chemical potentials.Another important conclusion is that the region in-between spinodal lines increases in chemical potential with increasing vector couplings.Matter inside the spinodal region corresponds to unstable matter which can only be reached in a non-equilibrium evolution of the system in the form of clusterized matter.Very different from the case without vector interactions, the increase of this region in chemical potential, due to finite vector interactions, indicates that it is possible to have clusterized chiral symmetric matter in a wider region of densities.
We also found that the region of negative entropy density is present on both sides of the first order phase transition line.The positive entropy density region and the negative entropy density region is, trivially separated by the s kIR = 0 line.However, this line behaves like an isentropic line: it passes through the first order line, touches one spinodal line, changes direction crossing the first order line again, touches the other spinodal and changes direction again.This leads us to conclude that the appearance of the negative entropy density region is a consequence of the displacement of the the s = 0 isentropic line from its T = 0 location.For a high enough    vector interaction the negative entropy density regions disappears leaving a physical phase diagram with a first order phase transition and CEP and without negative entropy.
Considering a difference of up and down quark chemical potentials δµ, so a finite isospin chemical potential, has a big effect on the chemical potential of the first order line but the location of the CEP is unchanged in temperature and marginally changed to smaller chemical potentials.Increasing at a finite δµ the coupling of the ρ3 0 vector field, G ρ , is equivalent to restore isospin symmetry while pushing the CEP to lower values of temperature, leading to a phase structure similar to the one with δµ = 0 with a CEP at smaller temperatures.
To better understand the origin of the unphysical negative entropy density region, as previously found by [48], a flow equation beyond the LPA could be derived and the phase diagram and entropy density calculated.A different regulator function could also influence the results.Due to the mathematical nature of the QM flow equation, we were only able to calculate the phase diagram down to low temperatures (T = 5 MeV) but not at zero temperature.Solving the T = 0 flow equation exactly could also provide some new analytical and numerical insights.Some efforts in this direction have been done in [66], where the authors try to solve the flow equation at T = 0 by executing a mathematical transformation to the differential equations in order to transform the rectangular initial condition on a circular one, due to the Fermi sphere.Another possible source for the appearance of the negative entropy density region is the fact that the UV potential is temperature and chemical potential independent.As future work we plan to explore how different, temperature and chemical potential dependent UV potentials affect the phase diagram and the negative density entropy region.
In [67] it was demonstrated how to ensure numerical stability during the integration of a generalized version of Eq. ( 28), through an optimal step-size.To derive such optimal step, for simplicity, it was considered that the function derivatives are calculated with low order finite-difference methods: forward difference for the RG time variable, t, and three-point rule for the σ−direction.It is supposed that the numerical stability condition derived within this simpler scheme, is also valid for the fourth-order Runge-Kutta method, used in the t variable and higher order finite-differences used for the σ derivatives.Following this approach the following conditions for the step-size ∆t was derived: Here, G and F are given by: As in [67], we do not consider these conditions for σ ∼ 0. Since these conditions only depend on the bosonic sector of the flow equation, the fact that one is dealing with effective finite chemical potentials does not change the conditions directly.The effect of finite chemical potential and finite vector mesons only change these conditions indirectly, since the potential and its derivatives will be different during the flow.
To solve the set of coupled differential equations, in such a way to get full access to the full effective potential, we employed the grid method.In this method, the field variable σ is discretized in an one-dimensional grid, and the first and second derivatives of the effective potential with respect to σ are calculated using finite differences.The five point midpoint rule was used except in the grid endpoints where the forward and backward rules were used.
One starts the calculation in the UV scale i.e., at k = Λ.At this momentum scale the effective potential and entropy density are calculated using the initial conditions provided in Eqs. ( 42) and (45), respectively.The needed derivatives with respect to σ are calculated for every σ-grid point, using finite differences.Next, an optimal step size in the renormalization group time, ∆t, is calculated using Eqs.( 46) and ( 47): the smaller ∆t is used.The flow equations are then solved using the fourth-order Runge-Kutta method, to provide the σ-dependent effective potential and entropy density in the next step t − ∆t, i.e., U k=Λ exp (t−∆t) (σ) and s k=Λ exp (t−∆t) (σ).This process is repeated until the infrared scale is reached at k = IR.After reaching the infrared scale, one is in possession of the U k=IR (σ) and s k=IR (σ) and can then calculate the minimum of the effective potential, in which all observables are defined.
When considering finite vector interactions, the self-consistent Eqs.(35) and (36) have to be solved at every σ-grid point, at every momentum scale k.Hence, for a given k, for every σ-grid point, a 2-dimensional root finding algorithm is used to find the values of g ω ω0,k and g ρ ρ3 0,k that fulfil this system of equations.To speed up the root finding process, the solutions at the momentum scale k are provided as guesses for the next momentum shell.Since we are using the forth-order Runge-Kutta method this process has to be performed four times to be able to calculate the effective potential and entropy density in a given momentum scale.
The computing time is then related to the σ grid size, the infrared cutoff, k IR , and the root fiding precision, when considering vector interactions.The complexity of the flow equations also dictates the computing time since the step-size ∆t dictates how fast one goes from the UV down to the IR and different values of temperature and chemical potential influence the overall magnitude of the adaptive step-size.
One very important observation is that the optimal step-size calculated using Eqs.( 48) and ( 49) does not depend on the chemical potential and can still be used in the calculation with finite vector interactions.
In order to arrive at the phase diagram, the flow equation was solved multiple times for different values of temperature and chemical potential.In order to speed-up calculations, the OpenMP interface was used to run the computer code in parallel.
Very easily one can solve the system of linear equations to get: When considering only one vector field i.e., if G ρ = 0 or G ω = 0, the temperature derivatives are much simpler:

Figure 1 :
Figure 1: Extrema of the infrared effective potential as a function of the quark chemical potential for T = 20MeV.The phase transition point is represented by the dot while the squares are the spinodal points.

Figure 2 :
Figure 2: First order phase transition of the QM model without vector interactions.The red lines are the spinodals, the black line is the first order phase transition line and the CEP is the black dot.The blue line corresponds to the s k=IR = 0 line and below this line entropy density is negative.

Figure 3 :
Figure 3: First order phase transition of the QM model, for increasing values of G ω and fixed G ρ = 0.The red, black and blue lines are the spinodals, first order phase transition and the s k=IR = 0 lines, respectively, for each value of G ω .The CEPs are represented by the black dots.Entropy density is negative below the s k=IR = 0 line.

Figure 4 :
Figure 4: First order phase transition of the QM model without vector interactions, for δµ = 0 and δµ = −30MeV.The red, black and blue lines are the spinodals, first order phase transitions and the s k=IR = 0 lines, respectively.The CEPs are represented by the black dots.Entropy density is negative below the s k=IR = 0 line.

Figure 5 :
Figure 5: First order phase transition of the QM model with δµ = −30 MeV, for increasing G ρ with fixed G ω = 0.The red, black and blue lines are the spinodals, first order phase transitions and the s k=IR = 0 lines, respectively.The CEPs are represented by the black dots.Entropy density is negative below the s k=IR = 0 line.
n o d a l l i n e s 1 s t -o r d e r l i n e C E P

Figure 6 :
Figure 6: First order phase transition of the QM model with δµ = −30 MeV, for G ω = G ρ = 0.008 MeV −1 .The red lines are the spinodals, the black line is the first order phase transition line and the black dot is the CEP.