What Is the Nature of the HESS J1731-347 Compact Object?

One open question of modern physics is it to find out if strongly interacting matter at high densities undergoes a phase transition to deconfined quarks and gluons. To tackle this problem, many fields of physics have been working together, including gravitational wave (GW) and multimessenger astrophysics, nuclear physics, and high-energy physics. While the latter probes finite temperature regimes, astrophysical observations of compact stars (CSs) mainly give access to vanishing temperature and high baryon density regimes, which cannot be probed with terrestrial experiments. In fact, recent detections of binary neutron star (NS) and NS–black hole mergers (Abbott et al. 2018, 2020) have opened a way to constrain zero-temperature NS matter properties during the inspiral phase, whereas the next generation of the GW telescopes is planned to reach enough sensitivity to relate it to the finite temperature equation of state (EOS) during the merger and postmerger phases (Raithel & Most 2023).

Considering the existing observational data of the heaviest known NSs (Antoniadis et al. 2013; Fonseca et al. 2021; Romani et al. 2021, 2022), the EOS of strongly interacting matter at densities above twice the normal saturation density (2n₀) is required to be stiff leading to typical NS radius measurements of ∼11–14 km at M ≳ 1.4 M_⊙. Therefore, despite being in agreement with theoretical calculations for the minimum mass of NSs, e.g., M = 0.88–1.28 M_⊙ (Lattimer & Prakash 2001; Strobel & Weigel 2001), the recently announced measurement of HESS J1731-347 with $M={0.77}_{-0.17}^{+0.20}\,{M}_{\odot }$ and $R={10.4}_{-0.78}^{+0.86}$ km (Doroshenko et al. 2022) challenges our understanding of the EOS at densities 1 − 2 n₀. In fact, simulations of supernova explosions predict the gravitational mass of the lightest possible NS to be M = 1.17 M_⊙, corresponding to a baryonic mass of M = 1.25–1.31 M_⊙ (Suwa et al. 2018).

The HESS J1731-347 measurement is also interesting as the first simultaneous measurement of the mass, radius, and surface temperature of a CS and opens the possibility to study its thermal evolution. The estimated $M={0.77}_{-0.17}^{+0.20}\,{M}_{\odot }$ relies on the fact that the object has a uniform-temperature carbon atmosphere and that the star is located at a distance of 2.5 kpc. Further studies are needed to reduce the uncertainties and understand the validity of the obtained results (Alford & Halpern 2023).

It was proposed that HESS J1731-347 could be a candidate for a strange quark star (QS) comprising all existing mass–radius measurements of CSs (Di Clemente et al. 2022; Das & Lopes 2023; Horvath et al. 2023). However, unpaired quark matter undergoes fast cooling, leading to difficulties in reproducing the surface temperature, presented by Doroshenko et al. (2022). In general, the cooling of a CS is defined by the internal composition, which is determined by the processes that operate inside the star (Potekhin et al. 2015). The fastest cooling direct Urca (DU) process corresponds to β- and inverse β-decay of neutron and d-quark in nuclear and quark matter, respectively. The nucleonic DU process has a threshold controlled by the number density of electrons according to the momentum conservation. As soon as the DU process is on, it leads to an intense neutrino emission and a consequent rapid drop in the surface temperature. On the other hand, the quark DU threshold condition is much easier to satisfy. In fact, such a fast cooling of strange stars contradicts the estimated redshifted surface temperature of HESS J1731-347 compact object ${T}_{s}^{\infty }={2.05}_{-0.06}^{+0.09}$ MK (Doroshenko et al. 2022) itself being rather high as for the age of ∼27 kyr (Beznogov & Yakovlev 2015). This difficulty could be overcome within the strange QS or hybrid star (HS) with an early quark deconfinement (Ivanytskyi & Blaschke 2022; Ivanytskyi et al. 2022; Brodie & Haber 2023) scenarios, where strong quark pairing leads to color superconductivity and strongly suppresses cooling of the quark core (Schaab et al. 1997; Weber 2005), while the baryonic matter (BM) provides a moderate cooling. The possible suppression of the nucleonic cooling is related to the superfluidity of neutrons and superconductivity of protons via the Cooper pairs breaking and formation (PBF; Yakovlev et al. 2001).

The questions at what density the deconfinement phase transition occurs and what the signals of the quark matter formation are still open. The elliptic flow in heavy-ion collisions and a combined analysis of multimessenger constraints (Annala et al. 2022; Huth et al. 2022) suggest that strongly interacting matter softens at a high density, which might correspond to a phase transition to quark–gluon plasma.

