Probing outbursts of the transient neutron star low mass X-ray binary Aql X-1 with NICER: a study of spectral evolution

X-ray observations of neutron star (NS) low mass X-ray binaries (LMXBs) are useful to probe physical processes close to the NS and to constrain source parameters. Aql X-1 is a transient NS LMXB which frequently undergoes outbursts provides an excellent opportunity to study source properties and accretion mechanism in strong gravity regime over a wide range of accretion rates. In this work, we systematically investigate the spectral evolution of Aql X-1 using NICER observations during the source outbursts in 2019 and 2020. The NICER observations cover the complete transition of the source from its canonical hard state to soft state and back. The spectra extracted from most observations can be explained by a partially Comptonised accretion disc. We find that the system can be described by an accretion disk with an inner temperature of $\approx0.7$ keV and a Comptonising medium of thermal electrons at $\approx2$ keV, while the photon index is strongly degenerate with the covering fraction of the medium. We also find evidence of Fe K$\alpha$ fluorescence emission in the spectra indicating reprocessing of the Comptonised photons. We observe an absorption column density higher than the Galactic column density for most of the observations indicating a significant local absorption. But for some of the observations in 2020 outburst, the local absorption is negligible.


INTRODUCTION
Low mass X-ray binaries (LMXB) consists of a neutron star (NS) or a black hole (BH) and a companion star with a mass ≲ 1 M ⊙ .In these binaries, the compact object accretes the matter from the companion star through the Roche-Lobe overflow.The NS LMXBs are broadly categorised into 'atoll' sources or 'Z' sources based on the shape of the track generated on their hardness intensity diagram (HID) or the colour-colour diagram (CCD) (see van der Klis 2004, for a review).Most of the NS LMXBs are transient systems which stay in the quiescent phase for long duration and undergo outbursts during which they are typically detected.These systems produce outbursts, due to the increased mass accretion rate, which are caused by the instabilities in the disc (Shakura & Sunyaev 1973;Done et al. 2007).During these outbursts, the luminosity of the system increases and it emits photons from the radio to the X-ray band (Done et al. 2007).
The X-ray emission from NS LMXBs has been modelled with a combination of spectral components.The primary component of the modelling is the accretion disc (Shakura & Sunyaev 1973), a blackbody component (from the accretion mound on the surface of the NS or from the boundary layer) and a non-thermal emission (either from the accretion column or the Comptonised emission).The Comptonised emission results from Compton up-scattering of the soft thermal photons from the accretion disc by a hot plasma of ★ E-mail: kputha@email.sc.edu electrons (often called corona) close to the NS.During the evolution of the source across various accretion states, there is a strong interplay between these components which decides the track followed by the source on the HID/CCD.Aquila X-1 (or Aql X-1, discovered by Kunte et al. 1973) is one of the most famous NS LMXBs which regularly goes into an outburst every year (although the outburst is not exactly periodic, Šimon 2002;Ootes et al. 2018).Due to the extensive observations of the source in various wavelengths, key parameters of the source are well known.Aql X-1 has a K-type donor and is placed at a distance of 6 ± 2 kpc (Mata Sánchez et al. 2017).The orbital period of the binary system is 18.95 d (Chevalier & Ilovaisky 1991;Welsh et al. 2000).The presence of the thermo-nuclear bursts in the X-ray light curves suggest the compact object is a NS (e.g.Galloway et al. 2008;Güver et al. 2022) and the detection of the millisecond pulsations in one such bursts (1.8 ms, Casella et al. 2008) confirms the nature of the accretor.Aql X-1 is one of the few NS LMXBs which shows a hysteresis in the HID (Maitra & Bailyn 2004;Güver et al. 2022) similar to BH LMXBs (e.g.Bhargava et al. 2019;Rawat et al. 2023, and references within).Using RXTE-PCA observations of Aql X-1, Maitra & Bailyn (2004) have studied the evolution of the source across the outburst in 2000 as it traverses various spectral states i.e. low hard state (LHS) and high soft state (HSS) and other intermediate states.The authors characterise the spectrum of the source as a combination of an accretion disc and a powerlaw component.Similar characterisation of the source evolution across different types of the outbursts (i.e.long high, medium low and short low) was conducted using monitoring observations of Aql X-1 with RXTE/PCA by Güngör et al. (2014).The long high outbursts are more frequent and during the majority of the outburst, the source stays in the HSS while tracing a hysteresis across the other states during the rise and the decay phase of the outburst (Maitra & Bailyn 2004).López-Navas et al. (2020) have undertaken preliminary investigation of the soft X-ray evolution of the source using Swift-XRT spectra but mainly focus on the correlation of the UV-optical properties and their correlation with the X-ray spectral evolution.Recently, Niwano et al. (2023) studied outbursts of the source from 2016-2019 using the MAXI instrument onboard ISS where they model the 2-20 keV spectrum as the emission from an absorbed accretion disc.
The soft X-ray evolution of NS LMXB holds the key to evolution of the thermal components (e.g.Bhargava et al. 2023).In case of Aql X-1, the soft X-ray spectrum is somewhat understudied.In this article, we investigate the properties of Aql X-1 in soft X-rays using extensive NICER monitoring during its outbursts in 2019 and 2020.We report the details of the observation and data reduction methods in section 2, spectral analysis methodology in section 3 and describe the results and our interpretation in section 4.

