{'fraction_non_alphanumeric': 0.03554856867567636, 'fraction_numerical': 0.020411268034179364, 'mean_word_length': 4.349247253492472, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Twin stars as probes of the nuclear equation of state: effects of rotation through the PSR J0952-0607 pulsar and constraints via the tidal deformability from the GW170817 event

In agreement with the constantly increasing gravitational wave events, new aspects of the internal structure of compact stars can be considered. A scenario in which a first order transition takes place inside these stars is of particular interest as it can lead, under conditions, to a third gravitationally stable branch (besides white dwarfs and neutron stars), the twin stars. The new branch yields stars with the same mass as normal compact stars but quite different radii. In the present work, we focus on hybrid stars undergone a hadron to quark phase transition near their core and how this new stable configuration arises. Emphasis is to be given on the aspects of the phase transition and its parametrization in two different ways, namely with Maxwell and Gibbs construction. We systematically study the gravitational mass, the radius, and the tidal deformability, and we compare them with the predictions of the recent observation by LIGO/VIRGO collaboration, the GW170817 event, along with the mass and radius limits, suggesting possible robust constraints. Moreover, we extent the study in order to include rotation effects on the twin stars configurations. The recent discovery of the fast rotating supermassive pulsar PSR J0952-0607 triggered the effort to constrain the equation of state and moreover to examine possible predictions related to the phase transition in dense nuclear matter. We pay special attention to relate the PSR J0952-0607 pulsar properties with the twin stars predictions and mainly to explore the possibility that the existence of such a massive object would rule out the existence of twin stars. Finally, we discuss the constraints on the radius and mass of the recently observed compact object within the supernova remnant HESS J1731-347. The estimations implies that this object is either the lightest neutron star known, or a star with a more exotic equation of state.


I. INTRODUCTION
Compact stars yield the most prominent natural laboratories for the study of exotic forms of matter [1][2][3][4]. Recently discovered pulsars alongside with gravitational waves detection, such as GW170817, have revealed new aspects of the internal structure of these stars, mainly in terms of their composition [5][6][7]. Whilst the equation of state (EoS) of nuclear matter is well established up to nuclear saturation density, one encounters the challenge of describing matter in fairly higher densities realized in the interior of these stars. At these densities the type of matter is yet to be determined and in turn, the construction of stellar models to agree with the aforementioned observations is still an open issue. Possible candidates are pure neutron stars composed by hadrons, strange quark stars composed of deconfined quarks, and hybrid stars composed by hadronic outer shells and cores of deconfined quarks. In the present article the latter case is claimed.
Stars of this branch are expected to have masses in the same range as normal neutron stars, yet fairly smaller radii. The existence of such stars is a strong indication that a hadronquark phase transition (HQPT) is a physical reality, a result of utmost importance, especially in the study of dense matter physics [8][9][10]. A study of the HQPT is presented here, examining the conditions under which the twin star configu-ration arises. The compatibility with the mass and radius constraints, as they are formed through up to date observations, is also considered.
The idea of a third family of compact stars and in particular, the connection with the possibility to be a signature of a strong phase transition in the interior of the star, was firstly introduced by Gerlach [11]. Later on, Kämpfer worked also on this issue [12,13]. Glendenning and Kettner introduced the term "twins" in their paper [14], while at the same time Schertler et al. [15] worked out that idea in detail. However, in all the previous studies, the maximum mass was approximately at the canonical binary pulsar mass 1.4 M . The revival of the idea of the twin stars started a few years later by Blaschke et al. [16,17]. Specifically, in the mentioned papers is suggested that high-mass twin stars, once detected by simultaneous mass and radius measurements, could provide the evidence for a strong first-order phase transition in cold matter, which then would imply the existence of at least one critical endpoint in the quantum chromodynamics phase diagram. Moreover, in the same work are also presented, for the first time, examples of EoSs that would not only provide twin solutions, but also fulfill the constraint on the maximum mass from the existence of pulsars as heavy as 2 M [16,17]. The previous idea was elaborated by Benic et al. [18] (see also Ref. [19]). A systematic Bayesian analysis of the new twin star EoS with observational constraints was presented in Ref. [20]. Finally, the analysis of the robustness of twin solutions against changing the Maxwell to a mixed phase construction, and the formation of structures in the mixed phase due to the interplay of the surface tension and Coulomb interaction effects, have been also considered respectively in Refs. [21,22]. Some recent works dedicated to the study of phase transition in the interior of the neutron stars and possible existence of twin stars, are included in Refs. .
Pulsar PSR J0952-0607 was firstly discovered by Bassa et al. [65] with a frequency of f = 709 Hz, making it the fastest known spinning pulsar in the disk of Milky Way. Very recently, Romani et al. [66] discovered that PSR J0952-0607 has a mass of M = 2.35 ± 0.17 M , the largest well measured mass found up to date. The above discoveries triggered the research into the constraints of the EoS of dense nuclear matter and is likely to revise many of the theoretical predictions concerning basic properties of neutron stars (for a recent study see Ref. [67]). One of them is the existence of twin stars, which is the main subject of the current study. In particular, the existence of a supermassive neutron star may rule out twin stars (see Ref. [27]). In this sense, the observation of such supermassive stars is decisive for checking the reliability of the corresponding theoretical predictions and therefore, a systematically and careful analysis should be consider. On the other hand, theoretical calculations for the maximum possible mass of neutron stars should be combined with recent observations concerning the radius and tidal deformability of stars with mass close to 1.4 M . In general, a nuclear model can be judged by its flexibility in being able to predict the maximum possible mass while correctly estimating the observational radial and tidal deformability values.
Moreover, one of the motivations of the present work is to examine in a more systematic way the applications of the two main formulation of the phase transition, that is the Maxwell (MC) and Gibbs (GC) constructions. These formulations are quite different, as the former imposes an energy jump between the two phases while the latter implements a smooth transition between the phases. It is worth mentioning that although the MC has been extensively applied in the literature, the GC have not received the attention it deserves (for a recent relevant study see Ref. [10]). In any case, the comparison of their predictions may offer useful physical insight.
Rotating neutron stars with exotic degrees of freedom in their cores have been already studied in limited papers by Banik et al. [59], Bhattacharyya et al. [42], Haensel et al. [62], and Bozzola et al. [63]. In the present study, in addition to the mentioned similar studies, the effects of rapid rotation on the twin stars scenario have been explored. The reason for this kind of study is twofold: on one hand we intend to examine to what extent the high rotation can differentiate the representation of the two stability branches of a static neutron star, and secondly to consider the possible cases where the existence of the twin stars may have arisen as a result of the high rotational frequency, where in the corresponding static case such a scenario is not foreseen. According to our knowledge this peculiar case has not been explored in details in the relevant bibliography, at least up to now. There is only an interesting mention in Ref. [59]. Nonetheless, future study is in progress, in order to further clarify the effects of rotation on twin stars problem.
Finally, we employ the constraints on mass and radius of the recently observed compact object within the supernova remnant HESS J1731-347 [68]. The corresponding estimations are M = 0.77 +0. 20 −0.17 M and R = 10.4 +0.86 −0.78 km, respectively, while the observations were carried out with the help of modelling of the X-ray spectrum and a robust distance estimate from Gaia observations [68]. According to the authors' guess, the above estimates imply that this object is either the lightest neutron star known, or a star with a more exotic equation of state. In any case, it is worth considering to what extent this compact object is compatible with the hybrid model and thus with the twin stars theory. The paper is organized as follows. In Sec. II we present the basic formalism of the hadron to quark phase transition, while in Sec. III we provide the tidal deformability. Sec. IV is dedicated to the presentation and discussion of the results of the present study and, finally, the concluding remarks are given in Sec. V.

A. Hybrid EoS and 1st-Order Phase Transition
The theoretical framework throughout this work is based on the two prescriptions for matching the low-intermediate density hadronic EoS with the one describing the free quark matter, the Maxwell and the Gibbs constructions. In particular, while both methods have a phenomenological origin, they mimic, in a way, the theoretical approach of phase transition. In this sense, they should be considered as a useful mathematical tool rather than an exact description of the phase transition from hadronic to free quark matter. Furthermore, although the first method (MC) appears to be more applicable, the recent observations of gravitational waves related to the binary neutron star merger GW170817 event, along with the flexibility of the second method (GC), have made the latter case a more compelling candidate [10]. Nevertheless, both methods have been used in the literature, leading under certain circumstances to the appearance of twin stars.
We note here that as the quark phase has a strong dependence on the speed of sound, for both configurations, in order to achieve the observed neutron star masses, at the selected transition densities, the speed of sound is set equal to the speed of light, having as a result the maximally stiff EoS.

Maxwell construction
The first case is the well-known MC case, which exhibits a sharp transition at the boundary and makes it hard for charged clusters of quarks to form in the hadronic matter. Of course, as it is expected, an abrupt phase change of this kind, yields a discontinuity in at least one physical quantity, with the most obvious one, the energy density, claiming the form [8][9][10] E(P ) = E HADRON (P ), P ≤ P tr E(P tr ) + ∆E + c −2 s (P − P tr ), P ≥ P tr . (1) In the above formula E(P ) denotes the energy density, P the pressure, c s = ∂P/∂E the speed of sound in units of speed of light, and ∆E the magnitude of the energy density jump at the transition point. During the quark phase the numerical value we assign for c s is equal to c s = 1, the maximum allowed value that is consistent with causality. That way we also ensure the stiffest EoS case and the greatest possible range of M-R relations. Moreover, P tr expresses the pressure that corresponds to the baryon density at phase transition point, n tr . It should be clear that the first line in Eq. (1) refers to the hadronic phase while the second one refers to the quark phase, seen as the 1st-order Taylor expansion of the energy density around the transition pressure plus the ∆E term. This is the so-called constant speed of sound parametrization (CSS) [8,10,[69][70][71]. The described process requires not just the pressure but also the rest of the intensive thermodynamic quantities, meaning temperature T (set equal to 0) and baryonic chemical potential µ B , of both phases to have the same value at the phase transition. The same is not true for the electric chemical potential µ Q , as there is a jump at the interface between the two phases [25]. Local charge neutrality conditions for both regions separately, must be imposed, in order to ensure β-equilibrium state.

Gibbs construction
In the case of a non-sharp HQPT, a finite region to embody the transition is implied. The mixed phase of this transition is composed, as its name suggests, by intermittent domains of pure hadronic and quark phases [10,25]. The Gibbs phase transition rule regarding the equality of the pressure of the two components (Hadron -Intermediate, Intermediate -Quark) is established here. Constraints requiring global charge neutrality and baryon conservation number are imposed throughout the whole process. This is quite important, as now both phases can be oppositely charged, as long as the mixed phase remains neutrally charged. Contrary to the MC case, where the pressure remains constant in the transition interval, in the GC case, the pressure increases with increasing baryon density, while also no discontinuities in the energy density appear, giving rise to the profile [10] −1 , and c s = 1 (maximally stiff EoS). The energy density is denoted by E(P ), the pressure by P , the speed of sound by c s , while a m , K m and the polytropic index Γ m are constants, with the former two evaluated, by requiring continuity of P and E at the transition points. As has already been explained, an energy density jump is not present here, opposed to the MC, but the increase in the quantity is assigned during the intermediate part of the formula. Subscript "css" denotes the corresponding quantity at the start of the quark phase.
The relation between the pressure and baryon density for the intermediate section in Eq. (2) is given through the wellknown, simple polytropic formula P (n) = K m n Γm . While the values of a m and K m depend on the transition point, the polytropic index Γ m in the GC is taken to be constant throughout the whole work and equal to 1.03. A value so close to unity is justified by taking into account the uncertainty related to how soft the EoS of the mixed phase is and also differentiating it as much as possible from the respective value of the MC, where according to the relevant discussion above, we have Γ m = 0 (see also Ref. [10]).

B. Seidov Criterion
The phase transition from the hadronic matter to the quark matter described in the previous section is set to produce an instability in the star when a certain criterion is met. This critical value, was derived by Seidov [72], based on the work of Lighthill [73], by performing a linear perturbative expansion in the size of the quark matter core and reads as In order to have a third family of compact objects appear in the M-R diagram, the aforementioned instability has to be satisfied. The theory dictates that as long as the mass M is an increasing function of the central pressure P c , the star will remain stable. As the central pressure increases, it reaches a certain point where its value becomes equal to the transition pressure P tr , leading to the formation of the quark matter core. Depending on whether the Seidov criterion for the energy density jump is surpassed or not, there are four possible outcomes, with only two of them producing a new type of compact object [3,8]. It should be noted here that the Seidov criterion is meaningful only in the presence of a sharp discontinuity in the energy density profile, which appears in the MC method as opposed to the GC method. In order to be able to compare the GC with the MC, where an energy jump appears in the form of ∆E cr , we define the corresponding energy increase in the GC as the quantity with E HADRON (P tr ) and P tr representing the respective values at the transition from the hadron phase to the mixed phase. The two possible outcomes for producing a new type of compact objects, are: • Case 1: If ∆E is taken to be much higher than the Seidov limit (∆E >> ∆E cr ), then the instability appears immediately after the formation of the quark matter core. The energy density jump that is present on the core bypasses the pressure that the quark matter exerts on the rest of the star. On the M-R diagram, the transition point takes the form of a cusp, meaning the sign of dM/dR is flipped. This case does not result to the formation of a stable hybrid star. In our work, just for reference, this particular scenario is shown for ∆E = ∆E cr + 100 MeV fm −3 and ∆E = ∆E cr + 200 MeV fm −3 of case n tr = 0.5 fm −3 in Figure 3.
• Case 2: The requirement of reaching the Seidov criterion is important here, since the exact value of ∆E should be higher that the value of Eq. (3) (∆E > ∆E cr ), but not high enough to reach Case 1 scenario. This case seems to be the most interesting and the one that will be one of the main focuses of our work, as it describes the appearance of a third family of compact objects. The cusp is also present here during the transition point, leading to an unstable disconnected area in the graph, which after a while flips around again, signalling the appearance of a stable hybrid star with a quark matter core. For completeness of the theory, we have also included the special case where ∆E = ∆E cr as seen in cases (a) of Figure 1 and Figure 3.

C. Violation of Seidov Criterion
At this section, we refer to the two remaining scenarios that describe the cases where the Seidov limit is not reached. They both include the hybrid branch connected to the nuclear star branch that may lead to the appearance of a third family of compact objects.
• Case 3: This is the first of the final two cases where the Seidov limit is not reached, known as both connected and disconnected hybrid branch (for more details see chapter 9 of Ref. [3] and Fig.3 of Ref. [8]). The reason is that, after the P c surpasses the transition value P tr and the quark matter core is formed, the energy density jump ∆E is not high enough to cause an instability at the moment (∆E < ∆E cr ), meaning that the quark matter is described by a curve on the map that is connected to the hadronic matter branch. The curve continues up to a point, after which an unstable region appears leading again to a stable hybrid branch, just like Case 2 scenario, albeit shorter in length. The existence of a hybrid star branch, as shown in Figure 12 for n tr = 0.38 fm −3 , proves to be an interesting topic of study. Additionally, since it is the only case that exhibits an intermediate instability appearing between the two quark matter sectors, it appears to be also a great candidate for future work.
• Case 4: Sharing a lot of similarities with the previous case, but overall a different final result, this case is also derived by not reaching the Seidov limit. The hybrid star branch connected to the nuclear branch is also present here with its length to be inversely proportional to the ratio of ∆E/E tr . It is shown that higher values of the above quantity, lead to a decrease in the size of the hybrid branch that spans the phase diagram. Immediately after a maximum mass value is reached, a sudden discontinuity appears that continues forth without the appearance of a third family of compact objects [8]. In our work no cases of this particular scenario have been derived.

III. TIDAL DEFORMABILITY
A very important source for the gravitational wave detectors are the gravitational waves from the late phase of the inspiral of a binary neutron star system, before the merger [74][75][76]. This kind of source leads to the measurement of various properties of neutron stars. In the inspiral phase the tidal effects can be detected [75].
The k 2 parameter, also known as tidal Love numbeer, depends on the equation of state and describes the response of a neutron star to the tidal field E ij [75]. The exact relation is given below where R is the neutron star radius and λ = 2R 5 k 2 /3G is the tidal deformability. The tidal Love number k 2 is given by [75,76] where β = GM/Rc 2 is the compactness of a neutron star. The paramter y R is determined by solving numerically the following differential equation with the initial condition y(0) = 2 [77]. F (r) and Q(r) are functionals of the energy density E(r), pressure P (r), and mass M (r) defined as [74] and Eq. (7) must be solved numerically and self consistently with the Tolman -Oppenheimer -Volkoff (TOV) equations under the following boundary conditions: y(0) = 2, P (0) = P c (P c denotes the central pressure), and M (0) = 0 [74,76]. From the numerical solution of TOV equations, the mass M and radius R of the neutron star can be computed, while the corresponding solution of the differential Eq. (7) provides the value of y R = y(R). The last parameter along with the quantity β are the basic ingredients of the tidal Love number k 2 .
The chirp mass M c of a binary neutron stars system is a well measured quantity by the gravitational wave detectors [7]. Its relation is given below where m 1 is the mass of the heavier component star and m 2 is the lighter's one. Hence, the binary mass ratio q = m 2 /m 1 lies within the range 0 < q ≤ 1.
Additionally, another quantity that is well measured is the effective tidal deformabilityΛ which is given by [7] Λ = 16 13 where Λ i is the dimensionless tidal deformability [7] The effective tidal deformabilityΛ is one of the main quantities that can be well measured by the detection of the corresponding gravitation waves.

A. Mass vs Radius diagram
In our study, as we have already mentioned, we used two constructions, a) the MC, and b) the GC, and two different EoSs, a) the MDI+APR1 EoS [78], and b) the GRDF-DD2 EoS [79]. Also, we focused on the following values of Especially, as we move to higher values of ∆E, the EoSs become softer (smaller higher masses and smaller radii). The purple horizontal shaded regions correspond to observational data from pulsars, while the light purple contour shaded region indicates the observation of the GW170817 event [7]. The dashed part of curves indicates their unstable region. Another remark is that between the two constructions, MC and GC, the first one provides branches of EoSs in a lower mass region than the last one. Therefore, the MC is more informative in our case of study, since the GW170817 event contains low values of component masses. Also, as one can observe, as we move to higher values of n tr the EoSs become softer.
Moreover, in Figures 1-4, the constraints on mass and radius of the recently observed remnant HESS J1731-347 have also been included [68]. Obviously, the prediction of the above constraints requires the use of a larger jump in the energy density for both constructions. In particular, we found that maybe this star is a hybrid star, belonging to the second stable branch. It is also remarkable that these predictions are compatible with the prediction of the maximum observed masses. It is worth mentioning that we do not claim that the remnant HESS J1731-347 is definitely a hybrid star. However, according to the present study, this case is more favorable than to be a neutron star. Further similar observations are needed in order to clarify this issue.
B. Tidal Deformability Figure 5 presents the tidal parameters k 2 and λ of a single neutron star as a function of its mass, for the MDI+APR1 EoS (MC) and for the cases with ∆E = ∆E cr + [100, 200] MeV fm −3 . The effect of the different ∆E and n tr can be observed, as they lead to distinct subbranches. In general, the higher value of ∆E lead to softer EoSs across the same bifurcation characterized by the value of transition density n tr .
In Figures 6-7 we present the Λ 1 − Λ 2 space for each EoS and configuration using the observational data of the GW170817 event (orange shaded region) [7]. To be more specific, we considered the three following combinations: Hybrid-Hybrid binary star system (HS-HS), Hybrid- for the majority of the cases of those hypothetical binary star systems, in some cases we modified and restricted the mass range for computational reasons.
In all cases, the HS-HS lead to smaller values of Λ, therefore to softer EoS, in accordance to the observational data of GW170817. The HS-NS case of a binary star system lies also inside the shaded region provided from LIGO, but in intermediate values of Λ. The NS-NS is the case where the two component stars correspond to the original EoS that we used in each case (the solid black line in M-R diagrams). As one can observe, the MC provides more cases for the HS-HS scenario compared to the GC, for both EoSs that we used in our study. In addition, the MDI+APR1 EoS for the NS-NS case lies inside the estimation of GW170817 nevertheless the combination that we considered, on the contrary of the GRDF-DD2 EoS, in which the NS-NS cases lie outside the estimated region.
In Figures 8-9 we show theΛ − q relation, by using the upper limit onΛ, provided by the GW170817 event [7]. The curves and colors are the similar to those of Λ 1 −Λ 2 diagrams. As one can observe, the HS-HS case for both EoSs and constructions lead to much lower values ofΛ. Hence, a possible lower limit onΛ would be very useful to restrict the lower values ofΛ, leading to constraints at least on the HS-HS case. Moreover, as we move to higher values of ∆E, all the curves for both HS-HS and HS-NS scenarios are shifted to lower values ofΛ, meaning that the increment of ∆E has as a result softer EoSs. This behavior was expected if we recall the effect of higher values of ∆E on the EoSs (see M-R diagrams). Between the two constructions, and for both EoSs, the MC provides higher number of HS-HS cases compared to the GC. Also, across the same EoS and n tr , the MC leads to a bit of lower values ofΛ compared to the GC. Therefore, for low mass events, such as the GW170817, MC is more suitable. The need for a lower limit onΛ that we described before, led us to the exploit of the constrained value of the dimen- sionless tidal deformability for a single 1.4M neutron star, derived by the study of the GW170817 event. In Figure 10 we show the relation between Λ 1.4 and ∆E. We notice that we used only those values of n tr that provide a separate branch from the original EoS. As one can observe, as we move from the ∆E cr (and ∆E G for the GC) to higher values of ∆E the variation between the marks decreases. In addition, all the marks that correspond to the GC (diamonds) predict higher values of Λ 1.4 compared to the MC, for both EoSs. We notice also that only for the ∆E = ∆E cr + 200 MeV fm −3 and ∆E = ∆E G + 200 MeV fm −3 there is a violation of the ac-cepted region. Moreover, in Figure 10(b), the square points indicate the ∆E = ∆E cr − 5 MeV fm −3 case under MC and for the GRDF-DD2 EoS. As one can observe, as we move to lower values compared to ∆E cr , the variation increases, which is in accordance with the behavior that we described before.
Furthermore, in order to shed more light on which specific cases should be excluded, we studied the relation between Λ 1.4 and n tr . In Figure 11 we present the aforementioned relation. As a first remark, the MC provides in all cases one more value of n tr compared to the GC. The GC shifts the curves (dashed) to higher values compared to the relevant curves of The left panel of the Figure 11 corresponds to the MDI+APR1 EoS. The curves that correspond to the ∆E cr and ∆E = ∆E cr + 100 MeV fm −3 lie inside the estimated region.
The same holds for ∆E G , and ∆E = ∆E G + 100 MeV fm −3 . On the other hand, the solid green curve, which corresponds to the MC with ∆E = ∆E cr + 200 MeV fm −3 is excluded. But if we apply the GC, the curve is shifted upwards (dashed green curve), with only a part being outside of the estimated region. Especially, this part is between n tr = 0.276 fm −3 and n tr = 0.353 fm −3 . Hence, not only the kind of construction we choose has a significant role, but the exact value of the transition density affects the final output. Therefore, a further understanding and possible constraints on the transition density n tr are necessary to shed more light on the twin star hypothesis.
The right panel of the Figure 11 corresponds to the GRDF-DD2 EoS. The curves that correspond to the ∆E cr and ∆E = ∆E cr + 100 MeV fm −3 lie inside the estimated region. Only the green solid curve which corresponds to the MC with ∆E = ∆E cr + 200 MeV fm −3 lies outside up to n tr = 0.316 fm −3 , meaning that above this value even this EoS could be acceptable. All the curves that correspond to the GC lie inside the estimated region. As we mentioned above, the construction and the transition density affect importantly the behavior of the curves. We notice that the purple line indicates the ∆E = ∆E cr − 5 MeV fm −3 case under MC. As we have already noticed, as we move to lower values of ∆E compared to ∆E cr , the EoS becomes stiffer, therefore in this diagram the curve is shifted sligtly to higher values of Λ 1.4 .
In Figure 12  corresponds to a prediction under the Seidov limit. It is interesting that this case predicts the existence of a second stable branch and consequently, the existence of twin star (and indeed within the predictions of GW170817 event [7]). However, this case cannot predict the existence of the compact object of the remnant HESS J1731-347 [68]. These results confirm those already presented in Figures 1-4, that is, hybrid EoSs with a large energy density gap are preferable for explaining the existence of the above compact object.

