THE THERMAL STATE OF KS 1731−260 AFTER 14.5 YEARS IN QUIESCENCE

Transient neutron star low-mass X-ray binaries (LMXBs) provide an opportunity to study the structural properties of neutron stars. Transient LMXBs typically alternate between periods of outburst (i.e., active accretion), and quiescence when accretion rates drop significantly. During outburst, the system's X-ray emission is dominated by the accretion disk and/or boundary layer. Accretion onto the neutron star surface can heat the stellar crust out of thermal equilibrium with the core (e.g., Ushomirsky & Rutledge 2001; Rutledge et al. 2002). Once the system returns to quiescence, thermal X-ray emission below a few keV originates from the surface of the neutron star (Brown et al. 1998) and the thermal relaxation of the crust can be observed directly. Direct observation of crustal cooling allows for the extraction of neutron star properties, such as thermal conductivity, crustal structure, and core temperature via cooling models (Shternin et al. 2007; Brown & Cumming 2009; Page & Reddy 2013; Deibel et al. 2015; Turlione et al. 2015; Ootes et al. 2016).

Crustal cooling of neutron star LMXBs has been observed in eight sources: KS 1731−260 (Wijnands et al. 2001, 2002; Cackett et al. 2006, 2010), MXB 1659−29 (Wijnands et al. 2003, 2004; Cackett et al. 2006, 2008, 2013), EXO 0748−676 (Degenaar et al. 2009, 2011b, 2014; Díaz Trigo et al. 2011), XTE J1701-462 (Fridriksson et al. 2010, 2011), IGR J17480−2446 (Degenaar et al. 2011a, 2013, 2015), MAXI J0556−332 (Homan et al. 2014), Swift J174805.3−244637 (Degenaar et al. 2015), and potentially Aql X-1 (Waterhouse et al. 2016). The cooling curves of these sources generally show a significant drop in temperature immediately following their return to quiescence, typically showing a close to exponential decay, with the curves then flattening as the cooling continues (see Figure 5 of Homan et al. 2014 for a comparison of several sources). However, the sources show a wide distribution of temperatures at the beginning of quiescence (∼100–300 eV), cooling timescale (e-folding time), and outburst timescales (a few months to nearly 24 years). While it was initially thought that long (quasi-persistent) outbursts were required to heat the crust out of thermal equilibrium, the observation of crustal cooling in IGR J17480-2446, Swift J174805.3−244637, and Aql X-1 demonstrates that it can occur with outbursts of only a few months.

KS 1731−260 was discovered by the Mir-Kvant instrument in 1989 August, and in previous data it was also observed in outburst in 1988 October (Sunyaev et al. 1989, 1990). It remained in outburst until early 2001 (Wijnands et al. 2001). KS 1731−260 was the first source for which crustal cooling was observed and has provided the longest cooling baseline for any source to date. A detailed history of the source can be found in Cackett et al. (2006). Since returning to quiescence, monitoring observations with Chandra and XMM-Newton have shown continued cooling up until 2009 (8 years into quiescence; Wijnands et al. 2001, 2002; Cackett et al. 2006, 2010). Physical models that track the thermal evolution of the neutron star reproduce the observed cooling in KS 1731−260 (e.g., Shternin et al. 2007; Brown & Cumming 2009; Page & Reddy 2013; Turlione et al. 2015; Ootes et al. 2016). Fits to the crustal cooling of KS 1731−260 consistently suggest that the crust has a high thermal conductivity, indicating a low impurity parameter (i.e., an ordered lattice crustal structure; discussed later). Moreover, in order to match the early evolution in the cooling curves, it has been found that an additional source of heating at shallow depths in the crust is often needed (e.g., Brown & Cumming 2009; Degenaar et al. 2011a; Deibel et al. 2015; Waterhouse et al. 2016). Such shallow heating is also required in order to provide the conditions for superburst ignition (Cumming et al. 2006; Gasques et al. 2007).

