Thermodynamically consistent accreted crust of neutron stars:
The role of proton shell effects

In Ref. [25] we started the derivation of the equations for the nHD crust from the two-parametric mEOS $\epsilon(n_{b},n_{N})$ and used the CLD model for numerical illustration. The CLD-based mEOS can be represented in this form by imposing beta-equilibrium, mechanical equilibrium, and chemical equilibrium conditions inside a spherical cell, which contains one nucleus (see the Supplementary Material of that work). Here, we take a step back and start from a three-parameter mEOS, $\epsilon(n_{b},n_{N},Z)$ (see Section III.1 for numerical implementation). This allows us to consider not only beta-equilibrated matter, but also a EOS for which $Z$ is constant in the inner crust. The latter model, referred to as “$Z$-fixed” EOS in what follows, is used as a simplified, yet adequate model mimicking the strong proton shell closure of nuclei with $Z=20$.

Let us introduce the chemical potentials, $\mu_{N}=\partial\epsilon/\partial n_{N}$ and $\mu_{b}=\partial\epsilon/\partial n_{b}$. The first one, $\mu_{N}$, describes the energy change resulting from the addition of a nucleus at fixed $n_{b}$ and $Z$ (alternatively, one can consider it as the creation of a nucleus from nucleons already available in the matter), while $\mu_{n}$ corresponds to the energy change due to an additional baryon at fixed $n_{N}$ and $Z$. Since the proton number density, $Zn_{N}$, remains also unchanged, the added baryon is a neutron and $\mu_{b}$ can be identified with the neutron chemical potential, $\mu_{n}$: $\mu_{b}=\mu_{n}$. Below, we will write $\mu_{n}$ instead of $\mu_{b}$ in all formulas.

In this section, we consider the two cases: (i) $Z$ is conserved during compression ($Z$-fixed EOS), and (ii) beta-equilibrated matter, for which (see Section III.1)

In both cases, the pressure can be written as

Here, the second equality can be derived straightforwardly by taking partial derivative at a fixed nucleon and baryon number, and using the definitions of the chemical potentials $\mu_{n}$ and $\mu_{N}$; see also Section II.1.

According to one of the Tolman-Oppenheimer-Volkoff equations [32]

Here the prime denotes the derivative with respect to the radial coordinate $r$. Combined with the nHD condition (1) and the Gibbs-Duhem relation in the form $dP=n_{b}d\mu_{n}+n_{N}d\mu_{N}$, which holds true for both the $Z$-fixed EOS and beta-equilibrated EOS, we arrive at the condition $\mu_{N}^{\infty}=\mathrm{constant}$. Thus,

Here $C$ is a constant that depends on the pressure $P_{\rm oi}$ at the outer-inner crust boundary. As we demonstrated in Ref. [25] (see also Section V), FAC EOS corresponds to a certain value of $C$. However, to determine this value, it is instructive to consider the whole nHD EOS family, i.e., a family of EOSs that are allowed by the nHD condition (1). The nHD EOS family can parametrized by the pressure $P_{\rm oi}$ (or, equally, by $C$).

The catalyzed crust corresponds to the global minimum of the energy density $\epsilon(n_{b},n_{N},Z)$ at fixed $n_{b}$, which is given by the beta-equilibrium condition (2) and the condition $\mu_{N}=\partial\epsilon/\partial n_{N}=0$. As pointed out in Ref. [25], it is a member of the beta-equilibrium nHD EOS family, corresponding to $C=0$.

For a realistic modeling of the accreted crust, one should utilize the mEOS, which takes into account nuclear shell effects. Obviously, the corresponding mEOS is more complicated than its CLD-based smooth analogue. In particular, the nuclear charge number $Z$ becomes discrete (integer), and the energy density dependence on $Z$ becomes rather complicated (see, e.g., Refs. [33, 31] and Section III.2). As a result, the beta-equilibrium condition (2) cannot be applied and should generally be replaced with some other requirement. Moreover, Eq. (5) can be violated if $Z$ varies in the inner crust due to nuclear reactions. To address these difficulties, we should modify our approach as discussed in Ref. [34].

Namely, we should replace Eq. (5) with its more general counterpart which, similarly to Eq. (5) follows from the Tolman-Oppenheimer-Volkoff equation (4) and nHD condition (1):

where, as we already indicated in the introduction, $m_{n}c^{2}$ is the neutron chemical potential at the top of the inner crust (located at $P=P_{\rm oi}$). This equation allows us to determine the chemical potential $\mu_{n}$ in the layer with the pressure $P$, if the EOS [e.g., the function $\epsilon(P)$] is known in the range from $P_{\rm oi}$ to $P$. In accordance with Refs. [35, 36, 37], our approach to constructing the inner crust model involves a step-by-step progression into deeper layers of the crust. At each layer, we determine the charge number $Z$ by minimizing the corresponding thermodynamic potential by means of permissible nuclear reactions (see Section II.3 for more details). Previously (e.g., in Refs. [35, 36, 37]), the redistribution of unbound neutrons was disregarded. Consequently, the reactions were assumed to occur at a constant pressure, and the potential to be minimized was associated with the Gibbs free energy (e.g., [35, 36]). However, as stipulated by Eq. (6), in the nHD crust, not only is the pressure $P$ fixed in a given layer, but also the chemical potential $\mu_{n}$. As shown in [34], the appropriate thermodynamic potential that should be minimized at fixed $P$ and $\mu_{n}$ is then

Here $V$ is the volume attached to nuclei, and $N_{b}=Vn_{b}$ is the total number of baryons in the volume.

If (pycnonuclear) fusion and fission reactions are not allowed (either too slow or energetically forbidden), the number of nuclei is conserved. In this case, to determine $Z$ one can minimize the potential $\Psi$ per one nucleus (note that it coincides with $\mu_{N}$ in a special case considered in this work when nuclei of only one species are present at any given pressure in the crust):

This minimization effectively replaces the beta-equilibrium condition. For a smooth CLD model with continuous $Z$, it reduces to the equation (2); see Section III.2 for details.

When the composition of a given layer is determined, we can apply Eq. (6) to consider the underlying layer. By repeating this procedure, we develop a step-by-step algorithm, which is formulated in the next subsection and applied in this work to construct the nHD crust.

Similar to the case of a smooth mEOS, we begin by constructing the nHD EOS family, which is parametrized by the pressure $P_{\mathrm{oi}}$ at the outer-inner crust interface. In subsequent sections, we discuss how to constrain the range of realistic $P_{\mathrm{oi}}$ and present numerical results for the CLD+sh model of Ref. [31].

