In this section, we present our results derived under the assumption that SQM represents the true ground state of matter. However, a key question, that requires a proper discussion, immediately arises: if SQM is more stable than nuclear matter why does normal matter persist? One explanation may be that nuclear matter is metastable [see Horvath et al. (1992), Olesen and Madsen (1994), Iida and Sato (1998), Bombaci et al. (2008), Bombaci et al. (2009), Ren and Zhang (2020) and references therein], separated from the favorable SQM state by a significant energy barrier. Thus, at low densities, the conversion to SQM may be suppressed due to quantum tunneling limitations. However, the extreme neutron star environment, characterized by high densities, may enhance the possibility of a SQM droplet appearing via quantum fluctuations. In addition, extreme events such as supernovae and neutron star mergers may also facilitate SQM formation. Given that the appearance of a SQM seed could trigger the conversion of the entire hadronic star into a strange star, it has long been hypothesized that all compact stars may be of SQM nature.

Proceeding with our calculations, we solved the system of the TOV equations (Oppenheimer and Volkoff, 1939), using a wide range of B ( 57 − 300 M e V ⋅ f m − 3 ) and Δ ( 50 − 300 M e V ) values to identify the M − R curves that are compatible with the light neutron star measurements for HESS J1731-347 and PSR J1231-1411 (at the 1 σ level). The latter was examined utilizing the results of both Salmi et al. (2024a) (Case I) and Qi et al. (2025) (Case II). Of course, we only considered B and Δ pairs that lie within the stability window defined by Equation 4 (fixing the strange quark mass at 95 MeV), while to investigate the possible effects of the a 4 parameter we employed two distinct values, a 4 = 1 and a 4 = 0.7 . The latter value was selected as it results to EOSs similar to the quark matter EOS from pQCD [e.g., see Alford et al. (2005a), Fraga et al. (2001)].

Notably, in the present study we work under the assumption of a universal EOS (a single EOS that accounts for all observations). Therefore, the derived EOSs should be compatible to all state-of-the-art multimessenger constraints on the mass and radius of compact stars. Thus, apart from HESS J1731-347 and PSR J1231-1411, we also utilized the corresponding data related to PSR J0952-0607 (Romani et al., 2022), PSR J0030+0451 (Miller et al., 2019) and PSR J0437-4715 (Choudhury et al., 2024). Lastly, to highlight the peculiar nature of the recent XTE J1814-338 measurement (Kini et al., 2024), we have included it in our analysis, although we did not attempt to interpret it simultaneously with the aforementioned constraints [for some recent works on its reconciliation see Zhou and Huang (2025), Yang et al. (2025), Laskos-Patkos and Moustakidis (2025), Lopes and Issifu (2025), M. Veselský et al. (2025)].

Figure 1 displays the ( B , Δ ) pairs that are compatible with the above measurements for two different values of a 4 , accounting for both estimations of PSR J1231-1411’s properties. The red, purple, blue, yellow and gray parameter spaces correspond to XTE J1814-338, PSR J1231-1411 (Case I), PSR J0952-0607, HESS J1731-347 and PSR J0030+0451 measurements, respectively. The last three areas extend all the way to the horizontal axis, though it is not visible in the graphs. The region between the black solid lines corresponds to PSR J0437-4715, while the one between the purple dashed lines to PSR J1231-1411 (Case II). Figure 1a illustrates the parameter spaces of the above objects for a 4 = 1 considering Case I, while Figure 1b shows the same graph for a 4 = 0.7 . As is evident, although a smaller a 4 value alters the CFL EOS, resulting in a modified stability window, it has only minimal impact on the parameter spaces that are compatible to different observations. Finally, in Figures 1c,d, Case II is examined for a 4 = 1 and a 4 = 0.7 respectively.

Figures 1a,b indicate that HESS J1731-347 and PSR J1231-1411 (Case I), can both be reconciled within the range of parameters that are compatible to the latter, while also respecting the maximum-mass constraints imposed by PSR J0952-0607. This region is characterized by small values for the bag constant, close to the minimum possible one, and large values for the color superconducting gap. As a result, the CFL EOS is extremely stiff, leading to very high maximum masses for all ( B , Δ ) pairs (above 3 M ⊙ ). The latter is demonstrated in Figure 2b, where the maximum masses for all ( B , Δ ) combinations of the purple region of Figure 1b are calculated (for the a 4 = 0.7 case).

