Crystallization of the outer crust of a non-accreting neutron star

The interior of a neutron star (NS) is usually assumed to be made of cold catalyzed matter. However, the outer layers are unlikely to remain in full equilibrium during the formation of the star and its cooling, especially after crystallization. We study the cooling and equilibrium composition of the outer layers of a NS down to crystallization. Here the impurity parameter, usually a free parameter in cooling simulations, is calculated self-consistently using a microscopic model for which a unified equation of state has recently been determined. We follow the evolution of the nuclear distributions of the multi-component Coulomb liquid plasma (MCP) fully self-consistently, adapting a general formalism originally developed for the description of supernova cores. We calculate the impurity parameter at the crystallization as determined in the one-component plasma (OCP) approximation. Our analysis shows that the sharp changes in composition obtained in the OCP approximation are smoothed out when a full nuclear distribution is allowed. The Coulomb coupling parameter at melting is found to be reasonably close to the canonical value of 175, except for specific pressures for which supercooling occurs in the OCP approximation. Our MCP treatment leads to non-monotonic variations of the impurity parameter with pressure. Its values can change by several orders of magnitude reaching about 50, suggesting that the crust may be composed of an alternation of pure (highly conductive) and impure (highly resistive) layers. The results presented here complement the recent unified equation of state obtained with the same model. Our self-consistent approach to hot MCP shows that the presence of impurities in the outer crust of a NS is non-negligible and may have a sizeable impact on transport properties. In turn, this may have important implications for the cooling of NS and their magneto-rotational evolution.


Introduction
Formed in the aftermath of gravitational core-collapse supernova explosions, neutron stars (NSs) are initially very hot. With temperatures exceeding 10 10 K, the outer part of the newly born NS is expected to be made of a dense Coulomb liquid containing various nuclear species in a charge compensating electron background. It is generally assumed that as the NS cools down by emitting neutrinos and photons, this multi-component plasma (MCP), which crystallizes at the temperature T m , remains in full thermodynamic equilibrium (with respect to all possible processes) until the ground state at T = 0 K is eventually reached. According to this so-called 'cold catalyzed matter' hypothesis, the outer crust of a mature NS is, thus, expected to be stratified into pure layers, each of which consists of a one-component Coulomb crystal (except, possibly, at the ⋆ The table of the impurity parameter at the crystallization temperature shown in Fig. 8 is 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/ boundaries between adjacent layers; see Chamel & Fantina (2016a) for a discussion).
However, if the interior of a NS cools down rapidly enough in comparison to the various reaction rates, the composition of the stellar material may be frozen at some finite temperature so that the ground state may never be attained, see. e.g. Goriely et al. (2011) (see also Haensel et al. (2007); Chamel & Haensel (2008)). Even in the simplified scenario of an adiabatic cooling process, the full thermodynamical equilibrium of the outer layers of the star is unlikely to be maintained after the crystallization, meaning that a more realistic picture of the outer crust of a NS is that of a multi-component Coulomb solid. With the crystallization temperature as low as ≈ 10 6 − 10 7 K (see Haensel et al. (2007)), the most probable ion species would presumably be close to or coincident with the one corresponding to the ground state. Consequently, the static properties of the frozen crust are not expected to be appreciably different from those of catalyzed crust. On the other hand, the co-existence of various nuclear species may have a dramatic impact on transport properties. However, the nuclear distributions in different crustal layers are to a large extent unknown. For this reason, NS cooling simulations have been generally performed using the groundstate composition. The presence of other nuclear species is taken into account by introducing an 'impurity factor', treated as a free parameter directly fitted to the cooling data. This parameter is important not only for thermal properties but for other transport properties as well, such as electrical conductivity (see e.g. Schmitt & Shternin (2018) for a recent review). The presence of impurities is, thus, also expected to have a strong impact on the magnetorotational evolution of NSs, see e.g. Pons et al. (2013) (see also Gourgouliatos & Esposito (2018) for a recent review).
In this paper, we study the composition and formation of the outer crust of a non-accreting unmagnetized NS. After determining the crystallization temperature in the one-component plasma (OCP) approximation, the nuclear distributions and the impurity parameter are calculated fully self-consistently, adapting a general formalism originally developed for the description of a hot dense MCP under conditions prevailing in supernova cores (Gulminelli & Raduta 2015;Grams et al. 2018). Our treatment of a OCP and a MCP plasma are presented in Sections 2 and 3, respectively. Results are discussed in Sect. 4 and conclusions are drawn in Sect. 5. In Appendix A, we derive the expression for the pressure of the MCP, while in Appendices B and C, we report for completeness the expressions used in this work for the free energy and pressure of the uniform electron gas and for the free energy of the Coulomb plasma of ions, respectively.

