A SURVEY OF CHEMICAL SEPARATION IN ACCRETING NEUTRON STARS

We adopt realistic mixtures that result from calculations of rp-process burning under different conditions. An important factor is whether the hydrogen and helium burning is thermally stable or unstable. Schatz et al. (1999) calculated the ashes of stable burning, and showed that the composition for accretion rates $\gtrsim {\dot{m}}_{{\rm{Edd}}}$ has a range of mean nuclear charge from $\langle Z\rangle \;\approx$ 26–50. We focus on the nuclear charge Z rather than mass A because that is the relevant quantity for freezing of a Coulomb liquid. The composition becomes heavier at larger accretion rates, up to a limit at $\approx 50{\dot{m}}_{{\rm{Edd}}}$ (where ${\dot{m}}_{{\rm{Edd}}}$ is the local Eddington accretion rate) beyond which the mixture reaches nuclear statistical equilibrium, driving the mixture back to iron group (this effect can be seen in Figures 7 and 10 of Schatz et al. 1999). Schatz et al. (2003) extended these calculations to lower accretion rates below ${\dot{m}}_{{\rm{Edd}}}$ and showed that the composition would be much lighter, with $\langle Z\rangle \approx 10$ at $\dot{m}\approx 0.1\ {\dot{m}}_{{\rm{Edd}}}$. Unstable burning gives a heavier composition than stable burning because it occurs at a significantly larger temperature, resulting in heavy ashes beyond the iron group with $A\approx 60$–100 (Schatz et al. 2003; Woosley et al. 2004).

Stevens et al. (2014) recently calculated the rp-process ashes for a large number of steady-state burning models with different accretion rates, base fluxes, and helium fraction in the accreted material. Here we study four models from that paper which have helium-rich accreting material with helium mass fraction Y = 0.55 and accretion rates $\dot{m}/{\dot{m}}_{{\rm{Edd}}}=0.1,0.5,1.0,2.0$ and then a series of models with Y = 0.2752 and $\dot{m}/{\dot{m}}_{{\rm{Edd}}}$ from 0.1 to 40. We also look at the composition used by Horowitz et al. (2007), taken from Gupta et al. (2008), which was the result of unstable H/He burning, and an additional composition resulting from unstable burning. These compositions span a nuclear charge range of $\langle Z\rangle =8$–35.

The rp-process ashes are produced well above the ocean floor, at densities $\sim {10}^{5}$–${10}^{6}\;{\rm{g}}\ {\mathrm{cm}}^{-3}$. We allow for beta decays and electron captures as the ashes move to higher densities by finding the Z for each mass chain A that is in beta equilibrium at an electron Fermi energy ${E}_{{\rm{F}}}=4\;{\rm{MeV}}$, which is a typical Fermi energy at the freezing depth (for a one-component plasma, the Fermi energy at the freezing depth is ${E}_{{\rm{F}}}=1.7\;{\rm{MeV}}\ {T}_{8}{(Z/26)}^{-5/3}$, where ${T}_{8}=T/{10}^{8}\;{\rm{K}}$). We have checked that our results are not sensitive to the exact choice of E_F.

Some of the compositions contain substantial amounts of the light elements helium and carbon, which are likely to burn before reaching the base of the ocean where chemical separation occurs. At lower accretion rates the validity of this assumption is much more uncertain and carbon could potentially reach the ocean basin possibly sourcing superbursts (Strohmayer & Brown 2002). Recent multi-zone modeling of superbursts found additional crustal heating (Q_b) to be necessary to transition from unstable to stable burning at low accretion rates (Keek & Heger 2011). This, in conjunction with rp-process ashes with large carbon mass fractions, makes the assumption of complete and stable burning of carbon even more uncertain at low accretion rates. The interaction of chemical separation with rp-process ashes of these low accretion rates will be discussed further in Section 3. Unless otherwise stated, we assume that light element burning converts the He and C into Mg (Z = 12), and so set the He and C number fractions to zero, and add the sum of the original He and C number fractions to the Mg number fraction.⁴

For each composition, we then use the semi-analytic approach of Medin & Cumming (2010) to find the composition of the liquid and solid that are in equilibrium with each other. Full details can be found in that paper, but for clarity we give a brief summary of the method here. We start with the linear mixing rule for the free energy of a mixture

$$\begin{eqnarray}&&{f}^{\mathrm{LM}}\quad =\quad \displaystyle \sum _{j}{x}_{j}\left[{f}^{\mathrm{OCP}}({{\rm{\Gamma }}}_{j})+\mathrm{ln}\left(\displaystyle \frac{{x}_{j}{Z}_{j}}{\langle Z\rangle }\right)\right],\end{eqnarray} \tag{ 1 }$$

