Phase transition and stiffer core fluid in neutron stars: Effects on stellar configurations, dynamical stability, and tidal deformability

TL;DR

This paper explores how phase transitions and stiffer core fluids in neutron stars affect their structure, stability, and tidal deformability, with implications for interpreting gravitational wave observations like GW170817.

Contribution

It introduces a model incorporating phase transitions and stiffer fluids in neutron star cores, analyzing their effects on observable properties and stability.

Findings

Phase transition influences neutron star mass and radius.

Stiffer core fluid affects tidal deformability and stability.

Results align with observational data from GW170817.

Abstract

In this work, we investigate the influence of the phase transition and a stiffer fluid in neutron stars' cores on the static equilibrium configuration, dynamical stability, and tidal deformability. For this aim, it is taken into account that the fluid in the core and the envelope follow the relativistic polytropic equation of state. We find that the phase transition and a stiffer fluid in the core will reflect in the total mass, radius, speed of sound, core radius, radial stability with a slow and rapid conversion at the interface, and tidal deformability. We also investigate the dimensionless tidal deformability $\Lambda_1$ and $\Lambda_2$ for a binary neutron stars system with chirp mass equal to GW$170817$. Finally, we contrast our results with observational data to show the role that phase transition and a stiffer core fluid could play in the study of neutron stars.

Tables

Table 1. Table 1 : The central energy density ρ c subscript 𝜌 𝑐 \rho_{c} , the total mass M / M ⊙ 𝑀 subscript 𝑀 direct-product M/M_{\odot} , the total radius R 𝑅 R , and the core radius R core subscript 𝑅 core R_{\rm core} where the null eigenfrequency of oscillation of the slow case ω 0 , s subscript 𝜔 0 𝑠 \omega_{0,s} and rapid case ω 0 , r subscript 𝜔 0 𝑟 \omega_{0,r} are derived. These parameter are found for ρ inn dis = 8 × 10 17 ​ [ kg / m 3 ] superscript subscript 𝜌 inn dis 8 superscript 10 17 delimited-[] kg superscript m 3 \rho_{\rm inn}^{\rm dis}=8\times 10^{17}\,[\rm kg/m^{3}] , Γ out = 2.4 subscript Γ out 2.4 \Gamma_{\rm out}=2.4 and different values of Γ inn subscript Γ inn \Gamma_{\rm inn} and λ 𝜆 \lambda .

		${\omega }_{0,s}=0$				${\omega }_{0,r}=0$
${\Gamma }_{\mathrm{inn}}$	$\lambda$	${\rho }_c\left[\mathrm{kg}/m^3\right]$	$M/M_{\odot }$	$R\left[\mathrm{km}\right]$	$R_{\mathrm{core}}\left[\mathrm{km}\right]$	${\rho }_c\left[\mathrm{kg}/m^3\right]$	$M/M_{\odot }$	$R\left[\mathrm{km}\right]$	$R_{\mathrm{core}}\left[\mathrm{km}\right]$
	$1.0$	$1.9776\times {10}^{18}$	$2.2391$	$12.254$	$12.254$	$1.9776\times {10}^{18}$	$2.2391$	$12.254$	$12.254$
	$0.9$	$2.5095\times {10}^{18}$	$2.0689$	$11.525$	$7.8405$	$2.1462\times {10}^{18}$	$2.0753$	$11.949$	$7.8218$
$2.4$	$0.8$	$3.0423\times {10}^{18}$	$1.8742$	$10.800$	$7.5842$	$2.5079\times {10}^{18}$	$1.8830$	$11.339$	$7.6380$
	$0.7$	$3.8103\times {10}^{18}$	$1.6619$	$9.8790$	$7.1888$	$3.1426\times {10}^{18}$	$1.6685$	$10.403$	$7.3066$
	$0.6$	$5.1074\times {10}^{18}$	$1.4405$	$8.7010$	$6.6230$	$4.3246\times {10}^{18}$	$1.4448$	$9.1068$	$6.7657$
	$0.5$	$7.0873\times {10}^{18}$	$1.2207$	$7.4132$	$5.9219$	$6.3371\times {10}^{18}$	$1.2227$	$7.6410$	$6.0283$
	$1.0$	$2.0514\times {10}^{18}$	$2.2966$	$12.050$	$12.050$	$2.0514\times {10}^{18}$	$2.2966$	$12.050$	$12.050$
	$0.9$	$2.5269\times {10}^{18}$	$2.1368$	$11.375$	$8.0505$	$2.2691\times {10}^{18}$	$2.1410$	$11.661$	$8.0615$
$2.6$	$0.8$	$3.0803\times {10}^{18}$	$1.9542$	$10.623$	$7.7841$	$2.6228\times {10}^{18}$	$1.9612$	$11.054$	$7.8585$
	$0.7$	$3.8087\times {10}^{18}$	$1.7576$	$9.7431$	$7.3940$	$3.2424\times {10}^{18}$	$1.7634$	$10.155$	$7.5138$
	$0.6$	$4.7597\times {10}^{18}$	$1.5543$	$8.7586$	$6.8949$	$4.2734\times {10}^{18}$	$1.5568$	$9.0081$	$6.9936$
	$0.5$	$6.1448\times {10}^{18}$	$1.3472$	$7.6503$	$6.2731$	$5.8217\times {10}^{18}$	$1.3477$	$7.7576$	$6.3259$

Equations

$$T_{\mu \nu }=\left(\rho +p\right)u_{\mu }{u}_{\nu }+p{g}_{\mu \nu },$$

