NuSTAR AND CHANDRA INSIGHT INTO THE NATURE OF THE 3–40 keV NUCLEAR EMISSION IN NGC 253

We began by extracting the NuSTAR spectra from the nuclear region of NGC 253 for each of the three observations. In order to encompass a large fraction of the source counts in the nuclear region, while minimizing contamination from nearby unrelated sources, we chose to extract the spectra using circular regions with radii of 20 arcsec (≈400 pc), somewhat larger than the NuSTAR half width at half maximum (i.e., ≈9 arcsec). This choice of aperture encompasses 32% of the NuSTAR PSF (half-power radius 29 arcsec). Background spectra were extracted from larger apertures (17–37 arcmin² depending on the observational epoch), which were defined by eye as source-free regions. We note that NGC 253 X-1, the bright off-nuclear source ≈30 arcsec to the south of the nucleus, is a ULX with L_2–10 keV ≈ (2–3) × 10³⁹ erg s⁻¹ (based on analysis below) and provides some non-negligible contribution to the nuclear region spectra (see Figure 2(a)). To estimate these contributions, we extracted spectra of this source using 20 arcsec radii circular apertures. The spectra of NGC 253 X-1 were fit using xspec v.12.7.1 (Arnaud 1996) and a broken power-law model bknpower with varying normalization. We found a best-fit low-energy slope of Γ₁ = 2.1 ± 0.1, break energy E_break = 6.8 ± 0.2 keV, and high-energy slope Γ₂ = 4.1 ± 0.1; no absorption was required in any of these fits. These values are similar to those of ULXs studied in detail by Gladstone et al. (2009), which have L_2–10 keV = (1–14) × 10³⁹ erg s⁻¹, and are fit well by broken power-law models with Γ₁ = 1–3, E_break = 3–8  keV, and Γ₂ = 2–7. Using a PSF model, we estimate that the contributions of NGC 253 X-1 to the nuclear region spectra are small (≲10% for all three epochs). From our best-fit spectral models for NGC 253 X-1, we constructed fake background spectra corresponding to its contributions to the nuclear region spectra. These were used when performing spectral fits to the nuclear region itself (see below). To test the robustness of this procedure, we measured the spectral contributions from the wings of the nuclear region PSFs to see if they influenced the extracted spectra of NGC 253 X-1. Once contributions from the nuclear region PSFs were subtracted, we found only very small differences (at the ≲7% level) in spectral parameters and normalizations derived for NGC 253 X-1. Such differences translate into ≲1% differences in the derived spectral properties of the nuclear region itself, well below the precision of our measurements. We therefore conclude that our procedure for estimating the contributions of NGC 253 X-1 to the nuclear region is robust.

In Figure 3(a), we present the NuSTAR spectra of the central nuclear region for the three epochs. We do not include the Chandra spectra of the nuclear point sources here, since a direct matching between the multiple components that lie within the NuSTAR beam is complex and unnecessary for understanding the hard X-ray properties of the nuclear region that we are interested in here. We know that the NuSTAR emission from the central nuclear region is a combination of hot plasma, as well as non-variable and variable power-law sources (i.e., X-ray binaries and a putative AGN). We therefore fit the spectra (with background and NGC 253 X-1 contributions subtracted) of all three epochs simultaneously using an apec plasma model to account for the hot gas plus two power-law components that account for the non-variable and variable power-law sources, the latter of which we expect to be primarily due to source B. We note that the X-ray emitting gas is likely to be complex and contain multiple components (e.g., Pietsch et al. 2001; Mitsuishi et al. 2013). However, since these components have their strongest influence at energies ≲1–2 keV, we chose to utilize a single hot plasma to account for the Fe-line emission. Such a component is expected to account for hot gas in the core of the starburst associated with SN remnants (see, e.g., Pietsch et al. 2001). Our adopted model provides a good fit to the data (see Figure 3(a)) by simply varying only the normalization of the power law associated with variable sources like source B. Our best-fit model parameters include an apec plasma with kT = (3.6 ± 0.3) keV, a non-variable unobscured power-law component with Γ_non-var = 2.1 ± 0.1, and an absorbed (N_H = [2.5 ± 0.4] × 10²³ cm⁻²) variable power law with Γ_var = 3.4 ± 0.2. The 7–20 keV flux of the nuclear region was factors of 1.69 ± 0.11 and 1.39 ± 0.08 (1σ errors) times higher in epochs 2 and 3, respectively, than its flux in epoch 1 ($f_{\rm 7\hbox{--}20\, keV}^{\rm epoch 1} = 7.5 \times 10^{-13}$ erg cm⁻² s⁻¹).