In this work, we examine the nature of HESS J1731-347 by considering how the simultaneously measured mass, radius, and surface temperature agree with the present theoretical understanding of the properties of strongly interacting matter and color superconductivity at high densities.

We also discuss an alternative origin of HESS J1731-347 as a dark matter (DM)–admixed NS, a scenario that gained a lot of attention recently (Goldman et al. 2013; Tolos et al. 2015; Ellis et al. 2018; Nelson et al. 2019b; Das et al. 2020; Ivanytskyi et al. 2020; Dengler et al. 2022; Di Giovanni et al. 2022; Karkevandi et al. 2022; Leung et al. 2022; Sagun et al. 2022). In fact, DM could be accumulated in the core of an NS leading to a decrease of the total gravitational mass, radius, and tidal deformability, which we will perceive as the effect similar to softening of the EOS (Giangrandi et al. 2023). At the same time, this scenario of asymmetric noninteracting DM agrees with cosmological and astrophysical observations, e.g., the Bullet Cluster (Randall et al. 2008), and provides a description of HESS J1731-347.

The article is organized as follows. In Sections 2.1, 2.2, and 2.3, we present the NS, strange QS, and HS scenarios, respectively. A detailed description of the thermal evolution of CSs is presented in Appendix A. The possibility of HESS J1731-347 to be a DM-admixed NS with a 1σ confidence interval (CI) is studied in Section 2.4 with a detailed explanation in Appendix B. Section 3 summarizes the results.

2.1. Neutron Star

Modeling the internal structure of HESS J1731-347 requires consistency with the properties of the nuclear matter ground state, chiral effective field theory (Tews et al. 2013), existing constraints on the mass–radius relation (Antoniadis et al. 2013; Miller et al. 2019; Riley et al. 2019; Raaijmakers et al. 2020; Fonseca et al. 2021; Miller et al. 2021; Riley et al. 2021; Romani et al. 2021, 2022) and tidal deformability of NSs (Abbott et al. 2018, 2020). The scenario of purely hadronic matter described with an EOS, which respects the above requirements, suggests a minimal assumption about the HESS J1731-347 nature. At the densities expected inside an NS of subsolar mass such a nuclear EOS can be strongly constrained by the microscopic Brueckner–Hartree–Fock calculations based on realistic nuclear potentials fitted to the nuclear scattering data (Yamamoto et al. 2016, 2017). However, the uncertainty of nuclear EOS obtained with these methods becomes important for modeling NSs heavier than 0.5 M_⊙. Therefore, due to the exploratory reasons in this article, we prefer to consider the possibilities of soft and stiff nuclear EOSs instead of relying on the results of the microscopic calculations. For this purpose, we utilize the set B of the induced surface tension (IST; Sagun et al. 2020) and BigApple (Fattoyev et al. 2020) EOSs, respectively. As is seen from Figure 1, while stiff hadronic EOS is completely outside the 2σ CI, the soft one is able to fit HESS J1731-347 constraint within the 1σ CI.

Zoom In Zoom Out Reset image size

The recently reported data on the thermal evolution HESS J1731-347 suggest its slow cooling (Doroshenko et al. 2022). To reproduce these data, we consider ¹ S₀ and ³ P₂ neutron superfluidity and ¹ S₀ proton superconductivity described by the SFB model (Schwenk et al. 2003), the phenomenological gap obtained from the fit of Cas A (Shternin et al. 2011), and the CCDK (Chen et al. 1993) model, respectively (see Appendix A).

As is seen from the left panel of Figure 2, the scenario of soft hadronic matter with paired nucleons and light-elements envelope (the green band with solid curves) is consistent with the observational data on the thermal evolution of HESS J1731-347.

Zoom In Zoom Out Reset image size

2.2. Quark Star

2.3. Hybrid Star

2.4. Dark Matter–Admixed Neutron Star

We model DM as a relativistic Fermi gas of noninteracting particles with a spin of one-half and mass m_DM accumulated inside NSs with a relative fraction f_DM. For more details about the DM EOS and two-fluid approach, see Appendix B. The dashed–dotted curves in Figure 1 show the effect of DM particles with m_DM = 2.8 GeV and f_DM = 4.75% accumulated in the core of NSs. To exclude the uncertainties of the underlying BM EOS from consideration the scan is performed for the soft IST EOS (blue curve; Sagun et al. 2020) and stiffer BigApple EOS (green curve; Fattoyev et al. 2020). The values of the DM particle's mass and relative fraction were chosen to provide an agreement with the 2σ constraints of HESS J1731-347 for both BM EOSs.