Main assumptions
In this study, we consider matter at densities high enough so that full ionization can be supposed, i.e. ρ 11AZ g cm −3 , which for iron, whose mass number A and charge number Z are A = 56 and Z = 26, yields ρ 10 4 g cm −3 . The nuclei are, thus, surrounded by a gas of highly-degenerate electrons, matter being electrically charge neutral. At finite temperatures, a free nucleon (neutron and proton) gas could also be present. However, this gas is expected to be very dilute at temperatures T < 3 × 10 9 K, which are of interest here (see, e.g. Haensel et al. (2007)). We shall, therefore, ignore the nucleon gas.
The properties of such dense matter in full (beta) equilibrium at temperature T and pressure P are determined by minimizing the Gibbs free energy under the constraint of baryon number conservation. In the OCP (single-nucleus) approximation, this procedure yields the mass number and charge number of the (unique) equilibrium nucleus (A, Z) at each temperature T in each layer at pressure P (see e.g. the pioneer works of Tondeur (1971); Baym et al. (1971) at T = 0 K). As a consequence, the baryon number density n B may vary discontinuously at the interface between two adjacent layers. These density jumps may be reduced (though not entirely removed) if one allows for the existence of multinary compounds (see Chamel & Fantina (2016a) for a recent discussion).
The total Gibbs free energy per nucleon g to be minimized is defined as where f is the total free energy per nucleon 1 and the baryon density n B is numerically calculated from the pressure P . The total free energy per ion reads In this expression, F e is the electron free energy, that accounts for the free (non-interacting) part, plus the corrections (exchange and correlation) in a uniform electron system. The term F i corresponds to the ion free energy including the Coulomb contribution, and is given by (see Chap. 2 in Haensel et al. (2007)) where M ′ (A, Z) is the ion mass (which coincides with the nuclear mass since atoms are fully ionized), c being the speed of light, F id i is the non-interacting ("ideal") contribution to the ion free energy, and F int i accounts for interactions. Specifically, F int where F ii includes all the Coulomb interactions (between ions, between electrons, and between ions and the uniform electron gas) and F pol ie represents the polarization correction that accounts for the deviation of the electron background from uniformity. For M ′ (A, Z), we make use of experimental masses, whenever available, from the 2016 Atomic Mass Evaluation (AME) (Wang et al. 2017), supplemented with the microscopic HFB-24 theoretical mass table based on the nuclear energydensity functional theory 2 . The underlying functional has been recently used to determine the ground-state composition and the equation of state in all regions of a non-accreting NS (Pearson et al. 2018). Usually, atomic masses are tabulated instead of the nuclear ones, which can be calculated as where M (A, Z)c 2 = ∆ǫ + Am u c 2 is the atomic mass (∆ǫ being the mass excess and m u being the atomic mass unit), m e is the electron mass, and B el is the binding energy of the atomic electrons (see Eq. (A4) in Lunney et al. (2003)) B el = 1.44381 × 10 −5 Z 2.39 + 1.55468 × 10 −12 Z 5.35 . (5) Similarly to the free energy, the total pressure can be written as where the ion pressure P i can be decomposed into a noninteracting ('ideal') part and a contribution due to the Coulomb interactions while P e is the pressure of the (uniform) electron background.

OCP in the liquid phase
At temperatures T > T m , ions form a Coulomb liquid. In this case, the non-interacting ('ideal') contribution to the ion free energy is given by (see Eq. (2.71) in Haensel et al. (2007)) where the ion density is the inverse of the Wigner-Seitz cell volume V , n N = 1/V , g s is the spin degeneracy, and λ is the de Broglie wavelength, k B being the Boltzmann constant and the Planck-Dirac constant. Baryon number conservation requires n B = An N . 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.

OCP in the solid phase
Below the crystallization temperature T m , we assume that ions arrange themselves in a perfect body-centred cubic (bcc) lattice (see, e.g. Chamel & Fantina (2016a)).
Since ions can still oscillate about their equilibrium positions, the 'ideal' part of the free energy, Eq. (8), is now replaced by the zero-point motion energy E zp with (an)harmonic corrections (see Sect. 2.3.3 in Haensel et al. (2007)). The ion free energy, Eq. (3), thus becomes where F ii,sol accounts for the Coulomb interactions (static lattice energy, plus thermal and anharmonic corrections), and F pol ie,sol includes the (electric charge) polarization corrections. The zero-point quantum vibration term is given by (Haensel et al. 2007) where u 1 ≡ (ω/ω p ) is a numerical constant (for a bcc crystal, u 1 = 0.511, see Table 2.4 in Haensel et al. (2007)) and the ion plasma frequency ω p is e being the elementary charge. The Coulomb interaction term is given by where the temperature-independent static lattice term reads (Haensel et al. 2007) with C M the Madelung constant (for a bcc lattice, C M = 0.895929, see Table 2.4 in Haensel et al. (2007)) and a N = (4πn N /3) −1/3 is the ion-sphere radius. As for the thermal corrections in the harmonic approximation, F th , and for the anharmonic corrections, F anharm , to the ion vibration, analytical representations have been derived in Baiko et al. (2001) and Potekhin & Chabrier (2010), respectively (see also Appendix C for the complete expressions used in this work). The last term in Eq. (14) accounts for the spin entropy. Although the spin degeneracy remains poorly known for several nuclei, this term has no direct effect on the determination of the melting temperature since it is the same in both the liquid and solid phases. However, the spin entropy might affect the determination of the equilibrium nucleus. Finally, the polarization correction, F pol ie,sol , is given by Eq. (42) in Potekhin & Chabrier (2000) (see also Appendix C).

