Crust cooling of the neutron star in Aql X-1: Different depth and magnitude of shallow heating during similar accretion outbursts

TL;DR

This study investigates the variability in shallow crust heating of the neutron star in Aql X-1 during similar accretion outbursts, revealing that the depth and magnitude of heating can differ significantly, affecting cooling behavior.

Contribution

The paper provides observational evidence that shallow heating varies in depth and strength between outbursts, challenging assumptions of uniform crustal heating mechanisms in neutron stars.

Findings

Shallow heating was stronger and deeper during the 2016 outburst compared to 2013.

Neutron star parameters remain constant, indicating other factors influence shallow heating.

Different outbursts with similar accretion profiles can produce diverse crust heating effects.

Abstract

The structure and composition of the crust of neutron stars plays an important role in their thermal and magnetic evolution, hence in setting their observational properties. One way to study the crust properties is to measure how it cools after it has been heated during an accretion outburst in a low-mass X-ray binary (LMXB). Such studies have shown that there is a tantalizing source of heat, of currently unknown origin, that is located in the outer layers of the crust and has a strength that varies between different sources and different outbursts. With the aim of understanding the mechanism behind this "shallow heating", we present Chandra and Swift observations of the neutron star LMXB Aql X-1, obtained after its bright 2016 outburst. We find that the neutron star temperature was initially much lower, and started to decrease at much later time, than observed after the 2013 outburst of the source, despite the fact that the properties of the two outbursts were very similar. Comparing our data to thermal evolution simulations, we infer that the depth and magnitude of shallow heating must have been much larger during the 2016 outburst than during the 2013 one. This implies that basic neutron star parameters that do not change between outbursts, do not play a strong role in shallow heating. Furthermore, it suggests that outbursts with a similar accretion morphology can give rise to very different shallow heating. We also discuss alternative explanations for the difference in quiescent evolution after the 2016 outburst.

Tables

Table 1. Table 1: Spectral analysis results for 2016–2017 quiescence observations of Aql X-1.

Instrument	Epoch	mean MJD	ObsID(s)	$k{T}^{\infty }$	$F_{\mathrm{th}}$	$F_X$	$L_X$
				(eV)	($\times {10}^{-12}\mathrm{erg}{\mathrm{cm}}^{-2}{s}^{-1})$		($\times {10}^{33}\mathrm{erg}{s}^{-1})$
Swift/XRT	2016 epoch 1	57677.7	34719016–20	$118.9\pm 4.5$	$1.15\pm 0.20$	$1.46\pm 0.30$	$4.4\pm 0.9$
Swift/XRT	2016 epoch 2	57692.9	34719021–26	$116.2\pm 2.8$	$0.97\pm 0.14$	$1.17\pm 0.19$	$3.5\pm 0.6$
Chandra/ACIS	2016 Nov 22	57714.8	18984	$117.5\pm 1.3$	$1.10\pm 0.13$	$1.50\pm 0.11$	$4.5\pm 0.3$
Chandra/ACIS	2017 Feb 20	57804.7	18985	$110.7\pm 1.1$	$0.81\pm 0.09$	$0.98\pm 0.06$	$2.9\pm 0.2$
Swift/XRT	2017 epoch 3	57859.1	33665090–99	$102.6\pm 2.5$	$0.62\pm 0.38$	$0.82\pm 0.08$	$2.5\pm 0.3$

Table 2. Table 2: Results of the crust cooling simulations.

Outburst	log $Y_C$	$Q_{sh}$	${\rho }_{\mathrm{sh},\min }$
		(MeV nucleon-1)	($g{\mathrm{cm}}^{-3}$)
2016	${6.6}_{\ast }^{+0.9}$	$9.2\pm 1.6$	${2.8}_{-0.2}^{+0.1}\times {10}^{10}$
2016 - cold core	${10.5}_{-0.8}^{+0.3}$	${5.3}_{-0.7}^{+5.4}$	${3.4}_{-0.6}^{+1.2}\times {10}^{10}$
2013	${8.8}_{-1.5}^{+1.1}$	${2.3}_{-0.3}^{+0.5}$	${0.4}_{\ast }^{+0.7}\times {10}^9$
2011	${8.3}_{-0.9}^{+0.7}$	${3.7}_{-0.9}^{+1.5}$	${0.4}_{\ast }^{+7.9}\times {10}^9$

Full text

††pagerange: Crust cooling of the neutron star in Aql X-1: Different depth and magnitude of shallow heating during similar accretion outbursts–LABEL:lastpage††pubyear: 0000

1 Introduction

The liquid, dense core of neutron stars is covered by a $\sim$1-km thick, solid crust. The structure and composition of the crust of neutron stars play a key role in the evolution of their magnetic field strength and interior temperature. As such, the crust properties are important for a variety of observational phenomena such as, for instance, pulsar glitches, thermonuclear X-ray bursts, and magnetar outbursts, as well as gravitational wave signals from neutron stars (e.g. Brown et al., 1998; Ushomirsky et al., 2000; Cumming & Bildsten, 2001; Horowitz & Berry, 2009; Pons et al., 2009; Page & Reddy, 2013). This provides a strong incentive to understand the detailed properties of the crust of neutron stars.

Neutron stars that are part of a binary system in which they can accrete gas from a companion star with a mass $\lesssim 1~{}\mathrm{M}_{\odot}$, are called low-mass X-ray binaries (LMXBs). Many neutron star LMXBs are transient and spend most of their time in quiescence. The gas that is transferred from the companion is then only accreted at high rates during weeks–years long outbursts. During such outbursts, LMXBs can exhibit different “spectral states” that are characterized by specific X-ray spectral and fast-variability properties (Hasinger & van der Klis, 1989), and likely reflect different accretion geometries (e.g. Done et al., 2007, for a review).