Here we present a new Chandra observation of KS 1731−260, taken 6 years after the previous Chandra observation and a total of approximately 14.5 years into quiescence. In addition to new data, recent updates to the crust cooling model used in Brown & Cumming (2009) (see Section 3.2) allow for the exploration of new parameters such as the influence of a nuclear pasta layer and additional shallow heating of the crust. Our findings suggests that cooling has likely halted in KS 1731−260 and that the neutron star crust has returned to thermal equilibrium with the core. In Section 2 we give an overview of data reduction and spectral analysis. In Section 3 we discuss empirical and physical models to fit the quiescent light curve. In Section 4 we discuss our findings and future work.

The new Chandra observation of KS 1731−260 was performed over 150 ks in three separate pointings. A 66 ks segment was performed on 2015 August 6/7 (ObsID: 16734). A 20 ks segment was performed on 2015 August 8 (ObsID: 17706), and a 64 ks segment was performed on 2015 August 9 (ObsID: 17707). As with the previous Chandra observations, KS 1731−260 was at the default aimpoint of the ACIS-S3 chip. Due to the time since the previous Chandra and XMM-Newton observations, we decided to reanalyze all data with the latest calibration files and software. Full details on the previous observations can be found in Wijnands et al. (2001, 2002) and Cackett et al. (2006, 2010); see also Table 1.

Table 1.  Neutron Star Atmosphere Fitting Parameters

Parameter	2428	013795201/301	3796	3797	0202680101	6279/5486	10037/10911	16734/17706/707
	(CXO)	(XMM)	(CXO)	(CXO)	(XMM)	(CXO)	(CXO)	(CXO)
MJD	51995.1	52165.7	52681.6	52859.5	53430.5	53512.9	54969.7	57242.1
${{{kT}}_{\mathrm{eff}}}^{\infty }$ (eV)	104.6 ± 1.3	89.5 ± 1.03	76.4 ± 1.8	73.8 ± 1.9	71.7 ± 1.4	70.3 ± 1.9	64.5 ± 1.8	64.4 ± 1.2
Fobs (10−15 erg cm−2 s−1)	41.6 ± ${}_{3.6}^{3.4}$	16.1 ± ${}_{0.97}^{1.10}$	5.98 ± ${}_{0.91}^{0.71}$	4.79 ± ${}_{0.63}^{0.58}$	3.97 ± ${}_{0.51}^{0.53}$	3.49 ± ${}_{0.41}^{0.61}$	1.99 ± ${}_{0.36}^{0.30}$	1.97 ± ${}_{0.22}^{0.23}$
Lbol (1033 erg s−1)	2.69	1.47	0.76	0.64	0.59	0.55	0.39	0.39

Note. Mass and radius are fixed to 1.4 M_⊙ and 10 km, respectively. Spectra were modeled using an absorbed neutron star atmosphere model (nsa) and a photoelectric absorption model (phabs). The distance to KS 1731−260 was set to 7 kpc, and the column density was tied between observations, giving a best-fitting value of N_H = 1.30 ± 0.06 × 10²² cm⁻². The effective temperature (${{kT}}_{\mathrm{eff}}^{\infty }$) is corrected for gravitational redshift (i.e., it is the effective temperature for an observer at infinity). The ObsID and observatory are indicated at the top of the table (CXO—Chandra and XMM—XMM-Newton). The observed flux was calculated over the 0.5–10 keV range. The bolometric luminosity was calculated over the 0.01–100 keV range. We do not include errors on L_bol, due to large systematic errors.

Download table as:  ASCIITypeset image

2.1. Chandra Data Reduction

2.2. XMM Data Reduction

2.3. Spectral Analysis

