A Long Look At MCG-5-23-16 With NuSTAR: I- Relativistic Reflection And Coronal Properties

TL;DR

This study analyzes NuSTAR observations of MCG-5-23-16, revealing stable relativistic reflection features, a correlation between cutoff energy and flux, and insights into the coronal plasma's composition and variability mechanisms.

Contribution

It provides the first detailed spectral analysis of MCG-5-23-16 over multiple years, highlighting stable reflection features and the relationship between coronal properties and X-ray emission.

Findings

Stable reflection spectrum over two years

Positive correlation between cutoff energy and flux

Coronal plasma contains mostly electrons, not pairs

Abstract

MCG-5-23-16 was targeted in early 2015 with a half mega-seconds observing campaign using NuSTAR. Here we present the spectral analysis of these datasets along with an earlier observation and study the relativistic reflection and the primary coronal source. The data show strong reflection features in the form of both narrow and broad iron lines plus a Compton reflection hump. A cutoff energy is significantly detected in all exposures. The shape of the reflection spectrum does not change in the two years spanned by the observations, suggesting a stable geometry. A strong positive correlation is found between the cutoff energy and both the hard X-ray flux and spectral index. The measurements imply that the coronal plasma is not at the runaway electron-positron pair limit, and instead contains mostly electrons. The observed variability in the coronal properties is driven by a variable optical depth. A constant heating to cooling ratio is measured implying that there is a feedback mechanism in which a significant fraction of the photons cooling the corona are due to reprocessed hard X-rays.

Figures

Tables

Table 1. Table 1 A summary of the observations used in this work.

Satellite	Obs. ID	UT Date	Exp. (ks)
NuSTAR	60001046002 (N2)	2013-06-03	160
	60001046004 (N4)	2015-02-15	210
	60001046006 (N6)	2015-02-21	98
	60001046008 (N8)	2015-03-13	220
Suzaku	700002010 (S0)	2005-12-07	191
	708021010 (S1)	2013-06-01	319
	708021020 (S2)	2013-06-05	221
Swift	0008042100[SW]	2015-01-28	$\sim 2$
	with SW = 2–11	to	each
		2015-03-13

Table 2. Table 2 Fit parameter for the relxill+xillver for the 1-10 keV fits for the Suzaku and NuSTAR datasets, and for 1-79 keV for the NuSTAR datasets. The 1 σ 𝜎 \sigma statistical uncertainty in the last significant figure is shown in brackets. N 𝑁 N is the normalization in units of log log {\rm log} (photons cm -2 s -1 ). Subscripts r 𝑟 r , p 𝑝 p and x 𝑥 x are for relxill , cutoffpl and xillver respectively. R 𝑅 R is the reflection fraction measured as the ratio of direct to reflected fluxes between 20–40 keV. R in subscript 𝑅 in R_{\rm in} is in units of gravitational radius r g subscript 𝑟 𝑔 r_{g} and the inclination θ 𝜃 \theta is in degrees. q 𝑞 q is the emissivity index, ξ 𝜉 \xi is the ionization parameter, Γ Γ \Gamma is the photon index and E c subscript 𝐸 𝑐 E_{c} is the cutoff energy. The observations in the first column are defined in Table 1 .

Obs	$N_H$	$\Gamma$	$R$	$\theta$	$R_{\mathrm{in}}$	$q$	$\log \left(\xi \right)$	$N_r$	$N_x$	$E_c$
Fits to individual spectra between 1–10 keV
S0	1.49(1)	1.84(3)	0.3(1)	57(7)		2(1)	2.6(1)	-3.25(1)	-3.38(6)	…
S1	1.41(1)	1.89(2)	0.2(1)	39(8)		3(1)	2.6(1)	-3.10(1)	-3.4(1)	…
S2	1.41(1)	1.90(1)	0.21(3)	47(7)	77(20)	$>6$	2.70(5)	-3.19(1)	-3.4(1)	…
N2	1.24(5)	1.91(1)	0.29(8)	34(7)	41(18)	5(2)	2.3(1)	-3.18(1)	-3.6(1)	…
N4	1.26(5)	1.83(2)	0.5(1)	64(6)			2.5(9)	-3.34(1)	-3.5(1)	…
N6	1.46(5)	1.88(2)	0.4(2)	52(13)		2(1)	2.3(8)	-3.29(1)	-3.6(5)	…
N8	1.34(5)	1.84(1)	0.29(7)	45(8)		3(2)	2.3(8)	-3.26(1)	-3.5(1)	…
Fits to individual spectra between 1–79 keV
N2	1.39(1)	1.87(1)	0.2(1)	51	$>50$	8(2)	2.5(5)	-3.23(1)	-3.45(1)	140(10)
N4	1.26(4)	1.77(1)	0.34(3)	…		1(1)	2.7(1)	-3.41(2)	-3.84(5)	120(8)
N6	1.41(5)	1.82(1)	0.25(7)	…		1(2)	2.5(1)	-3.33(1)	-3.7(1)	146(18)
N8	1.35(4)	1.77(1)	0.18(3)	…		0.2(7)	2.7(1)	-3.33(1)	-3.66(5)	124(7)
Joint fit to all the spectra between 1–79 keV
	$N_H$	$\Gamma$	$N_p$	$\theta$	$R_{\mathrm{in}}$	$q$	$\log \left(\xi \right)$	$N_r$	$N_x$	$E_c$
N2	1.40(1)	1.85(1)	-1.47(1)	51	47(23)	1.6(6)	2.66(3)	-3.80(1)	-3.60(2)	152(8)
N4	1.25(4)	1.78(1)	-1.68(1)	…	…	…	…	-4.21(6)	…	107(7)
N6	1.45(4)	1.83(1)	-1.59(1)	…	…	…	…	-4.12(7)	…	149(14)
N8	1.41(4)	1.79(1)	-1.60(1)	…	…	…	…	-4.11(1)	…	130(6)

Full text

1. Introduction

X-ray emission from AGN provides an excellent probe of the immediate vicinity of the central black hole. The spectra of Seyfert galaxies at hard energies ($>10$ keV) are characterized by a powerlaw that rolls over exponentially at energies around $\sim 100-200$ keV, along with a Compton reflection component. The powerlaw is the main X-ray source, produced by Comptonization of soft seed photons likely produced by the viscous dissipation of accretion energy (Haardt and Maraschi, 1991). Non-thermal Comptonization is not significant based on the detection of hard X-ray cutoffs and the non-detection of $\gamma$-rays (Zdziarski et al., 1995; Massaro et al., 2016).

