Issue 
A&A
Volume 633, January 2020



Article Number  A149  
Number of page(s)  13  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/201936359  
Published online  24 January 2020 
Crystallization of the outer crust of a nonaccreting neutron star^{⋆}
^{1}
Grand Accélérateur National d’Ions Lourds (GANIL), CEA/DRF – CNRS/IN2P3, Boulevard Henri Becquerel, 14076 Caen, France
email: anthea.fantina@ganil.fr
^{2}
Institut d’Astronomie et d’Astrophysique, Université Libre de Bruxelles, CP226, 1050 Brussels, Belgium
^{3}
LPC (CNRS/ENSICAEN/Université de Caen Normandie), UMR6534, 14050 Caen Cedex, France
Received:
23
July
2019
Accepted:
25
October
2019
Context. The interior of a neutron star is usually assumed to be made of cold catalyzed matter. However, the outer layers are unlikely to remain in full thermodynamic equilibrium during the formation of the star and its subsequent cooling, especially after crystallization occurs.
Aims. We study the cooling and the equilibrium composition of the outer layers of a nonaccreting neutron star down to crystallization. Here the impurity parameter, generally taken as a free parameter in cooling simulations, is calculated selfconsistently using a microscopic nuclear model for which a unified equation of state has recently been determined.
Methods. We follow the evolution of the nuclear distributions of the multicomponent Coulomb liquid plasma fully selfconsistently, adapting a general formalism originally developed for the description of supernova cores. We calculate the impurity parameter at the crystallization temperature as determined in the onecomponent plasma approximation.
Results. Our analysis shows that the sharp changes in composition obtained in the onecomponent plasma 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 values of the pressure for which supercooling occurs in the onecomponent plasma approximation. Our multicomponent treatment leads to nonmonotonic 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 within the same nuclear model.
Conclusions. Our selfconsistent approach to hot dense multicomponent plasma shows that the presence of impurities in the outer crust of a neutron star is nonnegligible and may have a sizeable impact on transport properties. In turn, this may have important implications not only for the cooling of neutron stars, but also for their magnetorotational evolution.
Key words: stars: neutron / dense matter / nuclear reactions, nucleosynthesis, abundances / plasmas
The table of the impurity parameter at the crystallization temperature shown in Fig. 8 is only available at the CDS via anonymous ftp to cdsarc.ustrasbg.fr (130.79.128.5) or via http://cdsarc.ustrasbg.fr/vizbin/cat/J/A+A/633/A149
© A. F. Fantina et al. 2020
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
Formed in the aftermath of gravitational corecollapse 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 multicomponent 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 socalled “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 onecomponent Coulomb crystal (except, possibly, at the 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 multicomponent 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 coexistence 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 nonaccreting unmagnetized NS. After determining the crystallization temperature in the onecomponent plasma (OCP) approximation, the nuclear distributions and the impurity parameter are calculated fully selfconsistently, 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 Sects. 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.
2. Onecomponent Coulomb plasma
2.1. 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 highlydegenerate 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 (singlenucleus) 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 (noninteracting) 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 Chapt. 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, is the noninteracting (“ideal”) contribution to the ion free energy, and accounts for interactions. Specifically, , where F_{ii} includes all the Coulomb interactions (between ions, between electrons, and between ions and the uniform electron gas) and 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 HFB24 theoretical mass table based on the nuclear energydensity functional theory^{2} (Goriely et al. 2013). The underlying functional has been recently used to determine the groundstate composition and the equation of state in all regions of a nonaccreting 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)
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.
2.2. OCP in the liquid phase
At temperatures T > T_{m}, ions form a Coulomb liquid. In this case, the noninteracting (“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 PlanckDirac 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.
2.3. OCP in the solid phase
Below the crystallization temperature T_{m}, we assume that ions arrange themselves in a perfect bodycentred 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 zeropoint 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 includes the (electric charge) polarization corrections. The zeropoint 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 temperatureindependent 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 ionsphere 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, , is given by Eq. (42) in Potekhin & Chabrier (2000) (see also Appendix C).
3. Multicomponent 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.
3.1. 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 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), , 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 is the ion density associated with the cell (j), with . 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 coexistence of various nuclear species. Accounting for nonlinear 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 centreofmass 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 singleion free energy of the MCP in the liquid phase, Eq. (3), becomes:
where M^{′(j)} = M′(A^{(j)}, Z^{(j)}) and , as given by Eqs. (3), (8), and (10) (in Eq. (8), n_{N} = 1/V^{(j)}). 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 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 ℱ_{e} of the free energy density does not depend on , the variation can be performed on the ion part only:
where the singleion canonical potential is given by:
In Eq. (23), the variations 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 :
with N^{(j)} = A^{(j)} − Z^{(j)}. Considering independent variations, the solutions are given by
with the normalization
The singleion grandcanonical potential reads:
where μ_{n} and μ_{p} can be identified with the neutron and proton chemical potentials, respectively. In the definitions above, the ion free energy contains the restmass energy, thus the chemical potentials include the restmass energies as well.
The origin of the rearrangement term, in Eq. (23) deserves a short discussion. Due to the uniformity of the electron background included in the expression for ℱ_{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) and Grams et al. (2018), it was pointed out that this introduces a selfconsistency problem. Indeed, the OCP ion free energy given by Eqs. (3), (8), and (10) depends on the local cell proton density because of the Coulomb interaction, and in turn this implies a dependence on the local density 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) (see Eqs. (15), (22), and (23) in Grams et al. 2018):
where in the last equality the following approximation has been made: , i.e. the pressure in each cell has been taken equal to its average value. This avoids the selfconsistency issue due to the dependence of the cell pressure on p_{j} (see Eq. (34) below and Appendix A). ℛ^{(j)}/V^{(j)} can be interpreted as the interaction part of the partial pressure of the (purephase) 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, ; 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 𝒢 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 allows for another representation for the electron chemical potential μ_{e}. We can introduce the Gibbs free energy of each cell:
where
() 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 𝒢_{e} is the Gibbs free energy density of electrons. Since 𝒢_{e} and ℱ_{e} depend solely on the electron density n_{e} (and are thus the same in each cell), the quantity on the righthandside 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 singleion grandcanonical 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, .
It is also interesting to express the singleion grandcanonical potential in terms of the thermodynamic quantities calculated in the OCP approximation. Introducing a OCP singleion grandcanonical potential as
where μ = μ_{n} is the baryonic chemical potential and is given by Eq. (3), can be equivalently written as
with the correction term given by
3.2. 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 singleion grandcanonical potential:
Similarly to the liquid state, the singleion grand canonical potential in the solid state can be written as
with
and 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 Eq. (49) should be numerically solved.
3.3. 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. For a system of m components, these conditions correspond to a set of m − 1 coupled equations (Medin & Cumming 2010):
where 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 Eqs. (29) and (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_{1},⋯,p_{m}) and the gradient operator ∇_{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, Eq. (52), for the unknown fractions 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 higher 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 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_{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, . This condition identifies the central zone of the spinodal region in a first order coexistence 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 socalled 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).
4. Numerical results
4.1. Method
We computed the finitetemperature composition of the outer crust of nonaccreting 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 & Fantina 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)
Starting from a highenough temperature for the plasma to be in a liquid phase, we first minimized the Gibbs free energy per baryon in the OCP approximation, (see Sect. 2.2), thus yielding and the corresponding neutron and proton chemical potentials, and ;

