© T. Carreau et al. 2020

1. Introduction

It is generally assumed that the composition of an isolated neutron star (NS) is of “cold catalyzed matter”, meaning that determined in the ground state at zero temperature (see, e.g., Haensel et al. 2007; Chamel & Haensel 2008; Blaschke & Chamel 2018). In this hypothesis, the crust of an NS is supposed to be made of pure layers, each consisting of a one-component Coulomb crystal. However NSs, being born from core-collapse supernova explosions, are initially hot, with temperatures exceeding 10¹⁰ K. At such temperatures, the crust of a (proto-)NS is expected to be made of a Coulomb liquid composed of different nuclear species in a charge-compensating electron background: see Oertel & Hempel (2017) for a review. As the NS crust cools down, it is generally supposed that this multicomponent plasma (MCP) remains in full thermodynamic equilibrium until the ground state is reached. However, it is unlikely for full equilibrium to be maintained, after crystallization occurs, until T = 0 K. Moreover, if the NS cools down rapidly enough, the composition of the crust could be frozen at finite temperature, T_f, above the crystallization temperature, T_m (see, e.g., Goriely et al. 2011). Therefore, a more realistic picture of the crust would be that of a multicomponent solid.

For the outer crust, the coexistence of different nuclear species is not expected to significantly impact the static properties of the crust. Indeed, because of the relatively low crystallization temperatures, T_m ≲ 10⁹ K, the most probable nucleus is very close to the ground-state one-component plasma (OCP) composition, and the contribution of other ions is typically very small. In the inner crust, the situation is a priori less obvious, and deviations from the ground-state composition may be larger, due to the higher crystallization temperature, 10⁹ ≲ T_m ≲ 10¹⁰ K (Haensel et al. 2007). Jones (1999, 2001) indeed suggested that thermal fluctuation of the charge and neutron numbers may be quite significant for mass densities of ρ_B ≳ 10¹³ g cm⁻³ near the crystallization temperature. The presence of amorphous and heterogeneous phases in the inner crust leads to a higher temperature-independent electrical resistivity and strong ohmic dissipation, and significant consequences on the magnetic field evolution were predicted by Jones (2004). More recently, Pons (2013) suggested that the increased resistivity due to the amorphous structure could reflect into observational timing properties of X-ray pulsars. More generally, the presence of impurities in the crust has notable effects on transport and magneto-rotational properties of the NS (see, e.g., Schmitt & Shternin 2018; Gourgouliatos & Esposito 2018 for recent reviews), which in turn affect the NS thermal evolution. For these reasons, although cooling simulations are usually carried out using the ground-state composition, the presence of various nuclear species is taken into account via an “impurity factor”, often taken as a free parameter adjusted on observational cooling data: see, for instance, Viganò (2013).

A microscopic calculation of the impurity parameter at the crystallization temperature for the outer crust of a non-accreting unmagnetized NS was recently performed in Fantina et al. (2020). In the latter work, the nuclear distributions of the multicomponent liquid plasma at the crystallization point was computed fully self-consistently, adapting a general formalism originally developed for the description of supernova cores (Gulminelli & Raduta 2015; Grams et al. 2018). The crystallization temperature was determined in the OCP approximation, using a microscopic nuclear mass model based on deformed Hartree–Fock–Bogoliubov calculations (HFB-24, Goriely et al. 2013). The study of Fantina et al. (2020), performed on the outer crust, has subsequently been extended in Carreau et al. (2020), who calculated the crystallization temperature and the associated composition in the inner crust using the compressible liquid-drop (CLD) model approach of Carreau et al. (2019), with parameters optimized on four different microscopic models: namely BSk22, BSk24, BSk25, and BSk26 (developed by Goriely et al. 2013). Shell effects, as calculated in Pearson et al. (2019) for the same functionals, were added to the CLD model. The use of such an approach instead of a fully microscopic one not only reduces the computational time, but more importantly allows us to quantitatively estimate the model dependence of the results. The outcomes of Carreau et al. (2020) suggest that, while shell effects are important at the lowest densities close to the outer crust, the highest source of uncertainties comes from the smooth part of the nuclear functional, specifically the surface tension at extreme isospin values.