The best spectra prior to NuSTAR were provided by CGRO-OSSE (Johnson et al., 1993), BeppoSAX (Boella et al., 1997) and INTEGRAL (Winkler et al., 2003). Modeling of the spectra obtained using these missions (Zdziarski et al., 1993; Petrucci et al., 2001; Perola et al., 2002; Molina et al., 2013; Malizia et al., 2014) provided some estimates of cutoff energies and reflection fractions, but offered only weak constraints on the physical models, mostly because the quality of the data was not high enough to break modeling degeneracies between the spectral index, the reflection strength and the cutoff energy. An additional difficulty was the requirement of simultaneous low energy ($<10$ keV) coverage that was not always available. Lubiński et al. (2016), for instance, compare results from different published studies with BeppoSAX, RXTE/HXTE, Swift-BAT and INTEGRAL, with lower ene,rgy coverage provided by several other instruments (mostly non-simultaneous), and found that estimates of the spectral index, the strength of reflection $R$ and the cutoff energy showed signi,ficant differences in values and correlations in the published work, even when analyzing the same datasets, mostly because of the low signal to noise ratio spectra above 30 keV.

The launch of NuSTAR (Harrison et al., 2013) offered a breakthrough for these studies. With its effective area and continuous energy coverage that includes the iron line at 6 keV and the Compton hump peaking at $\sim 30$ keV, NuSTAR provides a better handle on obtaining accurate estimates (or stringent limits) of the Comptonization parameters. Many estimates have been published so far for individual Seyfert and radio galaxies including IC 4329A (Brenneman et al., 2014), Ark 120 (Matt et al., 2014), MCG-6-30-15 (Marinucci et al., 2014b), SWIFT J2127.4+5654 (Marinucci et al., 2014a), NGC 4945 (Puccetti et al., 2014), MCG–5-23-16 (Baloković et al., 2015), NGC 2110 (Marinucci et al., 2015), NGC 5506 (Matt et al., 2015), NGC 4151 (Keck et al., 2015), NGC 7213 (Ursini et al., 2015b), 3C 273 (Madsen et al., 2015a), NGC 5548 (Ursini et al., 2015a), 3C 390.3 (Lohfink et al., 2015), 3C 382 (Ballantyne et al., 2014) and MRK 335 (Wilkins et al., 2015; Keek and Ballantyne, 2016). Cutoff energies in the range $100-500$ keV are measured, and although these values are outside the energy band of NuSTAR (3–79 keV), the fact that the spectra starts rolling over well below the cutoff energy, and the effect of the cutoff on the reflection spectrum that cannot be mimicked by other parameters, both allow accurate measurements of the cutoff energy (García et al., 2015).

There are many questions to address given these new observations. For instance, what is the geometry of the corona and what physical process controls the shape of the observed spectrum?. Population synthesis models for the X-ray background rule out significant emission above $\sim 300$ keV (Gilli et al., 2007; Ueda et al., 2014). Therefore, a cutoff should be observable in many sources. Additionally, depending on its compactness, the plasma cannot reach equilibrium at very high temperatures when photon-photon interactions become important. The production of electron-positron pairs acts as a thermostat, where increasing the power of the plasma produces more pairs rather than increases their temperature (Svensson, 1984; Zdziarski, 1985). Obtaining a measure of the electron temperature along with an estimate of the size of the corona is important in assessing to what extent these effects are important and help constrain the geometry of the corona.

Combining energy cutoff measurements from NuSTAR and size estimates from both spectral and timing information, Fabian et al. (2015) found that the implied electron temperatures are close to the boundary of the region in the compactness-temperature diagram which is forbidden due to runaway pair production. This suggests that pairs are an important ingredient in AGN coronal plasmas.

Theoretical models provide additional predictions that can be tested observationally. For example, in pair-free Compton-cooled coronae, an increase in cooling (keeping the power supplied to the electrons fixed) causes the temperature to drop, so $E_{c}$, the cutoff energy is anti-correlated with $\Gamma$, the photon index. On the other hand, in pair-dominated plasmas, Ghisellini and Haardt (1994) found that $E_{c}$ is positively correlated with the observed photon index $\Gamma$ for electron temperatures $T_{e}<m_{e}c^{2}$ (where $m_{e}$ is the electron mass and $c$ is the speed of light), while the reverse trend is predicted above it. Even below this limit, $E_{c}$ can remain constant for different values of $\Gamma$(e.g. Fig. 14b in Zdziarski et al., 2002). These arguments apply to a single source, and should be observable in a sample of sources, as has been explored with BeppoSAX and INTEGRAL (Lubiński et al., 2016), and is now being revisited with NuSTAR (Tortosa et al. in prep.).

In this work, we use a long-look observation of MCG–5-23-16 ($z=0.0085$; $M=10^{7.9}M_{\odot}$, Ponti et al. 2012) to attempt to address the question of cutoff variability and its relation to other parameters directly using a single bright object. We find that the cutoff energy is variable and shows interesting correlations with other parameters. The rest of this paper is organized as follows: Section 2 presents the data reduction and analysis, where we include data from NuSTAR, Swift (Gehrels et al., 2004) and Suzaku (Mitsuda et al., 2007) to obtain a complete spectral picture between 1-79 keV. Section 3 presents detailed spectral modeling, first for the iron line band and then for the whole observed band. The implications of the results are discussed in section 4. Results from the short time-scale variability will be published separately (Zoghbi et al. in prep.).

2. Observations & Data Reduction

2.1. Observations Log

Following the detection of both a cutoff energy (Baloković et al., 2015) and reverberation time delays (Zoghbi et al., 2014) in the first observation of 2013, MCG–5-23-16 was observed again with NuSTAR for a net exposure of $528$ ks in early 2015. The new observation was split into three exposures, with the first two separated by 2 days, and the third taken 22 days later. The observation ID’s, dates and exposures are shown in Table 1. Our analysis also includes data taken by Suzaku in 2013, simultaneously with the 2013 NuSTAR observation and also an earlier Suzaku observation taken in 2005 (see Table 1). Timing analysis of the 2013 Suzaku data was presented in Zoghbi et al. (2014), while the spectral analysis is presented here for the first time.