The gas that accretes on to the surface of neutron stars undergoes thermonuclear burning, which transfers light elements (i.e. the accreted H and/or He) into heavier ones (e.g. Schatz et al., 1999). This thermonuclear burning is often unstable and causes run-away energy production (e.g. Wallace & Woosley, 1981) that results in a thermonuclear X-ray burst (X-ray burst hereafter; e.g. Lewin et al., 1995; Galloway & Keek, 2017, for reviews).

Apart from causing thermonuclear burning on the surface, accretion fires up a series of nuclear reactions in the crust of neutron stars. This includes electron captures by nuclei in the outer layers of the crust, and density-driven fusion of nuclei at several hundreds meter depth. Based on theoretical calculations (e.g. Haensel & Zdunik, 1990; Yakovlev et al., 2006; Fantina et al., 2018), and laboratory data from nuclear experiments (e.g. Gupta et al., 2007; Estradé et al., 2011), the energy produced in these nuclear reactions is thought to be $\sim$2 MeV per accreted nucleon (e.g. Haensel & Zdunik, 2008). This energy release can significantly heat the crust and bring it out of thermal equilibrium with the core. Since most energy is generated in the nuclear fusion reactions, this process is referred to as “deep crustal heating” (Brown et al., 1998). The crust cools during quiescent phases, when the heat energy gained during outburst is thermally conducted throughout the neutron star (e.g. Colpi et al., 2001; Rutledge et al., 2002b; Brown & Cumming, 2009; Page & Reddy, 2012; Wijnands et al., 2013).

Over the past decade, evidence has been accumulating that during accretion phases in LMXBs, the crust of neutron stars are more strongly heated than predicted by nuclear heating models. This inference comes from observations of different types of phenomena. For instance, observations of very long X-ray bursts and mHz quasi-periodic oscillations (QPOs) point to higher crust temperatures than can be accounted for with standard nuclear heating (e.g., Cumming et al., 2006; Keek et al., 2008; Altamirano et al., 2012; in ’t Zand et al., 2012; Linares et al., 2012). Furthermore, detailed monitoring of the quiescent temperature evolution of neutron stars following accretion outbursts, has highlighted that there is something missing in our understanding of how neutron star crusts are heated (e.g. Brown & Cumming, 2009; Degenaar et al., 2011).

Ten neutron stars have been monitored after accretion outbursts with the aim to study how their accretion-heated crusts cool in quiescence (see Wijnands et al., 2017, for a review). The high temperatures observed within the first few hundred days of quiescence require an additional source of energy to heat the crust (e.g. Brown & Cumming, 2009; Degenaar et al., 2011). The strength of this “shallow heating”, if proportional to the mass-accretion rate, is inferred to be on the order of $\sim$1–2 MeV nucleon*-1* (e.g. Degenaar et al., 2011, 2014, 2015; Page & Reddy, 2013; Parikh et al., 2017b, 2019; Ootes et al., 2018), although one extreme case requires $\sim$15–17 MeV nucleon*-1* (e.g. Homan et al., 2014; Deibel et al., 2015; Parikh et al., 2017a). It is currently unclear what is causing this shallow heating.

It remains to be established to what extent shallow heating depends on neutron-star specific properties (e.g., spin, magnetic field strength, mass, radius, superfluid properties, age), and on the detailed properties of the accretion outburst (e.g., brightness, duration, accretion geometry). Studying multiple cooling curves of a single source may be a promising way to understand the origin of shallow heating; because the fundamental properties of a neutron star remain virtually unchanged, the effect of the outburst parameters on the crust heating and cooling can potentially be isolated (Waterhouse et al., 2016; Parikh et al., 2017a, 2019).

1.1 The frequently active neutron star LMXB Aql X-1

Aql X-1 is a transient neutron star LMXB that was discovered in the late sixties (Friedman et al., 1967) and has been seen active many times since (e.g. Kaluzienski et al., 1977; Kitamoto et al., 1993; Campana et al., 2013; Güngör et al., 2014; Ootes et al., 2018). Its accretion outbursts last $\sim$1–6 months, vary in brightness from L_{\mathrm{X}}$$\simeq$$10^{35} to $10^{38}~{}(D/5.0~{}\mathrm{kpc})^{2}~{}\mathrm{erg~{}s}^{-1}$, and recur on a timescale of $\sim$1 yr. The neutron star displays X-ray bursts (e.g. Koyama et al., 1981), and spins at a frequency of $\simeq$550 Hz (Zhang et al., 1998; Casella et al., 2008). The source is located at $\approx 5$ kpc (e.g. Rutledge et al., 2001; Galloway et al., 2008)111We note that Aql X-1 appears in Gaia DR2, but due to its faintness the parallax error is large and therefore the inferred distance estimate becomes strongly dependent on the assumed prior., and the binary companion is a K-type star that orbits the neutron star in $\approx$19 hr (Callanan et al., 1999; Mata Sánchez et al., 2017).

With its frequent outbursts of varying length and brightness, Aql X-1 could be a promising source to gain more insight into shallow crustal heating (Waterhouse et al., 2016). However, the origin of its quiescent X-ray emission has been debated. Firstly, apart from soft thermal emission from the neutron star surface, its quiescent X-ray spectrum often contains a hard emission tail. This component can be modeled as a $\Gamma\sim 1-2$ power law, and can contribute up to $\sim$80% of the total unabsorbed 0.5–10 keV flux (e.g. Verbunt et al., 1994; Rutledge et al., 2002a; Campana & Stella, 2003; Cackett et al., 2011; Campana et al., 2014). Such a power-law component is often seen in the spectra of quiescent LMXBs and could point to ongoing accretion (see e.g. Chakrabarty et al., 2014b; D’Angelo et al., 2015; Wijnands et al., 2015, for recent discussions). Secondly, days-long flares have been seen from Aql X-1 during which the quiescent X-ray emission increased by a factor of $\gtrsim$10 (Coti Zelati et al., 2014). It is commonly assumed that such flares are caused by a (short-lived) spurt of accretion (e.g. Degenaar & Wijnands, 2009; Fridriksson et al., 2011; Wijnands & Degenaar, 2013).

