Description of the first order phase transition region
of an equation of state for QCD with a critical point

We map the Ising model phase diagram (in terms of reduced temperature $r=(T-T_{c})/T_{c}$ and magnetic field $h$) onto the QCD one (in terms of temperature $T$ and baryonic chemical potential $\mu_{B}$), with the goal of introducing a critical point and first-order phase transition on top of the equation of state from lattice QCD, similarly to what was done in Refs. Parotto et al. (2020); Karthein et al. (2021); Kapusta and Welle (2022); Kahangirwe et al. (2024). Our main goal is to identify the spinodal lines in the QCD phase diagram as mapped from the Ising model, and provide a full description of the equation of state around the first-order phase transition.

When the transition is first order, the spinodal region is a distinctive feature of the QCD phase diagram, where quantities like entropy density, baryon density or energy density are multi-valued functions of the temperature at fixed chemical potential. Typically, only one of these branches is stable, while the others are metastable or unstable. When considering isothermal trajectories in the phase diagram, the spinodal points, or spinodal singularities, demarcate the limits of metastability Fisher (1967); Binder (1987). Furthermore, the spinodals are related to the Lee-Yang edge singularities in the Ising model, as studied in Ref. An et al. (2018). In the mean field Ising model, the spinodals are found at real values of the magnetic field $h$, i.e. they lie on the real $h-$axis. On the other hand, beyond mean field they are displaced from the real $h-$axis by an amount An et al. (2018):

$$\Delta\phi=\pi(\beta\delta-\frac{3}{2}),$$

(1)

where $\beta$ and $\delta$ are the critical exponents in the Ising model. When the critical exponents are different from their mean field values of $\beta=1/2,~{}\delta=3$, this phase is non-zero. As is the case in the 3D Ising model, the spinodal points are pushed into the complex plane, owing to the importance of fluctuations which are neglected in the mean field approach. The angle by which they are shifted into the complex plane is given by the difference between $(\beta\delta)_{3\mathrm{D}}=1.5648$ and the mean field value of 3/2, yielding $0.0648\pi$. Given the smallness in the difference of these critical exponents, we take this as a motivation to proceed with the mean field approximation. Later in this section, we shall introduce the Ginzburg criterion and discuss the range of applicability of mean field.

The mean field equation of state is given by the relationship between the free energy (Helmholtz $F$ or Gibbs $G$), magnetization ($M$), and temperature ($r=(T-T_{c})/T_{c}$). The external magnetic field ($h$), used as a control parameter, is given in terms of $M$ and $r$. In addition, the pressure ($P$) is equivalent to the Gibbs free energy up to a minus sign. The Helmholtz free energy, magnetic field, and pressure are defined as:

	$\displaystyle\centering F(M,r)\@add@centering$	$\displaystyle=h_{0}M_{0}\Big{(}\frac{1}{2}arM^{2}+\frac{1}{4}bM^{4}\Big{)}$		(2)
	$\displaystyle h(M,r)$	$\displaystyle=\Bigg{(}\frac{dF}{dM}\Bigg{)}_{r}=h_{0}(arM+bM^{3})$		(3)
	$\displaystyle P(M,r)$	$\displaystyle=-G(M,r)=Mh-F(M,r)$		(4)
		$\displaystyle=h_{0}M_{0}\Big{(}\frac{1}{2}arM^{2}+\frac{3}{4}bM^{4}\Big{)},$

where $a=3,~{}b=1$. This is consistent with the linear parametric model of the scaling EoS in terms of the parametric variables $(R,\theta)$ Bzdak et al. (2020); Schofield (1969); Guida and Zinn-Justin (1997). The normalization constants for the magnetization and the magnetic field are determined by imposing $M(r=-1,h=0^{+})=1$ and $M(r=0,h)\propto\mathrm{sgn}(h)|h|^{1/\delta}$ and found to be $h_{0}=\frac{1}{3\sqrt{3}},~{}M_{0}=\frac{1}{\sqrt{3}}$ Bzdak et al. (2020); Landau and Lifshitz (1980); Nonaka and Asakawa (2005). We note here that the choice of mean field Ising model allows us to forgo the formulation of the equation of state in terms of parametric variables in order to solve the system analytically with rational exponents. Furthermore, this gives us direct access to the Ising variables $(r,h)$, such that we have only a single mapping procedure to perform to QCD variables $(T,\mu_{B})$, unlike with the 3D Ising model scaling equation of state which relies on the $(R,\theta)$ parametrization Parotto et al. (2020); Karthein et al. (2021); Kahangirwe et al. (2024).