$$d{s}^2=-e^{\nu }d{t}^2+e^{\beta }d{r}^2+r^{2}d{\theta }^2+r^2{\sin }^2\theta d{\varphi }^2.$$

$$\frac{dm}{dr}=4\pi \rho r^2,$$

$$\frac{dp}{dr}=-(p+\rho )\left(4\pi rp+\frac{m}{r^2}\right)e^{\beta },$$

$$\frac{d\nu }{dr}=-\frac{2}{(p+\rho )}\frac{dp}{dr},$$

$$m(0)=0,p(0)=p_c,\rho (0)={\rho }_c,\mathrm{and}\nu (0)={\nu }_c.$$

$$p(R)=0.$$

$$e^{\nu (R)}=e^{-\beta (R)}=1-\frac{2M}{R},$$

$$\frac{d\xi }{dr}=\frac{\xi }{2}\frac{d\nu }{dr}-\frac{1}{r}\left(3\xi +\frac{\Delta p}{p\Gamma }\right),$$

$$\frac{d\Delta p}{dr}=(p+\rho ){\omega }^2\xi r{e}^{\beta -\nu }-4\xi \left(\frac{dp}{dr}\right)+{\left(\frac{dp}{dr}\right)}^2\frac{\xi r}{p+\rho }$$

$$-8\pi p(p+\rho )\xi r{e}^{\beta }-\left(\frac{1}{2}\frac{d\nu }{dr}+4\pi r{e}^{\beta }(p+\rho )\right)\Delta p,$$

$${\left(\Delta p\right)}_{\mathrm{center}}=-3{\left(\xi \Gamma p\right)}_{\mathrm{center}}.$$

$${\left(\Delta p\right)}_{\mathrm{surface}}=0.$$

$${\lambda }_1=-\frac{Q_{ij}}{{ϵ}_{ij}},$$

$$\Lambda =\frac{2}{3}\frac{k_2}{C^5},$$

$$k_2=\frac{8C^5}{5}(1-2C{)}^2[2+C(y_R-1)-y_R]$$

$$\times \{2C[6-3y_R+3C(5y_R-8)]+4C^3[13-11y_R$$

$$+C(3y_R-2)+2C^2(1+y_R)]+3(1-2C^2)$$

$$\times [2-y_R+2(y_R-1)]\ln (1-2C){\}}^{-1}.$$

$$r\frac{dy}{dr}+y^2+yF+r^{2}Q=0,$$

$$F=(\rho -p)\frac{r-4\pi r^3}{r-2m},$$

$$Q=4\pi e^{\beta }\left(5\rho +9p+\frac{\rho +p}{dp/d\rho }\right)-\frac{6e^{\beta }}{r^2}-{\left(\frac{d\nu }{dr}\right)}^2.$$

$${\xi }_{\mathrm{inn}}={\xi }_{\mathrm{out}}.$$

$${\left(\Delta p\right)}_{\mathrm{inn}}={\left(\Delta p\right)}_{\mathrm{out}}.$$

$${\left(\xi -\frac{\Delta p}{r{p}^{\prime }}\right)}_{\mathrm{inn}}={\left(\xi -\frac{\Delta p}{r{p}^{\prime }}\right)}_{\mathrm{out}},$$

$${\left(\Delta p\right)}_{\mathrm{inn}}={\left(\Delta p\right)}_{\mathrm{out}}.$$

$$y(r_{tr}+ϵ)=y(r_{tr}-ϵ)-\frac{4\pi r_{tr}^3\left[\rho (r_{tr}+ϵ)-\rho (r_{tr}-ϵ)\right]}{m(r_{tr})+4\pi r_{tr}^{3}p},$$

$$\rho$$

$$K_{\mathrm{inn}}=0.0195\times (1.67\times 10^{17}kg/{\mathrm{m}}^3{)}^{-1.34}.$$

$$K_{\mathrm{out}}=K_{\mathrm{inn}}\frac{{\left({\rho }_{\mathrm{inn}}^{\mathrm{dis}}-p_{\mathrm{inn}}^{\mathrm{dis}}/\left({\Gamma }_{\mathrm{inn}}-1\right)\right)}^{{\Gamma }_{\mathrm{inn}}}}{{\left({\rho }_{\mathrm{out}}^{\mathrm{dis}}-p_{\mathrm{out}}^{\mathrm{dis}}/\left({\Gamma }_{\mathrm{out}}-1\right)\right)}^{{\Gamma }_{\mathrm{out}}}}.$$

$${\rho }_{\mathrm{out}}^{\mathrm{dis}}=\lambda {\rho }_{\mathrm{inn}}^{\mathrm{dis}},$$

$$\mathcal{M}=\frac{(m_1{m}_2{)}^{3/5}}{(m_1+m_2{)}^{1/5}}.$$

Full text

††journal: Eur. Phys. J. C

1 Introduction

