NuSTAR Hard X-Ray Monitoring of Gravitationally Lensed Quasar RX J1131–1231

The X-ray emission from active galactic nuclei is believed to come from a combination of inverse Compton scattering of photons from the accretion disk and reprocessing of the direct X-ray emission by reflection. We present hard (10–80 keV) and soft (0.5–8 keV) X-ray monitoring of a gravitationally lensed quasar RX J1131−1231 (hereafter RXJ1131) with NuSTAR, Swift, and XMM-Newton between 2016 June 10 and 2020 November 30. Comparing the amplitude of quasar microlensing variability at the hard and soft bands allows a size comparison, where larger sources lead to smaller microlensing variability. During the period between 2018 June 6 and 2020 November 30, where both the hard and soft light curves are available, the hard and soft bands varied by factors of 3.7 and 5.5, respectively, with rms variability of 0.40 ± 0.05 and 0.57 ± 0.02. Both the variability amplitude and rms are moderately smaller for the hard X-ray emission, indicating that the hard X-ray emission is moderately larger than the soft X-ray emission region. We found the reflection fraction from seven joint hard and soft X-ray monitoring epochs is effectively consistent with a constant with low significance variability. After decomposing the total X-ray flux into direct and reprocessed components, we find a smaller variability amplitude for the reprocessed flux compared to the direct emission. The power-law cutoff energy is constrained at 96−24+47 keV, which positions the system in the allowable parameter space due to the pair production limit.