In order to obtain spectral coverage below 3 keV simultaneous with the new NuSTAR observations, we requested snapshots with Swift XRT. A total of nine exposures were taken while NuSTAR was observing the source, in addition to one that was simultaneous with the first NuSTAR observation. The IDs of these observations are also shown in Table 1.

2.2. Data Reduction

NuSTAR data were reduced and analyzed using the NuSTAR data analysis software, which is part of heasoft v6.19. The data were reduced by running the standard pipeline nupipeline. Source spectra were then extracted for modules A and B from circular regions 3 arcmin in radius centered on the source. Background spectra where extracted from source-free regions of the same size near the source. For the calibration files, we use caldb release v20160502. In the spectral analysis, spectra from modules A and B are fitted simultaneously, allowing for a multiplicative constant between them.

The XIS spectra from Suzaku were reduced also using the relevant software in heasoft v6.19. The initial reduction was done using aepipeline, using the caldb calibration release v20160607. Source spectra were extracted using xselect from circular regions 3 arcmin in radius centered on the source. Background spectra were extracted from a source-free region of the same size, away from the calibration source. The response files were generated using xisresp. Spectra from XIS0 and XIS3 were checked to be consistent and then combined to form the front-illuminated (FI) spectra. Comparing front- and back-illuminated spectra by fitting an absorbed power-law region between 3 and 5 keV gives an index that varies by $\sim 5\%$. Energies between 1.8 and 2.3 keV are ignored due to the calibration uncertainties associated with the CCD Si K edge. Data from other Suzaku instruments were not used because of low signal.

The Swift XRT observations were taken in the windowed timing mode except for SW2 and SW11 which were in imaging mode. The data were reduced using xrtpipeline. Source spectra were extracted from circular regions of 30 pixel radius (71 arcsec), and background spectra from similar regions away from the source. Observations SW2 and SW11 were taken in imaging mode and suffered some pileup. Therefore, the spectra were extracted from regions that excluded the central 3 pixels. We used the Swift   caldb release v20160121.

Spectral channels from NuSTAR, Suzaku XIS and Swift were grouped to have a minimum of one count per bin, and we use Poisson likelihood maximization in the modeling. Background spectra are handled using the W-statistics (Arnaud, 1996), where the background counts for each bin are considered fit parameters which can be solved for analytically as a function of other parameters. We use Poisson likelihood in order to be able to fit the swift low energy spectra and exploit the full resolution of NuSTAR at energies above 50 keV to constrain the high energy turnover. Using a Gaussian likelihood requires the channels to be grouped, which effectively removes some energy-dependent information at those energies.

The long term light curve from the Swift monitoring is shown in Figure 1. The flux changes by a factor of $\sim 3$ on a few days timescale, with the first 2013 (N2) observation having the highest flux. The NuSTAR observations sample the upper half of the flux variations observed with Swift. Given this variability and the simultaneity of observations, and for the purpose of constraining any variable column density at the host galaxy, we simultaneously fit observations that are taken within a day or so, or those having roughly the same 2–10 keV flux, modeling the NuSTAR and Suzaku data separately. Therefore, the following spectral groups are fitted with the same model, allowing for a multiplicative constant to account for cross calibration between instruments: [S0], [S1], [S2], [N2(A,B), SW8], [N4(A,B), SW3], [N6(A,B), SW6] and [N8(A,B), SW9]. We do not use SW2 as most of the counts are lost when correcting for pileup, removing its ability to constrain the column density below 3 keV. The cross calibration offset between the two NuSTAR modules is $<4\%$ is all cases, consistent with Madsen et al. (2015b).

3. Spectral Modeling

All the spectral modeling is performed in xspec v. 12.9.0n (Arnaud, 1996). The uncertainties quoted for the parameters are $1\sigma$ confidence, corresponding to $\Delta{\rm log}(Likelihood)=0.5$ (or $\Delta W=1$ where W is the W-statistic used in xspec), unless stated otherwise. The Galactic column in the direction of MCG–5-23-16 is $N_{\rm H}=9\times 10^{20}\,\rm{cm}^{-2}$ (Kalberla et al., 2005), and it is included in all subsequent fits using the model tbabs. We start this section by showing a model-independent representation of the data, then focus on the iron line first by using phenomenological models to track the long term variability, then use a physical model in section 3.1.2. We extend the analysis to the whole NuSTAR band in section 3.2.

Figure 2 shows the unfolded spectrum of MCG–5-23-16. Data from Suzaku observations S0 and S1 (see Table 1), along with the four NuSTAR observations, plotted after factoring out the effective area of the detectors. The spectra have been re-binned to a minimum signal to noise ratio of 6 for display clarity. The spectrum of MCG–5-23-16 is characterized by strong iron K emission and a broad excess at 30 keV, characteristic of a reflection spectrum. MCG–5-23-16 is seen through a Compton-thin absorber ($N_{\rm H}\sim 1.4\times 10^{22}\,\rm{cm}^{-2}$), and little emission from the nuclear regions escapes below 1 keV. The observed spectrum at these energies (below 1 keV), as revealed by XMM-Newton RGS spectra, is dominated by several emission lines superimposed on an unabsorbed scattered power-law continuum (Braito et al., 2007). They originate in a plasma in the Narrow Line Region. In the following Suzaku and Swift spectral modeling, we ignore spectral energies below 1 keV since they are not part of the nuclear emission.

3.1. The Fe Line Complex

The spectra of MCG–5-23-16 at the iron energies show both a narrow and a broad components (Reeves et al., 2007). The goal of this section is first to investigate whether the broad and narrow components are variable, and second to model the broad line with relativistic reflection models. All fits in this section are done in the 1–10 keV band.

3.1.1 Long Term Variability

We first attempt to check for column density variability using the Swift data alone. We fit an absorbed powerlaw to all the Swift spectra and track their changes. The column density showed some changes. However the strong degeneracy between column density and powerlaw index does not allow any firm conclusions to be drawn, and therefore analysis of all the spectral data is required.

To explore the changes in the iron line complex using a phenomenological approach, we show in Figure 3 the spectrum of the iron line from the highest and lowest flux NuSTAR observations (N2 and N4, respectively, which are separated by almost two years). The difference between the spectra is also shown. The spectrum in the iron line complex is clearly variable, and most of the variability is in the low energy wing of the line, not the core. The centroid energy for the N2 and N4 spectra are $6.30\pm 0.01$ and $6.36\pm 0.01$ keV, respectively, while the centroid from the difference spectrum is $6.17\pm 0.05$ keV. It appears therefore that the strongest changes in the iron line are in the broad component and not the narrow component.