To build the nHD EOS family, we apply an algorithm based on Eqs. (6) and (7). We stress that this algorithm is quite general and, in particular, applicable in the situation when both shell effects and odd-even staggering of nuclear energies are allowed for. Since the thermodynamic potential $\Psi$ should be minimized at fixed $P$ and $\mu_{n}$ we assume the mEOS to be parameterized by the pressure $P$, neutron chemical potential $\mu_{n}$, and nuclear charge number $Z$. In other words, the microphysical model should allow one to calculate the energy density $\epsilon$ and potential $\psi$ for a given $P$, $\mu_{n}$, and $Z$ (see Section III.2 for the realization of this form of mEOS based on the CLD+sh model, which is used in this work as a numerical example).

The algorithm contains the following stages (stages 2 and 3 are repeated at each step):

1. 1.

   Specify the initial conditions at the top of the inner crust.

   The EOS of the fully accreted outer crust is well studied, especially for a one-component ash discussed here (see, e.g., Refs. [35, 38]). This allows us to determine the charge number at the bottom of the outer crust, located at $P=P_{\mathrm{oi}}$. Combining with the neutron chemical potential at the top of the inner crust, which is equal to $\mu_{n,{\rm oi}}=m_{n}c^{2}$ by definition, we obtain the initial conditions (the $j=0$ step) for the construction of the nHD inner crust. Before turning to the next stage, we also analyze the reaction pathways at the oi interface according to stage 3 to determine the composition at the top of the inner crust.

2. 2.

   Advancing to $j$-th layer.

   Starting this step, we assume that the equation of state in the previous, $(j-1)$-th layer, has been determined and the respective pressure is $P_{j-1}$, the neutron chemical potential is $\mu_{n,j-1}$, and the charge number is $Z_{j-1}$. We increase the pressure by a small amount $\Delta P$, such that $P_{j}=P_{j-1}+\Delta P$. To guarantee the nHD equilibrium, we also increase the neutron chemical potential according to the formula (6):

   $\displaystyle\mu_{n,j}=\mu_{n,j-1}\,\cdot\exp\left[\frac{\Delta P}{\epsilon(P_{j-1})+P_{j-1}}\right]$

   (9)

   The composition will be adjusted to the equilibrium one at the next stage, and here we simply assign $Z_{j}=Z_{j-1}$.

3. 3.

   Choosing optimal $Z_{j}$

   As inputs for this stage, we have the pressure $P_{j}$ and neutron chemical potential $\mu_{n,j}$ at the layer $j$, but the nuclear charge $Z_{j}$ may differ from the optimal value. It can either be equal to $Z_{j-1}$ or correspond to the value obtained in the preceding iteration of this stage. In order to check whether $Z_{j}$ is further changing, we calculate the potential $\psi$ for $Z=Z_{j}$, $Z_{j}-1$, and $Z_{j}+1$. The latter two are nuclei that can be produced by electron capture and emission, respectively. By comparing the values $\psi(Z_{j})$, $\psi(Z_{j}-1)$, and $\psi(Z_{j}+1)$ we choose $Z$ corresponding to the minimal $\psi$ and assign $Z_{j}=Z$. If $Z_{j}$ is modified, we repeat this stage again. If it is not modified, $Z_{j}$ is considered optimal and we proceed to the next crustal layer by applying the algorithm of stage 2.

In principle, the fourth stage, which checks for pycnonuclear reactions at the layer $j$, can be added. However, for ⁵⁶Fe ash this stage is unnecessary, because the typical $Z\approx 20$, which is realized in the inner crust in this case, is too large, leading to extremely small rate of pycnonuclear reactions.

Obviously, this algorithm should be interrupted if the crust-core boundary is reached, i.e., when, for some $\mu_{n}$, the pressure at the crust matches that in the core.

In fact, the algorithm can be interrupted even earlier because of the disintegration of all nuclei at stage 3, before reaching the crust-core boundary. This disintegration could take place as a result of the onset of an instability at $P=P_{\rm inst}$, as described in Ref. [25] for the smooth CLD model. This is the most interesting case, because the instability is required for the formation of FAC.

In terms of the potential $\psi$ the criteria for the instability can be expressed as $\psi(Z_{\rm inst})=\psi(Z_{\rm inst}-1)>\psi(Z_{\rm inst}-2)>\ldots>\psi(Z=1)>0$, where $Z_{\rm inst}$ is the charge at the previous step of the algorithm (at the layer $P<P_{\mathrm{inst}}$). In principle, pycnonuclear fusion can become important once $Z$ has decreased to a low enough value during the disintegration process. In this case, further transformation of nuclei into neutrons proceeds via an unstoppable analogue of the superthreshold electron capture cascade (SEC, [39]). This modification does not affect the final outcome – all nuclei disintegrate in the considered layer. Specifically, fusion-produced nuclei will undergo beta captures and subsequent neutron emissions, leading to low-$Z$ nuclei, which, in turn, undergo pycnonuclear fusion. An overall result of this process is disintegration of one nucleus $(Z,A)$ per one pycnonuclear reaction and production of $A$ neutrons, which redistribute over the crust to keep fixed $\mu_{n}$ in the instability layer: $(Z,A)+(Z,A)\rightarrow(2Z,2A)\rightarrow(Z,A)+An$ ($A$ is the atomic mass number). It should be noted that in the traditional approach that neglects neutron redistribution, the SEC cascade increases the number of unbound neutrons, affecting $\mu_{n}$ and preventing complete disintegration in some layers (e.g., [40, 37]). In contrast, for the nHD crust, the SEC cascade becomes unstoppable: the nHD condition fixes the neutron chemical potential and thus the amount of free neutrons in the layer. The neutrons produced by disintegration of nuclei are removed from the layer by superfluid flow or diffusion.

To illustrate the onset of the instability, let us consider the “smooth” model, i.e., by treating $Z$ as a continuous variable. The schematic $\psi(Z)$ profiles for such model are shown in Figure 1. Panel (a) shows the typical profile for layers located at $P<P_{\rm inst}$. The function $\psi(Z)$ has a profound minimum, shown by the star, where nuclei are stable. However, there is also another extremum (local maximum) at a lower $Z$. With pressure increase, the minimum and the maximum become closer and closer, and at $P=P_{\mathrm{inst}}$, they merge, producing an inflection point, shown by the star in panel (b). At this point, nuclei become unstable and undergo a series of beta captures ($Z$ is lowered), until complete disintegration into neutrons. Clearly, at the instability point, we have

