1 Introduction↩︎

Neutron stars are extremely dense celestial objects formed from the gravitational collapse of massive stars. Their core densities can reach several times the nuclear saturation density, providing a unique natural laboratory for studying strongly interacting matter under extreme conditions [1], [2]. Consequently, precisely determining the Equation of State (EOS) of their internal matter is a key challenge at the intersection of nuclear physics and astrophysics, as it directly governs the macroscopic structure and properties of neutron stars [3]–[5].

In recent years, advances in multi-messenger astronomy have provided unprecedented opportunities to constrain the EOS, while also revealing its inherent complexity. On one hand, precise observations of massive pulsars, particularly PSR J0740+6620 (\(M \approx 2.08 M_{\odot}\)) [6], [7], require the EOS to be sufficiently "stiff" at high densities to support neutron stars exceeding two solar masses [8]–[10]. On the other hand, the tidal deformability parameter (\(\Lambda\)) inferred from the gravitational wave signal of the binary neutron star merger GW170817 demands that the EOS for a \(1.4 M_{\odot}\) neutron star be relatively "soft" at corresponding densities, with \(\Lambda_{1.4} \lesssim 800\) [5], [11]–[14]. Reconciling this tension between high-density stiffness and intermediate-density softness within a unified physical framework is a central problem in neutron star physics.

A possible solution is to introduce a phase transition from hadronic matter to quark matter in the core of the neutron star, forming what is known as a "hybrid star" [3]–[5]. In this study, we employ the Nambu-Jona-Lasinio (NJL) model, which effectively describes chiral symmetry breaking and restoration, to construct the quark matter EOS [15], [16]. By combining the NJL model with a well-established hadronic EOS, we construct a complete hybrid equation of state. We pay special attention to the effects of quark vector interactions within the model, as this interaction provides the necessary repulsion at high densities and is a key mechanism for stiffening the EOS to support massive neutron stars [10], [17]. The goal of this research is to systematically explore the parameter space of this hybrid star model to construct an EOS that is not only theoretically self-consistent but also capable of simultaneously passing tests from various astronomical observations, including mass, radius, and tidal deformability.

In the subsequent sections, we present a detailed description of the NJL quark model and the BSR6 hadronic model, along with the interpolation method used to ensure \(C^2\) continuity of the EOS. We then present the macroscopic properties of neutron stars predicted by our benchmark parameter set and compare these predictions in detail with astronomical observations. Following this, we systematically analyze two key parameters of the model—the vector coupling constant \(G_V\) and the phase transition endpoint \(BU\)—to reveal their distinct regulatory mechanisms governing the macroscopic properties of neutron stars. Finally, we summarize the findings of this paper.

2 NJL Model and Quark Matter Equation of State↩︎

The Nambu-Jona-Lasinio (NJL) model is an effective quantum field theory used to describe the strong interactions between quarks, particularly well-suited for studying the phase transition from hadronic to quark matter in nuclear matter and dense stars [15], [16]. It offers distinct advantages in describing phenomena such as spontaneous chiral symmetry breaking and restoration, the generation of dynamical quark masses, and the density-dependent nature of quark masses [15], [18]. This section will detail the fundamental structure of the two-flavor NJL model, its parameter choices, and the calculation of the quark matter EOS at zero temperature and finite chemical potential, thereby establishing a foundation for the subsequent study of neutron star structure.

2.1 Two-Flavor NJL Model Lagrangian↩︎

This research primarily focuses on the two-flavor (\(N_f=2\)) NJL model, which includes up (u) and down (d) quarks. Its Lagrangian can be written as: \[\mathcal{L} = \overline{\psi}(i\gamma^\mu \partial_\mu - \hat{m})\psi + \mathcal{L}_{\text{int}} \label{eq:lagrangian95total}\tag{1}\] Here, \(\psi\) represents the quark field, \(\gamma^\mu\) are the Dirac matrices, and \(\hat{m}\) is the current quark mass matrix. In standard configurations, the current masses for u and d quarks are considered equal, i.e., \(m_u = m_d\) [17].

The interaction term \(\mathcal{L}_{\text{int}}\) consists of a four-fermion contact interaction, with a structure designed to capture the key symmetries of Quantum Chromodynamics (QCD) in the low-energy regime. We primarily consider the following two important interaction channels [15], [17]:

1. Scalar-Pseudoscalar Channel: This term describes the attractive force responsible for spontaneous chiral symmetry breaking and is intimately related to the formation of mesons (such as \(\pi\) mesons). \[\mathcal{L}_{\sigma}^{(4)} = G_S[(\overline{\psi}\psi)^2 + (\overline{\psi}i\gamma_5\vec{\tau}\psi)^2] \label{eq:lagrangian95sigma}\tag{2}\] In this equation, \(G_S\) is the scalar coupling constant, and \(\vec{\tau}\) are the Pauli matrices acting in flavor space, representing the quark isospin degrees of freedom.

2. Vector Channel: This term accounts for the short-range repulsive force between quarks, which significantly affects the pressure of quark matter and the structure of compact stars [10], [17]. \[\mathcal{L}_{V}^{(4)} = -G_V(\overline{\psi}\gamma^\mu\psi)^2 \label{eq:lagrangian95vector}\tag{3}\] Here, \(G_V\) is the vector coupling constant. A positive \(G_V\) signifies a repulsive interaction.

Consequently, the total interaction Lagrangian is \(\mathcal{L}_{\text{int}} = \mathcal{L}_{\sigma}^{(4)} + \mathcal{L}_{V}^{(4)}\). This specific form of the NJL model is capable of capturing several non-perturbative features of low-energy QCD, including chiral symmetry breaking and the dynamical generation of quark masses [15].

The parameter selection for the NJL model in this study maintains physical consistency with the BSR6 model used for the hadronic phase. The BSR6 model predicts a nuclear matter saturation density of \(n_0 \approx 0.149 \text{ fm}^{-3}\), which we use as a reference benchmark. For this study, we establish a benchmark parameter set: the current quark mass is \(m_u = m_d = 5.50 \text{ MeV}\), the three-momentum cutoff is \(\Lambda = 660.0 \text{ MeV}\), the scalar coupling is constrained by \(G_S\Lambda^2 = 1.82952\), and the vector coupling ratio is \(G_V/G_S = 0.45\). The phase transition window for this benchmark is set by \(BL=1.0\) and \(BU=5.50\).

2.2 Chiral Symmetry Breaking and Quark Dynamical Mass↩︎

A central physical mechanism of the NJL model is the spontaneous breaking of chiral symmetry. This process allows quarks to acquire a substantial dynamical mass from their tiny current quark masses, which helps to explain why the quarks that form hadrons (like nucleons) appear to have a much larger effective mass [15], [16], [18]. Within the Mean-Field Approximation (MFA), the effective quark mass, \(M\) (also known as the constituent quark mass), is self-consistently determined by the following chiral (or "gap") equation [15]: \[M_f = m_f - 2G_S \langle \overline{\psi}_f \psi_f \rangle \label{eq:gap95equation}\tag{4}\] Here, \(m_f\) is the current quark mass for flavor \(f\), and \(\langle \overline{\psi}_f \psi_f \rangle\) is the expected value of the quark condensate for that flavor. Both the quark condensate \(\langle \overline{\psi}_f \psi_f \rangle\) and the quark number density \(\rho_f\) depend on the temperature \(T\) and the effective quark chemical potential \(\mu_f^*\) [15]. The effective chemical potential \(\mu_f^*\) accounts for the interaction between quarks and their environment, and it differs from the physical chemical potential \(\mu_f\), especially when vector interactions are present [17].