In the present work, we employed the same CLD model with parameters optimized on the same functionals as in Carreau et al. (2020), but we extended it by including a nuclear distribution in an MCP approach at equilibrium similar to that of Fantina et al. (2020). This also allows us to calculate the impurity parameter in the inner crust self-consistently, thus complementing the results obtained in Fantina et al. (2020) for the outer crust.

The formalism is described in Sect. 2. The numerical results are presented in Sect. 3; specifically, the composition of the inner crust is discussed in Sect. 3.1, and the impurity parameter in Sect. 3.2. Finally, we conclude in Sect. 4.

2. Model of the inner crust

2.1. MCP in nuclear statistical equilibrium

To model a full statistical equilibrium of ions in the inner crust, we extended the formalism of Fantina et al. (2020), allowing for the presence of dripped neutrons, which are supposed to constitute a homogeneous gas. The possible contribution of a free proton gas is expected to be small at the temperatures we considered, and thus it was neglected. This working hypothesis is a-posteriori confirmed by the calculation of the proton fugacity, z_p = exp[(μ_p − m_pc²)/(k_BT)], μ_p (m_p) being the proton chemical potential (mass), c the speed of light, and k_B the Boltzmann constant, which never exceeds −20 MeV in the density and temperature domain studied in this paper.

The NS crust at a given depth in the star is supposed to contain different ion species with mass and charge numbers (A^(j), Z^(j)) associated with different Wigner-Seitz cells of volume V^(j), such that p_j is the frequency of occurrence or probability of the component (j), with ∑_jp_j = 1. Thermodynamic quantities are defined in terms of the ion densities of the different species , which are related to the probabilities p_j through , where the bracket notation ⟨⟩ indicates ensemble averages.

The different (A^(j), Z^(j)) configurations are associated with different baryonic densities , such that the total baryonic density is (see Eq. (14)). Conversely, they share the same total pressure P imposed by the hydrostatic equilibrium and the same background densities of electrons, , and of free neutrons, . We also suppose that charge neutrality is realized in each cell, meaning that the proton density is the same in each cell (i.e. ) and equal to the electron density n_e (i.e. n_e = n_p = Z^(j)/V^(j)).

The free energy density of the multicomponent system is defined as:

(1)

where the free energy per ion of the component (j) accounts for the contribution of the ion, the dripped neutrons, and the electrons:

(2)

including their mutual interactions¹. For future convenience, the nuclear interactions between the ion and the neutron gas, and the Coulomb interactions between the ion and the electrons, are all included in the term . Therefore, the free neutron and electron components are simply given by:

(3)

where ℱ_g(e) is the free energy density of a uniform neutron (electron) gas at density n_g (n_e). The explicit expression of these terms is discussed in Sect. 2.4. The ion contribution, , can be written as:

(4)

The first term in Eq. (4), , noting m_n (m_p) the neutron (proton) mass, is given by:

(5)

where is the internal nuclear free energy and is the Coulomb interaction contribution. The explicit expressions of these terms, as well as of the last term in Eq. (4), δF^(j), accounting for the interaction between the ion and the surrounding (neutron) gas, depend on the adopted model and are discussed in Sects. 2.3 and 2.4. Finally, since in this work we are only interested in temperatures higher or equal to the melting temperature, where the MCP is expected to be in the liquid phase, the “ideal” contribution, , accounts for the translational center-of-mass motion:

(6)

where is the spin degeneracy. For this, we take for nuclei whose ground-state angular momentum is unknown, and the de Broglie wavelength of component (j) is given by:

(7)

ℏ being the Planck-Dirac constant, and the effective mass of the ion M^*(j) is defined as

(8)

The probabilities p_j and the densities are calculated so as to maximize the thermodynamic potential in the canonical ensemble. Because of the chosen free energy decomposition, we can observe that the electron and free neutron part of the free energy density, ℱ_e and ℱ_g, do not depend on , for instance,

(9)

where

(10)

Therefore, the variation can be performed on the ion part only:

(11)

where the single-ion canonical potential is given by:

(12)

In Eq. (11), the variations d are not independent, because of the normalization of probabilities and the baryonic number and charge conservation laws:

(13)