Figure 1 shows the Ising free energy for a given temperature $r$ as a function of the fixed magnetization from Eq. (2). The left-hand plot shows the free energy for temperatures above the critical point one, i.e. $r>0$. This corresponds to the crossover transition in the Ising model. In this case, we have a simple parabolic shape as both quadratic and quartic terms have positive coefficients. As the critical point is approached, the curvature at the minimum is reduced until it vanishes and the curve becomes completely flat. From this zero curvature feature at the critical point, the coexistence region emerges below the critical temperature, $T<T_{c}$ or $r<0$. The right-hand plot in Fig. 1 shows these features for $r<0$, corresponding to the first order phase transition. Here, the double-minimum character is apparent, due to the negative sign of the quadratic term in Eq. (2) when $r<0$. Furthermore, we can characterize the typical features of an isothermal curve passing through a first order phase transition: the coexistence and spinodal points. The coexistence points are the minima, $(\partial F/\partial M)_{r}=0$, separated by a free energy barrier, while the spinodals are the inflection points $(\partial^{2}F/\partial M^{2})_{r}=0$. Moreover, since $h=(\partial F/\partial M)_{r}$, it is clear that the coexistence points correspond to $h=0$. Thus, the spinodals are defined as $(\partial h/\partial M)_{r}=0$.

The two minima of equal magnitude of the free energy on the first order phase transition side imply that there are two phases which are equally probable and thus coexist. The value of magnetization $M$ and temperature $T$ of these minima define the coexistence line in the $(M,~{}T)$ plane. Phase coexistence arises at the point in the system evolution where one phase begins to change into the other. Because these phases occur together, this is known as the coexistence region. In a typical approach, one might employ a Maxwell construction between the coexistence points in order to work with a system in complete thermodynamic equilibrium. This would correspond to a straight line connecting the two coexistence points, along which the order parameter flips sign. However, we have further information in this mean field Ising approach: the spinodal points. As shown here, the spinodals are the inflection points, with $\partial^{2}F/\partial^{2}M=0$. By further studying the relationship between the magnetization and magnetic field, we shall discover the properties leading to the spinodal points.

Indeed, we shall further characterize the mean field equation of state by the behavior of the magnetization $M(h)$. In the case of the mean field Ising equation of state, as shown in Eq. (3), $h$ is a function that is defined in terms of $M$ and $r$. Thus, we must invert the equation $h(M,r)$, leading to three solutions for this cubic function. To do so, we write $h(M,r)=M(M+i\sqrt{3r})(M-i\sqrt{3r})$ yielding one real solution and an additional pair of complex solutions for $M(h,r)$:

$$\centering M(h,r)=\begin{cases}\dfrac{2^{1/3}}{\eta^{1/3}}r-\dfrac{\eta^{1/3}}{2^{1/3}}\\ -\dfrac{(1\pm i\sqrt{3})r}{2^{2/3}\eta^{1/3}}+\dfrac{(1\mp i\sqrt{3})\eta^{1/3}}{2^{4/3}},\end{cases}\@add@centering$$

(5)

where we define $\eta=\sqrt{h^{2}+4r^{3}}-h$. In order to access these features in the QCD phase diagram, we require a map between $(r,h)\rightarrow(T,\mu_{B})$. Therefore, it is necessary for us to invert $h(M,r)$ in favor of $M(h,r)$. Thus, we ensure that our method allows us to map $h$ as a variable onto the QCD phase diagram. In our approach, we choose to work within the $T-\mu_{B}$ phase diagram in order to utilize the first-principles input to the equation of state for strongly-interacting matter from lattice QCD Bellwied et al. (2015). The alternative choice would be the $T-n_{B}$ phase diagram Randrup (2010); Vovchenko et al. (2015). The analogous variables in the Ising system would, thus, be reduced temperature $r$ and magnetization $M$. Given this choice, $M$ is retained as a control variable allowing one to use the equation of state defined in Eq. (3), which provides a well-defined, single-valued function for $h(M,r)$. The behavior of this quantity can be inferred from Fig. 2. By rotating the figure, or equivalently switching the axes, one can see $h(M)_{r}$ is a single-valued function. To reiterate, we remain with the choice of $h$ as it is analogous to $\mu_{B}$ and allows us to utilize the constraints on the equation of state at vanishing $\mu_{B}$ from lattice QCD.

