The proto-neutron star inner crust in a multi-component plasma approach

Proto-neutron stars (PNS) are born hot, with temperatures exceeding a few times $10^{10}$ K. In these conditions, the PNS crust is expected to be made of a Coulomb liquid composed of an ensemble of different nuclear species. We perform a study of the beta-equilibrated PNS crust in the liquid phase in a self-consistent multi-component plasma (MCP) approach, thus allowing us to consistently calculate the impurity parameter, often taken as a free parameter in cooling simulations. We developed a self-consistent MCP approach at finite temperature using a compressible liquid-drop description of the ions, with surface parameters adjusted to reproduce experimental masses. The treatment of the ion centre-of-mass motion was included through a translational free-energy term accounting for in-medium effects. The results of self-consistent MCP calculations are systematically compared with those performed in a perturbative and in the one-component plasma treatment. We show that the inclusion of non-linear mixing terms arising from the ion centre-of-mass motion leads to a breakdown of the ensemble equivalence between the one-component and MCP approach. Our findings illustrate that the abundance of light nuclei becomes important, eventually dominating the distribution at higher density and temperature. This is reflected in the impurity parameter, which, in turn, may have a potential impact on NS cooling. For practical applications, we also provide a fitting formula for the impurity parameter in the PNS inner crust. Our results obtained within a self-consistent MCP approach show important differences in the prediction of the PNS composition with respect to those obtained with a one-component or a perturbative MCP approximation, particularly in the deeper region of the crust. This highlights the importance of a full, self-consistent MCP calculation for reliable predictions of the PNS crust composition.


Introduction
Formed from core-collapse supernova explosions, neutron stars (NSs) are initially born hot, with temperatures exceeding a few times 10 10 K ≈ 1 MeV (Prakash et al. 1997;Pons et al 1999;Haensel et al. 2007).At such temperatures, the crust of the proto-neutron star (PNS) is expected to be made of a Coulomb liquid composed of different nuclear species in a background of electrons and nucleons (see e.g.Oertel et al. (2017) for a review).It is generally assumed that, as the NS crust cools down, this multi-component plasma (MCP) remains in full thermodynamic equilibrium until the ground state at zero temperature is reached.Under this so-called 'cold catalysed matter' hypothesis, the NS crust is considered to be made of pure layers, each consisting of a one-component Coulomb crystal.However, if the NS cools down rapidly enough, the composition of the crust might be frozen at a temperature higher than the crystallisation temperature (see e.g.Goriely et al. (2011)), and the one-component picture becomes less reliable (see e.g.Fantina et al. (2020); ⋆ The tables of the impurity parameter shown in Fig. 15 are available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsweb.u-strasbg.fr/cgi-bin/qcat?J/A+A/ Carreau et al. (2020b)).Moreover, a correct description of the finite-temperature PNS crust requires the computation of a full nuclear ensemble.However, due to the complexity of calculating the nuclear distribution of MCP, most approaches to modelling the finite-temperature crust have been performed within the framework of the so-called one-component plasma (OCP) approximation.The latter assumes that, at a given thermodynamic condition, the ensemble of nuclei is represented by only one kind of nucleus, the one that is thermodynamically more favourable (see e.g.Onsi et al. (2008); Avancini et al. (2009Avancini et al. ( , 2017)); Fantina et al. (2020); Carreau et al. (2020a,b); Dinh Thi et al. (2023)).This approximation is justified at relatively low densities and temperatures, where the most probable nucleus is very close to the ground-state OCP composition and the contribution of other ions is typically very small, or when computing average thermodynamic quantities (Burrows & Lattimer 1984).On the other hand, the coexistence of different nuclear species can have an important impact on transport properties as well as on the magneto-rotational evolution of NSs (Jones 2004;Pons et al. 2013) (see also Schmitt & Shternin (2018); Gourgouliatos & Esposito (2018) for a review).For these reasons, (P)NS cooling simulations are usually performed using the ground-state OCP composition, but the presence of various nuclei is taken into ac-count via an 'impurity factor', often taken as a free parameter adjusted on observational cooling data (see e.g.Viganò et al. (2013)).
Self-consistent calculations of the impurity parameter for non-accreting unmagnetised NSs at relatively low temperatures, around the crystallisation point (k B T ≈ 0.01 − 1 MeV, k B being the Boltzmann constant), were recently presented by Fantina et al. (2020) and Carreau et al. (2020b).Specifically, Carreau et al. (2020b) investigated the MCP in the NS inner crust within a compressible liquid-drop (CLD) model with parameters optimised on microscopic energy-density functionals (Goriely et al. 2013).Carreau et al. (2020b) calculated the nuclear distributions using a perturbative implementation of the nuclear statistical equilibrium, as first proposed in Grams et al. (2018) for applications to supernova matter.In this perturbative treatment, for a given baryon density n B and temperature T , the beta-equilibrium composition is firstly found in the OCP approximation.This procedure yields the (OCP) neutron and proton chemical potentials, which are subsequently used in the computation of the ion abundances.Although not fully self-consistent, this approach has the advantage of leading to a very fast convergence, with reduced computational cost with respect to a full nuclear statistical equilibrium treatment (see e.g.Botvina & Mishustin (2010); Gulminelli & Raduta (2015); Raduta & Gulminelli (2019)).However, the applicability of such a perturbative approach requires a more thorough assessment at high densities, that is, in the deeper regions of the (P)NS crust, and at relatively high temperatures above the crystallisation temperature.
In the present work, we perform a study of the PNS crust in a self-consistent multi-component approach.To this aim, we extended and improved the previous work of Carreau et al. (2020b) in several aspects: (i) we extended the study of the (P)NS inner crust to higher temperatures above the crystallisation temperature, and to higher densities, until the crust-core transition; (ii) we improved the description of the ion centre-of-mass motion in the liquid crust by employing a more realistic prescription for the translational energy, where in-medium effects are accounted for (see the discussion in Dinh Thi et al. (2023) and references therein); and (iii) we went beyond the perturbative approach and performed a self-consistent calculation of the MCP in the whole PNS inner crust, which allowed us to calculate the impurity parameter in a self-consistent way.
The paper is organised as follows: The formalism is discussed in Sect.2: the MCP approach is introduced in Sect.2.1, the improved description of the translational free-energy term is presented in Sect.2.2, while the ion free energy in the OCP approximation is depicted in Sect.2.3, and the comparison between the MCP and OCP approaches is discussed in Sect.2.4.Numerical results are presented in Sect.3, where we describe both the implementation of the perturbative MCP approach in Sect.3.1 and the self-consistent MCP calculations in Sect.3.2; the outcomes regarding the impurity parameter as well as a fitting formula for its practical implementation in numerical simulations are presented in Sect.3.3.Finally, we draw conclusions in Sect. 4.