Multi-component plasma in nuclear statistical equilibrium
While matter at each pressure in the OCP can be described by identical Wigner-Seitz cells, centred on each ion, in the MCP, we expect that different configurations of the Wigner-Seitz cell are realized.

MCP in a liquid phase
Let us consider a very large volume containing different ion species (A (j) , Z (j) ) and, therefore, 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.
The different (A (j) , Z (j) ) configurations are associated with different baryonic densities n (j) B = A (j) /V (j) but share the same total pressure P (see Eq. (6)) imposed by the hydrostatic equilibrium. Moreover, it is supposed that charge neutrality is realized in each cell. This implies that the proton density n p associated with the different components is the same (and equivalent to the electron density n e ), i.e. n e = n p = Z (j) /V (j) .
The total free energy per ion of the system is given by: where the free energy per ion of the component (j), e , accounts for the contribution of the ion and the electrons, including their interactions. We make the hypothesis that this free energy depends only on the characteristics of the component (j), namely (A (j) , Z (j) , V (j) ), and on the global thermodynamic quantities, but it does not depend on the other components (j ′ ) = (j). This assumption is exact at the thermodynamic limit if the different components are associated with macroscopically separated domains. Even in the case of negligible interaction among different ion species, F (j) does not coincide with the free energy of a single Wigner-Seitz cell, as we shall later show. As discussed in Sect. 2.1, we neglect the effect of the nucleon gas.
We can also define the free energy density of the multicomponent system as: where n (j) N is the ion density associated with the cell (j), with j n (j) N A (j) = n B . The ion density is related to the probability p j through or equivalently Ensemble averages are given by: and similar relations hold for the other average quantities. Under the hypothesis of uncorrelated Wigner-Seitz cells (linear mixing approximation), the most probable values for A and Z correspond to those found in the OCP approximation in the same thermodynamic conditions and are denoted by A OCP and Z OCP , respectively. However, the average composition, A and Z , will generally be different due to the co-existence of various nuclear species. Accounting for non-linear mixing effects leads to larger deviations. It is important to note that a first deviation to the linear mixing rule appears due to the translational degree of freedom in the liquid phase (Gulminelli & Raduta (2015)). Indeed, the centre-of-mass position of each ion j of the MCP in the liquid phase is not confined to the single cell volume V (j) but can freely explore the whole volume, leading to Eq. (18) above 3 . Upon replacing this expression in Eq. (8), the single-ion free energy of the MCP in the liquid phase, Eq. (3), becomes: where The extra term on the right hand side of the previous equation is known as the mixing entropy term in the literature, see Medin & Cumming (2010).
Using standard methods in statistical mechanics and following Gulminelli & Raduta (2015); Grams et al. (2018), the probabilities p j and the densities n (j) N are calculated such as to maximize the thermodynamic potential in the canonical ensemble. In view of the chosen decomposition between F i and F e , we have Since the electron part F e of the free energy density does not depend on n (j) N , the variation can be performed on the ion part only: where the single-ion canonical potential is given by: (24) In Eq. (23), the variations dn (j) N are not independent because of the normalization of probabilities, and the baryonic number and charge conservation laws: These constraints are taken into account by introducing Lagrange multipliers (α, µ n , µ p ) leading to the following equations for the equilibrium densities n (j) with N (j) = A (j) − Z (j) . Considering independent variations, the solutions are given by with the normalization The single-ion grand-canonical potentialΩ 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 origin of the rearrangement term, N in Eq. (23) deserves a short discussion. Due to the uniformity of the electron background included in the expression for F e , charge conservation must be realized at the level of each cell: This is at variance with the baryonic density that can fluctuate from cell to cell. In Gulminelli & Raduta (2015); Grams et al. (2018), it was pointed out that this introduces a self-consistency problem. Indeed, the OCP ion free energy F OCP i given by Eqs. (3), (8), and (10) depends on the local cell proton density n (j) p = n p because of the Coulomb interaction, and in turn this implies a dependence on the local density n (j) N through Eq. (32). For this reason, a rearrangement term has to be added to guarantee the thermodynamic consistency of the model. The rearrangement term is calculated using Eqs. (18) and (32) where in the last equality the following approximation has been made: P the pressure in each cell has been taken equal to its average value. This avoids the self-consistency issue due to the dependence of the cell pressure on p j (see Eq. (34) below and Appendix A). R (j) /V (j) can be interpreted as the interaction part of the partial pressure of the (pure-phase) component (j), while, in the MCP, the total pressure reads: We can observe that the partial pressure of the MCP is modified with respect to the pressure defined in the OCP picture, P OCP i = −∂F OCP i /∂V ; in other words, the total pressure of the MCP cannot be calculated via a simple linear mixing rule employing the OCP pressures. The proof of Eq. (34) is given in Appendix A.
To evaluate the MCP composition, with the probability given by Eq. (29), we still have to evaluate the chemical potentials. To this aim, we exploit the thermodynamic relation: where G is the total Gibbs free energy density, n n is the neutron density, and µ e is the electron chemical potential. Using the chemical equilibrium condition µ n = µ p + µ e and the definition of the free energy density, Eq. (17), we get where g ( g e ) is the total (electron) Gibbs free energy per baryon of the MCP, and y p = Z / A is the average proton fraction of the mixture, and g = g i + g e , g i being the ion Gibbs free energy per baryon. The uniformity of the electron density over the different cells n (j) e = n e allows for another representation for the electron chemical potential µ e . We can introduce the Gibbs free energy of each cell: where e ) being the ion (electron) Gibbs free energy in the cell j. The following equalities then hold: where the quantities µ (j) coincide with the respective physical chemical potentials only in the OCP approximation. Dividing Eq. (42) by the cell volume yields where G e is the Gibbs free energy density of electrons. Since G e and F e depend solely on the electron density n e (and are thus the same in each cell), the quantity µ (j) e on the right-hand-side of Eq. (43) must, therefore, coincide with the electron chemical potential, which can be equivalently written as Using Eqs. (36), (37), and (44), we can finally express the single-ion grand-canonical potential in terms of the Gibbs free energies per particle as: In a perturbative treatment of nuclear statistical equilibrium, the average quantities can be replaced with the OCP solution, g ≈ g OCP liq . It is also interesting to express the single-ion grandcanonical potential in terms of the thermodynamic quantities calculated in the OCP approximation. Introducing a OCP single-ion grand-canonical potential as where µ = µ n is the baryonic chemical potential and F can be equivalently written asΩ with the correction term given by

