© The Authors 2023

1. Introduction

Formed from core-collapse supernova explosions, neutron stars (NSs) are initially born hot, with temperatures exceeding a few times 10¹⁰ K ≈ 1 MeV (Prakash et al. 1997; Pons et al. 1999; Haensel et al. 2007). At such temperatures, the crust of the proto-neutron star (PNS) is expected to be made of a Coulomb liquid composed of different nuclear species in a background of electrons and nucleons (see e.g. Oertel et al. 2017 for a review). It is generally assumed that, as the NS crust cools down, this multi-component plasma (MCP) remains in full thermodynamic equilibrium until the ground state at zero temperature is reached. Under this so-called ‘cold catalysed matter’ hypothesis, the NS crust is considered to be made of pure layers, each consisting of a one-component Coulomb crystal. However, if the NS cools down rapidly enough, the composition of the crust might be frozen at a temperature higher than the crystallisation temperature (see e.g. Goriely et al. 2011), and the one-component picture becomes less reliable (see e.g. Fantina et al. 2020; Carreau et al. 2020a). Moreover, a correct description of the finite-temperature PNS crust requires the computation of a full nuclear ensemble. However, due to the complexity of calculating the nuclear distribution of MCP, most approaches to modelling the finite-temperature crust have been performed within the framework of the so-called one-component plasma (OCP) approximation. The latter assumes that, at a given thermodynamic condition, the ensemble of nuclei is represented by only one kind of nucleus, the one that is thermodynamically more favourable (see e.g. Onsi et al. 2008; Avancini et al. 2009, 2017; Fantina et al. 2020; Carreau et al. 2020a,b; Dinh Thi et al. 2023). This approximation is justified at relatively low densities and temperatures, where the most probable nucleus is very close to the ground-state OCP composition and the contribution of other ions is typically very small, or when computing average thermodynamic quantities (Burrows & Lattimer 1984). On the other hand, the coexistence of different nuclear species can have an important impact on transport properties as well as on the magneto-rotational evolution of NSs (Jones 2004; Pons et al. 2013; see also Schmitt & Shternin 2018; Gourgouliatos et al. 2018 for a review). For these reasons, (P)NS cooling simulations are usually performed using the ground-state OCP composition, but the presence of various nuclei is taken into account via an ‘impurity factor’, often taken as a free parameter adjusted on observational cooling data (see e.g. Viganò et al. 2013).

Self-consistent calculations of the impurity parameter for non-accreting unmagnetised NSs at relatively low temperatures, around the crystallisation point (k_BT ≈ 0.01 − 1 MeV, k_B being the Boltzmann constant), were recently presented by Fantina et al. (2020) and Carreau et al. (2020a). Specifically, Carreau et al. (2020a) investigated the MCP in the NS inner crust within a compressible liquid-drop (CLD) model with parameters optimised on microscopic energy-density functionals (Goriely et al. 2013). Carreau et al. (2020a) calculated the nuclear distributions using a perturbative implementation of the nuclear statistical equilibrium, as first proposed in Grams et al. (2018) for applications to supernova matter. In this perturbative treatment, for a given baryon density n_B and temperature T, the beta-equilibrium composition is firstly found in the OCP approximation. This procedure yields the (OCP) neutron and proton chemical potentials, which are subsequently used in the computation of the ion abundances. Although not fully self-consistent, this approach has the advantage of leading to a very fast convergence, with reduced computational cost with respect to a full nuclear statistical equilibrium treatment (see e.g. Botvina & Mishustin 2010; Gulminelli & Raduta 2015; Raduta & Gulminelli 2019). However, the applicability of such a perturbative approach requires a more thorough assessment at high densities, that is, in the deeper regions of the (P)NS crust, and at relatively high temperatures above the crystallisation temperature.

In the present work, we perform a study of the PNS crust in a self-consistent multi-component approach. To this aim, we extended and improved the previous work of Carreau et al. (2020a) in several aspects: (i) we extended the study of the (P)NS inner crust to higher temperatures above the crystallisation temperature, and to higher densities, until the crust–core transition; (ii) we improved the description of the ion centre-of-mass motion in the liquid crust by employing a more realistic prescription for the translational energy, where in-medium effects are accounted for (see the discussion in Dinh Thi et al. 2023 and references therein); and (iii) we went beyond the perturbative approach and performed a self-consistent calculation of the MCP in the whole PNS inner crust, which allowed us to calculate the impurity parameter in a self-consistent way.

The paper is organised as follows: The formalism is discussed in Sect. 2: the MCP approach is introduced in Sect. 2.1, the improved description of the translational free-energy term is presented in Sect. 2.2, while the ion free energy in the OCP approximation is depicted in Sect. 2.3, and the comparison between the MCP and OCP approaches is discussed in Sect. 2.4. Numerical results are presented in Sect. 3, where we describe both the implementation of the perturbative MCP approach in Sect. 3.1 and the self-consistent MCP calculations in Sect. 3.2; the outcomes regarding the impurity parameter as well as a fitting formula for its practical implementation in numerical simulations are presented in Sect. 3.3. Finally, we draw conclusions in Sect. 4.

2. Model of the PNS inner crust

2.1. MCP in nuclear statistical equilibrium

To model the full statistical distribution of ions in the PNS liquid inner crust, we extended the formalism of Carreau et al. (2020a) for its applicability at higher densities and temperatures. As in Carreau et al. (2020a), we assume that beta-equilibrium holds, a condition that is met in late PNS cooling evolution stages (Prakash et al. 1997; Yakovlev & Pethick 2004; Page et al. 2006). For a given thermodynamic condition defined by the total baryonic number density n_B and the temperature T, the PNS inner crust is supposed to be composed of different spherical Wigner-Seitz (WS) cells of volume $V_{\mathrm{WS}}^{(j)}$, each containing an ion with mass and proton number (A^(j), Z^(j)) and internal density $n_{\mathrm{i}}^{(j)}=A^{(j)}/V_{\mathrm{N}}^{(j)}$ ($V_{\mathrm{N}}^{(j)}=4\pi r_{\mathrm{N}}^{(j)3}/3$ being the volume occupied by the cluster, with $r_{\mathrm{N}}^{(j)}$ the cluster radius), such that $u^{(j)}=V_{\mathrm{N}}^{(j)}/V_{\mathrm{WS}}^{(j)}$ is the volume fraction occupied by the cluster in the cell j. Each ion is surrounded by a uniform gas of electrons and neutrons of density n_e and n_gn, respectively. In addition, we suppose that charge neutrality is realised in every cell, meaning that the proton density is the same in each cell and equal to the electron density, that is, $n_{\mathrm{e}}=n_{\mathrm{p}}=Z^{(j)}/V_{\mathrm{WS}}^{(j)}$.

At finite temperature, a free proton gas could also be present and was included in Dinh Thi et al. (2023), where the PNS liquid crust was studied in the OCP approximation. However, we find that for the considered temperature regimes, k_BT ≲ 2.0 MeV, the proton-gas density remains very small, namely about a few times 10⁻³ fm⁻³ at most at the bottom of the crust, and its effects on the equation of state and on the crust composition are negligible. For these reasons, we ignore the presence of the free protons in the present work, as was done in Carreau et al. (2020a). Moreover, in this first step to improve our previous treatment (Carreau et al. 2020a), we neglect the possible presence of non-spherical configurations, which could exist in a narrow region at the bottom of the inner crust.

Denoting $n_{\mathrm{N}}^{(j)}$, the ion density of a species (A^(j), Z^(j)), the frequency of occurrence or probability of the component (j) is given by $n_{\mathrm{N}}^{(j)}=p_j/⟨V_{\mathrm{WS}}⟩$, with ∑_jp_j = 1 and the bracket notation ⟨⟩ indicates ensemble averages. The total free-energy density¹ of the system can therefore be written as

$$\begin{array}{c}\hfill {\mathcal{F}}^{\mathrm{MCP}}={\displaystyle \sum _j{n}_{\mathrm{N}}^{(j)}(F_{\mathrm{i}}^{(j)}-V_{\mathrm{N}}^{(j)}{\mathcal{F}}_{\mathrm{g}})+{\mathcal{F}}_{\mathrm{g}}+{\mathcal{F}}_{\mathrm{e}},}\end{array}$$(1)

where ℱ_g is the free-energy density of the uniform neutron gas, see Sect. 2.3 below, and ℱ_e is the free-energy density of the uniform electron gas, for which we use the same expressions as in Carreau et al. (2020a); see Eq. (2.65) in Haensel et al. (2007). The second term in Eq. (1), namely $-V_{\mathrm{N}}^{(j)}{\mathcal{F}}_{\mathrm{g}}$, accounts for the excluded volume, that is, the impenetrability of the different ions, while the cluster free energy $F_{\mathrm{i}}^{(j)}$ is given by

$$\begin{array}{c}\hfill F_{\mathrm{i}}^{(j)}=M_{\mathrm{i}}^{(j)}{c}^2+F_{\mathrm{bulk}}^{(j)}+F_{\mathrm{Coul}+\mathrm{surf}+\mathrm{curv}}^{(j)}+F_{\mathrm{trans}}^{\star ,(j),\mathrm{MCP}},\end{array}$$(2)

where $M_{\mathrm{i}}^{(j)}=(A^{(j)}-Z^{(j)})m_{\mathrm{n}}+Z^{(j)}{m}_{\mathrm{p}}$ is the bare mass of the cluster, with m_p (m_n) being the proton (neutron) mass, c the speed of light, $F_{\mathrm{bulk}}^{(j)}=\frac{A^{(j)}}{n_{\mathrm{i}}^{(j)}}{\mathcal{F}}_{\mathrm{B}}(n_{\mathrm{i}}^{(j)},1-2Z^{(j)}/A^{(j)},T)$ the cluster bulk free energy, and $F_{\mathrm{Coul}+\mathrm{surf}+\mathrm{curv}}^{(j)}=F_{\mathrm{Coul}}^{(j)}+F_{\mathrm{surf}+\mathrm{curv}}^{(j)}$ is the sum of the Coulomb, surface, and curvature energies. Finally, the last term in Eq. (2), $F_{\mathrm{trans}}^{\star ,(j),\mathrm{MCP}}$, represents the translational free energy of the cluster in the liquid phase,

