QCD equation of state matched to lattice data and exhibiting a critical point singularity

We construct a family of equations of state for QCD in the temperature range 30 MeV $\leq T\leq$ 800 MeV and in the chemical potential range $0\leq \mu_B \leq$ 450 MeV. These equations of state match available lattice QCD results up to $\mathcal{O}(\mu_B^4)$ and in each of them we place a critical point in the 3D Ising model universality class. The position of this critical point can be chosen in the range of chemical potentials covered by the second Beam Energy Scan at RHIC. We discuss possible choices for the free parameters, which arise from mapping the Ising model onto QCD. Our results for the pressure, entropy density, baryon density, energy density and speed of sound can be used as inputs in the hydrodynamical simulations of the fireball created in heavy ion collisions. We also show our result for the second cumulant of the baryon number in thermal equilibrium, displaying its divergence at the critical point. In the future, comparisons between RHIC data and the output of the hydrodynamic simulations, including calculations of fluctuation observables, built upon the model equations of state that we have constructed may be used to locate the critical point in the QCD phase diagram, if there is one to be found.


INTRODUCTION
The search for a possible QCD critical point is receiving increasing attention, which will culminate in the second Beam Energy Scan (BES-II) at the Relativistic Heavy Ion Collider (RHIC) at Brookhaven National Laboratory. The main goal of the BES-II program is to discover a critical point, or constrain its location, on the phase diagram of strongly interacting matter. One of the central questions that these experiments aim to answer is whether the continuous crossover [1] between quark-gluon plasma and hadronic matter that occurs as a function of decreasing T at µ B = 0 turns into a first order phase transition above some critical point at a nonzero µ B , corresponding to heavy ion collisions below some collision energy [2,3].
Lattice QCD simulations cannot currently be performed at finite density. For this reason, a first principle prediction of the existence and position of the critical point is still missing.
Several QCD-based models predict its location on the phase diagram, which naturally depends on the model parameters and approximations (for a review see e.g. [4]). This aspect makes the critical point search challenging, and is at the basis of the systematic scan in collision energies performed at RHIC. We anticipate that non-monotonous dependence of specific observables on collision energy will indicate the presence of the critical point as the freezeout point traverses the critical region [5,6]. As the BES-II approaches, it is therefore important to predict the effects of the critical point on several observables.
One of the main theoretical approaches to pursue this goal is represented by hydrodynamical simulations of the evolution of the fireball produced in heavy ion collisions (see e.g. [7] and references therein). While modifications of the hydrodynamical approach itself are needed in the vicinity of the critical point [8][9][10][11], the Equation of State (EoS) used as an input in these simulations needs to reflect all theoretical knowledge currently available as well as contain the singularity associated with the QCD critical point at a parametrically controllable location. Thus, the purpose of this manuscript is to produce a family of model equations of state for QCD, each of which contains a critical point somewhere in the region of the phase diagram covered by the BES-II at RHIC, and all of which respect what we know from lattice QCD calculations up to O(µ 4 B ). Previous hydrodynamical calculations at finite µ B have either incorporated models for the contribution from a critical point (for example, see Refs. [12][13][14][15]) or have used the Taylor-expanded equation of state from lattice QCD calculations to describe a crossover [16][17][18]. Our study is the first to incorporate both, and in addition correctly captures the universal physics near the hypothesized critical point.
The continuum limit for c 6 was published for the first time in [38], and later in [39]. In [40], a first determination of c 8 , at two values of the temperature and N t = 8 was presented. More recently, an estimate for c 8 as a function of the temperature for N t = 12 was presented in [41]. The advantage of the Taylor expansion method is that all the quantities are calculated at vanishing baryon chemical potential, where lattice QCD simulations do not suffer from the fermion sign problem. Moreover, the knowledge of the expansion coefficients can in principle provide information on the location of the critical point, under the assumption that such point is the closest singularity to µ B = 0 in the complex-µ B plane. Unfortunately, the fact that only a few coefficients are known makes this task extremely hard, leading to just an indication that the region corresponding to µ B 2 T in the phase diagram cannot contain the critical point [39]. In addition to this, the knowledge of the equation of state of QCD beyond the critical point (i.e. for µ B > µ BC ) could not come from a Taylor expansion, as a singularity cannot be reproduced in this method.
In this manuscript, we produce an equation of state which matches everything we know from lattice QCD simulations up to O(µ 4 B ) (for a recent review, see e.g. [42]) in the region where they are applicable, and which shows the correct singular behavior at and around the hypothesized critical point. The latter can be inferred from the fact that the critical point of QCD is expected to be in the same universality class as the one of the 3D Ising model [43][44][45][46][47]. Our approach is similar to the one presented in Refs. [48,49], but in our case the critical contribution is built on top of a first-principle result for the EoS, instead of relying on models such as the MIT bag or the quasiparticle model.
We will adopt the following strategy: i) Choose a location in the (µ B , T ) plane at which to put a critical point; ii) Make use of a suitable parametrization to describe the universal scaling behavior of the It is important to notice that our approach is based on the assumption that the critical point of QCD is the closest singularity to µ B = 0 on the real µ B axis. Only in this case we are allowed to merge the contributions of the lattice and 3D Ising approaches in the way detailed below. The result of this procedure will be an equation of state that meets our requirements and depends on the parameters of the non-universal map between Ising variables and QCD coordinates [50]. These parameters include the coordinates of the critical point. The ultimate goal of this project is to provide a family of model equations of state that can be used as inputs to future hydrodynamic calculations and calculations of fluctuation observables that can then be compared to experimental data from the BES-II program, resulting in constraints on the parameters in the map that we have constructed, in particular the parameters representing the location of the critical point. We will follow up on this discussion in the following sections. Note that we place a critical point in each of our equations of state entirely by construction; our calculation in isolation is therefore not a path toward determining whether the phase diagram of QCD features a critical point and where it lies. But, by comparing experimental data to predictions obtained by using our family of equations of state within future calculations of hydrodynamics and fluctuations progress toward this goal can be realized. The family of model equations of state that we construct can also be used directly to illustrate the divergence of quantities in the vicinity of the critical point, in particular the cumulants of the baryon number. We show in Section VI A our result for the second cumulant, which can be related to the variance of the net-proton distribution in heavy-ion collision experiments assuming thermal equilibrium [51]. In order to make contact with experiment, investigating out-of-equilibrium physics is important because of critical slowing down in the dynamics near a critical point and because the matter produced in a heavy-ion collision does not spend a long time near the critical point [52][53][54]. Simulations of hydrodynamics and fluctuations that are built upon the family of equations of state that we have constructed would be a good next step in this direction.