In Fig. 2, the behavior of $M(h)$ at fixed $r$ above and below the critical point is shown on the left- and right-hand plots, respectively. As can clearly be seen in the left-hand plot, only one solution is real on the crossover side where $r>0$, corresponding to the solution which is always real in Eq. (5). On the other hand, the interesting features of the first order phase transition are apparent for $r<0$, where the complex pair of solutions become real for certain values of $h$ at fixed $r$. Here we see the multi-valued behavior for this quantity that we expect from these three real solutions. Each of the three solutions for the magnetization from Eq. (5) is shown in a different color. The blue corresponds to the solution that is real across all temperatures, while the yellow and orange are the complex conjugate pair of solutions which become real when $r<0$. Both the coexistence points and spinodal points defined above as $h=0$ when $(\partial h/\partial M)_{r}>0$ and $(\partial h/\partial M)_{r}=0$, respectively, appear naturally here. By looking at this quantity in particular, we see that away from the coexistence points, i.e. away from $h=0$, the system is either in one phase or the other, given by the sign of the magnetization. On the other hand, as the system approaches $h=0$, moving through the coexistence points and then proceeding further to the spinodal points, matter becomes either super-heated or super-cooled. Between the coexistence and spinodal points we have this extreme environment for the phases where the system is in a metastable state. Then, beyond the spinodal points the system becomes unstable due to the negative slope, i.e. $(\partial M/\partial h)<0$. This is analogous to the thermodynamic instability $(\partial n_{B}/\partial\mu_{B})<0$ as encountered, for example, in the well-known liquid-gas phase transition.

The third solution exists only in this unstable region occurring between the two spinodal points along an isotherm. Therefore, we see here that the yellow and blue solutions describe the system up to metastability, and beyond the limit of metastability the final solution shown in orange is needed, corresponding uniquely to the instability $(\partial M/\partial h)<0$. Upon further consideration of this isotherm, we note that the Maxwell construction corresponds to proceeding from one phase to the other via the coexistence points at $h=0$. For example, in the case of negative magnetization, the system follows the bottom curve shown here in yellow until it reaches the coexistence point given by the open green circle, after which it proceeds through the phase transition along the $h=0$ line until the next coexistence point where it ends up on the solution shown here in blue. The phase has changed due to the change in sign of the magnetization, but this approach has completely disregarded the multi-valued nature of this quantity. We, thus, make use of the full information in the coexistence region by describing the metastable and unstable regions. For the rapidly expanding strongly-interacting system in heavy-ion collisions, it is possible to reach these meta- or un-stable regions.

We then turn to the behavior of the pressure in the Ising model. The quantity of interest for our purpose of studying these features in the QCD phase diagram is, in fact, the pressure. We, thus, utilize the expressions for the magnetization from Eq. (5) and obtain the pressure, Eq. (4), as a function of the external variables $r$ and $h$. Moreover, this allows us to track the spinodals in a straightforward way and study how they map onto the QCD phase diagram in temperature and baryon chemical potential. In this way, by utilizing the mean field approach, we preserve the information on multi-valued observables.

The pressure in the crossover ($r>0$) and first order ($r<0$) regions is shown in Fig. 3. Because the pressure relies on the form of $M(h)$, the presence of a single or three real solutions above or below the critical point is manifest in the pressure itself. In the case of the first order phase transition (Fig. 3, right), the coexistence point, where the pressure of the two phases is the same, corresponds to this location where the two solutions meet. The spinodal points, on the other hand, are located where the super-heated/super-cooled phases terminate, with the subsequent onset of the unstable region. The third solution describes this region between the spinodal points. Generally, the phase transition proceeds by either nucleation or spinodal decomposition, when the system is in the metastable or unstable state, respectively Binder (1987); Cahn and Hilliard (1959); Abyzov and Schmelzer (2007). Thus, between the spinodals the system likely undergoes spinodal decomposition Chomaz et al. (2004); Randrup (2009); Steinheimer and Randrup (2012); Kapusta et al. (2024), or possibly nucleation Cahn and Hilliard (1959); Mishustin (1999); Csernai and Kapusta (1992) but only when the nucleation rate is large enough compared to the expansion rate of the medium.

We study the pressure in the full Ising phase diagram in $h$ and $r$ in order to see its structure more clearly. The three solutions for the pressure are shown in terms of the Ising variables $(r,h)$ in the top panel of Fig. 4. This allows us to visualize the behavior of the pressure on the low ($r<0$) and high ($r>0$) temperature sheets. On the high temperature sheet, where the transition is a smooth crossover, only one solution exists which corresponds to the real root of $h(M,r)$, while on the low temperature sheet, where the transition becomes first order, the two complex roots also become real giving rise to three distinct planes for the pressure. In this figure, the three solutions for the pressure following from $P(h,r)$ that one obtains applying Eqs. (5) to (4) are shown in different colors with the same scheme applied as in the isotherms. The coexistence curve corresponds to the line along which the first two solutions, shown in yellow and blue, meet for $h=0,~{}r<0$, shown here as the dashed green line. On the other hand, the spinodals are located where the first two solutions meet the third one, shown here as the solid orange lines. Thus, as in the isotherm shown in Fig. 3, the spinodals are the lower edges of the triangle outlining the coexistence region. When crossing the $r=0$ line, i.e. the critical isotherm, the pressure collapses to one solution above $T_{c}$, i.e. on the high temperature sheet. Therefore, when moving along a trajectory in the $h-r$ plane for which $r$ is not fixed (which we will see later corresponds to a QCD isotherm), we will end up crossing the critical isotherm and going onto the low temperature sheet. It is important to note that there is no phase transition associated with crossing the critical isotherm.

