Comparison of Equations of State for Neutron Stars with First-Order Phase Transitions: A Qualitative Study

Anshuman Verma Indian Institue of Science Education and Research Bhopal Asim Kumar Saha Indian Institue of Science Education and Research Bhopal Ritam Mallick Indian Institue of Science Education and Research Bhopal

Abstract

The equation of state is fundamental in describing matter under the extreme conditions characteristic of neutron stars and is central to advancing our understanding of dense matter physics. A critical challenge, however, lies in accurately modelling first-order phase transitions while ensuring thermodynamic consistency and aligning with astrophysical observations. This study explores two frameworks for constructing EoSs with first-order phase transitions: the polytropic interpolation method and the randomized speed-of-sound interpolation approach. It is found that the mass-radius relation and pressure vs. energy density relation are blind towards the thermodynamic consistency check. The polytropic interpolation method can exhibit discontinuities in the chemical potential for first-order phase transition, raising concerns regarding potential causality violations and thermodynamic inconsistencies. In contrast, the speed of sound interpolation approach ensures continuity in the chemical potential, offering a more thermodynamically consistent and reliable framework. Moreover, the sound speed method effectively captures the softer segment of the mass-radius spectrum, a capability not achieved by the consistent piecewise-polytropic approach due to its monotonic stiffness constraints. The speed of sound definition involving number density and chemical potential reveals the thermodynamic inconsistency, making it a more consistent and robust definition. These findings underscore the importance of thermodynamic consistency in EoS construction and highlight the advantages of the randomized speed-of-sound method for modelling phase transitions in dense matter.

1 Introduction

Equation of state (EoS) is fundamental for understanding matter at extreme densities, particularly within neutron stars (Lattimer and Prakash, 2004; Lattimer, 2014; McLerran and Reddy, 2019; Bombaci, 2022; Gao et al., 2025). The EoS is the relation between the thermodynamic variables of the system, which describe matter properties. The EoS at extreme densities is still unknown, and both ab initio calculation and earth-based experiments are still being awaited. One, therefore, relies on model calculations and finally matches their predictions with experiments (Schoof et al., 2015; Piekarewicz, 2017; Mandal et al., 2020; Karthein, 2024). The only laboratory to test the models is the observation from neutron stars (NSs). In the last few decades, there has been tremendous work in this regard, both theoretically and observationally (Zdunik et al., 2008; Abbott et al., 2020; Komoltsev and Kurkela, 2022; Kumar et al., 2024). One of the most efficient ways to constrain the EoS is constructing an ensemble of EoS by randomizing the speed of sound and using thermodynamic bounds (Chu and Chen, 2014; Chatterjee et al., 2024; Thakur et al., 2024). Once the ensemble is constructed, the EoS ensemble is constrained using the astrophysical observational bounds (Haensel et al., 2002; Lattimer and Prakash, 2006; Miller, 2021). Initially, this was employed for smooth EoS (EoS does not have any first-order phase transition (FOPT)) and has been very successful in constraining the EoS to a large extent. Later, an ensemble of FOPT was also constructed using a piecewise-polytrope (PP) approach and then astrophysical observational bounds were employed to constrain the EoS (Vijaykumar et al., 2021; Mallick, 2022; Abbott et al., 2022; Gorda et al., 2023).

EoS with FOPT signifies a rapid, discontinuous change in the internal structure of dense matter, potentially leading to quark matter cores (Bombaci and Datta, 2000; Demorest et al., 2010). These transitions significantly impact macroscopic properties like neutron stars’ mass, radius, and frequencies (Klähn et al., 2007; Wang et al., 2012; Yang et al., 2017). The EoS also influences the NS stability and observable properties like star-quakes and magnetic field dynamics (Bhattacharyya et al., 2017). A PT may also result in abrupt changes in radius or mass, which can produce detectable signals in multi-messenger astrophysics (Alford et al., 2007; Zheng and Chen, 2016; Miller, 2021). Identifying signatures of FOPT is critical for advancing our understanding of compact objects and quantum chromodynamics (QCD) (Bissell et al., 2014; Lü et al., 2015; Abbott et al., 2020). However, astrophysical observation to date has not been able to clearly distinguish whether quark matter exists at the cores of neutron stars, let alone whether there is a smooth or FOPT (Xu et al., 2009; Acernese et al., 2014; Aasi et al., 2015; Bauswein and Stergioulas, 2019).