Zoom In Zoom Out Reset image size

We also tried accounting for the variability by allowing the variable power-law absorption column, photon index, and normalization components to vary in different combinations. Holding the normalization constant and varying the column density alone provided an acceptable fit and resulted in column densities ≳ 2.9 × 10²³ cm⁻² and Γ_var ≈ 3.7. Holding the normalization and varying only the photon index between epochs did not provide an acceptable fit to the data. Finally, varying the column density, photon index, and normalization components simultaneously provided a good fit to the data, and indicate that the variable component had high absorption columns (N_{H, var} ≳ 2.6 × 10²³ cm⁻²) and steep photon indices (Γ_var ≳ 3.6). None of the above variations improved the quality of the fits to the data over a simple variation in normalization between epochs. This analysis indicates that the 3–40 keV spectral slope is much steeper than that found for typical AGNs (Γ ≈ 1.5–2.2; e.g., Corral et al. 2011; Vasudevan et al. 2013).

As discussed in Section 2 above, the nuclear region of NGC 253 is resolved by Chandra into three bright point sources detected at 2–7 keV within a ≈60 arcsec² region. The NuSTAR PSF is too large to spatially distinguish these sources. However, the NuSTAR and Chandra spectral coverages overlap in the energy range of ≈3–7 keV, allowing us to measure indirectly the levels that each source contributes to the NuSTAR flux variations. With Chandra, we found that only the central source (labeled source "B" in Figures 2(b)–(e)) varied significantly over the three epochs. To estimate the fractional contributions that source B provided to the total NuSTAR extraction region spectra, we measured 4–7 keV count rates from source B and a larger nuclear region using circular apertures of 1.0 arcsec and 15 arcsec, respectively. We find that the 4–7 keV flux from the larger nuclear region was factors of 1.6 ± 0.1 and 1.2 ± 0.1 (1σ errors) higher in epochs 2 and 3, respectively, compared with epoch 1; this is consistent with the NuSTAR observations over the same energy range, which were 1.58  ±  0.07 and 1.33 ± 0.06 (1σ errors), respectively. Given that the shorter 20 ks (0.23 day) Chandra exposures did not cover the entire ≈3–3.5 day periods spanned by the NuSTAR observations (see Figure 1), this result indicates that the nuclear region was unlikely to vary substantially during each of the observational epochs. We verified this result by inspecting NuSTAR light curves of the nuclear region. These light curves were constructed by extracting event lists from a 10 arcsec radius circle centered on the nucleus and measuring the mean count rates of this region in 20 ks bins. We found no signatures of strong variability within each of the three observational epochs. A K-S test revealed that each of the nuclear region light curves were consistent with those expected from constant-intensity sources. We therefore conclude that the nuclear region did not vary substantially on timescales of ≈6–48 hr. Our Chandra observations indicate that source B itself was brighter in epochs 2 and 3 compared with epoch 1 by factors of 5.2 ± 1.0 and 2.6 ± 0.6 (1σ errors), respectively. Once we subtract the contributions of source B to the total Chandra fluxes, contained within the larger NuSTAR extraction region, we find no significant remaining variations in the 4–7 keV non-source B emission between all three epochs, indicating that nuclear region variability observed within the NuSTAR PSF, at least in the 4–7 keV band, are likely entirely attributable to variations in source B.