$$\begin{array}{c}\hfill F_{\mathrm{trans}}^{\star ,(j),\mathrm{MCP}}=k_{\mathrm{B}}T[ln\left(\frac{n_{\mathrm{N}}^{(j)}}{{\overline{u}}_{\mathrm{f}}}\frac{{({\lambda }_{\mathrm{i}}^{\star ,(j)})}^3}{g_{\mathrm{s}}^{(j)}}\right)-1],\end{array}$$(3)

where ${\lambda }_{\mathrm{i}}^{\star ,(j)}$ is the ion effective thermal wavelength, $g_{\mathrm{s}}^{(j)}$ is the spin degeneracy, and

$$\begin{array}{c}\hfill {\overline{u}}_{\mathrm{f}}=1-{\displaystyle \sum _j{n}_{\mathrm{N}}^{(j)}{V}_{\mathrm{N}}^{(j)}}\end{array}$$(4)

is the volume fraction available for the ion motion, ${\overline{u}}_{\mathrm{f}}=⟨V_{\mathrm{f}}⟩/⟨V_{\mathrm{WS}}⟩$, and the average free volume is given by $⟨V_{\mathrm{f}}⟩=\sum _j{p}_j{V}_{\mathrm{f}}^{(j)}$, with

$$\begin{array}{c}\hfill V_{\mathrm{f}}^{(j)}=\frac{4}{3}\pi (r_{\mathrm{WS}}^{(j)3}-r_{\mathrm{N}}^{(j)3}),\end{array}$$(5)

where $r_{\mathrm{WS}}^{(j)}$ is the WS cell radius. The translational term is discussed in more detail in Sect. 2.2, while the bulk and interaction terms are addressed in Sect. 2.3.

Following Gulminelli & Raduta (2015), Grams et al. (2018), Fantina et al. (2020), and Carreau et al. (2020a), the densities $n_{\mathrm{N}}^{(j)}$ were calculated so as to minimise the thermodynamic potential in the canonical ensemble. Because of the chosen free-energy decomposition in Eq. (1), the electron and neutron gas free-energy densities, ℱ_e and ℱ_g, do not depend on $n_{\mathrm{N}}^{(j)}$. Therefore, the variation with respect to the ion densities $n_{\mathrm{N}}^{(j)}$, $\mathrm{d}{\mathcal{F}}^{\mathrm{MCP}}\{n_{\mathrm{N}}^{(j)}\}/\mathrm{d}{n}_{\mathrm{N}}^{(j)}$, can be performed on the ion part only, that is, on the sum term in Eq. (1). However, the variations $\mathrm{d}{n}_{\mathrm{N}}^{(j)}$ are not independent because of the baryon number conservation and charge neutrality,

$$\begin{array}{cc}& n_{\mathrm{B}}=n_{\mathrm{gn}}+{\displaystyle \sum _j}{n}_{\mathrm{N}}^{(j)}(A^{(j)}-n_{\mathrm{gn}}{V}_{\mathrm{N}}^{(j)}),\hfill \end{array}$$(6)

$$\begin{array}{cc}& n_{\mathrm{e}}=n_{\mathrm{p}}={\displaystyle \sum _j}{n}_{\mathrm{N}}^{(j)}{Z}^{(j)}.\hfill \end{array}$$(7)

These constraints are taken into account by introducing two Lagrange multipliers, which can be shown to be directly connected to the chemical potentials, yielding the following equation for the equilibrium densities $n_{\mathrm{N}}^{(j)}$:

$$\begin{array}{cc}& {\displaystyle \sum _j(F_{\mathrm{i}}^{(j)}-V_{\mathrm{N}}^{(j)}{\mathcal{F}}_{\mathrm{g}}+\sum _{j^{\prime }}{n}_{\mathrm{N}}^{(j^{\prime })}\frac{\partial F_{\mathrm{i}}^{(j^{\prime })}}{\partial n_{\mathrm{N}}^{(j)}})\mathrm{d}{n}_{\mathrm{N}}^{(j)}}\hfill \\ \hfill & +{\mu }_{\mathrm{e}}{\displaystyle \sum _j{Z}^{(j)}\mathrm{d}{n}_{\mathrm{N}}^{(j)}}\hfill \\ \hfill & -{\mu }_{\mathrm{n}}{\displaystyle \sum _j(A^{(j)}-n_{\mathrm{gn}}{V}_{\mathrm{N}}^{(j)})\mathrm{d}{n}_{\mathrm{N}}^{(j)}=0.}\hfill \end{array}$$(8)

In Eq. (8), μ_e = dℱ_e/dn_e is the electron chemical potential and μ_n is given by

$$\begin{array}{c}\hfill {\mu }_{\mathrm{n}}\equiv \frac{\partial {\mathcal{F}}_{\mathrm{g}}}{\partial n_{\mathrm{gn}}}+\mathrm{\Delta }{\mu }_{\mathrm{n}}=\frac{\partial {\mathcal{F}}_{\mathrm{g}}}{\partial n_{\mathrm{gn}}}+\frac{{\sum }_j{n}_{\mathrm{N}}^{(j)}(\partial F_{\mathrm{i}}^{(j)}/\partial n_{\mathrm{gn}})}{1-{\sum }_j{n}_{\mathrm{N}}^{(j)}{V}_{\mathrm{N}}^{(j)}},\end{array}$$(9)

where one can recognise the uniform-matter chemical potential of the unbound neutrons, μ_gn ≡ ∂ℱ_g/∂n_gn, while the additional term Δμ_n accounts for the in-medium modification of the ion free energy arising from the relative motion of the ions with respect to the surrounding neutron fluid (see also Dinh Thi et al. 2023). This affects the global chemical potential of the system, breaking the identity μ_n = μ_gn, which is assumed in all OCP and MCP equation-of-state models based on the cluster degrees of freedom we are aware of; see Oertel et al. (2017) for a review. The presence of this in-medium term Δμ_n also implies that the pressure equilibrium condition inside each WS cell is modified with respect to the standard results, that is, the equivalence between the cluster and gas pressure (see e.g. Lattimer et al. 1985). In order to keep the same formal expression as in the literature – see Eq. (38) below – we include the in-medium modification in the definition of the gas pressure, which we define as P_g ≡ μ_nn_gn − ℱ_g.

The last term in the first row of Eq. (8), $\sum _{j\prime }{n}_{\mathrm{N}}^{(j\prime )}\frac{\partial F_{\mathrm{i}}^{(j\prime )}}{\partial n_{\mathrm{N}}^{(j)}}$, arises from the dependence on the ion density $n_{\mathrm{N}}^{(j)}$ of both the translational and the interaction part of the ion free energy. Indeed, because charge conservation holds both globally – see Eq. (7) – and at the level of each cell, $n_{\mathrm{e}}=n_{\mathrm{p}}=n_{\mathrm{p}}^{(j)}=Z^{(j)}/V_{\mathrm{WS}}^{(j)}$, all terms depending on the cell proton density lead to a dependence on the ion abundance $n_{\mathrm{N}}^{(j)}$. Using Eq. (3) for the translational free energy and noticing that the self-consistency is induced by the Coulomb part of the ion free energy, we obtain

$$\begin{array}{c}\hfill {\displaystyle \sum _{j^{\prime }}{n}_{\mathrm{N}}^{(j^{\prime })}\frac{\partial F_{\mathrm{i}}^{(j^{\prime })}}{\partial n_{\mathrm{N}}^{(j)}}=k_{\mathrm{B}}T+\sum _{j^{\prime }}{n}_{\mathrm{N}}^{(j^{\prime })}\frac{\partial F_{\mathrm{Coul}}^{(j^{\prime })}}{\partial n_{\mathrm{N}}^{(j)}}·}\end{array}$$(10)

The last term in Eq. (10) is the so-called rearrangement term, which was shown to be important to guarantee the thermodynamic consistency of the model and to recover the ensemble equivalence between the MCP and OCP approaches (see e.g. Grams et al. 2018; Fantina et al. 2020; Carreau et al. 2020a; Barros et al. 2020). Working out the rearrangement term, we find

$$\begin{array}{cc}\hfill {\mathcal{R}}^{(j)}& =Z^{(j)}{\displaystyle \sum _{j\prime }}{n}_{\mathrm{N}}^{(j\prime )}\frac{\partial F_{\mathrm{Coul}}^{(j^{\prime })}}{\partial n_{\mathrm{p}}}\hfill \\ \hfill & =V_{\mathrm{WS}}^{(j)}{\displaystyle \sum _{j^{\prime }}}{n}_{\mathrm{N}}^{(j^{\prime })}{V}_{\mathrm{WS}}^{(j^{\prime })}{P}_{\mathrm{int}}^{(j^{\prime })}\hfill \\ \hfill & =V_{\mathrm{WS}}^{(j)}{\overline{P}}_{\mathrm{int}},\hfill \end{array}$$(11)

where $P_{\text{int}}^{(j^{\prime })}$ is the pressure due to the Coulomb interaction (see also the discussion in Appendix A),

$$\begin{array}{c}\hfill P_{\mathrm{int}}^{(j^{\prime })}=\frac{n_{\mathrm{p}}^2}{Z^{(j^{\prime })}}\frac{\partial F_{\mathrm{Coul}}^{(j^{\prime })}}{\partial n_{\mathrm{p}}}·\end{array}$$(12)

We also note that the rearrangement term, Eq. (11), is not exactly equivalent to that derived in Carreau et al. (2020a). Indeed, in the latter work, the approximation that $F_{\mathrm{i}}^{(j)}$ was only dependent on $n_{\mathrm{N}}^{(j)}$ was made, while in reality a dependence on the other cells arises because of the Coulomb interaction, and therefore $\partial F_{\mathrm{i}}^{(j\prime )}/\partial n_{\mathrm{N}}^{(j)}\ne 0$.