MCP in nuclear statistical equilibrium
To model the full statistical distribution of ions in the PNS liquid inner crust, we extended the formalism of Carreau et al. (2020b) for its applicability at higher densities and temperatures.As in Carreau et al. (2020b), we assume that beta-equilibrium holds, a condition that is met in late PNS cooling evolution stages (Prakash et al. 1997;Yakovlev & Pethick 2004;Page et al. 2006).For a given thermodynamic condition defined by the total baryonic number density n B and the temperature T , the PNS inner crust is supposed to be composed of different spherical Wigner-Seitz (WS) cells of volume V ( j)  WS , each containing an ion with mass and proton number (A ( j) , Z ( j) ) and internal density n ( j)  i = A ( j) /V ( j) N (V ( j) N = 4πr ( j)3 N /3 being the volume occupied by the cluster, with r ( j) N the cluster radius), such that u

WS
is the volume fraction occupied by the cluster in the cell j.Each ion is surrounded by a uniform gas of electrons and neutrons of density n e and n gn , respectively.In addition, we suppose that charge neutrality is realised in every cell, meaning that the proton density is the same in each cell and equal to the electron density, that is, n e = n p = Z ( j) /V ( j) WS .At finite temperature, a free proton gas could also be present and was included in Dinh Thi et al. (2023), where the PNS liquid crust was studied in the OCP approximation.However, we find that for the considered temperature regimes, k B T ≲ 2.0 MeV, the proton-gas density remains very small, namely about a few times 10 −3 fm −3 at most at the bottom of the crust, and its effects on the equation of state and on the crust composition are negligible.For these reasons, we ignore the presence of the free protons in the present work, as was done in Carreau et al. (2020b).Moreover, in this first step to improve our previous treatment (Carreau et al. 2020b), we neglect the possible presence of non-spherical configurations, which could exist in a narrow region at the bottom of the inner crust.
Denoting n ( j) N , the ion density of a species (A ( j) , Z ( j) ), the frequency of occurrence or probability of the component ( j) is given by n ( j) N = p j /⟨V WS ⟩, with j p j = 1 and the bracket notation ⟨⟩ indicates ensemble averages.The total free-energy density 1 of the system can therefore be written as where F g is the free-energy density of the uniform neutron gas, see Sect.2.3 below, and F e is the free-energy density of the uniform electron gas, for which we use the same expressions as in Carreau et al. (2020b); see Eq. (2.65) in Haensel et al. (2007).
The second term in Eq. ( 1), namely −V ( j) N F g , accounts for the excluded volume, that is, the impenetrability of the different ions, while the cluster free energy F ( j)  i is given by where j) m p is the bare mass of the cluster, with m p (m n ) being the proton (neutron) mass, c the speed of light, j) , T ) the cluster bulk free energy, and surf+curv is the sum of the Coulomb, surface, and curvature energies.Finally, the last term in Eq. (2), F ⋆,( j),MCP trans , represents the translational free energy of the cluster in the liquid phase, where λ ⋆,( j) i is the ion effective thermal wavelength, g ( j)  s is the spin degeneracy, and is the volume fraction available for the ion motion, ūf = ⟨V f ⟩/⟨V WS ⟩, and the average free volume is given by ⟨V f ⟩ = j p j V ( j) f , with where r ( j) WS is the WS cell radius.The translational term is discussed in more detail in Sect.2.2, while the bulk and interaction terms are addressed in Sect.2.3.
Following Gulminelli & Raduta (2015), Grams et al. (2018), Fantina et al. (2020), andCarreau et al. (2020b), the densities n ( j)  N were calculated so as to minimise the thermodynamic potential in the canonical ensemble.Because of the chosen free-energy decomposition in Eq. ( 1), the electron and neutron gas free-energy densities, F e and F g , do not depend on n ( j) N .Therefore, the variation with respect to the ion densities n ( j) N , dF MCP {n ( j)  N }/dn ( j) N , can be performed on the ion part only, that is, on the sum term in Eq. ( 1).However, the variations dn ( j)  N are not independent because of the baryon number conservation and charge neutrality, These constraints are taken into account by introducing two Lagrange multipliers, which can be shown to be directly connected to the chemical potentials, yielding the following equation for the equilibrium densities n ( j) N : In Eq. ( 8), µ e = dF e /dn e is the electron chemical potential and µ n is given by where one can recognise the uniform-matter chemical potential of the unbound neutrons, µ gn ≡ ∂F g /∂n gn , while the additional term ∆µ n accounts for the in-medium modification of the ion free energy arising from the relative motion of the ions with respect to the surrounding neutron fluid (see also Dinh Thi et al. (2023)).This affects the global chemical potential of the system, breaking the identity µ n = µ gn , which is assumed in all OCP and MCP equation-of-state models based on the cluster degrees of freedom we are aware of; see Oertel et al. (2017) for a review.The presence of this in-medium term ∆µ n also implies that the pressure equilibrium condition inside each WS cell is modified with respect to the standard results, that is, the equivalence between the cluster and gas pressure (see e.g.Lattimer et al. (1985)).In order to keep the same formal expression as in the literature -see Eq. ( 38) below-we include the in-medium modification in the definition of the gas pressure, which we define as P g ≡ µ n n gn − F g .
The last term in the first row of Eq. ( 8), , arises from the dependence on the ion density n ( j) N of both the translational and the interaction part of the ion free energy.Indeed, because charge conservation holds both globally -see Eq. ( 7)and at the level of each cell, n e = n p = n ( j)  p = Z ( j) /V ( j) WS , all terms depending on the cell proton density lead to a dependence on the ion abundance n ( j)  N .Using Eq. ( 3) for the translational free energy and noticing that the self-consistency is induced by the Coulomb part of the ion free energy, we obtain The last term in Eq. ( 10) is the so-called rearrangement term, which was shown to be important to guarantee the thermodynamic consistency of the model and to recover the ensemble equivalence between the MCP and OCP approaches (see e.g.Grams et al. (2018); Fantina et al. (2020); Carreau et al. (2020b); Barros et al. (2020)).Working out the rearrangement term, we find where P ( j ′ ) int is the pressure due to the Coulomb interaction (see also the discussion in Appendix A), We also note that the rearrangement term, Eq. ( 11), is not exactly equivalent to that derived in Carreau et al. (2020b).Indeed, in the latter work, the approximation that F ( j) i was only dependent on n ( j)  N was made, while in reality a dependence on the other cells arises because of the Coulomb interaction, and therefore ∂F ( j ′ ) i /∂n ( j) N 0 .Considering independent variations of n ( j)  N and introducing the single-ion canonical potential, with Eq. ( 8) becomes The ion densities can be obtained as with We note that, with respect to the derivation in Carreau et al. (2020b), here a self-consistent problem arises.Indeed, the terms ∆µ n , ūf , and Pint entering Ω( j) i on the right-hand side of Eq. ( 16) via Eqs.( 9), (4), and (11), respectively, themselves depend on the ion densities n ( j) N .Therefore, for each given baryon density and temperature, we solve the following system of equations simultaneously: together with the definition of the free volume fraction, Eq. ( 4), and the baryon number and charge conservation equations, Eqs. ( 6) and ( 7), with n ( j) N being given by Eq. ( 16).

