Probing hadron-quark phase transition in twin stars using 𝑓 -modes

Although it is conjectured that a phase transition from hadronic to deconfined quark matter is possible in the ultrahigh density environment in Neutron Stars, the nature of such a transition is still unknown. Depending on whether there is a sharp or slow phase transition, one may expect a third family of stable compact stars or “twin stars" to appear, with the same mass but different radii compared to Neutron stars. The possibility of identifying twin stars using astrophysical observations has been a subject of interest, which has gained further momentum with the recent detection of gravitational waves from binary neutron stars. In this work, we investigate for the first time the prospect of probing the nature of hadron-quark phase transition with future detection of gravitational waves from unstable fundamental (f) mode oscillations in Neutron Stars. By employing a recently developed model that parametrizes the nature of the hadron-quark phase transition via “pasta phases", we calculate f-mode characteristics within a full general relativistic formalism. We then recover the stellar properties from the detected mode parameters using Universal Relations in GW asteroseismology. Our investigations suggest that the detection of gravitational waves emanating from the f-modes with the third-generation gravitational wave detectors offers a promising scenario for confirming the existence of the twin stars. We also estimate the various uncertainties associated with the determination of the mode parameters and conclude that these uncertainties make the situation more challenging to identify the nature of the hadron-quark phase transition.


