$\Delta$-admixed neutron stars: spinodal instabilities and dUrca processes

Within the covariant density functional theory of nuclear matter we build equations of state of $\Delta$-admixed compact stars. Uncertainties in the interaction of $\Delta(1232)$ resonance states with nuclear matter, due to lack of experimental data, are accounted for by varying the coupling constants to scalar and vector mesonic fields. We find that, over a wide range of the parameter space allowed by nuclear physics experiments and astrophysical observations, cold catalyzed star matter exhibits a first order phase transition which persists also at finite temperature and out of $\beta$-equilibrium in the neutrino-transparent matter. Compact stars featuring such a phase transition in the outer core have small radii and, implicitly, tidal deformabilities. The parameter space is identified where simultaneously $\Delta$-admixed compact stars obey the astrophysical constraint on maximum mass and allow for dUrca processes, which is otherwise forbidden.


INTRODUCTION
With densities exceeding several times the nuclear saturation density, n s ≈ 0.16 fm −3 , the core of neutron stars (NS) has been acknowledged since long ago to represent a unique environment for studies of compressed baryonic matter. In the innermost shells of massive stars several non-nucleonic degrees of freedom (d.o.f.) -such as hyperons, kaon and pion condensates, ∆-resonances and quark gluon plasma -have been conjectured [1] to exist in addition or instead of the nucleonic ones.
The high-precision mass measurements, during this decade, of several massive pulsars in binary systems with white dwarfs [2][3][4][5][6] re-opened the issue of dense matter hyperonisation. A large variety of theoretical frameworks has been employed, such as relativistic density functional theory (DFT), quark-meson coupling model, auxiliary field diffusion Monte Carlo approach, cluster variational method and Brueckner-Hartree-Fock theory. Within the most frequently employed DFT class of models reconciliation among two solar mass NS and hyperonic d.o.f. has been possible either employing a sufficiently stiff nucleonic equation of state (EoS) or going beyond the SU(6) symmetry ansatz to fix the vector meson couplings, for a detailed discussion see [7]. The relativistic quark model and DFT have been also used to address nucleation of ∆(1232)-resonances in stellar matter, see [8][9][10][11][12][13][14][15][16]. With masses lying between those of Σ and Ξ hyperons and an attractive potential in nuclear matter, ∆-resonances are expected to be populated based on the same generic energetic arguments with which hyperonisation is advocated. Similarly to the appearance of any other new hadronic species above its production threshold ∆s will soften the EoS and, thus, modify the mass-radius diagram. Refs. [10,11] have shown that nucleation of ∆ in purely nucleonic NS decreases the maximum mass by up to ≈ 0.5M ⊙ and the radii of intermediate mass NS by up to ≈ 2 km. The magnitude of these effects is nevertheless very sensitive to the underlying nucleonic EoS, ∆s effective mass and strength of ∆ − N interaction. Ref. [9,[12][13][14][15][16] have shown that in hypernuclear stars the only significant effect induced by the nucleation of ∆s consists in the reduction of radii. High precision determination of NS radii, expected from the currently operating NICER mission [17] and also from future x-ray observatories like Athena x-ray telescope [18] and eXTP [19], could thus contribute to shed light on the ∆ d.o.f. in NS and, implicitly, on ∆ − N interaction. Complementary information can, in principle, be provided also by NS cooling history, considering that, because of its negative charge, ∆ − can shift the threshold of nucleonic dUrca to lower densities and, possibly, open new dUrca processes. The first effect is particularly interesting as hyperonisation has been shown to not alter the threshold of nucleonic dUrca [20].
A sufficiently attractive ∆N interaction potential is expected not only to modify NS global properties but also the thermodynamic stability of dense matter. The occurrence of thermodynamic instabilities related to the onset of ∆s has been so far discussed in connection with hot and dense hadronic matter produced in high-energy heavy-ion collisions [21] and, more recently, cold hypernuclear stars [22]. [22] shown that such instabilities manifest over a baryonic density domain which corresponds to the inner core, 4n s n B 7n s , for coupling constants spanning wide domains of values in agreement with available experimental constraints. The limited number of EoS discussed in [22] nevertheless suggests that when the two solar mass constraint is imposed ∆-driven instabilities are ruled out.
The first aim of the present work is to investigate whether ∆-driven instabilities occur also in NS built upon EoS which fulfill the 2M ⊙ constraint [3]. To this aim a systematic investigation of the parameter space is performed. As [22] we use the covariant density functional theory but rely on a different nucleonic EoS and, for the sake of clarity, disregard hyperonic d.o.f. Muons are also disregarded. The second motivation of this study is to identify the parameter space of the ∆ − N interaction where compact stars that obey the two solar mass limit allow for dUrca processes. To better address this issue we select a density-dependent nucleonic EoS which does not allow for n → p + e − +ν e and p + e → n + ν e dUrca processes in stable stars.
The paper is organized as follows. In Section we present the EoS models. In Section the formalism for the analysis of thermodynamic instabilities in multicomponent systems is revisited. Properties of ∆-admixed compact star are discussed in Section . Finally, the conclusions are drawn in Section .

EQUATIONS OF STATE
The phenomenological EoS considered in our study are obtained by extending the density-dependent DDME2 [23] nucleonic parameterisation such as to additionally account for ∆(1232) resonance states of the baryon J 3/2 decuplet. Interaction among different baryonic species is realized by the exchange of scalar-isoscalar (σ), vectorisoscalar (ω) and vector-isovector (ρ) mesons. Coupling constants of ∆s to mesonic fields are expressed in terms of coupling constants of nucleons to mesonic fields x m,∆ = g m,∆ /g m;N , where m = σ, ω, ρ; as common in literature, the assumption on the same density dependence for ∆-and nucleon-meson couplings is made.
Our choice of the nucleonic DDME2 is motivated by the following factors: i) the parameters of isospin symmetric nuclear matter around saturation density are in good agreement with present experimental constraints; ii) the values of the symmetry energy at saturation J = 32.3 MeV as well as its slope L = 51.2 MeV and curvature K sym = −87.1 MeV lie within the domains deduced from nuclear data and neutron star observations; for constraints on L, see [24,25]; for constraints on K sym , see [26][27][28]; the relatively low value of L explains the good agreement of the energy per baryon with ab initio calculations of low density pure neutron matter [29,30] (see fig. 12 in [31]); iii) the maximum mass of purely nucleonic NS, M (N ) max = 2.48M ⊙ [31], built upon DDME2 fulfills the 2M ⊙ observational constraint on the lower limit of maximum masses [3]; iv) the radius of the canonical 1.4M ⊙ NS, R (symmetric 90% credible interval) for a low spin prior extracted from the GW170817 event [34], vi) the agreement with above cited astrophysical data is maintained upon inclusion of hyperons [13,31,35,36]; vii) the nucleonic dUrca process is not allowed in stable stars; this provides us with the perfect framework to assess the impact the nucleation of ∆s has on chemical composition.
In vacuum ∆(1232)s are broad resonances which decay into nucleons with emission of a pion. In medium they are considered to be stabilized by the Pauli blocking of the final nucleon states. Other effects of the medium regard the narrowing of the quasi-particle width and shift of quasi-particle energy to larger values [37,38]. Within the DFT and relativistic quark model it was assumed that ∆s preserve their vacuum masses and have vanishing widths [8][9][10][11][12][13][14][15]. In the present work we employ the same hypothesis.
The information about ∆-nucleon interaction is scarce. Data extracted from pion-nucleus scattering and pion photo-production [39], electron scattering on nuclei [40] and electromagnetic excitations of the ∆-baryons [41] have been reviewed by Drago et al. [9] and Kolomeitsev et al. [12] with the aim of constraining the values of the coupling constants at saturation. Their conclusions may be summarized as follows: i) the potential of the ∆s in the nuclear medium is slightly more attractive than the nucleon potential −30 MeV + U this translates in values of x σ∆ slightly larger than 1, ii) 0 x σ∆ −x ω∆ 0.2, and iii) no experimental constraints exist for the value of x ρ∆ .
The still remaining uncertainties are customarily accounted for by allowing for variation of coupling constants within large domains of values. Previous works [9,[12][13][14][15] have considered the ranges 0.8 ≤ x σ∆ ≤ 1.8, 0.6 ≤ x ω∆ ≤ 1.6 and 0.5 ≤ x ρ∆ ≤ 3. We shall hereafter adopt the same strategy and the following domains The core EoS calculated within DFT is smoothly merged with the crust EoS. For the latter we employ the Negele and Vautherin [42] and Haensel-Zdunik-Dobaczewski [43] EoS.

THERMODYNAMIC INSTABILITIES
Spinodal instabilities manifest as convexity anomalies of thermodynamic potentials expressed in terms of extensive variables. Mathematically they are signaled by negative eigenvalues of the curvature matrix ., N ; f stands for the free energy density; n i represents the number density of each conserved species i, whose total number is N . The eigenvectors associated to negative eigenvalues correspond to the directions in the density space along which density fluctuations are spontaneously and exponentially amplified in order to achieve phase separation. Each negative eigenvalue corresponds to a phase transition whose order parameter is given by the eigenvector. The maximum number of phase transitions in a multi-component system is equal to the dimension of the curvature matrix which, in its turn, is equal to the number of conserved species. The curvature matrix analysis has been often employed in mean-field studies of baryonic matter. For the case of liquid-gas phase transition taking place at sub-saturation densities, see [44][45][46][47]; for strangeness-induced instabilities in dense strange hadronic matter, see [48][49][50][51][52].
Spinodal instabilities of multi-component systems whose direction of phase separation is dominated by one of the conserved densities, n j , are revealed by convex intruders in the free energy density Legendre conjugated with respect to the remaining (N − 1) chemical poten- .,N ;i =j µ i n i or, alternatively, back-bending behavior of the chemical potential µ j = ∂f /∂n j ni;i=1,...,N ;i =j as a function of the conjugate species density n j . For a detailed discussion, see [53].
In the most simple case, that we employ here, in which muons and hyperons are disregarded ∆-admixed stellar matter is made of neutrons, protons, the four ∆isobars (∆ − , ∆ 0 , ∆ + , ∆ ++ ) and electrons. The net charge neutrality condition links the baryon charge density, n Q = n p + n ∆ + + 2n ∆ ++ − n ∆ − , to the electron density, n e = n Q . In the neutrino-transparent matter regime -characterized by vanishing lepton chemical potential, µ L = 0, -the electron density equals the lepton density, n e = n L . This means that the thermodynamic potentials may be expressed in terms of the total baryon number density, n B , and charge (or lepton) density, n Q(L) .
In the following spinodal instabilities will be identified via the negative eigenvalues of C B,L = ∂ 2 e(n B , n L )/∂n B ∂n L , where e stands for the total energy density, the relevant thermodynamic potential at zero temperature. Despite the fact that, being interested in neutron stars, we impose β-equilibrium, the stability analyses are performed in two dimensions, which means that we allow density fluctuations to drive matter out of µ L = 0. Phase transitions will be discussed in the hybrid ensemble (n B , µ L ) [53] for µ L = 0, which corresponds to β-equilibrated matter. Extension to finite temperatures and non-vanishing values of µ L will allow us to check the persistence of ∆-driven instabilities in hot matter out of β-equilibrium.

RESULTS
We now turn to investigate chemical composition and thermodynamic instabilities of cold catalyzed ∆-admixed stellar matter, for different values of the coupling constants of ∆ to mesonic fields. The onset of ∆ isobars is discussed in connection with density thresholds of dUrca processes involving nucleons, n → p + e − +ν e , as well as the two most abundant ∆s: ∆ − → n + e − +ν e , ∆ 0 → p + e − +ν e . Then we address the compatibility of ∆-driven phase transitions and dUrca reactions with compact stars EoS that obey the 2M ⊙ constraint [3]. Finally the relevance of thermodynamic instabilities in hot stellar matter beyond β-equilibrium is studied by  Nucleation of ∆-isobars obviously depends on their interaction potential in nuclear matter, whereσ,ω,ρ stand for the mean field expectation values of the mesonic field and t 3 represents the third component of the isospin, with the convention t 3∆ ++ = 3/2; the dependence on baryonic particle number densities has been omitted for the potential, coupling constants and mesonic fields. Eq. (1) suggests that, for a given nucleonic EoS, more attractive potentials are provided by large values of x σ∆ and low values of x ω∆ . It also shows that in neutron rich matter, as it is the case of NS, the potential felt by positively charged isobars is larger in absolute values than that felt by negative isobars sincē ρ < 0. At variance with this, in symmetric nuclear matter (SNM), all isobars experience the same potential. Quantitative information is provided in Fig. 1 for the extreme cases of pure neutron matter (PNM) and SNM. Different combinations of (x σ∆ , x ω∆ ) are considered; for all cases we assume x ρ∆ = 1. It comes out that, depending on the coupling constants, U determine the onset densities of each species. Here µ B and µ Q stand for baryon and charge chemical potentials and can be expressed in terms of chemical potentials of nucleons, leptons and electrons as µ B = µ n and µ Q = µ p − µ n = µ L − µ e . The left panels of Fig. 2 illustrate the role played by each coupling constant when the other two are kept constant. For all considered sets of (x σ∆ , x ω∆ , x ρ∆ ) the first particle that nucleates is ∆ − . Despite the low attractive potential its population is favored by charge conservation. The other three particles nucleate only for certain values of the coupling constants; the second most favored particle is ∆ 0 , while the less favored one is ∆ ++ . We also note that: i) the onset densities decrease (increase) with the increase of x σ∆ (x ω∆ ), reflecting thus the evolution of U (N ) ∆ with the two coupling constants; ii) x ρ∆ has almost no influence on the onset densities; iii) the onset density of ∆ − is significantly lower than those of the other three particles, which do not differ much one from the other.
By partially replacing the electrons, which regulate the relative abundances of neutrons and protons, ∆ − will implicitly modify the neutron and proton abundances. If protons will become sufficiently abundant such as the triangle inequalities of Fermi momenta are satisfied [54], nucleonic dUrca -forbidden in DDME2 -will become energetically allowed. If, additionally, ∆s are abundant enough to satisfy, along with electrons and nucleons, the appropriate triangle inequalities, also dUrca processes involving these particles will open up.
The right panels of Fig. 2 illustrate the threshold densities of nucleonic dUrca as well as those of two dUrca processes involving ∆s as a function coupling constants. The same cases as in the left panels are considered. With the exception of x σ∆ = 1.1 and x ω∆ ≥ 1.1, all considered sets of parameters allow for nucleonic dUrca, otherwise forbidden in NS built upon DDME2. ∆ − → n + e − +ν e is energetically allowed roughly over the same parameter ranges which allow for nucleonic dUrca; its threshold density is only slightly higher than the one of nucleonic dUrca. ∆ 0 → p + e − +ν e opens up at higher densities and operates over a reduced range of parameters. Similarly to the onset densities, density thresholds of dUrca processes show strong (poor) sensitivity to x σ∆ and x ω∆ (x ρ∆ ).