If we consider the behavior in the complex $h$-plane on each of the low and high temperature sheets separately, we can understand the presence of the spinodal lines in the phase diagram. In the middle panel of Fig. 4, we show the pressure in the complex $h$ plane for $r<0$. We see that the Lee-Yang edge singularities are located along the real $h$-axis in this case. Conversely, for the high temperature $r>0$ sheet, these manifestations corresponding to the spinodals on the other sheet lie completely along the imaginary axis (bottom panel of Fig. 4). Beyond mean field, these spinodal singularities no longer lie along the real or imaginary axes, but rather move into the complex plane.

In the present study, we are concerned with characterizing the thermodynamical aspects of the first order phase transition via the equation of state. Moreover, because the equation of state serves as crucial input to the hydrodynamic simulations of heavy-ion collisions, we require full information on the baryon density and energy density as functions of the temperature and chemical potential. In the first order region, the baryon and energy densities become multi-valued solutions of $T$ and $\mu_{B}$, where the phases coexist. Therefore, the simulations need the ability to distinguish the different branches of the baryon and energy density in to order establish where precisely the system lies in terms of $T$ and $\mu_{B}$ within the coexistence region. This is exactly the behavior that we capture by utilizing the mean field Ising equation of state, as shown by the magnetization and pressure in the Ising model in Figs. 2 and 3. By having explicit control over the three real solutions appearing for $r<0$, we are able to characterize the nature of the phase coexistence. Thus, by translating these features into QCD variables, this equation of state framework allows the hydrodynamical simulations to track the phases in the first order region. This is particularly necessary when modeling low energy collisions where a first order transition may be present. However, it is important to note that the spinodals are valid in the regime of metastability which has a finite lifetime, and as such, they do not belong to a system in thermal equilibrium Langer (2000); Binder (1987). They are, nonetheless, manifest from the mathematics as properties of the equation of state, belonging to the broader class of Lee-Yang edge singularities Lee and Yang (1952). Furthermore, spinodal decomposition (together with nucleation) helps to describe the dynamics of real physical systems near a first order phase transition. From an applied perspective, spinodal decomposition is ubiquitous in industrial applications Jones (2002). The importance of understanding the spinodals is, therefore, clear, and we proceed to produce and study an equation of state for QCD that contains these features of a first order phase transition. The scope of this work is to study how these features map onto the QCD phase diagram, and as such we do not consider dynamical aspects owing to the finite lifetime of metastability.

Before mapping the features of the Ising model onto QCD, we first discuss the range of validity of the mean field approximation which we are utilizing in this study. By applying a mean field approach, one inherently assumes that the average behavior is sufficient to describe the system and that fluctuations can be neglected. The well-known Ginzburg criterion establishes the regime of validity for mean field where the fluctuations must be smaller than the magnitude of the order parameter itself Ginzburg (1960). This idea was taken a step further by Binder Binder (1987) who showed that, in the case of the spinodal singularity, the interaction range, $R$, must be sufficiently large for systems with dimensionality below the upper critical dimension, to apply mean field treatment:

$$1\ll R^{d}(h_{c}-h)^{(6-d)/4}.$$

(6)

Indeed, the real space Ginzburg criterion as shown by Als-Nielsen and Birgeneau Als-Nielsen and Birgeneau (1977) can be written as:

$$\chi_{M}\ll VM^{2},$$

(7)

where $\chi_{M}=\partial h/\partial M$ is the magnetic susceptibility and $V$ is the volume, maximally determined by the size of the correlation length $\xi^{d}$. By utilizing the well-known critical scaling Zinn-Justin (2002) for the susceptibility $\chi_{M}\sim r^{-\gamma}$, order parameter, $M\sim r^{\beta}$, and correlation length, $\xi\sim\xi_{0}r^{-\nu}$, the Ginzburg criterion becomes:

$$r^{-\gamma+d\nu-2\beta}\ll\xi_{0}^{d}.$$

(8)

Inside this region, the mean field approximation begins to breakdown as fluctuations become important. This is consistent with what is found in Ref. An et al. (2018) in terms of the behavior of the scaling-invariant variable $w=ht^{-\beta\delta}$. In the case of the 3D Ising model, when considering the current level of precision of the determination of the critical exponents, the left-hand side becomes $r^{0}$.

