Diffusive nuclear burning in cooling simulations and application to new temperature data of the Cassiopeia A neutron star

TL;DR

This paper models diffusive nuclear burning in neutron star envelopes, deriving temperature relations, and applies these models to interpret long-term cooling data of the Cassiopeia A neutron star, highlighting the importance of envelope composition evolution.

Contribution

It introduces analytic temperature relations incorporating DNB effects into neutron star cooling models and applies them to observational data of Cassiopeia A.

Findings

DNB prevents stable hydrogen buildup in envelopes.

DNB makes helium envelopes more heat-transparent at high temperatures.

Envelope composition evolution affects long-term neutron star cooling.

Abstract

A critical relation in the study of neutron star cooling is the one between surface temperature and interior temperature. This relation is determined by the composition of the neutron star envelope and can be affected by the process of diffusive nuclear burning (DNB), which occurs when elements diffuse to depths where the density and temperature are sufficiently high to ignite nuclear burning. We calculate models of H-He and He-C envelopes that include DNB and obtain analytic temperature relations that can be used in neutron star cooling simulations. We find that DNB can lead to a rapidly changing envelope composition and prevents the build-up of thermally stable hydrogen columns y$_H$ > 10$^{7}$ g cm$^{-2}$, while DNB can make helium envelopes more transparent to heat flux for surface temperatures $T_s$ > 2 $\times 10^6$ K. We perform neutron star cooling simulations in which we evolve temperature and envelope composition, with the latter due to DNB and accretion from the interstellar medium. We find that a time-dependent envelope composition can be relevant for understanding the long-term cooling behaviour of isolated neutron stars. We also report on the latest Chandra observations of the young neutron star in the Cassiopeia A supernova remnant; the resulting 13 temperature measurements over more than 18 years yield a ten-year cooling rate of $\approx$ 2%. Finally, we fit the observed cooling trend of the Cassiopeia A neutron star with a model that includes DNB in the envelope.

Figures

Tables

Table 1. Table 1: Surface temperature T s , 6 subscript 𝑇 s 6 T_{\rm s,6} ( 10 6 ​  K superscript 10 6  K 10^{6}\mbox{ K} ), absorbed 0.5–7 keV flux F − 13 abs superscript subscript 𝐹 13 abs F_{\rm-13}^{\rm abs} ( 10 − 13 ​  erg cm -2  s -1 superscript 10 13  erg cm -2  s -1 10^{-13}\mbox{ erg cm${}^{-2}$ s${}^{-1}$} ), and absorption column N H , 22 subscript 𝑁 H 22 N_{\rm H,22} ( 10 22 ​  cm -2 superscript 10 22  cm -2 10^{22}\mbox{ cm${}^{-2}$} ) determined from model fits to Chandra  ACIS-S GRADED spectra of the neutron star in the Cassiopeia A supernova remnant. Three spectral fit results are shown, two with constant N H subscript 𝑁 H N_{\rm H} and one with changing N H subscript 𝑁 H N_{\rm H} . Each set of 13 temperatures and fluxes are fit to a linear decline, with decline rate and fit statistic as shown. For merged ObsIDs, the MJD listed is that of the first ObsID. Number in parentheses is 1 ​ σ 1 𝜎 1\sigma uncertainty in last digit.

		$N_{H,22}$	1.73		1.67		see below
		$M$ ($M_{\mathrm{Sun}}$)	1.65		1.65		1.65
		$R$ (km)	10.3		12.9		12.9
		${\chi }^2$/dof	1562/1449		1522/1449		1498/1436
ObsID	Date	MJD	$T_{s,6}$	$F_{-13}^{\mathrm{abs}}$	$T_{s,6}$	$F_{-13}^{\mathrm{abs}}$	$T_{s,6}$	$F_{-13}^{\mathrm{abs}}$	$N_{H,22}$
114	2000 Jan 30	51573.4	2.127(10)	7.4(2)	1.867(8)	7.3(1)	1.877(10)	7.3(2)	1.73(3)
1952	2002 Feb 6	52311.3	2.126(10)	7.4(2)	1.867(8)	7.3(1)	1.876(10)	7.4(2)	1.72(3)
5196	2004 Feb 8	53043.7	2.107(10)	7.1(1)	1.851(8)	7.0(2)	1.849(10)	7.0(2)	1.66(3)
9117/9773	2007 Dec 5/8	54439.9	2.098(9)	7.0(2)	1.841(8)	6.9(2)	1.849(11)	6.9(2)	1.72(4)
10935/12020	2009 Nov 2/3	55137.9	2.090(10)	6.9(2)	1.834(8)	6.8(2)	1.839(11)	6.8(1)	1.70(4)
10936/13177	2010 Oct 31/Nov 2	55500.2	2.081(10)	6.8(2)	1.826(8)	6.7(2)	1.819(11)	6.6(1)	1.63(4)
14229	2012 May 15	56062.4	2.046(9)	6.4(2)	1.797(8)	6.3(1)	1.804(11)	6.4(1)	1.71(4)
14480	2013 May 20	56432.6	2.064(9)	6.6(2)	1.813(8)	6.5(1)	1.817(10)	6.5(1)	1.69(4)
14481	2014 May 12	56789.1	2.047(9)	6.3(2)	1.799(7)	6.2(1)	1.806(10)	6.3(1)	1.71(4)
14482	2015 Apr 30	57142.5	2.066(9)	6.7(1)	1.817(8)	6.6(1)	1.801(10)	6.6(1)	1.58(4)
19903/18344	2016 Oct 20/21	57681.2	2.061(9)	6.6(2)	1.810(8)	6.5(2)	1.795(11)	6.5(1)	1.58(4)
19604	2017 May 16	57889.7	2.052(9)	6.5(2)	1.804(7)	6.4(2)	1.802(10)	6.4(2)	1.66(4)
19605	2018 May 15	58253.7	2.046(9)	6.4(2)	1.800(7)	6.4(1)	1.783(10)	6.4(1)	1.56(4)
		10-year
		decline rate	$2.2\pm 0.2\%$	$7.5\pm 1.0\%$	$2.1\pm 0.2\%$	$7.3\pm 0.9\%$	$2.7\pm 0.3\%$	$7.7\pm 0.9\%$
		${\chi }^2$/dof	14.8/11	11.0/11	15.7/11	15.8/11	6.2/11	10.9/11

Table 2. Table 2: Best fit parameters for the H-He mixture with ρ b subscript 𝜌 𝑏 \rho_{b} = 10 8 g cm -3 (Equations  10 ).

$p_1$	$p_2$	$p_3$	$p_4$	$p_5$	$p_6$	$p_7$	$p_8$
0.8254	2.086	18.70	-1.025	3.346	1.589	0.03848	1.572

Table 3. Table 3: Best fit parameters for the He-C mixture with ρ b subscript 𝜌 𝑏 \rho_{b} = 10 10 g cm -3 (Equations  11 ).

$p_1$	$p_2$	$p_3$	$p_4$	$p_5$	$p_6$	$p_7$	$p_8$
10.927	0.977	2.925	3.521	0.0346	1.650	963.474	-19.318
$p_9$	$p_{10}$	$p_{11}$	$p_{12}$	$p_{13}$	$p_{14}$	$p_{15}$	$p_{16}$
-0.846	2.931	1.190	-0.175	3.110 $\times$ 103	10.812	2.181 $\times$ 10-9	2.314

Table 4. Table 4: Surface temperature T s , 6 subscript 𝑇 s 6 T_{\rm s,6} ( 10 6 ​  K superscript 10 6  K 10^{6}\mbox{ K} ), absorbed 0.5–7 keV flux F − 13 abs superscript subscript 𝐹 13 abs F_{\rm-13}^{\rm abs} ( 10 − 13 ​  erg cm -2  s -1 superscript 10 13  erg cm -2  s -1 10^{-13}\mbox{ erg cm${}^{-2}$ s${}^{-1}$} ), and absorption column N H , 22 subscript 𝑁 H 22 N_{\rm H,22} ( 10 22 ​  cm -2 superscript 10 22  cm -2 10^{22}\mbox{ cm${}^{-2}$} ) determined from model fits using CALDB 4.7.8 (in comparison to Table  1 which results from using CALDB 4.8.1) to Chandra  ACIS-S GRADED spectra of the neutron star in the Cassiopeia A supernova remnant. Three spectral fit results are shown, two with constant N H subscript 𝑁 H N_{\rm H} and one with changing N H subscript 𝑁 H N_{\rm H} . Each set of 13 temperatures and fluxes are fit to a linear decline, with decline rate and fit statistic as shown. For merged ObsIDs, the MJD listed is that of the first ObsID. Number in parentheses is 1 ​ σ 1 𝜎 1\sigma uncertainty in last digit.

		$N_{H,22}$	1.73		1.68		see below
		$M$ ($M_{\mathrm{Sun}}$)	1.65		1.65		1.65
		$R$ (km)	10.3		13.0		13.0
		${\chi }^2$/dof	1550/1450		1512/1450		1486/1437