INTRODUCTION
Since their discovery, the interior composition of neutron stars (NS) remains one of the most intriguing unanswered questions in astrophysics.Although terrestrial experiments provide hints about the behaviour of matter at high densities, the densities in the core of neutron stars surpass those of laboratories by several orders of magnitude.This makes compact objects the ideal environment to probe matter under extreme conditions of low temperature and high densities not accessible to terrestrial experiments.Quantum chromodynamics (QCD) predicts the appearance of strange matter at high densities via a first-order phase transition from hadronic matter.It is therefore conjectured that deconfined quark matter may appear as a stable component in the high-density environment of the neutron star inner core (Christian 2023; Kumar et al. 2023a;Takátsy et al. 2023;Shirke et al. 2023;Blaschke et al. 2020a).The phase diagram of QCD corresponding to dense nuclear matter in the hot and low density regime has been explored by laboratory experiments, namely relativistic nuclear collisions in the Large Hadron Collider (LHC) as well as by numerical solutions of Lattice QCD.The result is a smooth phase transition or crossover.Both the intermediate and lower temperature and density regions are not very well studied experimentally or theoretically, therefore the nature ★ E-mail: bikramp@iucaa.in(KTS) of the phase transition is quite uncertain (Kumar et al. 2023a).In particular, compact star matter located at the low-temperature and high-density region of the QCD phase diagram can potentially feature a phase transition from hadronic matter to deconfined quark matter.This phase transition is conjectured to be a first order which might or might not feature geometrical shapes at the interface.Alternatively, there are proposals that state that the hadron-quark phase transition ought to be a crossover, see for instance Baym et al. (2018); Sotani & Kojo (2023).Within this work, the first possibility namely the QCD first order phase transition is addressed.
Although equations of QCD remain unsolved at energy scales relevant for describing NS matter, it is possible to model the NS interior via equations of state (EOS) based on microscopic or phenomenological density functional theories.Through the EOS, the NS interior composition can then be connected to global NS properties, such as its mass or radius in electromagnetic observations or tidal deformation of the components of a binary system during a merger from gravitational wave emission.It is expected that the appearance of strange degrees of freedom would result in a softer EOS and correspondingly to a lower maximum NS mass.This can then be compared to mass measurements from NSs in binary systems to test whether a particular EOS can support the observed maximally massive NS.Moreover, neutrinos from core-collapse supernova, from accreting compact stars, and gravitational wave signals from their binary mergers produced at the postmerger phase can potentially identify and constrain the properties of a first order phase transition of hadronic matter into quark matter (Khosravi Largani et al. 2023b,a;Bauswein et al. 2019).
Further, depending on the nature of the hadron-quark phase transition, distinct branches in the NS mass-radius relation may appear (Glendenning 1992;Blaschke et al. 2020a;Zacchi et al. 2017;Espino & Paschalidis 2022).In recent years, a significant amount of investigations have been carried out connecting the astrophysical observations and twin stars (Christian & Schaffner-Bielich 2022, 2020;Tsaloukidis et al. 2023;Montaña et al. 2019).The detection of compact star mass twins, i.e., two stars of about the same mass but different radii, would be a smoking gun evidence for a strong first order phase transition.On the contrary, a smooth phase transition that might be caused by the appearance of non-homogeneous matter or pasta phases at the quark-hadron interface might not result in a third branch of compact stars in the mass-radius diagram.Studies carried out in (Ayriyan et al. 2018;Blaschke & Alvarez-Castillo 2020;Blaschke et al. 2020b;Maslov et al. 2019) reveal an effective assessment of the amount of quark-hadron pasta that could support the neutron star while preserving the third, disconnected compact star branch.
Recently there has been significant interest in trying to identify whether hadron-quark transition occurs via either a rapid or slow process (Pereira et al. 2018) at the boundary between pure phases (accompanied by an energy density jump sometimes resulting in compact stars twins (Landry & Chakravarti 2022;Alvarez-Castillo 2021)).In this work we also consider a phase transition via a mixed phase (Abgaryan et al. 2018;Ayriyan et al. 2018;Ayriyan & Grigorian 2018;Ayriyan et al. 2021).Several phenomenological interpolation schemes were proposed to mimic the hadron-quark mixed phase in compact stars via geometrical structure of different shapes denoted as "pasta phases" (Abgaryan et al. 2018;Ayriyan et al. 2021), and their effect on properties of hybrid stars were investigated.(Pereira et al. 2018) studied whether the nature of the phase transition can be assessed from future detections of global NS parameters such as mass, radius and tidal deformability.Landry and Chakraborty (Landry & Chakravarti 2022) explored the possibility of distinguishing neutron stars from compact twins from future detections of tidal deformability of binary NS mergers.
The LIGO-VIRGO collaboration has, until now, directly witnessed gravitational waves (GW) from binary neutron stars (BNS), and the discovery of the BNS event GW170817 along with the electromagnetic counterpart opens an exciting new chapter in multi messenger astronomy (Abbott et al. 2018(Abbott et al. , 2017a(Abbott et al. ,b, 2019a)).In addition to the binary system, isolated NSs can also emit GWs by non-axisymmetric deformations.The parameters of the quasi-normal modes (QNM) are dependent on the stellar interior.Therefore, knowing how QNM parameters behave in relation to the interior of NS can aid in revealing the interior of NS from QNM observations in the future.The properties of the different oscillation modes of hybrid stars with phase transitions have been the subject of extensive research (Flores & Lugones 2014;Ranea-Sandoval et al. 2018;Rodríguez et al. 2021;Ranea-Sandoval et al. 2022a;Zhao et al. 2022;Zhao & Lattimer 2022;Constantinou et al. 2021;Jaikumar et al. 2021;Kumar et al. 2023c;Constantinou et al. 2023;Kumar et al. 2023b;Ranea-Sandoval et al. 2023).It has been noted that the quadrupolar  -mode, one of the QNMs, is among the most promising sources of GW emission and is also accessible with the improving sensitivity of the current LIGO-VIRGO detectors or with the next-generation GW detectors, such as the Einstein Telescope (ET) and Cosmic Explorer (CE) (Ho et al. 2020;Pradhan et al. 2023a).Alternatively, asteroseismology is another approach, where EOS independent or universal relations (URs) (Yagi & Yunes 2013;Largani et al. 2022) are employed to infer the NS interior from the detection of QNM parameters (Völkel et al. 2021;Völkel & Krüger 2022;Pradhan et al. 2023a;Ranea-Sandoval et al. 2022b).According to a recent study from Ranea-Sandoval et al. (2023), the inclusion of a slow phase transition ( the conversion speed is slower compared to the radial perturbation time scale ) appears to violate the conventional universal relations involving  -mode parameters and other NS observables.Additionally, the recent study by Laskos-Patkos & Moustakidis (2023) examines the r-mode instability windows and spin-down evolution of twin stars.It suggests that hybrid equations of state, depending on the transition density, could explain the observed spin and temperature evolution, and the future detection of r-mode GW could serve as an indicator of the presence of twin stars.
In this work, we investigate the possibility of distinguishing between a sharp and slow hadron-quark phase transition via hypothetical future observations of gravitational waves from  -mode oscillations in NSs.Using the methodology developed in Abgaryan et al. (2018), we parametrize the transition from hadronic to quark matter EOS beyond the Maxwell point using a single parameter, the pressure increment at the critical chemical potential.We study the effect of the nature of the phase transition on  -mode characteristics.Further, we examine whether detection of gravitational waves from current or future  -mode observations can constrain the nature of phase transition convincingly.
The paper is organized in the following way.In section 2, we discuss the model applied to describe the NS interior as well as its global structure.The results of our study are presented in section 3 and we summarize our conclusions in section 4.

EOS Model
The pressure density relationship, referred to as the equation of state (EOS), is vital in connecting the microscopic behavior of NS matter (NSM) to the NS observables.As suggested in the literature, a deconfined quark phase could appear in high-density NSM.The nature of the phase transition from hadronic to deconfined quark phase is still a matter of debate, and it is unclear whether there exists a jump in the thermodynamic variables or whether the transition is a smooth crossover.A significant number of works have described the transition either by the Maxwell or Gibbs construction.Furthermore, consideration of finite surface tension effects can lead to the existence of a mixed phase with different geometric shapes (referred to as the "pasta" phase) (Ravenhall et al. 1983;Voskresensky et al. 2003;Yasutake et al. 2014).The description of the crossover-like transition from hadronic to the deconfined quark phase was first discussed in Masuda et al. (2013) and later adequately in (Alvarez-Castillo & Blaschke 2017;Alvarez-Castillo et al. 2017;Abgaryan et al. 2018).
For this study, we consider the four-parameter realization of (Alvarez-Castillo & Blaschke 2017; Paschalidis et al. 2018) abbreviated as "ACB4" to describe the NSM at densities higher than saturation density  0 = 0.15fm −3 , for the reasons outlined below.The ACB4 EOS model satisfies the Seidov condition at the phase transition to produce high-mass twin stars (Seidov 1971).The EOS model mimics the pasta phases and is relatively simpler to construct as required for the NS physics.The pressure () as a function of density () can be expressed as; where each density region is described by the pre-factor   and polytropic index Γ  .It is convenient and thermodynamically consistent to express the pressure as a function of chemical potential () as, where  0, is the effective mass of the constituent: nucleon effective mass for the hadronic region and quark effective mass in the quark phase.Further, using the thermodynamic relations, the number density and energy density () can be obtained as a function of  following the methodology given in (Alvarez-Castillo & Blaschke 2017;Alvarez-Castillo et al. 2017).The EOS model describes the hadron quark phase transition by the Maxwell construction with a constant transition pressure   =  2 at a critical chemical potential (  ), making the pressure for both phases equal, i.e.,   (  ) =   (  ) =   .The second region with polytrope Γ 2 = 0 describes the first-order phase transition.In the literature, there are two different main approaches (i) replaced interpolation method (RIM) and (ii) mixed interpolation method (MIM) that have been used to mimic the pasta phase and construct the mixed phase Abgaryan et al. (2018).We consider the RIM approach and replace the pressure in the relevant region of the hadronic and quark phase near the Maxwell critical point (  ,   ) by a parabolic polynomial function defined as (Ayriyan et al. 2018;Ayriyan & Grigorian 2018); where Δ = Δ *   represents the additional pressure of the mixed phase with   = (  ) being the critical pressure of the Maxwell construction.The mixed phase polynomial (3) connects the EOS smoothly from the hadronic phase to the mixed phase at   and mixed phase to quark phase at   .The unknowns  1 ,  2 ,   ,   are determined by the continuity of thermodynamic quantities  and  at   and   given as; We tabulate the EOS model parameters corresponding to the ACB4 realization in table 1.We display the EOSs resulting from the described model in fig.1a characterized by Δ arbitrarily ranging from 0 up to 8%, where Δ = 0 represents the first-order Maxwell construction.We observe that for Δ values above 5% the third family branch disappears implying that all the compact stars lie in a connected, second family branch.This result clearly shows how robust the existence of the third family of compact stars is within this EoS framework against the appearance of pasta phases.

Macroscopic Observables
For a given EOS  = (), the compact star (CS) configurations are obtained by integrating the general relativistic hydrostatic equilibrium Tolman-Oppenheimer-Volkoff (TOV) equations (Tolman 1939;Oppenheimer & Volkoff 1939), . (5) Integrating TOV Eqs (5) from the center ( = 0, with a central pressure  =   ) to the surface of the star with the boundary condition that the pressure vanishes at the surface, i.e, () = 0, one can obtain the radius ( ) of the star.The mass enclosed within  is the stellar mass () i.e.,  = ().
In a binary system, the CSs get tidally deformed due to the tidal field of the companions, and the deformation can be measured from the GW observation of the binary system represented by the tidal deformability.The tidal deformability can be used to constrain the CS EOS as it strongly depends upon the CS EOS.The dimensionless tidal deformability (Λ) is expressed in terms of the love number  2 in (6).Furthermore, one needs to solve an additional set of differential equations along with the TOV equations Hinderer (2008) to obtain the love number  2 , where  is the compactness i.e., the ratio of mass and radius  = /.We display the hybrid EoSs in fig.1a and the corresponding  −  and  − Λ relations in figs.1b and 2 respectively.The second and third family branches are connected using the dotted lines, showing the unstable configurations for    ⩽ 0. On increasing the value of Δ, the unstable region gets shortened, and beyond Δ = 5%, the second and third families merge to form a single branch.The discontinuity in  and Λ as a function of  indicates the presence of twin stars.Hence, the simultaneous inference of  −  or  − Λ from observations can reveal the existence of discontinuity in  −  or  − Λ relations and can provide hints to the presence of twin stars or even the nature of phase transition.

Solving for 𝑓 -mode Characteristics
Among the different quasi-normal modes of the NS, the nonradial fundamental mode (  -modes) is strongly coupled with the NS fluid (Thorne & Campolattaro 1967) and the dominant source of GW emission.In the literature, a significant amount of effort has been spent in developing the methodology to solve for mode characteristics including the resonance matching method (Chandrasekhar & Ferrari 1991), direct integration method (Lindblom & Detweiler 1983; Detweiler   & Lindblom 1985), method of continued fraction (Leins et al. 1993;Sotani et al. 2001) and WKB approximation (Andersson & Kokkotas 1996).Though several works used the Cowling approximation method to find mode frequency ignoring the metric perturbation, the importance of inclusion of linearized general relativistic treatment over the Cowling approximation has been discussed in several recent works (Yoshida & Kojima 1997;Chirenti et al. 2015;Pradhan et al. 2022) which concluded that the Cowling approximation could overestimate the  -mode frequency upto ∼ 30% compared to the frequency obtained within general relativistic treatment.
In this work, we obtain the mode parameters by solving the perturbations in full general relativistic treatment.We use the direct integration method developed in (Detweiler & Lindblom 1985;Sotani et al. 2001;Pradhan et al. 2022) to find the NS  -mode frequency.Shortly, the coupled perturbation equations for perturbed metric and fluid variables are integrated within the NS interior with appropriate boundary and junction (for hybrid stars with density discontinuity) conditions (Sotani et al. 2001).Further, the fluid variables are set to zero outside the star, and then Zerilli's wave equation (Zerilli 1970) is integrated too far from the star.Further, a search is carried out for the complex  -mode frequency ( = 2  +    ) corresponding to only outgoing wave solution to the Zerilli's equation at infinity.The real part of  represents the  -mode angular frequency, and the imaginary part represents the damping time.For finding the mode characteristics, we use the numerical methods developed in our previous work (Pradhan et al. 2022) along with including jump conditions for hybrid stars from Sotani et al. (2001).

𝑓 -mode parameters and Universal Relations
We display in fig.3a  transition.However, the frequency jump disappears beyond Δ = 5%.Therefore, simultaneous observations of  -mode frequency (or damping time) with mass can lead us to comment on the presence of the jump or the nature of phase transition.The simultaneous mass and mode frequency measurements can be obtained from the detection of GW events from a binary system (Williams et al. 2022).We discuss different possibilities and prospects of binary systems involving twin stars in section 3.5.However, in the case of isolated stars, we may not have the privilege of measuring the mass.The QNMs are the only source of GW emission, implying the detectable parameters are  and .Hence, one has to rely on asteroseismology to infer the stellar properties or, inversely, to infer the interior of the compact star.
The idea of asteroseismology involving the inference of the stellar observables from detection of NS quasi-normal modes (QNMs) was first addressed in Andersson & Kokkotas (1996, 1998), where it was shown that there exist EOS-independent relations or universal relations (URs) among the NS mode parameters and the NS observables like , .Though the initial empirical relations involve the mean density of the star, later works have shown that the empirical relations may depend upon the EOS model considered (Pradhan & Chatterjee 2021;Pradhan et al. 2022;Lopez et al. 2022).However, the URs involving stellar compactness (/) are EOS independent and can be used to reconstruct NS observables from the detection of QNMs (Tsui & Leung 2005;Lioutas & Stergioulas 2018).We obtain the URs considering a wide range of EOS models considered in Pradhan et al. (2023b) 7) and ( 8).We display URs involving compactness,  -mode characteristics, and fit relations in fig. 4 and tabulate the fit parameters of eqs.( 7) and ( 8) in table 2.
and for the UR involving mass-scaled damping time (/) or Im() and compactness (/) can be given as, We present the dependence of scaled complex QNM frequency (scaled with NS mass) as a function of compactness along with the Universal Relations (URs) from this work and previous studies in fig. 4. We display the URs proposed in the literature, notably those from (Tsui & Leung 2005)  Λ with radius scaled tidal deformability parameter Λ 10km and can be written as, where,  = log [Λ 10 ] with  10 = /10km.Hence, there will be 2 URs for eq.( 9): one for  -mode angular frequency Re() and the other for the damping time, and similarly 2 URs for eq.( 10).The URs for  Thus, the twin stars discussed in this study should not be confused with SSHS.