Another interesting point from Figures 1a,b is related to the fact that the parameter spaces which are compatible with PSR J0437-4715 and PSR J1231-1411 (Case I) do not overlap at the 1 σ level, implying that they cannot be simultaneously explained within this model. However, this tension is mild, as we have checked that a narrow common parameter space exists, when considering the 2 σ estimation for PSR J0437-4715 (see Figure 2a, where we have plotted the M − R diagrams for representative EOSs which are compatible with PSR J1231-1411).

Considering Case II, Figures 1c,d show that there is a small range of parameters which is compatible with all measurements at the 1 σ level (of course except XTE J1814-338). These areas are defined by the lower boundary from PSR J1231-1411 (purple dashed line) and the upper limit from PSR J0030+0451 (gray region). Notably, for PSR J0030+0451 we employed the measurement of Miller et al. (2019). The aforementioned common areas do not change significantly when adopting the measurements of Riley et al. (2019), and they would increase when using the lower (regarding the radius) of the three estimations provided by Vinciguerra et al. (2024).

The parameter space in which all measurements can be explained by pure CFL matter (Case II) is shown in Figure 2c, alongside the corresponding maximum mass predictions (for the a 4 = 0.7 case). The resulting maximum masses span the range of approximately 2.25 − 3.25 M ⊙ , which is consistent with current observations of heavy pulsars. This maximum mass range, along with that obtained for the parameter space compatible with PSR J1231-1411 (Case I) for a 4 = 0.7 , points toward the possible existence of quark stars with masses exceeding 3 M ⊙ . Notably, such masses reside within the so-called lower mass gap (Shao, 2022; Samsing and Hotokezaka, 2021), in which the nature of compact objects remains undetermined. Interestingly, evidence for compact objects within the mass gap has emerged from multiple recent studies, based on observations of non-interacting binary systems (Thompson et al., 2019; Jayasinghe et al., 2021), radio pulsar surveys (Barr et al., 2024), and gravitational wave detections from binary mergers (Abbott R. et al., 2020; Abbot et al., 2023; Abac et al., 2024; Abbot et al., 2017; Abbot B. P. et al., 2020).

A final remark needs to be made with regards to the possible explanation of the thermal evolution for the CCO in the HESS J1731-347 remnant. Notably, the CCO has a rather high temperature for its estimated age which suggests slow cooling (Di Clemente et al., 2024), analogous to those of purely hadronic stars. Previous works (Di Clemente et al., 2024; Horvath et al., 2023) have qualitatively suggested that superconductivity may suppress the rapid cooling processes that are expected in unpaired SQM and therefore the temperature of the CCO could be explained. Interestingly, according to Horvath et al. (2023) this could occur only under the consideration of vanishingly small gap values Δ < 1 MeV. This is related to the fact that for large values of pairing gaps the suppression of standard cooling processes in the stellar core is dramatic. Therefore, there is a potential tension between the necessity for large gaps, dictated by the large maximum mass constraints, and the low gaps suggested by qualitative analyses on the cooling of the CCO in HESS J1731-347. In addition, it remains an open question whether the low central density expected in the low mass CCO is sufficient for CFL matter to appear for vanishingly small pairing gaps [considering the conditions that need to be met for CFL matter to be energetically favorable (Alford et al., 2001; Alford et al., 2005b)]. In any case, a future precise analysis on the cooling of CFL stars incorporating all of the aforementioned issues, including also potential uncertainties on the estimated temperature and age, is essential to safely conclude about whether the CCO in the HESS J1731-347 could be in fact a CFL strange star.

In the previous section, we considered parametrizations that support the absolute stability of CFL quark matter. At this point, we aim to examine the scenario of explaining all current multimessenger constraints by considering a first-order phase transition, from hadronic to CFL quark matter, in the stellar interior.

Notably, the simultaneous reconciliation of both HESS J1731-347 and PSR J1231-1431 (in Case I) measurements would be rather difficult assuming a purely hadronic EOS. In particular, the hadronic model should be rather soft, at low densities, to support the low radius associated with HESS J1731-347 in the sub-solar mass region, and then it should rapidly stiffen to achieve the higher radius of PSR J1231-1431 at slightly larger masses. Interestingly, the low-density domain of the nuclear EOS can be effectively constrained through parity-violating electron scattering experiments, which aim to measure the neutron skin thickness of different nuclei. In particular, the PREX-II collaboration provided a measurement for the neutron skin thickness of lead, which pointed to stiff EOS behavior at low densities (Adhikari et al., 2021; Reed et al., 2021). Then, the subsequent CREX experiment extracted a puzzling value for the neutron skin thickness of calcium (Adhikari et al., 2022), supporting softer models. More precisely, the simultaneous explanation of both CREX and PREX-II was not possible with the use of traditional energy density functionals (Reinhard et al., 2023; Tagami et al., 2022; Miyatsu et al., 2023; Kumar et al., 2023; Burgio et al., 2024). However, recent attempts have proposed sophisticated modifications to the Lagrangians which are used to describe nuclear matter, and they achieved an explanation of both experimental values (Reed et al., 2024; Kumar et al., 2024). Nonetheless, the predicted EOSs are characterized as extremely stiff. Notably, the latter issue was at some extent resolved in Salinas and Piekarewicz (2024).