Despite indications that there is at least some level of quiescent accretion activity in Aql X-1, Waterhouse et al. (2016) argued that the crust of the neutron star may be so hot that its cooling drives the long-term quiescent flux evolution (because the injection of heat from continued low-level accretion is then very small compared to the heat content of the crust; see also e.g. Fridriksson et al., 2011; Turlione et al., 2015; Homan et al., 2014; Bahramian et al., 2015; Deibel et al., 2015; Parikh et al., 2017b). Waterhouse et al. (2016) studied X-ray data from the Neil Gehrels Swift observatory (Swift; Gehrels et al., 2004) after three different outbursts of Aql X-1 (2011, 2013 and 2015). The quiescent spectra were dominated by a thermal emission component and the inferred temperatures were compared to a neutron star thermal evolution code (DStar; Brown & Cumming, 2009). It was shown that the measured temperatures and decay trends could naturally be explained within the crust cooling paradigm (Waterhouse et al., 2016).

The possibility of observing crust cooling in Aql X-1 was further investigated by Ootes et al. (2018), who tracked the evolution of the neutron star temperature over its 1996–2015 outburst history, using the thermal evolution code NSCool (Page, 2016), and allowing it to evolve both in quiescence and in outburst (as implemented in NSCool by Ootes et al., 2016). The data could be reproduced if the depth and magnitude of shallow heating were allowed to vary between different outbursts, although no clear correlation was found between the shallow heating parameters and the outburst properties. Moreover, it was found that the crust requires $\sim$1500 days to fully cool, which is longer than the recurrence time of the outbursts in Aql X-1 (Ootes et al., 2018).

A new bright outburst from Aql X-1 was detected in 2016 July (e.g. Sanna et al., 2016). Here we report on Chandra and Swift observations of Aql X-1 obtained after this bright outburst ceased. This study was aimed as a further test of the crust cooling hypothesis for the quiescent X-ray emission of Aql X-1, and to constrain the properties of shallow heating based on that assumption.

2 Observations and data reduction

2.1 Chandra observations

As part of our Chandra ToO program, Aql X-1 was observed at $\sim$50, 150, and 250 days after the end of its 2016 outburst. For all these observations the source was placed on the ACIS-S chip, using a 1/8 sub-array and the faint, timed data mode.

The outburst ceased in 2016 September and our first Chandra observation, lasting $\sim$14 ks, was obtained on November 22 starting at 15:57 ut. The second observation was performed on 2017 February 20, when it was observed for $\sim$18 ks starting at 12:43 ut. The third and final Chandra observation was carried out on 2017 May 30, starting at 17:57 ut, and had an exposure time of $\sim$23 ks. Unfortunately, Aql X-1 had just gone into outburst at this time (Vlasyuk & Spiridonova, 2017), after a relatively short quiescent phase of $\sim$235 days. Since the X-rays in this third Chandra observation will thus track the accretion emission and not that of the cooling neutron star, it cannot be used for our purposes and we will not discuss it further.

The Chandra data were reduced within ciao (v. 4.9). A circular region with a radius of $1.5^{\prime\prime}$ was used to extract source events and a $15^{\prime\prime}$-radius circular region was used for the background. We used dmextract to extract count rates and light curves. During the first and second Chandra observation, Aql X-1 was detected at a net 0.3–10 keV count rate of $(9.45\pm 0.26)\times 10^{-2}$ and $(5.92\pm 0.21)\times 10^{-2}~{}\mathrm{c~{}s}^{-1}$, respectively. No prominent variability is seen in the source light curves. We used specextract to extract spectra and to create response files. All spectral data were grouped to contain a minimum of 20 photons per bin using grppha.

2.2 Swift observations

2.2.1 Outburst and quiescent XRT light curve

Aql X-1 was monitored with Swift during its 2016 outburst and the subsequent quiescent phase, with the exclusion of a $\sim$3 month period during which the source was too close to the Sun. We used the Swift/XRT data performed in quiescence to obtain additional constraints on the temperature evolution of the neutron star. To determine which observations were suited for our quiescent analysis, we first produced an XRT count rate light curve. This was obtained using the online XRT repository222www.swift.ac.uk/user_objects/ (Evans et al., 2007, 2009), and is shown in Figure 1.

One of our prime aims was to compare the quiescent evolution of Aql X-1 after its bright 2016 outburst with that observed after its bright 2011 and 2013 outbursts (studied by Waterhouse et al., 2016; Ootes et al., 2018). Figure 1 illustrates the striking resemblance between these three outbursts in terms of duration, peak flux, and overall shape. This is further illustrated by the Swift/BAT (Krimm et al., 2013) and MAXI (Matsuoka et al., 2009) monitoring light curves, which are shown in Figure 2. This suggest that the spectral state evolution was very similar during the three outbursts (see also Díaz Trigo et al., 2018).

The XRT light curve of the 2016 outburst shows a transition from a rapid decay that takes place over $\sim$2 weeks time (days $\sim$44–57 in Figure 1) to a much slower decay that continues until the new outburst commences. Such a prominent change in decay rate has been seen more often in neutron stars that were monitored for their crust cooling, and has been interpreted as the transition from the outburst decay to quiescence (e.g. Fridriksson et al., 2010; Homan et al., 2014; Parikh et al., 2017a, b). As in Waterhouse et al. (2016) and Ootes et al. (2018), we fitted the two different decay parts to exponential functions to estimate the onset of quiescence, $t_{0}$, as the intercept of the two. This yielded $t_{0}=$ MJD 57664 (2016 October 3), with exponential decay time scales of $\sim$1.7 and 145 days for the rapid and slow decay, respectively. This is similar to the results obtained for the 2013 outburst (Waterhouse et al., 2016).333No Swift/XRT observations were available during the decay of the 2011 outburst due to Sun constraints (see Figure 1).