Considering independent variations of $n_{\mathrm{N}}^{(j)}$ and introducing the single-ion canonical potential,

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\mathrm{i}}^{(j)}=F_{\mathrm{i}}^{\star ,(j)}-V_{\mathrm{N}}^{(j)}{\mathcal{F}}_{\mathrm{g}}+{\mathcal{R}}^{(j)}-k_{\mathrm{B}}Tln{\overline{u}}_{\mathrm{f}},\end{array}$$(13)

with

$$\begin{array}{c}\hfill F_{\mathrm{i}}^{\star ,(j)}=F_{\mathrm{i}}^{(j)}-F_{\mathrm{trans}}^{\star ,(j),\mathrm{MCP}}+k_{\mathrm{B}}Tln\left(\frac{{({\lambda }_{\mathrm{i}}^{\star ,(j)})}^3}{g_{\mathrm{s}}^{(j)}}\right),\end{array}$$(14)

Eq. (8) becomes

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\mathrm{i}}^{(j)}+k_{\mathrm{B}}Tln{n}_{\mathrm{N}}^{(j)}+{\mu }_{\mathrm{e}}{Z}^{(j)}-{\mu }_{\mathrm{n}}(A^{(j)}-n_{\mathrm{gn}}{V}_{\mathrm{N}}^{(j)})=0.\end{array}$$(15)

The ion densities can be obtained as

$$\begin{array}{c}\hfill n_{\mathrm{N}}^{(j)}=exp(-\frac{{\stackrel{\sim }{\mathrm{\Omega }}}_{\mathrm{i}}^{(j)}}{k_{\mathrm{B}}T}),\end{array}$$(16)

with

$$\begin{array}{c}\hfill {\stackrel{\sim }{\mathrm{\Omega }}}_{\mathrm{i}}^{(j)}={\mathrm{\Omega }}_{\mathrm{i}}^{(j)}+{\mu }_{\mathrm{e}}{Z}^{(j)}-{\mu }_{\mathrm{n}}(A^{(j)}-n_{\mathrm{gn}}{V}_{\mathrm{N}}^{(j)}).\end{array}$$(17)

We note that, with respect to the derivation in Carreau et al. (2020a), here a self-consistent problem arises. Indeed, the terms ${\overline{u}}_{\mathrm{f}}$, Δμ_n, and ${\overline{P}}_{\mathrm{int}}$ entering ${\stackrel{\sim }{\mathrm{\Omega }}}_{\mathrm{i}}^{(j)}$ on the right-hand side of Eq. (16) via Eqs. (4), (9) and (11), respectively, themselves depend on the ion densities $n_{\mathrm{N}}^{(j)}$. Therefore, for each given baryon density and temperature, we solve the following system of equations simultaneously:

$$\begin{array}{cc}& \mathrm{\Delta }{\mu }_{\mathrm{n}}=\frac{{\sum }_j{n}_{\mathrm{N}}^{(j)}(\partial F_{\mathrm{i}}^{(j)}/\partial n_{\mathrm{gn}})}{1-{\sum }_j{n}_{\mathrm{N}}^{(j)}{V}_{\mathrm{N}}^{(j)}},\hfill \end{array}$$(18)

$$\begin{array}{cc}& {\overline{P}}_{\mathrm{int}}={\displaystyle \sum _{j\prime }}{n}_{\mathrm{N}}^{(j\prime )}{V}_{\mathrm{WS}}^{(j^{\prime })}{P}_{\mathrm{int}}^{(j^{\prime })},\hfill \end{array}$$(19)

together with the definition of the free volume fraction, Eq. (4), and the baryon number and charge conservation equations, Eqs. (6) and (7), with $n_{\mathrm{N}}^{(j)}$ being given by Eq. (16).

2.2. Translational free energy of the ions in the MCP

The translational free-energy term appearing in Eq. (2), $F_{\mathrm{trans}}^{\star ,(j),\mathrm{MCP}}$, accounts for the centre-of-mass motion of each ion of the liquid MCP. Contrarily to the OCP description of the crust (see Dinh Thi et al. 2023 for a detailed discussion of the translational term and its impact on the OCP inner-crust composition and equation of state), in the MCP picture, the centre-of-mass position of each ion is not confined to the single WS cell; the cluster can explore the whole volume, albeit with an effective mass that accounts for the fact that only bound nucleons can be in relative motion with respect to the unbound neutron gas. This leads to Eq. (3) for the translational free energy, where we set the spin degeneracy $g_{\mathrm{s}}^{(j)}=1$, because the most abundant nuclei are expected to be even-even. The effective thermal wavelength of the ion in Eq. (3) is given by

$$\begin{array}{c}\hfill {\lambda }_{\mathrm{i}}^{\star ,(j)}=\sqrt{\frac{2\pi {ħ}^2}{M_{\mathrm{i}}^{\star ,(j)}{k}_{\mathrm{B}}T}},\end{array}$$(20)

with ℏ being the Planck constant, and the effective mass reads

$$\begin{array}{c}\hfill M_{\mathrm{i}}^{\star ,(j)}=M_{\mathrm{i}}^{(j)}(1-{\gamma }^{(j)}),\end{array}$$(21)

with ${\gamma }^{(j)}=n_{\mathrm{gn}}/n_{\mathrm{i}}^{(j)}$ being the ratio between the neutron-gas density and the cluster density. We note that Eq. (3) differs from the translational free energy used in Carreau et al. (2020a) in two respects: (i) here we use a different expression for the effective mass, that is, the one derived in the hydrodynamic approach considering that only neutrons in the continuum participate in the (collective) flow (see also Sect. 2 in Dinh Thi et al. 2023); and (ii) the centre-of-mass position of each ion of the MCP cannot freely explore the whole volume, which is consistent with the excluded-volume approach of Eq. (1). In accordance with the approach of Hempel & Schaffner-Bielich (2010), the latter condition gives rise to the extra ${\overline{u}}_{\mathrm{f}}$ term in Eq. (3).

Fantina et al. (2020) noted that, because of the different expressions of the translational energy in the OCP and MCP approaches, a first deviation from the linear mixing rule appears. As a consequence, the total free energy of the ions is not just the sum of the OCP ion free energies, but an extra term arises, known in the literature as the mixing entropy (see Eq. (21) in Fantina et al. 2020; see also Medin & Cumming 2010). This is still the case here, but the additional (non-ideal) ‘mixing entropy’ term now reflects the in-medium (excluded-volume) effects included in the translational energy $F_{\mathrm{trans}}^{\star ,(j),\mathrm{MCP}}$,

$$\begin{array}{c}\hfill F_{\mathrm{i}}^{(j)}=F_{\mathrm{i}}^{(j),\mathrm{OCP}}+k_{\mathrm{B}}Tln\left[p_j\frac{Z^{(j)}}{⟨Z⟩}\frac{{(1-{(u^{(j)})}^{1/3})}^3}{{\overline{u}}_{\mathrm{f}}}\right]·\end{array}$$(22)

The computation of the free energy of the ions in the MCP therefore still requires knowledge of the ion free energy in the OCP phase, which is addressed in Sect. 2.3.

2.3. Ion free energy in the OCP approximation

We define in this section the different terms of the ion free energy in the OCP approximation, $F_{\mathrm{i}}^{(j),\mathrm{OCP}}$; see Eq. (22)². A detailed discussion on these terms can also be found in Dinh Thi et al. (2023); here we present the main expressions, for completeness.

In the CLD approach, the ion free energy is decomposed in the bulk, Coulomb, finite-size, and translational contributions (see also Eq. (2)),

$$\begin{array}{c}\hfill F_{\mathrm{i}}^{(j),\mathrm{OCP}}=F_{\mathrm{i}}=M_{\mathrm{i}}{c}^2+F_{\mathrm{bulk}}+F_{\mathrm{Coul}+\mathrm{surf}+\mathrm{curv}}+F_{\mathrm{trans}}^{\star }.\end{array}$$(23)

The bulk free energy of the cluster, $F_{\mathrm{bulk}}=\frac{A}{n_{\mathrm{i}}}{\mathcal{F}}_{\mathrm{B}}(n_{\mathrm{i}},1-2Z/A,T)$, is expressed in terms of the free-energy density of homogeneous nuclear matter ℱ_B(n, δ, T) at a total baryonic density n = n_n + n_p (n_n and n_p being the neutron and proton densities, respectively), isospin asymmetry δ = (n_n − n_p)/n, and temperature T. We note that the same expression of the bulk energy is used to compute the free energy of the free neutron gas in Eq. (1), ℱ_g = ℱ_B(n = n_gn, δ_g = 1, T)+n_gnm_nc². To determine ℱ_B, we used the self-consistent mean-field thermodynamics (Lattimer et al. 1985; Ducoin et al. 2007; see Sect. 3.1 in Dinh Thi et al. 2023 for details), where the free-energy density is decomposed into a ‘potential’ and a ‘kinetic’ part,

$$\begin{array}{c}\hfill {\mathcal{F}}_{\mathrm{B}}(n,\delta ,T)={\mathcal{F}}_{\mathrm{kin}}(n,\delta ,T)+{\mathcal{V}}_{\mathrm{MM}}(n,\delta ).\end{array}$$(24)

The ‘kinetic’ term ℱ_kin in Eq. (24) is given by

$$\begin{array}{c}\hfill {\mathcal{F}}_{\mathrm{kin}}={\displaystyle \sum _{q=\mathrm{n},\mathrm{p}}}[\frac{-2k_{\mathrm{B}}T}{{\lambda }_q^{3/2}}{F}_{3/2}\left(\frac{{\stackrel{\sim }{\mu }}_q}{k_{\mathrm{B}}T}\right)+n_q{\stackrel{\sim }{\mu }}_q],\end{array}$$(25)