(2)
For the same nucleus , we calculated the Gibbs free energy per baryon of the solid phase, (see Sect. 2.3), and we checked whether crystallization had occurred for the OCP, that is, whether ; see Eq. (54);

(3)
Starting from the OCP solution, that is, from and , 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 selfconsistent 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 selfconsistent 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 HFB24 nuclear mass model (Goriely et al. 2013). 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, Eqs. (B.5), (C.3), and (C.19), but we dropped the electron correlation energy. For the solid phase, we included the zeropoint vibration energy, Eq. (12), as well as the thermal harmonic correction, Eq. (C.14), and the anharmonic corrections, Eq. (C.17).
4.2. 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 in Eq. (54) by the average Gibbs energy per baryon in the liquid phase ⟨g⟩. Indeed, we verified that the relative differences between 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 with the model of Haensel & Pichon (1994) for the outer crust.
Fig. 1. Crystallization temperature for the onecomponent plasma (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 dotdotdashed line) corrections. The inset shows a zoom in the highpressure regime. See text for details. 
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 dotdotdashed line) corrections.
4.3. 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 coexistence of different nuclear species is 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.
Fig. 2. Variation with pressure P of the average (solid lines) and most probable (dotted lines) values of the charge number Z and mass number A of ions in a multicomponent liquid plasma in full equilibrium for two selected temperatures: T = 10^{9} K (panel a) and T = 2 × 10^{9} K (panel b). For comparison, results obtained in the onecomponent plasma approximation (dashed lines) are also shown. See text for details. 
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.
Fig. 4. Coulomb parameter at melting, Γ_{m}, as a function of pressure. The dashed horizontal line indicates the value of Γ_{m} = 175. 
Fig. 5. Gibbs free energy per baryon with respect to the neutron mass energy as a function of temperature for the liquid (dashed lines) and solid (solid lines) phase of different nuclei for different pressures (labelled in units of MeV fm^{−3}). Panels a and b: ^{66}Ni (black thin lines) and ^{86}Kr (red thick lines); panels c and d: ^{80}Ni (black thin lines) and ^{124}Mo (red thick lines). See text for details. 
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).
Fig. 6. Normalized probability distribution p(Z) as a function of Z for two selected temperatures and pressures. Panel a: T = 0.4 × 10^{9} K, P = 4.2 × 10^{−7} MeV fm^{−3} (black circles) and P = 4.4 × 10^{−7} MeV fm^{−3} (blue squares). Panel b: T = 1.3 × 10^{9} K, P = 1.19 × 10^{−4} MeV fm^{−3} (black circles) and P = 1.22 × 10^{−4} MeV fm^{−3} (blue squares). See text for details. 
4.4. 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.
Fig. 7. Normalized probability distribution of ions p(Z) as a function of the charge number Z, for pressures P = 10^{−5} MeV fm^{−3} (panels a and c) and P = 2 × 10^{−4} MeV fm^{−3} (panels b and d). Upper (lower) panels: distributions of the multicomponent plasma in the liquid phase around (above) the crystallization temperature. Arrows indicate the charge number Z^{OCP} predicted by the onecomponent plasma approximation. See text for details. 
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.
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. 
5. Conclusions
In this work, we studied the cooling and the equilibrium composition of the outer layers of a nonaccreting unmagnetized NS down to crystallization. To this end, we took into account the coexistence of different nuclear species in a selfconsistent nuclear statistical equilibrium treatment using the latest experimental atomic mass data supplemented with the microscopic nuclear mass table HFB24. 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 nonmonotonic 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 energydensity functional BSk24 for which a unified equation of state of nonaccreting NSs has 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.
The mass table is available on the BRUSLIB online database http://www.astro.ulb.ac.be/bruslib/ (Xu et al. 2013).
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.
We note that with respect to the expression for given by Eq. (2.65) in Haensel et al. (2007) there are two differences: (i) the term −8x^{3} of Eq. (B.4) is not present in Eq. (2.65) in Haensel et al. (2007) because the latter equation includes the restmass energy while our Eq. (B.2) does not; (ii) the finitetemperature 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).
Note that in the second line of Eq. (16) in Potekhin & Chabrier (2000), B_{2}ln(1 + Γ/B_{1}) should be replaced by B_{2}ln(1 + Γ/B_{2}). 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 respect to Potekhin & Chabrier (2010), we have indicated d_{1} instead of b_{1}, and f_{n} instead of a_{n} to avoid conflicting notation for the numerical coefficients with previous expressions of the thermal (harmonic) term.
http://www.ioffe.ru/astro/EIP/index.html. We have employed here the routines for unmagnetized plasmas.
Acknowledgments
The research leading to these results has received funding from the CNRS PICS07889; this work was also partially supported by the PHAROS European Cooperation in Science and Technology (COST) action CA16214. The work of N.C. was supported by Fonds de la Recherche Scientifique (Belgium) under grant IISN 4.4502.19. The authors would like to thank A. Y. Potekhin for valuable discussions.
References
 Baiko, D. A., Potekhin, A. Y., & Yakovlev, D. G. 2001, Phys. Rev. E, 64, 057402 [NASA ADS] [CrossRef] [Google Scholar]
 Baym, G., Pethick, C., & Sutherland, P. 1971, ApJ, 170, 299 [Google Scholar]
 Chamel, N., & Fantina, A. F. 2016a, Phys. Rev. C, 94, 065802 [NASA ADS] [CrossRef] [Google Scholar]
 Chamel, N., & Fantina, A. F. 2016b, Phys. Rev. D, 93, 063001 [NASA ADS] [CrossRef] [Google Scholar]
 Chamel, N., & Haensel, P. 2008, Liv. Rev. Relat., 11, 10 [Google Scholar]
 Chamel, N., Fantina, A. F., & Zdunik, J.L., & Haensel, P., 2015, Phys. Rev. C, 91, 055803 [NASA ADS] [CrossRef] [Google Scholar]
 Farouki, R. T., & Hamaguchi, S. 1993, Phys. Rev. E, 47, 4330 [NASA ADS] [CrossRef] [Google Scholar]
 Goriely, S., Chamel, N., Janka, H.T., & Pearson, J. M. 2011, A&A, 531, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Goriely, S., Chamel, N., & Pearson, J. M. 2013, Phys. Rev. C, 88, 024308 [NASA ADS] [CrossRef] [Google Scholar]
 Gourgouliatos, K. N., & Esposito, P. 2018, in The Physics and Astrophysics of Neutron Stars, eds. L. Rezzolla, P. Pizzochero, D. I. Jones, N. Rea, & I. Vidaña (Berlin: Springer), Astrophys. Space Sci. Lib., 457, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Grams, G., Giraud, S., Fantina, A. F., & Gulminelli, F. 2018, Phys. Rev. C, 97, 035807 [CrossRef] [Google Scholar]
 Gulminelli, F., & Raduta, Ad. R. 2015, Phys. Rev. C, 92, 055803 [NASA ADS] [CrossRef] [Google Scholar]
 Haensel, P., & Pichon, B. 1994, A&A, 283, 313 [NASA ADS] [Google Scholar]
 Haensel, P., Potekhin, A. Y., & Yakovlev, D. G. 2007, Neutron Stars 1. Equation of state and structure (New York: Springer) [Google Scholar]
 Lattimer, J. M. 1996, in The Nuclear Equation of State and Supernovae, in Nuclear Equation of State, eds. A. Ansari, & L. Satpathy (Singapore: World Scientific), 83 [Google Scholar]
 Lunney, D., Pearson, J. M., & Thibault, C. 2003, Rev. Mod. Phys., 75, 1021 [NASA ADS] [CrossRef] [Google Scholar]
 Medin, Z., & Cumming, A. 2010, Phys. Rev. E, 81, 036107 [NASA ADS] [CrossRef] [Google Scholar]
 Meisel, Z., Deibel, A., Keek, L., Shternin, P., & Elfritz, J. 2018, J. Phys. G, 45, 093001 [Google Scholar]
 Pearson, J. M., Goriely, S., & Chamel, N. 2011, Phys. Rev. C, 83, 065810 [NASA ADS] [CrossRef] [Google Scholar]
 Pearson, J. M., Chamel, N., Potekhin, A. Y., et al. 2018, MNRAS, 481, 2994 [NASA ADS] [Google Scholar]
 Pons, J. A., Viganò, D., & Rea, N. 2013, Nat. Phys., 9, 431 [NASA ADS] [CrossRef] [Google Scholar]
 Potekhin, A. Y., & Chabrier, G. 2000, Phys. Rev. E, 62, 8554 [NASA ADS] [CrossRef] [Google Scholar]
 Potekhin, A. Y., & Chabrier, G. 2010, Contrib. Plasma Phys., 50, 82 [NASA ADS] [CrossRef] [Google Scholar]
 Salpeter, E. E. 1961, ApJ, 134, 669 [NASA ADS] [CrossRef] [Google Scholar]
 Schmitt, A., & Shternin, P. 2018, in The Physics and Astrophysics of Neutron Stars, eds. L. Rezzolla, P. Pizzochero, D. I. Jones, N. Rea, & I. Vidaña (Berlin: Springer), Astrophys. Space Sci. Lib., 457, 455 [NASA ADS] [CrossRef] [Google Scholar]
 Stolzmann, W., & Blöcker, T. 1996, A&A, 314, 1024 [NASA ADS] [Google Scholar]
 Tondeur, F. 1971, A&A, 14, 451 [NASA ADS] [Google Scholar]
 Wang, M., Audi, G., Kondev, F. G., Huang, W. J., Naimi, S., & Xu, X. 2017, Chin. Phys. C, 41, 030003 [Google Scholar]
 Weiss, A., Hillebrandt, W., Thomas, H. C., & Ritter, H. 2004, Cox and Giuli’s Principles of Stellar Structure, 2nd edn. (Cambridge: Cambridge Scientific Publishers) [Google Scholar]
 Xu, Y., Goriely, S., Jorissen, A., Chen, G. L., & Arnould, M. 2013, A&A, 549, A106 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Appendix A: Pressure of the multicomponent 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):
