The Study of X-Ray Flux Variability of M87

A radio galaxy is an active galactic nucleus (AGN) with a strong jet. Some radio galaxies have been detected in the TeV gamma-ray band, indicating that particle acceleration up to TeV energy occurs. The mechanism of particle acceleration up to TeV in the AGN is still under debate. For example, Axford et al. (1977) suggested a shock-in-jet model. Intraday TeV gamma-ray variability found that IC 310 cannot be explained by the shock-in-jet model (Aleksić et al. 2014); thus, other models have been suggested, e.g., the minijet model (Giannios et al. 2009), the cloud-jet interaction model (Barkov et al. 2012), and the magnetospheric model (Blandford & Znajek 1977). Apart from the shock-in-jet model in the core jet region, particle acceleration up to TeV energy in the outer jet region has been studied. A recent discovery of extended TeV gamma-ray emission from the kpc-scale jet of Cen A showed particle acceleration up to TeV in the kpc-scale jet region (H.E.S.S., 2020).

M87 is a nearby FR-i radio galaxy at z ∼ 0.004 (∼16 Mpc) (Jordan et al. 2005), and jet components have been resolved through radio, optical, and X-ray observations. Therefore, it is ideal to study the location of particle acceleration regions up to TeV, and many X-ray and gamma-ray observations have been performed. In 2017, the Event Horizon Telescope team observed black hole shadow, and central black hole mass was estimated to be M_BH = (6.5 ± 0.7) × 10⁹ M_⊙ (The Event Horizon Telescope Collaboration 2019). This mass is more massive than other known black holes (Impey et al. 2006). Gamma-ray emission has been detected in the GeV and TeV bands (Aharonian et al. 2003; Abdo et al. 2009). Daily TeV gamma-ray variability was found by H.E.S.S (Aharonian 2006) and MAGIC (Albert et al. 2008), indicating that the emission region was as compact as the Schwarzschild radius of the M87 supermassive black hole. Aharonian (2006) suggested that daily TeV gamma-ray variability was due to particle acceleration at HST-1, ∼100 pc from the core. Moreover, from later X-ray observations of the core and HST-1 flux variability, the core was suggested for the TeV emission region (Acciari et al. 2010). Therefore, whether HST-1 is a site of particle acceleration up to TeV energy or not is under debate. In this article, to study the location of particle acceleration up to TeV, we searched for short-term X-ray variability from each jet component of M87 from archival X-ray data of Chandra, Suzaku, and NuSTAR. If we find a rapid variability with a steep X-ray spectrum, synchrotron emission by TeV electrons is suggested.

In this study, we assumed that the distance to M87 is 16 Mpc. In Section 2, we describe observations and data analysis. In Section 3, we describe the results of time variability and spectral analysis. Finally, in Section 4, we discuss the origin of intraday X-ray time variability.

To search for X-ray intraday variability, we chose X-ray data that have a longer exposure time of more than 10 ks. For Chandra data, we chose data taken directly by Advanced CCD Imaging Spectrometer (ACIS) and High Resolution Camera (HRC); imaging observations without grating. For ACIS data, we chose data taken with a short frame time of 0.4 s (subarray mode) to minimize pileup effects. As a result, 15 Chandra data (13 ACIS and 2 HRC), 6 NuSTAR data, and 1 Suzaku data were found. There are 6 NuSTAR observations with more than 10 ks, and Chandra observed simultaneously with each NuSTAR observation. Since NuSTAR cannot resolve each jet component, Chandra data was useful and thus these data are analyzed, even though their exposure time is short. Tables 1, 2, and 3 summarize the list of observation data analyzed in this article.

We analyzed Chandra data with CIAO (version 4.12), and Suzaku and NuSTAR data with HEASoft v6.28. CALDB (version 20180925) were used as the calibration data for Suzaku and NuSTAR data analysis. For Chandra data, we checked pileup effects in the 3 × 3 pixels around the core in each observation, following the method of Davis (2001), and we confirmed that a pileup fraction was ≤5% in the core region for all observations analyzed here, except the observation (ID 21076) in the bright state in 2018, therefore the effects are negligible. For this observation (ID 21076), the pileup fraction was around 7.5%. Therefore, this is not so significant for light-curve analysis, but could affect the spectral fitting. and thus we note this issue in Section 3.

