1 Introduction

The direct multimessenger detection from binary black holes merger carried out by the LIGO-Virgo scientific network [1,2,3,4,5] 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., [6,7,8]). Another important event of GWs comes from a merger of a pair of neutron stars (NSs) [9], known as event GW170817, 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 [10], 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 [11,12,13,14,15]. Additional constraints are feasible because of the implication of tidal deformation [16,17,18].

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 [19,20,21,22,23,24,25,26,27,28,29,30,31,32] or the existence of a phase transition between layers with different mechanical properties [33,34,35,36,37,38,39,40,41,42,43,44,45,46].

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 [34,35,36,37,38,39,40,41,42,43,44,45], as well as the possibility of arising of the so-called gravitational pulsation mode (g-mode) [33,34,35].

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 [35, 36]. 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 [34], a sharp phase transition [37, 38], the mixed-phase [39,40,41,42,43], and electric charge [44]. 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 [36], how these are affected in the face of a magnetic field [45], and their influences on non-radial oscillations [35]. 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: Sect. 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 Sect. 4 we show the numerical results for equilibrium configurations, radial stability, and tidal deformability of NSs with two-phase. Finally, we conclude in Sect. 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

$$\begin{aligned} T_{\mu \nu }=\left( \rho +p\right) \,u_{\mu }\,u_{\nu }+p\,g_{\mu \nu }, \end{aligned}$$
(1)

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

$$\begin{aligned} ds^2=-e^{\nu }dt^2+e^{\beta }dr^2+r^2d\theta ^2+r^2\sin ^2\theta \,d\phi ^2. \end{aligned}$$
(2)

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

$$\begin{aligned}{} & {} \frac{dm}{dr}=4\pi \rho r^2, \end{aligned}$$
(3)
$$\begin{aligned}{} & {} \frac{dp}{dr}=-(p+\rho )\left( 4\pi rp+\frac{m}{r^2}\right) e^{\beta },\end{aligned}$$
(4)
$$\begin{aligned}{} & {} \frac{d\nu }{dr}=-\frac{2}{(p+\rho )}\frac{dp}{dr}, \end{aligned}$$
(5)

where the parameter m represents the mass inside the sphere of radius r. Equation (4) is known as the hydrostatic equilibrium equation for a spherically symmetric static astrophysical object, also called as Tolman-Oppenheimer-Volkoff equation [47, 48].

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

$$\begin{aligned} m(0)=0,\quad p(0)=p_c,\quad \rho (0)=\rho _c,\quad \textrm{and}\quad \nu (0)=\nu _c.\nonumber \\ \end{aligned}$$
(6)

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

$$\begin{aligned} p(R)=0. \end{aligned}$$
(7)

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:

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

with M being the total mass of the star.

2.2 Radial oscillations equations

The radial pulsation equation is obtained by Chandrasekhar [49] 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 \(\varDelta r/r\) and \(\varDelta p\) [24]; with \(\varDelta r\) and \(\varDelta p\) representing respectively the relative radial displacement and Lagrangian perturbations of pressure. Thus, the system of equations is, \(\xi = \varDelta r/r\):

$$\begin{aligned}{} & {} \frac{d\xi }{dr}=\frac{\xi }{2}\frac{d\nu }{dr}-\frac{1}{r}\left( 3\xi +\frac{\varDelta p}{p\varGamma }\right) , \end{aligned}$$
(9)
$$\begin{aligned}{} & {} \frac{d\varDelta p}{dr}= (p+\rho )\omega ^2\xi re^{\beta -\nu }-4\xi \left( \frac{dp}{dr}\right) +\left( \frac{dp}{dr}\right) ^2\frac{\xi r}{p+\rho }\nonumber \\{} & {} \quad -8\pi p\,(p+\rho )\xi r e^{\beta } -\left( \frac{1}{2}\frac{d\nu }{dr}+4\pi re^\beta (p+\rho )\right) \varDelta p,\nonumber \\ \end{aligned}$$
(10)

