X-Ray Spectral Analysis of the Jet Termination Shock in Pictor A on Subarcsecond Scales with Chandra

Relativistic jets launched from high-accretion rate active galactic nuclei (AGNs), such as quasars and high-excitation radio galaxies, terminate by forming powerful shock waves, observed as prominent hot spots at the edges of extended radio cocoons/lobes inflated by the jets in the ambient medium (Blandford & Rees 1974; Scheuer 1974). In more detail, a light but high-power relativistic jet, when interacting with much denser interstellar/intergalactic medium, forms a double-shock structure: the nonrelativistic forward shock propagates within the surrounding gas, compressing and heating the thermal plasma (see, e.g., Carilli & Barthel 1996; O’Sullivan et al. 2018), while the relativistic reverse shock converts the bulk kinetic energy of the outflow to the internal energy of jet particles (e.g., Meisenheimer et al. 1989; Kino & Takahara 2004). Magnetic field amplification and acceleration of some fraction of the jet particles to high, and even ultrahigh energies, is expected to take place at the front of the reverse shock as well, although the exact acceleration processes, or the efficiency of the magnetic amplification, are still under the debate (e.g., Stawarz et al. 2007; Fan et al. 2008; Araudo et al. 2016, 2018; Matthews et al. 2019).

Hot spots in cosmologically distant radio quasars and high-power FR II radio galaxies are typically of the size of a few/several kiloparsecs, and so in order to study them properly, one needs instruments with at least arcsecond resolution. A considerable effort was made to resolve such structures at radio and infrared/optical frequencies, where hot spots shine through the synchrotron emission downstream of the reverse shock (e.g., Prieto et al. 2002; Brunetti et al. 2003; Mack et al. 2009; Perlman et al. 2010; Orienti et al. 2012, 2017, 2020; Pyrzas et al. 2015; Dabbech et al. 2018; Migliori et al. 2020; Sunada et al. 2022a). Hot spots are also the sources of nonthermal X-ray photons, as established by numerous Chandra observations (Hardcastle et al. 2004; Kataoka 2005; Tavecchio et al. 2005; Harris & Krawczynski 2006; Massaro et al. 2011, 2015; Mingo et al. 2017). The origin of the X-ray hot spots’ emission is, in many cases, unclear: while in some sources the X-ray spectrum seems to fall into the extrapolation of the radio-to-optical synchrotron continuum, in other sources the X-ray excess suggests an additional emission component, typically ascribed to inverse-Comptonization of cosmic microwave background photons, or of the hot spot’s own synchrotron photons, by lower-energy electrons.

Among the other targets, the western (W) hot spot in the radio galaxy Pictor A is exceptionally well suited for deep observational studies, due to the combination of its relatively large angular size, very large angular separation from the bright galactic nucleus, and its high surface brightness. As such, it has been subjected to a number of multiwavelength programs, including the radio domain with the Very Large Array (VLA; Perley et al. 1997), the mid-infrared range with the IRAC camera on board the Spitzer Space Telescope (Werner et al. 2012), the Wide-field Infrared Survey Explorer (Isobe et al. 2017), and the SPIRE camera on board the Herschel Space Observatory (Isobe et al. 2020), at optical wavelengths with the Faint Object Camera on the Hubble Space Telescope (HST; Thomson et al. 1995), in X-rays with the Advanced CCD Imaging Spectrometer (ACIS) on board the Chandra X-ray Observatory (Wilson et al. 2001; Hardcastle et al. 2016; Thimmappa 2020), as well as the EPIC MOS1 camera of XMM-Newton (Migliori et al. 2007), and lastly in hard X-rays with NuSTAR (Sunada et al. 2022b). The hot spot has also been the target of high-resolution radio imaging by the Very Long Baseline Array (VLBA; Tingay et al. 2008).