Translational free energy of the ions in the MCP
The translational free-energy term appearing in Eq. ( 2), F ⋆,( j),MCP trans , accounts for the centre-of-mass motion of each ion of the liquid MCP.Contrarily to the OCP description of the crust (see Dinh Thi et al. (2023) for a detailed discussion of the translational term and its impact on the OCP inner-crust composition and equation of state), in the MCP picture, the centre-of-mass position of each ion is not confined to the single WS cell; the cluster can explore the whole volume, albeit with an effective mass that accounts for the fact that only bound nucleons can be in relative motion with respect to the unbound neutron gas.This leads to Eq. (3) for the translational free energy, where we set the spin degeneracy g ( j)  s = 1, because the most abundant nuclei are expected to be even-even.The effective thermal wavelength of the ion in Eq. ( 3) is given by with ℏ being the Planck constant, and the effective mass reads with γ ( j) = n gn /n ( j)  i being the ratio between the neutron-gas density and the cluster density.We note that Eq. (3) differs from the translational free energy used in Carreau et al. (2020b) in two respects: (i) here we use a different expression for the effective mass, that is, the one derived in the hydrodynamic approach considering that only neutrons in the continuum participate in the (collective) flow (see also Sect. 2 in Dinh Thi et al. ( 2023)); and (ii) the centre-of-mass position of each ion of the MCP cannot freely explore the whole volume, which is consistent with the excluded-volume approach of Eq. ( 1).In accordance with the approach of Hempel & Schaffner-Bielich (2010), the latter condition gives rise to the extra ūf term in Eq. (3).Fantina et al. (2020) noted that, because of the different expressions of the translational energy in the OCP and MCP approaches, a first deviation from the linear mixing rule appears.As a consequence, the total free energy of the ions is not just the sum of the OCP ion free energies, but an extra term arises, known in the literature as the mixing entropy (see Eq. ( 21) in Fantina et al. (2020); see also Medin & Cumming (2010)).This is still the case here, but the additional (non-ideal) 'mixing entropy' term now reflects the in-medium (excluded-volume) effects included in the translational energy F ⋆,( j),MCP trans , The computation of the free energy of the ions in the MCP therefore still requires knowledge of the ion free energy in the OCP phase, which is addressed in Sect.2.3.

Ion free energy in the OCP approximation
We define in this section the different terms of the ion free energy in the OCP approximation, F ( j),OCP i ; see Eq. ( 22)2 .A detailed discussion on these terms can also be found in Dinh Thi et al. (2023); here we present the main expressions, for completeness.
In the CLD approach, the ion free energy is decomposed in the bulk, Coulomb, finite-size, and translational contributions (see also Eq. ( 2)), The bulk free energy of the cluster, ), is expressed in terms of the free-energy density of homogeneous nuclear matter F B (n, δ, T ) at a total baryonic density n = n n + n p (n n and n p being the neutron and proton densities, respectively), isospin asymmetry δ = (n n − n p )/n, and temperature T .We note that the same expression of the bulk energy is used to compute the free energy of the free neutron gas in Eq. ( 1), To determine F B , we used the self-consistent mean-field thermodynamics (Lattimer et al. 1985;Ducoin et al. 2007) (see Sect. 3.1 in Dinh Thi et al. (2023) for details), where the free-energy density is decomposed into a 'potential' and a 'kinetic' part, The 'kinetic' term F kin in Eq. ( 24) is given by where q = n, p labels neutrons and protons, is the thermal wavelength of the nucleon with m ⋆ q being the density-dependent nucleon effective mass, F 3/2 denotes the Fermi-Dirac integral, and the effective chemical potential μq is related to the thermodynamical chemical potential µ q through μq = µ q − U q , with U q the mean-field potential.Regarding the 'potential' energy, we employed the meta-model approach of Margueron et al. (2018a,b), where V MM (n, δ) is expressed as a Taylor expansion truncated at order N (here, N = 4) around the saturation point (n = n sat , δ = 0), where x = (n − n sat )/(3n sat ) and u , with b = 10 ln(2) being a parameter governing the function behaviour at low densities.The factor u N k (x) was introduced to ensure the convergence at the zero-density limit, while the parameters v is k and v iv k are linear combinations of the so-called nuclear matter empirical parameters {n sat , E sat,sym , L sym , K sat,sym , Q sat,sym , Z sat,sym } (see Margueron et al. (2018a) for details).In this work, we employ the BSk24 empirical parameters (Goriely et al. 2013), as an illustrative example; indeed, this functional is in very good agreement with all current astrophysical data as well as nuclear physics experiments and theory (Dinh Thi et al. 2021).
The Coulomb and finite-size contributions, F Coul+surf+curv = V WS (F Coul + F surf+curv ), account for the total interface energy.As we consider here only spherical clusters, the Coulomb energy density is given by with e being the elementary charge and I = 1−2Z/A.For the surface and curvature contributions, we employed the same expression as in Lattimer & Swesty (1991), Maruyama et al. (2005), and Newton et al. (2013), that is where σ s and σ c are the surface and curvature tensions (Lattimer & Swesty 1991), The expressions of the surface and curvature tensions at zero temperature are taken from Ravenhall et al. (1983), who proposed a parameterization based on the Thomas-Fermi calculations at extreme isospin asymmetries, where y p = (1 − I)/2 and the surface parameters (σ 0 , σ 0,c , b s , β, p) were optimised for each set of bulk parameters and effective mass to reproduce the experimental nuclear masses in the 2016 Atomic Mass Evaluation (AME) table (Wang et al. 2017).The temperature dependence of the surface terms is encoded in the function h(T ) (Lattimer & Swesty 1991) where T c is the critical temperature (see Eq. (2.31) in Lattimer & Swesty (1991)).
Finally, the translational energy in the OCP approximation is given by (Dinh Thi et al. 2023) where the reduced ('free') WS volume reads (Lattimer et al. 1985) We remark that the definition of the reduced volume, Eq. ( 34), differs from that introduced in Sect.2.1 for the MCP translational free energy, Eq. ( 5).Indeed, in the MCP picture, the WS approximation is relaxed and the integral over the cluster centre of mass is extended to the whole (macroscopic) volume, except that occupied by the ensemble of the different ions, which leads to Eqs. ( 3) and (5).