ObsID	Date	MJD	$T_{s,6}$	$F_{-13}^{\mathrm{abs}}$	$T_{s,6}$	$F_{-13}^{\mathrm{abs}}$	$T_{s,6}$	$F_{-13}^{\mathrm{abs}}$	$N_{H,22}$
114	2000 Jan 30	51573.4	2.127(10)	7.4(2)	1.862(8)	7.3(2)	1.870(10)	7.3(2)	1.73(3)
1952	2002 Feb 6	52311.3	2.125(10)	7.4(2)	1.861(8)	7.3(2)	1.869(10)	7.3(2)	1.73(3)
5196	2004 Feb 8	53043.7	2.106(10)	7.1(1)	1.845(8)	7.0(2)	1.842(10)	7.0(1)	1.66(3)
9117/9773	2007 Dec 5/8	54439.9	2.097(9)	7.0(2)	1.835(8)	6.9(2)	1.842(10)	6.9(2)	1.72(4)
10935/12020	2009 Nov 2/3	55137.9	2.089(10)	6.9(2)	1.828(8)	6.8(2)	1.832(11)	6.8(2)	1.71(4)
10936/13177	2010 Oct 31/Nov 2	55500.2	2.080(10)	6.7(1)	1.820(8)	6.6(2)	1.812(11)	6.6(1)	1.64(4)
14229	2012 May 15	56062.4	2.044(9)	6.4(2)	1.791(7)	6.3(1)	1.797(11)	6.3(2)	1.72(4)
14480	2013 May 20	56432.6	2.062(9)	6.5(1)	1.806(8)	6.4(2)	1.810(10)	6.5(1)	1.70(4)
14481	2014 May 12	56789.1	2.045(9)	6.3(1)	1.792(7)	6.2(2)	1.799(10)	6.2(2)	1.72(4)
14482	2015 Apr 30	57142.5	2.064(9)	6.6(2)	1.809(7)	6.6(2)	1.795(10)	6.5(2)	1.60(4)
19903/18344	2016 Oct 20/21	57681.2	2.059(9)	6.6(2)	1.803(8)	6.5(1)	1.789(11)	6.4(1)	1.60(4)
19604	2017 May 16	57889.7	2.050(9)	6.5(2)	1.798(7)	6.4(2)	1.795(10)	6.4(1)	1.67(4)
19605	2018 May 15	58253.7	2.050(9)	6.5(2)	1.798(7)	6.5(2)	1.777(10)	6.4(1)	1.53(4)
		10-year
		decline rate	$2.2\pm 0.2\%$	$7.3\pm 1.0\%$	$2.1\pm 0.2\%$	$7.2\pm 1.0\%$	$2.6\pm 0.3\%$	$7.6\pm 0.9\%$
		${\chi }^2$/dof	15.8/11	12.9/11	16.8/11	15.4/11	6.3/11	11.1/11

Equations

$$\frac{\partial T}{\partial z}=-\frac{3\kappa \rho }{16T^3}{T}_s^4,$$

$$J_{\text{light}}=-\mathcal{D}\frac{d{n}_{\text{light}}}{dz}+n_{\text{light}}{w}_{\text{light}}=n_{\text{light}}{v}_{\text{light}},$$

$$\mathcal{D}\approx 10^{-3}\frac{A_{\text{light}}^{0.1}{T}_6^{1.3}}{Z_{\text{light}}^{1.3}{Z}_{\text{heavy}}^{0.3}{\rho }_5^{0.6}}{\mathrm{cm}}^2{\mathrm{s}}^{-1},$$

$$\frac{{\partial }^2{f}_{\text{light}}}{\partial z^2}+\frac{\partial f_{\text{light}}}{\partial z}\frac{\partial \text{ln}{n}_{\text{light},0}}{\partial z}=\frac{f_{\text{light}}}{\mathcal{D}{\tau }_{\text{light}}},$$

$${\zeta }_{\mathrm{DNB}}=\frac{y_{\text{light}}}{{\tau }_{\text{light,col}}}=\int dz\frac{A_{\text{light}}{n}_{\text{light}}{m}_p}{{\tau }_{\text{light}}(n_{\text{light}},n_{\text{sub}},T)},$$

$${\dot{y}}_{\text{acc}}\approx {\dot{M}}_{\text{acc}}\frac{g_s}{4\pi GM}.$$

$$\Delta y_{\text{light}}=\frac{d{y}_{\text{light}}}{dt}dt=({\dot{y}}_{\text{acc}}-{\zeta }_{\mathrm{DNB}})dt,$$

$$\dot{M}=\frac{4\pi {\lambda }_a{G}^2{M}^2{\rho }_{\mathrm{ISM}}}{{\left(v^2+c_s^2\right)}^{3/2}},$$

$$\begin{split}T_{b}(Y,\rho*)=10^{7}K\times\Biggl{(}f_{4}(Y)+[f_{1}(Y)-f_{4}(Y)]\\ \times\left[f_{5}(Y,\rho*)+\left(\frac{\rho*}{f_{2}(Y)}\right)^{f_{3}(Y)}\right]^{p_{12}}\Biggr{)}\end{split}$$

$$\begin{array}{rl}{f}_1(Y)& =p_1{Y}^{p_2}\sqrt{1+p_3{Y}^{p_4}}\\ f_2(Y)& =p_5{Y}^{p_6}\sqrt{1+p_7{Y}^{p_8}}\\ f_2(Y)& =\frac{p_9{Y}^{p_{1}0}}{{\left(1-p_{11}Y+p_{12}{Y}^2\right)}^2}\\ f_3(Y)& =p_{13}{Y}^{-p_{14}}\\ f_5(Y)& =1\end{array}$$

$$\begin{array}{rl}{f}_1(Y)& =p_1{Y}^{p_2{\text{log}}_{10}Y+p_3}\\ f_2(Y)& =p_7{Y}^{p_8({\text{log}}_{10}Y{)}^2+p_9}\\ f_3(Y)& =p_{10}\sqrt{\frac{Y}{Y^2+p_{11}^2}}\\ f_4(Y)& =p_4{Y}^{p_5{\text{log}}_{10}Y+p_6}\\ f_5(Y,\rho \ast )& =p_{13}{Y}^{p_{14}}+p_{15}\rho {\ast }^{p_{16}}\end{array}$$

Full text

††pubyear: 2018††pagerange: Diffusive nuclear burning in cooling simulations and application to new temperature data of the Cassiopeia A neutron star–B

1 Introduction

The physical properties, such as composition and structure, of neutron stars (NSs) are still uncertain. Valuable information about the neutron star interior can be obtained by comparing observations of thermal radiation of cooling NSs with theoretical models of how a NS cools over time. Isolated NSs cool via neutrino and photon emission after they are formed hot in a supernova explosion (see review by Potekhin et al. 2015). Additionally, cooling models can be used to study the cooling of accretion-heated NSs in binary systems (for a recent review see Wijnands et al. 2017).

Observations of several isolated neutron stars suggest that the envelopes of younger pulsars ($<$10${}^{4}-$105 yr) may consist of light elements such as hydrogen or helium, whereas the envelopes of older pulsars ($>$105 yr) may be composed of heavier elements (see, e.g., Yakovlev et al. 2001; Potekhin 2014). Evidence for the presence of carbon in a NS envelope was found for the first time using spectral observations of the young NS in Cassiopeia A (Ho & Heinke, 2009). Since then, several other NS spectra are found to be fit by carbon atmosphere spectra (see De Luca 2017, and references therein). These observations indicate a possible evolution of the composition of the envelope through nuclear burning, after the initial composition is set by fallback material onto the NS after a supernova explosion. The NS in Cassiopeia A (we hereafter refer to the NS as Cas A) is an intriguing object not only because of the chemical composition of its envelope, but also because surface temperatures derived from X-ray observations suggest that the source might be cooling rapidly (Heinke & Ho 2010; c.f. Posselt et al. 2013). Constraining cooling models with multiple temperature observations can provide fundamental knowledge on the properties of the NS interior, such as critical temperatures for superconductivity and superfluidity (e.g., Shternin et al. 2011; Page et al. 2011; Ho et al. 2015). Other possibilities include fast or rotationally induced neutrino cooling (Negreiros et al., 2013; Taranto et al., 2016), magnetic field decay (Bonanno et al., 2014), slow thermal relaxation (Blaschke et al., 2012, 2013), stellar fluid oscillations (Yang et al., 2011), and transition to axions or quark matter (Noda et al., 2013; Sedrakian, 2013; Leinson, 2014; Hamaguchi et al., 2018).

A major source of uncertainty in cooling models is the chemical composition of the thin outer envelope which extends from just below the surface down to a boundary mass density of $\rho_{b}\sim 10^{8}-10^{10}$ g cm*-3*. The envelope acts as a thermal insulator due to its relatively poor thermal conductivity and thus sets a strong temperature gradient between the surface temperature Ts and the interior temperature at the bottom of the envelope Tb [$\equiv T(\rho_{b})$]. The composition affects the thermal conductivity and thus how transparent the envelope is to the heat flux. This means that a measured surface temperature can correspond to different interior temperatures depending on the composition of the envelope.