The direct multimessenger detection from binary black holes merger carried out by the LIGO-Virgo scientific network abbott2016 ; abbott2016_2 ; abbott2017 ; abbott2017_1 ; abbott2017_2 has marked the starting of the era of Gravitational Waves (GWs) astronomy. The detection of GWs has opened a new window to explore the cosmos and supplied some astrophysics and fundamental physics implications (check, e.g., abbott2017_3 ; abbott2016_3 ; yunes2016 ). Another important event of GWs comes from a merger of a pair of neutron stars (NSs) abbott2017_NS , known as event GW$170817$, which was also reported by the LIGO-Virgo scientific network. This new signal opened the GWs multi-messenger astronomy, being the first detection with electromagnetic counterpart abbott2017_4 , and providing a set of valuable information about the properties of NSs and their equation of state (EOS). After the first detection of GWs from the NSs binary system, many important efforts have been realized to constrain, e.g., their radii and EOS malik_2018 ; most_2018 ; annala_2018 ; paschalidis_2018 ; lugones2019 . Additional constraints are feasible because of the implication of tidal deformation hinderer_2008 ; damour_2009 ; hinderer_2010 .

It is known that the matter density that makes up NSs reach densities up to a few times the nuclear saturation density, however, until these days, detailed information about the characteristics and nature of their deep interiors is still lacking. Future multi-messenger signatures hold the promise of identifying the specific internal aspect of NSs. Theoretically, asteroseismology is widely employed to analyze the internal structure of compact stars -the name used for white dwarfs, neutron stars, hybrid stars, or strange quark stars- to investigate the thermodynamic properties inside these objects. Through this diagnostic technique, analyzing the frequency modes can obtain a solid way to learn more about the physics inside compact stars. For example, if inside these stars a single component fluid is present chanmugam1977 ; Alcock_1986 ; vath_chanmugam1992 ; bombaci1996 ; gondek1997 ; gondek1999 ; kokkotas2001 ; lugones2010 ; doneva2012 ; arbanil_malheiro2015 ; arbanil_malheiro_2016 ; flores2017 ; sagun_2020 ; annala_2022 or the existence of a phase transition between layers with different mechanical properties pereira_flores2018 ; tonetto_2020 ; miniutti_2003 ; mishustin_2003 ; sahu2002 ; gupta2001 ; flores_lenzi2012 ; brillante_2014 ; clemente_2020 ; jimenez_2021 ; parisi_2021 ; sun_2021 ; mariani_lugones2019 ; ivanytskyi_2022 .

In literature, unlike the study of one-phase static compact stars, two-phase stars and the impact of the phase transition on the properties of these stars were not widely investigated. There are studies analyzing how density jumps affect the static stellar equilibrium configuration and radial frequency of oscillations sahu2002 ; gupta2001 ; flores_lenzi2012 ; brillante_2014 ; clemente_2020 ; jimenez_2021 ; sun_2021 ; parisi_2021 ; mariani_lugones2019 ; pereira_flores2018 ; tonetto_2020 ; mishustin_2003 , as well as the possibility of arising of the so-called gravitational pulsation mode ($g$-mode) miniutti_2003 ; mishustin_2003 ; tonetto_2020 .

As regards the radial perturbations of compact stars with a sharp interface, the set of equations must be solved by taking into account the additional boundary conditions at the phase-splitting interface. Around this point, there are two types of physical behavior due to radial perturbations: the slow and rapid phase conversion pereira_flores2018 ; tonetto_2020 . In the case of slow conversion, there is no change of matter over the pulsating interface. On the contrary, the rapid conversion case involves a flow of mass from one phase to the other, and vice-versa, through the moving phase boundary. In recent years, the impact of the phase transition on the radial oscillations of compact stars has been reported in different articles. In the case of slow transition, for example, the authors concentrate on investigating the effects of the core formation mishustin_2003 , a sharp phase transition clemente_2020 ; parisi_2021 , the mixed-phase sahu2002 ; gupta2001 ; flores_lenzi2012 ; sun_2021 ; jimenez_2021 , and electric charge brillante_2014 . On the other hand, among those reported considering both phase transitions in compact stars, we find: the ones that analyze the fast and slow conversion in the context of general relativity pereira_flores2018 , how these are affected in the face of a magnetic field mariani_lugones2019 , and their influences on non-radial oscillations tonetto_2020 . In the rapid phase transition case, unlike the slow phase transition, in a sequence of equilibrium configuration, the maximum mass peak marks the beginning of radial instability. Indeed, in slow transitions, after this turning point, it is possible to find additional stable equilibrium configurations. Therefore, in a sequence of equilibrium configurations with increasing central energy density, some stars with the same mass but different radii are obtained. These stars are known as twin stars.

In the aforementioned articles, different models of equations of state are studied from the perspective of observational deformability data from the event GW170917. Some of these works study this phenomenon against the possibility of the existence of phase transitions inside the compact star, some taking the aspect of an analysis of the stability of these stars, and calculating radial oscillations. However, these same articles often make these calculations assuming only slow transitions, which allow the appearance of stable regions, in the mass-radius diagram, after the maximum mass. In this work, we present a detailed study of the influence of the phase transition and a stiffer fluid in NSs core on the equilibrium configuration, radial stability, and tidal deformability. In this sense, we analyze how the radius, mass, speed of sound, core radius, radial frequency of oscillation, and tidal deformation change when a phase transition and stiffer fluid in NSs core are considered. In the analysis of the radial stability of NSs, we will focus on the slow and rapid phase conversions. We also contrast our results with observational data to see the role that phase transition and a stiffer core fluid could play in the study of NSs.

The present article is arranged as follows: Section 2 presents the equilibrium equations and radial stability equations; moreover, this section is also devoted to presenting the junction conditions at the interface of the two-phase, which are required to investigate the slow and rapid phase conversions. Section 3 presents the EOSs employed for NSs, as well as the numerical method used to solve the complete set of equations required to investigate the equilibrium and radial stability. In Section 4 we show the numerical results for equilibrium configurations, radial stability, and tidal deformability of NSs with two-phase. Finally, we conclude in Section 5. Throughout the paper, we work with geometric units, i.e., $c=1=G$, and the metric signature $+2$.