For a two-flavor NJL model, the expression for the quark condensate typically involves an integral over momentum space. Due to the non-renormalizable nature of the model, a cutoff parameter \(\Lambda\) (such as a three-momentum cutoff or Proper-Time Regularization) is necessary to handle divergent integrals [15]. At finite temperature and chemical potential, the gap equation includes Fermi-Dirac distribution functions to account for medium effects: \[M_f = m_f + 4N_c G_S \int \frac{d^3 p}{(2\pi)^3} \frac{M_f}{E_p} (1 - n_p(E_p, \mu_f^*) - \bar{n}_p(E_p, \mu_f^*)) \label{eq:gap95equation95finite95T}\tag{5}\] where \(E_p = \sqrt{\vec{p}^2 + M_f^2}\) is the quark energy, and \(n_p\) and \(\bar{n}_p\) are the Fermi-Dirac distribution functions for quarks and antiquarks, respectively [15]. At zero temperature, these distribution functions simplify to step functions.

The introduction of the vector interaction term \(\mathcal{L}_{V}^{(4)} = -G_V(\overline{\psi}\gamma^\mu\psi)^2\) establishes a clear link between the physical quark chemical potential \(\mu_f\) and its effective chemical potential \(\mu_f^*\). The effective chemical potential \(\mu_f^*\) incorporates the mean-field interaction effects between quarks and the vector meson field, with the following expression [15], [17]: \[\mu_f^* = \mu_f - 2G_V \rho_f \label{eq:effective95chemical95potential}\tag{6}\] Here, \(\mu_f\) is the physical chemical potential for quark flavor \(f\), \(G_V\) is the vector coupling constant, and \(\rho_f = \langle \psi_f^\dagger \psi_f \rangle\) is the corresponding quark number density. The number density \(\rho_f\) for quarks at zero temperature is given by: \[\rho_f = \frac{N_c}{\pi^2} \int_0^{p_{F,f}} p^2 dp = \frac{N_c p_{F,f}^3}{3\pi^2} \quad \text{where } p_{F,f} = \sqrt{(\mu_f^*)^2 - M_f^2} \label{eq:quark95density95expression}\tag{7}\] Here, \(N_c=3\) is the number of colors, and \(p_{F,f}\) is the Fermi momentum. When \(\mu_f^* < M_f\), the Fermi momentum is 0, and the number density \(\rho_f\) is also 0. This relationship highlights that the presence of vector interactions means the actual energy state of quarks in the medium (described by the effective chemical potential) differs from the externally applied physical chemical potential. This interaction, which is generally repulsive, "offsets" a portion of the physical chemical potential, requiring a higher physical chemical potential to reach the same effective chemical potential state for the quarks.

2.3 Charge Neutrality and Beta-Equilibrium Conditions↩︎

Matter within a neutron star must satisfy specific equilibrium conditions to remain stable under its extreme conditions. For quark matter, two fundamental conservation laws are charge neutrality and beta-equilibrium [5], [15]. These conditions impose strict constraints on the relationship between quark flavors and chemical potentials, which in turn profoundly impacts the quark matter EOS and the macroscopic properties of neutron stars.

1. Charge Neutrality Condition: Given the nature of the strong interaction, quark matter must maintain overall charge neutrality to prevent the accumulation of immense Coulomb energy [15]. This requires that the total charge density from quarks and leptons (such as electrons) must be zero. For a two-flavor (u and d) quark system, the charge neutrality condition is expressed as [5]: \[\frac{2}{3}\rho_u - \frac{1}{3}\rho_d - \rho_e = 0 \label{eq:charge95neutrality}\tag{8}\] where \(\rho_u\) and \(\rho_d\) are the number densities of up and down quarks, respectively, and \(\rho_e\) is the number density of electrons. The contribution of electrons, as leptons, cannot be neglected. The electron number density at zero temperature is determined by its chemical potential \(\mu_e\) [5]: \[\rho_e(\mu_e) = \frac{\mu_e^3}{3\pi^2} \label{eq:electron95density}\tag{9}\]

2. Beta-Equilibrium Condition: After a neutron star is formed, its internal matter achieves thermodynamic equilibrium through weak interaction processes, referred to as beta-equilibrium. These processes involve the interconversion of quarks, as well as quarks and leptons. For two-flavor (u, d) quark matter, the primary weak interaction processes are [5]: \[\begin{align} d &\leftrightarrow u + e^- + \bar{\nu}_e \tag{10} \\ u + e^- &\leftrightarrow d + \nu_e \tag{11} \end{align}\] Assuming neutrinos (\(\nu_e\)) are not trapped and can freely escape the star (which is the case for an old, cooled neutron star), these weak interactions lead to the following relationship between the chemical potentials of quarks and electrons [5]: \[\mu_d = \mu_u + \mu_e \label{eq:beta95equilibrium}\tag{12}\] When the charge neutrality condition (8 ) is combined with the beta-equilibrium condition (12 ), only one independent variable remains among the quark chemical potentials. We typically choose the up-quark chemical potential \(\mu_u\) as this variable and derive the expressions for \(\mu_d\) and \(\mu_e\) from these relationships. For instance, the down-quark chemical potential is \(\mu_d = \mu_u + \mu_e\), where the value of \(\mu_e\) is determined self-consistently by the charge neutrality condition [5], [15].

When calculating the quark matter EOS, these conditions must be satisfied simultaneously. This means the number densities of quarks (\(\rho_u, \rho_d\)) and electrons (\(\rho_e\)) are interconnected, collectively determining the system’s pressure and energy density. This self-consistent calculation is a critical step in understanding the complex phase structure of matter within neutron stars.

2.4 Thermodynamic Potential and Equation of State Calculation↩︎

At zero temperature (\(T=0\)) and finite chemical potential, the macroscopic properties of quark matter are determined by its Grand Canonical Potential, also known as the thermodynamic potential \(\Omega\) [15], [17]. By integrating over the quark energy spectrum, we can derive the expression for this potential. In the mean-field approximation, the total thermodynamic potential of the system includes contributions from both quarks and electrons [15].

The general form of the total thermodynamic potential \(\Omega(T, \mu; M, \tilde{\mu})\) is given by [15]: \[\Omega(T,\mu;M,\tilde{\mu}) = \Omega_M(T,\tilde{\mu}) + \frac{(M-m)^2}{4G_S} - \frac{(\mu-\tilde{\mu})^2}{4G_V} + \text{const.} \label{eq:omega95total}\tag{13}\] Here, \(\Omega_M(T=0, \tilde{\mu})\) represents the contribution from a free Fermi gas (quarks and antiquarks) at zero temperature. For two-flavor quarks (\(N_f=2\)), its expression is: \[\Omega_M(T=0, \tilde{\mu}) = -\frac{N_c}{24\pi^2} \sum_{f=u,d} \left[ \tilde{\mu}_f \sqrt{\tilde{\mu}_f^2 - M_f^2} (2\tilde{\mu}_f^2 - 5M_f^2) + 3M_f^4 \ln\left(\frac{\tilde{\mu}_f + \sqrt{\tilde{\mu}_f^2 - M_f^2}}{M_f}\right) \right] \label{eq:omega95M95quark95expression}\tag{14}\] where the summation \(\sum_{f=u,d}\) runs over up and down quarks, and \(N_c=3\) is the number of colors.

The self-consistent equations are then obtained by minimizing the thermodynamic potential with respect to its auxiliary variables (such as \(M\) and \(\tilde{\mu}\)) [15]. Once a stable self-consistent solution is found, fundamental EOS quantities like pressure \(P\) and energy density \(\epsilon\) can be derived using standard thermodynamic relations [15], [17]: \[\begin{align} P &= -\Omega \tag{15} \\ \epsilon &= \sum_f \mu_f \rho_f - P \tag{16} \end{align}\] Here, \(\mu_f\) and \(\rho_f\) are the chemical potential and particle number density for quark flavor \(f\). The summation \(\sum_f\) includes all existing quark flavors, which in this study are primarily the u and d quarks. For neutron star matter, the contribution of leptons (e.g., electrons) must also be considered, so the energy density is more accurately expressed as \(\epsilon = \sum_i \mu_i \rho_i - P\), where \(i\) runs over all constituent particles (quarks and leptons). The baryon number density \(\rho_B\) is given by the derivative of pressure with respect to the baryon chemical potential: \(\rho_B = \frac{\partial P}{\partial \mu_B}\).

