Unified equations of state for cold non-accreting neutron stars with Brussels-Montreal functionals. I. Role of symmetry energy

TL;DR

This paper develops unified equations of state for cold neutron stars using nuclear energy-density functionals, highlighting how uncertainties in symmetry energy influence their structure and composition.

Contribution

It introduces a comprehensive, thermodynamically consistent framework for modeling neutron stars with multiple functionals, emphasizing the role of symmetry energy.

Findings

Symmetry energy uncertainties significantly affect neutron star composition.

Analytic fits provided for practical astrophysical modeling.

Different functionals yield consistent pressure-density relations.

Abstract

The theory of the nuclear energy-density functional is used to provide a unified and thermodynamically consistent treatment of all regions of cold non-accreting neutron stars. In order to assess the impact of our lack of complete knowledge of the density dependence of the symmetry energy on the constitution and the global structure of neutron stars, we employ four different functionals. All of them were precision fitted to essentially all the nuclear-mass data with the Hartree-Fock-Bogoliubov method and two different neutron-matter equations of state based on realistic nuclear forces. For each functional, we calculate the composition, the pressure-density relation, and the chemical potentials throughout the star. We show that uncertainties in the symmetry energy can significantly affect the theoretical results for the composition and global structure of neutron stars. To facilitate astrophysical applications, we construct analytic fits to our numerical results.

Tables

Table 1. Table 1: J 𝐽 J , L 𝐿 L , K v subscript 𝐾 𝑣 K_{v} and K sym subscript 𝐾 sym K_{\mathrm{sym}} coefficients for the functionals of this paper (defined in Eq. ( 2 )).

	BSk22	BSk24	BSk25	BSk26
$J$ [MeV]	32.0	30.0	29.0	30.0
$L$ [MeV]	68.5	46.4	36.9	37.5
$K_v$ [MeV]	245.9	245.5	236.0	240.8
$K_{\mathrm{sym}}$ [MeV]	13.0	-37.6	-28.5	-135.6

Table 2. Table 2: Root-mean-square ( σ 𝜎 \sigma ) deviations between the experimental nuclear mass data (AME 2016, Wang et al. 2017 ) and the predictions for the models based on the functionals of this paper. The first line refers to the 2408 measured masses M 𝑀 M with Z 𝑍 Z and N ≥ 𝑁 absent N\geq 8, and the second to the masses M nr subscript 𝑀 nr M_{\mathrm{nr}} of the subset of 286 neutron-rich nuclei (neutron separation energy ≤ \leq 5.0 MeV).

	HFB-22	HFB-24	HFB-25	HFB-26
$\sigma \left(M\right)$ [MeV]	0.648	0.565	0.556	0.580
$\sigma \left(M_{\mathrm{nr}}\right)$ [MeV]	0.904	0.781	0.829	0.811

Table 3. Table 3: Composition and EoS of the outer crust of a cold non-accreting neutron star as obtained using experimental data from the 2016 Atomic Mass Evaluation (Wang et al., 2017 ) (above the dotted line) supplemented with the nuclear mass model HFB-22 (Goriely et al., 2013 ) . We have also made use of the very recent measurements of copper isotopes (Welker et al., 2017 ) . In the table are listed: the atomic number Z 𝑍 Z , the neutron number N 𝑁 N , and the mass number A 𝐴 A of each nuclide, the minimum and maximum baryon number densities n ¯ min subscript ¯ 𝑛 min \bar{n}_{\mathrm{min}} and n ¯ max subscript ¯ 𝑛 max \bar{n}_{\mathrm{max}} at which the nuclide is present, the corresponding maximum pressure P max subscript 𝑃 max P_{\mathrm{max}} , and the neutron, proton, and electron chemical potentials at density n ¯ max subscript ¯ 𝑛 max \bar{n}_{\mathrm{max}} . Note that the neutron and proton chemical potentials are given with the rest mass subtracted out. The baryon number densities are measured in units of fm -3 , pressures are given in units of MeV fm -3 and chemical potentials are given in units of MeV. See text for details.

$Z$	$N$	$A$	${\overline{n}}_{\min }$	${\overline{n}}_{\max }$	$P_{\max }$	${\mu }_n-M_n{c}^2$	${\mu }_p-M_p{c}^2$	${\mu }_e$
26	30	56	$-$	$4.93\times {10}^{-9}$	$3.36\times {10}^{-10}$	$-8.96$	$-8.62$	0.95
28	34	62	$5.08\times {10}^{-9}$	$1.63\times {10}^{-7}$	$4.34\times {10}^{-8}$	$-8.25$	$-9.56$	2.61
28	36	64	$1.68\times {10}^{-7}$	$8.01\times {10}^{-7}$	$3.56\times {10}^{-7}$	$-7.53$	$-10.57$	4.33
28	38	66	$8.28\times {10}^{-7}$	$8.79\times {10}^{-7}$	$3.87\times {10}^{-7}$	$-7.49$	$-10.62$	4.42
36	50	86	$8.98\times {10}^{-7}$	$1.87\times {10}^{-6}$	$1.04\times {10}^{-6}$	$-7.00$	$-11.36$	5.65
34	50	84	$1.94\times {10}^{-6}$	$6.83\times {10}^{-6}$	$5.62\times {10}^{-6}$	$-5.87$	$-13.16$	8.58
32	50	82	$7.09\times {10}^{-6}$	$1.67\times {10}^{-5}$	$1.78\times {10}^{-5}$	$-4.81$	$-14.94$	11.43
30	50	80	$1.74\times {10}^{-5}$	$2.97\times {10}^{-5}$	$3.63\times {10}^{-5}$	$-4.01$	$-16.37$	13.65
28	48	76	$3.03\times {10}^{-5}$	$3.72\times {10}^{-5}$	$4.79\times {10}^{-5}$	$-3.67$	$-16.99$	14.62
28	50	78	$3.82\times {10}^{-5}$	$4.55\times {10}^{-5}$	$6.07\times {10}^{-5}$	$-3.36$	$-17.58$	15.51
28	52	80	$4.68\times {10}^{-5}$	$7.79\times {10}^{-5}$	$1.20\times {10}^{-4}$	$-2.40$	$-19.49$	18.39
42	82	124	$8.16\times {10}^{-5}$	$1.12\times {10}^{-4}$	$1.84\times {10}^{-4}$	$-1.73$	$-20.96$	20.52
40	82	122	$1.16\times {10}^{-4}$	$1.34\times {10}^{-4}$	$2.23\times {10}^{-4}$	$-1.42$	$-21.64$	21.52
39	82	121	$1.36\times {10}^{-4}$	$1.59\times {10}^{-4}$	$2.76\times {10}^{-4}$	$-1.06$	$-22.46$	22.69
38	84	122	$1.65\times {10}^{-4}$	$2.16\times {10}^{-4}$	$3.97\times {10}^{-4}$	$-0.42$	$-23.96$	24.84
38	86	124	$2.20\times {10}^{-4}$	$2.25\times {10}^{-4}$	$4.09\times {10}^{-4}$	$-0.37$	$-24.10$	25.02
38	88	126	$2.29\times {10}^{-4}$	$2.59\times {10}^{-4}$	$4.83\times {10}^{-4}$	$-0.06$	$-24.86$	26.09
38	90	128	$2.64\times {10}^{-4}$	$2.69\times {10}^{-4}$	$4.99\times {10}^{-4}$	0.00	$-25.01$	26.30

Table 4. Table 4: Same as in Table  3 , for the nuclear mass model HFB-24.

$Z$	$N$	$A$	${\overline{n}}_{\min }$	${\overline{n}}_{\max }$	$P_{\max }$	${\mu }_n-M_n{c}^2$	${\mu }_p-M_p{c}^2$	${\mu }_e$
26	30	56	$-$	$4.93\times {10}^{-9}$	$3.36\times {10}^{-10}$	$-8.96$	$-8.62$	0.95
28	34	62	$5.08\times {10}^{-9}$	$1.63\times {10}^{-7}$	$4.34\times {10}^{-8}$	$-8.25$	$-9.56$	2.61
28	36	64	$1.68\times {10}^{-7}$	$8.01\times {10}^{-7}$	$3.56\times {10}^{-7}$	$-7.53$	$-10.57$	4.33
28	38	66	$8.28\times {10}^{-7}$	$8.79\times {10}^{-7}$	$3.87\times {10}^{-7}$	$-7.49$	$-10.62$	4.42
36	50	86	$8.98\times {10}^{-7}$	$1.87\times {10}^{-6}$	$1.04\times {10}^{-6}$	$-7.00$	$-11.36$	5.65
34	50	84	$1.94\times {10}^{-6}$	$6.83\times {10}^{-6}$	$5.62\times {10}^{-6}$	$-5.87$	$-13.16$	8.58
32	50	82	$7.09\times {10}^{-6}$	$1.67\times {10}^{-5}$	$1.78\times {10}^{-5}$	$-4.81$	$-14.94$	11.43
30	50	80	$1.74\times {10}^{-5}$	$3.47\times {10}^{-5}$	$4.45\times {10}^{-5}$	$-3.76$	$-16.83$	14.36
28	50	78	$3.62\times {10}^{-5}$	$6.63\times {10}^{-5}$	$1.00\times {10}^{-4}$	$-2.65$	$-18.93$	17.57
28	52	80	$6.81\times {10}^{-5}$	$7.54\times {10}^{-5}$	$1.15\times {10}^{-4}$	$-2.44$	$-19.34$	18.18
42	82	124	$7.89\times {10}^{-5}$	$1.22\times {10}^{-4}$	$2.07\times {10}^{-4}$	$-1.52$	$-21.36$	21.13
40	82	122	$1.27\times {10}^{-4}$	$1.58\times {10}^{-4}$	$2.79\times {10}^{-4}$	$-1.01$	$-22.47$	22.75
39	82	121	$1.61\times {10}^{-4}$	$1.64\times {10}^{-4}$	$2.88\times {10}^{-4}$	$-0.95$	$-22.59$	22.93
38	82	120	$1.68\times {10}^{-4}$	$1.95\times {10}^{-4}$	$3.54\times {10}^{-4}$	$-0.59$	$-23.44$	24.14
38	84	122	$1.99\times {10}^{-4}$	$2.39\times {10}^{-4}$	$4.54\times {10}^{-4}$	$-0.13$	$-24.53$	25.69
38	86	124	$2.44\times {10}^{-4}$	$2.56\times {10}^{-4}$	$4.87\times {10}^{-4}$	0.00	$-24.85$	26.14

Table 5. Table 5: Same as in Table  3 , for the nuclear mass model HFB-25.

$Z$	$N$	$A$	${\overline{n}}_{\min }$	${\overline{n}}_{\max }$	$P_{\max }$	${\mu }_n-M_n{c}^2$	${\mu }_p-M_p{c}^2$	${\mu }_e$
26	30	56	$-$	$4.93\times {10}^{-9}$	$3.36\times {10}^{-10}$	$-8.96$	$-8.62$	0.95
28	34	62	$5.08\times {10}^{-9}$	$1.63\times {10}^{-7}$	$4.34\times {10}^{-8}$	$-8.25$	$-9.56$	2.61
28	36	64	$1.68\times {10}^{-7}$	$8.01\times {10}^{-7}$	$3.56\times {10}^{-7}$	$-7.53$	$-10.57$	4.33
28	38	66	$8.28\times {10}^{-7}$	$8.79\times {10}^{-7}$	$3.87\times {10}^{-7}$	$-7.49$	$-10.62$	4.42
36	50	86	$8.98\times {10}^{-7}$	$1.87\times {10}^{-6}$	$1.04\times {10}^{-6}$	$-7.00$	$-11.36$	5.65
34	50	84	$1.94\times {10}^{-6}$	$6.83\times {10}^{-6}$	$5.62\times {10}^{-6}$	$-5.87$	$-13.16$	8.58
32	50	82	$7.09\times {10}^{-6}$	$1.67\times {10}^{-5}$	$1.78\times {10}^{-5}$	$-4.81$	$-14.94$	11.43
30	50	80	$1.74\times {10}^{-5}$	$3.28\times {10}^{-5}$	$4.13\times {10}^{-5}$	$-3.85$	$-16.66$	14.10
28	50	78	$3.42\times {10}^{-5}$	$7.46\times {10}^{-5}$	$1.17\times {10}^{-4}$	$-2.41$	$-19.40$	18.28
44	82	126	$7.78\times {10}^{-5}$	$7.84\times {10}^{-5}$	$1.19\times {10}^{-4}$	$-2.39$	$-19.49$	18.39
42	82	124	$8.09\times {10}^{-5}$	$1.29\times {10}^{-4}$	$2.23\times {10}^{-4}$	$-1.39$	$-21.61$	21.51
40	82	122	$1.34\times {10}^{-4}$	$1.68\times {10}^{-4}$	$3.02\times {10}^{-4}$	$-0.86$	$-22.78$	23.21
39	82	121	$1.71\times {10}^{-4}$	$1.71\times {10}^{-4}$	$3.04\times {10}^{-4}$	$-0.85$	$-22.80$	23.24
38	82	120	$1.75\times {10}^{-4}$	$2.12\times {10}^{-4}$	$3.94\times {10}^{-4}$	$-0.38$	$-23.89$	24.80
38	84	122	$2.16\times {10}^{-4}$	$2.50\times {10}^{-4}$	$4.83\times {10}^{-4}$	0.00	$-24.79$	26.08

Table 6. Table 6: Same as in Table  3 , for the nuclear mass model HFB-26.

$Z$	$N$	$A$	${\overline{n}}_{\min }$	${\overline{n}}_{\max }$	$P_{\max }$	${\mu }_n-M_n{c}^2$	${\mu }_p-M_p{c}^2$	${\mu }_e$
26	30	56	$-$	$4.93\times {10}^{-9}$	$3.36\times {10}^{-10}$	$-8.96$	$-8.62$	0.95
28	34	62	$5.08\times {10}^{-9}$	$1.63\times {10}^{-7}$	$4.34\times {10}^{-8}$	$-8.25$	$-9.56$	2.61
28	36	64	$1.68\times {10}^{-7}$	$8.01\times {10}^{-7}$	$3.56\times {10}^{-7}$	$-7.53$	$-10.57$	4.33
28	38	66	$8.28\times {10}^{-7}$	$8.79\times {10}^{-7}$	$3.87\times {10}^{-7}$	$-7.49$	$-10.62$	4.42
36	50	86	$8.98\times {10}^{-7}$	$1.87\times {10}^{-6}$	$1.04\times {10}^{-6}$	$-7.00$	$-11.36$	5.65
34	50	84	$1.94\times {10}^{-6}$	$6.83\times {10}^{-6}$	$5.62\times {10}^{-6}$	$-5.87$	$-13.16$	8.58
32	50	82	$7.09\times {10}^{-6}$	$1.67\times {10}^{-5}$	$1.78\times {10}^{-5}$	$-4.81$	$-14.94$	11.43
30	50	80	$1.74\times {10}^{-5}$	$3.56\times {10}^{-5}$	$4.62\times {10}^{-5}$	$-3.71$	$-16.91$	14.49
28	50	78	$3.72\times {10}^{-5}$	$5.91\times {10}^{-5}$	$8.59\times {10}^{-5}$	$-2.88$	$-18.49$	16.91
28	52	80	$6.07\times {10}^{-5}$	$7.55\times {10}^{-5}$	$1.15\times {10}^{-4}$	$-2.44$	$-19.35$	18.20
42	82	124	$7.91\times {10}^{-5}$	$1.21\times {10}^{-4}$	$2.03\times {10}^{-4}$	$-1.55$	$-21.29$	21.04
40	82	122	$1.25\times {10}^{-4}$	$1.51\times {10}^{-4}$	$2.62\times {10}^{-4}$	$-1.13$	$-22.23$	22.39
40	84	124	$1.53\times {10}^{-4}$	$1.72\times {10}^{-4}$	$3.07\times {10}^{-4}$	$-0.85$	$-22.86$	23.30
38	82	120	$1.76\times {10}^{-4}$	$1.79\times {10}^{-4}$	$3.16\times {10}^{-4}$	$-0.80$	$-22.97$	23.46
38	84	122	$1.83\times {10}^{-4}$	$2.28\times {10}^{-4}$	$4.25\times {10}^{-4}$	$-0.26$	$-24.24$	25.27
38	86	124	$2.32\times {10}^{-4}$	$2.51\times {10}^{-4}$	$4.73\times {10}^{-4}$	$-0.06$	$-24.73$	25.96
38	88	126	$2.55\times {10}^{-4}$	$2.61\times {10}^{-4}$	$4.90\times {10}^{-4}$	0.00	$-24.88$	26.18

Table 7. Table 7: Parameters of the neutron drip point for the functionals of this paper.

	BSk22	BSk24	BSk25	BSk26
$J$ [MeV]	32	30	29	30
${\overline{n}}_{\mathrm{nd}}$ [$\times {10}^{-4}$ fm-3]	$2.69$	$2.56$	$2.50$	$2.61$
$Z$	38	38	38	38
$N$	90	86	84	88

Table 8. Table 8: Principal values of Z eq subscript 𝑍 eq Z_{\mathrm{eq}} in the inner crust below proton drip for BSk22.

$Z_{\mathrm{eq}}$	${\overline{n}}_{\min }$ [fm-3]	${\overline{n}}_{\max }$ [fm-3]
40	2.69$\times {10}^{-4}$	0.0340
20	0.0350	0.0533
40	0.0544	0.591

Table 9. Table 9: Principal value of Z eq subscript 𝑍 eq Z_{\mathrm{eq}} in the inner crust below proton drip for BSk24 (corrected as per the Erratum).

$Z_{\mathrm{eq}}$	${\overline{n}}_{\min }$ [fm-3]	${\overline{n}}_{\max }$ [fm-3]
40	2.56$\times {10}^{-4}$	0.0715

Table 10. Table 10: Principal values of Z eq subscript 𝑍 eq Z_{\mathrm{eq}} in the inner crust below proton drip for BSk25 (corrected as per the Erratum).

$Z_{\mathrm{eq}}$	${\overline{n}}_{\min }$ [fm-3]	${\overline{n}}_{\max }$ [fm-3]
40	2.50$\times {10}^{-4}$	2.70$\times {10}^{-4}$
50	2.70$\times {10}^{-4}$	0.0138
58	0.0139	0.0474
92	0.0478	0.07663
138	0.0768	0.0806

Table 11. Table 11: Principal value of Z eq subscript 𝑍 eq Z_{\mathrm{eq}} in the inner crust below proton drip for BSk26.