To investigate this in a more systematic way, we fit the 1–10 keV band of the spectral groups using a powerlaw, a narrow and a broad Gaussian line. The width of the narrow line is fixed at the instrument resolution as neither NuSTAR nor Suzaku data are able to resolve it. Constraints to the column density for the NuSTAR data are provided by including Swift data, as discussed in section 2. The resulting parameter changes are shown in Figure 4.

NuSTAR observation N2 was simultaneous with Suzaku S2, but there are clear systematic differences between them, both in the absolute flux values and in the photon indices. These are due to absolute calibration uncertainties between the two detectors and the uncertain cross calibration between the front- and back-illuminated detectors in Suzaku. The photon index differences are the main cause of the large difference in $N_{\rm H}$ between S2 and N2 (which are degenerate in the fits). Nonetheless, it is clear that most of the variability is in the powerlaw continuum flux and its photon index. The energy of the narrow component is consistent with a constant. In fact, a model where the line energy is fixed to the average is statistically preferred over a model where it is free to change.

We assess the significance of the parameter changes using the sample-corrected Akaike Information Criterion (cAIC; Akaike, 1974; Burnham et al., 2011), where we compare models in which the parameter of interest is either fixed or allowed to change between observations. For the narrow line, we find that fixing the line flux between observation gives a lower cAIC than allowing it to change, with $\Delta\rm{cAIC}=-2$ from each observation. This implies that the model with a fixed line flux is preferred. For the line energy, $\Delta\rm{cAIC}<0$ for the Suzaku-only data and $\Delta\rm{cAIC}>10$ when the NuSTAR data are included. In other words, the line appears to change only between different instruments but is constant within the same instrument. This suggests that the line flux is physically constant and the observed changes are caused by instrument absolute calibration. For the column density, there appears to be significant changes between observations, even when using the same instrument ($\Delta\rm{cAIC}>15$ corresponding to a $>3.3\sigma$ significance). Changes in the parameters of the broad component are significant at the $3\sigma$ level when the narrow component is assumed constant. Degeneracies between the two reflection components are addressed in section 3.1.2 where we use full physical models.

3.1.2 Relativistic Model

The narrow Gaussian line is due to distant reflection from the Broad Line Region or the torus. To model it more physically, we replace the narrow Gaussian with the reflection model xillver (García et al., 2014). xillver includes self-consistently emission from Fe K$\beta$ and allows for the reflection to be ionized. Using the xillver+powerlaw model accounts for most of the residuals around 6.4 keV, including the Fe K$\beta$ line present in MCG–5-23-16, also seen in previous observations (Reeves et al., 2007). The xillver model provides a better fit compared to pexmon (Nandra et al., 2007), for example, with $\Delta W\sim 28$ per observation for the same number of degrees of freedom for Suzaku spectra where the K$\beta$ line is most clearly seen. This corresponds to $\Delta\rm{cAIC}=17$ and a significance of $>3.5\sigma$. The best fit ionization parameter of the reflector in this case is ${\rm log}(\xi)=0$, consistent with neutral reflection111$\xi$ in subsequent discussions is in units of erg s*-1* cm*-1*.

Continuing the analysis of the spectra in the 1–10 keV band, we next model the broad component of the iron line with a full relativistic reflection model. We use the relxill model (version v0.4a) (Dauser et al., 2010; García et al., 2014). The reflection spectrum is a result of a hard X-ray source illuminating a constant density disk. The observer sees emission from both the illuminating source and the reflector. The reflection spectrum is convolved with a relativistic kernel to model the strong gravity effects of the black hole. As we will show, the inner radius of the disk is $>10r_{g}$, and because the inner radius is degenerate with the black hole spin, we fix the latter at the maximum. The fact that the inner radius is $>10r_{g}$, the exact fixed spin value has little effect on the fits. We use a single powerlaw emissivity profile, and assume that the directly observed component has the same spectral index as that illuminating the disk. The high energy cutoff is fixed at 300 keV, and it is modeled when fitting the whole NuSTAR band in section 3.2. The abundance of the inner and outer reflectors are assumed the same. We fit all seven spectra independently, and the results are summarized in Table 2.

The reflection parameters are generally better constrained in the Suzaku spectra. This is because the resolution of NuSTAR in the iron K band does not allow the narrow and broad components of the line to be unambiguously separated. One effect of this is that the values of the iron abundance are not consistent between Suzaku and NuSTAR. The three Suzaku spectra suggested a value of about 1, while it is was around 0.5 (the model minimum) in the NuSTAR spectra. We therefore fixed the NuSTAR value at the Suzaku value. The uncertainties in Table 2 are calculated from Monte Carlo Markov Chains (MCMC) as the $1\sigma$ standard deviation of the chains. We found when exploring the likelihood space that it is multi-modal, with multiple parameter combinations having close likelihood values near the maximum. MCMC is therefore well suited to explore this multi-modality. All the chains reported here were generated using the affine invariant ensemble sampler (Goodman and Weare, 2010). All the chains were run several times and long enough to ensure convergence. The convergence is assessed with both the autocorrelation of the chains and the stability of the chain variances.

Table 2 quotes only the average values. There are however two general solutions in the modeling with different values for the ionization parameter $\xi$ of the relativistic reflection component (${\rm log}(\xi)\sim 0$ and ${\rm log}(\xi)\sim 2.7$). This is illustrated in the top panel of Figure 5, where we plot the confidence contours for the best fit parameters for (the log of) the ionization parameter versus the photon index for all the observations. Although the best fit value for the ionization parameter is ${\rm log}(\xi)\sim 2.3-2.7$, lower values ($\leq 1$) are also possible with a slightly higher value to $\Gamma$. Two parameters which are of interest are $R_{\rm in}$ and the emissivity index, which locate the emitting region relative to the central object. These two parameters are highly correlated in the modeling as shown by the middle panel of Figure 5. Although the best fit values are at $R_{\rm in}\sim 20$ gravitational radii ($r_{g}=GM/c^{2}$) and $q\sim 3$, lower and higher values for $R_{\rm in}$ are also supported by the data.