2 General relativistic formulations

2.1 Equilibrium equations

We take into account that the unperturbed neutron star is made up of layers of effective perfect fluids, whose energy-momentum tensors can be expressed as

[TABLE]

with $\rho$, $p$, and $u_{\mu}$ representing respectively the energy density, the fluid pressure, and the four-velocity.

To analyze the effect of phase transition in the dense matter on the equilibrium and radial stability of neutron stars, we set the space-time metric, in Schwarzschild coordinates, as

[TABLE]

The potential metric functions $\nu=\nu(r)$ and $\beta=\beta(r)$ depend on the radial coordinate $r$ only.

For the energy-momentum tensor (Eq. (1)) and line element (Eq. (2)) adopted, with the potential metric $e^{-\beta}=\left(1-2m/r\right)$, we derive the set of stellar structure equations

[TABLE]

where the parameter $m$ represents the mass inside the sphere of radius $r$. Eq. (4) is known as the hydrostatic equilibrium equation for a spherically symmetric static astrophysical object, also called as Tolman-Oppenheimer-Volkoff equation tolman ; oppievolkoff .

The stellar structure equations (3)-(5) are integrated from the center toward the star’s surface. At the center $(r=0)$ the integration starts

[TABLE]

The surface of the star $(r=R)$ is determined by

[TABLE]

At this point, the interior solution connects smoothly with the exterior Schwarzschild vacuum solution. This indicates that at the star’s surface the interior and exterior potential metrics are related through of the form:

[TABLE]

with $M$ being the total mass of the star.

2.2 Radial oscillations equations

The radial pulsation equation is obtained by Chandrasekhar chandrasekhar_rp making perturbation in the fluid and space-time variables. The perturbed quantities are placed into Einstein’s field equation and in the linearized form of the conservation of stress-energy tensor.

The solution of the radial oscillation equation provides information about the eigenfrequency of oscillations $\omega$. Intending to set this equation in a more appropriate form for numerical integration, we place it into two first-order equations for the variables $\Delta r/r$ and $\Delta p$ gondek1999 ; with $\Delta r$ and $\Delta p$ representing respectively the relative radial displacement and Lagrangian perturbations of pressure. Thus, the system of equations is, $\xi=\Delta r/r$:

[TABLE]

where $\Gamma=\left(1+\frac{\rho}{p}\right)\frac{dp}{d\rho}$. The variables $\xi$ and $\Delta p$ have a time dependence of the form $e^{i\omega t}$, with $\omega$ being the eigenfrequency.

To solve the differential equations (9) and (10), boundary conditions in the center and on the star’s surface are required. Moreover, to find regular solutions in the center of the star, the second term of the right-hand side of Eq. (9) must vanish in $r\rightarrow 0$. In this way, it is considered

[TABLE]

At this point, for normalized eigenfunctions, we regard $\xi(r=0)=1$. On the other hand, as established above, the surface of the star is determined when $p(R)=0$. It implies

[TABLE]

2.3 Tidal deformability

Tidal effects are very common in the context of NSs binary systems. In fact, the gravitational field generated by one star in a binary system can result in deformation in its companion. The parameter of tidal deformability is the measure of the deformation in compact stars due to an external field. From a mathematical point of view, this parameter can be expressed in terms of the fraction,

[TABLE]

where $Q_{ij}$ is the quadrupole moment perturbed by an external tidal field $\epsilon_{ij}$ hinderer_2008 ; damour_2009 ; hinderer_2010 . The tidal deformability parameter $\lambda_{1}$ is connected with the Love number $k_{2}$ through the relation $k_{2}=\frac{3}{2}\lambda_{1}R^{-5}$. Moreover, the dimensionless tidal deformability $\Lambda$ can be written in terms of the Love number $k_{2}$ as

[TABLE]

with $C=M/R$ being the compactness parameter. $k_{2}$ can be expressed in terms of parameter $y_{R}\equiv y(r=R)$ as follows

[TABLE]

The parameter $y(r)$ is calculated along the whole of the star -from the center to the surface of the star- integrating the equation

[TABLE]

together with the set of equations (3)-(5), considering at the center of the star $y(r=0)=2$. The functions $F=F(r)$ and $Q=Q(r)$ are represented by the relations:

[TABLE]

2.4 Junction conditions at the interface

In the last years, compact stars with two different phases have been considered a real possibility. However, there are still a bunch of open questions such as, for instance, the density at which a hadron-quark phase transition occurs and some discussions about the kind of phase transition depending on the surface tension between the phases pereira_flores2018 ; tonetto_2020 ; Ilda1998 ; mariani_lugones2019 ; T-endo2011 . In our case, a first-order phase transition is considered, which results in the presence of a finite energy density discontinuity. Moreover, some approaches are regarded to investigate the radial stability and deformability of stars.

2.4.1 Radial oscillations

The phase transitions can be classified as slow or rapid depending of the time scale of the reaction of the matter in the neighborhood of the hadron-quark interface haensel1989 ; pereira_flores2018 . A scenario of slow phase transition appears when the rate of reaction transforming one phase into another is much greater than those of the radial perturbations. In such circumstances, there is no flow of matter across the surface splitting the two phases. Such a condition implies that $\xi$ must always be continuous at the interface across the interface, i.e.,

