Dissecting a Disk-instability Outburst in a Symbiotic Star: NuSTAR and Swift Observations of T Coronae Borealis during the Rise to the "Superactive" State

We observed T CrB with the NuSTAR satellite on 2015 September 23 for a total of 79.8 ks (hereafter NuSTAR₂₀₁₅), well into the rising phase of the latest optical brightening event. We used the NUPRODUCTS v0.3.0 software to extract both spectral and light-curve products, resulting in a spectrum and light curve for each FPMA and FPMB units. Source events were extracted from a 70'' circular region centered on α = 15^h59^m30160 and δ = +25°55'1259, while background events were extracted from an equally sized circular region off-source. The resulting spectra were binned to have a minimum of 25 counts per bin.

The Swift satellite observed T CrB twice during the NuSTAR₂₀₁₅ observation (ObsId 00081659001 and 00081659002, hereafter Swift${}_{2015}^{1}$) and a week later (ObsID 00045776003 and 00045776004, hereafter Swift${}_{2015}^{2}$). The total duration of the combined exposures of Swift${}_{2015}^{1}$ was 9853 s and of Swift${}_{2015}^{2}$ is 2953 s. Data products were created by first combining the event files from the two observations and then extracting source events from a circular region of radius 20 pixels centered on the SIMBAD coordinates of T CrB, and background events from an off-source circular region with a radius of 80 pixels. The ancillary response files were created using the xrtmkarf tool.

During these same observations, we also obtained UVOT images in the U (λ3465 Å, FWHM = 785 Å) and UVW2 (λ1938 Å, FWHM = 657 Å) filters. T CrB was bright enough to saturate the UVOT detector, resulting in severe coincidence losses that meant that source brightness could not be estimated using the standard photometry tools in HEAsoft. Instead, we utilized the readout streak in the images to estimate the source brightness (Page et al. 2013).

The NuSTAR₂₀₁₅ spectra are well suited to looking for evidence of reflection of X-rays from the WD surface. We explored possible models for T CrB by jointly fitting the two NuSTAR₂₀₁₅ modules and the Swift${}_{2015}^{1}$/XRT data (Figure 2). The short exposure time and low number of counts detected during the Swift${}_{2015}^{2}$ observation did not provide a significant improvement to the joint fit. We fit the NuSTAR₂₀₁₅ spectra in the range 3–78 keV, and the Swift${}_{2015}^{1}$/XRT spectrum in the range 0.3–10 keV. For models that include reflection we used the reflect model in XSpec, which calculates the reflection of X-rays from neutral material using the method of Magdziarz & Zdziarski (1995). The angle between the line of sight and the reflecting surface normal (θ) is a key parameter in the reflection model, but Magdziarz & Zdziarski (1995) showed that cos(θ) = 0.45 provided a reasonable approximation for a wide range of θ (18°–76°; see their Figure 5). For the specific case of the WD surface reflecting X-rays from an equatorial boundary layer, θ will be close to 90 if the binary is very close to exactly face on (binary inclination i ∼ 0), while for an edge-on system, what we need is the average over azimuth changing from θ = 90° to 0° back to 90°, with real systems being intermediate cases. Done & Osborne (1997) evaluated this average for SS Cyg (assumed i ∼ 37°) and estimated the average θ to be close to 60° (so cos(θ) ∼ 0.5); for T CrB, with a higher binary inclination, this will be somewhat lower. However, given the relatively slow dependence of reflection on cos(θ) in this regime, we have fixed it to the Xspec default value of 0.45 in our analysis. To use the reflection model, we must extend the energy response over which XSpec evaluates the model. For T CrB, we need to do this well beyond the NuSTAR₂₀₁₅ response. We chose 200 keV as the high energy cutoff, and model in 400 linear bins using the command energies extend,200.0,400,lin in XSpec.

Zoom In Zoom Out Reset image size

Given our knowledge of the current superactive state, which is most likely due to an enhancement of the accretion rate through the disk and an increase of the optical depth of the boundary layer (Luna et al. 2018a), we considered that both complex absorption and reflection could be at play. Our preferred spectral model that fits the NuSTAR₂₀₁₅ + Swift${}_{2015}^{1}$ spectrum includes both a partial covering absorber and reflection, i.e., TBabs × (partcov × TBabs) × (reflect × mkcflow+ Gauss) in XSpec. We prefer this model over the others described below both on physical grounds and because of the better distribution of fit residuals revealed with the runs test statistic (or Wald–Wolfowitz test; a non-parametric test used to check the randomness of the residuals, where a run is defined as a succession of one or more identical values of the data that are followed and preceded by a different value). This model yields a shock maximum temperature of kT_max = 40 ± 3 keV and a reflection fraction of 1.1 ± 0.2 (model "Both" in Table 1).