As can be seen, the presence of a dense DM core leads to a strong reduction of the total gravitational mass and radius of stars, which resembles a softening of the BM EOS. This degeneracy between the effect of DM and possible change of the strongly interacting matter properties at high density was studied in Ivanytskyi et al. (2020), Sagun et al. (2022), and Giangrandi et al. (2023). In fact, the scan over mass and fraction of DM in Figure 3 shows that the same effect could be seen for heavier DM particles and smaller fractions. The color map represents the total maximum gravitational mass of DM-admixed NSs for the BigApple EOS. The overlap between the two scans yields the allowed region of mass and fraction of DM (above the red curve in Figure 3) that reproduces M and R measurement of HESS J1731-347 (2σ CI), free from the BM EOS uncertainties.

Zoom In Zoom Out Reset image size

Moreover, the DM-admixed NS scenario remains the thermal evolution unaffected, as DM, within the considered candidate, does not take part in cooling. Effectively, the thermal evolution of an NS admixed with asymmetric DM interacting only through gravity with BM is equivalent to cooling of a pure NS of a smaller mass (A. Ávila et al. 2023, in preparation). Thus, this scenario is consistent with the HESS J1731-347 data.

We analyzed four different scenarios of the possible internal composition of a central compact object within the supernova remnant HESS J1731-347 in light of the recent measurement presented by Doroshenko et al. (2022). We discuss that the soft nucleonic EOS is able to simultaneously reproduce the mass and radius measurements within the 1σ CI, and give a good description of the surface temperature. This scenario also includes the possibility of the soft hadronic EOS at 1–2n₀, while at n ≥ 2n₀ the EOS could attain an extra stiffening due to, e.g., density dependent repulsion between the constituents. Such behavior is also consistent with the proton flow constraint (Danielewicz et al. 2002). From the analysis of various models of pairing gaps, we conclude that a combination of n and p singlet pairing with n triplet pairing described within the SFB and CCDK models, and the phenomenological gap obtained from the fit of Cas A, respectively, with the light-elements envelope provide the best fit of the data.

HESS J1731-347 could also be explained as an HS with an early deconfinement phase transition that occurs below twice nuclear saturation density. Thus, HS would contain a big quark–gluon plasma core, which could potentially lead to rapid cooling related to an operating quark DU process. However, a self-consistent derivation of the quark pairing gap and RDF model formulation makes us conclude that quarks exist in the 2SC phase that suppresses rapid cooling and provides an agreement with the surface temperature measurement.

Moreover, the strange QS scenario can also reproduce mass, radius, and surface temperature very well, as, similar to the HS scenario, paired quarks suppress neutrino emission, while, in comparison to the HS case, the photon emission from the surface will be even lower due to the smaller star's radius. However, the three abovementioned scenarios are in tension with the recent supernova simulations that predict the lowest compact star to be 1.17M_⊙.

As an alternative scenario, we considered HESS J1731-347 to be an NS admixed with DM, which results in the effective softening of the EOS and creation of more compact configurations. This scenario leaves the thermal evolution unaffected, as asymmetric noninteracting DM interacts only through gravity with BM. Based on the performed scan over model parameters we found that fermionic DM particles with mass above 1.15 GeV and fraction above 4.2% provide a full agreement with HESS J1731-347 2σ CI measurement for both stiff and soft baryonic EOSs. The analysis was made for two different EOSs that cover the parameter range to exclude the BM uncertainties from consideration. The performed scan over mass and fraction of DM shows that the same effect could be seen by increasing the DM particle's mass and decreasing its fraction.

We argue that in comparison to GW170817, GW190425, NICER, and heaviest CSs measurements that probe the properties of strongly interacting matter at high densities, HESS J1731-347 provides an important piece of information in the range of 1–2 nuclear saturation density.

Future observations of the HESS J1731-347 object are required, as well as the study of the impact of different effects, e.g., existence of the possible hot/cold spots on the surface of the star, atmosphere composition, distance to the object, etc. While the low mass and radius measurement is confirmed, it will put the most stringent constraint on the strongly interacting matter at 1–2n₀ density range and possible DM-rich environment around the star.