MCP versus OCP
In order to compare the MCP and OCP approaches, we evaluated the most probable configuration in the MCP.The latter corresponds to the maximum probability, p max j , or equivalently the minimum grand-canonical potential, Ω( j),min i .The minimisation was performed with respect to the three variables associated to the cluster, that is (r ( j)  N , I ( j) , n ( j) i ), yielding the following system of equations: This system of variational equations can be compared to that obtained in the OCP approach of Carreau et al. (2020b) and Dinh Thi et al. (2023), without the free proton gas, It is easy to show that Eqs. ( 35)-(37) for the most probable ion in the MCP are equivalent to the first three OCP variational equations Eqs. ( 38)-( 40) if the following conditions are satisfied: (i) F ⋆,( j) i is the same as in the OCP, that is, F ⋆,( j) i = F ( j),OCP i = F i (see Sect. 2.3; this is the case if the non-linear mixing term, which arises due to the translational motion (Gulminelli & Raduta 2015;Fantina et al. 2020), see Sect.2.2, is negligible); and (ii) the gas densities, n gn and n e , as well as ∆µ n and Pint , are identical to the OCP values.This latter condition is expected to be satisfied if the ion distributions are strongly peaked on a unique ion species, such that the averages in Eqs. ( 18)-( 19) can be approximated to a single term corresponding to the OCP solution.The validity of these conditions, which are necessary for the OCP approximation to give a satisfactory description of the finite-temperature configuration, are discussed in Sect.3.1.
The additional OCP condition in Eq. ( 41) is identical to the well-known Baym virial theorem (Baym et al. 1971), which holds for a nucleus in vacuum.This condition arises from the fact that, for a given thermodynamic condition, there is systematically one more independent variable in the OCP than in the MCP.Indeed, the WS cell volume in the MCP is simply defined by charge conservation, while in the OCP it corresponds to an additional variable which has to be variationally determined.
In a CLD picture, for a given mass number A and atomic number Z, a nuclear species is also characterised by its radius, r N , or equivalently its density n i = 3A/(4πr 3 N ), which can in principle fluctuate from cell to cell in the MCP equilibrium.Changing variables from (r N , I, n i ) to (A, Z, n i ), the ion-radius distribution of a given (A, Z) nucleus is expressed as where N is a normalization.This distribution will be peaked at a value r N corresponding to ∂ Ω( j) i /∂n i | A ( j) ,Z ( j) = 0, yielding the following pressure equilibrium condition: Equation ( 38) is not necessarily compatible with the so-called virial condition, or radius-minimisation ansatz given by As a consequence, even if conditions (i) and (ii) above are met, the minimisation equation for the OCP radius, Eq. ( 41), is not necessarily fulfilled in the MCP.This first important difference between the OCP and MCP approaches is presented in Fig. 1, where we plot the normalised distribution of the cluster radius for the most probable (A, Z) nucleus obtained at different thermodynamic conditions (orange solid lines).This distribution is compared with the solution obtained from the condition of pressure equilibrium in Eq. ( 43) (green dash-dotted lines) and that resulting from the minimisation of the free energy per nucleon with respect to r N in Eq. ( 44) (red solid lines).To obtain these distributions, at each thermodynamic condition defined by (n B , T ), we considered independent variations of the three variables (A ( j) , Z ( j) , r ( j) N ), and solved the system of equations given by Eqs. ( 6), ( 7), ( 4), (18), and ( 19), together with the equation for the ion density n ( j) N , Eq. ( 16) (see also Sect.3.2).From Fig. 1, one can also observe that the solution of the pressure equilibrium equation always coincides with the peak of the MCP distribution.This is true not only for the radius of the most probable cluster, but also for the r ( j) N of any (A ( j) , Z ( j) ) cluster.On the other hand, using Eq. ( 44) to determine the cluster radius fails in reproducing the most probable cluster radius obtained within the MCP approach.In fact, it only produces the correct peak for the most probable cluster if the OCP chemical potentials are used in the absence of the non-linear mixing term, as it was analytically demonstrated in Grams et al. (2018).
Computation of the crust observables considering the cluster distribution in the three-dimensional variable space (A, Z, r N ) is a heavy task from a numerical point of view.In the literature, most MCP codes designed to determine the composition and  To limit the numerical cost, in the following sections we only consider a two-dimensional cluster distribution p AZ , and fix for each (A, Z) the cluster radius r N from the pressure equilibrium condition using Eq. ( 43).Indeed, we verified that, in all thermodynamic conditions explored in the present work, the radius distribution is always symmetric and relatively narrow.

Crust composition at finite temperature in the MCP
Using the formalism described in Sect.2.1, we calculated the composition and equation of state of the inner crust of a PNS, employing the BSk24 (Goriely et al. 2013) empirical parameters.We performed our analysis for temperatures in the range from 1 MeV to 2 MeV, where the crust is expected to be in the liquid phase (Carreau et al. 2020a) and the beta-equilibrium condition be realised.