where ${f}^{\mathrm{OCP}}$ is the free energy of a one-component plasma, and the logarithmic term is the entropy of mixing. The Coulomb parameter ${{\rm{\Gamma }}}_{j}$ is given by ${{\rm{\Gamma }}}_{j}={Z}_{j}^{5/3}{{\rm{\Gamma }}}_{e}$ with ${{\rm{\Gamma }}}_{e}={e}^{2}/{a}_{e}{k}_{B}T$, where a_e is the mean electron spacing ${a}_{e}={[3/(4\pi {n}_{e})]}^{1/3}$ at the local electron density n_e. While the linear mixing rule is adequate to describe the free energy of the liquid phase, the solid phase requires a correction term be included, and so we write the liquid and solid free energies as

$$\begin{eqnarray}&&{f}_{l}\quad =\quad {f}_{l}^{\mathrm{LM}}\phantom{\rule{1cm}{0ex}}{f}_{s}\quad =\quad {f}_{s}^{\mathrm{LM}}+{\rm{\Delta }}{f}_{s}.\end{eqnarray} \tag{ 2 }$$

Extending the work of Ogata et al. (1993) on the three-component plasma to an arbitrary number of m components, Medin & Cumming (2010) write ${\rm{\Delta }}{f}_{s}$ as a sum over pairwise interactions,

$$\begin{eqnarray}&&{\rm{\Delta }}{f}_{s}\quad =\quad \displaystyle \sum _{j=1}^{m-1}\;\displaystyle \sum _{k=j+1}^{m}{{\rm{\Gamma }}}_{j}{x}_{j}{x}_{k}{\rm{\Delta }}g\left(\displaystyle \frac{{x}_{k}}{{x}_{j}+{x}_{k}},\displaystyle \frac{{Z}_{k}}{{Z}_{j}}\right),\end{eqnarray} \tag{ 3 }$$

where the function ${\rm{\Delta }}g$ is taken from Ogata et al. (1993).

Once the free energy is obtained, we then look for an $m-1$-dimensional plane in the m-dimensional space of composition that is tangent to the free energy surfaces for the liquid and solid. This m-dimensional version of the usual tangent construction then allows one to decompose any composition that lies between the liquid and solid tangent points into coexisting phases of that liquid and solid.

We calculate the liquid and solid compositions that are in equilibrium for a mixture of 50% liquid and 50% solid. We calculate the chemical separation using the 17 most abundant species from each rp-process mixture. The choice of 17 species is for convenience, since it matches the number of species in the composition used by Horowitz et al. (2007) and Medin & Cumming (2010). However, we have checked that using a smaller or larger number of species in the calculation of chemical separation does not significantly change the results.

The results for all mixtures studied in this paper ($Y=0.55,0.2752$ and two X-ray burst compositions) are summarized in Table 1 which gives the mean charge $\langle Z\rangle$ and impurity parameter in the initial mixture, liquid, and solid for each case. The mean charge is plotted in Figure 1. In agreement with the results of Horowitz et al. (2007) for a single mixture, we find the general result that $\langle Z\rangle$ increases for the solid and decreases for the liquid relative to the initial composition.⁵ This implies that the equilibrium liquid and solid phases will be preferentially enriched in lighter and heavier nuclei, respectively, relative to the initial mixture.

Zoom In Zoom Out Reset image size

Table 1.  Properties of the Liquid–Solid Equilibrium Compositions