[TABLE]

Moreover, this also leads to a continuity pressure at the interface of the two phases. At this point, it assures that

[TABLE]

In the case of a rapid phase transition, the rate of reaction transforming one phase into another is much lower than those of the radial perturbations. In this scenario, a change of mass flow through the interface occurs. The diminution of mass on one side should be equal to the increase of mass on the other side. This condition, together with the demand for the continuity of pressure, lead to

[TABLE]

where the prime stands the operation with respect to the radial derivative, and

[TABLE]

2.4.2 Tidal deformability

At the interface of the two layers, we can clearly note that exists a singularity in Eq. (18) due to the speed of sound ($dp/d\rho$) since, at this point, we have the same value of the fluid pressure for two different energy densities. In Ref. damour_2009 , the authors discussed this problem in the context of the surface vacuum discontinuity for incompressible stars. After in Refs postnikov_2010 ; takatsy_2020 ; zhang_2002 , this approach was extended for the case of first-order transitions inside hybrid-stars, where authors concluded that at this point the function $y(r)$ must follow the next condition:

[TABLE]

with $r_{tr}$ and $\epsilon$ being the radial position where phase-transition occurs inside of the star and $\epsilon$ an infinitesimal parameter, respectively. Obviously, the region $r<r_{tr}$ represents the core of the star, and the region $r>r_{tr}$ depicts the envelope of the star.

3 Equation of state and numerical method

3.1 Equation of state

To describe the matter that makes up the compact object, in the two-phase configurations, the relativistic polytropic equation of state tooper1965 is adopted. Then, the energy density and fluid pressure of each phase are respectively connected through the relations:

[TABLE]

These two relations bear parameters from the inner and outer regions denoted by the sub-indexes “${\rm inn}$” and “${\rm out}$”, respectively; namely, $K_{\rm inn}$ and $K_{\rm out}$ are the polytropic constants, $\Gamma_{\rm inn}$ and $\Gamma_{\rm out}$ being the polytropic exponents, and $p^{\rm dis}_{\rm inn}$ and $p^{\rm dis}_{\rm out}$ represent the phase transition pressure. Following doneva2012 , we set the inner polytropic constant value:

[TABLE]

Since the fluid pressure should be continuous along the star, at the phase transition point, where the inner and outer phase transition pressures meet the condition $p_{\rm out}^{\rm dis}=p_{\rm inn}^{\rm dis}$, the outer polytropic constant takes the form:

[TABLE]

At this point, the inner and outer phase transition energy densities are related by:

[TABLE]

where $\lambda$ is known as the density jump parameter. Some examples of the EOS with a sharp density jump are presented in Fig. 1.

In the next section, we consider the value of $\Gamma_{\rm inn}=2.4$, $\Gamma_{\rm inn}=2.6$ and $\Gamma_{\rm out}=2.4$ and the phase transition parameter $\lambda$ between the values $0.5$ and $1.0$, in order to analyze both the role of a stiffer core fluid and the effects of a more abrupt phase transition in the star equilibrium configuration, respectively. The values of $\Gamma_{\rm inn},\Gamma_{\rm out}$ and $\lambda$ chosen allow us to obtain comparable results with the data found through observation, which are reported by the LIGO-Virgo network in abbott_2018a_tidal and by NICER in riley2019 ; miller2019 ; riley2021 ; miller2021 ; cromartie2020 ; antoniadis2013 ; demorest2010 .

3.2 Numerical method

The effects of the phase transition on the equilibrium configurations and tidal deformations are investigated through the numerical solution of the system of equations, boundary conditions, and junction conditions established in the section 2, for each $\rho_{\rm inn}^{\rm dis}$, $\lambda$, $\Gamma_{\rm inn}$, and $\Gamma_{\rm out}$. This system of equations is integrated from the center toward the star’s surface.

The analysis of the radial stability starts by solving the stellar structure equations by using the Runge-Kutta fourth-order method in order to determine the radial pulsation coefficients. After, we begin at the star’s core with the solution of Eqs. (9)-(10) for a test value of $\omega^{2}$. These equations are numerically integrated outwards until the interface is found, where the junction conditions are employed (for the slow case Eqs. (19) and (20) and for the rapid case Eqs. (21) and (22)) to find the values of $\xi$ and $\Delta p$ at the other side of the interface. Then, the numerical integration continues towards the surface of the stars attempting to reach the conditions (7) and (12). Whether, after each integration, the equality (12) is not fulfilled, $\omega^{2}$ is corrected until satisfying this equality in the next integration. The parameters $\omega^{2}$ which satisfy the oscillation are called eigenvalues of the radial pulsation equation and $\omega$ of eigenfrequencies (review pereira_flores2018 ).

We are interested in analyzing the radial stability of stars, we only analyze the lowest eigenvalue, i.e., $\omega_{0}^{2}$. When $\omega_{0}^{2}>0$, the star is stable against small radial perturbations. $\omega_{0}$ is known as the eigenfrequency of the fundamental mode.

4 Results

4.1 Equilibrium configuration of neutron stars

The compact star total mass sequence, normalized to the Sun’s mass $M_{\odot}$, as a function of the central energy density is plotted in Fig. 2 for two values of $\Gamma_{\rm inn}$, four different values of $\lambda$, $\rho_{\rm inn}^{\rm dis}=8\times 10^{17}\,[\rm kg/m^{3}]$, and $\Gamma_{\rm out}=2.4$. The central energy density goes from $10^{18}$ to $10^{20}\,[\rm kg/m^{3}]$. On the panel, the total mass increases with the central energy density until reaches the maximum mass of the sequence, after this point, $M/M_{\odot}$ decreases monotonically with the increment of $\rho_{c}$.