Model construction of an ensemble of EoS is a novel idea; however, one should be careful and check its thermodynamic consistency. As these are not microscopic models, one should check them rigorously so that they do not violate thermodynamical consistency. More specifically, FOPT does have a discontinuity in energy and number density; however, pressure and chemical potential must be continuous. There lies a pitfall as the construction of NSs mainly involves energy density, number density and pressure. The chemical potential is not directly related to the Tolmann-Volkoff-Oppendeimer (TOV) (Tolman, 1939). The PP approach (Gorda et al., 2023) does not effectively spell out the chemical potential continuity and the mass-radius curves are blind to it. For a correctly modelled FOPT EoS, the pressure and chemical potential must be continuous at the transition point. This ensures mechanical and thermodynamic equilibrium, ensuring Euler theorem and the Gibbs phase rule (Radhakrishnan, 1990; Janakiraman, 2018). The first-order transition is characterized by a discontinuous jump in density (or volume per particle) as the system shifts from a lower-density phase (e.g., hadronic matter) to a higher-density phase (e.g., quark matter) or vice-versa (Bauswein et al., 2019; Zhou et al., 2024).

The speed of sound interpolation (SoSI) is not usually used to generate FOPT and is restricted only to generate smooth EoS; it can also be employed to generate FOPT. The main aim of this paper is to generate thermodynamic consistent EoS using the two formalisms: the SoSI and PP. More importantly, the condition that needs to be fulfilled to construct a correct thermodynamic consistent EoS using the PP approach is derived. The structure of the paper is organized as follows: Section 2.3 addresses the construction of the EoS ensemble, ensuring thermodynamic consistency. Section 3 presents the results, and Section 4 summarises our findings and concludes the study.

2 Construction Methods

The construction methods addressed in this section are:

•

Piecewise-polytropic (PP) method (Gorda et al. (2023)).
•

Piecewise-polytropic consistent (PP-consistent) method.
•

Speed of sound interpolation (SoSI) method (Annala et al. (2020)).

The PP method suffers a chemical potential jump across the transition, which leads to inconsistencies, as discussed in section 3. We propose the PP-consistent method to address the inconsistency, which ensures a continuous chemical potential across the jump.

2.1 Piecewise Polytropic (PP) Construction

An ensemble of EoS using a PP approach is constructed up to a baryon density of 10 $n_{0}$ (where $n_{0}=0.16\,\text{fm}^{-3}$ , nuclear saturation density).The construction of the EoS involves three key steps:

1.

For densities $n\leq 0.57n_{0}$ , we employ the BPS crust EoS from Baym et al. (1971). From $n=0.57n_{0}$ up to $n_{\text{CET}}=1.1n_{0}$ , we interpolate between the soft and stiff EoSs within the chiral CET band of Hebeler et al. (2013).
2.

For densities $n_{\text{CET}}<n\leq n_{\text{PT}}$ , where $n_{PT}\in[1.1,10]n_{0}$ ; a single polytropic form is used:

$p(n)=p_{\text{CET}}(n_{\text{CET}})\left(\frac{n}{n_{\text{CET}}}\right)^{\Gamma_{1}}$
3.

For $n_{\text{PT}}+\Delta n<n\leq 10n_{0}$ ; where $\Delta n/n_{PT}\in[0,9.09]$ , a second polytropic EoS is applied:

$p(n)=p(n_{\text{PT}})\left(\frac{n}{n_{\text{PT}}+\Delta n}\right)^{\Gamma_{2}}$

The phase transition (PT) occurs between the two polytropic segments. It is characterized by two key parameters: the onset of the phase transition $n_{PT}$ and the density jump $\Delta n$ .

2.2 Speed-of-Sound Interpolation (SoSI) Construction

Refer to caption — Figure 1: Diagram illustrating the speed of sound randomization in the $c_{s}^{2}-\mu$ space. A first-order phase transition happens at $\mu_{PT}$ , and the pQCD region begins at $\mu_{end}=2.6$ GeV.

The speed of sound ( $c_{s}^{2}$ ) is parametrised as a function of chemical potential ( $\mu$ ) (Annala et al., 2020), which allows us to write the number density as

n(\mu)=n_{\scriptscriptstyle CET}\exp\left[\int_{\mu_{CET}}^{\mu}\frac{d\mu^{\prime}}{\mu^{\prime}c_{s}^{2}(\mu^{\prime})}\right]

(1)

where $n_{CET}$ is fixed by matching the lower-density CET EoS to the initials of the interpolated EoS, and $\mu_{CET}$ is the chemical potential corresponding to $n_{CET}$ . Pressure can then readily be calculated as

p(\mu)=p_{\scriptscriptstyle CET}\!+n_{\scriptscriptstyle CET}\int_{\mu_{CET}}^{\mu}d\mu^{\prime}\exp\left[\int_{\mu_{CET}}^{\mu^{\prime}}\frac{d\mu^{\prime\prime}}{\mu^{\prime\prime}c_{s}^{2}(\mu^{\prime\prime})}\right]

(2)