We fit all available Chandra and XMM-Newton spectra since the end of outburst using XSPEC (ver. 12.9.0; Arnaud 1996) following a procedure similar to Cackett et al. (2006, 2010). We modeled the spectra using a neutron star atmosphere model (nsa; Zavlin et al. 1996), modified by photoelectric absorption within our galaxy (phabs). The nsa model has been used previously in studies of KS 1731−260, and other models, such as nsatmos (Heinke et al. 2006), provide consistent results within 1σ. A comparison of atmosphere models can be found in Section 3.4 of Heinke et al. (2006). For our analysis and modeling, we fix the neutron star radius to canonical values of 10 km and the mass to 1.4 M_⊙. We set the distance to the source at 7 kpc (Muno et al. 2000), resulting in a normalization parameter = 2.041 × 10⁻⁸ pc⁻². The normalization is also a fixed parameter. We show the X-ray spectrum, and best-fitting model, for the new Chandra observation in Figure 1, compared to the 2009 Chandra observation. The latest spectrum is still thermal, with no need for an additional power-law component, and shows no significant change between 2009 and 2015 (see Figure 1).

Zoom In Zoom Out Reset image size

All spectra were fit simultaneously with absorption column density, N_H, tied between all observations and the effective temperature set as a free parameter. Due to the close proximity of some observations, several spectra had their parameters tied together. We tied the parameters of the spectra for XMM-Newton observations taken around MJD 52165.7 (ObsIDs 012795201/301) and Chandra observations taken around MJD 53512.9 (ObsIDs 6279/5468), MJD 54969.7 (ObsIDs 10037/10911), and MJD 57242.1 (ObsIDs 16734/17706/17707). Using GRPPHA, we binned the spectra to have one count per bin. We used C-statistic to fit the binned spectra. C-statistic was used rather than χ² statistics due to the low number of total counts in some spectra.

The spectral fitting results are shown in Table 1. The effective temperature is corrected for gravitational redshift (i.e., it is the effective temperature for an observer at infinity). The new Chandra observation gives ${{kT}}_{\mathrm{eff}}^{\infty }=64.4\pm 1.2\,\mathrm{eV}$; this temperature is consistent with the previous 2009 Chandra observation (${{kT}}_{\mathrm{eff}}^{\infty }=64.5\pm 1.8$ eV) within 1σ, suggesting that the crust of KS 1731−260 has returned to thermal equilibirum with the core. The progression of effective temperature over time (cooling curve) is shown in Figure 2.

Zoom In Zoom Out Reset image size

3.1. Empirical Models

3.2. Physical Models

In order to fit the quiescent crustal cooling, we model the thermal evolution of the neutron star using the open-source code dStar ⁹ (Brown 2015). This code uses the same microphysics and integration scheme discussed in Brown & Cumming (2009). In addition, the code now has a more flexible interface that allows the distribution of heat sources and impurities with depth specified by the user; this allows us to model the effect of an additional heat source in the shallow ocean and an insulating pasta layer in the deep crust. Fixed parameters in dStar include the crust composition, where the composition of Haensel & Zdunik (1990) is used, the atmosphere model with the column depth of the light-element layer is set to 10⁴ g cm⁻², crust-core transition density is 8.13 × 10¹³ g cm⁻³, and the radius and mass are set to 10 km and 1.4 M_⊙ to be self-consistent with the spectral fits. Our dStar models use the superfluid critical temperature in the crust from Schwenk et al. (2003).

Prior to modeling the cooling, we simulate 12.5 years of constant accretion at the rate of $\dot{M}={10}^{17}$ g s⁻¹, which is consistent with the time-averaged accretion rate found by Galloway et al. (2008). For the initial run of models we varied three parameters: core temperature (T_c), the impurity parameter (Q_imp ≡ ${n}_{\mathrm{ion}}^{-1}{\sum }_{i}$n_i(Z_i − $\langle Z\rangle$)²) of the crust, and additional shallow heating of the crust (Q_sh). The impurity parameter measures the distribution of nuclide charge numbers. This influences the conductivity and structure of the crust (i.e., low Q_imp means high thermal conductivity and well-structured lattice). The shallow heating is an additional few MeVs of crustal heating caused by an unknown source (e.g., Cumming et al. 2006; Brown & Cumming 2009; Medin & Cumming 2014; Deibel et al. 2015). In the model, we set the depth of the shallow heating to occur between pressure values of ${P}_{{Q}_{\mathrm{sh}},\min }=1\times {10}^{27}$ g cm⁻¹ s⁻² and ${P}_{{Q}_{\mathrm{sh}},\max }\,=1\times {10}^{28}$ g cm⁻¹ s⁻². The pressure values were chosen so that the heat source was shallower than the depth corresponding to the thermal time of the first observation, but otherwise it was arbitrary. The cooling curve is not that sensitive to these boundaries, since the first observation is 65 days after the end of the outburst. The influence of the T_c, Q_imp, and Q_sh parameters on the shape of the cooling curve can be seen in panels (a)–(c) of Figure 3.

