Probing the nuclear equation of state from the existence of a $\sim 2.6~M_{\odot}$ neutron star: the GW190814 puzzle

On August 14, 2019, the LIGO/Virgo collaboration observed a compact object with mass $\sim 2.59_{-0.09}^{+0.08}~M_{\odot}$, as a component of a system where the main companion was a black hole with mass $\sim 23~M_{\odot}$. A scientific debate initiated concerning the identification of the low mass component, as it falls into the neutron star - black hole mass gap. The understanding of the nature of GW190814 event will offer rich information concerning open issues, the speed of sound and the possible phase transition into other degrees of freedom. In the present work, we made an effort to probe the nuclear equation of state along with the GW190814 event. Firstly, we examine possible constraints on the nuclear equation of state inferred from the consideration that the low mass companion is a slow or rapidly rotating neutron star. In this case, the role of the upper bounds on the speed of sound is revealed, in connection with the dense nuclear matter properties. Secondly, we systematically study the tidal deformability of a possible high mass candidate existing as an individual star or as a component one in a binary neutron star system. As the tidal deformability and radius are quantities very sensitive on the neutron star equation of state, they are excellent counters on dense matter properties. We conjecture that similar isolated neutron stars or systems may exist in the Universe and their possible future observation will shed light on the maximum neutron star mass problem.


Introduction
In Ref. [1] the authors reported the observation of a compact binary coalescence involving a 22.2 − 24.3 M black hole and a compact object with a mass of a 2.50 − 2.67 M (all measurements quoted at the 90% credible level). The announcement of the GW190814 event [1] triggered various theoretical studies concerning the equation of state (EoS) of dense nuclear matter, in order to explain the possibility of the second partner to be a very massive neutron star. It is worth pointing out that the authors in Ref. [1] did not exclude the possibility that the second partner to be a neutron star or an exotic compact object, i.e. quark star, boson star or gravastar. The consideration of a neutron star as the second partner has been studied in recent Refs. [2][3][4][5][6][7][8][9][10][11][12][13][14]. On the other hand, the case that the second partner is a quark or hybrid star has been explored in Refs. [15][16][17][18][19]. Finally, some modified theories of gravity have also been applied as possible solutions to the problem [20][21][22].
In the present work we concentrate our study on the case where the EoS of neutron star matter is pure hadronic and the hydrostatic equilibrium is described by the general relativistic equations (Tolman-Oppenheimer-Volkoff (TOV) equations) [23][24][25][26][27]. Moreover, we focus on two possible cases, that the second partner to be a slow rotating neutron star and, in the extreme case, to be a very rapidly rotating one (even close to the Kepler limit). Firstly, we consider that the dense nuclear matter properties are described by the MDI-APR [28] (MDI: momentum dependent interaction, APR: Akmal, Pandharipande and Ravenhall) nuclear model. This model has recently been applied successfully in similar studies and, according to our opinion, is a robust guide for neutron star studies. However, since the behavior of dense nuclear matter at high densities remains uncertain, the parametrization of the EoS via the speed of sound is almost inevitable, at least in the framework of hadronic EoSs. To be more specific, we construct large number of EoSs where for low densities (concerning the neutron star crust) we employ well established results, but for the EoS concerning the core, we apply a twofold consideration. For the outer part we employ the EoS predicted by the MDI-APR model and for the inner one, a parametrization based on the speed of sound upper limits is applied. In this study, the transition density and the upper limit of the speed of sound are the two free parameters.
In the first part of our study we concentrate on the effect of the speed of sound and transition density on the bulk neutron star (non-rotating and rapidly rotating) properties including the maximum mass, the Kepler frequency, the kerr parameter, and the maximum central density. We explore under which circumstances the prediction of the mass range 2.5 − 2.67 M of the second partner is possible. More importantly, we provide the constraints which are inferred by the above consideration.
In the second part of the paper we systematically study the tidal deformability of neutron stars by employing the large set of EoSs. We mainly focus in the case of high mass candidates existing as an individual star or as a partner in binary neutron stars system. Until now, there are no observations of an individual, or as a partner of binary system, very massive neutron star (close to 2.5 M ). However, we consider that it is worth to examine this possibility, by focusing on the predictions of the tidal deformability and the radius, quantities that are very sensitive on the neutron star EoS. These quantities are excellent counters on dense matter properties. In the present work, we provide predictions about both the individual and averaged tidal deformability of a hypothetical binary neutron star system where the most massive partner has a mass in the region 2.5 − 2.67 M .
The article is organized as follows: In Section 2, we present the MDI-APR nuclear model along with quantities at the rotating configuration while in Section 3, we present the speed of of sound parametrization of the EoS. In Section 4 we provide the basic formalism for the tidal deformability. The results and the discussion are provided in Section 5 while Section 6 includes the concluding remarks of the present study. Finally, Section 7 contains information about the rotating configuration code.