An important concept in the NJL model is the "Bag Constant," \(B\), which can be physically interpreted as the difference in vacuum energy density resulting from chiral symmetry breaking [3], [15]. In EOS calculations, the bag constant is often treated as a phenomenological parameter that sets the energy density of quark matter at zero pressure and indirectly accounts for quark confinement effects. This study will rely on the widely accepted NJL model parameter settings and calculation methods found in existing literature to ensure the reliability and comparability of our results [10], [17].

3 Hadronic and Hybrid Equation of State Construction↩︎

A central challenge in understanding the internal structure and evolution of neutron stars is the accurate construction of an Equation of State (EOS) that describes matter under extreme conditions. The density within a neutron star varies drastically, from the relatively low densities in its crust to ultra-high densities in the core that can far exceed the nuclear saturation density. No single physical model can comprehensively cover such a wide range. Therefore, this study employs a method of layered construction and smooth interpolation to combine hadronic and quark EOSs, with the aim of creating a thermodynamically consistent hybrid EOS across the entire density range.

3.1 Hadronic Equation of State Selection↩︎

The outer regions of a neutron star are composed of nuclear matter, and its EOS determines the physical properties of the stellar crust and outer core. Based on the density, we have carefully selected and processed the hadronic EOS:

1. Crust EOS: In the outermost layer, or the crust, the matter density is relatively low, and atomic nuclei and electrons form a lattice structure [3]. This study uses the Baym-Pethick-Sutherland (BPS) model [19] to describe the EOS in this region. The BPS model considers the lattice effects of the nuclei and the contribution of free electrons, providing a basis for calculating pressure and energy density at low densities.

2. Outer Core Hadronic EOS: As the depth increases, the density of matter gradually rises, reaching and surpassing the nuclear saturation density, forming the liquid core of the neutron star. To accurately describe the nuclear matter in this region, we selected the BSR6 model, which is based on the relativistic mean-field (RMF) theory [20]. The BSR6 model describes the properties of nuclear matter by introducing effective interactions between nucleons and mesons (such as \(\sigma\), \(\omega\), and \(\rho\) mesons). It also includes self-interaction and mixed-interaction terms for the meson fields to provide a more comprehensive description of the complex behavior of nuclear matter [20], [21]. Its nuclear properties parameters, such as the saturation density \(n_0=0.149 \text{ fm}^{-3}\), binding energy \(E/A=-16.1 \text{ MeV}\), incompressibility \(K=235.8 \text{ MeV}\), and symmetry energy \(J=35.6 \text{ MeV}\), are all consistent with experimental data and theoretical constraints [20]. It is important to note that RMF models do not typically include explicit hyperon degrees of freedom. This is because the interactions between hyperons and nucleons, as well as hyperon-hyperon interactions, are still subject to considerable uncertainty, and most models predict that hyperons begin to appear at densities around \(n_B \sim 2-3n_0\), a range that coincides with the hadron-quark crossover region. Therefore, their impact on the EOS must be considered as part of the hybrid state construction.

In RMF theory, the Lagrangian density for hadronic matter typically includes a nucleonic component, meson self-interaction components, and mixed-interaction terms between mesons [21], [22]. For the BSR6 model, the Lagrangian density can be expressed as: \[\mathcal{L}=\mathcal{L}_{NM}+\mathcal{L}_{\sigma}+\mathcal{L}_{\omega}+\mathcal{L}_{\rho}+\mathcal{L}_{\sigma\omega\rho} \label{eq:hadronic95lagrangian}\tag{17}\] where:

• \(\mathcal{L}_{NM}\) is the nucleonic part of the Lagrangian, which describes the free motion of nucleons (neutrons \(n\) and protons \(p\)) and their coupling to the meson fields: \[\mathcal{L}_{NM}=\sum_{H=n,p}\overline{\psi}_{H}[i\gamma^{\mu}\partial_{\mu}-(M-g_{\sigma}\sigma)-(g_{\omega}\gamma^{\mu}\omega_{\mu}+\frac{1}{2}g_{\rho}\gamma^{\mu}\vec{\tau}\cdot\vec{\rho}_{\mu})]\psi_{H}\] Here, \(\psi_H\) represents the nucleon field, \(M\) is the nucleon mass, \(g_{\sigma}, g_{\omega}, g_{\rho}\) are the coupling constants for the nucleons to the \(\sigma, \omega, \rho\) meson fields, respectively, and \(\vec{\tau}\) is the isospin matrix.

• \(\mathcal{L}_{\sigma}\), \(\mathcal{L}_{\omega}\), and \(\mathcal{L}_{\rho}\) describe the dynamics and self-interaction terms of the \(\sigma\) meson (scalar-isoscalar), \(\omega\) meson (vector-isoscalar), and \(\rho\) meson (vector-isovector) fields. For example, the \(\sigma\) meson term includes a mass term and nonlinear self-coupling terms: \[\mathcal{L}_{\sigma}=\frac{1}{2}(\partial^{\mu}\sigma\partial_{\mu}\sigma-m_{\sigma}^{2}\sigma^{2})-\frac{\kappa_{3}}{6M}g_{\sigma}^{3}m_{\sigma}^{2}\sigma^{3}-\frac{\kappa_{4}}{24M^{2}}g_{\sigma}^{4}m_{\sigma}^{2}\sigma^{4}\] where \(m_{\sigma}\) is the \(\sigma\) meson mass, and \(\kappa_3, \kappa_4\) are nonlinear coupling coefficients. For the \(\omega\) meson term \(\mathcal{L}_{\omega}\), its form typically includes: \[\mathcal{L}_{\omega} = -\frac{1}{4}F_{\omega}^{\mu\nu}F_{\omega, \mu\nu} + \frac{1}{2}m_{\omega}^{2}\omega^{\mu}\omega_{\mu} + \frac{\zeta_0}{4!} (g_{\omega}\omega^{\mu}\omega_{\mu})^2\] Here, \(F_{\omega}^{\mu\nu} = \partial^{\mu}\omega^{\nu} - \partial^{\nu}\omega^{\mu}\) is the field strength tensor for the \(\omega\) meson, \(m_{\omega}\) is the \(\omega\) meson mass, and \(\zeta_0\) is its nonlinear self-coupling coefficient. For the \(\rho\) meson term \(\mathcal{L}_{\rho}\), its form typically includes: \[\mathcal{L}_{\rho} = -\frac{1}{4}\vec{F}_{\rho}^{\mu\nu}\cdot\vec{F}_{\rho, \mu\nu} + \frac{1}{2}m_{\rho}^{2}\vec{\rho}^{\mu}\cdot\vec{\rho}_{\mu} + \frac{\xi_0}{4!} (g_{\rho}\vec{\rho}^{\mu}\cdot\vec{\rho}_{\mu})^2\] Here, \(\vec{F}_{\rho}^{\mu\nu} = \partial^{\mu}\vec{\rho}^{\nu} - \partial^{\nu}\vec{\rho}^{\mu} - g_{\rho}(\vec{\rho}^{\mu}\times\vec{\rho}^{\nu})\) is the field strength tensor for the \(\rho\) meson, \(m_{\rho}\) is the \(\rho\) meson mass, and \(\xi_0\) is its nonlinear self-coupling coefficient. These self-interaction terms are vital for describing the saturation behavior of mesons at high densities, particularly for accurately characterizing the properties of asymmetric nuclear matter [21]–[23].