A first comparison of the quiescent evolution after the 2011, 2013 and 2016 outbursts is provided by the XRT light curves shown in Figure 1. It is of note that monitoring after the 2013 outburst stopped $\sim 100$ days into quiescence due to Sun-constraints, and for the same reason the decay and first $\sim 100$ days of quiescence were missed for the 2011 outburst (Waterhouse et al., 2016). Nevertheless, we can see from Figure 1 that the XRT count rate light curve after the 2016 outburst deviates from that obtained after the 2013 outburst, i.e. up to $\sim 100$ days into quiescence. This provides a first hint that the early temperature evolution was not the same after (and perhaps during) the 2013 and 2016 outbursts.

For the 2011 outburst we can only compare the late-time quiescent evolution, i.e. at $\gtrsim$100 days into quiescence. The first series of data points obtained after the 2011 outburst have no direct overlap with the 2016 data set, but the count rate is higher than the those measured shortly after the 2016 outburst. At $\sim$150 days in quiescence, when we can first directly compare the two data sets, the absolute count rates and evolution largely overlap (see Figure 1).

2.2.2 Quiescent XRT spectra

To obtain temperature measurements of Aql X-1 in quiescence in addition to those provided by our two Chandra observations, we used 11 XRT observations performed between 2016 October 10 and November 7 (obsID 00034719016–26) and 10 observations performed between 2017 March 16 and May 10 (obsID 00033665090–99). The XRT light curve shown in Figure 1 suggests that there were no accretion flares around the time of our Chandra observations.

Given the short exposures and low count rates (see Figure 1), the number of counts obtained for individual XRT observations was low ($\sim$10–30 counts). Therefore, we combined subsequent observations with similar count rates for our spectral extraction. To facilitate fits with a two-component model (see Section 3.1), we aimed for a total of $\sim$100–200 counts per composite spectrum. This resulted in 3 different Swift/XRT epochs (see Table 1).

All XRT observations that we used for our spectral analysis were obtained in PC mode. Reduction of these Swift data was performed within heasoft (v. 6.23). As an initial step, all observations were reprocessed with the xrtpipeline, using standard quality cuts. We next used XSelect to extract accumulated spectra. To obtain source spectra we used a circular extraction region with a radius of $30^{\prime\prime}$. Three circular regions with $30^{\prime\prime}$-radii, placed away from the source, were used to obtain background spectra. After summing the exposuremaps for the different observations in a particular epoch with XImage, an arf was created using xrtmkarf. The appropriate response matrix file (v. 14) was obtained from the calibration data base. All obtained spectra were grouped to contain a minimum of 15 photons per bin using grppha.

2.2.3 Spectral fitting

All spectral fits were performed within XSpec (v. 12.10; Arnaud, 1996). Based on previous quiescent studies of Aql X-1, we fitted the quiescent spectral data to a two-component model comprised of a neutron star atmosphere component (nsatmos; Heinke et al., 2007), and a power-law component (pegpwrlw). We modeled the interstellar absorption with tbabs, using vern cross-sections and wilm abundances (Verner et al., 1996; Wilms et al., 2000).

As in Waterhouse et al. (2016), we fixed the nsatmos-model parameters $M=1.6~{}\mathrm{M}_{\odot}$, $R=11$ km, $D=5$ kpc, and $N=1$ (where the latter implies that the neutron star surface is uniformly radiating), and only left the neutron star temperature ($kT_{\mathrm{eff}}$) free to vary.444The temperature measured in the neutron star frame was converted to that of an observer in infinity via $kT^{\infty}_{\mathrm{eff}}=kT_{\mathrm{eff}}(1+z)^{-1}$, where $(1+z)=1/\sqrt{1-2GM/(Rc^{2})}$ is the gravitational redshift. Since we perform our spectral fits for $M=1.6~{}\mathrm{M}_{\odot}$ and $R=11$ km, we used $(1+z)=1.33$. For the pegpwrlw model, the energy boundaries were set such that the model normalization represents the unabsorbed power-law flux in the 0.5–10 keV band. Furthermore, as the power-law index was poorly constrained in our data set, we fixed this parameter in all our spectral fits (see Section 3.1).

To obtain our final results, all quiescent Chandra and Swift spectra were fitted simultaneously with the hydrogen column density ($N_{\mathrm{H}}$) tied between all epochs. The total and thermal 0.5–10 keV unabsorbed fluxes and errors were determined with the cflux model. All errors quoted in this work reflect 1$\sigma$ confidence levels.

3 Analysis and results

3.1 Quiescent spectral evolution after the 2016 outburst

We first investigated the two individual Chandra observations. For both observations one-component fits with an nastmos model leave strong residuals at energies above 3 keV, leading to unacceptable fits ($\chi_{\nu}^{2}$/dof $=2.55/46$ for observation 1 and $\chi_{\nu}^{2}$/dof $=1.72/41$ for observation 2). Adding a power-law component with $\Gamma=1.7$ (see next paragraph), significantly improved the fits ($\chi_{\nu}^{2}$/dof $=1.07/45$ for observation 1 and $\chi_{\nu}^{2}$/dof $=1.00/40$ for observation 2), demonstrating that an additional spectral component is statistically required to describe the data (with f-test probabilities of $3\times 10^{-10}$ and $2\times 10^{-6}$, respectively).