Several models of the heat blanketing envelope have been developed after the work by Gudmundsson et al. (1983), who calculated the thermal structure and conductivity of an iron envelope. Later models went beyond the single element envelope and considered the envelope as multiple shells of light elements (H, He, C) separated by an abrupt boundary (Potekhin et al., 1997, 2003). These models have subsequently been used in many interior studies. Most recently, Beznogov et al. (2016) developed new models for binary ion mixtures (H-He, He-C and C-Fe) where ion diffusion is taken into account but the possible effect of nuclear burning due to diffusion is not included.

Numerical calculations of the neutron star envelope including the effects of diffusive nuclear burning (DNB) have been performed by Chang & Bildsten (2003, 2004); Chang et al. (2010). The key idea of DNB is that the tail of the light element layer may diffuse down to depths where the density and temperature is large enough to ignite burning. The depletion of the light elements at this depth can subsequently drive a light element current from the surface to the burning layer. In Chang & Bildsten (2003) and Chang & Bildsten (2004), a H-C envelope was modelled, and it was found that diffusive nuclear burning can readily consume all of the initial hydrogen in the NS envelope within $\sim$105 yr. Motivated by the detection of carbon on the surface of Cas A, Chang et al. (2010) calculated models for DNB in He-C envelopes and found that any primordial helium is likely depleted during the hot ($T_{b}$> 108 K), early evolution of the NS.

Typically, NS cooling models calculate the NS structure and thermal evolution from the centre out to the bottom of the envelope, where they are coupled to theoretical envelope temperature relations (between Tb and Ts) based on the envelope models above. Diffusive nuclear burning can change the composition of the envelope and hence its thermal conductivity. Therefore, it can affect the temperature relations between Ts and Tb. We use the work of Chang & Bildsten (2003, 2004); Chang et al. (2010) and calculate models including diffusive nuclear burning for H-He and He-C envelopes and provide analytic envelope relations (based on those presented in Beznogov et al. 2016) that can be used in cooling models. We use these temperature relations to investigate the effect of diffusive nuclear burning and time-variable envelopes on neutron star cooling. We also report on four new Chandra observations of the NS in Cassiopeia A resulting in the longest spanning ( $>$18 years) dataset currently available. We then use the new temperature relations including DNB to model the observed cooling of Cas A.

2 Theory of diffusive nuclear burning

2.1 The physics of DNB

We summarize the key physics input for diffusive nuclear burning as calculated in Chang et al. (2010). For a full description of the physics we refer the reader to this work. The neutron star envelope, at $\rho$ < 10${}^{8}-$1010 g cm*-3*, has a thickness ($\sim$ a few hundred metres) that is much smaller than the stellar radius and can be well approximated as a plane-parallel layer with a constant downward gravitational acceleration $g$. The thermal profile is obtained by solving the heat diffusion equation for a constant flux in radiative equilibrium

[TABLE]

where $z$ is the depth, $T_{s}$ is the effective (non-redshifted) surface temperature and $\kappa$ is the opacity. In the outermost layer of the envelope the radiative opacity is dominated by free-free absorption and Thomson scattering. We use the radiative opacities of Potekhin & Yakovlev (2001) and the analytic formulae for the conductive opacity given by Potekhin et al. (1999).

We consider the case where a trace ion species, hereafter referred to as ‘light elements’, is immersed in a fixed background ‘heavy’ ion species. The equation for hydrostatic equilibrium, which is dominated by electric and gravitational forces (see Chang et al. 2010) is used to obtain the number density profiles of ion species in the absence of nuclear burning. We require overall charge neutrality such that, $n_{e}~{}=~{}Z_{\text{light}}n_{\text{light}}+Z_{\text{heavy}}n_{\text{heavy}}$, where $n_{e}$ is the electron number density, and $Z_{\text{light}}$, $n_{\text{light}}$, $Z_{\text{heavy}}$ and $n_{\text{heavy}}$ are the charge number and number density of the light (trace) and heavy (background) element species, respectively. The diffusive tail of the light element layer may extend deep down into the underlying heavier element layer. When lighter elements diffuse to depths where they are consumed by nuclear burning, the depletion of light elements at this depth induces a diffusive light element current from the surface to the burning layer. Because the timescale at which the local column density changes is much longer than the local nuclear timescale, a steady state approximation is used for the light element current induced by DNB ($\partial n_{\text{light}}/\partial t=0$). The light element current is defined as

[TABLE]

where $v_{\text{light}}$ and $w_{\text{light}}$ are the light element relative velocity and drift velocity, respectively. The diffusion coefficient $\mathcal{D}$ is calculated as in Chang & Bildsten (2004) (similar to Brown et al. 2002)

[TABLE]

where $T_{6}$ = $T/$106 K and $\rho_{5}$ = $\rho$/105 g cm*-3*. The atomic mass and charge number are denoted by $A$ and $Z$, respectively. From the light element continuity equation and the definition of the current, one can derive the general case of the steady-state diffusion equation (which is derived in detail in Chang & Bildsten 2004)

[TABLE]

where $n_{\text{light},0}$ is the number density of the trace ion in the absence of nuclear burning obtained by solving the hydrostatic balance equations numerically, $f_{\text{light}}$ is the correction factor to the number density due to nuclear burning (i.e., $f_{\text{light}}\equiv n_{\text{light}}/n_{\text{light},0}$) and $\tau_{\text{light}}$ is the local lifetime against nuclear capture set by the local nuclear burning rate. With increasing depth and temperature, the ion scale height (=$[d$ln n_{\text{light,0}}/dz$$]^{-1}) rapidly becomes larger than the nuclear scale height (=$\sqrt{\mathcal{D}\tau_{\text{light}}})$. Therefore, at the burning layer where the second term can effectively be dropped, Equation 4 represents a diffusion equation for $f_{\text{light}}$ with a nuclear-driven source.

In our models, hydrogen can be consumed steadily by the proton capture reactions of the CNO-cycle and pp chain reactions. Helium is consumed by proton- and alpha capture reactions and the triple alpha process, although the reaction rates of the triple alpha process were found to be negligible in the parameter space of interest (see Chang et al. 2010) and thus they are ignored. Nuclear reaction rates and experimental reaction cross section values are obtained from the NACRE111http://pntpm.ulb.ac.be/Nacre/nacre.htm (Angulo et al., 1999) and REACLIB222http://nucastro.org/reaclib.html database compilations.

2.2 Model description

We model the neutron star envelope for 1H-4He and 4He-12C mixtures following the work of Beznogov et al. (2016), but including the physics of DNB. These envelope compositions were chosen to be able to compare our solutions with those found in Beznogov et al. (2016), but in reality other ions may be present in neutron star envelopes. In particular, one might consider a H-He-C envelope, with an intervening helium layer between the hydrogen and carbon components. The effects of such an envelope are strongly dependent on the assumed size of each layer. For simplicity in the present work, we limit our study of DNB to two component envelopes. We use a surface gravity of $g_{s}$ = 2.43 $\times$ 1014 cm s*-2* in our calculations, but the resulting temperature relations at other surface gravities can be approximated using the temperature scaling relation by Gudmundsson et al. (1983).

We start our numerical calculations with a given surface temperature $T_{s}$ and fraction $x_{\text{light}}$ of the light element species present at the photosphere. The partial pressure $p_{\text{light}}$ and column density $y_{\text{light}}$ of the ion components at the photosphere is set by $x_{\text{light}}$ as ycol = $y_{\text{light}}+y_{\text{bkg}}$ = $x_{\text{light}}\times y_{col}+x_{\text{bkg}}\times y_{col}$. We then integrate the system of equations in Section 2.1 to obtain profiles for the internal temperature, ion densities and nuclear burning rates. A given photospheric light element fraction can correspond to different values of the total light element column density depending on the surface temperature. Therefore, we perform a large number of numerical calculations for varying $T_{s}$ and $x_{\text{light}}$ to obtain sufficient coverage in the resulting $T_{s}$-$T_{b}$-$y_{\text{light}}$ parameter space that is used for fitting analytic functions.

We note that in the analytic temperature relations calculated by Beznogov et al. (2016), the amount of light elements is characterized by their effective transition density $\rho*$. This is an artificially chosen density at which a total mass of light elements is contained in the outer shell $\rho<\rho*$. Note that the light element column density $y_{\text{light}}$ can be converted to $\rho*$ making use of the relevant equation of state.

2.3 Parameter ranges