Spinodal instabilities in cold catalyzed matter
We now turn to investigate the occurrence of ∆-driven instabilities in cold catalyzed matter. With the purpose to limit the parameter space and the motivation that x ρ∆ does not impact the onset of ∆s, we fix x ρ∆ = 1.
We generate EoS of ∆-admixed cold β-equilibrated matter in the neutrino-transparent regime for 0.9 ≤ x σ∆ ≤ 1.5 and x σ∆ − 0.2 ≤ x ω∆ ≤ x σ∆ + 0.2. Thermodynamic instabilities are sought after by calculating the eigenvalues of the curvature matrix C B;L = ∂ 2 e(n B , n L )/∂n B ∂n L . Three situations are encountered: a) no instability, b) one domain of instabilities, c) two domains of instabilities. The sets of (x σ∆ , x ω∆ ) corresponding to the latter two situations are illustrated in Fig. 3  with red stars and, respectively, blue solid dots. Instabilities occur for a narrow domain of 0.15 ≤ x σ∆ −x ω∆ ≤ 0.2 and, for the lowest values of x σ∆ , 0.9 ≤ x σ∆ ≤ 1.1, two instability domains are present. For all matter states featuring spinodal instability only one negative eigenvalue was found for C B;L .
Further insight into the thermodynamic behavior is provided in Fig. 4 in terms of µ B (n B ) (top panels); the chemical composition as a function of baryonic particle number density is represented as well (bottom panels). Two sets (x σ∆ , x ω∆ ) are considered, as mentioned in left and right panels. For the lower value of x σ∆ two instability domains are obtained, signaled by backbendings of µ B (n B ). Inspection of the bottom panel indicates that the first instability is due to the onset of ∆ − and the second one to the onset of the other three isobars. For the higher value of x σ∆ only the low density instability related to nucleation of ∆ − is present. The fact thatdespite their earlier onset -∆ 0 , ∆ + and ∆ ++ do not lead to instabilities is due to the less steep increase of their abundances, caused by the higher value of x ω∆ . µ B (n B , µ L = 0) considered in Fig. 4 corresponds to hybrid ensembles which allow multi-component systems be treated as uni-component ones and have the Gibbs construction reduced to the Maxwell construction [53]. The instabilities induced by ∆ − occur over n s n B 2n s , which corresponds to the outer core of NS. As such it is expected to impact the radii of all mass NS and have little or no effect on the maximum mass. The instabilities induced by ∆ 0 , ∆ + , ∆ ++ occur over 2.5n s n B 4n s , which corresponds to more central shells in the core. In addition to radii also the maximum NS mass may be affected. However some of the coupling constant sets might not be compatible with the 2M ⊙ constraint [3]. To check the issue we plot in Fig. 3 the sets of (x σ∆ , x ω∆ ) for which the maximum NS mass equals 2M ⊙ and 2.2M ⊙ . Only the coupling constants sets lying at the r.h.s. of these curves provide maximum NS masses larger than the value to which the contour corresponds. This means that the occurrence of thermodynamic instabilities associated to ∆ 0 , ∆ + , ∆ ++ is ruled out by the astrophysical limit on pulsar masses. This is the case of EoS discussed by [22] which, indeed, provide low mass NS. The explanation of the low maximum NS masses obtained for this domain of the parameter space consists in the fact that, because of vanishing nucleon effective masses, the baryonic density can not exceed 3 − 4n s . In conclusion the only phase transition compatible with 2M ⊙ is the low density one.
Also shown in Fig. 3 is the domain of parameter sets for which the density threshold of the first energetically allowed dUrca process is smaller than 0.8 fm −3 . This arbitrarily chosen value roughly coincide with the central baryonic density of the maximum mass configuration provided by DDME2 when only nucleonic d.o.f. are accounted for.
The phase diagram of (n, p, ∆, e) matter The phase diagram of (n, p, ∆, e) matter corresponding to the coupling constants set considered in the right panels of Fig. 4 is plotted in Fig. 5. In addition to T = 0 also the domains of phase coexistence between pure (n, p, e)matter and ∆-resonance admixed matter for T = 5 and 10 MeV are represented. As in Fig. 4 the open dots mark, for each temperature and µ L = 0, the coexisting phases. It comes out that, contrary to the liquid-gas phase transition taking place in sub-saturated nuclear matter, the order parameter of the ∆-driven phase transition has only a tiny component along n Q . Other features are nevertheless similar to those of the liquid-gas phase transition: the coexistence domain gets quenched as the temperature increases and second order phase transitions occur as well. First order phase transitions in hot and dense matter have been shown to impact star stability with respect to collapse to a black hole during the core collapse [55]. Second order phase transitions have been proven to drastically reduce the neutrino mean free path and, thus, slow down the proto-NS cooling [49]. Though addressed in connection with a strangeness driven phase transition associated to hyperonisation the two effects are expected to hold also in the case of the presently discussed phase transition and, thus, lead to observable effects. 11 11.5 12 12. 5 Fig. 6. The vertical arrow marks the limits obtained for a 1.4M⊙ NS, 70 < Λ1.4 < 580, as derived from the observation of GW170817 by the LVC collaboration [56].