2.1. Light Curve

To study a short-term variability, we derived an X-ray light curve in each observation. For Chandra, we derived light curves of each component in the jet: core, HST-1, knot D, and knot A (see Figure 1). The extraction region of each component is shown in Figure 2. The position at 08 from the core is defined as HST-1. In addition, knots D and A are defined in order from the core, and their positions are set at 29 and 124 from the core, respectively. In each observation, the distance from the core for each component was fixed, and the coordinate of the core was adjusted. Each region is defined as a circle with a radius of 20, but a fan-shaped region is excluded when the neighboring component exists in the circular region. Most Chandra data were taken after 2016 when the core was brighter than HST-1. Therefore, HST-1 signals were contaminated by the core signals. To evaluate this effect, Figure 1(b) shows a projected count rate profile along the R.A. direction. The core and HST-1 positions are at 27 and 36 pixels, respectively. The core emission is dominant, and a significant fraction of signals at HST-1 come from the core. Moreover, the effect of HST-1 on the core is insignificant and can be ignored. Therefore, it is difficult to study the HST-1 variability but we will be able to detect the HST-1 variability if there is an independent variation with respect to the core. Chandra light curves were created by the dmextract command in CIAO in the energy range of 2.0–10.0 keV (HRC data was not filtered) with a time bin of 3600 ks. To evaluate the variability of the background signal, light curves of the background region are created in the same way as source light curves, but the position is at 50 from the core in the direction opposite to the jet.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For NuSTAR, we derived light curves of FPMA and FPMB separately in each observation. The source area is defined as a circle with a radius of 30''. A light curve was created with the lcurve command in an energy range of 8.0–30.0 keV with a time bin of 3600 s. As shown in Section 2.2.1, the AGN emission is dominant above 8.0 keV. The light curves of the background region were made with the same radius as the source region at 32 from M87.

The Suzaku satellite observed M87 for two days from 2006 November 29 to December 1, (MJD: 54068–54070). In this study, we used data obtained via X-ray Imaging Spectrometer (XIS). XIS consists of three available CCDs in this observation, XIS0, XIS1, and XIS3. XIS0 and XIS3 are front-illuminated CCDs, whereas XIS1 is a back-illuminated CCD. Therefore, the sensitivity of the low energy range of XIS1 is superior to XIS0 and XIS3 (Koyama et al. 2007). In each instrument, a circle with a radius of 3.0 arcmin was defined as the source region. The light curves of each CCD are created with the lcurve command in an energy range of 3.0–10.0 keV with a time bin of 8400 s. Circles with a radius of 30 at 80 from the core were defined as background regions.

2.2. Spectral Analysis

2.2.1. NuSTAR and Chandra

There are six observations taken by both NuSTAR and Chandra simultaneously. Following the method of Wong et al. (2017), we performed simultaneous fitting of NuSTAR and Chandra spectra to obtain the photon index and flux of each component. In the simultaneous fitting, one NuSTAR spectrum and three Chandra spectra of the core, HST-1, and jet fitted with the model of hot gas emission were considered for the NuSTAR spectrum, and each power-law model was considered for the NuSTAR spectrum and one Chandra spectrum of each component. Power-law parameters were common between NuSTAR and Chandra.

NuSTAR spectra were extracted in a circular region with a radius of 30'', and background regions were taken at 80''–130'' from M87. Spectra were binned so that 1 bin contained more than 20 photons. For Chandra spectra of each component, the source region of core and HST-1 were defined as a circle with a position specified in Figure 1(a) and a radius of 04. For the background of core and HST-1 spectra, a semicircular band region at 2''–4'' was specified in the direction opposite to the jet (Figure 1(a)). The jet area, including knots D and A, is defined as shown in Figure 1(a), and the background region was specified as 195 ×3'' boxes on both sides of the jet. We created Chandra spectra with the specextract command and binned the spectra so that 1 bin contained more than 20 photons.