MCP in a solid phase
In the solid state, the equilibrium distribution of ions is given by an equation similar to Eq. (29), but using an appropriate expression for the single-ion grand-canonical potential: Similarly to the liquid state, the single-ion grand canonical potentialΩ (j) i,sol in the solid state can be written as with Ω and δΩ (j) sol is the deviation from linear mixing in the solid phase, see Medin & Cumming (2010). This term contains the rearrangement that can be analytically worked out, but it also couples the probabilities of the m components, and the set of Eqs. (49) should be numerically solved.

Thermodynamic conditions for crystallization
From the thermodynamical point of view, the crystallization temperature at each pressure is, thus, obtained from the Gibbs conditions of phase equilibrium for all ion species.
where F MCP liq(sol) = F i is the ion part of the free energy per ion in the MCP liquid (solid) phase and the partial derivatives should be computed at the equilibrium solutions of each phase given by Eq. (29) and Eq. (49), respectively. Equation (52) has to be supplemented with the extra condition ensuring that the two phases share the same thermodynamic potential: with p p p = (p 1 , . . . , p m ) and the gradient operator ∇ p ∇ p ∇ p has components ∂/∂p j . If the complete set of equations is satisfied by the equilibrium solid and liquid solutions, Eqs. (29) and (49), this means that the two phases can coexist at equilibrium, and crystallization occurs.
An alternative procedure consists of directly solving the Gibbs equilibrium conditions, Eqs. (52), for the unknown fractions p p p sol = (p 1,sol , . . . , p m,sol ) in the solid phase, together with the condition (53). Both procedures are numerically costly. Moreover, they suppose that the crystallization occurs at the thermodynamical transition point. In the case of NS cooling, time scales are such that it is not clear whether nuclear statistical equilibrium is maintained until the transition point, see Goriely et al. (2011). Depending on the dynamics of the process, the ion distribution could be frozen at temperatures larger than the crystallization temperature. In view of these uncertainties, we do not solve the full equations of phase equilibrium. Rather, we consider the much simpler crystallization condition of a OCP: where g OCP sol(liq) is the OCP solution for the Gibbs free energy per baryon in the solid (liquid) phase.
We can see from Eq. (36) that in the case of MCP, the local condition g liq = g sol , where both terms are calculated at the composition p p p liq = (p 1,liq , . . . , p m,liq ) obtained from the MCP equilibrium in the liquid phase, is equivalent to the local equality of the free energy per ion, F MCP sol (p p p liq ) = F MCP liq (p p p liq ). This condition identifies the central zone of the spinodal region in a first order co-existence zone. Since the spinodal is always included inside the binodal, we can expect that our simplified condition for the crystallization transition, Eq. (54), will yield a lower limit estimate for the crystallization temperature.
In our approach for the MCP, it is also possible to calculate the so-called impurity parameter of the solid crust, defined as where p(Z (j) ) is the normalized probability distribution (integrated over all N (j) ) of the element Z (j) and it is assumed that the most abundant species (contributing the most to Z ) form a crystalline structure. This quantity, which also represents the variance of the ionic charge distributions, is important for the calculation of transport coefficients hence also for NS cooling simulations (see, e.g. the discussion in Sect. 9 in Chamel & Haensel (2008) and in Sect. 7 in Meisel et al. (2018) for a review).

