X-ray AGN Emission of J033226#
The following example shows how we can use Lightning to fit the X-ray to IR SED of a bright AGN, J033226.49-274035.5, in the Chandra Deep Field South.
Data#
Typically it will only be necessary to fit with 2-4 band X-ray photometry (or even one band, e.g. if we
assume negligible intrinsic absorption). Here, to show what’s possible in principle with the Lightning X-ray
implementation, we’re going to fit 15-band X-ray photometry. To do so, we downloaded the level 3 spectral file
and response from the Chandra Source Catalog for the single deepest (\(\approx\) 163 ks) observation of
the source (ObsID 5019). We then subtracted the background and grouped the net counts from the spectrum into
15 log-spaced bins from \(0.5-6.0~\rm keV\) using Sherpa. [2] We calculated the uncertainty
\(\sigma N\) on the net counts \(N\) as
where
and
for source and background counts \(S\) and \(B\), and \({\tt BKG\_SCALE}\) is the background exposure scale. [3]
We retrieved optical-IR photometry from the Guo et al. (2013) CANDELS photometric catalogs, including IR photometry from the Barro et al. (2019) catalogs. We added calibration uncertainties (typically \(<5\%\)) to each band as specified in Table 1 of Doore et al. (2022).
The relevant data files can be found in examples/AGN_J033226/ as xray/J033226_photometry.fits,
xray/J033226_xray_photometry.fits, xray/J033226_summed.arf for the UV-to-IR photometry, X-ray data,
and ARF, respectively.
Configuration#
For this example we’ll configure Lightning to use the SKIRTOR AGN Emission (UV-to-IR) templates and fit both with and without an AGN X-ray Emission model to see how adding X-ray data to the fits affects the results.
Common Settings#
For both runs, we change the below lines of the configuration files in examples/AGN_J033226/noxray/
and examples/AGN_J033226/xray/ from the defaults.
First we change the output filename:
71OUTPUT_FILENAME: 'J033226_noxray_output' ,$
We add \(10\%\) model uncertainty:
91MODEL_UNC: 0.10 ,$
The default initializations on the SFH are then widened, as we expect this \(z \sim 1\) AGN host to have a large SFR:
189PSI: {Prior: ['uniform', 'uniform', 'uniform', 'uniform', 'uniform'] ,$
190 Prior_arg: [[0.0d, 0.0d, 0.0d, 0.0d, 0.0d] ,$
191 [1.d3, 1.d3, 1.d3, 1.d3, 1.d3]] ,$
192 Initialization_range: [[0.0d, 0.0d, 0.0d, 0.0d, 0.0d] ,$
193 [1.d2, 1.d2, 1.d2, 1.d2, 1.d2]] $
194 } ,$
We then limit the mass abundance of PAHs in the dust model
385QPAH: {Prior: 'fixed' ,$
386 Prior_arg: [4.7d-3] ,$
387 Initialization_range: [4.7d-3, 4.58d-2] $
388 } ,$
And turn on the UV-IR AGN model:
486AGN_MODEL: 'SKIRTOR' ,$
Since J033226 is classified as a type 1 AGN from optical spectra, we limit the inclination of the models to type 1 AGN accordingly:
516AGN_COSI: {Prior: 'uniform' ,$
517 Prior_arg: [0.75d, 1.0d] ,$
518 Initialization_range: [0.75d, 1.0d] $
519 } ,$
We also set the prior on the AGN normalization to log-uniform on \(10^{11}-10^{13}\ \rm L_{\odot}\) and initialize in a narrow range of luminosities, as initializing the AGN component far from the solution has a negative effect on convergence.
LOG_L_AGN: {Prior: 'uniform' ,$
Prior_arg: [11.d0, 13.d0] ,$
Initialization_range: [12.d, 12.5d] $
} ,$
Finally, we increase the number of MCMC trials per walker to \(4 \times 10^4\):
539NTRIALS: 4e4 ,$
and change the proposal width parameter to \(a = 1.8\) to increase the MCMC acceptance rate:
572AFFINE_A: 1.8d ,$
X-ray Fit#
For the fit with the X-ray model, we additionally modify the following settings in
examples/AGN_J033226/xray/lightning_configure.pro.
We turn on the X-ray emission module,
415XRAY_EMISSION: 1 ,$
set the module to use the errors we calculated on the net counts,
427XRAY_UNC: 'USER' ,$
modify the default initialization range for the column density,
439NH: {Prior: 'uniform' ,$
440 Prior_arg: [1.d-4, 1.d5] ,$
441 Initialization_range: [1.d0, 1.d2] $
442 } ,$
and narrow the initialization range for the SMBH mass:
459AGN_MASS: {Prior: 'uniform' ,$
460 Prior_arg: [1.d5, 1.d10] ,$
461 Initialization_range: [1.d6, 1.d8] $
462 } ,$
Finally, we generate high resolution models for all final MCMC chain elements so that we can use them to calculate other properties later.
654HIGH_RES_MODEL_FRACTION: 1.0 ,$
Running Lightning#
Note
The IDL code snippets below are also available in batch file format, as examples/AGN_J033226/J033226_batch.pro.
At this point, we are ready to run Lightning with each configuration. We can do this with:
cd, !lightning_dir + 'examples/AGN_J033226/'
restore, !lightning_dir + 'lightning.sav'
lightning, 'noxray/J033226_lightning_input.fits'
lightning, 'xray/J033226_lightning_input.fits'
With the MCMC configuration we’ve selected for the fit with no X-ray model (since we’re aiming for a comprehensive sampling of the posterior), this may take around 40-50 minutes on a moderately powerful laptop CPU (we ran it on a ca. 2015 2.9 GHz Intel Core i5). For the fit with the X-ray data, this will take about an hour, due to the added complexity of the X-ray model and the additional photometry. We note that the fits could also be run simultaneously in separate IDL sessions.
Analysis#
Once the fits finish, Lightning will automatically create post-processed files for us containing the models and the sampled posterior distributions.
noxray_res = mrdfits('noxray/lightning_output/J033226_noxray_output.fits.gz', 1)
xray_res = mrdfits('xray/lightning_output/J033226_xray_output.fits.gz', 1)
Convergence#
With the results loaded, our first step should be to check on the status of our fits:
print, '//Convergence for the fit without the X-ray model//'
print, 'Mean acceptance fraction: ', strtrim(mean(noxray_res.ACCEPTANCE_FRAC), 2)
print, 'Convergence flag: ', strtrim(noxray_res.CONVERGENCE_FLAG, 2)
print, 'Short chain flag: ', strtrim(noxray_res.SHORT_CHAIN_FLAG, 2)
print, 'Number of "stranded" walkers: ', strtrim(total(noxray_res.STRANDED_FLAG), 2)
which outputs
//Convergence for the fit without the X-ray model//
Mean acceptance fraction: 0.25888400
Convergence flag: 0
Short chain flag: 0
Number of "stranded" walkers: 2.00000
It seems that we have a reasonable acceptance fraction (i.e., \(>20 \%\)), and that we can be confident in the solution and the number of samples in the posterior, since the convergence flag is 0. We see however that 2/75 of the walkers in the ensemble didn’t reach the same solution as the others; this tends to happen when they’re initialized too close to the boundary set by the priors. We could adjust the settings and retry, but we have enough samples in the posterior as-is, so we’ll leave it.
Now we check the X-ray fit:
print, '//Convergence for the fit with the X-ray model//'
print, 'Mean acceptance fraction: ', strtrim(mean(xray_res.ACCEPTANCE_FRAC), 2)
print, 'Convergence flag: ', strtrim(xray_res.CONVERGENCE_FLAG, 2)
print, 'Short chain flag: ', strtrim(xray_res.SHORT_CHAIN_FLAG, 2)
print, 'Number of "stranded" walkers: ', strtrim(total(xray_res.STRANDED_FLAG), 2)
//Convergence for the fit with the X-ray model//
Mean acceptance fraction: 0.22930733
Convergence flag: 1
Short chain flag: 0
Number of "stranded" walkers: 2.00000
We can see that Lightning is concerned about whether the MCMC chains converged. If we do
print, 'Number of walkers with low acceptance fractions: ', strtrim(total(xray_res.ACCEPTANCE_FLAG), 2)
print, "Number of parameters with high autocorrelation times: ", strtrim(total(xray_res.AUTOCORR_FLAG), 2)
we see
Number of walkers with low acceptance fractions: 0.00000
Number of parameters with high autocorrelation times: 3.00000
The convergence flag is set because a three parameters had high autocorrelation times.
Let’s check what parameter this is and how it compares to our factor of 50 threshold for
considering convergence (i.e., NTRIALS/AUTOCORR_TIME should be ~50).
print, '//Autocorrelation Time Checks//'
print, 'Parameter with high autocorrelation time: ', $
strtrim((xray_res.parameter_names)[where(xray_res.autocorr_flag eq 1, /null)], 2)
print, 'Ratio of Ntrials to autocorr_time: ', $
strtrim(4e4/((xray_res.autocorr_time)[where(xray_res.autocorr_flag eq 1, /null)]), 2)
//Autocorrelation Time Checks//
Parameter with high autocorrelation time: TAUV AGN_LOGMDOT AGN_COSI
Ratio of Ntrials to autocorr_time: 44.566111 47.838389 38.319755
So as we can see, the “problem” parameters are the optical depth, accretion rate, and AGN inclination. Additionally, we can see that the optical depth and accretion rate are just barely below our default threshold of 50, while the inclination is well below. Since we want a value of ~50, this is close enough for the optical depth and accretion rate. However, for the inclination, this is likely due to the strong prior we introduced, which is limiting the model to an inclination range it does not prefer. Therefore, we will still use this solution and note that our inclination estimates may not be the most independently sampled.
Figures#
Now, we want to look at the SED fit. We’ll plot the UV-IR SED component in the “traditional” way,
but since we’ve got so much X-ray photometry we’ll plot the X-ray component in the way you might expect
from XSpec or Sherpa, just for fun. We have prepared a function for this: J033226_spectrum_plots.pro.
p = J033226_spectrum_plots(xray_res)
The plot of the best-fitting models looks like so:
We can see that the joint X-ray-IR fit works quite well in this case. Remember that though we’ve plotted the X-ray and optical-IR components on separate axes, they are fit jointly: the likelihood in the MCMC is the combined likelihood of both components.
Now, let’s compare the posteriors. However, since the two fits have different parameters, we’ll do some additional
postprocessing to generate the posterior on the UV-IR AGN luminosity for the fit with the X-ray model. The
script calc_integrated_AGN_luminosity.pro in the lightning/examples/AGN_J033226/ directory
generates this posterior.
restore, 'xray/lightning_output/lightning_configure.sav'
; The ARF is not strictly necessary, but LIGHTNING_MODELS expects it
arf = mrdfits('xray/J033226_summed.arf', 1)
AGN_model_lum = calc_integrated_AGN_luminosity(xray_res, config, arf)
help, AGN_model_lum
which shows the following:
** Structure <62b85a08>, 9 tags, length=56024, data length=56024, refs=1:
SED_ID STRING 'J033226.49-274035.5'
REDSHIFT DOUBLE 1.0310000
NH DOUBLE Array[1000]
AGN_MASS DOUBLE Array[1000]
AGN_LOGMDOT DOUBLE Array[1000]
TAU97 DOUBLE Array[1000]
AGN_COSI DOUBLE Array[1000]
L2500 DOUBLE Array[1000]
LBOL_AGN_MODEL DOUBLE Array[1000]
The agn_model_lum.LBOL_AGN_MODEL is the posterior on the integrated optical-IR luminosity for the X-ray fit.
With this in hand we can compare the posterior distributions of the AGN model parameters. Again, we have a plotting
function for this, in posterior_comparison.pro:
p = posterior_comparison(xray_res, noxray_res, AGN_model_lum)
Which generates the following plot:
Footnotes