Considering that current experimental constraints might point to stiff behavior for the hadronic EOS at low density, in the present study, we are going to employ a nuclear model which is sufficiently stiff, so that it crosses the PSR J1231-1411 contour (Case I), and, therefore, potentially incompatible to the HESS J1731-347 constraints. In particular, we are going to use a widely employed Skyrme model, namely Ska (Köhler, 1976; Typel et al., 2022; Gulminelli and Raduta, 2015; Danielewicz and Lee, 2009; Baym et al., 1971). Thus, in our attempt for explaining all current astronomical measurements, by considering a first-order phase transition to CFL quark matter, we will rely on the hybrid branch to account for the existence of the central compact object in the HESS J1731-347 remnant. It is important to comment that, the employed nuclear model is of nucleonic composition. In principle, as density increases, hyperonic degrees of freedom may appear, altering the properties of the EOS. Given the uncertainties related to hyperon-hyperon and hyperon-nucleon interactions, in the present study we work under a simplified framework that neglects their existence.

In Figure 3a, we depicted the mass-radius dependence for hybrid EOSs constructed by varying the values of Δ and B (fixing a 4 = 1 ). Notably, the pairs of Δ and B are selected so that the phase transition density remains fixed, to allow for the hadronic branch to account for PSR J1231-1411. As one can observe, the hybrid branch crosses the HESS J1731-347 inferred region. However, the EOS softening that is induced due to the phase transition plays a critical role on the determination of the resulting maximum mass. As a consequence, the hybrid branch never reaches the 2 M ⊙ and, therefore, it fails to reproduce one of the most robust constraints derived by compact object observations. In Figure 3b we attempted to resolve this issue by lowering the value of the parameter a 4 . The motivation to do that relies on the knowledge that reducing a 4 reduces the value of pressure for a given baryon chemical potential and, as a consequence, for the same parameter values of B and Δ , shifts the onset of the phase transition to higher densities. Therefore, to recover the desired phase transition onset (i.e., when the hadronic branch crosses the PSR J1231-1411 contour) we need to reduce the value of the bag constant (Baym et al., 2018), leading to EOS stiffening. While the constructed EOSs achieve higher masses (still not the conservative 2 M ⊙ limit), the induced stiffening does not allow for an explanation of the HESS J1731-347 measurement. Finally, we tried to stiffen the quark EOS by adding an ad hoc vector interaction between quarks. We refer to this attempt as ad hoc since the resulting models do not follow from an initial unified description, but are produced by simply adding a vector mean-field term in the CFL MIT bag model expressed in the canonical ensemble [following Kanakis-Pegios et al. (2024), Yang et al. (2021), Yang et al. (2023), the strength of the vector repulsion is quantified via the parameter G v = ( g / m v ) 2 , where g is the coupling constant and m v the mass of the mediating boson]. For more details on the inclusion of this vector term the reader is referred to Kanakis-Pegios et al. (2024). In any case, the results depicted in Figure 3c indicate that such an approach also does not suffice to reconcile all of the considered astronomical constraints. However, while the CFL MIT bag model might fail in reproducing all measurements (when assuming a phase transition), we can not exclude that a different, comprehensive analysis (e.g., via a relativistic mean-field type model), that includes such a stiffening term, may enable the simultaneous explanation of all observations.

Up to this moment, we have only considered Case I for PSR J1231-1411. However, it is clear that the consideration of Case II would not alter any conclusion about our inability of reproducing all observational constraints when we combine the CFL model with a stiff hadronic EOS (similar to the employed one). Nevertheless, future refinement of theoretical models or nuclear experiments that point towards a softer nuclear model may alter the current picture.