The work of E.G., C.P., and V.S. was supported by national funds from FCT - Fundação para a Ciência e a Tecnologia, I.P., within the Projects No. UIDB/04564/2020, UIDP/04564/2020, EXPL/FIS-AST/0735/2021. E.G. also acknowledges the support from Project No. PRT/BD/152267/2021. C.P. is supported by Project No. PTDC/FIS-AST/28920/2017. The work of O.I. was supported by the program Excellence Initiative–Research University of the University of Wrocław of the Ministry of Education and Science. R.N. acknowledges financial support from CAPES, CNPq, and FAPERJ. This work is part of the project INCT-FNA Proc. No. 464898/2014-5 as well as FAPERJ JCNE Proc. No. E-26/201.432/2021.

Born in supernova explosions or through the coalescence of light CSs, NSs cool down through a combination of neutrino emission from their interior and thermal radiation from the surface. The former process is directly determined by the internal composition of a star. Low- and medium-mass stars usually cool down through so-called slow and intermediate cooling processes, i.e., modified Urca, bremsstrahlung, and PBF (Page et al. 2006).

Once the relative abundances of involved baryonic and leptonic species are high enough, the DU process starts to operate leading to rapid cooling. It corresponds to the β- and inverse β-decay that operate after the triangle inequality of Fermi momenta p_F,i + p_F,j ≥ p_F,k is satisfied. Accounting for charge neutrality and the relation between the Fermi momenta and the number density of each particle, we obtain the DU threshold corresponding to the minimal proton fraction of ∼11% of the total baryon density (Lattimer et al. 1991). In quark matter the threshold for the DU reactions, $d\to u+{e}^{-}+{\bar{\nu }}_{e}$ and u + e⁻ → d + ν_e , is much easier to satisfy leading to an emissivity of ∼ T⁶. For comparison, the modified Urca processes give only ∼ T⁸ (Iwamoto 1980, 1982).

At vanishing temperature, an attraction between nucleons or quarks leads to the existence of pairs, where the particle excitations are gapped and the cooling mechanism is drastically suppressed. At temperatures below the critical temperature of nuclear superfluidity, T ≪ T_c , the neutrino emission is suppressed by a Boltzman factor ${e}^{\tfrac{-{\rm{\Delta }}}{T}}$, where Δ is the energy gap. At T_c , the effect of PBF results in neutrino emissivity of ∼ T⁷ (Page et al. 2006).

In quark matter depending on the abundance of strange quarks, it is possible to distinguish the color-flavor-locked (CFL) phase, in which the quarks form Cooper pairs, whose color properties are correlated with their flavor properties in a one-to-one correspondence between three color pairs and three flavor pairs, and the two-color superconducting (2SC) phase, characterized by the absence of the s-quark and the appearance of u-d diquark condensate in a selected direction in color space (Alford 2009). In fact, the quark pairing could be more diverse, e.g., gapless 2SC, crystalline CSC, gapless CFL (gCFL), etc. (Buballa 2005). However, they are out of the scope of this article.

We performed calculations of the thermal evolution of NSs modeled within the IST EOS as this model provides an agreement with HESS J1731-347 mass–radius measurements within 1σ CI. The IST EOS is formulated in terms of nucleons characterized by an effective hard-core radius yielding a short-range repulsion between them. The latter was fixed from the fit of heavy-ion collision data (Sagun et al. 2018), while the IST contribution was implemented by accounting for an interparticle interaction at high density. The model was applied to describe the nuclear liquid-gas phase transition and its critical point (Sagun et al. 2017), proton flow constraint (Ivanytskyi et al. 2018), and further generalized to describe NSs showing a big application range of the unified IST approach (Sagun et al. 2019a, 2019b, 2019c). The considered parameterization gives the values of the symmetry energy E_sym = 30.0 MeV, symmetry energy slope L = 93.2 MeV and nuclear incompressibility factor K₀ = 201.0 MeV at the normal nuclear density (Sagun et al. 2020). For realistic modeling of the outer layers, the IST EOS is supplemented by the Haensel–Zdunik (HZ) EOS for the outer crust and the Negele–Vautherin (NV) EOS for the inner crust (Negele & Vautherin 1973; Haensel & Zdunik 1990).

The pairing of nucleons in simulations of the thermal evolution of NSs depends on the pairing channel and the considered gap model. By adopting the thermal evolution code described in Sales et al. (2020), we found that the best agreement with HESS J1731-347 data is obtained for ¹ S₀ neutron and proton pairing, in the inner crust and core, respectively, described by the SFB (Schwenk et al. 2003) and CCDK (Chen et al. 1993) models, as well as ³ P₂ pairing of neutrons (Shternin et al. 2011). These results are very much in line with predictions of other works (Negreiros et al. 2012; Dexheimer et al. 2013; Negreiros et al. 2018; Lyra et al. 2023). The results are presented in the left panel of Figure 2.