RESULTS
As discussed, GW asteroseismology aims to recover the stellar properties from the detected mode parameters using the URs.We discuss the impact of uncertainties associated with the URs on the reconstruction of stellar properties from  -mode parameters in section 3.1.Additionally, the detection of GWs from  -modes itself can be associated with uncertainty in the mode parameters, which will reflect on the stellar properties, which is discussed in section 3.2.Furthermore, log Table 4. Fit parameters for real (  ) and imaginary (  ) of the complex frequency , related to  = log [Λ 10 ] through the URs eq. ( 10).The coefficient corresponding to   and   are denoted with subscript  and  respectively.
we discuss the inverse problem to constrain the compact star EOS from  -mode GW observation in section 3.3.

Effect of uncertainty in the URs
As mentioned in the section 2.4, the idea behind constructing the URs is to reconstruct the NS observables like ,  from detecting mode parameters.For an observed  -mode frequency , the URs eq. ( 7) will result in one curve in the  and / plane corresponding to Re() and similarly for the observed damping time , the UR eq. ( 8) will result in another curve in the  and / plane corresponding to the Im().The intersection of the two resulting curves will give us the value of  and /, hence .As one of our main concerns is to probe the nature of the phase transition, which can be decided from the value of Δ or focusing on the particular  −  region of the fig.1b: mainly in the mass range 1.8 ⊙ to 2.0  ⊙ for the EOS model considered here.Hence, if the observed  and  along with their uncertainties of one EOS Δ ≠ 0% do not overlap with the unstable region of the  −  curve of the Δ = 0, then that could lead us to confirm that the  -mode observations can differentiate the nature of Δ.So, assuming that the  -mode parameters are observed precisely for a few randomly selected stars, the  and  can be estimated using URs, as explained before.The consideration of the uncertainties associated with the UR parameters will further result in uncertainties of the recovered stellar observables such as  and .
Under the assumed scenario of the precise measurement of the mode parameters, the mass and radius recovered for different EOSs sing URs eqs.( 7) and ( 8) with different Δ along with their uncertainties are displayed in the fig. 5.It is clear that the uncertainty region recovered for  and  for Δ = 8% are distinguishable from that of recovered for Δ = 0%, particularly in the region where the  −  sequence of Δ = 0% have an unstable region (connecting the second and third family stable branches) within it (see fig. 5a).This indicates that observing  -mode parameters can lead to commenting about the value of Δ and hence the nature of phase transition.However, if one considers the EOSs with Δ = 5%, the recovered regions for  and  overlap with - recovered from EOS with Δ = 0% and indicated no distinguishability.Hence although the observations of  -modes and the use of URs can differentiate the - region for EOSs with Δ = 8% from EOS with Δ = 0%, it fails to differentiate the  −  region between Δ = 0% and Δ = 5%.
The analysis can be performed using the URs from eqs. ( 9) and ( 10).The methodology of recovering stellar properties using the URs eqs.( 9) and ( 10), have been demonstrated in Ranea-Sandoval et al. (2022a) for axial w modes, which can be extended to the mode.Performing similar tests as mentioned above but with the URs eqs.( 9) and ( 10), the recovered  and  regions for EOSs Δ = 0% and Δ = 8% are displayed in fig.6.We notice that using any pair of URs results in similar conclusions regarding the distinguishability of the value of Δ.However, in the case of using URs eqs.( 9) and ( 10), the URs eq. ( 9) provides the measurement of (,Λ) and use of eq. ( 10) recovers (,Λ) hence, during the recovery, it is needed to be checked that the resulting Λ with different URs should be same or should be within minimal measurement uncertainty.In principle, to get the joint posterior of (, ), one has to eliminate the Λ from (,Λ) and (,Λ) obtained using eq.( 9) and eq.( 10) respectively.