where $p_{CET}$ corresponds to the pressure value at the matching point of the CET EoS. Our study here includes 5 randomised segments of ( $c_{s,i}^{2},\mu_{i}$ ), where $\mu_{i}\in[\mu_{\scriptscriptstyle CET},2.6$ $GeV]$ . The upper limit of $\mu$ is chosen by Fraga et al. (Fraga et al., 2014), such that the uncertainty in the pQCD regime is roughly the same as the uncertainty at CET. A piecewise linear function is employed to connect the points $\left\{\mu_{i},c_{s,i}^{2}\right\}$ as:

c_{s}^{2}(\mu)=\frac{(\mu_{i+1}-\mu)c_{s,i}^{2}+(\mu-\mu_{i})c_{s,i+1}^{2}}{\mu_{i+1}-\mu_{i}}

(3)

where $c_{s,i}^{2}\in[0,1]$ . This function aids us to carry out the integrals of (1) and (2).

SoSI construction method inherently permits FOPT; however, they are statistically suppressed. To address this, we modify the construction by dividing it into two segments: the first and second branches (see Fig 1). The first branch, which mimics the hadronic branch, is extended up to a chemical potential of $\mu_{\scriptscriptstyle PT}$ ; where $\mu_{\scriptscriptstyle PT}$ is chosen randomly in the range between $\mu_{\scriptscriptstyle CET}$ and onset of pQCD potential. In order to induce a FOPT between the two branches, the first branch’s endpoint and the second branch’s starting point are particularly significant. The integrals for the thermodynamic quantities at the endpoint of the first branch are expressed as:

	$\displaystyle{n(\mu_{\scriptscriptstyle{PT}})}$	$\displaystyle={n_{\scriptscriptstyle CET}\exp\left[\int_{\mu_{\scriptscriptstyle CET}}^{\mu_{\scriptscriptstyle PT}}\frac{d\mu^{\prime}}{\mu^{\prime}c_{s}^{2}(\mu^{\prime})}\right]}$		(4)
	$\displaystyle{p(\mu_{\scriptscriptstyle{PT}})}$	$\displaystyle={p_{\scriptscriptstyle CET}\!+n_{\scriptscriptstyle CET}\int_{\mu_{\scriptscriptstyle CET}}^{\mu_{\scriptscriptstyle PT}}d\mu^{\prime}\exp\left[\int_{\mu_{\scriptscriptstyle CET}}^{\mu^{\prime}_{\scriptscriptstyle PT}}\frac{d\mu^{\prime\prime}}{\mu^{\prime\prime}c_{s}^{2}(\mu^{\prime\prime})}\right]}$		(5)

To exhibit FOPT, the number density at the starting point of the second branch is induced with a jump of $\Delta n$ , which allows us to write the integrals of (1) and (2) for the thermodynamic quantities of the second branch as:

	$\displaystyle{n(\mu)}$	$\displaystyle={(n_{\scriptscriptstyle PT}+\Delta n)\exp\left[\int_{\mu_{PT}}^{\mu}\frac{d\mu^{\prime}}{\mu^{\prime}c_{s}^{2}(\mu^{\prime})}\right]}$		(6)
	$\displaystyle{p(\mu)}$	$\displaystyle={p_{\scriptscriptstyle PT}\!+(n_{\scriptscriptstyle PT}+\Delta n)\int_{\mu_{PT}}^{\mu}d\mu^{\prime}\exp\left[\int_{\mu_{PT}}^{\mu^{\prime}}\frac{d\mu^{\prime\prime}}{\mu^{\prime\prime}c_{s}^{2}(\mu^{\prime\prime})}\right]}$		(7)

where $p_{\scriptscriptstyle PT}=p(\mu_{\scriptscriptstyle PT})$ and $n_{\scriptscriptstyle PT}=n(\mu_{\scriptscriptstyle PT})$ . The value of $\Delta n$ is chosen to span the space similar to that of the polytropic construction. The parameters defining FOPT for the SoSI construction are $\mu_{PT}$ and $\Delta n$ , while for the PP construction, they are $n_{PT}$ and $\Delta n$ . The stiffness of the two branches in the PP construction is determined by the parameters $\Gamma_{1}$ and $\Gamma_{2}$ , while the SoSI construction is governed by the speed of sound segments, which are randomized. Adhering to the astrophysical constraints of Demorest et al. (2010), we ensure that every EoS is capable of generating M ${}_{\scriptscriptstyle TOV}>2.0M_{\odot}$ as well as respecting the tidal deformability limit of $\tilde{\Lambda}<720$ .

This study analyses the EoS derived from polytrope-based construction and the SoSI method. In order to facilitate a valid comparison between the two methodologies, we ensure that the transition onset density for the SoSI method is maintained within the interval $[1.1,10]n_{0}$ , which aligns with the range employed in the polytropic construction method. This alignment is achieved by imposing a constraint on the value of the integral in equation 4, keeping it below the threshold of $10n_{0}$ .