$Z_{\mathrm{eq}}$	${\overline{n}}_{\min }$ [fm-3]	${\overline{n}}_{\max }$ [fm-3]
40	2.61$\times {10}^{-4}$	0.0730

Table 12. Table 12: Onset of proton drip in the inner crust (corrected as per the Erratum).

	BSk22	BSk24	BSk25	BSk26
$\overline{n}$ [fm-3]	0.059	0.073	0.081	0.074

Table 13. Table 13: Comparison of inner-crust and outer-crust codes at the neutron drip-point density n ¯ nd subscript ¯ 𝑛 nd \bar{n}_{\mathrm{nd}} of the functional in question; results for latter code in parentheses (except for the neutron-drip density). e eq subscript 𝑒 eq e_{\mathrm{eq}} is the internal energy per nucleon (with the neutron rest-mass energy M n ​ c 2 subscript 𝑀 n superscript 𝑐 2 M_{\mathrm{n}}c^{2} subtracted), and P 𝑃 P the pressure.

Force	${\overline{n}}_{\mathrm{nd}}$ (fm-3)	$Z$	$N$	$e_{\mathrm{eq}}$ (MeV)	$P$ (MeV fm-3)
BSk22	$2.69\times {10}^{-4}$	40 (38)	98 (90)	$-1.78$ ($-1.85$)	$5.10\times {10}^{-4}$ ($4.99\times {10}^{-4}$)
BSk24	$2.56\times {10}^{-4}$	40 (38)	94 (86)	$-1.83$ ($-1.90$)	$4.95\times {10}^{-4}$ ($4.87\times {10}^{-4}$)
BSk25	$2.50\times {10}^{-4}$	40 (38)	93 (84)	$-1.85$ ($-1.93$)	$4.85\times {10}^{-4}$ ($4.83\times {10}^{-4}$)
BSk26	$2.61\times {10}^{-4}$	40 (38)	95 (88)	$-1.80$ ($-1.87$)	$5.04\times {10}^{-4}$ ($4.90\times {10}^{-4}$)

Table 14. Table 14: Parameters relating to the crust-core transition.

Functional	$n_{\mathrm{cc}}$ (fm-3)	$Y_{e,\mathrm{cc}}$	$P_{cc}$ (MeV fm-3)	$n_{\mathrm{cc}}^{\prime }$ (fm-3)
BSk22	0.0716068	0.028294	0.290934	0.0705000
BSk24	0.0807555	0.033671	0.267902	0.0790740
BSk25	0.0855534	0.035829	0.210878	0.0852758
BSk26	0.0849477	0.035721	0.363049	0.0835000

Table 15. Table 15: Highest densities in inner crust beyond which instabilities set in.

Functional	$\overline{n}$ (fm-3)	$Z_{\mathrm{eq}}$	$e_{\mathrm{eq}}$ (MeV)	$\Lambda$	$P$ (MeV fm-3)
BSk22	0.0635018	28.5	6.75925	0.02798	0.23010
BSk24	0.0749994	64.6	8.23139	0.02504	0.24594
BSk25	0.0832615	139.1	8.88852	0.01830	0.22250
BSk26	0.0786984	56.5	8.56597	0.02355	0.32862

Table 16. Table 16: Hydrostatic stability limits

EoS	$M_{\mathrm{NS}}^{\mathrm{stab}}$	$R_{\mathrm{NS}}^{\mathrm{stab}}$	$n_c^{\mathrm{stab}}$	${\rho }_c^{\mathrm{stab}}$
	$\left(M_{\odot }\right)$	(km)	(fm-3)	(g cm-3)
BSk22	$2.264$	11.20	0.967	$2.26\times {10}^{15}$
BSk24	$2.279$	11.08	0.973	$2.26\times {10}^{15}$
BSk25	$2.224$	11.05	0.987	$2.26\times {10}^{15}$
BSk26	$2.170$	10.20	1.123	$2.67\times {10}^{15}$

Table 17. Table 17: Radius, central density of baryons and central mass-energy density of a neutron star with M NS = 1.4 ​ M ⊙ subscript 𝑀 NS 1.4 subscript 𝑀 direct-product M_{\mathrm{NS}}=1.4\,M_{\odot} .

EoS	$R_{\mathrm{NS}}^{1.4}$ (km)	$n_c^{1.4}$ (fm-3)	${\rho }_c^{1.4}$ (g cm-3)
BSk22	13.04	0.385	$6.92\times {10}^{14}$
BSk24	12.57	0.408	$7.31\times {10}^{14}$
BSk25	12.37	0.416	$7.46\times {10}^{14}$
BSk26	11.77	0.506	$9.19\times {10}^{14}$

Table 18. Table 18: Causality limits

EoS	$n^{\mathrm{caus}}$ (fm-3)	${\rho }^{\mathrm{caus}}$ (g cm-3)
BSk22	1.095	$2.74\times {10}^{15}$
BSk24	1.088	$2.69\times {10}^{15}$
BSk25	1.378	$3.81\times {10}^{15}$
BSk26	0.982	$2.17\times {10}^{15}$

Table 19. Table 19: Direct Urca thresholds

EoS	$n_{\mathrm{DU}}$ (fm-3)	${\rho }_{\mathrm{DU}}$ (g cm-3)	$M_{\mathrm{DU}}/M_{\odot }$
BSk22	0.333	$5.88\times {10}^{14}$	$1.151$
BSk24	0.453	$8.25\times {10}^{14}$	$1.595$
BSk25	0.469	$8.56\times {10}^{14}$	$1.612$
BSk26	1.458	$4.19\times {10}^{15}$	$\left(2.115\right)$

Table 20. Table 20: Parameters of Eq. ( 75 ).

$i$	$p_i$
	BSk22	BSk24	BSk25	BSk26
1	7.02$\times {10}^8$	6.59$\times {10}^8$	6.411$\times {10}^8$	1.3995$\times {10}^8$
2	1.133$\times {10}^{11}$	9.49$\times {10}^{10}$	8.76$\times {10}^{10}$	1.066$\times {10}^9$
3	6.19$\times {10}^7$	6.95$\times {10}^7$	7.40$\times {10}^7$	7.472$\times {10}^8$
4	4.54$\times {10}^6$	5.63$\times {10}^6$	6.31$\times {10}^6$	3.599$\times {10}^7$
5	5.46$\times {10}^5$	6.51$\times {10}^5$	7.13$\times {10}^5$	2.906$\times {10}^6$
6	15.24	19.37	22.11	0.3051
7	0.0683	0.1028	0.1217	$-0.54068$
8	8.86	4.09	2.54	388.5
9	4611	6726	8317	776.5
10	48.07	29.57	25.63	15.435
11	2.697	2.6728	2.507	2.2483
12	81.7	19.51	7.92	0.3029
13	7.05	4.39	3.92	18.66
14	1.50	1.75	2.06	0.569

Table 21. Table 21: Parameters of Eq. ( 78 ).

$i$	$p_i$
	BSk22	BSk24	BSk25	BSk26
1	6.682	6.795	7.210	3.672
2	5.651	5.552	5.196	7.844
3	0.00459	0.00435	0.00328	0.00876
4	0.14359	0.13963	0.12516	0.22604
5	2.681	3.636	4.624	3.129
6	11.972	11.943	12.16	11.939
7	13.993	13.848	9.348	13.738
8	1.2904	1.3031	1.6624	1.3389
9	2.665	3.644	4.660	3.112
10	$-$27.787	$-$30.840	$-$28.232	$-$23.031
11	2.0140	2.2322	2.0638	1.6264
12	4.09	4.65	5.27	4.83
13	14.135	14.290	14.365	14.272
14	28.03	30.08	29.10	23.28
15	$-$1.921	$-$2.080	$-$2.130	$-$1.542
16	1.08	1.10	0.865	2.10
17	14.89	14.71	14.66	15.31
18	0.098	0.099	0.069	0.083
19	11.67	11.66	11.65	11.66
20	4.75	5.00	6.30	6.16
21	$-$0.037	$-$0.095	$-$0.172	$-$0.042
22	14.10	14.15	14.18	14.18
23	11.9	9.1	8.6	14.8

Table 22. Table 22: Parameters of Eq. ( 79 ).

EoS	$p_1$	$p_2$ [fm3]	$p_3$	$p_4$ [fm-3]
BSk22	14.5	220	$6\times {10}^4$	0.066
BSk24	$-68$	1770	$8\times {10}^5$	0.075
BSk25	$-445$	7380	$4\times {10}^6$	0.0806
BSk26	$-40.7$	1237	0	0

Table 23. Table 23: Parameters of Eq. ( 80 ).

EoS	$p_1$	$p_2$	$p_3$	$p_4$
BSk22	$6.08\times {10}^{-4}$	$3.42\times {10}^{-3}$	133.6	15.89
BSk24	$5.91\times {10}^{-4}$	$2.88\times {10}^{-3}$	196.3	12.14
BSk25	$5.82\times {10}^{-4}$	$2.67\times {10}^{-3}$	247.2	10.15
BSk26	$5.99\times {10}^{-4}$	$1.66\times {10}^{-3}$	137.5	8.91

Table 24. Table 24: Parameters of Eqs. ( 81 ) and ( 82 ).

EoS	$p_1$	$p_2$	$p_3$
BSk22	0.057	18.0	$-0.97$
BSk24	0.069	15.5	$-0.80$
BSk25	0.089	7.3	$-0.59$
BSk26	0.060	28.7	$-0.88$

Table 25. Table 25: Parameters of Eq. ( 85 ).

	BSk22	BSk24	BSk25	BSk26
$p_1$	173.2	195.2	203.7	216.4
$p_2$	204.8	233.8	251.1	323.0
$p_3$	57.3	96.8	132.2	$1.482\times {10}^4$
$p_4$	61.4	86.8	107.7	$1.241\times {10}^4$
$p_5$	1.86	2.41	3.37	1.99
$p_6$	0.080	0.116	0.154	0.062
$p_7$	18	29	56	10.4

Table 26. Table 26: Parameters of Eq. ( 91 ).

	BSk22	BSk24	BSk25	BSk26
$p_1$	$-0.0024$	0.0021	0.0164	$-0.0064$
$p_2$	0.456	0.581	0.410	0.690
$p_3$	12.349	9.874	4.471	0.6074
$p_4$	16.15	16.20	3.89	28.11
$p_5$	48.72	39.20	17.64	2.455
$p_6$	0.311	$-0.148$	$-0.517$	0.344
$p_7$	3.674	4	4	4.6

Table 27. Table 27: Parameters of Eq. ( 92 ).

EoS	$p_1$	$p_2$	$p_3$
BSk22	9.41	7.62	259
BSk24	4.027	18.31	106
BSk25	1.172	26.97	61.9
BSk26	8.042	$-6.86$	11.1

Table 28. Table 28: Parameters of Eq. ( 94 ).

EoS	$p_1$	$p_2$	$p_3$	$p_4$	$p_5$	$p_6$
BSk22	6824	2.4322	9.250	0.10698	2.248	0.0440
BSk24	2807	2.1785	8.055	0.18142	1.915	0.0716
BSk25	1434	1.9741	7.932	0.24071	2.009	0.1132
BSk26	3302	2.2416	6.268	0.14566	0.770	0.0400

Table 29. Table 29: Parameters of Eq. ( 95 ).

EoS	$p_1$	$p_2$	$p_3$	$p_4$	$p_5$	$p_6$	$p_7$
	for $C_p$:
BSk22	5.50	$3.04\times {10}^5$	4.22	2.75	0.5	0.462	0.486
BSk24	5.527	$5.29\times {10}^4$	4.00	1.574	1	0.631	0.720
BSk25	5.71	455	1.93	202	3	0.197	0.438
BSk26	5.555	$-3.90$	0.315	$-60.4$	2	0.294	0.444
	for $C_n$:
BSk22	5.73	$4.16\times {10}^5$	4.43	2.0	0.5	0.459	0.478
BSk24	5.706	$1.18\times {10}^4$	3.60	0.063	1	0.516	0.640
BSk25	3.77	45.1	0.855	$-139$	3	0.0626	0.031
BSk26	5.775	$-2.51$	0.250	$-48.4$	2	0.283	0.431
	for $a_p$:
BSk22	0.4376	3.23	0.877	$-731.55$	2.568	$-29.8$	3.153
BSk24	0.431	4.89	1.005	$-3729$	3.356	$-6.7\times {10}^8$	9.70
BSk25	0.434	7.86	1.159	$-2199$	3.355	$-18.5$	3.52
BSk26	0.434	4.09	0.970	$-2523$	3.244	$-1.07\times {10}^{12}$	12.7
	for $a_n$:
BSk22	0.656	5.50	0.671	$-1.39\times {10}^5$	4.643	$-0.463$	1.43
BSk24	0.631	5.27	0.729	$-2701$	3.377	$-335$	4.13
BSk25	0.636	12.6	0.972	$-1.70\times {10}^6$	6.34	$-42.6$	3.86
BSk26	0.633	1.94	0.511	$-136.6$	2.175	$-3.65$	2.442
	for $x_p$:
BSk22	0.1023	1.292	0.5114	343	2.605	1.12	2.59
BSk24	0.1035	1.944	0.5717	608	3.143	0.0225	1.26
BSk25	0.09975	2.747	0.6316	1748	3.906	$5.94\times {10}^{-4}$	0.192
BSk26	0.1025	1.383	0.526	6.69	1.238	45.8	4.25
	for $x_n$:
BSk22	0.1134	1.620	0.5163	132	2.41	229	4.89
BSk24	0.1137	2.562	0.5932	3333	4.63	160	4.66
BSk25	0.1116	3.132	0.6207	$-2.2$	2.78	$9.8\times {10}^4$	7.96
BSk26	0.1113	1.571	0.5196	2.49	0.966	71	4.47

Equations

$$S(n)=e_{\mathrm{NeuM}}(n)-e_{\mathrm{SNM}}(n),$$

$$\begin{gathered} e(n,\eta )=a_v+\left(J+\frac{1}{3}Lϵ\right){\eta }^2+\frac{1}{18}(K_v+{\eta }^2{K}_{\mathrm{sym}}){ϵ}^2 \\ +\cdots . \end{gathered}$$

$$e_{\mathrm{eq}}=M^{\prime }(A,Z)c^2/A+{\mathcal{E}}_{\mathrm{e}}(n_{\mathrm{e}})/\bar{n}+E_L(A,Z)/A-M_{\mathrm{n}}{c}^2.$$

$$\rho =\bar{n}\left(e_{\mathrm{eq}}/c^2+M_{\mathrm{n}}\right);$$

$$\begin{gathered} M^{\prime }(A,Z)c^2=M(A,Z)c^2+1.44381\times 10^{-5}{Z}^{2.39} \\ +1.55468\times 10^{-12}{Z}^{5.35} \end{gathered}$$

$$g=e_{\mathrm{eq}}+M_{\mathrm{n}}{c}^2+\frac{P}{\bar{n}}.$$

$${\mu }_{\mathrm{n}}=g,$$

$${\mu }_{\mathrm{p}}={\mu }_{\mathrm{n}}-{\mu }_{\mathrm{e}},$$

$${\mu }_{\mathrm{n}}=g=M_{\mathrm{n}}{c}^2$$

$${\mathcal{E}}(\bm{r})={\mathcal{E}}_{\mathrm{nuc}}(\bm{r})+{\mathcal{E}}_{\mathrm{ee,ep}}^{\mathrm{c}}(\bm{r})+{\mathcal{E}}_{\mathrm{e}}(\bm{r})\\ +\Big{[}\bar{n}_{\mathrm{n}}M_{\mathrm{n}}+\bar{n}_{\mathrm{p}}(M_{\mathrm{p}}+m_{\mathrm{e}})\Big{]}c^{2}\quad,$$

$${\mathcal{E}}_{\mathrm{nuc}}={\mathcal{E}}_{\mathrm{Sky}}+{\mathcal{E}}_{\mathrm{pp}}^{\mathrm{c}},$$

$$\mathcal{E}(\mathbf{r})={\mathcal{E}}_{\mathrm{Sky}}(\mathbf{r})+{\mathcal{E}}_{\mathrm{c}}(\mathbf{r})+{\mathcal{E}}_{\mathrm{e}}(\mathbf{r})+\bar{n}(M_{\mathrm{n}}{c}^2-Y_{\mathrm{p}}{Q}_{\mathrm{n},\beta }),$$

$$e=\frac{1}{A}{\int }_{V_{\mathrm{cell}}}\mathcal{E}(\mathbf{r}){\mathrm{d}}^{3}r=e_{\mathrm{Sky}}+e_{\mathrm{c}}+e_{\mathrm{e}}-Y_{\mathrm{p}}{Q}_{\mathrm{n},\beta },$$

$$n_q(r)=n_{\mathrm{B}q}+n_{\Lambda q}{f}_q(r),$$

$$f_{q}(r)=\frac{1}{1+\exp\left[\Big{(}\frac{C_{q}-R}{r-R}\Big{)}^{2}-1\right]\exp\Big{(}\frac{r-C_{q}}{a_{q}}\Big{)}}$$

$$A=V_{\mathrm{cell}}\bar{n}=\frac{4\pi }{3}{R}^3\bar{n},$$

$$Z=4\pi {\int }_0^R{r}^2{n}_{\mathrm{p}}(r)\mathrm{d}r$$

$$N=4\pi {\int }_0^R{r}^2{n}_{\mathrm{n}}(r)\mathrm{d}r.$$

$$Z_{\mathrm{cl}}=4\pi n_{\mathrm{\Lambda p}}{\int }_0^R{r}^2{f}_{\mathrm{p}}(r)\mathrm{d}r$$

$$N_{\mathrm{cl}}=4\pi n_{\mathrm{\Lambda n}}{\int }_0^R{r}^2{f}_{\mathrm{n}}(r)\mathrm{d}r.$$

$$e_{\mathrm{Sky}}^{\mathrm{ETF}}=\frac{4\pi }{A}{\int }_0^R{r}^2{\mathcal{E}}_{\mathrm{Sky}}^{\mathrm{ETF}}(r)\mathrm{d}r.$$

$$e_{\mathrm{Sky}}=e_{\mathrm{Sky}}^{\mathrm{ETF}}+\frac{1}{A}\left(E_{\mathrm{p}}^{\mathrm{sc},\mathrm{pair}}+E_{\mathrm{pair}}\right),$$

$${\mathcal{E}}^{\mathrm{c}}=\frac{e^2}{2}(n_{\mathrm{p}}(\mathbf{r})-n_{\mathrm{e}})\int \frac{n_{\mathrm{p}}({\mathbf{r}}^{\prime })-n_{\mathrm{e}}}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}{\mathrm{d}}^3{\mathbf{r}}^{\prime };$$