Perturbative MCP calculations
We begin this section by discussing the results obtained by calculating the MCP distribution in a perturbative approach, as first proposed in Grams et al. (2018).For each (n B , T ), we first calculated the OCP solution, solving the system of equilibrium equations Eqs. ( 38)-( 41).This yields the OCP composition (i.e. A OCP , Z OCP , and the nuclear radius r OCP N ), the OCP chemical potentials, and the neutron and electron densities, n gn = n OCP gn , .With the latter five quantities as input, and the MCP calculations were performed, giving the normalised probabilities p j ≡ p AZ .
As discussed in Sect.2.4, the OCP and MCP results should coincide if, in addition, the non-linear mixing term coming from the translational free energy is set to zero.This is illustrated in Fig. 2, where the cluster radius r N (blue curves), proton number Z (red curves), and mass number A (black curves) predicted by the OCP (solid lines) are plotted as a function of the input baryon density, n B , together with the MCP average (dashed lines) and most probable quantities (diamonds).The latter values correspond to those maximising p AZ .For all the considered densities (from the neutron-drip up to the crust-core transition density) and temperatures, we can see that the most probable A, Z, and r N coincide with the OCP solutions, and follow the average values very closely, implying that the distributions are centred around the OCP predictions.This is further shown in Fig. 3, where we display the joint distributions of the cluster proton number Z and mass number A for the same temperatures as in Fig. 2 and for three selected densities, n B = 2 × 10 −3 fm −3 (green contours), n B = 10 −2 fm −3 (orange contours), and n B = 2 × 10 −2 fm −3 (blue contours).The distributions are indeed Gaussian-like, peaked at the OCP solutions (black stars).We can also observe that, while the proton number Z is relatively constant, the most probable mass number A increases with density, as already noted in Carreau et al. (2020b) and Dinh Thi et al. (2023).Moreover, the width of the distributions gets broader with temperature and density, with lighter clusters being populated at higher T and n B (see panel c), underlying the importance of considering a full nuclear ensemble instead of a single-nucleus (OCP) approach.
We now turn to a discussion of the influence of the non-linear mixing term, which comes from the translational motion.In Fig. 4, we show the evolution with the baryon density of the cluster variables A, Z, and r N , for the same conditions as in Fig. 2, but with the translational free energy included in both the OCP and the MCP calculations, Eqs. ( 33) and (3), respectively.We can see that the different description of the centre-of-mass motion induces a discrepancy between the predictions of the two approaches.Indeed, in the OCP approximation, the cluster moves in the reduced ('free') volume V OCP f -see Eqs. ( 33)-( 34)-associated with the single WS cell determined from the variational procedure, while in the MCP, clusters of different species are considered to move in the same macroscopic volume, leading to an average 'free' volume ⟨V f ⟩; see Eqs. (3)-( 5).This correlation between the different ion species breaks the linear mixing rule.At lower temperatures (see panel (a) of Fig. 4), the deviation between the MCP and OCP predictions is negligible at densities below about 0.05 fm −3 , implying that the influence of the non-linear mixing term is not very important.On the other hand, from panels (b) and (c) of Fig. 4 we can see that, with increasing temperature, the discrepancy between the two approaches starts to become significant at progressively lower densities.This can be explained by the fact that the contribution of the translational free energy becomes more important at higher T and n B , as dis- cussed in Dinh Thi et al. (2023).In particular, we can see that the distributions in the (perturbative) MCP are shifted to larger clusters with respect to the OCP solution, and that the crust-core transition occurs earlier.This is because Eq. ( 41), which results in small clusters in the OCP approximation (see Dinh Thi et al. ( 2023)), no longer holds in the MCP, as discussed in Sect.2.4.In addition, the contribution of the translational free-energy term in the OCP approximation is more important than that in the MCP.Indeed, in the former, the free volume term is a variable entering the minimisation procedure, whereas ⟨V f ⟩ in the MCP is the common volume for all nuclear species.As a result, it does not affect the nuclear distribution; see Eqs. ( 13)-( 14).The discontinuous behaviour of the most probable (A, Z) cluster in Fig. 4 is due to the fact that, as the temperature increases, the contribution of very light clusters, such as neutronrich helium isotopes, becomes increasingly favoured.As shown in Fig. 5, the distributions become double-peaked when the temperature overcomes k B T ≈ 1 MeV, and a density domain exists where the helium abundance overcomes that of the heavy cluster in the iron region predicted by the OCP.Still, it is also visible from Fig. 5 that the light-cluster contribution is limited to a small number of helium isotopes, while a large variety of nuclear species have a comparable probability around the (A OCP , Z OCP ) value.Therefore, if we consider the probabilities of each element Z by integrating the distributions over the isotopic content (i.e. over A), the (perturbative) MCP predictions appear to be in closer agreement with the OCP approximation, at least at relatively low densities, as can be seen from Fig. 6.
The results presented in this section confirm that the nonlinear mixing term induced by the translational energy leads to a breakdown of the ensemble equivalence between the OCP and MCP predictions.
The perturbative MCP approach employed here has the clear advantage of computing the full nuclear distributions at a reduced computational cost.However, this approach is not fully self-consistent, because the gas variables and chemical potentials resulting from the OCP solutions might not exactly satisfy the constraints of baryon number conservation and charge neutrality in the MCP, Eqs. ( 6) and ( 7), respectively.For this reason, we computed the full self-consistent MCP, which is discussed in the following section.

