A strongly changing accretion morphology during the outburst decay of the neutron star X-ray binary 4U 1608-52

It is commonly assumed that the properties and geometry of the accretion flow in transient low-mass X-ray binaries (LMXBs) significantly change when the X-ray luminosity decays below $\sim 10^{-2}$ of the Eddington limit ($L_{\rm Edd}$). However, there are few observational cases where the evolution of the accretion flow is tracked in a single X-ray binary over a wide dynamic range. In this work, we use NuSTAR and NICER observations obtained during the 2018 accretion outburst of the neutron star LMXB 4U 1608-52, to study changes in the reflection spectrum. We find that the broad Fe-K$\alpha$ line and Compton hump, clearly seen during the peak of the outburst when the X-ray luminosity is $\sim 10^{37}$ erg/s ($\sim 0.05$ $L_{\rm Edd}$), disappear during the decay of the outburst when the source luminosity drops to $\sim 4.5 \times 10^{35}$ erg/s ($\sim 0.002$ $L_{\rm Edd}$). We show that this non-detection of the reflection features cannot be explained by the lower signal-to-noise at lower flux, but is instead caused by physical changes in the accretion flow. Simulating synthetic NuSTAR observations on a grid of inner disk radius, disk ionisation, and reflection fraction, we find that the disappearance of the reflection features can be explained by either increased disk ionisation ($\log \xi \geq 4.1$) or a much decreased reflection fraction. A changing disk truncation alone, however, cannot account for the lack of reprocessed Fe-K$\alpha$ emission. The required increase in ionisation parameter could occur if the inner accretion flow evaporates from a thin disk into a geometrically thicker flow, such as the commonly assumed formation of an radiatively inefficient accretion flow at lower mass accretion rates.