At $P>P_{\mathrm{inst}}$ (panel c) $\psi$ monotonically increases with $Z$, hence no beta-stable solutions [see Eq. (10)] are available.

In a more realistic model, the shell effects introduce many local minima in the $\psi(Z)$ curve (see Figure 4). However, for sufficiently large $P$, the general slope of the function $\psi(Z)$ becomes strong enough to smooth out the local minima, leading to the instability.

As long as all nuclei in the considered volume disintegrate, the construction of the nHD accreted crust cannot be continued unambiguously to higher pressures. However, reaching $P=P_{\mathrm{inst}}$ does not necessarily mark the end of the crust. At larger pressures, $P>P_{\mathrm{inst}}$, the crust can be continued by “relic” layers that are not replaced by the accreted material in the FAC regime. These layers are formed at the initial stages of accretion and can be composed of the spherical nuclei or more complicated nuclear shapes referred to as “pasta” (see, e.g., [41, 42, 43, 44, 45, 46, 47, 48, 49, 50] for a recent discussion of pasta in pristine crust). Clearly, the unambiguous determination of the composition of these layers requires consideration of their formation, which is beyond the scope of the present paper.

When applying the CLD model without shell effects, it is natural to assume that the nuclear charge $Z$ evolves continuously in the inner crust. Of course, these assumptions will be relaxed in the next section, where the shell effects are taken into consideration. Here, we apply the CLD model suggested in Ref. [31].

The CLD model of Ref. [31] starts from an explicit expression for the free energy (see section 2 of that reference). In the zero temperature limit, which is adopted here, it reduces to the energy density $\epsilon(n_{ni},n_{pi},n_{no},n_{N},n_{e},w)$, where $n_{ni}$ and $n_{pi}$ are, respectively, the neutron and proton number densities inside nuclei; $n_{no}$ is the number density of unbound neutrons; $n_{N}$ and $n_{e}$ are, respectively, the number densities of nuclei and electrons; and $w=V_{p}n_{N}$ is the fraction of volume occupied by nuclei ($V_{p}$ is the volume of a single nucleus).

Following Ref. [25], it is useful to introduce the “total” number densities $n_{ni}^{\rm(tot)}=n_{ni}w$, $n_{pi}^{\rm(tot)}=n_{pi}w$, and $n_{no}^{\rm(tot)}=n_{no}(1-w)$, instead of $n_{ni}$, $n_{pi}$, and $n_{no}$. For example, $n_{ni}^{\rm(tot)}$ can be interpreted as the total number of neutrons in nuclei divided by the total volume (not the volume occupied by nuclei). Using these variables, the differential $d\epsilon$ can be expressed as

Introducing the baryon number density $n_{b}=n_{ni}^{\rm(tot)}+n_{pi}^{\rm(tot)}+n_{no}^{\rm(tot)}$, the neutron chemical potentials inside $\mu_{ni}=\partial\epsilon/\partial n_{ni}^{\rm(tot)}$ and outside $\mu_{no}=\partial\epsilon/\partial n_{no}^{\rm(tot)}$ nuclei, the proton chemical potential inside nuclei, $\mu_{pi}=\partial\epsilon/\partial n_{pi}^{\rm(tot)}$, as well as the chemical potentials of nuclei $\mu_{N}=\partial\epsilon/\partial n_{N}$, and electrons, $\mu_{e}=\partial\epsilon/\partial n_{e}$, Eq. (12) can be rewritten as

where we make use of the quasineutrality condition, $n_{e}=n_{pi}^{\rm(tot)}$. It is worth noting, that $\mu_{pi}$ and $\mu_{ni}$ are not equal to their bulk counterparts. In particular, both these quantities include contributions from the surface energy, and $\mu_{pi}$ additionally includes a correction associated with Coulomb energy.

To obtain EOS in the form of Section II.1, we should keep $n_{N}$, $n_{b}$, and $Z=n_{pi}^{\rm(tot)}/n_{N}$ fixed, while determining the remaining parameters $w$ and $n_{no}^{\rm(tot)}$ by minimizing $\epsilon$

These equations have a natural physical meaning: the first one represents the mechanical equilibrium of nucleus and unbound neutrons, while the second one represents the diffusion equilibrium for neutrons outside and inside nuclei (an analogue of the equilibrium with respect to neutron emission and captures). Since $\mu_{no}=\mu_{ni}$, below we refer to both of these quantities as the neutron chemical potential and denote it simply as $\mu_{n}$.

Let us consider the beta equilibrium condition (2). It is equivalent to $\partial\epsilon/\partial n_{pi}^{\rm(tot)}=0$, leading to

The equilibrium (catalyzed) crust can be obtained by imposing an additional condition,

This condition allows one to find the optimal (equilibrium) number density of nuclei, $n_{N}$.

For the accreted crust, which is considered here, additional nuclei are provided by accretion, leading to a nonequilibrium $n_{N}$ and nonequilibrium crust [25, 27]. As shown in Section II.1, for the nHD inner crust, Eq. (17) should be replaced by Eq. (5).

For the sake of completeness, let us point out that the minimization of the $\psi$ potential with respect to $Z$ can be shown to be equivalent to the beta-equilibrium condition

Summarizing, the equations (14)–(16) and (5) constitute a complete system of equations that allows one to calculate all the thermodynamic quantities if the parameter $C\equiv\mu_{N}/\mu_{n}$ and one of the thermodynamic quantities (e.g., $n_{b}$ or pressure $P$) are given.

Following Ref. [31], we derive the CLD+sh model by constraining the nuclear charge $Z$ in the smooth CLD model to integer values and incorporating precalculated shell energies. Consequently, the total energy density is expressed in the form

where we employ the notations of Section III.1; $\epsilon^{\rm CLD}$ is the energy density as it is given by the CLD model, and $\Delta\epsilon^{\rm shell}$ is a correction to the energy density, associated with the shell effects.

It is important to emphasize that in Refs. [54, 55, 33, 31] the shell corrections were primarily applied to determine $Z$ for the ground state composition, i.e., $Z$ that minimizes energy density at a fixed $n_{b}$. However, shell corrections appear to have been overlooked when calculating the “secondary” thermodynamic quantities, such as pressure. As a result, the pressure, as it was calculated in [54, 53], may not be fully thermodynamically consistent. While this inconsistency is likely insignificant for determining the catalyzed EOS, it could be important in the context of the energy release in the accreted crust discussed here. Therefore, we endeavor to avoid this inconsistency.