• \(\mathcal{L}_{\sigma\omega\rho}\) describes the mixed interaction terms between the meson fields, which are crucial for precisely characterizing nuclear matter properties (especially the density dependence of the symmetry energy) [22]: \[\begin{align} \mathcal{L}_{\sigma\omega\rho} &= \frac{\eta_{1}}{2M}g_{\sigma}m_{\omega}^{2}\sigma\omega^{\mu}\omega_{\mu} + \frac{\eta_{2}}{4M^{2}}g_{\sigma}^{2}m_{\omega}^{2}\sigma^{2}\omega^{\mu}\omega_{\mu} + \frac{\eta_{3}}{2M}g_{\sigma}m_{\rho}^{2}\sigma\rho^{\mu}\rho_{\mu} \\ &+ \frac{\eta_{4}}{4M^{2}}g_{\sigma}^{2}m_{\rho}^{2}\sigma^{2}\rho^{\mu}\rho_{\mu} + \frac{\eta_{5}}{4M^{2}}g_{\omega}^{2}m_{\rho}^{2}\omega^{\mu}\omega_{\mu}\rho^{\mu}\rho_{\mu} \end{align}\] These \(\eta\) coefficients are phenomenologically determined in the RMF model by fitting the ground-state properties of finite nuclei and the nuclear matter parameters at saturation density to optimize the model.

3. Crust and Core Stitching: To ensure that the entire hadronic EOS is smooth and thermodynamically consistent across the full density range, we use a quintic polynomial interpolation to smoothly connect the crust EOS with the core hadronic EOS. This method ensures the continuity of the pressure \(P\), energy density \(\epsilon\), and their first and second derivatives with respect to pressure (\(\partial\epsilon/\partial P\) and \(\partial^2\epsilon/\partial P^2\)) at the connection point, thereby preventing unphysical jumps or abrupt changes [5], [21]. The stitched hadronic EOS then serves as the baseline for the low-density region in the subsequent construction of the hybrid EOS.

3.2 Hadron-Quark Hybrid Equation of State Construction↩︎

The extreme density environment in the core of a neutron star may induce a deconfinement phase transition of hadronic matter, leading to the formation of quark matter composed of free quarks and gluons. This gives rise to the possibility of "hybrid stars" [3]–[5]. This study employs a "three-window" approach [3], [10] to construct the hadron-quark hybrid EOS. This method allows for a smooth transition, or "crossover," region between the hadronic and quark phases.

This approach avoids potential issues like pressure discontinuities or unphysical sound speeds associated with first-order phase transitions and better represents a continuous evolution from hadrons to quarks [3], [10]. We divide the EOS into three regions based on the baryon chemical potential \(\mu_B\):

1. Low-Density Hadronic Region (\(\mu_B < \mu_{BL}\)): In this region, matter exists entirely in the hadronic phase, which is described by the stitched hadronic EOS (including the crust and core parts) as selected in the previous subsection. \(\mu_{BL}\) is the specified upper boundary chemical potential for the hadronic phase, typically corresponding to a baryon number density of \(n_B \sim (1-2)n_0\). Above this density, the reliability of traditional hadronic models may diminish [3], [10].

2. High-Density Quark Region (\(\mu_B > \mu_{BU}\)): In this region, matter is assumed to be fully deconfined quark matter. We use the two-flavor NJL model as described in Section 2, which incorporates quark vector interactions and chiral symmetry breaking effects. \(\mu_{BU}\) is the defined lower boundary chemical potential for the quark phase, generally corresponding to a baryon number density of \(n_B \sim (4-7)n_0\). Below this density, quark confinement effects become significant, limiting the applicability of a quark model [3], [10].

3. Intermediate Transition Region (\(\mu_{BL} \le \mu_B \le \mu_{BU}\)): This is a hadron-quark mixed phase or a continuous crossover region, which is difficult to calculate precisely from first principles. This study uses a phenomenological interpolation method to describe the EOS in this region. We use a quintic polynomial \(\mathcal{P}(\mu_B) = \sum_{m=0}^{5}C_{m}\mu_{B}^{m}\) to connect the pressure-baryon chemical potential relationship of the hadronic and quark phases (as detailed in Appendix 8) [3], [5], [10]. The polynomial coefficients \(C_m\) are determined by applying the following boundary conditions, which ensure the thermodynamic consistency and smoothness of the entire hybrid EOS:

   • At \(\mu_B = \mu_{BL}\), the pressure \(P(\mu_B)\) and its first two derivatives with respect to \(\mu_B\) (i.e., the baryon number density \(\rho_B(\mu_B)\) and the baryon number susceptibility \(\partial \rho_B / \partial \mu_B\)) must match those of the hadronic EOS.

   • At \(\mu_B = \mu_{BU}\), the pressure \(P(\mu_B)\) and its first two derivatives with respect to \(\mu_B\) must match those of the quark EOS.

   These conditions ensure the continuity of pressure, number density, and number susceptibility in the transition region, thereby avoiding unphysical jumps or unstable areas [21].

This piecewise construction and interpolation method yields a hybrid EOS, \(P(\mu_B)\), that smoothly transitions across the entire density range and satisfies fundamental physical constraints [3], [10]:

• Pressure Continuity: The pressure \(P\) must be a continuous function of the baryon chemical potential \(\mu_B\).

• Thermodynamic Stability: The baryon number density \(\rho_B = \partial P / \partial \mu_B\) must be a monotonically increasing function of \(\mu_B\), i.e., \(\partial^2 P / \partial \mu_B^2 > 0\), to ensure the system’s stability against density fluctuations.

• Causality: The speed of sound squared, \(v_s^2 = \partial P / \partial \epsilon\), in the medium must be less than or equal to the speed of light squared, \(c^2\), i.e., \(v_s^2/c^2 \le 1\). This is a fundamental physical requirement that signals cannot propagate faster than light.

A hybrid EOS that satisfies these constraints will provide a reliable physical input for subsequent calculations of neutron star structure and properties.

3.3 Speed of Sound Squared \(v_s^2\) and Equation of State Analysis↩︎

The speed of sound is a key physical quantity that describes how quickly a disturbance propagates through a medium and reflects the matter’s responsiveness to small perturbations. Inside a compact star, the speed of sound squared, \(v_s^2\), is a crucial feature of the EOS, defined as the derivative of the pressure \(P\) with respect to the energy density \(\epsilon\): \[v_s^2 = \frac{dP}{d\epsilon} \label{eq:speed95of95sound95new}\tag{18}\] The magnitude of the speed of sound squared provides an intuitive measure of the EOS’s "stiffness": a higher speed of sound indicates a "stiffer" EOS, meaning the matter is less compressible under changes in pressure. Within a relativistic framework, the speed of physical signals cannot exceed the speed of light. Therefore, the speed of sound squared must be less than or equal to the speed of light squared, \(c^2\), i.e., \(v_s^2/c^2 \le 1\). This causality condition is a basic criterion for verifying the physical validity of any EOS [24]. The behavior of the sound speed in the hadron-quark phase transition region is particularly key to understanding the process of matter’s structural transformation.

The parameters \(BL\) and \(BU\), which appear frequently in this paper and its figures, are baryon number density coefficients that define the boundaries of the hadron-quark mixed phase region:

• \(BL\) (Lower Boundary Coefficient): This coefficient denotes the baryon number density at which the hadron-quark mixed phase begins. Specifically, when the baryon number density reaches \(BL \times n_0\) (where \(n_0\) is the nuclear saturation density), hadronic matter starts to transition into quark matter.

• \(BU\) (Upper Boundary Coefficient): This coefficient denotes the baryon number density at which the hadron-quark mixed phase ends. Specifically, when the baryon number density reaches \(BU \times n_0\), the matter is considered to have completely transformed into quark matter.