II. SCALING EOS IN THE 3D ISING MODEL
The first ingredient for this work is a parametrization of the Ising model equation of state in the vicinity of the critical point, which corresponds to a map between two variables (R, θ) to Ising variables (r, h), where r is the reduced temperature r = (T − T c )/T c and h is the magnetic field. The map needs to accommodate the correct behavior of the order parameter M (the magnetization) as a function of r and h themselves. The following form for the parametrization meets the requirements [48,49,55,56]: where Starting from this parametrization, it is possible to define the Gibbs free energy density: where F (M, r) is the free energy density, defined as: where α 0.11 is another critical exponent of the 3D Ising model (also, the relation 2 − α = β(δ + 1) holds). The function g(θ) is fixed by noticing that h = (∂F/∂M ) h and solving the following differential equation: which results in: with: Now everything is determined, and one can build an expression for the pressure in the 3D Ising model in the scaling regime, noticing that the Gibbs free energy density equals the pressure up to a minus sign: G = −P , hence: Notice that this pressure is dimensionless. Besides, this expression is completely analytic in (R, θ) in the whole range of parameter values. However, the map (R, θ) −→ (r, h) is not globally invertible.

III. NON-UNIVERSAL MAP FROM ISING TO QCD
The next step is to build a map from Ising variables to QCD coordinates, so that Eq.
(10) that we derived for the pressure becomes useful for our purpose. We want to map the phase diagram of the 3D Ising model onto the one of QCD, so that the critical point of the Ising model r = h = 0 corresponds to the one of QCD, and that the lines of first order phase transition and crossover in the Ising model are mapped onto those of QCD.
The simplest way to do so is through a linear map as follows [57]: which can be visualized in Fig At this point, we have a double map between coordinates: where the second step is globally invertible. We will now apply the thermodynamics we developed in the previous section for the Ising model, making use of the additional variables (R, θ), to the QCD phase diagram, in a parametrized form given by Eqs. (11), (12).
In order to do this analytically, we would need the map (R, θ) −→ (T, µ B ), which unfortunately cannot be globally inverted. Therefore it is necessary to solve the following relations numerically: for each value of (T, µ B ) needed in the QCD phase diagram. We proceed in the following way: we choose a range of interest for T and µ B , and given a choice of the parameters in the Ising-QCD map, we solve Eqs. (14) and (15) numerically for a two-dimensional grid in T and µ B in the desired range, thus providing a discrete inverse map (T, µ B ) −→ (R, θ).
With this solution, although not analytic, it is possible to transport the thermodynamics of the Ising model (written in terms of (R, θ)), into the QCD phase diagram, given a choice of parameters for the map. IV.