The radio structure of the W hot spot at GHz frequencies with the 15 VLA resolution is complex, including the main compact knot at the westernmost edge of the system, and the diffuse plateau region extending to the east/southeast (Perley et al. 1997). With subarcsecond VLA resolution (reaching 017), the main knot remains unresolved, while the plateau region reveals distinct filaments. The 74 MHz–5 GHz spectral index of the main knot is α ≃ 0.6–0.7, and the degree of polarization reaches 70%; the upstream filaments seem to be characterized by a steeper spectrum (Δα ≳ 0.1) and decreased polarization level (down to 10%–30%). The projected magnetic field aligns with the levels of constant radio brightness, such that if the main compact knot denotes the position of the terminal reverse shock and its near downstream, the magnetic field configuration corresponds to that of a perpendicular shock. On the optical HST image with ≃01 resolution (Thomson et al. 1995), the main knot is decomposed into a system of highly polarized (≳50%) wisps elongated perpendicular to the jet axis.

Such complexity can hardly be followed at X-ray frequencies even with Chandra's superb resolution. However, after a careful ACIS image deconvolution with subpixel resolution, presented in Thimmappa (2020), the W hot spot could, in fact, be resolved into (i) a disk-like feature perpendicular to the jet axis, located ≃15 to the southeast of the intensity peak of the main radio knot, but coinciding with the peak of the hot spot’s optical emission, and (ii) an elongated feature aligned with the jet axis, and located even further upstream, i.e., within the region of the radio plateau. The disk-like feature could be traced for ∼4″ in its longitudinal direction, but is resolved in its transverse direction only on subpixel scale.

The overall interpretation of the observed multiwavelength morphology of the W hot spot therefore emerges, in which the perpendicular disk-like structure at the position of the hot spot’s optical and X-ray intensity peaks corresponds to the very front of the reverse shock, where the most efficient particle acceleration is expected to take place. The radio intensity peak located further away, on the other hand, marks the downstream of the reverse shock, where the radiative cooling of the plasma convected away from the shock front prevents the production of high-energy optical and X-ray synchrotron photons. Finally, the nature of the X-ray jet-like feature upstream of the shock, as well as optical and radio filaments within the extended plateau region, remain unclear, although such structures may be related to a network of weaker oblique shocks formed around the head of the jet by the plasma backflowing from the downstream of the reverse shock (see, e.g., Saxton et al. 2002; Mizuta et al. 2010).

The good agreement between the optical and X-ray maps, along with the general X-ray spectral properties of the hot spot, as well as hints of the X-ray time variability of the target, all imply in accord with a synchrotron origin of the observed X-ray photons (Hardcastle et al. 2016; Thimmappa 2020; Sunada et al. 2022b). For this reason, the X-ray spectral properties of the hot spot are crucial for a proper understanding of particle acceleration processes taking place at mildly relativistic perpendicular shocks in general. And indeed, the very presence of X-ray synchrotron photons means that such shocks are able to accelerate electrons up to energies E_e ∼ 10⁸ m_e c², assuming the hot-spot magnetic field of the order of B ∼ 0.1–1 mG (see the discussion in Thimmappa 2020; Sunada et al. 2022b).

However, for an X-ray spectral analysis with any of the available X-ray instruments, the source extraction region has to be relatively large, in order to maximize the photon statistics for a given point-spread function (PSF) and the source intrinsic extent. Such an extended region unavoidably includes therefore various subcomponents of the system, and hence the resulting spectral constraints do not correspond to the reverse shock exclusively, but instead to a superposition of the reverse shock, its downstream region, and also of the upstream filaments/jet-like features. Here, we propose a novel, alternative method for constraining the shape of the X-ray continuum emission at the very position of the reverse shock, with subarcsecond resolution. The method is based on hardness map analysis, for separately deconvolved soft and hard maps; this novelty resolves the problem of artifact features appearing on X-ray hardness maps due to the energy-dependent Chandra PSF.

Throughout this paper, we assume a Lambda cold dark matter cosmology with H₀ = 70 km s⁻¹, Ω_m = 0.3, and Ω_Λ = 0.7. The Pictor A redshift z = 0.035 (Eracleous & Halpern 2004), corresponds therefore to the luminosity distance of 154 Mpc, and the conversion angular scale 0.7 kpc arcsec⁻¹. The photon index Γ is defined here as F_ε ∝ ε^−Γ for the photon flux spectral density F_ε and the photon energy ε; the spectral index is α = Γ − 1.