The MDI-APR model and the rapidly rotating neutron star
The structure of the EoS and the properties of neutron stars are studied under the MDI model. In this model, the energy per particle is available through the form [29,30] where u = n/n s , with n s denoting the saturation density (n s = 0.16 fm −3 ), I = (n n − n p )/n, X 0 = x 0 + 1/2, and X 3 = x 3 + 1/2. The parameters for the symmetric (SNM) and asymmetric nuclear matter are denoted as A, B, σ, C 1 , C 2 , and B , and x 0 , x 3 , Z 1 , and Z 2 , respectively. The finite range parameters are Λ 1 = 1.5k 0 F and Λ 2 = 3k 0 F with k 0 F being the Fermi momentum at the saturation density. The model is flexible enough to reproduce various values of the slope parameter L and symmetry energy S 2 , defined as [28] L = 3n s dS 2 (n) dn n s and S 2 (n) = 1 2 and as a consequence the parametrization of its stiffness. The MDI model, that combines both density and momentum dependent interaction among the nucleons, is suitable for studying neutron star matter at zero (present study) as well as at finite temperature. In particular, although it was introduced by Gale et al. [31,32] to examine the influence of momentum dependent interactions on the momentum flow of heavy-ion collisions, the model has been modified, elaborated, and applied also in the study of the properties of nuclear matter at neutron stars. The advantages of the MDI model are: (a) reproduces with high accuracy the properties of SNM at the saturation density, including isovector quantities, (b) reproduces the microscopic properties of the Chiral model for pure neutron matter and the results of state-of-the-art calculations of Akmal et al. [33], (c) predicts maximum neutron star mass higher than the observed ones [34][35][36], and (d) maintains the causal behavior of the EoS even at densities higher than the ones at the maximum mass configuration.
In this work we apply the EoS produced in Ref. [28], where for the construction of the EoS, the MDI model and data from Akmal et al. [33] had been used (for more details see Ref. [28]). This EoS, not only has the mentioned advantages, but also reproduces the mass of the second component of GW190814 event.
In addition, as a possible scenario is the rotation, we apply rotating configuration in the EoS. In fact, we are interested about the Kepler frequency and the maximum mass of the neutron star at this configuration. This frequency is considered as the one where the star would shed matter at its equator and consequently is the maximum one (mass-shedding limit). An interesting quantity, which connects the gravitational mass with the angular momentum of the star, is the kerr parameter defined as where M and J are the gravitational mass and angular momentum, respectively. For the construction of the rotating equilibrium model we used the RNS code [37]. Furthermore, in Ref. [38] had been found an analytical relation which connects the kerr parameter with the gravitational mass of the non-rotating case (TOV), expressed as which is used by the authors to imply constraints on the possible kerr parameter of the second component, as well as the upper limit of a neutron star mass [13].