Table 1.  NuSTAR₂₀₁₅ + Swift${}_{2015}^{1}$ and Suzaku₂₀₀₆ Spectral Fits

	NuSTAR + Swift			Suzaku
	Complex Absorption Only	Reflection Only	Both	Complex Absorption Only	Both
NH (ISM) (1022 cm−2)	0.05	0.05	0.05	0.05	0.05
NH (int.) (1022 cm−2)	32 ± 5	32 ± 1	⋯	27 ± 2	22 ± 3
NH (P.C.) (1022 cm−2)	162 ± 80	⋯	43 ± 3	49 ± 5	32 ± 3
Covering Fraction	0.4 ± 0.1	⋯	0.93 ± 0.06	0.71 ± 0.05	0.8 ± 0.1
kTmax (keV)	42 ± 5	39 ± 2	40 ± 3	49 ± 3	43 ± 3
${\dot{M}}_{\mathrm{thin}}$ (10−9 M⊙ yr−1)	0.9 ± 0.1	0.48 ± 0.02	0.55 ± 0.02	0.93 ± 0.07	0.77 ± 0.05
EW${}_{\mathrm{FeK}\alpha }$ (eV)		191 ± 2			138 ± 2
Z/Z⊙	1.00 ± 0.01	1.01 ± 0.01	1.04 ± 0.09	1.04 ± 0.08	1.03 ± 0.09
Reflection Fraction	⋯	1.8 ± 0.3	1.1 ± 0.2	⋯	1.1
${\chi }_{\nu }^{2}$	1.1	1.06	1.01	1.01	1.06
d.o.f	1090	1091	1090	1322	1321
Flux	0.95 ± 0.02	0.50 ± 0.01	0.58 ± 0.01	1.00 ± 0.05	1.11 ± 0.01
Luminosity	7.4	3.89	4.5	7.8	8.7
Runs Statistica	−2.50	−1.22	−1.16	⋯	⋯

Note. Unabsorbed X-ray flux and luminosity, in units of 10⁻¹⁰ erg s⁻¹ cm⁻²  and 10³³ erg s⁻¹, respectively, are calculated in the 0.3–50 keV energy band. Elemental abundances are quoted in units of abundances from (Wilms et al. 2000). Luminosity and $\dot{M}$ are determined assuming a distance of 800 pc. Statistical errors are calculated at the 90% confidence level. Systematic errors between Suzaku and NuSTAR are of about 5%–10% depending on the energy range (Madsen et al. 2017).

^aThe Runs test evaluates if the fit residuals are randomly distributed above and below zero. We tested if the hypothesis that the residuals are randomly distributed can be rejected at the 5% significance level, with higher values of abs(Runs) indicating the hypothesis can be rejected.

Download table as:  ASCIITypeset image

Using the more up-to-date absorption model TBabs (Wilms et al. 2000), we also fit the spectra using the model that Luna et al. (2008) fit to the Suzaku₂₀₀₆ spectrum, i.e., TBabs × (partcov × TBabs) × (mkcflow + Gauss). In this model, X-rays from the accretion disk boundary layer are attenuated by two high column density absorbers, one of which partially covers the source. This complex absorbing system was required to explain the lack of X-rays below 2 keV, and to reproduce the high energy curvature observed in the XIS and HXD Suzaku₂₀₀₆ spectra. We find a statistically acceptable fit to the NuSTAR₂₀₁₅ + Swift${}_{2015}^{1}$ data with the Suzaku₂₀₀₆ model; the best-fit parameters are shown in the first column ("Complex absorption only") of Table 1. This model yields a very high column density for the partial covering absorber of 1.62 × 10²⁴ cm⁻². The high column densities on the order of 10²⁴ cm⁻² and 10²³ cm⁻² of the partial and full covering absorbers, respectively, imply that, for the photons created behind the absorber, only those above 10 keV are being transmitted. This is the energy range precisely where reflected X-rays from the WD surface are expected to be observed. Furthermore, with a column density in excess of 10²⁴ cm⁻², the partial covering medium is itself optically thick enough to be a reflecting surface for X-rays. The results of X-ray spectral fitting using models without an explicit reflection component thus suggest that significant reflection was indeed present in the spectrum.

