Inner crust of a neutron star at crystallization in a multi-component approach

The possible presence of amorphous and heterogeneous phases in the inner crust of a neutron star is expected to reduce the electrical conductivity of the crust, with potentially important consequences on the magneto-thermal evolution of the star. In cooling simulations, the disorder is quantified by an impurity parameter which is often taken as a free parameter. We aim to give a quantitative prediction of the impurity parameter as a function of the density in the crust,performing microscopic calculations including up-to-date microphysics of the crust. A multi-component approach is developed at finite temperature using a compressible liquid drop description of the ions with an improved energy functional based on recent microscopic nuclear models and optimized on extended Thomas-Fermi calculations. Thermodynamic consistency is ensured by adding a rearrangement term and deviations from the linear mixing rule are included in the liquid phase. The impurity parameter is consistently calculated at the crystallization temperature as determined in the one-component plasma approximation for the different functionals. Our calculations show that at the crystallization temperature the composition of the inner crust is dominated by nuclei with charge number around $Z \approx 40$, while the range of the $Z$ distribution varies from about 20 near the neutron drip to about 40 closer to the crust-core transition. This reflects on the behavior of the impurity parameter that monotonically increases with density up to around 40 in the deeper regions of the inner crust. Our study shows that the contribution of impurities is non-negligible, thus potentially having an impact on the transport properties in the neutron-star crust. The obtained values of the impurity parameter represent a lower limit; larger values are expected in the presence of non-spherical geometries and/or fast cooling dynamics.