• By adjusting the combination of \(BL\) and \(BU\), we can simulate equations of state with different phase transition starting and ending densities, and subsequently explore their impact on the macroscopic properties of neutron stars.

Figure 1 (a) shows the trend of the speed of sound squared \(v_s^2/c^2\) as a function of energy density \(\epsilon\) for a typical parameter set (Set 1) in this study. The "Hadronic" (black solid line) curve shows that the sound speed of hadronic matter gradually increases with energy density, exhibiting a "stiffer" characteristic. In contrast, the "Quark (Set 1)" (blue dotted line) curve has a relatively lower and flatter sound speed, reflecting the "softer" compressibility of quark matter. The core of this study, the "Hybrid (Set 1)" (cyan dot-dashed line) curve, clearly demonstrates a continuous transition from the hadronic to the quark phase. In the transition region (around \(500 \text{ MeV/fm}^3\)), the speed of sound squared shows a distinct peak followed by a rapid drop after a brief rise. This "softening" is a typical feature of a hadron-quark phase transition, reflecting the sharp change in the system’s compressibility during the structural rearrangement of matter. Notably, most EOS curves constructed in this study satisfy the causality condition \(v_s^2/c^2 \le 1\). However, as will be discussed in Section 6, certain parameter choices for the hadron-quark transition window can lead to a violation of this condition, providing an important physical constraint on the model parameters.

3.4 Pressure-Energy Density Relation Analysis↩︎

The Equation of State (EOS) is a fundamental physical quantity that describes the constitutive relationship between a material’s pressure (\(P\)) and its energy density (\(\epsilon\)). It is this relationship that fundamentally determines a neutron star’s macroscopic structure and properties. By analyzing the \(P-\epsilon\) relation plot, we can gain an intuitive understanding of the "stiffness" or "softness" of matter in different density regions, which is crucial for predicting the neutron star’s mass, radius, and its response to external perturbations [24], [25].

Figure 1 (b) provides a comprehensive illustration of the pressure-energy density relationship for the hybrid EOS constructed with a typical parameter set (Set 1). This plot clearly depicts the evolution of the matter phase inside a neutron star, from low to high density. The "Hadronic" (black solid line) curve represents hadronic matter, where pressure rises rapidly with energy density, indicating a relatively stiff EOS. In contrast, the "Quark (Set 1)" (blue dotted line) curve shows a lower pressure at the same energy density, reflecting the relatively "softer" nature of quark matter. The central "Hybrid (Set 1)" (cyan dot-dashed line) curve represents our constructed hybrid EOS. In the low-energy density region, this curve precisely overlaps with the hadronic curve, indicating that matter is in the hadronic phase. As the energy density increases, the hybrid curve gradually deviates from the hadronic curve and smoothly transitions towards the quark curve at high energies. This smooth connection visually confirms the effectiveness of our quintic polynomial interpolation method in constructing a thermodynamically self-consistent hadron-quark crossover, providing theoretical support for the possible existence of a continuous phase transition inside neutron stars [3], [10].

4 Neutron Star Structure Calculation↩︎

The macroscopic properties of a neutron star-such as mass, radius, and deformability under an external gravitational field, characterized by the tidal deformability parameter-are directly governed by the Equation of State (EOS) of its internal matter. To derive these macroscopic characteristics from the microscopic EOS, we must solve the stellar structure equations—typically the Tolman–Oppenheimer–Volkoff (TOV) equations—within the framework of general relativity and account for physical processes such as tidal deformation.

4.1 Tolman-Oppenheimer-Volkoff (TOV) Equations and Their Numerical Solution↩︎

For a static, spherically symmetric neutron star, its internal structure is described by the Tolman-Oppenheimer-Volkoff (TOV) equations[26], [27]. This set of equations is a simplified form of Einstein’s field equations for a spherically symmetric fluid distribution, which precisely describes how the internal pressure \(P(r)\) and enclosed mass \(M(r)\) change with the radial coordinate \(r\): \[\begin{align} \frac{dP}{dr} &= -\frac{G(\epsilon + P/c^2)(M(r) + 4\pi r^3 P/c^2)}{r(r - 2GM(r)/c^2)} \tag{19} \\ \frac{dM}{dr} &= 4\pi r^2 \epsilon \tag{20} \end{align}\] Here, \(G\) is the gravitational constant, \(c\) is the speed of light in a vacuum, and \(\epsilon\) is the energy density, which is closely related to the pressure \(P\) through the EOS, \(\epsilon(P)\), constructed previously.

To numerically solve the TOV equations, we typically use the Runge-Kutta method for integration. The process begins at the center of the star (\(r=0\)), where an initial central pressure \(P_c\) and an initial mass of \(M(0)=0\) are set. As the radial distance \(r\) increases, the pressure \(P(r)\) gradually decreases. The integration stops when the pressure falls to a preset surface threshold. At this point, the radial distance \(r\) is the star’s radius \(R\), and \(M(R)\) is the total mass \(M\) of the neutron star corresponding to the central pressure \(P_c\). By systematically scanning a range of different central pressures \(P_c\), we can calculate multiple sets of \((R, M)\) values, which allows us to plot the complete mass-radius (\(M-R\)) relationship curve for the neutron star. This curve provides a direct visual representation of how the EOS influences the star’s macroscopic structure and can be directly compared with astronomical observational data [24], [25].

Figure 2 (a) depicts the \(M-R\) relationship for the neutron star predicted by our typical hybrid EOS (Set 1), compared with the hadronic model and the latest astronomical data. The "Hadronic" (black solid line) curve predicts a relatively large maximum mass and radius. The "Set 1" (blue solid line) hybrid EOS curve overlaps with the hadronic curve at low masses but shows a more compact structure at the high-mass end.

We compare these theoretical predictions with three representative observational constraints:

• PSR J0740+6620: A known massive pulsar with a mass of approximately \(M = 2.08 \pm 0.07 M_\odot\) and a radius of about \(12.49^{+1.28}_{-0.88}\) km [6], [7].

• PSR J0030+0451: NICER observations provide a joint constraint on its mass and radius, \(M = 1.34^{+0.15}_{-0.16} M_\odot\) and \(R = 12.71^{+1.14}_{-1.19}\) km [28].

• PSR J0437-4715: Observations of this pulsar also constrain its mass and radius to approximately \(M = 1.44 \pm 0.07 M_\odot\) and about \(11.36^{+0.95}_{-0.63}\) km [29].

As the figure shows, the M-R curve for the "Set 1" hybrid star is consistent with the observational constraints from all three pulsars, demonstrating the compatibility of our constructed hybrid EOS with current astronomical data.

4.2 Quantitative Analysis of Key Macroscopic Neutron Star Properties↩︎

To more quantitatively evaluate the physical effectiveness of the constructed EOS, we have extracted and compared the key macroscopic properties of neutron stars predicted by the hadronic model and our typical hybrid model (Set 1), as shown in Table 1.

The table clearly shows that, compared to the hadronic model, the hybrid star (Set 1) with a quark core exhibits significantly different macroscopic features. Its maximum mass (\(2.20 M_\odot\)) is slightly lower but still well above the observational lower limit. More importantly, its radius (\(R_{1.4} = 12.27\) km) and tidal deformability parameter (\(\Lambda_{1.4} = 352.97\)) at \(1.4 M_\odot\) are both significantly smaller than those of the hadronic model. This result quantitatively demonstrates that the hadron-quark phase transition makes the neutron star more compact, which is crucial for simultaneously satisfying the diverse constraints from different astronomical observations.

4.3 Tidal Deformability \(\boldsymbol{\Lambda}\) and Neutron Star Macroscopic Properties↩︎

The tidal deformability parameter \(\boldsymbol{\Lambda}\) is a key physical quantity that measures how much a neutron star deforms under an external gravitational field (for example, during a binary neutron star merger event) [30], [31]. This parameter, by influencing the phase evolution of the gravitational wave signal, provides a unique and powerful way to precisely constrain the EOS of dense matter [11], [32], [33].