In general, proceeding within the CLD approach, the shell correction $\Delta\epsilon^{\rm shell}$ to the energy density should be considered as a function of five independent thermodynamic variables: $n_{b}$, $n_{N}$, $Z$, $n^{\mathrm{(tot)}}_{ni}$, and $w$. However, in this section, for the sake of simplicity, we postulate that the shell energy corrections can be expressed in the form

This assumption is made taking into account that the shell corrections are computed on top of the CLD model, which is minimized over internal CLD variables ($n^{\mathrm{(tot)}}_{no}$ and $w$) according to Eqs. (14) and (15). As a result, we arrive at the three-parameter mEOS, with the energy density $\epsilon(n_{b},n_{N},Z)$, which should be used to calculate the remaining thermodynamic quantities consistently. In particular, the pressure can be determined as

where the chemical potentials are given by

Here, as in Section II.1, the baryon chemical potential $\mu_{b}$ corresponds to the energy change due to addition of a baryon at fixed proton number density, $Zn_{N}$. Since $Zn_{N}$ is fixed, the added baryon is a neutron and $\mu_{b}$ can (again) be identified with the neutron chemical potential, $\mu_{b}=\mu_{n}$. In Eqs. (22) and (23) $\mu_{n}^{\mathrm{CLD}}$ and $\mu_{N}^{\mathrm{CLD}}$ are corresponding chemical potentials obtained from the CLD model

while corrections, associated with the shell effects are

Technically, to derive mEOS in the form $\epsilon(P,\mu_{n},Z)$, utilized in Section II.2, we employ numerical solvers to obtain the parameters $n_{b}$, $n_{N}$ for given $P$, $\mu_{n}$ and $Z$, in accordance with Eqs. (21) and (22).

In this section, we discuss the proton shell models used in this study. We recall that calculating shell effects in the inner crust is a model-dependent problem. Given the uncertainties outlined below, we consider the results obtained with shell models employed in this paper as a reasonable first step and a proof-of-principle calculation, providing a foundation for subsequent calculations with more refined models.

All our shell models are primarily based on the supplementary tables presented in Ref. [33, 53].¹¹1The tables were downloaded from https://doi.org/10.1093/mnras/stz800 on December 04, 2019. However, adapting these tables to our problem is not straightforward, and thus, we need to delve into certain technical details and assumptions to clarify their application.

Firstly, it is important to note that the tables presented as supplementary data in Ref. [53] may not be directly applicable to our case. The reported values were obtained by minimizing the ETF energy density over $n_{N}$ at fixed $Z$ and $n_{b}$, which is not appropriate for accreted crust (see Section III.1). Nevertheless, given that the resulting EOS closely resembles the catalyzed case, it seems reasonable that actual shell energies can be well approximated by assuming that $E^{\rm shell}$ does not explicitly depend on $n_{b}$, but primarily depends on the cell size, i.e., $n_{N}$. Consequently, the shell corrections are expressed in the form

where $E^{\rm shell}(n_{N},Z)$ can be obtained from the tables of Ref. [53]. Note that, in this approximation, shell corrections to the neutron chemical potential vanish, $\mu_{n}^{\mathrm{shell}}=0$ [see Equation (26)].

Secondly, it is worth mentioning that for nucleon interaction potential BSK24, applied in this work, the tables in Ref. [53] only provide shell energies for even $Z$. However, to consider beta-capture and beta-emission reactions, we need information about energies of odd-$Z$ nuclei. To analyze the role of the pairing effects, we consider two models for odd-$Z$ nuclei. In the first model, referred to as “Shell”, the shell energies for odd $Z$ are assumed to vanish. In the second model, referred to as “Shell+Pairing”, we set the shell energies for odd-$Z$ nuclei to be equal to the pairing term. The latter is estimated within the qualitative model suggested in [56] and also used by Mackie and Baym [57]:

Here $A(n_{N})$ represents the dependence of the mass number $A$ on the number density of nuclei $n_{N}$, which, for simplicity, was adopted from the results of $Z$-fixed calculations (see below).

Thirdly, we assume that shell corrections become negligible above the proton drip density. This implies that we set $E^{\rm shell}(n_{N},Z)=0$, if $n_{b}$ in the respective line of the table [53] exceeds the proton drip density $n_{b}^{p,\mathrm{drip}}$. For simplicity, we assume the proton drip density to be the same for all $Z$, $n_{b}^{p,\mathrm{drip}}=0.073$ fm^-3 [33]. For some $Z$, we apply a smooth suppression of the function $E^{\rm shell}(n_{N},Z)$ as we approach the proton drip to avoid unrealistically sharp jumps, which could impact the pressure in our thermodynamically consistent approach [see Eqs. (21)–(23)]. For numerical applications, we also eliminate some outliers from the tables in Ref. [53] and fit the remaining data on the shell energies as functions of $n_{N}$ separately for each $Z$.

Finally, it is worth noting that the shell energy tables in [53] also do not contain data for $Z<18$. For such nuclei, the shell energy was assumed to be zero, which most certainly has a negligible effect on our calculations because formation of nuclei with $Z<18$ is blocked by high potential $\psi$ for $Z=18$ nuclei (see Figure 4).

We also employ the third (simplified) model as a sensitivity test for our results. In this model we assume that the charge number $Z$ remains fixed, $Z=20$, up to the proton drip, effectively mimicking very strong shell and pairing effects. While, strictly speaking, the dependence of the shell energy on $n_{N}$ influences the results by affecting the pressure, we simplify this model by assuming $E^{\rm shell}(n_{N},Z=20)=0$. The results for the $Z$-fixed model are qualitatively aligned with both the Shell and Shell+Pairing models. The quantitative differences can be used to estimate uncertainties associated with shell effects (see Sections IV.2 and VII for details).

Let us begin the discussion of the results for the smooth CLD model by examining the profiles of the potential $\psi$ as a function of $Z$ for five values of pressure. The panels (from top to bottom) correspond to $P=0.085$, $0.22$, $0.25$, $0.27$, and $0.28$ MeV fm^-3, respectively. In each panel, three members of the nHD EOS family are shown: $P_{\mathrm{oi}}=P_{\mathrm{nd}}^{\mathrm{(cat)}}$ (solid line), $P_{\mathrm{oi}}=1.005\,P_{\mathrm{nd}}^{\mathrm{(cat)}}$ (dotted line), and $P_{\mathrm{oi}}=P_{\mathrm{oi}}^{\mathrm{(min)}}\approx 1.011\,P_{\mathrm{nd}}^{\mathrm{(cat)}}$ (dashed line). Here $P_{\mathrm{nd}}^{\mathrm{(cat)}}$ is the pressure at the oi interface for catalyzed crust, which coincides with the neutron drip pressure (see [58] for discussion). Following [27], we indicate this by the subscript “nd”.