The spectral model for simultaneous fit consists of two temperature thermal plasma and three power-law models with Galactic absorption: wabs∗(vapec+vapec+pegpwrlw+pegpwrlw+pegpwrlw) in XSPEC. The normalization of vapec was set to 0 for Chandra spectra, and vapec parameters for NuSTAR were fixed as described later. A column density of Galactic absorption was fixed to 1.94 × 10²⁰ cm⁻² (Kalberla et al. 2005). Three pegpwrlw models represent emission from the core, HST-1, and jet. Thus, for each Chandra spectrum of the core, HST-1, and jet, only one pegpwrlw normalization was left free and other normalizations were set to 0. For these analyses, photon index and normalization were unified between NuSTAR and Chandra. Considering cross-calibration is uncertain between FPMA and FPMB of NuSTAR, we added a constant model. The constant was fixed to 1 for FPMA and Chandra spectra, whereas the constant was left to be free for FRMB. The vapec parameters for NuSTAR spectra were determined as follows. Since the extraction regions of the NuSTAR spectra are as wide as 30'', we must consider X-rays from the thermal plasma of clusters of galaxies. However, since thermal emission is dominant in the soft X-ray band, it is difficult to estimate them only by the NuSTAR spectrum. Therefore, first, we estimate the thermal emission using Chandra data. Among observations of M87 by Chandra/ACIS, a full-array observation (2000-07-29; ID352) was selected. We extracted a spectrum within 30'' of M87, excluding the core and jet region. For the background region, we specified a ring shape with an inner and outer radius of 80'' and 130'', respectively. We fitted this spectrum with two apec models, following Russell et al. (2015), but instead of apec we used vapec. A redshift was fixed to 0.0043 (NASA/IPAC 1999). The metal abundances of elements other than He, C, and N were free, but treated as common parameters between two vapec models. He, C, and N abundances were fixed to 1 solar. This model fits the spectrum well with a reduced chi-square of 1.28. As a result, temperatures were determined as kT = 1.66, 0.98 keV, and normalizations were determined as 5.10 × 10⁻³ and 8.68 × 10⁻⁴. Thus, the obtained parameters were used as fixed parameters in the simultaneous fitting of NuSTAR and Chandra spectra. An example of simultaneous fittings is shown in Figure 3. The spectra were fitted well.

Zoom In Zoom Out Reset image size

For reference, we performed spectral fitting for only the NuSTAR spectrum with wabs∗(vapec+vapec+pegpwrlw) and obtained a photon index and flux in 10.0–40.0 keV.

2.3. Suzaku

Spectra of Suzaku/XIS was extracted for the source region with a radius of 3' (Figure 4(a)). Response files were made by the xisarfgen and xisrmfgen commands.

Zoom In Zoom Out Reset image size

We fitted the spectra in 1.0–10.0 keV with wabs∗(vapec+vapec+pegpwrlw). The metal abundances of H, C, and N were fixed to 1, and the abundance parameters of other metals are unified between the two vapec models. Figure 4(b) shows a result of spectral fitting. As a result, temperatures became 2.45 and 1.53 keV. The spectra were well fitted with a reduced χ square of ${\chi }^{2}/{\rm{d}}.{\rm{o}}.{\rm{f}}=1.25$. The photon index became $\gamma ={2.38}_{-0.04}^{+0.07}$ (2.0–10.0 keV), and the power-law flux became ${2.71}_{-0.18}^{+0.18}\times {10}^{-11}$ erg cm⁻² s⁻¹.

To evaluate a short-term time variability, we calculated a rms from X-ray light curves of each observation. The rms is defined as the following equation (Shirai et al. 2008):

$\begin{eqnarray}&&\mathrm{RMS}=\sqrt{\tfrac{{\sum }_{i}[{({d}_{i}-\bar{d})}^{2}-{e}_{i}^{2}]}{N{\bar{d}}^{2}}}\end{eqnarray} \tag{ 1 }$

where d_i is a count rate of ith data, N is the number of data points, $\bar{d}$ is an average, and e_i is a statistical error of d_i . Table 1 summarizes the rms values of each observation. The rms of observations with ${\sum }_{i}[{({d}_{i}-\bar{d})}^{2}-{e}_{i}^{2}]\lt 0$ is not shown, indicating no variability. As a result, some data showed significant variability. Figure 5 shows light curves only for observations with a nonzero rms with 2σ significance.

Zoom In Zoom Out Reset image size