Method
We computed the finite-temperature composition of the outer crust of non-accreting unmagnetized NSs, both in the OCP approximation and in the MCP, thus including a distribution of nuclei in nuclear statistical equilibrium, as well as the crystallization temperature for a OCP. We started our calculations at P = 10 −9 MeV fm −3 , which also ensures that the atoms are completely ionized, and repeated the process until the neutron drip sets in, the condition for which is µ n = g = m n c 2 , m n being the neutron mass (see, e.g. Chamel et al. (2015); Pearson et al. (2018) for a recent discussion on the neutron drip). For each value of the pressure, which we increased in steps of ∆P = 0.003P , we determined the composition as follows: (1)  , we performed the calculation of the MCP in the liquid phase, as described in Sect. 3. We went beyond the perturbative approach and computed a self-consistent calculation of the MCP, updating the neutron and proton chemical potentials at each iteration. We found that convergence is reached only after a few additional iterations since the chemical potentials of the OCP are already very close to the self-consistent MCP solution 4 ; (4) We repeated the first three steps, decreasing the temperature until the crystallization temperature, T m , was reached for the OCP.
Step (3) allowed us to calculate the average A and Z , as well as the impurity parameter at T m .
To reduce the computational time, we first estimated the crystallization temperature for the OCP from Eq. (2.28) in Haensel et al. (2007), assuming the Coulomb parameter at melting Γ m = 175 and (A, Z) to be the same as in cold catalyzed matter, for which the composition had been already calculated in Pearson et al. (2018) with the same functional 5 . The density n B in Eq. (56) was estimated from the zerotemperature equation of state (see Table 4 in Pearson et al. (2018)) using a linear interpolation of the pressure. All the results presented in this Section were obtained making use of the experimental masses from AME2016 (Wang et al. 2017) complemented with the HFB-24 nuclear mass model . Unless explicitly stated, we included the following corrections to the free energy: in both the liquid and solid phases, we included the electron exchange and polarization corrections, Eq. (B.5) and Eqs. (C.3) and (C.19), but we dropped the electron correlation energy. For the solid phase, we included the zero-point vibration energy, Eq. (12), as well as the thermal harmonic correction, Eq. (C.14), and the anharmonic corrections, Eq. (C.17).

Crystallization temperature
In Fig. 1 (black solid line), we show the crystallization temperature for the outer crust obtained in the OCP approximation, see Eq. (54). We do not expect that the obtained values of T m will be substantially affected if we replace g OCP liq in Eq. (54) by the average Gibbs energy per baryon in the liquid phase g . Indeed, we verified that the relative differences between g OCP liq and g lie below 0.5%, except at the interface between the outer crust and the inner crust where the deviations become very large. This may be attributed to the neglect of a free nucleon gas which becomes questionable near the neutron drip and at the relative high crystallization temperature (above 2 × 10 9 K). For the considered mass model, the crystallization temperature varies between ≈ 10 8 K and ≈ 2.8 × 10 9 K in the outer crust. These values are in agreement with those presented in the left panel of Fig. 3.17 in Haensel et al. (2007) and obtained 4 The criterion for converge is determined by requiring the difference in the average Gibbs energy per baryon between two consecutive iterations to be below 10 −9 MeV. 5 We actually started the calculations from a value of temperature slightly higher than that given by Eq. (56), thus ensuring that the OCP is in the liquid phase. with the model of Haensel and Pichon (1994) for the outer crust.
The results are quite sensitive to the (even small) corrections included in the free energy. While the inclusion of the exchange correction to the electron energy, Eq. (B.5), has a negligible impact on the determination of T m , including the polarization correction, Eqs. (C.3) and (C.19), changes the crystallization temperature of about a few %, and up to about 40% − 50% around P ≈ 1.25 × 10 −4 MeV fm −3 , where the curve becomes very steep and the composition changes from the liquid 80 Ni to the solid 124 Mo (also see the discussion in Sect. 4.3 and Fig. 5). Concerning the anharmonic correction to the ion vibrations, its inclusion lowers the crystallization temperature in almost all the explored pressure interval, reducing T m up to ≈ 10%. This is shown in Fig. 1, where we plot the crystallization temperature for the OCP with all corrections included (black solid line) or without taking into account either the exchange (red dotted line), the polarization (blue dashed line), or the anharmonic (green dot-dot-dashed line) corrections.