In the smooth CLD model, the nHD EOS for the entire crust can be specified by setting the pressure $P_{\mathrm{oi}}$ at the outer-inner crust interface. Consequently, the solid line, corresponding to $P_{\mathrm{oi}}=P_{\mathrm{nd}}^{\mathrm{(cat)}}$, represents the $\psi$ profiles for catalyzed crust. At each pressure the profiles exhibit profound minima, as clearly visible in Figure 2. These minima correspond to the catalyzed crust composition ($\psi=\mu_{N}=0$ at the minima; see Section II.1), indicating the absence of disintegration instability. It is important to note, however, that all panels, except the top one, correspond to pressures where the pasta phases are more energetically favorable in the CLD model, as reported in Ref. [47]. Hence, nuclei with large $Z$ ($Z>50$) are likely absent in the catalyzed crust, being shown here only for illustrative purposes.

The dotted line represents a slightly higher value of $P_{\mathrm{oi}}=1.005\,P_{\mathrm{nd}}^{\mathrm{(cat)}}$, which, however, is not sufficient to induce instability. The $\psi$ profiles still exhibit minima at all pressures, although these minima are less pronounced than for solid curves ($P_{\mathrm{oi}}=P_{\mathrm{nd}}^{\mathrm{(cat)}}$). Nonetheless, they still correspond to a stable composition. Therefore, the FAC model would require a larger $P_{\mathrm{oi}}$ to reach instability.

The dashed line corresponds to $P_{\mathrm{oi}}=P_{\mathrm{oi}}^{\mathrm{(min)}}\approx 1.011\,P_{\mathrm{nd}}^{\mathrm{(cat)}}$, the lowest $P_{\mathrm{oi}}$ value that leads to disintegration instability in the nHD inner crust. For this $P_{\mathrm{oi}}$, the minima in the $\psi$ profiles become progressively shallower with increasing pressure and eventually disappear, transforming into an inflection point in the panel corresponding to $P=0.27$ MeV fm^-3 (the second panel from the bottom). At this pressure, the nuclei disintegrate, allowing for the existence of stationary accreted crust. However, it is worth noting that the neutron chemical potential at this point does not match the core $\mu_{n}$ (at the same pressure), and thus, the crust must be extended to higher pressures. These layers can be called “relic” [27], since the accreted nuclei do not penetrate into these layers (they disintegrate earlier). The relic layers can be filled with spherical nuclei formed during the initial stages of accretion. Indeed, at a higher pressure, $P=0.28$ MeV fm^-3 (the bottom panel), a minimum in the $\psi(Z)$ profile for the dotted curve reappears, suggesting that stable relic crust can exist. Alternatively, the relic layers could be filled with pasta, or a layer of relic spherical nuclei could be followed by a pasta layer.²²2To avoid any confusion, by “spherical nuclei” we mean, following, e.g., Ref. [33] nuclei whose energy is calculated within the ETFSI approach using a spherically symmetric single-particle effective potential. A detailed study of the structure of relic layers is beyond the scope of the present work.

For a smooth CLD model, the nHD solutions with $P_{\mathrm{oi}}>P_{\mathrm{oi}}^{\mathrm{(min)}}$ contain an unstable region where $\Psi/N_{\mathrm{N}}$ is a monotonically increasing function of $Z$. This region is located between the accreted and relic parts of the crust, indicating that the crustal model with $P_{\mathrm{oi}}>P_{\mathrm{oi}}^{\mathrm{(min)}}$ is thermodynamically inconsistent and cannot be applied to describe neutron stars [25, 27]. Therefore, the nHD model with $P_{\mathrm{oi}}=P_{\mathrm{oi}}^{\mathrm{(min)}}$ is the only possible stationary crust model, and it is referred to as simply the “FAC model” in what follows.

Figure 3 compares two EOSs: the catalyzed crust and the FAC model. The density and composition of catalyzed crust are shown by solid lines for $P<P_{\mathrm{pasta}}\approx 0.13$ MeV fm ^-3 and dotted lines for $P\geq P_{\mathrm{pasta}}$. The dotted lines are used to indicate that the region $P\geq P_{\mathrm{pasta}}$ may actually contain strongly nonspherical (pasta-like) structures [47], although the curves were calculated assuming spherical nuclei.

The density and composition of the FAC model are shown in Figure 3 by thick dashes. The instability takes place at $P_{\mathrm{inst}}\approx 0.27$ MeV fm^-3, marked with asterisks in each panel of Figure 3. As discussed above, $P_{\mathrm{inst}}$ is lower than the crust-core boundary and the crust should continue with relic layers. In Figure 3, the density and composition of relic layers are represented by thin dash-dotted lines, assuming that these layers consist of spherical nuclei. It is worth noting that for the smooth CLD model, the composition of these layers can be unambiguously determined for a given $P_{\mathrm{oi}}$, as the nHD equilibrium of relic layers allows us to apply formulas from Section II.1. However, if shell effects are included, determining the composition of relic layers generally requires consideration of their formation history from the beginning of accretion process.

As seen from Figure 3, both density and composition profiles of the FAC state closely resemble those of the catalyzed crust. The fact that the instability occurs at a higher density than the crust-pasta transition, as predicted by [47], suggests that the thickness of the pasta layers may be affected by accretion. However, we defer a more detailed discussion of these effects and the role of the pasta in the accreted crust to subsequent studies.

The upper panel in Figure 3 illustrates the difference in $\mu_{n}(P)$ between the catalyzed crust and FAC model. It is evident that the difference is quite small ($\lesssim 20$ keV) and diminishes with increasing pressure. As in the other panels, the dotted part of the curve is calculated assuming that the catalyzed crust is composed of spherical nuclei at such densities.