Intraday variability is seen in Suzaku observation in 2006 (Figure 5(f)) and Chandra observation in 2017 (ID:18612, Figure 5(b)). For the Suzaku light curve, count rate changes by approximately 30% with a timescale of 0.3 days. There is a difference in count rates between XIS0 and XIS3, which is attributable to the off-axis observation (JAXA/ISAS & NASA/GFSC 2015). The Chandra light curve of the core in 2017 showed approximately 40% variability with a timescale of 0.5 days. However, since the extraction region of Chandra data was as small as 20, the fluctuation of the satellite attitude may cause variability of the count rate. Therefore, we made light curves with different extraction radii of 04, 10, 20, and 40 to evaluate the systematic effects. From this analysis, rms became 0.10 ± 0.01, 0.09 ± 0.01, 0.08 ± 0.01, and 0.06 ± 0.01 for 04, 10, 20, and 40, respectively. The rms decreased as the radius increased, but variability was significant, even for a large radius. Figure 5(c) shows a significant variability of HST-1 in 2017. However, the variability was quite similar to that of the core, and thus the variability was due to the core signals. The light curves of knot A in 2016 and 2017 showed an insignificant variability pattern (Figure 5(a) and (d)). The NuSTAR light curve in 2016 showed an insignificant variability pattern because two instruments showed a different variability pattern (Figure 5(e)).

Therefore, variability was detected from the Suzaku X-ray observation data in 2006 and Chandra observation of the core data in 2017.

For Suzaku and Chandra data with significant variability, we analyzed the X-ray spectra as described in Section 2 to obtain a photon index and flux of the power-law component. As a result, a power-law photon index became ${2.38}_{-0.04}^{+0.07}$ and ${1.96}_{-0.04}^{+0.05}$, and a flux in 2.0–10.0 keV became ${2.71}_{-0.18}^{+0.18}\times {10}^{-11}$ erg cm⁻² s⁻¹ and ${1.86}_{-0.16}^{+0.17}\times {10}^{-12}$ erg cm⁻² s⁻¹ for Suzaku in 2006 and Chandra core in 2017, respectively. These fluxes were consistent with a long-term light curve of core and HST-1 (Yang et al. 2019). The photon index was different between 2006 and 2017, inferring a different emission mechanism.

We did not find any short-term significant variability for NuSTAR data, but found that one of the six observations showed a very different X-ray flux. As described in Section 2.2.1, we fitted Chandra and NuSTAR spectra simultaneously or only the NuSTAR spectrum of each observation, and obtained a power-law photon index and flux. Figure 6 shows the 10.0–40.0 keV X-ray light curves obtained from the NuSTAR spectrum and the 2.0–10.0 keV light curves for the core, HST-1, and jet obtained with a simultaneous fitting of Chandra and NuSTAR spectra. Remarkably, the flux in 2018 was 2.5 times higher than other observations. Chandra core flux was also higher in the same observation, but the fluxes of HST-1 and jet were almost constant. Wong et al. (2017) reported a flux of ${7.7}_{-1.0}^{+1.1}\times {10}^{-13}$ erg cm⁻² s⁻¹ and a photon index of ${2.12}_{-0.13}^{+0.12}$ (20.0–40.0 keV) for the NuSTAR data in 2017. Our result of the flux of ${9.1}_{-1.3}^{+1.4}\times {10}^{-13}$ erg cm⁻² s⁻¹ and the photon index of ${2.20}_{-0.11}^{+0.11}$ (20.0–40.0 keV) were consistent with those of Wong et al. (2017). As described in Section 2, the Chandra core in the brightest state in 2018 was affected by pileup. The pileup fraction is around 7.5%, and thus the flux could be reduced by at most 7.5%. In addition, this could make the photon index smaller by ∼0.2, based on the Chandra ABC guide to Pileup. ¹ However, these effects are within the statistical errors, and thus our result and discussion do not change significantly.

Zoom In Zoom Out Reset image size