$\dot{m}/{\dot{m}}_{{\rm{Edd}}}$	$\langle Z{\rangle }_{i}$	$\langle Z{\rangle }_{l}$	$\langle Z{\rangle }_{s}$	Qi	Ql	Qs	${Y}_{e,i}$	${Y}_{e,l}$	${Y}_{e,s}$	${Y}_{e,l}-{Y}_{e,i}$	${{\rm{\Gamma }}}_{i}$	${{\rm{\Gamma }}}_{l}$	${{\rm{\Gamma }}}_{s}$
Y = 0.2752
0.1	11.4	11.6	11.2	8.3	13.5	3.1	0.4600	0.4592	0.4609	−0.0008	286.5	299.5	273.5
0.2	16.4	13.2	19.7	41.2	23.1	38.4	0.4483	0.4581	0.4419	0.0098	509.3	348.4	670.2
0.3	18.8	12.8	24.8	48.1	19.1	4.6	0.4460	0.4655	0.4367	0.0194	556.7	288.8	824.5
0.4	21.6	16.9	26.3	51.2	50.1	8.2	0.4434	0.4582	0.4343	0.0148	351.5	240.6	462.2
0.5	22.5	18.3	26.8	49.8	55.3	8.7	0.4405	0.4535	0.4319	0.0130	318.4	232.8	404.0
0.6	23.1	19.1	27.1	49.7	58.3	8.8	0.4388	0.4507	0.4307	0.0119	305.3	229.4	381.2
0.7	23.4	19.4	27.5	50.6	60.4	7.9	0.4380	0.4496	0.4302	0.0116	303.1	228.4	377.8
0.8	23.7	19.6	27.7	51.8	62.3	8.0	0.4376	0.4490	0.4300	0.0114	301.4	227.4	375.3
0.9	23.8	19.8	27.9	52.7	64.0	8.2	0.4372	0.4484	0.4297	0.0111	298.5	226.2	370.9
1.0	24.0	20.0	28.0	52.8	65.1	8.5	0.4368	0.4475	0.4295	0.0107	294.6	225.1	364.0
1.26	24.4	20.5	28.3	53.1	66.9	9.3	0.4363	0.4462	0.4294	0.0099	289.0	224.3	353.6
1.58	24.7	21.0	28.5	52.6	67.6	9.4	0.4363	0.4455	0.4297	0.0092	285.3	224.3	346.2
2.0	25.1	21.6	28.6	51.1	67.2	10.8	0.4363	0.4443	0.4304	0.0081	278.5	224.2	332.7
2.51	25.5	22.3	28.7	48.4	65.1	11.3	0.4362	0.4433	0.4309	0.0070	272.3	224.2	320.3
3.16	26.0	23.2	28.7	44.1	61.2	12.0	0.4357	0.4414	0.4312	0.0057	260.0	221.3	298.6
4.0	26.7	24.5	28.8	38.1	54.0	13.0	0.4351	0.4394	0.4316	0.0042	245.6	217.6	273.6
5.0	27.5	25.9	29.2	34.1	48.6	14.1	0.4334	0.4357	0.4315	0.0022	234.1	214.2	253.9
6.31	28.7	27.4	30.0	32.4	45.6	16.0	0.4311	0.4320	0.4304	0.0008	253.8	247.0	260.5
8.0	30.1	29.1	31.0	32.1	42.3	20.1	0.4300	0.4306	0.4294	0.0006	232.2	222.0	242.4
10.0	31.3	30.5	32.1	35.5	42.0	27.7	0.4292	0.4300	0.4285	0.0007	241.2	232.1	250.3
12.59	32.5	31.6	33.4	43.5	47.3	37.9	0.4278	0.4288	0.4268	0.0010	252.1	240.9	263.4
15.85	33.8	32.8	34.7	52.1	53.3	49.0	0.4263	0.4275	0.4253	0.0011	254.9	243.4	266.4
20.0	34.9	33.9	36.0	45.9	47.5	42.0	0.4242	0.4252	0.4233	0.0009	254.0	242.1	265.8
25.0	36.6	35.2	37.9	42.8	44.4	37.4	0.4214	0.4223	0.4206	0.0009	251.2	236.1	266.3
30.0	25.4	22.3	28.6	48.4	65.3	11.5	0.4366	0.4434	0.4314	0.0068	270.1	223.0	317.1
40.0	26.0	23.4	28.6	43.3	60.7	12.3	0.4360	0.4413	0.4318	0.0053	255.2	219.2	291.2
Y = 0.55
0.1	11.5	11.3	11.6	1.4	1.8	1.0	0.4620	0.4623	0.4618	0.0003	199.1	195.2	203.0
0.5	13.0	13.4	12.5	18.8	25.5	11.7	0.4569	0.4543	0.4597	−0.0025	402.7	434.3	371.0
1.0	13.8	12.2	15.4	28.0	15.9	34.9	0.4628	0.4739	0.4543	0.0111	555.1	442.7	667.6
2.0	14.4	12.3	16.6	31.3	11.8	41.5	0.4661	0.4833	0.4541	0.0172	590.9	434.1	747.7
X-ray Bursts
Horowitza	29.3	27.7	30.9	38.9	51.4	21.0	⋯	⋯	⋯	⋯	236.0	216.9	255.1
XRB	34.5	28.6	40.4	127.5	146.5	39.1	0.4251	0.4324	0.4201	0.0073	325.2	246.7	403.6

Note.

^aThe nuclear abundances in terms of atomic mass number (A) were not given for this mixture and so corresponding Y_e values could not be calculated.

Download table as:  ASCIITypeset image