Self-consistent MCP calculations
In this section, we present the results obtained by performing fully self-consistent MCP calculations.At each thermodynamic condition defined by (n B , T ), we solved the system of equations given by Eqs. ( 6), ( 7), ( 4), (18), and ( 19), together with the equation for the ion density n ( j) N , Eq. ( 16).The translational free energy is always included in the computation of the ion free energy.
With Fig. 7, we begin by showing the evolution with the baryon number density of the neutron chemical potential µ n presented in Eq. ( 9) (panel a), the electron chemical potential potential µ e = dF e /dn e (panel b), the average free-volume fraction ūf in Eq. ( 4) (panel c), and the MCP interaction pressure Pint presented in Eq. ( 19) (panel d).All these quantities enter the computation of the ion abundances via Eq.( 17).For illustrative purposes, we only plot the results at one selected temperature, k B T = 1 MeV.The solutions from the MCP calculations are shown by solid lines, while the dotted lines correspond to the values calculated at the OCP composition, that is, µ OCP n , µ OCP e , ūOCP f , and POCP int .As we can see from panels (a) and (b), the chemical potentials in the MCP are very similar to those in the OCP approximation; indeed, the two curves are almost indistinguishable.However, from the inset in panel (a) we can notice that the neutron chemical potential in the MCP is actually slightly lower than the OCP counterpart, suggesting that the corresponding gas density is not the same in the two approaches, and in particular that it is smaller in the MCP.The same trend is observed for the electron chemical potential (see panel (b) in Fig. 7), except at densities above 0.06 fm −3 , where µ e calculated in the MCP is slightly above that obtained in the OCP approach, and therefore the corresponding density is higher in the MCP.The lower density of the unbound nucleons in the MCP was also observed by Burrows & Lattimer (1984) and by Gulminelli & Raduta (2015) (see their Fig.13).Although the difference between the MCP and OCP chemical potentials is numerically small, given that these latter enter in the calculation of the ion abundances through the exponential, Eq. ( 16), this difference could still lead to a significant deviation between the MCP and OCP results, as already noted by Gulminelli & Raduta (2015).From panels (c) and (d) in Fig. 7  almost identical to those calculated in the OCP approximation 3 .We therefore expect that, at low density and temperature, the results shown in Fig. 4 are a good approximation.On the other hand, at higher densities, the absolute value of the interaction pressure, which directly enters into the computation of the rearrangement term, tends to zero in the MCP approach, while | POCP int | increases.As the latter is an input in the perturbative MCP, we might also expect the discrepancies between the perturbative and self-consistent MCP approaches to become more pronounced at higher densities.
Indeed, this is shown in Fig. 8, where the average (dashdotted violet lines) and the most probable (violet squares) values of Z (top panel) and A (lower panel) in the self-consistent MCP calculations are plotted as a function of n B at k B T = 1 MeV, together with the OCP predictions (black solid lines) and the average (orange dashed lines) and most probable (orange dots) values obtained within the perturbative MCP procedure.Below   8. Average (dash-dotted violet lines, labelled 'avg-µ sc ') and most probable (violet squares, labelled 'MP-µ sc ') values of the cluster proton number Z (top panel) and mass number A (lower panel) in the selfconsistent MCP calculations as a function of the total baryon density n B at k B T = 1 MeV.For comparison, the OCP solutions (black solid lines), as well as the average (orange dashed lines labelled 'avg-µ OCP ') and most probable (orange dots labelled 'MP-µ OCP ') values in the perturbative MCP are also plotted.about 0.02 fm −3 , the results of the three different approaches almost coincide.However, as the density increases, the outcomes of the different methods start to diverge, and therefore the perturbative MCP and OCP predictions become less reliable.On the one hand, these results confirm the validity of these approximations at relatively low densities, while on the other, they highlight the importance of full MCP calculations in the deeper region of the PNS inner crust.
From Fig. 8, it is also interesting to see that, in the selfconsistent MCP calculations, lighter nuclei dominate at high density, which is in agreement with previous works based on the nuclear statistical equilibrium (see e.g.Souza et al. (2009); Gulminelli & Raduta (2015)).Conversely, in the OCP and perturbative MCP approximations, heavier clusters survive until the crust-core transition.
In order to further investigate this point, we plot in Fig. 9 the normalised distributions of the cluster proton number Z (left panel), A p AZ , and mass number A (right panel), Z p AZ , at different densities in the inner crust, n B ∈ [0.001, 0.04] fm −3 , at k B T = 1 MeV, for both the self-consistent MCP calculations (violet areas) and the perturbative MCP ones (orange areas).As one may expect, at lower densities, the distributions obtained in the two approaches are almost identical; they are narrow and peaked at the OCP predictions (green triangles).Therefore, employing a perturbative MCP where the gas variables are fixed from the OCP solution is a very good approximation in these thermodynamic conditions.As the density increases, at n B ≈ 0.02 − 0.03 fm −3 , the self-consistent MCP distributions are displaced towards lower values of A and Z, and even exhibit a double-peak structure, meaning that lighter clusters coex-ist with heavier ones and their contribution becomes important.In these regimes, the OCP and perturbative MCP predictions significantly overestimate the cluster mass and proton number.Eventually, at the bottom of the crust, nuclei with Z < 5 and A < 50 may even dominate, with helium being the most probable element, although the tails of the distributions extend up to Z ≃ 50 and A ≃ 400, implying that heavier nuclei may still be present near the crust-core transition.However, it is to be noted that, while the Z = 2 element appears to be the most abundant, 4 He, which is the most bound helium isotope in the vacuum, is not the most abundant cluster (see also Fig. 8).This is because in the NS inner crust, the medium is very dense and (helium) clusters are immersed in continuum neutron states.Indeed, as already noted by Gulminelli & Raduta (2015), alpha particles are abundant for matter close to isospin symmetry, while more neutron-rich (hydrogen and helium) isotopes prevail in neutronrich matter.
The findings presented above were obtained at a given temperature, namely k B T = 1 MeV.In order to study the dependence of the results on the temperature, we ran the calculations for different temperatures up to k B T = 2 MeV.In Fig. 10, we plot the distributions of Z (left panel) and A (right panel), defined as in Fig. 9, at n B = 0.001 fm −3 in the aforementioned temperature range.As in Fig. 9, the results from both the self-consistent MCP calculations (violet areas) and the perturbative MCP ones (orange areas) are displayed and compared with the OCP solutions (green triangles).At lower temperatures, the peaks of the distribution coincide with or are very close to the OCP predictions.However, as the temperature increases, more nuclear species are populated, and therefore the distributions flatten and the OCP prediction tends to overestimate the cluster size.At higher temperature, k B T ≥ 1.8 MeV, small clusters start to appear.At 2 MeV, the double-peak structure already observed at high density is clearly visible in the self-consistent MCP calculations, while it is absent in the perturbative MCP results.Nevertheless, the deviation between the predictions of the two approaches is not as significant as in Fig. 9, which suggests that the effect of density overcomes that of temperature as far as the cluster distributions are concerned.
A complementary comparison of the distributions obtained in the perturbative and self-consistent MCP is given in Fig. 11, where the normalised distributions A p AZ are plotted as a function of Z for k B T = 2 MeV and for three selected densities: n B = 0.01 fm −3 (panel a), n B = 0.02 fm −3 (panel b), and n B = 0.03 fm −3 (panel c).These conditions correspond to those for which a double-peaked distribution is observed in the perturbative MCP calculations (see Fig. 5).For comparison, the OCP results are also shown (green triangles).At all three densities considered, the distributions obtained with the self-consistent MCP are shifted towards lower Z with respect to those of the perturbative MCP calculations, as already observed in Fig. 9 at k B T = 1 MeV.Indeed, at k B T = 2 MeV, in the self-consistent MCP, light nuclei are always more abundant than heavier ones.Moreover, the distributions become peaked around Z ≈ 2 (see panels (b) and (c)), while they remain relatively broad and peaked towards Z ≈ 30 in the perturbative MCP.
The results presented in Figs.8-11 show that even for temperatures as low as 1 MeV (and more so at higher temperatures), the OCP and the perturbative MCP approximations are no longer reliable in the deepest region of the crust.Therefore, for studies requiring accurate knowledge of the crust composition, a full MCP calculation is needed.However, the effect is less important as far as more global properties such as the equation of state are concerned.This is illustrated in Fig. 12, where we display the neutron and electron gas density over the baryon number density, n gn /n B and n e /n B , respectively, and in Fig. 13, where we plot the total pressure versus the mass-energy density in the PNS inner crust.Results are shown for three different temperatures, k B T = 1 MeV (red lines), 1.5 MeV (blue lines), and 2 MeV (black lines), for both the self-consistent MCP (solid lines) and the OCP approximation (dotted lines).From Fig. 12, we can see that the MCP neutron gas densities are slightly smaller than the OCP ones in all the considered cases, although differences of up to about 20% are observed for k B T = 2 MeV at the highest densities in the crust.As for the electron density, that calculated within the MCP is generally slightly higher than that calculated in the OCP approximation, except possibly at k B T = 1 MeV around a few 10 −2 fm −3 , in accordance with Fig. 7. From Fig. 13, we note that, concerning the equation of state, at k B T = 1 MeV, the pressure-density relations provided within the self-consistent MCP approach and the OCP approximation are very similar at all densities, except very close to the crust-core transition.At higher temperatures, deviations between the two calculations start to emerge, particularly at high density.However, these discrepancies amount to approximately 10% at most in the vicinity of the crust-core transition.