Effect of Observational Uncertainties
Section 3.1 demonstrates the  -mode asteroseismology problem dealing with the uncertainty on the URs under the assumption that the mode parameters are measured precisely.However, detecting mode will always result in some uncertainties in the mode parameters, which we address in this section while discussing the fact that whether or not future observation can help us distinguish the nature of phase transition.The  -mode GW signal for a source at a distance , with frequency (  ) and damping time scale (), can be modeled as a damped sinusoidal (Kokkotas et al. 2001;Ho et al. 2020), where, One needs the value of the  gw for proceeding further.There are different phenomena have been used to stimulate the excitation of the NS oscillation modes: like mini collapse (Lin et al. 2011), newborn NSs (Ferrari et al. 2003), star quakes (Keer & Jones 2014;Mock & Joss 1998), magnetars (Abbott et al. 2019b;Abbott et al. 2022b), the pre-merger (Steinhoff et al. 2016;Andersson & Ho 2018) and post-merger stages of a NS in binary (Shibata 1994;Stergioulas et al. 2011;Bauswein & Janka 2012).Even though the connection between pulsar glitches and the mode excitation is not clear, the assumption of  -mode excitation with an energy similar to that of typical pulsar glitches has been widely considered while discussing NS seismology (Andersson & Comer 2001;Andersson 2021) and the detectability (Ho et al. 2020;Abbott et al. 2019b;Abbott et al. 2021Abbott et al. , 2022a,b) ,b) of transient  -mode GW signal.A significant number of explorations have been made regarding the relationship between glitches and  -mode excitations in recent works (Keer & Jones 2014;Yim & Jones 2023;Lopez et al. 2022).Recent works have also been performed on  -mode GW searches with the assumption of mode excitation with energy typical to that of pulsar glitches by LIGO-VIRGO-KAGRA (LVK) collaboration (Abbott et al. 2019b;Abbott et al. 2022aAbbott et al. ,b, 2021)).Hence to consider the observational uncertainties in this work, we consider the  -mode excitation in the isolated glitching pulsars with energy same as the typical energy of pulsar glitches.Now assuming that GW energy  gw is supplied by the energy of the glitch; one can have (Ho et al. 2020), where  and  are the moments of inertia and spin frequency, respectively, whereas Δ  is the relative change in spin frequency of a glitch event.
To demonstrate the asteroseismology problem, we consider an mode GW event from Vela pulsar with an energy corresponding to the strongest glitch of the Vela pulsar and assign a random mass 1.75 ⊙ .Then the other stellar properties required for modeling the  9) and ( 10).
GW signal as per eqs.( 11) to ( 13), such as  , ,  corresponding to a 1.75 ⊙ star is assigned from the hybrid EOS model with Δ = 0%.Then we perform the parameter estimation (PE) of the GW parameters in a Bayesian framework using the nested sampling algorithm dynesty (Speagle 2020), as implemented on the GW inference package bilby (Ashton et al. 2019).We keep a log uniform prior in  gw ∈ logU[10 −25 , 10 −4 ], uniform prior in  ∈  [500, 4000]Hz, uniform prior in  ∈ U[0.05, 0.5]s.We keep the distance  and sky positions fixed at their observed values.We consider two GW network configurations: first, 2 LIGO detectors H1, L1 operating at O5 sensitivity Abbott et al. (2020) 1 as anticipated for the 5th observation run (A+), and then consider the next generation Einstein 1 https://dcc.ligo.org/LIGO-T2000012/publictelescope (ET) with ET-D sensitivity (Hild et al. 2011) 2 .We display the joint distribution of (  , ) (marginalised over  gw ) recovered with A+ and ET in fig.7a.Further, we reconstruct the (, ) from the recovered posterior (  , ) using the URs eqs.( 7) and ( 8) and display the recovered (, ) in fig.7b.It is clear that in both A+ and ET configuration, with a 90% credible interval (CI), the frequency  can be estimated within ⩽ 1%.However, for the same source, the damping time  can be recovered within ∼ 5% in ET compared to the ∼ 20% error recovered with A+ configuration.For the injected scenario in A+, the  and  with a 90% credible interval are recovered within ∼20% and ∼ 6%, respectively.In ET, the mass and radius, with a 90% credible interval, are recovered within ∼6% and ∼ 2%, respectively.It is not always that the posteriors of recovered ( − ) is unique as shown in fig.7b.The URs that have given in eqs.( 7) and ( 8) are not linear, and that can result in more than one solution in the (, ) plane for a given pair of (  , ).In the case of multiple solutions, one can further combine different types of URs wisely to find out the correct solution for (, ) better representing the observational data (  , ).Furthermore, we only use the URs eqs.( 7) and ( 8) to recover the stellar observables from the posterior of (  , ) and do not consider the URs eqs.( 9) and ( 10), as one needs to perform additional works on eliminating the Λ to get the joint posterior of (, ).Here, we also ignore the uncertainties on the URs as we focus on the observational uncertainties.