For the same initial mixture used in the molecular dynamics simulations of Horowitz et al. (2007), we find good agreement with their results (as did Medin & Cumming 2010). Horowitz et al. (2007) found the average nuclear charge, $\langle Z\rangle$, of the initial and the liquid and solid equilibrium phases to be $\langle Z{\rangle }_{i}=29.30$, $\langle Z{\rangle }_{l}=28.04$, and $\langle Z{\rangle }_{s}=30.48$ with equilibrium impurity parameters of Q_i = 38.9, Q_l = 52.7, Q_s = 22.3, and Γ values of ${{\rm{\Gamma }}}_{i}=247$, ${{\rm{\Gamma }}}_{l}=233$, ${{\rm{\Gamma }}}_{s}=261$. Our semi-analytic method yields $\langle Z\rangle$ values of $\langle Z{\rangle }_{i}=29.3$, $\langle Z{\rangle }_{l}=27.47$, and $\langle Z{\rangle }_{s}=30.9$ with impurity parameters of Q_i = 38.9, Q_l = 51.4, Q_s = 21.0 and Γ values of ${{\rm{\Gamma }}}_{i}=236$, ${{\rm{\Gamma }}}_{l}=217$, ${{\rm{\Gamma }}}_{s}=255$ (listed as "Horowitz" in Table 1).

The initial number fractions x_i and relative number fractions of solid to liquid ${\text{}}{x}_{s}/{\text{}}{x}_{l}$ for five different models are shown in Figure 2 for helium abundance Y = 0.2752. Lighter nuclei, relative to the average Z of the composition, are preferentially retained in the liquid (${x}_{s}/{x}_{l}\lt 1$), whereas heavier nuclei are found in the solid (${x}_{s}/{x}_{l}\gt 1$). This extends the result of Horowitz et al. (2007) to a wider range of compositions. The average atomic number ($\langle Z{\rangle }_{i}$) of the initial mixture (vertical dashed lines) provides a quick but imperfect means of determining the "cross-over point" below which nuclei are preferentially retained in the liquid and above which tend to crystallize in the solid. This method tends to underestimate the point for lighter compositions ($\langle Z{\rangle }_{i}\;\lt 25$) and overestimate it for heavier compositions ($\langle Z{\rangle }_{i}\;\gt 25$).

Zoom In Zoom Out Reset image size

The top panels of Figure 2, corresponding to an accretion rate $\dot{m}$ = 0.1 ${\dot{m}}_{{\rm{Edd}}}$, show that for light compositions the general rule of light nuclei being retained in the liquid and heavy nuclei crystallizing in the solid no longer applies. This observation for the light compositions is further verified by equivalent plots for Y = 0.55 data shown in Figure 3. In this light composition regime ($\langle Z{\rangle }_{i}\lt 15$) the correlation of phase with nuclear charge is no longer as strong as for heavier mixtures. Instead, the lightest compositions corresponding to the top panel of Figure 2 and the top two panels of Figure 3 show that the solid phase is preferentially populated by elements with $Z\sim \langle Z{\rangle }_{i}$, which includes the most abundant element of the composition (Z = 12). These nuclei are not necessarily all heavy relative to $\langle Z{\rangle }_{i}$, allowing for the possibility of phase separation producing heavier equilibrium liquid than the corresponding solid as is the case for the $\dot{m}=0.5$ ${\dot{m}}_{{\rm{Edd}}}$, Y = 0.55 composition. As the compositions become heavier (i.e., the bottom two panels of Figure 3) the relative number fraction plots (${\text{}}{x}_{s}/{\text{}}{x}_{l}$) begin their convergence to the general trend seen in the heavy composition regime.

Zoom In Zoom Out Reset image size

If we drop the assumption that He and C burn to Mg before freezing begins, we see a substantial change in the results. The presence of carbon in these light compositions appears to exert a strong influence on the phase separation trend of the ${\text{}}{x}_{s}/{\text{}}{x}_{l}$ plots. In the light composition regime the rp-process ashes produce substantial amounts of carbon (${x}_{i}^{C}=0.42-0.78$) which when included with any residual He gives less regular behavior in the trend of ${\text{}}{x}_{s}/{\text{}}{x}_{l}$ and tends to magnify ${\text{}}{x}_{s}/{\text{}}{x}_{l}$ for certain elements. As seen in the bottom panel of Figure 4, for these light mixtures carbon's presence significantly lightens the composition which in turn drastically increases the relative number fractions for specific elements. This effect tends to decrease as the compositions become heavier and the contributions of C and He are no longer as significant. An exception to this was found for the relatively heavy X-ray burst composition (the green data of Figures 1 and 5) with its broad range of nuclear species spanning Z = 2–44. The top panel of Figure 4 shows that including C and He in the composition with initial number fractions (${x}_{i}^{C},{x}_{i}^{\mathrm{He}}=0.044,0.096$) causes irregularities to develop in the ${\text{}}{x}_{s}/{\text{}}{x}_{l}$ trend even though the corresponding $\langle Z{\rangle }_{i}$ of the composition does not change substantially and is well outside the domain where we expect these irregularities to develop (i.e., outside the light mixture regime). Although these results are interesting it is unlikely that these light compositions with high carbon content reflect the actual composition undergoing phase separation at the base of a neutron star's ocean. These models do indicate, however, that even slight variations in the initial composition can elicit substantial changes to the relative number fractions of certain elements.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 5 shows the impurity parameter for the initial, solid, and liquid compositions for all mixtures (red: Y = 0.55; and black: Y = 0.2752). Also plotted are the impurity parameters for the 17-species composition of Horowitz et al. (2007) (blue) and the additional X-ray burst composition (green). Across all data sets the impurity parameter of the solid phase is always lower than its initial or liquid counterpart. Comparing Figure 5 with Figure 1, we see that the compositions with the largest difference in impurity parameter between liquid and solid, ${\text{}}{Q}_{l}\quad -\quad {\text{}}{Q}_{s}$, correspond to the compositions with the largest difference in $\langle Z\rangle$ between liquid and solid. In general, a larger enrichment of light elements in the liquid leads to a purer solid phase.