A. Strategy
The strategy we wish to pursue in order to produce an equation of state for QCD which meets the requirements stated in Section I is the following. Starting from the Taylor expansion coefficients up to O(µ 4 B ) in Eq. (I), available from lattice QCD simulations, we re-write them as a sum of an "Ising" contribution coming from the critical point of QCD, and a "Non-Ising" contribution, which would contain the regular part as well as any other possible criticality present in the region of interest: where f (T, µ B ) is a regular function of the temperature and chemical potential, with dimension of energy to the fourth power. Away from the critical regime, f just reshuffles the regular terms and can be chosen arbitrarily. Near the critical point, the choice for f is almost arbitrary, with the only requirement being that it must not add any leading singular behavior. In general, though, any term in f beyond a constant introduces sub-leading behavior in the vicinity of the critical point. For this reason, the simplest choice is to take f to be a constant, with the appropriate dimension. This also ensures that no sub-leading behavior is introduced near the critical point, which cannot be predicted through universality. Note that Eq. (16) is to be understood as a definition for the c Non-Ising n coefficients.
Once these coefficients are obtained, we will build a Taylor expansion in µ B analogous to the lattice one, using the "Non-Ising" coefficients. The latter have the advantage that the critical behavior coming from the critical point has been removed, so that the expansion can be pushed to larger values of µ B . This provides an expression for the "Non-Ising" pressure over a broad region of the QCD phase diagram. The assumption here is that the Ising critical point contribution to the Taylor coefficients from lattice QCD can be reproduced upon imposing the correct scaling behavior in the vicinity of the critical point.
Once this expansion is carried out, the full pressure is then reconstructed simply by adding the critical contribution at any (T, µ B ) to the Taylor expanded "Non-Ising" one: Note that in Eq. (17), the critical pressure is obtained from Eq. (10) with the use of the relation in Eq. (13) and the multiplication by the regular function f (T, µ B ) in Eq. (16): which is extremely easy to calculate using the above relations. We will hereafter consider the following choice for the function f (T, µ B ): The prescription we follow in order to enforce the matching of our EoS to lattice QCD at vanishing chemical potential should not be understood as a way to force a critical point in the QCD phase diagram without direct relation to the thermodynamics. Although we can, for any choice of the parameters, ensure that the coefficients at µ B = 0 coincide with lattice, this does not guarantee that the resulting EoS is thermodynamically stable. An essential step of our procedure is to check that thermodynamic inequalities are satisfied by the EoS we generate. Thanks to this, given a choice of parameters in the Ising-to-QCD map, we can ensure that the resulting EoS realizes a merging with lattice QCD without leading to pathological behavior of the thermodynamic quantities. In Section VI B (see Fig. 13) we show an example of such an analysis using our EoS when varying two of the parameters in the Ising-to-QCD map, holding the others fixed.