Can we differentiate the nature of Δ𝑝 from future 𝑓 -mode GW observations?
To check whether the future observations of  -modes can help us distinguish the nature of Δ from the  −  plane, the following tests are performed (as the uncertainties on the measurements of (, ) in A+ is quite high, we investigate this with ET): • Few  -mode GW events are considered with random masses with assumption of the particular EOS is Δ = 0%.A random glitch energy is assigned to each event such that the Signal to Noise Ratio (SNR) in ET is ⩾ 103 .We perform a parameter estimation using Bilby to get the posteriors of (  , ).From the recovered (  , ), we reconstruct the (, ) using the URs eqs.( 7) and ( 8).
• We repeat the above exercise considering the EOS with Δ = 8% by choosing a few  − configurations where there is an unstable region in the  −  plane for the EOS Δ = 0%.
From fig. 8, one can conclude that there are overlapped regions of recovered  −  for different EOSs.Additionally, the uncertainty on the  −  recovered for the EOS model with Δ = 8% overlaps with the particularly unstable region connecting the second and third family of the stellar configurations in the  −  relation of the Δ = 0% EOS model, and one can barely distinguish the value of Δ based upon the recovered  − .However, it is clear that under the assumption of Δ = 0%, i.e., in the case of the presence of twin stars, the  −  of the star from the second family and its twin companion from the third family is clearly distinguishable, which can lead to confirming the presence of twin stars.We have also considered the scenario distinguishing the EOS models with Δ < 8% other than 8% from Δ = 0%.However, one can notice that even the consideration of the extreme cases of Δ (0% and 8%) in this work results in overlapping regions in the recovered  −  and further consideration of lower Δ values makes the distinguishability of the nature of Δ more challenging.7) and ( 8) from a posterior of recovered (  , ) for a random choice of masses and glitches injected in ET.The injections are shown with empty squares (in red for Δ  = 0% and black for Δ  = 8%).The uncertainties are shown in blue (green) for EOS with Δ  = 0% (Δ  = 8%).