In the framework of general relativity, for a static, spherically symmetric, non-rotating star, when it is perturbed by an external quadrupolar tidal field \(\mathcal{E}_{ij}\), the star induces its own quadrupolar moment \(Q_{ij}\). The Love number \(k_l\) (typically referring to the quadrupolar deformation, \(k_2\)) is defined as the dimensionless proportionality constant that links this induced quadrupolar moment to the external tidal field [34]–[36]. The tidal deformability parameter \(\Lambda\) further relates the Love number to the star’s mass \(M\) and radius \(R\) with the following specific relationship [36], [37]: \[\Lambda = \frac{2}{3} k_2 \left(\frac{c^2 R}{G M}\right)^5 \label{eq:lambda95k295relation95new}\tag{21}\] where \(R\) and \(M\) are the neutron star’s radius and mass, respectively, and \(k_2\) is the quadrupolar Love number. Calculating \(k_2\) is a critical step in determining \(\Lambda\).

The calculation of \(k_2\) involves solving the linearized Einstein equations for the star when it is subjected to a small external tidal field. First, the unperturbed spherically symmetric stellar metric (obtained from the solution of the TOV equations) is linearly perturbed. In the Regge-Wheeler gauge, for a static, even-parity (electric-type) quadrupolar (\(l=2\)) perturbation, the metric perturbation can be described by a master function \(H(r)\) [30], [38]. The master function \(H(r)\) satisfies a second-order ordinary differential equation: \[r H^{\prime\prime}(r) - H^{\prime}(r) F(r) - H(r) Q(r) = 0 \label{eq:H95equation}\tag{22}\] Here, the prime denotes a derivative with respect to the radial coordinate \(r\), and the coefficients \(F(r)\) and \(Q(r)\) are complex functions determined by the background spacetime (i.e., the solution of the TOV equations) (their explicit forms are given in Appendix 9), depending on the pressure, energy density, and enclosed mass at \(r\).

Solving this equation typically involves a numerical integration method. Starting from the center of the star (\(r \rightarrow 0\)), the initial condition for the master function \(H(r)\) is determined by the regularity requirement at the origin, which for \(l=2\) is \(H(r) = a_0 r^2\)[30], [31], where \(a_0\) is an arbitrary constant. The integration proceeds outward to the stellar surface at \(r=R\). By matching the internal numerical solution \(H(R)\) and its radial derivative \(H'(R)\) with the external analytical solution at the stellar surface \(r=R\), we can ultimately calculate the Love number \(k_2\). Specifically, \(k_2\) can be expressed as a function that depends on the star’s compactness \(C=GM/Rc^2\) and the surface logarithmic derivative \(y=RH'(R)/H(R)\) [30], [37].

The Love number \(k_2\) is highly sensitive to the "stiffness" of the EOS: typically, a stiffer EOS leads to a larger Love number, and vice versa [37], [39].

Figure 2 (b) shows the relationship between the tidal deformability parameter and neutron star mass. The gravitational wave event GW170817 placed a strong constraint of \(\Lambda_{1.4} \lesssim 800\) on the tidal deformability parameter of a \(1.4 M_\odot\) neutron star [5], [14], providing crucial astronomical data for testing high-density matter EOSs. The gray vertical dashed line at \(1.4 M_\odot\) and the red horizontal dashed line at \(\Lambda=800\) in the figure mark this critical constraint region.

The figure clearly shows that the \(\Lambda_{1.4}\) value predicted by the hadronic EOS is 814.78, which is slightly above the observational upper limit. However, the hybrid EOS (Set 1), which includes a hadron-quark phase transition, predicts a \(\Lambda_{1.4}\) value of only 352.97, well within the observationally allowed range. This compelling comparison suggests that a hadronic EOS may struggle to satisfy all observational constraints. Introducing a quark core appears to be an effective way to resolve this tension. This finding provides powerful evidence for using gravitational wave data to constrain the phase transition behavior of extreme matter.

5 Impact of the Vector Coupling Constant \(G_V\) on Neutron Star Macroscopic Properties↩︎

In the NJL model, the vector coupling constant \(G_V\) describes the repulsive interaction between quarks mediated by vector meson exchange. This interaction is particularly important at high densities, as it can effectively increase the "stiffness" of the equation of state, thereby having a decisive impact on macroscopic properties such as the maximum mass of a neutron star [10], [17]. To systematically investigate the effect of this key parameter, this section fixes all other model parameters at their benchmark values and varies only the value of \(G_V\), constructing three different sets of hybrid equations of state for a detailed comparative analysis.

5.1 Parameter Sets and Comparison of Key Properties↩︎

We select three representative values of \(G_V\) for comparison: \(G_V/G_S = 0.40, 0.45,\) and \(0.50\). The set with \(G_V/G_S=0.45\) serves as our benchmark. All parameter sets use the same phase transition window, with \(BL=1.0\) and \(BU=5.5\). Table 2 quantitatively summarizes the key macroscopic properties of neutron stars for these three parameter sets.

The table clearly shows a strong positive correlation between the value of \(G_V\) and the macroscopic properties of the neutron star. As the \(G_V/G_S\) ratio increases, the maximum mass (\(M_{max}\)), the corresponding radius (\(R_{M_{max}}\)), the radius at \(1.4M_\odot\) (\(R_{1.4}\)), and the tidal deformability parameter (\(\Lambda_{1.4}\)) all systematically increase. This directly confirms the stiffening effect of the repulsive vector interaction on the equation of state and its role in supporting the star’s structure.

5.2 Graphical Analysis of the EOS and Macroscopic Properties↩︎

To better visualize the impact of \(G_V\), Figure 3 provides a side-by-side comparison of the three parameter sets at both the equation of state level and the macroscopic properties level.

A detailed analysis of Figure 3 reveals the following:

• Speed of Sound (Panel a) clearly illustrates the mechanism by which \(G_V\) affects the EOS stiffness. In the high-energy density region after the phase transition (\(\epsilon \gtrsim 700 \text{ MeV/fm}^3\)), the magnitude of the speed of sound squared \(v_s^2/c^2\) is directly proportional to the value of \(G_V\): the larger the \(G_V\) (the \(G_V/G_S=0.50\) set), the higher the sound speed and the stiffer the EOS. This demonstrates the critical role of vector repulsion in stiffening the matter at extreme densities.

• Pressure-Energy Density Relation (Panel b) corroborates the sound speed analysis from another perspective. In the low-density region, the three hybrid curves coincide with the hadronic curve. After the phase transition point, the curves diverge, with the curve corresponding to the largest \(G_V\) (\(G_V/G_S=0.50\)) exhibiting a higher pressure and a steeper slope at the same energy density, indicating a stiffer character.

• Mass-Radius Relation (Panel c) shows how the "stiffness" of the EOS directly determines the load-bearing capacity of the neutron star. As shown, a stiffer EOS (larger \(G_V\)) can support a larger maximum mass. As \(G_V/G_S\) increases from 0.40 to 0.50, the maximum mass grows from \(2.16 M_\odot\) to \(2.25 M_\odot\). The M-R curves for all three parameter sets successfully pass through the observational constraint regions of the three key pulsars, indicating good observational compatibility of the model.

• Tidal Deformability vs. Mass (Panel d) shows that the stiffening of the EOS also leads to an increase in the tidal deformability parameter. At \(1.4 M_\odot\), as \(G_V/G_S\) increases, the \(\Lambda_{1.4}\) value rises from 348.83 to 385.32. However, the most critical finding is that even for the set with the largest \(G_V\) and the stiffest EOS, the \(\Lambda_{1.4}\) value remains far below the upper limit of \(\Lambda_{1.4} \lesssim 800\) set by GW170817.