$$\Lambda =\frac{1}{V_{\mathrm{cell}}}\int {\mathrm{d}}^3\mathbf{r}{\left(\frac{n(\mathbf{r})}{\bar{n}}-1\right)}^2,$$

$$t(r=0)>U_{\mathrm{p}}(r=R)-U_{\mathrm{p}}(r=0),$$

$$t(r=0)=\frac{{\hslash }^2}{2M_{\mathrm{p}}^{\ast }(r=0)}{\left[3{\pi }^2{n}_{\mathrm{p}}(r=0)\right]}^{2/3},$$

$$P=P_{\mathrm{nuc}}+P_{\mathrm{e}},$$

$$C_1^n=-\frac{1}{4}{t}_0\left(\frac{1}{2}+x_0\right)-\frac{1}{24}{t}_3\left(\frac{1}{2}+x_3\right)n_{B0}^{\alpha };$$

$$P={\bar{n}}^2{\left(\frac{\partial e}{\partial \bar{n}}\right)}_{T,Y_{\mathrm{p}}},$$

$$Y_{\mathrm{e}}+Y_{\mu }=Y_{\mathrm{p}}.$$

$$e=\frac{1}{n}\left({\mathcal{E}}_{\mathrm{Sky}}+{\mathcal{E}}_{\mathrm{lept}}\right)-Y_{\mathrm{p}}{Q}_{\mathrm{n},\beta }+Y_{\mu }(m_{\mu }-m_{\mathrm{e}})c^2,$$

Full text

††pagerange: Unified equations of state for cold non-accreting neutron stars with Brussels-Montreal functionals. I. Role of symmetry energy–C.5

1 Introduction

Three distinct regions can be recognized in a neutron star below its thin atmosphere: a locally homogeneous core and two concentric shells characterized by different inhomogeneous phases (e.g., Haensel P., Potekhin A. Y. & Yakovlev, 2007; Chamel & Haensel, 2008). The outermost of the shells, the ‘outer crust’, consists of a lattice of nuclei and electrons that globally is electrically neutral. In the absence of accreted material, the surface of this crust is expected to be made of 56Fe. However, on moving towards the interior, more and more neutron rich nuclei appear (Rüster, Hempel & Schaffner-Bielich, 2006; Roca-Maza & Piekarewicz, 2008; Pearson, Goriely & Chamel, 2011; Kreim et al., 2013; Wolf et al., 2013; Chamel et al., 2015a; Utama, Piekarewicz & Prosper, 2016; Chamel et al., 2017), until at a mean baryon number density $\bar{n}$ of around $2.6\times 10^{-4}$ nucleons fm*-3* (a mass-energy density $\rho$ of around $4.3\times 10^{11}$ g cm${}^{-3})$ neutron drip sets in (see, e.g., Chamel et al. 2015b for a recent discussion). This marks the transition to the ‘inner crust’, an inhomogeneous assembly of neutron-proton clusters and unbound neutrons, neutralized by electrons (proton drip can also set in at higher densities). By the point where $\bar{n}$ has risen to about half the equilibrium density $n_{0}$ of infinite (homogeneous) nuclear matter (INM) that is charge-symmetric (we denote by SNM this special case of INM), the inhomogeneities have been smoothed out and we enter the core of the star. The homogeneous medium of which the core is comprised is known as ‘neutron-star matter’ (NM), and is made up primarily of neutrons, with an admixture of protons neutralized by electrons and, at densities above $n\simeq 0.12$ fm-3*, by muons. At higher densities, other particles such as hyperons might appear (Haensel et al., 2007; Weber et al., 2007), but we do not consider this possibility here.

In this paper we continue our project of developing a unified treatment of neutron stars within the framework of the picture of ‘cold catalysed matter’ (CCM), by which we mean that thermal, nuclear and beta equilibrium prevail at a temperature $T$ low enough that thermal effects are negligible for the composition and pressure. The equilibrium conditions can reasonably be expected to be valid in any neutron star that is not accreting from a neighbour, but will otherwise fail because of the relative slowness with which accreted matter acquires nuclear equilibrium. For densities greater than around $10^{-10}$ nucleons fm*-3* (106 g cm${}^{-3})$ all the atoms of the outer crust are completely ionized and the electrons form a degenerate gas. Since nuclear degeneracy likewise holds everywhere it follows that the CCM picture is valid throughout the star except in a thin layer of $\rho\lesssim 10^{6}$ g cm*-3*, where the atomic ionization and electron degeneracy can be incomplete. The unified treatment of all three regions of neutron stars that we present in this paper therefore does not include these outermost parts, which do not in any case involve any new nuclear physics, and have been extensively discussed by Haensel et al. (2007).

Our calculations of the degenerate equation of state (EoS) and the composition of the three regions are microscopic, and the unifying feature to which we have alluded lies in the fact that in each region we use the energy-density functional theory with the same functional. For this we shall take one or other of the functionals that we have developed in the last few years not only for the study of neutron-star structure but also for the general purpose of providing a unified treatment of a wide variety of phenomena associated with the birth and death of neutron stars, such as supernova-core collapse and neutron-star mergers, along with the r-process of nucleosynthesis (both in the neutrino-driven wind and via the decompression of N*M). These functionals are based on generalized Skyrme-type forces and contact pairing forces, the formalism for which is presented in the Appendix of Chamel, Goriely & Pearson (2009). The first set of such functionals (Goriely, Chamel & Pearson, 2010) that we devised, labelled BSk19, BSk20 and BSk21, has already been applied to a unified treatment of neutron-star structure by Pearson, Goriely & Chamel (2011); Chamel et al. (2011); Pearson et al. (2012); Potekhin et al. (2013); Fantina et al. (2013).

The parameters of these functionals were determined primarily by fitting to essentially all the data of the 2003 Atomic Mass Evaluation, AME 2003 (Audi, Wapstra & Thibault, 2003), with the nuclear masses being calculated using the Hartree-Fock-Bogoliubov (HFB) method. To do this it was necessary to add to the HFB energy phenomenological Wigner terms and correction terms for the spurious collective energy. Note that our HFB code takes account of axial deformations (Samyn et al., 2002). For each functional complete HFB nuclear-mass tables, running from one drip line to the other, and labelled HFB-19, HFB-20, and HFB-21, respectively, were constructed (Goriely et al., 2010). The supplementary terms that we had to add to the HFB energy to calculate nuclear masses do not enter into our calculation of the inner crust and core, the functional alone being sufficient, while for the outer crust the entire nuclear dependence is subsumed into the nuclear mass. Thus this latter region is calculated directly from the appropriate mass table (or from experimental mass data, where available).

Since all the astrophysical applications that we envisage involve a long-range extrapolation from experimentally accessible environments to highly neutron-rich environments our functionals incorporate as much well established theoretical knowledge of neutron-rich systems as possible. The most significant such constraint that we impose is to require that our nuclear-mass fits be consistent, up to the densities prevailing in neutron-star cores, with the EoS of homogeneous pure neutron matter (NeuM), as calculated by many-body theory from realistic two- and three-nucleon forces. Several such EoSs have have been published, differing considerably in their stiffness at supersaturation densities. Functional BSk19 was fitted to the EoS of Friedman & Pandharipande (1981) (FP), BSk20 to the somewhat stiffer EoS of Akmal, Pandharipande & Ravenhall (1998), which they label as ‘A18 + $\delta\,v$ + UIX∗’ and which we refer to as APR, while functional BSk21 was fitted to a still stiffer EoS of NeuM, the one labelled ‘V18’ in Li & Schulze (2008), and which we refer to as LS2.

When we came to apply our three EoSs, BSk19-21, to neutron-star structure we found that the choice of the EoS of NeuM to which we fit our functionals has a significant impact on the maximum possible mass of neutron stars. In particular, we found that if we constrained to the EoS of FP then it was impossible to support the heaviest neutron star that had been observed (Chamel et al., 2011). On the other hand, there was no such problem with the APR and LS2 constraints.

The present paper stems in part from the realization that fitting NeuM and nuclear masses (along with some other properties of real nuclei such as charge radii) does not exhaust completely the degrees of freedom allowed by our generalized Skyrme functionals. In particular, it does not tie down the symmetry energy of our functionals, defined here as the difference between the energy per nucleon in NeuM and the energy per nucleon in SNM,

[TABLE]

essentially because the fit to nuclear masses does not determine $e_{\mathrm{SNM}}$ uniquely. We are thus left with a certain flexibility in the choice of the symmetry coefficient $J$, defined as follows in an expansion of the energy per nucleon of INM of density $n=n_{0}(1+\epsilon)$ and charge asymmetry $\eta=(n_{\mathrm{n}}-n_{\mathrm{p}})/n$ (with $n_{\mathrm{n}}$ and $n_{\mathrm{p}}$ the neutron and proton number densities respectively) about the equilibrium density $n=n_{0}$ and $\eta=0$:

[TABLE]

(see Eqs. (5) – (7) of Goriely, Chamel & Pearson 2013).

Now the functionals BSk19-21 had all been constrained to $J$ = 30 MeV. Thus to investigate the degree of freedom remaining on $J$ we constructed a new series of functionals constrained to different values of $J$, but fitted to the same EoS of NeuM, for which we took LS2 (Goriely et al., 2013). In the meantime the 2012 Atomic Mass Evaluation, AME 2012 (Audi et al., 2012), had become available, and we fitted the new functionals to the 2353 measured masses of nuclei having $N$ and $Z\geq$ 8 given there. In this way we generated four new functionals, BSk22, BSk23, BSk24 and BSk25 fitted to $J$ = 32, 31, 30 and 29 MeV, respectively; we were unable to find acceptable mass fits outside this range.

In addition to the new functionals and the new mass data that have become available, another recent development has been a significant improvement in our code for computing the inner crust, with the inclusion of pairing and the consideration of proton drip for the first time in any of our work (see Section 2.2). This makes it worthwhile to revisit our earlier work with functionals BSk19-21 relating to the NeuM constraint, and so we constructed a fifth new functional, BSk26, constrained like BSk24 to $J$ = 30 MeV, but fitted to the softer APR EoS of NeuM. On the other hand, we do not consider functional BSk23 in the present paper, since it was fitted to $J=31$ MeV, and thus is intermediate between BSk22 and BSk24.

Then to assess the impact of the symmetry energy we compare functionals BSk22, BSk24 and BSk25, while to assess the impact of the NeuM constraint we compare BSk24 and BSk26. In Table 1 we list for these four functionals the values of $J$ along with the corresponding values of the symmetry-slope coefficient $L$, as well as the incompressibility $K_{v}$ and the symmetry-incompressibility $K_{\mathrm{sym}}$. A comparison of BSk22, BSk24 and BSk25 shows the correlation between $J$ and $L$ that has been noted many times over the last forty years: see, for example, Farine, Pearson & Rouben (1978); Lattimer (2012); Lattimer & Lim (2013); Baldo & Burgio (2016). This correlation arises primarily from the fits to nuclear masses, an increase in $J$ being offset by the $L$ term in the coefficient of $\eta^{2}$ in Eq. (2) for subnuclear densities ($\epsilon<0$), which are dominant for finite nuclei. However, $J$ cannot determine $L$ uniquely since the EoS must play a role: while BSk24 and BSk26 have the same $J$, the former has a larger value of $L$, as might be expected for a functional fitted to an EoS of NeuM that is stiffer at higher densities, $\epsilon>0$. These questions are discussed in greater detail in Section IIIC1 of Goriely et al. (2013).

The quality of the mass fits of the four functionals is shown in Table 2. Note particularly that we calculate the deviations with respect to the data of the 2016 AME (Wang et al., 2017) and not the data to which our functionals were fitted, those of the 2012 AME (Audi et al., 2012). For this reason the deviations we show in Table 2 are slightly different from those of Table IV of Goriely et al. (2013). We see that BSk24 gives the best fit to the masses of the neutron-rich nuclei, which are of greater relevance for our purpose. It also is apparent that BSk22 gives a significantly worse fit than do the other functionals.

Fig. 1 shows the EoSs of completely degenerate NeuM for our four functionals; it will be seen that while the fits of BSk24 and BSk25 to the LS2 EoS are excellent, that of BSk22 is less so at high densities, being somewhat stiffer than LS2. The fit of BSk26 to the APR EoS is very good.

In Figs. 2 and 3 we plot for our four functionals the symmetry energy defined in Eq. (1). Comparing BSk24 and BSk26, we see how at high densities the symmetry energy of the latter increases much less steeply, as is appropriate, given the much softer EoS of NeuM to which it has been fitted. On the other hand, at subnuclear densities the two functionals have virtually identical symmetry energy (this is particularly apparent in fig. 2b of Goriely et al. 2013), showing that the constraining EoS of NeuM is irrelevant at these densities, and that it is the symmetry coefficient $J$ that is the determining factor. Comparing now BSk22, BSk24 and BSk25 to investigate the role of $J$ we see that while there is a correlation between $S(n)$ and $J$ at super-nuclear densities (at least up to around 4$n_{0}$), these quantities are anticorrelated at lower densities ($n<n_{0}$). Moreover, while the functionals BSk24 and BSk25 behave very similarly in NeuM their symmetry energies start to diverge at higher densities. This can only be a result of differing $e_{\mathrm{SNM}}$, a difference that is easily understood given that the same mass data are being fitted with different values of $J$. It will be seen in the course of this paper how far these symmetry properties of INM are reflected in the EoSs for the different functionals in the various regions of the neutron star.

While the LS2 EoS of NeuM to which we fit functionals BSk22-25 can be regarded as typically hard, and the APR EoS to which we fit BSk26 as typically soft, the question arises as to what extent the ‘real’ EoS of NeuM could lie outside these limits. Certainly, an EoS much softer than APR, such as FP, can be ruled out as long as we neglect the possibility of neutron stars with exotic cores (Chamel et al., 2013). Concerning the possibility of an EoS that is stiffer than LS2 we recall the case of the EoS labelled BOB by Li & Schulze (2008) and LS3 by Goriely et al. (2013) and Goriely, Chamel & Pearson (2016), but we ruled it out of our considerations as probably being too soft at low densities to satisfy the constraints imposed by chiral effective field theory (Tews et al., 2013): see the discussion in Section IVA of Goriely et al. (2016). At the same time we recall from the discussion in Section IIIA of Goriely et al. (2013) that the functionals of this paper are consistent with the quantum Monte Carlo calculations of Gandolfi, Carlson & Reddy (2012). In the time since our functionals (Goriely et al., 2013) were constructed in 2013, a number of new calculations of NeuM using chiral effective field theory (Roggero, Mukherjee & Pederiva, 2014; Rrapaj, Roggero & Holt, 2016; Drischler et al., 2016; Tews et al., 2016) have appeared. However, these calculations were restricted to densities below 0.3 fm*-3*; two of them at least (Roggero et al., 2014; Rrapaj et al., 2016) appear to favour LS2 over APR in this low-density regime, where LS2 is softer than APR. At higher densities we are not aware of any later calculations indicating that the real EoS of NeuM is stiffer than LS2. This is corroborated by recent gravitational-wave observations (see Section 4.1).

Concerning the value of the symmetry coefficient $J$, we have been unable to fit, with the same accuracy, our generalized Skyrme density functionals to the mass data with values of $J$ lying outside the range 32 MeV $\geq J\geq$29 MeV; indeed, $J$ = 32 MeV might be too high. On the other hand, many relativistic mean-field (RMF) functionals are characterized with significantly higher values of $J$. However, the high values quoted do not come from mass fits at all, but rather from measurements on giant dipole resonances, pygmy resonances or neutron-skin thicknesses. The point is that so far it has not been possible to fit masses in the RMF framework with a precision at all comparable to what has been achieved with Skyrme functionals, and as a result the quality of the fits becomes relatively insensitive to $J$. It should nevertheless be noted that the best RMF mass fits, those of Ref. Peña-Arteaga, Goriely & Chamel (2016), for which rms deviations of about 1.2 MeV were found for 2353 masses, have $J$ = 30 or 31 MeV.

Our methods of calculation for the three different regions of neutron stars are presented in Section 2, where we summarize much of our work from earlier papers and describe our current refinements, particularly with regards to proton drip. Section 3 presents our numerical results for the three different regions of neutron stars in both graphical and tabular form (extensive tables relating to the inner crust will be found in the supplementary material). In Section 4 we apply our EoSs to some gross properties of neutron stars, namely the mass-radius relation, the maximum mass and the direct Urca process. Our conclusions are summarized in Section 5. In Appendix A we derive a result used to calculate the chemical potentials of neutrons and protons, while our treatment of leptonic gases appears in Appendix B. Finally, Appendix C presents what is in many ways the central feature of this paper: in order to facilitate the application of our results to the modelling of neutron-star structure they are fitted by analytic parametrizations.

2 Methods of calculation

2.1 Outer crust

As discussed in the Introduction, we calculate the outer crust only for densities $\rho\gtrsim 10^{6}$ g cm*-3* ($\bar{n}\gtrsim 10^{-9}$ nucleons fm*-3*). We also suppose that the temperature is not only low enough for complete degeneracy to hold, but is also lower than the crystallization temperature, so that the nuclei of the outer crust can be supposed to be arranged in a regular crystal lattice. It is easy to check (see, e.g., Chapter 2 of Haensel et al. 2007) that for 56Fe this is true at $\rho\gtrsim 3\times 10^{6}\,(T/10^{8}\mbox{~{}K})^{3}$ g cm*-3*. For simplicity, we suppose the crystal lattice to be made of only one type of ion ($A,Z$) with proton number $Z$ and mass number $A$, for which the most stable crystal structure is the body-centred cubic (bcc) lattice (see Chamel & Fantina 2016a for a recent discussion on multinary compounds in the outer crust of neutron stars).

We calculate the EoS and the composition of the outer crust mainly as described in detail by Pearson et al. (2011). To recall the essentials, in a region of mean density $\bar{n}$ let the equilibrium nucleus have mass number $A$ and atomic number $Z$. Then the energy per nucleon in that region is

[TABLE]

Here $M^{\prime}(A,Z)$ is the atomic mass for the element ($A,Z$) with the binding energy of the atomic electrons subtracted, $\mathcal{E}_{\mathrm{e}}(n_{\mathrm{e}})$ is the energy density of the electrons (without the electron mass), $E_{L}(A,Z)$ is the lattice energy per nucleus, and the constant term in the neutron mass $M_{\mathrm{n}}$ is subtracted out for convenience, as in our treatment of the inner crust and the core. Note particularly that this missing neutron mass has to be restored if we want the corresponding total mass-energy density, thus