This leads to the conclusion that the size of the system determined by the size of the correlation length must be much greater than 1. This large (diverging) correlation length indicates the presence of the critical point singularity. If the correlation length is taken to be in units of fm, a correlation length of several fm has been shown to be large enough to see the critical effects in a system with finite size and lifetime Berdnikov and Rajagopal (2000). Thus, in the regime of large correlation lengths, corrections to mean field may become relevant. On the other hand, this formulation of the Ginzburg criterion shows that mean field is an appropriate description when the correlation length is small. Additionally, we will compare the overall equation of state between the 3D and mean field Ising models in Sec. IV.

In order to determine the location of the spinodal points within the QCD phase diagram after the mapping procedure, we derive the dependence of the spinodal location on the mapping parameters. Firstly, we begin with the spinodal definition in the mean field Ising model, consistent with what has been shown in Ref. An et al. (2018)

$$h_{\mathrm{sp}}=\pm\frac{2r^{3/2}}{3\sqrt{3}},$$

(10)

which is obtained by solving the equation $(\partial h/\partial M)_{r}=0$. Here, we note that we only consider the spinodals in their region of validity, i.e. only for $r<0$, which leads to the absorption of the “$i$” appearing in the more general formulation of Lee-Yang edge singularities in Ref. An et al. (2018). Now, plugging in the mapping equations $h(T,\mu_{B}),~{}r(T,\mu_{B})$ with orthogonal angles $\alpha_{1},~{}\alpha_{2}$ leads to the following expression:

		$\displaystyle\frac{(\mu_{B}-\mu_{B,c})\sin\left(\alpha_{1}\right)+(T-T_{c})\cos\left(\alpha_{1}\right)}{T_{c}w}=$		(11)
		$\displaystyle\pm\frac{2}{3\sqrt{3}}\Bigg{(}\frac{\left(-(\mu_{B}-\mu_{B,c})\cos(\alpha_{1})+(T-T_{c})\sin\alpha_{1}\right)}{\rho T_{c}w}\Bigg{)}^{3/2}.$

We choose orthogonal angles to eliminate the additional dependence on the angle $\alpha_{2}$ and study the effect of the angle $\alpha_{1}$. We leave the study of relaxing the orthogonality condition and exploring the effect of $\alpha_{2}$ to future work. In order to study the spinodal location in $\mu_{B}$ in a more closed-form expression, we can consider the two limiting cases of the angle $\alpha_{1}=0^{\mathrm{o}},~{}90^{\mathrm{o}}$, which correspond to directly mapping the Ising phase transition line $h=0$, i.e. the $r$-axis, along lines parallel to the QCD $\mu_{B}$-axis or $T$-axis, respectively. The former case is more appropriate for QCD at small values of the chemical potential since we know that the chiral phase transition line is very flat in this region, giving rise to a small curvature or angle Bellwied et al. (2015); Borsanyi et al. (2020). On the other hand, at large values of the chemical potential, the transition line should approach the $\mu_{B}$-axis orthogonally Halasz et al. (1998), corresponding to the latter case. It is also important to note that the case of $\alpha_{1}=90^{\mathrm{o}}$ corresponds to mapping analogous quantities onto one another, i.e. $h\rightarrow\mu_{B},~{}r\rightarrow T$. Conversely, for vanishing $\alpha_{1}$ it is the opposite, $h\rightarrow T,~{}r\rightarrow\mu_{B}$. Thus, when $\alpha_{1}=0$, fixing the QCD temperature $T$ corresponds to constant magnetic field lines rather than constant temperature in Ising variables. In other words, an isotherm in the Ising model now corresponds to a trajectory at fixed baryon chemical potential $\mu_{B}$. This also applies to the case of small $\alpha_{1}$ and will be important for understanding the behavior of the pressure mapped to $T$ and $\mu_{B}$.

We begin with $\alpha_{1}=0^{\mathrm{o}}$ and the positive solution for $h_{\mathrm{sp}}$ from Eq. (10). In this case, the location of the spinodal in $\mu_{B}$ corresponds to:

$$\mu_{B,\mathrm{sp}}^{\alpha_{1}=0}(T)=\mu_{B,c}+\frac{3}{2^{2/3}}\rho w^{1/3}T_{c}^{1/3}(T-T_{c})^{2/3}.$$

(12)

In fact, this expression is the same regardless of which $h_{\mathrm{sp}}$ solution we begin with. Thus, this already shows that we will only find a single spinodal point in $\mu_{B}$ for this choice of angle. In order to further illuminate this situation for the choice of small angle, it is important to note that the critical isotherm ($r=0$) rather than the coexistence line ($h=0$) is being crossed along the QCD isothermal trajectories. In other words, if we considered constant $\mu_{B}$ trajectories we would find the two solutions we expect corresponding to each of the $h_{\mathrm{sp}}$ solutions.

On the other hand, if we consider the case with $\alpha_{1}=90^{\mathrm{o}}$, we find the spinodal locations to be:

$$\mu_{B,\mathrm{sp}}^{\alpha_{1}=\pi/2}(T)=\mu_{B,c}\pm\frac{2\sqrt{3}}{9}\rho^{-3/2}w^{-1/2}T_{c}^{-1/2}(T_{c}-T)^{3/2}.$$

(13)

Here, the sign corresponds to the sign of the $h_{\mathrm{sp}}$ solution. Thus, we are able to determine the location of both spinodals for a given QCD isothermal trajectory in the case of $\alpha_{1}=90^{\mathrm{o}}$.

These limiting cases help to understand what we can expect to see for the spinodals in terms of $T$ and $\mu_{B}$. We see that because the spinodals in the Ising model are defined by $(\partial h/\partial M)_{r}=0$, at fixed Ising temperature $r$, the way $r$ maps to $T$ is very important for isothermal trajectories. By mapping the Ising coexistence line tangential to the chiral crossover line, there will generally be a mixing between the Ising variables as mapped onto $T$ and $\mu_{B}$, since the angle $\alpha_{1}\neq 90^{\mathrm{o}}$. Thus, isotherms from the Ising phase diagram become straight lines in the QCD phase diagram characterized by some value of $\mu_{B}/T$.

Next, we study isothermal curves for the Ising pressure as mapped onto the QCD phase diagram. As can be seen from the derivation of the spinodal location in $\mu_{B}$ when the angle $\alpha_{1}$ is small, a QCD isotherm will only capture one of the spinodals while the other lies in a region of the phase diagram not reached by that isotherm. The behavior of the Ising pressure in its own phase diagram, as shown in Fig. 4, becomes important for understanding this effect. As can be seen in this 3D representation, for any Ising isotherm $r<0$, both edges of the spinodal lines will be reached by such an isothermal curve. On the other hand, once the Ising variables are traded for the QCD ones, an isotherm for $T<T_{c}$ may not necessarily capture both sides of the triangular structure for all mapping parameters. This is especially true for small angles, where the mapping is such that $h\rightarrow T$ and $r\rightarrow\mu_{B}$, as discussed previously. Furthermore, as the mapping in Eq. (9) shows, $h$ and $r$ can mix in the expressions for $T$ and $\mu_{B}$. This means that these triangular structures in the pressure can become distorted along the isotherms. In order to elucidate this isothermal behavior in QCD variables, Fig. 9 shows the spinodal and coexistence curves as mapped onto the QCD phase diagram. From these figures, it is clear that for smaller angles, or equivalently smaller values of $\mu_{B,c}$, a line of constant $T$ cannot capture both spinodal lines. This is due to the fact that the spinodal line above the phase transition line, which represents super-heated hadronic matter, opens up into the phase diagram towards higher temperatures. These results will be discussed in detail in Sec. III.3.

We attempted to obtain access to the spinodal region for the same choice of mapping parameters from Ref. Parotto et al. (2020), which determine the location of the critical point ($\mu_{B,c}$ from which we obtain $T_{c}$ and $\alpha_{1}$) as well as the size and shape of the critical region ($w,~{}\rho,~{}\alpha_{2}$). The location in this example choice is $\mu_{B,c}=350$ MeV with scaling parameters $w=1,~{}\rho=2$ and orthogonal Ising axes such that $\alpha_{2}=\alpha_{1}+90^{\mathrm{o}}$. In Ref. Parotto et al. (2020), the chiral phase transition line was utilized up to $\mathcal{O}(\mu_{B}^{2})$, which results in a very small mapping angle of $\alpha_{1}=3.8^{\mathrm{o}}$. In this work, we utilize updated results from the Wuppertal-Budapest collaboration on the chiral phase transition line up to $\mathcal{O}(\mu_{B}^{4})$ from lattice QCD Borsanyi et al. (2020):

$$\frac{T_{c}(\mu_{B}/T_{c})}{T_{c}(\mu_{B}=0)}=1-\kappa_{2}\Big{(}\frac{\mu_{B}}{T_{c}(\mu_{B})}\Big{)}^{2}-\kappa_{4}\Big{(}\frac{\mu_{B}}{T_{c}(\mu_{B})}\Big{)}^{4}.$$

(14)

In Fig. 6 we show the results for this curve with all errors included as a band in blue. As the chemical potential increases, the lattice results become less constraining. As such, we choose a parametrization of the chiral phase transition line within the range of values estimated by the lattice calculations, shown as the black, solid line in Fig. 6. This choice is motivated by obtaining a curve that bends steeply enough to yield larger angles $\alpha_{1}$ as $\mu_{B}$ increases. Any such parametrization would be acceptable, as long as the curve terminates at baryon chemical potentials larger than that of the proton mass where we know confined hadronic matter exists. Additionally, the transition line should bend down to touch the $\mu_{B}$-axis and become orthogonal at $T=0$ due to the third law of thermodynamics Halasz et al. (1998). As in Eq. (14), we use a fourth order polynomial in order to parametrize our fit curve. As such, we do not expect to capture the singularity at $T=0$, where $\partial T/\partial\mu_{B}$ becomes infinite.