(14)

(15)

The correction factor on the right hand side of Eq. (14) accounts for the excluded volume, meaning the gas cannot occupy the nucleus volume. In the same equation, n_B is the total baryonic density, and is the average density of the ion (j). This latter can be calculated by imposing equilibrium with the nucleon gas via:

(16)

where is the pressure of the neutron gas. This expression is explicitly demonstrated in Sect. 2.4 (see Eq. (38)).

The constraints Eqs. (13)–(15) are taken into account by introducing Lagrange multipliers (α, μ_n, μ_p) leading to the following equations for the equilibrium densities :

(17)

with . Considering independent variations, the equilibrium distributions are given by

(18)

with the normalization

(19)

The single-ion grand-canonical potential reads:

(20)

where μ_n and μ_p can be identified with the neutron and proton chemical potentials, respectively. In the definitions above, the ion free energy contains the rest-mass energy, thus the chemical potentials include the rest-mass energies as well.

The calculation of the grand-canonical potential, , requires the evaluation of the chemical potentials μ_n, μ_p, as well as of the rearrangement term (last term in Eq. (12)):

(21)

These terms are worked out in Sects. 2.2 and 2.5, respectively.

Once the abundances of the different ions are calculated via Eq. (18) at the crystallization temperature, it is also possible to calculate the impurity parameter of the solid crust, which represents the variance of the ionic charge distributions and is defined as (see, e.g., the discussion in Sect. 7 in Meisel et al. 2018 for a review)

(22)

where p(Z^(j)) is the normalized probability distribution (integrated over all N^(j)) of the element Z^(j).

2.2. Evaluation of the chemical potentials

In a given thermodynamic condition expressed by a temperature T and a pressure P, the proton and neutron chemical potentials can be determined using the thermodynamic relation ℱ + P = μ_nn_n + μ_pn_p + μ_en_e, giving, together with the beta-equilibrium condition μ_n = μ_e + μ_p (μ_e being the electron chemical potential),

(23)

where the baryon and proton densities, n_B and n_p, are given by Eqs. (14) and (15), respectively, the free energy density ℱ is given by Eq. (1), and ℱ_e () is the free energy density (pressure) of the electron gas at density n_p = n_e. With this prescription, the equilibrium probabilities can only be determined by the solution of a complex nonlinear system of coupled equations which is a challenging numerical task. Such a complete nuclear statistical equilibrium formalism has been adopted by different authors (see Oertel & Hempel 2017; Burgio & Fantina 2018 for a review); however, simplified nuclear functionals were adopted, the density (instead of the pressure) was imposed, and the rearrangement term was neglected.

In the outer crust regime, it was found by Fantina et al. (2020) that a perturbative implementation of the nuclear statistical equilibrium as proposed by Grams et al. (2018) leads to a very fast convergence, with reduced computational cost and increased numerical precision. We therefore adopt this same prescription in the inner crust. In the perturbative treatment, the equilibrium problem is solved in the OCP approximation, as detailed in Sect. 2.3 below. This gives a first guess for the chemical potentials as:

(24)

where ℱ^OCP is the equilibrium free energy density in the OCP approximation, and , are the baryon and proton densities that, in the OCP approximation, lead to the pressure P (see Sect. 2.3). Similarly, the electron quantities ℱ_e and P_e are calculated at . With this guess, the ion abundances are readily calculated via Eq. (18), and again using Eq. (23) we can get an improved estimation of the chemical potentials as

(25)

(26)

where y_p = ⟨Z⟩/(n_B⟨V⟩) is the average proton fraction of the mixture, with ⟨Z⟩ = ∑_jp_jZ^(j). The problem can thus be solved by iteration. It turns out that the difference between the initial guess, Eq. (24), and the result of the first iteration, Eqs. (25)–(26), is so small for all pressures and temperatures considered in this work, that the simple OCP estimation, Eq. (24), can be kept.

The full MCP calculation becomes therefore computationally equivalent to the much simpler OCP one, with the additional advantage that the MCP results can be compared to the more standard OCP ones with no extra computational cost.

2.3. The OCP approximation