Zoom In Zoom Out Reset image size

We calculate the thermal evolution of the crust for a wide range in all parameters, creating a grid of cooling curve models. This is the first time a full exploration of parameter space has been conducted. To find the best-fit parameters, we search the grid to find the model that fits the data with the lowest χ² value. The best-fitting parameters are T_c = 9.35 ± 0.25 × 10⁷ K, ${Q}_{\mathrm{imp}}={4.4}_{-0.5}^{+2.2}$, and Q_sh = 1.36 ± 0.18 MeV/nucleon. The uncertainties quoted here, and throughout the paper, are at the 1σ level. The best-fitting model is shown as a solid black line in Figure 5. This model gives a reduced-χ² = 2.00 (dof = 5). Our core temperature value is approximately two times greater than the core temperature values found in Cackett et al. (2010) (T_c = 4.6 × 10⁷ K) and Brown & Cumming (2009) (T_c = 5.4 × 10⁷). We find our impurity parameter value is higher than but similar to Cackett et al. (2010) (Q_imp = 4.0) and approximately three times greater than Brown & Cumming (2009) (Q_imp = 1.5). The slightly different core temperature is not surprising since it is strongly constrained by this latest data point, as it is the surface temperature once the crust and core are in thermal equilibrium. Moreover, our choice of light-element column depth is smaller than in Brown & Cumming (2009) (who use 10⁹ g cm⁻²), the effect of which leads to a higher core temperature because the envelope is less opaque (Cumming et al. 2016). The 1σ, 2σ, and 3σ confidence regions of the parameter space can be seen in Figure 4.

Zoom In Zoom Out Reset image size

We ran a second set of physical models that included a low thermal conductivity layer, which is consistent with nuclear pasta. A nuclear pasta phase is expected to exist at the base of neutron star crusts (e.g., see Horowitz et al. 2015 and references therein). Pons et al. (2013) suggested that the irregular shapes in the pasta would lead to low conductivity. Horowitz et al. (2015) showed that there could be defects in the pasta and suggested that they could act as scattering sites that would lower the thermal conductivity. We kept T_c, Q_imp, and Q_sh of the crust as free parameters. We added Q_imp of the pasta layer and the density of the transition to the pasta phase, ρ_pasta, as variable parameters. The pasta layer is insulating. The presence of a pasta phase requires a greater temperature gradient to carry a thermal flux. As the transition density becomes smaller, the pasta layer becomes thicker, which causes the temperature of the crust to increase, and the crust cools more slowly. Additionally, the pasta layer becomes more insulating as the impurity parameter is increased, resulting in a similar effect as the decreased transition density. The influence of the pasta parameters on the cooling curve can be seen in panels (d)–(e) of Figure 3. Again, we create a large grid of models and find the model with the best (lowest) χ² value. The best-fitting parameters are T_c = 9.34 ± 0.21 × 10⁷ K, Q_imp =2.1 ± 1.0, Q_sh = 1.43 ± 0.15 MeV/nucleon, ${Q}_{\mathrm{imp},\mathrm{pasta}}\,=12.4\pm 5.1$, and ρ_pasta = 2.7 ± 0.8 × 10¹³ g cm⁻³. The best-fitting pasta model is shown as a dashed red line in Figure 5. This model gives a reduced-χ² = 3.232 (dof = 3). The 1, 2, and 3σ confidence regions of the parameter space can be seen in Figure 6.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We have analyzed the latest Chandra observation of the neutron star LMXB KS 1731−260, 14.5 years into quiescence and 6 years after the previous Chandra observation. From fitting the X-ray spectrum, we find an effective temperature of ${{{kT}}_{\mathrm{eff}}}^{\infty }$ = 64.4 ± 1.2 eV, which is consistent (within 1σ) of the 2009 Chandra observation which measured ${{{kT}}_{\mathrm{eff}}}^{\infty }$ = 64.5 ± 1.8 eV, implying that the crust has likely thermally relaxed and may have returned to equilibrium with the core. A broken power law fits the cooling curve better than a simple power law, at approximately the 3σ level. Alternatively, if we exclude the latest data point and re-fit a simple power-law decay, we find that the newest observation is 3.5σ away from an extrapolation of the best-fitting model, further suggesting the crust has stopped cooling. We caution, however, that further observations are needed to confirm that this is not a deviation from the cooling curve caused by low-level accretion, as has been seen in, e.g., XTE J1701−462 (Fridriksson et al. 2011). The lack of a power-law component in this latest Chandra observation of KS 1731−260 suggests this may not be happening here.