As for the smooth CLD model, let us begin the discussion of the results by examining the profiles of the potential $\psi$, which are shown in Figure 4. In contrast to Figure 2, the lines in each panel correspond to different shell models (see Section III.3): $Z$-fixed (solid line), Shell (long dashes), and Shell+Pairing (short dashes). We have chosen $P_{\mathrm{oi}}$ to be equal to $P_{\mathrm{oi}}^{\mathrm{(min)}}$, which is specific for each model ($P^{\mathrm{(min)}}_{\mathrm{oi}}\approx 0.432$ keV fm^-3 for Shell and Shell+Pairing models; $P^{\mathrm{(min)}}_{\mathrm{oi}}\approx 0.424$ keV fm^-3 for the $Z$-fixed model). The panels, from top to bottom, correspond to $P=0.085$, 0.16, 0.22, 0.25 MeV fm^-3 and $P=P^{\mathrm{inst}}$. The charge of the nuclei is denoted by a cross (for each model in each panel). As expected, for Shell and Shell+Pairing models, it corresponds to a local minimum of the potential $\psi$ in each panel. The exception is the bottom panel, where the minimum disappears and becomes an unstable inflection point. For the $Z$-fixed model, the charge is artificially fixed at $Z=20$ up to the proton drip (three upper panels). That is, $Z=20$ does not correspond to the minimum of the $\psi$ potential.

In two upper panels ($P=0.085$ and $0.16$ MeV fm^-3) a profound minimum of the $\psi$ potential is formed by the shell effects for both Shell and Shell+Pairing models. This guarantees the conservation of the nuclear charge number since beta-capture and beta-emission reactions are not energetically favourable. For the $Z$-fixed model $Z$ is fixed, thus for all considered models, $Z$ is equal to 20 at these pressures.

In the third panel ($P=0.22$ MeV fm^-3) the charge number starts to differ. For the $Z$-fixed model, it is kept equal to 20 by construction, while for the Shell model it evolves to $Z=35$. This is because the shell effects for low-$Z$ nuclei become small (according to the applied model), and the shell structure cannot prevent beta-reactions driven by the general trend of $\psi$ provided by the CLD part of the model. The proton drip has not yet been reached, and we retain pairing effects in the Shell+Pairing model. They prevent a strong increase in $Z$; however, $Z$ is increased to $Z=22$ due to a pair of beta emissions, which occur between $P=0.16$ MeV fm^-3 and $P=0.22$ MeV fm^-3, when combination of the shell and pairing effects remove a $\psi$-potential barrier for the transition to $Z=21$ nuclei via beta emission.

The fourth panel, $P=0.25$ MeV fm^-3, corresponds to the matter after the proton drip. Within our approximation, the shell and pairing corrections vanish, and $\psi(Z)$ becomes a smooth function of $Z$, as given by the CLD model. The composition is driven to the minima located at $Z=32$, 34, and 35 for Shell, Shell+Pairing, and $Z$-fixed models, respectively. The $\psi$ profiles become rather close for all the considered models. Note, however, that they do not coincide exactly because the $\psi$ potential depends not only on $P$, but also on $\mu_{n}$, which is specific for each of the considered shell models due to differences in $P_{\mathrm{oi}}^{\mathrm{(min)}}$ and in the respective EOSs. This feature leads to differences in nuclear evolution after the proton drip (see Figure 5 below).

In the bottom panel, where $P=P^{\mathrm{inst}}$, there are no minima for all models; the “optimal” $Z$ correspond to the inflection point. Nuclei become unstable and disintegrate through a sequence of beta captures, driving them to lower $Z$ and $\psi$.

Let us now turn to discussing the evolution of nuclei in the inner crust. For $Z$-fixed and Shell+Pairing models, $Z$ at the oi interface starts from the value $Z=20$ for all the considered members of the nHD family of crust models, corresponding to different $P_{\mathrm{oi}}$. These nuclei are formed by reactions in the outer crust and at the oi interface. Subsequent compression up to $P\lesssim 0.15$ MeV fm^-3 does not lead to beta-reactions and the charge number remains to be $Z=20$. No heat is released. The evolution consists of continuous growth of nuclear mass numbers $A$ and $A_{c}$ with pressure, $P$. Neutrons required for this growth are provided by diffusion/superfluid flow from the higher-density regions. The evolution for the Shell model is generally the same, but $Z=21$ nuclei are formed at the oi interface and converted into $Z=20$ nuclei by electron capture at a slightly higher pressure. This is a clear artifact of neglecting the pairing effects in that model. We expect that it does not significantly affect our results. In particular, the electron capture does not lead to any energy release, because it takes place exactly at the threshold.

Nuclear evolution in the inner layers of the inner crust is more interesting (see Figure 5). Panels in that figure, from top to bottom, present the accumulated heat $Q_{\mathrm{i}}$ in the inner crust (per accreted baryon, not including the heat released at the oi interface), as well as $A_{c}$, $A$, and $Z$ as functions of $P$. Each panel contains two types of curves for each model (Z-fixed, Shell, and Shell+Pairing): the thin curve shows the evolution for $P_{\mathrm{oi}}=P^{\mathrm{(min)}}_{\mathrm{oi}}$, while the thick curve is for $P_{\mathrm{oi}}=P^{\mathrm{(cat)}}_{\mathrm{nd}}$. In the lower panels, all curves coincide at low pressure. In the upper panel, the energy release for $P_{\mathrm{oi}}=P^{\mathrm{(cat)}}_{\mathrm{nd}}$ is divided by a factor of 5 to fit the scale of the plot. The instability points are shown by asterisks (except in the upper panel).

Let us first consider the case $P_{\mathrm{oi}}=P^{\mathrm{(cat)}}_{\mathrm{nd}}$. Then, for Shell and Shell+Pairing models, the instability occurs for $Z=20$ nuclei before the proton drip. This is due to the high value of $P_{\mathrm{oi}}$, for which the general trend in the $\psi(Z)$ dependence becomes so strong that the shell effects fail to prevent beta captures. For the $Z$-fixed model, the instability occurs exactly at the proton drip point, simply because we do not allow $Z$ to vary at lower $P$ in this model. All the heat $Q_{\mathrm{i}}$ is released at the instability point, resulting in a jump in the $Q_{\mathrm{i}}(P)$ dependence, as shown in the upper panel.

For $P_{\mathrm{oi}}=P^{\mathrm{(min)}}_{\mathrm{oi}}$, the nuclear evolution for the considered shell models is somewhat different. However, the general trend is a growth of $Z$ up to $33-35$ with subsequent decrease before the onset of instability. The increase in $Z$ is accompanied by heat release; however, the released heat at the instability is dominant for all models.