In the OCP approximation, the equilibrium configuration of inhomogeneous dense matter in the inner crust in full thermodynamic equilibrium is obtained by minimizing the free energy density in a Wigner-Seitz cell of volume V with the constraint of a given baryon density, n_B. (see Lattimer & Swesty 1991; Gulminelli & Raduta 2015; Carreau et al. 2020).

Similarly to the general MCP case of Sect. 2.1, we write:

(27)

where the variational variables are the mass number A and isospin ratio I = 1−2Z/A of the ion, its internal density n₀, the proton density in the cell n_p = n_e, and the density of the homogeneous gas of dripped neutrons n_g. As in Sect. 2.1, see Eq. (4), we include the interactions of the nucleus with the neutrons and electrons in the term F_i:

(28)

with

(29)

In the OCP approximation, the translational motion is limited to the single Wigner-Seitz cell (n_N = 1/V):

(30)

where the de Broglie wavelength λ is given by the same expression as in Eq. (7), with M^⋆(j) = M^⋆. The interacting part of the ion free energy can be decomposed as:

(31)

Analytical formulae were derived by Potekhin & Chabrier (2000) for these two terms; see their Eqs. (16) and (19), respectively. For this study, only the first term is included; indeed, the polarization correction is found to have no effect in the density and temperature regime studied in the present paper and is therefore neglected. In addition, the nuclear finite-size correction is also included. The latter is derived from the Gauss theorem and reads:

(32)

with r₀ = (4πn₀/3)^−1/3 and e as the elementary charge.

Finally, the interaction between the ion and the surrounding neutron gas is handled in the excluded volume approximation:

(33)

The equilibrium configuration is obtained by minimizing Eq. (27) with respect to the variational variables using the baryon density constraint limited to a single cell,

(34)

This leads to the following system of coupled differential equations²:

(35)

(36)

(37)

(38)

where the gas pressure is given by , and the baryon chemical potential μ_B results in:

(39)

In our parameterization (see Sect. 2.4), the in-medium modification of the nuclear energy arising from the external gas is governed by a single parameter, p, which does not depend on the external neutron density, but only on the isospin asymmetry, I. Therefore, and the baryon chemical potential can be identified with the chemical potential of the gas: μ_B = μ_g ≡ dℱ_g/dn_g.

At each value of the baryon density and temperature T ≥ T_m above the crystallization point, the system of coupled differential equations, Eqs. (35)–(38), is numerically solved as in Carreau et al. (2019). This procedure leads to the determination of the favored liquid composition (A, I, n₀, n_p, n_g)|_OCP and to the evaluation of the total free energy and pressure, ℱ^OCP and P, as well as of the electron component, , . These quantities allow one to compute the chemical potentials of the MCP using Eq. (24).

2.4. The free energy functional

The free energy functional for an isolated nucleus in the vacuum is modeled using the CLD model of Carreau et al. (2020), which we briefly outline here.

The nuclear free energy at temperature T of a nucleus of mass number A, isospin asymmetry I, and average density n₀, is decomposed into a bulk, surface, and Coulomb part as:

(40)

where f_b(n_B, δ, T) represents the free energy per baryon of bulk nuclear matter, with n_B = n_p + n_n, δ = (n_n − n_p)/n_B, and n_p(n_n) is the homogeneous proton (neutron) density. Assuming spherical nuclei, we write the Coulomb energy as:

(41)

with the surface and curvature free energies as in Newton et al. (2013), Lattimer & Swesty (1991):

(42)

with α = 5.5, and an isospin-dependent surface tension given by:

(43)

The surface and curvature parameters σ₀, b_s, p, σ_0, c, and β are optimized on extended Thomas-Fermi (ETF) mass tables built with the same nuclear functional adopted for the bulk term, see Carreau et al. (2020) for details. The same functional is also used to compute the free energy density of the neutron gas,

(44)

In principle, a shell and pairing correction should be added to the nuclear free energy expression, Eq. (40). However, it was shown in Carreau et al. (2020) that these corrections rapidly fade away with the temperature, and we thus consider that they can be neglected in the temperature range we explore in this work.