C. Rotational frequency at 709 Hz and the PSR J0952-0607 pulsar
The simultaneous observation of mass and frequency of the PSR J0952-0607 pulsar [66] opens a new window in studying the rotating configuration for neutron stars and the possible existence of rotating twin stars. In Figures 13  and 15 we present three cases of transition density, n tr = [0.2, 0.25, 0.3] fm −3 , for both the non-rotating and rotating at 709 Hz configuration, with MC and GC respectively. Furthermore, we extended the study in a region of the energy jump, starting from the ∆E cr for MC and ∆E G for GC and reaching values up to +150 MeV fm −3 . In particular, while in Figure 13(a) there are numerous cases that exist either in the mass limit of the rotating neutron star or in the non-rotating one, there are only three cases where both configurations meet the mass limits simultaneously [66,[81][82][83]. Specifically, (a) the n tr = 0.25 fm −3 and ∆E = ∆E cr + 150 MeV fm −3 , (b) the n tr = 0.3 fm −3 and ∆E = ∆E cr + 60 MeV fm −3 , and (c) the n tr = 0.3 fm −3 and ∆E = ∆E cr + 90 MeV fm −3 . In addition, the mentioned cases also fulfil the non-rotating ra-  [68], GW170817 event [7], PSR J1614-2230 [81], PSR J0348+0432 [82], PSR J0740+6620 [83], and PSR J0952-0607 [66] pulsar observations for possible maximum mass. The inset plot indicates the case with ntr = 0.38 fm −3 , in which the square point represents the phase transition. dius limit [84], as shown in Figure 13   The dash-dotted lines represent the non-rotating configuration, while the solid lines correspond to the 709 Hz configuration. To guide the eye, the unstable region is marked with the dashed lines. The shaded regions from bottom to top represent the PSR J1614-2230 [81], PSR J0348+0432 [82], PSR J0740+6620 [83], and PSR J0952-0607 [66] pulsar observations for possible maximum mass. GW170817 event [84] for the 1.4 M configuration. In the latter, we also present, for completeness and comparison, the corresponding radius at 709 Hz. However, in Figure 13(a), there are numerous cases that can describe solely the 709 Hz configuration. It has to be noted that the mentioned cases are eligible for the description of the PSR J0952-0607 [66] pulsar, since we do not know the corresponding non-rotating limit for the mass. On the contrary, in Figure 15(a) there are only a few cases that represent twin stars, which nearly all of them (except one) meet the current limits of the gravitational mass only at the non-rotating configuration. In addition, their radius fulfils the limits for the 1.4 M configuration, imposed from the GW170817 event, indicated in Figure 15(b). The exception case can only describe the 709 Hz. Moreover, we also present the corresponding Figures 14 and 16 where we highlight the transition density over the energy jump.
The GC configuration presented an effect where while the non-rotating configuration is not able to produce twin stars, that is not the case at 709 Hz [59]. Specifically, Figure 17 displays the representative cases: (a) n tr = 0.25 fm −3 and ∆E = ∆E G + 120 MeV fm −3 and (b) n tr = 0.3 fm −3 and ∆E = ∆E G + 60 MeV fm −3 . The last case, n tr = 0.25 fm −3 and ∆E = ∆E G + 90 MeV fm −3 , is not presented for the clearness of the figure. Moreover, for clarity, we indicate in the inner figure, the region of an EoS at 709 Hz, where the twin star branch appears. Also, we have marked the twin star that corresponds to the first maximum mass configuration along with the relevant radii. This effect added the reported cases in the list where the rotation meets the mass limit of the PSR J0952-0607 pulsar.
The MC and GC configurations are both suitable methods to describe twin stars, where the second one, using similar parametrization, is more elaborated. The key differences between them are that (a) the MC provides higher gravitational masses than the GC, and (b) the MC rises fairly more twin stars than the GC. A possible explanation for both of them may lie to the existence of the sharp transition and the energy jump applicable on the MC, where in the GC, both of them are extinct. However, the GC can provide twin stars in the rotating configuration without their appearance in the non-rotating scenario. This peculiar effect may have its origin at the smooth phase transition that the GC uses, alongside with the transition density, the rotational frequency, and the hadronic EoS. Nevertheless, a thorough study for the latter effect should be done and will be the subject of a forthcoming paper. The analysis of Figures 13-17 enrich the scenario of twin stars existence as they can describe the observed neutron star masses both at non-rotating and rotating configurations.
For the numerical integration of the equilibrium equations at the rotating configuration, we used the publicly available numerical code nrotstar from the C++ Lorene/Nrotstar library [85].

V. CONCLUDING REMARKS
The existence of two stable branches in a neutron star configuration, which lead to the concept of twin stars, may reveal the scenario of phase transition in the interior of dense nuclear matter. Therefore, it is extremely important to study the above hypothesis both from theoretical and and observant point of view. In the present work a systematic study, concerning the twin star hypothesis, has been carried out and the main conclusions could be summarized as follows: 1. There could be hybrid or twin star branches that can describe the current observations of binary neutron star systems, such as the GW170817 event. Not only the transition density has an important role on the possibility of the existence of a twin star branch, but also the ∆E affects this scenario. Especially, as we move to higher values of ∆E, the twin star branch is more easily to exist even for the lowest values of n tr that we used in our study.
2. Across the same EoS, the kind of construction that we use affects the behavior of the curves. In general, the MC offers more hybrid and twin star branches in lower neutron star masses, compared to the GC. Moreover, the MC leads to a softer behavior of the EoS. Therefore, binary neutron star systems that contain lower component masses (otherwise, lower chirp mass), suit better to the MC. Also, the GC shifts the curves to higher values of Λ, as one can see from the Λ 1.4 − n tr diagram. This fact has as a consequence that some cases which were to be excluded by using the MC, could lie (by using the GC) as a whole or partly inside the estimated area.
3. In our study we used the observational value ofΛ using as a reference the GW170817 event. This upper limit was not helpful to constrain further the cases that we studied. On the contrary, because of the fact that the HS-HS binary system scenario has a very soft behavior on the EoS, a possible lower limit onΛ could be very helpful to shed more light on the problem, from another point of view. Moreover, a further constrain on the radius (e.g. the radius of a 1.4M ) could be informative in order to impose constraints on the very soft cases.
4. The observation of the PSR J0952-0607 pulsar, where both mass and frequency are evident, does not exclude the existence of twin stars. In fact, both MC and GC are able to provide numerous configurations that can describe either solely the rotating configuration, solely the non-rotating configuration, or both of them simultaneously. The latter, since we do not have the corresponding non-rotating mass of the PSR J0952-0607 pulsar, provide us with strong evidence that twin stars can exist as both static and rotating.
5. The rotational frequency, in some cases of the GC, is able to introduce twin stars from normal stars. In particular, as the rotating twin star loses its rotational frequency, the neutron star moves to higher energies and stable configurations, where the EoS is more stiff and extincts the phenomenon of twin stars. However, an analysis using fixed baryon mass sequences would provide more information for the effect. Nevertheless, we speculate that this effect due to the stiffening of the EoS, comes as a consequence to the GC, where a smooth phase transition is employed, in contrast to the MC, where an energy jump appears.
6. The recent observation of a compact object with the lowest observable mass M = 0.77 +0.20 −0.17 M , within the supernova remnant HESS J1731-347, opens a new window to impose constraints on the EoS of dense nuclear matter. In the present study, our predictions reinforce the estimation of the authors in Ref. [68] that the observed compact object is rather a hybrid star with exotic core than the lightest neutron star ever known. This result strengthens the hypothesis of hybrid stars and therefore of the existence of twin stars.
The above conclusions can be summarized as follows: Further systematic theoretical research is required to clarify the role of the method that describes the phase transition in dense nuclear matter, as well as its particular characteristics (density transition, energy, etc.). In addition, more relevant observations from binary neutron star merger (and not only) are necessary to be able to check the plausibility of the theoretical predictions. Possibly, in this case we will be able to confirm, with enough confidence, the existence of twin stars but even more importantly, to confirm, not only qualitatively but also quantitatively, the phenomenon of phase transition in dense nuclear matter.