Inverse Problem: constraining the EOS parameters
Future observations of  -mode GW events can be further used to constrain the interior of the compact stars by using the observed  -mode parameters (Völkel et al. 2021;Völkel & Krüger 2022;Pradhan et al. 2023a).As the mode parameters can be observed more precisely with the next generation GW detectors, the posterior can then be translated to get a better constraint on the EOS model parameters.Furthermore, depending upon the nature and posterior of Δ, further inferences can be drawn regarding the nature of the preferred phase transition or even the existence of the pasta phase.

Sensitivity to Twin star detection for low mass twins
Depending upon the onset of phase transition, an early phase transition from hadron to quark matter can result in the existence of low-mass twin stars.In our previous discussions, particularly in sections 3.1 and 3.2, we explored various perspectives on the detectability of high-mass twin stars.Before we conclude, we also examine our methodology in the context of low-mass twin stars using the ACB5 parameterized EOS model (see detailed description in (Paschalidis et al. 2018)).We present the Equation of State (EOS) and the corresponding Mass-Radius ( − ) relations for the ACB5 model in figs.9a and 9b respectively.Notably, compared to the ACB4 EOS model, the ACB5 EOS model indicates a phase transition occurring at a lower density, with the transition onset at around 1.4 ⊙ .
For the ACB5 EOS model, the second and third families of compact stars merge into a single branch for Δ ⩾ 2%.Furthermore, we illustrate the variation of tidal deformability with mass for the ACB5 model in fig.10.Interestingly, the differences in compact star radii (Δ) and tidal deformability (ΔΛ) between hadronic neutron stars (NS) and their twin companions are significantly reduced compared to the high-mass twins.In an optimistic scenario involving the detection of binaries with next-generation gravitational wave (GW) detectors, it becomes possible to measure ΔΛ within approximately 15% (Landry & Chakravarti 2022).This indicates that hybrid stars can be distinguished using these advanced GW detectors, where the tidal deformability Λ differs by more than 15% compared to their hadronic companions.We display the  -mode characteristics for the ACB5 Model in fig.11.To distinguish the low-mass twins, one needs a precise measurement of mode parameters simultaneously with the mass measurement.Though GW from the binary system simultaneously measures mass and mode frequency, the more significant errors in the  -mode frequency Williams et al. (2022) make the distinguishability of low-mass twins more challenging.For the ACB5 EOS model with Δ = 0%, the  -mode frequency of the third family twin star differs by 5% compared to its second family pair at the onset  ∼ 1.4 ⊙ .Hence, a future detection of the  -mode frequency of the pulsars with precisely known mass (from other observations, such as radio) within ⩽ 5% may put insights into the existence of low-mass twins.Furthermore, from detecting  -mode GWs from the pulsars (nearby pulsars are more likely) with unknown macroscopic properties, the conclusion can be made following section 3.2: the detection of mode GWs with ET can measure the compact star radius  up to ∼ 2% (at a 90% CI, i.e.,  ,90 ).To completely distinguish the twins, the radii of the twins should be separated at least by 2 ,90 such that the two posterior distributions corresponding to the measurement of the radii of twins do not overlap with each other4 .Twins separated by radius Δ greater than 2 ,90 can be distinguishable.Hence, from the detection of  -mode GWs from two compact objects having the same mass with 2 ,90 Δ ⩾ 1, one might clearly distinguish the twins.Looking at fig. 7 of section 3.2 we have 2 ,90 ∼ 0.4 km for the  -mode detection with ET, which concludes that twins having separation ⩾ 0.4km can be distinguished from the future  -mode observations.
The investigations of Alvarez-Castillo (2021) demonstrated the possibilities of the existence of twin stars at different onset masses, and the results can be used to comment on the detection of twins at other mass regions using  -mode GW observations.Our analysis suggests that confirming the presence of low-mass twins and identifying the earliest phase transition may be less feasible even with  -mode GW detection compared to high-mass twins.However, various population studies have indicated that a compact star mass distribution peak occurs around ∼ 1.4 ⊙ (de Sá et al. 2023) and given the statistically significant number of observations near this mass range, combining data from different observations can offer valuable insights into the existence of low-mass twin stars (if they indeed exist).

Twin Stars in Binary System
The binary neutron star system provides an excellent scenario for measuring NS mass and frequency.However,  -mode parameters can be constrained for GW events in binary systems having high SNR by the next generation GW detectors (Williams et al. 2022), while the leading order parameters are chirp mass M  (see eq. ( 14)) and tidal parameter Λ (see eq. ( 15)) can be well constrained for events detected by both current and next generation GW detectors.As the measurement of M  and Λ can be used to comment on the presence of strong or crossover phase transition, additional information regarding the nature of the phase transition can be addressed by the simultaneous measurement of mass and  -mode frequency.Considering a series of detections, Landry & Chakravarti (2022) recently discussed the prospects of the detections of twin stars using the jump in the  − Λ plane, and the Λ − M  behavior has been used to comment on the presence of the twin branch, where . In eqs.( 14) to ( 16),   , Λ  ,   are the mass, quadruple tidal deformability and  -mode angular frequency of the   ℎ companion, respectively.  =   /( 1 +  2 ).
The  -mode frequency may be detected within a few 100 Hz from a BNS event in the era of XG GW detectors Williams et al. (2022) through the dynamical tidal parameter   (see eq. ( 16)).Simultaneously, the leading order tidal parameters would be more precisely measured.Looking at fig. 12a, it is clear that there is a jump in Λ − M  plane for binaries with at least one twin star, compared to the binaries with only stars from the second family branch.Similar behavior can be seen in fig.13a for the dynamical parameter   involving the  -mode frequency.However, if one considers an EOS with a crossover-like phase transition and with no discontinuous twin branch (say the EOS with Δ = 5%), there is no discontinuous jump in Λ or   as seen for the EOS with Δ = 0% (see fig. 13b).So the presence of the jump in Λ − M  or   − M  can help us indicate the presence of a twin branch or the presence of a strong phase transition.