2.3 Piecewise Polytrope Consistent (PP-consistent) method

To confirm the continuity of the chemical potential within the framework of the polytropic construction method, as elaborated in 2.1, we define the endpoint thermodynamic variables of the hadronic branch as

$\displaystyle p^{\scriptscriptstyle{h}}_{\scriptscriptstyle{PT}}(n^{\scriptscriptstyle{h}}_{\scriptscriptstyle{PT}})$	$\displaystyle=p_{\scriptscriptstyle{CET}}\left(\frac{n^{\scriptscriptstyle{h}}_{\scriptscriptstyle{PT}}}{n_{\scriptscriptstyle{CET}}}\right)^{\Gamma_{1}},$
$\displaystyle\epsilon^{\scriptscriptstyle{h}}_{\scriptscriptstyle{PT}}(n^{\scriptscriptstyle{h}}_{\scriptscriptstyle{PT}})$	$\displaystyle=m_{\scriptscriptstyle{b}}n_{\scriptscriptstyle{PT}}^{\scriptscriptstyle{h}}+\frac{p^{\scriptscriptstyle{h}}_{\scriptscriptstyle{PT}}}{\Gamma_{1}-1},$
$\displaystyle\mu^{\scriptscriptstyle{h}}_{\scriptscriptstyle{PT}}$	$\displaystyle=m_{\scriptscriptstyle{b}}+\frac{p^{\scriptscriptstyle{h}}_{\scriptscriptstyle{PT}}}{n^{\scriptscriptstyle{h}}_{\scriptscriptstyle{PT}}}\left(\frac{\Gamma_{1}}{\Gamma_{1}-1}\right)$	(8)

After undergoing FOPT, at the onset of the quark branch, the thermodynamic parameters are expressed as:

	$\displaystyle\epsilon^{\scriptscriptstyle{q}}_{\scriptscriptstyle{PT}}(n^{\scriptscriptstyle{q}}_{\scriptscriptstyle{PT}})$	$\displaystyle=m_{\scriptscriptstyle{q}}n_{\scriptscriptstyle{PT}}^{\scriptscriptstyle{q}}+\frac{p^{\scriptscriptstyle{q}}_{\scriptscriptstyle{PT}}}{\Gamma_{2}-1},$
	$\displaystyle\mu^{\scriptscriptstyle{q}}_{\scriptscriptstyle{PT}}$	$\displaystyle=m_{\scriptscriptstyle{q}}+\frac{p^{\scriptscriptstyle{q}}_{\scriptscriptstyle{PT}}}{n^{\scriptscriptstyle{q}}_{\scriptscriptstyle{PT}}}\left(\frac{\Gamma_{2}}{\Gamma_{2}-1}\right)$		(9)

where the pressure at this junction is maintained as $p_{\scriptscriptstyle{PT}}^{\scriptscriptstyle{q}}=p_{\scriptscriptstyle{PT}}^{\scriptscriptstyle{h}}$ . The number density adapts according to $n_{\scriptscriptstyle{PT}}^{\scriptscriptstyle{q}}=n_{\scriptscriptstyle{PT}}^{\scriptscriptstyle{h}}+\Delta n$ . By equating the chemical potentials given by equations 8 and 9, the thermodynamic consistency condition reduces to

\frac{(m_{q}-m_{b})n_{\scriptscriptstyle{PT}}^{\scriptscriptstyle{q}}}{p_{\scriptscriptstyle{PT}}}\frac{(\Gamma_{1}-1)}{\Gamma_{1}}+\frac{\Gamma_{1}-\Gamma_{2}}{\Gamma_{1}(\Gamma_{2}-1)}=\frac{\Delta n}{n_{\scriptscriptstyle{PT}}^{\scriptscriptstyle{h}}}

(10)

This relationship is essential to ensure that the chemical potential remains continuous during the transitional phase. By imposing this condition, one can operate within the framework of the piecewise polytropic formalism, guaranteeing thermodynamic consistency across the transition.

3 Results and Discussion