Given our parametrized fit curve, we find a different angle $\alpha_{1}$ for $\mu_{B,c}=350$ MeV than Ref. Parotto et al. (2020). In addition, the value of $T_{c}$ for a given choice of the critical chemical potential, $\mu_{B,c}$, will also be adjusted. In this case, our updated parameters which are fixed by the placement of the critical point along the chiral phase transition line are $\alpha_{1}=4.6^{\rm{o}}$ and $T_{c}=146$ MeV for the choice $\mu_{B,c}=350$ MeV.

Fig. 7 shows a QCD isotherm for the pressure in the case of such a shallow mapping of $\alpha_{1}=4.6^{\rm{o}}$ along the chiral phase transition line. First, the mapping of the features of the first order phase transition from the Ising model to QCD should be understood. In this figure, the three solutions for the pressure, coming from the three solutions for $M(h)$ as shown in Eq. (5), are shown in different colors along with the value of temperature to which this isotherm corresponds. The coexistence region begins where the first two solutions, shown here in yellow and blue, meet at $\mu_{B}>\mu_{B,c}$, while the spinodal points are located where the first two solutions meet the third one shown in orange (see Fig. 3). In this case of a small angle, we do not capture the point where the second and third solutions meet, since this is at a different temperature (see discussion in Sec. III.1). Here, we chose an isotherm that is close to the critical point, $T=0.95T_{c}=140$ MeV, in order to depict the behavior of the spinodals, particularly to point out that the shape of the isotherms becomes distorted (compare to Fig. 3). In fact, the rightmost spinodal becomes stretched out, moving through the phase diagram to another value of $T$ and $\mu_{B}$.

After showing the isothermal behavior for small angles $\alpha_{1}$, we explore other choices for the location of the critical point, which consequently change $\alpha_{1}$. The angle of the mapping tangent to the chiral phase transition line increases with increasing $\mu_{B}$ as the parabola approaches its terminus at $T=0$. Results for the isotherms for several different options for the location of the critical point at $\mu_{B,c}=550,\,750,\,900$ MeV are shown in Fig. 8. As discussed in Sec. III, when choosing the critical point along the chiral phase transition line we obtain a constraint on the angle $\alpha_{1}$. These options for placing the critical point lead to mapping angles of $\alpha_{1}=9.2^{\rm{o}},\,16.6^{\rm{o}},\,28.3^{\rm{o}}$, respectively. We show that the spinodal features seen in the isothermal trajectories are still distorted for an angle of $9.2^{\rm{o}}$, at least at this temperature relatively close to the critical point. On the other hand, for the subsequent choices corresponding to $\alpha_{1}=16.6^{\rm{o}},28.3^{\rm{o}}$, we capture the typical isothermal behavior. In particular, we can now identify the metastable phases beyond the coexistence point. We see that both these choices for the placement of the critical point give rise to super-heated and super-cooled phases, where the blue and yellow solutions track past the coexistence point and begin to describe a metastable state. In between the spinodal points lies the instability and the unique unstable solution shown in orange. From this, we demonstrate the strong influence that the angle $\alpha_{1}$ has on the mapping of the spinodal points onto the QCD phase diagram. This is, however, not the only parameter affecting the spinodals in the QCD phase diagram. We show the effect of additional parameters in Appendix A.

Beyond a single isothermal curve for the pressure as a function of baryon chemical potential, we determine the spinodal curves in the QCD phase diagram. From Eq. (11) we obtain the $T$- and $\mu_{B}$-dependence of the spinodals in order to determine the extent of the spinodal region in the phase diagram for a given choice of mapping parameters. We choose the same three locations for the critical point as those shown in Fig. 8, which are $\mu_{B,c}=550,750,900$ MeV.

The results for the spinodal lines in the QCD phase diagram are shown in Fig. 9 for $\mu_{B,c}=550$ MeV (top), $\mu_{B,c}=750$ MeV (middle), $\mu_{B,c}=900$ MeV (bottom), all with the choice of scaling parameters $w=1,~{}\rho=2$. In each of these plots, our parametrization of the chiral phase transition line is shown along which sits the critical point, given as the filled circle, for each of the three choices. We see that, for smaller values of $\mu_{B,c}$ and consequently of the angle $\alpha_{1}$, the metastable region covers more of the phase diagram. We also show the coexistence line as mapped onto the phase diagram given our choice of the mapping shown in Eq. (9). In this case, the transition line ($h=0$) from the Ising model is a straight line tangent to the chiral phase transition line ($T_{c}(\mu_{B}/T_{c})$) due to the fact that we employ this linear map. The coexistence line should rather curve down toward the $\mu_{B}$-axis, and when $T=0$ it should meet the axis at $90^{\rm{o}}$, as discussed in Sec. III.2. For the largest angle choice, we see that the linear approximation follows closely the chiral phase transition line, as the tangent line becomes a better approximation as the angle increases. As we are focused mainly on the spinodal implementation in the QCD phase diagram, we leave the study of an alternative mapping procedure to future work.