where q = n, p labels neutrons and protons, ${\lambda }_q={\left(\frac{2\pi {ħ}^2}{k_{\mathrm{B}}T{m}_q^{\star }}\right)}^{1/2}$ is the thermal wavelength of the nucleon with $m_q^{\star }$ being the density-dependent nucleon effective mass, F_3/2 denotes the Fermi-Dirac integral, and the effective chemical potential ${\stackrel{\sim }{\mu }}_q$ is related to the thermodynamical chemical potential μ_q through ${\stackrel{\sim }{\mu }}_q={\mu }_q-U_q$, with U_q the mean-field potential. Regarding the ‘potential’ energy, we employed the meta-model approach of Margueron et al. (2018a,b), where 𝒱_MM(n, δ) is expressed as a Taylor expansion truncated at order N (here, N = 4) around the saturation point (n = n_sat, δ = 0),

$$\begin{array}{c}\hfill {\mathcal{V}}_{\mathrm{MM}}^{N=4}(n,\delta )={\displaystyle \sum _{k=0}^4}\frac{n}{k!}(v_k^{\mathrm{is}}+v_k^{\mathrm{iv}}{\delta }^2)x^k{u}_k^{N=4}(x),\end{array}$$(26)

where x = (n − n_sat)/(3n_sat) and $u_k^N(x)=1-{(-3x)}^{N+1-k}$exp(−b(1 + 3x)), with b = 10ln(2) being a parameter governing the function behaviour at low densities. The factor $u_k^N(x)$ was introduced to ensure the convergence at the zero-density limit, while the parameters $v_k^{\mathrm{is}}$ and $v_k^{\mathrm{iv}}$ are linear combinations of the so-called nuclear matter empirical parameters {n_sat, E_sat, sym, L_sym, K_sat, sym, Q_sat, sym, Z_sat, sym} (see Margueron et al. 2018a for details). In this work, we employ the BSk24 empirical parameters (Goriely et al. 2013), as an illustrative example; indeed, this functional is in very good agreement with all current astrophysical data as well as nuclear physics experiments and theory (Dinh Thi et al. 2021).

The Coulomb and finite-size contributions, F_{Coul + surf + curv} = V_WS(ℱ_Coul + ℱ_{surf + curv}), account for the total interface energy. As we consider here only spherical clusters, the Coulomb energy density is given by

$$\begin{array}{c}\hfill {\mathcal{F}}_{\mathrm{Coul}}=\frac{2}{5}\pi {(e{n}_{\mathrm{i}}{r}_{\mathrm{N}})}^{2}u{\left(\frac{1-I}{2}\right)}^2[u+2(1-\frac{3}{2}{u}^{1/3})],\end{array}$$(27)

with e being the elementary charge and I = 1 − 2Z/A. For the surface and curvature contributions, we employed the same expression as in Lattimer & Swesty (1991), Maruyama et al. (2005), and Newton et al. (2013), that is

$$\begin{array}{c}\hfill {\mathcal{F}}_{\mathrm{surf}}+{\mathcal{F}}_{\mathrm{curv}}=\frac{3u}{r_{\mathrm{N}}}({\sigma }_{\mathrm{s}}(I,T)+\frac{2{\sigma }_{\mathrm{c}}(I,T)}{r_{\mathrm{N}}}),\end{array}$$(28)

where σ_s and σ_c are the surface and curvature tensions (Lattimer & Swesty 1991),

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{s}(\mathrm{c})}(I,T)={\sigma }_{\mathrm{s}(\mathrm{c})}(I,T=0)h(T).\end{array}$$(29)

The expressions of the surface and curvature tensions at zero temperature are taken from Ravenhall et al. (1983), who proposed a parameterization based on the Thomas-Fermi calculations at extreme isospin asymmetries,

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{s}}(I,T=0)& ={\sigma }_0\frac{2^{p+1}+b_{\mathrm{s}}}{y_{\mathrm{p}}^{-p}+b_{\mathrm{s}}+{(1-y_{\mathrm{p}})}^{-p}},\hfill \end{array}$$(30)

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{c}}(I,T=0)& =5.5{\sigma }_{\mathrm{s}}(I,T=0)\frac{{\sigma }_{0,\mathrm{c}}}{{\sigma }_0}(\beta -y_{\mathrm{p}}),\hfill \end{array}$$(31)

where y_p = (1 − I)/2 and the surface parameters (σ₀, σ_0, c, b_s, β, p) were optimised for each set of bulk parameters and effective mass to reproduce the experimental nuclear masses in the 2016 Atomic Mass Evaluation (AME) table (Wang et al. 2017). The temperature dependence of the surface terms is encoded in the function h(T) (Lattimer & Swesty 1991)