INTRODUCTION
Quasars spectral energy distributions range from radio to X-ray bands (e.g., Elvis et al. 1994;Shang et al. 2011;Gupta et al. 2020), with jetted ones extending to the gamma-ray regime (e.g., Abdollahi et al. 2020;Yang et al. 2022).This diversity of features comes from the aggregated emission of many different physical components in the central engine (e.g., Antonucci 1993;Urry & Padovani 1995).In particular, the X-ray spectrum of quasars is composed of a direct X-ray continuum, a soft excess, and a Compton hump (e.g., Fabian et al. 1995;Nandra et al. 1997a;Fukazawa et al. 2011).The direct X-ray flux is generated by inverse Compton scattering, where ultra-violet photons from the accretion disk interact with relativistic electrons from the corona, causing the photons to gain energy and emit a power-law X-ray continuum (Haardt & Maraschi 1991, 1993).This X-ray continuum will irradiate the cold, optically thick material, and Compton scattering and photoelectric absorption will produce a reflection spectrum with associated metal fluorescent lines.The most prominent features of the reflection component are the FeKα line (6.4 keV) due to the high cosmic iron abundance and fluorescent yield and further enhanced by the fact that the reflection hump peaks at ∼ 20 keV (see Reynolds & Nowak (2003) for a review).
The spatial origin of the reflection region has been debated in the literature (e.g., Andonie et al. 2022;Masterson & Reynolds 2022), on whether the reflection occurs in a region immediately around the supermassive black hole, accretion disk, disk wind, or torus.The FeKα line profile provides a first clue to the reflection location.Since a narrow feature is commonly observed in the X-ray spectrum of AGN (e.g., Kaspi et al. 2001;Yaqoob et al. 2001;Zhou & Wang 2005), remote reflections, e.g., off the outer regions of the accretion disk or torus, should contribute to a portion of the reflection flux.For some nearby Compton-thick AGN, spatially resolved reflection component has been measured in a range of a few tens to several hundreds of parsecs around the central engine (e.g., Fabbiano et al. 2017;Kawamuro et al. 2019;Ma et al. 2020;Yi et al. 2021).In addition, a broad and skewed relativistic line component is observed in many high signal-to-noise ratio spectra of AGN (e.g., MCG-6-30-5, Tanaka et al. (1995);NGC 3516, Nandra et al. (1997b)), suggesting that the reflection in a region immediately around the supermassive black hole is present in many AGN (Fabian et al. 1989;Laor 1991).Light bending models were introduced (Miniutti & Fabian 2004), where the main X-ray emission is from a very compact source immediately above the black hole.Because of the strong gravitational light bending, the majority of the X-ray emission will return and irradiate the accretion disk and only a smaller fraction of the direct continuum will escape.If the source has vertical motion, it will result in a time varying partition between the direct and returned emission, and the fractional change of the direct component will be much larger than that of the reflected due to the different strengths of the light bending effect.These models require a very compact power-law source, which is confirmed by recent quasar microlensing (e.g., Morgan et al. 2008Morgan et al. , 2012;;Dai et al. 2010;Mosquera et al. 2013;Chartas et al. 2016;Guerras et al. 2017) and X-ray reverberation mapping methods (e.g.Reis & Miller (2013)).Separating the broad iron line from the underlying continuum may be complicated sometimes by the spectral curvature introduced by the partially ionized, and possibly non-uniform, absorbers (e.g., Reeves et al. 2018).It is also certainly possible that the reflection component is produced by a combination of heterogeneous sources.
A key to resolving this issue is to measure the relative size of the reflection component with respect to the (power law emitting) X-ray continuum to constrain the majority of the reflection region.Quasar microlensing has become a powerful tool to measure different components of AGN central engines that cannot be spatially resolved by current telescopes, with equivalent nano-arcsec resolutions (Kochanek et al. 2007;Moustakas et al. 2019).Quasar microlensing (Irwin et al. 1989;Vanderriest et al. 1989) is induced by stars in the lensing galaxy close to the gravitationally lensed images yielding extra magnifications or demagnifications to those images.As the source, lens, and observers move relatively, quasar microlensing will produce uncorrelated variability between the lensed images.Furthermore, over the relatively large energy band considered in our study, there is the possibility of two components that could in principle vary independently.An important scale in these studies is the Einstein ring size in the source plane for a typical microlens, where sources smaller than this scale are sensitive to microlensing effects, and for a 1M ⊙ star, the typical Einstein ring size is slightly smaller than the estimated broad line region sizes (e.g., Mosquera & Kochanek 2011).Thus, quasar microlensing has been used to constrain quasar central engine components from the broad line to the X-ray regions.Recent studies also show that planet mass objects in the lens galaxy will contribute significantly to the X-ray microlensing signal, because their Einstein ring scales match the compact emission regions immediately around the supermassive black hole (Dai & Guerras 2018;Bhatiani et al. 2019;Guerras et al. 2020).Since quasar microlensing has successfully constrained the relativistic Fe Kα region to be 2-8 r g in a sample of gravitationally lensed quasars (Dai et al. 2019), this implies that the reflection hump should have a similar size.Microlensing size measurements of the reflection region will provide an independent test to the current paradigm of relativistic reflection model, where the corona is compact with 10 r g and the reflection region is dominated by relativistic effects including light bending.
RX J1131−1231 (hereafter RXJ1131) is a quadruply lensed quasar system with the lens and source redshift at z l = 0.295 and z s = 0.658, respectively (Sluse et al. 2003), where the central black hole mass of the lensed quasar is estimated to be M BH = (1.3 ± 0.3) × 10 8 M ⊙ (Dai et al. 2010).The system shows large quasar microlensing variability based on 38 Chandra monitoring observations over a decade, where the image resolved flux ratios exhibit large variations (Chartas et al. 2017).It is also the brightest radio-quiet lensed quasar in the X-ray band (Jackson et al. 2015).These properties make RXJ1131 an ideal target to monitor with NuSTAR to constrain the relative size of the reflection hump compared to the direct X-ray continuum, and here we present the analysis of these NuSTAR monitoring observations paired with soft X-ray observations by Swift or XMM-Newton.The paper is organized as follows.We describe the observations and data analysis procedures in Section 2, spectral analysis in Section 3, and present the analysis results and discussion in Section 4. We assumed the set of cosmological parameters Ω m = 0.3, Ω Λ = 0.7, and H 0 = 70 km s −1 Mpc −1 in this paper.