The confidence contours of the inclination $\theta$ and the normalization of the distant reflector $N_{\rm x}$ are shown in the bottom panel of Figure 5. The best fit inclination value in this case is $\theta\sim 50^{\circ}$, but higher inclinations with stronger distant reflection are also supported by the data. The data in this case supports either a low inclination disk ($\theta<60^{\circ}$) with a relatively weak distant reflector or a highly inclined disk ($\theta\sim 88^{\circ}$) and stronger reflection. Given that this object is seen through a Compton-thin rather than thick absorber, suggests an intermediate inclination, making the first solution more physically plausible. We note here that parameters reported in the analysis of the first NuSTAR observation (Baloković et al., 2015) are consistent with one of the local minima in the fit, with ${\rm log}(\xi)\sim 0$. Our best fit has a value of ${\rm log}(\xi)\sim 2.7$. Other parameters change accordingly, with reflection parameters not differing significantly (apart from the normalization).

3.2. Full Band Relativistic Model

Here, we extend the analysis to higher energies (up to 79 keV) and focus on the full NuSTAR data. Column density constraints are provided by Swift XRT for observations N4, N6 and N8, while for N2 we use the Suzaku observation S2 as it is of higher quality and it is simultaneous with N2. We model the spectrum with a model similar to that in section 3.1.2, except that we fit for a powerlaw explicitly so we can track flux variations of individual components and allow for the cutoff energy to be a free parameter. The model has an xspec form: tbabsztbabs (relxill+cutoffpl+xillver). Allowing all parameters to be free showed the same strong correlations between the inclination and distant reflection flux discussed in section 3.1.2, with the highly inclined disk being the best fit. Because this is physically unlikely, we fix the inclination at $\theta=51^{\circ}$, the weighted average from modeling the Suzaku spectra in section 3.1.2. We found that the exact value for the inclination does not affect the following results significantly. The main effect is that a lower (higher) value causes the distant reflector to be weaker (stronger), as already suggested by the bottom panel of Figure 5.

The best fit parameters from modeling the whole NuSTAR band are shown in Table 2. The best fit model and residuals are shown in Figure 6. The residuals in the middle set of panels in Figure 6 have been binned so that the signal to noise ratio of every spectral bin is at least 6. This is done so the residual plot is meaningful. The residuals at the high energy part of the spectra are shown in the bottom set of panels, plotted as the residuals of integrated (unbinned) data to the integrated model. This is a convenient way of plotting to account for the fact that bins at these energies have small number of counts. The bottom set of panels show that the deviations between the model and the data are comparable to the deviations between the two NuSTAR modules due to counting noise or cross calibration uncertainties.

The goodness-of-fit statistic is estimated using Monte Carlo simulations. We start from the best fit parameters and generate a large number of parameters drawn from the MCMC chains. For each parameter set, a spectrum is faked using fakeit in xspec, taking into account counting noise. The faked spectra are then refitted with the model and a distribution of fit statistics is produced from the resulting fits. The fraction of simulated data that have a fit statistic that is at least as good as the observed value are 0.96, 0.59, 0.21 and 0.4 for observations N2, N4, N6 and N8, respectively. These goodness parameters correspond to the probability of rejecting the null hypothesis corresponding to the best fit model. These fits are very good given the high quality data. The high value for the N2-S2 combination is due to cross calibration uncertainties between NuSTAR and Suzaku and between the front- and back-illuminated detectors in Suzaku.

Most of the parameters remained similar to those in section 3.1 and Figure 5. $R_{\rm in}$, $q$ and the parameters of xillver were all consistent with being constant between observations. Therefore, and in order to obtain further constraints on the variable parameters, we fit all four NuSTAR observations (and the matching Suzaku and Swift datasets) together, allowing only parameters that showed variability in the individual fits to vary. These parameters are: $N_{\rm H}$, $\Gamma$, $E_{c}$ and the normalizations of relxill ($N_{r}$) and cutoffpl ($N_{p}$). The best fit parameters in this case are presented in Table 2.

The column density appears to change between observations in a way that is not directly related to the observed flux. The variability is however only marginally significant. The 99.5% confidence limits on $N_{\rm H}$ are consistent with a constant column. The remaining parameters change significantly between observations. The results of their variability is summarized in Figure 7, which shows the best fit parameter confidence contours for $\Gamma$, reflection fraction222What we refer to as reflection fraction here is reflection strength in the nomenclature of Dauser et al. (2016). $R$ and the 2–10 keV flux of the powerlaw component ($F_{p}$), plotted against the high energy cutoff. Although flux from both the powerlaw and reflection components change, their ratio remains relatively constant. The difference between high- and low-flux model spectra in this case resembles the shape of the relativistic reflection (plus a powerlaw), and it directly explains the observed difference spectrum between N2 and N4 shown in Figure 3.

The photon index is correlated with the high energy cutoff $E_{c}$ within individual spectra. This is a consequence of the model parameterization when the data above $\sim 50$ keV have a relatively low signal to noise ratio. Additionally, both $\Gamma$ and the continuum flux appear to be correlated with the cutoff energy when all four observations are considered. The reflection fraction is independent of cutoff energy. To quantify these correlations, we use the MCMC chains already calculated to calculate the Pearson correlation coefficient $r$. We find correlation coefficients of 0.90(6), 0.0(3) and 0.7(1) for the relations between $E_{c}$ and $\Gamma$, $R$ and $F_{p}$, respectively. The number in bracket is the uncertainty in the last significant digit taken as the $1\sigma$ sample standard deviation.

We emphasize that $R$ in Figure 7 is the reflection fraction from the relativistic reflection. The reflection fraction from the distant reflector (not plotted), increases with $E_{c}$, driven by the $F_{p}-E_{c}$ correlation and the fact that the flux from the distant reflector is constant. Figure 7 also shows one measurement from a BeppoSAX taken in 1998 (Perola et al., 2002). That measurement appears to follow the same trends we observe, albeit with larger uncertainties.

We note here that the two, low and high $\xi$, solutions found when fitting data below 10 keV (section 3.1) no longer give a comparable goodness of fit. The higher $\xi$ solution gives a better fit by $\Delta W>30$, corresponding to a significance of $>3.7\sigma$ per observation. This solution is therefore preferred over that reported in Baloković et al. (2015) for observation N2. Forcing the $\xi\sim 0$ solution gives, in addition to a worse fit, lower cutoff energies but the correlations in Figure 7 hold.

We also note the low $N_{\rm H}$ value for N4. It is unlikely that the column density changes significantly in days time-scale. We therefore tested tying the $N_{\rm H}$ values between observations. We found a slightly worse fit ($\Delta W\sim 14$ or a significance of $\sim 2\sigma$), with no large changes in $E_{c}$. The reason is that NuSTAR data quality is much better than Swift, so forcing a new $N_{\rm H}$ does not affect the NuSTAR data significantly but makes the Swift fit slightly worse.