We find that the degree of chemical separation depends on the fractional spread in Z in the mixture, which we measure with the parameter $\bar{\sigma }={Q}_{{\rm{imp}}}^{1/2}/\langle Z{\rangle }_{i}$. The bottom panel of Figure 6 shows how this parameter depends on the $\langle {Z}_{i}\rangle$ for our mixtures. In general, heavier mixtures tend to have a smaller fractional spread in Z and vice versa, with the exception of the lightest mixtures. The top three panels of Figure 6 show that the Coulomb coupling parameter of the initial mixture at equilibrium (${{\rm{\Gamma }}}_{i}$), the contrast in Y_e between liquid and initial mixture (${Y}_{e,l}-{Y}_{e-,i}$), and the heavy element enrichment in the solid ($\langle Z{\rangle }_{s}-\langle Z{\rangle }_{i}$) all increase with $\bar{\sigma }$. This means that mixtures with a larger $\bar{\sigma }$ have a lower melting temperature and a greater degree of chemical separation. The lightest mixtures deviate from these trends because they are in the light mixture regime where the solid forms from the most abundant element. In that case, there is no longer a clean separation in Z between elements that go into the liquid and elements that go into the solid, and the degree of chemical separation is much smaller. For example, the models with $\langle Z{\rangle }_{i}\leqslant 13$ in Table 1 have much smaller values of ${Y}_{e,l}-{Y}_{e,i}$. These lightest mixtures do still lie on the same trend of ${{\rm{\Gamma }}}_{i}$ against $\bar{\sigma }$ however, so that $\bar{\sigma }$ is a good predictor of ${{\rm{\Gamma }}}_{i}$ no matter what the value of $\langle Z{\rangle }_{i}$.

Zoom In Zoom Out Reset image size

Medin & Cumming (2010) compared their results for the 17-species mixture with a two-component model in which the chemical separation for any given element was calculated by approximating the mixture as consisting of that element plus the most abundant element (which was Se for the 17-species mixture). The pattern of enrichment of heavy elements in the solid and light elements in the liquid was well-reproduced by this simpler model, although with differences in the detailed values of ${\text{}}{x}_{s}/{\text{}}{x}_{l}$.

Here, we investigate how well the two-component model reproduces our results for a wider range of compositions, and extend it to include three components. Figure 7 compares the results from the 17-species calculation to the two- and three-component approximation for three steady-state models summarized in Table 2. For the two-component model, we follow Medin & Cumming (2010) and calculate the equilibrium solid–liquid number fractions for each element assuming the plasma is composed of only two ion species, the element itself and the most abundant element in the mixture. The initial composition of the mixture is chosen such that the ratio of the abundances of the two elements is the same as the 17-component system. In the three-component approximation, we calculate ${\text{}}{x}_{s}/{\text{}}{x}_{l}$ of each element assuming the plasma is composed of the two most abundant elements in addition to the element in question. The number fractions are renormalized for each element so that its number fraction relative to the two most abundant elements agrees with the 17-component mixture.

Zoom In Zoom Out Reset image size

Table 2.  Initial Number Fractions (x_i) for the Two Most Abundant Species

$\dot{m}/{\dot{m}}_{{\rm{Edd}}}$	Y	Z1	x1	Z2	x2
2	0.55	12	0.49	8	0.12
0.5	0.2752	28	0.41	12	0.27
10	0.2752	28	0.31	30	0.30

Download table as:  ASCIITypeset image