where M^{′(j)} is the ion mass and 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
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), , we have
The ionic pressure becomes:
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 is the interaction part of the pressure as calculated in the (pure phase) OCP approximation:
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 calculated as in Eq. (A.10).
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 restmass energy), is the exchange part, and accounts for the electroncorrelation free energy. The last term is the restmass 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 , 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. Chapt. 24 of Weiss et al. 2004 and Sect. 2 in Lattimer 1996)^{7}
where λ_{e} = ℏ/(m_{e}c) is the electron Compton wavelength,
is the relativity parameter, p_{F} being the Fermi momentum, and
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
and C_{exc} = −0.7046. As for the correlation energy, since it is expected to be negligible, especially in the relativistic regime (see, e.g. the discussion in Pearson et al. 2011 and in Sect. 2.4.3 in Haensel et al. 2007), we neglect it here.
The pressure can be similarly decomposed as
where the kinetic term reads (Weiss et al. 2004)
The exchange term can be written as (Stolzmann & Blöcker 1996)
where is given by Eq. (B.5) and the Gibbs free energy density is expressible as (see Eqs. (49)–(51) in Stolzmann & Blöcker 1996)
with
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. (3), with the “ideal” and interaction parts given by Eqs. (8) and (10), respectively. The analytical representation of the total Coulomb contribution, F_{ii, liq}, has been derived by Potekhin & Chabrier (2000)^{8}:
where A_{1}, A_{2}, , B_{1}, B_{2}, B_{3}, and B_{4} are numerical constants, and Γ is the Coulomb parameter,
a_{N} = (4π/3 n_{e}/Z)^{−1/3} being the interion spacing.
As for the polarization correction to the free energy, , an analytical fit is given by Eq. (19) in Potekhin & Chabrier (2000):
where r_{s} ≡ a_{e}/a_{0} is the density parameter with a_{e} = (4πn_{e}/3)^{−1/3} the electronsphere radius and a_{0} = ℏ^{2}/(m_{e}e^{2}) the Bohr radius, Γ_{e} is the coupling parameter for nondegenerate electrons,
and
and
C.2. Coulomb crystal
For a Coulomb crystal, the free energy in the OCP is given by Eq. (11), with Eqs. (12) and (14). Analytical expressions for the thermal contribution due to the ion vibrations around their equilibrium position in the harmonic approximation and the anharmonic correction have been derived by Baiko et al. (2001) and Potekhin & Chabrier (2010), respectively. The analytical fitting formula for the thermal (harmonic) contribution, F_{th}, can be found in Baiko et al. (2001, see their Eq. (13)),
where θ ≡ ℏω_{p}/(k_{B}T) = T_{p}/T, ω_{p} being the ion plasma frequency, Eq. (13), and
with α_{n}, a_{n}, and b_{n} numerical constants (see Table II in Baiko et al. 2001). The anharmonic correction, F_{anharm}, is only known for a bcc lattice. Analytical expressions have been derived in Potekhin & Chabrier (2010); see their Eq. (8):
where
with c_{1}, d_{1}, and f_{n} numerical constants^{9}. In Eq. (C.17) (Eq. (8) of Potekhin & Chabrier 2010), the anharmonic correction for a classical Coulomb crystal derived in Farouki & Hamaguchi (1993), Eq. (C.18), has been modified by the inclusion of two additional terms reproducing the zerotemperature and classical limits. This expression is valid for any value of θ and ensures that the anharmonic corrections to the heat capacity and entropy do not exceed the dominant (harmoniclattice) contribution (Potekhin & Chabrier 2010).
The polarization correction in the solid phase has been analytically fitted in Potekhin & Chabrier (2000) as
where
with α the fine structure constant. The parameters s and b_{1}–b_{4}, that depend on Z only, are given by Potekhin & Chabrier (2000)
For a bcc lattice, a_{1} = 1.1866, a_{2} = 0.684, a_{3} = 17.9, a_{4} = 41.5, and q = 0.205 (see Table III in Potekhin & Chabrier 2000).
In the limit of low temperature, θ ≡ T_{p}/T ≫ 1, for which Q(θ)→qθ, the polarization correction to the free energy density reduces to (see Appendix B in Pearson et al. 2018)
where
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 ThomasFermi 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}.
All Figures
Fig. 1. Crystallization temperature for the onecomponent plasma (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 dotdotdashed line) corrections. The inset shows a zoom in the highpressure regime. See text for details. 