X-RAY OBSERVATIONS AND DATA ANALYSIS
We have performed six pairs of joint observations with the Nuclear Spectroscopic Telescope Array (NuSTAR, Harrison et al. (2013)) and the Neil Gehrels Swift Observatory (Swift, (Burrows et al. 2005)) to monitor the hard (10-80 keV) and soft (0.5-8 keV) X-ray variability of RXJ1131 between 6 June 2018 and 30 November 2020.The Swift observations averaged an exposure time of 1.49 ks, and the NuSTAR observations 30 ks.We also analyzed a pair of archival joint observations from the X-ray Multi-Mirror Mission (XMM-Newton, Jansen et al. ( 2001)) with 97.7 ks exposure and NuSTAR with 77.4 ks exposure on 6 June 2018.We additionally analyzed 32 Swift observations between 20 June 2016 and 30 November 2020 besides the six paired with NuSTAR observations.The log of observations is listed in Table 1.We reduced the data using HEAsoft version 6.25 (Nasa High Energy Astrophysics Science Archive Research Center (Heasarc) 2014).The NuSTAR observations from Focal Plane Modules A and B (FPMA, FPMB) were processed with NUSTARDAS using CALDB version 1.0.2.We reprocessed the data using the NuSTAR data analysis pipeline nupipeline.We used the SCIENCE observation mode and did not filter for SAA after experimentation with different filtering methods showed consistent results with unfiltered data.The Stage 2 event files were used to further select source and background regions in accordance with the recommendations from the analysis guide.The circular source region was centered on the coordinates of RXJ1131 with a radius of 61.5".The background region had a radius of 196.8" and overlapped two quadrants.Default grade selection of grades 0-26 was used.Because of the slight differences in the FPMA and FPMB telescopes, the spectra were not co-added and instead were simultaneously fit in the spectral analysis.Ancillary response files were generated for each spectrum as part of Stage 3 processing with nupipeline through the numkarf subroutine.
The Swift data reduction was performed using XSelect.We used the response matrix file (RMF) for photon counting mode for default grades 0-12.Source and background region files were created.The circular source region was centered on the coordinates of RXJ1131 with a radius of 47.2", and the background region was placed significantly away from the source and had a radius of 295".The ancillary response file (ARF) was made for each observation using xrtmkarf.We binned the spectra by 15 counts to standardize the SNR for each bin.Some observations did not have enough counts to support this binning scheme; they were binned by 7 counts (see note in Table 1).We performed the spectral fitting in XSPEC1 .
The XMM-Newton EPIC data for observation 0820830101 were reduced using epchain in SAS (xmmsas-20201028-0905-19.0.0).Source and background photons with standard grades are extracted from circular regions of 50" radius centered on and away from the location of the source respectively.The response and area files were produced using the rmfgen and arfgen tools respectively.

Swift Spectra
We fit the 38 Swift spectra in the energy range 0.5 − 8 keV, ignoring the rest as recommended by the manual.We used χ 2 statistics for spectra binned by a minimum of 15 counts and C-statistics for those binned by 7 counts.
Because of the low SNR of these Swift spectra, with a median net count of 185, we can only fit these spectra with a simple powerlaw (zpowerlw) modified by Galactic (tbabs) and intrinsic absorption (ztbabs).We froze the Galactic neutral hydrogen column densities at 2.68 × 10 20 atoms/cm2 , which we obtained from the HEASARC online tool 2 .We found during the fitting process that the intrinsic absorption component is consistent with zero and therefore we dropped it in the subsequent analysis.This is consistent with the results from much higher SNR Chandra spectra (Chartas et al. 2017), which also constrained the intrinsic absorption to be zero.Several example spectra are shown in Figure 1.
Thus for each Swift spectrum, we are able to constrain the photon index and normalization.The measured photon indices range from 1.21 to 2.71 keV with a median of Γ = 1.69, which is consistent with the median of the Chandra spectra (Chartas et al. 2017).The uncertainties of the photon index are typically between 10 to 20 percent set by the typical SNR of these spectra.Finally, we set the photon index to the median value and refit the model to measure the unabsorbed model flux in the soft X-ray band from the 38 Swift spectra (Table 1).