Assuming the NuSTAR variability is primarily due to variations in source B, we can estimate its spectra by differencing the total nuclear region spectra between bright and faint epochs. Since the variations in the total NuSTAR spectra appear to be consistent with only changes in spectral normalization and not shape (see above), we assume that S_i, the total nuclear region spectrum in epoch i, is described by

$$\begin{equation} S_i \equiv S_{\rm const} + c_i S_B, \end{equation} \tag{ 1 }$$

where the constant S_const represents all the non-variable emission in the nuclear region (e.g., due to sources A and C, as well as diffuse emission), the constant S_B represents the spectrum of source B, and c_i represents the variable normalization of source B at epoch i. Indeed, Figure 3(b) shows the difference spectra, which are well fit by an absorbed broken power-law model with best-fit N_H = (1.6 ± 0.5) × 10²³ cm⁻², low-energy slope of Γ₁ = 2.4 ± 0.5, break energy E_break = 7.9 ± 0.9 keV, and high-energy slope Γ₂ = 3.9 ± 0.4 for both cases, with only the normalization varying. Fitting the difference spectra of epoch 2 minus epoch 1 and epoch 3 minus epoch 1 separately yields similar best-fit parameters and no improvement in spectral fit, implying the spectral shape of source B did not vary significantly between epochs 2 and 3. Equation (1) then implies

$$\begin{equation} c_iS_B = \frac{S_i - S_1}{1 - c_1/c_i}. \end{equation} \tag{ 2 }$$

From the Chandra observations, we know that c₁/c₂ ≈ 0.19 and c₁/c₃ ≈ 0.38. We note that no Fe line was necessary in our fits to source B implying that variable Fe emission is not seen in this source. This result alone likely rules out the possibility that source B is a more luminous version of the variable Fe reflection nebulae seen in the Galactic center, since we would expect corresponding variability in the Fe emission line (e.g., Ponti et al. 2010).

Integration of our spectral model and Equation (2) implies the 2–10 keV fluxes of source B were ≈7.4 and 3.0 × 10⁻¹³ erg cm⁻² s⁻¹, for observations 2 and 3, respectively. These fluxes correspond to observed 2–10 keV luminosities of L_2–10 keV ≈ 1.4 × 10³⁹ and 5.6 × 10³⁸ erg s⁻¹, respectively. The unabsorbed, intrinsic 2–10 keV luminosities are therefore L_2–10 keV ≈ 5.1 and 2.6 × 10³⁹ erg s⁻¹, respectively. We note that the central black hole of NGC 253 has been estimated to be ≈5 × 10⁶ M_☉ (Rodríguez-Rico et al. 2006). Such a black hole would have an Eddington luminosity of $L_{\rm 2\hbox{--}10\, keV}^{\rm Edd} \approx 3 \times 10^{43}$ erg s⁻¹. For this approximation, we assumed a 2–10 keV bolometric correction of ≈22.4, which corresponds to the median bolometric correction for local AGNs with L_2–10 keV ≈ 10⁴¹–10⁴⁶ erg s⁻¹ (Vasudevan & Fabian 2007). If source B were powered by the black hole, then it would be accreting at ∼10⁻⁴ Eddington. AGNs with these levels of L_2–10 keV/L_Edd typically have spectral slopes of Γ = 1.4 ± 0.4 (e.g., Shemmer et al. 2006; Younes et al. 2011), much shallower than the spectrum measured for source B. This suggests that source B is unlikely to be powered by an LLAGN. Instead, the measured luminosities, lack of variability on the ≈5–48 hr timescales, and X-ray spectral shape of source B closely resemble the properties of binaries in the ultraluminous state, which include ULXs that are likely to be stellar-mass black holes accreting above the Eddington limit (e.g., Roberts 2007; Heil et al. 2009; Gladstone et al. 2009). It is therefore likely that source B is a ULX and not an AGN. Given that the number of ULXs in galaxies has been observed to correlate with galaxy-wide SFR (e.g., Mineo et al. 2012), it would not be surprising to find a ULX associated with the nuclear starburst of NGC 253. Indeed, multiple bright point sources and transient ULXs have also been observed in the nuclear starburst in M82 (e.g., Feng & Kaaret 2007). Such a feature may be ubiquitous among starburst galaxies. At distances ≳5–10 Mpc, even Chandra may have difficulty in distinguishing between such sources.