4. Discussion

4.1. X-ray reflection and the inner disk

The most detailed analysis of the broad component of the iron line in MCG–5-23-16 prior to this work was presented in Reeves et al. (2007). The results presented here regarding the reflection spectrum are consistent with that analysis. Taking the best fit parameters for the relativistic reflection suggest a truncated disk at $\sim 40r_{g}$ with a emissivity index close to a standard non-relativistic value of 3. There is, however, a strong degeneracy between the inner radius and emissivity index parameters, such that the data also supports (at the 99% confidence level) a disk that extends down to $R_{\rm in}<6r_{g}$, with a flatter emissivity index ($q<2$). X-ray reverberation lags detected in this source, taken at face value suggest that the latter solution is preferred (Zoghbi et al., 2014).

The reflection component, producing the narrow core of the iron line and a considerable fraction of the reflection hump at 30 keV, appears to be constant across the two year period spanned by the data, unless the inclination of the disk changes, which is unlikely. This is not surprising if emission comes from material far from the central source, thereby smoothing out any variability. A consequence of this observation is that the reflection fraction from this component varies, even when the flux and index of the illuminating source are the only parameters that change, not the geometry of the system. A correlation of the distant reflection strength $R$ with the flux in a single object is therefore naturally explained (e.g. Malzac and Petrucci, 2002).

The flux from the inner reflector on the other hand, tracks closely variations of the illuminating source. The strength of the reflection (measured as the ratio of fluxes between 20 and 40 keV) remains constant, a result also seen in other objects (e.g. Lohfink et al., 2016). Combining this with the fact that the column density changes very little, and also the lack of significant changes in the relativistic reflection parameter, gives a picture in which the flux changes seen in this source (e.g. Figure 1) are driven by intrinsic flux fluctuations in the primary source, which are closely matched, on days to months timescales, by variations in the flux of the relativistic reflection component. This seems to be the long term extension of the relativistic reverberation signatures seen on short time-scales in this object (Zoghbi et al., 2014).

4.2. Plasma properties

Modeling the reflection spectrum properly allows us to extract information about the Comptonization process. We find that most of the variability between observations is due to changes in the flux of the primary powerlaw component. Its photon index and cutoff energy are also found to be significantly variable, while the relativistic reflection component remains constant in shape and constant in flux relative to the primary component. The primary continuum changes in a structured manner, as indicated by the high correlation coefficients between the photon index and flux with the cutoff energy.

It is known that in the simple cutoff powerlaw model, the flux, index and cutoff energy can be correlated by construction. This is apparent in the elongated contours shown in Figure 7 from individual fits. It should be noted however, that the correlations between observations are robust, in a sense that if the parameters from one observation are fixed at the best fit parameters of another observation, the fit significantly worsens. Another way to see this is to observe that in the $\Gamma-E_{c}$ and $F_{p}-E_{c}$ plots in Figure 7, the elongated contours are not parallel to the observed correlations, indicating that although the parameters might be correlated within a single spectrum, their correlations between observations are robust.

4.2.1 $\Gamma-E_{c}$ Correlation

Previous results on possible $\Gamma-E_{c}$ correlations in samples of objects were not conclusive. Piro et al. (1999) first noted a possible $\Gamma-E_{c}$ relation using two BeppoSAX observations of NGC 4151. Petrucci et al. (2001) found a weak relation with six Seyfert Galaxies, and a relatively stronger relation was found by Perola et al. (2002) using a slightly larger sample. Using INTEGRAL data, Molina et al. (2013) reported a weak relation while Malizia et al. (2014) reported no relation. It appears therefore that as far as a sample of AGN is concerned, there is at most a weak relation between $\Gamma$ and $E_{c}$.

The data for individual objects is less clear. Mostly because of the difficulty of obtaining high quality data in single epoch observations, with low energy coverage. We note that from the INTEGRAL study of Molina et al. (2013), who analyzed separate observations of individual objects, in almost all cases of objects with multiple observations, flatter spectra are accompanied by small cutoff energies. The uncertainties in the parameters are, however, large. Using NuSTAR, Ballantyne et al. (2014) analyzed two observations of the radio galaxy 3C 382 in two flux intervals. The low flux observation had a flatter spectrum and a higher cutoff energy compared to the higher flux observation (i.e. opposite the trend seen in INTEGRAL data and seen here in MCG–5-23-16). Also using NuSTAR data, Keek and Ballantyne (2016) found a positive correlation between $\Gamma$ and $E_{c}$ in Mrk 335, although we note that both the photon indices and cutoff energies ($<50$ keV) found there are small, differing substantially from other studies using the same datasets (Parker et al., 2014; Wilkins et al., 2015). The result we found here suggests a strong positive correlation between the photon index $\Gamma$ similar to Mrk 335.

4.2.2 $Luminosity-E_{c}$ Correlation

Variability of the cutoff energy (or electron temperature) with flux or luminosity has not been explored in detail extensively in AGN, unlike black hole binaries. As we pointed out in section 4.2.1, results from AGN are not conclusive yet, where both a correlation and an anti-correlation of $E_{c}$ and flux have been reported for 3C 382 and Mrk 335 respectively (Ballantyne et al., 2014; Keek and Ballantyne, 2016).

For X-ray binaries, Cyg X-1 showed an increase in cutoff energy when the luminosity in the hard X-rays drops during the hard state (Gierlinski et al., 1997; Ibragimov et al., 2005). A strong anti-correlation of the cutoff energy and luminosity was also observed in GX 339-4 when the luminosity was above $\sim 10\%$ the Eddington luminosity, and it remained constant below that (Miyakawa et al., 2008). A similar result was found by Motta et al. (2009), who additionally observed the reverse trend when the source softened before transiting to the soft state. Similar behavior is seen in other objects, including V404 Cyg in its recent outburst as observed with Fermi (Jenke et al., 2016), and possibly also Cyg X-1 in the soft state from NuSTAR observations (Walton et al., 2016). We note that when $E_{c}$ is correlated with luminosity, so is $\Gamma$, and when the spectra soften during the hard intermediate state, both $\Gamma$ and $E_{c}$ reverse their dependence on luminosity. The result we find for MCG–5-23-16 appears to match the behavior of black hole binaries not in the hard state where there $E_{c}$ is anti-correlated with $L$, but during the intermediate state. A correlation between $\Gamma$ and the flux is well established in bright AGN and Galactic black holes (e.g. Sobolewska and Papadakis, 2009; Yang et al., 2015).