In the text 
Fig. 2. Variation with pressure P of the average (solid lines) and most probable (dotted lines) values of the charge number Z and mass number A of ions in a multicomponent liquid plasma in full equilibrium for two selected temperatures: T = 10^{9} K (panel a) and T = 2 × 10^{9} K (panel b). For comparison, results obtained in the onecomponent plasma approximation (dashed lines) are also shown. See text for details. 

In the text 
Fig. 3. Same as Fig. 2 at the crystallization temperature T_{m}. 

In the text 
Fig. 4. Coulomb parameter at melting, Γ_{m}, as a function of pressure. The dashed horizontal line indicates the value of Γ_{m} = 175. 

In the text 
Fig. 5. Gibbs free energy per baryon with respect to the neutron mass energy as a function of temperature for the liquid (dashed lines) and solid (solid lines) phase of different nuclei for different pressures (labelled in units of MeV fm^{−3}). Panels a and b: ^{66}Ni (black thin lines) and ^{86}Kr (red thick lines); panels c and d: ^{80}Ni (black thin lines) and ^{124}Mo (red thick lines). See text for details. 

In the text 
Fig. 6. Normalized probability distribution p(Z) as a function of Z for two selected temperatures and pressures. Panel a: T = 0.4 × 10^{9} K, P = 4.2 × 10^{−7} MeV fm^{−3} (black circles) and P = 4.4 × 10^{−7} MeV fm^{−3} (blue squares). Panel b: T = 1.3 × 10^{9} K, P = 1.19 × 10^{−4} MeV fm^{−3} (black circles) and P = 1.22 × 10^{−4} MeV fm^{−3} (blue squares). See text for details. 

In the text 
Fig. 7. Normalized probability distribution of ions p(Z) as a function of the charge number Z, for pressures P = 10^{−5} MeV fm^{−3} (panels a and c) and P = 2 × 10^{−4} MeV fm^{−3} (panels b and d). Upper (lower) panels: distributions of the multicomponent plasma in the liquid phase around (above) the crystallization temperature. Arrows indicate the charge number Z^{OCP} predicted by the onecomponent plasma approximation. See text for details. 

In the text 
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. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.