A SIGNATURE OF CHEMICAL SEPARATION IN THE COOLING LIGHT CURVES OF TRANSIENTLY ACCRETING NEUTRON STARS

The observation of surface cooling of accreting neutron stars on timescales of days to years after they go into quiescence has opened up a new window on the physics of neutron star crusts. Six neutron stars in low mass X-ray binaries have now been observed to cool after extended accretion outbursts that were long enough to heat the crust significantly out of thermal equilibrium with the core (Cackett et al. 2010, 2013; Fridriksson et al. 2011; Degenaar et al. 2011, 2013a, 2013b). The subsequent thermal relaxation of the crust depends on its physical properties such as its thickness, thermal conductivity, and heat capacity, with deeper regions being probed at successively later times in the cooling light curve (Eichler & Cheng 1989; Brown & Cumming 2009, hereafter BC09). Shternin et al. (2007) and BC09 showed that the observed cooling curves of KS 1731–260 and MXB 1659–29 imply that the inner crust has a thermal conductivity corresponding to an impurity parameter of order unity. Page & Reddy (2013) find that the cooling curve of XTE J1701–462 is compatible with similar crust microphysics to KS 1731–260 and MXB 1659–29.

Chemical separation occurs when, as a material freezes, the equilibrium compositions of the liquid and solid phases are different. This important process can drive sedimentation or mixing in white dwarf (Althaus et al. 2012) and giant planet (e.g., Wilson & Militzer 2010) interiors and Earth's core (Hirose et al. 2013). Horowitz et al. (2007) carried out molecular dynamics simulations of the freezing of the multicomponent plasma expected in the outer layers of accreting neutron stars, and found chemical separation occurred, with lighter nuclear species being preferentially retained in the liquid phase. In a previous paper (Medin & Cumming 2011, hereafter MC11) we considered the effect of chemical separation while accretion is ongoing. In that case, as matter is driven to higher pressure and crosses the freezing depth, the light elements released into the neutron star ocean will drive convection, mixing them throughout the ocean and raising the light element fraction. At the same time, because the convection is driven by composition gradients in an otherwise thermally stable layer, the convective heat flux is inward, potentially heating the layers deep in the ocean.

In this Letter, we consider the cooling of the outer layers of the star in quiescence, during which the liquid layers freeze into a solid and chemical separation occurs, enriching the ocean further in light elements. We show that the resulting compositionally driven convection can significantly modify the expected cooling curves of accreting transients. We first discuss the expected size of the compositionally driven heat flux and show that it is easily comparable to the cooling flux in the ocean (Section 2). We then present simulations of ocean cooling including heat transport by convection, showing that the light curve is significantly modified (Section 3). We conclude with a discussion of the theoretical uncertainties and implications for observed sources (Section 4).

The outermost layers of an accreting neutron star form a gaseous atmosphere and liquid ocean (Bildsten & Cutler 1995), in regions where the Coulomb coupling parameter Γ ≡ 〈Z^5/3〉e²/ak_BT < 175 (Potekhin & Chabrier 2000), where 〈Z^5/3〉 is the number average over the mixture of nuclei with charge Z_i, a is the electron sphere radius, and T is the temperature. A useful measure of the depth is the column depth y, related to the pressure by P = gy where g is the gravity. The base of the ocean lies at

$$\begin{equation} y_{b,14} = 1.8 \left(\frac{T_{b,8}}{3}\right)^4 \left(\frac{\left\langle Z_b^{5/3} \right\rangle }{100}\right)^{-4}\left(\frac{g_{14}}{2}\right)^{-1}, \end{equation} \tag{ 1 }$$

where y₁₄ = y/10¹⁴ g cm⁻², T₈ = T/10⁸ K, g₁₄ = g/10¹⁴ cm s⁻², and the subscript "b" signifies that the quantities are taken at the base of the ocean. Note that in Equation (1) we assume $\langle Z_b^{5/3}\rangle \sim 100$, appropriate for an oxygen-enriched ocean; a heavy-element-rich one can be much shallower, y_b ∼ 10¹² g cm⁻².

BC09 showed that when accretion ends, the cooling of the star proceeds by the temperature profile relaxing from the outside-in. At any given time, lower density layers have cooled and adopted a constant flux temperature profile, whereas deeper layers have yet to thermally relax. The transition occurs at a depth given by setting the thermal timescale

$$\begin{equation} \tau \simeq \frac{\rho c_P H_P^2}{2K}=19\, {\rm days}\, y_{14}\left(\frac{T_8}{3}\right)^{-1}\left(\frac{Y_e}{0.4}\right)^{2}\left(\frac{g_{14}}{2}\right)^{-1} \end{equation} \tag{ 2 }$$