[TABLE]

and likewise in Eq. (6) below.

For the first of the quantities in Eq. (3) we have, in units of MeV,

[TABLE]

(see Eq. (A4) of Lunney, Pearson & Thibault 2003), in which $M(A,Z)$ is just the usual atomic mass. We use for this latter quantity the HFB-22, HFB-24 HFB-25 or HFB-26 tables available on the bruslib database (Xu et al.., 2013), except when experimental values are available, for which we use the 2016 Atomic Mass Evaluation (Wang et al., 2017), supplemented by the very recent measurements of copper isotopes (Welker et al., 2017).

For the energy density $\mathcal{E}_{\mathrm{e}}(n_{\mathrm{e}})$ of the electrons ($n_{\mathrm{e}}=Z\bar{n}/A$), we now, unlike Pearson et al. (2011), use complete expressions for the electron exchange and screening (also sometimes referred to as ‘polarization’) corrections to the electron energy density, as discussed in Appendix B. The electron-correlation energy is included as in Pearson et al. (2011). The lattice energy $E_{L}(A,Z)$ is calculated as in Pearson et al. (2011), with quantum zero-point and finite nuclear-size corrections.

The EoS and the composition are now determined by minimizing, at fixed pressure $P$, the Gibbs free energy per nucleon

[TABLE]

The pressure is determined from the energy $e_{\mathrm{eq}}$, as described by Pearson et al. (2011) and in Appendix B.

We calculate the neutron chemical potential using the identity

[TABLE]

valid only because of the beta equilibrium that holds throughout the neutron star (see Baym, Pethick & Sutherland 1971a and Appendix A). Because of this beta equilibrium we can then calculate the proton chemical potential through

[TABLE]

where $\mu_{\mathrm{e}}$ is the electron chemical potential (see Appendix B).

The minimization at constant pressure of the Gibbs free energy per nucleon $g$ is repeated with increasing value of the pressure until neutron drip sets in, the condition for which is

[TABLE]

(this is equivalent to Eq. (2.7.3) of Shapiro & Teukolsky 1983).

2.2 Inner crust

Since the theoretical nuclear masses used in Section 2.1 for the outer crust were derived from our functionals by application of the HFB method, it might seem appropriate to use this same method for the inner-crust calculations also. As in the pioneer HF calculations of Negele & Vautherin (1973), such calculations are generally performed within the framework of the spherical Wigner-Seitz (WS) approximation (Montani, May & Müther, 2004; Baldo, Saperstein & Tolokonnikov, 2007; Grill, Margueron & Sandulescu, 2011). But an inevitable consequence of the WS approximation is to introduce shell effects in the spectrum of unbound neutron states, which dominate the properties of the inner crust. Such shell effects are to a large extent spurious, since in reality the unbound neutron states form a quasi-continuum, as shown by band-structure calculations (Chamel, 2005, 2006, 2012). This difficulty is analysed in detail by Chamel et al. (2007); Margueron, van Giai & Sandulescu (2007), the latter reference showing that the error thereby introduced in the energy per nucleon cannot easily be reduced below 50 keV, which is incompatible with a reliable calculation of the composition of the inner crust. Pastore et al. (2017) have recently proposed a new method to reduce the errors to a few keV by considering a supercell with a fixed radius of 80 fm. But they also stressed that the WS treatment becomes unreliable at densities above 0.02 fm*-3*. In the last few years 3D calculations (mostly at finite temperatures and fixed proton fractions) have been carried out (Magierski & Heenen, 2002; Gögelein & Müther, 2007; Newton & Stone, 2009; Pais & Stone, 2012; Schuetrumpf et al., 2013, 2014; Pais, Newton & Stone, 2014; Schuetrumpf et al., 2015; Sagert et al., 2016). However, the use of a cubic box with periodic boundary conditions can still lead to spurious shell effects (see, for example, Section IIC2 in Newton & Stone 2009) that can only be removed by choosing a large enough box (Sébille, Figerou & de la Mota, 2009, 2011) or imposing Bloch boundary conditions (Chamel, 2012; Schuetrumpf & Nazarewicz, 2015); moreover such calculations require computer times that are quite impractical for the applications that we are undertaking here.

In view of these problems it is not surprising that a more popular approach to the calculation of the inner crust has been to use the much simpler compressible liquid-drop model; a typical such calculation is that of Douchin & Haensel (2001). Within each WS cell this method makes a clear separation of INM into two distinct homogeneous phases, the densities of which are free parameters of the model. The bulk properties of the two phases are calculated microscopically using the adopted functional, as are the surface properties (preferably including curvature corrections) of the interface between them. A more accurate treatment of spatial inhomogeneities is to employ semi-classical methods such as the Thomas-Fermi (TF) approximation, as for instance in Oyamatsu & Iida (1997); Okamoto et al. (2013); Sharma et al. (2015), or its extension to second order in $\hbar^{2}$ by Martin & Urban (2015) (also private communication from N. Martin). However, none of these calculations makes any shell correction.

The shortcomings of all these different methods are avoided in our application of the ETFSI (fourth-order extended Thomas-Fermi plus Strutinsky integral) method to the calculation of the EoS of the inner crust (Dutta, Onsi & Pearson, 2004; Onsi et al., 2008; Pearson et al., 2012). In these papers we extend a method that was originally developed as a fast approximation to a Skyrme-HF-BCS treatment of finite nuclei (Dutta et al., 1986): it was based on a full ETF treatment of the kinetic-energy and spin current densities, with neutron and proton shell corrections added perturbatively. When compared to exact HF-BCS calculations performed with the same Skyrme and pairing forces it was found to overbind by 170 keV/nucleon at most (for 16O), and only a few tens of keV/nucleon for heavier elements. For comparison, zero-order TF calculations with parametrized distributions (whose parameters are fully determined by $Z$, $A$, and nuclear-matter properties) lead to errors of order 0.5-1 MeV/nucleon (Papakonstantinou et al., 2013). The extension to second order reduces the errors below 300 keV/nucleon (Aymard, Gulminelli & Margueron, 2014).

However, for most purposes, what counts is not the absolute binding energy but the differences in energy between different configurations, e.g., different nuclei, or WS cells of different composition. Viewed in this light the ETFSI method emerges as a much better approximation to the Skyrme-HF-BCS method for finite nuclei. Referring to the 69 spherical nuclei listed in Table 2 of Dutta et al. (1986), we took the differences of 48 pairs (in most cases giving thereby nucleon separation or beta-decay energies), and found that the greatest discrepancy between the ETFSI and HF results was 5 keV per nucleon, and generally was much less.

When the ETFSI method is applied, as here, to the EoS problem, it allows, like the TF method, for a continuous variation of the density of nuclear matter within each spherical WS cell, without any artificial separation into two distinct phases. However it is expected to provide a much better description of nuclear clusters than does the TF method, since the semi-classical expressions for the kinetic-energy and spin current densities include density-gradient terms up to the fourth order. Moreover, our application of the ETFSI method to the EoS problem (Dutta et al., 2004; Onsi et al., 2008; Pearson et al., 2012) includes proton shell corrections. On the other hand, we do not calculate the neutron shell corrections since they are known to be much smaller than the proton shell corrections when there are unbound neutrons (Oyamatsu & Yamada, 1994; Chamel, 2006; Chamel et al., 2007); the error thereby introduced is negligible compared to the spuriously large neutron shell corrections that arise in current implementations of the HF equations, as noted above. Also, while our earlier ETFSI calculations of the EoS (Dutta et al., 2004; Onsi et al., 2008; Pearson et al., 2012) did not include pairing, Pearson et al. (2015) showed how to include, for zero temperature, a BCS treatment of proton pairing, and we follow that procedure here. It is to be noted that the errors of a few keV per nucleon incurred by adopting the BCS approximation (instead of solving the full HFB equations) lie within the errors of the ETFSI approach (Pastore, Shelley & Diget, 2016).

To recall the main features of our method, we write the total energy density at any point in the inner crust as

[TABLE]

where $M_{\mathrm{n}},M_{\mathrm{p}}$ and $m_{\mathrm{e}}$ denote the masses of the neutron, proton and electron, respectively. The first term here can be broken up into a term corresponding to the generalized Skyrme force (in which we include the nucleonic kinetic energy) and a proton-proton Coulomb term,

[TABLE]

while the second term represents the total direct part of the Coulomb interactions of the electrons, i.e., the $e-e$ and $e-p$ interactions. The third term, which is purely electronic, contains both kinetic-energy and Coulomb-exchange terms: see Appendix B. Bracketing together the two Coulomb terms, ${\mathcal{E}}_{\mathrm{pp}}^{\mathrm{c}}+{\mathcal{E}}_{\mathrm{ee,ep}}^{\mathrm{c}}={\mathcal{E}}_{\mathrm{c}}$, we can rewrite Eq. (10) as

[TABLE]

where $Y_{\mathrm{p}}=Z/A$ with $Z$ and $A$ the total number of protons and nucleons in the WS cell respectively, and $Q_{\mathrm{n,\beta}}$ is the $\beta$-decay energy of the neutron (= 0.782 MeV). Then for the total energy per nucleon we have, after integrating over the WS cell of volume $V_{\mathrm{cell}}$,

[TABLE]

in which we have dropped for convenience the constant term $M_{\mathrm{n}}c^{2}$, as we do for the outer crust and the core.

The following salient points are to be noted for the calculation of the first two terms. With spherical symmetry being assumed the neutron and proton density distributions within the WS cells are parametrized as the sum of a constant ‘background’ term and a ‘cluster’ term according to

[TABLE]

in which, with $q=n$ or $p$,

[TABLE]

(if $n_{\Lambda q}$ is negative the cluster becomes a bubble). This ‘damped’ modification of the usual simple Fermi profile is introduced in order to satisfy the requirement of the usual implementation of the ETF method that the first three derivatives of the density vanish at the cell surface (see section II of Pearson et al. 2012, where other relevant details will be found). At the same time, a smooth matching of the nucleonic distributions between adjacent cells is ensured.

With $Z$ protons and $N$ neutrons in the cell, the cell radius $R$ is determined by the mean density $\bar{n}$ through

[TABLE]

where $A=Z+N$. The eight independent cell parameters that are left, four for each charge type, will then be constrained by

[TABLE]

It is convenient now to define a proton-cluster number and a neutron-cluster number through

[TABLE]

Using the density distribution (14), one can calculate (Chamel et al., 2009) an ETF energy density, ${\mathcal{E}}_{\mathrm{Sky}}^{\mathrm{ETF}}(r)$ for the Skyrme force, whence for the ETF energy per nucleon we have

[TABLE]

Since we now include pairing as well as shell corrections we have

[TABLE]

in which $E^{\mathrm{sc,pair}}_{\mathrm{p}}$ is the Strutinsky-integral shell correction, as modified by pairing, and $E_{\mathrm{pair}}$ is the BCS energy (see Eqs. (5) and (6), respectively, of Pearson et al. 2015).

For the direct part of the Coulomb term $e_{\mathrm{c}}$ in Eq. (13) we have

[TABLE]

there is no protonic Coulomb exchange term for the functionals of this paper (see Sect II of Goriely et al. 2013). Both the direct and exchange parts of the electronic term $e_{\mathrm{c}}$ are discussed in Appendix B.

For a given mean density $\bar{n}$ nucleons per unit volume the total energy $e$ per nucleon has to be minimized with respect to eight parameters: $N$, $Z$ and six of the eight geometric parameters defined in Eqs. (14) and (15). In making this minimization care has to be taken that in allowing for bubble configurations we do not encounter negative densities at any point for either charge type of nucleon: see Onsi, Przysiezniak & Pearson (1997).

Insofar as shell corrections are included for protons, $Z$ has to be discretized to integral values in the minimization, but $N$ is treated as a continuous variable. Even though the total number of neutrons in the crustal layer is integral, the notion of a fractional number of neutrons per WS cell corresponds to the physical reality, since the neutrons are delocalized. Then with $Z$ being held fixed at different integral values automatic minimization is performed with respect to $N$ and six geometric parameters.

Since our 2012 paper (Pearson et al., 2012) we have made significant improvements to the code that calculates the inner crust. It will be recalled that at higher densities close to the point of transition to a homogeneous medium our search routine broke down and was unable to find any solution, when minimizing with respect to all seven parameters. We have now rectified this problem with our code, and can in fact find solutions at densities known to correspond to homogeneous nuclear matter. To quantify this capability it is convenient to use as a global measure of the departure from homogeneity the ‘inhomogeneity factor’

[TABLE]

where the integration goes over the cell. At the interface between the inner and outer crusts we find values of around 400, but at densities corresponding to homogeneous nuclear matter we find values of $\Lambda$ of the order of 10*-7*, or even less, values that are appropriately small, suggesting that our code is working correctly even at densities higher than for which it was intended.

Proton drip. Now that our inner-crust code works properly right into the homogeneous region we can address the problem of proton drip, which we neglected in our earlier work, but which becomes significant at higher densities. At all values of $\bar{n}$ in the inner crust, there is a non-vanishing proton density at the surface of the WS cell, but it rises sharply when $\bar{n}$ is sufficiently high for the motion of the highest-energy proton to become infinite in the classicle single-particle approximation, that is when

[TABLE]

where $t(r=0)$ is the maximum proton kinetic energy at the centre of the cell and $U_{\mathrm{p}}(r)$ is the value of the proton single particle effective potential at the point $r$. For the former quantity we can safely use the TF expression at the centre of the cell, thus

[TABLE]

where $M^{*}_{\mathrm{p}}(r=0)$ is the effective proton mass at the centre of the cell (see the Appendix of Chamel et al. 2009). Once this proton drip has begun the unbound proton single particle states will form a quasi-continuum, and proton shell effects will largely vanish, exactly as do neutron shell effects at all densities in the inner crust, i.e., beyond the neutron-drip point. For values of $\bar{n}$ below the proton drip point, as defined by the condition (23), the shell correction $E^{\mathrm{sc,pair}}_{\mathrm{p}}$ of Eq. (20) is kept in its entirety, along with the BCS energy $E_{\mathrm{pair}}$; otherwise both terms are dropped completely. This sharp cutoff of the shell correction probably exaggerates the actual situation, but it is not devoid of all physical sense, and we see no simple way of smoothing it that would not introduce an arbitrary element. Simultaneously dropping $E_{\mathrm{pair}}$ is more contentious, and is done for simplicity; nevertheless, since pairing contributes to the shell correction $E^{\mathrm{sc,pair}}_{\mathrm{p}}$ it would not be entirely logical to drop the latter and not $E_{\mathrm{pair}}$. In any case, $E_{\mathrm{pair}}$ will be quite small, since the gas of free protons will be much more dilute than the system of bound protons, with the result that the pairing field is much weaker.

Once protons begin to unbind and we drop shell corrections, $Z$ is no longer restricted to integral values and in principle should be treated, like $N$, as a continuous variable. However, since it is more reliable to minimize with respect to seven than eight variables, we discretize $Z$ in intervals of 0.1.

Pressure. The pressure in the inner crust at any mean density $\bar{n}$ is calculated as described in Appendix B of Pearson et al. (2012). Briefly, we can write the pressure as a sum of nucleonic and electronic terms,

[TABLE]

in which both the direct and exchange contributions to the electronic term are found in Appendix B of the present paper, while the nuclear term, which contains contributions from both neutrons and protons, is given by the simple expression (B28) of Pearson et al. (2012). We point out here that there is a misprint in Eq. (B30b) of this latter reference, which should read

[TABLE]

the calculations of Pearson et al. (2012) were made with this correct expression. Note that the second term of Eq. (B27) of Pearson et al. (2012) vanishes here because we drop the Coulomb-exchange terms for protons.

This procedure is more reliable and much faster than applying numerical differentiation to the identity

[TABLE]

in which it is understood that electrical neutrality is maintained.

Chemical potentials. The nucleonic chemical potentials $\mu_{\mathrm{n}}$ and $\mu_{\mathrm{p}}$ are calculated as in the outer crust, i.e., by using Eqs. (7) and (8). Note that these equations can only be applied when beta equilibrium holds.

Non-spherical configurations. It should be stressed that in both this and our previous papers we restrict ourselves to WS cells that are spherically symmetric. However, many, but not all, calculations, beginning with the work of Ravenhall, Pethick & Wilson (1983) and Hashimoto, Seki & Yamada (1984) (see also Section 3.3 of Chamel & Haensel 2008 and Section 3.4.2 of Haensel et al. 2007) suggest that as the interface with the core is approached there might be a slight energetic preference for non-spherical shapes, referred to collectively as ‘pasta’. Of especial relevance to our own work is that of Martin & Urban (2015), in which ETF calculations without shell corrections were performed over a wide range of densities with a number of Skyrme functionals, including BSk22 and BSk24. These authors report a transition to non-spherical shapes at densities slightly below 0.06 fm*-3* in the inner crust, but their conclusions are sensitive to the approximations made since the energy differences involved are very small. In particular, the parametrization that they adopted for the nucleon distributions does not satisfy the necessary boundary condition that the first three derivatives vanish at the cell surface (see the discussion after Eq. (15), and also in Section 3.5). Moreover, their ETF calculations are limited to second order (Martin & Urban, 2015), and it is known from work on fission barriers that the fourth-order terms can make departure from sphericity energetically expensive: see section 4.3 of Brack, Guet & Håkansson (1985). We have confirmed that the fourth-order terms become increasingly positive as the deformation increases, using a code based on Tondeur et al. (1987) that aims to extend the present work to deformed WS cells but which has not yet been fully developed.

Adding further doubt to the existence of pasta, Viñas et al. (2017) have recently performed self-consistent TF calculations in which the Euler-Lagrange equations were solved without any parametrization of the density distributions, and have found no pasta for several Skyrme functionals. In particular, their calculations predict that spherical clusters remain stable throughout the inner crust for SLy4, in disagreement with the results obtained by Martin & Urban (2015). However, our code for deformed WS cells shows that the second-order ETF terms, unlike the fourth-order terms, can reduce the deformation energy. Thus it is not clear that the conclusion of Viñas et al. (2017) would survive a full ETF calculation.

2.3 Core

Because of translational invariance a considerable simplification occurs in the N*M that defines the core, even though muons may now be present. To take account of the latter we denote by $Y_{\mu}$ the number of muons per nucleon, where, because of overall charge neutrality,

[TABLE]

Then in place of Eq. (13) we now have for the total energy per nucleon

[TABLE]