OBSERVATIONS AND DATA REDUCTION
NICER (Gendreau et al. 2016) has regularly observed Aql X-1 since its launch in 2017.Since then Aql X-1 has undergone several outbursts.For this analysis, we consider the outbursts that happened from 2019 August 2 -2019 September 21 (outburst A) and from 2020 August 15 -2020 October 23 (outburst B).These observations are considered for the analysis as they are representative of the complete outburst (from the rise to the decay).The details of the observations are reported in tables A1 and A2.In the draft, we refer to the observations with A/BXXX, where A/B denote which outburst they are part of and the XXX indicates the last three digits of the observation id (ObsID) as mentioned in tables A1 and A2.
The data is processed through a standard filtering pipeline, using nicerl2 i.e, excluding the detectors 14, 34 and excluding the time intervals corresponding to South Atlantic Anomaly (SAA) passages, the Earth elevation angle > 15 • , the elevation angle with respect to the bright Earth limb ≲ 30 • , and considering the intervals where the undershoot and overshoot count rates per module are 0-500 and 0-30 respectively.We used HEASOFT V6.31, NICERDAS V10a and CALDB version xti20221001.
Using nicerl3-lc, light curves are extracted in the energy ranges 0.5-12.0keV, 0.5-2.0keV, and 2.0-5.0 keV.For further analysis, we ignore the observations with the net exposure ≲ 200 s.We depict the light curves in 0.5-12 keV in the top panels of figure 1 and the hardness ratio of computed using 2-5 keV and 0.5-2 keV light curves in the bottom panels of figure 1.We also construct a HID of the source for both the outbursts in figure 2. Some of the observations analysed here show bursts in their light curve.Since we want to focus on the persistent emission of the source, we identify the time stamps of the bursts from the light curve using xselect and create Good Time Intervals (GTIs) excluding the bursts 1 .Information about these bursts are listed in table A3.To exclude the burst, we remove the interval during which the source exceeded the mean count rate before and after the burst, without assuming any decay timescale.These GTIs are then passed to nicerl3-spec routine while extracting the spectrum.
1 https://swift.gsfc.nasa.gov/analysis/threads/batlightcurvethread.html, is used to construct GTI The burst intervals have been excluded from the light curves plotted in the above-mentioned figures and the further spectral analysis.