We fit the numerical results of the models described above to obtain analytic $T_{s}-T_{b}$ relations that can be coupled to thermal evolution codes for the neutron star interior. Note that the $T_{s}-T_{b}$ relations are actually relations for $T_{s}$, $T_{b}$ and $y_{\text{light}}$ (or equivalently, $\rho*$) for a given surface gravity, i.e. Tb($T_{s}$, $y_{\text{light}}$). We are limited in our parameter ranges by both physical and numerical constraints. Physical constraints arise for example from the fact that the light element column density cannot be larger than the envelope-crust boundary column density $y_{b}$ or smaller than the photospheric density. For the H-He envelopes, we let $\rho_{b}$ = 108 g cm*-3* and increase $\rho^{*}$ up to $\sim$107 g cm*-3*. For He-C envelopes we use $\rho_{b}$ = 1010 g cm*-3* and let $\rho^{*}$ vary up to $\sim$108 g cm*-3*. The surface temperature range for all compositions is $0.3\times 10^{6}$ K$<T_{s}<3\times 10^{6}$ K, such that it encompasses the observed temperatures for both isolated neutron stars and transient neutron stars in quiescence.

3 Burning rates

The reaction rates due to DNB depend on the composition of the envelope, the amount of burning material available and the temperature in the burning region. Using the nuclear reactions that consume hydrogen and helium discussed in Section 2, the total burning rate is

[TABLE]

where $A_{\text{light}}$ is the atomic number of the light element species, $\tau_{\text{light,col}}$ is the characteristic time for the total integrated light element column to be consumed and $\tau_{\text{light}}$ is the local lifetime of the light element species due to DNB. The local number densities of the $\alpha$- or $p$-capturing substrate are denoted by $n_{\text{sub}}$.

In Figure 1 we plot the resulting burning rates for several values of ylight for both the H-He and He-C envelope models. For hydrogen burning (dashed lines in Figure 1), we see that the burning rates are more sensitive to the amount of material available for burning (as characterised by the hydrogen column density) than to the temperature (characterised by the surface temperature). The strong dependence on column density is not surprising, as the burning region is located in the exponentially suppressed diffusive hydrogen tail. Interestingly, Figure 1 shows that DNB can alter the hydrogen column density significantly for the entire range of plotted Ts. Especially for intermediate and large hydrogen column densities (y${}_{\rm H}\gtrsim$ 105 g cm*-2*) the hydrogen column density lifetime is $<<$ 104 yr. It is clear that the helium burning rates drop off rapidly with decreasing temperatures, which is expected as the dominant burning processes are the highly temperature sensitive $\alpha$-capturing reactions. We find that the helium burning rate effectively switches off for T${}_{s}~{}\lesssim$ 106 K. At temperatures T${}_{s}~{}\gtrsim$ 106 K the burning rate is sufficiently large to change the helium column density within years.

Using the burning rates calculated above, we illustrate the evolution of the light element column over time for a fixed surface temperature in Figure 2 and 3 for hydrogen and helium, respectively. We use an initial hydrogen column density of ylight = 106 g cm*-2* and vary the surface temperature between $0.4\times 10^{6}<$ $T_{\rm s}<3.2\times 10^{6}$ K. For all temperatures, the sharpest drop in hydrogen column density occurs within the first $\sim$ 10 years due to the strong density dependence of the burning rate. After that time, the column decreases more gradually. For helium burning we use an initial helium column density of ylight = 1010 g cm*-2* and vary the surface temperature between $2\times 10^{6}<$ $T_{\rm s}<3.2\times 10^{6}$ K. It is clear that the helium column change is more strongly coupled to the temperature and that for $T_{\rm s}\lesssim$ 106 K the column does not change significantly within 100 years. At higher temperatures, Figure 3 illustrates that DNB can even substantially alter the helium column on relatively short timescales of $<10-100$ years.

4 Static envelope temperature relations

We present analytic temperature relations that include diffusive nuclear burning in the envelope. The fit functions and parameters for each of the mixtures are given in Appendix A. In this section, we review their main characteristics and implications.

The temperature relations are most sensitive to the composition in a small region in the envelope which is commonly referred to as the sensitivity strip. The location of the sensitivity strip in the envelope is determined by where the radiative $\kappa_{r}$ and conductive $\kappa_{c}$ opacity are comparable, i.e., $\kappa_{r}\approx\kappa_{c}$. Even a small amount of light elements in this region can have a large effect on the resulting temperature profile in the envelope. The composition outside of the sensitivity strip barely affects the temperature relation. Therefore, in order to understand the effect that DNB has on the temperature relations, $T_{\rm b}(T_{\rm s},y_{\rm light})$, we will consider for which values of $y_{\rm light}$ the composition in the sensitivity strip is altered by DNB.

In Figure 4, we show the $T_{\rm s}$-$y_{\rm H}$ parameter space for a H-He envelope; it is important to note that what is illustrated in Figures 4 and 5 is light element column density, not total column density. For small hydrogen columns and high surface temperatures (lower right blue region), the bulk of the hydrogen remains close to the surface and does not reach the sensitivity strip; thus the temperature relations resemble those of a pure helium envelope and are only a function of $T_{\rm s}$. The region enclosed by solid lines indicates temperatures and hydrogen columns which affect the composition of the sensitivity strip and thus the resulting $T_{\rm b}$. For large hydrogen columns (red region), hydrogen makes up most of the composition in the sensitivity strip, and increasing the column further does not change the composition inside the sensitivity strip. In this case, the temperature relations resemble a pure hydrogen envelope and are unaffected by the size of the hydrogen column. Hydrogen DNB occurs in the diffusive tail and the reaction rates are highly sensitive to the amount of material available at the burning depth (see Section 3). Therefore, DNB does not impact the composition distribution in the sensitivity strip for small hydrogen columns. For large H columns (y${}_{\rm H}\gtrsim 10^{7}$ g cm*-2*), the time for hydrogen to diffuse down to the burning depth is shorter than the rate at which it is consumed. In other words, the burning reactions are not limited by the diffusion time but by the nuclear reaction rate. Thus, hydrogen at the burning depth is consumed rapidly, which prevents the build up of large thermally stable hydrogen columns ($y_{\rm H}~{}\gtrsim$ 107 g cm*-2*).

In Figure 5, we show $T_{\rm s}$-$y_{\rm He}$ for the He-C envelope model. Diffusive helium burning takes place deep in the envelope (108 g cm*-2* $\lesssim$ y $\lesssim$ 3 $\times$ 1012 g cm*-2*) and reaction rates are highly sensitive to the temperature (see Section 3). At $T_{\rm s}<$ 1.25 $\times$ 106 MK (denoted by the dashed line in Figure 5), burning rates are too low to affect the envelope composition and temperature relations. Above this surface temperature, DNB prevents the helium column from extending deeper into the envelope beyond the burning depth. As a result, more helium resides at lower densities in the sensitivity strip compared to the case when DNB is not taken into account. Thus a He-C envelope is more transparent to heat flux with DNB for sufficiently high temperatures.

Figures 4 and 5 emphasise that the effect of hydrogen on the conductivity of the envelope is very different than that of helium, as first noted by Beznogov et al. (2016). Typically, increasing the amount of light elements leads to a higher thermal conductivity of the envelope, such that for a given $T_{\rm s}$ the corresponding $T_{\rm b}$ is lower. An exception to this is the case of hydrogen, which has the opposite behaviour, as its larger radiative opacities results in a lower thermal conductivity. Therefore, increasing the amount of hydrogen, leads to a larger $T_{\rm b}$.

The resulting $T_{\rm s}-T_{\rm b}$ relations for a H-He envelope including DNB are shown in Figure 6. For comparison, the temperature relations where DNB is neglected are also plotted. It is clear that including DNB has a negligible effect on the thermal conductivity of the H-He envelope as the resulting temperature relations are indistinguishable. However, including DNB in the H-He envelope leads to an upper limit on the size of the hydrogen column, making the range in $T_{\rm b}$ for a given $T_{\rm s}$ smaller than when DNB is ignored (the difference is indicated by the shaded region in Figure 6).

For a He-C envelope, the temperature relations are shown in Figure 7. As expected, these relations are unaffected by DNB at low temperatures. At higher temperatures ($T_{\rm s}>2\times 10^{6}$ K; see Section 3) where helium burning is significant, DNB changes the temperature relations. In particular, for a given surface temperature, the boundary temperature is lower with DNB; equivalently, the same boundary temperature corresponds to a higher surface temperature. The He-C envelope is a better conductor and more transparent to the heat flux at large temperatures when DNB is taken into account since helium burning at large depths implies a greater helium abundance at shallower depths (in the sensitivity strip) for a given helium column (see above).

5 Time variable envelopes

Envelope temperature relations, such as the ones presented in Section 4, are used to calculate the NS temperature evolution while assuming that the envelope composition does not change (see, e.g., Potekhin et al. 2015 for a review of NS cooling). The burning rates in Section 3 indicate that modelling the envelope with a fixed composition over the course of the NS lifetime can be inaccurate. In this section, we explore time-variable envelope compositions due to DNB and their effect on NS temperature evolution. Instead of calculating a cooling curve with a fixed envelope composition, we use a fixed initial composition and evolve both heat diffusion in the star and envelope composition.