Concerning the nuclear models, we use the same functionals as in Carreau et al. (2020), namely the recent functionals of the BSk family BSk22, BSk24, BSk25, and BSk26 introduced by Goriely et al. (2013). These realistic microscopic models span a relatively large range in the symmetry energy parameters consistent with existing experimental constraints, thus covering the most important part of the present equation-of-state uncertainty (Pearson et al. 2014, 2018), meaning that the spread of the predictions of those models can be taken as a reasonable estimation of the model dependence of our results.

Finally, the free energy density ℱ_e and pressure P_e of the electron gas are calculated within a relativistic Sommerfeld expansion. The complete expressions can be found in Haensel et al. (2007): see their Eqs. (2.65) and (2.67), respectively. Exchange and correlation contributions are found to be very small in the ranges of density and temperature explored in this work and can be safely neglected.

2.5. Evaluation of the rearrangement term

The computation of the equilibrium distributions, Eq. (18), associated with a thermodynamic condition characterized by a temperature T, and chemical potentials μ_n, μ_p, requires the evaluation of the rearrangement term entering Eq. (12):

(45)

As already discussed in Fantina et al. (2020), the rearrangement term arises from the self-consistency induced by the Coulomb part of the ion free energy. This stems from the fact that, due to the strong incompressibility of the electrons, we have imposed charge conservation at the level of each cell:

(46)

This is at variance with the baryonic density, which can fluctuate from cell to cell (see Eq. (14)). As a consequence of that, any component of the free energy density that depends on the local cell proton density leads to a dependence on the local density through Eq. (46). Within the functional described in Sect. 2.4, this is only the case for the Coulomb interaction . The rearrangement term of component (j) thus reduces to:

(47)

where we used Eq. (46), implying the relation . Following Grams et al. (2018), to avoid the complication of a self-consistent resolution of Eq. (18), we looked for an approximation of Eq. (47) using the requirement that the most probable ion in the MCP mixture should coincide with the OCP result, if nonlinear mixing terms in the MCP are omitted. This condition is a direct consequence of the principle of ensemble equivalence in the thermodynamic limit (see Gulminelli & Raduta 2015).

To look for the extremum of Eq. (18), one has to consider that in the MCP, both n_p and n_g are imposed once the thermodynamic condition is specified. Therefore, these densities no longer act as constraints and should not be varied. The variation of Eq. (18) with respect to the ion variables A, I, n₀ thus gives:

(48)

(49)

(50)

where the partial derivatives are calculated at the values corresponding to the equilibrium OCP solution, and is given by Eq. (29) using Eq. (30), meaning by the OCP functional, that is nonlinear mixing terms are excluded.

By comparing Eqs. (48)–(50) to the OCP ones, Eqs. (35)–(38), and using P_g = n_gμ_n − ℱ_g, we can deduce that ℛ^(j) should not depend on , that is ℛ^(j) = ℛ^(j)(A^(j), I^(j)), and that at the OCP solution, we should have

(51)

This is satisfied if ℛ^(j) linearly depends on Z^(j) = A^(j)(1 − I^(j))/2. Our final expression for the rearrangement term is therefore:

(52)

where the quantity in the parentheses is calculated at the OCP solution.

3. Numerical results

We computed the finite-temperature composition of the inner crust of non-accreting unmagnetized NSs within our MCP approach, thus including a distribution of nuclei in nuclear statistical equilibrium. For the considered BSk functionals, the recent calculations of Pearson et al. (2020) show that nonspherical pasta structures are expected to be present at the highest densities above n_B ≈ 0.05 fm⁻³, close to the crust-core transition point. Since we only consider spherical nuclei in the present study, we limited our calculation to the density domain extending from the neutron-drip point to n_B = 0.04 fm⁻³.

All the results presented in this section were obtained using the BSk CLD models, with the surface and curvature parameters fitted to the corresponding ETF calculations and crust-core transition densities (see Table 1 of Carreau et al. 2020 for the explicit parameter values). In Sect. 3.1, the results for the inner-crust composition at a finite temperature are shown for the BSk24 CLD model, as an illustrative example, while the impurity parameter is presented in Sect. 3.2 for all four of the considered CLD models based on the BSk22, BSk24, BSk25, and BSk26 functionals.