Equilibrium composition of the MCP
The average and most probable values for the mass and charge numbers of ions in a MCP in full equilibrium are plotted in Fig. 2 as a function of the pressure P for two different temperatures: T = 10 9 K (left panel) and T = 2 × 10 9 K (right panel). At these temperatures and pressures, the MCP is in a liquid state. Results obtained in the OCP approximation are also shown for comparison. As expected, the discontinuous changes of composition with pressure found in the OCP approximation are smoothed out when the co-existence of different nuclear species are taken into account. Moreover, the most probable ions are found to coincide with the OCP predictions except for a few values of the pressures, e.g. P ∼ 8 × 10 −7 MeV fm −3 for T = 2 × 10 9 K. This shows that the linear mixing rule is generally a very good approximation in the liquid phase.
The equilibrium composition of the MCP at the crystallization is shown in Fig. 3. The average values for the mass and charge numbers, A and Z , follow the OCP values closely, with two noticeable exceptions around P 1 ≈ 4.2 × 10 −7 MeV fm −3 and P 2 ≈ 1.2 × 10 −4 MeV fm −3 . The deviations appear more clearly as spikes in the pressure variations of the Coulomb coupling parameter at melting, Γ m , displayed in Fig. 4, as calculated using Eq. (C.2) with T = T m (solid line) and Γ m = 175 (horizontal dashed line, see Haensel et al. (2007)). The two pressures P 1 and P 2 signal changes of compositions associated with supercooling in the OCP approximation: the liquid phase of the newly formed ionic species turns out to be unstable, the equilibrium state of those species corresponding to the solid phase for the same pressure. This is illustrated in Fig. 5, where the variations with pressure of the Gibbs free energy per baryon (with respect to the neutron mass) around P 1 and P 2 are plotted. As shown in panel (a) for P P 1 the OCP made of 66 Ni crystallizes when the temperature decreases to T m ≈ 3.3 × 10 8 K, before 66 Ni could be converted into 86 Kr (this would occur at the lower temperature ≈ 3 × 10 8 K if 66 Ni remained liquid). On the contrary, for a slightly higher pressure, the composition of the liquid changes from 66 Ni to 86 Kr at T ≈ 3.9 × 10 8 K before 66 Ni crystallizes, as shown in panel (b). However, the liquid made of 86 Kr at this temperature is supercooled, the solid phase of 86 Kr having a lower Gibbs free energy per baryon. A similar behaviour can be inferred around P 2 for 80 Ni and 124 Mo, as shown in panels (c) and (d) of Fig. 5. However, such supercooling instabilities are the direct consequence of the OCP and are, therefore, spurious. They would disappear in the MCP approach. Except for the two pressures P 1 and P 2 , the Coulomb coupling parameter at melting varies from ≈ 155 to ≈ 180 over almost all the explored range of pressures, in fairly good agreement with the canonical value Γ m = 175. The crystallization temperature can thus be well estimated by Eq. (56). The abrupt changes in composition found in the OCP approximation at pressures P 1 and P 2 disappear in the MCP approach. This is best seen in Fig. 6, where the normalized probability distribution p(Z) is plotted for temperatures close to the crystallization temperature and for two different pressures around P 1 (left panel) and P 2 (right panel). The distribution exhibits a bimodal character around 66 Ni and 86 Kr in the former case, and around 80 Ni and 124 Mo in the latter case, leading to a gradual change of the most probable nuclide from one to the other as the pressure is increased. For this reason, the change in the most probable nucleus in the MCP is shifted to a slightly higher pressure with respect to the OCP case, see Fig. 3. Moreover, despite the apparent discontinuity in the most probable nucleus, the composition actually varies very smoothly, as can be seen from the average values of the mass and charge numbers.
To better assess the validity of the OCP approximation, we plot in Fig. 7 the normalized probability distribution p(Z) as a function of Z, for two different pressures, P = 10 −5 MeV fm −3 (left panels) and P = 2 × 10 −4 MeV fm −3 (right panels), and for two different temperatures, T ≈ T m (upper panels) and T m < T = 5 × 10 9 K (lower panels). The charge numbers Z OCP predicted in the OCP approximation are indicated by arrows. As can be seen, Z OCP coincides with the most probable Z , thus indicating that deviations from the linear mixing rule are negligibly small. At relatively low pressure and temperature (panel (a)), the OCP treatment is a very good approximation since the distribution is very peaked around the most thermodynamically favoured nuclide. With increasing pressure and temperature, the broadening of the distribution makes the OCP approximation less accurate. In particular, panel (d) shows that for some pressure and temperature the distribution may even become bimodal (a similar situation is also displayed in Fig. 6).