We therefore also tested simpler models in which the high energy X-ray spectral shape was due only to reflection of X-rays near the WD surface, with no complex absorption. We modeled this in Xspec as TBabs × (reflect × mkcflow+ Gauss). In this scenario, the cooling flow component is reflected by material near the WD surface, which modifies its shape above 10 keV. Finally, the Gaussian component models the 6.4 keV Fe Kα line, which is also due to reflection but is not included in the reflect model. The resulting fit is also statistically acceptable, with ${\chi }_{\nu }^{2}$/dof = 1.06/1091. The best-fit maximum temperature of the cooling flow was kT = 39 ± 2 keV, and the accretion rate through the optically thin X-ray plasma was lower than the value returned from the model "Complex Absorption Only," ${\dot{M}}_{\mathrm{thin}}=4.8\pm 0.2\times {10}^{-10}\ {M}_{\odot }$ yr⁻¹. The absorber had a column density of 3.2 × 10²³ cm⁻², and the reflection fraction was 1.8 ± 0.3. We disfavor this model based on the results from the runs test when compared with the model that includes both reflection and absorption (see Table 1).

The Suzaku₂₀₀₆ data alone cannot probe the presence of reflection due to the systematic uncertainties in the HXD background spectrum, which are energy-dependent. We can, however, apply the knowledge about reflection we have derived from the NuSTAR₂₀₁₅ data to the Suzaku₂₀₀₆ spectrum. Fixing the reflection fraction to 1.1 (as we found for the NuSTAR₂₀₁₅ data), we find a maximum temperature for the cooling flow in 2006 (kT = 43 ± 3 keV) that is very similar to that in 2015. The normalization for this model was slightly higher than for the NuSTAR₂₀₁₅ observation, suggesting a higher accretion rate through the optically thin portion of the boundary layer of ${\dot{M}}_{\mathrm{thin}}$ = 7.7 ± 0.5 × 10⁻¹⁰ M_⊙ yr⁻¹.

To summarize our spectral results: (i) the distribution of fit residuals, as measured by the runs test statistic, indicates that the model that includes reflection and complex absorption is preferable over other models; (ii) both the NuSTAR₂₀₁₅ and Suzaku₂₀₀₆ spectra can be described by highly absorbed cooling flow models. There is clear evidence of reflection in the NuSTAR₂₀₁₅ spectrum, and including this source of emission in the Suzaku₂₀₀₆ data resulted in acceptable model fits. The maximum temperature of the cooling flow was similar in 2006 and 2015, but the Suzaku₂₀₀₆ data implies a slightly higher accretion rate through the optically thin part of the boundary layer in 2006.

First, we used the readout streak method on the Swift UVOT data to estimate count rates through the U and UVW2 filters (corrected for coincidence losses) of 0.220 and 0.070 c s⁻¹, or Vega system magnitudes of 10.1 and 10.33, respectively. We used the UVW2 values to estimate the broadband UV flux of T CrB at the time of the NuSTAR₂₀₁₅ observation, which we take as a proxy of its disk luminosity.

Previously, Selvelli et al. (1992) analyzed the International Ultraviolet Explorer (IUE) observations of T CrB from 1979 January 5 through 1990 February 9. They estimated the 1250–3200 Å flux (hereafter F_UV) for each pair of SWP and LWP/LWR spectra, assuming a reddening of $E(B-V)$ = 0.15, and showed that T CrB was highly variable in the UV during the IUE era: F_UV ranged from 6.9 × 10⁻¹¹ erg cm⁻² s⁻¹ (on 1989 August 1) to 1.33 × 10⁻⁹ erg cm⁻² s⁻¹ (on 1983 May 1; see their Table 1). They measured an increase in F_UV from 1.57 × 10⁻¹⁰ to 3.02 × 10⁻⁹ erg cm⁻² s⁻¹ between 1989 March 2 and April 8, so T CrB is demonstrably capable of a factor of almost two change in approximately one month.

These IUE observations provided low-resolution, wide-band spectroscopy covering the 1250–3200 Å range, and the Swift UVOT data provide filter photometry in a small part of the IUE band. We downloaded archival IUE data from MAST,⁸ already reduced and fluxed, and combined short and long-wavelength spectra that were taken on the same day or on nearby dates. We folded them through the UVOT filter effective area curve to predict the magnitude that would have resulted, had Swift UVOT observations taken place during the IUE era. We also made our own estimates of F_UV during the 1980s, by dereddening the IUE spectra using $E(B-V)$ = 0.15 and integrating over the 1250–3200 Å range. Our values are similar to those found by Selvelli et al. (1992), although often ∼10% lower. This may be due to different reduction and/or calibration of the IUE data, or due to numerical differences in, e.g., the dereddening process. The only other caution regarding the results of Selvelli et al. (1992) is that they used an assumed distance to T CrB of 1300 pc, so their luminosity values need to be reduced by a factor of approximately 2.6 for a direct comparison with ours.