Impurity parameter
Within our MCP approach, it is also possible to self-consistently calculate the impurity parameter, which is defined as the variance of the Z distribution, The latter is plotted in Fig. 14 for three temperatures, k B T = 1.0 MeV (red lines), k B T = 1.5 MeV (blue lines), and k B T = 2.0 MeV (black lines).Results obtained with the self-consistent (perturbative) MCP calculations as a function of the baryon density n B in the crust are shown as solid (dash-dotted) lines.At the lower temperature of k B T = 1.0 MeV (red lines), the predictions from the two treatments coincide at lower densities, until n B ≈ 0.01 fm −3 .With increasing density, the impurity parameter is, at first, larger in the self-consistent MCP approach with respect to the perturbative one.However, at variance with the perturbative MCP predictions, Q imp does not increase monotonically in the self-consistent calculations, but peaks around n B ≈ 0.025 fm −3 , reaching Q imp ≈ 100 for the considered BSk24 model, and subsequently decreases at higher densities.A similar behaviour is observed for higher temperatures, although the discrepancy between the two treatments and the peak in the impurity parameter appear at smaller densities.This can be understood from the Z distributions shown in the left panel of Figs.9-10 and in Fig. 11.Indeed, while the presence of the second peak at small Z for moderate densities increases the variance of Z, the transition to light nuclei at high density overall decreases the value of the impurity parameter.These findings also show that the impurity parameter calculated within the perturbative MCP approach is underestimated at low densities, and is severely overestimated at high densities.
In the cooling process of a NS, it is reasonable to assume that the crust composition is frozen after the solidification of the crust.Neutron absorption or β-decays might still occur at lower temperatures, but pycnonuclear reactions that involve overcoming a Coulomb barrier might be considerably inhibited even above crystallisation.A realistic estimate of the temperature at which the ion distribution is frozen and the impurity parameter is settled would require a comparison between the cooling time and the different reaction rates, which is beyond the scope of the present study.As a reasonable estimate, the impurity parameter of the solid crust can be computed from that obtained at the crystallisation temperature of a pure Coulomb plasma (Haensel et al. 2007):  We can see that a high impurity parameter 10 ≲ Q imp ≲ 100 should be expected in the whole inner crust, which could have important consequences for the magneto-thermal evolution of Xray pulsars (see Pons et al. (2013); Newton (2013); Horowitz et al. (2015)).
For comparison, we also plot Q imp for four different selected temperatures: k B T = 0.5 MeV (magenta points), k B T = 1.0 MeV (red points), k B T = 1.5 MeV (blue points), and k B T = 2.0 MeV (black points).We can see that the general behaviour of Q imp is similar for all temperatures; that is, all curves show a peak and a subsequent drop.Indeed, even at a relatively low temperature of k B T = 0.5 MeV, Q imp is still reduced at high densities, n B ≈ 0.06 fm −3 , indicating that the charge distribution is dominated by small Z.Moreover, as expected, a higher Q imp is obtained if the composition is frozen at higher temperature, but the sharp drop occurs at lower density because of the appearance of light nuclei (see also .However, our results should be considered with care at high density, because non-spherical pasta structures could appear, and are not considered in the present paper.Such structures are also expected to be associated with high impurity parameter values Q imp ∼ 30; see e.g.Schneider et al. (2016) and Caplan et al. (2021).Also, the CLD description can be questioned for light clusters and a more microscopic treatment of the in-medium modifications might be needed; see Pais et al. (2018) and Röpke (2020).
The results for the impurity parameters shown in Fig. 15 are given in Table 1 and are available in tabular format at the CDS.For practical applications to numerical simulations, we also provide a fitting formula: where x B ≡ n B /fm −3 , x T ≡ k B T/MeV, a 0 = 3.328, and a k j are given in Table 2, and the function ∆Q reads with f 0 = 65.55,f 1 = 9.046 × 10 −2 , f 2 = −6.324,f 3 = −1.087× 10 −2 , and f 4 = 4.802.The variable x BT in Eq. ( 51) depends on both n B and T and is defined as x BT (n B , T ) ≡ 5 + log 10 x B 2 x T , (52) while f (n B , T ) is given by The impurity parameter calculated at different temperatures with the fitting expression Eq. ( 50) is shown by solid lines in panel (b) of Fig. 15, while the absolute error with respect to the computed values of Q imp in the MCP approach is displayed in panel (c).We can see that the error remains relatively small, |∆Q imp | ≈ 5 at most, in the whole inner-crust density range.Indeed, even if the fit tends to break down at higher densities and higher temperatures (see black and blue lines in panel (b) of Fig. 15 for n B ≳ 0.05 fm −3 ), the error is very small because of the low values of Q imp ≈ 0.1 − 0.2.

Conclusions
In this work, we studied the properties of the beta-equilibrated inner crust of a PNS in the liquid phase within a multicomponent approach.To this aim, the formalism of Fantina et al. (2020) and Carreau et al. (2020b) was extended for its applicability over a wider range of densities and temperatures in the inner crust, from the neutron drip up to the crust-core transition.a k j j = 1 j = 2 j = 3 k = 0 −2.285 × 10 1 2.205 × 10 1 −5.370 k = 1 1.818 × 10 3 −1.087× 10 3 2.000 × 10 2 k = 2 −9.456 × 10 4 6.849 × 10 4 −2.673 × 10 4 k = 3 1.730 × 10 6 −5.766 × 10 5 4.969 × 10 5 k = 4 −1.161 × 10 7 6.356 × 10 6 −1.666 × 10 7 The cluster energetics was described within a CLD model approach, in which the surface parameters were optimised to reproduce the experimental masses from the AME2016 (Wang et al. 2017), and the nuclear matter properties were calculated within finite-temperature mean-field thermodynamics.We performed the calculations employing both a perturbative MCP treatment, where the chemical potentials and gas densities were taken from the OCP calculations, and a full self-consistent MCP approach.
We find that, if non-linear mixing terms arising from the centre-of-mass motion of the ion are neglected, the most probable clusters in the MCP (perturbative) approach coincide with the OCP predictions.However, as already discussed in Dinh Thi et al. (2023), this translational free-energy term should be taken into account in the calculations of the (liquid) inner-crust composition, leading to a breakdown in the ensemble equivalence.Our outcomes also show that the OCP and the perturbative MCP treatments are a good approximation at relatively low densities and temperatures, especially as far as global properties such as the equation of state are concerned.However, we believe that a full self-consistent MCP is needed for reliable predictions of the PNS crust composition, particularly in the deeper region of the crust.
Moreover, our findings reveal that, with increasing density and temperature, the abundance of light nuclei becomes important, and eventually dominates the whole distribution.This result has an important effect on the impurity parameter, and therefore potentially on the NS cooling.For applications to numerical simulations, we also provide both values of the impurity parameter calculated within our fully self-consistent MCP approach in tabular form and a fitting formula able to reproduce the numerical results with good accuracy.
In our study, we only considered spherical clusters.However, non-spherical structures, so-called pasta phases, may appear at the bottom of the inner crust (see Pelicer et al. (2021) for recent calculations of pasta phases in a MCP approach at finite temperature within a relativistic mean-field functional).We defer an investigation of the co-existence of different (non-spherical) structures within our approach to future studies.

Fig. 1 .
Fig. 1.Normalised distribution of the cluster radius, p AZ (r N ) (orange solid line), in a full MCP calculation.For comparison, the solutions from two different equilibrium prescriptions for the most probable cluster are displayed.Results are presented for two different temperatures: k B T = 1 MeV (top panels) and k B T = 1.5 MeV (bottom panels), and two different densities: n B = 10 −3 fm −3 (left panels) and n B = 10 −2 fm −3 (right panels).See text for details.

Fig. 2 .
Fig. 2. Cluster radius r N (blue), proton number Z (red), and mass number A (black) as a function of the baryon density n B at different temperatures: k B T = 1.0 MeV (panel a), k B T = 1.5 MeV (panel b), and k B T = 2.0 MeV (panel c).OCP results are shown by solid lines, while dashed lines and symbols correspond to the average (avg) and the most probable (MP) value, respectively, of the cluster distribution in a perturbative MCP calculation.The translational free energy is not included in the calculations.

Fig. 3 .
Fig. 3. Joint distributions of the cluster proton number Z and mass number A for the same temperatures as in Fig. 2 and for three selected input baryon densities: n B = 2 × 10 −3 fm −3 (green contours), n B = 10 −2 fm −3 (orange contours), and n B = 2 × 10 −2 fm −3 (blue contours), in a perturbative MCP calculation.The black stars indicate the OCP solution.The translational free energy is not included in the calculations.

Fig. 4 .
Fig.4.Same as in Fig.2but with the translational free energy, F ⋆ trans , included in both the OCP and perturbative MCP calculations.

Fig. 5 .
Fig. 5. Normalised probability p AZ in a perturbative MCP calculation as a function of the cluster proton number Z and cluster proton fraction Z/A.Results were obtained at the selected temperature of k B T = 2.0 MeV for three different baryon densities: n B = 0.01 fm −3 (top panel), n B = 0.02 fm −3 (middle panel), and n B = 0.03 fm −3 (bottom panel).The green arrow in each panel indicates the corresponding OCP solution.The translational free energy is included in both MCP and OCP calculations.

Fig. 6 .
Fig. 6.Normalised distributions P(Z) = A p AZ , where p AZ are taken from Fig. 5, as a function of Z in a perturbative MCP calculation.Calculations are performed for k B T = 2 MeV and n B = 0.01 fm −3 (green dash-dotted line), n B = 0.02 fm −3 (red solid line), and n B = 0.03 fm −3 (black dashed line).Arrows indicate the OCP solutions.

Fig. 7 .
Fig. 7. Self-consistent MCP solution (solid lines) for the neutron chemical potential µ n (panel a), electron chemical potential µ e (panel b), average free volume fraction ūf (panel c), and MCP interaction pressure Pint (panel d) as a function of the total baryon density n B at k B T = 1 MeV.The dotted lines show the OCP predictions.The tick marks of the x and y axes in the insets of panels (a) and (b) are spaced by 2 × 10 −4 fm −3 and 0.1 MeV, respectively.

Fig.
Fig.8.Average (dash-dotted violet lines, labelled 'avg-µ sc ') and most probable (violet squares, labelled 'MP-µ sc ') values of the cluster proton number Z (top panel) and mass number A (lower panel) in the selfconsistent MCP calculations as a function of the total baryon density n B at k B T = 1 MeV.For comparison, the OCP solutions (black solid lines), as well as the average (orange dashed lines labelled 'avg-µ OCP ') and most probable (orange dots labelled 'MP-µ OCP ') values in the perturbative MCP are also plotted.

Fig. 9 .Fig. 10 .
Fig. 9. Normalised probability distributions of the cluster proton number Z (left panel), A p AZ , and mass number A (right panel), Z p AZ , at k B T = 1 MeV for different baryon densities, n B ∈ [0.001, 0.04] fm −3 .The violet distributions are obtained from the self-consistent MCP calculations, while the orange ones correspond to the perturbative MCP.The OCP solutions are marked by the green triangles.

Fig. 11 .
Fig. 11.Normalised distributions P(Z) = A p AZ as a function of Z in a fully self-consistent ('SC') and a perturbative ('P') MCP calculation.Calculations are performed for k B T = 2 MeV and n B = 0.01 fm −3 (panel a), n B = 0.02 fm −3 (panel b), and n B = 0.03 fm −3 (panel c).Green triangles indicate the OCP solutions.

Fig. 13 .Fig. 14 .Fig. 15 .
Fig.13.Total pressure P vs. mass-energy density ρ B resulting from the self-consistent MCP calculations (solid lines) at k B T = 1 MeV (red lines), 1.5 MeV (blue lines), and 2 MeV (black lines).For comparison, the OCP results are plotted with dotted lines.In the inset, the tick marks on the x and y axis are spaced by 10 12 g/cm 3 and 2 × 10 30 dyn/cm 2 , respectively.
Ratio of the neutron gas density (upper panel) and of the electron gas density (lower panel) over the baryon density, n gn /n B and n e /n B , respectively, as a function of the total baryon density n B resulting from the self-consistent MCP calculations (solid lines) at k B T = 1 MeV (red lines), 1.5 MeV (blue lines), and 2 MeV (black lines).For comparison, the OCP results are plotted with dotted lines.

Table 1 .
Impurity parameter Q imp in the inner crust for different temperatures.The results are obtained with the fully self-consistent MCP calculation.

Table 2 .
Fitting parameters a k j in Eq. (50) for the impurity parameter.