Leaving the power-law index free in the Chandra fits resulted in $\Gamma=1.0\pm 0.6$ for the first spectrum, while this parameter was completely unconstrained for the second. Studying the decay of the bright 2010 outburst of Aql X-1 with Chandra and XMM-Newton, Campana et al. (2014) showed that the power-law index was consistent with being constant at a value of $\Gamma=1.7\pm 0.1$. It is not clear if the power-law index should be the same during the outburst decay and in quiescence, nor if the power-law index should be constant in quiescence. Various studies of Aql X-1 in quiescence report a range of power-law indices (\Gamma$$\sim0.8–2.7; e.g. Cackett et al., 2011; Marino et al., 2018). In the Swift quiescent study of Waterhouse et al. (2016), all spectra were fitted with $\Gamma=1.7$ fixed, based on the results of Campana et al. (2014). To allow for a direct comparison with the outbursts studied in Waterhouse et al. (2016), we choose to fix the power-law index to the same value. To test if this choice affects our conclusions, we also carried out fits for $\Gamma=1.0$ (i.e. the value obtained for our first Chandra spectrum).

In the $\sim 3$ months separating our two Chandra observations, the temperature of the neutron star atmosphere decreased from $117.5\pm 1.3$ eV to $110.7\pm 1.1$ eV, and the corresponding thermal flux decreased by $\sim$30% (see Table 1). The flux in the power-law spectral component also decreased in strength during this time. As a result, the fractional contribution of this hard emission tail to the overall unabsorbed 0.5–10 keV flux is similar for the two observations; $\sim$20% in the first and $\sim$14% in the second. Figure 3 shows the two Chandra spectra.

When using an absorbed nsatmos model to describe the XRT spectra, we find excess emission above the model fit at energies $\gtrsim$2 keV, most notably for the second and third Swift epoch. This is shown in Figure 4 and suggests that a hard emission tail is present in the XRT spectra as well. We therefore included a power-law component, with a fixed index of $\Gamma=1.7$ and variable normalization, when fitting the XRT data.

Table 1 summarizes the results of fitting the three Swift/XRT data sets together with the two Chandra spectra. The temperature of the neutron star atmosphere is observed to decrease from $118.9\pm 4.5$ eV in 2016 October to $102.6\pm 2.5$ eV in 2017 March–May. The 0.5–10 keV flux of the neutron star atmosphere is observed to decrease accordingly. All spectra are dominated by the neutron star atmosphere component, which has a fractional contribution to the total 0.5–10 keV flux of $\sim$70% to 90%.

Looking at the temperature evolution in more detail, we see that during the first $\sim$50 days (covered by the first two Swift epochs and the first Chandra observation) the neutron star temperature does not strongly change, while there is clear decrease in the $\sim$50 days separating the two Chandra observations. As can be seen in Table 1 and Figure 5, the temperature had decreased further by the time of the last Swift epoch.

In Figure 5 we compare the temperature evolution inferred from our Chandra and Swift analysis for the 2016 outburst with that observed by Swift/XRT after the 2011 and 2013 outbursts (from Waterhouse et al., 2016), which were analysed in the same way as we do here. The comparison with the 2011 (cyan squares) and 2013 (red filled circles) outbursts is of particular interest because these had such similar properties as the 2016 outburst that we study here (see Sections 2.2.1 and 3.2). It is immediately clear that at early times the absolute temperatures after the 2016 outburst are lower, and the temperature evolution is flatter, than seen after the 2013 outburst. There is thus clearly a difference in quiescent behavior after the two outbursts, as is also illustrated by the Swift/XRT quiescent light curves (Section 2.2.1). It is worth emphasizing that the entire shape of the cooling curve is different, i.e. there is not simply a systematic shift in temperature between the different years.

Modeling all spectra with a power-law component of $\Gamma=1.0$ yields a systematic upward shift in temperature by $\sim$1–2 eV, which is well within the typical $1\sigma$ errors that we obtain (see Table 1). Moreover, our main conclusions are based on the relative temperature evolution (i.e. the shape of the cooling curve), in particular in comparison with the 2013 data, which are unaffected by this systematic shift, hence our choice of $\Gamma$. Having measured the temperature evolution of Aql X-1 following the end of its 2016 outburst, we proceed by modeling these data with a thermal evolution model to put physical constraints on the properties of the crust.

3.2 Properties of the 2016 outburst

Since the energy released in the crust of neutron stars due to nuclear reactions is proportional to the mass-accretion rate (e.g. Haensel & Zdunik, 2008; Steiner, 2012), the outburst properties have to be taken into account in crust cooling simulations (e.g. Brown & Cumming, 2009; Ootes et al., 2016). Ootes et al. (2018) calculated the long-term (1996–2015) bolometric flux light curve of Aql X-1 from RXTE, MAXI, and Swift data by using instrument-specific count-rate conversion factors for hard and soft spectral states (see their Table 1). Applying the same method to calculate the energetics of the 2016 outburst yields an average mass-accretion rate of $3.2\times 10^{-9}~{}\mathrm{M_{\odot}~{}yr}^{-1}$. This is $\sim$20–30% higher than the values obtained for the 2011 and 2013 outbursts ($2.7\times 10^{-9}$ and $2.3\times 10^{-9}~{}\mathrm{M_{\odot}~{}yr}^{-1}$, respectively; Ootes et al., 2018).

Another ingredient in crust cooling simulations that can, in principle, be constrained from observations is the amount of light elements present in the accreted envelope. This determines how the observed effective temperature maps on to the actual interior temperature. As mass is accreted, the amount of light elements in the envelope increases but it suddenly drops when an X-ray burst occurs (since the light elements are then fused into heavier elements). We did not find any reports of X-ray bursts during the 2016 outburst in the literature, nor are there any X-ray bursts detected in the Swift/XRT observations. However, since Aql X-1 is known to display X-ray bursts regularly (e.g. Galloway et al., 2008), it is very likely that these were simply missed due to limited sampling and observing time. This implies that we cannot determine the last instance at which an X-ray burst occurred during the 2016 outburst, which would provide an upper limit on the He column depth accumulated before the start of quiescence (because it is plausible that a later X-ray burst was missed, which would lower the amount of light elements). We therefore left this parameter free in our crust cooling simulations (Section 3.3).

3.3 Crust cooling simulations