INTRODUCTION
Low-mass X-ray binaries (LMXBs), in which a neutron star or a black hole attracts gas from a low-mass companion star, are prime tools to study the physics of accretion. The X-ray luminosity (LX) of LMXBs scales with the rate at which matter is accreted by the compact object. The population of transient sources alternates outbursts of accretion with quiescent episodes and can therefore be studied over a wide range of X-ray luminosities; from the Eddington limit, L Edd 10 38 (M/M ) erg s −1 , down to a very small fraction of it (LX∼10 −8 L Edd ).
It is observationally established that at high X-ray luminosities e-mail: A.J.vandenEijnden@uva.nl ( 10 −2 L Edd ), the accreted gas spirals in a thin Keplerian disk that typically extends close to the compact object (e.g., Done et al. 2007, for a review). At very low accretion rates ( 10 −4 − 10 −5 L Edd ), on the other hand, the accretion disk is likely truncated further away from the black hole or the neutron star (e.g.,, Zdziarski et al. 1999;Esin et al. 2001;Kalemci et al. 2004;Tomsick et al. 2004). Nevertheless, it is currently unclear at what accretion luminosity the disk begins to move away, how this truncation proceeds with changing luminosity, and if this process differs for neutron stars and black holes. Standard accretion theory suggests that with decreasing accretion rates, an increasing part of the inner disk evaporates into a radiatively inefficient hot flow (e.g., Narayan & Yi 1994;Esin et al. 1997). This classical disk truncation mechanism is expected to be more efficient for black holes, because for neutron stars the soft photons from the stellar surface can cool the hot flow and make it condense back on the disk (Narayan & Yi 1995). In neutron star LMXBs, a second mechanism that can lead to disk truncation may be at play; the stellar magnetic field may push the inner disk out (e.g., Ibragimov & Poutanen 2009;Miller et al. 2011;Degenaar et al. 2014Degenaar et al. , 2017van den Eijnden et al. 2017). One may expect this effect to become larger with decreasing accretion rate, and hence lower pressure exerted by the disk. However, there is little observational data to test this idea in a single source (van den Eijnden et al. 2018). The interaction between the stellar magnetic field and accretion flow plays an important role in a number of key processes, including the production of outflows (e.g., Romanova et al. 2009;Deller et al. 2015;van den Eijnden et al. 2018), and the spin evolution of neutron stars (e.g., Ghosh & Lamb 1978;Alpar et al. 1982;Wijnands & van der Klis 1998;D'Angelo 2017;Bhattacharyya & Chakrabarty 2017). Yet, while the interaction between the accretion disk and the magnetosphere is clearly important for the evolution and appearance of neutron stars, it remains poorly understood.
X-ray reflection modelling is a powerful tool to constrain the geometry of accretion flows in LMXBs. This X-ray emissionoriginating from close to the accretor and then reprocessed or 'reflected' toward the observer by the accretion disk -manifests itself most prominently as an Fe-Kα emission line at ∼6-7 keV and the Compton hump at ∼20-40 keV. The shape of these features is modified by Doppler and gravitational redshift effects as the gas in the disk moves in high-velocity Keplerian orbits inside the gravitational well of the compact accretor. In particular, reflection modelling allows for a measure of the inner radial extent of the accretion disk, Rin. The reflection spectrum thus encodes information about the morphology of the accretion flow (e.g., Fabian & Ross 2010;Dauser et al. 2016, for reviews). In addition, the reflected spectrum is affected by the internal properties of the disk, such as level of ionisation and composition.
Comparing the reflection spectrum at different fractions of the Eddington accretion rate in the same source can reveal changes in the accretion flow structure and properties as the X-ray luminosity decays. An example of a source where the reflection spectrum has been tracked over a range in X-ray luminosity is the black hole LMXB GX 339−4 (e.g., García et al. 2015). In this source, the relation between X-ray luminosity and disk inner radius has been determined in detail, showing a general trend of decreasing inner radius with Eddington fraction (Miller et al. 2006;Reis et al. 2008;Tomsick et al. 2009;Petrucci et al. 2014;García et al. 2015;Wang-Ji et al. 2018;García et al. 2019). This large set of measurements is possible due to high duty cycle of GX 339−4, yielding over twenty measurements of the reflection spectrum in the past fifteen years.
However, such systematic studies have to date been rare, especially in neutron stars, given the deep exposures required to detect reflection or otherwise obtain constraining upper limits at low luminosities (e.g., LX< 10 −2 L Edd ). In addition, the latter luminosity regime is often traversed quickly and unpredictably, complicating the scheduling of observational campaigns in sources that show fewer outbursts. Finally, for neutron stars, a mere measurement of the inner radius does not necessarily reveal the driving factor of the disk truncation: the formation of a radiatively inefficient hot flow, truncation by the magnetic field, or possibly both (Cackett et al. 2008;Degenaar et al. 2017). So while NuSTAR reflection spectra have been measured for many individual neutron star LMXBs (see Ludlam et al. 2019, for a recent overview), differences between their neutron star properties mean variations in the reflection spectrum can not be directly attributed to accretion rate. For such an inference, a reflection study of a single source at different accretion rates is essential.
Given these challenges for neutron star LMXBs, 4U 1608−52 (Grindlay & Gursky 1976;Tananbaum et al. 1976) is a promising source to track the reflection spectrum and accretion geometry at different outburst stages, down to LX< 10 −2 L Edd . Firstly, a NuS-TAR observation performed during its 2014 outburst, when it was accreting at ∼ 10 −2 L Edd (∼ 10 36 erg s −1 ), revealed very strong reflection features that allowed to measure an inner disk radius of 8.5 ± 1.5 rg(rg= GM/c 2 Degenaar et al. 2015). The strength of its reflection is a very important advantage over other frequently active neutron star LMXBs, such as Aql X-1 and SAX J1808.4-3658, for which the reflection spectrum is much weaker (e.g., Cackett et al. 2008;Ludlam et al. 2017). Secondly, 4U 1608−52 is relatively close (∼ 3 − 4.5 kpc; Galloway et al. 2008;Poutanen et al. 2014), so that a large flux is obtained even when it accretes at a low Eddington luminosity. Thirdly, 4U 1608−52 is among the most active transient neutron star LMXBs, going into outburst about once every 2-4 yr (e.g., Lochner & Roussel-Dupre 1994;Šimon 2004), easing the scheduling of observations. Finally, the outburst profiles are often similar, making it feasible to plan observations at a specific X-ray luminosity.
In this work, we report on our efforts to use reflection spectroscopy to constrain the changing accreting morphology in the neutron star LMXB 4U 1608−52 in its 2018 outburst, using the NuSTAR, NICER, and Swift X-ray observatories.

Set up of the observation program
The MAXI transient alert system (Negoro et al. 2016) signalled activity from 4U 1608−52 on 19 June 2018 June 19 (MJD 58288; nova ID 8288849999). When MAXI monitoring indicated that 4U 1608−52 had evolved beyond its outburst peak, we started to monitor the source with the X-ray telescope (XRT) aboard the Neil Gehrels Swift observatory (Swift; Gehrels et al. 2004). This Swift monitoring was set up to facilitate triggering a NuSTAR observation (PI Degenaar), aiming to catch the source at a luminosity of ∼ 10 35 erg s −1 (∼ 10 −3 L Edd ), i.e. a factor ∼10 lower than the 2014 observation of the source. In addition to our Swift and NuS-TAR program, we used observations obtained with NICER (Gendreau et al. 2016) to investigate if, and how, the profile of the Fe-Kα line changed along the course of the 2018 outburst.

Spectral modelling
All spectral fits reported in this work were performed using XSPEC (v. 12.10.1;Arnaud 1996). To account for interstellar absorption along our line of sight, we used the model TBABS, with VERN cross-sections (Verner et al. 1996) and WILM abundances (Wilms et al. 2000), in our spectral modelling. Reported errors for spectral parameters indicate 1σ confidence levels. To calculate luminosities from measured fluxes, we adopted a distance of 3.6 kpc.

Swift
To put the NuSTAR and NICER observations in perspective, we use 22 Swift/XRT observations (obsID 00010741001-26) to illustrate the source's evolution along the 2018 outburst. These observations were obtained between 3 July and 24 October 2018 (MJD 58302 and 58415) and had typical exposure times of ∼1 ks. We followed standard data reduction procedures, using XSELECT to extract spectra and XRTMKARF to create arf files, while the rmf file was taken from the CALDB. All XRT spectra were fitted to a simple absorbed powerlaw model (TBABS*PEGPWLW in XSPEC). In these fits, we fixed the hydrogen column density at NH = 1.58 × 10 22 cm −2 (as determined from a simultaneous fit to all XRT spectra where NH was tied). The unabsorbed 0.5-10 keV flux was subsequently determined using the convolution model CFLUX. While in the main paper we only show the XRT flux light curve (cf. Figure 1), details on these spectral fits can be found in Appendix A.

NuSTAR
4U 1608−52 was observed for ∼53 ks with NuSTAR on July 12-14 2018 (obsID 80401321002), providing data in the 3-79 keV energy band. We processed the data with NUPIPELINE, incorporated in HEASOFT (v. 6.23), using SAAMODE=OPTIMIZED. We extracted light curves and spectra for the two co-aligned grazing incidence hard X-ray imaging telescopes, with focal plane modules (FPM) A and B, using NUPRODUCTS. A circular region with a radius of 60 was used to extract source events and a void region of the same size, placed on the same detector, was used to obtain background events.
All spectral data were grouped into bins with a minimum of 20 photons using GRPPHA. 4U 1608−52 was detected significantly above the background in the entire NuSTAR bandpass. Furthermore, when fitting the data of the two mirrors simultaneously with a freely-floating constant factor in between, we found the difference between the FPMA and B to be <1% for this data set. Therefore, we summed the two spectra using ADDASCASPEC and in the following report only on the results obtained using this summed spectrum. In this work, we compare this 2018 data set to the NuSTAR observation performed in 2014. For details on the 2014 observation and data reduction, we refer to Degenaar et al. (2015).

NICER
We used 12 NICER observations (obsIDs 1050070125-1050070136) obtained between 21 June and 9 July 2018, and having typical exposure times of 0.1-2.5 ks. No Type-I X-ray bursts occurred during the NICER monitoring. The data were reduced using NICERDAS 2018-10-07 V005 with the standard NICER filtering. Additionally, we created GTIs to select a SUN ANGLE > 60 • , COR SAX > 4, and KP < 5 via the NIMAKETIME tool. The latter two corrections define the cut-off magnetic rigidity and the space weather conditions. We applied the GTIs to the data using the NIEXTRACT-EVENTS command. A light curve was then extracted by loading the cleaned event files into XSELECT. We separate the NICER monitoring of the 2018 outburst into five intervals shown in Figure 2. A spectrum was were extracted from each interval.
Because NICER does not have imaging capabilities, backgrounds were generated from periodic observations of RXTE blank sky field locations (Jahoda et al. 2006). We extracted a timeaveraged background from 51 observations of blank sky field 5 that were reduced in the same manner as 4U 1608−52. The background spectrum was subtracted from the source spectra prior to normalizing the data to the Crab Nebula (see Ludlam et al. 2018, for more details). This accounts for instrumental uncertainties within the data as NICER's calibration continues to be refined. We note  . NICER count rate light curve of the 2018 outburst of 4U 1608−52 (0.5-6.8 keV, binned by 128 s per observation). The dashed red lines indicates the different phases of the outburst for which we extracted X-ray spectra (presented in Figure 4). that the source was significantly detected above the background in all used observations.

RESULTS
The Swift/XRT light curve in Figure 1 (top) shows how 4U 1608−52 initially decayed steadily over a timescale of ∼ 10 days, after which the 0.5-10 keV flux rose by a factor of a few, remaining steady for ∼20 days until finally returning to quiescence ∼2 months after the outburst peak. We note that the 0.5-10 keV flux detected in the last set of Swift observations is around the quiescent level reported from Chandra observations (Marino et al. 2018). The 2018 NuSTAR observation was obtained exactly in the valley of the initial decay of the outburst (red bar in Figure 1 top). At this time, the 0.5-10 keV unabsorbed flux measured by Swift (Obsid 00010741006) was (2.9 ± 0.15) × 10 −10 erg cm −2 s −1 (see Appendix A for details on the Swift spectral fits). This translates into a luminosity of (4.5 ± 0.2) × 10 35 erg s −1 which is a factor ∼10 lower, as intended, than during the NuSTAR observation that was obtained during the 2014 outburst (Degenaar et al. 2015).
4U 1608−52 shows regular outbursts, that broadly fall into two categories: relatively faint and short outbursts, and significantly brighter and longer ones. As the two NuSTAR observations were taken in two outbursts, we compare their profiles in Figure 3. While slight differences exist in peak count rate and duration, the comparison with a 2016 outburst (shown in grey) shows that both outbursts clearly fall in the faint/short category. Therefore, in the remainder of this paper, we will assume the accretion geometry follows a similar evolution throughout both outburst decays. The comparison of the 2018 NICER and the 2014 NuSTAR observations support this assumption, as shown in the next Section.

The disappearance of reflection during the decay
During the 2014 outburst peak, NuSTAR revealed strong reflection features in 4U 1608−52, detecting both a broad Fe-Kα line and a Compton hump (Degenaar et al. 2015). With NICER monitoring, we can assess whether such features were also initially present during the 2018 outburst. For this purpose, we initially performed a quick-look analysis of both the NICER spectra of the first four epochs (rise, peak, decay start, middle; see Figure 2) and the archival NuSTAR spectrum using a simple continuum model, TBABS×(BBODYRAD+NTHCOMP). To highlight the reflection features and only fit the continuum, the 5-7 keV and 15-30 keV ranges were ignored. In the left panel of Figure 4, we show the datato-model ratio for these spectral fits in the 1-79 keV energy range.
We can see clearly that during the rise, peak, and early/middle decay of the outburst, a strong, broad Fe-Kα line is present with a stable shape similar to the line profile seen in the 2014 NuSTAR observation.
Repeating the same quick-look exercise with the spectrum of the final NICER epoch, and the 2018 NuSTAR observation (and using a more simple TBABS×NTHCOMP model due to the lower flux), we can see that both the Fe-Kα line and the Compton hump disappear. We show this evolution of the reflection spectrum in the right panel of Figure 4.
The disappearance of the Fe-Kα line, combined with its remarkable prior stability, as observed by NICER is suggestive of a distinct change in the accretion flow geometry and/or intrinsic properties on a days time scale. We choose, however, not to fit physical reflection models to the NICER data, for two reasons: firstly, NICER spectral calibration continues, at the time of writing, to be Table 1. Fits to the 2018 NuSTAR observation of 4U 1608−52. The continuum model consists of a TBABS×(BBODYRAD+ NTHCOMP) model, where the blackbody temperature kT BB is linked between the thermal and comptonized component. In the reflection model, the NTHCOMP component is replaced by a RELXILL model, which contains both the Comptonized continuum and the broad-band reflection spectrum. In all reflection modelling, we fix the emissivity index, outer disk radius and cutoff energy to their default values, and set the spin parameter to a = 0.29, corresponding to the known spin frequency. All uncertainties are given at the 1σ level.

Parameters 2018 Continuum
Reflection (866) refined. Secondly, and importantly, the soft NICER bandpass does not include the Compton hump clearly seen in the 2014 NuSTAR spectrum (Figure 4, left panel). As this difference in coverage limits the ability to constrain the parameters of the reflection spectrum, we instead focus on the deep 2018 NuSTAR observation in the remainder of this work. For our detailed analysis, we start by fitting the 2018 NuS-TAR spectrum with a simple phenomenological continuum model, TBABS*(BBODYRAD+NTHCOMP), where we tie the temperature of the NTHCOMP thermal input spectrum to that of the BBODYRAD component. We fit the data in the full 3-79 keV NuSTAR bandpass. This simple model, not containing any reflection features, provides an adequate description of the data, with a reduced χ 2 of 0.89 for 870 degrees of freedom. In Figure 5, we show the spectrum and this continuum model, while the fitted spectral parameters are listed in Table 1. The data-to-model ratio shown in the bottom panel of Figure 5 confirms the lack of detected reflection features seen in the final NICER epoch.
While no reflection features appear to be present, we also attempt a fit with a full reflection model, TBABS*(BBODYRAD+RELXILL). This model is based on that used by Degenaar et al. (2015) to fit the 2014 NuSTAR observation. The RELXILL model (version v1.0.4) includes both a power law continuum component and a fully relativistic, broad-band reflection model. The parameters that set the reflected spectrum are the emissivity profile, the neutron star spin, the system's inclination, inner and outer disk radius, ionisation parameter of the disk, iron abundance, cutoff energy, and reflection fraction (defined as the ratio of reflected and direct flux between 20 and 40 keV). The illuminating flux is characterised by the power law index. Following Degenaar et al. (2015), we set the spin parameter to 0.29, and freeze the outer disk radius to 400 gravitational radii rg.  Additionally, we freeze the emissivity profile and the cutoff energy to their default values. Fitting this model does not yield a better fit than the continuum model, with a similar reduced χ 2 of 0.90 for 866 degrees of freedom (i.e. four additional free parameters). Unsurprisingly, given the lack of reflection features, the RELXILL parameters are also poorly constrained. All parameters and fluxes with their uncertainties can be found in Table 1.

Could reflection be hiding in the noise?
No reflection features are detected during the last NICER epoch and the 2018 NuSTAR observation. However, the lower flux raises the question whether the lack of detected reflection might arise from a lower signal-to-noise instead of physical changes in the source. To assess this question, we simulated synthetic NuSTAR observations: as a first step, we fit the TBABS*(BBODYRAD+RELXILL) model to the 2014 NuSTAR spectrum, reproducing the model fit that was reported in Degenaar et al. (2015). We then scale the normalisations to the flux during the 2018 NuSTAR observations, and use the XSPEC FAKEIT command to simulate a spectrum with the same exposure time as the 2018 observation. We then fit this simulated spectrum with both the reflection and the continuum models, and compare those fits to test whether the 2014-level reflection would have been detected at the lower flux and 2018 exposure time.
The fit with the (correct) reflection model yields χ 2 ν = 1343.6/1283 = 1.05, while the continuum model returns a worse fit with χ 2 ν = 1550.8/1287 = 1.20. An f-test comparison of these models yields a chance probability of such an improvement by the reflection model of p = 10 −38 . However, this does not take the change in continuum shape at lower flux into account. Therefore, we repeated this test by simulating a second fake spectrum from the reflection model, now using the continuum parameters (power law index, RELXILL normalisation, and blackbody temperature) as found for the 2018 NuSTAR data. For this second synthetic observation, the reflection and continuum models yielded χ 2 ν = 1250.3/1283 = 0.97 and χ 2 ν = 1305.4/1287 = 1.01 respectively. While the difference between the two fits is smaller, the (correct) reflection model still fits significantly better at an f-test probability of 2 × 10 −11 . In other words, both simulations show that the 2018 NuSTAR observational setup should have detected reflection features as present during the 2014 (and 2018, as shown by NICER) outburst peaks.
Therefore, we conclude that the disappearance of the reflection features at lower flux is not only due to the decrease in signalto-noise, but that physical changes in the accretion flow must also contribute.

Constraining the accretion geometry in the absence of reflection
Without a detected reflection signal during the 2018 NuSTAR observation, we cannot directly fit reflection parameters to determine how the accretion flow properties have changed. However, we can make use of the fact that reflection features are not detected in another way: we can instead simulate synthetic spectra with the same observational setup and flux as the 2018 observation, using different combinations of reflection parameters. For each set of reflection parameters, we can then assess whether a reflection signal would have been detected and hence constrain the changes in the accretion flow during the outburst decay, despite the absence of fittable reflection features.
Considering the relevant parameters of the RELXILL model, only three parameters can be expected to change during an outburst decay: the inner disk radius, the disk ionisation, and the reflection fraction. The remaining parameters are either continuum parameters with measured changes (power law index and normalisation), flux-independent (spin, inclination, and iron abundance), or frozen to their default values in our modelling (emissivity profile, cutoff energy). Therefore, we set up a grid in inner disk radius, ionisation, and reflection fraction, and for each combination simulate a NuS-TAR spectra from the TBABS*(BBODYRAD+RELXILL) model with the continuum parameters from the 2018 NuSTAR spectral shape. For each simulated spectrum at each gridpoint, we then fitted three models: the full (underlying) reflection model, a continuum-only model TBABS*(BBODYRAD+NTHCOMP), and a continuum model with an additional Gaussian Fe-Kα line. The energy of the latter Gaussian was constrained to lie between 6.4 and 6.96 keV, and its equivalent width provides a diagnostic of the strength of the Fe-Kα line in the simulated data. A second diagnostic of the strength of reflection is given by the reduced χ 2 of the continuum-only model, which should approach ∼ 1.0 as the reflection becomes undetectable.
We performed our simulations for the following parameter grid: inner disk radius between Rin = 6 and Rin = 120 rg, ionisation parameter between log ξ = 0.5 and log ξ = 4.7, and reflection fraction between F ref = 0.0 and F ref = 3.0. We first ran a lowresolution simulation to determine regions of interest in the probed parameter space, using 11 points in inner disk radius, 6 in ionisation parameter, and 5 in reflection fraction. Then, after concluding that the inner disk radius had little effect, while the reflection fraction and both high and low ionisation parameters did (as discussed below), we reran the simulations twice: both times with a lower resolution in inner disk radius (6 points) and a higher resolution in reflection fraction (9 points). These two added simulations used a higher resolution in ionisation parameters (6 points) to zoom in on high (log ξ = 3.2 − 4.7) and low (log ξ = 0.5 − 1.5) values, respectively.
We show representative results of this grid search in Figures 6  and 7: each panel shows one of the two aforementioned diagnostics of reflection strength -reduced χ 2 of the continuum-only model and Gaussian iron line equivalent width -at a given value of one parameter, as a function of the two other parameters. In Figure 6, we show the χ 2 ν of the continuum-only model at low (log ξ = 1.2; left) and high (log ξ = 4.5; right) ionisation. From this figure, we can draw two conclusions: first, the lack of vertical changes reveals that the value of the inner disk radius does not greatly affect the detectability of the reflection signatures. Secondly, the reflection signatures are always weak at high values of the ionisation parame-ter, while for lower ionisation they are weak only for low reflection fractions.
The same results can be inferred from Figure 7. Here, the panels show the two reflection strength diagnostics for approximately the same disk inner radius. The continuum-only model fits best, again implying weak to no reflection, at low reflection fraction or high ionisation. The Gaussian Fe-Kα line's equivalent width (bottom panel) confirms this trend, as it decreases towards those same regions of parameter space. Finally, we stress that the χ 2 ν values of the continuum model shown in Figures 4 and 5 nowhere reach below the χ 2 ν found when fitting the real 2018 observation with the continuum model (i.e. χ 2 ν = 0.89). Finally, we show a one-dimensional representation of our simulation output in Figure 8. Here, we show the reduced χ 2 for the simulated data fitted with the continuum-only model, as a function of disk ionisation. The lines of different colour and width represent different reflection fractions. Given the lack of any systematic effect as the inner disk radius varies (see above), we pick a representative inner radius of 63 rg and plot the results for that value. Figure 8 confirms the inferences made above: the reflection signatures always become undetectable towards very high ionisation ( log ξ 4.1-4.3, depending on the reflection fraction). In addition, the non-detection of reflection features can be explained by a low reflection fraction (F ref ≈ 0) independently of ionisation.

DISCUSSION
We have presented a systematic analysis of the X-ray reflection spectrum during the outburst decay of the NS LMXB 4U 1608−52. Comparing Swift, NICER, and NuSTAR spectra obtained over two similar outbursts, we find the following results: (i) Reflection features, such as a broad Fe-Kα line and Compton hump, are systematically present at high mass accretion rates without changing shape, but disappear completely within days towards the end of the outburst. (ii) This non-detection of reflection features at low flux is not due to a lower signal-to-noise, but can instead be explained by either a high ionisation (log ξ 4.1) or low reflection fraction (F ref ≈ 0) (iii) the lack of reflection features does not result from changes in inner disk radius alone, although we stress that we cannot directly constrain the value of this radius. Here, we will discuss these results in a physical context and compare it to other LMXB systems.
From our simulations, we conclude that the disappearance of measurable reflection features can be explained by either a very low reflection fraction, or a high disk ionisation, approaching the grid edge of the RELXILL model. The reflection fraction in the RELXILL model is defined as the ratio of the reflected and direct flux in the 20-40 keV range (Dauser et al. 2016). While for some flavours of reflection models, reflection fractions have been calculated self-consistently assuming given illuminating geometries, in the RELXILL model used in this work no specific geometry is assumed (Dauser et al. 2016). Therefore, while it can explain the 2018 NuSTAR spectrum, a very low reflection fraction does not directly map onto a physical picture.
Therefore, in the remainder of the discussion, we focus more on the second explanation of a high disk ionisation, as this provides more direct physical insight into the evolution of the system. The ionisation parameter scales with the ionising X-ray flux divided by particle density, i.e. ξ = LX,ionising/(nr 2 ). Comparing the two NuSTAR spectra of 4U 1608−52, we find that the spectrum is harder when no reflection is detected. In other words, a higher fraction of flux comes out at high (ionising) energies. However, the total flux also decays, which counteracts this effect. More likely, an increase in ionisation parameter could be explained by a decrease in particle number density as the source decays to lower X-ray luminosities.
What might explain such a decrease in the number density of particles in the disk? As X-ray binaries decay towards quiescence, their accretion flows are expected to change morphology. For instance, the inner accretion flow might evaporate into a radiatively inefficient hot flow (for instance an advection-dominated accretion flow; Narayan & Yi 1994, 1995. During the same transition, the inner radius of the thin disk might move outwards, truncating the Keplerian disk (e.g Zdziarski et al. 1999;Esin et al. 2001;Kalemci et al. 2004;Tomsick et al. 2004). As the inner accretion flow evaporates, the density of the gas will decrease. This transition could therefore explain the increase in ionisation parameter and resulting disappearance of the reflection signatures observed in 4U 1608−52. The NICER observations at higher mass accretion rates show a remarkably stable Fe-Kα line profile, which suggest that in this scenario, the disk initially remains stable. Subsequently, the disk changes morphology within days, as it evaporates into a radiatively inefficient hot flow below ∼ 10 36 erg s −1 (i.e. below ∼ 10 −2 L Edd ). Due to gaps in the X-ray monitoring and required grouping of the NICER observations, the exact transition X-ray luminosity cannot be measured precisely.
The possible transition to a radiatively inefficient hot flow before the NuSTAR observation in the 2018 outburst can be compared to the behaviour of other neutron star LMXBs. For instance, the continuum blackbody and comptonization parameters during the 2018 NuSTAR observation (e.g., Table 1) are very similar to those observed for Aql X-1 with Suzaku at a similar luminosity level (Sakurai et al. 2014). At this luminosity, no reflection features are detected in Aql X-1 either, which might similarly arise from the evaporation into a hot flow. It should be noted, however, that reflection features in Aql X-1 are also weak at higher luminosity and that the disk is likely always truncated away from the neutron star (King et al. 2016;Ludlam et al. 2017). Therefore, the inner disk geometry might not be exactly comparable and the reflection spectrum could simply be undetectable due to the low flux. This comparison with Aql X-1 highlights the exceptional combination of strong reflection and low distance presented by 4U 1608−52, which allows us to conclude that the lack of reflection does not result only from limited signal-to-noise.
An additional comparison can be made with very-faint X-ray binaries (VFXBs), which never reach high accretion luminosities, as they either persistently accrete between ∼ 10 34 to 10 36 erg s −1 or exhibit outbursts that peak below 10 36 erg s −1 . Of the few VFXBs with detailed spectral observations, most do not show reflection features (e.g., Armas Padilla et al. 2011Padilla et al. , 2013aLotti et al. 2016;Sanna et al. 2018a,b). Moreover, the transient black hole LMXB Swift J1357.2-0933, which does not become brighter than ∼ 10 35 erg s −1 during its accretion outbursts, does not show an Fe-Kα line either (Armas Padilla et al. 2014;Beri et al. 2019). Being in the same X-ray luminosity range as 4U 1608−52 during the 2018 NuSTAR observation, the lack of reflection in these VFXBs might also be explained by a radiatively inefficient accretion flow with a high ionisation parameter.
Interestingly, weak reflection features were detected with NuSTAR, XMM-Newton and Chandra in the (persistently accreting) neutron star VFXB IGR J17062−6143. Reflection modelling of those data implied that the disk could be truncated at ∼ 100 rg van den Eijnden et al. 2018). In this source, and generally in neutron star X-ray binaries showing truncated disks, it is difficult to distinguish between two possible explanations for the disk truncation: a relatively strong magnetic field, or the evaporation of the inner thin disk. However, our results for 4U 1608−52 might provide some insight: in IGR J17062−6143, reflection is detected despite a strongly truncated accretion flow. Based on our results on 4U 1608−52, we can envision three scenarios: firstly, that the disk in IGR J17062−6143 is also evaporated but has not reached high enough ionisation levels to yield reflection features undetectable; secondly that the disk truncation is instead  Figure 6, now showing the reduced χ 2 as a function of reflection fraction and disk ionisation parameter at a given inner disk radius. Reflection features become undetectable (i.e. χ 2 ν ≈ 1.0) at low reflection fraction and very high disk ionisation. Bottom: the equivalent width of a Gaussian iron line added to the simple continuum model, fitted to the simulated spectra. The axes are the same as the top two panels. A low equivalent width corresponds to undetectable reflection features, confirming the trends visible in the reduced χ 2 diagnostic. caused by a different mechanism, such as the neutron star magnetosphere truncated an otherwise relatively thin disk. The presence of X-ray pulsations, making IGR J17062−6143 an accreting millisecond X-ray pulsar (Strohmayer & Keek 2017;Strohmayer et al. 2018), plausibly fits within this picture, as matter is channeled to the magnetic poles; or thirdly, that the disappearance of reflection features in 4U 1608−52 is caused by a low reflection fraction instead, which doesn't occur in IGR J1706-2−6143.
Here, we briefly consider the effect of the continuum shape on the detectability of reflection features. Analysing four Suzaku observations of 4U 1608−52 throughout its 2010 outburst, Armas Padilla et al. (2017) find no evidence for reflection at any moment, placing limits on the equivalent width of the Fe-Kα line of < 14 eV. These spectra were taken in the range of ∼ (1−20)×10 36 erg s −1 , a factor two to forty brighter than the 2018 NuSTAR spectrum. Armas Padilla et al. (2017) note that this lack of reflection could result from either an intrinsic lack of reflection features, or burying by the continuum. We find that at lower luminosity it can be a combination of both: while reflection should have been seen in the 2018 NuSTAR observation if it were present (Section 3.2), it becomes undetectable at lower ionisation parameter for the 2018 continuum shape. However, the comparison between NuSTAR and Suzaku is more accurate using the 2014 observation, which overlaps in flux with the Suzaku range (Degenaar et al. 2015). Given the strong reflection in the NuSTAR spectrum, and in the NICER monitoring around the 2018 outburst peak, a low signal-to-noise seems the most likely explanation of the lack of reflection detected by Armas Padilla et al. (2017).
If reflection had been detected in 4U 1608−52 during the 2018 NuSTAR observation, we could have measured the disk inner radius and compared it with the earlier NuSTAR measurement. However, in the absence of reflection, we do not obtain any constraints on the inner disk radius. While it could very well change in tandem with the changing ionisation state, our simulations show that varying the inner disk radius alone cannot explain the disappearance of reflection features. Therefore, we cannot place any constraints on the evolution or value of rgin the 2018 observation.
Finally, in the few black hole LMXBs observed at similarly low X-ray luminosities, the inner disk appears to be truncated: in GX 339−4, the inner disk radius appears to move from a few gravitational radii at high accretion rate to 35 rg at ∼ 10 −3 L Edd (Tomsick et al. 2009;García et al. 2019). A NuSTAR observation of GRS 1739−278 places its inner disk at ∼ 15 − 35 rg when accreting at ∼ 10 −2 L Edd (Fürst et al. 2016). Interestingly, in these sources, the Fe-Kα line is clearly detected. In the scenario proposed above, where we assume ionisation changes caused by a transition to a low-density, radiatively inefficient accretion flow drive the iron line disappearance, instead of low reflection fraction, this FeKαline is not expected. Therefore there might be a difference in the evaporation of the disk as the thin disk recedes in these two black hole systems, when compared to 4U 1608−52. Inner disk radius = 63.0R g Figure 8. The reduced χ 2 of the continuum model fitted to the synthetic NuSTAR spectra, as a function of disk ionisation parameter, shown for low (left) and high (right) ionisation. Lines of different colors and thickness indicate different reflection fractions. As expected, for zero reflection fraction in the synthetic data, the continuum model always provide a good description. However, for non-zero reflection fractions, the reduced χ 2 only approaches 1 for the highest ionisation states (log ξ 4.1).