The EoS constructed in two different ways are considerably different. Figure 2 presents a comparative analysis of the EoS and mass-radius diagram obtained from two different methods. In the left (top) panel, the plot illustrates the relationship between pressure ( $P$ ) and central energy density ( $\epsilon$ ) on a logarithmic scale for EoS models created using the PP-consistent method (green contour) and the SoSI method (blue contour). Both approaches comply with causality constraints. The PP-consistent method employs predefined functional forms, resulting in a more constrained FOPT EoS. This leads to a narrower range of possible behaviours, with controlled variations at intermediate densities and convergence toward the pQCD limit at high densities. Conversely, the SoSI method (blue contour) offers greater flexibility, allowing a broader spectrum of phase transition scenarios. The wider spread of EoS curves at lower energy densities (up to approximately 1000 MeV/fm³) highlights the uncertainties associated with nuclear matter properties. However, this variability diminishes at higher densities as the EoS approaches the pQCD limit. This broader distribution at low densities reflects the sensitivity of the speed of sound to phase transition characteristics, capturing a wider range of potential EoS behaviours. The left (bottom) panel displays the corresponding mass-radius (M-R) curves obtained by solving the TOV equations. Both methods satisfy key astrophysical constraints, including a maximum neutron star mass exceeding $2M_{\odot}$ and a tidal deformability limit of $\widetilde{\Lambda}<720$ . The PP-consistent method (green contour) produces more tightly clustered M-R curves, indicating a more constrained range of neutron star properties. In contrast, the SoSI method (blue contour) exhibits a broader spread, reflecting more significant variability in phase-transition behaviour (Annala et al., 2020).

The right panel of Fig. 2 presents the PDF plot of the M-R curves for the SoSI construction method (top) and the PP-consistent method (bottom). While the SoSI construction covers a larger area in the M-R space, it is mainly concentrated in a more confined region. A similar pattern is observed for the PP-construction method, with the concentration increases in a confined region having a substantial overlap in coverage with the SoSI ensemble.

For early phase transitions, generally happening at a density less than 2 $n_{0}$ , a substantial increase in stiffness beyond the transition is necessary to meet the constraints set by observations (simultaneous satisfaction of the $2M_{\odot}$ bound and also maintaining the limit of tidal deformability $\Lambda<720$ ). This requirement imposes stringent limits on the EoS and the corresponding M-R diagram. As both the PP construction methods are restricted compared to the SoSI method, they could not generate consistent EoS which satisfies astrophysical bounds. It is most interesting to note that although PP (inconsistent) was able to generate at least a few EoS satisfying the astrophysical bound, the PP-consistent method was not able to generate any (detailed in the Appendix section B). For the PP-consistent approach, all the EoS undergo transitions at densities greater than $2n_{0}$ . The limitation of the PP-consistent method arose because of the consistency condition (eqn 10). For the PP method, $\Gamma_{2}$ is considered a free parameter, offering more flexibility in stiffness adjustment. However, within the PP-consistent methodology, $\Gamma_{2}$ is not independently variable, as its value is directly influenced by the parameters $\Gamma_{1}$ , $\Delta n$ , and $n^{h}_{\scriptscriptstyle{PT}}$ as specified by equation 10. On the other hand, the SoSI method does not have such restriction and allows for significant variation in stiffness even after the early PT to satisfy astrophysical bounds.

A particularly noteworthy observation is that the PP-based EoS cannot cover the softer region of the M-R spectrum. This arises from the functional nature of polytropes, which exhibit monotonically increasing stiffness. Consequently, a single polytrope cannot simultaneously produce the required soft segment to describe low-mass stars; thereafter, the stiff segment is needed to satisfy astrophysical constraints (Ferrari et al., 2010; Kurkela et al., 2014).

As already stated, it is difficult to infer the consistency of the FOPT EoS from $P$ vs. $\epsilon$ curve or the M-R sequence curve. One of the most essential thermodynamic criteria, even for FOPT, is that the pressure and chemical potential should be smooth without any discontinuity. The left panel of Figure 3 illustrates the relationship between $\mu$ and $P$ for three distinct EoS models. The EoS constructed using the PP method (depicted in red) exhibits a discontinuity in the $\mu$ - $P$ curve, which stems from the specific formulation of this approach. In contrast, the SoSI-based EoS (blue) and PP-consistent method (green) produce a smooth, continuous curve.

A critical parameter for identifying phase transitions in dense matter is $c_{s}^{2}$ . This can be defined in two ways: in terms of pressure and energy density ( $c_{s}^{2}(P,\epsilon)$ ) or using number density and chemical potential ( $c_{s}^{2}(n,\mu)$ ) (Abgaryan et al., 2018).

c_{s}^{2}=\left(\frac{\partial(log_{e}\mu)}{\partial(log_{e}n)}\right)_{s}=\left(\frac{\partial P}{\partial\epsilon}\right)_{s}

(11)

Both definitions should maintain the causality for an EoS to be physically consistent. In the case of the PP-based model, the discontinuity in the $\mu$ - $P$ relationship manifests itself as a discrepancy between these two relations. As confirmed by fig 3 right, the SoSI method EoS maintains causality and generates the same curve with two different definitions. Whereas the PP-consistent EoS maintains causality and generates the same curve for two different definitions, the same is not true for inconsistent PP EoS. The speed-of-sound curve differs, and the definition involving chemical potential and number density violates causality. The figure also hints at one of the pitfalls of the consistency check of the EoS, as the speed of sound definition with pressure and energy density does not violate causality and seems to be thermodynamically consistent.