Speed of sound formalism and stiffness of equation of state
Another scenario that is followed to reproduce the mass of the second component is the possible stiffness of the EoS. This is achievable by studying the upper and lower limit on the speed of sound, as well as the possible transition density. In this consideration, we have parametrized the EoS, according to Refs. [39][40][41][42][43][44][45][46], as where P and E denote the pressure and energy density, respectively, and E tr is the transition energy density. For the construction of the EoSs, we adopted the following: (a) in region E ≤ E c−edge , we used the equation of Feynman et al. [47] and also of Baym et al. [48] for the crust and low densities of neutron star, (b) in the intermediate region, E c−edge ≤ E ≤ E tr , we employed a specific EoS based on the MDI model and data from Akmal et al. [33], and (c) for E tr ≥ E region, the EoS is maximally stiff with the speed of sound, defined as v s = c (∂P/∂E ) S (where S is the entropy) fixed in the present work in the range [c/ √ 3, c]. The implementation of speed of sound values between the limited ones will lead to results well constrained by the two mentioned limits. Although the energy densities below the E c−edge have negligible effects on the maximum mass configuration, we used them in calculations for the accurate estimation of the tidal deformability. The cases which took effect in this study can be divided into two categories based on the fiducial baryon transition density, n tr = pn s , and the speed of sound as: (a) the ones where p takes the values [1.5, 2, 3, 4, 5], while the speed of sound is parametrized in the two limiting cases, (v s /c) 2 = 1/3 and (v s /c) 2 = 1 and (b) the ones where p takes the values [1.5, 2], while the speed of sound is parametrized in the range (v s /c) 2 = [1/3, 1]. The predicted EoSs are functional of n tr and (v s /c) and implemented to study the possibly existence of a neutron star with ∼ 2.6 M , either non-rotating or a rotating one.
In approach followed in Eq. (5), while the continuity on the EoS is well ensured, the continuity in the speed of sound at the transition density, due to its artificial character, is not. Therefore, in order to ensure the continuity and a smooth phase transition, we employ a method presented in Ref. [49]. We proceeded with the matching of the EoSs on the transition density by considering that, above this value, the speed of sound is parametrized as follows (for more details see Ref. [49]) where the parameters c 1 , c 2 , and w are fit to the speed of sound and its derivative at n tr , and also to the demands v s (n tr ) = [c/ √ 3, c] [39] according to the value of α. Using Eq. (6), the EoS for n ≥ n tr can be constructed with the help of the following recipe [49] The treatment with both discontinuity and continuity in the speed of sound is presented in Table V of Ref. [39]. The outline was that the two approaches converge and consequently the effects of the discontinuity are negligible.