For most of the compositions the two-component approximation tends to reproduce the general number fraction trend of the 17-component plasma with discrepancies further reduced in the three-component model, as illustrated by the two examples in the upper panels of Figure 7. However, we noticed that compositions with a large $\bar{\sigma }$ show significant differences. In particular, the two- and three-component models show the opposite trend of ${\text{}}{x}_{s}/{\text{}}{x}_{l}$ with Z compared to the 17-component calculation. This can be seen in the lower panel of Figure 7. The difference arises because compositions with large $\bar{\sigma }$ mark the transition between the light and heavy element regimes. In this case, whereas the full composition is close to being in the heavy element regime, the two- and three-component approximations remain firmly in the light elements regime and therefore show the opposite trend.

Our results have implications for the depth of the neutron star ocean, the composition of the outer parts of the neutron star crust, and for driving convection in the neutron star ocean. One caveat is that we have calculated the liquid–solid equilibrium for a 50/50 mixture of liquid and solid, and so have not followed the complete freezing of the mixture. The depth of the ocean floor likely does not correspond to the value of ${\rm{\Gamma }}$ for a 50/50 mixture, nor is the composition at the top of the crust the same as the 50/50 equilibrium solid composition. In fact, in a steady state, the average composition entering the crust must be the same as the rp-process ashes entering the ocean, and therefore the ${Q}_{{\rm{imp}}}$ in the outer crust would be the same on average as the rp-process ashes.

It is unlikely, however, that a single-phase solid with the same composition as the rp-process ashes can form. For a general multicomponent mixture, multiple solid phases are likely to form instead, with the liquid adopting the composition at the eutectic point. Horowitz et al. (2007) found this for the particular mixture they studied, and a simpler example is shown in Figure 1 of Medin & Cumming (2011) for a two-component mixture of O and Se. The outer crust is therefore likely to consist of multiple solid phases. Diffusion could result in the different solid phases merging in the outer crust. Hughto et al. (2011) find a diffusion coefficient of $D\sim {10}^{-6}{\omega }_{p}{a}^{2}\sim {10}^{-8}\ {\rho }_{9}^{-1/6}\mathrm{at}\ {\rm{\Gamma }}$ ∼180, just after the solid forms. Since the ocean may take ∼100 years to accrete, the diffusion length is $\sim 3\ {\rm{cm}}$ close to the top of the crust. This may be enough to merge different solid phases, although the diffusion coefficient decreases sharply with density (Hughto et al. 2011).

We are not able to confidently calculate beyond the 50/50 equilibrium mixture and the possible multiphase solid composition at present, without further comparisons with molecular dynamics simulations. As discussed by Medin & Cumming (2010), as the solid fraction increases the resulting compositions depend more sensitively on the particular form of ${\rm{\Delta }}{f}_{s}$ chosen (Equation (3)) because ${\rm{\Delta }}{f}_{s}$ dominates the free energy at large ${\rm{\Gamma }}$. It is also not clear if a steady state is reached in the outer crust, as the rp-process ashes are sensitive to variations in accretion rate onto the star. Given these uncertainties, we will assume in exploring the implications of our results that the results we have obtained for 50/50 mixtures give a good estimate of the typical mixtures that make up the outer crust. Further work is needed to constrain the nucleation rate of the equilibrium solid at the ocean–crust boundary. The actual equilibrium mass fraction could therefore differ substantially from the 50/50 mass fraction assumed here.

The equilibrium liquid–solid compositions were also determined for the lightest mixtures (i.e., low accretion rates $\dot{m}\leqslant 0.5{\dot{m}}_{{\rm{Edd}}}$) assuming no burning of carbon. We find large variability in the equilibrium ${{\rm{\Gamma }}}_{i}$ values for these mixtures (${{\rm{\Gamma }}}_{i}=265-626$) which roughly equates to an order of magnitude density difference in the approximate location of the ocean–crust boundary. We also find that for the three lightest compositions, the corresponding ${{\rm{\Gamma }}}_{i}$ values increase with $\langle Z{\rangle }_{i}$. This result runs contrary to one-component plasmas which crystallize more readily with increase in nuclear charge. This discrepancy is likely due to the increase in the contribution of the entropy of the mixing term in Equation (1) as compositions become more heterogeneous with an increase in accretion rate. Research on the variability of the ocean–crust boundary with composition could offer opportunities to further constrain properties of cooling neutron stars post-outburst or accretion regimes of superbursting sources.

The impurity parameter ${Q}_{{\rm{imp}}}$ sets the electron-impurity scattering rate in the crust. Whether this determines the electrical and thermal conductivities depends on how it compares to the electron–phonon scattering rate (Flowers & Itoh 1976). At a density $\rho$, impurity scattering will dominate if ${Q}_{{\rm{imp}}}\gt {Q}_{{\rm{crit}}}$, where