Based on the properties we infer for the 2016 outburst (Section 3.2), and taking into account the outburst history of Aql X-1, we modeled the temperatures obtained in the subsequent quiescent phase with the thermal evolution code NSCool (Page, 2016; Ootes et al., 2016, 2018). To be consistent with our spectral fits, we assumed $M=1.6~{}\mathrm{M}_{\odot}$ and $R=11$ km in all thermal evolution simulations. The impurity parameter of the crust ($Q_{\mathrm{imp}}$) was assumed to be 1, since Ootes et al. (2018) did not find any evidence of a higher impurity parameter when modeling the long-term outburst and quiescence data of Aql X-1. This is not surprising, since at the high crust temperatures found for Aql X-1, electron–impurity scatterings do not influence the thermal conductivity (see for instance figure 9 of Page & Reddy, 2012).

Modeling the outburst (and cooling) history is particularly important for Aql X-1 because the recurrence time of the outbursts is too short to allow the crust to fully cool (Ootes et al., 2018). In other words, the observed crust cooling behavior seen after an outburst depends on the crust temperature profile when the outburst commences, and in case of Aql X-1 this depends on the accretion history. In first instance, we ran a fit keeping all parameters for the previous outbursts fixed at the values found by Ootes et al. (2018), including the core temperature in the neutron star frame ($T_{0}=8.9\times 10^{7}$ K). In this simulation therefore only the envelope composition ($Y_{\mathrm{C}}$), shallow heating strength ($Q_{\mathrm{sh}}$), and shallow heating depth ($\rho_{\mathrm{sh,min}}$) were free to vary for the 2016 outburst.

The results of our simulations are listed in Table 2 and shown as the black solid curve in Figure 6. We find that reproducing the data requires a large amount of shallow heat during the 2016 outburst, located relatively deep in the crust. The magnitude of shallow heating is a factor $\sim$2.5–4 higher, and the depth is a factor $\sim 6$ larger, than obtained for the 2011 and 2013 outbursts of Aql X-1 (Ootes et al., 2018), which are listed in Table 2 for comparison.555It is of note that due to the lack of data within the first $\sim 100$ days of quiescence (see Figure 1), the shallow heating parameters for the 2011 outburst are less well constrained than for the 2013 and 2016 outbursts. The fitted light envelope composition is not strongly constrained, but the obtained values are similar for the three outbursts.

The main observational features that determine our derived shallow heating parameters are the flatness of the cooling curve in the first 50 days of quiescence, and the temperature drop that is observed after this. Indeed the reason that, in the model for the 2016 outburst, the shallow heating must be located very deep, results from the fact that the cooling starts relatively late.666For comparison, in MAXI J0556–332 for instance, the strong temperature drop that constrains the depth of the shallow heating is observed at $\sim$10–20 days after the outburst (Parikh et al., 2017a). Moreover, the inferred strength of the shallow heating is driven by the relatively strong temperature drop observed at this time. We note that although the inferred strength of the shallow heating is higher for the 2016 outburst than for the 2013 one, the temperature observed shortly after the outburst is lower for 2016, because the shallow heating is located much deeper (see Table 2).

Seeking alternative ways to explain the sudden temperature drop at $\sim 100-150$ days into quiescence, we considered the possibility that Aql X-1 might have a very cold core. As pointed out by Ootes et al. (2018), we likely have never observed Aql X-1 at its true base level because the crust does not have time to fully cool in between outbursts. We therefore performed another run in which we fit for the core temperature in addition to the envelope composition and shallow heating parameters. The results are listed in Table 2 and shown as the red dashed curve in Figure 6. Since a lower core temperature results in a stronger temperature gradient, hence stronger cooling, the amount of required shallow heating is reduced. However, we still need a significant amount of shallow heating ($\sim 5$ MeV nucleon*-1*) to explain the observed quiescent temperature evolution. The requirement for having this heat located relatively deep in the crust also remains, since it is determined by the late time of the temperature drop.

We emphasize that the need for strong shallow heating, relatively deep in the crust, is driven by the shape of the cooling curve and not the absolute temperatures. As stated in Section 3.1, the temperatures inferred from spectral fitting shift $\sim 1-2$ eV upward if a power law with $\Gamma=1.0$, instead of $\Gamma=1.7$, is used. This small systematic increase in temperature does not affect our conclusions about the depth and strength of shallow heating.

4 Discussion

We report on Chandra and Swift observations obtained after the bright 2016 outburst of Aql X-1 ended. The aim was to i) further test the hypothesis that the bulk of the quiescent flux evolution of Aql X-1 is driven by cooling of an accretion-heated crust, and ii) based on this assumption, better constrain the properties of the neutron star crust. In particular, we are interested in gaining more insight into the nature of the puzzling source of shallow heating that has been inferred for several neutron star LMXBs, including Aql X-1, in the crust-cooling hypothesis (e.g. Page & Reddy, 2013; Degenaar et al., 2014, 2015; Parikh et al., 2017b, 2019; Ootes et al., 2016, 2018). If crust cooling can be studied in Aql X-1, its frequent outbursts provide the opportunity to study the properties of this shallow heating after several different outbursts, thereby allowing to break degeneracies with neutron star specific parameters that may be involved in shallow heating and do not change between different outbursts (e.g. mass, spin, magnetic field strength; Waterhouse et al., 2016).