Introduction
It is generally assumed that the composition of an isolated neutron star (NS) is that of "cold catalyzed matter", i.e. 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 a NS is supposed to be made of pure layers, each consisting of a one-component Coulomb crystal (except possibly at the interface between different layers where multinary compounds may exist; see Chamel & Fantina (2016) for a discussion). However NSs, being born from core-collapse supernova explosions, are initially hot, with temperatures exceeding 10 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 compen- The tables of the impurity parameter at the crystallization temperature shown in Fig. 6, as well as tables of the impurity parameter for different values of density and temperature relevant for the inner crust, are available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsweb.ustrasbg.fr/cgi-bin/qcat?J/A+A/ sating electron background, see Oertel et al. (2017) for a review. As the NS crust cools down, it is generally supposed that this multi-component plasma (MCP) remains in full thermodynamic equilibrium until the ground state is reached. However, it is unlikely that full equilibrium is 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 some 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 multi-component solid.
For the outer crust, the co-existence 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 9 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 larger crystallization temperature, 10 9 T m 10 10 K (Haensel et al. 2007). Indeed, Jones (1999Jones ( , 2001 suggested that thermal fluctu-Article number, page 1 of 10 arXiv:2005.02089v1 [astro-ph.HE] 5 May 2020 ation of the charge and neutron numbers may be quite significant for mass densities ρ B 10 13 g cm −3 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 important consequences on the magnetic field evolution were predicted by Jones (2004). More recently, Pons et al. (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ò et al. (2013).
A microscopic calculation of the impurity parameter at the crystallization temperature for the outer crust of a nonaccreting unmagnetized NS has been recently performed in Fantina et al. (2020). In the latter work, the nuclear distributions of the multi-component liquid plasma at crystallization has been 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 been subsequently 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 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 employ the same CLD model with parameters optimized on the same functionals as in Carreau et al. (2020), but we extend it by including a nuclear distribution in a 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 selfconsistently, 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.

MCP in nuclear statistical equilibrium
To model a full statistical equilibrium of ions in the inner crust, we extend 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 it will be neglected. This working hypothesis is a-posteriori confirmed by the calculation of the proton fugacity, z p = exp[(µ p − m p c 2 )/(k B T )], µ 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 number (A (j) , Z (j) ), associated to different Wigner-Seitz cells of volume V (j) , such that p j is the frequency of occurrence or probability of the component (j), with j p j = 1. Thermodynamic quantities are defined in terms of the ion densities of the different species n (j) N , which are related to the probabilities p j through n (j) N = p j / V , where the bracket notation indicates ensemble averages.
The different (A (j) , Z (j) ) configurations are associated with different baryonic densities n (j) B , such that the total baryonic density is n B = j p j n (j) B (see Eq. (14)). Conversely, they share the same total pressure P imposed by the hydrostatic equilibrium and the same background densities of electrons, n (j) e = n e , and of free neutrons, n (j) g = n g . We also suppose that charge neutrality is realized in each cell, meaning that the proton density is the same in each cell, i.e. n (j) p = n p , and equal to the electron density n e , i.e. n e = n p = Z (j) /V (j) .
The free energy density of the multi-component system is defined as: where the free energy per ion of the component (j) accounts for the contribution of the ion, the dripped neutrons and the electrons, including their mutual interactions 1 . 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 F (j) i . Therefore, the free neutron and electron components are simply given by: where F 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, F i , can be written as: The first term in Eq. (4), F (j),0 i , noting m n (m p ) the neutron (proton) mass, is given by where F (j),nuc i is the internal nuclear free energy and F (j),int i 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 will be 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, F (j),id i , accounts for the translational centerof-mass motion: where g (j) s is the spin degeneracy, which we take g (j) s = 1 for nuclei whose ground-state angular momentum is unknown, and the de Broglie wavelength of component (j) is given by being the Planck-Dirac constant, and the effective mass of the ion M * (j) is defined as The probabilities p j and the densities n (j) N are calculated such 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, F e and F g , do not depend on n (j) where Therefore, the variation can be performed on the ion part only: where the single-ion canonical potential is given by: In Eq. (11), the variations dn (j) N are not independent because of the normalization of probabilities, and the baryonic number and charge conservation laws: The correction factor on the right hand side of Eq. (14) accounts for the excluded volume, i.e. the gas cannot occupy the nucleus volume. In the same equation, n B is the total baryonic density and n (j) 0 is the average density of the ion (j). This latter can be calculated by imposing equilibrium with the nucleon gas via: where P g = n 2 g d(F g /n g )/dn g 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 n (j) with . Considering independent variations, the equilibrium distributions are given by with the normalization The single-ion grand-canonical potentialΩ (j) N reads: 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,Ω i , requires the evaluation of the chemical potentials µ n , µ p , as well as of the rearrangement term (last term in Eq. (12)), These terms will be worked out in Sects. 2.2 and 2.5, respectively. Once the abundancies 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) where p(Z (j) ) is the normalized probability distribution (integrated over all N (j) ) of the element Z (j) .

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 F + P = µ n n n + µ p n p + µ e n e , giving, together with the beta equilibrium condition µ n = µ e + µ p (µ e being the electron chemical potential), where the baryon and proton densities, n B and n p , are given by Eqs. (14) and (15), respectively, the free energy density F is given by Eq.
(1), and F e (P e = n 2 e d(F e /n e )/dn 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 non-linear 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 et al. (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: where F OCP is the equilibrium energy density in the OCP approximation, and n OCP B , n OCP p are the baryon and proton densities that, in the OCP approximation, lead to the pressure P (see Sect. 2.3). Similarly, the electron quantities F e and P e are calculated at n e = n OCP p . With this guess, the ion abundancies are readily calculated via Eq. (18), and using again Eq. (23) we can get an improved estimation of the chemical potentials as where y p = Z /(n B V ) is the average proton fraction of the mixture, with Z = j p j Z (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.

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: where the variational variables are the mass number A and isospin ratio I = 1−2Z/A of the ion, its internal density n 0 , 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 : with In the OCP approximation, the translational motion is limited to the single Wigner-Seitz cell (n N = 1/V ): 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: Analytical formulae have been 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 with r 0 = (4πn 0 /3) −1/3 . Finally, the interaction between the ion and the surrounding neutron gas is treated in the excluded volume approximation: 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, This leads to the following system of coupled differential equations 2 : where the gas pressure is given by P g = n 2 g d(F g /n g )/dn g = n g µ B − F g , and the baryon chemical potential µ B results: In our parametrization (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, ∂F 0 i /∂n g = 0 and the baryon chemical potential can be identified with the chemical potential of the gas, µ B = µ g ≡ dF g /dn g .
At each value of the baryon density n OCP B 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 0 , n p , n g )| OCP and to the evaluation of the total free energy and pressure, F OCP and P , as well as of the electron component, F e (n OCP e ), P e (n OCP e ). These quantities allow one to compute the chemical potentials of the MCP using Eq. (24).

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), that we briefly recall.
The nuclear free energy F nuc i at temperature T of a nucleus of mass number A, isospin asymmetry I, and average 2 These equations are equivalent to Eqs. (8)-(11) in Carreau et al. (2020). Indeed, the notation Fi in Carreau et al. (2020) is equivalent to the notation F 0 i used in the present paper. We note however that there is a misprint in Eq. (9) in Carreau et al. (2020) (although the calculations have been done correctly); indeed, the term ∆mn,pc 2 should not appear in their Eq. (9). density n 0 , is decomposed into a bulk, surface and Coulomb part as: 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: with e the elementary charge, and the surface and curvature free energies as in Newton et al. (2013); Lattimer & Swesty (1991): with α = 5.5, and an isospin dependent surface tension given by: The surface and curvature parameters σ 0 , 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, 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 will 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 relative large range in the symmetry energy parameters consistent with existing experimental constraints thus covering the most important part of the present EoS uncertainty (Pearson et al. 2014(Pearson et al. , 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 F 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.

Evaluation of the rearrangement term
The computation of the equilibrium distributions, Eq. (18), associated to a thermodynamic condition characterized by a temperature T and chemical potentials µ n , µ p , requires the evaluation of the rearrangement term entering Eq. (12), 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: This is at variance with the baryonic density that 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 n (j) p = n p leads to a dependence on the local density n where we have used Eq. (46) implying the relation ∂n p /∂n (j) N = Z (j) . Following Grams et al. (2018), to avoid the complication of a self-consistent resolution of Eq. (18), we look 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 non-linear 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 do not act anymore as constraints and should not be varied. The variation of Eq. (18) with respect to the ion variables A, I, n 0 thus gives: where the partial derivatives are calculated at the values corresponding to the equilibrium OCP solution and F 0 i is given by Eq. (29) using Eq. (30), i.e. by the OCP functional, that is non-linear mixing terms are excluded.  to the OCP ones, Eqs. (35)-(38), and using P g = n g µ n − F g , we can deduce that R (j) should not depend on n (j) 0 , i.e. R (j) = R (j) (A (j) , I (j) ), and that at the OCP solution we should have: This is satisfied if R (j) linearly depends on Z (j) = A (j) (1 − I (j) )/2. Our final expression for the rearrangement term is therefore: where the quantity in the parenthesis is calculated at the OCP solution.

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 non-spherical pasta structures are expected to be present at the highest densities above n B ≈ 0.05 fm −3 , 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 −3 . All the results presented in this Section were obtained using the BSk CLD models, with the surface and curvature parameters fitted to the corresponding Extended Thomas-Fermi 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 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 the four considered CLD models based on the BSk22, BSk24, BSk25, and BSk26 functionals.
Our calculations of the liquid MCP are performed at the crystallization temperature T m and, for comparison, at 10 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 9 K and ≈ 8 × 10 9 K for the considered CLD models (see their Figs. 5 and 7, panel (a)). Therefore, we have chosen T f = 10 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.

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 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 very closely the OCP ones. This means that the deviations from the linear mixing rule in the liquid phase are small, as already noticed 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.
To evaluate the width of the distribution, we show in Fig. 2 the normalized probability distribution p(Z) for T = 10 10 K and T = T m and for two selected densities in the inner crust: n B = 5 × 10 −4 fm −3 (panel (a)) and n B = 10 −2 fm −3 (panel (b)). The peaks of the distributions, i.e. 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 mark with vertical lines 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 higher density. Indeed, without taking into account this term, the distribution is systematically and considerably shifted towards 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 in- ner 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.
To better assess the evolution of the nuclear distribution, both in charge and mass number with density and temperature, we show in Fig. 4 the normalized probability distribution p(Z, N ) for two selected densities in the inner crust: n B = 5 × 10 −4 fm −3 (panels (a) and (b)) and n B = 10 −2 fm −3 (panels (c) and (d)), both at T = 10 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 ions 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).

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, that 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 interval Q imp ≈ 0, 1 − 2. Going 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 much beyond the purpose 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.
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, as well as tables of the impurity parameter on a density-temperature grid for the four considered models, 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 sub-saturation density (Pearson et al. 2018), we can take the spread of Q imp values as 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 that the OCP approach is a good approximation, as it 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 show 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, with the model associated to the lower (larger) symmetry energy coefficient at saturation having the larger (lower) Q imp . The same trend is observed at T = 10 10 K (dashed lines in Fig. 6), although the hierarchy of the models is not preserved.

Conclusions
In this paper, we have presented a multi-component approach for 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, non-linear 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 consider four different realistic microscopic nuclear functionals of the BSk family, which cover the present uncertainty in the nuclear modeling below the nuclear matter saturation density. The impurity parameter is seen to increase with the density, and values in the interval Q imp ≈ 20 − 40 are reached at the highest densities considered in this study, namely n B = 0.04 fm −3 . 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 non-spherical pasta phases, which have not been considered in the present study.