The results obtained in this study are summarized below. First, the Suzaku observation in 2006 and the Chandra core observations in 2017 showed significant intraday X-ray variability. In the Suzaku light curve, the variability amplitude was approximately 30% with a timescale of about 0.3 days. For Chandra, the amplitude was approximately 40% with a timescale of 0.5 days. In 2006, the HST-1 was much brighter than the core in the X-ray band (Yang et al. 2019). Moreover, in 2017 the core was brighter than the HST-1. The spectrum obtained by Suzaku in 2006 showed a steep photon index of ${2.38}_{-0.04}^{+0.07}$, and the core obtained by Chandra and NuSTAR in 2017 had a hard photon index of ${1.96}_{-0.04}^{+0.05}$. The minimum size ${R}_{\min }$ of the emission region was estimated by the X-ray variability timescale Δt as follows:

$\begin{eqnarray}&&{R}_{\min }=\displaystyle \frac{c{\rm{\Delta }}t{\delta }_{D}}{1+z}\end{eqnarray} \tag{ 2 }$

where δ_D is a Doppler factor and z is a redshift. We found Δ t ∼ 0.3 day in 2006, and thus we obtain R_min ≃ 1.0δ_D  × 10¹⁵ cm. Here, a central black hole mass was assumed to be M_BH =(6.5 ± 0.7) × 10⁹ M_⊙ (The Event Horizon Telescope Collaboration 2019). The Schwarzschild radius is R_s  ∼ 2.0 × 10¹⁵ cm. Thus, the size of the emission region was as compact as about ∼R_s . Suzaku/XIS cannot resolve each component of the M87 jet, but knots D and A were quite fainter than the core and HST-1. In addition, X-ray spectra of knots D and A were reported to be soft with a photon index of 2.4–2.6, ≃2.2, respectively (Sun et al. 2018), and thus it is difficult for the knots D and A to explain a rapid variability above 3.0 keV.

From the result of Chandra on 2006 November 13, a count rate ratio of HST-1 to the core was approximately 4.45 (Yang et al. 2019). If the variability was only due to the core, a count rate of the core must vary with an amplitude of ∼260%. Such an extremely large amplitude of variability has never been observed for the M87 core, and thus it is unlikely.

Hubble Space Telescope observed M87 on 2006 May 23, and November 28. The ultraviolet flux of the core and HST-1 increased by a factor of 1.7 and 1.4, respectively, from May to November (Madrid 2009). Notably, Suzaku observation was performed in 2006 November. Radio observations were performed by VLBI in 2006 May, 2007 January, and 2007 April, and flux ratio of HST-1 in January was higher by a factor of 1.5 than that in 2006 May, but it decreased in 2007 April (Chang et al. 2010). Moreover, other components of the M87 did not show variability, indicating that HST-1 was active in 2006 and likely the origin of intraday X-ray variability. The X-ray spectrum in 2006 was steep with a photon index of 2.38. De Jong et al. (2015) showed that X-ray emission observed by Suzaku could be explained by a high-energy tail of synchrotron, based on the spectral energy distribution (SED) modeling, and thus concluded that HST-1 is likely the origin of TeV gamma-rays. These results implied that X-ray intraday variability in 2006 observed by Suzaku/XIS can be explained by synchrotron radiation from HST-1.

Assuming that X-ray decay was due to electron cooling, we could estimate a magnetic field strength. The cooling timescale was expressed by the following equation from Tashiro et al. (1995):

$\begin{eqnarray}&&{\tau }_{s}=3.2\times {10}^{4}\times {B}^{-3/2}\times {E}_{\mathrm{ph}}^{-1/2}\times {\delta }^{-1/2}({\rm{s}})\end{eqnarray} \tag{ 3 }$

where τ_s is the decay timescale, B is the strength of the magnetic field, E_ph is the photon energy in the observer frame, and δ is the Doppler factor. Using a decay timescale of 0.3 days, and E_ph = 6500 eV, and a Doppler factor δ, the magnetic field became B = 1.94δ^1/3 mG. This value was consistent with that in Park et al. (2019), which was based on the study of Faraday rotation. However, it was quite smaller than 0.6 G obtained by de Jong et al. (2015), which was based on the SED modeling of HST-1. This inconsistency might be due to the SED modeling in de Jong et al. (2015) using TeV gamma-ray flare data between 2003 and 2006. If the magnetic field was increased in only the TeV gamma-ray flaring period, by shock compression, our result could be explained. Using the estimated magnetic field, electron energy emitting synchrotron X-ray emission should be 8.7 ${\delta }^{-1/6}$ TeV, and thus particle acceleration up to TeV was indicated to occur in the HST-1.