SPECTRAL ANALYSIS
To investigate the spectral evolution of Aql X-1, we characterise the spectrum as seen by NICER in 0.5-10 keV for all the observations.For modelling the spectrum, we use the  2 -statistic and report the 1 uncertainty in the parameters throughout the article.We use the SCORPEON model (Markwardt et al. 2023) to estimate the background, and fix the background model parameters to the default values.In case of the fainter observations (count rate in 0.5-12 keV ≲ 100 cts/s, B149-B159), we restrict the spectrum to the energy ranges where the source counts are above the expected background level.A systematic error of 0.8% is used during the spectrum modelling in the operational energy range of 0.5-10 keV.We have used the optimal binning by the spectral response (Kaastra & Bleeker 2016) for the spectral modelling.We have depicted spectrum of observation A116 in the top panel of the figure 3.
To test the nature of the spectrum, we model the spectra from the individual observations with, a single additive component (viz.diskbb, bbodyrad or nthcomp), along with an absorption model to correct for the interstellar+intrinsic absorption (modelled with tbabs).The abundances from Wilms et al. (2000) and the crosssections from Verner et al. (1996) are used to model the spectrum.The single-component model is found to be insufficient to explain the spectrum with count rate >100 counts.For example, ObsID A121 the model tbabs*(diskbb) result in  2 /degrees of freedom (DoF) of 8852.37/158, but on inclusion of Comptonisation model thcomp (Zdziarski et al. 2020) result in  2 /DoF of 144.22/156.The drastic change in  2 for a minimal change in the DoF indicates a strong presences of a Comptonized emission.Thus we model all the spectra at least as a combination of a thermalised accretion disc and a Comptonised emission i.e, tbabs*(thcomp⊗diskbb).Across all ObsIDs (except the fainter observations, i.e.B151-B159), we identify residuals around 0.75 keV which can be modelled by a Gaussian.This feature has also been seen in other sources with NICER indicating that it is an instrumental feature (e.g.Zhang et al. 2023).For example, in A112, adding a Gaussian component reduces  2 /DoF from 191.97/143 to 153/140, resulting in a f-test statistic of 11.88 and probability of 5.55×10 −7 .Thus we model the observations by including the Gaussian in the model (tbabs*(gauss+thcomp⊗diskbb); Model A).The residuals of some spectra indicate the presence of a weak Fe K emission line which is modelled by a Gaussian, i.e Tbabs*(gauss+gauss+thcomp⊗diskbb) (Model B).For example, the addition of Gaussian for A116 (which shows a structured residual near 6.5 keV with Model A) results in a f-test value of 22.55 with probability 2.97×10 −9 , whereas adding a Gaussian to A112 (with no significant residual near 6.5 keV with Model A) results in an f-test value of 2.84 with probability 0.061.Some observations at the end of outburst B, i.e, from B151-B153 do not show evidence for the instrumental Gaussian component and are modelled by tbabs*(thcomp⊗diskbb); (Model C) and fainter observations B154-B159 are modelled by tbabs*(diskbb) (Model D) .We depict the spectrum and its spectral decomposition of A116 in the first panel of figure 3 and the residuals to the spectral fits in the bottom panels of figure 3. The second panel shows the residuals with both Gaussians, third panel shows the residual without the Gaussian at 6.5 keV and last panel displays the residual with out Gaussian at 0.75 keV.
The convolution model thcomp up-scatters a fraction of the seed Hardness ratio (2 5) keV/(0.5 2) keV  Hardness ratio is defined as the ratio between the 2-5 and 0.5-2 keV light curve is plotted against the intensity in 0.5-12 keV energy range.The inset plot shows a better view of the clustered points in the HSS.The marker scheme is same as Figure 1.photons (  ) from the disc and returns a combination of the unscattered disc emission and Comptonised emission (of a photon index Γ and electron temperature   ).Since convolution models requires computation of the model beyond the spectral fitting range, we used a custom energy grid from 0.01-100 keV using 1000 logarithmic bins.
If we keep all parameters free across the different observations, we find that some of the observations indicated a low Γ (∼1) while some required high Γ (∼1.8), while they lie at a similar HID position.These observations also show low and high  respectively.To probe if there is a degeneracy between these parameters, we generate  2 contour maps by sampling a uniform grid of Γ from 1.001 to 2 and  from 0 to 1 in the steps of 25.We find that for all observations these two parameters are degenerate but the degeneracy curve depends on the source position on the HID.We depict the contour plot for some of the observations from different parts of the HID in figure 4.
Due to a degeneracy in the parameter space, we cannot probe both parameters simultaneously.To understand the parameter evolution and observe the trends in different parameters as a function of time, we fix the  to a range of values 0.1, 0.3, 0.5, 0.75, and 1.0, and estimated the spectral parameters for all the observations.The evolution of the parameters is shown in figure 5. We note that not all observations result in acceptable fits (indicating a  2 /DoF > 1.5) on fixing  and the parameters of unacceptable fits are not shown in figure 5. We report the parameters for  = 0.75 in tables A4 and A5 as most of the observations have a  2 /DoF < 1.5.
We compute the unabsorbed flux in the energy range of 0.5-12.0keV using the convolution model cflux.We use the same energy grid as we did for the thcomp.For computing the flux, we freeze all the parameters and then add the cflux component, fitting the model.Then refit the model by thawing the non-normalisation parameters (e.g. in ) while keeping the Gaussian and f frozen to get Figure 4:  2 contours of observations from different parts of HID, displaying the degeneracy between the parameters Γ and f are depicted in different colours.The minimum of  2 for each observation are indicated by different markers (and the symbol is included in the legend).We depict the contours corresponding to 68% (solid), 90% (dashed), 99% (dashed-dot) and 99.99%(dotted) confidence levels.
the flux value.The flux remains the same irrespective of f.We report the unabsorbed flux in Tables A4 and A5.