B. Taylor coefficients in the Ising model
The other quantities we need to calculate from the parametrization of the Ising model thermodynamics are the contributions to the expansion coefficients of the pressure, which are simply the derivatives of the latter with respect to the baryonic chemical potential at fixed temperature: Unfortunately, the expression for the critical pressure is given in terms of the additional variables (R, θ) and not as a function of (r, h) or (T, µ B ). In order to obtain the derivatives we need, we will have to use the rules for the derivative of the inverse and for the multivariate chain rule in order to be able to express everything analytically as a function of (R, θ), and convert to QCD coordinates only at the end.
We have to calculate expressions such as: which we will have to re-write as (n = 1 as an example): where: Since we do not have explicit expressions for the dependence of (R, θ) on (r, h), we need to proceed in the following way: i. Use the rule for the derivative of the inverse, so that we can express derivatives of (R, θ) wrt (r, h) as combinations of derivatives of (r, h) wrt (R, θ); ii. Use the rule for derivatives of a function holding another function constant.

Derivatives of the inverse
Naming Q n the n th derivative of an invertible function y = y(x): we can exploit the recursive relationship: where the P n are polynomials in {Q k } and P 1 = 1. P indicates derivation with respect to x.
For example, one can find:

Derivatives with functions held constant
We will have to use the following: where, in our case (x 1 , x 2 ) = (R, θ) and (y 1 , y 2 ) = (r, h).
For example, one has: .
The sequential application of Eq. (23) gives the correct expression for higher order derivatives. The explicit expressions increase in complexity extremely fast when higher order derivatives are considered, but they remain completely analytic in terms of the variables (R, θ), and allow us to have any of these derivatives defined at any point in the QCD phase diagram, provided a choice of parameters for the transformation map is given, the only step to be performed numerically being the solution of Eqs. (14), (15).

C. Critical pressure
Our procedure reduces to the use of Eqs. (10) and (20), because the dependence on (R, θ) is well defined and the numerical inversion allows us to transport any quantity to any point in the QCD phase diagram. A remark is in order at this point.
Because of the charge conjugation symmetry, in QCD the partition function needs to be an even function of the baryon chemical potential:

V. METHODOLOGY AND PARAMETER CHOICE
With the prescription exposed in Section IV C, we now have a well defined procedure to produce an expression for the pressure that meets all our requirements, needing only to make a choice for the parameters in the Ising-to-QCD map.