where the leptonic term ${\mathcal{E}}_{\mathrm{lept}}={\mathcal{E}}_{\mathrm{e}}+{\mathcal{E}}_{\mu}$ now includes the muonic contribution ${\mathcal{E}}_{\mu}$, and we have again dropped the neutron mass. Note that because of exact neutrality at all points there is no direct Coulomb term ${\mathcal{E}}_{\mathrm{c}}$, while the exchange terms are included in ${\mathcal{E}}_{\mathrm{lept}}$ (see Appendix B), there being no protonic Coulomb exchange term for the functionals of this paper.

A great simplification in the $\mathcal{E}_{\mathrm{Sky}}$ term of Eq. (29) arises, since in the zero-temperature approximation considered here it is given by an analytic expression. Assuming the system to be unpolarized, i.e., if time-reversal invariance holds, we have from Eq. (A13) of Chamel et al. (2009)

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

The leptonic term in Eq. (29) is ${\mathcal{E}}_{\mathrm{lept}}={\mathcal{E}}^{\mathrm{kin}}_{\mathrm{e}}+{\mathcal{E}}^{\mathrm{kin}}_{\mu}+{\mathcal{E}}^{\mathrm{ex}}_{\mathrm{e}}+{\mathcal{E}}^{\mathrm{ex}}_{\mu}$, in which all quantities are calculated as described in Appendix B.

For the total pressure we can now write

[TABLE]

For the first term here we have

[TABLE]

and for the second

[TABLE]

where the first two terms are given by Eq. (60) and the last two terms by Eq. (62).

Because there is no Coulomb term coupling leptons with protons the nucleonic chemical potentials in the core will depend only on $\mathcal{E}_{\mathrm{Sky}}$. Given that we are dealing with a homogeneous system, the neutron chemical potential becomes

[TABLE]

(These equations can be obtained as the homogeneous limit of Eq. (B20) of Pearson et al. 2012.) In these expressions there appears the following partial derivative:

[TABLE]

For the leptonic chemical potentials in the core see Appendix B.

The equilibrium composition at any given nucleonic density $n$ is characterized by the parameters $Y_{\mathrm{e}}$ and $Y_{\mu}$, but rather than minimize numerically the total energy per nucleon $e$, as given by Eq. (29), with respect to these two parameters we solve the equations of beta equilibrium

[TABLE]

where $Y_{\mathrm{p}}$ is given by the neutrality condition (28).

3 Results

3.1 Outer crust

As explained in Section 2.1, the EoS and the composition of the outer crust are determined by minimizing the Gibbs free energy per nucleon $g$ at constant pressure $P$. We begin this process at $P=9\times 10^{-12}$ MeV fm*-3*, thereby ensuring a density greater than 106 g cm*-3*, a sufficient condition for complete ionization and degeneracy of the electron gas. The process is then repeated with $P$ increasing in steps of $\Delta P=0.003P$ until neutron drip sets in, as given by Eq. (9). This marks the passage to the inner crust.

The numerical results are summarized in Tables 3-6, where we list the nuclei that minimize $g$, the minimum and maximum densities $\bar{n}_{\mathrm{min}}$ and $\bar{n}_{\mathrm{max}}$ at which that nucleus is present, the pressure $P_{\mathrm{max}}$ and the chemical potentials corresponding to $\bar{n}_{\mathrm{max}}$.

The variation of the energy per nucleon $e_{\mathrm{eq}}$, as given by Eq. (3), with the mean number density $\bar{n}$ is shown in Fig. 4. The four models coincide up to $\bar{n}\approx 3\times 10^{-5}$ fm*-3*, since in this low-density regime we use the masses of experimentally known nuclei, while for higher densities the differences are only just discernible in the figure. The lower panel shows the deviations of the calculated energies $e_{\mathrm{eq}}$ from the analytic parametrization (75). Given that this continuous analytic function was fitted to the calculated energies over the entire density range encountered in neutron stars, the deviations are remarkably small.

The variation of the pressure $P$ with the mean number density $\bar{n}$ is shown in Fig. 5. Again, the four models coincide up to $\bar{n}\approx 3\times 10^{-5}$ fm*-3*, while for higher densities the differences are barely discernible in the figure. The steps at the transition from one nucleus to another are, however, clearly visible, the pressure remaining constant across the transition. Similar steps exist for the chemical potential, because at thermodynamic equilibrium the chemical potential, like the pressure, should be the same on both sides of the interface between two phases in contact, where the density has a discontinuity. (No such steps are visible in Fig. 4 since the energy per nucleon continues to rise across the transition.) The lower panel of Fig. 5 shows the deviations of the calculated pressures from the analytic parametrization (78). The smallness of these deviations is made explicit in Fig. 6 for the case of the functional BSk24. The continuous parametrizing function clearly fails to fit the steps, but otherwise the agreement is excellent, especially in view of the fact that the function was fitted to the calculated pressures over the entire density range encountered in neutron stars.

A more accurate (but discontinuous) parametrization that fits even the steps in the $P-\bar{n}$ curves is given by the analytical thermodynamic model of a fully ionized Coulomb plasma (e.g., Potekhin & Chabrier 2010) with the composition determined by Tables 3 – 6.

Turning now to the composition, we plot in Fig. 7 the variation of the proton fraction $Y_{\mathrm{p}}=Z/A$ as a function of the mean density $\bar{n}$, and see again that the differences between the four functionals are very small. The dominant role of nuclear-structure effects is evident in this figure, but certain weak trends can be discerned in the differences between the four functionals. Comparing BSk22, BSk24 and BSk25 shows generally, but not everywhere, a tendency for higher $J$ to be associated with nuclei that are slightly more neutron-rich. This could be a consequence of the anticorrelation between symmetry energy and $J$ at subnuclear densities, apparent in Fig. 3. On the other hand, we see that BSk24 and BSk26, which have the same value of $J$, follow each other very closely, suggesting that the NeuM constraint is of even less significance than the constraint to $J$.

Fig. 7 also makes apparent a general feature found in all calculations on the outer crust, regardless of the nuclear mass model that is adopted: as we pass from one layer to another with $\bar{n}$ increasing, $Y_{\mathrm{p}}$ always decreases. This is a consequence of the fact that in the outer crust the pressure is determined mainly by the electrons, there being no free neutrons, and the lattice and ion-electron contributions being small. Then since mechanical equilibrium requires that the pressure be continuous across the interface between two layers it follows that the electron density must be likewise almost continuous, and with it the mean proton density $\bar{n}_{\mathrm{p}}$. Then, in an obvious notation

[TABLE]

Thus if $\bar{n}_{2}>\bar{n}_{1}$ it follows that $Z_{2}/A_{2}<Z_{1}/A_{1}$.

The variation of $Z$ and $N$ is plotted as a function of $\bar{n}$ in Fig. 8. The composition is identical for the four models up to $\bar{n}\approx 3\times 10^{-5}$ fm*-3*, since in the low-density regime we use the masses of experimentally known nuclei. At higher densities, the sequence of nuclei in Tables 3–6 becomes model dependent, but it is similar for the different models, apart from some missing nuclides. For example, HFB-22 (HFB-25) predicts the presence of the nucleus 76Ni (126Ru), unlike the other models. Without the inclusion of the mass measurements of Welker et al. (2017), 76Ni is replaced by the odd-$A$ nucleus 79Cu from densities $\bar{n}_{\mathrm{min}}=2.78\times 10^{-5}$ fm*-3* to $\bar{n}_{\mathrm{max}}=4.21\times 10^{-5}$ fm*-3* ($P_{\mathrm{max}}=5.63\times 10^{-5}$ MeV fm*-3*). On the other hand, HFB-25 (HFB-22) does not support the presence of 80Ni (120Sr), unlike the other models. Comparing HFB-24 and HFB-26, both having the same $J$ but different $L$, HFB-26 predicts the presence of 124Zr, unlike HFB-24, and a different Sr isotope at the neutron drip. The discrepancies in the predicted nuclei arise from the uncertainties in the masses of neutron-rich nuclei. Indeed, the masses of 76Ni and 120Sr calculated with HFB-22 and HFB-24 differ by 650 keV and 3.58 MeV respectively, while for 80Ni (126Ru) the difference in the theoretical mass calculated with HFB-25 and HFB-24 amounts to 750 keV (260 keV). For 124Zr, the masses calculated with HFB-26 and HFB-24 differ by 810 keV.

However, as can be seen from the last lines of Tables 3-6, all four functionals predict the same element (Sr) at their respective drip points, although not the same isotopes. For convenience, we summarize the relevant parameters of the neutron-drip points in Table 7. We see that both the drip density $\bar{n}_{\mathrm{nd}}$ and the neutron number $N$ increase as $J$ increases. The values for $\bar{n}_{\mathrm{nd}}$ are slightly different from those given by Fantina et al. (2016a) because of the present inclusion of various corrections, as discussed in Section 2.1. Goriely et al. (2013) have discussed how the masses of drip-line nuclei are correlated with $J$. Since mass is the only nuclear quantity on which the composition of the crust depends, one can expect some correlations between $J$ and the composition. This correlation may be masked by the noise arising from the different errors with which the different functionals fit the data and the numerical errors with which masses are calculated with a given functional. Nevertheless, it is possible to establish some correlations between the properties of the crust at the neutron-drip transition and $J$, (or, equivalently, $L$) as discussed by Fantina et al. (2016a).

Concerning the corrections, the screening correction to the electron energy density may, or may not, change the composition in the outer crust (see Chamel & Fantina 2016b; Fantina et al. 2016b for a discussion). Indeed, as shown by Chamel & Fantina (2016b), for some models, e.g., HFB-22, the differences in the Gibbs free energies for some nuclei may be so small that including or not even small corrections leads to different equilibrium nuclei. Nevertheless, the current uncertainties in nuclear masses are generally more important than the effects induced by the different corrections.

Indeed, within the CCM hypothesis, a small mass difference can lead to a significant difference in the equilibrium nuclide at any depth, as exemplified by the measurement of the mass of 82Zn a few years ago (Wolf et al., 2013). It is thus perhaps remarkable that there is so much agreement between the different mass models. For example, all our HFB models predict the presence of 78Ni, 124Mo, 122Zr, and 120Sr in the deepest layers of the outer crust. These nuclides were also found with our earlier mass models HFB-19, HFB-20, and HFB-21 (Wolf et al., 2013) as well as with the HFB model using the Gogny D1M effective interaction (see Pearson et al., 2011). The existence of these nuclei in the crust is supported by other models. These include the RMF model using the TMA parametrization (see Rüster et al., 2006), the model of Duflo & Zucker (1995) (see Roca-Maza & Piekarewicz, 2008), the finite-range droplet model of Möller et al. (1995) (see Kreim et al., 2013), and the recent BNN model of Utama et al. (2016). Descending through the outer crust to the neutron-drip point, nuclei become more and more exotic (see Fig. 7). For this reason, the composition of the deepest layers remains uncertain. The neutron enrichment may lead to a profound change in the nuclear structure that can hardly be anticipated by phenomenological approaches or models that were adjusted to stable nuclei only. It is therefore not surprising that such models generally agree in predicting the occurrence at the bottom of the outer crust of nuclei with the standard neutron magic number 82, but with a wide spread in the proton number ranging from 32 to 40 (see Rüster et al., 2006). On the other hand, it is noteworthy that all the aforementioned HFB models, which were constrained to reproduce realistic EoSs of NeuM, consistently predict the existence of strontium isotopes at the bottom of the outer crust. More importantly, although some mass models such as that of Duflo & Zucker (1995), the finite-range droplet and BNN yield a very good fit to known masses, they are not applicable beyond the neutron-drip point and therefore do not permit a unified treatment of neutron-star interiors.

3.2 Inner crust

The complete inner-crust results for each of our four functionals are presented in the supplementary material.

Proton fraction. Fig. 9 shows the variation of the equilibrium value of the proton fraction $Y_{\mathrm{p}}=Z_{\mathrm{eq}}/A$ as a function of $\bar{n}$ for our four functionals. Comparing BSk22, BSk24 and BSk25, we see that $Y_{\mathrm{p}}$ has a modest dependence on the symmetry coefficient $J$, increasing $J$ being associated with a decreasing proton fraction, i.e., increasing asymmetry. This is in accordance with the anti-correlation shown in Fig. 3 between $J$ and the symmetry energy $S(\bar{n})$ of INM at subnuclear densities $\bar{n}<n_{0}$; it reflects the fact that higher symmetry energies imply lower asymmetries at beta equilibrium, i.e., higher values of $Y_{\mathrm{p}}$. On the other hand, BSk26 and BSk24, which have the same value of $J$, show similar values of $Y_{\mathrm{p}}$ over the entire inner crust. Thus the choice of the constraining EoS of NeuM has no significant impact on the proton fraction $Y_{\mathrm{p}}$ in the inner crust: if masses are fitted $Y_{\mathrm{p}}$ is determined entirely by the symmetry coefficient $J$.

The analytic function (80) has been fitted to our computed values of $Y_{\mathrm{p}}$ for all four functionals; the lower panel of Fig. 9 shows the deviations between this fit and the computed data.

Energy per nucleon. In Fig. 10 we show the variation of the equilibrium energy per nucleon $e_{\mathrm{eq}}$ as a function of $\bar{n}$ for our four functionals. Comparing functionals BSk22, BSk24 and BSk25, we see from Fig. 3 that $e_{\mathrm{eq}}$ is correlated with $J$ in the same way as is the symmetry energy $S(n)$ at subnuclear densities. That is, the higher the symmetry energy (or equivalently, the lower the asymmetry) the higher $e_{\mathrm{eq}}$. In other words, an increase in the symmetry energy is not entirely offset by the reduction of the asymmetry, i.e., by the increase in $Y_{\mathrm{p}}$, at beta equilibrium (this can be shown analytically). Nevertheless, as core densities are approached in the inner crust the gradient of $e_{\mathrm{eq}}$ is greater the lower the value of $e_{\mathrm{eq}}$ ; this may be reflecting the fact that all three functionals tend to the same EoS of NeuM.

On the other hand, the values for $e_{\mathrm{eq}}$ in the inner crust show very little difference between functionals BSk24 and BSk26: once again, it is $J$ rather than the constraining EoS of NeuM that is the relevant factor. Nevertheless, it will be seen that the rate of variation of $e_{\mathrm{eq}}$ with $\bar{n}$ is somewhat steeper for BSk26 than for BSk24. This will have consequences for the pressure, as we will discuss below.

Our computed values of $e_{\mathrm{eq}}$ in the inner crust have been fitted by the same analytic function (75) with which we fit the two other regions of the neutron star; the lower panel of Fig. 10 shows the deviation between this fit and the computed data.

Pressure. The variation of the equilibrium pressure $P$ with $\bar{n}$ is shown in Fig. 11 for our four functionals. Comparing the first three functionals, we see that the higher $J$ the stiffer is the EoS, in the sense that the pressure rises more rapidly with $\bar{n}$; the differences in the pressure at any given density are simply reflecting the differences in the gradients of the $e_{\mathrm{eq}}$ curves of Fig. 10.

On the other hand, although $J$ has the same value for BSk24 and BSk26, the latter shows throughout most of the inner crust a consistently higher pressure than the former; it is only in the homogeneous core that BSk26 begins to exert a lower pressure than BSk24. This may cause some surprise, given that BSk26 has been fitted to a softer EoS for NeuM than has BSk24, but it can be traced to the purely nuclear properties of INM, as follows. From Eq. (2) we have

[TABLE]

This shows that constraining BSk26 at higher densities, $n>n_{0}$, to a softer EoS of NeuM than BSk24 can only be achieved by having a lower value of the factor $K_{v}+\eta^{2}K_{\mathrm{sym}}$, given that $L$ has roughly the same value for the two functionals. Since our functionals are all fitted to the experimental value of the incompressibility $K_{v}$, the coefficient $K_{\mathrm{sym}}$ will have to be much lower for BSk26 than for BSk24. Table 1 shows that this is indeed the case. It follows that at sufficiently low subnuclear densities, $n<n_{0}$, $P_{\mathrm{INM}}$ will be higher for BSk26 than for BSk24 (we find that this happens for $n/n_{0}<0.72)$.

Our computed values of $P$ in the inner crust have been fitted by the same analytic function (78) with which we fit the two other regions of the neutron star; the lower panel of Fig. 11 shows the deviations between this fit and the computed data.

Equilibrium value of $Z$. Fig. 12 shows how the value of $Z_{\mathrm{eq}}$ varies throughout the inner crust for functionals BSk22, BSk24, BSk25 and BSk26; a zoom of this figure in the free-proton region is shown in Fig. 13 for functionals BSk22, BSk24 and BSk25, and in Fig. 14 for functionals BSk24 and BSk26. While there is no evidence of any shell structure in Figs. 9 – 11 for $Y_{\mathrm{p}},e_{\mathrm{eq}}$ and $P$, respectively, we see from Fig. 12 that before proton drip sets in the value of $Z_{\mathrm{eq}}$ expressed as a function of $\bar{n}$ has a strong shell structure, which differs significantly according to the value of $J$. Furthermore, this $J$-dependent shell structure is superimposed on smoothly varying ETF estimates of the optimal values of $Z$ (dotted lines in Fig. 12) that themselves are strongly $J$-dependent.

Exceptional stability is seen in Fig. 12 for $Z$ = 20 and 40 in the case of functional BSk22, 40 for BSk24 and BSk26, and 40 (at the neutron drip line), 50, 58 and 92 for BSk25 (see Tables 8-11). These ‘magic numbers’ are reflected as local minima or at least as kinks in the plots of the energy per nucleon $e$ as a function of $Z$ that we show in Fig. 15 for two different mean densities $\bar{n}$; note that at each $Z$ the energy is minimized with respect to the neutron number $N$. The disappearance of the familiar magic numbers 28 and 82 associated with bound finite nuclei and the appearance of the unfamiliar numbers 58 and 92 is related to the very strong quenching of the proton spin-orbit splitting that we find in the presence of a large neutron excess. Experimental evidence for such a quenching has been found in neutron-rich bound finite nuclei (Schiffer et al., 2004).

While there are sharp variations in the equilibrium value of $Z$ as a function of density, the fact that the proton fraction $Y_{\mathrm{p}}$ varies smoothly without any apparent shell effects (see Fig. 9) means that the shell effects seen in $Z_{\mathrm{eq}}$ are manifesting themselves primarily as sharp variations in the size of the WS cells.

We stress that all these calculations assumed complete degeneracy, i.e., effectively zero temperature. However, the calculations of Onsi et al. (2008) suggest that the shell effects observed here should survive at realistic crust temperatures of 0.01 MeV.