Tidal deformability
It has mentioned that the gravitational waves emitted from the final stages of an inspiraling binary neutron star system are one of the most important sources for the terrestrial gravitational waves detectors [50][51][52][53][54][55][56][57][58]. In such case, properties like the mass of the component stars can be measured. As Flanagan and Hinderer [52] articulated, the tidal effects can be measurable during this final stage of the inspiral.
The response of a neutron star to the presence of the tidal field, is described by a dimensionless tidal parameter, the tidal Love number k 2 . This parameter depends on the neutron star structure; hence its mass and EoS. The tidal Love number k 2 is the coefficient of proportionality between the induced quadrupole moment Q ij and the applied tidal field E ij [52,59], given below where R is the neutron star radius and λ = 2R 5 k 2 /3G is the other tidal parameter that we use in our study, the so-called tidal deformability. The tidal Love number k 2 is given by [52,53] where β = GM/Rc 2 is the compactness parameter of a neutron star. The quantity y R is determined by solving the following differential equation with the initial condition y(0) = 2 [55]. F(r) and Q(r) are functionals of E (r), P(r) and M(r) defined as [50,55] and The numerical solution requires that the Eq. (12) must be integrated self consistently with the TOV equations using the boundary conditions y(0) = 2, P(0) = P c and M(0) = 0 [50,53]. From the solution of TOV equations the mass M and radius R of the neutron star can be extracted, while the corresponding solution of the differential Eq. (12) provides the value of y R = y(R). This parameter along with the quantity β are the basic ingredients of the tidal Love number k 2 .
One parameter that is well constrained by the gravitational waves detectors is the chirp mass M c , which is a combination of the component masses of a binary neutron star system [60,61] 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 is within the range 0 ≤ q ≤ 1.
In addition, another binary parameter that can be measured from the analysis of the gravitational wave signal is the effective tidal deformability [60,61] Λ = 16 13 (12q where the key quantity q characterizes the mass asymmetry, and Λ i is the dimensionless deformability defined as [60,61] By combining Eq. (17) with Eq. (11), one can find that Λ i depends both on star's compactness and the value of y(R).
We notice that Λ i depends directly on the stiffness of the EoS through the compactness β and indirectly through the speed of sound which appears in Eq. (14). The dependence of Λ i on the behavior of y(R) can lead to useful estimations or constraints on the tidal deformability itself. To be more specific, the applied EoS affects also the behavior of Λ regarding the neutron star's mass M and radius R . In our study we use the secondary very massive component of GW190814 system (see Ref. [1]) to examine the tidal deformability and the behavior of dense nuclear matter in the extreme scenario of such a massive neutron star.

Results and Discussion
The merger of a very massive black hole (∼ 23 M ) with a ∼ 2.6 M compact object has recently announced by the LIGO/Virgo collaboration, as the GW190814 event. The scenarios that follow the second merger component are that of (a) the lightest black hole, (b) the most compact neutron star, (c) a rapidly rotating neutron star, and (d) an exotic compact object. In the present work we study only the second and third case scenarios, that is either a compact non-rotating neutron star or a rapidly rotating one.

Slow/rapid rotation: Implications to neutron star properties
In Fig. 1 we display the gravitational mass as a function of the kerr parameter for the MDI-APR EoS. In addition, we note the universal relation (4) [13]. We note that the MDI-APR EoS is in the range of the described limits for the gravitational mass and kerr parameter, being a suitable hadronic EoS to simulate the compact object of ∼ 2.6 M .
In addition, taking into consideration the limiting case that the compact object is rotating at its mass-shedding limit, then constraints on the maximum value of the kerr parameter, the corresponding equatorial radius, and the central energy density are possible. In particular, firstly we employ the relation found in Ref. [28] [13] is shown. In addition, the region for the kerr parameter K max = [0.67, 0.69] from Ref. [28], if the low mass component was rotating at its mass-shedding limit, is presented with the vertical shaded region (right). The markers point the maximum mass configuration at the mass-shedding limit.
for the observable gravitational mass. For the mass of the second component, the kerr parameter lies in the range K max = [0.67, 0.69], which is also noted in Fig. 1. Secondly, using the derived relation from the recent Ref. [62], which connects the maximum value of kerr parameter with the one of compactness parameter, as it is possible to extract a specific range for the corresponding equatorial radius. In this case, the corresponding equatorial radius lies in the range R max = [14.77, 14.87] km. We concentrate now on the microscopic properties of the neutron star, the speed of sound and the transition density. In Fig. 2 we display the gravitational mass and kerr parameter as a function of the transition density for the two limiting cases of the speed of sound based on Ref. [39]. From Fig. 2(a), the intersection of GW190814 mass area with the extracted curves provide us two regions of the possible transition density with respect to the applied speed of sound. By employing the formula from Ref. [39] M max M = α 1 coth α 2 n tr n s we were able to restrict the transition density with respect to the speed of sound. More precisely, we took under consideration two possible cases: (a) non-rotating and (b) maximally-rotating neutron star. In the first case, as Fig. 2(a) shows, the possible transition density region is restricted between the two limiting cases of the speed of sound. However, as the lower limit in the speed of sound is not able to represent the gravitational mass of the low mass component, we found the lower possible speed of sound value that reproduces this mass at the specific value of the transition density. Consequently, the transition density is constrained in the range n tr = [1.5, 3.2] n s and the speed of sound in the range (v s /c) 2 = [0. 45,1]. In the second case, the possible transition density can be constrained both from the gravitational mass and the kerr parameter. From Fig. 2(a) it is clear that the transition density can take all the values in the area under consideration. However, from Fig. 2(b) the transition density is constrained from the lower limit of the speed of sound. According to the derived formula from Ref. [39] K max = α 3 coth α 4 n tr n s  M.R. Hadronic EoSs [28] Cook et al. [63] Salgado et al. [64] (vs/c) 2 = 1/3 -N.R. Approximate Central Baryon Density n c /n s Figure 3. Gravitational mass as a function of the central energy/baryon density at the maximum mass configuration both at non-rotating and maximally-rotating case. Circles correspond to 23 hadronic EoSs [28] at the non-rotating case (N.R.), squares to the corresponding maximally-rotating (M.R.) one, stars to data of Cook et al. [63], and triangles to data of Salgado et al.

Isolated non-rotating neutron star
Firstly, we concentrated our study of tidal deformability on the isolated non-rotating neutron star case, by using two transition densities n tr = [1.5, 2]n s and eight values of speed of sound bounds (v s /c) 2 = [1/3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1]. The values of transition density were taken to be close to the constraints of Ref. [65]. The numerical solution of TOV equations' system, by using the previous bounds for sound speed, provided the mass-radius diagram presented in Fig. 4. In Fig. 4 one can observe two main branches, related to the transition density; the solid (dashed) curves correspond to the n tr = 1.5n s (n tr = 2n s ) case. In each branch, there are bifurcations in the families of EoSs, in analog to each speed of sound boundary condition. The higher the speed of sound, the lighter the representing color of the curves in the figure. The purple solid horizontal line, with the shaded region, indicates the estimation of the recently observed massive compact object of Ref. [1]. As Fig. 4 shows, the branch of EoSs with n tr = 1.5n s provides stiffer EoSs compare to the n tr = 2n s branch. The EoSs of the n tr = 1.5n s case are more likely to provide such a massive non-rotating neutron star, than the n tr = 2n s case in which three EoSs of the total sum lie outside of the shaded region. Especially, between the same kind of transition density n tr the EoSs with higher speed of sound bounds lead to higher values of neutron star mass and radius, hence a high bound of the speed of sound (even more close to the causality as the transition density is getting higher) is needed for the description of such a massive compact object.
From the observation of Fig. 4 a trend across the maximum masses contained in each branch of EoSs, characterized by the speed of sound bound, seems to be inherent. Therefore, we constructed the appropriate diagram of Fig. 5. The cross (star) marks represent the maximum masses of n tr = 1.5n s (n tr = 2n s ) case. As the speed of sound bound is getting higher, the marks' color lightens. Similarly to the previous figure, the purple solid horizontal line, with the shaded region, indicates the estimation of the recently observed massive compact object of Ref. [1]. The red (green) curves represent the following fit formula for the n tr = 1.5n s (n tr = 2n s ) case, given below where w = (v s /c) 2 . The coefficients are given in Table 1. By using the mass estimation of the secondary component of GW190814 system, in combination with the fitting formulas mentioned above, we obtained estimations on the speed of sound values for each transition density n tr scenario. In particular, for a non-rotating massive neutron star with M = 2.59 M the value of the speed of sound must be (a) (v s /c) 2 = 0.485 (n tr = 1.5n s ), and (b) (v s /c) 2 = 0.659 (n tr = 2n s ).
The exact values' interval is given in Table 1. We observe that for higher values of transition density n tr the fitted curve and marks are shifted downwards; hence the higher the point of the transition in density, the smaller the provided maximum mass. The higher values of speed of sound are more suitable to describe such massive neutron stars, until a specific boundary value of transition density n tr in which even the causality would not be suitable. Therefore, a very massive non-rotating neutron star favors higher values of speed of sound than the v s = c/ √ 3 limit. We notice that a lower bound on the transition density n tr is needed to be able in the description of the observed neutron star mergers [65]. Therefore there is a contradiction since the transition density n tr must be above a specific lower limit and not big enough to predict very massive masses. This kind of remark arises in the speed of sound value, respectively.
Moving on to the tidal parameters, we investigated the tidal Love number k 2 and the tidal deformability λ. In Fig. 6 we display the two tidal parameters for the single neutron star case that we examined. In both diagrams, the vertical purple shaded region and line correspond to the GW190814 system's secondary component compact object. There are two main families of EoSs, distinguished by the transition density n tr . In general, the EoSs with higher values of speed of sound bounds lead to larger values on both tidal parameters. Therefore, a neutron star with a higher speed of sound more easily deformable, rather than a more compact star (smaller tidal deformation) with lower speed of sound. As Fig. 6 shows, the EoSs with smaller transition density n tr and higher (v s /c) 2 values are more likely to predict a very massive neutron star of M = 2.59 M . We postulate that a further study with higher transition density n tr would lead to smaller values of tidal parameters, therefore to more compact stars and more difficult to be deformed. In this case, a very high value of speed of sound, even close to the causal limit, would be necessary to predict such a massive non-rotating neutron star.

A very massive neutron star component
Regarding the binary neutron star system case, we considered the scenario of a very heavy component neutron star, in agreement with the recent observation of GW190814 event [1]. Especially, we consider a heavy mass of m 1 = 2.59 M and we let the second star to fluctuate within the range m 2 ∈ (1, 2.59) M . By subtracting the component masses m 1 , m 2 in Eq. (15) we obtain the corresponding range for the values of M c . Since the masses are defined, from the Eqs. (16) and (17), the effective tidal deformabilityΛ can be determined. Fig. 7(a) shows the effective tidal deformabilityΛ as a function of the chirp mass M c , for all the possible binary neutron star systems with such a massive neutron star component. We notice that from the total sum of EoSs that we studied in the single neutron star case above, in the binary case we used only those who can predict a neutron star with 2.59 M mass. In general, there are two families of EoSs, distinguished by the transition density n tr . Inside each family of EoSs, the EoSs with higher speed of sound value provide higher values ofΛ. We notice that for a binary system with m 1 = 2.59 M and m 2 = 1.4 M the chirp mass is M c = 1.642 M . In addition, binary neutron star systems with both heavy components, therefore higher M c , lead to smaller values ofΛ. In such case, possible limits on the lower bound ofΛ may be more suitable to extract constraints on the EoS.
In the same way, Fig. 7(b) shows the dependence ofΛ to the corresponding binary mass ratio q. One can observe a similar behavior of the curves; two main families and the EoSs with higher speed of sound provide higher values ofΛ. As the binary neutron star systems are more symmetric (q → 1), the binary tidal deformability is getting smaller. The highest values ofΛ correspond to the most asymmetric binary neutron star systems. We notice that for a binary system with m 1 = 2.59 M and m 2 = 1.4 M the asymmetry ratio is q = 0.541. Beyond the general behavior ofΛ that we studied above, it is in our interest to examine the radius and possible constraints that can be derived from it. Following the previous steps, we focus on the R 1.4 case of a m 2 = 1.4 M secondary component neutron star. In Fig. 8 we display this dependence; the EoSs are in five main groups, characterized by the transition density n tr . We notice that we expanded our study to transition densities n tr = [1. 25, 1.75, 2.25] to be more accurate in calculations and study in more detail the curves' behavior. The higher speed of sound values correspond to lighter marks' color. In analog to the remarks of the previous Fig. 7, the high speed of sound bounds lead to higherΛ and R 1.4 . Moreover, we applied a fitting expression to the (v s /c) 2 = [0.8, 0.9, 1] cases. The expression was taken to be in the kind of the proposed relations of Refs. [66,67] where the coefficients for each case are given in Table 2.
By applying an upper limit on R 1.4 one can obtain an upper limit onΛ for each case. We adopted the general limit of Ref. [68] that led us to the constraints of Table 2

Concluding remarks
The GW190814 puzzle and its nature through the nuclear EoS has been addressed in this study. In particular, an effort to explain the existence of a ∼ 2.6 M neutron star, which falls into the neutron star -black hole mass gap, had been made both for non-rotating and maximally-rotating neutron stars. For this reason, the MDI-APR EoS and its parametrization for various values of sound speed and transition density in the ranges (v s /c) 2 = [1/3, 1] and n tr = [1.5, 5] n s , respectively, have been studied.
Firstly, we compare the MDI-APR EoS with the applicable range of the proposed kerr parameter from Ref. [13], where the authors restraint it using an analytical relation connecting this property with the maximum mass of a non-rotating neutron star. The results shown that the MDI-APR EoS lies in the range of the gravitational mass of the low mass component, as well as the one for the kerr parameter. In fact, the MDI-APR EoS proposes that the ∼ 2.6 M compact star is a rapidly rotating neutron star, close to its mass-shedding limit.
By considering that the ∼ 2.6 M compact star had been rotating at its mass-shedding limit, possible constraints can be extracted for the corresponding equatorial radius. The kerr parameter at the mass-shedding limit can be expressed as a relation with the gravitational mass, and hereafter a region of K max = [0.67, 0.69] is extracted. This region also includes the upper limit of the relevant region from Ref. [13] in a narrow range. In addition, using a relation that connects the kerr parameter and the compactness parameter at its mass-shedding limit, a possible tight region for the equatorial radius of the star is implied as R max = [14.77, 14.87] km.
The upper limit on the central energy/baryon density is a very interesting property, as it is connected with the evolution of the neutron star and the possible appearance of a phase transition. From this analysis, the central energy density must be lower than the values in the range ε c /c 2 = [2.53, 2.89] 10 15 gr cm −3 , while for the central baryon density, the corresponding range is n c = [7.27, 8.09] n s . The latter can inform us about the stability of the neutron star, as a neutron star with higher values of central energy/baryon density cannot exist, as well as the appearance of the back-bending process.
The transition density, along with the speed of sound, can infer various structures for the EoS. Assuming that the neutron stars is a non-rotating one, the transition density is constrained in the region n tr = [1.5, 3.2] n s while the corresponding value of sound speed must be in the range (v s /c) 2 = [0. 45,1].
Otherwise, if the neutron star is considered as a maximally-rotating one, although the speed of sound implies no constraints, the transition density must be higher than 1.6 n s .
In the case of non-rotating neutron star, the construction of the M-R diagram showed at first glance the cases that can describe the extreme scenario of our study. Moving to a more detailed diagram of the mass vs. the speed of sound bounds, it was feasible to extract stringent constraints on the speed of sound bounds for each case of transition density n tr . For n tr = 1.5n s this bound is (v s /c) 2 ∈ [0.448, 0.52] while for n tr = 2n s is (v s /c) 2 ∈ [0.597, 0.72]. We observe that the first lower bound is in agreement with the bound extracted above, which is a good validation of our result. We postulate that for higher transition densities n tr it is more difficult to achieve such a massive non-rotating neutron star. As the transition density n tr grows the speed of sound would need to be even close to the causal limit.
The study of the tidal parameters for a single non-rotating neutron star allowed us to examine the behavior of EoSs in each case. This lead to the general conclusion that the lower transition densities n tr lead to higher tidal parameters. Therefore the transition density n tr = 2n s corresponds to a more compact and less deformable neutron star. Among the same kind of transition density n tr , the EoSs with higher speed of sound values provide higher tidal parameters. Hence, in a second level across the same kind of n tr EoSs, the higher speed of sound bound signifies that the tidal deformation is higher and the star is less compact.
Concerning the binary neutron star system case, the adoption of a very massive component with m 1 = 2.59 M allowed us to investigate a variety of possible binary neutron star systems with such a heavy component. We notice that as the binary neutron star system consists of both heavy component stars, therefore high chirp mass M c , the effective tidal deformabilityΛ is taking smaller values. Hence, the binary deformation is smaller in such systems. Similarly, the same behavior was noticed in theΛ − q diagram in which the increasing binary mass symmetric ratio q leads to smaller values ofΛ. In the case that the second component has a mass m 2 = 1.4 M , the chirp mass of the system M c and the ratio q are M c = 1.642 M and q = 0.541 respectively.
Lastly, we considered the case of a binary neutron star system with m 1 = 2.59 M with a secondary component m 2 = 1.4 M . This selection permitted the study of the radius R 1.4 and the extraction of possible constraints. In general, the transition density n tr = 1.5n s provides higher values of R 1.4 andΛ than the n tr = 2n s case. The examination of other transition densities n tr permitted us to confirm this behavior. In addition, the high values of speed of sound (v s /c) 2 exhibits a similar behavior; high speed of sound bounds provide higher values on both referred parameters. The adoption of an upper limit on the radius, allowed us to extract some upper limits onΛ for each case of sound speed. Nevertheless the value of n tr , as the speed of sound bound is getting higher the upper limit onΛ performs an analogous course. We conclude that the existence of such a massive non-rotating neutron star would require a significant differentiation from all the so far known cases, consisting in any case a unique and very interesting challenge for physics.

Materials and Methods
The numerical integration of the equilibrium equations for neutron stars is being under the publicly available RNS code [37] by Stergioulas and Friedman [69]. This code is developed based on the method of Komatsu, Eriguchi, and Hachisu (KEH) [70], while modifications where introduced by Cook, Shapiro, and Teukolsky [71]. The input of the code is the EoS in a tabulated form which includes the energy density, the pressure, the enthalpy, and the baryon density.
Funding: This research received no external funding.