The W hot spot of Pictor A was observed on-axis with the ACIS (Garmire et al. 2003) on board the Chandra X-ray Observatory (Weisskopf et al. 2000) on 2002 September 17 (ID 3090) and 2002 September 22 (ID 4369). A combination of relatively long uninterrupted exposures for both pointings, totaling an observing time of 95.5 ks, and a small off-axis angle , makes them ideal data set for our high-resolution study.

The observational data were reprocessed using the chandra_repro script as per CIAO v4.14; Fruscione et al. 2006) analysis threads, ⁴ using CALDB v4.9.7. Pixel randomization and readout streaks were removed from the data during processing. Point sources in the vicinity of the hot spot were detected with the wavdetect tool using the minimum PSF method, and removed. For our analysis, we selected photons in the range of 0.5–7.0 keV. Photon counts and spectra were extracted for the source and background regions from individual event files using the specextract script. Spectral fitting was done with the Sherpa package (Freeman et al. 2001).

The total number of counts for the hot spot, ∼10,000 for both exposures together (see Table 1), places us in the regime where calibration uncertainties dominate over statistical uncertainties (Drake et al. 2006). Methods to account for calibration uncertainties in the analysis of Chandra data have been discussed by Lee et al. (2011) and Xu et al. (2014). The moderate count rate of ≃0.1 s⁻¹ for the hot spot located at the center of the S3 chip, implies only small chances of a pileup in the detector (Davis 2001); we have verified this during the spectral analysis, but nonetheless have included the pileup model when performing MARX simulations anyway (see the following section).

3.1. Spectral Modeling

A composite hot-spot spectrum was extracted for each ObsID from a circular region (position: R.A. = 5:19:26.2993, decl. = –45:45:54.377) with a radius of 20 px (≃10″, for the conversion scale 0492 px⁻¹), and the background set as a concentric annulus of 30–60 px radius (see Thimmappa 2020, Figure 1 therein). The background-subtracted hot-spot spectra were next fitted within the soft (0.5–2.0 keV) and hard (2.5–7.0 keV) bands separately, assuming a power-law model modified by the Galactic column density N_H,Gal = 3.6 × 10²⁰ cm⁻² (HI4PI Collaboration et al. 2016). The results of spectral fitting are summarized in Table 1, and the fitted spectra are shown in Figure 1.

Zoom In Zoom Out Reset image size

The power-law models with photon indices Γ ≃ 1.9 in the soft band, and significantly larger Γ ∼ 2.7 in the hard band, provide a reasonable description of the source composite spectra, sufficient in particular for the purpose of the PSF modeling. We note that analogous fits with the Galactic absorption set free returned similar results, only with slightly decreased values of N_H,Gal and Γ. Finally, including the jdpileup model in the fitting procedure does not affect the best-fit values of the model parameters, as the fraction of piled-up events that result in a good grade turns out to be very low.

3.2. PSF Modeling

To model the Chandra PSF at the position of the W hot spot, we used the ChaRT online tool (Carter et al. 2003) ⁵ and the MARX software (Davis et al. 2012). ⁶ For both ObsIDs 3090 and 4369, the centroid coordinates of the selected source region were taken as the position of a point source. The source spectra for ChaRT were the respective power-law models in the 0.5–2.0 and 2.5–7.0 keV bands, as described in Section 3.1. Since each particular realization of the PSF is different due to random photon fluctuations, in each case a collection of 50 event files was made, with 50 iterations using ChaRT by tracing rays through the Chandra X-ray optics. The rays were projected onto the detector through MARX simulation, taking into account all the relevant detector effects, including pileup and energy-dependent subpixel event repositioning. The PSF images were created with the size of 32 × 32 pix², and binned with 0.5 px resolution. An example of the simulated PSF images for ObsID 3090 in the soft and hard bands is presented in Figure 2.

Zoom In Zoom Out Reset image size