Joint Soft and Hard X-ray Spectral Analysis
We analyzed the six Swift/NuSTAR and one XMM-Newton/NuSTAR spectral pairs jointly to constrain the reflection fraction of RXJ1131 over the soft and hard X-ray bands.The spectral analysis was performed in 0.5 -8.0, 10.0 -80.0, and 0.5 -8.0 keV bands for Swift, NuSTAR, and XMM-Newton, respectively.We first focused on the high SNR XMM-Newton/NuSTAR spectral pair to guide our spectral analysis process.We modeled the XMM-Newton/NuSTAR pair as an exponentially cut off power law reflected from neutral material (pexrav) or ionized material (pexriv) (Magdziarz & Zdziarski 1995) modified by Galactic absorption (tbabs).We binned the NuSTAR spectra by SNR of 1 which yielded the tightest constraints on model parameters based on our experiments of a range of binning schemes.The XMM spectrum was binned to have a minimum of 25 counts per bin with a minimum SNR of 6.We froze the model parameters for redshift, Galactic nH, metal and iron abundances (at solar abundances), and inclination (at 30 • , consistent with Seyfert 1 AGN).The values of the fit parameters were tied across all three spectra.The neutral reflection model also included a Gaussian centered at 0.7 keV to account for the soft excess component (Reis et al. 2014).The final model fit for photon index, reflection fraction, normalization, and cutoff energy is listed in Table 2, and the spectra and best-fit pexrav model were shown in Figure 2.Both the neutral and ionized reflection models yield acceptable fits to the joint XMM-Newton/NuSTAR spectrum in terms of reduced χ 2 , with the reflection fraction and cutoff energy constrained at r ref l = 0.22 +0.12 −0.10 , 0.34 +0.11 −0.10 and E cut = 96 +47 −24 , 107 +157 −27 , respectively, for the neutral and ionized reflection models.Both models have similar uncertainties on the reflection fraction; however, the cutoff energy is less constrained in the ionized model on the upper range and the ionization parameter ξ is poorly constrained for the ionized model.
Next, for the six lower SNR Swift/NuSTAR pairs, we used the characteristic photon index from the Swift observations and the cutoff energy from the XMM-Newton/NuSTAR joint fit to constrain the reflection fraction.We used the simpler neutral reflection model and fixed most of the parameters as with the XMM-Newton/NuSTAR paired observation, but only varying the reflection fraction and the normalization of each pair.As before, the parameters for all spectra in an epoch were tied together.Our best constraints were derived from binning NuSTAR spectra to a minimum SNR of 2 and Swift spectra a minimum of 15 counts per bin.Our results are listed in Table 3 and the temporal evolution of the reflection fraction is plotted in Figure 5. Model fits for the six Swift/NuSTAR pairs are plotted in Figure 3.We compare the hard (10.0-80.0keV) and soft (0.5-8.0 keV) band variability of RXJ1131.In the period where both hard and soft light curves are available between 5 June 2018 and 29 November 2020, where the hard band flux is measured by NuSTAR and soft by XMM-Newton and Swift, we measure a factor of 3.7 change in the hard band and 5.5 in the soft band, showing smaller hard X-ray variability.We next calculated the fractional root mean square (rms) amplitude F var (Vaughan et al. 2003) of the hard and soft band light curves, yielding 0.40 ± 0.05 and 0.57 ± 0.02 for the hard and soft bands respectively, where the soft band F var is larger by 3.2σ.
The total hard and soft X-ray flux is a combination of direct and reflected emission components.Our spectral analysis of 7 joint hard and soft X-ray spectra has constrained the average reflection scales as r ref l = 0.26, and the evolution of the reflection fraction is plotted in Figure 5 (left).We evaluate the variability of the reflection fraction over the seven epochs by fitting a constant model, yielding a residual χ 2 = 7.6 with 6 degrees of freedom, consistent with a constant model with 26% probability.For the last data point, there is an indication of an increase of the reflection fraction.The average reflection fraction for the first six observations is r ref l = 0.20, and the reflection fraction in the last observation, 0.65 +0.22  −0.2 , is 2.3σ away from this mean.We decompose the total X-ray flux into direct and reprocessed fluxes using the reflection fraction measurements for the seven epochs (Figure 5, right).The total, direct, and reprocessed emission curves change by factors of 4.1 ± 0.8, 10.0 ± 5.8, and 3.4 ± 3.2 respectively, by comparing the highest over lowest flux measurements.The large uncertainties for the direct and reprocessed components are caused by the measurement uncertainties of the reflection fraction.Examining the reprocessed emission curve shows that the excess variance is dominated by measurement errors, so the fractional variance is undefined, and a constant model fit to the reprocessed flux gives a residual χ 2 of 3.47, or a 75% probability for the constant model.However, if we discard the measurement uncertainties, the fractional variance of the flux estimates is 0.56 and 0.34 for the direct and reprocessed emission, respectively, suggesting a smaller variability amplitude for the reprocessed emission.
Assuming quasar microlensing is contributing significantly to the total variability, the smaller variability amplitude and fractional rms in the hard band measured by NuSTAR indicate a larger hard X-ray emitting region than soft X-ray emitting region.However, variability amplitude and fractional rms are only moderately smaller, qualitatively suggesting the hard X-ray emission region is only moderately larger than the soft size of 10r g measured by Dai et al. (2010).Although statistically insignificant, the reprocessed light curve suggests smaller variability amplitude as well.This is consistent with a picture that the reflection is occurring in a region in the accretion disk with a range larger than the soft X-ray size in RXJ1131.Remote reflection in the outer region of the accretion disk or torus as the majority of the reflection region in RXJ1131 can still be viable because of large measurement uncertainties.Our result can be consistent with remote reflection contributing to a fraction of the reprocessed flux, e.g.10-20% as suggested by other studies (Walton et al. 2014;Zoghbi et al. 2017).We defer the quantitative quasar microlensing analysis to a future paper with the on-going extended monitoring data.
We have constrained the cutoff energy as 96 +47 −24 keV in the neutral reflection model, consistent with the survey results in Ricci et al. (2018) which show characteristic cutoff energies for AGN of the order 10 2 keV.The ionized model yields similar cutoff energy values albeit with larger uncertainties in the upper bound.We place our measurement in context with other quasars in the compactness-temperature diagram.As photon-photon collisions occur in the corona at very high energies, the photons decay into electron-positron pairs and then further annihilate into photons.Therefore pair production can become a runaway process in a sufficiently hot and compact corona, and the mechanism can acts as a natural thermostat.This theory is demonstrated by many quasars sitting along the pair production balance line in the compactness versus temperature (l − θ) diagram, which is translated into observables in the L X − E cut diagram (Bertola et al. 2022).We show RXJ1131 in the L X − E cut plane compared to other quasars collected by Bertola et al. (2022) and references therein in Figure 6, where the X-ray luminosity has been corrected by a total magnification factor µ AGN = 48.2(Paraficz et al. 2018), and RXJ1131 lies on the allowable parameter space without violating the electron pair production limit.This lens model is based on the point-like (quasar) image position from the CO emission measured by ALMA, and the model is consistent with those based on HST or 2.2µm quasar positions.Since the cutoff energy has a relatively low energy of ∼100 keV, whether RXJ1131 lies in the allowable parameter space is almost independent of the magnification correction value.(Bertola et al. 2022).The green and blue curves represent the boundary line for runaway pair production presented as a polynomial in Ricci et al. (2018).We assume an X-ray source size of 10rg.Arrows in place of error bars represent a lower limit on the cutoff energy.