We also analyze the effect of different envelope compositions: a hydrogen-rich envelope that contains the fraction of light elements η = ΔM/M = 10⁻⁷ (depicted by solid curves in Figure 2) and one containing more heavy elements (dashed curves in Figure 2). Here ΔM is the mass of light elements in the upper envelope.

The cooling simulations presented on the right panel of Figure 2 were performed for HSs described within the DD2npY-T EOS (hadron phase) and RDF approach (quark phase). Modeling of the thermal evolution of HSs incorporated 2SC pairing between quarks obtained in a self-consistent calculation within the RDF model. Moreover, for the inner and outer crusts, we adopted the same HZ and NV EOSs (Negele & Vautherin 1973; Haensel & Zdunik 1990) as for the NS case.

As the 2SC pairing yields a very good agreement with the HESS J1731-347 measured surface temperature, we do not see the necessity to include the CFL phase, as it will cause even stronger neutrino suppression providing an equivalently good data fit.

The observational data were taken from Beznogov & Yakovlev (2015). We consider 2σ error bars for the available data, otherwise a factor of 0.5 and 2 in both the temperature and the age, excluding the upper limits. The sources are 0—Cas A NS, 1—PSR J0205+6449 (in 3C58), 2—PSR B0531+21 (Crab), 3—PSR J1119-6127, 4—RX J0822-4300 (in PupA), 5—PSR J1357-6429, 6—PSR B1706-44, 7—PSR B0833-45 (Vela), 9—PSR J0538+2817, 10—PSR B2334+61, 11—PSR B0656+14, 12—PSR B0633+1748 (Geminga), 13—PSR J1741-2054, 14—RX J1856.4-3754, 15—PSR J0357+3205 (Morla), 16—PSR B1055-52, 17—PSR J2043+2740, and 18—RX J0720.4-3125. The object 8—XMMU J1731-347 in Beznogov & Yakovlev (2015) was substituted by the HESS J1731-347 (Doroshenko et al. 2022). Updated data show a slight increase in the surface temperature.

In this work, we consider the DM component as a relativistic Fermi gas composed of noninteracting massive particles possessing a spin of one-half. The corresponding EOS has been extensively studied in the literature, e.g., by Nelson et al. (2019a), Ivanytskyi et al. (2020), and Sagun et al. (2022).

Following the constraint from the Bullet Cluster (Clowe et al. 2006) on the negligible cross section between BM and DM, we assume two components to interact only through gravity. As a result, stress-energy tensors of two fluids (i = DM, BM) are conserved separately leading to two coupled Tolman–Oppenheimer–Volkoff (TOV) equations (Oppenheimer & Volkoff 1939; Tolman 1939)

$$\begin{eqnarray}&&\displaystyle \frac{{{dp}}_{i}}{{dr}}=-\displaystyle \frac{({\epsilon }_{i}+{p}_{i})({M}_{\mathrm{tot}}+4\pi {r}^{3}{p}_{\mathrm{tot}})}{{r}^{2}\left(1-2{M}_{\mathrm{tot}}/r\right)},\end{eqnarray} \tag{ B1 }$$

where M_tot = M_DM + M_BM and p_tot = p_DM + p_BM are the total gravitational mass and pressure, respectively. Since it is a two-fluid system we define the value of the central density for each component. After the integration of the TOV equations, we get the gravitational masses of each of the components. Using these gravitational masses, the DM fraction can be expressed as ${f}_{\mathrm{DM}}=\tfrac{{M}_{\mathrm{DM}}}{{M}_{\mathrm{tot}}}$. By adjusting the central energy densities of each component we are able to obtain different scenarios of admixed stars, and, in particular, stars with different DM fractions. As was shown by Ivanytskyi et al. (2020), the chemical potentials of two components are related to each other as

$$\begin{eqnarray}&&\displaystyle \frac{d\mathrm{ln}{\mu }_{B}}{{dr}}=\displaystyle \frac{d\mathrm{ln}{\mu }_{\mathrm{DM}}}{{dr}}=-\displaystyle \frac{{M}_{\mathrm{tot}}+4\pi {r}^{3}{p}_{\mathrm{tot}}}{{r}^{2}(1-2{M}_{\mathrm{tot}}/r)}.\end{eqnarray} \tag{ B2 }$$