DISCUSSIONS
In this work, we investigated for the first time whether future detection of gravitational waves from  -mode oscillations could allow us to probe the nature of the hadron-quark phase transition in NSs.The nature of the phase transition is intimately related to the stiffness of dense matter.The effect of the sharp first-order phase transitions it to soften matter just as the appearance of pasta phases in our study does.On the other hand, crossover phase transitions might either soften or stiffen matter at the phase transition, see (Sotani & Kojo 2023) for examples of the latter.The general tendency of the  -modes curves as a function of stellar mass in the case of stiffer matter is to lower their values with respect to the pure hadronic stars case as opposed to increased values for the case of softer matter, as can be seen in figure 4 of (Sotani & Kojo 2023) for stiffening matter and in the left panel of figure 3 of this work for softening matter.In our study we employed a recently developed phenomenological interpolation scheme from Abgaryan et al. (2018) to mimic the thermodynamic behaviour of the mixed phase via pasta phases, in which the nature of the phase transition is described in terms of a single parameter, the pressure increment (Δ) at critical chemical potential.Within this EOS scheme, we calculated the global NS properties such as mass, radius and tidal deformability as well as  -mode characteristics in full general relativistic framework.We systematically investigated how these observable NS properties are affected by the variation in parameter Δ, i.e., in going from a sharp (Δ = 0) to a smooth (Δ = 8%) phase transition.
For the case of isolated NSs, we then used Universal Relations (URs) to recover stellar properties such as mass and radius considering future detection of  -modes excited by NS glitches.We further analysed whether, given the uncertainties in the URs, the measurement of masses and radii can allow us to comment on the nature of the hadron-quark phase transition.We concluded that if the  -mode characteristics are precisely known, for a given UR one may be able to distinguish between a strong phase transition (Δ = 0) and a smooth one (Δ = 8%), but the distinguishability depends on the value of Δ.We further considered uncertainties in the  -mode observations and found that this would further decrease the distinguishability between the mixed-phase scenarios.However, in all the cases, twin stars (corresponding to the strong phase transition Δ = 0) can be clearly distinguished from a normal NS.
We also examined whether GW observations from  -modes can be used to comment on the presence or distinguishability of the low mass twin stars.Contrary to high mass twins, we discover that applying the same methodology it becomes more challenging to detect low mass twin stars when the  -mode detections are made from compact stars of unknown mass.We go over various astrophysical observational  scenarios and estimate the amount of observational precision that would be necessary to distinguish between twins at different twin star onset masses.
For the case of binary NSs, we probed whether future observations of dynamical  -modes can lead to constraints on the EOS.We found that although one may be able to distinguish NSs from twins using the detection of the dynamical tidal parameter, the constraints from the tidal parameter in the case of binary NSs provide better evidence for the existence of twin stars, which supports the scenario of a sharp hadron-quark phase transition.
Previous studies Ayriyan & Grigorian (2018); Pereira et al. (2022); Landry & Chakravarti (2022) attempted to differentiate between the strong hadron-quark phase transition or smooth crossover by investigating their effects on the mass, radius or tidal deformability.(Pereira et al. 2022) concluded that even considering the most optimistic case for future generation GW detectors, distinguishing a sharp phase transition from a mixed state may be observationally challenging.Landry and Chakravarti Landry & Chakravarti (2022) analyzed the number of binary NS observations required to infer the existence of twin stars from the measurement of tidal deformation.Recently, in the work of Suleiman & Read (2024), the broadening of a few relevant URs due to the consideration of different physical constraints and the impact on future NS measurements is discussed.The effect of twin stars or hybrid stars on the URs involving the moment of inertia or relevant for GW analysis from binary systems is subject to future investigation.It is of great interest to the GW community to calculate the quasi-normal modes of hybrid stars relevant for third-generation detectors (Hild et al. 2011;Hall 2022) or the planned GW mission NEMO (Ackley et al. 2020).Interestingly, Sotani & Kojo (2023) find that the fundamental  -mode frequencies with Quark-hadron crossover EOS basically are smaller and the 1st pressure p1-mode frequencies with QHC EOS are larger than those with hadronic EOS and, moreover, are able to distinguish between these two possibilities using the so-called universal relations.It was suggested that oscillation modes in NSs could provide smoking-gun evidence for the nature of the mixed phase.A few recent studies also attempted to study the effects of the mixed phase on NS oscillation modes such as g-modes Constantinou et al. (2023) and w-modes (Ranea-Sandoval et al. 2022b,a), however, they did not comment on their detectability from GW observations.The results of our investigation demonstrate that for future  -mode detections from isolated NSs, whether or not one may be able to differentiate between sharp and smooth phase transitions depends on the sharpness of the phase transition as well as uncertainties in the observations and the universal relations.Improved universal relations and high-precision measurements may provide hope to be able to comment conclusively about the existence of a mixed phase in the NS interior.9) are also displayed.(fig.A1b) Same as fig.A1a, but with Λ 10 as the independent variable and the 'Fit' relations correspond to eq. ( 10).

