Role of the Symmetry Energy on the Structure of Neutron Stars with Unified Equations of State
Nicolas Chamel, John Michael Pearson, Alexander Y. Potekhin, Anthea F. Fantina, Camille Ducoin, Anup Dutta, Stéphane Goriely, Loïc Perot
aa r X i v : . [ a s t r o - ph . H E ] A p r Role of the Symmetry Energy on the Structure of NeutronStars with Unified Equations of State
Nicolas Chamel , John Michael Pearson , Alexander Y. Potekhin , Anthea F.Fantina , Camille Ducoin , Anup Dutta , St´ephane Goriely and Lo¨ıc Perot a) Corresponding author and speaker: [email protected] Institut d’Astronomie et d’Astrophysique, CP-226, Universit´e Libre de Bruxelles, 1050 Brussels, Belgium D´ept. de Physique, Universit´e de Montr´eal, Montr´eal (Qu´ebec), H3C 3J7 Canada Io ff e Institute, Politekhnicheskaya 26, 194021 St. Petersburg, Russia Grand Acc´el´erateur National d’Ions Lourds (GANIL), CEA / DRF - CNRS / IN2P3, Boulevard Henri Becquerel, 14076Caen, France Institut de Physique Nucl´eaire de Lyon, CNRS / IN2P3, Universit´e Claude Bernard Lyon 1, Villeurbanne, France School of Physics, Devi Ahilya University, Indore, 452001 India
Abstract.
The role of the symmetry energy on the internal constitution and the global structure of a cold nonaccreted neutronstar is studied using a set of unified equations of state. Based on the nuclear energy-density functional theory, these equationsof state provide a thermodynamically consistent treatment of all regions of the star and were calculated using the four di ff erentBrussels-Montreal functionals BSk22, BSk24, BSk25 and BSk26. Our predictions are compared to various constraints inferredfrom astrophysical observations including the recent detection of the gravitational wave signal GW170817 from a binary neutron-star merger. INTRODUCTION
Formed in the aftermath of gravitational core-collapse supernova explosions, neutron stars are initially fully fluidand very hot with internal temperatures as high as 10 − K [1]. In such a furnace, the compressed stellarmaterial is likely to undergo all kinds of nuclear and electroweak processes eventually reaching full thermodynamicalequilibrium. The proto-neutron stars cool rapidly by emitting neutrinos. After a few months, the surface of the star -possibly surrounded by a very thin atmospheric plasma layer of light elements - still remains hot enough to be liquid.However, the layers beneath crystallize thus forming a solid crust. After about a hundred thousand years, heat from thestellar interior di ff uses to the surface and is dissipated in the form of radiation. Assuming that full thermodynamicalequilibrium is maintained throughout the di ff erent cooling stages, the interior of an old neutron star consists of “coldcatalyzed matter”, i.e., electrically charge-neutral matter in its absolute ground state with density exceeding that foundinside the heaviest atomic nuclei [2, 3]. This hypothesis can reasonably be expected to be valid in any neutron starthat has not accreted material from a stellar companion, but will otherwise fail because of the relative slowness withwhich accreted matter reaches nuclear equilibrium. Although properties of hot nuclear matter o ff equilibrium can beexplored in terrestrial laboratories up to a few times saturation density n ≃ .
16 nucleons fm − by analysing heavy-ion collisions, the extreme conditions prevailing in the very dense and highly neutron-rich interior of a cold neutronstar cannot be experimentally reproduced.Neutron stars thus o ff er the unique opportunity to probe novel phases of matter [4]. In particular, the interior ofa neutron star can be decomposed into three distinct regions below its thin atmosphere and liquid “ocean”: an outercrust made of a body-centered cubic lattice of exotic nuclei in a charge neutralizing background of highly degenerateelectrons, an inner crust consisting of neutron-proton clusters immersed in a neutron sea (possibly enriched withprotons at su ffi ciently high densities), and a liquid core of nucleons and leptons. Other particles such as hyperonsmight be present in the most massive neutron stars [1, 4], but we do not consider this possibility here. Because theore of a neutron star contains asymmetric nuclear matter at densities up to ∼ n , astrophysical observations mayshed light on the longstanding issue of the density dependence of the symmetry energy, defined here as the di ff erencebetween the energy per nucleon of pure neutron matter (NeuM) and that of symmetric nuclear matter (SNM). Tothis end, we have developed a series of four unified equations of state (EoS) in the framework of the nuclear energy-density functional (EDF) theory [5]. These EoS provide a thermodynamically consistent description of all regions ofa neutron star. The underlying EDFs were precision fitted to essentially all experimental atomic mass data and weresimultaneously adjusted to realistic NeuM EoS as obtained from many-body calculations. They mainly di ff er in theirpredictions for the symmetry energy. These EoS are used to study the role of the symmetry energy on the structure ofa cold nonaccreted neutron star. NUCLEAR ENERGY-DENSITY FUNCTIONAL THEORYFormalism
The nuclear energy-density functional (EDF) theory aims at providing a universal description of various nuclearsystems, from atomic nuclei to extreme astrophysical environments such as neutron stars and supernova cores (see,e.g. Ref. [6] for a review).In this theory, the energy E of a nuclear system (assumed here to be invariant under time reversal) is expressedas a universal functional of the so-called normal and abnormal density matrices defined by [7] n q ( r , σ ; r ′ , σ ′ ) = < Ψ | c q ( r ′ , σ ′ ) † c q ( r , σ ) | Ψ > , (1) e n q ( r , σ ; r ′ , σ ′ ) = − σ ′ < Ψ | c q ( r ′ , − σ ′ ) c q ( r , σ ) | Ψ > , (2)respectively where | Ψ > is the many-body wave function, c q ( r , σ ) † and c q ( r , σ ) are the creation and destruction opera-tors for nucleons of charge type q ( q = n , p for neutron, proton respectively) at position r with spin coordinate σ = ± E under the constraint of fixed numbers of neutrons and pro-tons. The minimization is usually carried out by expressing the density matrices in terms of auxiliary independenttwo-component quasi-particle wavefunctions ψ ( q )1 k ( r , σ ) and ψ ( q )2 k ( r , σ ), as n q ( r , σ ; r ′ , σ ′ ) = X k ( q ) ψ ( q )2 k ( r , σ ) ψ ( q )2 k ( r ′ , σ ′ ) ∗ (3)and e n q ( r , σ ; r ′ , σ ′ ) = − X k ( q ) ψ ( q )2 k ( r , σ ) ψ ( q )1 k ( r ′ , σ ′ ) ∗ = − X k ( q ) ψ ( q )1 k ( r , σ ) ψ ( q )2 k ( r ′ , σ ′ ) ∗ , (4)where the index k represents the set of suitable quantum numbers, the symbol ∗ denotes complex conjugation, and thequasiparticle states are subject to the following conditions [7]: X σ Z d r n q ( r , σ ; r , σ ) e n q ( r , σ ; r , σ ) − e n q ( r , σ ; r , σ ) n q ( r , σ ; r , σ ) = , (5) X σ Z d r n q ( r , σ ; r , σ ) n q ( r , σ ; r , σ ) + e n q ( r , σ ; r , σ ) e n q ( r , σ ; r , σ ) = n q ( r , σ ; r , σ ) . (6)Although the EDF theory can potentially describe the exact ground state of the system, the correspondingEDF is unknown. Phenomenological functionals have been traditionally constructed from density-dependent e ff ec-tive nucleon-nucleon interactions in the framework of the self-consistent mean-field methods [8]. However, it shouldbe stressed that the resulting EDFs actually include beyond mean-field e ff ects (with respect to “bare” nucleon-nucleoninteractions) through the fitting to nuclear data. Although such a formulation imposes stringent restrictions on the formof the EDF, it guarantees the cancellation of self-interaction errors [9] (nonetheless, the EDFs may still be contami-nated by many-body self-interactions errors, see, e.g. Ref. [6] and references therein). The Brussels-Montreal EDFsonsidered here are based on generalized Skyrme e ff ective interactions of the form [10] v ( r i , r j ) = t (1 + x P σ ) δ ( r i j ) + t (1 + x P σ ) 1 ~ h p i j δ ( r i j ) + δ ( r i j ) p i j i + t (1 + x P σ ) 1 ~ p i j · δ ( r i j ) p i j + t (1 + x P σ ) n ( r ) α δ ( r i j ) + t (1 + x P σ ) 1 ~ h p i j n ( r ) β δ ( r i j ) + δ ( r i j ) n ( r ) β p i j i + t (1 + x P σ ) 1 ~ p i j · n ( r ) γ δ ( r i j ) p i j + i ~ W ( ˆ σ i + ˆ σ j ) · p i j × δ ( r i j ) p i j , (7)where r i j = r i − r j , r = ( r i + r j ) / p i j = − i ~ ( ∇ i − ∇ j ) / ˆ σ i and ˆ σ j are Pauli spin matrices, P σ is the two-body spin-exchange operator, and n ( r ) denotes the average nucleon number density. The t and t termswere originally introduced to remove spurious spin and spin-isospin instabilities [11]. Nuclear pairing is treated usinga di ff erent e ff ective interaction of the form v ( r i , r j ) =
12 (1 − P σ ) f ± q v π q [ n n ( r ) , n p ( r )] δ ( r i j ) , (8)where n n ( r ) and n p ( r ) denote the average neutron and proton number densities respectively. The pairing strength v π q [ n n ( r ) , n p ( r )] given by [12] v π q [ n n , n p ] = − π ~ k F q M q " ε ( q )F ∆ q + Λ ε Λ ε ( q )F − , (9) Λ ( x ) = log(16 x ) + √ + x − (cid:16) + √ + x (cid:17) − , (10)where M q is the nucleon mass, k F q = (3 π n q ) / is the Fermi wave number, ε ( q )F = ~ k q / (2 M q ), and ε Λ is a cuto ff inthe single-particle energy above the Fermi level, was determined in such a way as to reproduce the S pairing gaps ∆ q of NeuM and SNM as obtained from extended Brueckner-Hartree-Fock calculations including medium-polarizatione ff ects [13] (see Refs. [14–16] for details). But because of Coulomb e ff ects and a possible charge-symmetry breakingof nuclear forces, the proton pairing strength may be slightly di ff erent from the neutron pairing strength. This fine-tuning was made possible by introducing the parameters f ± q . Our neglect of the time-odd fields implicit in our use ofthe equal-filling approximation was compensated phenomenologically by allowing the pairing force to also dependon whether there is an even ( + ) or odd ( − ) number of nucleons of the charge type in question. By definition, f + n = E can be expressed as E = Z d r E ( r ) , (11)where the local energy density E ( r ) depends on the local normal and abnormal nucleon number densities, respectively, n q ( r ) = X σ = ± n q ( r , σ ; r , σ ) , (12) e n q ( r ) = X σ = ± e n q ( r , σ ; r , σ ) , (13)the kinetic-energy density τ q ( r ) = X σ = ± Z d r ′ δ ( r − r ′ ) ∇ · ∇ ′ n q ( r , σ ; r ′ , σ ) , (14)nd the spin-current vector density J q ( r ) = − i X σ,σ ′ = ± Z d r ′ δ ( r − r ′ ) ∇ n q ( r , σ ; r ′ , σ ′ ) × ˆ σ σ ′ σ = i X σ,σ ′ = ± Z d r ′ δ ( r − r ′ ) ∇ ′ n q ( r , σ ; r ′ , σ ′ ) × ˆ σ σ ′ σ . (15)Note that all the terms in J and J q from the energy density are dropped as discussed in Ref. [17]. Minimizing theenergy leads to the self-consistent Hartree-Fock-Bogoliubov (HFB) equations [7] X σ ′ = ± h q ( r ) σσ ′ ψ ( q )1 k ( r , σ ′ ) + ∆ q ( r ) ψ ( q )2 k ( r , σ ) = ( E k + µ q ) ψ ( q )1 k ( r , σ ) , ∆ q ( r ) ψ ( q )1 k ( r , σ ) − X σ ′ = ± h q ( r ) σσ ′ ψ ( q )2 k ( r , σ ′ ) = ( E k − µ q ) ψ ( q )2 k ( r , σ ) , (16)where E k denotes the quasiparticle energy, µ q the nucleon chemical potential, h q ( r ) σσ ′ = −∇ · δ E δτ q ( r ) ∇ δ σσ ′ + δ E δ n q ( r ) δ σσ ′ − i δ E δ J q ( r ) · ∇ × ˆ σ σσ ′ (17)represents the single-particle Hamiltonian, ∆ q ( r ) = δ E δ e n q ( r ) = v π q [ n n ( r ) , n p ( r )] e n q ( r ) (18)is a potential responsible for pairing. Expressions for these fields can be found in Refs. [10, 12]. Depending on thechoice of boundary conditions, Equation (16) can describe the atomic nuclei present in the outer crust of a neutronstar, the inhomogeneous nuclear matter constituting the inner crust or the homogeneous neutron-proton mixture inthe core. However, a reliable description of all these di ff erent nuclear phases requires a suitable adjustment of theparameters of the EDF. Brussels-Montreal energy-density functionals
The parameters of the Brussels-Montreal EDF series BSk22-BSk26 [18] were determined primarily by fitting to the2353 measured masses of atomic nuclei having Z , N ≥ E coll = E crankrot (cid:26) b tanh( c | β | ) + d | β | exp {− l ( | β | − β ) } (cid:27) , (19)in which E crankrot denotes the cranking-model value of the rotational correction and β the quadrupole deformation,while all other parameters were fitted freely. This correction was shown to be in good agreement with calculationsusing 5D collective Hamiltonian [17]. To the HFB energy, a phenomenological Wigner correction was also added E W = V W exp − λ N − ZA + V ′ W | N − Z | exp − AA . (20)This term contributes significantly only for light nuclei ( A < A ) or nuclei with N close to Z (see, e.g., Ref. [14] for adiscussion of the physical interpretation of these terms). The spurious centre-of-mass energy was removed followingan essentially exact procedure [21]. Moreover, a correction for the finite size of the proton was made to both thecharge radius and the energy [14]. Finally, Coulomb exchange for protons was dropped, thus simulating neglectede ff ects such as Coulomb correlations, charge-symmetry breaking, and vacuum polarization [22].o ensure reliable extrapolations to the highly neutron-rich and very dense interiors of neutron stars, the Brussels-Montreal EDFs were further constrained to reproduce the EoS of homogeneous NeuM, as calculated by many-bodytheory. Although the EoS is fairly well determined at subsaturation densities, it remains highly uncertain at super-saturation densities prevailing in the core of the most massive neutron stars. Two di ff erent EoS were considered: therather sti ff EoS labelled as ‘V18’ by Li and Schulze [23] and the softer EoS labelled as ‘A18 + δ v + UIX ∗ ’ by Akmal,Pandharipande, and Ravenhall [24]. The fit to nuclear masses along with these constraints do not lead to a uniquedetermination of the EDF. Expanding the energy per nucleon of infinite nuclear matter (INM) of density n = n (1 + ǫ )and charge asymmetry η = ( n n − n p ) / n about the equilibrium density n = n and η = e ( n , η ) = a v + J + L ǫ ! η +
118 ( K v + η K sym ) ǫ + · · · (21)the incompressibility coe ffi cient K v was further restricted to lie in the experimental range K v = ±
10 MeV [25].To achieve a good fit to nuclear masses with a root-mean-square deviation as low as 0 . − . ffi cient J from 29 to 32 MeV. The EDFs BSk22, BSk23,BSk24 and BSk25 were thus all fitted to the NeuM EoS of Ref. [23] while having J =
32, 31, 30 and 29 MeV,respectively. To assess the role of the NeuM EoS, the EDF BSk26 was fitted to the softer EoS of Ref. [24] with J =
30 MeV. Nuclear-matter properties for the EDFs are summarized in Table 1. The intermediate EDF BSk23will not be further considered here. As shown in Fig. 1, the variation of the symmetry energy S ( n ) with density n aspredicted by the Brussels-Montreal EDFs are compatible with experimental constraints from transport-model analysesof midperipheral heavy-ion collisions of Sn isotopes [26], from the analyses of isobaric-analog states and neutron-skin data [27], and from the electric dipole polarizability of Pb [28]. The symmetry energy is defined here as thedi ff erence between the energy per nucleon in NeuM and the energy per nucleon in SNM, S ( n ) = e NeuM ( n ) − e SNM ( n ) . (22)The EDFs mainly di ff er in their predictions for the symmetry energy at densities n > n , as shown in Fig. 2. TABLE 1.
Nuclear-matter properties for the Brussels-Montreal func-tionals.
BSk22 BSk23 BSk24 BSk25 BSk26 J [MeV] 32.0 31.0 30.0 29.0 30.0 L [MeV] 68.5 57.8 46.4 36.9 37.5 K v [MeV] 245.9 245.7 245.5 236.0 240.8 K sym [MeV] 13.0 -11.3 -37.6 -28.5 -135.6 UNIFIED EQUATIONS OF STATE OF COLD CATALYZED MATTERMethods
For densities greater than ∼ − nucleons fm − ( ∼ g cm − ) atoms are completely ionised by the pressure andelectrons form a degenerate Fermi gas. Since nuclear degeneracy likewise holds everywhere it follows that the coldcatalyzed matter hypothesis is valid throughout the star except in a thin layer at densities ρ . g cm − , wherethe atomic ionization and electron degeneracy can be incomplete. The properties of this outermost region have beenextensively discussed in Ref. [1]. Therefore, we calculated the outer crust only for densities ρ & g cm − mainlyas described in Ref. [29] except that we used complete expressions for the electron exchange and charge polarizatione ff ects (see Ref. [5] for further details). The equilibrium properties were determined by minimizing the Gibbs freeenergy per nucleon g at each given pressure P assuming that each crustal layer consists of a perfect body-centeredcubic crystal made of a single nuclear species (binary compounds may be present, but only at the interface betweenpure adjacent crustal layers, see Ref. [30]). For the masses of nuclei, we made use of experimental data from the 2016Atomic Mass Evaluation [31], supplemented by the very recent measurements of copper isotopes [32]. For nucleiwith experimentally unknown masses, we used the HFB-22, HFB-24, HFB-25 or HFB-26 tables [33]. Above somepressure P drip such that g = M n c , some neutrons drip out of nuclei marking the transition to the inner crust [34]. IGURE 1.
Variation of the symmetry energy S ( n ) with density n for the Brussels-Montreal functionals. The shaded areas areexperimental constraints: from heavy-ion collisions [26] (blue), from isobaric-analog states and neutron skins [27] (yellow), andfrom the electric dipole polarizability of Pb [28] (green).
Full HFB calculations of the inner crust properties would be computationally extremely costly due to the necessity toaccount for unbound states. For this reason, we adopted the fourth-order extended Thomas-Fermi method with protonshell and pairing corrections added perturbatively using the Strutinsky integral theorem [35, 36]. All details can befound in Ref. [5]. We recall only the main features here. Nuclear clusters were supposed to be spherical. To furtherspeed up the computations, the Wigner-Seitz (WS) cell was approximated by a sphere of radius R . Moreover, thenucleon density distributions in the WS cell were parameterized as the sum of a constant “background” term and a“cluster” term according to n q ( r ) = n Bq + n Λ q f q ( r ) , (23) f q ( r ) = + exp ((cid:18) C q − Rr − R (cid:19) − ) exp (cid:18) r − C q a q (cid:19) (24)(if n Λ q is negative the cluster becomes a bubble), where q = n , p for neutrons or protons respectively, while n B , q , n Λ , q , C q , and a q , are free parameters. This profile guarantees that the first three derivatives of the density vanish atthe cell surface, as required by the usual implementation of the ETF method (see Section II of Ref. [35], where otherrelevant details will be found). At high enough densities, free protons may appear. The onset of proton emission wasdetermined by the condition ~ M ∗ p ( r = h π n p ( r = i / > U p ( r = R ) − U p ( r = , (25) .2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 n [fm -3 ] S ( n ) [ M e V ] BSk22BSk24BSk25BSk26
FIGURE 2.
Variation of the symmetry energy S ( n ) at densities n > n for the Brussels-Montreal functionals. where M ∗ p ( r =
0) is the e ff ective proton mass at the center of the cell and U p ( r ) is the value of the proton single-particlepotential at the point r . In the free-proton region, the proton shell and pairing corrections were dropped entirely. Forconvenience, we minimized the energy per nucleon at fixed average baryon number density. As shown in Ref. [35],density discontinuities are negligibly small in the inner crust so that minimizing the energy per nucleon or the Gibbsfree energy per nucleon yields practically the same results without having recourse to a Maxwell construction. Thepressure P at any mean density ¯ n was calculated semi-analytically as described in Appendix B of [35]. Calculationsin the inner crust were performed using the same EDF as that underlying the HFB nuclear mass model used in theouter crust. Calculations in the core of a neutron star are comparatively much simpler since the energy density andthe pressure are given by analytic expressions (see Ref. [5]). The equilibrium composition (allowing for the presenceof electrons and muons) was determined by solving the equations of beta equilibrium under the condition of chargeneutrality using the same EDF as in the crust. Results
The composition of the outer crust is completely determined by experimentally measured masses of known nuclei upto ¯ n ≈ × − fm − . The use of the HFB models at higher densities introduces some uncertainties. Nevertheless,the di ff erences between models remain small. Indeed, the sequence of nuclei is essentially governed by nuclear-structure e ff ects (large regions containing nuclei with neutron magic numbers N =
50 and N = Y p = Z / A and the pressure P as a function of the mean density ¯ n . Transitions between adjacent layers are accompaniedby density discontinuities, which will be reduced (though not entirely removed) if binary compounds are allowed [30].In all cases, matter becomes progressively more neutron rich with increasing depth, a consequence of mechanicalequilibrium and the fact that the pressure arises mainly from the degenerate electron gas. Moreover, nuclei are less andess bound (see, e.g. Ref. [4]) until neutrons are emitted delimiting the boundary with the inner crust. The compositionbeyond this point is found to be more sensitive to the symmetry energy. As shown in Fig. 3, the inner crust tends to bemore neutron-rich with increasing J reflecting the anticorrelation between J and the symmetry energy S ( n ) of INMat subnuclear densities n < n . Comparing BSk24 and BSk26, which were fitted to the same value of J , shows thatthe choice of the constraining EoS of NeuM has no significant impact on the crustal composition. As can be seen inFig. 4, the equilibrium proton number Z eq is governed to a large extent by shell e ff ects, as in the outer crust. Apartfrom the well-known magic numbers 20 and 50 from atomic nuclei, the numbers 40, 58, 92, and 138 appear to bealso favored in this extremely neutron-rich environment. This change of nuclear structure mainly arises from the verystrong quenching of the spin-orbit coupling. Evidence for such a quenching have been also reported in exotic nuclei[37]. Not only does the shell structure strongly depend on the symmetry energy, but also the smooth ETF part of theenergy per nucleon. In the free proton region where the shell correction is dropped entirely, Z eq is thus found to beanticorrelated with J . Given that the proton fraction Y p varies smoothly it follows that the mean number of neutrons N eq per WS cell must display the same shell e ff ects as does the proton number Z eq . The variation of the pressure P with ¯ n in the inner crust, shown in Fig. 5 for our four EDFs, directly reflects the slope of the symmetry energy. Indeed,the pressure (essentially given by free neutrons) is approximately given by P ≈ n " K v ǫ n + η dSdn , (26)with η ∼
1. Since the EDFs were all fitted to the experimental value of the incompressibility K v , the di ff erences inthe pressure originate from the di ff erences in the slope of the symmetry energy. The coincidence of the slopes of thesymmetry energy at density ∼ .
44 fm − for BSk22, BSk24 and BSk25 leads to a crossing of the pressure curves.Because the slope of the symmetry energy for BSk26 is greater than that for BSk24 at any density below n , theBSk26 EDF yields a sti ff er EoS than BSk24 in the inner crust. For the same reason, the EoS obtained with BSk26 issofter in the core, while BSk22 leads to the sti ff est EoS. Comparing only the slopes of the symmetry energy betweenthe di ff erent EDFs, one could have expected an even softer EoS in the core for BSk26. However, the softer symmetryenergy is to a large extent compensated by the higher neutron abundance such that the pressure (26) in the core is notmuch di ff erent. A remarkable feature we found is that the proton fraction Y p ( n ) increases with density in the core forall EDFs, in contrast to what happens in both the outer and inner crusts.Complete numerical results and analytic fits applicable to the entire star can be found in Ref. [5]. All these fitshave been implemented in Fortran subroutines, freely available at URL http: // ff e.ru / astro / NSG / BSk / . GLOBAL STRUCTURE OF NEUTRON STARS
The global structure of a non-rotating neutron star was computed from the Tolman-Oppenheimer-Volko ff (TOV)equations [38, 39] d P ( r )d r = − G E ( r ) M ( r ) c r " + P ( r ) E ( r ) + π P ( r ) r c M ( r ) − G M ( r ) c r − , (27)where G is the gravitational constant, and M ( r ) = π c Z r E ( r ′ ) r ′ d r ′ . (28)Here E ( r ) is the mass-energy density of matter at the radial coordinate r . The gravitational mass of the star is given by M NS = M ( R NS ), where R NS is the radial coordinate at which the radiative surface is reached.The mass-radius relation for the EoSs BSk22, BSk24, BSk25, and BSk26 are shown in Figure 6. As expected,the maximum mass is mainly determined by the NeuM constraint; BSk26 being fitted to a softer NeuM EoS than theother EDFs thus yields the lowest maximum mass. All our models are compatible with the measured masses of themost massive neutron stars known [40–42]. The figure also shows some constraints inferred from the analysis of thegravitational-wave signal GW170817 from a binary neutron-star merger. In particular, all models are compatible withthe maximum-mass limits derived in Ref. [43] under the assumption that the binary merger product collapsed to ablack hole with a mass very close to the mass-shedding limit. Quasi-universal relations between the maximum massof nonrotating stellar models and the maximum mass supported through uniform rotation were also used. Also shownin Figure 6 are the tighter limits of Refs. [44, 45]. All our models are consistent with the constraint of Ref. [45], butSk22 and BSk24 appear to be disfavored by that of Ref. [44]. However, these two constraints may be prone to modeluncertainties as they rely on numerical simulations with further assumptions on the gamma-ray burst.Figure 6 shows that for a given constraining EoS of NeuM the radius is correlated with the slope L of thesymmetry energy, as previously found, see e.g. Refs. [46, 47]. Lower bounds on the neutron-star radii have beenobtained in Refs. [48, 49] from the assumption of delayed collapse using an empirical relation for the threshold binarymass for prompt collapse. All our models satisfy these robust constraints, especially the most stringent one fromRef. [49]. Accounting for the observational estimate of the tidal deformability of a neutron star leads to a narrowerrange of allowed radii [50–54]. The radii predicted by all our models fall below the upper limit, but BSk26 appearto agree marginally with the lower limit obtained in Ref. [54] at 2 σ confidence level. Additional constraints can beobtained from the requirement that the extremely powerful direct Urca processes of neutrino emission operate ina small number of neutron stars (see, e.g., [55, 56]). The direct Urca processes are allowed in su ffi ciently massiveneutron stars. The low value of the threshold mass M DU ≃ . M ⊙ given by the BSk22 EDF implies that theseprocesses would occur in the majority of neutron stars, thus leading to their rapid cooling, which does not agreewith observations (see, e.g., Ref. [57]) unless these processes are suppressed by nuclear pairing in low-mass neutronstars. With threshold masses M DU ≃ . M ⊙ and M DU ≃ . M ⊙ , the BSk24 and BSk25 EDFs respectively, areconsistent with neutron-star cooling observations. On the contrary, the direct Urca reactions are forbidden for allhydrostatically stable neutron stars by the BSk26 EDF, which must therefore be ruled out. CONCLUSIONS
We have computed a series of four di ff erent unified EoS for the cold dense matter constituting the interior of non ac-creted neutron stars within the framework of the nuclear EDF theory using the accurately calibrated Brussels-MontrealEDFs BSk22, BSk24, BSk25 and BSk26 [18]. These EDFs were precision fitted to essentially all the available atomicmass data (with Z , N ≥
8) by using the HFB method. In addition, these functionals were constrained to fit, up to thedensities prevailing in neutron-star cores, two di ff erent realistic EoS of NeuM: the sti ff EoS ‘V18’ of Ref.[23] forBSk22, BSk24 and BSk25, and the softer EoS ‘A18 + δ v + UIX ∗ ’ of Ref.[24] for BSk26. The BSk22, BSk24 andBSk25 EDFs mainly di ff er in their values of the symmetry-energy coe ffi cient J , fixed to 32, 30 and 29 MeV, respec-tively. Good-quality mass fits could not be achieved outside this range, the optimum fit being obtained for J = ff erent methods to implement the EDF equations in the di ff erent regions. Whereas the HFBmethod was used in the outer crust (except when the appropriate atomic mass data were available) and in the core, thecomputationally much faster ETFSI method was followed in the inner crust. The composition and the EoS of the innercrust and core of a neutron star are found to be very sensitive to the symmetry energy S (¯ n ) at densities ¯ n below andabove normal density n respectively. Full numerical results as well as analytic fits valid over the whole star for the allneutron-star masses are given in Ref. [5]. While the maximum neutron-star mass is mainly determined by the NeuMconstraint, the radius is found to be correlated with the slope L of the symmetry energy for a given NeuM EoS. Allour EoS are compatible with the robust radius constraints of Refs. [48, 49] inferred from the analysis of GW170817.However, the tighter constraints obtained from the estimated tidal deformability of neutron stars disfavor BSk26. Thissame model also predicts the absence of the direct Urca processes in stable neutron stars, which is di ffi cult to reconcilewith observations of some neutron stars that apparently cool via such processes. ACKNOWLEDGMENTS
The work of N. C. was partially supported by the Fonds de la Recherche Scientifique - FNRS (Belgium) under grantn ◦ CDR-J.0187.16 and CDR-J.0115.18. J. M. P. acknowledges the partial support of the NSERC (Canada). The workof A. Y. P. was supported by the RFBR grant 16-29-13009-ofi-m. S. G. acknowledges the support of Fonds de laRecherche Scienti ffi que-FNRS (Belgium). This work was also partially supported by the European Cooperation inScience and Technology (COST) Actions MP1304 and CA16214. REFERENCES [1] P. Haensel, A. Y. Potekhin, and D. G. Yakovlev,
Neutron Stars 1: Equation of state and structure , Springer(2007).2] B. K. Harrison, and J. A. Wheeler, in
Onzi`eme Conseil de Physique Solvay , Stoops, Bruxelles, Belgium(1958).[3] B. K. Harrison, K. S. Thorne, M. Wakano, and J. A. Wheeler,
Gravitation Theory and Gravitational Collapse ,The University of Chicago Press (1965).[4] D. Blaschke and N. Chamel, Phases of dense matter in compact stars. In: L. Rezzolla, P. Pizzochero, D.Jones, N. Rea, I. Vida˜na (eds) The Physics and Astrophysics of Neutron Stars. Astrophysics and SpaceScience Library , pp.337-400 (Springer, 2018).[5] J. M. Pearson, N. Chamel, A. Y. Potekhin, A. F. Fantina, C. Ducoin, A. K. Dutta, S. Goriely, Mon. Not. Roy.Astron. Soc. , 2994 (2018); Erratum: Mon. Not. Roy. Astron. Soc. , 768 (2019).[6] T. Duguet, Lect. Notes Phys. (Springer-Verlag Berlin Heidelberg, 2014), pp 293-350.[7] J. Dobaczewski, H. Flocard, J. Treiner, Nucl. Phys.
A422 , 103 (1984).[8] M. Bender, P.-H. Heenen and P.-G. Reinhard, Rev. Mod. Phys. , 121 (2003).[9] N. Chamel, Phys. Rev. C , 061307(R) (2010).[10] N. Chamel, S. Goriely, J. M. Pearson, Phys. Rev. C , 065804 (2009).[11] N. Chamel, S. Goriely, Phys. Rev. C , 045804 (2010).[12] N. Chamel, Phys. Rev. C , 014313 (2010).[13] L.G. Cao, U. Lombardo, and P. Schuck, Phys. Rev. C , 064301 (2006).[14] N. Chamel, S. Goriely, and J. M. Pearson, Nucl. Phys. A812 , 72 (2008).[15] S. Goriely, N. Chamel, and J. M. Pearson, Phys. Rev. Lett. , 152503 (2009).[16] S. Goriely, N. Chamel, and J. M. Pearson, Eur. Phys.
A42 , 547 (2009).[17] S. Goriely, N. Chamel, and J. M. Pearson, Phys. Rev. C , 035804 (2010).[18] S. Goriely, N. Chamel, J. M. Pearson, Phys. Rev. C , 024308 (2013).[19] G. Audi, M. Wang, A. H. Wapstra, F. G. Kondev, M. MacCormick, X. Xu, B. Pfei ff er, Chin. Phys. C ,1287 (2012).[20] M. Samyn, S. Goriely, P.-H. Heenen, J. M. Pearson, and F. Tondeur, Nucl. Phys. A700 , 142 (2002).[21] S. Goriely, M. Samyn, M. Bender, and J. M. Pearson, Phys. Rev. C , 054325 (2003).[22] S. Goriely and J. M. Pearson, Phys. Rev. C , 031301(R) (2008).[23] Z. H. Li and H.-J. Schulze, Phys. Rev. C , 028801 (2008).[24] A. Akmal, V. R. Pandharipande, and D. G. Ravenhall, Phys. Rev. C , 1804 (1998).[25] G. Col`o, N. V. Giai, J. Meyer, K. Bennaceur, and P. Bonche, Phys. Rev. C , 024307 (2004).[26] M. B. Tsang et al., Phys. Rev. Lett. , 122701 (2009).[27] P. Danielewicz and J. Lee, Nucl. Phys. A922 , 1 (2014).[28] Zhen Zhang and Lie-Wen Chen, Phys. Rev. C , 031301(R) (2015).[29] J. M. Pearson, S. Goriely, N. Chamel, Phys. Rev. C , 065810 (2011).[30] N. Chamel and A. F. Fantina, Phys. Rev. C , 065802 (2016).[31] M. Wang, G. Audi, F. G. Kondev, W. J. Huang, S. Naimi, X. Xu, Chin. Phys. C , 030003 (2017).[32] A. Welker et al., Phys. Rev. Lett. , 192502 (2017).[33] http: // / bruslib / [34] N. Chamel, A. F. Fantina, J.-L. Zdunik, P. Haensel, Phys. Rev. C , 055803 (2015).[35] J. M. Pearson, N. Chamel, S. Goriely, C. Ducoin, Phys. Rev. C , 065803 (2012).[36] J. M. Pearson, N. Chamel, A. Pastore, S. Goriely, Phys. Rev. C , 018801 (2015).[37] J. P. Schi ff er et al., Phys. Rev. Lett. , 162501 (2004).[38] R. C. Tolman, Phys. Rev. , 364 (1939).[39] J. R. Oppenheimer and G. M. Volko ff , Phys. Rev. , 374 (1939).[40] P. B. Demorest, T. Pennucci, S. M. Ransom, M. S. E. Roberts, J. W. T. Hessels, Nature , 1081 (2010).[41] J. Antoniadis, P. C. C. Freire, N. Wex, T. M. Tauris, R. S. Lynch, M. H. van Kerkwijk, M. Kramer, C. Bassa,V. S. Dhillon, T. Driebe, J. W. T. Hessels, V. M. Kaspi, V. I. Kondratiev, N. Langer, T. R. Marsh, M. A.McLaughlin, T. T. Pennucci, S. M. Ransom, I. H. Stairs, J. van Leeuwen, J. P. W. Verbiest, D. G. Whelan,Science , 1233232 (2013).[42] E. Fonseca et al, Astrophys. J. , 167 (2016).[43] L. Rezzolla, E. R. Most, L. R. Weih, Astrophys. J. , L25 (2018).[44] M. Shibata et al., Phys. Rev. D , 123012 (2017).[45] M. Ruiz, S. L. Shapiro, A. Tsokaros, Phys. Rev. D , 021501(R) (2018).[46] M. Fortin et al., Phys. Rev. C , 035804 (2016).[47] J. Margueron, R. Ho ff mann Casali, F. Gulminelli, Phys. Rev. C , 025806 (2018).[48] A. Bauswein, O. Just, H.-T. Janka, N. Stergioulas, Astrophys. J. , L34 (2017).49] S. K¨oppel, L. Bovard, L. Rezzolla, Astrophys. J. , L16 (2019).[50] B. P. Abbott et al., Phys. Rev. Lett. , 161101 (2018).[51] E. Annala, T. Gorda, A. Kurkela, A. Vuorinen, Phys. Rev. Lett. , 172703 (2018).[52] F. J. Fattoyev, J. Piekarewicz, C. J. Horowitz, Phys. Rev. Lett. , 172702 (2018).[53] S. De, D. Finstad, J. M. Lattimer, D. A. Brown, E. Berger, C. M. Biwer, Phys. Rev. Lett. , 091102 (2018).[54] E.R. Most et al., Phys. Rev. Lett. , 261103 (2018).[55] D. G. Yakovlev, K. P. Levenfish, A. Y. Potekhin, O. Y. Gnedin, G. Chabrier, Astron. Astrophys. , 169(2004).[56] E. F. Brown, A. Cumming, F. J. Fattoyev, C. J. Horowitz, D. Page, S. Reddy, Phys. Rev. Lett. , 182701(2018).[57] A. Y. Potekhin and G. Chabrier, Astron. Astrophys. , A74 (2018). IGURE 3.
Upper panel: Curves show the computed equilibrium proton fraction Y p = Z eq / A in the inner crust as a function ofmean baryon density ¯ n for our four functionals; arrows indicate onset of proton drip. The inset shows the cross-over between BSk24and BSk26 just before proton drip. Lower panel: Deviations between the computed data and the fitted analytic function ( ∆ Y p = fit − data). IGURE 4.
Equilibrium value Z eq of number of protons in inner crust as a function of mean baryon density ¯ n for our fourfunctionals. For clarity only every second point is shown. Arrows indicate onset of proton drip. The dotted curves relate to thevalues of Z eq calculated in the ETF approximation. IGURE 5.
Upper panel: Curves show the computed pressure P in the inner crust as a function of mean baryon density ¯ n for ourfour functionals; arrows indicate onset of proton drip. Lower panel: Fractional deviations between the computed data and the fittedanalytic function ( ∆ P = fit − data). IGURE 6.