title: Relativistic accretion disk reflection in AGN X-ray spectra at z=0.5--4: a study of four \textit{Chandra} deep fields
authors: Linda Baronchelli Kirpal Nandra Johannes Buchner
text: 
  
  firstpage–lastpage
  
   The polar-generalized normal distribution
      
   Masoud  Faridi,
  Majid Jafari Khaledi^a
  
  masoudfaridi@modares.ac.ir
  
  ^a Department of Statistics, Tarbiat Modares University, P.O. Box 14115-134, Tehran, Iran
  
  jafari-m@modares.ac.ir
  
      
  ============================================================================================================================================================================================
  
  
  
  
  
  We confirm that the spectra are best fit by a model containing two Compton reflection components, one from distant material, and the other displaying relativistic broadening, most likely from the inner accretion disk. The degree of relativistic broadening indicates a preference for high black hole spin, but the reflection is weaker than that expected for a flat disk illuminated by a point source. 
  We investigate the Compton reflection signatures as a function of luminosity, redshift and obscuration, confirming an X-ray Baldwin effect for both the narrow and broad components of the iron line. 
  Anti-correlations are also seen with redshift and obscuring column density, but are difficult to disentangle from the Baldwin effect. 
  Our methodology is able to extract information from multiple spectra with low signal-to-noise ratio, and can be applied to future data sets such as eROSITA. We show using simulations, however, that it is necessary to apply an appropriate signal-to-noise ratio cut to the samples to ensure the spectra add useful information.  
  
  
  
  
  
  galaxies: active – galaxies: Seyfert – galaxies: high-redshift – X-rays: galaxies
  
  
  
  
  
  
  
  
  § INTRODUCTION
  
  Active galactic nuclei (AGN) are powered by matter falling onto a supermassive black hole (SMBH) <cit.>, 
  via an accretion disk <cit.>. 
  
  In addition to the direct disk emission AGN emit up to 20% of their bolometric luminosity in the X-ray waveband <cit.>. The X-ray emission is produced via Compton upscattering by a hot corona <cit.> which illuminates the innermost regions of the accretion disk. The X-ray emission can be backscattered and induce fluorescence in the inner disk <cit.>, which leaves imprints due to the large velocities and gravitational field. Hence, the analysis of the X-ray spectrum from AGN allows us to probe general relativity and investigate the behaviour of matter in extreme gravitational fields, only a few gravitational radii from SMBHs <cit.>.
  
  In many of these systems a strong iron (Fe) Kα fluorescent line is observed <cit.>. The Fe Kα feature peaks around 6.4 keV in the rest frame. Part of the emission comes from material at scales of several parsec, most likely the torus envisaged in orientation-dependent unification schemes <cit.> and hence the line is relatively narrow, with velocities of a few ∼ 100 km s^-1 <cit.>.
  As discussed above, the remainder of the iron Kα feature is emitted in a region of the accretion disk in the proximity of the SMBH, and it is broadened and skewed by relativistic effects, e.g. gravitational redshift and relativistic Doppler shifts <cit.>.
  The shape and width of the line and the amount of broadening can help us gain information on the geometry of the system and the spin of the SMBH <cit.>. Constraining the distribution of black hole spins in the whole AGN population would help determine the nature of the accretion history of the SMBH, for example by allowing us to distinguish whether the SMBH grew mostly through mergers, continuous, or "chaotic" accretion <cit.>.
  
  Strong, relativistically broadened Fe Kα lines are observed in a number of individual, bright Seyfert galaxies in the nearby Universe, e.g. MCG-6-30-15 <cit.> and NGC 3516 <cit.>. Indeed such emission is found to be common in the nearby AGN population <cit.>, though the evidence for such a component is not universal <cit.>. 
  This is probably due to fact that the strength of 
  the broad reflection component is weaker than that expected from a standard accretion disc illuminated by a point source <cit.>, meaning that very high signal-to-noise ratio (S/N) is needed to detect it convincingly <cit.>.
  
  The fact that the reflection is weaker than expected, and differs in strength from source to source, can nevertheless reveal important information about the system such as the geometry. 
  A well known effect of this type is the X-ray Baldwin or Iwasawa-Taniguchi effect <cit.> where the equivalent width of the iron line decreases with luminosity. This effect has been observed for the narrow <cit.> 
  and broad <cit.> components of the line. 
  
  Clearly it is then important to characterize the iron line emission of a representative sample of AGN, including more typical objects at higher redshift. There are a few observations of relativistic broadened Fe Kα lines at high-redshift, sometimes thanks to studies of lensed Quasars <cit.>. The ubiquity or otherwise of the broad features in samples of typical AGN beyond the local Universe is hard to establish, however, given the faintness of the targets. The deepest X-ray surveys offer the opportunity to investigate this issue, and several attempts have been made to use deep surveys to infer the properties of a population of AGN and to verify the ubiquity of feature like the broadened Fe Kα line in a population of AGN  <cit.>.
  
  Most past studies have relied on stacking of large samples of low count X-ray spectra <cit.>. By stacking, one may be able to infer the population properties in cases where fitting individual spectra would not yield meaningful results. Evidence for broadening of the Fe Kα line has been reported in several stacking studies, such as  <cit.> and <cit.>. <cit.>, however, combined the rest-frame spectra of 600 XMM-Newton type-1 AGN without finding compelling evidence for a relativistic broadening of the Fe Kα line.
  Some of this work showed, however, that stacking the spectra might induce artificial broadening of the Fe Kα in samples with a wide redshift distribution <cit.>. This problem was addressed by <cit.> and <cit.> who compared the stacked spectra to simulations. While this approach can increase confidence in the existence of the broadened features, the simulations require an assumption about the true underlying spectrum. 
  
  In our previous work <cit.> we employed an alternative technique to characterise the and the iron line and reflection properties of a large sample of AGN, whose X-ray spectra individually have low S/N ratio.
  Instead of fitting a stacked spectrum, we used the Bayesian X-ray Analysis (BXA) software <cit.> to fit the individual spectra of low S/N sources. We then combined the Bayesian evidence for the putative broad reflection component to establish whether or not it was present in the whole sample. 
  The analysis was performed on a sample of 199 hard X-ray selected sources from the Chandra Deep Field South 4Ms exposure, and revealed strong evidence for a relativistically broadened X-ray reflection component from the accretion disk. The properties of the reflection also implied a preference for a maximally spinning SMBH, as compared to a non-rotating Schwarzschild black hole.
  
  In this work, we expand the study of <cit.> to a total of four Chandra fields.
  Our aim is to investigate further the prevalence of the Fe Kα and reflection features within typical AGN up to z=4, and characterize their properties. With our expanded sample, we aim to confirm our previous results, and investigate the dependence of the reflection strength with other parameters such as the luminosity, redshift and obscuration. 
  
  The paper is structured as follows: in Sect. <ref>, we describe the data used in this work, and our methods of spectral and statistical analysis. The results of the work are reported in Sect. <ref> and interpreted and discussed in Sect. <ref>.  Sect. <ref> summarises our results and presents our primary conclusions. 
  
  Throughout this work, we adopt Ω_m= 0.272, Ω_Λ= 0.728, and H_0 = 70.4 km s^-1 Mpc^-1 <cit.>.
  
  
  
  
  § SAMPLE AND METHOD
  
  
  
  
  
  
  
  
  
   §.§ Data
  
  
  In this work we analyse four of the deepest fields observed by Chandra, the 7Ms exposure of the Chandra Deep Field South (CDFS), the Chandra Deep Field North (CDFN), AEGIS and the COSMOS fields. The long exposures in these fields ensure that meaningful spectral information can be extracted for typical AGN at moderate redshifts (z=0-4) which dominate the accretion history of the Universe <cit.>. In our analysis, we focus on the properties of the iron Kα emission line and hard X-ray reflection continuum from Compton thick structures surrounding the AGN, such as the accretion disk and molecular torus. We hence use a hard X-ray selection (>2 keV) for all fields to ensure that meaningful spectral constraints can be obtained using these features. Here we summarise the data and sample selection in each field: 
  
   
  
    §.§.§ The Chandra Deep Field South
  
  
  The CDFS <cit.> is, with a nominal total exposure time of ∼7 Ms, the deepest of all the Chandra surveys, and indeed the deepest X-ray survey of all. While it reaches extremely faint fluxes, it covers a relatively small area of ∼0.13 deg^2 
  The CDF-S 7Ms is a collection of observations performed over multiple epochs between Oct 14, 1999 and Mar 24, 2016. 
  
  All 102 observation used ACIS-I, which offers spectral imaging over an approximataely 17 × 17 arcmin  field of view 
  and is often used for surveys. 
  
  We limit our CDFS sample to the 199 hard X-ray (2-7 keV) selected AGN at redshift z<4 previously studied in <cit.> and <cit.>. These were selected form the source catalog of <cit.>, which was based on the 4Ms Chandra exposure. In the current analysis, however, we extract the spectra of the 4Ms sources from the deeper 7Ms exposure. 
  The 7Ms spectra of these 199 sources contain a total of 313914 source counts in the observed 1–8 keV energy band, 
  which is most relevant for our analysis. All of these sources have a redshift measurement, the majority of which are spectroscopic redshifts for which we adopt a single value. For the remaining 38% of these sources (76/199) photometric redshifts and their probability distributions from <cit.> were used. These are specially tailored for AGN, and are based on the methods of <cit.> and <cit.>.
  
  
  
    §.§.§ The Chandra Deep Field North
  
  
  The CDFN is a field which is has received a Chandra exposure of 2 Ms over a sky area of ∼0.12 deg^2
  
  It comprises 20 different pointings taken between November 1999 and February 2002. We use the source catalog from <cit.> to obtain the redshift values and to exclude stars. 
  We analyse a sub-sample of 376 sources with redshift information from the 411 hard selected (2–7 keV band) sources in <cit.>. Of these sources, 159 have spectroscopic redshift <cit.>. 
  The spectra contain a total of 
  127563 source counts in the 1–8 keV energy band. The PDFs of the photometric redshift for the CDFN are not available, thus we used the preferred redshift adopted in <cit.>. The fact that the photo-z PDFs are not available for the CDF-N means that, in some cases, an inaccurate redshift will be adopted in the spectral fit. 
  On average, this will have 
  the effect of reducing the significance when comparing the true underlying model with any other model, so should be conservative with respect to the significance
  of the results presented below. 
  Because a single value for the redshift is used, it will also result in an overly narrow posterior distributions for those sources.
  This would also tend to lead to an overestimate of the intrinsic scatter and an 
  underestimate of the statistical uncertainty of the mean values derived for the sample.
  
  
  
  
    §.§.§ AEGIS-X
  
  
  The AEGIS-X Deep survey <cit.> is the result of deep Chandra imaging of the central region of the Extended Groth Strip. The survey encompasses an area of approximately 0.29 deg^2 with a nominal exposure time of 800 ks. AEGIS-X is currently the third deepest Chandra blank field survey after the Chandra Deep Fields (CDF). While being shallower than the CDFs by a factor of ∼2–3 it covers an area ∼3 times larger. We use the source catalog from <cit.>, selecting as a parent sample the sources detected in the 2-7 keV band, comprising 572 sources. After removing the sources identified as stars in <cit.>, we select a sub-sample of 540 sources with redshift information from the original 572 sources in the hard selected sample, with a total of 119462 source counts in the 1–8 keV energy range. Of these 540 sources, 
  202 (∼ 37%)  have a spectroscopic redshift. For the remainder, <cit.> provide photometric redshifts tailored for AGN, and their probability distribution functions, which we use in the spectral fitting.
  
  
  
    §.§.§ COSMOS
  
  
  The Chandra COSMOS Legacy survey spans an area of 2.2 deg^2 on the sky. The central 1.5 deg^2 has a nominal exposure of ∼ 160 ks while the surrounding regions are nominally exposed with ∼ 80 ks depth <cit.>. This makes it the largest area survey in our complication, but also the shallowest. We study the hard band selected sample from <cit.> after removing the sources identified as stars in <cit.>, which comprises 1122 objects with a total of 95014 source counts in the 1-8 keV energy range. Of these sources, 534 have spectroscopic redshift, while the remaining 53% have accurate photometric redshifts and photo-z probability distributions from the work of <cit.>.
  
  
  
   §.§ Combined sample
  
  
  
  
  
  The parent sample for this study, combining all four fields, comprises a total of 2237 sources (See Table <ref>). 
  Figure <ref> shows the distribution of counts and S/N of our sample. 
  We calculate the S/N using the formalism of <cit.> (see also <cit.>), which takes into account the Poisson nature of both the source and background count measurements. In particular we use the function  from the python library  developed by <cit.>. The sample spans a wide range of both total counts and signal-to-noise ratio. Numerically it is dominated by sources from the COSMOS survey, which has the largest area, but the objects from the deepest field, the CDF-S, have the highest number of counts and S/N overall. The right hand panel of  Fig. <ref> shows a vertical line at S/N=7. We use sub-samples cut at this S/N in the subsequent analysis, as discussed below, and the number of sources and source counts obtained after applying this S/N cut are shown in Table <ref>. 
  
  The luminosities L, defined in the 2–10 keV energy range, and redshifts of the selected sources are presented in Figure <ref>, which also shows the S/N split sub-samples. The luminosities were calculated from the spectral fits described below and corrected for galactic and intrinsic absorption. 
  The sample covers a broad redshift range up to z∼ 4, and the bulk of the sample covers the luminosity range log(L/erg s^-1)=42-45[The log here represents the logarithm in base 10 (log_10).], with just a very few low redshift sources fainter than this. The sample becomes increasingly incomplete at luminosities below log(L/erg s^-1)=43 above a redshift of about 1.   
  
  
  
   §.§ Spectral Extraction
  
  
  
  We extracted the spectra of the sources using the software package ACIS Extract (AE) <cit.>.
  AE was developed to automate as much as possible the analysis of X-ray data taken with the ACIS instrument of Chandra. It is well suited for our application as, given a source catalog, spectra can be extracted from multiple observations of the same field. As input files, AE requires the Level 2 event list of the observations, the exposure maps and aspect histograms corresponding to the field of view of the event data, the aspect solution file covering the time range of the observations and the mask file of the observations. We produced these files with a combination of the ciao 4.10 <cit.> and ftools 6.25 <cit.> software. 
  
  The process of extracting spectra using AE can be summarized in four steps. First, the AE tool ae_make_catalog builds extraction regions sized to encompass 90% of the local point spread function (PSF) but small enough to avoid overlaps in crowded regions.
  Secondly, the tool ae_standard_extraction extracts source and background spectra of the sources in the catalog. Note that for data taken at -110C on certain CCDs the event file will not be corrected for charge transfer inefficiency (CTI) and AE will be set to use the ciao tool mkrmf instead of the default mkacisrmf to build the RMF files. The next step uses ae_adjust_backscal_range to analyze the source's existing background extraction region and to choose a target background scaling range individually for each source. The process of extracting background and choosing a scaling range has to be repeated until the scaling range is stable. Finally, AE merges the observations combining the extraction from all the ObsIds and performing the photometry.
  
  We produced the appropriate input files for AE following the method presented in <cit.>. The final data products are the source and background spectra together with RMF and ARF files.
  
  
  
  
  
  
  
  
   §.§ Model fitting and model comparison
  
  
  
  Four physically justified models are considered to represent different scenarios for the major gas structures surrounding the central black hole. The first model is an absorbed power-law,  in Xspec terminology, which describes the emission from an X-ray corona behind a screen of obscuring gas. In the second model we add a non-relativistic reflection component to the simple absorbed power-law, , to represent reflection from distant material such as the obscuring torus. We chose the  model <cit.> to describe the reflection component since it combines a exponentially cutoff power-law emission reflected by neutral material (; <cit.>) with self-consistently generated Fe Kα, Fe Kβ, Ni Kα and Fe Kα Compton shoulder emission <cit.>. The third and the fourth model add a relativistically broadened reflection component to the second model, as expected from an accretion disk. This last component is modelled by convolving a narrow reflection spectrum with a convolution model to represent the expected Doppler and gravitational shifts expected from an accretion disk. Specifically, we use the the  model <cit.>. The  allows the BH spin to be a free parameter. However, to avoid having unnecessarily many free parameters in the broadened reflection model we constrain our analysis by fixing the spin parameter to the two special cases of spin a=0 (Schwarzschild metric, model three) and maximally spinning a=0.998 (model four). The parameter priors are chosen to be consistent with <cit.>, except for the parameter describing the inclination angle of the broad component, that in the current work is chosen to be uniformly distributed in cosine space (see Table <ref>). In Table <ref> we list the chosen model parameters for the components , ,  and . For the parameters allowed to vary, we list the minimum and maximum value of the prior distributions. The normalization parameters and the N_H are chosen to be uniform in logarithmic space, the inclination angle in  is uniform in cosine space, while all the other parameters are uniform in linear space. When only the photometric redshift of the source is available, we use the probability distribution function (PDF) produced using the SED fitting procedures and templates from <cit.> and <cit.> to take into account of the uncertainty of the photometric redshift estimation. The exception is the CDFN for which the photometric redshift PDFs are not available, and we use a single value for the redshift.
  
  We fit a background model simultaneously to the data, following <cit.>, instead of subtracting the background. We fit the Chandra background model provided in BXA to all the background spectra and then we include the background model with best fit parameters frozen in the model of the source spectra <cit.>.
  As discussed in the introduction, we do not stack the spectra to determine the average properties of the sample, but instead fit each source individually and then infer the properties of the sample by combining the information from these individual fits. 
  
  
  
  
  
  
  
  
   §.§ Bayesian X-ray Analysis
  
  
  The Bayesian X-ray Analysis (BXA) <cit.> package is a Bayesian framework to determine the best fit parameters and their posterior distribution for X-ray spectra. 
  BXA applies  <cit.>, a python wrapping of an implementation of the nested sampling algorithm (Multinest) <cit.> combined with Sherpa or Xspec to compute the  Bayesian evidence Z for X-ray data and hence parameter constraints. The Bayesian evidence Z is the integral of the likelihood over the prior and can be interpreted as the probability P(D|M) of the model M given the data D marginalized over the model parameters θ.
  
  
      Z=P(D|M)=∫ L(θ)P(θ|M)dθ,
  
  MultiNest <cit.> provides an efficient approximation to this integral. This algorithm samples a number of live points in the parameter prior space evaluating the likelihood of the model for every point. At every step, the live point with lowest likelihood will be replaced by a newly sampled point until the algorithm converges to the highest likelihood value. MultiNest is particularly specialized to deal efficiently with multi-modal distributions by using a recursive clustering algorithm and proposal regions in the shape of ellipsoids.
  
  To perform model comparison, we assume that all sources are described by the same model. We calculate the total evidence for that model by adding the log(Z) <cit.> values from the individual fits. This allows one to compare the total evidence of the sample for different models. The difference in the logarithmic evidence then corresponds to a Bayes factor (BF), which can be used to discriminate between the models. A commonly used way to interpret the BF values is the Jeffrey scale, which strengthens the choice of one model over the other approximately every time that the logarithm of the BF increases by one in natural logarithmic units <cit.>.  However, Bayes factors are continuous quantities and such discretisations should not be over-interpreted <cit.>. 
  In <cit.>, we verified with simulations that the Bayes factors scatter around one (Δlog Z≈ 0) for low signal-to-noise data.
  
  
  
  § RESULTS
  
  
  
  
  
  
  
  
  
  
   §.§ Initial fitting and S/N effects
  
  
  
  We first fit the total sample of 2237 sources with BXA to calculate the Bayesian evidence Z for the four models described in Section <ref>. We fit each source with each of the models individually and use the Bayes factor (BF) method to compare the models and determine which one provides the better fit. The single sources are generally too faint to significantly favour one model over the others. Instead, as discussed above, we combine the Bayesian evidence from the individual source fits to obtain the evidence for the full sample <cit.>. 
  
  As can be seen in Table <ref>, the Bayesian evidence for the total sample shows a strong preference for the model including both narrow and broad reflection components, with a maximally spinning black hole preferred over a Schwarzschild solution. This preference is also found in the field-by-field subsamples, with the exception of the COSMOS field. Here the preferred model is that without the blurred reflection component. We discuss this result further below, but note here that there are a number of differences between the COSMOS sample and the remainder of the fields. One such difference is that, as shown in Fig <ref>, the COSMOS sources typically have lower signal-to-noise ratio than the sources in the other fields. This raises the possibility that, below a certain S/N threshold, the spectra become insensitive to the broad reflection component and do not add information about the presence or properties of that component. 
  
  This hypothesis is borne out by the data. Fig. <ref> (top) shows the cumulative evidence as a function of signal-to-noise ratio for the sample. The evidence rises rapidly when adding sources with high S/N, but then flattens off and eventually becomes approximately horizontal, showing that the lowest quality spectra are not adding additional information/evidence. To estimate at which S/N the source will not add further information to the total evidence of the sample, we perform a set of simulations using the  tool of Xspec <cit.>.
  
  We simulate a sample of 300 sources using the ancillary files (ARF and RMF), the background spectra and the redshifts from the AEGIS sample. As a first step to simulate a sample of spectra, we have to choose an input model that will define the spectral shape of the simulations. Since we are interested in studying how our method would perform on a sample of relativistic broadened spectra at different S/N, we chose to simulate the spectra using the fourth model.  We simulate the sources so that the total number of net counts is comparable with the one of CDFS 7Ms sample.
  
  We fit the simulated sample with the four models described in Section <ref>. In the bottom panel of Fig. <ref>, we show the cumulative evidence ratios log(Z_blur)-log(Z_pow) and log(Z_pex)-log(Z_pow) that represent the cumulative distributions of the broad and narrow models normalized by the evidence of the simple powerlaw. It can be seen that the two curves start to flatten below S/N∼ 7, very similar to what is seen with the real data. We thus conclude that below this S/N value, little additional information is being added about the properties of the reflection, and henceforth restrict our analysis to a subsample with S/N>7. The number of sources meeting this S/N criterion are shown in Table <ref>.
  
  
  
   §.§ Compton reflection properties of the sample
  
  
  
  
  The results of the model comparison for the S/N-censored sample are also given are given in Table <ref>. The results are, in fact, rather similar to the full sample, as the number of very low S/N ratio sources is small. More specifically, we find that the BF method selects the model including both broad and narrow reflection components to be the best fitting model in the CDFS 7Ms, CDFN and AEGIS fields. 
  For the COSMOS field, the BF method selects the model with only distant reflection as the best-fitting model over the more complex model with both narrow and broad reflection. 
  Considering the total sample, comprising all four fields, the evidence for the broad reflection is very strong. The difference in the logarithmic evidence can be interpreted like a probability difference. Thus, we can see from Table <ref> that the model containing a blurred component is selected to be 10^26.8 more probable than the scenario with only a narrow reflection component to describe the sample with S/N larger than 7.
  The total evidence also shows that the model with a maximally spinning (Kerr) black hole has a probability 10^6.6 of being preferred over a non-rotating (Schwarzschild) solution, again seen also in the individual fields with the exception of COSMOS. 
  
  The fact that the COSMOS field shows a preference for narrow reflection only in the evidence comparison shown in Table <ref> does not necessarily imply that broad reflection is not present in the COSMOS source population. 
  This is because we test for the presence of the broad reflection in addition to narrower reflection from more distant material, e.g. the torus. As the model with broad reflection has more free parameters, this additional complexity is penalised in the evidence comparison. We therefore performed an additional test by fitting the spectra also with a model with a blurred reflection component but no additional narrow reflection component (see Table <ref>). This then tests whether there is evidence for broad reflection as an alternative to the narrow reflection. According to this test the preferred model is that with a blurred reflection component with a non-rotating (Schwarzschild) black hole. The preference for this model over the narrow reflection or maximally-spinning black hole is, however, marginal. We conclude that, while the COSMOS data are of sufficient quality to confirm the presence of Compton reflection in the spectra, they are not able to distinguish the properties of the reflection e.g. whether it is broad or narrow, or if broad the value of the black hole spin implied. Based on the fields with higher S/N ratio spectra, however, it seems most likely that both components are present also in the COSMOS data. 
  
  The average strength of both the narrow and broad reflection components, as measured by the R parameter, is an important diagnostic of the system, as it depends on the geometry and, in the case of the broad reflection, potentially also on relativistic effects close to the black hole <cit.>. The task of calculating the mean and sigma of the underlying parent population that describes the R value is not trivial, however, since the posterior distributions for this parameter are not always well described by a normal distribution. To address this, we us a Hierarchical Bayesian model (HBM), described in detail in Appendix <ref>. For a sample it models the intrinsic log(R) distribution of the sample as a Gaussian. Taking into account the posterior uncertainties on each individual object, the HBM fit returns mean and standard deviation σ of the distribution. The results are shown in Fig <ref>. Applying the HBM to sources with S/N>7 from all fields, we find a mean of log(R_pex) = -0.53 (thus, R_ pex = 10^-0.53=0.30) with spread σ_log(R_pex) = 0.2 for the narrow reflection component and mean of log(R_blur) = -0.57 (thus, R_ blur = 10^-0.57 = 0.27) with spread of σ_log(R_blur) = 0.14 for the blurred reflection component. The population mean values obtained with the HBM for the reflection fraction are similar for the narrow and broad components. The strength of the blurred reflection component in particular is smaller than would be expected from a flat disk illuminated by a point source, as has been found previously <cit.>. Both R values also show a significant spread of ∼ 0.2 dex, indicating there is considerable diversity in reflection strength within the population. 
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
   §.§ Dependence on other parameters
  
  
  Our expanded sample compared to that of <cit.> gives a more robust detection of the reflection components and thus offers the opportunity to investigate any dependence of the reflection on other parameters. To this end, we computed the fraction of sources with S/N>7 showing Z_blur > Z_pex. This is shown in Figure <ref> as function of luminosity, redshift and N_H bins. Errors are calculated assuming that the fractions follow a binomial distribution an the bin size is chosen to have the same interval size for the parameter on the x-axis.
  
  In Fig. <ref>, we see that the fraction of sources showing larger evidence for broad reflection decreases with increasing luminosity (see Figure <ref>, left panel). In the plot, we remove the 12 sources with luminosity less than L<10^41 erg/s. The excluded sources are the 12 with the lowest luminosity in the sample. They cover a 3 dex luminosity interval and hence the lowest luminosity bin in the Baldwin plot would be severely underpopulated with large error bars. In addition there is a minor concern that at these very low luminosities, sources may be partially contaminated by non-AGN emission. While this is likely to have a negligible effect on the overall sample properties, the reflection properties determined in this lowest luminosity bin may be somewhat unreliable. Moreover, including those sources does not change the conclusions from Figure <ref>.
  We also notice an anti-correlation with increasing redshift and a clear increase of the fraction of broadened sources with increasing N_H (see Figure <ref>, middle and right panel). 
  
  The anti-correlation of the equivalent width (EW) of the Fe Kα line and the luminosity of AGN is well-known characteristic <cit.>. This anti-correlation, called the Iwasawa-Taniguchi or X-ray Baldwin effect (Baldwin effect hereafter), has been seen both for the narrow core of the Fe Kα, and has also been claimed for the broad component of the line <cit.>. The anti-correlation between luminosity and fraction of sources selected as broad (Figure <ref>, left panel) could be a consequence of the Baldwin effect for the broadened component of the Fe Kα line.
  
  To study this phenomenon in more detail, we explore the relationship between the reflection fraction R (see Table <ref>) of both the broad and narrow reflection components in the most complex model and the luminosity of the sample sources (see Figure <ref>, left). The mean R values and their intrinsic dispersion were calculated using the HBM, as for the mean values for the whole sample, and exclude from the analysis the 12 sources with L<10^41 erg/s to avoid contamination from star forming galaxies.
  
  In Figure <ref>, we see a clear effect that both the narrow and broad Compton reflection fractions decrease significantly as a function of luminosity, from R∼ 0.5 at the lowest luminosities to R∼ 0.1 at the highest. Thus we confirm the existence of the Baldwin effect for both the broad and narrow components of the line.
  
  Fig. <ref> also shows the dependence of the reflection strengths with redshift and N_H. The same trends shown in Figure <ref> are seen, with a reduction in the R values with redshift, and a weak increase seen with obscuration. Both of these trends might be wholly or partially a consequence of the Baldwin effect, given the usual correlation between luminosity and redshift seen in flux-limited samples (Fig. <ref>, and the anti-correlation seen between luminosity and obscured fraction <cit.>.
  
  We test this hypothesis by splitting the sample (see Figure <ref>) into low-luminosity (L < 10^43.8 erg/s, blue in Figure <ref>) and high-luminosity (L > 10^43.8 erg/s, red in Figure <ref>) sub-samples. We chose this luminosity threshold because it splits the sample almost in half, with 1042 sources in the low-luminosity sample and 1120 in the high-luminosity one. The solid curves show the behaviour of the reflection fraction of the broad disk reflection component while the dashed lines show the R from the narrow torus reflection. We notice that the curves of R as function of redshift and N_H seem to flatten for higher luminosities, thus the trends of R vs. redshift and N_H might be indeed be mirroring the dependency of R with luminosity. Even in Figure <ref>, the stronger trend we can observe is the one with luminosity. In fact, the R values of broad and narrow reflection components at luminosities of L > 10^43.8 erg/s are consistently lower than the R values at L < 10^43.8 erg/s, as the X-ray Baldwin effect would predict.
  
  
  
  
     
  
  
  
  
  
  
  § DISCUSSION
  
  
  
  In this work we have used the deepest X-ray fields performed by Chandra to place constraints on 
  X-ray Compton reflection in a sample of AGN typical of the overall population, covering a luminosity range log L=41-45 erg/s out to z∼ 4. We build on the previous work by <cit.>, adopting the same Bayesian framework BXA to fit the spectra, determine parameters, and compare models. We confirm strong evidence for Compton reflection, and by implication also iron Kα emission, both from distant material, most likely the torus envisaged in orientation-dependent unification schemes, and relativistically broadened reflection modelled as arising from the inner accretion disk. 
  
  This relativistic reflection is expected to be ubiquitous in the standard scenario where the SMBH is surrounded by an accretion disk and a hot corona of electrons. As such, perhaps the most important result of our study is to confirm that paradigm in the general population of AGN, which are responsible for the bulk of black hole growth in the Universe, and the majority of the X-ray background radiation <cit.>. We observe that the model 
  with maximally spinning Kerr BH is preferred over a model with Schwarzschild BH, reproducing the result from <cit.>. This confirms that a portion of the reflection comes from very close to the black hole, indeed perhaps from within the 6 R_ g innermost stable circular orbit of a non-rotating black hole. In turn this implies both that a relatively cool accretion disk extends within this radius, and that the X-ray emission comes from the innermost regions, and is compact enough that a substantial proportion of the disk illumination is at these small radii.
  The result is supported by studies of microlensed quasars at high redshift. For example, <cit.> show evidence of high spin (a>0.8) in a sample of five lensed quasars at high redshift (z > 1.2) and an ultra-compact X-ray emitting region with size <10 R_ g.
  
  This provides one possible explanation for one of the more puzzling aspects of our analysis, that being the average strength of the relativistic reflection component. We find a value for the average reflection fraction of R_ blur∼ 0.3, contrasting with the R_ blur∼ 1 expected for a flat disk illuminated by a point X-ray source. This relatively weak blurred reflection, at least on average, is in agreement with observations of local AGN <cit.>, although some nearby AGN also show anomalously strong reflection <cit.>.  
  Both of these facts can be explained by strong relativistic effects and in particular light bending close to the central black hole. This can result in the reflection being either weaker, or stronger than that expected for a flat disc depending on geometrical considerations <cit.>. A reduced reflection strength would also be expected if the X-ray emission is beamed away from the disk <cit.>. 
  
  Apparently weak reflection may also be explained by ionization of the disk. Some reduction of the line flux is expected for moderate ionization due to resonant trapping, and at very high ionization parameters no line is produced at all once iron becomes fully ionized <cit.>. In this case the reflection continuum would also become hard to distinguish from the primary continuum. 
   
  There may also be geometric effects. The strength of the reflection component is maximised for our assumed geometry of a semi-infinite slab illuminated by a point source. If the real accretion disk-corona geometry is different to this then the reflection is reduced. For example, if the corona of hot electrons were at a height h comparable than the ISCO radius r_ms, fewer photons from the corona would intercept the disk. The divergence from the "lamp post" supported is corroborated by studies of microlensed quasars <cit.>, and furthermore suggest a very compact primary X-ray source, consistent with this idea. This effect would be further exacerbated if the accretion disk is truncated before the last stable orbit, although in this case the relativistic signatures would be less prominent, contrary to the strong evidence for their presence found in this work. 
  
  Our expanded sample compared to that of <cit.> has enabled an analysis of the dependence of the strength of reflection from both the torus and the accretion disk with luminosity (see Figure <ref>, left). We confirm an X-ray "Baldwin effect" in which the strength of the reflection component, and by implication the equivalent with of the iron Kα line, reduces with luminosity. 
  
  This effect is quite well established for the narrow core of the iron Kα line. The most common interpretation is that of the "receding torus" in which the covering fraction reduces with luminosity (e.g. Simpson 2005). This is in agreement with the observation of a higher fraction of optical type 1 galaxies, and lower prevalence of X-ray absorption, at high luminosities <cit.>. 
  
  The possible presence of a Baldwin effect for the broad part of the emission line was suggested by <cit.>. This should not have the same physical origin as the narrow-line X-ray Baldwin effect if the broad line comes from the accretion disk, rather than the torus. The near-absence of reflection in the highest luminosity objects may in part explain why the average reflection fraction in our sample is so low, and hence might be due to the same effects e.g. geometry, special or general relativistic beaming and/or disk ionization. A priori it is difficult to see where there should be a strong relationship between the disk-corona geometry and the luminosity, disfavouring this interpretation. On the other hand, if photons are beamed away from the disk this would result in an enhanced luminosity for a given object, along with weaker reflection, as observed. At higher luminosities, the disk may also be more highly ionized suppressing the iron Kα line and reflection continuum, as discussed above, and providing a natural explanation for the Baldwin effect. 
  
  Looking next to the apparent dependence of R with redshift, this seems most likely to be a consequence of the Baldwin effect, given the very strong luminosity-redshift correlation in our flux-limited samples. We therefore do not speculate further on the possibility of evolution of the disk-corona system over cosmic time, but if this can be confirmed with better data and samples it would certainly be an intriguing phenomenon. 
  
  We observe also that the fraction of sources selected as broad and the R value of both disk and torus reflection shows a mild increase with N_H (see Figures <ref> and <ref>, right panel).
  At face value this is the opposite of what might be expected, because in standard orientation-dependent unification schemes, absorbed sources should be seen at high inclination, whence the observed reflection signatures are weaker. The observed increase could also partially be a consequence of the Baldwin effect. As discussed above lower luminosity AGN are more likely to be absorbed, and when we restrict the luminosity range in our analysis the effect does indeed seem weaker. 
  
  A further effect is the possible degeneracy of the broad iron line with complex absortion in sources with high N_H. If the absorption is in fact more complicate than the simple model assumed here, then for moderately high values around log N_H=23, mismodelling could introduce apparent curvature in the continuum in the 5-6 keV range, mimicking a broad red wing to the line. Unfortunately, the quality of the spectra in our sample prevents us from constraining the disk inclination parameter for most of the sources. Thus we can neither confirm nor rule out the degeneracy of the obscuration with the inclination of the disk on this basis.
  
  One caveat to the above discussion is that is clearly challenging with spectra of the quality used in this work to decompose the reflection into its broad and narrow components. It can be noted from Figure <ref> (left), for example, that the luminosity dependence of the R value for the broad disk reflection and for the narrow reflection from the torus have a very similar behavior. Since the two components arise from very different regions around the SMBH, this may indicates that the component are not well decoupled in the model, and that there is considerable co-variance between them.
  
  Properly decoupling and measuring the properties of X-ray reflection in individual AGN at high redshift requires an X-ray telescope with significantly higher throughput then the current generation of instrument. Once launched, Athena <cit.> will provide this capability. Deep field osbervations with the Athena Wide Field Imager <cit.> will yield one or two orders of magnitude more photons per unit exposure than Chandra,
  giving high quality spectra for individual objects like those in our sample. Follow-up observations of selected objects with the Athena X-ray Integral Field Unit <cit.> of brighter examples found with the WFI will enable the first high resolution spectra of such objects. 
  
  On the other hand, with this project we have confirmed the potential of of X-ray spectroscopy combined with Bayesian inference to reveal information about the population properties of AGN, even with individual spectra of low signal-to-noise ratio. There are some limitations, in that we have also demonstrated that adding spectra with very low signal-to-noise ratio at some points fails to add further information. Once appropriate signal-to-noise ratio cuts are applied, however, each new objects added to the analysis is able to strengthen our inferences about the accretion processes in SMBHs and the gas structures around them. The recent launch of the instrument eROSITA <cit.> aboard the SRG satellite present a particularly exciting opportunity to take this forward. eROSITA will detect millions of AGN <cit.> spread over the full sky and filling out an extremely broad luminosity-redshift plane. This will open up new and exciting possibilities for the application of our methods.
  
  
  
  § SUMMARY AND CONCLUSIONS
  
  
  We present an analysis of the X-ray spectra of sources taken from the four Chandra deep fields CDFS 7Ms, CDFN, AEGIS and COSMOS to determine the Compton reflection and iron Kα line properties of typical AGN outside the nearby Universe. To this purpose, we fit all the spectra individually using BXA, rather than stacking them. We fit four models of increasing complexity, starting from a simple absorbed power-law and adding to this model a narrow reflection component () and subsequently a further relativistically blurred disk reflection () with dimensionless spin parameterfixed at two values, a=0 and a=0.998. We perform simulations to determine from which signal-to-noise ratio, S/N, a source adds information to the total Bayesian evidence of the sample, finding that below a S/N of 7 the sources are too faint to add any new information, so we restrict the bulk of our analysis to sources above this limit. The outputs of BXA are the best fit parameters for and the Bayesian evidence of the model, thus we can use the latter to calculate the Bayes factor (BF) for model comparison. Based on the fits to the individual spectra, we adopt a hierarchical Bayesian model (see Appendix <ref>) to determine the sample properties. 
  
  
  
  Our main findings are:
  
      
    * When considering the sample as a whole, the Bayesian evidence comparison shows a preference for a model containing both narrow and broad, relativistic Compton reflection. This is in agreement with the result of <cit.> using a smaller sample of spectra from the CDFS 4Ms.
      
    * As in <cit.>, we find that the broad disk reflection model with a maximally spinning BH is preferred over one with spin=0.   
      
    * The HBM shows that on average both narrow and broad Compton reflection is relatively weak. We find a mean of log(R_pex) = -0.53 with spread σ_log(R_pex) = 0.2 for the narrow reflection component and mean of log(R_blur) = -0.57 with spread of σ_log(R_blur) = 0.14 for the blurred reflection component. This implies a departure from the simple "lamp post" geometry assumed in our reflection models, and in the case of the broad reflection possibly light bending or ionization effects. 
      
    * We investigate the presence of an X-ray Baldwin effect in our sample, confirming a decrease in the reflection strength for both the distant and blurred components, and by inference both the narrow and broad components of the iron Kα line. The former may be explained by a "receding torus" model, whereas as the latter implies a dependence of the inner disk geometry and or ionization with luminosity. 
      
    * We also find anti-correlations of the reflection fraction of the disk and torus with redshift, and a weak positive correlation with N_H. Both may, however, be artifacts of the Baldwin effect. 
      
    * With this analysis, we confirm the power of Bayesian statistics to infer important physical characteristics and features of AGN using a sample of relatively low S/N X-ray spectra, a technique which can be applied powerfully to the upcoming eROSITA survey. To measure the properties of X-ray reflection in individual high-redshift AGN we will need instrument with significantly higher effective area, such Athena.
   
  
  
  
  
  § ACKNOWLEDGEMENTS
  
  
  We thank the referee, Francisco Carrerra, for a number of helpful suggestions which considerably improved this manuscript.
  
  This research made use of Astropy, a community-developed core Python package for Astronomy <cit.> and the NASA's Astrophysics Data System. This research made use of APLpy, an open-source plotting package for Python <cit.>. We also used extensively the Python package Matplotlib <cit.>.
  
  JB acknowledges support from the CONICYT-Chile grants Basal-CATA PFB-06/2007 & AFB-170002, FONDECYT Postdoctorados 3160439 and the Ministry of Economy, Development, and Tourism's Millennium Science Initiative through grant IC120009, awarded to The Millennium Institute of Astrophysics, MAS. 
  This research was supported by the DFG cluster of excellence ,,Origin and Structure of the Universe”.
  
  
  
  § DATA AVAILABILITY
  
  The data underlying this article will be shared on a reasonable request to the corresponding author. Part of the data (the posterior distribution of log(R)) can be found here: <https://github.com/blinda/HierarchicalBayes>.
  
  
  
  
  
  
  mnras
   
  
  
  
  
  
  
  
  
  
  
  
  § HIERARCHICAL BAYESIAN MODEL TO INFER THE INTRINSIC R DISTRIBUTION
  
  
  
  
  
  
  
  
  
  
  In this work, we have derived posterior distributions for each model parameter over a large sample (N>1000 objects). The parameter of foremost interest in this work is R, which we derived from fitting a model to X-ray spectra of AGN using BXA. 
  
  The posterior distributions are diverse and often wide, as illustrated in Figure <ref>. The question is now how to combine these uncertain posterior probability distributions to infer the intrinsic distribution.
  
  Hierarchical Bayesian modelling (HBM) can use the information we have on the single object (the posterior distribution) to infer the features of the parent population from which the objects are drafted <cit.>. In typical HBMs, the constraints per-object (e.g., its R posterior) is solved simultaneously with the parameters of the population distribution (e.g., the mean and standard deviation R distribution of the population). This is a complex, very high-dimensional (>1000-dimensional) problem. In practice, we address this problem in two simplified ways. They are both based on the same probabilistic model laid out below, but make use of the R posteriors already derived under flat priors. In one approach (Section <ref>), we use the posterior distributions as weights in numerical integration. In another approach, we solve the high-dimensional problem directly with Hamiltonian Monte Carlo (Section <ref>), but employ auxiliary distributions that approximate the posterior. The former approach has the benefit of not assuming a shape of the distributions. The latter approach avoids numerical sampling issues. In practice, the two methods of Section <ref> and <ref> show consistent results, which gives confidence in the method. The presented tools for Hierarchical Bayesian modeling are thus powerful and robust for inferring the intrinsic distribution given a large number of uncertain measurements, including upper limits.
  
  
  
  
  
  To derive the probabilistic model, we first assume that the parent distribution of the log(R) parameter is distributed as a Gaussian N(log R|μ, σ) with unknown mean μ and standard deviation σ. Adopting a different parent distribution shape, such as a skewed normal or a Beta distribution did not change the results significantly. For each object, we also have an a priori (before considering the data) unknown parameter R. We have already constrained this, as encoded in the posterior distributions P(log R| D_i) for each object. We can thus write the combined likelihood for a single object as:
  
  
      ∫ P(log R| D_i) N(log R|μ, σ) dlogR.
  
  
  Because the same parent distribution should hold for all objects, we multiply their probabilities and find the HBM likelihood:
  
  
      ℒ= ∏_i ∫ P(log R| D_i) N(log R|μ, σ) dlogR.
  
  
  Reusing the derived per-object posterior works here because we have adopted wide priors that are uniform over the integration variable in Eq.<ref> (log R).
  
  After adopting priors on μ (uniform) and σ (log-uniform), this forms a 2+N-dimensional Bayesian inference problem. To derive posterior distributions on μ and σ, we use two techniques explained in the following sections.
  
  
  
   §.§ Hierarchical Bayesian Model inference with Stan
  
  
  
  One way to solve Eq.<ref> is to fit for all N+2 parameters simultaneously. This requires advanced Hamiltonian Monte Carlo techniques which rely on likelihood gradients to navigate the search space.
  One issue is that we do not want to refit the spectra in this process.
  Therefore, to still allow each per-object R to vary according to its spectral constraints, we adopt an analytic approximation to its posterior.
  
  We first fit the posterior distributions of every object with a beta distribution. The free parameters are the shape parameters α and β and the location and width of the distributions. Since the parameter range from -2<log(R)<-1 does not have much physical sense and reflection fractions below 0.1 are virtually indistinguishable, we constrain the fit distributions to lie between -1 and 1.
  The posterior shapes differ depending on whether the parameter is well constrained, not constrained or an upper/lower limit (see Figure <ref>, left panels).
  Thus, we choose the beta distribution because it is flexible enough to fit reliably distributions with different shapes. The bottom panels of Figure <ref> show our best-fit beta approximations. 
  
  Next, we implement with the Hamiltonian Monte Carlo framework Stan[See <https://pystan.readthedocs.io/en/latest/>.] <cit.> a model that reads the parameters of all Beta distributions (vectors of α, β, location and scale). The model (Stan code in Listing <ref>) has free R parameters, which both follow these distributions and a parent normal distribution. The MCMC algorithm then simultaneously determines the posterior of the parameters (mean and sigma) of the normal parent distribution and that of the R values.
  
  For the sample with S/N>7 we obtain a mean of log(R_pex) = -0.53 with spread σ_log(R_pex) = 0.2 for the narrow reflection component and mean of log(R_blur) = -0.57 with spread of σ_log(R_blur) = 0.14 for the blurred reflection component (see Figure <ref>). For comparison, if we average the means of every single posterior distribution of the R values we obtain a mean of log(R_pex) = -0.64 with spread of σ_log(R_pex) = 0.38 for the narrow reflection component and mean of log(R_blur) = -0.6 with spread of σ_log(R_blur) = 0.26 for the blurred reflection component (see Figure <ref>).
  
  This method takes into account the large parameter uncertainties and upper limits. Since many of the posterior distributions for the R value have the shape of an upper limit (see Figure <ref>, middle), the mean of the population we obtain with a HBM is much smaller than the mean we would obtain by simply averaging the mean of every posterior distribution (see Figure <ref>, for the simple mean and Figure <ref> for the mean and sigma obtained with a HBM method). 
  
  The values presented in Figures <ref> and <ref> were calculated by applying the HBM to the subsamples of objects in 6 bins of luminosity, column density and redshift.
  
  [caption=Stan definition of a HBM where the input data has the shape of a beta distribution with parameters a, b, loc and scale and the model to be fit is a normal with parameters mu and sigma.]
  data 
      int<lower=0> N;
      vector[N] a;
      vector[N] b;
      vector[N] loc;
      vector[N] scale;
  
  parameters 
      real<lower=-1, upper=1> mu;
      real<lower=-2, upper=2> logsigma;
      vector<lower=0, upper=1>[N] u;
  
  transformed parameters 
      vector[N] x;
      real<lower=0> sigma;
      x = u .* scale + loc;
      sigma = pow(10, logsigma);
  
  model 
      u   beta(a, b);
      x   normal(mu, sigma);
  
  
  
  
  
   §.§ Numerical Hierarchical Bayesian Model inference
  
  
  
  
  
  Another approach is to use importance sampling to numerically simplify the problem to 2 parameters. In practice, we already have posterior samples R_i,j for each object i that approximate the (sometimes complex) posterior distributions. Therefore, we can write <ref>, dropping constant factors, with an importance sampling estimate:
  
      ℒ(μ, σ) ≈∏_i ∑_j N(log R_i,j|μ, σ).
  
  When using too few posterior samples, this approach can induce numerical noise into the population posterior. Care has to be taken when this approach is used for multi-dimensional integrations <cit.>.
  Akin to cross-validation, this could be further improved by using sub-samples of the posterior samples in Eq.<ref>, and averaging the estimators. However, by varying the number of posterior samples used from hundreds to thousands, we verified that for our problem this Monte Carlo one-dimensional integration is stable.
  
  The two-dimensional log-likelihood defined in Eq.<ref> is Monte Carlo sampled using [See <https://johannesbuchner.github.io/UltraNest/index.html>], a python nested sampling package developed in <cit.>. Figure <ref> shows the posterior of the mean and sigma of the normal distribution.
  For the total sample with S/N>7 we obtain a mean of log(R_pex) = -0.54 with spread σ_log(R_pex) = 0.18 for the narrow reflection component and mean of log(R_blur) = -0.59 with spread of σ_log(R_blur) = 0.14 for the blurred reflection component.
  
  
  
  
  
  
  
  