There are two other neutron stars for which crust cooling has been studied after different outbursts; MAXI J0556–332 (Homan et al., 2014; Parikh et al., 2017a) and MXB 1659–298 (e.g. Wijnands et al., 2004; Cackett et al., 2013b; Parikh et al., 2019). For MAXI J0556–332, which exhibited 3 outbursts of different duration and peak intensity, there is clear evidence that the depth and magnitude of shallow heating varies between the different outbursts. This rules out that basic neutron star parameters play an important role in regulating shallow heating (Parikh et al., 2017a). MXB 1659–298, on the other hand, showed remarkably consistent heating parameters for 2 outbursts that were of similar brightness but different duration (2.5 and 1.5 yr; Parikh et al., 2019). In the case of Aql X-1, it was previously noted by Waterhouse et al. (2016) that the 2011 and 2013 outburst properties were very similar, and that the quiescent Swift/XRT count rate light curves and temperature evolution also gave consistent results between the two outbursts. Indeed, the detailed study of Ootes et al. (2018) yielded similar shallow heating parameters for the 2011 and 2013 outbursts (see also Table 2).777A note of caution is that there is only short time overlap in the data; for the 2013 outburst there is only quiescent coverage up to $\sim 100$ days into quiescence, whereas monitoring after the 2011 outburst did not start until this point (see Figure 1 and 5).

Based on X-ray burst studies of the persistently accreting neutron star LMXB 4U 1820–30, it was previously suggested that differences in the accretion geometry could perhaps lead to different levels of shallow heating (in ’t Zand et al., 2012). Aql X-1 is a promising target to test this idea, because it exhibits different classes of outbursts that have different spectral state evolution and hence likely different accretion geometries (e.g. Maitra & Bailyn, 2008; Asai et al., 2012). In modeling multiple outbursts of Aql X-1, Ootes et al. (2018) tested if shallow heating could be operational only in a particular spectral state (hard or soft), but ruled out such a simple connection. The data that we present here provides another opportunity to test any possible link between shallow heating and spectral states. The MAXI and Swift (BAT and XRT) light curves of the 2016 outburst of Aql X-1 are remarkably similar to that of the bright 2011 and 2013 outbursts (see Figures 1 and 2). This makes it likely that the spectral state behavior (i.e. the accretion geometry) was very similar during the three outbursts (which was also noted by Díaz Trigo et al., 2018). If the global outburst properties are the main driver in determining the shallow heating, we would therefore expect to see similar crust cooling behavior after the three outbursts.

The observed early thermal evolution after the 2016 outburst was, however, markedly different than seen after the 2013 one. The neutron star temperature was much lower right after the outburst ended, and it did not strongly evolve until $\sim$100 days into quiescence, when a decrease in temperature is observed. By using thermal evolution simulations, we determined that, within our current understanding of heating and cooling models, the late temperature drop seen after the 2016 outburst can only be achieved if the shallow heating was much stronger, and located much deeper in the crust, than for the 2013 outburst. Our results imply that the spectral state behavior during an outburst (i.e. the accretion geometry) cannot play a strong role in setting the properties of shallow heating.

4.1 Implications for our understanding of shallow heating

Whereas the concept of shallow heating has been known for over a decade now (e.g. Brown & Cumming, 2009), its physical origin is still not understood. Nevertheless, ongoing crust cooling studies are providing us with pieces of information regarding the mechanism of shallow heating. In particular, the studies of crust cooling of 3 sources that exhibited multiple outbursts (MAXI J0556–332, MXB 1659–28 and Aql X-1) has established that shallow heating does not strongly depend on i) neutron star parameters that do not change between outbursts (this work and Parikh et al., 2017a), and ii) the accretion geometry (this work and Ootes et al., 2018). Therefore, there must be other factors that drive the shallow heating.

One framework in which different shallow heating for similar types of outbursts can be accommodated, is that of chemical convection. In this model, the mixture of elements left in the neutron star envelope after thermonuclear burning organizes itself into layers of light and heavy elements (Horowitz et al., 2007). This chemical separation may drive a convective heat flux that can potentially heat the outer layers of the crust (Medin & Cumming, 2011, 2014, 2015). Depending on the last instance at which an X-ray burst occurred before accretion switched off, the envelope composition can differ, even if the outbursts are very similar (Brown et al., 2002). Chemical convection will act differently depending on the envelope composition and could therefore lead to different shallow heating. However, it can likely generate only a few tenths of an MeV of energy per accreted nucleon (Medin & Cumming, 2015), and therefore cannot account for the strong shallow heating ($\sim 5-10$ MeV nucleon*-1*, depending on the core temperature) that we infer for the 2016 outburst of Aql X-1. Similar conclusions were drawn for the strong shallow heating inferred for the main outburst of MAXI J0556–332 (Deibel et al., 2015; Parikh et al., 2017a).

One source of energy that can potentially provide the strong shallow heating inferred for MAXI J0556–332 and Aql X-1, is the orbital energy of the accreted material (Inogamov & Sunyaev, 2010). It was previously suggested, for MAXI J0556–332, that dissipation of accretion-generated oscillations (“g modes”) in the liquid part of the neutron star crust (the “ocean”) could inject heat in the crust (Deibel et al., 2015; Deibel, 2016). However, since the observed properties of 2016 outburst of Aql X-1 were so similar to that of its 2011 and 2013 outbursts, it seems difficult to understand why the heat injected in the crust would not be similar. This is, at least, under the assumption that the X-ray properties are a reasonable proxy for the mass-accretion rate onto the neutron star, which may not necessarily be true (e.g. van der Klis, 2001).

Our results for the 2016 outburst of Aql X-1 do not only stand out from its 2013 (and 2011) outburst, but also from the other crust cooling sources. Indeed, the magnitude of the shallow heating that we require is much higher ($\sim 5-10$ MeV nucleon*-1*) than inferred for most other crust cooling sources, which typically require 1–2 MeV nucleon*-1* (e.g. Page & Reddy, 2013; Degenaar et al., 2014, 2015; Parikh et al., 2017b, 2019; Ootes et al., 2016, 2018), except for MAXI J0556–332 during its main outburst, for which $\sim 15-17$ MeV nucleon*-1* was inferred (Deibel et al., 2015; Parikh et al., 2017a). The depth of shallow heating that we infer for Aql X-1, $\sim 10^{10}$ g cm*-3*, is also higher than the typical values inferred for other crust cooling sources (using NSCool); these are often on the order of $10^{8}-10^{9}$ g cm*-3* (e.g. for KS 1731–260 and 1RXS J180408.9–342058; Ootes et al., 2016; Parikh et al., 2017b), although sometimes the constraints are poor and much higher densities of $\sim 10^{10}$ g cm*-3* are allowed (e.g. for MXB 1659–29 and the second outburst of MAXI J0556–332; Parikh et al., 2017a, 2019).