In the panel, it is noted the diminution of the total mass with the density jump parameter. This is associated with the fact that the fluid pressure decays abruptly due to the presence of a phase transition, being this declines greater for lower $\lambda$ (see Fig. 1). The change of the mass with $\Gamma_{\rm inn}$ is also observed in Fig. 2. For a fixed $\rho_{c}$ and $\lambda$, for a greater $\Gamma_{\rm inn}$, a larger total mass is derived. This point can be understood since for larger interior polytropic exponents larger central pressures are obtained. Thus, larger central pressure supports more mass against gravitational collapse. In addition, it is important to say that, for larger $\Gamma_{\rm inn}$, neutron stars with more compact cores are obtained (review, e.g., alz-poli-qbh ; alz-2eos-qbh ).

Fig. 3 shows the mass as a function of the total radius for two values of $\Gamma_{\rm inn}$, for $\rho_{\rm inn}^{\rm dis}=8\times 10^{17}[\rm kg/m^{3}]$, few values of $\lambda$, and $\Gamma_{\rm out}=2.4$. As in the case of Fig. 2, in this figure is considered $\rho_{c}$ between $10^{18}$ and $10^{20}\,[\rm kg/m^{3}]$. The panel shows an increment of $M/M_{\odot}$ with the diminution of $R$ until to find a $M_{\rm max}/M_{\odot}$. After this point, the curves turn anti-clockwise, to $M(R)$ starts to decrease with $R$ until to reach $R_{\rm min}$. From here on, the mass decays with the increment of the total radius.

In Fig. 3, for some range of central energy density, we also find an increment of the total radius with the diminution of the density jump parameter. This is due to the fact that in these compact stars, the pressure of the fluid decays slower with the increase of the radial coordinate, thus obtaining larger radii. In addition, for some range of $\rho_{c}$ and $\lambda$, we can also observe decrements of the total radius with the increment of $\Gamma_{\rm inn}$. Despite having an increase of the central pressure with $\Gamma_{\rm inn}$, it decays faster with the growth of the radial coordinate. In this way, compact objects have a smaller total radius. Finally, in Fig. 3, we also see that the radius of the stars with maximum mass decreases with the jump in phase transition density. This is due to the fact that the central pressure of the star decreases with $\lambda$, thus, the pressure decays faster with the radial coordinate.

The total mass versus the speed of sound is plotted in Fig. 4, for $\rho_{\rm inn}^{\rm dis}=8\times 10^{17}[\rm kg/m^{3}]$, some values of $\lambda$, and $\Gamma_{\rm out}=2.4$. In this figure, we only present equilibrium configurations with the speed of sound lower than the speed of light $c_{s}^{2}=1.0$. As can be seen, in particular, at stars that present first-order phase transition with low central energy densities (stars with low total masses), the speed of sound never exceeds the conformal limit $c_{s}^{2}=1/3$. However, in stars with larger total masses, the speed of sound exceeds the conformal limit value but is far to attain the speed of light. These results are in concordance with those reported in the article annala_2022 .

4.2 Radial stability of neutron stars

Fig. 5 shows the behavior of the slow and rapid eigenfrequency of the fundamental oscillations as a function of the total mass and against the central energy density for some different values of $\lambda$, $\rho_{\rm inn}^{\rm dis}=8\times 10^{17}\,[\rm kg/m^{3}]$, $\Gamma_{\rm out}=2.4$, and $\Gamma_{\rm inn}=2.4$ on the left panels, and $\Gamma_{\rm inn}=2.6$ on the right panels. From the figure, for $\lambda=1$, it can be noted that the maximum total mass is found at the zero eigenfrequencies of oscillation. This case represents the usual study of radial oscillation of neutron stars in the absence of phase transition. In turn, for $\lambda<1$, as well as in pereira_flores2018 , we note that the total mass at the null eigenfrequency of oscillation depends on the type of the phase transition. At this point, when $\Gamma_{\rm inn}=\Gamma_{\rm out}=2.4$, the difference of the mass attained in the rapid and slow case is almost $0.246\%$ for $\lambda=0.9$, $0.470\%$ for $\lambda=0.8$, $0.397\%$ for $\lambda=0.7$, $0.299\%$ for $\lambda=0.6$, and $0.164\%$ for $\lambda=0.5$. In the case $\Gamma_{\rm inn}=2.6$ and $\Gamma_{\rm out}=2.4$, the difference of the mass is around $0.197\%$ for $\lambda=0.9$, $0.358\%$ for $\lambda=0.8$, $0.330\%$ for $\lambda=0.7$, $0.161\%$ for $\lambda=0.6$, and $0.0371\%$ for $\lambda=0.5$ (review Table 1).

In Fig. 5, in the slow case, the total mass at the zero eigenfrequencies of oscillation is derived at $\rho_{c}$ larger than the one employed to obtain the maximum mass value; i.e., twins stars are derived, stable stars with the same total mass but with both different central energy densities and total radii. However, in the rapid case, the maximum mass and the null eigenfrequency of oscillation are obtained by using the same value of $\rho_{c}$. It indicates that in a sequence of static equilibrium configurations, the maximum mass point marks the beginning of the instability against small radial perturbations. This characteristic in each phase transition could be useful to differentiate them.