In order to illustrate the size of the PSFs in both bands, we calculated the enclosed count fraction (ECF) for all the simulated PSFs, i.e., the fraction of counts that would be detected within a certain circular aperture for a particular realization of the PSF. The resulting ECFs are presented in Figure 3 for the soft and hard bands (left and right panels, respectively), ObsID 3090. As shown, for the soft band, the 2σ radius is typically ≃6 px, while in the hard band the corresponding 2σ radius has a wider spread, ranging from ≃6 px for some realizations of the PSF, up to even ≃15 px.

Zoom In Zoom Out Reset image size

Note that, since the region encompassing the hot-spot structure is relatively compact, one should not expect any significant change of the PSF across the field subjected to the image deconvolution, as described in the next section. The spectral information provided for the PSF modeling, on the other hand, corresponds to the composite radiative output of the entire structure, while below we argue for the presence of significant spectral changes on subpixel scale within the brightest segments of the hot spot. This inherent inconsistency does not however affect the main results of the analysis.

3.3. Image Deconvolution

We used LRDA, which is implemented in the CIAO tool arestore, to remove the PSF blurring, and in this way to restore the intrinsic surface brightness distribution of the hot spot. This method does not affect the number of counts on the image, but only their distribution.

The algorithm requires an image form of the PSF, which is provided by our ChaRT and MARX simulations as described in Section 3.2 above, and exposure-corrected maps of the source (for more details see Marchenko et al. 2017; Thimmappa 2020). Here, we perform the deconvolution separately for the soft and hard bands, in each case for 50 random realizations of the simulated PSF; those 50 deconvolved images were then averaged to a single image using the dmimgcalc tool. The resulting images are shown in Figure 4. The two main features of the hot spot observed by Thimmappa (2020), namely, the disk-like feature perpendicular to the jet axis, as well as a weaker jet-like feature extending to the southeast along the jet axis, are present in both soft and hard maps for both ObsIDs, although the jet-like feature is much less prominent on the hard maps.

Zoom In Zoom Out Reset image size

3.4. Hardness Ratios

Based on the deconvolved soft and hard images of the W hot spot in Pictor A with 0.5 px resolution, we performed a spatially resolved hardness ratio (HR) analysis, in order to investigate the spectral structure of the system on subarcsecond scales, free, as much as possible, from the PSF blurring.

HR analysis of Chandra data has been widely applied to various classes of astrophysical sources before (e.g., Bałucińska-Church & Ostrowski 2005; Siemiginowska 2007; Nandra et al. 2015; Haggard et al. 2019), in particular being considered as a useful tool that allows constraints on spectra of unresolved weak sources, for which the standard spectral modeling approach is not possible due to low numbers of counts. Spatially resolved HR analysis for extended sources, however, remains largely unexplored, because of artifact features appearing on the HR maps, in relation to (i) the energy dependence of Chandra's PSF, and also (ii) random fluctuations of photons, relevant especially in the low surface brightness regime. Our approach resolves the aforementioned problems, since (i) the HR analysis is based on the separately deconvolved soft and hard maps, and (ii) we remove the effect of random fluctuations by averaging over 50 realizations of each modeled PSF. In particular, based on the soft (S) and hard (H) deconvolved images, we produce spectral maps defined as HR = (H − S)/(H + S). Next, we average the HR maps corresponding to ObsIDs 3090 and 4369, obtaining at the end the final distribution of the X-ray HR for the W hot spot of Pictor A, shown in the left panel of Figure 5.

Zoom In Zoom Out Reset image size

An HR variance map, i.e., the values of the variance at a given position (x,y) on the map, was generated based on the same N = 100 deconvolved HR images (50 for ObsID 3090 and another 50 for ObsID 4369), HR_i (x, y), simply as

where denotes the averaged HR image. The square-root of this variance, , corresponds therefore to the statistical uncertainty in the derived values of the HR at a given position (x,y) on the map. This uncertainty is shown in the right panel of Figure 5.