$$\begin{eqnarray}&&{Q}_{{\rm{crit}}}\quad =\quad 32\ \displaystyle \frac{{T}_{8}^{2}}{{\rho }_{11}^{5/6}}\left(\displaystyle \frac{Z}{20}\right){\left(\displaystyle \frac{{Y}_{e}}{0.4}\right)}^{-4/3}\end{eqnarray} \tag{ 4 }$$

(Cumming et al. 2004), and ${T}_{8}=T/{10}^{8}\ {\rm{K}}$. As density increases, impurity scattering is more likely to dominate, as the increasing Debye temperature reduces the phonon-scattering contribution.

Constraints on the crust thermal conductivity, and therefore ${Q}_{{\rm{imp}}}$, have come from modeling the observed cooling of neutron star transients in quiescence (Shternin et al. 2007; Brown & Cumming 2009; Page & Reddy 2012; Turlione et al. 2015). Brown & Cumming (2009) assumed a constant ${Q}_{{\rm{imp}}}$ throughout the entire crust and showed that fits to the transient sources MXB 1659-29 and KS 1731-260 required ${Q}_{{\rm{imp}}}\;\leqslant 10$, even when uncertainties in the neutron star surface gravity, and therefore crust thickness, were taken into account. This upper limit is indicated by the horizontal dashed line in Figure 5.

Brown & Cumming (2009) noted that their constraint on ${Q}_{{\rm{imp}}}$ applies to the inner crust, since electron–phonon scattering sets the conductivity in the outer crust. The value of ${Q}_{{\rm{imp}}}$ in the inner crust likely does not reflect the ${Q}_{{\rm{imp}}}$ of the freshly made solid at the top of the crust. Steiner (2012) modeled the evolution of multicomponent mixtures in the crust as they are compressed and undergo electron captures, neutron emissions, and captures near neutron drip. Those models showed a reduction in ${Q}_{{\rm{imp}}}$ near neutron drip, where neutrons are able to redistribute among nuclei (see also Gupta et al. 2008). The simplification of the mixture near neutron drip means that a relatively large ${Q}_{{\rm{imp}}}$ at the top of the crust could be consistent with the constraint ${Q}_{{\rm{imp}}}\lt 10$ from Brown & Cumming (2009). Our calculations of the mixture entering the outer crust can be used as inputs to calculations of the nuclear evolution through neutron drip.

Page & Reddy (2013) were able to derive constraints on the outer crust ${Q}_{{\rm{imp}}}$. They fit the neutron star transient XTE J1701-462 using different values of ${Q}_{{\rm{imp}}}$ for the outer and inner crust and found best fitting models with ${Q}_{{\rm{imp}}}\approx 15-30$ at densities $\rho \lt {10}^{12}\ {\rm{g}}\ {\mathrm{cm}}^{-3}$, transitioning to ${Q}_{{\rm{imp}}}\approx 3-4\ \mathrm{for}\quad \rho \gt {10}^{13}\;{\rm{g}}\;{\mathrm{cm}}^{-3}$. Their constraint for the outer crust is represented by a blue shaded band in Figure 5. Our results in this paper show that chemical separation on freezing could explain the inferred values of ${Q}_{{\rm{imp}}}$ for the outer crust. Figure 5 shows that the initial values of ${Q}_{{\rm{imp}}}$ are larger than 30 for compositions with $\langle Z\rangle \;\gtrsim$ 16, but that most of the corresponding solid compositions have ${Q}_{{\rm{imp}}}\lt 30$. Further work on the constraints on ${Q}_{{\rm{imp}}}$ in the outer crust from other sources would be interesting.

We can use our results to estimate the strength of compositionally driven convection in the neutron star ocean (Medin & Cumming 2011). The rate of buoyancy production at the ocean–crust boundary depends mostly on the difference in ${Y}_{e}$ between the rp-process ash mixture and the equilibrium liquid composition. Electrons dominate the pressure in the ocean (by a factor of $\sim {E}_{{\rm{F}}}/{k}_{{\rm{B}}}T$ ∼100), so that the compositional buoyancy is dominated by the gradient of ${Y}_{e}$. In that case, we can simplify the expression for the Brunt–Väisälä frequency (Bildsten & Cumming 1998; Medin & Cumming 2015) as

$$\begin{eqnarray}&&\displaystyle \frac{{N}^{2}H}{g}\approx \displaystyle \frac{{\chi }_{T}}{{\chi }_{\rho }}\left({{\rm{\nabla }}}_{{ad}}-{\left.\displaystyle \frac{d\mathrm{ln}T}{d\mathrm{ln}P}\right|}_{\star }\right)-\displaystyle \frac{{\chi }_{{Y}_{e}}}{{\chi }_{\rho }}{\left.\displaystyle \frac{d\mathrm{ln}{Y}_{e}}{d\mathrm{ln}P}\right|}_{\star },\end{eqnarray} \tag{ 5 }$$