In summary, the vector coupling constant \(G_V\) is a key parameter for controlling the properties of the hybrid star model, especially its maximum mass. By fine-tuning \(G_V\), we can precisely match the latest astronomical observations for maximum mass while ensuring the model satisfies constraints on radius and tidal deformability.

6 Impact of the Phase Transition Endpoint \(BU\) on Neutron Star Macroscopic Properties↩︎

Besides the vector coupling constant \(G_V\), which determines the "stiffness" of quark matter, the extent of the hadron-quark mixed phase is also a key factor affecting the macroscopic properties of neutron stars. In this section, we fix all NJL model parameters, including \(G_V\), to their benchmark values and vary only the phase transition endpoint coefficient \(BU\). The parameter \(BU\) defines the baryon number density at which matter fully transitions to the quark phase (\(n_B = BU \times n_0\)), so changing \(BU\) is equivalent to altering the width of the phase transition region. We construct three sets of equations of state (Set 1, Set 2, and Set 3) to systematically investigate the impact of the transition width on neutron star structure.

6.1 Parameter Sets and Comparison of Key Properties↩︎

We construct three sets of equations of state with \(BU = 4.5, 5.5,\) and \(6.5\) to systematically investigate the impact of the transition width on neutron star structure. All parameter sets use the same phase transition starting point, \(BL=1.0\). Table 3 quantitatively summarizes the key macroscopic properties of neutron stars for these three parameter sets.

The table shows that changes in the phase transition endpoint \(BU\) have a relatively small impact on the maximum mass, which remains around \(2.20 M_\odot\). However, this parameter significantly influences the radius and tidal deformability, revealing a notable trend. As the value of \(BU\) increases (i.e., the transition region widens):

• The radius (\(R_{1.4}\)) exhibits a clear and systematic increase, rising from \(12.21\) km to \(12.35\) km. This indicates that a wider transition window leads to a less compact star for intermediate masses.

• Similarly, the tidal deformability parameter (\(\Lambda_{1.4}\)) shows a systematic increase, rising from \(349.06\) to \(378.73\).

This positive correlation between the radius and tidal deformability for the BU parameter is a key finding, highlighting the complex influence of the crossover EOS on stellar structure.

6.2 Graphical Analysis of the EOS and Macroscopic Properties↩︎

To better visualize the impact of \(BU\), Figure 4 provides a side-by-side comparison of the three parameter sets at both the equation of state level and the macroscopic properties level.

From Figures 4 (a) and (b), we can clearly uncover the mechanism by which \(BU\) affects the "softness" of the EOS. Increasing the value of \(BU\) means widening the interval for the quintic polynomial interpolation. To smoothly connect the fixed starting point of the hadronic phase with the fixed form of the quark phase, a wider interpolation interval results in a lower and broader sound speed peak in the transition region. This is shown in panel (a), where the curve with the largest \(BU\) (Set 3, \(BU=6.5\)) has the lowest sound speed peak. Crucially, the narrowest transition window (Set 1, \(BU=4.5\)) leads to a very sharp peak that significantly exceeds the causality limit (\(v_s^2/c^2 > 1\)), rendering this parameter set physically invalid. The other two sets remain causal, thus establishing a lower bound for the width of the transition region.

The effect of a wider transition region on the star’s structure is clearly visible in Figure 4 (c). The curve with a larger \(BU\) value corresponds to a larger radius for a given intermediate mass. For instance, at \(1.4 M_\odot\), the radius for the \(BU=6.5\) set (approximately 12.35 km) is noticeably larger than for the \(BU=4.5\) set (approximately 12.21 km). It is important to note that despite this significant change in radius, the maximum masses of all three curves are very similar and all satisfy observational constraints.

The impact of \(BU\) on tidal deformability is shown in Figure 4 (d). Consistent with the trend in radius, a wider transition region (larger \(BU\)) leads to a larger tidal deformability parameter. The curve for \(BU=6.5\) is the highest in the intermediate mass range, indicating the largest tidal response, while the \(BU=4.5\) curve is the lowest. Nevertheless, all calculated \(\Lambda_{1.4}\) values for these parameter sets remain well below the upper limit from GW170817.

In summary, the phase transition endpoint parameter \(BU\) is an effective tool for tuning the radius and tidal deformability of neutron stars. It primarily adjusts the properties in the intermediate-density range without significantly affecting the maximum mass. The observation that a wider transition (larger BU) produces both a larger radius and a larger tidal deformability highlights the complex influence of the crossover EOS on stellar structure.

7 Conclusion↩︎

In this paper, we have constructed and systematically studied an Equation of State (EOS) for hadron-quark hybrid stars. The model describes hadronic matter using the BSR6 parameter set of relativistic mean-field theory and high-density quark matter using a two-flavor Nambu-Jona-Lasinio (NJL) model that includes both scalar and vector interactions. A core element of our work was the use of a quintic polynomial interpolation method ensuring \(C^2\) continuity to build a smooth transition region from the hadronic to the quark phase. This approach not only ensures the thermodynamic self-consistency and causality of the EOS over the entire density range but also provides a flexible framework for exploring the physical properties of the phase transition region.

Our primary goal was to test whether this hybrid star model could resolve a central tension in multi-messenger astronomy: the existence of massive pulsars requires the EOS to be sufficiently "stiff" at high densities, while constraints on tidal deformability from gravitational wave events require the EOS to be relatively "soft" at intermediate densities. Through a systematic scan and analysis of the model’s key parameters—the quark matter vector coupling constant \(G_V\) and the phase transition endpoint coefficient \(BU\)—we have reached the following core conclusions:

1. We successfully constructed a hybrid EOS that simultaneously satisfies all key current astronomical observational constraints. Taking the benchmark parameter set (Set 1) as an example, its predicted maximum neutron star mass is \(M_{max} = 2.20 M_\odot\), which fully satisfies the observational lower limit from massive pulsars like PSR J0740+6620 [6], [7]. At the same time, its predicted radius for a \(1.4 M_\odot\) neutron star is \(R_{1.4} = 12.27\) km, which is highly compatible with NICER satellite observations of PSR J0030+0451 and PSR J0437-4715 [28], [29]. More importantly, the model’s tidal deformability parameter is \(\Lambda_{1.4} = 352.97\), significantly below the strict upper limit of \(\Lambda_{1.4} \lesssim 800\) from the gravitational wave event GW170817 [13]. This result strongly demonstrates that a hybrid star with a hadron-quark phase transition is a self-consistent and highly promising physical picture for explaining neutron star observations.

2. We revealed that the vector coupling constant \(G_V\) is the key physical mechanism for controlling the maximum mass of neutron stars. The analysis in Section 6 shows that the vector repulsion provided by \(G_V\) is the main source of stiffening in the high-density EOS. As shown in Table 2, As shown in Table 2, increasing the \(G_V/G_S\) ratio from 0.40 to 0.50 directly raises the maximum mass from \(2.16 M_\odot\) to \(2.25 M_\odot\). This strong positive correlation indicates that \(G_V\) is the core parameter ensuring the model can support massive neutron stars. Our analysis also shows that the explored range of \(G_V\) is fully consistent with the causality condition.

3. We found that the phase transition endpoint \(BU\) is a critical parameter for adjusting the properties of intermediate-mass stars and is constrained by causality. The analysis in Section 6 revealed that widening the transition region (increasing BU from 4.5 to 6.5) makes the EOS softer, which leads to both a larger radius (\(R_{1.4}\) increases) and a larger tidal deformability (\(\Lambda_{1.4}\) increases). Most importantly, a transition region that is too narrow (e.g., \(BU=4.5\)) results in a violation of causality, establishing a physical lower bound on the extent of the crossover. This demonstrates how causality considerations can effectively constrain the nature of the hadron-quark transition.