The main structure of the hot spot that is prominent on the total intensity map is characterized by the values −1 < HR < 0. This structure is surrounded by a soft halo with HR = −1, meaning simply no hard photons; outside of the soft halo, where the X-ray flux also drops in the soft band, the HR values fluctuate around 0. This assures the reality of the spectral map produced, as no artifact features are present in the regions with background-level flux, and all the physically meaningful HR values are concentrated exclusively in the high-flux region. Moreover—and this is the major finding of the analysis—there is a clear systematic HR gradient across the main disk-like feature, ranging from approximately −0.4 down to −0.9 and below (see the left panel of Figure 5) across the main disk-like feature, i.e., from the upstream (southeast) to the downstream (northwest) of the shock. The HR uncertainty in that region is on average, ±0.2 (see the right panel in Figure 5), so that the HR gradient is statistically significant.

In addition to the statistics, however, a careful investigation of the systematic uncertainty is also required. We have therefore produced hardness maps of other astrophysical sources appearing pointlike for Chandra, in the analogous way as described above. For a fair comparison with the Pictor A hot spot, we selected sources that are unresolved and were observed by Chandra around 2002 (in order to avoid complications related to the CCD degradation), were not variable during the Chandra exposure, were free of pileup, and had comparable photon statistics to those of the analyzed Pictor A pointings. The best targets fulfilling such criteria, were the BL Lac object AO 0235+16 (z = 0.94), and quasar 4C+13.85 (z = 0.673). In both cases, we found no evidence of any substructure introduced by the HR map. Thus, we do not believe that the gradient seen above is a systematic effect of our method. The corresponding maps for the two targets are presented in the Appendix.

Assuming a single power-law emission model, a given value of the HR corresponds to a particular set of values for the photon index Γ and Galactic column density N_H,Gal (assuming zero intrinsic absorption). In Figure 6, left panel, we plot this dependence, adopting N_H,Gal = 3.6 × 10²⁰ cm⁻². With such, the value HR = −0.4 gives the photon index Γ ≃ 1.2, while, for example, HR = −0.9 gives Γ ≃ 2.8. The resulting 0.5–7.0 keV photon index map of the W hot spot in Pictor A, corresponding to the 0.5 px resolution HR map discussed above, is shown in the right panel of Figure 6. The uncertainties in the exact N_H,Gal value, even if at the level of ∼50%, are in this context much less relevant than the statistical HR mean uncertainty of ≃0.2, following from the square-root variance mapping of the disk feature. This statistical uncertainty would in particular translate into a wider range of the allowed photon indices, roughly speaking Γ ≤ 1.6 for the upstream edge, and Γ ≥ 1.9 for the downstream region. In the case of the synchrotron origin of the detected X-ray photons, those values of photon indices would then correspond to the index of the electron energy distribution , ranging from ≤2.2 up to >2.8.

Zoom In Zoom Out Reset image size

It is, however, important to emphasize at this point, that the broadband spectrum of ultrarelativistic electrons accelerated at the shock front, may be much more complex than a single power law. A single power-law model is used here rather for illustrative purposes, to provide a basic insight into the slope of the high-energy segment of the electron distribution, and the amount of spectral steepening observed across the shock front in the Pictor A hot spot.

The gradient in the HR values across the terminal reverse shock we have found has several important implications for understanding particle acceleration at relativistic shocks in general. First, the fact that the hardest X-ray spectra we see are concentrated at the upstream edge of the X-ray intensity peak, means that the efficient electron acceleration—forming flat electron energy distributions with indices s ≤ 2.2 (when approximated by a single power law) and electron energies of the order 10–100 TeV—takes place at the very front of the mildly relativistic shock with perpendicular magnetic field configuration, and not in the far downstream, for example, where compact radio knots have been found in VLBA observations (Tingay et al. 2008). Second, the HR gradient suggests that high-energy electrons advected from the shock front cool radiatively (leading to a steepening of their energy distribution and the corresponding X-ray spectrum). This confirms the origin of the offset between the X-ray hot spot and the VLA hot spot: the propagation length of the ultrarelativistic electrons that produce kiloelectronvolt photons is of the order of a parsec for the expected hot-spot magnetic field of 0.1–1 mG, and at most 100 parsecs for unrealistically low magnetic field intensity of a few microgauss (Thimmappa 2020), while it is much longer for radio-emitting electrons. By the time the jet has traveled ≃1 kpc between the X-ray hot spot and the compact radio knots, there are essentially no X-ray emitting electrons left.