where we have dropped the ion terms. In this limit, the expression for the convective heat flux in the ocean from Medin & Cumming (2011) (see their Equation (43)) becomes

$$\begin{eqnarray}&&{F}_{{\rm{conv}}}\quad =\quad \displaystyle \frac{{c}_{P}T\dot{m}}{{\chi }_{T}}{\chi }_{{Y}_{e}}\left(\displaystyle \frac{{Y}_{e,l}-{Y}_{e,i}}{{Y}_{e,l}}\right).\end{eqnarray} \tag{ 6 }$$

This equation gives the expected heat flux for steady accretion; in a time-dependent situation such as cooling after an accretion outburst, the $\dot{{m}}$ factor should be replaced by $\partial {y}_{m}/\partial$ ${\rm{t}}$, the rate of change of the column depth of the ocean floor due to temperature changes (Medin & Cumming 2015). Since the pressure is dominated by relativistic degenerate electrons, ${\chi }_{{Y}_{e}}\approx 4/3$, so that the composition enters Equation (6) through the factor $({Y}_{e}-{Y}_{e,i})/{Y}_{e}$.

The values of ${Y}_{e}$ for the initial, liquid, and solid mixtures are given in Table 1. The largest values of ${Y}_{e,l}-{Y}_{e,i}\mathrm{are}$ ≈0.02. This is in the range considered by Medin & Cumming (2011), who looked at two-component mixtures only for simplicity. In that paper, the steady state was defined as the point where the composition crystallizing at the bottom of the ocean matched the composition accreting from the top. This agreement suggests that their results represent a realistic range of outcomes for compositionally driven convection despite their two-component approximation. It is likely that the values ${Y}_{e,l}-{Y}_{e,i}$ tabulated in this paper provide lower bounds on the actual convective heat flux of that mixture. This is because unlike Medin & Cumming (2011), the results calculated here are for a system where 50% of the mixture crystallizes into solid. Generally, this equilibrium solid state is substantially different (heavier) from the initial mixture accreting at the top of the ocean. Therefore we expect further phase separation and enrichment of the ocean in light elements to arrive closer to the steady state of Medin & Cumming (2011). Their O–Se mixture, which was chosen to have a similar $\langle Z\rangle$ to the mixture studied by Horowitz et al. (2007), had 2% oxygen by mass at the top of the ocean, so that ${Y}_{e,i}\approx {Y}_{e,{\rm{Se}}}=34/79=0.430$. The steady-state composition at the base of the ocean was 37% oxygen, giving ${Y}_{e}-{Y}_{e,i}=0.026$. Their Fe–Se mixture had ${X}_{\mathrm{Fe}}=0.23$ at the top of the ocean and ${X}_{\mathrm{Fe}}=0.37$ at the base, giving ${Y}_{e}-{Y}_{e,i}=0.005$. Interestingly, some compositions that we study in this paper have ${Y}_{e}-{Y}_{e,i}\lt 0$, so that chemical separation would not lead to convection in those cases. Also we note that the lightest mixtures with $\langle Z\rangle \ \leqslant$ 13 have smaller values of ${Y}_{e,l}-{Y}_{e,i}$ by about a factor of ten compared to the heavier mixtures. Therefore, light element oceans should have much less convective driving.

Our results rely on an extension of the analytic fits made by Ogata et al. (1993) to Monte Carlo simulations of three-component plasmas to multicomponent plasmas. Even for the three-component case, Hughto et al. (2012) found systematically lower melting temperatures (higher ${\rm{\Gamma }}$) in their MD simulations compared to the semi-analytic model. They noted that this discrepancy seemed to grow with impurity parameter ${Q}_{{\rm{imp}}}$, perhaps suggesting a problem with the form of the deviation from linear mixing ${\rm{\Delta }}{f}_{s}$ (Equation (3)) assumed in the semi-analytic model. Figure 6 shows that ${{\rm{\Gamma }}}_{i}$ increases with ${Q}_{{\rm{imp}}}$ at fixed $\langle Z\rangle$, which may explain the trend with ${Q}_{{\rm{imp}}}$ seen by Hughto et al. (2011), since ${\rm{\Delta }}{f}_{s}$ becomes more and more important at larger values of ${\rm{\Gamma }}$. It is important to carry out further comparisons with molecular dynamics simulations (Horowitz et al. 2007; Hughto et al. 2012) to check and refine these assumptions about the functional form of the free energy as well as to investigate the formation of multiple solid phases. We hope that our results will give a useful starting point for such comparisons.