4. A key finding of this study is the differential control mechanism that the model’s internal parameters exert on macroscopic properties. \(G_V\) primarily controls the EOS stiffness at high densities, thus determining the maximum mass; while \(BU\) primarily controls the EOS softness at intermediate densities—thus tuning the radius and tidal deformability—and is itself constrained by the causality requirement. This approximate decoupling of controls is a core strength of this hybrid star model. It means we are no longer in a situation of "robbing Peter to pay Paul" but can synergistically adjust different physical parameters to separately satisfy the seemingly contradictory constraints from different types of astronomical observations (massive pulsars vs. gravitational wave merger events).

8 Solving for the Interpolation Coefficients of the Hybrid EOS↩︎

In Section 3 of this thesis, a quintic polynomial is employed to construct a smooth crossover region from the hadronic phase to the quark phase. This polynomial describes the relationship between pressure \(P\) and baryon chemical potential \(\mu_B\) within the transition interval [\(\mu_{BL}, \mu_{BU}\)]: \[P(\mu_B) = \sum_{m=0}^{5} C_{m} \mu_{B}^{m} = C_0 + C_1 \mu_B + C_2 \mu_B^2 + C_3 \mu_B^3 + C_4 \mu_B^4 + C_5 \mu_B^5\] To ensure that the entire hybrid equation of state is \(C^2\) continuous—meaning that the pressure, baryon number density (\(\rho_B = dP/d\mu_B\)), and its first derivative (\(d\rho_B/d\mu_B = d^2P/d\mu_B^2\)) are all continuous at the boundaries—we impose six boundary conditions.

These six boundary conditions form a system of linear equations for the six unknown coefficients (\(C_0, \dots, C_5\)), specified as follows:

1. \(P(\mu_{BL}) = P_H(\mu_{BL})\)

2. \(P'(\mu_{BL}) = \rho_{B,H}(\mu_{BL})\)

3. \(P''(\mu_{BL}) = \chi_{B,H}(\mu_{BL})\)

4. \(P(\mu_{BU}) = P_Q(\mu_{BU})\)

5. \(P'(\mu_{BU}) = \rho_{B,Q}(\mu_{BU})\)

6. \(P''(\mu_{BU}) = \chi_{B,Q}(\mu_{BU})\)

Here, the subscripts \(H\) and \(Q\) denote quantities from the hadronic and quark phases, respectively, and \(\chi_B = d\rho_B/d\mu_B\) is the baryon number susceptibility.

This system of equations can be written in matrix form as \(\mathbf{A}\mathbf{x} = \mathbf{b}\), where \(\mathbf{x} = [C_0, C_1, C_2, C_3, C_4, C_5]^T\) is the vector of coefficients to be solved for. The matrix \(\mathbf{A}\) and vector \(\mathbf{b}\) are given by: \[\mathbf{A} = \begin{pmatrix} 1 & \mu_{BL} & \mu_{BL}^2 & \mu_{BL}^3 & \mu_{BL}^4 & \mu_{BL}^5 \\ 0 & 1 & 2\mu_{BL} & 3\mu_{BL}^2 & 4\mu_{BL}^3 & 5\mu_{BL}^4 \\ 0 & 0 & 2 & 6\mu_{BL} & 12\mu_{BL}^2 & 20\mu_{BL}^3 \\ 1 & \mu_{BU} & \mu_{BU}^2 & \mu_{BU}^3 & \mu_{BU}^4 & \mu_{BU}^5 \\ 0 & 1 & 2\mu_{BU} & 3\mu_{BU}^2 & 4\mu_{BU}^3 & 5\mu_{BU}^4 \\ 0 & 0 & 2 & 6\mu_{BU} & 12\mu_{BU}^2 & 20\mu_{BU}^3 \end{pmatrix}\] \[\mathbf{b} = \begin{pmatrix} P_H(\mu_{BL}) \\ \rho_{B,H}(\mu_{BL}) \\ \chi_{B,H}(\mu_{BL}) \\ P_Q(\mu_{BU}) \\ \rho_{B,Q}(\mu_{BU}) \\ \chi_{B,Q}(\mu_{BU}) \end{pmatrix}\] By solving this linear system, the polynomial coefficients \(\mathbf{x} = \mathbf{A}^{-1}\mathbf{b}\) can be uniquely determined, thereby constructing a thermodynamically consistent and smooth crossover equation of state.

For the baseline parameter set (Set 1) used in this work, the chemical potential boundaries for the transition region were found through optimization to be \(\mu_{BL} \approx 980.12 \text{ MeV}\) and \(\mu_{BU} \approx 1612.56 \text{ MeV}\). Solving the system of equations above yields the following approximate coefficient values:

• \(C_0 \approx 2.088 \times 10^{10}\)

• \(C_1 \approx -5.841 \times 10^7\)

• \(C_2 \approx 5.256 \times 10^4\)

• \(C_3 \approx -12.975\)

• \(C_4 \approx -3.720 \times 10^{-3}\)

• \(C_5 \approx 1.728 \times 10^{-6}\)

These coefficients ultimately define the crucial part of the hybrid equation of state used in this paper.

9 Explicit Forms of the Tidal Perturbation Functions↩︎

In Section 4.3 of this thesis, the master function \(H(r)\), which describes the quadrupole deformation of a neutron star under an external tidal field, is obtained by solving a second-order ordinary differential equation (Eq. 22 ). This equation arises from perturbing the background spacetime of a spherically symmetric star within the framework of general relativity.

Specifically, this result is derived by linearizing the Einstein field equations under the Regge-Wheeler gauge, considering static, even-parity, quadrupole (\(l=2\)) perturbations. The coefficients \(F(r)\) and \(Q(r)\) in the equation are entirely determined by the unperturbed background spacetime (the solution to the TOV equations) and the equation of state (EOS) of the matter. Their explicit expressions are as follows [30], [37]:

First, the background spacetime is described by the spherically symmetric Schwarzschild metric, which takes the form: \[ds^2 = -e^{2\nu(r)}dt^2 + e^{2\lambda(r)}dr^2 + r^2(d\theta^2 + \sin^2\theta d\phi^2)\] where the metric function \(e^{2\lambda(r)} = (1-2GM(r)/r)^{-1}\). The terms \(M(r)\) and \(P(r)\) are the radial mass and pressure profiles, respectively, obtained by integrating the TOV equations (Eqs. 19 and 20 ). The derivative of the other metric function, \(\nu(r)\), is given by: \[\nu'(r) = \frac{dP/dr}{\epsilon(r) + P(r)} = \frac{G(M(r) + 4\pi r^3 P(r))}{r(r - 2GM(r))}\] Here, \(\epsilon(r)\) is the energy density corresponding to \(P(r)\).

Using these background quantities, the coefficients \(F(r)\) and \(Q(r)\) from Eq. 22 can be explicitly expressed as: \[\begin{align} F(r) &= \left(1-\frac{2GM(r)}{r}\right)^{-1} \left[1 - 4\pi G r^2 (\epsilon(r) - P(r))\right] \\ Q(r) &= \left(1-\frac{2GM(r)}{r}\right)^{-1} \left[ 4\pi G \left(5\epsilon(r) + 9P(r) + \frac{\epsilon(r) + P(r)}{v_s^2(r)}\right) - \frac{6}{r^2} - (\nu'(r))^2 \right] \end{align}\] where \(v_s^2 = dP/d\epsilon\) is the local speed of sound squared, which directly reflects the stiffness of the equation of state. These expressions connect the microscopic properties of the EOS (via \(\epsilon, P, v_s^2\)) with the macroscopic gravitational structure of the star (via \(M(r), \nu(r)\)), and are crucial for calculating the tidal deformability.

References↩︎

---

1. Email: zhuchunran@outlook.com↩︎

2. Corresponding author email: blli@usst.edu.cn↩︎