The discrepancy in the M-R diagram (not having stars with consistent PP EoS for transition density less than twice the saturation density) results from this discrepancy. The PP method is a more restricted method, only using a single polytrope before and after the PT. If more than one polytrope is used, one can have more freedom in obtaining more EoS, which satisfies both thermodynamic consistency and astrophysical bound. This reinforces the conclusion that the SoSI method offers a more comprehensive and reliable framework for modelling first-order phase transitions in dense neutron star matter.

4 Summary and Conclusion

The study examines the construction of an ensemble of EoS involving FOPT. Two separate methods are used for the construction: the SoSI and the PP method. One of the paper’s key findings is that the $P$ vs. $\epsilon$ diagram or the $M-R$ curve are blind towards the consistency check of the EoS. A consistent EoS, even with FOPT, should not have any pressure or chemical potential discontinuity. Even if they follow the low and high-density limits and satisfy the astrophysical bounds of neutron stars and GW170817 tidal deformability, the EoS can be thermodynamically inconsistent.

The study identifies the consistency condition for the PP ensemble of EoS. However, even with consistent EoS, the PP-consistent EoS and SoSI method produces very different EoS ensembles. The SoSI method has more flexibility and can generate all possible EoS models, whereas even PP-consistent EoS is largely restricted. The PP-consistent method has a single polytrope on either side of the PT discontinuity and, therefore, becomes very restricted in adhering to the consistency of the astrophysical constraints.

One of the essential aspects of the study is that early phase transitions at densities below $2n_{0}$ are associated with larger neutron star radii but necessitate significant EoS stiffening to support $2M_{\odot}$ stars. Many early transition models fail to satisfy tidal deformability constraints, particularly the limit $\widetilde{\Lambda}<720$ . This limitation is highlighted in the PP method, which restricts their ability to achieve the required balance of softness and stiffness to meet astrophysical criteria.

The inconsistency of the PP model is also revealed in the definition of speed of sound. The speed of sound can be expressed as a function of pressure and energy density or as a function of number density and chemical potential. For inconsistent PP EoS, the definition of $c_{s}^{2}$ with energy density and pressure maintains causality, whereas it violates causality for the definition involving number density and the chemical potential. For consistent EoS, both definitions give the same $c_{s}^{2}$ maintaining causality. Therefore, one of the key points of the study is that the definition of $c_{s}^{2}$ involving $n$ and $\mu$ is more robust and should be rigorously used.

The main research contribution is the qualitative comparison of two widely used methods for constructing EoS with phase transitions. By demonstrating the critical role of chemical potential continuity and thermodynamic consistency, this work identifies the SoSI method as a more physically reliable approach for modelling dense matter in neutron stars. Future directions include expanding the analysis to include multi-messenger data to improve the robustness of these models further and provide deeper insights into the nature of dense matter and phase transitions.

Acknowledgements

The authors thank the Indian Institute of Science Education and Research Bhopal for providing all the research and infrastructure facilities. AV acknowledges the Prime Minister’s Research Fellowship (PMRF), Ministry of Education, Government of India, for a graduate fellowship, and RM acknowledges the Science and Engineering Research Board (SERB), Govt. of India, for financial support in the form of a Core Research Grant (CRG/2022/000663).

References