Global properties of NS featuring a ∆-driven phase transition
In the following we discuss some properties of ∆admixed NS with focus on modifications induced by the phase transition. Equilibrium configurations of spherically-symmetric non-rotating NS are obtained by solving the Tolman-Oppenheimer-Volkoff set of differential equations. The tidal deformabilities are calculated following [57,58].
Figs. 6 and 7 illustrate the mass-radius diagram and, respectively, the tidal deformabilities as a function of gravitational mass. Results corresponding to several ∆admixed NS are confronted with those of purely nucleonic stars. As previously discussed, nucleation of ∆s entails a reduction of the maximum mass of up to a few tens of solar mass [10,11] and a diminish of NS radii by up to 1-2 km [9,[12][13][14][15]. Related to the latter effect a significant reduction of tidal deformabilities is also obtained. The amplitude of these modifications increases with the amount of heavy baryons, that is with the increase of x σ∆ and the decrease of x ω∆ . For the two considered sets of coupling constants for which matter manifests instabilities the Maxwell construction causes an extra reduction of the order of several percent of both radii and tidal deformabilities.

CONCLUSIONS
Inspired by previous findings of Lavagno and Pigato [22] regarding spinodal instabilities in cold catalyzed NS matter with admixture of ∆(1232)-resonances we have investigated in a systematic way the parameter space of the ∆ − N interaction and identified the domains where the onset of ∆ generates thermodynamic instabilities. EoS which fulfill the two solar mass astrophysics constraint manifest spinodal instabilities over n s n B 2n s , which corresponds to the outer shells of the core. When mean-field anomalous behavior is cured by Maxwell construction an extra diminish, of the order of several percent, is obtained for both NS radii and tidal deformabilities. ∆ − -driven instabilities persist also out of βequilibrium in the neutrino-transparent matter and at temperatures of a few tens MeV, which suggests thatif present -the thermodynamic instabilities may affect also other astrophysical phenomena where dense matter is populated.
Finally we have shown that, for a large domain of the parameter space, nucleation of ∆s opens-up the nucleonic dUrca process which is otherwise forbidden. The consequences on thermal evolution of isolated NS and Xray transients are beyond the aim of this work and will addressed elsewhere. Nevertheless it is straightforward to anticipate that ∆-admixed massive NS, whose inner core accommodates unpaired nucleons, might meet the conditions required to describe the data of the faintest Xray transients, XRT SAX J1808.4-3658 and 1H 1905+000 [59]. The situation is interesting the more as to date none of the covariant density functional models with coupling constants fitted on experimental and observational data has been able to simultaneously describe the whole set of thermal data from isolated NS and accreting NS in quiescence, even if the pairing gaps in different channels were allowed to vary over wide ranges and hyperonic degrees of freedom have been accounted for [60].