Although the shallow heating depth for the 2016 outburst is high, it is not excessive. For example, for the main outburst of MAXI J0556–332, which also required the very strong shallow heating, the inferred depth was $\sim 5\times 10^{9}$ g cm*-3* Parikh et al. (2017a). We also note that this depth is still well below that of deep crustal heating, which occurs at densities of $10^{12}-10^{13}$ g cm*-3* (e.g. Haensel & Zdunik, 2008). Steiner (2012) proposed that uncertainties in the nucleon symmetry energy could allow for much stronger deep crustal heating ($\sim 5$ MeV nucleon*-1*) than previously assumed ($\sim 2$ MeV nucleon*-1*), although more detailed calculations are needed to confirm this (Fantina et al., 2018). Moreover, it is unlikely that a different amount of deep shallow heating can account for our results of Aql X-1, as it acts at so much higher densities and would thus manifest itself at much later time in the cooling curve (not probed by our observations; Brown & Cumming, 2009).

We can speculate that how, or when, shallow heating is operating, somehow depends on the initial temperature of the crust when an outburst commences. The outburst that preceded the 2016 activity of Aql X-1 was shorter and fainter than the outbursts that occurred before the 2011 and 2013 ones (see Ootes et al., 2018). Since the crust cooling time of this neutron star is shorter than the quiescent time between two outbursts, the initial temperature of the crust can significantly differ between outbursts, and perhaps this can explain the different levels of shallow heating inferred. A similar argument might apply for MAXI J0556–332; when its second and third outburst occurred, the crust had not yet relaxed from the strong heating of its first (main) outburst (Parikh et al., 2017a). Perhaps the much higher crust temperature at the start of the second and third outburst is related to the fact that the shallow heating inferred for these outbursts was very different from that obtained for the first outburst ($\lesssim 2$ and $\sim 17$ MeV nucleon*-1*, respectively; Parikh et al., 2017a). For MXB 1659–29, on the other hand, the crust temperature at the start of its last two outbursts was likely similar (Parikh et al., 2017b), and perhaps that can explain why there was no apparent difference in shallow heating for this source.888Although the quiescent phases preceding the 1999–2001 and 2015–2017 outbursts differed by a few years, the quiescent temperature evolves only slowly after a decade. It is therefore likely that the crust temperature was similar before both outbursts.

4.2 Other explanations for the quiescent evolution of Aql X-1

Other than crust cooling, the quiescent flux evolution of Aql X-1 could possibly be powered by ongoing low-level accretion (e.g. Kuulkers et al., 2009; Coti Zelati et al., 2014). Nevertheless, in this interpretation it may not obvious either why the quiescent evolution of Aql X-1 after its 2016 outburst should be significantly different from that observed after the 2013 outbursts. The similarity in outburst properties suggests that a similar part of the disk was involved in the outburst, and that a comparable amount of mass drained from the disk on to the neutron star. Therefore, one might expect a similar quiescent evolution, opposed to what is seen.

In assessing the low-level accretion scenario for Aql X-1, it is worth noting that the fractional contribution of the hard emission component to the 0.5–10 keV unabsorbed flux is generally low; $\lesssim$20% in that data that we presented here and in Waterhouse et al. (2016). Wijnands et al. (2015) proposed that if the power-law and thermal component both contribute $\sim$50% to the total quiescent flux, the emission is likely powered by accretion, while a lower power-law contribution may point to a different origin. Indeed, the most proximate neutron star LMXB, Cen X-4, exhibits equal flux contributions of its two quiescent emission components (e.g. Cackett et al., 2010), and there is strong evidence that both are powered by low-level accretion onto the neutron star surface (Bernardini et al., 2013; Cackett et al., 2013a; Chakrabarty et al., 2014a; D’Angelo et al., 2015). For Aql X-1, the power-law contribution is generally lower than $\lesssim$50% (e.g. Marino et al., 2018). Nevertheless, the fractional contribution of the power-law remained approximately constant between our two Chandra observations, i.e. the neutron star atmosphere and the power-law component decreased in tandem. This suggests that the two components are connected, possibly arising from the same emission process, which could be more naturally explained by continued accretion (e.g. Cackett et al., 2010).

Other than turning to a different explanation for the quiescent emission altogether, it is also possible that our results on Aql X-1 are exposing that we are still missing something in our understanding of heating and cooling of neutron star crusts. For instance, as noted earlier, a different envelope composition would in principle produce a systematic temperature shift, because it yields a different mapping between the surface and interior temperature of the neutron star, and hence not explain the different shape of the 2016 cooling curve of Aql X-1. However, it was very recently shown that the envelope composition can potentially significantly change during quiescence as a result of diffuse nuclear burning, a process in which elements diffuse to such depths where the density and temperature are sufficiently high to ignite nuclear burning (Wijngaarden et al., 2019). The impact of this process on crust cooling studies, and if it can potentially explain the observed behavior of Aql X-1, needs to be explored.

Acknowledgements

ND is supported by a Vidi grant awarded by the Netherlands organization for scientific research (NWO). LO, RW and AP are supported by an NWO top grant, module 1, awarded to RW. DP is partially supported by the Consejo Nacional de Ciencia y Tecnologíaa with a CB-2014-1 grant #240512. DA acknowledges support from the Royal Society. JMM and JH acknowledge support from Chandra grant GO7-18031B. This work made use of data supplied by the UK Swift Science Data Centre at the University of Leicester. We also made use of MAXI data provided by RIKEN, JAXA and the MAXI team, and public light curves from the Swift/BAT transient project.