Table 1 presents the central energy densities and total masses where the zero eigenfrequencies of oscillations for the slow and the rapid phase transition. These parameters are derived for $\rho_{\rm inn}^{\rm dis}=8\times 10^{17}\,[\rm kg/m^{3}]$, $\Gamma_{\rm out}=2.4$ and some values of $\Gamma_{\rm inn}$ and $\lambda$. In the table, for a fixed $\Gamma_{\rm inn}$, at the null eigenfrequency of oscillation for the slow and rapid case, we note that the total mass, the total radius, and the core radius decrease with the density jump parameter. This could be understood since $p_{c}$ decays with $\lambda$, in this way, the total pressure diminishes faster with the growth of the radial coordinate. On the other hand, when $\Gamma_{\rm inn}$ is increased, i.e., when a stiffer core fluid is considered, stars with a core radius $R_{\rm core}$ closer to the total radius $R$ are found.

4.3 Tidal deformability in the light of GW$170817$

Tidal deformability against the total mass of stable NS is plotted at the top panel of Fig. 6 for $\rho_{\rm inn}^{\rm dis}=8\times 10^{17}\,[\rm kg/m^{3}]$, $\Gamma_{\rm out}=2.4$, different values of $\lambda$ and $\Gamma_{\rm inn}$. On the panel, it is also presented the tidal deformability constrained by the event GW$170817$ for a star of $1.4M_{\odot}$ to be $70\leq\Lambda_{1.4M_{\odot}}\leq 580$ abbott_2018a_tidal , at 90% confidence level, for low-spin priors.

For a fixed $\lambda$ and masses range, at the top of Fig. 6, we note an increment of the tidal deformability when a stiffer fluid is considered in the core. On the other hand, by setting the parameter $\Gamma_{\rm inn}$ for different values of the phase transition parameter $\lambda$, we can notice that the importance of the parameter $\Lambda$ changes in a more relevant way. The lower the value of $\lambda$, the less the pressure of transition (see Fig. 1) and the smaller the value of the deformability for the same mass. From these results, the effect of phase transition and stiffer fluid in the core are noticeable in the tidal deformability. Nonetheless, between these two factors, the phase transition parameter affects the NS properties more.

At the bottom panel of Fig. 6, the curves $\Lambda_{1}-\Lambda_{2}$ are presented for a binary NS system with chirp mass equal to GW$170817$. Since the total mass is associated with the dimensionless tidal deformability (top panel of Fig. 6), the curves $\Lambda_{1}-\Lambda_{2}$ are obtained once chosen a value of $m_{1}$ and calculating $m_{2}$ for the fixed value of the chirp mass ${\cal M}=1.188\,M_{\odot}$ abbott2017_NS , which is theoretically calculated by the relation:

[TABLE]

The values considered for $m_{1}$ and $m_{2}$ run from $1.36M_{\odot}\leq m_{1}\leq 1.60M_{\odot}$ and $1.17M_{\odot}\leq m_{2}\leq 1.36M_{\odot}$, respectively.

At the bottom panel of Fig. 6, we investigate the effects of the phase transition and a stiffer fluid in the core (for $\rho_{\rm inn}^{\rm dis}=8\times 10^{17}\,[\rm kg/m^{3}]$) of two neutron stars in a binary system. In this case, we note that both the phase transition and a stiffer equation of state could play an important role in the detection of these compact objects. From the results, we note that there is an interval for the density jump parameter $\lambda$, $0.7\leq\lambda\leq 0.8$, which plays the compact stars inside of $50\%$ and $90\%$ regions. The compact stars with $\lambda\geq 0.9$ are outside the $90\%$ region and do not appear in the panel. This is realized with the aim that the curves shown inside $50\%$ and $90\%$ can be seen clearly. Furthermore, we observe that for smaller values of $\lambda$ the curve change to smaller values of dimensionless deformability. On the other hand, it can be noted that a stiffer fluid in the core produces larger values of deformability.

4.4 Change of stellar physical parameters with phase transition energy density and its comparison with observational data

In Fig. 7, the mass-radius curves are compared with the observational data considering $\Gamma_{\rm inn}=\Gamma_{\rm out}=2.4$, some values of $\lambda$ and two $\rho_{\rm inn}^{\rm dis}$. On the top and bottom panel are employed $\rho_{\rm inn}^{\rm dis}=7.0\times 10^{17}[\rm kg/m^{3}]$ and $\rho_{\rm inn}^{\rm dis}=9.0\times 10^{17}[\rm kg/m^{3}]$, respectively. In this figure is used $10^{18}\leq\rho_{c}\leq 10^{20}[\rm kg/m^{3}]$. The observation data correspond to the NICER constraints obtained from the pulsars PSR J$0030+0451$ riley2019 ; miller2019 and PSR J$0740+6620$ riley2021 ; miller2021 . The corresponding bands of the pulsars PSR J$0740+6620$ cromartie2020 , PSR J$0348+0432$ antoniadis2013 and PSR J$1614+2230$ demorest2010 are also presented. From the figure, we observe that the change of $\rho_{\rm inn}^{\rm dis}$ affects the stellar structure configuration; being this change more noticeable in the range of low central energy densities. In this interval, for larger $\rho_{\rm inn}^{\rm dis}$, greater total mass and smaller total radius are found. From these results, we note that the change of internal phase transition energy density allows us to obtain some results more accurate and closer to empirical evidence of the neutron stars PSR J$0030+0451$. In addition, we see that with the increment of $\rho_{\rm inn}^{\rm dis}$, grow the possibility of having equilibrium solutions with sharper phase transition (smaller density jump parameter $\lambda$) within PSR J$0030+0451$. On the other hand, from Figs. 3 and 7, we observe that in the range of larger total mass, a stiffer fluid in the core could help to reach empirical evidence of neutron stars PSR J$0740+6620$, PSR J$0348+0432$, and PSR J$1614+2230$.