Our calculations of the liquid MCP were performed at the crystallization temperature T_m and, for comparison, at 10¹⁰ K = T >  T_m. The reason of this choice stems from the fact that, depending on the NS cooling timescales, the composition may be already frozen at some temperature T_f >  T_m (see e.g., Goriely et al. 2011). In Carreau et al. (2020), the crystallization temperature of the inner crust was estimated to lie between ≈2.5 × 10⁹ K and ≈8 × 10⁹ K for the considered CLD models (see their Figs. 5 and 7, panel (a)). Therefore, we chose T_f = 10¹⁰ K as an illustrative example. Indeed, a more realistic estimate of T_f would require dynamical simulations, which are beyond the scope of this paper.

3.1. Equilibrium composition of the MCP

The average and most probable mass and charge number in the MCP are displayed in Fig. 1 as a function of the baryon density in the inner crust for T = 10¹⁰ K (panel a) and T = T_m (panel b). For comparison, the results obtained in the OCP approximation are also shown (dotted lines). We can see that the average and most probable values in the MCP approach follow the OCP ones very closely. This means that the deviations from the linear mixing rule in the liquid phase are small, as already noted in Fantina et al. (2020) for the outer crust. While the mass numbers increase with density, the charge number is almost constant, Z ≈ 40. The latter value is very close to that obtained at zero temperature (see also, the dotted curve in Fig. 6, panel (b), in Carreau et al. 2020 and Fig. 12 in Pearson et al. 2018), suggesting that the presence of Z ≈ 40 ions in the inner crust is a robust result.

Fig. 1. Variation with baryon density nB of the average (solid lines) and most probable (dashed lines) values of the charge number Z (blue lines), cluster mass number A (orange lines), and total mass number Acell (red lines) in the inner crust at two selected temperatures: T = 1010 K (panel a), and T = Tm (panel b). Results obtained in the one-component plasma (OCP) approximation are also shown (dotted lines).

To evaluate the width of the distribution, in Fig. 2 we show the normalized probability distribution p(Z) for T = 10¹⁰ K and T = T_m and for two selected densities in the inner crust: n_B = 5 × 10⁻⁴ fm⁻³ (panel a) and n_B = 10⁻² fm⁻³ (panel b). The peaks of the distributions, meaning the most probable Z, coincide with the charge numbers predicted in the OCP approximation (shown by the associated arrows), thus indicating that the linear mixing rule is a good approximation. To assess the importance of the rearrangement term, Eq. (47), we use vertical lines to mark the average values of the charge number, ⟨Z⟩, obtained when this term is not included in the calculations. We observe that the effect of the rearrangement term is significant, particularly at a higher density. Without taking this term into account, the distribution is systematically and considerably shifted toward a lower Z, proving that the rearrangement term is actually needed to satisfy the thermodynamic consistency. We can also notice that, as expected, the distributions become broader with increasing temperature and density, thus making the OCP approximation less reliable. The flattening of the distribution is more clearly visible in Fig. 3, where the normalized probability distribution p(Z) at the crystallization temperature is plotted for different increasing baryon densities in the inner crust. While the average value of Z is centered around 40 throughout the inner crust, the range of Z of the distribution varies from ≈20 closer to the neutron drip up to ≈40 near the crust-core transition.

Fig. 2. Normalized probability distribution p(Z) for nB = 5 × 10−4 fm−3 (panel a) and nB = 10−2 fm−3 (panel b) at two selected temperatures: T = Tm (orange squares), and T = 1010 K (blue circles). Arrows indicate the OCP solutions. Vertical dashed lines correspond to the value of ⟨Z⟩ obtained without considering the rearrangement term (see text for details).

Fig. 3. Normalized probability distribution p(Z) with increasing baryon density nB in the inner crust at the crystallization temperature Tm.

To better assess the evolution of the nuclear distribution, both in charge and mass number with density and temperature, in Fig. 4 we show the normalized probability distribution p(Z, N) for two selected densities in the inner crust: n_B = 5 × 10⁻⁴ fm⁻³ (panels (a) and (b)) and n_B = 10⁻² fm⁻³ (panels (c) and (d)), both at T = 10¹⁰ K (panels b and d) and T = T_m (panels (a) and (c)). As expected, going from lower to higher densities (upper to lower panels), we observe that the ion species become more neutron rich and that the distribution, both in Z and N, broadens when going from lower to higher temperatures (left to right panels).

