Baryonic models of ultra-low-mass compact stars for the central compact object in HESS J1731-347

The recent attempt on mass and radius inference of the central compact object within the supernova remnant HESS J1731-347 suggests for this object an unusually low mass of $M = 0.77^{+0.20}_{-0.17}\,M_{\odot}$ and a small radius of $R = 10.4^{+0.86}_{-0.78}$\,km. We explore the ways such a result can be accommodated within models of dense matter with heavy baryonic degrees of freedom which are constrained by the multi-messenger observations. We find that to do so using only purely nucleonic models, one needs to assume a rather small value of the slope of symmetry energy $L_{\rm sym}$. Once heavy baryons are included higher values of the slope $L_{\rm sym}$ become acceptable at the cost of a slightly reduced maximum mass of static configuration. These two scenarios are distinguished by the particle composition and will undergo different cooling scenarios. In addition, we show that the universalities of the $I$-Love-$Q$ relations for static configurations can be extended to very low masses without loss in their accuracy.


Introduction
Over the last decade, there has been exciting progress in gathering new astrophysical information on compact star (CS) parameters such as mass, radius, and tidal deformability. This progress boosts the search for the ultimate state of dense matter and, simultaneously, narrows down the admissible theories of dense matter.
The most massive neutron star yet observed, PSR J0740-+6620, with a mass of 2.08 +0.07 −0.07 M (68.3% CI) [1,2], excludes equation of state (EoS) models which cannot support star with masses that reach that limit. Recently the companion of "black widow" pulsar PSR J0952-0607 was reported to have a mass of 2.35 +0. 17 −0.17 M (68.3% CI) [3], which makes it the fastest and heaviest known neutron star to date. The data from the NICER instrument and pulse modelling of X-ray emission from hot spots on the surfaces of CSs have provided simultaneous mass and radius measurements for PSR J0030+0451 [4,5] and PSR J0740+6620 [6,7]. The detection of gravitational wave (GW) event GW170817 [8] involving a merger of two CSs has placed constraints on the tidal deformability of this system,Λ ≤ 720 [9][10][11][12][13]. While the associated electromagnetic counterparts, AT2017gfo and GRB170817A [14,15], have provided an upper limit on the maximum static mass of neutron stars M ∼ 2.3 M [16][17][18][19][20][21]. The X-ray observations of the temperature of the neutron star in the Cassiopeia A supernova remnant and modeling of its thermal evolution give insights into the superfluid properties of the core matter (see Ref. [22] and references therein).
Email addresses: jiajieli@swu.edu.cn (Jia Jie Li), sedrakian@fias.uni-frankfurt.de (Armen Sedrakian) Very recently, an estimate of the mass and radius of the central compact object (CCO) within the supernova remnant HESS J1731-347 was reported [23] using modeling of its Xray spectrum and a distance estimate from Gaia observations. The mass and the radius of this object are M = 0.77 +0. 20 −0.17 M and R = 10.4 +0. 86 −0.78 km (68.3% CI), respectively, when only parallax priors and X-ray data are considered. Taking in addition the full distance priors and EoS constraint priors into account, the improved estimates are M = 0.83 +0. 17 −0.13 M and R = 11.25 +0. 53 −0.37 km (68.3% CI). This estimation makes the CCO in HESS J1731-347 the lightest neutron star known to date [23,24], and potentially a candidate for an exotic object -"strange star" [23,25,26], as stars made of ordinary matter in supernova explosions have a lower limit on the mass which is above 1.17 M [27]. Previously, the confirmed lightest neutron star with M = 1.174 +0.004 −0.004 M is the companion star in the double neutron star system with the primary being the pulsar J0453+1559 which has an estimated mass [28,29]. Note that these two mass estimates are compatible at better than 2σ accuracy.
Ref. [23] based their discussion on EoS models and associated mass-radius (M-R) curves which were based on the χEFT EoS up to the density of 1.5 ρ sat (where ρ sat is the nuclear saturation density) that are extrapolated to higher densities using the constant speed-of-sound EoS. It was shown that a large number of EoS models meets all current observational constraints. An alternative scheme of analysis is based on the covariant density functional (CDF) approach [30,31], which is a powerful tool that allows one to model finite nuclei and nuclear matter across a wide range of densities and temperatures relevant to neutron stars and their coalescences. In this work, we test how well the CDF approach to baryonic matter can accommodate the pu-tative mass and radius values inferred for the CCO in HESS J1731-347 by Ref. [23].