In Fig. 8, the top and bottom panels, respectively, present the tidal deformability against the total mass and $\Lambda_{1}$-$\Lambda_{2}$ curves for a binary NS system with chirp mass equal to GW$170817$ considering the relation (29) where $m_{1}$ and $m_{2}$ goes from $1.36M_{\odot}\leq m_{1}\leq 1.60M_{\odot}$ and $1.17M_{\odot}\leq m_{2}\leq 1.36M_{\odot}$. The inner phase transition energy density considered on the left and right panels are $7.0\times 10^{17}[\rm kg/m^{3}]$ and $9.0\times 10^{17}[\rm kg/m^{3}]$, respectively. In all panels of Fig. 8, only stable equilibrium configurations for the rapid transition case, employing $\Gamma_{\rm inn}=\Gamma_{\rm out}=2.4$, are considered. As stated above, the observation data correspond to the event GW$170817$.

In Fig. 8, on the top panels, it can be seen that all curves tidal deformability ($\Lambda$)-total mass ($M/M_{\odot}$) decay with the increment of $\rho_{\rm inn}^{\rm dis}$. From this, we understand that this phenomenon allows us to have equilibrium configurations with a lower density jump parameter $\lambda$ (a sharper phase transition) within the observational data of the GW$170817$ event. From the bottom panels, we note that some equilibrium solutions located outside the range of observational data derived from the GW$170817$ event fall within this interval when we increase the phase transition energy density $\rho_{\rm inn}^{\rm dis}$. In addition, from Figs. 8 and 6, we can also say that having a stiffer fluid in neutron stars’ cores helps us to have static equilibrium configurations within the range of observational data.

5 Conclusions

In this work, we investigated the influence of the phase transition on the equilibrium, radial stability, and tidal deformability of NSs with a stiffer fluid in the core. In the core and the envelope of the star, the relativistic polytropic equation of state is considered. The spherical equilibrium configurations are connected smoothly with the Schwarzschild exterior spacetime. We examined the change of the mass, radius, speed of sound, core radius, the eigenfrequency of the fundamental mode of the star with a slow and rapid phase conversion at the interface, and tidal deformability for different density jump parameters $\lambda$, phase transition energy densities $\rho_{\rm inn}^{\rm dis}$, interior polytropic exponents $\Gamma_{\rm inn}$, and exterior polytropic exponent $\Gamma_{\rm out}=2.4$.

As well as in the study of NSs developed in sotani_2001 , which employs a non-relativistic polytropic equation of state, we note that some aspects of the static equilibrium configurations -such as the mass and radius- are affected by the phase transition, stiffer fluid in the core (which change with $\Gamma_{\rm inn}$), and phase transition energy density.

For the values $\Gamma_{\rm inn}$ and $\lambda$ employed, in the slow case, the zero eigenfrequencies of the fundamental mode are attained beyond the maximum mass points and, in the rapid case, the maximum masses points mark the beginning of the radial instability thus indicating that the regions constituted by stable and unstable stars can be recognized by the conditions $dM/d\rho_{c}>0$ and $dM/d\rho_{c}<0$, respectively. These results are in concordance with those one reported in the works pereira_flores2018 ; tonetto_2020 ; mariani_lugones2019 .

The change of the tidal deformability for a NS ($\Lambda$) and a binary NS system ($\Lambda_{1}$ and $\Lambda_{2}$), with equal chirp mass as GW$170817$ event, as a function of $\Gamma_{\rm inn}$ and $\lambda$, has been analyzed. We obtained a dependence of the dimensionless tidal deformability with these two factors in the aforementioned frameworks. For NS configurations, for some interval of masses, we noted that $\Lambda$ grows and decreases with the increment of $\Gamma_{\rm inn}$ and diminution of $\lambda$. In turn, for a binary NS scenario, we showed that the phase transition and stiffer fluid in the NS core could also play an important role in the detection of NSs. These results are in agreement with those ones published in parisi_2021 .

We also investigated the dependence of some physical parameters of NSs with the phase transition energy density. At the interval of low central energy densities, we found that $\rho_{\rm inn}^{\rm dis}$ can also be important in the study of NSs since their physical parameters could be significantly affected by the value of phase transition energy density.

Finally, we noted that a change in jump density parameter $\lambda$, phase transition energy density $\rho_{\rm inn}^{\rm dis}$, and stiffer core fluid $\Gamma_{\rm inn}$ could lead to the possibility that some equations of state that are outside of the observational data, can be inside this framework for accurate values of $\lambda$, $\rho_{\rm inn}^{\rm dis}$ and $\Gamma_{\rm inn}$.

Acknowledgements.

JDVA would like to thank the Universidad Privada del Norte and Universidad Nacional Mayor de San Marcos for funding - RR Nº$\,005753$-$2021$-R$/$UNMSM under project number B$211\-31781$. CHL is thankful to the Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP) under thematic project $2017/05660-0$, Grant No. $2020/05238-9$ and Power Data Tecnologia Ltda for providing a technological environment for data processing.