Impurity parameter
We show, in Fig. 8, the impurity parameter, Eq. (55), as a function of pressure at the crystallization temperature T m (solid line). These data are available in tabular format at the CDS. Since the impurity parameter represents the variance of the charge distribution, low values of Q imp (say below 1) indicate that the distribution is quite peaked and, therefore, the OCP treatment is a good approximation. This is in accordance with Fig. 7, where it can be seen that low values of Q imp correspond to pressures for which A and Z are very close or nearly coincide with A OCP and Z OCP , respectively. On the contrary, appreciable deviations from the OCP predictions translate into large values for Q imp , reaching, for P ≈ 1.2 − 1.3 × 10 −4 MeV fm −3 , about 50 at crystallization (see also panel (b) of Fig. 6). The variations of Q imp with pressure suggest that the outer crust may actually consist of an alternation of pure (highly conductive) and impure (highly resistive) layers.
This calculation was performed based on the hypothesis that the statistical equilibrium is maintained during the cooling process down to the crystallization temperature. However, if the interior of a NS cools down rapidly enough in comparison to the various reaction rates, the composition may be frozen at some finite temperature T f > T m , see e.g. Goriely et al. (2011) (see also Haensel et al. (2007); Chamel & Haensel (2008)). A realistic calculation of T f requires dynamical simulations and is left for future works. For comparison, we show in Fig. 8 the impurity parameter, assuming that the composition is frozen at a fixed temperature of T = 10 9 K (dashed line). The most prominent deviations are seen in the shallowest layers of the crust, where the differences between T f and T m are the largest.