(Henyey & L'Ecuyer 1969; Equation (7) of BC09) equal to the current time. In Equation (2), ρ is the mass density, c_P is the specific heat, K is the thermal conductivity, and H_P = −dr/dln y = y/ρ is the pressure scale height. In the ocean, the pressure is dominated by the relativistic degenerate electrons with $E_F/k_BT=400 y_{14}^{1/4}(3/T_8)(g_{14}/2)^{1/4}$, giving the scalings K∝y^1/4T, H_P∝y^1/4, and c_P ≈ 3k_b/Am_p and T constant (e.g., MC11 and references therein).

Equations (1) and (2) show that tens of days after the onset of quiescence, the base of the ocean will start to cool and solidify. As new crust is formed, the bottom of the ocean becomes enriched in light elements because of chemical separation (Horowitz et al. 2007; Medin & Cumming 2010). In MC11, we estimated the timescales for particle nucleation, growth, and sedimentation of solid particles, finding that these "microscopic" timescales were much shorter than "macroscopic" timescales such as the time to accrete the ocean, or more relevant here, the time for the ocean to cool. The picture then is that fluid elements with a light composition are deposited at the base of the ocean and will rise upward, driving convective mixing. Moreover, the convection occurs in a medium that is otherwise thermally stratified. Therefore the convective heat flux, which is proportional to the excess of the temperature gradient ∇ = dln T/dln P in the star compared to the adiabatic gradient ∇_ad ≃ 0.4, or F_conv∝∇ − ∇_ad, is negative; the compositionally driven convection transports heat inward.

We can estimate the expected size of F_conv by comparing the heat and composition fluxes. In mixing length theory, the heat flux is F_conv ≈ ρv_convc_PT(∇ − ∇_ad), where v_conv is the convective velocity; whereas the composition flux is F_X ≈ ρv_convX∇_X, where ∇_X is the composition gradient in the star ∇_X = dln X/dln P (MC11), and X is the light element mass fraction. The convective velocity is given by the superadiabaticity,

$$\begin{equation} v_{\rm conv}^2\approx g l^2 [\chi _T(\nabla -\nabla _{\rm ad})+\chi _X\nabla X], \end{equation} \tag{ 3 }$$

where l is the mixing length, χ_T = ∂ln P/∂ln T|_{ρ, X}, and χ_X = ∂ln P/∂ln X|_{T, ρ}. In the ocean, the convection is extremely efficient and therefore close to marginal stability, χ_T(∇_ad − ∇) ≈ χ_X∇X. The heat and composition fluxes are then related by

$$\begin{equation} F_{\rm conv} = - c_PT \left(\frac{\chi _X}{\chi _T} \right) \frac{F_X}{X}. \end{equation} \tag{ 4 }$$

During cooling, the composition of the ocean is changing with time as the ocean is enriched in light elements. For efficient convection, the light element fraction is fairly constant over the ocean (∇_X ≈ χ_T∇_ad/χ_X∝k_BT/E_F ≪ 1), and so ∂X/∂t is roughly constant throughout the ocean. The composition flux is then F_X ≈ y(∂X/∂t), and the heat flux is

$$\begin{equation} F_{\rm conv} = -y c_PT \frac{\chi _X}{\chi _T}\frac{\partial \ln X}{\partial t}. \end{equation} \tag{ 5 }$$

Taking χ_X ≈ 0.1 and χ_T ≈ 10k_BT/〈Z〉E_F (MC11), we find a flux that increases steeply inward, F_conv∝y^5/4.

If cooling is unaffected by compositionally driven convection, we estimate the timescale on which X is changing, ∂t/∂ln X, to be a typical thermal time at the ocean floor: in the standard cooling model the bulk of the ocean freezes on a thermal time; as it does it releases nearly all of its light elements into the remaining ocean and approximately doubles the light element content there. Equation (5) becomes

$$\begin{equation} F_{\rm conv}\approx {-} 10^{25} {\rm \, erg\, cm^{-2} s^{-1}} y_{14}^{5/4} \left(\frac{\partial t/\partial \ln X}{10\, {\rm days}}\right)^{-1}, \end{equation} \tag{ 6 }$$

where we have assumed Y_e = 0.4 and g₁₄ = 2 (and y_{b, 14} = 1 and T₈ = 3 in Equation (2)). The flux from compositionally driven convection in Equation (6) easily outweighs the cooling flux, $\sigma T_{\rm eff}^4 = 10^{20} {\rm \,erg\, cm^{-2} s^{-1}} (T_{\rm eff}/100 {\rm \, eV})^4$. This suggests that the cooling of an ocean with chemical separation included should be significantly different than without, and motivates our numerical calculations that will be presented in the next section. As we shall see below, the freezing of the ocean is strongly regulated by the convective heat flux, and in turn keeps the heat flux at a much lower level than is suggested by Equation (6).

We solve for the thermal relaxation of the neutron star ocean and crust by solving the thermal diffusion equation following BC09, but including convective heat fluxes modeled using mixing length theory and assuming the convection is efficient as described in MC11, and in addition by following the composition profile as chemical separation and mixing occur. For the examples shown here, we take the ocean composition to be a mixture of oxygen and selenium, for which the phase diagram is shown in the lower panel of Figure 1 in MC11. A full description of our numerical code will be presented elsewhere (Z. Medin & A. Cumming 2014, in preparation); here we focus on the resulting light curves and describe the influence of compositionally driven convection on the evolution.

Figure 1 shows an example light curve with and without compositionally driven convection included. Figures 2 and 3 show the temperature and composition profiles, respectively, at different times as the ocean and crust cool. The temperature profile at the end of the accretion outburst is similar to that assumed by BC09, with an inward directed heat flux, and (for Figures 2 and 3) an initial base composition $X_b^{\rm init} = 0.37$, corresponding to the steady state (MC11). The kink in the temperature profile at 10⁸ g cm⁻² is due to the jump in composition between the ocean and overlying H/He layer (cf. Paper I). The neutron star mass and radius are 1.62 M_☉ and R = 11.2 km, giving a redshift factor of 1 + z_surf = 1.32.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We find that the evolution proceeds in four stages. During stage 1, the base of the ocean has not yet started to cool and so the evolution is the same with or without convection included. The light curve is a power-law with slope given by Equation (8) of BC09, with the modification that ∂ln τ/∂ln y = 1 in the ocean rather than 3/4 as in the outer crust.

In stage 2 (20 days ≲ t_∞ ≲ 80 days in Figures 1–3), the cooling wave has reached the bottom of the ocean, and new crust begins to form, driving convection. Because convection transports heat inward from the ocean to the ocean–crust boundary, cooling at the boundary is delayed at the expense of more rapid cooling in the ocean. As a result, the temperature profile in the ocean becomes very steep and the light curve drops faster than without convection. During this stage, the convective heat flux is strong and the boundary cools very slowly; the transition depth and temperature y_b and T_b remain close to their values at the onset of quiescence. In this way, the convection acts analogously to the latent heat.³ The cusp in the temperature profile is due to the jump in conductive flux that must occur to balance the large inward convective flux in the ocean.

The temperature profile cannot steepen indefinitely, because eventually the temperature gradient at the base of the ocean ∇_b approaches ∇_L ≃ 0.25, the liquidus temperature gradient (i.e., how the melting temperature varies with pressure; see MC11). When ∇_b = ∇_L, multiple depths at the bottom of the ocean freeze simultaneously. This rapid freezing quickly suppresses itself, however, as strong compositionally driven convection heats the base of the ocean and melts the top of the crust, mixing a heavy-element fluid into the bottom of the ocean and thereby stabilizing the ocean against further convection. Cooling resumes, but convection remains off, such that the steep temperature profile at the ocean base that had been supported by convection quickly flattens due to heat conduction. Once the heavy-element-enriched fluid at the base of the ocean completely solidifies, compositionally driven convection resumes and the temperature profile steepens again.

In this manner ∇_b oscillates around ∇_L. We refer to this phase as stage 3 (80 days ≲ t_∞ ≲ 2100 days in Figures 1–3). In this stage convection is sporadic in the ocean and can no longer prevent the ocean base from cooling. As a result, the ocean boundary moves outward and the cooling wave then continues its inward motion through the crust. We find generally that ∂ln T_b/∂t ≫ ∂ln y_b/∂t; i.e., that the ocean–crust boundary cools rapidly but moves outward slowly. This is due to the compensating effect of light element enrichment on the freezing depth. During stage 3, the surface temperature T_eff still drops with time, but at a rate similar to without convection; and now ∇_b = ∇_L is constant while T_b and y_b drop and X_b increases.

During stage 4 (t_∞ ≳ 2100 days in Figures 1–3), the crust is thermally relaxed, the ocean cools too slowly for convection to support the steep gradient ∇_b = ∇_L, and the light curve returns to the shape it would have if there was no ocean convection. Note, however, that due to light element enrichment the asymptotic value of T_eff is slightly higher in the case with convection than without (the ocean thermal conductivity K∝〈Z〉⁻¹, and so for a given base temperature the outward flux is greater for a lower 〈Z〉). Alternatively, the ocean experiences an abrupt transition to stage 4 when the base is saturated with light elements (X = 1) and chemical separation halts; this happens at t_∞ ≃ 2000 days in the $X_b^{\rm init}=0.7$ case of Figure 1.

Previous calculations of the thermal relaxation of accreting neutron stars in quiescence have assumed that the ocean cools passively. In contrast, we have shown here that mixing in the ocean driven by chemical separation at the base leads to a significantly different evolution, changing the expected cooling curve. The early time (1–100 days) cooling curve of quiescent neutron stars potentially offers a remarkable new probe of the freezing and chemical separation of a classical plasma in "real time." The timing of the rapid drop in flux is sensitive to the composition and temperature of the ocean at the end of the outburst. The magnitude of the effect depends on the composition of the ocean, in particular the fraction of light elements and contrast in the atomic mass of the light and heavy nuclei.

Detecting the signature of convection will require better sampling of the early phase of the cooling curve. We can generally fit currently available light curves equally well with and without convection. These fits will be presented in a companion paper (Z. Medin & A. Cumming 2014, in preparation); two examples are shown in Figure 4. XTE J1701–462 (Fridriksson et al. 2011) is interesting because its high temperature means that the ocean remains liquid for hundreds of days. Page & Reddy (2013) fit the light curve with a standard cooling model (neglecting the two data points marked in Figure 4 that are argued to be contaminated with residual accretion). Including convection we find that we can match the drop in luminosity at 100–200 days. A general result is that convection lessens, but does not remove, the need for a shallow heat source in the ocean during accretion (BC09; Paper I; Page & Reddy 2013), because light element enrichment increases the thermal conductivity and reduces the temperature gradient, making the ocean hotter.

Zoom In Zoom Out Reset image size

Convection could help with recent observations of two classical transients (with shorter 2–3 month outbursts). In IGR J17480–2446 (Degenaar et al. 2013a), the flux remains elevated above the value observed before the outburst. Convection allows $X_b^{\rm init}$ to change from one accretion episode to the next, modifying the ocean conductivity and thereby changing the late time temperature (compare the solid curve with the hatched region in Figure 4). XTE J1709–267 showed a rapid decrease in temperature during a single 8 hr XMM observation (Degenaar et al. 2013b). If due to crust cooling, Degenaar et al. (2013b) suggested that a strong heat source must be operating at low densities within the ocean during the outburst, necessary for significant thermal relaxation after a short outburst of only 2–3 months. We can reproduce the rapid drop in temperature with convection if stage 2 of the light curve (Section 3) occurs during the observation.

There remains much to be explored theoretically. We have included only two species in our calculations, oxygen and selenium, which approximates the rp-process ashes used by Horowitz et al. (2007). The phase diagram for multicomponent mixtures is complex but can be calculated (Horowitz et al. 2007; Medin & Cumming 2010) and should be included. We have assumed that solid particles form at a single depth. However, electron capture reactions may occur in the ocean (for example, ⁵⁶Fe captures at a density of 1.5 × 10⁹ g cm⁻³; Haensel & Zdunik 1990), lowering the 〈Z〉 at that depth, and potentially leading to formation of solid particles pre-electron capture above the post-electron capture liquid layers. It will be important to include carbon burning in the models. Enrichment of the ocean with carbon remains a major issue for superburst models (Schatz et al. 2003). Chemical separation during the cooling phase will significantly enrich the ocean in light elements immediately following an outburst, much more efficiently than gravitational sedimentation during quiescence. This may have implications for the puzzling superburst observed immediately before the onset of an accretion outburst in EXO 1745–248 (Altamirano et al. 2012).

We thank Chuck Horowitz, Nathalie Degenaar, and Chris Fontes for useful discussions. A.C. is supported by an NSERC Discovery Grant and is an associate member of the CIFAR Cosmology and Gravity program. We are grateful for the support of an International Team on Neutron Star Crusts by ISSI in Bern. Z.M. was supported by a LANL Director's Postdoctoral Fellowship. This research was carried out in part under the auspices of the National Nuclear Security Administration of the U.S. Department of Energy at Los Alamos National Laboratory and supported by Contract No. DE-AC52-06NA25396.