RESULTS AND DISCUSSION
We investigated the evolution of the continuum properties of the LMXB Aql X-1 in the soft X-rays using the NICER observations of the outbursts of the source in 2019 and 2020 (Outburst A and B respectively).The outbursts covered the complete evolution of the source and revealed an interesting trend in the spectral properties of the source.

Light curves and HID
The NICER observations of Aql X-1 in 2019 and 2020 indicate an evolution of the source which had been extensively seen in the literature (Maitra & Bailyn 2004;Güngör et al. 2014;Niwano et al. 2023).The outbursts of the source investigated here (see figure 1) are the 'long high' kind with the profile resembling typical fast rise exponential decay (FRED, Chen et al. 1997).Such outbursts have been seen often in Aql X-1, indicating that these outbursts are typical of the source.
Aql X-1 is known to show a Q-shaped hysteresis in its HID (using the hardness of 9.7-16/6.4-9.7 keV from RXTE/PCA, e.g.Maitra & Bailyn 2004) as the source evolves in the outburst.In the case of the hardness defined as 9.7-16/6.4-9.7 keV, the count rate and the hardness typically show an anti-correlation (for a higher count rate, a lower hardness was observed and vice versa).The track changes its shape when a softer definition of the hardness was used (3.8-6.8/2-3.8keV using NICER, Güver et al. 2022, and 2-5/0.5-2keV figure 2), with a change in the trend between the count rate and hardness.For the softer definition of hardness, the count rate was proportional to hardness for the majority of the observations.The trend between the hardness and intensity depends mainly on the spectral shape of the source.A positive trend indicates that with increasing count rate, the harder band has more photons than the softer band and vice versa.In case of NICER observations, the harder band was 2-5 keV in the current work and 3.8-6.8keV in Güver et al. ( 2022  band was 0.5-2 keV and 2-3.8 keV respectively.Since in both of these cases, we see a positive trend between hardness and intensity, it indicated that the count rate was increasing by a larger amount in the bands 2-6.8 keV as compared to 0.5-2 keV.And the negative trend for a similar change in the count rate in 9.7-16/6.4-9.7 keV from RXTE/PCA observations indicated in 6.4-9.7 keV band the count rates were increasing more as compared to 9.7-16 keV.Thus, in a model-independent way, we ascertain that during an outburst, there was a strong spectral change, with ≈ 2-6 keV count rate changing with the total count rate more strongly as compared to the energy bands around it.The Q-shape hysteresis in Maitra & Bailyn (2004) shows that for a similar count rate, during the rise phase and the decay phase, the source could show different hardness and thus have a different spectral description.But the HID using a softer hardness ratio didn't depict the hysteresis and we observe that the path taken during the rise part of the outburst is similar to the decay part.Although the source showed jumps in hardness at similar counts (e.g.B114, B116, B150 had a higher hardness as compared to B113, B117 and B149) this is different from the traditional hysteresis where the transition from the LHS to HSS and HSS to LHS takes different paths.It was interesting to note that some of these observations have also indicated bursts (table A3).

Spectral and parameter evolution
The HID of the source indicates a clear evolution of the source throughout its outburst as it traverses through various states.In particular, the source stays for a significant amount of time in the outburst in the HSS, where the source spectrum is adequately described as a combination of a thermalised accretion disc and Compton upscattering of a fraction of the accretion disc photons.
Due to a limited energy range covered by NICER, the spectral model is found to have strong degeneracies.Mainly, we found that the different combinations of photon index (Γ) and covering fraction (f ) of the Comptonising medium result in similar  2 statistic (see figure 4 for a  2 contours for a sample of observations).Although there was a degeneracy between these two parameters, for individual observations the loci of the contours follow different tracks.This indicates, that although these parameters are degenerate, they have a systematic evolution as the source undergoes an outburst.Due to the degeneracy, we investigated the evolution of the parameters by fixing f for the duration of the outburst at different values.
The choice of the f did not affect the estimated value of the absorption column density ( H ) which has a typical value of 0.48-0.51for almost all observations.This value is higher than the line of sight column density (= 0.31 × 10 22 cm −2 , HI4PI Collaboration et al. 2016) indicating a local absorbing column.The  H observed is similar to the value reported in the literature (Gusinskaia et al. 2020;Li et al. 2021).Some of the observations (B113 -B117 and B150), however, suggest that the absorption column density is similar to the line of sight value (also seen in Aql X-1 in HXMT observations in a different outburst; Güngör et al. 2020).Notably, these observations also lie on a separate part of HID (see figure 2), occupying a typically harder region for a similar count rate in rest of the observations.This behaviour suggests that the local absorption column is perhaps caused due to a clumpy medium which disappears for a small part of the outburst.And since the local absorbing column variation is independent of assumed covering fraction of the corona, we can safely assume that they are separate entities in the system.
The measurement of the temperature of the inner edge of the accretion disc ( in ) depends on the assumption of f.A lower covering fraction requires a hotter accretion disc for a similar fit (see figure 5).At the same time, for a lower f, the normalisation of the accretion disc component (  ) is also lower.Barring these differences in the value offset, the trends for both parameters of the accretion disc are similar for different f.The trends for the accretion disc parameters also conform with each other in different outbursts, indicating a similar behaviour across different epochs.The evolution of  in is similar to count rate and hardness evolution i.e. for the observations in HSS, the temperature remains the same.Additionally for the observations in HSS, the disc normalization is observed to be constant indicating a stable disc during the HSS (e.g.Done et al. 2007).The RXTE observations have also built a similar picture of the source in the HSS (Maitra & Bailyn 2004) but the actual values differ due to the energy range considered and the model assumption used to describe the spectrum.As the source moves away from the HSS, the disc temperature reduces and the normalisation (and thus the inner radius) increases.
As expected from the  2 contours, Γ has a strong dependence on the assumed covering fraction.But for a given covering fraction, the Γ follows the count rate evolution inversely.For a low covering fraction, the Γ is consistent with 1.0 (see panels (d) and (i) in figure 5, where these points are indicated by triangles), which corresponds to an infinite optical depth of the corona.Electron temperature (  ) on the other hand is insensitive to the assumption of a fixed covering fraction (typically varying between 1.5-2 keV) but has a stronger deviation for the observations not in the HSS.In the HSS, the parameters of the corona are also roughly constant with minor fluctuations in the electron temperature which is also seen as the jitters in the source position on the HID.But as the source moves away from the HSS, the spectrum steepens with a simultaneous reduction in the accretion disc temperature and the normalization.Some of the observations indicate residuals near 6.5 keV after modelling the spectrum with Model A. Thus we modelled these residuals with a Gaussian component to account for the Fe K emission line at 6.5 keV (Reynolds & Nowak 2003).For example, for observation A116, the emission line can be clearly seen in the residuals plotted in the third panel of 3. The observations depicting significant Fe K emission are marked with crosses in figures 1, 2, and 5.In all cases of detection, the line width is ≈ 0.5 keV and line normalization is ≈ 5 × 10 −3 photons cm −2 s −1 .

CONCLUSIONS
We investigate the spectral properties of Aql X-1 as observed with NICER in 0.5-10 keV and infer the spectral components present in various canonical states of the source.The extensive NICER observations probed the evolution of the source using the sensitive soft X-ray coverage.The HID constructed using a softer definition of the X-ray colour revealed a track distinct from the hysteresis track typically observed for the source using a hard colour.We find that the spectrum of the source in the HSS can be well described as a partially Comptonised accretion disk.Some of the observations also exhibit reprocessing of the Comptonised emission as an iron line emission.Our analysis is limited by the degeneracy in the parameters of the Comptonising medium, likely due to limited energy range of the instrument.We find evidence for a variation in the local absorption column which also manifests as a separate position on the HID.
Division of the Smithsonian Astrophysical Observatory.KGP thanks Dr Aru Beri for the help in learning HEASoft.

APPENDIX A: TABLES
This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 1 :
Figure 1: The plots have been updated Outburst of Aql X-1 in 2019 (left panels) and 2020 (right panels).The top panels (both left and right) show the evolution of the count rate in 0.5-12 keV while the bottom panels show the evolution of the hardness ratio in 2-5/0.5-2keV energy bands.Each point corresponds to a single observation.The shape of the points have been assigned according to the spectral decomposition,i.e,Model A = tbabs*(gauss+thcomp⊗diskbb) by a circle, Model B = tbabs*(gauss+gauss+thcomp⊗diskbb) by a cross, Model C = tbabs*(thcomp⊗diskbb) by a diamond, and Model D = tbabs*(diskbb) by a square (see the text for a detailed description of the models).

Figure 2 :
Figure 2: The plot has been updated Hardness intensity diagram of 2019 (left panel) and 2020 (right panel) outbursts of Aql X-1.Hardness ratio is defined as the ratio between the 2-5 and 0.5-2 keV light curve is plotted against the intensity in 0.5-12 keV energy range.The inset plot shows a better view of the clustered points in the HSS.The marker scheme is same as Figure1.

Figure 3 :
Figure 3: Count spectra and spectral decomposition (top panel), and residuals ( = (Data -Model)/Error) without both Gaussians (second panel), without Gaussian for iron line (third panel) and with both Gaussians (bottom panel) for A116 observation with the model B having a fixed covering fraction f =0.75.

Figure 5 :
Figure 5: Spectral parameter evolution of Outburst A (Right panels) and B (left panels).Only the parameters of the observations with  2 /DoF< 1.5.Different colours are used to show the evolution of parameters for different f.The marker shapes are kept identical to figure 1.The inverted triangles in panels (d) and (i) indicate the observations which peg of Γ=1.The inset plot in panel (h) shows the disc norm evolution of the fainter observations in outburst B.The parameter values for the outbursts for  = 0.75 are noted in tables A4 and A5.

Table A1 :
NICER observation log for outburst A. Exposures are net exposure after the standard filtering.

Table A2 :
NICER observation log for outburst B. Exposures are net exposure after the standard filtering.

Table A3 :
Log of all the burst intervals observed in the NICER monitoring of 2019 and 2020 outburst of the source.The time mentioned here is in spacecraft time.The mean count rate is the mean of the count rate before the outburst.The pean count rate is the peak count rate of the burst.For the observation 3050340116, the decay of a burst is observed, so the peak count rate is uncertain.

Table A4 :
Outburst A (2019) spectral parameter table for f =0.75 assuming the models: Model A: tbabs*(gauss+thcomp⊗diskbb) and Model B: tbabs*(gauss+gauss+thcomp⊗diskbb).Errors listed here are 1 errors and they are not listed for the observations with  2 /DoF > 1.5 and for parameter with error of the order 0.001.