Proton drip. Unlike the calculations of Baym, Bethe & Pethick (1971b), based on the compressible liquid-drop model, we find that proton drip can occur in the inner crust for all four of our functionals, Eq. (23) becoming satisfied at the densities indicated in Table 12. It will be seen from Table 1 that these drip densities are tightly correlated with $L$, more so than with $J$, despite the loose correlation that exists between these coefficients (see Roggero et al. 2018).

To understand why we should find proton drip while Baym et al. (1971b) do not, it should be noted that rather than have a continuously varying density within the WS cell, as we do, the model of Baym et al. (1971b) adopts a simpler, but less realistic, picture of just two distinct homogeneous phases. The tendency for ‘protons to uncluster’ was observed by Buchler & Barkat (1971); Barkat, Buchler & Ingber (1972) from TF calculations. More recently, the appearance of free protons was also found in the ETF calculations of Martin & Urban (2015).

In the region of free protons the optimal value $Z_{\mathrm{eq}}$ of $Z$ is less well determined than in the bound-proton region, as a comparison of Fig. 16 with Fig. 15 shows. At any given density we take for our final value of $Z_{\mathrm{eq}}$, as displayed in Figs. 13 and 14, the arithmetic mean of the lowest and highest values of $Z$ for which the minimum value (to six significant figures) of the energy per nucleon $e$ is found.

The curves in the free-proton region of the upper panels of Figs. 13 and 14 represent the fit of the computed values of $Z_{\mathrm{eq}}$ to Eq. (79) and of $Z_{\mathrm{cl}}$ to Eq. (81). The lower panels of these figures show the deviations between these fits and the computed data.

Consulting the supplementary material shows that over the range of uncertainty in the value of $Z_{\mathrm{eq}}$, the proton fraction $Y_{\mathrm{p}}\equiv Z/A$ remains sensibly constant. This translates the uncertainty in $Z_{\mathrm{eq}}$ into an uncertainty in the cell size, a circumstance propitious to the development of non-spherical cell shapes, such as pasta. We return to this question in Section 3.5.

The absence of shell effects above the proton-drip density is apparent in Figs. 13 and 14, and reflects our decision to drop the proton shell corrections, along with the pairing term, in this density domain. The rationale for this was discussed in Section 2.2, and is altogether consistent with our treatment of the neutronic shell and pairing corrections. In any case, referring to the supplementary material shows that there is a general tendency for the Strutinsky-integral shell correction $E^{\mathrm{sc,pair}}_{\mathrm{p}}$ to decrease as the density increases to the point where proton drip sets in.

Consulting now Table 12 for the specific cases of functionals BSk22, BSk24 and BSk25, we see that the proton drip density increases as the symmetry energy $J$ decreases. Fig. 13 likewise shows that $Z$ is consistently higher the lower $J$. From Fig. 9 it is seen that the anti-correlation found between $Y_{\mathrm{p}}$ and $J$ at lower densities is maintained in the free-proton region.

Open symbols in Fig. 12 and its zooms (Figs. 13 and  14) denote the cluster component $Z_{\mathrm{cl}}$ of $Z_{\mathrm{eq}}$, given by Eq. (18a). These data points have been fitted by the analytic function (81). The difference

[TABLE]

represents the number of what can be described as free protons, i.e., those associated with the background term in Eq. (14). It should be realized that this number is non-vanishing, although small, even at densities below proton drip. In any case, the free proton fraction $Y_{\mathrm{pf}}=Z_{f}/A$ rises very rapidly with $\bar{n}$ above proton drip. The curves in Fig. 17 represent analytic fits to the computed $Y_{\mathrm{pf}}$, in which we combine Eq. (42) with Eqs. (79) and (81), and use Eqs. (79) and (80) for $A=Z_{\mathrm{eq}}/Y_{\mathrm{p}}$.

Equilibrium value of $N$. Given that the proton fraction $Y_{\mathrm{p}}$ varies smoothly it follows that the number of neutrons $N_{\mathrm{eq}}$ per WS cell in the inner crust at equilibrium must display the same shell effects as does the proton number $Z_{\mathrm{eq}}$. This is confirmed in Figs. 18 and 19. The discontinuities seen in $N_{\mathrm{eq}}$ are a direct consequence of the proton shell effects, given that $Y_{\mathrm{p}}$ varies continuously. The analytic expressions for $N_{\mathrm{eq}}$ represented by the upper set of curves in Figs. 18 and 19 correspond to the quantities appearing in Eq. (84).

We can also define $N_{\mathrm{cl}}$, the cluster component of $N_{\mathrm{eq}}$, given by Eq. (18b), and then

[TABLE]

the number of free neutrons in the WS cell. With the ratio $Y_{\mathrm{nf}}=N_{f}/A$ parametrized according to Eq. (85), $N_{\mathrm{cl}}$ becomes parametrized through Eq. (86), which defines the lower set of curves in Figs. 18 and 19.

The computed ratio $Y_{\mathrm{nf}}=N_{f}/A$ is plotted separately in Fig. 20; this figure also shows the ratio $Y_{\mathrm{n}}=N_{\mathrm{eq}}/A=1-Y_{\mathrm{p}}$.

Neutron chemical potential. The neutron chemical potential $\mu_{\mathrm{n}}$ is calculated (for configurations in beta equilibrium) using Eq. (7), and its variation over the inner crust is shown in Fig. 21 for our four functionals (with the rest-mass energy of the neutron subtracted). For all four functionals the trend is for $\mu_{\mathrm{n}}$ to increase monotonically with $\bar{n}$, following the growing neutron excess.

The analytic function (92) has been fitted to our computed values of $\mu_{\mathrm{n}}$; the lower panel of Fig. 21 shows the deviations between this fit and the computed data.

Proton chemical potential. This is calculated (for configurations in beta equilibrium) from the calculated values of $\mu_{\mathrm{n}}$ and $\mu_{\mathrm{e}}$, using Eq. (8), and its variation over the inner crust is shown in Fig. 22 for all four of our functionals (with the rest-mass energy of the proton subtracted). The general trend for $\mu_{\mathrm{p}}$ is the opposite to that for $\mu_{\mathrm{n}}$: it becomes more and more negative with increasing $\bar{n}$, following the growing neutron excess. The $J$-dependence also is reversed, with higher $J$ being associated with higher $\mu_{\mathrm{p}}$, simply because $Y_{\mathrm{p}}$ and hence $\mu_{\mathrm{e}}$ are lower for higher values of $J$.

3.3 Matching between inner and outer crust.

Our inner-crust code, as used here, is in principle applicable to the outer crust, and it is thus meaningful to compare this code with the HFB code that we used for the outer-crust calculation of Section 3.1. In Table 13 we make this comparison at the drip-point density for the functional in question; the results for the outer-crust code are shown in parentheses. The values for $\bar{n}_{\mathrm{nd}}$ and $P$ from our outer-crust code are slightly different from those previously obtained by Fantina et al. (2016a) without all the corrections considered here.

It will be seen that there is a slight disagreement between our inner- and outer-crust codes concerning the values of $Z$ and $N$ at the drip point: for all four functionals we find $Z$ = 40 rather than 38, and somewhat bigger discrepancies for $N$. This can be attributed to the several approximations made in our ETFSI method for calculating the internal energy per nucleon $e$ in the inner crust; we see from Table 13 that the inner-crust code (ETFSI) underbinds with respect to the outer-crust code (HFB) by around 4 %; this rises to about 5% if we take the same values of $Z$ and $N$ in the two codes. (The corresponding figures for the pressure are 2% and 5%.)

We recall that the approximations made in the ETFSI method, but not in the HFB method adopted in our outer-crust calculations, are as follows. i) The kinetic energy and spin currents are calculated with the semiclassical ETF method. ii) Proton shell corrections are put in perturbatively, and neutron shell corrections (shown to be much smaller than proton shell corrections as soon as neutron drip sets in (Chamel, 2006; Chamel et al., 2007), but obviously not zero, in the outer crust) are neglected completely. This source of discrepancy between the two codes will be maximal at the interface between the inner and outer crust, since the neglected neutron shell effects will decrease as the density increases. iii) Rather than allowing arbitrary density variations when minimizing the total energy, the density is parametrized according to Eqs. (14) and (15). iv) Proton pairing is treated approximately, while neutron pairing is neglected completely. v) The collective and Wigner correction terms that are added to the HFB energy in the outer-crust calculations are neglected in the inner crust.

On the other hand, we have checked that errors arising from the assumption of sphericity in the inner-crust code are negligible in this region of the nuclear chart.

We must stress that even if there is a disagreement of around 4% between the two codes the energy differences between adjacent values of $Z$ and $N$ calculated by the inner-crust code are much more precise: the perceived regularities suggest at least 5-figure, and possibly 6-figure accuracy: see, for example, Fig. 15.

3.4 Core

Our results for the equilibrating values in the core of the number of protons $Y_{\mathrm{p}}$ and of the number of electrons $Y_{\mathrm{e}}$ per nucleon, calculated as described in Section 2.3, are shown for our four functionals in Fig. 23; the difference between the two sets of curves represents the number of muons $Y_{\mu}=Y_{\mathrm{p}}-Y_{\mathrm{e}}$ per nucleon. This figure also indicates the densities corresponding to the breakdown of causality and the maximum neutron-star mass for the different functionals (see Section 4). The lower panel indicates the deviations between our calculated results and the analytic fits to these results given by Eq. (91). The most striking feature of this figure is the fact that in the core both $Y_{\mathrm{p}}$ and $Y_{\mathrm{e}}$ increase with increasing density, in contrast to what happens in both the outer and inner crusts (see Figs. 7 and 9).

Comparing functionals BSk22, BSk24 and BSk25 to look for a possible $J$-dependence, we see from Fig. 2 that at all core densities there is a close correlation of $Y_{\mathrm{p}}$ (and $Y_{\mathrm{e}}$) with the symmetry energy $S(n)$, which is only partially correlated with $J$. There is, however, a tendency for the three functionals (especially BSk22 and BSk25) to converge at high densities, reflecting the fact that all three have been constrained to the same EoS of NeuM. In fact, the BSk26 curve in Fig. 23 shows that the NeuM constraint has much more influence than the choice of $J$: the softer EoS of NeuM favours higher neutron excesses, and thus lower values of $Y_{\mathrm{p}}$ (and $Y_{\mathrm{e}}$).

The energy per nucleon $e_{\mathrm{eq}}$ and the pressure $P$ at equilibrium in the core are shown for our four functionals in Figs. 24 and 25, respectively. In the lower panels of these two figures we see the deviations between the calculated values and the analytic fits of Eqs. (75) and  (78), respectively. We stress once again that these analytic fits to $e_{\mathrm{eq}}$ and $P$ are valid over the entire star. Comparing these two figures with Figs. 10 and 11, respectively, we see that BSk22, BSk24 and BSk25 behave much more similarly in the core than in the inner crust; this simply reflects the fact that all three have been constrained to the same EoS of NeuM. Not surprisingly, in both Figs. 24 and 25, BSk26 is seen to be somewhat softer than the other three functionals, but one might have expected a stronger dependence on the NeuM constraint. The point is that the softer EoS of NeuM is partially offset by the higher asymmetries.

Figs. 26 and 27 show for our four functionals the variation over the core of the neutron chemical potential $\mu_{\mathrm{n}}$ and the proton chemical potential $\mu_{\mathrm{p}}$, respectively. It will be seen that the slope of $\mu_{\mathrm{p}}$ turns positive in this region of the star. A comparison of BSk22, BSk24 and BSk25 shows that the choice of symmetry coefficient $J$ has little systematic impact on the values of the chemical potentials in the core. On the other hand, both $\mu_{\mathrm{n}}$ and $\mu_{\mathrm{p}}$ tend to be lower for BSk26, the functional fitted to the softer EoS of NeuM. In the lower panel of Fig. 26 we see the deviations between the calculated values of $\mu_{\mathrm{n}}$ and the analytic fit of Eq. (94). The deviation shown in the lower panel of Fig. 27 for $\mu_{\mathrm{p}}$ is calculated from the analytic expression (94) for neutrons using the beta-equilibrium condition (8), as explained in Appendix C.

3.5 Crust-core interface

In Table 14 we tabulate for each of our four functionals the values for the density of the crust-core transition $n_{\mathrm{cc}}$, along with the corresponding values of $Y_{\mathrm{e,cc}}$ and the pressure $P_{cc}$, as calculated by the method of Ducoin, Chomaz & Gulminelli (2007). This is the same method that was used by Pearson et al. (2012), (where $n_{\mathrm{cc}}$ was denoted by $n_{\mathrm{trans}}^{N*M}$), and it consists of a determination of the conditions for homogeneous NM to be unstable against breakup into finite-size clusters (rather than the conditions for equilibrium between infinite homogeneous phases). The last column of Table 14 shows $n_{\mathrm{cc}}^{\prime}$, the lowest density for which we have found no inhomogeneous solutions with our inner-crust code, by which we mean solutions with values of the inhomogeneity parameter $\Lambda$, defined in Eq. (22), greater than 10-7* (typically, this value of $\Lambda$ is found by our code for densities known to correspond to NM). For all four functionals we find $n_{\mathrm{cc}}^{\prime}$ to be slightly smaller than $n_{\mathrm{cc}}$, which means that NM remains stable against breakup down to slightly lower densities than predicted by the method of Ducoin et al. (2007). In this respect it should be pointed out that this method is perturbative, assuming small-amplitude fluctuations of the density. On the other hand it is possible that we are missing some inhomogeneous solutions because of the discrete grid that we take for the proton number $Z$ (see Section 2.2), although it is unlikely that this would happen for all four functionals.

For all four functionals we begin to find, as the density approaches $n_{\mathrm{cc}}^{\prime}$, equilibrium solutions that are still sensibly inhomogeneous, with values of $\Lambda$ of the order of 10*-2*), but that are mechanically unstable, in the sense that the calculated pressure decreases with an increase in $\bar{n}$. The situation is resumed in Table 15, where we show for each functional the highest density beyond which instability sets in. Instability of this sort is typical of what happens at a first-order phase change when minimizing the Helmholtz free energy, rather than the Gibbs energy. The actual physical situation could be restored by a Maxwell construction, but this is beyond the scope of the present paper. Thus all our inner-crust results for densities higher than those indicated in Table 15 are excluded from the supplementary tables and from the results shown in Section 3.2.

It is nevertheless of interest to speculate on the nature of these instabilities in our inner-crust results. Phase changes occur throughout the inner crust, at every point where the number of protons $Z$ changes, but only at higher densities near the transition to homogeneous matter do we see any instabilities. Thus the phase changes associated with changes in $Z$ must either be of higher order, or else, if of first order, too weak to be visible in our calculations. This suggests that the instabilities that we do see in our calculations involve something more significant than a change in $Z$, and an obvious possibility is that changes to non-spherical pasta phases are being signalled, in the sense that this is what we would see if our code allowed it. This interpretation is strengthened by the fact that some of our solutions in the unstable region correspond to spherical bubbles, which means that non-spherical configurations such as lasagna or spaghetti might be favoured at slightly lower densities, since spherical bubbles are generally the densest types of pastas to appear (Ravenhall et al., 1983; Hashimoto et al., 1984). However, with the WS cells of our code being limited to spherical shapes we are unable to investigate this question any further here.

It is, of course, possible that the apparent bubbles we find may simply be a result of numerical error, since it will be seen from the supplementary material that in this region the pressure can vary significantly as $Z$ varies while the energy per nucleon $e$ remains more or less constant. Thus a slight error in the calculation of the latter could lead to a much larger error in the pressure. Another possibility is that the apparent instabilities and bubbles result from our restricted parametrization of the density profile. Nevertheless, the fact that we find bubble-like solutions emphasizes the need for a generalization of the present work to include the possibility of non-spherical pasta phases within the full fourth-order ETFSI framework.

4 Gross properties of neutron stars

4.1 Mass-radius relation

We have compared the mass $M_{\mathrm{NS}}$ and radius $R_{\mathrm{NS}}$ of a neutron star, calculated using the tabulated EoSs and their analytical representations (78). We integrated the Tolman-Oppenheimer-Volkoff (TOV) equation from the centre, with the central mass-energy density $\rho_{\mathrm{c}}$ as a free parameter, outward to $\rho=10^{6}$ g cm*-3*, using the fourth-order Runge-Kutta method with an adaptive step and controlled accuracy. Since the adaptive step does not generally coincide with the interval between the density points that are given in the EoS tables, we used points obtained by interpolation in the tables. Two different methods of interpolation were used, linear and cubic, the two methods agreeing to 5 significant figures in $M_{\mathrm{NS}}$. The star was assumed spherically symmetric, nonrotating and nonmagnetized. The outermost layer of $\rho<10^{6}$ g cm*-3*, which is excluded from consideration, is unimportant for the gross properties of a neutron star, because its thickness does not exceed a few meters, and its mass is only $\sim 10^{-12}\,M_{\odot}$.

Figure 28 shows the mass-radius relation for the EoSs BSk22, BSk24, BSk25, and BSk26. The neutron-star configurations obtained with the original EoSs and with their analytical representations are drawn as symbols and lines, respectively. The maximum neutron-star masses $M_{\mathrm{NS}}^{\mathrm{stab}}$, the corresponding radii $R_{\mathrm{NS}}^{\mathrm{stab}}$, central number densities of baryons $n_{\mathrm{c}}^{\mathrm{stab}}$, and central mass densities $\rho_{\mathrm{c}}^{\mathrm{stab}}$ are listed in Table 16. The values of $M_{\mathrm{NS}}^{\mathrm{stab}}$, calculated using the tables and the fits, differ by less than 0.06% for each of the four EoSs. For configurations with $n_{\mathrm{c}}>n_{\mathrm{c}}^{\mathrm{stab}}$ ($\rho_{\mathrm{c}}>\rho_{\mathrm{c}}^{\mathrm{stab}}$), the condition of hydrostatic stability $\mathrm{d}M_{\mathrm{NS}}/\mathrm{d}\rho_{\mathrm{c}}>0$ is violated; the unstable configurations are shown by the parts of the curves to the left of the maxima in Fig. 28, marked by the empty circles.

Comparing the data listed in Table 16 for BSk22, BSk24 and BSk25, we see that a decrease of $J$ is accompanied by an increase of the number density at the centre of the most massive stable configuration, $n_{\mathrm{c}}^{\mathrm{stab}}$, although the corresponding mass-energy density $\rho_{\mathrm{c}}^{\mathrm{stab}}$ is almost unaffected. This can be traced back to the correlation between $J$ and the behavior of the symmetry energy at suprasaturation densities, as shown in Fig. 2. Now looking at the last line of the table, we see that the NeuM constraint has a bigger influence on the most massive stable configuration, as expected: the softer EoS BSk26 corresponds to the smaller $M_{\mathrm{NS}}^{\mathrm{stab}}$ and larger number density and mass-energy density. Very recently, gravitational-wave observations have been interpreted as indicating that the maximum mass of non-rotating neutron stars cannot exceed $2.16^{+0.17}_{-0.15}$ solar masses (Rezzolla, Most & Weih, 2018), supporting our assumption that the real EoS cannot be much stiffer than LS2.