Lattimer and Prakash (2004) J. M. Lattimer and M. Prakash, Science 304, 536 (2004), arXiv:astro-ph/0405262 [astro-ph] .
Lattimer (2014) J. M. Lattimer, General Relativity and Gravitation 46, 1713 (2014).
McLerran and Reddy (2019) L. McLerran and S. Reddy, Phys. Rev. Lett. 122, 122701 (2019), arXiv:1811.12503 [nucl-th] .
Bombaci (2022) I. Bombaci, in Astrophysics and Space Science Library, Astrophysics and Space Science Library, Vol. 465, edited by S. Bhattacharyya, A. Papitto, and D. Bhattacharya (2022) pp. 281–317.
Gao et al. (2025) Q. Gao, N. Dai, Y. Gong, et al., ApJ 978, 8 (2025), arXiv:2306.06645 [gr-qc] .
Schoof et al. (2015) T. Schoof, S. Groth, J. Vorberger, and M. Bonitz, Phys. Rev. Lett. 115, 130402 (2015).
Piekarewicz (2017) J. Piekarewicz, in Handbook of Supernovae, edited by A. W. Alsabti and P. Murdin (2017) p. 1075.
Mandal et al. (2020) R. Mandal, P. J. Bhuyan, P. Chaudhuri, et al., Nature Communications 11, 2581 (2020).
Karthein (2024) J. Karthein, in European Physical Journal Web of Conferences, European Physical Journal Web of Conferences, Vol. 296 (2024) p. 01018.
Zdunik et al. (2008) J. L. Zdunik, M. Bejger, P. Haensel, and E. Gourgoulhon, A&A 479, 515 (2008), arXiv:0707.3691 [astro-ph] .
Abbott et al. (2020) R. Abbott, T. D. Abbott, S. Abraham, et al., ApJ 896, L44 (2020), arXiv:2006.12611 [astro-ph.HE] .
Komoltsev and Kurkela (2022) O. Komoltsev and A. Kurkela, Phys. Rev. Lett. 128, 202701 (2022).
Kumar et al. (2024) R. Kumar, V. Dexheimer, J. Jahan, et al., Living Reviews in Relativity 27, 3 (2024), arXiv:2303.17021 [nucl-th] .
Chu and Chen (2014) P.-C. Chu and L.-W. Chen, ApJ 780, 135 (2014), arXiv:1212.1388 [astro-ph.SR] .
Chatterjee et al. (2024) S. Chatterjee, H. Sudhakaran, and R. Mallick, The European Physical Journal C 84, 1291 (2024).
Thakur et al. (2024) P. Thakur, S. Chatterjee, K. K. Nath, and R. Mallick, Phys. Rev. D 110, 103045 (2024).
Haensel et al. (2002) P. Haensel, J. L. Zdunik, and F. Douchin, A&A 385, 301 (2002), arXiv:astro-ph/0201434 [astro-ph] .
Lattimer and Prakash (2006) J. M. Lattimer and M. Prakash, Nucl. Phys. A 777, 479 (2006).
Miller (2021) M. C. Miller, Timing Neutron Stars: Pulsations, Oscillations and Explosions , 1 (2021).
Vijaykumar et al. (2021) A. Vijaykumar, S. J. Kapadia, and P. Ajith, Phys. Rev. Lett. 126, 141104 (2021), arXiv:2003.12832 [gr-qc] .
Mallick (2022) R. Mallick, in European Physical Journal Web of Conferences, European Physical Journal Web of Conferences, Vol. 274 (2022) p. 07002, arXiv:2212.00462 [astro-ph.HE] .
Abbott et al. (2022) R. Abbott, T. D. Abbott, F. Acernese, et al., ApJ 932, 133 (2022), arXiv:2112.10990 [gr-qc] .
Gorda et al. (2023) T. Gorda, K. Hebeler, A. Kurkela, et al., The Astrophysical Journal 955, 100 (2023).
Bombaci and Datta (2000) I. Bombaci and B. Datta, ApJ 530, L69 (2000), arXiv:astro-ph/0001478 [astro-ph] .
Demorest et al. (2010) P. B. Demorest, T. Pennucci, S. M. Ransom, et al., Nature 467, 1081 (2010).
Klähn et al. (2007) T. Klähn, D. Blaschke, F. Sandin, et al., Physics Letters B 654, 170 (2007), arXiv:nucl-th/0609067 [nucl-th] .
Wang et al. (2012) J. Wang, C. M. Zhang, and H. K. Chang, A&A 540, A100 (2012), arXiv:1101.3042 [astro-ph.HE] .
Yang et al. (2017) Y.-Y. Yang, C.-M. Zhang, D. Li, et al., ApJ 835, 185 (2017).
Bhattacharyya et al. (2017) S. Bhattacharyya, I. Bombaci, D. Logoteta, and A. V. Thampan, ApJ 848, 65 (2017), arXiv:1709.02415 [astro-ph.HE] .
Alford et al. (2007) M. Alford, D. Blaschke, A. Drago, et al., Nature 445, 7 (2007), arXiv:astro-ph/0606524 [astro-ph] .
Zheng and Chen (2016) H. Zheng and L.-W. Chen, ApJ 831, 127 (2016), arXiv:1603.07518 [nucl-th] .
Bissell et al. (2014) M. L. Bissell, J. Papuga, H. Naïdja, et al., Phys. Rev. Lett. 113, 052502 (2014).
Lü et al. (2015) H.-J. Lü, B. Zhang, W.-H. Lei, et al., ApJ 805, 89 (2015), arXiv:1501.02589 [astro-ph.HE] .
Xu et al. (2009) J. Xu, L.-W. Chen, B.-A. Li, and H.-R. Ma, ApJ 697, 1549 (2009), arXiv:0901.2309 [astro-ph.SR] .
Acernese et al. (2014) F. Acernese, M. Agathos, K. Agatsuma, et al., Classical and Quantum Gravity 32, 024001 (2014).
Aasi et al. (2015) J. Aasi, B. Abbott, R. Abbott, et al., Classical and quantum gravity 32, 074001 (2015).
Bauswein and Stergioulas (2019) A. Bauswein and N. Stergioulas, Journal of Physics G: Nuclear and Particle Physics 46, 113002 (2019).
Tolman (1939) R. C. Tolman, Physical Review 55, 364 (1939).
Radhakrishnan (1990) T. P. Radhakrishnan, Journal of Mathematical Chemistry 5, 381 (1990).
Janakiraman (2018) D. Janakiraman, Journal of Chemical Education 95, 2086 (2018), https://doi.org/10.1021/acs.jchemed.8b00377 .
Bauswein et al. (2019) A. Bauswein, N.-U. F. Bastian, D. B. Blaschke, et al., Phys. Rev. Lett. 122, 061102 (2019).
Zhou et al. (2024) W. Zhou, H. Shen, J. Hu, and Y. Zhang, Phys. Rev. D 110, 043017 (2024).
Annala et al. (2020) E. Annala, T. Gorda, A. Kurkela, et al., Nature Physics 16, 907 (2020).
Baym et al. (1971) G. Baym, C. Pethick, and P. Sutherland, ApJ 170, 299 (1971).
Hebeler et al. (2013) K. Hebeler, J. M. Lattimer, C. J. Pethick, and A. Schwenk, The Astrophysical Journal 773, 11 (2013).
Fraga et al. (2014) E. S. Fraga, A. Kurkela, and A. Vuorinen, The Astrophysical Journal 781, L25 (2014).
Ferrari et al. (2010) L. Ferrari, P. C. R. Rossi, and M. Malheiro, International Journal of Modern Physics D 19, 1569 (2010).
Kurkela et al. (2014) A. Kurkela, E. S. Fraga, J. Schaffner-Bielich, and A. Vuorinen, The Astrophysical Journal 789, 127 (2014).
Abgaryan et al. (2018) V. Abgaryan, D. Alvarez-Castillo, A. Ayriyan, et al., Universe 4, 94 (2018), arXiv:1807.08034 [astro-ph.HE] .
Mindel (1962) J. Mindel, J. Chem. Educ. 39, 512 (1962).
Blundell and Blundell (2009) S. J. Blundell and K. M. Blundell, Concepts in Thermal Physics (Oxford University Press, 2009).