$$\begin{array}{c}\hfill h(T)=\{\begin{array}{cc}0\hfill & \mathrm{if}T>T_{\mathrm{c}}\hfill \\ {[1-{\left(\frac{T}{T_{\mathrm{c}}}\right)}^2]}^2\hfill & \mathrm{if}T\le T_{\mathrm{c}}\hfill \end{array},\end{array}$$(32)

where T_c is the critical temperature (see Eq. (2.31) in Lattimer & Swesty 1991).

Finally, the translational energy in the OCP approximation is given by (Dinh Thi et al. 2023)

$$\begin{array}{c}\hfill F_{\mathrm{trans}}^{\star }=k_{\mathrm{B}}Tln\left(\frac{1}{V_{\mathrm{f}}^{\mathrm{OCP}}}\frac{{\lambda }_{\mathrm{i}}^{\star 3}}{g_{\mathrm{s}}}\right)-k_{\mathrm{B}}T,\end{array}$$(33)

where the reduced (‘free’) WS volume reads (Lattimer et al. 1985)

$$\begin{array}{c}\hfill V_{\mathrm{f}}^{\mathrm{OCP}}=\frac{4}{3}\pi {(r_{\mathrm{WS}}-r_{\mathrm{N}})}^3.\end{array}$$(34)

We remark that the definition of the reduced volume, Eq. (34), differs from that introduced in Sect. 2.1 for the MCP translational free energy, Eq. (5). Indeed, in the MCP picture, the WS approximation is relaxed and the integral over the cluster centre of mass is extended to the whole (macroscopic) volume, except that occupied by the ensemble of the different ions, which leads to Eqs. (3) and (5).

2.4. MCP versus OCP

In order to compare the MCP and OCP approaches, we evaluated the most probable configuration in the MCP. The latter corresponds to the maximum probability, $p_j^{\text{max}}$, or equivalently the minimum grand-canonical potential, ${\stackrel{\sim }{\mathrm{\Omega }}}_{\mathrm{i}}^{(j),\min }$. The minimisation was performed with respect to the three variables associated to the cluster, that is ($r_{\mathrm{N}}^{(j)},I^{(j)},n_{\mathrm{i}}^{(j)})$, yielding the following system of equations:

$$\begin{array}{cc}& {\left(n_{\mathrm{i}}^{(j)}\right)}^2\frac{\partial (F_{\mathrm{i}}^{\star ,(j)}/A^{(j)})}{\partial n_{\mathrm{i}}^{(j)}}-\frac{r_{\mathrm{N}}^{(j)}{n}_{\mathrm{i}}^{(j)}}{3}\frac{\partial (F_{\mathrm{i}}^{(j),\star }/A^{(j)})}{\partial r_{\mathrm{N}}^{(j)}}=P_{\mathrm{g}},\hfill \end{array}$$(35)

$$\begin{array}{cc}& \frac{r_{\mathrm{N}}^{(j)}}{3A^{(j)}}\frac{\partial F_{\mathrm{i}}^{\star ,(j)}}{\partial r_{\mathrm{N}}^{(j)}}+(1-I^{(j)})\frac{\partial (F_{\mathrm{i}}^{\star ,(j)}/A^{(j)})}{\partial I^{(j)}}={\mu }_{\mathrm{n}}-\frac{P_{\mathrm{g}}}{n_{\mathrm{i}}^{(j)}},\hfill \end{array}$$(36)

$$\begin{array}{cc}& 2\frac{\partial (F_{\mathrm{i}}^{\star ,(j)}/A^{(j)})}{\partial I^{(j)}}-\frac{{\overline{P}}_{\mathrm{int}}}{n_{\mathrm{p}}}={\mu }_{\mathrm{e}}.\hfill \end{array}$$(37)

This system of variational equations can be compared to that obtained in the OCP approach of Carreau et al. (2020a) and Dinh Thi et al. (2023), without the free proton gas,

$$\begin{array}{cc}& n_{\mathrm{i}}^2\frac{\partial }{\partial n_{\mathrm{i}}}\left(\frac{F_{\mathrm{i}}}{A}\right)=P_{\mathrm{g}},\hfill \end{array}$$(38)

$$\begin{array}{cc}& \frac{F_{\mathrm{i}}}{A}+(1-I)\frac{\partial }{\partial I}\left(\frac{F_{\mathrm{i}}}{A}\right)={\mu }_{\mathrm{n}}-\frac{P_{\mathrm{g}}}{n_{\mathrm{i}}},\hfill \end{array}$$(39)

$$\begin{array}{cc}& 2[\frac{\partial }{\partial I}\left(\frac{F_{\mathrm{i}}}{A}\right)-\frac{n_{\mathrm{p}}}{1-I}\frac{\partial }{\partial n_{\mathrm{p}}}\left(\frac{F_{\mathrm{i}}}{A}\right)]={\mu }_{\mathrm{e}},\hfill \end{array}$$(40)

$$\begin{array}{cc}& \frac{\partial (F_{\mathrm{i}}/A)}{\partial r_{\mathrm{N}}}=0.\hfill \end{array}$$(41)

It is easy to show that Eqs. (35)–(37) for the most probable ion in the MCP are equivalent to the first three OCP variational equations Eqs. (38)–(40) if the following conditions are satisfied: (i) $F_{\mathrm{i}}^{\star ,(j)}$ is the same as in the OCP, that is, $F_{\mathrm{i}}^{\star ,(j)}=F_{\mathrm{i}}^{(j),\mathrm{OCP}}=F_{\mathrm{i}}$ (see Sect. 2.3; this is the case if the non-linear mixing term, which arises due to the translational motion (Gulminelli & Raduta 2015; Fantina et al. 2020), see Sect. 2.2, is negligible); and (ii) the gas densities, n_gn and n_e, as well as Δμ_n and ${\overline{P}}_{\mathrm{int}}$, are identical to the OCP values. This latter condition is expected to be satisfied if the ion distributions are strongly peaked on a unique ion species, such that the averages in Eqs. (18)–(19) can be approximated to a single term corresponding to the OCP solution. The validity of these conditions, which are necessary for the OCP approximation to give a satisfactory description of the finite-temperature configuration, are discussed in Sect. 3.1.

The additional OCP condition in Eq. (41) is identical to the well-known Baym virial theorem (Baym et al. 1971), which holds for a nucleus in vacuum. This condition arises from the fact that, for a given thermodynamic condition, there is systematically one more independent variable in the OCP than in the MCP. Indeed, the WS cell volume in the MCP is simply defined by charge conservation, while in the OCP it corresponds to an additional variable which has to be variationally determined.

In a CLD picture, for a given mass number A and atomic number Z, a nuclear species is also characterised by its radius, r_N, or equivalently its density $n_{\text{i}}=3A/(4\pi r_{\text{N}}^3)$, which can in principle fluctuate from cell to cell in the MCP equilibrium. Changing variables from (r_N, I, n_i) to (A, Z, n_i), the ion-radius distribution of a given (A, Z) nucleus is expressed as

$$\begin{array}{c}\hfill p_{\mathit{AZ}}(r_{\mathrm{N}})=\mathcal{N}{n}_{\mathrm{N}}(A^{(j)}=A,Z^{(j)}=Z,r_{\mathrm{N}}),\end{array}$$(42)

where 𝒩 is a normalization. This distribution will be peaked at a value r_N corresponding to $\partial {\stackrel{\sim }{\mathrm{\Omega }}}_{\mathrm{i}}^{(j)}/\partial n_{\mathrm{i}}{|}_{A^{(j)},Z^{(j)}}=0$, yielding the following pressure equilibrium condition:

$$\begin{array}{c}\hfill P_{\mathrm{cl}}^{\star }\equiv {\frac{n_{\mathrm{i}}^2}{A}\frac{\partial F_{\mathrm{i}}^{\star }}{\partial n_{\mathrm{i}}}|}_{A,Z}=P_{\mathrm{g}}.\end{array}$$(43)

Equation (38) is not necessarily compatible with the so-called virial condition, or radius-minimisation ansatz given by

$$\begin{array}{c}\hfill {\frac{\partial (F_{\mathrm{i}}^{\star }/A)}{\partial r_{\mathrm{N}}}|}_{A,Z}=0.\end{array}$$(44)

As a consequence, even if conditions (i) and (ii) above are met, the minimisation equation for the OCP radius, Eq. (41), is not necessarily fulfilled in the MCP.

This first important difference between the OCP and MCP approaches is presented in Fig. 1, where we plot the normalised distribution of the cluster radius for the most probable (A, Z) nucleus obtained at different thermodynamic conditions (orange solid lines). This distribution is compared with the solution obtained from the condition of pressure equilibrium in Eq. (43) (green dash-dotted lines) and that resulting from the minimisation of the free energy per nucleon with respect to r_N in Eq. (44) (red solid lines). To obtain these distributions, at each thermodynamic condition defined by (n_B, T), we considered independent variations of the three variables $(A^{(j)},Z^{(j)},r_{\mathrm{N}}^{(j)})$, and solved the system of equations given by Eqs. (4), (6), (7), (18), and (19), together with the equation for the ion density $n_{\mathrm{N}}^{(j)}$, Eq. (16) (see also Sect. 3.2). From Fig. 1, one can also observe that the solution of the pressure equilibrium equation always coincides with the peak of the MCP distribution. This is true not only for the radius of the most probable cluster, but also for the $r_{\mathrm{N}}^{(j)}$ of any (A^(j), Z^(j)) cluster. On the other hand, using Eq. (44) to determine the cluster radius fails in reproducing the most probable cluster radius obtained within the MCP approach. In fact, it only produces the correct peak for the most probable cluster if the OCP chemical potentials are used in the absence of the non-linear mixing term, as it was analytically demonstrated in Grams et al. (2018).

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

Computation of the crust observables considering the cluster distribution in the three-dimensional variable space (A, Z, r_N) is a heavy task from a numerical point of view. In the literature, most MCP codes designed to determine the composition and equation of state for core-collapse simulations or for PNS modelling either ignore the cluster density degree of freedom (see e.g. Hempel & Schaffner-Bielich 2010; Raduta & Gulminelli 2019) or impose the virial condition, Eq. (44) (see e.g. Furusawa et al. 2017; Barros et al. 2020; Pelicer et al. 2021). To limit the numerical cost, in the following sections we only consider a two-dimensional cluster distribution p_AZ, and fix for each (A, Z) the cluster radius r_N from the pressure equilibrium condition using Eq. (43). Indeed, we verified that, in all thermodynamic conditions explored in the present work, the radius distribution is always symmetric and relatively narrow.

3. Crust composition at finite temperature in the MCP

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

3.1. Perturbative MCP calculations

We begin this section by discussing the results obtained by calculating the MCP distribution in a perturbative approach, as first proposed in Grams et al. (2018). For each (n_B, T), we first calculated the OCP solution, solving the system of equilibrium equations Eqs. (38)–(41). This yields the OCP composition (i.e. A^OCP, Z^OCP, and the nuclear radius $r_{\mathrm{N}}^{\mathrm{OCP}}$), the OCP chemical potentials, and the neutron and electron densities, $n_{\mathrm{gn}}=n_{\mathrm{gn}}^{\mathrm{OCP}}$, $n_{\mathrm{e}}=n_{\mathrm{e}}^{\mathrm{OCP}}$. With the latter five quantities as input, and

$$\begin{array}{cc}& \mathrm{\Delta }{\mu }_{\mathrm{n}}\equiv \mathrm{\Delta }{\mu }_{\mathrm{n}}^{\mathrm{OCP}}={\left(\frac{\partial F_{\mathrm{i}}/\partial n_{\mathrm{gn}}}{V_{\mathrm{WS}}-V_{\mathrm{N}}}\right)}_{\mathrm{OCP}},\hfill \end{array}$$(45)

$$\begin{array}{cc}& {\overline{P}}_{\mathrm{int}}\equiv {\overline{P}}_{\mathrm{int}}^{\mathrm{OCP}}={\left(\frac{n_{\mathrm{p}}^2}{Z}\frac{\partial F_{\mathrm{Coul}}}{\partial n_{\mathrm{p}}}\right)}_{\mathrm{OCP}},\hfill \end{array}$$(46)

$$\begin{array}{cc}& {\overline{u}}_{\mathrm{f}}\equiv {\overline{u}}_{\mathrm{f}}^{\mathrm{OCP}}={\left(\frac{V_{\mathrm{f}}}{V_{\mathrm{WS}}}\right)}_{\mathrm{OCP}},\hfill \end{array}$$(47)

the MCP calculations were performed, giving the normalised probabilities p_j ≡ p_AZ.

As discussed in Sect. 2.4, the OCP and MCP results should coincide if, in addition, the non-linear mixing term coming from the translational free energy is set to zero. This is illustrated in Fig. 2, where the cluster radius r_N (blue curves), proton number Z (red curves), and mass number A (black curves) predicted by the OCP (solid lines) are plotted as a function of the input baryon density, n_B, together with the MCP average (dashed lines) and most probable quantities (diamonds). The latter values correspond to those maximising p_AZ. For all the considered densities (from the neutron-drip up to the crust–core transition density) and temperatures, we can see that the most probable A, Z, and r_N coincide with the OCP solutions, and follow the average values very closely, implying that the distributions are centred around the OCP predictions.

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

This is further shown in Fig. 3, where we display the joint distributions of the cluster proton number Z and mass number A for the same temperatures as in Fig. 2 and for three selected densities, n_B = 2 × 10⁻³ fm⁻³ (green contours), n_B = 10⁻² fm⁻³ (orange contours), and n_B = 2 × 10⁻² fm⁻³ (blue contours). The distributions are indeed Gaussian-like, peaked at the OCP solutions (black stars). We can also observe that, while the proton number Z is relatively constant, the most probable mass number A increases with density, as already noted in Carreau et al. (2020a) and Dinh Thi et al. (2023). Moreover, the width of the distributions gets broader with temperature and density, with lighter clusters being populated at higher T and n_B (see panel c), underlying the importance of considering a full nuclear ensemble instead of a single-nucleus (OCP) approach.

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

We now turn to a discussion of the influence of the non-linear mixing term, which comes from the translational motion. In Fig. 4, we show the evolution with the baryon density of the cluster variables A, Z, and r_N, for the same conditions as in Fig. 2, but with the translational free energy included in both the OCP and the MCP calculations, Eqs. (33) and (3), respectively. We can see that the different description of the centre-of-mass motion induces a discrepancy between the predictions of the two approaches. Indeed, in the OCP approximation, the cluster moves in the reduced (‘free’) volume $V_{\mathrm{f}}^{\mathrm{OCP}}$ – see Eqs. (33) and (34) – associated with the single WS cell determined from the variational procedure, while in the MCP, clusters of different species are considered to move in the same macroscopic volume, leading to an average ‘free’ volume ⟨V_f⟩; see Eqs. (3)–(5). This correlation between the different ion species breaks the linear mixing rule. At lower temperatures (see Fig. 4a), the deviation between the MCP and OCP predictions is negligible at densities below about 0.05 fm⁻³, implying that the influence of the non-linear mixing term is not very important. On the other hand, from Figs. 4b,c we can see that, with increasing temperature, the discrepancy between the two approaches starts to become significant at progressively lower densities. This can be explained by the fact that the contribution of the translational free energy becomes more important at higher T and n_B, as discussed in Dinh Thi et al. (2023). In particular, we can see that the distributions in the (perturbative) MCP are shifted to larger clusters with respect to the OCP solution, and that the crust–core transition occurs earlier. This is because Eq. (41), which results in small clusters in the OCP approximation (see Dinh Thi et al. 2023), no longer holds in the MCP, as discussed in Sect. 2.4. In addition, the contribution of the translational free-energy term in the OCP approximation is more important than that in the MCP. Indeed, in the former, the free volume term is a variable entering the minimisation procedure, whereas ⟨V_f⟩ in the MCP is the common volume for all nuclear species. As a result, it does not affect the nuclear distribution; see Eqs. (13) and (14).

Fig. 4. Same as in Fig. 2 but with the translational free energy, $F_{\mathrm{trans}}^{\star }$, included in both the OCP and perturbative MCP calculations.

The discontinuous behaviour of the most probable (A, Z) cluster in Fig. 4 is due to the fact that, as the temperature increases, the contribution of very light clusters, such as neutron-rich helium isotopes, becomes increasingly favoured. As shown in Fig. 5, the distributions become double-peaked when the temperature overcomes k_BT ≈ 1 MeV, and a density domain exists where the helium abundance overcomes that of the heavy cluster in the iron region predicted by the OCP. Still, it is also visible from Fig. 5 that the light-cluster contribution is limited to a small number of helium isotopes, while a large variety of nuclear species have a comparable probability around the (A^OCP, Z^OCP) value. Therefore, if we consider the probabilities of each element Z by integrating the distributions over the isotopic content (i.e. over A), the (perturbative) MCP predictions appear to be in closer agreement with the OCP approximation, at least at relatively low densities, as can be seen from Fig. 6.

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

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

The results presented in this section confirm that the non-linear mixing term induced by the translational energy leads to a breakdown of the ensemble equivalence between the OCP and MCP predictions.

The perturbative MCP approach employed here has the clear advantage of computing the full nuclear distributions at a reduced computational cost. However, this approach is not fully self-consistent, because the gas variables and chemical potentials resulting from the OCP solutions might not exactly satisfy the constraints of baryon number conservation and charge neutrality in the MCP, Eqs. (6) and (7), respectively. For this reason, we computed the full self-consistent MCP, which is discussed in the following section.

3.2. Self-consistent MCP calculations

In this section, we present the results obtained by performing fully self-consistent MCP calculations. At each thermodynamic condition defined by (n_B, T), we solved the system of equations given by Eqs. (4), (6), (7), (18), and (19), together with the equation for the ion density $n_{\mathrm{N}}^{(j)}$, Eq. (16). The translational free energy is always included in the computation of the ion free energy.

With Fig. 7, we begin by showing the evolution with the baryon number density of the neutron chemical potential μ_n presented in Eq. (9) (panel a), the electron chemical potential potential μ_e = dℱ_e/dn_e (panel b), the average free-volume fraction ${\overline{u}}_{\mathrm{f}}$ in Eq. (4) (panel c), and the MCP interaction pressure ${\overline{P}}_{\mathrm{int}}$ presented in Eq. (19) (panel d). All these quantities enter the computation of the ion abundances via Eq. (17). For illustrative purposes, we only plot the results at one selected temperature, k_BT = 1 MeV. The solutions from the MCP calculations are shown by solid lines, while the dotted lines correspond to the values calculated at the OCP composition, that is, ${\mu }_{\mathrm{n}}^{\mathrm{OCP}},{\mu }_{\mathrm{e}}^{\mathrm{OCP}},{\overline{u}}_{\mathrm{f}}^{\mathrm{OCP}}$, and ${\overline{P}}_{\mathrm{int}}^{\mathrm{OCP}}$. As we can see from panels a and b, the chemical potentials in the MCP are very similar to those in the OCP approximation; indeed, the two curves are almost indistinguishable. However, from the inset in panel a we can notice that the neutron chemical potential in the MCP is actually slightly lower than the OCP counterpart, suggesting that the corresponding gas density is not the same in the two approaches, and in particular that it is smaller in the MCP. The same trend is observed for the electron chemical potential (see Fig. 7b), except at densities above 0.06 fm⁻³, where μ_e calculated in the MCP is slightly above that obtained in the OCP approach, and therefore the corresponding density is higher in the MCP. The lower density of the unbound nucleons in the MCP was also observed by Burrows & Lattimer (1984) and by Gulminelli & Raduta (2015; see their Fig. 13). Although the difference between the MCP and OCP chemical potentials is numerically small, given that these latter enter in the calculation of the ion abundances through the exponential, Eq. (16), this difference could still lead to a significant deviation between the MCP and OCP results, as already noted by Gulminelli & Raduta (2015). From Figs. 7c,d, one can see that at relatively low baryon density, n_B ≲ 0.02 fm⁻³, ${\overline{u}}_{\mathrm{f}}$ and ${\overline{P}}_{\mathrm{int}}$ computed within the MCP are almost identical to those calculated in the OCP approximation³. We therefore expect that, at low density and temperature, the results shown in Fig. 4 are a good approximation. On the other hand, at higher densities, the absolute value of the interaction pressure, which directly enters into the computation of the rearrangement term, tends to zero in the MCP approach, while $|{\overline{P}}_{\mathrm{int}}^{\mathrm{OCP}}|$ increases. As the latter is an input in the perturbative MCP, we might also expect the discrepancies between the perturbative and self-consistent MCP approaches to become more pronounced at higher densities.

Fig. 7. Self-consistent MCP solution (solid lines) for the neutron chemical potential μn (panel a), electron chemical potential μe (panel b), average free volume fraction ${\overline{u}}_{\mathrm{f}}$ (panel c), and MCP interaction pressure ${\overline{P}}_{\mathrm{int}}$ (panel d) as a function of the total baryon density nB at kBT = 1 MeV. The dotted lines show the OCP predictions. The tick marks of the x and y axes in the insets of panels a and b are spaced by 2 × 10−4 fm−3 and 0.1 MeV, respectively.

Indeed, this is shown in Fig. 8, where the average (dash-dotted violet lines) and the most probable (violet squares) values of Z (top panel) and A (lower panel) in the self-consistent MCP calculations are plotted as a function of n_B at k_BT = 1 MeV, together with the OCP predictions (black solid lines) and the average (orange dashed lines) and most probable (orange dots) values obtained within the perturbative MCP procedure. Below about 0.02 fm⁻³, the results of the three different approaches almost coincide. However, as the density increases, the outcomes of the different methods start to diverge, and therefore the perturbative MCP and OCP predictions become less reliable. On the one hand, these results confirm the validity of these approximations at relatively low densities, while on the other, they highlight the importance of full MCP calculations in the deeper region of the PNS inner crust.

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

From Fig. 8, it is also interesting to see that, in the self-consistent MCP calculations, lighter nuclei dominate at high density, which is in agreement with previous works based on the nuclear statistical equilibrium (see e.g. Souza et al. 2009; Gulminelli & Raduta 2015). Conversely, in the OCP and perturbative MCP approximations, heavier clusters survive until the crust–core transition.

In order to further investigate this point, we plot in Fig. 9 the normalised distributions of the cluster proton number Z (left panel), ∑_Ap_AZ, and mass number A (right panel), ∑_Zp_AZ, at different densities in the inner crust, n_B ∈ [0.001, 0.04] fm⁻³, at k_BT = 1 MeV, for both the self-consistent MCP calculations (violet areas) and the perturbative MCP ones (orange areas). As one may expect, at lower densities, the distributions obtained in the two approaches are almost identical; they are narrow and peaked at the OCP predictions (green triangles). Therefore, employing a perturbative MCP where the gas variables are fixed from the OCP solution is a very good approximation in these thermodynamic conditions. As the density increases, at n_B ≈ 0.02 − 0.03 fm⁻³, the self-consistent MCP distributions are displaced towards lower values of A and Z, and even exhibit a double-peak structure, meaning that lighter clusters coexist with heavier ones and their contribution becomes important. In these regimes, the OCP and perturbative MCP predictions significantly overestimate the cluster mass and proton number. Eventually, at the bottom of the crust, nuclei with Z < 5 and A < 50 may even dominate, with helium being the most probable element, although the tails of the distributions extend up to Z ≃ 50 and A ≃ 400, implying that heavier nuclei may still be present near the crust–core transition. However, it is to be noted that, while the Z = 2 element appears to be the most abundant, ⁴He, which is the most bound helium isotope in the vacuum, is not the most abundant cluster (see also Fig. 8). This is because in the NS inner crust, the medium is very dense and (helium) clusters are immersed in continuum neutron states. Indeed, as already noted by Gulminelli & Raduta (2015), alpha particles are abundant for matter close to isospin symmetry, while more neutron-rich (hydrogen and helium) isotopes prevail in neutron-rich matter.

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

The findings presented above were obtained at a given temperature, namely k_BT = 1 MeV. In order to study the dependence of the results on the temperature, we ran the calculations for different temperatures up to k_BT = 2 MeV. In Fig. 10, we plot the distributions of Z (left panel) and A (right panel), defined as in Fig. 9, at n_B = 0.001 fm⁻³ in the aforementioned temperature range. As in Fig. 9, the results from both the self-consistent MCP calculations (violet areas) and the perturbative MCP ones (orange areas) are displayed and compared with the OCP solutions (green triangles). At lower temperatures, the peaks of the distribution coincide with or are very close to the OCP predictions. However, as the temperature increases, more nuclear species are populated, and therefore the distributions flatten and the OCP prediction tends to overestimate the cluster size. At higher temperature, k_BT ≥ 1.8 MeV, small clusters start to appear. At 2 MeV, the double-peak structure already observed at high density is clearly visible in the self-consistent MCP calculations, while it is absent in the perturbative MCP results. Nevertheless, the deviation between the predictions of the two approaches is not as significant as in Fig. 9, which suggests that the effect of density overcomes that of temperature as far as the cluster distributions are concerned.

Fig. 10. Same as Fig. 9, but at a fixed baryon density nB = 0.001 fm−3 and for different temperatures, kBT ∈ [1.0, 2.0] MeV.

A complementary comparison of the distributions obtained in the perturbative and self-consistent MCP is given in Fig. 11, where the normalised distributions ∑_Ap_AZ are plotted as a function of Z for k_BT = 2 MeV and for three selected densities: n_B = 0.01 fm⁻³ (panel a), n_B = 0.02 fm⁻³ (panel b), and n_B = 0.03 fm⁻³ (panel c). These conditions correspond to those for which a double-peaked distribution is observed in the perturbative MCP calculations (see Fig. 5). For comparison, the OCP results are also shown (green triangles). At all three densities considered, the distributions obtained with the self-consistent MCP are shifted towards lower Z with respect to those of the perturbative MCP calculations, as already observed in Fig. 9 at k_BT = 1 MeV. Indeed, at k_BT = 2 MeV, in the self-consistent MCP, light nuclei are always more abundant than heavier ones. Moreover, the distributions become peaked around Z ≈ 2 (see panels b and c), while they remain relatively broad and peaked towards Z ≈ 30 in the perturbative MCP.

Fig. 11. Normalised distributions P(Z) = ∑ApAZ as a function of Z in a fully self-consistent (‘SC’) and a perturbative (‘P’) MCP calculation. Calculations are performed for kBT = 2 MeV and nB = 0.01 fm−3 (panel a), nB = 0.02 fm−3 (panel b), and nB = 0.03 fm−3 (panel c). Green triangles indicate the OCP solutions.

The results presented in Figs. 8–11 show that even for temperatures as low as 1 MeV (and more so at higher temperatures), the OCP and the perturbative MCP approximations are no longer reliable in the deepest region of the crust. Therefore, for studies requiring accurate knowledge of the crust composition, a full MCP calculation is needed. However, the effect is less important as far as more global properties such as the equation of state are concerned. This is illustrated in Fig. 12, where we display the neutron and electron gas density over the baryon number density, n_gn/n_B and n_e/n_B, respectively, and in Fig. 13, where we plot the total pressure versus the mass-energy density in the PNS inner crust. Results are shown for three different temperatures, k_BT = 1 MeV (red lines), 1.5 MeV (blue lines), and 2 MeV (black lines), for both the self-consistent MCP (solid lines) and the OCP approximation (dotted lines). From Fig. 12, we can see that the MCP neutron gas densities are slightly smaller than the OCP ones in all the considered cases, although differences of up to about 20% are observed for k_BT = 2 MeV at the highest densities in the crust. As for the electron density, that calculated within the MCP is generally slightly higher than that calculated in the OCP approximation, except possibly at k_BT = 1 MeV around a few 10⁻² fm⁻³, in accordance with Fig. 7. From Fig. 13, we note that, concerning the equation of state, at k_BT = 1 MeV, the pressure-density relations provided within the self-consistent MCP approach and the OCP approximation are very similar at all densities, except very close to the crust–core transition. At higher temperatures, deviations between the two calculations start to emerge, particularly at high density. However, these discrepancies amount to approximately 10% at most in the vicinity of the crust–core transition.

Fig. 12. Ratio of the neutron gas density (upper panel) and of the electron gas density (lower panel) over the baryon density, ngn/nB and ne/nB, respectively, as a function of the total baryon density nB resulting from the self-consistent MCP calculations (solid lines) at kBT = 1 MeV (red lines), 1.5 MeV (blue lines), and 2 MeV (black lines). For comparison, the OCP results are plotted with dotted lines.

Fig. 13. Total pressure P vs. mass-energy density ρB resulting from the self-consistent MCP calculations (solid lines) at kBT = 1 MeV (red lines), 1.5 MeV (blue lines), and 2 MeV (black lines). For comparison, the OCP results are plotted with dotted lines. In the inset, the tick marks on the x and y axis are spaced by 1012 g cm−3 and 2 × 1030 dyn cm−2, respectively.

3.3. Impurity parameter

Within our MCP approach, it is also possible to self-consistently calculate the impurity parameter, which is defined as the variance of the Z distribution,

$$\begin{array}{c}\hfill Q_{\mathrm{imp}}=⟨Z^2⟩-{⟨Z⟩}^2.\end{array}$$(48)

The latter is plotted in Fig. 14 for three temperatures, k_BT = 1.0 MeV (red lines), k_BT = 1.5 MeV (blue lines), and k_BT = 2.0 MeV (black lines). Results obtained with the self-consistent (perturbative) MCP calculations as a function of the baryon density n_B in the crust are shown as solid (dash-dotted) lines. At the lower temperature of k_BT = 1.0 MeV (red lines), the predictions from the two treatments coincide at lower densities, until n_B ≈ 0.01 fm⁻³. With increasing density, the impurity parameter is, at first, larger in the self-consistent MCP approach with respect to the perturbative one. However, at variance with the perturbative MCP predictions, Q_imp does not increase monotonically in the self-consistent calculations, but peaks around n_B ≈ 0.025 fm⁻³, reaching Q_imp ≈ 100 for the considered BSk24 model, and subsequently decreases at higher densities. A similar behaviour is observed for higher temperatures, although the discrepancy between the two treatments and the peak in the impurity parameter appear at smaller densities. This can be understood from the Z distributions shown in the left panel of Figs. 9 and 10 and in Fig. 11. Indeed, while the presence of the second peak at small Z for moderate densities increases the variance of Z, the transition to light nuclei at high density overall decreases the value of the impurity parameter. These findings also show that the impurity parameter calculated within the perturbative MCP approach is underestimated at low densities, and is severely overestimated at high densities.

Fig. 14. Impurity parameter Qimp as a function of the baryon density nB in the inner crust for kBT = 1.0 MeV (red lines), kBT = 1.5 MeV (blue lines), and kBT = 2.0 MeV (black lines) in the full MCP calculation (solid lines). For comparison, results obtained with the perturbative MCP (dash-dotted lines) are also shown.

In the cooling process of a NS, it is reasonable to assume that the crust composition is frozen after the solidification of the crust. Neutron absorption or β-decays might still occur at lower temperatures, but pycnonuclear reactions that involve overcoming a Coulomb barrier might be considerably inhibited even above crystallisation. A realistic estimate of the temperature at which the ion distribution is frozen and the impurity parameter is settled would require a comparison between the cooling time and the different reaction rates, which is beyond the scope of the present study. As a reasonable estimate, the impurity parameter of the solid crust can be computed from that obtained at the crystallisation temperature of a pure Coulomb plasma (Haensel et al. 2007):

$$\begin{array}{c}\hfill T_{\mathrm{m}}\approx \frac{Z^2{e}^2}{k_{\mathrm{B}}{r}_{\mathrm{WS}}{\mathrm{\Gamma }}_{\mathrm{m}}},\end{array}$$(49)

where Z and r_WS correspond to the ground-state composition at zero temperature, and Γ_m ≈ 175 is the Coulomb coupling parameter at the melting point. It is interesting to observe that this simple expression was shown in Carreau et al. (2020b; see their Fig. 4) to give a good order-of-magnitude estimate of the temperature where the OCP free-energy densities in the liquid equate to those of the solid phase.

Values of k_BT_m from Eq. (49) are displayed in Fig. 15a: k_BT_m increases monotonically from ∼0.3 MeV at n_B = 10⁻³ fm⁻³ to ∼1.0 MeV at n_B = 0.07 fm⁻³ (see also Fig. 4 in Carreau et al. 2020b). The corresponding evolution of the impurity parameter in the fully self-consistent MCP calculations at the crystallisation point as a function of the baryon density in the inner crust is shown with green points in Fig. 15b. We can see that a high impurity parameter 10 ≲ Q_imp ≲ 100 should be expected in the whole inner crust, which could have important consequences for the magneto-thermal evolution of X-ray pulsars (see Pons et al. 2013; Newton 2013; Horowitz et al. 2015).

Fig. 15. Impurity parameter at different temperatures. Panel a: crystallisation temperature kBTm calculated using Eq. (49) as a function of the baryon density nB. Panel b: impurity parameter Qimp in the fully self-consistent MCP calculation as a function of nB evaluated at the crystallisation temperature kBTm (green points) and at four fixed temperatures: kBT = 0.5 MeV (magenta points), kBT = 1.0 MeV (red points), kBT = 1.5 MeV (blue points), and kBT = 2.0 MeV (black points). Solid lines are obtained from the fit given by Eq. (50). Panel c: absolute error in the fitting formula (|ΔQimp| = |calculated Qimp – fit|).

For comparison, we also plot Q_imp for four different selected temperatures: k_BT = 0.5 MeV (magenta points), k_BT = 1.0 MeV (red points), k_BT = 1.5 MeV (blue points), and k_BT = 2.0 MeV (black points). We can see that the general behaviour of Q_imp is similar for all temperatures; that is, all curves show a peak and a subsequent drop. Indeed, even at a relatively low temperature of k_BT = 0.5 MeV, Q_imp is still reduced at high densities, n_B ≈ 0.06 fm⁻³, indicating that the charge distribution is dominated by small Z. Moreover, as expected, a higher Q_imp is obtained if the composition is frozen at higher temperature, but the sharp drop occurs at lower density because of the appearance of light nuclei (see also Figs. 9–11). However, our results should be considered with care at high density, because non-spherical pasta structures could appear, and are not considered in the present paper. Such structures are also expected to be associated with high impurity parameter values Q_imp ∼ 30; see e.g. Schneider et al. (2016) and Caplan et al. (2021). Also, the CLD description can be questioned for light clusters and a more microscopic treatment of the in-medium modifications might be needed; see Pais et al. (2018) and Röpke (2020).

The results for the impurity parameters shown in Fig. 15 are given in Table 1 and are available in tabular format at the CDS. For practical applications to numerical simulations, we also provide a fitting formula:

$$\begin{array}{c}\hfill Q_{\mathrm{imp}}=exp{\displaystyle (a_0+\sum _{k=0}^4{x}_{\mathrm{B}}^k\sum _{j=1}^3{a}_{\mathit{kj}}{x}_{\mathrm{T}}^j)+\mathrm{\Delta }Q,}\end{array}$$(50)

where x_B ≡ n_B/fm⁻³, x_T ≡ k_BT/MeV, a₀ = 3.328, and a_kj are given in Table 2, and the function ΔQ reads

$$\begin{array}{c}\hfill \mathrm{\Delta }Q=\frac{f_0{x}_{\mathrm{T}}{(x_{\mathrm{BT}}(n_{\mathrm{B}},T))}^{f_1}}{[{(x_{\mathrm{BT}}(n_{\mathrm{B}},T)/f_2+f_3)}^4+f_4]f(n_{\mathrm{B}},T)},\end{array}$$(51)

with f₀ = 65.55, f₁ = 9.046 × 10⁻², f₂ = −6.324, f₃ = −1.087 × 10⁻², and f₄ = 4.802. The variable x_BT in Eq. (51) depends on both n_B and T and is defined as

$$\begin{array}{c}\hfill x_{\mathrm{BT}}(n_{\mathrm{B}},T)\equiv {(5+{log}_{10}{x}_{\mathrm{B}})}^2{x}_{\mathrm{T}},\end{array}$$(52)

while f(n_B, T) is given by

$$\begin{array}{c}\hfill f(n_{\mathrm{B}},T)=1+x_{\mathrm{B}}{x}_{\mathrm{T}}^{3/2}exp(\frac{200}{3}{x}_{\mathrm{B}}-\frac{1}{4})·\end{array}$$(53)

The impurity parameter calculated at different temperatures with the fitting expression Eq. (50) is shown by solid lines in Fig. 15b, while the absolute error with respect to the computed values of Q_imp in the MCP approach is displayed in panel c. We can see that the error remains relatively small, |ΔQ_imp|≈5 at most, in the whole inner-crust density range. Indeed, even if the fit tends to break down at higher densities and higher temperatures (see black and blue lines in Fig. 15b for n_B ≳ 0.05 fm⁻³), the error is very small because of the low values of Q_imp ≈ 0.1 − 0.2.

4. Conclusions

In this work, we studied the properties of the beta-equilibrated inner crust of a PNS in the liquid phase within a multi-component approach. To this aim, the formalism of Fantina et al. (2020) and Carreau et al. (2020a) was extended for its applicability over a wider range of densities and temperatures in the inner crust, from the neutron drip up to the crust–core transition. The cluster energetics was described within a CLD model approach, in which the surface parameters were optimised to reproduce the experimental masses from the AME2016 (Wang et al. 2017), and the nuclear matter properties were calculated within finite-temperature mean-field thermodynamics. We performed the calculations employing both a perturbative MCP treatment, where the chemical potentials and gas densities were taken from the OCP calculations, and a full self-consistent MCP approach.

We find that, if non-linear mixing terms arising from the centre-of-mass motion of the ion are neglected, the most probable clusters in the MCP (perturbative) approach coincide with the OCP predictions. However, as already discussed in Dinh Thi et al. (2023), this translational free-energy term should be taken into account in the calculations of the (liquid) inner-crust composition, leading to a breakdown in the ensemble equivalence. Our outcomes also show that the OCP and the perturbative MCP treatments are a good approximation at relatively low densities and temperatures, especially as far as global properties such as the equation of state are concerned. However, we believe that a full self-consistent MCP is needed for reliable predictions of the PNS crust composition, particularly in the deeper region of the crust.

Moreover, our findings reveal that, with increasing density and temperature, the abundance of light nuclei becomes important, and eventually dominates the whole distribution. This result has an important effect on the impurity parameter, and therefore potentially on the NS cooling. For applications to numerical simulations, we also provide both values of the impurity parameter calculated within our fully self-consistent MCP approach in tabular form and a fitting formula able to reproduce the numerical results with good accuracy.

In our study, we only considered spherical clusters. However, non-spherical structures, so-called pasta phases, may appear at the bottom of the inner crust (see Pelicer et al. 2021 for recent calculations of pasta phases in a MCP approach at finite temperature within a relativistic mean-field functional). We defer an investigation of the co-existence of different (non-spherical) structures within our approach to future studies.

Acknowledgments

This work has been partially supported by the IN2P3 Master Project NewMAC, the ANR project ‘Gravitational waves from hot neutron stars and properties of ultra-dense matter’ (GW-HNS, ANR-22-CE31-0001-01), and the CNRS International Research Project (IRP) ‘Origine des éléments lourds dans l’univers: Astres Compacts et Nucléosynthèse (ACNu)’.

References

Appendix A: Pressure in the MCP

The pressure in each WS cell can be written as

$$\begin{array}{c}\hfill P^{(j)}=P_i^{(j)}+P_{\mathrm{gn}}+P_e,\end{array}$$(A.1)

where the gas and electron pressures are the same in all cells since both gases are uniformly distributed over the whole volume, and $P_i^{(j)}$ is the ion pressure,

$$\begin{array}{c}\hfill P_i^{(j)}=P_{\mathrm{trans}}^{\star ,(j)}+P_{\mathrm{int}}^{(j)}.\end{array}$$(A.2)

The first term in Eq. (A.2), $P_{\mathrm{trans}}^{\star ,(j)}$, accounts for the contribution from the translational centre-of-ion motion:

$$\begin{array}{c}\hfill P_{\mathrm{trans}}^{\star ,(j)}=P_{\mathrm{trans}}^{\star }=\frac{k_{\mathrm{B}}T}{{\overline{u}}_{\mathrm{f}}⟨V_{\mathrm{WS}}⟩},\end{array}$$(A.3)

while the second term, $P_{\mathrm{int}}^{(j)}$, is the interaction (Coulomb) term of the pressure as calculated in the (pure phase) OCP approximation (see also Eq. (12)).

To derive the total pressure in the MCP, we used the definition of the pressure in the canonical ensemble,

$$\begin{array}{c}\hfill P^{\mathrm{MCP}}=-{\frac{\partial F}{\partial V}|}_{p_j,T}=P_{\mathrm{gn}}+P_e+P_i,\end{array}$$(A.4)

where P_gn = μ_gnn_gn − ℱ_g (P_e = μ_en_e − ℱ_e) is the neutron gas⁴ (electron) pressure, P_i is the ion pressure, p_j is the set of probabilities of the different species, with p_j being the probability of each component j. However, as already pointed out by Fantina et al. (2020), the total pressure in MCP is not equal to the sum of the pressure of the OCP phases, that is, P^MCP ≠ ∑_jp_jP^(j), but it is rather given by

$$\begin{array}{c}\hfill P^{\mathrm{MCP}}=P_{\mathrm{gn}}+P_e+P_{\mathrm{trans}}^{\star }+\frac{n_p^2}{⟨Z⟩}{\displaystyle \sum _j}{p}_j\frac{\partial F_{\mathrm{Coul}}^{(j)}}{\partial n_p}.\end{array}$$(A.5)

All Tables

All Figures

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

In the text

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

In the text

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

In the text

	Fig. 4. Same as in Fig. 2 but with the translational free energy, $F_{\mathrm{trans}}^{\star }$, included in both the OCP and perturbative MCP calculations.
In the text

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

In the text

	Fig. 6. Normalised distributions P(Z) = ∑ApAZ, where pAZ are taken from Fig. 5, as a function of Z in a perturbative MCP calculation. Calculations are performed for kBT = 2 MeV and nB = 0.01 fm−3 (green dash-dotted line), nB = 0.02 fm−3 (red solid line), and nB = 0.03 fm−3 (black dashed line). Arrows indicate the OCP solutions.
In the text

Fig. 7. Self-consistent MCP solution (solid lines) for the neutron chemical potential μn (panel a), electron chemical potential μe (panel b), average free volume fraction ${\overline{u}}_{\mathrm{f}}$ (panel c), and MCP interaction pressure ${\overline{P}}_{\mathrm{int}}$ (panel d) as a function of the total baryon density nB at kBT = 1 MeV. The dotted lines show the OCP predictions. The tick marks of the x and y axes in the insets of panels a and b are spaced by 2 × 10−4 fm−3 and 0.1 MeV, respectively.

In the text

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

In the text

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

In the text

	Fig. 10. Same as Fig. 9, but at a fixed baryon density nB = 0.001 fm−3 and for different temperatures, kBT ∈ [1.0, 2.0] MeV.
In the text

	Fig. 11. Normalised distributions P(Z) = ∑ApAZ as a function of Z in a fully self-consistent (‘SC’) and a perturbative (‘P’) MCP calculation. Calculations are performed for kBT = 2 MeV and nB = 0.01 fm−3 (panel a), nB = 0.02 fm−3 (panel b), and nB = 0.03 fm−3 (panel c). Green triangles indicate the OCP solutions.
In the text

Fig. 12. Ratio of the neutron gas density (upper panel) and of the electron gas density (lower panel) over the baryon density, ngn/nB and ne/nB, respectively, as a function of the total baryon density nB resulting from the self-consistent MCP calculations (solid lines) at kBT = 1 MeV (red lines), 1.5 MeV (blue lines), and 2 MeV (black lines). For comparison, the OCP results are plotted with dotted lines.

In the text

	Fig. 13. Total pressure P vs. mass-energy density ρB resulting from the self-consistent MCP calculations (solid lines) at kBT = 1 MeV (red lines), 1.5 MeV (blue lines), and 2 MeV (black lines). For comparison, the OCP results are plotted with dotted lines. In the inset, the tick marks on the x and y axis are spaced by 1012 g cm−3 and 2 × 1030 dyn cm−2, respectively.
In the text

	Fig. 14. Impurity parameter Qimp as a function of the baryon density nB in the inner crust for kBT = 1.0 MeV (red lines), kBT = 1.5 MeV (blue lines), and kBT = 2.0 MeV (black lines) in the full MCP calculation (solid lines). For comparison, results obtained with the perturbative MCP (dash-dotted lines) are also shown.
In the text

Fig. 15. Impurity parameter at different temperatures. Panel a: crystallisation temperature kBTm calculated using Eq. (49) as a function of the baryon density nB. Panel b: impurity parameter Qimp in the fully self-consistent MCP calculation as a function of nB evaluated at the crystallisation temperature kBTm (green points) and at four fixed temperatures: kBT = 0.5 MeV (magenta points), kBT = 1.0 MeV (red points), kBT = 1.5 MeV (blue points), and kBT = 2.0 MeV (black points). Solid lines are obtained from the fit given by Eq. (50). Panel c: absolute error in the fitting formula (|ΔQimp| = |calculated Qimp – fit|).

In the text