Figure 1 .
Figure 1.We display (a) hybrid EOSs corresponding to ACB4 parametrization with different Δ  values (b) the corresponding  −  relations.The dotted lines in 1b and 2 show the unstable region with   ⩽ 0. Horizontal bands correspond to masses  = 2.072 +0.067 −0.066  ⊙ of PSR J0740+6620 and  = 2.01 +0.04 −0.04  ⊙ of PSR J0348+0432.In fig.1b, the 90% contours for PSR J0740+6620 corresponding toMiller et al., 2021(Miller et al. 2021) is shown in grey color while the contour of  −  measurement for PSR J0740+6620 corresponding toRiley et al., 2021(Riley et al. 2021) is shown in brown.The  −  estimates of the two companion neutron stars of the merger event GW170817 are shown by the shaded area labeled with GW170817 M1 (M2) in (b).
, a few hybrid EOS models from Ayriyan et al. (2021); Paschalidis et al. (2018); Alvarez-Castillo (2021), and the EOSs shown in figs.1a and 9a.All the hybrid EOS models considered in this work are displayed in fig.B1 of appendix B The URs involving stellar compactness ( = /) and mass-scaled  -mode angular frequency () are given in eqs.(

Figure 3 .Figure 4 .
Figure 3.We show the variation of  -mode (a) frequency and (b) damping time as a function of stellar mass corresponding to stable  −  configurations of the EOSs presented in fig.1a.

Figure 7 .
Figure 7. (a) Marginalised corner distribution of  and  recovered with different GW network configurations (blue for A+ and orange for ET).Red lines mark the injected values.The injected values for  = 1603.2Hz and  = 0.23 s corresponds to  = 1.75 ⊙ and EOS with Δ  = 0%.(b) Figure showing the joint distribution of  and  obtained by using the URs eqs.(7) and (8) from the recovered (  , ) posterior.The injected values of  = 1.75 ⊙ and  = 14.10 km are shown using the red lines.In both figures, the title shows the median and symmetric 90% credible intervals of the parameters.

Figure 13 .
Figure 13.Same as figure fig.12 but for the dynamical tidal parameter   .

Figure A2 .Figure A3 .Figure B2 .
Figure A2.(Upper panel)  − Love relation resulting from consideration of a wide range of hybrid EOSs along with the different  − Love relations from the literature includes URs from Chan et al. (Chan et al. 2014), Sotani & Kumar (Sotani & Kumar 2021), and Pradhan et al., (Pradhan et al. 2023a).(lower panel) The relative error on the Re(  ) for the corresponding fit relation has also been displayed.Additionally, different URs and the resulting error due to different URs on the   of a representative hybrid EOS ACB4,Δ  = 5%, are also shown.

Table 1 .
The values of the polytropic index Γ  , the pre factor   and the

Table 2 .
Fit parameters for the URs (7) and (8) obtained in this work.We provide the quadratic fit for Re(M) and a fit relation (8) for Im(M).
and fig.3bthe fundamental mode (  -mode) frequency (  ) and damping time () respectively as a function of mass (of the stable branch) corresponding to the EOSs in fig.1a.From fig.3, one can conclude that there is a sudden increase in the frequency (decrease in the damping time) at the onset of the

Table 3 .
Ranea-Sandoval et al. (2022b)18).From fig.4, it is evident that different URs are in good agreement.Additionally, we show the errors on the  -mode characteristics resulting from Fit parameters for real (  ) and imaginary (  ) of the complex frequency , related to  = log  1.4⊙ Λ through the URs eq.(9).The coefficients corresponding to   and   are denoted with subscript  and , respectively.theproposedURsfor the hybrid Equation of State (EOS) considered, as well as for a representative hybrid EOS, ACB4 with Δ = 5%.We observe better agreement for Re() among the proposed URs and those given inTsui & Leung (2005).The deviation of the fit relation ofLioutas & Stergioulas (2018)from General Relativity (GR) can be attributed to the approximations considered for solving the  -mode damping times in their work.For the Im(), we notice a comparatively lower error with the proposed URs compared to the relation given inTsui & Leung (2005).For the URs involving the love number, e.g., the  − Love relation discussed in the appendix A, although the different URs introduce different uncertainties in different ranges of Λ, we notice that the errors are minimal, irrespective of the UR used.As we focus on the distinguishability of hybrid starsKumar (2021), involving the tidal deformability (Λ) scaled  -mode parameters and later considered inRanea-Sandoval et al. (2022b)for -modes.The new URs involving Λ can be written as, and to avoid clutter/confusion in fig. 5 and fig.6, we only show the estimations of Mass-Radius ( − ) recovered with the newly obtained URs.The uncertainties in the  −  estimation from future  -mode observations discussed later in section 3.2 result from the measurement uncertainties on the mode frequency and damping time.Therefore, fixing the URs to a different one will not change the uncertainties in the  −  measurement, as discussed in fig.8.However, a change in UR might lead to a bias or shift in the  −  measurement.Recently, a new class of URs has been proposed in Sotani & log  , Λ =  , +  ,  +  2 (9) where,  = log Λ 1.4 ⊙ with  1.4 ⊙ = /1.4⊙ ,   = Re() = 2f and   = Im() = 1/ .Similarly, other URs exist for  , , Λ as a function of ln  1.4 ⊙ Λ are displayed in fig.A1a of appendix A and the fit parameters are tabulated in table 3. We display the URs among  , with ln [Λ 10km ] fig.A1b of appendix A and tabulate the fit parameters in table 4. The universal relations (URs) discussed in this study incorporate a broad spectrum of hybrid equations of state (EOSs) derived from several sources, including Ayriyan et al.