Figure 1 .
Figure 1.Example Swift spectra.Top left and top right: Typical SNR for the Swift spectra (19 and 18, respectively).Bottom left: An example of low SNR spectrum (4.8) excluded from analysis.Bottom right: An example of high SNR spectrum (24).

Figure 2 .
Figure 2. Top: Combined XMM-Newton (low energy) / NuSTAR (high energy) fit.The model is represented by the solid red (purple) line overlaying the XMM-Newton (NuSTAR) data.NuSTAR data is shown in two data sets coming from the two focal plane modules FPMA and FPMB.Bottom: Residuals of the fit.

Figure 3 .
Figure 3. Model fits for the six Swift/NuSTAR paired observations.Soft X-ray Swift spectra are in black, hard X-ray NuSTAR spectra are in gray, and the model is in red.

Figure 4 .
Figure 4. Unabsorbed flux light curve of all spectra from Nustar, Swift, XMM, and Chandra, where the Chandra flux is from Chartas et al. (2017); all other data was described in this paper.The dotted lines are not physical and only used for the purpose of better reading the individual data points.

Figure 5 .
Figure 5. Left: Reflection fraction curve.The first point is from the archival XMM-Newton-NuSTAR pair; the rest are paired Swift-NuSTAR observations.Right: Change in total, direct and reprocessed emission over time.

Figure 6 .
Figure 6.X-ray luminosity versus cutoff energy from a selection of quasars(Bertola et al. 2022).The green and blue curves represent the boundary line for runaway pair production presented as a polynomial inRicci et al. (2018).We assume an X-ray source size of 10rg.Arrows in place of error bars represent a lower limit on the cutoff energy.

Table 1 .
Observation Log a Archival observation reanalyzed in this paper.

Table 3 .
Reflection fraction Fit Results