4.2.3 Physical Interpretation

Before discussing the details of the plasma physics, it is worth mentioning the possibility that the changes in $E_{c}$ may not be due to changes in the intrinsic election temperature, but rather due to changes in the gravitational redshift of a constant spectrum (e.g. Niedźwiecki et al., 2016). Such changes in the gravitational redshift of the emitted photons (due to changes in the size of the corona for instance) could artificially introduce variations in $E_{c}$ without the plasma properties changing. This, however, also produces changes in the reflection fraction due to the focusing of light rays into and away from the disk. The fact that the observed reflection fraction is constant suggests that the geometry does not change significantly, strengthening the interpretation in which the $E_{c}$ variations are intrinsic to the plasma. Also, the inner radius we measure is relatively large, so GR effects are present but not extreme, and therefore the discussion of relativistic modeling in Niedźwiecki et al. 2016 are not applicable in this work.

In the simplest considerations of a pure thermal plasma (Sunyaev and Titarchuk, 1980) (no electron-positron pairs), an increase in the soft flux impinging on the corona leads to softer Comptonization spectra and lower temperatures as the electrons are cooled efficiently (e.g. Zdziarski et al., 2002). This picture cannot be applied directly here, first because pairs are not included and their effect could be important (Fabian et al., 2015), and second because the correlation we measure is in the hard flux (emitted by the corona) rather than the soft flux impinging on it. Therefore, we would like first to assess the importance of pair production given the measurements we have.

We start by comparing the observations to predictions of pair-dominated plasmas. The temperature of a plasma cannot be arbitrarily high for a given size. The key parameter that is often used is the compactness $l=4\pi(m_{p}/m_{e})(r_{g}/r)(L/L_{\rm edd})$, which measures the luminosity to source size ratio. As the compactness increases, photon-photon interactions become important, and any extra heating goes into producing electron-positron pairs rather than heating the plasma, causing the temperature to saturate (Guilbert et al., 1983; Zdziarski, 1985).

Following Zdziarski (1985) and Ghisellini and Haardt (1994), we calculate the spectral index $\alpha$ ($\alpha=\Gamma-1$) and cutoff energy predicted from models in thermal and pair equilibria, for different values of the compactness $l_{h}$ and $l_{h}/l_{s}$, where $l_{h}$ is the compactness of the Comptonization plasma (i.e. the power heating the corona), and $l_{s}$ is the compactness of the soft source producing the seed photons. We use the model eqpair (Coppi, 1999) to generate spectra for a grid of parameters for $l_{h}$ (between 10 and $5\times 10^{4}$) and $l_{h}/l_{s}$ (between 1 and 100). We assume the soft source is a blackbody with a temperature of 10 eV and the plasma is spherical and contains no background electron plasma ($\tau_{p}=0$). The generated spectra are then fitted with a cutoff powerlaw model to simulate the fitting procedure and obtain the energy spectral index $\alpha$ and $E_{c}$. The results are shown in the left panel of Figure 8. Contours of $\alpha$ and $E_{c}$ are shown, and the $1\sigma$ measured values of $\alpha$ and $E_{c}$ are shown with the green boxes. This plot shows that, first, the inferred compactness ratio $l_{h}/l_{s}\sim 5-6$, which is not atypical of AGN. Secondly and most importantly, it shows that, if the plasma is dominated by pairs, then the inferred $l_{h}$ is large (for reference, a source radiating at the Eddington limit that is 1 gravitational radius in size has $l\sim 10^{5}$). Therefore, based just on the the measured $\alpha$ and $E_{c}$, if the plasma is dominated by pairs, the source has to be very compact suggesting a small size and/or high luminosity.

Further information is provided by including the source size measurement we have from the reflection spectrum and the observed luminosity. The results in this case are shown in the right panel of Figure 8, showing the temperature-$l$ relations at the pair limit from the modeling of Stern et al. (1995) for three geometries of the corona. For a given compactness and a geometry, a source cannot have a temperature above the lines shown. Below the lines, the effect of pairs decreases and the plasma contains only electrons. As we have discussed in section 3, the reflection spectrum constrains the inner radius of the disk to be $R_{\rm in}=47\pm 23\,r_{g}$. If we assume the corona is of the same size (otherwise, relativistic features in the reflection will be washed out), we obtain the green circles shown in Figure 8. We used the cutoff powerlaw flux in the range 0.1–200 keV to measure the luminosity assuming a standard cosmology ($\Omega_{m}=0.3,\Lambda=0.7$) with $H_{0}=68\,km\,s^{-1}\,Mpc^{-1}$, and a black hole mass of $M=10^{7.9}\,M_{\odot}$ (Ponti et al., 2012)333We follow a similar procedure to Fabian et al. (2015).. The red arrows show the upper limits on the observed $l_{h}$ obtained by setting $R_{\rm in}=1.23r_{g}$, the innermost stable circular orbit for a maximally spinning black hole. The electron temperature is estimated as $kT_{e}\sim E_{c}/2$ (Petrucci et al., 2001).

We can see that all the measurements fall below the pair limit lines for the three geometries, indicating that the plasma in MCG–5-23-16 is not dominated by pairs and consists mostly of electrons. This is a robust statement given the small uncertainties in the cutoff measurements. Using the nthcomp (Zdziarski et al., 1996) model (as implemented in the reflection model relxillCp which calculates the reflection spectrum when illuminated by nthcomp; Garcia et al. in prep.) instead of the cutoffpl shifts the electron temperatures up only by $\sim 10\%$, and our conclusion about the plasma content is not altered. The points can of course be shifted to the right (i.e increasing $l_{h}$) if the black hole mass is erroneous. We find that in order for the plasma to be in pair balance for the slab geometry, the black hole mass needs to be smaller by $\sim$ two orders of magnitude (i.e. $M\sim 10^{6}\,M_{\odot}$). We note that the black hole mass in Ponti et al. (2012) is uncertain. A mass estimate using the fundamental plane of black hole activity (Gültekin et al., 2009) using the radio fluxes from (Mundell et al., 2009) and the X-ray fluxes from our observations suggest a lower mass of $10^{7}\,M_{\odot}$. Other estimates suggest similar smaller values (e.g. Zhou and Zhao, 2010), but not low enough to alter our conclusions about the plasma content. Note also that the correlation between $E_{c}$ and $l_{h}$ is weaker than the $E_{c}-F_{p}$ correlation in Fig. 7. That is because we use the wider energy range to calculated the flux and also because of the uncertainties in the radius measurements.