There are no major uncertainties in estimating equivalent Swift UVOT magnitudes from IUE spectra, apart from the above mentioned minor issues of exact reduction procedure and calibration. Going in the other direction, from Swift UVOT photometry to F_UV over the total IUE range, on the other hand, involves extrapolation of the source spectrum, which we do not directly observe, and hence is highly uncertain. It would be wrong to assume a typical accretion disk spectrum (${F}_{\lambda }\propto {\lambda }^{-2.33}$ power law) in the extrapolation; Selvelli et al. (1992) showed that T CrB had a highly variable UV spectrum that was usually much flatter than the theoretical accretion disk, except on one occasion when it had a power-law index of 2.2. However, we can convert Swift UVOT magnitudes to F_UV by assuming the observed IUE spectra, collectively, are representative of the UV spectra of T CrB at the respective UV flux levels. Indeed, the equivalent Swift UVOT (UVW2) magnitudes for IUE data show a tight correlation with F_UV (see Figure 3). During the Swift₂₀₁₅ observations, the UVW2 magnitudes were estimated to be 10.33–10.20, which corresponds to UV (1250–3200 Å) fluxes of 20.4–23.0 (in units of 10⁻¹⁰ erg cm⁻² s⁻¹) or luminosities of (1.5–1.7) × 10³⁵ erg s⁻¹. These values should be interpreted with caution because given that T CrB is so bright in the UV during the current high state, we are extrapolating the IUE-based UVW2 versus F_UV relationship. We use this estimate as our best (though imperfect) proxy for the luminosity of the Keplerian disk.

Zoom In Zoom Out Reset image size

A major question for accreting WDs, particularly those accreting at rates of around 10⁻⁹ M_⊙ yr⁻¹, is the accretion rates at which the boundary layer transitions (or begins to transition) from optically thin to thick (${\dot{M}}_{\mathrm{trans}}^{1}$) regimes and vice versa (${\dot{M}}_{\mathrm{trans}}^{2}$). These transition accretion rates are a function of M_WD and are closely tied to the assumed structure of the boundary layer. Observational constraints on them can thus provide insight into boundary layer structure. Given that the timescale for the transition in one direction has been observed to be different from the timescale of the transition in the opposite direction in dwarf novae, it seems reasonable to suspect that the accretion rates for these transitions could also be different. However, only a few theoretical predictions of these thresholds exist. For the purpose of comparing observations with theory, we return here to the simple picture in which all the empirical indicators of a boundary layer transition from optically thin to thick for a particular system appear at a single accretion rate ${\dot{M}}_{\mathrm{trans}}^{1}$, and all the empirical indicators of a boundary layer transition from optically thick to thin appear at a distinct, single accretion rate ${\dot{M}}_{\mathrm{trans}}^{2}$. Narayan & Popham (1993) concluded that when $\dot{M}$ = 3 × 10⁻¹⁰ M_⊙ yr⁻¹ and M_WD = 1 M_⊙, the boundary layer is generally expected to be optically thin, hot, and a source of hard X-rays, with a transition from this regime starting at ${\dot{M}}_{\mathrm{trans}}^{1}$ ∼ 10⁻⁹ M_⊙ yr⁻¹ (assuming a WD rotational speed of Ω = 0.5). At $\dot{M}\gtrsim {\dot{M}}_{\mathrm{trans}}^{1}$, Narayan & Popham (1993) found the boundary layer to be optically thick.

Later, Popham & Narayan (1995) studied the optically thick solutions for the boundary layers in CVs and derived values of ${\dot{M}}_{\mathrm{trans}}^{2}$ during the decay of an outburst for different values of M_WD. By assuming that the transition would occur at τ_* = 1 (which includes opacity from free–free absorption from a fully ionized gas and electron scattering), the authors found that in a nonrotating WD, the transition would occur at a rate of 7.5 × 10⁻⁷ M_⊙ yr⁻¹ (for a 1 M_⊙ WD). Taking a definition of ${\dot{M}}_{\mathrm{trans}}^{2}$ with a slightly smaller optical depth, ${\tau }_{* }$ = 0.8, yielded ${\dot{M}}_{\mathrm{trans}}^{2}({\tau }_{* }=0.8)$ = 4.6 × 10⁻⁸ M_⊙ yr⁻¹ for the same WD mass. Suleimanov et al. (2014) found that these thresholds might be overestimated due to an underestimation of the Rosseland opacity which might have previously been underestimated by two orders of magnitude.