5.1 Cooling model

We model the thermal evolution of the NS using the general relativistic cooling code NSCool333An older version of this code is available at http://www.astroscu.unam.mx/neutrones/NSCool/ (Page, 2016) for a NS with mass $M$ = 1.4 $M_{\rm Sun}$ and radius $R$ = 11.5 km, which is consistent with a surface gravity g_{\mathrm{s}}=1.8\times 10^{14}\mbox{ cm s{}^{-2}}. We use the nuclear A18+$\delta\nu$+UIX* equation of state (Akmal et al., 1998). Neutron superfluidity is taken into account using the ‘SFB’ 1S0 gap model Schwenk et al. (2003), and proton superconductivity is modelled using the ‘CCDK’ ${}^{1}S_{0}$ gap model (Chen et al., 1993). Possible pairing of neutrons in the spin triplet-state is not considered in this section but will be explored in the next one. While uncertainties in the superfluid properties of neutrons and superconducting properties of protons affect the cooling curves, we do not expect this to significantly change our results, as we use the same assumptions for all cooling curves and are only interested in the differences resulting from varying envelope compositions. The energy transport and conservation equations are solved up to the outer boundary defined as the bottom of the envelope: \rho_{\rm b}=10^{8}\mbox{ g cm{}^{-3}} and \rho_{\rm b}=10^{10}\mbox{ g cm{}^{-3}} for a H-He and He-C envelope, respectively. We use the $T_{\rm s}-T_{\rm b}$ models including DNB as described in Section 4 to relate the surface temperature and the boundary temperature.

5.2 Variable envelope composition

We use fits of the numerical burning rates to calculate the evolution of the light element column density $y_{\rm light}(t)$. In addition to solving equations describing the interior thermal profile, we evaluate the change in $y_{\rm light}$ at each time step. We model two competing processes that affect the evolution of column density. On the one hand, the amount of light elements decreases due to DNB when there is enough material available to burn and the temperature is sufficiently large. On the other hand, the amount of light elements increases due to accretion from the interstellar medium (ISM) or from a companion star. A mass accretion rate $\dot{M}_{\text{acc}}$ is equivalent to an increasing light element column density $\dot{y}_{\text{acc}}$ via

[TABLE]

The total change in y${}_{\text{light}}$ is obtained by combining the decreasing effect from DNB and the increasing effect from accretion (from the ISM). Thus the change in y${}_{\text{light}}$ in a time interval $dt$ is

[TABLE]

where $\zeta_{\rm DNB}$ is calculated as in Equation 5 for a grid of $T_{\rm s}-T_{\rm b}-y_{\rm light}$ and interpolated between gridpoints. Note that our steady state approximation is not strictly valid at early, hot stages immediately after NS birth, when the NS cools faster than diffusion through the envelope. A calculation of neutron star interior temperature evolution and diffusion in the envelope is needed for a self-consistent solution, which is beyond the scope of this work.

5.3 ISM accretion rates

Whether the NS is accreting from the ISM and at what rate depends on the local density of the ISM and properties of the NS (e.g., its mass, radius and magnetosphere; see, e.g. Blaes & Madau 1993). Furthermore, when considering central compact objects (CCOs), the NS is within a supernova remnant, for which the distribution and composition differ from the ISM (see, e.g. Vink 2012; Lopez & Fesen 2018).

Following classical accretion physics (see, e.g. Shapiro & Teukolsky (1983)), the accretion rate onto the NS surface can be described using the Bondi formula

[TABLE]

where $\rho_{\rm ISM}$ is the local interstellar mass density, $v$ is relative speed between the NS and ISM, and $c_{s}$ is the speed of sound. The parameter $\lambda_{a}$ can vary between 0.25 and 1 (Hunt, 1971). Here we take $c_{s}$ = 10 km s*-1* and $\lambda_{a}$ = 1.

For simplicity, we consider three constant values of mass accretion rate for our models, i.e., we do not account for potential short time scale variability in the accretion rate. In the present work, we consider only the effect of a light element column density on the resulting cooling curve and ignore the effect of heat generated by nuclear reactions (see Section 7). We can obtain the highest accretion rate considered here (\dot{M}_{\rm acc}=10^{-15}\mbox{ M_{\rm Sun}$yr${}^{-1}}) by using Equation 8 for an isolated $1.4\,M_{\rm Sun}$ NS with v=20\mbox{ km s{}^{-1}} and ISM number density n=1\mbox{ cm{}^{-3}}. Lower accretion rate values ($10^{-18}$ and 10^{-21}\mbox{ M_{\rm Sun}$yr${}^{-1}}) are possible for higher velocities (e.g., 200 or 1000\mbox{ km s{}^{-1}}) and/or lower ISM densities. These low accretion rates are further motivated by possible suppression of the Bondi accretion rate such that NSs in this case accrete at lower rates (see, e.g., Illarionov & Sunyaev 1975; Perna et al. 2003; Krumholz et al. 2005).

5.4 Cooling curves

For H-He envelopes, we show the cooling curves and hydrogen column evolution in Figure 8 for an initial hydrogen column density of $y_{\rm H}$ = $10^{6}$ g cm*-2*; note that the column evolution curves here and in Figures 9 and 10 are smoothed (averaged over short time intervals) in order to remove short timescale fluctuations due to the finite time steps of the numerical simulations. The first conclusion that is clear from this plot is that the cooling curve is almost unaffected by changes in hydrogen column density, as we already found in Section 4. Nevertheless, it is interesting to see how the hydrogen column changes over time for varying ISM accretion rates. Note that for both the H-He and He-C models, we have artificially set the lower limit for the light element column at 100 g cm*-2*. Changes in the light element column below this limit will not affect the cooling curve, as they occur at low densities outside of the sensitivity strip. If DNB is the only mechanism for depleting surface hydrogen, the column density initially drops sharply due to high temperatures, until hydrogen reaction rates become limited by the amount of material available in the burning region. In other words, when the temperature has dropped significantly, the burning is limited by the time it takes for hydrogen to diffuse to the burning depth. Thus for most of such a NS’s life (when $T_{\rm s}\lesssim 10^{6}\mbox{ K}$), DNB in the NS envelope is not limited by the nuclear reaction rate, but by the diffusion time. For the two scenarios with higher accretion rates we find that, after a period of accreting and burning at similar rates ($\gtrsim$ 105 years), the accretion rate dominates, and the hydrogen column grows.

In Figure 9, we show cooling curves and helium column density evolution for an He-C envelope with an initial helium column density of $y_{\rm He}$ = $10^{10}$ g cm*-2* for different accretion rates. In all cases where $\dot{M}_{\rm acc}$ is less than the maximum accretion rate, more than 99$\%$ of the initial helium is consumed at the hot early stages and the helium column density drops to $y_{\rm He}$ $<$ 2 $\times$ 106 g cm*-2* within 1 year. During this time, helium is more rapidly burned than it is accreted. Depending on the accretion rate, the temperature has dropped enough at $t\sim 1-100$ years for the burning rate to become comparable to the accretion rate. After this time (when $\dot{y}_{\rm acc}>\zeta_{\rm DNB}$), the helium column gradually increases as matter accumulates from the ISM, and the temperature is too low for DNB to sufficiently counteract. The scenario with no accretion results in an effectively constant low column density after the initial helium column is depleted. The resulting cooling curve at times $>$ 1 year is the same as a cooling curve with no helium column.The cooling curves for the scenarios with larger accretion rates show short timescale variations, corresponding to variations in helium column. Although the short timescale changes in column density can be large, the effect on the surface temperature is small ($\sim$ a few percent), as the $T_{\rm s}-T_{\rm b}$ relations are only sensitive to changes in the column density that affect the composition in the sensitivity strip.

5.5 Implications for photospheric composition

In the discussion above, we focus on the effect of DNB on the composition of the envelope and the resulting temperature relations between the surface and the interior. We now consider the implications of the envelope evolution for the composition at the photosphere, which is uncertain for isolated NSs and is important for interpreting observations of NS surface radiation.

In Figures 8 and 9, the (artificially set) minimum light element column of $y_{\text{light}}$ = 100 g cm*-2* is large enough such that the photosphere, and corresponding spectrum, will always be dominated by hydrogen and helium, respectively. Due to this large lower limit, these figures do not directly address whether or not heavier elements can be present in the photosphere. Nevertheless Figure 8 shows that, even at the low accretion rate of 10${}^{-21}M_{\rm Sun}$ yr*-1*, the chosen minimum column is never reached, as the hydrogen column plateaus and increases after a few thousand years. Thus for a H-He envelope, the photosphere will be dominated by hydrogen even for very low accretion rates. On the other hand, Figure 9 shows that, for accretion rates $<10^{-15}M_{\rm Sun}$ yr*-1*, there will be a period after $\sim$1 month, where the spectrum of the He-C envelope can be dominated by carbon.