As noted in Section 2 (see also Figure 2(b)), the Chandra position of source B is ∼1 arcsec from the radio sources TH2 and TH4 (from Ulvestad & Antonucci 1997) and the galactic dynamical center (within the 1.2 arcsec 3σ uncertainty radius; Müller-Sánchez et al. 2010). A similar hard source located ≈0.4 arcsec from the dynamical center was previously reported by Müller-Sánchez et al. (2010), based on the Chandra observation from 2003. In Figure 2(e), we show the nuclear region image of the first 20 ks of the 2003 observation (to be equivalent to the depth of our 2012 epochs), with the average locations of the 2012 sources circled. It appears that the position of the 2003 hard near-nuclear source (labeled "N-2003" in Figure 2(e)), located at ($\alpha,\delta)_{\rm J2000} = 00^{\rm h} 47^{\rm m} 33\buildrel{\mathrm{s}}\over{.}12,\,-25^\circ 17^\prime17{\buildrel{\prime\prime}\over{.}} 87$, is offset from source B by ∼1 arcsec (∼19 pc) in the direction of TH2/TH4, and the dynamical center. We note that three other archival Chandra exposures of NGC 253 (ObsID: 790, 969, and 383) were inspected for evidence of N-2003 or source B. Although not formally detected, N-2003 is visually apparent in a 4–7 keV image from the shallower ≈14 ks archival observation (ObsID: 969) that was conducted in 1999 December. The moderately deep 45 ks observation (ObsID: 790), also conducted in 1999 December, had an aim point displaced ≈5 arcmin from the nuclear region, which effectively blended the PSFs making it impossible to spatially measure whether N-2003 or source B were present. Finally, the ≈2 ks exposure conducted in 2000 August (ObsID: 383) was too shallow to detect either N-2003 or source B. We therefore choose to restrict further comparisons of our 2012 observations to the 2003 exposure.

Given the brightness of source B and N-2003 (≈50–300 net counts in the 2–7 keV band) and small off-axis angles with respect to the Chandra aim points (≲0.5–1 arcmin), we expect the 90% confidence positional errors related to PSF centroiding to be ≈0.1–0.2 arcsec (based on Equation (13) of Kim et al. 2007). These positional errors are much smaller than the observed offsets and strongly indicates that the two sources are unrelated. Although the Chandra astrometric frames were aligned to the 2003 observation (see Section 2 for details), some non-negligible error related to image registration is expected. To test further whether the offsets between N-2003 and source B are statistically significant, given our best image registrations, we matched Chandra sources detected in the 2003 observation to counterparts detected in each of the three 2012 observations, excluding source B and N-2003. For our matching, we required that sources in each catalog have ≳40 net counts in the 2–7 keV band and be located within ≈3 arcmin of the nucleus, so that the signal-to-noise ratios and Chandra PSFs are comparable to those of N-2003 and source B. Two sources were considered to match if they were separated by <2 arcsec. In Figure 4(a), we show the right ascension (α) and declination (δ) offsets for all 20 matches and highlight the offsets between the N-2003 and source B in epochs 2 and 3 (source B was not detected in epoch 1). We find that the offsets between N-2003 and source B are 1.10 and 0.99 arcsec for the epoch 2 and 3 positions, respectively; both are much larger than the maximum offset of the 20 matches (0.6 arcsec). In Figure 4(b), we show the cumulative number of matches as a function of offset with a best-fit error function, which describes well the cumulative offset distribution. The data in Figure 4(b) have been normalized by the error function maximum to allow for direct estimates of the confidence intervals. Our error function suggests that the source B offsets are both significant at the ≳99.99% confidence levels. Taken together, we infer that source B is unrelated to N-2003.