where \(\varGamma =\left( 1+\frac{\rho }{p}\right) \frac{dp}{d\rho }\). The variables \(\xi \) and \(\varDelta 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

$$\begin{aligned} \left( \varDelta p\right) _{\textrm{center}}=-3\left( \xi \varGamma p\right) _{\textrm{center}}. \end{aligned}$$
(11)

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

$$\begin{aligned} \left( \varDelta p\right) _{\textrm{surface}}=0. \end{aligned}$$
(12)

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,

$$\begin{aligned} \lambda _1=-\frac{Q_{ij}}{\epsilon _{ij}}, \end{aligned}$$
(13)

where \(Q_{ij}\) is the quadrupole moment perturbed by an external tidal field \(\epsilon _{ij}\) [16,17,18]. The tidal deformability parameter \(\lambda _1\) is connected with the Love number \(k_2\) through the relation \(k_2=\frac{3}{2}\lambda _1R^{-5}\). Moreover, the dimensionless tidal deformability \(\varLambda \) can be written in terms of the Love number \(k_2\) as

$$\begin{aligned} \varLambda = \frac{2}{3}\frac{k_2}{C^5}, \end{aligned}$$
(14)

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

$$\begin{aligned}{} & {} {k_2}=\frac{8C^5}{5}(1-2C)^2[2+C(y_R-1)-y_R] \nonumber \\{} & {} \quad \times \{2C [6-3y_R+3C(5y_R-8)]+4C^3[13-11y_R\nonumber \\{} & {} \quad +C(3y_R-2)+2C^2(1+y_R)]+3(1-2C^2)\nonumber \\{} & {} \quad \times [2-y_R+2(y_R-1)]\ln (1-2C)\}^{-1}. \end{aligned}$$
(15)

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

$$\begin{aligned} r \frac{dy}{dr} + y^2 + y F + r^2 Q = 0, \end{aligned}$$
(16)

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:

$$\begin{aligned}{} & {} F = (\rho -p)\frac{r-4\pi r^3}{r-2m}, \end{aligned}$$
(17)
$$\begin{aligned}{} & {} Q = 4\pi e^{\beta } \left( 5 \rho + 9 p + \frac{\rho +p}{dp/d\rho }\right) - \frac{6e^{\beta }}{r^2}-\left( \frac{d\nu }{dr} \right) ^2.\nonumber \\ \end{aligned}$$
(18)

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 [35, 36, 45, 50, 51]. 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 [36, 52]. 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.,

$$\begin{aligned} \xi _{\textrm{inn}}=\xi _{\textrm{out}}. \end{aligned}$$
(19)

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

$$\begin{aligned} \left( \varDelta p\right) _{\textrm{inn}}=\left( \varDelta p\right) _{\textrm{out}}. \end{aligned}$$
(20)

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

$$\begin{aligned} \left( \xi -\frac{\varDelta p}{rp'}\right) _{\textrm{inn}}=\left( \xi -\frac{\varDelta p}{rp'}\right) _{\textrm{out}}, \end{aligned}$$
(21)

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

$$\begin{aligned} \left( \varDelta p\right) _{\textrm{inn}}=\left( \varDelta p\right) _{\textrm{out}}. \end{aligned}$$
(22)

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. [17], the authors discussed this problem in the context of the surface vacuum discontinuity for incompressible stars. After in Refs. [53,54,55], 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:

$$\begin{aligned} y(r_{tr}+\epsilon ) = y(r_{tr}-\epsilon ) \!-\! \frac{4\pi r_{tr}^3\left[ \rho (r_{tr}+\epsilon ) \!-\! \rho (r_{tr}-\epsilon )\right] }{m(r_{tr})+4\pi r_{tr}^3 p},\nonumber \\ \end{aligned}$$
(23)

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 [56] is adopted. Then, the energy density and fluid pressure of each phase are respectively connected through the relations:

$$\begin{aligned} \rho= & {} \left( \frac{p}{K_{\textrm{inn}}}\right) ^{1/\varGamma _\textrm{inn}}+\frac{p}{\varGamma _{\textrm{inn}}-1},\quad p_{\textrm{inn}}^{\textrm{dis}}\le p\le p_c, \end{aligned}$$
(24)
$$\begin{aligned} \rho= & {} \left( \frac{p}{K_{\textrm{out}}}\right) ^{1/\varGamma _\textrm{out}}+\frac{p}{\varGamma _{\textrm{out}}-1}, \quad 0\le p\le p_\textrm{out}^{\textrm{dis}}. \end{aligned}$$
(25)

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

$$\begin{aligned} K_{\textrm{inn}}=0.0195\times (1.67\times 10^{17}\,\mathrm{kg/m^3})^{-1.34}. \end{aligned}$$
(26)
Fig. 1
figure 1

Some equations of state for two density jump parameters \(\lambda \). The fluid pressure p and energy density \(\rho \) are normalized by the nuclear density \(\rho _{\textrm{nuclear}}=2.68\times 10^{17}\,[\mathrm kg/m^3]\). On the top and bottom panel are respectively adopted \(\rho _{\textrm{inn}}^{\textrm{dis}}=7\times 10^{17}\,[\mathrm kg/m^3]\) and \(8\times 10^{17}\,[\mathrm kg/m^3]\)

Full size image

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_{\textrm{out}}^{\textrm{dis}}=p_{\textrm{inn}}^\textrm{dis}\), the outer polytropic constant takes the form:

$$\begin{aligned} K_{\textrm{out}}=K_{\textrm{inn}}\frac{\left( \rho _{\textrm{inn}}^{\textrm{dis}}-p_\textrm{inn}^{\textrm{dis}}/\left( \varGamma _{\textrm{inn}}-1\right) \right) ^{\varGamma _\textrm{inn}}}{\left( \rho _{\textrm{out}}^{\textrm{dis}}-p_{\textrm{out}}^\textrm{dis}/\left( \varGamma _{\textrm{out}}-1\right) \right) ^{\varGamma _{\textrm{out}}}}. \end{aligned}$$
(27)

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

$$\begin{aligned} \rho _{\textrm{out}}^{\textrm{dis}}=\lambda \,\rho _{\textrm{inn}}^{\textrm{dis}}, \end{aligned}$$
(28)

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 \(\varGamma _{\textrm{inn}}=2.4\), \(\varGamma _{\textrm{inn}}=2.6\) and \(\varGamma _{\textrm{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 \(\varGamma _{\textrm{inn}}, \varGamma _{\textrm{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 [57] and by NICER in [58,59,60,61,62,63,64].

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 Sect. 2, for each \(\rho _{\textrm{inn}}^{\textrm{dis}}\), \(\lambda \), \(\varGamma _{\textrm{inn}}\), and \(\varGamma _{\textrm{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 \(\varDelta 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 [36]).

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

Fig. 2
figure 2

Mass against the central energy density for different values of \(\lambda \), \(\rho _{\textrm{inn}}^{\textrm{dis}}=8\times 10^{17}\,[\mathrm kg/m^3]\), \(\varGamma _{\textrm{out}}=2.4\), and two different values of \(\varGamma _\textrm{inn}\)

Full size image

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 \(\varGamma _{\textrm{inn}}\), four different values of \(\lambda \), \(\rho _{\textrm{inn}}^\textrm{dis}=8\times 10^{17}\,[\mathrm kg/m^3]\), and \(\varGamma _{\textrm{out}}=2.4\). The central energy density goes from \(10^{18}\) to \(10^{20}\,[\mathrm 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 \(\varGamma _{\textrm{inn}}\) is also observed in Fig. 2. For a fixed \(\rho _c\) and \(\lambda \), for a greater \(\varGamma _{\textrm{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 \(\varGamma _{\textrm{inn}}\), neutron stars with more compact cores are obtained (review, e.g., [65, 66]).

Fig. 3
figure 3

Total mass as a function of the radius for four values of \(\lambda \), \(\rho _{\textrm{inn}}^{\textrm{dis}}=8\times 10^{17}\,[\mathrm kg/m^3]\), \(\varGamma _{\textrm{out}}=2.4\), and two different values of \(\varGamma _\textrm{inn}\)

Full size image

Figure 3 shows the mass as a function of the total radius for two values of \(\varGamma _{\textrm{inn}}\), for \(\rho _{\textrm{inn}}^\textrm{dis}=8\times 10^{17}\,[\mathrm kg/m^3]\), few values of \(\lambda \), and \(\varGamma _{\textrm{out}}=2.4\). As in the case of Fig. 2, in this figure is considered \(\rho _c\) between \(10^{18}\) and \(10^{20}\,[\mathrm kg/m^3]\). The panel shows an increment of \(M/M_{\odot }\) with the diminution of R until to find a \(M_{\textrm{max}}/M_{\odot }\). After this point, the curves turn anti-clockwise, to M(R) starts to decrease with R until to reach \(R_{\textrm{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 \(\varGamma _{\textrm{inn}}\). Despite having an increase of the central pressure with \(\varGamma _{\textrm{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.

Fig. 4
figure 4

Total mass against the speed of sound for four values of \(\lambda \), \(\rho _{\textrm{inn}}^{\textrm{dis}}=8\times 10^{17}\,[\mathrm kg/m^3]\), \(\varGamma _{\textrm{out}}=2.4\), and two different values of \(\varGamma _\textrm{inn}\)

Full size image

The total mass versus the speed of sound is plotted in Fig. 4, for \(\rho _{\textrm{inn}}^{\textrm{dis}}=8\times 10^{17}\,[\mathrm kg/m^3]\), some values of \(\lambda \), and \(\varGamma _{\textrm{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 [32].

4.2 Radial stability of neutron stars

Fig. 5
figure 5

The slow (s) and rapid (r) eigenfrequency of the fundamental mode \(\omega _0\) as a function of the total mass \(M/M_{\odot }\) and against the central energy density for \(\rho _{\textrm{inn}}^\textrm{dis}=8\times 10^{17}\,[\mathrm kg/m^3]\), \(\varGamma _{\textrm{out}}=2.4\), four values of \(\lambda \), and two values of \(\varGamma _{\textrm{inn}}\). On the left panels, it is used \(\varGamma _{\textrm{inn}}=2.4\) and on the right panels, it is employed \(\varGamma _{\textrm{inn}}=2.6\)

Full size image

Figure 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 _{\textrm{inn}}^{\textrm{dis}}=8\times 10^{17}\,[\mathrm kg/m^3]\), \(\varGamma _{\textrm{out}}=2.4\), and \(\varGamma _{\textrm{inn}}=2.4\) on the left panels, and \(\varGamma _{\textrm{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 [36], we note that the total mass at the null eigenfrequency of oscillation depends on the type of the phase transition. At this point, when \(\varGamma _{\textrm{inn}}=\varGamma _\textrm{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 \(\varGamma _{\textrm{inn}}=2.6\) and \(\varGamma _{\textrm{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 The central energy density \(\rho _c\), the total mass \(M/M_{\odot }\), the total radius R, and the core radius \(R_{\textrm{core}}\) where the null eigenfrequency of oscillation of the slow case \(\omega _{0,s}\) and rapid case \(\omega _{0,r}\) are derived. These parameter are found for \(\rho _{\textrm{inn}}^{\textrm{dis}}=8\times 10^{17}\,[\mathrm kg/m^3]\), \(\varGamma _{\textrm{out}}=2.4\) and different values of \(\varGamma _{\textrm{inn}}\) and \(\lambda \)
Full size table

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 _{\textrm{inn}}^{\textrm{dis}}=8\times 10^{17}\,[\mathrm kg/m^3]\), \(\varGamma _{\textrm{out}}=2.4\) and some values of \(\varGamma _{\textrm{inn}}\) and \(\lambda \). In the table, for a fixed \(\varGamma _{\textrm{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 \(\varGamma _{\textrm{inn}}\) is increased, i.e., when a stiffer core fluid is considered, stars with a core radius \(R_\textrm{core}\) closer to the total radius R are found.

Fig. 6
figure 6

Top: dimensionless tidal deformability \(\varLambda \) as a function of the total mass in solar masses. The vertical dash-dot-dot line marks the tidal deformability of event GW170817 estimated in [57]. Bottom: the dimensionless tidal deformability \(\varLambda _1\) and \(\varLambda _2\) for a binary NS system with masses \(m_1\) and \(m_2\) and the same chirp mass as the event GW170817 [9]. We only take into account the combination with \(m_1>m_2\). The diagonal short dot line denotes the \(\varLambda _1=\varLambda _2\) limit. The solid top yellow line indicates the \(90\%\) credibility level and the bottom yellow solid line is the \(50\%\) level established by LIGO-Virgo scientific network in the low-spin prior scenario. Only stable equilibrium configurations with slow conversions at the interface are shown

Full size image
Fig. 7
figure 7

Comparison of mass–radius curves with observational data using different values of \(\lambda \), \(\varGamma _{\textrm{inn}}=\varGamma _\textrm{out}=2.4\) and two different values of \(\rho _{\textrm{inn}}^{\textrm{dis}}\). On the top and bottom panel are respectively used \(\rho _{\textrm{inn}}^\textrm{dis}=7.0\times 10^{17}\, [\mathrm kg/m^3]\) and \(\rho _{\textrm{inn}}^\textrm{dis}=9.0\times 10^{17}\, [\mathrm kg/m^3]\)

Full size image

4.3 Tidal deformability in the light of GW170817

Tidal deformability against the total mass of stable NS is plotted at the top panel of Fig. 6 for \(\rho _{\textrm{inn}}^\textrm{dis}=8\times 10^{17}\,[\mathrm kg/m^3]\), \(\varGamma _{\textrm{out}}=2.4\), different values of \(\lambda \) and \(\varGamma _{\textrm{inn}}\). On the panel, it is also presented the tidal deformability constrained by the event GW170817 for a star of \(1.4M_{\odot }\) to be \(70\le \varLambda _{1.4M_{\odot }}\le 580\) [57], 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 \(\varGamma _{\textrm{inn}}\) for different values of the phase transition parameter \(\lambda \), we can notice that the importance of the parameter \(\varLambda \) 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 \(\varLambda _1-\varLambda _2\) are presented for a binary NS system with chirp mass equal to GW170817. Since the total mass is associated with the dimensionless tidal deformability (top panel of Fig. 6), the curves \(\varLambda _1-\varLambda _2\) are obtained once chosen a value of \(m_1\) and calculating \(m_2\) for the fixed value of the chirp mass \({{\mathscr {M}}}=1.188\,M_{\odot }\) [9], which is theoretically calculated by the relation:

$$\begin{aligned} {{\mathscr {M}}}=\frac{(m_1\,m_2)^{3/5}}{(m_1+m_2)^{1/5}}. \end{aligned}$$
(29)

The values considered for \(m_1\) and \(m_2\) run from \(1.36M_{\odot }\le m_1\le 1.60M_{\odot }\) and \(1.17M_{\odot }\le m_2\le 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 _{\textrm{inn}}^{\textrm{dis}}=8\times 10^{17}\,[\mathrm 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\le \lambda \le 0.8\), which plays the compact stars inside of \(50\%\) and \(90\%\) regions. The compact stars with \(\lambda \ge 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 \(\varGamma _\textrm{inn}=\varGamma _{\textrm{out}}=2.4\), some values of \(\lambda \) and two \(\rho _{\textrm{inn}}^{\textrm{dis}}\). On the top and bottom panel are employed \(\rho _{\textrm{inn}}^{\textrm{dis}}=7.0\times 10^{17}\,[\mathrm kg/m^3]\) and \(\rho _{\textrm{inn}}^{\textrm{dis}}=9.0\times 10^{17}\,[\mathrm kg/m^3]\), respectively. In this figure is used \(10^{18}\le \rho _c\le 10^{20}\,[\mathrm kg/m^3]\). The observation data correspond to the NICER constraints obtained from the pulsars PSR J\(0030+0451\) [58, 59] and PSR J\(0740+6620\) [60, 61]. The corresponding bands of the pulsars PSR J\(0740+6620\) [62], PSR J\(0348+0432\) [63] and PSR J\(1614+2230\) [64] are also presented. From the figure, we observe that the change of \(\rho _{\textrm{inn}}^{\textrm{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 _\textrm{inn}^{\textrm{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 _{\textrm{inn}}^{\textrm{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 \(\varLambda _1\)-\(\varLambda _2\) curves for a binary NS system with chirp mass equal to GW170817 considering the relation (29) where \(m_1\) and \(m_2\) goes from \(1.36M_{\odot }\le m_1\le 1.60M_{\odot }\) and \(1.17M_{\odot }\le m_2\le 1.36M_{\odot }\). The inner phase transition energy density considered on the left and right panels are \(7.0\times 10^{17}\,[\mathrm kg/m^3]\) and \(9.0\times 10^{17}\,[\mathrm kg/m^3]\), respectively. In all panels of Fig. 8, only stable equilibrium configurations for the rapid transition case, employing \(\varGamma _\textrm{inn}=\varGamma _{\textrm{out}}=2.4\), are considered. As stated above, the observation data correspond to the event GW170817.

In Fig. 8, on the top panels, it can be seen that all curves tidal deformability (\(\varLambda \))-total mass (\(M/M_{\odot }\)) decay with the increment of \(\rho _{\textrm{inn}}^\textrm{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 GW170817 event. From the bottom panels, we note that some equilibrium solutions located outside the range of observational data derived from the GW170817 event fall within this interval when we increase the phase transition energy density \(\rho _\textrm{inn}^{\textrm{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.

Fig. 8
figure 8

Top: Dimensionless tidal deformability \(\varLambda \) against the total mass in Sun masses. The vertical line of the dash-dot-dot plots the tidal deformability of the GW170817 event measured at [57]. Bottom: The dimensionless tidal deformability \(\varLambda _1\) and \(\varLambda _2\) for a binary NS system with masses \(m_1\) and \(m_2\) and a chirp mass of \(1.4M_{\odot }\) considering the combination with \(m_1>m_2\). The diagonal short dot line marks \(\varLambda _1=\varLambda _2\) limit. The solid top and bottom yellow lines denote respectively \(90\%\) and \(50\%\) level established by LIGO-Virgo scientific network in the low-spin prior scenario. On the left and right panels are used \(\rho _{\textrm{inn}}^\textrm{dis}=7.0\times 10^{17}\,[\mathrm kg/m^3]\) and \(\rho _{\textrm{inn}}^\textrm{dis}=9.0\times 10^{17}\,[\mathrm kg/m^3]\), respectively. Only stable equilibrium configurations with rapid conversions at the interface considering \(\varGamma _{\textrm{inn}}=\varGamma _{\textrm{out}}=2.4\) are presented

Full size image

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 _{\textrm{inn}}^{\textrm{dis}}\), interior polytropic exponents \(\varGamma _{\textrm{inn}}\), and exterior polytropic exponent \(\varGamma _\textrm{out}=2.4\).

As well as in the study of NSs developed in [67], 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 \(\varGamma _{\textrm{inn}}\)), and phase transition energy density.

For the values \(\varGamma _{\textrm{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 [35, 36, 45].

The change of the tidal deformability for a NS (\(\varLambda \)) and a binary NS system (\(\varLambda _1\) and \(\varLambda _2\)), with equal chirp mass as GW170817 event, as a function of \(\varGamma _{\textrm{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 \(\varLambda \) grows and decreases with the increment of \(\varGamma _{\textrm{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 [38].

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 _{\textrm{inn}}^{\textrm{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 _{\textrm{inn}}^{\textrm{dis}}\), and stiffer core fluid \(\varGamma _{\textrm{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 _{\textrm{inn}}^{\textrm{dis}}\) and \(\varGamma _{\textrm{inn}}\).