As accreted material from the ISM consists mainly of hydrogen, and thus accretion will mostly make the photosphere (and envelope) more hydrogen rich, it is interesting to consider the evolution of a H-C envelope. Therefore, we compute a H-C envelope model with $\rho_{b}$ = 1010 g cm*-3* and $g_{s}$ = 2.43 $\times$ 1014 cm s*-2*, for which the possible nuclear reactions are proton-proton captures and proton captures onto carbon (see Chang & Bildsten 2004). In Figure 10 we show the column evolution for an initial hydrogen column of $y_{\text{H}}$ = 106 g cm*-2* using a lower limit of 0.1 g cm*-2*. Note that this lower limit corresponds to a hydrogen fraction at the photosphere of $\sim$5$\%$ and thus a carbon dominated spectrum. Only slightly larger hydrogen columns are enough to result in a hydrogen spectrum, as a hydrogen column of $\sim 1.5$ g cm*-2* corresponds to a hydrogen fraction of $\sim$90$\%$ at the photosphere. Figure 10 shows that due to the highly efficient capture onto carbon, the initial hydrogen column is rapidly depleted within an hour for accretion rates $<10^{-15}M_{\rm Sun}$ yr*-1*. For these accretion rates, a period can exist when the spectrum is that of a carbon atmosphere. An accretion rate of $10^{-20}M_{\rm Sun}$ yr*-1* leads to an optically thick hydrogen atmosphere after $\sim$1000 year. For higher accretion rates, the hydrogen column will effectively always be large enough to maintain a hydrogen atmosphere. On the other hand, accretion rates $<10^{-21}M_{\rm Sun}$ yr*-1* are too low to build up an optically thick hydrogen column and will result in a carbon atmosphere after the initial hydrogen column is consumed. These results are relevant for the understanding of the neutron star in Cas A, of which the spectra can be well described by a carbon atmosphere (Ho & Heinke, 2009), as well as other NSs for which a carbon atmosphere is suggested (see Section 7). We will consider temperature relations for H-C envelopes and their application to these neutron stars in future work.

6 Application to Cassiopeia A

6.1 Chandra data analysis

We process all Chandra ACIS-S GRADED data using Chandra Interactive Analysis of Observations (CIAO) 4.10 and Chandra Calibration Database (CALDB) 4.8.1, following the procedure described in Heinke & Ho (2010). Spectra are extracted using specextract from a 4 pixel circular region centered on the source, while background is taken from a 5$-$8 pixel annulus centered on the source. Observations taken close in time (see Table 1) are merged using combine_spectra and dmgroup with a minimum of 200 counts per energy bin.

Spectra are fit with Xspec 12.10.0c (Arnaud, 1996) and a four component model consisting of pileup, tbabs, spexpcut, and nsx. For pileup (Davis, 2001), we set the grade migration parameter to 0.27 for observations through 2004, which have a frame time of 3.24 s, and to 0.24 for observations starting from 2007, which have a frame time of 3.04 s; the maximum number of photons is set to 5, and the point-spread-function fraction is set to 0.95 (see Heinke & Ho 2010, for details). To model photoelectric absorption, we use tbabs with abundances from Wilms et al. (2000). To model interstellar dust scattering, we use spexpcut with exponent index $\alpha=-2$ and characteristic energy E_{\rm cut}(\mbox{keV})=[0.49N_{\rm H}(10^{22}\mbox{ cm{}^{-2}})]^{1/2} (Predehl et al., 2003), where $N_{\rm H}$ is the absorption column determined from tbabs. For nsx (Ho & Heinke, 2009), we assume a partially ionized carbon atmosphere and fix the distance to 3.4 kpc (Reed et al., 1995) and normalization to 1.

6.2 Spectral fits

To determine the neutron star mass $M$ and radius $R$ for use in spectral fitting, we first allow fit parameters $M$, $R$, $N_{\rm H}$, and $T_{\rm s}$ to vary but tie them between observations, so that each has the same value across all observations. The resulting fit produces best-fit values of $M=1.44\pm 0.15\,M_{\rm Sun}$, $R=13.78\pm 0.36\mbox{ km}$, N_{\rm H}=(1.623\pm 0.030)\times 10^{22}\mbox{ cm{}^{-2}}, and $T_{\rm s}=(1.728\pm 0.030)\times 10^{6}\mbox{ K}$, with $\chi^{2}=1646$ for 1458 degrees of freedom (dof). If we instead allow $T_{\rm s}$ to have different values between observations, then the best-fit values are $M=1.65\pm 0.16\,M_{\rm Sun}$, $R=12.94\pm 0.34\mbox{ km}$, and N_{\rm H}=(1.674\pm 0.031)\times 10^{22}\mbox{ cm{}^{-2}}, with $\chi^{2}/\mbox{dof}=1522/1446$. If we allow both $N_{\rm H}$ and $T_{\rm s}$ to have different values between observations, then the best-fit values are $M=1.46\pm 0.16\,M_{\rm Sun}$ and $R=13.69\pm 0.36\mbox{ km}$, with $\chi^{2}/\mbox{dof}=1497/1434$. A f-test between the first and second fits ($\Delta\chi^{2}=124$) and between the second and third fits ($\Delta\chi^{2}=25$) yield probabilities of $9.8\times 10^{-19}$ and 0.024, respectively, which suggest that allowing $T_{\rm s}$, and possibly $N_{\rm H}$, to change between observations is warranted. Therefore in the spectral fits performed below, we untie $T_{\rm s}$ and fix the neutron star mass and radius to $M=1.65\,M_{\rm Sun}$ and $R=12.9\mbox{ km}$, respectively.

We show the $M$-$R$ confidence contour in Fig. 11, assuming a distance of 3.4 kpc. The best-fit values from our previous work (Heinke & Ho, 2010) and used in subsequent works (Shternin et al., 2011; Elshamouty et al., 2013; Ho et al., 2015) are $M=1.65\,M_{\rm Sun}$ and $R=10.3\mbox{ km}$. It is important to note that, while the $M$-$R$ constraints we obtain are somewhat constraining, there are caveats. First, there are calibration issues with the data, which are discussed in Section 7 (see also Posselt et al. 2013). Second, the magnetic field strength of Cas A is assumed here to be low enough as to not affect its atmospheric emission. This may not necessarily be the case, as three other CCOs have surface magnetic fields $\sim 10^{10}-10^{11}\mbox{ G}$ (Gotthelf et al., 2013; Ho, 2013; De Luca, 2017), which is somewhat larger than accommodated for in our atmosphere model (c.f. Potekhin et al. 2014). Finally, it is possible that the NS surface is not uniform in temperature, and hot spots could dominate the observed emission (see, e.g., Gotthelf et al. 2010; Bogdanov 2014), although pulsations from this emission are not detected in Cas A (Murray et al., 2002; Halpern & Gotthelf, 2010).

Next, we perform two fits of all spectra with the same set of models but fixing the parameters $M=1.65\,M_{\rm Sun}$ and $R=12.9\mbox{ km}$. The first fit allows $T_{\rm s}$ to vary between observations but holds the absorption parameter fixed at N_{\rm H}=1.67\times 10^{22}\mbox{ cm{}^{-2}}. This fit yields $\chi^{2}/\mbox{dof}=1522/1449$, and the temperature evolution is shown in Fig. 12 and Table 1. The second fit allows both $T_{\rm s}$ and $N_{\rm H}$ to change between observations. This fit yields $\chi^{2}/\mbox{dof}=1498/1436$, and temperature and absorption column evolutions are shown in Figs. 13 and 14, respectively (see also Table 1). A f-test between these two fits ($\Delta\chi^{2}=24$) yields a probability of 0.043, so allowing $N_{\rm H}$ to vary is possibly justified (see Section 7). For comparison to our previously published results, we also perform a fit using fixed (M,R,N_{\rm H})=(1.65\,M_{\rm Sun},10.3\mbox{ km},1.73\times 10^{22}\mbox{ cm{}^{-2}}). Results are shown in Fig. 12 and Table 1 and have $\chi^{2}/\mbox{dof}=1562/1449$. The approximately constant offset in temperatures (e.g., $\sim 2.13/1.87=1.14$) between the two fixed-$N_{\rm H}$ fit results plotted in Fig. 12 is primarily due to the difference in assumed $R$ and its anti-correlation with $T_{\rm s}$, i.e., for constant luminosity, $T_{\rm s}\propto R^{-1/2}$, which yields $(12.9\mbox{ km}/10.3\mbox{ km})^{1/2}=1.12$, as well as a difference in gravitational redshift.