We constrain crust properties using a physical model of the crust's thermal relaxation. This model calculates the thermal evolution of crust, and we initially allow the core temperature, impurity parameter, and presence of additional shallow heat sources in the crust as free parameters to fit the cooling curve. We found best-fit values for the core temperature and impurity parameter (T_c = 9.35 ± 0.25 × 10⁷ K and ${Q}_{\mathrm{imp}}={4.4}_{-0.5}^{+2.2}$). The low impurity parameter value is consistent with previous discussions that the crust of KS 1731−260 has high thermal conductivity (e.g., Wijnands et al. 2002; Cackett et al. 2006; Shternin et al. 2007; Brown & Cumming 2009; Page & Reddy 2013). The core temperature of the neutron star is most sensitive to cooling at late times. This new late time observation (∼5300 days after the end of outburst), along with the full exploration of parameter space, places significantly improved constraints on the core temperature. In turn, constraints on the core temperature are important for placing a lower limit on the core's heat capacity (Cumming et al. 2016). Moreover, we find that a significant additional shallow heat source in the crust is preferred, with a best-fitting value of Q_sh = 1.36 ± 0.18 MeV/nucleon. This amount of extra heating is consistent with that found in other neutron star LMXBs (e.g., EXO 0748−676 and Aql X-1; Degenaar et al. 2014; Waterhouse et al. 2016), but is significantly less than the amount of shallow heating required for MAXI J0556−332 (Q_sh ≈ 4–10 MeV/nucleon; Deibel et al. 2015).

dStar also allows testing for the presence of a layer of low thermal conductivity material close to the crust-core transition, possibly a nuclear pasta phase. However, we find that the inclusion of such a layer did not have significant influence on the best-fitting model, with the best fit coming from the lowest pasta density we allowed from the model (see contours in Figure 6). This is consistent with Horowitz et al. (2015) that KS 1731−260 can be modeled equally as well with or without a nuclear pasta layer. A parameter we did not vary in our pasta model was the superfluid critical temperature in the crust. The influence of this parameter can be seen in Deibel et al. (2016).