Mildly relativistic magnetized shocks in electron–ion plasmas—meaning shock bulk Lorentz factors and magnetization parameters, defined as a ratio of the upstream Poynting flux to the kinetic energy flux, 10⁻³ < σ ≤ 0.1, matching the conditions expected to hold in the western hot spot of Pictor A—have been investigated with 2D kinetic particle-in-cell simulations by Sironi & Spitkovsky (2011), and more recently by Ligorini et al. (2021a, 2021b). These studies do show some energization of electrons, due to the resonant interactions with large-amplitude longitudinal Langmuir waves, combined with shock-surfing acceleration (Lyubarsky 2006; Hoshino 2008), however with a rather low efficiency when compared to ultrarelativistic shocks (i.e., shocks with γ_sh ≫ 10). As a consequence, the downstream electron spectra observed in such simulations (i) are basically thermal with little or no nonthermal power-law components, (ii) have total energy density much below that of the ions, at the level of about 10%, and (iii) have limited maximum energies E_e /m_e c² <(m_p /m_e )γ_sh < 10⁴. This is in contrast to the observational findings presented here, and elsewhere in the literature, regarding the ion-electron energy equipartition (see Stawarz et al. 2007), electron energies of the order of 10–100 TeV (see Sunada et al. 2022b), and flat slopes of the electron energy distribution (this work). Together, these observational findings indicate that electron acceleration is both fast and efficient at the jet termination shocks of luminous radio galaxies and quasars.

This research has made use of data obtained from the Chandra Data Archive. This work was supported by the Polish NSC grant 2016/22/E/ST9/00061 (R.T. and Ł.S.) and NASA award 80NSSC20K0645 (R.T., and J.N). The authors thank the anonymous referee and O. Kobzar for valuable comments and suggestions on the manuscript.

BL Lac object AO 0235+16 (z = 0.94), has been observed on 2000 20 August by Chandra on the ACIS-S3 chip (ObsID 884) with a 30.625 ks exposure time. The source spectrum was extracted from a circular region (position: R.A. = 2:38:38.9560, decl. = 16:36:59.440) with a radius of 3 px (≃15), and a 5–10 px annulus background. The background-subtracted source spectra were fitted within the soft (0.5–2.0 keV) and hard (2.5–7.0 keV) bands, assuming single power-law models modified by the Galactic column density N_H,Gal = 6.79 × 10²⁰ cm⁻² (HI4PI Collaboration et al. 2016).

Based on those spectra, we performed PSF simulations, and next produced deconvolved and HR images, all as presented in Figure 7. As shown, there is no substructure on the HR map: a point source at the position of the blazar, characterized by HR ≃ −0.5, is surrounded by the background with HR values of either −1 or 0.

Zoom In Zoom Out Reset image size

The radio-loud quasar B2251+134 (=4C+13.85, z = 0.673) has been observed (ObsID 2146) on 2000 January 18 with 25.8 ks exposure time. The source spectrum was extracted from a circular region (position: R.A. = 22:54:20.9771, decl. = 13:41:48.802) with a radius of 7 px (≃35), and the background annulus of 10-15 px radii.

We have modeled the Chandra spectra of B2251+134 in the soft (0.5–2.0 keV) and hard (2.5–7.0 keV) bands with single power-law models modified by the Galactic column density N_H,Gal = 4.67 × 10²⁰ cm⁻² (HI4PI Collaboration et al. 2016), this time however allowing for the intrinsic absorption in addition to the Galactic one. The intrinsic column density was kept as a free parameter when fitting the soft spectrum; the resulting best-fit value was then frozen when fitting the hard segment of the source spectrum. The results of the following image deconvolution are presented in Figure 8. Again, what we see is a well-defined point source in the center with HR ∼ −0.8, surrounded by the HR = −1 or 0 background.

Zoom In Zoom Out Reset image size