Finally, we perform a linear fit to each set of declining temperatures $T_{\rm s}$. For (M,R,N_{\rm H})=(1.65\,M_{\rm Sun},10.3\mbox{ km},1.73\times 10^{22}\mbox{ cm{}^{-2}}), the slope is -13\pm 1\mbox{ K d{}^{-1}} ($\approx-2.2\pm 0.2\%$ per 10 years) with $\chi^{2}/\mbox{dof}=14.8/11$. For (M,R,N_{\rm H})=(1.65\,M_{\rm Sun},12.9\mbox{ km},1.67\times 10^{22}\mbox{ cm{}^{-2}}), the slope is -11\pm 1\mbox{ K d{}^{-1}} ($\approx-2.1\pm 0.2\%$ per 10 years) with $\chi^{2}/\mbox{dof}=15.7/11$. For $(M,R)=(1.65\,M_{\rm Sun},12.9\mbox{ km})$ and changing $N_{\rm H}$, the slope is -14\pm 1\mbox{ K d{}^{-1}} ($\approx-2.7\pm 0.3\%$ per 10 years) with $\chi^{2}/\mbox{dof}=6.2/11$. We discuss uncertainties beyond statistical ones in our measured temperature decline rates and compare to rates measured in previous works in Section 7.

6.3 Cas A cooling curve

The observed temperature drop has been interpreted as a result of neutron superfluidity and proton superconductivity in Cas A (Page et al. 2011; Shternin et al. 2011; Yakovlev et al. 2011; Ho et al. 2015, see Section 1 for references on alternative interpretations). With the addition of 4 new data points presented above, we now have a set of 13 temperatures spanning over 18 years. We apply the cooling model described in Section 5.1 to the full set of data. We use the temperature data corresponding to the spectral fits with (M,R,N_{\rm H})=(1.65\,M_{\rm Sun},12.9\mbox{ km},1.67\times 10^{22}\mbox{ cm{}^{-2}}). The cooling model makes use of a time-variable He-C envelope model with an initial He column of 106 g cm*-2* and assumes no accretion from the ISM. As we found in Section 5, the cooling curve is not sensitive to the initial helium column density, as practically all helium will be depleted within 10 years. Therefore, the cooling curve at the time of the Cas A observations (at age $\approx 320-340\mbox{ yr}$) is the same as a cooling curve using a pure carbon envelope.

In order to describe the rapid temperature decline, we follow the previous works (e.g., Page et al. 2011; Shternin et al. 2011) by invoking superfluidity of neutrons in the inner crust through pairing in the ${}^{1}S_{0}$ state and in the core through pairing in the ${}^{3}P_{2}$ state. Within the minimal cooling paradigm (Gusakov et al., 2004; Page et al., 2004), this means that stronger neutrino emission is only allowed as a result of neutron or proton pairing. The exact properties, such as the critical temperature $T_{\rm c}$, of neutron ${}^{3}P_{2}$ superfluidity in the core are uncertain and can be constrained by comparing the cooling curve to the observed temperatures (Page et al., 2011; Shternin et al., 2011; Ho et al., 2015). In our model here, we use superfluidity and superconductivity properties similar to those in Page et al. (2011); specifically, we use a neutron ${}^{3}P_{2}$ superfluid model with a maximum critical temperature $T_{\rm c}\sim 5.9\times 10^{8}\mbox{ K}$ for $M=1.7\,M_{\rm Sun}$ (with the APR EOS, resulting in $R$ = 11.4 km) and the ‘CCDK’ proton ${}^{1}S_{0}$ superconducting model described in Page et al. (2004). Note that our considered theoretical mass and radius using the APR EOS (1.7 Msun and 11.4 km) is within the 90 percent confidence contour of mass and radius obtained from the spectral analysis (see Figure 11). The aim of our calculations here is not to find a best fitting model but to show that models fitting the new data are in line with previous results for superfluid properties and can also be explained using the envelope $T_{\rm b}-T_{\rm s}$ relations that include DNB.

To compare the model to the observations, we use the date for the birth of the Cas A NS of $1681\pm 19$ (Fesen et al., 2006). The resulting cooling curve is shown in Figure 15 and can adequately describe the observed temperatures. It is important to note that several combinations of superfluid and neutron star properties can describe the observed cooling, so the model presented here is not a unique solution and is merely for illustrative purposes.

6.4 Implications of DNB on Cas A photosphere

In Section 5.5, we show that DNB rapidly consumes any initial hydrogen from H-C envelopes for low to moderate accretion rates, thus allowing for an atmosphere composed of carbon to be visible. Effectively, hydrogen atmosphere formation by accretion in the dense post-supernova environment is delayed by DNB to a later, lower density epoch. For the relatively young Cas A NS, this is likely the case. In order for Cas A to have a carbon atmosphere, it must accrete at an average rate lower than 10^{-20}\mbox{ M_{\rm Sun}$yr${}^{-1}}; for larger accretion rates, DNB cannot consume hydrogen fast enough to prevent formation of a hydrogen atmosphere. While DNB can explain a carbon atmosphere, it is unlikely to be responsible for the observed temperature variations in Cas A, as the light element burning rate is negligible when the atmosphere is dominated by carbon.

7 Discussion

In this present work, we consider neutron star envelopes consisting of either H-He or He-C and including the physics of diffusive nuclear burning, in order to investigate how DNB impacts the thermal evolution of NSs. DNB is relevant for the thermal evolution in both a static and time-dependent manner: static because DNB affects the interior composition profile of the envelope and thus can alter the $T_{\rm s}-T_{\rm b}$ relations, and time-dependent because it can alter the composition of the envelope over time. In Section 3, we calculate the burning rates and find that DNB can alter the composition on short time-scales ($\lesssim 10-100$ yr). For H-He, the hydrogen burning rate is strongly coupled to the amount of hydrogen available (characterised by its column density) whereas for the He-C envelope the burning rate is highly sensitive to the temperature and is only significant for $T_{\rm s}>10^{6}$ K.

In Section 4, we obtained static analytic relations for the surface temperature $T_{\rm s}$ and the temperature at the bottom of the envelope $T_{\rm b}$. These temperature relations (in Appendix A) can be used to model the thermal evolution of isolated and accreting neutron stars. We find the following main differences in the $T_{\rm s}-T_{\rm b}$ relations due to DNB:

• •

For a H-He envelope, we find that when DNB is taken into account, there is a clear upper limit on the size of the hydrogen column. Larger hydrogen columns ($y_{H}\gtrsim 10^{7}$ g cm*-2*) are not thermally stable, as hydrogen is consumed rapidly when DNB is nuclear-limited. For temperatures above $T_{\rm s}\sim$ $10^{6}$ K, this means the spread in $T_{\rm b}$ is smaller than without DNB (see Figure 4).

• •

For He-C envelopes, DNB is highly sensitive to the temperature, and DNB makes the envelope more transparent to the heat flux at $T_{\rm s}>2\times 10^{{}^{6}}$ K for the helium column sizes that affect the sensitivity strip (see Figure 5).

We then combine these results in Section 5 and use the DNB burning rates (Section 3) and the temperature relations (Section 4) to calculate NS cooling curves with envelope compositions that evolve over time through nuclear burning and accretion. As the H-He temperature relations only span a small range of $T_{\rm b}$, the changing envelope composition does not impact the cooling curves significantly. We find that if DNB is the only burning mechanism and no accretion takes place, an initial hydrogen column of $10^{6}$ g cm*-2* is depleted after $\sim 5\times 10^{3}$ years because the burning rate drops as the hydrogen column decreases. For the He-C envelope all helium will be depleted during the hot early cooling phase (within 1 year) if no accretion takes place. Note that we have only considered illustrative cases and have not included heat generated by the nuclear reactions nor the possibility of sudden accretion events.

For several CCOs, the spectra could be well described by a carbon atmosphere (see, e.g., Klochkov et al. 2016; Suleimanov et al. 2017; Doroshenko et al. 2018). We consider the photosphere of a H-C envelope to investigate the possibility of a carbon atmosphere on isolated neutron stars. Low levels of accretion lead to an optically thick hydrogen atmosphere, with a hydrogen column of $1.5$ g cm*-2* corresponding to a hydrogen abundance in the atmosphere of 90$\%$. The composition evolution of the H-C envelope in Figure 10 for accretion rates $\leq$ $10^{-20}M_{\rm Sun}$ yr*-1* indicate that it is possible for Cas A, at an age of 340 years, to not yet accrete enough hydrogen to form an optically thick hydrogen atmosphere. For older NSs, the possibility of a carbon atmosphere becomes less likely, as even a little bit of accretion leads to a hydrogen dominated atmosphere.

Note that we have only considered constant accretion rates here, whereas scenarios with changing accretion rates are more realistic for CCOs, as their environment is highly variable (see, e.g., Vink 2012; Lopez & Fesen 2018). For example, it is possible that the NS initially accretes at very low rates after the nearby environment has been cleared out by the SN shock. Later on, a reverse shock could introduce a higher density of the surrounding medium, and thus a higher accretion rate onto the NS. Early fallback of material could also lead to varying accretion rates. As the variability in the accretion rate and the composition of the accreted material in SNRs are highly model dependent, we have limited our models to illustrative cases with constant rates and compositions.