Assuming that the WD mass can be estimated from the maximum cooling flow temperature in the X-ray spectral model that we refer to as "Both," we derive a mass for the WD in T CrB of at least 1.15 ± 0.03 M_⊙; in this case the theoretical models from Popham & Narayan (1995) predict an ${\dot{M}}_{\mathrm{trans}}^{2}$ of about 4.6 × 10⁻⁸ M_⊙ yr⁻¹ (taking τ = 0.8). The observations with XMM-Newton in 2017 (XMM-Newton₂₀₁₇) indicated an accretion rate of $\dot{M}\sim 6.6$ × 10⁻⁸ at that time (Luna et al. 2018a). Observations during the end of the current state will therefore be crucial to test the theoretical prediction for ${\dot{M}}_{\mathrm{trans}}^{2}$.

Observations of six dwarf novae compiled by Fertig et al. (2011) during quiescence and outburst showed that the observationally derived ${\dot{M}}_{\mathrm{trans}}^{1}$ are much lower than the theoretical values predicted by Popham & Narayan (1995). Moreover, the ${\dot{M}}_{\mathrm{trans}}^{1}$ for these six dwarf novae are different from each other, ranging from ∼1.6 × 10⁻¹⁰ M_⊙ yr⁻¹ in the case of SS Cyg to ∼1.2 × 10⁻¹¹ M_⊙ yr⁻¹ in the case of U Gem. Note, however, that Fertig et al. (2011) implicitly assumed a single ${\dot{M}}_{\mathrm{trans}}^{1}$, and took the highest recorded luminosity of the optically thin X-rays as corresponding to that ${\dot{M}}_{\mathrm{trans}}^{1}$. If the emergence of an optically thick boundary layer and the decrease in luminosity of the optically thin X-ray component happens at different transition $\dot{M}$, this procedure does not provide an unambiguous result. The behavior of U Gem is unusual for a dwarf nova, and is a case in point: it has simultaneously exhibited a detectable soft component from the optically thick boundary layer and an enhanced optically thin X-ray component during outburst (see, e.g., Mattei et al. 2000). During the 1993 December/1994 January outburst observed with the Extreme-Ultraviolet Explorer (EUVE) three times, the soft component persisted to an estimated accretion rate of 3 × 10⁻¹⁰ M_⊙ yr⁻¹ (scaled to a distance of 100 pc) so that could be taken as an indicator of ${\dot{M}}_{\mathrm{trans}}^{2}$ (Long et al. 1996). On the other hand, the hard X-ray component does not disappear in U Gem even when the accretion rate through the boundary layer exceed 5 × 10⁻⁹ M_⊙ yr⁻¹, if all outbursts of U Gem are similar to each other.

We observed an increase in accretion rate of several orders of magnitude between 2004 and 2018 (Luna et al. 2018a). Other evidence, primarily from IUE spectra obtained between 1979 and 1989, suggests that the accretion rate is generally variable by at least an order of magnitude (Selvelli et al. 1992). Observations obtained with Galaxy Evolution Explorer using the far ultraviolet (FUV) (λ1528 Å; Δλ442 Å) filter in 2006 July, show that the luminosity in the FUV filter was about L_FUV = 2.2 × 10³³ (d/800 pc)² (deriving the flux in the FUV band using gPhoton; Million et al. 2016). If this luminosity is roughly half of the available accretion luminosity, then the accretion rate at that time in 2006 was 2.2 × 10⁻¹⁰ M_⊙ yr⁻¹ (and with L_BL/L${}_{\mathrm{disk}}\approx 1$, the boundary layer was optically thin). Theoretical nova outburst models, such those presented by Yaron et al. (2005), predict that for a WD with a mass of about 1.0 M_⊙, accreting at ∼10⁻⁸ M_⊙ yr⁻¹, a nova outburst is produced every ≳2 × 10³ yr (see their Tables 2 and 3). The observed recurrence time of ∼80 yr in T CrB requires either a higher WD mass and average accretion rate, or a highly variable accretion rate. The episode of increased $\dot{M}$ reported in Luna et al. (2018a) might not have been unique, and similar events might have passed unnoticed.