Fig. 4. Normalized probability distribution of nuclei p(N, Z) for four chosen thermodynamic conditions. Panel a: nB = 5 × 10−4 fm−3, T = Tm; panel b: nB = 5 × 10−4 fm−3, T = 1010 K; panel c: nB = 10−2 fm−3, T = Tm; panel d: nB = 10−2 fm−3, T = 1010 K. In each panel, the OCP solution coincides with the intersection of the black lines.

3.2. Impurity parameter

The impurity parameter at the crystallization temperature is shown in the whole crust in Fig. 5, as calculated with the BSk24 CLD model (solid line). The black dot marks the neutron-drip point. We can see that the impurity parameter in the inner crust is higher than in the outer crust, meaning that the distribution is less peaked, thus the OCP approximation is less reliable than in the outer crust. Indeed, larger values of Q_imp indicate more appreciable deviations from the OCP predictions. For comparison, we also plot the impurity parameter, taken from Fantina et al. (2020), calculated in the outer crust with the HFB-24 model (dashed line). The latter calculations show more prominent variations of Q_imp with respect to the CLD model calculations. This is due to the natural inclusion of shell effects in the fully microscopic calculations, which exhibit bimodal distributions around values of pressure corresponding to the simultaneous presence of the two characteristic elements of adjacent layers (see Fig. 6 in Fantina et al. 2020). These strong fluctuations are naturally smoothed out in the CLD model, because the nuclear functional varies continuously with A and Z, and so does the probability. However, we can observe that the CLD calculation nicely interpolates the microscopic results, with an average impurity factor steadily increasing with the density and lying in the Q_imp ≈ 0, 1−2 interval. In the inner crust, neutrons drip out of the finite ion volume, and the associated shell effects naturally disappear (Chamel 2006). The inclusion of proton shell effects in the inner crust would require a formidable numerical investment, which is far beyond the scope of the paper. Moreover, it was suggested in Carreau et al. (2020) that these effects are small at the higher melting temperature of the inner crust, and that their effect on the observables is smaller than the uncertainty brought by our imperfect knowledge of the smooth part of the energy functional. For this reason, shell effects were completely neglected in our study. We expect that a fully microscopic calculation would still present oscillations in Q_imp beyond the drip point, that these oscillations should be progressively damped going deeper in the star, and that our calculation can be taken as a smooth interpolation of that oscillating behavior.

Fig. 5. Variation with baryon density nB of the impurity parameter Qimp in the crust at the crystallization temperature T = Tm. In the outer crust regime, the solid (dashed) line represents the BSk24 CLD (HFB-24) prediction. In the inner crust regime, the impurity parameter is calculated in the CLD approximation. Points indicate the neutron-drip transition.

To have a quantitative prediction of the impurity factor, the problem of model dependence has to be addressed. Apart from the modeling of finite temperature shell effects discussed above, the main source of uncertainty of the calculation comes from the choice of the nuclear functional. We show, in Fig. 6, the impurity parameter (Eq. (22)) as a function of the baryon density in the inner crust, at the crystallization temperature T_m (solid line), for the four considered BSk CLD models. These data are available in tabular format at the CDS. Considering that the chosen models are believed to cover the main uncertainty on the nuclear equation of state at subsaturation density (Pearson et al. 2018), we can take the spread of Q_imp values obtained by the four calculations as a reasonable estimation of the uncertainty on the impurity parameter. Since this latter represents the variance of the charge distribution, low values of Q_imp indicate that the distribution is quite peaked, and thus the OCP approach is a good approximation, which can also be seen from Figs. 1 and 2, panel (a). The monotonic increase of the impurity parameter with density is also in accordance with Fig. 3, which clearly shows the growth of the width of the charge distribution with increasing density. While at lower densities all the models predict similar values of the impurity parameter (≲5), at higher densities the spread among the models becomes larger, and the model associated with the lower (larger) symmetry energy coefficient at saturation has the larger (lower) Q_imp. The same trend is observed at T = 10¹⁰ K (dashed lines in Fig. 6), although the hierarchy of the models is not preserved.