For the first choice of $\mu_{B,c}$, the high-temperature spinodal never reaches the $T=0$ axis, but rather moves away from it at large chemical potential. At this point, we compare our result for the spinodal lines in the phase diagram to those obtained in the holographic model in Refs. Grefa et al. (2021); Hippert et al. (2023). The holographic model provides an excellent description of lattice QCD thermodynamics, it naturally incorporates the ideal fluidity of the Quark-Gluon Plasma and it predicts a critical point on the QCD phase diagram. Starting from the critical point, at chemical potentials $\mu_{B}>\mu_{B,c}$, the first-order phase transition (coexistence) line departs, surrounded by the two spinodals. Entropy density, energy density and net-baryon density become multi-valued solutions of the temperature at fixed chemical potential in this region, and the high-temperature spinodal line runs almost parallel to the $T=0$ axis at large $\mu_{B}$, similarly to what happens in the middle panel of Fig. 9 in our case. Additionally, critical exponents were extracted from this holographic approach and were found to be that of mean field DeWolfe et al. (2011), thus further linking these two sets of results. From this we also see that, if the parameters of the mapping can be constrained by experimental efforts in the search for the QCD critical point, we will be able to determine the size of the spinodal region in the QCD phase diagram.

Since QCD is symmetric for the exchange of matter with anti-matter, the charge conjugation symmetry must be obeyed in our equation of state across the $\mu_{B}=0$ line. Furthermore, as can be seen from the isothermal trajectories, without symmetrization the density would be negative for the solution shown in Fig. 8 in yellow which has a negative slope. The symmetrization proceeds by utilizing the information at $-\mu_{B}$ such that:

$$\begin{split}P_{\text{symm}}(T,\mu_{B})&=\frac{1}{2}(P(T,\mu_{B})+P(T,-\mu_{B})),\\ n_{B\text{,symm}}(T,\mu_{B})&=\frac{1}{2}(n_{B}(T,\mu_{B})-n_{B}(T,-\mu_{B})).\end{split}$$

(15)

In order to show the effect of the symmetrization and calculate further characteristic signatures of the first order phase transition, we will proceed with one of the choices of the placement of the critical point, $\mu_{B,c}=550$ MeV. Given the results for the spinodal curves in the phase diagram shown in Fig. 9, we choose an isotherm that is sufficiently close to the critical point such that both spinodal curves are captured by that isotherm. With this in mind, we choose an isotherm of $T=0.98T_{c}=120$ MeV.

In Fig. 10 we show the isothermal pressure curve for the choice motivated here of $\mu_{B,c}=550$ MeV, $T=120$ MeV, before and after the symmetrization process. It is evident that after symmetrization the sign of the slope changes for the left and middle solutions shown in yellow and orange. This is an important development as a result of the symmetrization because the slope of the curve $P(\mu_{B})$ corresponds directly to the baryon density $n_{B}=(\partial P/\partial\mu_{B})_{T}$. The baryon density, thus, becomes positive, as we would expect for our system in a matter-dominant regime with positive baryon chemical potential.

After the symmetrization is complete, we can study the behavior of the baryon density itself in the first order phase transition region for $\mu_{B}>\mu_{B,c}$. Figure 11 (left panel) shows the symmetrized baryon density from the Ising model as a function of the chemical potential, as well as the behavior of the pressure as a function of the baryon density after symmetrization (right panel). Important features are present in the baryon density, which were not previously available in the critical equation of state mapped to QCD with the 3D Ising model Parotto et al. (2020); Karthein et al. (2021). Since the spinodal features lie in the complex plane for the 3D Ising model as discussed in Sec. II, the density would simply jump from the coexistence point given by the open green marker on the blue curve to the next coexistence point on the yellow curve during the cooling process. On the other hand, in our mean field Ising approach, we can access the full curve by proceeding through the spinodal points as well. In this case, both the coexistence points along the first order phase transition line and the spinodal points reaching out from them are now clearly present for each of the solutions of the critical baryon density. Here, we show the three different solutions given by the yellow, orange, and blue curves with the coexistence and spinodal points highlighted in open and filled circles, respectively. Similarly for the behavior of the pressure as a function of the baryon density, we achieve a result that shows the progression through the mixed phase. The orange curve shown here depicts a mechanical instability where $\partial P/\partial n_{B}<0$, further affirming that this solution represents the unstable phase.