4.2.4 Cutoff Variability

The additional information provided by the four measurements of the plasma properties provide further constraints. The cutoff energy in pair-dominated plasmas scales inversely with compactness (or with luminosity when the source size is constant, as is the case here where we measure a constant reflection fraction and reflection parameters): $E_{c}\propto l_{h}^{-1}$ (Svensson, 1994). This comes from the fact that increasing $l_{h}$ produces more pairs, and for balance to hold, the same energy is now distributed to more particles so the energy per particle drops. This trend is not what we observe in MCG–5-23-16. Instead we find that the cutoff energy is higher for higher fluxes, providing further evidence that the plasma is not dominated by pairs. In electron plasmas on the other hand, the electron temperature is not expected to depend on $l_{h}$ for a given ratio $l_{h}/l_{s}$ and optical depth $\tau$. This comes from the fact that although electrons gain energy when $l_{h}$ increases, the constant $l_{h}/l_{s}$ means cooling also increases so to keep the temperature constant. The fact that the cutoff energy (or electron temperature) varies with on $l_{h}$ observationally means that either $l_{h}/l_{s}$ or $\tau$ is variable. The first is ruled out by virtue of Fig. 8-left. We conclude, therefore, that the optical depth $\tau$ varies in a way to produce the $E_{c}-$flux relation in Figure 7.

We now turn to the $E_{c}-\Gamma$ correlation. In pair-dominated plasma, and for temperatures below the rest energy of the electron (as observed here), $E_{c}$ is expected to be either positively correlated or independent of $\Gamma$, depending on the temperature and the compactness $l_{h}$ (Ghisellini and Haardt, 1994; Zdziarski et al., 2002). In pair-free plasma on the other hand, an $E_{c}-\Gamma$ anti-correlation is expected for a fixed $l_{h}$. However, as we have already established observationally, $l_{h}$ is not constant, which implied a variable optical depth $\tau$. To investigate the variable $\tau$ possibility further in pair-free plasma, we employ a similar method to that used to produce Figure 8 and discussed in 4.2.3, but now we use the model compps (Poutanen and Svensson, 1996), and use a grid of $T_{e}$ and $\tau$. The results are shown in Figure 9 for a slab geometry (geometry parameter in compps set to 1). We find that the optical depth $\tau$ varies significantly with $T_{e}$, being lower for higher temperature. The exact values of $\tau$ depend on the assumed geometry (e.g. $\tau$ varies between 1.3–1.9 if we assumed a spherical geometry instead of the slab geometry shown in Figure 9).

The conclusion here is that, in order to explain the $E_{c}-\Gamma$ and $E_{c}-F_{p}$ correlations with, the plasma need to be pair-free and its optical depth need to change as shown in Figure 9. This could possibly be accompanied by geometry changes, but it has to be small enough for the inner radius and $R$ measurements to remain constant within the observational uncertainty.

One additional observation that can be noted from Fig. 8-left, is that the ratio of heating to cooling in the corona changes very little between observations. The observed changing coronal flux therefore suggests that the photon flux cooling the corona changes too, and in the same direction. This could be achieved if the UV photons of the disk vary with the X-rays. The UV flux from the source measured with the Swift UVOT camera however, shows little variations compared to the X-rays (the fractional RMS variation in X-rays is $24\pm 2\%$ while in the UV, it is $10\pm 4\%$ in the W2 filter and no more than $5\%$ in other bands red-ward of $2000{\rm\AA}$). The soft flux reaching the corona could in principle change if only the disk temperature changes, so the flux in the UVOT filters is not affected. This however would suggest that the inner radius changes, which not observed. The constant heating over cooling we find implies that there is a feedback between the hot corona and the disk photons cooling it (e.g. Haardt and Maraschi, 1991), suggesting that a significant part of the photons cooling the corona are due to reprocessing in the disk.

Comparing our results with Fabian et al. (2015) indicates that AGN coronae are not always pair-dominated, or that some sources are in that regime and others are not. Mrk 335 appears to show a behavior similar to that reported here for MCG–5-23-16 (Keek and Ballantyne, 2016), so it too, is unlikely to be pair-dominated. The presence or absence of reflection close to the black hole is unlikely to be be the reason of the difference, nor is the Eddington ratio (Mrk 335 is accreting close to the Eddington limit while MCG–5-23-16 accretes at the few percent level). Studies of other sources with NuSTAR in the near future will help address the issue.

5. Conclusion

We use data from the longest NuSTAR observing campaign of a Seyfert galaxy to study the properties of the plasma responsible for the hard X-ray emission. The sensitivity of NuSTAR allows us to constrain the plasma properties and probe its variability. Our main results are as follows:

• •

The inner radius of the disk and its emissivity remain constant between observing epochs, suggesting a constant geometry. Most of the spectral variability is due to changes in the flux and spectral index of the primary X-ray source. Flux from the relativistic reflection follows the flux from the direct component.

• •

The measured cutoff energies (and inferred electron temperatures) are not high enough for the plasma to be dominated by electron-positron pairs, unless the black hole mass is two orders of magnitude lower. This means that the plasma contains mostly electrons.

• •

We find that the cutoff energy is strongly correlated with both the source flux and the spectral index. The former correlation is another indication that the plasma is not dominated by pairs. The two correlations are driven by changes in the optical depth of the plasma.

• •

A constant heating to cooling ratio is inferred for the plasma. This, along with the constant UV flux observed, suggest a feedback mechanism in which most of the photons cooling the hot corona are due to reprocessing in a cold disk.

Acknowledgment

We thank the referee A. Zdziarski for the useful comments and suggestions that helped with the interpretation of the data. This work has been partly supported by NASA grant NNX14AF89G. This work made use of data from the NuSTAR mission, a project led by the California Institute of Technology, managed by the Jet Propulsion Laboratory, and funded by the National Aeronautics and Space Administration. We thank the NuSTAR Operations, Software and Calibration teams for support with the execution and analysis of these observations. This research has made use of the NuSTAR Data Analysis Software (NuSTARDAS) jointly developed by the ASI Science Data Center (ASDC, Italy) and the California Institute of Technology (USA).