Our modeling with dStar assumes a constant accretion rate throughout the 12.5-year outburst of KS 1731−260. However, a recent investigation by Ootes et al. (2016) has shown that variations in the accretion rate throughout the outburst influence the cooling curve. Especially important are variations at the end of the outburst, which can strongly influence the early part of the cooling curve and hence have a significant impact on the implied amount of extra heating at shallow depths. Their modeling of KS 1731−260 (without including the latest observation) finds the need for a 1.4 MeV/nucleon shallow heat source when variations in accretion rate are included and only 0.6 MeV/nucleon when a constant accretion rate is assumed. (Although, note that Ootes et al. 2016 do not optimize their fits to get the best-fit parameters.) This value for Q_sh is significantly less than the value implied from our modeling, suggesting some dependence on the input parameters that we do not fit for, such as the depth of the shallow heating and crust composition included in the cooling code. For instance, the minimum shallow heating density we use is ${\rho }_{\mathrm{sh},\min }=2.12\times {10}^{9}$ g cm⁻³, while Ootes et al. (2016) use ${\rho }_{\mathrm{sh},\min }=4\times {10}^{8}$ g cm⁻³. The choice of ${\rho }_{\mathrm{sh},\min }$ was used in previous calculations, and we have tested that our results are insensitive to the precise value.

KS 1731−260 has shown superbursts (Kuulkers et al. 2002), type-I X-ray bursts due to thermonuclear burning of carbon, leading to bursts that last for hours rather than tens of seconds. The conditions required for carbon to ignite places constraints on the thermal properties of the crust (Cumming et al. 2006), independent of the constraints placed by crustal cooling. The presence of superbursts requires that the crust has a temperature of ∼6 × 10⁸ K at a depth of ∼10¹² g cm⁻² (Cumming et al. 2006). Without any additional heat source, the crust is too cold (Cumming et al. 2006; Gasques et al. 2007). dStar gives the crust temperature profile at all times; thus we can study the temperature in the crust at the time of the superburst. As discussed in Ootes et al. (2016), the accretion rate during the outburst varies and is significantly higher in the earlier part than at the end. The persistent flux before the superburst (taken from Kuulkers et al. 2002) is about 1.66 times higher than the average mass accretion rate we assume. We therefore recalculate the thermal evolution of the crust during the outburst assuming this higher mass accretion rate, and a correspondingly higher value of Q_sh, and look at the crust temperature profile at a time of 2894 days into the outburst (when the superburst took place). Doing this gives a temperature in the outer part of the crust (close to 10¹² g cm⁻²) of ∼5 × 10⁸ K, which is broadly consistent with the required 6 × 10⁸ K for superburst ignition, within the uncertainties and assumptions made.

KS 1731−260 provides the rare opportunity to place two independent mass–radius constraints on the same source. During outburst, KS 1731−260 displayed photospheric radius expansion (PRE) bursts (Muno et al. 2000). PRE bursts are thought to reach the Eddington luminosity; thus fitting the blackbody emission from the burst can lead to both mass and radius constraints (both the Eddington luminosity and observed emitted radius depend on M and R), assuming the distance to the source is known (Özel 2006, though see, e.g., Poutanen et al. 2014 for an opposing view). Özel et al. (2012) applied this technique to KS 1731−260, implying R ≤ 12.5 km and M ≤ 2.1 M_⊙ (95% confidence level). In future work, we will explore constraints on M and R through fitting crustal cooling models to KS 1731−260.

R.L.M. and E.M.C. gratefully acknowledge support for this work that was provided by the National Aeronautics and Space Administration through Chandra Award Number GO5-16043X issued by the Chandra X-ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of the National Aeronautics Space Administration under contract NAS8-03060. E.F.B. and A.D. are supported by the National Science Foundation under Grant No. AST-1516969. A.C. is supported by an NSERC Discovery Grant. N.D. is supported by an NWO/Vidi grant and an EU Marie Curie Intra-European fellowship (contract no. FP-PEOPLE-2013-IEF-627148). R.W. acknowledges support from a NWO top grant, module 1. D.P. is partially supported by the Consejo Nacional de Ciencia y Tecnología with a CB-2014-1 grant #240512. This work was enabled in part by the National Science Foundation under Grant No. PHY-1430152 (JINA Center for the Evolution of the Elements). We also thank the International Space Science Institute, Bern, for hosting the Neutron Star Crust Team, where many productive discussions on this topic took place.