Interestingly, we can attempt to reconcile all of the considered astronomical constraints by allowing the hadronic branch to reach the two solar masses and then induce an extremely strong phase transition to make the mass-radius diagram drop to cross the HESS J1731-347 contour. Notably, hybrid stars may be stable even at a descending branch of a mass-radius diagram, assuming that the phase conversion is sufficiently slow (Pereira et al., 2018). The characterization slow refers to the magnitude of the phase conversion timescale τ conv , compared to the period τ osc of radial oscillations that the stellar configurations may undergo. In particular, if τ conv ≫ τ osc , then hybrid stars could be stable even at regions where the well-known turning point criterion would rule them out (Pereira et al., 2018; Rather et al., 2024). The idea of explaining the HESS J1731-347 CCO as a slow stable hybrid star is not new. More precisely, Mariani et al. (2024) established the plausibility of this scenario in their recent work. In the present study, we aim to examine if the CFL MIT bag model can lead to such an explanation. Specifically, the authors of Mariani et al. (2024) considered a constant speed of sound parametrization where there are three free parameters, namely the transition pressure, the speed of sound and the energy density jump. As a consequence, one can construct a sufficiently stiff EOS, by controlling the speed of sound, which would be characterized by a large energy density jump at any chosen transition density. In contrast, in the current framework we have two free parameters (we are keeping the value of a 4 = 0.8 fixed to account for pQCD constraints; see below) and their interplay has to self-consistently provide an appropriate transition density (so that the hadronic branch reaches the 2 M ⊙ ), and a large enough energy density discontinuity to trigger an almost vertical drop of the mass-radius curve. Therefore, it is particularly interesting to examine which pairs of the phenomenological parameters B and Δ could lead to a simultaneous explanation of all current astronomical constraints.

In Figure 4a, we depict the mass-radius dependence for hybrid EOSs constructed by varying the color superconducting gap in the range [0,600] MeV. The bag constant was derived by considering that the phase transition occurs when the hadronic mass reaches the 2 M ⊙ . Notably, all of the hybrid branches, in Figure 4a, are plotted until the last slow stable configuration or terminal configuration as named in Mariani et al. (2024). As one can observe, the reconciliation of the HESS J1731-347 constraints is possible, but it would require extremely large values for both B and Δ (see the legend of Figure 4), at least compared to those found usually in the literature. This is, however, somewhat expected as the majority of studies focuses on phase transitions of moderate energy density jumps and low to moderate transition densities. Hence, the extreme conditions required for the explanation of HESS J1731-347 require extreme values for the phenomenological parameters of CFL quark matter. In any case, it is worth noting that recent and robust constraints have been imposed on the color superconduncting gap in Abbott et al. (2025). While a value of around ∼ 400 MeV is potentially compatible to their results, larger values may be disfavored [see Figure 2 in Abbott et al. (2025)].

Notably, the central baryon chemical potential range that is being considered, in order to achieve the reconciliation of the HESS J1731-347 constraints, reaches such high values that it may cross the regime where the results of pQCD, for the EOS of strongly interacting quark matter, are potentially reliable. Typically, pQCD is considered to be credible at densities around ∼ 40 n 0 (where n 0 is the nuclear saturation density), which correspond to values of baryon chemical potential that most likely exceed 2 GeV (Abbott et al., 2025; Kurkela et al., 2010; Fraga et al., 2014; Kurkela et al., 2014; Annala et al., 2018; Annala et al., 2020; Kurkela et al., 2024). For that matter, we wanted to test whether the constructed EOSs are compatible to pQCD results (for values beyond 2 GeV), at least for values of baryon chemical potential equal to those found at the center of the last slow stable configurations. In Figure 4b, we plotted the quark EOSs, compared to O ( a s 2 ) results of pQCD (Kurkela et al., 2010) [employing the fit to the numerical data provided by Fraga et al. (2014)]. Notably, each curve depicts the pressure (divided by the Stefan-Boltzmann pressure P S B = 3 ( μ B / 3 ) 4 / ( 4 π 2 ℏ 3 ) ) as a function of the chemical potential, neglecting the contribution of color superconductivity [as such a contribution is not encapsulated in the considered pQCD calculations (Fraga et al., 2014)]. Finally, the ’x’ appearing on each curve denotes the terminal configuration for each EOS. As is evident, for values of Δ up to 400 MeV the derived EOSs are compatible to the results of pQCD. However, this is not the case for larger Δ values. The latter is of utmost importance as it places strong constraints on the possible range of parameters that may enable the explanation of the HESS J1731-347 constraints in a slow phase conversion scenario.

A final remark is appropriate regarding the fact that while the existence of slow stable hybrid stars is theoretically intriguing it is not clear how such objects are born. In that sense, future work on possible formation scenarios of such objects would be of utmost importance. Some recent progress on the astrophysical paths that may lead to the existence of twin stars (but for a hybrid branch respecting the turning point criterion) has been made in the work of Naseri et al. (2024). In that direction, we expect that future research will hopefully shed light on the possible existence of slow stable hybrid star branches.