In Table 17 we list the radii $R_{\mathrm{NS}}^{\mathrm{1.4}}$, central number densities $n_{\mathrm{c}}^{\mathrm{1.4}}$ and central mass densities $\rho_{\mathrm{c}}^{\mathrm{1.4}}$ of configurations with $M_{\mathrm{NS}}=1.4\,M_{\odot}$. All the radii except that of functional BSk26 agree very well with the recent estimate by Margueron, Hoffmann Casali & Gulminelli (2018) of 12.7$\pm 0.4$ km. The first three lines of the table show that for a given constraining EoS of NeuM decreasing $J$ (and thus decreasing $L$) is accompanied with decreasing $R_{\mathrm{NS}}$ and increasing central density. This reflects the correlation between $L$ and $R_{\mathrm{NS}}$ known from previous studies (e.g., Fortin et al., 2016; Margueron et al., 2018). Comparing BSk24 and BSk26 in order to assess the role of the EoS of NeuM shows a bigger difference in the radius, a difference large enough to be observable (Psaltis, Özel & Chakrabarty, 2014). As expected, this pattern is found not only for the radii at a fixed mass, but also for the radii of the most massive stable configurations, $R_{\mathrm{NS}}^{\mathrm{stab}}$ (Table 16). All models are consistent with the radius constraint inferred from the recent detection of gravitational waves from the binary neutron star merger GW170817 (Annala et al., 2018; De et al., 2018; Fattoyev et al., 2018).

4.2 Causality limits

Table 18 lists the largest values of baryon number density ($n^{\mathrm{caus}}$) and mass-energy density ($\rho^{\mathrm{caus}}$), for which the condition $\mathrm{d}P/\mathrm{d}\rho<c^{2}$ is satisfied, i.e., for which the speed of sound is smaller than the speed of light. At higher densities the EoS becomes superluminal and causality breaks down; for this reason $n^{\mathrm{caus}}$ and $\rho^{\mathrm{caus}}$ are often named causality limits (see, e.g., Section 5.15 of Haensel et al. 2007).

The crosses in Fig. 28 correspond to $\rho_{\mathrm{c}}=\rho^{\mathrm{caus}}$. The line segments to the left of the crosses correspond to configurations where the central part of the star has $\rho>\rho^{\mathrm{caus}}$. We see that for the functionals BSk22, BSk24, and BSk25 these segments correspond to configurations that are unstable anyway, since $n_{\mathrm{c}}^{\mathrm{stab}}<n^{\mathrm{caus}}$, so there is no breakdown of causality in any situation that would otherwise be physically meaningful. On the other hand, for functional BSk26 we have $n_{\mathrm{c}}^{\mathrm{stab}}>n_{\mathrm{c}}^{\mathrm{caus}}$, which means that over the range of central densities $n_{\mathrm{c}}$ satisfying $n_{\mathrm{c}}^{\mathrm{stab}}>n_{\mathrm{c}}>n^{\mathrm{caus}}$ there will be stable solutions for which causality breaks down in the central region of the star. The minimum mass for which this can happen, corresponding to $n_{\mathrm{c}}=n^{\mathrm{caus}}$, is $M_{\mathrm{NS}}=2.15\,M_{\odot}$ (see the BSk26 panel in Fig. 28). The most massive possible star (for BSk26) has $M_{\mathrm{NS}}=2.17\,M_{\odot}$, and the central region over which causality fails has a maximum radius of 3.6 km and a maximum mass of $0.23\,M_{\odot}$.

Such a breakdown of causality for model BSk26 is, of course, quite unphysical, and it is a consequence of its failure in the realistic APR EoS. The problem is discussed in detail by Akmal et al. (1998), where it is shown that the reduction in the total stellar mass resulting from a restoration of causality is quite small. It should be safe to suppose that the correct limiting mass for functional BSk26 lies somewhere between 2.15 and 2.17 $M_{\odot}$.

4.3 Direct Urca process

Number fractions of the electrons and muons in the core of a neutron star are important in the neutron-star cooling theory, because they determine whether the extremely powerful direct Urca processes of neutrino emission operate or not (e.g., Haensel 1995 and references therein). For strongly degenerate particles, the energy-momentum conservation law requires the condition $n_{\mathrm{e}}^{1/3}+n_{\mathrm{p}}^{1/3}>n_{\mathrm{n}}^{1/3}$ to be fulfilled, in order for these processes to work. For the theoretical models where $Y_{\mathrm{e}}$ monotonically increases with the increase of density, the direct Urca processes can work in the central regions of sufficiently massive neutron stars. Using Eqs. (88) – (91), we obtain the corresponding threshold values of central number densities $n_{\mathrm{DU}}$. Using the results of Sec. C.1, we find the respective mass-energy density values $\rho_{\mathrm{DU}}$. Using Eq. (78) and solving the TOV equation with these central mass densities, we find the minimum mass $M_{\mathrm{DU}}$ for a neutron star to cool rapidly via the direct Urca processes. The values of $n_{\mathrm{DU}}$, $\rho_{\mathrm{DU}}$, and $M_{\mathrm{DU}}$ for the considered EoS models are listed in Table 19. In the case of BSk26, the $M_{\mathrm{DU}}$ value is enclosed in brackets, which indicates that in this case $n_{\mathrm{DU}}>n_{\mathrm{c}}^{\mathrm{stab}}$, so that the corresponding spherically symmetric configuration belongs to the unstable branch in Fig. 28. Thus the direct Urca processes cannot work in stable neutron stars described by the BSk26 functional.

Observations indicate that the direct Urca processes operate in a relatively small number of neutron stars. None of the around forty isolated neutron stars with measured ages and thermal luminosities are as cool as required by the “rapid cooling” models (see, for example, Potekhin & Chabrier 2018 and references therein). Most of them can be explained within the so called “minimal cooling paradigm”, not involving the direct Urca process (Page et al., 2004; Gusakov et al., 2004). However, a few objects show thermal luminosities between the minimal and rapid cooling predictions, while most of the magnetars are considerably hotter than predicted by both models. These exceptions may be explained by internal heating, which can be provided by a number of physical mechanisms (see, e.g., Potekhin, Pons, & Page 2015, for review and references). Thermal luminosities of accreting neutron stars in soft X-ray transients in quiescence can be explained only if the direct Urca processes are forbidden in the hottest of them (Aql X-1, 4U 1608$-$52) but allowed in the coldest one (SAX J1808.4$-$3658) (Yakovlev et al., 2004). Recently, observations of thermal relaxation of the neutron star in the transient system MXB 1659$-$29 have delivered an evidence of the direct Urca processes operating in its core (Brown et al., 2018).

The low value of $M_{\mathrm{DU}}$ given by the BSk22 functional implies that the direct Urca processes would operate in most neutron stars, which is hardly compatible with the fact that at least a large fraction of them are well described by the minimal cooling model. Thus the BSk22 model, which in any case gives the worst atomic mass fit of our four functionals, is disfavored by neutron-star observations. On the other hand, the presence of the direct Urca processes in some neutron stars rules out the BSk26 functional. Thus, among our four functionals, only BSk24 and BSk25 provide $n_{\mathrm{DU}}$ values compatible with observations.

5 Concluding remarks

We have presented a microscopic treatment of the nuclear physics of the outer crust, the inner crust and the core of neutron stars. Our treatment is unified in the sense that the EoS (i.e., the pressure as a function of the density), the composition and the chemical potentials of all three regions are calculated with the same energy-density functional. We have performed calculations of the entire neutron star (assumed to be completely degenerate and non-accreting) with four different functionals, BSk22, BSk24, BSk25 and BSk26 (Goriely et al., 2013), which are based on generalized Skyrme-type forces and contact pairing forces. These functionals were precision fitted to essentially all the available atomic mass data (with $Z,N\geq 8$) by using the HFB method. In addition, these functionals were constrained to fit, up to the densities prevailing in neutron-star cores, the EoS of homogeneous pure NeuM. Since this latter EoS is by no means uniquely determined by our present knowledge of nuclear physics at high densities we considered two quite different EoSs, the hard EoS ‘V18’ of Li & Schulze (2008) that we label LS2 here, and the soft EoS ‘A18 + $\delta\,v$ + UIX∗’ of Akmal et al. (1998) that we label APR: functionals BSk22, BSk24 and BSk25 were all fitted to LS2, while BSk26 was fitted to APR. We have argued that the real EoS of NeuM can probably not be much stiffer than LS2, and certainly not much softer than APR.

The NeuM constraint, even combined with the atomic mass fit, does not completely determine the symmetry energy, allowing us thereby some flexibility on the symmetry coefficient $J$. Accordingly, BSk22, BSk24 and BSk25 were fitted to values of $J=$ 32, 30 and 29 MeV, respectively. The quality of the mass fits deteriorates rapidly outside this range and indeed the quality of the fits indicates a slight preference for $J$ = 30 MeV, i.e., BSk24. We did not exploit this element of flexibility in $J$ under the constraint of the APR EoS of NeuM, but instead imposed the unique value of 30 MeV, constructing thereby functional BSk26. Thus to study the role of the NeuM constraint one has to compare BSk24 with BSk26, while for the role of the symmetry coefficient it is BSk22, BSk24 and BSk25 that have to be compared.

To calculate the properties of neutron stars from the given functionals we use the HFB method in the outer crust (except when the appropriate atomic mass data are available), while for the locally homogeneous core the calculation is essentially exact. For the inner crust we use the ETFSI approximation to the HF method, with pairing handled in the BCS approximation. The first, semi-classical, stage of our implementation of the ETFSI method is based on the picture of spherically symmetric WS cells, and adopts a parametrization of the nucleonic density distributions that respects the boundary conditions necessary both for the validity of the fourth-order ETF formalism with which we calculate the energy and for continuity between adjacent cells. The second stage calculates proton shell effects perturbatively, but neglects the much smaller neutron shell effects, and thereby avoids the inevitable problems associated with the continuum in the WS approach.

Concerning our results for the composition, the proton fraction $Y_{\mathrm{p}}$ is nowhere in the star very sensitive to $J$ (see Fig. 30), but in the core it depends strongly on the NeuM constraint. While the value of $J$ influences only weakly $Y_{\mathrm{p}}$, the actual values of $N$ and $Z$ in the inner crust depend strongly on it, doing so both through the semi-classical ETF part of the energy and also through the strong proton shell effects: the magic numbers that emerge depend strongly on $J$. We stress here that the range of values of $J$ that we consider is the widest possible consistent with a good fit to the atomic mass data: this latter condition constitutes a very stringent constraint.

Turning to our results for the EoS, it might appear from Fig. 29 that the dependence of the EoS on both $J$ and the NeuM constraint is rather modest. However, we have shown in Section 4 that these apparently small differences have a significant and potentially observable impact on the global structure of neutron stars. All our results for the EoS, the composition and the chemical potentials have, for the convenience of modellers, been fitted by analytic expressions for all four functionals.

A limitation of our calculations is the restriction to spherical cells, which makes it impossible to study the question of nuclear pasta. The very existence of these non-spherical cell shapes, which at the most can be favoured only in a narrow shell close to the crust-core interface, is controversial because of the tiny energy differences between spherical and non-spherical cell shapes, but for precisely this reason the impact of pasta on the EoS is very weak. Nevertheless, pasta strongly influences other properties of neutron stars, and it is important to know whether or not it exists. A more detailed study, in which the limitation of spherical cell shapes is removed, is desirable.

After the construction of the functionals of this paper (Goriely et al., 2013) a new set was published, BSk30, BSk31 and BSk32 (Goriely et al., 2016). These functionals were characterized by an improved treatment of the pairing, with self-energy effects being included. In this way more realistic gaps for INM were obtained, an essential condition for the study of superfluidity in neutron stars, but to maintain the high quality of the mass fits of our older functionals it was necessary to add a phenomenological pairing term dependent on the density gradient. The three new functionals were fitted to $J$ = 30, 31 and 32 MeV, respectively, all being fitted to the LS2 EoS of NeuM. So far, no functional with the new pairing has been mass fitted under the softer APR NeuM constraint, which means that the new family is not suitable for the sort of study that we have undertaken in the present paper.

However, we can still assess the impact that the improved treatment of the pairing functional would have on the conclusions that we have drawn from the older family of functionals considered here. The optimal mass fit with the new pairing is found with BSk31, for which $J$ = 31 MeV. We therefore compare it with functional BSk23 (Goriely et al., 2013), which also has $J$ = 31 MeV and is constrained to the LS2 EoS of NeuM, but belongs to the older family of functionals in that it has the same form of pairing. Preliminary calculations show that the effect of changing the pairing functional on the EoS and on the composition is much smaller than the effect of changing $J$ from 31 to 29 MeV, or of changing the constraining EoS of NeuM from LS2 to APR. Nevertheless, given the importance of the improved treatment of pairing for reliable calculations of neutron-superfluid properties, we intend to publish a more complete study of the new functionals in a later paper.

In the meantime, we believe that the results we have presented here for functionals BSk22, BSk24 and BSk25 span the current range of uncertainty associated with the gaps in our knowledge of the symmetry coefficient $J$ (and of the symmetry-slope coefficient $L$, which is roughly correlated to $J$ by the fits to nuclear masses). Comparing, on the other hand, functionals BSk24 and BSk26 gives some insight into the impact of the uncertainty in the EoS of NeuM. Not surprisingly, the influence of the uncertainty in $J$ is greatest in the crust and diminishes as we move to the core of the star, while the contrary is the case for the uncertainty in the EoS. Thus we can conclude with some confidence that we have tied down the nuclear uncertainties (aside from those related to pairing) in the crust of the star. This confidence decreases somewhat as we approach the centre of the core, but it is unlikely that the actual properties lie far outside the limits that we have considered here.

Acknowledgements

We wish to thank M. Onsi for her invaluable contributions to earlier stages of this project, and the referee for valuable remarks. J. M. P. acknowledges the partial support of the NSERC (Canada). The work of A.Y.P. was partially supported by the RFBR grant 16-29-13009-ofi-m. The work of N. C. was partially supported by the Fonds de la Recherche Scientifique - FNRS (Belgium) under grant n∘ CDR-J.0187.16 and CDR-J.0115.18. S. G. acknowledges the support of Fonds de la Recherche Scientiffique-FNRS (Belgium). This work was also partially supported by the European Cooperation in Science and Technology (COST) Actions MP1304 and CA16214.

Appendix A Neutron chemical potential: proof of Eq. (7)

Although the bulk of this paper has been devoted to the inhomogeneities of the outer and inner crusts, we must recognize that this feature holds only at the microscopic level, the isolated nuclei of the former containing no more than a few hundred nucleons, and the WS cells of the latter at most a few thousand. Thermodynamically the crust, like the core, can be regarded as homogeneous, in the sense that over an appreciable range of sizes the mass, energy and entropy (for $T>0$) of a piece of crust (or of the core) are proportional to its volume. It is then easy to show that the total Gibbs energy $G$ of the piece can be written as

[TABLE]

where $\mu_{i}$ is the chemical potential of component $i$ and $N_{i}$ is the number of particles of component $i$ in the piece (see, for example, Section 11.3 of Adkins 1983).

Discounting the possible presence of hyperons, we have just four components in any star: neutrons, protons, electrons and muons (it makes no difference in the following whether or not the nucleons form clusters or nuclei). Then given charge neutrality it follows from Eq. (44) that the Gibbs energy per nucleon is

[TABLE]

If there is beta equilibrium, as in a neutron star, then Eqs. (39a) and (39b) will hold, and the last two terms of Eq. (45) vanish, leaving us with Eq. (7). The same conclusion follows if there are no muons, since in that case $Y_{\mathrm{e}}=Y_{\mathrm{p}}$. It is worth noting that Eq. (7) holds for all temperatures, including $T=0$.

Appendix B Treatment of lepton gas

The electron gas in the inner crust and the core, and the muon gas in the latter region, are assumed to be completely uniform and degenerate, with special-relativistic effects treated in all generality. Then following, for example, Chapter 24 of Weiss et al. (2004), we define

[TABLE]

where $\lambda=\hbar/(mc)$ is the Compton wavelength of the lepton (electron or muon), $m$ being the corresponding mass, and $n$ is the lepton number density. Defining also

[TABLE]

we can write the kinetic-energy density of the leptons as

[TABLE]

which reduces in the high-density extreme-relativistic limit to

[TABLE]

Only in the outer crust will deviations from uniformity be non-negligible, and in that region we add to the energy density the screening correction given by the interpolation formula of Potekhin & Chabrier (2000), instead of the Thomas-Fermi expression (Salpeter, 1961) used by Pearson et al. (2011) (see Chamel & Fantina 2016b for a recent discussion on this approximation). Then at temperatures $T$ much lower than the plasma temperature

[TABLE]

where

[TABLE]

is the ion-plasma frequency (the ion mass coincides with the nuclear mass since atoms are fully ionized), the energy density of the screening correction is given by Potekhin & Chabrier (2000)

[TABLE]

Here $x_{\mathrm{e}}$ is the relativity factor (46) for electrons, and

[TABLE]

in which $a_{\mathrm{n}}=(4\pi/3\ n_{\mathrm{e}}/Z)^{-1/3}$ is the inter-ion spacing and $\alpha=e^{2}/(\hbar c)$ is the fine-structure constant. Also, the parameters $s$ and $b_{1}$–$b_{4}$, which depend only on $Z$, are given by Potekhin & Chabrier (2000) as

[TABLE]

where, for a bcc lattice, $a_{1}=1.1866$, $a_{2}=0.684$, $a_{3}=17.9$, $a_{4}=41.5$, and $q=0.205$ (Potekhin & Chabrier, 2000).

Note that the expressions suggested by Potekhin & Chabrier (2000) provide good approximations to the energy and pressure at arbitrary temperature, but they produce an unphysically slow decrease of the electron-screening contribution to the heat capacity with decreasing $T$ below $T_{\mathrm{pl}}$, i.e., for a quantum crystal. The problem was recognized and rectified by Potekhin & Chabrier (2010). However, the expressions of Potekhin & Chabrier (2000) and Potekhin & Chabrier (2010) have the same limit at $T/T_{\mathrm{pl}}\to 0$, Eq. (52), which we use here.