Appendix A Thermodynamic Consistency

The chemical potential $\mu$ , is related to the Gibbs free energy per particle $G$ as,

G=F+PV

where $F$ is Helmholtz free energy, $P$ is the pressure, V is the volume. For a system with particle number $N$ , it can also be written as,

G=\mu N

which, by the Euler theorem, shows that the Gibbs free energy per particle is equivalent to the chemical potential $\mu$ Mindel (1962). The partial derivative of the chemical potential with respect to pressure at constant temperature is,

\left(\frac{\partial\mu}{\partial P}\right)_{T}=\frac{V}{N}=v

(A1)

where $v$ is the volume per particle. Since $v$ is positive by definition, $\mu$ must be a strictly increasing function of pressure at a constant temperature. So, for a stable thermodynamic system, the Gibbs free energy is a strictly concave function of pressure. This concavity ensures that there is thermodynamic stability in equilibrium. In the context of a first-order phase transition, this concavity implies that the chemical potentials of two homogeneous phases near the transition point intersect. At the crossing point, the phase with the lower Gibbs free energy (or chemical potential) is the more stable, representing the minimum Gibbs potential for the system (Radhakrishnan, 1990; Blundell and Blundell, 2009; Janakiraman, 2018).

Appendix B M-R curve for transition below 2 $n_{0}$

Figure 4 demonstrates the region within the M-R space that is spanned by the SoSI method, compared to the region spanned by the PP (inconsistent) method, for PT occurring at densities below 2 $n_{0}$ . The M-R contour generated from the sequences satisfies the astrophysical bounds of M $>2$ M_⊙ and the tidal deformability bound of $\tilde{\Lambda}<720$ . Notably, the SoSI method encompasses an area similar to that associated with transitions above $2n_{0}$ . In contrast, the inconsistent PP method exhibits a more confined contour due to its inherently polytropic characteristics. This results in an even more limited region within the space. Also, the PP-consistent method fails to generate any EoS when PT happens below $2n_{0}$ .

Appendix C Comparison with PP method

In this section, we compare the EoS ensemble obtained using SoSI, PP-consistent and PP (inconsistent) methods. We also show the M-R region, encompassed by three different construction methods. In the left panel, the plot illustrates the relationship between pressure ( $P$ ) and central energy density ( $\epsilon$ ) on a logarithmic scale. The SoSI method (blue contour) offers greater flexibility, allowing a broader spectrum of phase transition scenarios. The wider spread of EoS curves at lower energy densities (up to approximately 1000 MeV/fm³) highlights the uncertainties associated with nuclear matter properties. This broader distribution at low densities reflects SoSI capturing a wider range of potential EoS behaviours. The polytropic method (PP and PP-consistent) employs predefined functional forms, producing more tightly clustered M-R curves and indicating a more constrained range of neutron star properties.