Conclusions
In this work, we studied the cooling and the equilibrium composition of the outer layers of a non-accreting unmagnetized NS down to crystallization. To this end, we took into account the co-existence of different nuclear species in a self-consistent nuclear statistical equilibrium treatment using the latest experimental atomic mass data supplemented with the microscopic nuclear mass table HFB-24. We calculated the crystallization temperature in the OCP approximation for the range of pressures relevant for the outer crust, starting from P = 10 −9 MeV fm −3 . We found that the crystallization temperature varies from ≈ 10 8 K to ≈ 2.8 × 10 9 K. The corresponding Coulomb coupling parameter at melting is found to be reasonably close to the canonical value of 175, except for specific values of the pressure for which supercooling occurs.
As for the composition, the discontinuous behaviour with pressure observed in the OCP approximation is smoothed out when matter is modelled according to a MCP approach. However, the average and most probable values for the mass and charge numbers follow the OCP predictions closely at the crystallization temperature, except when supercooling occurs in the OCP approximation. This confirms that the linear mixing rule usually adopted in the description of the liquid phase is generally a very good approximation, as long as the thermodynamical equilibrium is maintained during the NS cooling, down to the crystallization temperature.
Within our approach for the MCP, we also consistently calculated the impurity parameter in the range of pressure of interest for the outer crust. The non-monotonic variations of Q imp , whose values can change by several orders of magnitude, amounting up to about 50 at crystallization, suggests that the crust may be composed of an alternation of pure (highly conductive) and impure (highly resistive) layers. In the scenario where a NS cools down sufficiently rapidly and the composition is frozen at some finite temperature T f higher than the crystallization temperature T m , the impurity parameter may be significantly larger than that obtained at T m , especially in the shallowest layer of the crust where the deviations between T f and T m are expected to be the largest. Therefore, the results that we obtained for Q imp at crystallization can be considered a lower limit. The precise determination of T f (hence, of the impurity parameter as well) would require dynamical simulations with a nuclear reaction network and is left for future studies. The results we obtained are based on the same nuclear energy-density functional BSk24 for which unified equations of state of non-accreting NSs have been recently calculated and can be directly implemented in NS cooling simulations.
In this work, we applied our treatment for the MCP in the outer layers of the NS, however, a similar approach can be also employed in the deeper layers of the NS with a proper account of the free nucleon gas. This improvement deserves further investigation in view of the significance of the presence of impurities for the evolution of NSs.   Fig. 8. Impurity parameter Q imp as a function of pressure inside the outer crust assuming that the composition is frozen at the crystallization temperature T m (black solid line) and at a fixed temperature of T = 10 9 K (red dashed line). See text for details.
Appendix A: Pressure of the multi-component plasma The pressure of the MCP is more easily worked out if we consider the canonical ensemble, where P i (P e ) is the ion (electron) pressure, T is the temperature and {p j } is the set of probabilities of the different ion species, p j being the probability of the component j characterized by an ion with mass (charge) A (j) (Z (j) ). In Eq. (A.1), the total free energy F = F e + F i (electron plus ion part), volume V , baryonic number A, and charge Z, are calculated per ion 6 . Specifically, the volume entering Eq. (A.1) is the average volume per ion: The ion and electron parts of the free energy are given by where F (j) is the free energy per ion of the component (j) as given for the liquid phase by Eq. (21): is the ion mass and n (j) N is given in Eq. (18). Omitting for simplicity the constant variables (p j , T ) in the derivatives, Eq. (A.1) can be written as Since the electron density n e is the same in each cell, the derivative of the electron free energy yields directly the electron pressure − ∂F e ∂V = n 2 e ∂f e ∂n e ≡ P e . (A.6) As for the ion contribution, we consider separately the ideal part (second term in Eq. (A.4)) and the interaction part (last term in Eq. (A.4)). Using the definition of the partial density, Eq. (18), n (j) The ionic pressure becomes: As in the main text, also in the Appendices we use capital letters for the energy per ion, e.g F for the free energy per ion, small letters for the (free) energy per baryon, e.g. f , and the notation F for the free energy density.
where we have used n B = A/V , n B being the baryon density, with A = j p j A (j) . Making use of the charge conservation, n e = n B Z/A, and considering that in the canonical ensemble the derivatives are evaluated for fixed numbers of particles, where P (j),int i is the interaction part of the pressure as calculated in the (pure phase) OCP approximation: (A.10) In the case of a MCP, we can still define the partial pressure of the (pure) (j) component as but the total pressure in a MCP is not just the sum of the pressures of the (pure) OCP phases. Rather, it is given by with P Appendix B: Free energy and pressure of the electron gas For completeness, we give the expressions for the free energy and pressure of the (uniform) electron gas at finite temperature. The former can be written as where the first term denotes the kinetic ('ideal') contribution (without the rest-mass energy), F exc e is the exchange part, and F corr e accounts for the electron-correlation free energy. The last term is the rest-mass energy, m e being the electron mass. We note that the correction due to the polarization is not included here since it is explicitly included in F ie , and accounted for in the ion free energy, Eq. (3).
The kinetic free energy density, without the rest mass energy, is given by (see, e.g. Chap. 24 of Weiss et al. (2004) and Sect. 2 in Lattimer (1996) We note that with respect to the expression for F kin e given by Eq. (2.65) in Haensel et al. (2007) there are two differences: (i) the term −8x 3 is not present in Eq. (2.65) in Haensel et al. (2007) because the latter equation includes the rest-mass energy while our Eq. (B.2) does not; (ii) the finite-temperature corrections are not the same. This second discrepancy comes from a different expansion of the integrals at finite temperature. Therefore, also the temperature corrections in the pressure are different, see our Eq. (B.11) and Eq. (2.67) in Haensel et al. (2007). where λ e = /(m e c) is the electron Compton wavelength, x r = p F m e c = (3π 2 n e ) 1/3 m e c (B.3) is the relativity parameter, p F being the Fermi momentum, and g(x) = −8x 3 r + 3x r (1 + 2x 2 r ) 1 + x 2 r − 3 sinh −1 (x r ) . (B.4) The exchange correction to the free energy density for a strongly degenerate electron system is given by (see Eq. (2.151) in Haensel et al. (2007); see also Stolzmann & Blöcker (1996)) where the kinetic term reads (Weiss et al. 2004) P kin e = m e c 2 24π 2 λ 3 e x r 1 + x 2 r (2x 2 r − 3) +3 ln(x r + 1 + x 2 r ) + 4π 2 (k B T ) 2 (m e c 2 ) 2 x r 1 + x 2 r . (B.11) The exchange term can be written as (Stolzmann & Blöcker 1996)  The correlation correction to the pressure being negligible, as for the free energy density, we neglect it here.

Appendix C: Free energy of the Coulomb plasma of ions
For the completeness and reproducibility of the results, here we report the expressions for the free energy of the Coulomb plasma of ions that we have used in this work.

C.1. Coulomb liquid
In the liquid phase, the ion free energy in the OCP approximation is given by Eq.
8 Note that in the second line of Eq. (16) in Potekhin & Chabrier (2000), B2 ln(1 + Γ/B1) should be replaced by B2 ln(1 + Γ/B2). The correct expression is given by Eq. (2.87) in Haensel et al. (2007) and implemented in the FITION9 routine available on the Ioffe website http://www.ioffe.ru/astro/EIP/index.html. with T p = ω p /k B . Note that for finite values of Z, and assuming Γ p ≫ 1, the electron polarization correction to the energy for a bcc lattice at zero temperature can be approximately expressed as (Chamel & Fantina 2016b) where the Thomas-Fermi correction is given by Salpeter (1961)  with E L the static lattice term given by Eq. (15). The corresponding pressure terms can be derived from the thermodynamic definition, Eq. (A.1). The routines that compute the analytical representations of both the free energy and pressure of Eqs. (C.1), (C.3), (C.14), and (C.19) are available on the Ioffe Institute website 10 .