Zoom In Zoom Out Reset image size

The above analysis shows convincingly that N-2003 is unrelated to source B implying that at least one of these sources is not associated with the central black hole. Our spectral constraints, presented in Section 3.2, show that source B is most likely to be a ULX. On the other hand, N-2003 is offset from the dynamical center by only ≈0.4 arcsec (i.e., within the 1.2 arcsec 3σ uncertainty of the dynamical center), indicating that it is a better candidate for a "true" nuclear X-ray point source that may be an AGN (see Figure 2(e)). In Figure 5, we show the 3–8 keV Chandra spectra of N-2003 and source B (when at its peak in epoch 2). We note that the morphology of the diffuse 3–8 keV Chandra emission in the nuclear region is complex and appears brightest in the immediate vicinity of N-2003 and source B (see Figures 2(b)–(e)). As such, detailed modeling of the background spectrum associated with the diffuse emission in this region is difficult. To mitigate this limitation, we restricted our Chandra spectral analysis to energies above 3 keV to exclude strong continuum contributions from the hot interstellar medium (ISM) and added an Fe component at 6.7 keV to account for the gas emission line (see Section 3.1 for motivation). Given that the NuSTAR difference spectra of source B do not show a strong Fe line (see Figure 3(b)), we do not expect that this component is intrinsic to source B; however, we are less certain about the nature of the Fe line associated with N-2003. The best-fitting spectral model for source B, presented in Section 3.1, renormalized to the 3–5 keV flux of N-2003, is obviously far too steep to fit the higher-energy 5–8 keV spectrum of N-2003, supporting the conclusion that source B and N-2003 are different, despite having similar 2–10 keV fluxes. As a caveat, we note that X-ray binaries accreting at near-Eddington rates can have variable spectral properties without substantial changes in luminosity (e.g., GRS 1915+105; Vierdayanti et al. 2010), implying that the harder spectrum of N-2003 does not rule out a ULX origin similar in nature to source B. The Chandra spectrum of N-2003 can be fit well using an absorbed power-law model with N_H = (2.8 ± 0.6) × 10²³ cm⁻² and Γ = 1.9 ± 0.6, and an Fe line fixed at 6.7 keV. If we used only the Chandra data to fit the spectrum of source B, we find N_H = (1.9 ± 0.8) × 10²³ cm⁻² and Γ = 3.1 ± 1.0, consistent with the values found from our NuSTAR analysis in Section 3.2. Although the fitting parameters of N-2003 and source B are formally consistent, the spectrum of N-2003 is harder than that of source B (see bottom panel of Figure 5), and more consistent with an AGN (Γ ≈ 1.5–2.2; see further discussion above).

Zoom In Zoom Out Reset image size

Due to its close proximity to the dynamical center and its harder spectrum, it is plausible that N-2003 is an AGN that was in a low state during our 2012 monitoring campaign. However, we cannot currently rule out a ULX nature for N-2003. If N-2003 were an AGN, the Γ ≈ 1.9 power-law spectrum would extend well into the NuSTAR bandpass out to beyond ≈40 keV. The observed 2–7 keV flux of N-2003 is 1.1 times that of source B in epoch 2. However, if N-2003 were an AGN, then we predict that it would have had a 10–20 keV flux ≳2.5 times higher than that of source B at its peak in epoch 2. Therefore, future monitoring of NGC 253 with both Chandra and NuSTAR would be able to resolve the degeneracy between the AGN and ULX nature of N-2003 if caught in a high state.