In all the envelope models currently used in cooling simulations, it is assumed that the flux generated by nuclear reactions in the envelope is negligible compared to the radiative flux from the neutron star interior. However, we find that the assumption of constant flux through the envelope can be invalid in relevant regions of parameter space. The assumption of negligible heat production in the envelope is certainly incorrect for the hot, initial post-supernova stage, when surface temperatures are $\gg 10^{6}\mbox{ K}$. However, as the energy is produced at low densities ($\rho<10^{8}-10^{10}$ g cm*-3*) and the interior is still hot, this heat is radiated from the surface and will not affect cooling at later times. Including nuclear heating at early times will likely lead to an even more rapid consumption of light elements than calculated here. At later times, heat generated by DNB is negligible for H-He envelopes, while heat created by DNB can be important for high temperatures (T $\gtrsim$ 2 $\times$ 106 K) and large helium columns (y${}_{\text{He}}\gtrsim 10^{11}$ g cm*-2*) in the case of He-C envelopes. Interestingly, for H-C envelopes the luminosity generated by DNB can be larger than the interior luminosity for hydrogen columns greater than y${}_{\text{H}}$ $>10^{6}$ g cm*-2* and surface temperatures higher than $T_{\rm s}\gtrsim 10^{6}$ K. We will examine the effect of heat generated by nuclear reactions in future work.

In Section 6, we present new observations of the NS in the Cassiopeia A supernova remnant and apply our model to this data. We performed spectral fits of 13 Chandra ACIS-S GRADED observations of Cas A taken over more than 18 years (see Table 1). For constant interstellar absorption (N_{\rm H}\approx 1.7\times 10^{22}\mbox{ cm{}^{-2}}), these data indicate that the NS temperature is decreasing at a ten-year rate of $2.1\pm 0.2\%$ ($1\sigma$) and a decreasing ten-year flux rate of $\approx 7\%$. However, it is known that there is a contaminant building up on the Chandra ACIS-S detector that can affect GRADED mode data if uncorrected (Posselt et al., 2013; Plucinsky et al., 2016, 2018). Elshamouty et al. (2013) analyse data from different detectors and operating modes and infer a ten-year temperature decline rate of $3.5\pm 0.4\%$ for ACIS-S GRADED data and a ten-year decline rate of $1.0\pm 0.7\%$ for HRC-S data, suggesting that at least one dataset is affected by systematic errors. Posselt et al. (2013) analyse an alternative set of 2 Chandra ACIS-S subarray observations taken over 6 years in FAINT mode, to minimize the effects of pile-up and contaminant; they find a lower ten-year rate of flux decline and no apparent temperature decline ($10\pm 5\%$ and $1.3\pm 1.0\%$, respectively, where errors here are at 90% confidence level). An important result of our work here is that the temperature decline rate ($2.1\pm 0.2\%$ for constant $N_{\rm H}$ and $2.7\pm 0.3\%$ for varying $N_{\rm H}$) is significantly lower than values previously reported based on a smaller set of ACIS-S GRADED observations ($3.9\pm 0.7\%$ per 10 years from Heinke & Ho 2010 and $3.5\pm 0.4\%$ per 10 years from Elshamouty et al. 2013). Therefore, the rates of temperature decline measured using HRC-S data and ACIS-S subarray and GRADED data all agree at $\sim 90\%$ confidence level.

Very recently, Posselt & Pavlov (2018) follow up their previous analysis by including a new ACIS-S subarray observation, which expands this dataset to 3 measurements in 9 years, and use CALDB 4.7.3444CALDB 4.7.3 includes contaminant model N0010 for ACIS-S. Similarly, CALDB 4.7.8 (used in Appendix B) includes N0010. On the other hand, the most recent calibration CALDB 4.8.1 (used in Section 6) includes the latest contaminant model N0012 for ACIS-S (see http://cxc.cfa.harvard.edu/caldb4/downloads/Release_notes/CALDB_v4.8.1.html). (also testing versions up to 4.7.7); their analysis yields 3$\sigma$ upper limits on the ten-year rate of temperature decline of 2.4% assuming a constant $N_{\rm H}$ and 3.3% assuming a varying $N_{\rm H}$. In our new analysis using 13 measurements in 18 years, as well as the latest Chandra calibration and contamination models (CALDB 4.8.1; see footnote 4 and Appendix B), the temperature decline we measured using ACIS-S GRADED mode is in agreement with the upper limits found by Posselt & Pavlov (2018) using the subarray observations.

As noted above, our spectral fits suggest $N_{\rm H}$ may be variable (see Table 1 and Fig. 14). Alp et al. (2018) study the effect on $N_{\rm H}$ of material expelled by a supernova using three-dimensional supernova simulations. They find that different supernova models do not manifest measurable differences in $N_{\rm H}$ at the current age of the Cassiopeia A neutron star. They also analyse the same three Chandra observations (with $N_{\rm H}$ tied across observations) as those used by Posselt & Pavlov (2018) and measure a ten-year rate of temperature decline of $\sim$1$\pm$1$\%$ and a supernova ejecta contribution to $N_{\rm H}$ that is indistinguishable from that due to the intervening interstellar medium. While motion of intervening supernova remnant material (or the NS) as the cause of varying $N_{\rm H}$ cannot be ruled out, a varying $N_{\rm H}$ likely reflects incomplete modelling of the detector contaminant (see, e.g., Plucinsky et al. 2018). We note however that in the case of our spectral fits allowing for a varying $N_{\rm H}$, the ten-year rates of temperature and flux decline are actually larger, i.e., $2.7\pm 0.3\%$ and $\approx 8\%$, respectively (see Fig. 13), than the case when $N_{\rm H}$ is constant. Recent analysis by Plucinsky et al. (2018) of the ACIS contaminant suggests the build up on the detector is not increasing as much as predicted by model N0010, and this effect is taken into account in the revised model N0012 (see footnote 4). Future refinements of the ACIS-S calibration may aid attempts to more accurately measure the cooling rate of the NS in Cassiopeia A.

Acknowledgements

COH thanks B. Posselt for discussions. The authors are grateful to the anonymous referee for comments which improved the manuscript. WCGH acknowledges support through grant ST$/$R00045X$/$1 from Science and Technology Facilities Council in the UK. COH is supported by NSERC Discovery Grant RGPIN-2016-04602 and a Discovery Accelerator Supplement. MB is fully and DP partially supported by the Consejo Nacional de Ciencia y Tecnología with a CB-2014-1 grant $\#$240512. This research made use of NumPy (Van Der Walt et al., 2011), SciPy (Jones et al., 2001) and matplotlib, a Python library for publication quality graphics (Hunter, 2007) .

Appendix A Analytic Ts-Tb relations

In this section, we present accurate fits to the computed $T_{\rm s}$-$T_{\rm b}$-$\rho^{*}$ data which can be used in cooling simulations. All the fits are obtained for an envelope with surface gravity $g_{s,0}$ = 2.4271 $\times$ 1014 cm s*-2* but can be scaled for any $g_{s}$ using $Y=(T_{\rm s}/1{\rm MK})(g_{s,0}/g_{s})^{1/4}$ (Gudmundsson et al., 1983). Firstly, we obtain $T_{\rm b}$($T_{\rm s}$, $\rho^{*}$) fits and use an adapted version of the shape of the analytic functions (Equation 9) presented by Beznogov et al. (2016).

[TABLE]

A.1 H-He mixture

The analytic results for the H-He mixture are applicable to an envelope with a lower boundary of $\rho_{b}$ = 108 g cm*-3*. While the H-He temperature relations can be fit with a simple function needing less fit parameters, we fit the same function as Beznogov et al. (2016) for consistency. The best fit parameters correspond to a maximum relative error of 0.0054 and a rms of the relative error of 0.0028.

[TABLE]

A.2 He-C mixture

The analytic results for the He-C mixture for an envelope with a lower boundary of $\rho_{b}$ = 1010 g cm*-3*. Equation 9 consists of the functions given in 11. Here, the functions $f_{1}-f_{4}$ are identical to those in Beznogov et al. (2016). The function $f_{5}$ is an additional term to accurately fit the DNB models. For the He-C mixture, we obtain a maximum relative error of 0.019 which is slightly larger than the maximum relative error of 0.017 that Beznogov et al. (2016) obtain for the C-Fe envelope. The root mean square of the relative error is 0.0037.

[TABLE]

Appendix B Analysis of Cassiopeia A data using previous Chandra contamination model N0010

Here we show results of our analysis of the same data as that in Section 6 but using a previous CALDB (4.7.8 in CIAO 4.9), with ACIS-S contamination model N0010, and Xspec (12.9.1). The results presented in this Appendix are to illustrate the small changes that the ACIS-S contamination model N0012 (in CALDB 4.8.1) introduces, as well as to facilitate comparisons with the results of Posselt & Pavlov (2018) who use CALDB 4.7.3 and N0010. The analysis follows the same procedure as that described in Section 6, and the final fit results are given in Table 4. Comparison with the fit results in Table 1 show differences in temperature and flux well within uncertainties, keeping in mind that some of the small changes are due to differing best-fit radius (13.0 km versus 12.9 km) and hence gravitational redshift.