CDF for baryonic matter
The theoretical framework implemented in this work has been presented in great detail in Ref. [31]. We provide below the most relevant points of our approach for the sake of completeness. The baryonic matter is described within a CDF approach with density-dependent baryon-meson couplings [31,32]. The interaction part of the Lagrangian is given by where ψ stands for the Dirac spinor with index B labeling the spin-1/2 baryonic octet, namely, nucleons N ∈ {n, p} and hyperons Y ∈ {Λ, Ξ 0,− , Σ +,0,− }, while ψ ν for the Rarita-Schwinger spinor [33] with index D referring to the spin-3/2 resonance quartet ∆ ∈ {∆ ++,+,0,− }. The baryons interact via exchanges of σ, ω and ρ mesons, which comprise the minimal set necessary for a quantitative description of nuclear phenomena [34]. In addition, the two hidden-strangeness mesons (σ * , φ) describe interactions between hyperons [35][36][37].
In the nucleonic sector, these are three meson-nucleon (mN) coupling constants (g σN , g ωN , g ρN ) at saturation density ρ sat , and four parameters that control their density dependence. These seven parameters allow one to establish the correspondence between this type of a CDF and the purely phenomenological expansion of the energy density of nuclear matter [38,39] in the vicinity of ρ sat , with respect to the number density ρ and isospin asymmetry δ = (ρ n − ρ p )/ρ where ρ n(p) is the neutron (proton) density. The coefficients of this double expansion are referred to commonly as saturation energy E sat , incompressibility K sat , and the skewness Q sat in isoscalar channel, and the symmetry energy J sym , and its slope parameter L sym in isovector channel.
Among the above-mentioned coefficients, of particular interest are the quantities that arise at a higher order of the expansion, specifically, Q sat and L sym . Their values are weakly constrained by the conventional fitting protocol used in constructing the density functionals. However, they play complementary roles as the value of Q sat controls the high-density behavior of the nucleonic EoS, and thus the maximum mass of nucleonic CSs, while the value of L sym determines the intermediatedensity behavior of the nucleonic EoS, and thus the radii of low-mass stars.
In hyperonic sector, the vector meson-hyperon (mY) couplings are given by either the SU(6) spin-flavor symmetric quark model or the SU(3) flavor symmetric model [31,[35][36][37], depending on the stiffness of the utilized nucleonic EoS model. On the other hand the scalar meson-hyperon couplings are determined by fitting to certain preselected properties of hypernuclear systems. To be specific, as it is common, we determine the coupling constants, g σY , using the following hyperon potentials in the symmetric nuclear matter at ρ sat [40,41]: In addition, we adopt the ΛΛ bond energy at ρ sat /5, which reproduces the most accurate experimental value to date [42], to fix the value of the coupling g σ * Λ . The coupling of remaining hyperons Ξ and Σ to the σ * is determined by the relation g σ * Y /g φY = g σ * Λ /g φΛ [37]. Finally, let us turn to the ∆-resonance sector. The information on the meson-∆ (m∆) couplings is scarce, as no consensus has been reached yet on the magnitude of the ∆ potential. In this work, we limit ourselves to the case where and R ∆σ = g σ∆ /g σN is determined by fitting to the following representative values of ∆-potentials in nuclear matter at ρ sat , where U N is the nucleon isoscalar potential. This means the potential of ∆'s in the nuclear medium is more attractive than the nucleon potential, consistent with the studies of the scattering of electrons and pions off nuclei and photoabsorption with phenomenological models [43]. The resultant values for R ∆σ , respectively, are, R ∆σ 1.10, 1.16, 1.23.
Note that within the range (R ∆σ − R ∆ω ) 0.15, ∆-admixed matter undergoes spinodal instability [44]. Our choice avoids such instabilities with the most attractive case U (N) ∆ = 5/3 U N being at the limit of the stability. This completes our discussion of the parameters entering our CDF for baryonic matter and they are now fully determined. Note that we assume that the hyperon and ∆ potentials scale with a density as the nucleonic one, therefore their high-density behavior is inferred from that of the nucleons.
With this input, we compute the EoS of the stellar matter by implementing the additional conditions of weak equilibrium and charge neutrality that prevail in CSs. We further match smoothly our EoS for the core to that of the crust EoS given in Refs. [45,46] at the crust-core transition density ∼ ρ sat /2.
In the present work, as in Refs. [39,47], we map the nucleonic EoS given by the Lagrangian (1) for each set of parameters Q sat and L sym . For our analysis below we adopt the lower-order coefficients in the isoscalar channel, i.e., E sat = −16.14, and K sat = 251.15 MeV, as those inferred from the DDME2 parametrization [32] which were adjusted to the properties of finite nuclei. For the isovector channel, instead fixing J sym = 32.31 MeV the DDME2 value, we hold the value of the symmetry energy at a lower density ρ c = 0.11 fm −3 at the constant value E sym (ρ c ) = 27.09 MeV [39]. The density ρ c is distinguished by the fact that a large variety of nuclear models predict almost identical values of E sym (ρ c ). In the following discussion, we will use the pair of parameters (Q sat , L sym ) to identify any specific nucleonic CDF. To appreciate the quality of current functional and the range of variation of parameters we show in Fig. 1 the energy per particle of symmetric nucleonic matter (SNM) and pure neutron matter (PNM) as a function of density for six representative (Q sat , L sym ) pairs, while keeping all other coefficients at their default values of DDME2 parametrization [32,39]. We further restrict the set of EoS by choosing only those which reproduce the combined result for PNM derived from several recent many-body calculations with χEFT interactions for densities up to ∼ 1.3 ρ sat [48]. This results in the range 30 L sym 80 MeV, which is independent of the values of Q sat . In addition, we have checked that all the alternative parametrizations can reasonably reproduce the binding energies and charge radii of several closed-shell nuclei with ∼ 2% relative deviation.

CSs with only nucleons
We start our discussion with the CS models containing only nucleons. We limit the discussion to three sets of representative EoS models by taking three values of Q sat = −600, −200, and 600 MeV, and a range of L sym which allow an overlap of our models with the range limited by the analysis of the CCO in HESS J1731-347 [23]. Large values of L sym correspond to a stiffer EoS close to nuclear saturation density, as shown in Fig. 1, and lead to larger radii for low-mass stars. Larger values of Q sat imply EoS which are stiffer at high density and, thus, predict larger maximum masses for static nucleonic stars [39,49]. For Q sat = −600 MeV the maximum mass is about 2.0 M , which matches the mass of PSR J0740+6620 [1,2]; for Q sat = −200 MeV the maximum mass is consistent with the (approximate) theoretical upper limit on the maximum mass of static CSs ∼ 2.3 M inferred from analysis of the GW170817 event and the corresponding electromagnetic counterparts in a scenario where a supramassive remnant formed after a merger collapses into a black hole [16][17][18][19][20][21]; finally, for Q sat = 600 MeV the maximum mass is slightly higher than 2.5 M , which would be  [9,50], the ellipses indicate the regions compatible with the inferences from NICER observations [4][5][6][7], the contours show the M-R constraints for the CCO in HESS J1731-347 [23]. See text for details.
compatible with the mass of the secondary in the GW190814 event [50] and its interpretation as a nucleonic CS [51][52][53].
The M-R relations for our nucleonic EoS models are shown in Fig. 2, along with the current astrophysical observational bands derived from X-ray or gravitational wave observations. These include: (i) the ellipses obtained by the NICER collaboration which simultaneously determined the radius and mass of PSR J0030+0451 and J0740+6620 via X-ray pulse-profile modeling (at 68.3% and 95.4% CIs) [4][5][6][7]; (ii) the M-R contours for the CCO in HESS J1731-347, where solid lines correspond to the case that only parallax priors and X-ray data are considered, and the dashed lines correspond to the joint fit including all prior information (at 68.3% and 95.4% CIs) [23]; (iii) the regions for each of the two CSs that merged in the GW170817 event (at 90% CIs) [10]; and (iv) the mass of the secondary component of the GW190814 event (at 90% CI) [50].
As seen from Fig. 2, these observational ellipses can be accounted for by appropriate choices of the parameters Q sat and L sym at the 95.4% CI. The choice of these parameters gives sufficient flexibility to allow for intermediate-density soft EoS to account for ellipses for CCO in HESS J1731-347 and GW170817 and sufficiently stiff EoS at high density to account for the ellipses of PSR J0740+6620. To give more specifical examples, a 95.4% CI consistency is achieved by choosing L sym 60 MeV and Q sat ∼ −600 MeV, L sym 40 MeV and Q sat ∼ −200 MeV, and L sym 30 MeV and Q sat ∼ 600 MeV. Among these, the model with (Q sat , L sym ) = (−600, 30) MeV intersects also the 68.3% CI ellipse. Thus, we conclude that the nucleonic models are consistent with the current information available from multimessenger astrophysics within 95.4% CI, in particular, the small radius of the CCO in HESS J1731-347, if fairly low values of L sym are adopted. These low values for L sym lie outside the 1σ confidence placed by the analysis of neutron skin in the PREX-II experiment for 208 Pb [54,55], but are compatible with the analysis including also dipole polarizabilities in a set of finite nuclei [56], the analysis of the CREX experiment for 48 Ca [57,58], as well as with recent χEFT predictions [59]. A large sample of nucleonic EoS based on CDF theory with nonlinear couplings was considered in Ref. [24] and was found to be inconsistent with the 68.3% CI ellipses of HESS J1731-347. This is in contrast to our finding, even though our EoS are consistent with the χEFT computations of PNM (see Fig. 1) which were used to tune the CDFs in Ref. [24]. The difference in the conclusions arises from the fact that our collection includes EoS with small values of the slope of symmetry energy, e.g., (Q sat , L sym ) (−600, 30) MeV, which are still consistent with the χEFT band. This is in contrast, the collection of Ref. [24] which is limited to larger values of L sym 40 MeV. Anticipating the discussion in the next section, we point out here that if the ∆-resonance threshold is low enough, their nucleation relieves further the tension between the observations and CDF-based models.

CSs with heavy baryons
Next, we consider models which allow for the nucleation of heavy baryons in dense matter. The hyperonic EoS models that support massive two-solar mass CSs require the nucleonic sector to be a priori stiff, which results in a large radius of a canonical-mass star R 1.4 13 km [36,37,[60][61][62][63][64]. Only during recent years, hyperon-∆ admixed CDF-based models were developed which were compatible with the nuclear data, radius measurements of CSs, and the tidal deformability inferred from the analysis of the GW170817 event [31,47,[65][66][67]. The appearance of ∆-resonances, in parallel with hyperons, does not affect the maximum mass of a static CS, but reduces the radius of the star by tens of percent, depending on the value of ∆potential in nuclear matter U ∆ . Therefore, it is convenient for further discussion to denote underlying EoS model by the potential of ∆'s, whereby the parameters specifying the nucleonic and hyperonic sectors are fixed.
In Fig. 3 we present the M-R relations for purely nucleonic, hyperonic, and hyperon-∆ admixed stellar matter for several parameter values, along with current astrophysical observational constraints. In panel (a) the EoS models are constructed from nucleonic model with Q sat = 600 MeV, SU(6) symmetric parameterization for the hyperonic sector which results in a maximum mass M max 2.0 M , and a SU(3) symmetric parameterization which yields M max 2.3 M . To assess the range of variations in the sequences arising from the uncertainties in ∆ potential U ∆ and the value of L sym in nucleonic sector which are crucial for the radius of the star, we consider for each model U ∆ /U N = 1, 4/3, and 5/3 for L sym = 30, 60 MeV. In panel (b) we present results for alternative EoS models which are constructed from nucleonic model has Q sat = −200 MeV in combination with a SU(3) symmetric model for the hyperonic sector, which gives M max 2.0 M . In Fig. 3 (a), we show the results for the nucleonic model with (Q sat , L sym ) = (600, 30) MeV which are marginally compatible with the CCO in HESS J1731-347 constraint at 95.4% CI limit. Allowing for hyperons and ∆-resonances and varying the value of U ∆ we observe the reduction of the stellar radius with increasing value of U ∆ . The M-R track in this case could lie within the 68.3% CI ellipse. For a (isovector) stiffer nucleonic EoS model with (Q sat , L sym ) = (600, 60) MeV which produces M-R tracks outside the HESS J1731-347 region, the onset of hyperons and resonances leads to M-R tracks that are within 95.4% CI region if a sufficiently attractive ∆ potential is chosen, e.g., U ∆ /U N = 5/3. In Fig. 3 (b), for EoS models that are constructed from a (isoscalar) softer nucleonic model with  be noted that if HESS J1731-347 is a baryonic star it contains likely only non-strange particles, due to its low mass. Indeed, as shown in Fig. 3, the onset of hyperons occurs in stars with larger mass with M 1.2 M .  Fig. 4 (a), and the deviation in radius for a fixed mass is found less than 3%. For star with a mass ∼ 0.77 M , the inner core density is about 2 ρ sat , around where significant difference in particle composition is observed. Thus, the knowledge of the mass and radius alone does not allow one to distinguish the internal composition of CSs. In fact, even under the restrictive assumption of purely nucleonic stars, current empirical knowledge allows for multiple solutions for the composition of β-equilibrium matter, i.e., the composition of CSs cannot be determined unequivocally [68].

Composition of the CCO in HESS J1731-347
However, as well known, the particle content of dense matter is an important factor in the cooling of CSs due to new channels of neutrino emissions in the presence of heavy baryons. Because of same values of L sym the nucleonic and heavy baryon admixed models lead to a similar composition at low density ρ 0.25 fm −3 in Figs. 4 (b) and (c). We note that the requirement of a small radius of CCO in HESS J1731-347, and equivalently a low value of coefficient L sym preclude the onset of the direct Urca process in such models. In contrast, in the ∆admixed hyperonic matter a trace fraction of hyperons and ∆'s will open up a host of direct Urca processes involving these particles. According to Fig. 4 (c) the ∆-admixed hyperonic EoS allows the direct Urca processes in a large region of the density,  Fig. 4 (c) as well. The onset of hyperons then will affect the cooling behaviour of stars more massive than HESS J1731-347, see Ref. [31] and references therein.

Tidal deformability and I-Love-Q relations
The gross quantities of a CS such as the mass, radius, deformability, moment of inertia, quadrupole moment, etc., sensitively depend on the microscopic EoS. In Fig. 5 we show the mass vs dimensionless tidal deformability relations for our models that satisfy all relevant astrophysical constraints at the 95.4% CI. All these models satisfy the constraint placed for a 1.  HESS J1731-347 are indicated as well. Since the low-mass stars have very large tidal deformability Λ, spanning the range from 10 3 to 10 5 , it is clear that they could be good targets for GW observatories if involved in a merger process. Various approximately universal relations connecting different CS properties have been established and intensively studied in recent years. Of particular interest is the I-Love-Q relations that first discovered in Ref. [69] which connect the dimensionless moment of inertiaĪ, tidal deformability Λ, and the spin-induced quadrupole momentQ of neutron stars in slow rotation approximation. These relations can be numerically described with a polynomial [69,70]: ln y = a 0 + a 1 ln x + a 2 (ln x) 2 + a 3 (ln x) 3 + a 4 (ln x) 4 , (7) where pairs (x, y) represent (Λ,Ī), (Λ,Q) and (Q,Ī). These relations are commonly studied for stars with mass M 1.0 M [70].
Here we briefly investigate whether these relations hold for stars with heavy baryons, in particular, the hyperon-∆ admixed CSs, especially in the low-mass region extending down to about 0.2 M . Our results for I-Love and Q-Love relations, derived from a collection of our models that satisfy all relevant astrophysical constraints at the 95.4% CI, are shown in Fig. 6, together with the fits according to Eq. (7), where the bottom panels present the fractional differences between the data and the fits. It is seen in Fig. 6 that the absolute fractional differences are 1% for all two relations. Clearly, the universality holds for the third pairĪ-Q as well. The best-fit coefficients with Eq. (7) for the three relations are summarized in Table 1, which are consistent with those found in Ref. [70]. Thus, our results extend the previously reported universal I-Love-Q relations for CSs into the low-mass domain.

Summary
In this work, we have studied the problem of the nature of the CCO in HESS J1731-347 by using a CDF approach for dense matter with and without heavy baryons in high-density regions. While some authors suggested the "strange star" scenario for this object [25,26], here we explore a "less strange" alternative. The density functional we used was parameterized and varied via two parameters -the skewness coefficient Q sat of symmetric nuclear matter and the slope coefficient L sym of the symmetry energy. The hyperon potentials were tuned to the most plausible potentials extracted from hypernuclear data. The ∆-potential was taken to be in the range 1 ≤ U ∆ /U N ≤ 5/3,a reasonable range consistent with the current data.
We found that purely nucleonic models for the EoS can accommodate the estimate for mass and radius of the CCO in HESS J1731-347, but only if the slope of symmetry energy coefficient L sym is fairly small, i.e., L sym 60 MeV for models with Q sat ∼ −600 MeV, and L sym 30 MeV for models with Q sat ∼ 600 MeV. These low values for L sym are compatible with the analysis of the PREX-II and CREX experiments of neutron skins, dipole polarizabilities, as well as recent χEFT computations. Allowing for hyperon and ∆-resonance populations in dense matter, our EoS models allow a broader range of values of L sym . We note that these models do allow for maximum masses in the range 2.0 M M max 2.3 M . Another important point of our analysis is that our EoS models with or without heavy baryons which predict very similar M-R relations have different compositions with different sets of neutrino emission processes and, therefore, cooling behavior. The core matter of nucleonic CS does not allow for the direct Urca process, as the low values of coefficient L sym used in our 6 models result in a low proton fraction. Note that the DDME2 parametrization, i.e., without variations of the L sym parameter does not allow direct Urca as well. In the case of the ∆-admixed hyperonic matter, direct Urca processes involving ∆'s and hyperons are allowed for a large region of the density relevant for high mass CSs. Finally, we showed that the presence of heavy baryons in CSs does not alter the universalities of the I-Love-Q relations for static configurations. We further extended this relation to the low-mass domain not been studied previously, but which is relevant for the CCO in HESS J1731-347. These relations were found accurate up to 1% for mass range from 0.2-2.5 M .