In addition to a correction for non-uniformity of the electron gas in the outer crust we include a correction for Coulomb exchange in all regions of the star. The general expression that we take for the exchange energy density is Eq. (9) of Salpeter (1961), which can be simplified (Engel, 2002), without approximation, to

[TABLE]

We use this expression as it stands for the electrons in the outer crust and the muons in the core, while for the electrons of the inner crust and the core it suffices to take the extreme relativistic limit

[TABLE]

There is no direct Coulomb term in the core, since strict neutrality prevails throughout, but there is in the outer and inner crusts. However, in the outer crust these direct terms constitute the lattice energy $E_{L}(A,Z)$ of Eq. (3), while in the inner crust they are represented by the ${\mathcal{E}}_{\mathrm{ee,ep}}^{\mathrm{c}}(\bm{r})$ term of Eq. (10). Thus we do not have to consider any further the direct Coulomb interactions of the electrons, either with themselves or with protons.

As in Pearson et al. (2011), we also include in the outer crust a correction for the electron-correlation energy, using the high-density limit

[TABLE]

where

[TABLE]

(see, e.g., Engel 2002). This correction is quite negligible in the inner crust and core.

Thus, for the total energy density of the electrons, excluding their rest mass, we write

[TABLE]

in which all quantities are as already defined. A similar expression holds for muons.

Pressure. Defining

[TABLE]

the lepton pressure corresponding to the energy density (48) is

[TABLE]

This reduces for the extremely relativistic electrons of the inner crust and core to

[TABLE]

The exchange pressure corresponding to the general expression (54) is

[TABLE]

For the extremely relativistic electrons this reduces to

[TABLE]

In the outer crust there is also an electron-screening contribution to the pressure given by

[TABLE]

with

[TABLE]

The pressure corresponding to the electron-correlation energy, Eq. (56), is given by

[TABLE]

Thus, for the total pressure exerted by the electrons we have

[TABLE]

in which all quantities are as already defined, with a similar expression for muons.

Chemical potentials. For the chemical potential of the electrons we have

[TABLE]

where $e_{\mathrm{e}}={\mathcal{E}}_{\mathrm{e}}/n_{\mathrm{e}}$ is the total energy per electron, excluding the rest mass. With ${\mathcal{E}}_{\mathrm{e}}$ given by Eq. (58), we have

[TABLE]

For the first term here Eqs. (48) and (60) lead simply to the familiar Fermi energy of a free fermion gas with the rest mass subtracted,

[TABLE]

For the second term in Eq. (69) we find from Eqs. (54) and (62)

[TABLE]

which reduces in the extreme relativistic limit to

[TABLE]

For the screening term in Eq. (69) we find from Eqs. (52) and (64)

[TABLE]

while from Eqs. (56) and (66) we find for the last term in Eq. (69)

[TABLE]

All these expressions for the electron chemical potential hold in all regions of the neutron star, and equally well for the muon chemical potential. Note, however, that since we neglect the last two terms in Eq. (69) everywhere except in the outer crust they have no relevance for muons.

Appendix C Analytic representations

For each of our four functionals BSk22, BSk24, BSk25, and BSk26, we construct parametrizations in terms of analytic functions for the quantities of astrophysical interest that we have calculated, following the approach previously developed by Haensel & Potekhin (2004); Potekhin et al. (2013). In the case of the energy per nucleon and pressure we were able to find unified fits to a single continuous analytic function of the number density $\bar{n}$ over a broad density range covering the entire neutron star, i.e., the outer and inner crusts and the core (in the core we have $\bar{n}\equiv n$). These unified fits smear away all density discontinuities between layers of different composition, and can be useful for hydrodynamic modeling. Other quantities that are required in modeling neutron-star structure and evolution, such as particle fractions, are given by separate parametrizations for the inner crust and the core. All parametrizations presented in this Appendix have been implemented in Fortran subroutines, freely available at the Ioffe Institute website.111http://www.ioffe.ru/astro/NSG/BSk/

C.1 Energy as a function of density

For each of the four functionals, the equilibrium energy per baryon that we have calculated is fitted to a single continuous analytic function of the number density $\bar{n}$ over all three regions of the star. We recall that the equilibrium energy $e_{\mathrm{eq}}$ per nucleon that we have shown in the figures and tables always has the neutron mass subtracted out: see Eqs. (3), (13) and (29). Our fitting function consists of a constant term and three parts corresponding respectively to low, moderate, and high densities, matched together by appropriate weight functions $w_{1,2}(\bar{n})$, thus

[TABLE]

where

[TABLE]

The requirement that 56Fe be fitted for the vanishingly small values of $\bar{n}$ in the outer crust fixes the constant term at $e_{\mathrm{gr}}=-9.1536$ MeV. The numerical coefficients $p_{i}$ are given in Table 20, with $e_{\mathrm{eq}}$ measured in MeV and $\bar{n}$ in fm*-3*. The fit is valid over the range $10^{-9}$ fm${}^{-3}\lesssim\bar{n}\lesssim 3$ fm*-3*, the typical error of Eq. (75) over this range being $\lesssim 0.5$% for all three functionals; the maximum errors ($<3$%) occur at the phase boundaries, because the dependence $e_{\mathrm{eq}}(\bar{n})$ is fitted by a continuous function across the discontinuities of the argument $\bar{n}$.

For some purposes it is more convenient to deal with the total mass-energy density $\rho$ than with the energy per nucleon $e_{\mathrm{eq}}$. The function $\rho(\bar{n})$ is then equally well fitted by the function (75) through the identity (4).

The analytic fit (75) can be also used for the calculation of the baryon number density $\bar{n}$ and energy $e_{\mathrm{eq}}$ as functions of mass density $\rho$. For this purpose, we employ an iterative procedure, using the secant method for the logarithms (see, for example, Section 9.2 of Press et al. 1997),

[TABLE]

For each $i$th approximation to the number density, $\bar{n}_{i}$, Eq. (75) provides the mass-energy density estimate $\rho_{i}=\bar{n}_{i}\,(e_{\mathrm{eq,i}}+M_{\mathrm{n}}c^{2})/c^{2}$, which is used to correct $\bar{n}$ at the next iteration. Two starting values, $n_{0}$ and $n_{1}$, are required. For the first of these we take $\bar{n}_{0}=\rho c^{2}/e_{0}$, which, with $e_{0}=e_{\mathrm{gr}}+M_{\mathrm{n}}c^{2}$ = 930.4118 MeV, is exact at vanishingly small $\rho$, while for the second we take $\bar{n}_{1}=\bar{n}_{0}\rho/\rho_{0}$. The procedure rapidly converges: a fractional accuracy better than $10^{-6}$ is reached in 1 – 2 iterations for crustal densities and in no more than 5 iterations for the core (within the stability and causality limits listed in Tables 16 and 18).

C.2 Pressure as a function of density

For all four functionals we fit the calculated pressure, like the energy per nucleon, to a single continuous analytic function of the number density $\bar{n}$ that covers the entire star. It is convenient to work through the mass-energy density $\rho$, expressing it in terms of $\bar{n}$ through Eqs. (75) and (4). Our fitting function for the pressure then has the same form as in Potekhin et al. (2013),

[TABLE]

where $\xi\equiv\log_{10}(\rho/\textrm{g cm}^{-3})$. Setting $K$ = 0 gives the pressure in units of dyn cm*-2*, setting $K$ = -33.2047 gives it in units of MeV$\cdot$fm*-3*, the units of the figures and tables of this paper. The parameters $p_{i}$ are given in Table 21. The typical fit error of $P$ is $\approx 1$% for $6\lesssim\xi\lesssim 16$. The maximum errors of $\lesssim 4$% are reached at phase boundaries, because the dependence $P(\rho)$ is fitted by a continuous function across the discontinuities of the argument $\rho$.

C.3 Particle numbers

C.3.1 Particle numbers in the inner crust

As we see from Fig. 12, the equilibrium values $Z_{\mathrm{eq}}$ of the number of protons in a WS cell are discontinuous at certain densities. At the values of $\bar{n}$ below the proton drip density, $Z_{\mathrm{eq}}$ is a constant integer between the discontinuities, and is presented by Tables 8 – 11. In the ETF domain beyond proton drip, however, $Z_{\mathrm{eq}}$ is a smooth continuous function. The values of this function at any given $\bar{n}$ are subject to uncertainty, because the minimum of energy as a function of $Z$ is very shallow (see Fig. 16). Within this uncertainty, $Z_{\mathrm{eq}}$ can be roughly expressed as

[TABLE]

with the parameters $p_{i}$ given in Table 22. The differences between this fit and the data are shown in the lower panel of Figs. 13, 14.

Jumps of $Z$ from one integer to another throughout the inner crust are accompanied by simultaneous jumps of the number of neutrons $N$ in such a way that the proton fraction $Y_{\mathrm{p}}=Z/(Z+N)$ is continuous (see Fig. 9). It can be fitted as

[TABLE]

where the parameters $p_{i}$ are given in Table 23, assuming that $\bar{n}$ is measured in fm*-3*. The differences between the fit and the data are shown in the lower panel of Fig. 9.

It is also of interest to parametrize the cluster component $Z_{\mathrm{cl}}$ of $Z_{\mathrm{eq}}$, as given by Eq. (18a). At all densities $Z_{\mathrm{cl}}$ can be approximated analytically by

[TABLE]

where the parameter $p_{3}$, which is negative, is given in the last column of Table 24, and

[TABLE]

In this last equation the parameters $p_{1}$ and $p_{2}$ are given in Table 24, while

[TABLE]

which runs from 0 to 1 in the inner crust, $\bar{n}_{\mathrm{nd}}$ representing the neutron-drip density and $n_{\mathrm{cc}}$ the density at the core-crust interface. Actually, $\tilde{Z}$ itself, as given by Eq. (82), is a good approximation to $Z_{\mathrm{f}}$ for $\tilde{Z}\ll Z_{\mathrm{eq}}$, but for higher values of $\tilde{Z}$ the correction (81) is needed to ensure that $Z_{\mathrm{cl}}$ does not become negative.

With the number of free protons being defined by Eq. (42), the number fraction of free protons relative to all nucleons is given by $Y_{\mathrm{pf}}=Z_{\mathrm{f}}/A=Y_{\mathrm{p}}Z_{\mathrm{f}}/Z_{\mathrm{eq}}$ and thus can be represented analytically using the analytic fits already made. It is shown in the upper panel of Fig. 17, and the differences between the fit and the data are shown in the lower panel.

The equilibrium values $N_{\mathrm{eq}}$ of the number of neutrons in a WS cell can be expressed by the identity

[TABLE]

Now $Z_{\mathrm{eq}}$ has been parametrized in the proton-drip region by Eq. (79), while at lower densities in the inner crust it takes well defined integer values. Moreover, $Y_{\mathrm{p}}$ is parametrized throughout the inner crust by Eq. (80). Thus Eq. (84) suffices to parametrize $N_{\mathrm{eq}}$ at all densities in the inner crust; it defines the curve showing $N_{\mathrm{eq}}$ in Fig. 18 and 19. The deviations between this parametrization and the data are shown in the lower panels of Figs. 18 and 19.

The fraction of neutrons in the crust is $Y_{\mathrm{n}}=1-Y_{\mathrm{p}}$. In the inner crust, $N_{f}$ neutrons per unit WS cell are considered as free, where $N_{f}$ is given by Eqs. (18b) and (43). The number fraction of the free neutrons relative to the total number of baryons can be approximated as

[TABLE]

where $x$ is the same as in Eq. (83), and the parameters $p_{i}$ are given in Table 25. The differences between the fit and the data are shown in the lower panel of Fig. 20.

With $Z_{\mathrm{eq}},Y_{\mathrm{p}}$ and $Y_{\mathrm{nf}}$ already parametrized, $N_{\mathrm{cl}}$ can now be parametrized through the identity

[TABLE]

it defines the curves showing $N_{\mathrm{cl}}$ in Figs. 18 and 19. The deviations between this parametrization and the data are shown in the lower panels of Figs. 18 and 19.

C.3.2 Particle numbers in the core

At $n>n_{\mathrm{cc}}$, the nuclear clusters disappear, so that all baryons become essentially free: $Y_{\mathrm{nf}}=Y_{\mathrm{n}}$, $Y_{\mathrm{pf}}=Y_{\mathrm{p}}$. When the density increases still further, the chemical potential of electrons continues to increase. Eventually it exceeds the muon rest energy $m_{\mu}c^{2}=105.6584$ MeV. Then free muons are at equilibrium with the electrons, whence their chemical potentials are the same,

[TABLE]

as follows from Eqs. (39a) and  (39b). Thus, neglecting exchange, it follows that the Fermi energies (B) of the electrons and muons are equal (with the rest masses included), whence for a given number density of electrons $n_{\mathrm{e}}$, the number density of muons is given by the relation

[TABLE]

where

[TABLE]

is the relativity factor at the Fermi surface of the electrons (muons). (Actually, the validity of Eq. (88) depends on the muon gas being degenerate, but since the temperature of the star is not identically zero it follows that at the threshold where the muons just start to appear this condition will not be satisfied. However, the density range over which we do not have complete degeneracy for the muons is insignificant.)

It follows from Eqs. (88) and (89) that

[TABLE]

provided that $1+x_{\mathrm{e}}^{2}>(m_{\mu}/m_{\mathrm{e}})^{2}$ (otherwise $n_{\mu}=0$). Therefore, it is sufficient to have a fit to $Y_{\mathrm{e}}=n_{\mathrm{e}}/n$ in order to evaluate $Y_{\mu}=n_{\mu}/n$ and $Y_{\mathrm{p}}=Y_{\mathrm{e}}+Y_{\mu}$. We represent $Y_{\mathrm{e}}$ by

[TABLE]

with parameters listed in Table 26 (for $n$ measured in fm*-3*). Fig. 23 shows electron and muon number fractions as functions of $n$ in the core of a neutron star and the differences between the fit (91) and the numerical data.

C.4 Chemical potentials of nucleons

In the outer crust, the nucleons are bound in the nuclei. In the inner crust and in the core, there are free nucleons, and their chemical potentials can be of interest for some astrophysical problems. In this section we present fits to these chemical potentials.

C.4.1 Chemical potential of nucleons in the inner crust

In the inner crust, the chemical potential of free neutrons without the rest mass, $\mu_{\mathrm{n}}^{*}=\mu_{\mathrm{n}}-M_{\mathrm{n}}c^{2}$, is fitted by

[TABLE]

[TABLE]

Here, $n_{\mathrm{nf}}=Y_{\mathrm{nf}}\,\bar{n}$ is the number density of free neutrons, $Y_{\mathrm{nf}}$ is represented by Eq. (85), $\mu_{\mathrm{id}}^{*}$ is the Fermi energy of free neutrons (without the rest energy) in the model of a gas of nonrelativistic fermions, and the denominator in Eq. (92) for $\mu_{\mathrm{n}}^{*}$ is a correction factor to this model. Parameters $p_{i}$ are given in Table 27, assuming that $\bar{n}$ is measured in fm*-3*. Deviations of this parametrization from the computed data are shown in the lower panel of Fig. 21.

The proton chemical potential is then parametrized by using the beta-equilibrium condition (8), with $\mu_{\mathrm{n}}$ given by the parametrization (92), while we approximate $\mu_{\mathrm{e}}$ by the Fermi energy $\mu_{\mathrm{e}}^{\mathrm{kin}}+m_{\mathrm{e}}c^{2}$, where $\mu_{\mathrm{e}}^{\mathrm{kin}}$ is given by Eq. (B) and Eq. (80) for $Y_{\mathrm{p}}$; the neglected terms in $\mu_{\mathrm{e}}$ have a negligible impact. Deviations of this parametrization from the computed data are shown in the lower panel of Fig. 22.

C.4.2 Chemical potential of nucleons in core

In the core, the model of free fermion gas ($\mu_{\mathrm{id}}^{*}$) would be an inadequate starting approximation, because of the great role of strong interactions. We parametrize the neutron chemical potential without the rest mass in the core by

[TABLE]

Parameters $p_{i}$ are given in Table 28, assuming that $n$ is measured in fm*-3* and $\mu_{\mathrm{n}}^{*}$ in MeV. Deviations of this parametrization from the computed data are shown in the lower panel of Fig. 26.

The proton chemical potential $\mu_{\mathrm{p}}$ is now parametrized in the core the same way as it was in the inner crust. Deviations of this parametrization from the computed data are shown in the lower panel of Fig. 27.

C.5 Geometrical parameters and sizes of nuclei in the inner crust

In applications one sometimes needs more detailed information on microscopic distributions of nucleons in the inner crust than is given simply by the numbers of free and bound nucleons considered in Section C.3.1. For example, cross sections of scattering of electrons on nuclei in the crust depend on the charge distribution in a nucleus (e.g., Kaminker et al. 1999; Gnedin, Yakovlev & Potekhin 2001). For the use in such applications, we construct analytic approximations for the geometrical parameters $C_{q}$ and $a_{q}$, which enter Eq. (15) and are tabulated as functions of $\bar{n}$. The general form of these approximations is

[TABLE]

where $y=C_{\mathrm{p}}$, $a_{\mathrm{p}}$, $C_{\mathrm{n}}$, or $a_{\mathrm{n}}$. Equation (95) is similar to Eq. (17) of Potekhin et al. (2013), but more complicated, mainly because of the last term, which appears since $Z_{\mathrm{eq}}$ varies with increasing density in the crust, unlike constant $Z_{\mathrm{eq}}=40$ in Potekhin et al. (2013). As in Potekhin et al. (2013), the same form of parametrization (95) also provides approximations to the dimensionless nuclear-size parameters $x_{\mathrm{p}}$ and $x_{\mathrm{n}}$, which enter the expressions for electrical and thermal conductivities (Gnedin et al., 2001) and are related to the root-mean-square radii of the proton and neutron distributions in the clusters, $R_{\mathrm{rms,p}}$ and $R_{\mathrm{rms,n}}$, via $x_{\mathrm{p,n}}=\sqrt{5/3}\,R_{\mathrm{rms,p,n}}/R$, where $R$ denotes the cell radius. The parameters $p_{i}$ of approximation (95) for the parameters $a_{q}$, $C_{q}$, and $x_{q}$ ($q=p,n$) are listed in Table 29 (for $\bar{n}$ measured in fm*-3*). Figure 31 demonstrates the behavior of $C_{q}$, $a_{q}$, and $R_{\mathrm{rms},q}$ as functions of density (the top panels) and the fractional errors of the corresponding fits (95) (middle and bottom panels). The errors of the fits are within (1 – 2)%.