Let us discuss the details of nuclear evolution, starting with the Shell model. The first electron emission takes place at $P\approx 0.17$ MeV fm^-3, when the shell effects for nuclei with $Z=20$ become small and cannot prevent beta capture. The absence of pairing correction in this model allows for the formation of $Z=21$ nuclei. Subsequent electron captures are associated with a further decrease of the shell corrections, leading to gradual disappearance of local minima of the function $\psi(Z)$. Finally, at $P\approx 0.20$ MeV fm^-3, the shell corrections vanish for low-$Z$ nuclei, and the nuclear charge arrives at a minimum of $\psi(Z)$, determined by the smooth CLD model. The negative shell energy corrections for high-$Z$ nuclei constitute the global minimum of $\psi(Z)$; however, this minimum remains unattainable. This pattern stays the same with a subsequent increase in pressure, and the proton drip does not leave any imprints on the evolution, since the local shape of $\psi(Z)$ is already driven by the CLD model. The decrease of $Z$ before the instability onset is associated with changes in the $\psi(Z)$ shape with growing $P$.

Within the Shell+Pairing model, the decrease of shell corrections at $P\approx 0.17$ MeV fm^-3 does not affect the evolution because the local minimum at $Z=20$ is well defined by the pairing correction. Subsequent compression leads to a pair of electron emissions and the associated energy release from the second capture. However, the pairing correction adopted in the Shell+Pairing model is strong enough to form a local minimum at $Z=22$ and prevent subsequent electron emissions until the proton drip takes place. At the proton drip, we suppress the pairing corrections along with shell corrections, and the $\psi(Z)$ profile becomes determined by the CLD model. The nuclear charge reaches a minimum of $\psi(Z)$ potential, located at $Z=34$. During subsequent compression, $Z$ follows the position of the minimum (see the second from the bottom panel in Figure 4).

The $Z$-fixed model is applied as a sensitivity test for shell effects. It prevents any evolution of the nuclear charge $Z$ until the proton drip point. For higher pressure, as in the Shell and Shell+Pairing models, the $\psi(Z)$ profile becomes determined by the CLD model, which drives the subsequent evolution.

Figure 6 presents the dependence of $P_{\mathrm{inst}}$ on $P_{\mathrm{oi}}$ for the considered models. One can observe that $P_{\mathrm{inst}}$ decreases with an increase in $P_{\mathrm{oi}}$ for all models.

Let us begin the discussion of this dependence with the $Z$-fixed model. As seen in Figure 5, for $P_{\mathrm{oi}}=P_{\mathrm{oi}}^{\mathrm{(min)}}$, the instability occurs after the proton drip. An increase in $P_{\mathrm{oi}}$ shifts the instability point to lower pressure, closer to the proton drip, located at $P\approx 0.231$ MeV fm^-3. Finally, at $P_{\mathrm{oi}}\approx 0.425$ keV fm^-3, the instability starts at the proton drip. In the $Z$-fixed model, we do not allow $Z$ to evolve before the proton drip. As a result, the instability cannot occur at a lower pressure. Consequently, for $P_{\mathrm{oi}}>0.425$ keV fm^-3, the instability always takes place at the proton drip point (some fluctuations in Figure 6 are numerical noise). This behavior is an artifact of the simplified $Z$-fixed approach, prompting us to switch to more realistic models.

For the Shell model, as we discussed above, even for $P_{\mathrm{oi}}=P_{\mathrm{oi}}^{\mathrm{(min)}}$, the nuclear charge starts to evolve before the proton drip point, and the proton drip does not affect the evolution. This statement holds true for higher $P_{\mathrm{oi}}$ values, and $P_{\mathrm{inst}}(P_{\mathrm{oi}})$ is a monotonically decreasing function with no peculiarities at the proton drip point. The break at $P_{\mathrm{oi}}\approx 0.44$ keV fm^-3 is associated with a change in the nuclear evolution behavior: for lower $P_{\mathrm{oi}}$, the nuclear charge starts to grow (above some pressure), and the instability occurs for $Z>20$, while for higher $P_{\mathrm{oi}}$, the instability occurs for $Z=20$ nuclei without any charge evolution before it (see the thick line for $P_{\mathrm{oi}}=P_{\mathrm{nd}}^{\mathrm{(cat)}}$ and the thin line for $P_{\mathrm{oi}}=P_{\mathrm{oi}}^{\mathrm{(min)}}$ in Figure 5).

The behavior of $P_{\mathrm{inst}}(P_{\mathrm{oi}})$ is a bit more complicated for the Shell+Pairing model. As with the Shell model, for $P_{\mathrm{oi}}=P_{\mathrm{oi}}^{\mathrm{(min)}}$, the nuclear charge starts to evolve before the proton drip, but this evolution is not too strong; pairing effects prevent $Z$ from growing above $Z=22$ (see Figure 5). At the proton drip point, we remove shell and pairing corrections, and the nuclear charge increases.

Increasing $P_{\mathrm{oi}}$ results in a decrease in $P_{\mathrm{inst}}$, while qualitatively nuclear evolution remains unchanged: The nuclear charge increases at the proton drip point and (generally) evolves a little further before the instability onset. However, at $P_{\mathrm{oi}}\approx 0.433$ keV fm^-3, the instability point reaches the proton drip, where instability leads to the disintegration of $Z=22$ nuclei. Subsequent growth of $P_{\mathrm{oi}}$ up to $0.44$ keV fm^-3 does not affect $P_{\mathrm{inst}}$, because the pairing effects prevent disintegration of $Z=22$ nuclei before the proton drip. However, as with the Shell model, for $P_{\mathrm{oi}}>0.44$ keV fm^-3, the instability occurs before the proton drip in the form of disintegration of $Z=20$ nuclei, without an increase of nuclear charge before the instability. In this high-$P_{\mathrm{oi}}$ region, the value of $P_{\mathrm{inst}}$ is slightly higher for the Shell+Pairing model than for the Shell model because the pairing correction provides additional stabilization for $Z=20$ nuclei.

For each shell model, one can determine $Q_{\mathrm{inst}}^{\mathrm{(min)}}$ – the minimum possible $Q_{\mathrm{inst}}$. According to Figure 7, this corresponds to the highest possible pressure, $P^{\mathrm{(max)}}_{\mathrm{inst}}$, and numerically $Q_{\mathrm{inst}}^{\mathrm{(min)}}$ is the same for all considered shell models, $Q_{\mathrm{inst}}^{\mathrm{(min)}}\approx 0.05$ MeV per accreted nucleon.

This latter statement should hold true for any shell model added on top of a given smooth CLD model, provided that shell effects are negligible at high pressures. The proof of this crucial statement is, in fact, straightforward: if shell effects are neglected at high pressures, the nuclear physics there is specified by the smooth model⁴⁴4It can be either a CLD model (as in the present work) or, for example, a more sophisticated smooth model based on the ETF calculations [33, 53]. . As previously stated, in the smooth model the amount of energy release is determined by $P_{\mathrm{inst}}$. Thus, the function $Q_{\mathrm{inst}}(P_{\mathrm{inst}})$ and its minimum, $Q_{\mathrm{inst}}^{\mathrm{(min)}}$, are not affected by the shell effects.