Hardcastle et al. (2015) mentioned that the hotspot of radio galaxy Pictor A, which was 160 kpc from the core, showed yearly variability. However, the hotspot of Pictor A was further away from the core than HST-1, indicating that particle acceleration could occur far away at sub-kpc or kpc from the central black hole. Such particle acceleration could be in the recollimation zone of the jet (Stawarz et al. 2006).

What about an X-ray variability of the core with a timescale of Δt ∼ 0.5 days in 2017? VERITAS observations in 2008 showed that the X-ray flux of the core simultaneously varied with TeV gamma-ray flux, suggesting that the core X-ray emission originated from the jet as well as TeV gamma-ray (Acciari et al. 2008). Therefore, the 2017 variability could be due to jet emission at the core. The power-law photon index of the core X-ray spectrum in 2017 is ${1.96}_{-0.04}^{+0.05}$, which could be the inverse Compton scattering component rather than synchrotron radiation if it comes from the jet. Other possibilities of the origin of harder X-ray emission besides the jet include an advection-dominated accretion flow (ADAF), which had been highlighted in previous studies (e.g., Harris et al. 2011 and Yuan & Narayan (2014). Let us consider whether the observed variability timescale of 0.5 day in 2017 can be explained by ADAF. It is believed that X-ray emission of low luminosity AGN like M81 comes from ADAF. The mass of the supermassive black hole in M81 was ${7.0}_{-1}^{+2}\times {10}^{7}$ M_⊙ (Schödel et al. 2007), and the X-ray variability timescale was confirmed to be about one day (La Parola et al. 2004). Therefore, in the case of M81, the size of the radiation region is 25 R_s , which is much larger than ∼R_s of M87 estimated in this study. The origin of X-ray emission of M81 was investigated by Quataert et al. (1999) who modeled the multiwavelength SED. They concluded that the ADAF model is promising. The size of 25 R_s for the emission region was consistent with the ADAF view. In other words, the intraday variability of the M87 core leads to the emission region size of ∼R_s and thus ADAF origin is unlikely. Accordingly, it is inferred that the intraday X-ray variability of the core in 2017 is due to jet emission. Intraday variability of the core implies a smaller emission region than the photon ring (6 R_s ) (The Event Horizon Telescope Collaboration 2019). Future observations by the Event Horizon Telescope may reveal a radio variability of the photon ring region.

In the six NuSTAR and Chandra observations in 2017–2019, the flux of the core was high by a factor of 2.5 in only one observation of 2018 April. The core emission in this period is considered as jet inverse Compton or ADAF radiation rather than jet synchrotron radiation, as discussed above. Since the variability timescale is constrained to be longer than 1 day, we cannot discuss the origin of variability any further.

We have intraday X-ray variability of M87 from Chandra, NuSTAR, and Suzaku archival data. As a result, Suzaku observation in 2006 and Chandra core observation in 2017 showed significant intraday variability with a timescale of 0.3 days and 0.5 days, respectively. A much brighter X-ray flux of HST-1 than the core in 2006 indicates that the variability comes from HST-1.

The X-ray spectrum of M87 in 2006 was steep with a photon index of ${2.38}_{-0.04}^{+0.07}$, but the core in 2017 was hard with a photon index of ${1.96}_{-0.04}^{+0.05}$. The steep spectrum of HST-1 seems to be a high-energy tail of synchrotron emission. Assuming the decay timescale corresponds to a synchrotron-cooling timescale, the magnetic field strength is calculated to 1.94δ^1/3 mG, yielding the electron energy of $1.1\times {10}^{8}{\gamma }^{-1}{\delta }^{-1/3}$ TeV, which indicates that particle acceleration up to TeV occurs at subkiloparsec scale from the central black hole.

The core hard spectrum in 2017 indicates the possibility of ADAF emission or inverse Compton in the jet, but we obtain minimum emission size as compact as the Schwarzschild radius from the variability timescale, and thus ADAF origin is not preferred. From these results, we suggest that particle acceleration occurs in both the core and subkiloparsec jets in M87.

This research has made use of data and software provided by the High Energy Astrophysics Science Archive Research Center (HEASARC), which is a service of the Astrophysics Science Division at NASA/GSFC.