A. Acceptable parameter values
The choice of parameters is a key part of the procedure, because it can provide physical information, e.g. indications on the location of the critical point through comparison with experimental data. In general, the linear map we introduced has a total of six parameters.
As we will detail below, most of them are not arbitrary.
Some indication or constraint on this choice comes from our current knowledge of the QCD phase diagram. For example, since the chiral/deconfinement transition temperature is T 155 MeV at µ B = 0 [1,[58][59][60][61], and the curvature of the transition line is negative [16,62,63], we can safely expect the temperature of the critical point to be T C 155 MeV.
Other works have also shown that the presence of the critical point in the region µ B 2T appears to be strongly disfavored [39].
From the fact that the curvature of the transition line is negative and extremely small, we can easily argue that the angle α 1 in the map needs to be positive, and very small as well. The choice of the second angle is rather arbitrary, because no argument of symmetry can be made to guide the choice. For simplicity, we will consider parameter sets in which the r and h axes are orthogonal.
For the chemical potential at the critical point µ BC , although not required by any physical argument, we will restrict ourselves to values that are within reach of the BES-II program, The choice of the scale factors w and ρ is definitely less intuitive. When we consider the contribution from the pressure and its derivatives at µ B = 0, the effect of changing the scale factor w, keeping µ BC fixed, is the same as moving away or towards the critical point.
In particular, since ∂ µ B ∼ 1/w, reducing w would result in a larger "Ising" contribution to the pressure and its derivatives at µ B = 0, and a larger critical region. This can be seen in Fig. 3: a smaller value of w increases the value of the "Ising" contribution to the pressure. Moreover, the pressure grows faster both in the T and µ B directions. As for the other scale parameter ρ, its role is to govern the behavior in the pressure and its derivatives when moving away from the critical point in the temperature direction. Since we do not know what the scaling should look like in the temperature direction relative to the one in the chemical potential direction, the choice of ρ remains mainly arbitrary.

Reducing the number of parameters
In practice, although the most general linear map between Ising variables and QCD coordinates requires the use of six parameters, it is possible to impose some constraint in the choice by making use of additional arguments for the location of the critical point.
For example, the curvature of the transition line at µ B = 0 has been estimated in lattice simulations [16,62,63]. The shape of such transition line can be approximated with a parabola: where T 0 and κ are the transition temperature and curvature of the transition line at µ B = 0, respectively. The number of parameters is thus reduced to four, the angle α 1 also being fixed by: In the following, remembering that the aim of the EoS is to be employed in hydrodynamic simulations for heavy-ion collisions in the BES-II program, we will consider a choice of the baryonic chemical potential which is µ BC = 350 MeV, resulting in: In addition, the axes are chosen to be orthogonal, as we already mentioned, so that α 2 93.85 • . Finally, the scaling parameters are initially chosen as: Later we will explore different choices for w and ρ, trying to reduce their acceptable range on the basis of physical conditions for the thermodynamic quantities. For completeness, in Appendix B we show the results for the EoS obtained with other allowed parameter choices.

B. Lattice results
Recalling Eq. (16), we can see that the other defining ingredients for our procedure, besides the calculation of the Ising model thermodynamics and its "translation" to QCD, are the Taylor coefficients from lattice QCD. In the following we will use data from the Wuppertal-Budapest Collaboration [20,64] for the pressure and its derivatives at µ B = 0.
Before actually using the lattice results, we need to address a couple of issues: i. The range of temperatures of the available lattice results is not sufficient to provide an equation of state as needed in hydrodynamical simulations, namely for temperature values 30 MeV T 800 MeV; ii. The dependence of such quantities on the temperature needs to be smooth enough, such that when we take derivatives of the thermodynamic quantities wrt to T and µ B (to calculate e.g. entropy density and baryon density) they do not present an unphysical wiggly behavior.
To solve these issues, we took the following steps: i. Generate data for temperatures below the reach of lattice (T ≥ 135 MeV in this case) using the HRG model; ii. Provide a parametrization of the temperature dependence of the pressure and its derivatives in the desired temperature range.
The parametrizations of χ 0 (T ) and χ 4 (T ) were performed through a ratio of 5 th order polynomials in the inverse temperature: while for χ 2 (T ), a different expression was used: The values of the parameters are given in the following tables: In Fig. 4 we can see the comparison between the lattice data (and the extension with the HRG model) and the resulting parametrization. The HRG model employed to calculate the pressure does not contain any interaction, and makes use of the most up to date particle list available from the Particle Data Group [66] (list PDG2016+ in [67]).
The smooth curves obtained from the parametrization will be the c LAT

VI. RESULTS
At this point, we have all the ingredients in Eq. (17): which is now straightforward. However, although the overall behavior is correct, at low temperatures and in particular in regions where the ratio µ B /T is very large, the pressure becomes negative, and so do other observables as well. This is due to the fact that, given our way to smoothly merge the pressure coming from the procedure we developed so far with the one from the HRG model.
The smooth merging can be obtained through a hyperbolic tangent as: where T (µ B ) works as the "switching temperature", and ∆T is roughly the size of the "overlap region" where both pressures contribute to the sum. The dependence on the baryon chemical potential of the "switching temperature" is chosen to be parabolic, and parallel to the chiral transition line we assumed in Eq. (26):  and speed of sound normalized by the correct power of the temperature: lation for the speed of sound directly from the definition; however, it is possible to re-write the expression for this observable in terms of derivatives of the pressure with respect to the temperature or the chemical potential only [68]: In Figs We note that, at large values of µ B , small wiggles appear in the thermodynamic observables, in particular the speed of sound. This is due to the truncation in the Taylor expansion of the non-Ising contribution to the pressure. In general these wiggles are only appearing at the edge of the allowed µ B region and only for some parameter choices (see Appendix B for other choices in which they are not present). Besides, their magnitude is very small so we are confident that they will not affect the results of the hydrodynamic simulations.
However, in order to extend our EoS to larger values of µ B , we need to devote future work to incorporate higher order terms in the Taylor expansion into the construction of model equations of state in order to improve their behavior at higher µ B .
Although not very evident from the pressure, the critical point manifests itself clearly in first order derivatives (entropy, baryon and energy density), where the discontinuity due to a first order phase transition is clearly visible at µ B > µ BC . Furthermore, the speed of sound shows a clear dip at the critical point, as well as a (less evident) discontinuity at µ B > µ BC . Figure 11 shows the trajectories at constant S/n B in the QCD phase diagram.
In addition to the thermodynamic quantities mentioned above, it is possible to calculate observables that are more sensitive to critical fluctuations: in Fig. 12 we show the second cumulant of the baryon number for the choice of parameters in Section V A 1.
In heavy-ion collision experiments, it is possible to measure related quantities, constructed from moments of the event-by-event distribution of the measured number of protons in a given acceptance. The critical contribution to χ B 2 diverges at the critical point like ξ 2 , where ξ is the correlation length of the order parameter fluctuations. Higher, non-Gaussian, moments of the proton number distribution (or the baryon number distribution) receive larger contributions from critical fluctuations, with the third and fourth cumulants diverging like ξ 9/2 and ξ 7 respectively [69], making these observables more sensitive to the presence of a critical point [70][71][72][73], and motivating their experimental measurement in the RHIC Beam Energy Scan [74][75][76][77][78], as well as the future investigations of higher moments of the n B distribution in this model.

B. Exploration of parameter space
By requiring thermodynamic stability, i.e. positivity of pressure, entropy density, baryon density, energy density and speed of sound, and causality, i.e c 2 s < 1, over the whole phase diagram, it is possible to reduce the range of acceptable parameters in the non-universal Ising → QCD map. By keeping the location of the critical point fixed (µ BC = 350 MeV, T C 143 MeV), as well as the orientation of the axes (α 1 3.85 • , α 2 − α 1 = 90 • ), we investigated the role of the scaling parameters w, ρ. In Fig. 13, we can see in red the points corresponding to pathological parameter choices, while the blue dots correspond to acceptable ones. We notice that, while most commonly specific parameter choices are unacceptable because of the negativity of n B , for very low w (w = 0.25) we observe violation of causality as well (c 2 s > 1).  Higher order cumulants, which are expected to display a more pronounced divergence at the critical point, will be investigated in future work. We also leave to future work the inclusion of strangeness and electric charge chemical potentials, which play an important role in heavy-ion collision physics. The main reason for this is that the same universality arguments we employ in this work cannot be used to answer the question how the critical behavior in the (T, µ B ) plane would vary with the two additional chemical potentials. Because of this, such an exercise would require a great deal of additional modeling, which is beyond the scope of this work. The following relations will be used in Eq. (23): and hence:

Second order
The relationships start to become more complicated at the second order.
Derivatives of G(R, θ) wrt (R, θ): Derivatives of G(R, θ) wrt (T, µ B ) are constructed as: where the terms with second derivatives of (r, h) wrt (T, µ B ) have been dropped, since the transformation between the two sets is linear.
Derivatives of G(R, θ) wrt (r, h): where the expression for the derivatives of (r, h) wrt (R, θ) are already quite long.

APPENDIX B
In this appendix we present a few more results for the EoS, obtained with different parameter choices in the allowed ranges discussed in the main text. We mainly aim at exploring the effect of a different size for the critical region and a different location for the critical point. Figure 14 shows our EoS for the same location of the CP and angles presented in the main text, but for w = 4 and ρ = 1. This choice was made to show that the size of the critical region can be made arbitrarily small by increasing w. This is evident by the fact that the discontinuity in entropy, energy and baryonic density is much reduced in this case, as well as the dip in the speed of sound and peak in χ 2 . This is evident by the fact that the discontinuity in entropy, energy and baryonic density is larger in this case, as well as the dip in the speed of sound and peak in χ 2 .
In Figure 16 we change the location of the CP to µ BC =400 MeV, the angles become α 1 4.40 • , α 2 − α 1 = 90 • , and we set w = 2 and ρ = 2. When the data from the second beam energy scan at RHIC become available, it will be possible to perform a Bayesian analysis to systematically scan the parameter region (in a similar fashion as presented in this appendix, but for many more parameter sets) and test the results of hydrodynamic simulations with this family of EoSs as input against the data. This will hopefully help to constrain the size of the critical region, and the location of the critical point.
T and µ B . All panels correspond to the same location of the CP and angles presented in the main text, but with w = 4 and ρ = 1.

APPENDIX C
A. Treatment of the phase coexistence region The first order phase transition appearing at chemical potentials larger than the critical one, leads to a phase coexistence region in the plane of pressure vs net-baryon density. The knowledge of the EoS in such a region, even though it does not correspond to a stable state of the system, can be needed in hydrodynamic simulations [79][80][81][82][83].
Due to the fact that our treatment is carried out in the temperature vs chemical potential plane, our EoS does not contain information about the coexistence region. In order to use the standard Maxwell construction, we would need a thermodynamic potential defined even in the un-physical region, but this is unfortunately unavailable in our approach. Thus, we take a phenomenological approach to reconstruct the EoS in this region. In Fig. 17 the isothermal curves in the pressure vs baryon density plane are shown, with the same choice of parameters as in Section V A 1: the solid blue line corresponds to a temperature T > T C and shows only a slight inflection, while the solid black lines correspond to temperatures T < T C , and display a clear discontinuity; the dashed red line corresponds to T T C , and shows the insurgence of the discontinuous behavior.
Since the thermodynamic conditions for the coexistence of two phases (thermal, chemical and mechanical equilibrium) are all fulfilled by construction in our procedure, what is needed is a continuous curve connecting the two branches of each discontinuous isotherm. The phenomenological approach we follow here is to construct a form that is similar to many other analytic theories (e.g. real gases) and presents the typical double lobe structure.
We perform a polynomial fit to the two branches of each isotherm and the middle point in the discontinuity. In the left panel of Fig. 18 the original, discontinuous isotherm at T = 142MeV is shown (black solid), along with the points used to generate the curve in the coexistence region (dark pink) and the resulting function inside the coexistence region (pink dashed).
The results of this treatment are shown in the right panel of Fig. 18: the gap in the (P, n B ) phase diagram corresponding to the coexistence region between two phases is now filled by a phenomenological dependence of P (n B ) along the isotherms. This will allow one to use our EoS in hydrodynamic codes inside the phase coexistence region emerging from the presence of a first order chiral/deconfiment phase transition.