In principle, $Q_{\mathrm{inst}}^{\mathrm{(min)}}$ can depend on the applied smooth model. The detailed analysis of this dependence is left for future work. Here we would like to point out that for another Skyrme-type nuclear potential (Sly4 [59]), the smooth CLD model leads to a very similar value for the minimal heat release ($0.04$ MeV per accreted nucleon, see supplementary material in Ref. [27]). This suggests that $Q_{\mathrm{inst}}^{\mathrm{(min)}}\sim(0.04\div 0.05)$ MeV per accreted nucleon can be considered a rather model-independent estimate.

In Ref. [27], we constrained $P_{\mathrm{oi}}$ by imposing a condition that $Q^{\mathrm{AMT}}_{\mathrm{i}}(P_{\mathrm{oi}})>0$. The corresponding lower bound for the pressure was denoted as $P_{\mathrm{oi}}^{(0)}$. Here we suggest applying a tighter constraint: $Q^{\mathrm{AMT}}_{\mathrm{i}}(P_{\mathrm{oi}})>Q_{\mathrm{inst}}^{\mathrm{(min)}}$, where $Q_{\mathrm{inst}}^{\mathrm{(min)}}$ represents the minimal energy release at the instability. As discussed in Section VI.1, for BSK24 and SLy4 models $Q_{\mathrm{inst}}^{\mathrm{(min)}}$ can be estimated as $Q_{\mathrm{inst}}^{\mathrm{(min)}}\sim(0.04\div 0.05)$ MeV per accreted nucleon, provided that shell corrections are negligible in the deepest layers of the inner crust. By using the dependence $Q^{\mathrm{AMT}}_{\mathrm{i}}(P_{\mathrm{oi}})$ calculated in Ref. [27] (see inset in Figure 8, $Q_{\mathrm{inst}}^{\mathrm{(min)}}\sim(0.04\div 0.05)$ MeV is filled in gray in the inset) for the BSK24 HFB mass tables [61], we obtain $P_{\mathrm{oi}}\geq 0.46$ keV fm^-3. However, this constraint is only slightly stronger (by $\sim 2\%$) than the originally suggested lower bound $P_{\mathrm{oi}}^{(0)}$ [27]. Note, however, that a detailed analysis within the advanced shell model will likely lead to the same result.

As discussed in Section VI.1, $Q_{\mathrm{inst}}^{\mathrm{(min)}}$ is determined by the CLD model and thus should not depend on the initial composition. Thus, the approach discussed in this subsection is likely applicable to accreted crust with a multicomponent composition of initial ashes. Note, however, that in the latter case, noticeable fraction of the heat can be released in the intermediate layers of the inner crust [29, 30], so that the constraint should be applied to the residual part of the heat release, which also can be estimated within the AMT approach [29, 30].

Here we extend the approach from Section VII.1 in order to heuristically predict the dependence of $P_{\mathrm{inst}}(P_{\mathrm{oi}})$ for the advanced shell model. As before, we employ the energy conservation law. In addition, we make use of the assumption of constancy of $Z$ for the majority of the inner crust, up to the innermost regions, where the advanced shell model aligns with the respective CLD+sh model. This assumption is supported by our calculations (see Sections IV.2 and VI) and can be considered as an applicability condition for the approach proposed in this subsection.

As pointed out above, $Q^{\mathrm{AMT}}_{\mathrm{i}}(P_{\mathrm{oi}})$ calculated following Ref. [27], should be equally valid for the advanced shell model, i.e., $Q_{\rm i}(P_{\rm oi})=Q^{\mathrm{AMT}}_{\mathrm{i}}(P_{\mathrm{oi}})$. At the same time, one can calculate $Q_{\mathrm{i}}(P_{\mathrm{inst}})$ for a certain CLD+sh model (see Section VI). If all heat release in the inner crust is concentrated in its deepest layers, it is reasonable to assume that $Q_{\mathrm{i}}(P_{\mathrm{inst}})$ is fully determined by the nuclear physics in the deepest layers of the inner crust (see SectionVI), and thus $Q_{\mathrm{i}}(P_{\mathrm{inst}})$ as calculated for the CLD+sh model, should also be valid for the advanced shell model.

Using the (known) functions $Q_{\mathrm{i}}(P_{\mathrm{inst}})$ and $Q_{\mathrm{i}}(P_{\mathrm{oi}})=Q^{\mathrm{AMT}}_{\mathrm{i}}(P_{\mathrm{oi}})$, it is straightforward to find $P_{\rm inst}$ as a function of $P_{\rm oi}$. This prediction does not require any additional calculation for the advanced shell model; it indicates that the function $P_{\rm inst}(P_{\rm oi})$ does not depend on the form of the shell corrections in the interiors of the inner crust, provided that the shell effects are strong enough to prevent the evolution of $Z$ (and thus heat release) up to the deepest layers of the inner crust.

To illustrate this approach, we combined the function $Q_{\mathrm{i}}(P_{\mathrm{inst}})$ calculated in this work (Figure 7) with the function $Q_{\mathrm{i}}^{\mathrm{AMT}}(P_{\mathrm{oi}})$ calculated for the BSK24 HFB mass tables [61] in Ref. [27] (see inset in Figure 8). The resulting heuristic prediction for the corrected dependence $P_{\mathrm{inst}}(P_{\mathrm{oi}})$ is shown in Figure 8.

Comparing Figures 6 and 8, one can see that for low $P_{\mathrm{oi}}$ the corrected predictions for the $Z$-fixed model are closer to the Shell and Shell+Pairing models than in the case of our original calculations in Section IV.2. This is because the nuclear physics input for CLD+sh versions of these models differs even in the outermost layers of the inner crust. Specifically, as discussed in Section III.2, we neglect shell corrections for the $Z$-fixed model. As a result, $Q_{\mathrm{i}}(P_{\mathrm{oi}})$ calculated within the $Z$-fixed model differs from calculations for Shell and Shell+Pairing models simply because the mass of $Z=20$ nuclei predicted by the $Z$-fixed model differs from the mass of the same nuclei, predicted by Shell and Shell+Pairing models. The corrected results rely on the more accurate nuclear physics in the shallow inner crust region and should be considered more reliable.