Fig. 6. Variation with baryon density nB of the impurity parameter Qimp in the inner-crust regime at two selected temperatures: T = 1010 K, and T = Tm, based on BSk22 (red lines), BSk24 (black lines), BSk25 (blue lines), and BSk26 (green lines) CLD calculations.

4. Conclusions

In this paper, we present a multicomponent approach to the modeling of the crust of isolated, unmagnetized NSs. Completing the work of Fantina et al. (2020) on the outer crust, the composition of the inner crust is evaluated with an extended nuclear statistical equilibrium based on a CLD model for the nuclear part of the ion energetics. To achieve thermodynamical consistency, a rearrangement term is explicitly worked out. This term has an important effect on the distributions and it is necessary to recover the correct limit at zero temperature. Since NSs are born hot, the equilibrium composition of a mature NS can be determined assuming a liquid phase for the MCP, at the lowest temperature at which strong and weak equilibrium are attained. In the absence of a dynamical estimation of the associated reaction rates, we consider the lowest temperature limit as given by the OCP crystallization temperature T_m. We show that at that temperature, the OCP approximation gives a very good estimation of the average composition of the inner crust, with nonlinear mixing terms playing a very small role in the liquid phase. However, an important contribution of impurities is obtained, favoring the picture of a temperature-independent high resistivity in the inner crust for all T <  T_m.

In order to reach quantitative predictions for the associated impurity parameter, we considered four different realistic microscopic nuclear functionals of the BSk family, which cover the present uncertainty in the nuclear modeling below the saturation density of nuclear matter. The impurity parameter is seen to increase with the density, and values in the Q_imp ≈ 20−40 interval are reached at the highest densities considered in this study, namely n_B = 0.04 fm⁻³. Higher values of the impurity parameter might be expected in the deepest region of the inner crust, close to the core-crust transition, due to the presence of nonspherical pasta phases, which have not been considered in the present study.

Acknowledgments

This work has been partially supported by the IN2P3 Master Project MAC and the CNRS PICS07889. Discussions with N. Chamel are gratefully acknowledged.

References

All Figures

Fig. 1. Variation with baryon density nB of the average (solid lines) and most probable (dashed lines) values of the charge number Z (blue lines), cluster mass number A (orange lines), and total mass number Acell (red lines) in the inner crust at two selected temperatures: T = 1010 K (panel a), and T = Tm (panel b). Results obtained in the one-component plasma (OCP) approximation are also shown (dotted lines).

In the text

	Fig. 2. Normalized probability distribution p(Z) for nB = 5 × 10−4 fm−3 (panel a) and nB = 10−2 fm−3 (panel b) at two selected temperatures: T = Tm (orange squares), and T = 1010 K (blue circles). Arrows indicate the OCP solutions. Vertical dashed lines correspond to the value of ⟨Z⟩ obtained without considering the rearrangement term (see text for details).
In the text

	Fig. 3. Normalized probability distribution p(Z) with increasing baryon density nB in the inner crust at the crystallization temperature Tm.
In the text

	Fig. 4. Normalized probability distribution of nuclei p(N, Z) for four chosen thermodynamic conditions. Panel a: nB = 5 × 10−4 fm−3, T = Tm; panel b: nB = 5 × 10−4 fm−3, T = 1010 K; panel c: nB = 10−2 fm−3, T = Tm; panel d: nB = 10−2 fm−3, T = 1010 K. In each panel, the OCP solution coincides with the intersection of the black lines.
In the text

	Fig. 5. Variation with baryon density nB of the impurity parameter Qimp in the crust at the crystallization temperature T = Tm. In the outer crust regime, the solid (dashed) line represents the BSk24 CLD (HFB-24) prediction. In the inner crust regime, the impurity parameter is calculated in the CLD approximation. Points indicate the neutron-drip transition.
In the text

	Fig. 6. Variation with baryon density nB of the impurity parameter Qimp in the inner-crust regime at two selected temperatures: T = 1010 K, and T = Tm, based on BSk22 (red lines), BSk24 (black lines), BSk25 (blue lines), and BSk26 (green lines) CLD calculations.
In the text