11institutetext: Department of Astronomy, University of Geneva, ch. d’Écogia
16, CH-1290 Versoix, Switzerland
11email<EMAIL_ADDRESS>22institutetext: Max-Planck-Institut für
extraterrestrische Physik, Gießenbachstraße 1, D-85748 Garching, Germany

# Reconstructing AGN X-ray spectral parameter distributions with Bayesian
methods II: Population inference

Lingsong Ge 11 Stéphane Paltani 11 Dominique Eckert 11 Mara Salvato 22

(Received xx, 2021; accepted xx, 202x)

We present a new Bayesian method to reconstruct the parent distributions of
X-ray spectral parameters of active galactic nuclei (AGN) in large surveys,
using the probability distribution function (PDF) of posteriors obtained by
fitting a consistent physical model to each object with a Bayesian method,
presented in our companion paper. The PDFs are often broadly distributed and
some present systematic biases, hence naive point estimators or even some
standard parametric modelling are not sufficient to reconstruct the parent
population without obvious bias. Our method uses a redistribution matrix
computed from a large realistic simulation with the same selection as in the
actual sample to redistribute the stacked PDF, and then forward-fit a non-
parametric model to it in a Bayesian way, so that the biases in the PDFs are
properly taken into account. In this way, we are able to get accurate
uncertainties on the model parameters. We apply our spectral fitting and
population inference methods to the _XMM_ -COSMOS survey as a pilot study. For
the 819 AGN detected in the COSMOS field, 663 (81%) of which having
spectroscopic redshifts (spec-z) and the others having high quality
photometric redshifts (photo-z), we find prominent bi-modality with widely
separated peaks in the distribution of the absorbing hydrogen column density
($N_{\mathrm{H}}$), and an indication that absorbed AGN have harder photon
indices. A clear decreasing trend of the absorbed AGN fraction vs the
intrinsic 2–10 keV luminosity is presented, but there is no clear evolution
trend of the absorbed fraction with redshift. Our method is designed to be
readily applicable to large AGN samples like XXL, and eventually eROSITA.

###### Key Words.:

galaxies: active – X-ray: galaxies – Methods: data analysis – Methods:
statistical

##  1 Introduction

Active galactic nuclei (AGN) represent the innermost regions of active
galaxies, where supermassive black holes (SMBHs) are accreting large amounts
of material from their surroundings, converting a large fraction of the
accretion energy ($\approx$5–42% depending on the spin of the SMBH, Kerr 1963;
Shapiro & Teukolsky 1983) into radiation throughout most of the
electromagnetic spectrum. AGN are among the most luminous steady sources in
the Universe, which makes them detectable out to high redshifts (Vito et al.
2017). All galaxies are thought to host a SMBH in their center (Magorrian et
al. 1998) and the majority of galaxies are believed to have gone through an
AGN phase in the past (Marconi et al. 2004; Merloni & Heinz 2008). For this
reason, determining the history of accretion onto SMBHs and their properties
is key to understanding galaxy evolution. The evolution of AGN activity also
allows us to probe the origin of the tight relationship between SMBH mass and
host galaxy properties (e.g. stellar mass) (Magorrian et al. 1998; Ferrarese &
Merritt 2000; Marconi & Hunt 2003; Häring & Rix 2004; Gültekin et al. 2009;
Kormendy & Ho 2013), and their similar growth history (Hopkins et al. 2006;
Silverman et al. 2008; Aird et al. 2015).

To see the full picture of the evolution of SMBHs and its impact on the host
galaxies, a complete demographic study of AGN is required. However, this is a
challenging task as most of the AGN ($\approx 90$%; Hickox & Alexander 2018)
are obscured in the optical and ultraviolet (UV) by the intervening gas and
dust, which is thought to consist mostly in a dusty torus (Krolik & Begelman
1988). A typical AGN consists in particular of the accretion disk, the broad-
line region (BLR) and the dusty torus. The BLR is a dust-free region close to
the central SMBH on sub-parsec scale, and produces broad emission lines (full
width at half maximum FWHM $\geq 2000$ km s-1) in the optical/UV band; they
are used to classify AGN as type-1 or type-2 depending on their presence or
absence in the optical spectrum (Antonucci 1993), and thus provide a tracer of
obscuration. The dusty torus is a geometrically and optically thick toroidal
structure well outside of the accretion disk and the BLR on a scale of several
parsecs, which largely obscures the direct observation of emission from the
central region of AGN when we observe through it (Krolik & Begelman 1988; Urry
& Padovani 1995). Therefore, optical surveys are heavily biased towards bright
type-1 AGN and very inefficient at detecting obscured AGN, even though the AGN
spectral energy distribution usually peaks in the optical/UV through thermal
emission from the accretion disk due to the gas viscosity (Shakura & Sunyaev
1973). Alternatively, mid-infrared (mid-IR) AGN surveys are less affected by
the obscuring torus compared to the optical, because they are sensitive to the
reprocessed emission from the torus that is observed in this band (Sanders et
al. 1989). However, in cases where the host galaxies are also dusty, the AGN
signal is diluted, restricting the mid-IR selected sample to luminous sources
(Donley et al. 2008).

Among the many different AGN detection methods, X-ray observations are one of
the most complete and least biased ways to select AGN (Brandt & Alexander
2015), thanks to the ubiquitous X-ray emission in AGN (Avni & Tananbaum 1986;
Brandt et al. 2000; Gibson et al. 2008), which always outshines that of the
host galaxy (Maccacaro et al. 1988), and to the strong penetrating power of
X-rays through gas and dust (Wilms et al. 2000). Both wide and deep X-ray
surveys have been conducted to study the evolution and population properties
of AGN up to redshift larger than 3 (e.g., Giacconi et al. 2002; Hasinger et
al. 2007; Elvis et al. 2009; Pierre et al. 2016; Luo et al. 2017, and many
others), taking advantage of the efficiency and completeness of X-ray
selection, especially in the hard X-ray band, i.e. above 2 keV.

The intrinsic X-ray emission of AGN is thought to be produced by the inverse
Comptonization of disk photons in a hot and optically thin corona in the
vicinity of the SMBH (Blandford et al. 1990; Zdziarski et al. 1995, 1996;
Krolik 1999), which appears in the X-ray spectrum as a power law, namely
$N(E)\propto E^{-\Gamma}$, with $\Gamma$ being the photon index. This
parameter is an important indicator of the AGN growth rate, since previous
studies have shown that the 2–10 keV intrinsic photon index is strongly
correlated with the Eddington ratio (Pounds et al. 1995; Risaliti et al. 2009;
Jin et al. 2012; Brightman et al. 2013), the ratio of the AGN bolometric
luminosity to the Eddington luminosity ($L_{\mathrm{Edd}}$), the maximum
theoretical luminosity of an accreting object. This correlation may also be
used to estimate the SMBH mass through the dependence of $L_{\mathrm{Edd}}$
(Shemmer et al. 2008).

A reliable measurement of the intrinsic photon index requires a careful
correction for the X-ray absorption due to the hydrogen column density
$N_{\mathrm{H}}$ along the line of sight, since these two parameters are
highly correlated (see figure 12 in Ge et al. 2021, submitted; hereafter Paper
I). Therefore, measuring the $N_{\mathrm{H}}$ distribution of the AGN
population is essential to constrain the distribution of $\Gamma$ and infer
the accretion history of SMBHs, as well as other demographic properties, such
as the evolution of the intrinsic X-ray luminosities ($L_{\mathrm{X}}$) (Ueda
et al. 2014; Aird et al. 2015; Buchner et al. 2015). In addition, it is also
important to constrain the structure and understand the physical properties of
the nuclear regions of AGN. For example, the unified AGN model (Antonucci
1993; Urry & Padovani 1995) successfully explains the different types of AGN
classified at optical wavelengths based on the orientation of the torus with
respect to the line of sight. However, more recent surveys found clear
evidence showing that the absorbed AGN fraction, defined as the fraction of
sources with $N_{\mathrm{H}}$ greater than $10^{22}$ cm-2, is decreasing with
increasing intrinsic X-ray luminosity (Ueda et al. 2003; Hasinger 2008; Brusa
et al. 2010; Merloni et al. 2014; Buchner et al. 2015), which points towards
differences between absorbed and non-absorbed AGN. To explore the physics
behind this dependency, Ricci et al. (2017b) analyzed a hard X-ray-selected
(14–195 keV) sample from the _Swift_ /BAT survey, showing that the absorbed
AGN fraction is intrinsically dependent on $L_{\mathrm{Edd}}$, which suggests
that the obscuration of AGN is mostly driven by the mass-normalized accretion
rate.

Besides the primary power law, X-ray spectra of AGN also exhibit a number of
secondary components, which makes the total spectrum very complex, including
an excess of soft X-ray emission below 1–2 keV (Arnaud et al. 1985; Singh et
al. 1985), a reflected component from the inner region of the AGN by the
circumnuclear material (Pounds et al. 1990; Nandra & Pounds 1994), a
scattering component by warm, fully-ionized gas relatively far away from the
center (Bianchi et al. 2006; Guainazzi & Bianchi 2007), and thermal emission
from diffuse plasma in the host galaxy (Iwasawa et al. 2011). Therefore,
extracting the main spectral parameters (e.g. $N_{\mathrm{H}}$, $\Gamma$ and
$L_{\mathrm{X}}$) without strong bias is not trivial considering the
complexity of the spectral components of AGN X-ray spectra, especially when
considering large samples of sources with low signal-to-noise ratios (S/N),
which is typical of large surveys. This point is addressed in Paper I, which
presents our method for AGN X-ray spectral fitting validated using
simulations.

In this paper, we develop a new Bayesian method to reconstruct the
distributions of spectral parameters in the parent AGN population based on the
PDFs that are obtained from the spectral fitting step. Since the PDFs from
fits to low S/N data can be biased due to the inhomogeneous sensitivity of the
likelihood dependent on the parameters (e.g., $N_{\mathrm{H}}$), we introduce
a redistribution matrix to calibrate this bias in the PDF and perform a
reliable population inference (see Sect. 4.2). We validate our method using
extensive simulations (Sect. 4.3) and apply it to the _XMM_ -COSMOS survey as
a pilot study, including reconstructing the distributions of $N_{\mathrm{H}}$
and $\Gamma$, as well as the absorbed fraction versus $L_{\mathrm{X}}$ and
redshift (Sect. 5). Our data reduction and source detection of the _XMM_
-COSMOS survey are introduced before presenting the method (Sect. 2).
Throughout the paper we assume the WMAP9 $\Lambda$CDM cosmology: $H_{0}=69.3$
km/s/Mpc, $\Omega_{\mathrm{\Lambda}}=0.72$, $\Omega_{\mathrm{m}}=0.28$,
$\sigma_{\mathrm{8}}=0.82$ (Hinshaw et al. 2013).

##  2 The _XMM_ -COSMOS survey

The _XMM_ -COSMOS survey centered at RA +150.11916667 (10:00:28.600), DEC =
+2.20583333 (+02:12:21.00) (Hasinger et al. 2007; Scoville et al. 2007)
consists of $\sim$60 ks observations of the 2 deg2 COSMOS field with XMM-
Newton. 1848 point-like sources are soft and/or hard X-ray detected
(Cappelluti et al. 2009). It is an excellent test bed for our new method,
because it has been studied extensively in various multi-wavelength surveys
(Brusa et al. 2010), thanks to which almost all sources have either
spectroscopic redshifts (spec-z) or high-quality photometric redshifts
(photo-z; Salvato et al. 2009, 2011; Laigle et al. 2016; Weaver et al. 2021).

###  2.1 Data preparation

We retrieved Observation Data Files (ODF) in the COSMOS field from the XMM-
Newton Science Archive, and found 56 pointings (see the detailed observation
log in Appendix. Reconstructing AGN X-ray spectral parameter distributions
with Bayesian methods II: Population inference). Starting from the ODFs, we
reduce the observation data using the data analysis pipeline described in
detail in section 2.2 of Ghirardini et al. (2019), which is based on the XMM-
Newton Science Analysis Software v13.5 (XMMSAS). emchain and epchain are used
to extract calibrated event files. The mos-filter and pn-filter executables
are run to filter out time periods affected by soft-proton flares. We also
measure count rates in unexposed corners of both instruments to compare them
with those inside the field of view (FOV) of the telescope, in order to
estimate the contamination of residual soft protons in the spectrum.

Count images are extracted from the three EPIC detectors in 9 bins: [0.6,
0.9], [0.9, 1.3], [1.3, 1.9], [1.9, 2.8], [2.8, 4.2], [4.2, 6.3], [6.3, 7.2],
[7.2, 9.2], [9.2, 12] keV. The 0.4–0.6 keV band is ignored because of the
unstable energy gain of the EPIC cameras at low energy, which renders the
energy calibration uncertain. The sizes of the bins are selected to be close
to a uniform distribution in log space, except for two of them, which are set
to 1.3–1.9 keV and 7.2–9.2 keV because of the Al K${\alpha}$ and Si
K${\alpha}$ fluorescence emission lines at 1.5 and 1.7 keV for both pn and
MOS, and those of Cu K${\alpha}$, Ni K${\alpha}$ and Zn K${\alpha}$ around 8
keV for pn (See XMM-Newton Users Handbook 111https://xmm-
tools.cosmos.esa.int/external/xmm_user_support/documentation/uhb/epicintbkgd.html).
We use eexpmap to compute exposure maps for the three detectors independently,
taking into account vignetting. For source detection, we extract images and
exposure maps in the 2–7 keV band, which allows us to be sensitive to absorbed
sources that are undetected at low energy. At the same time, restricting the
upper bound of the band to 7 keV enhances the sensitivity of the map, as the
ratio of effective area to background is unfavorable beyond 7 keV. The total
EPIC image is created by summing up count maps of all detectors, and its
corresponding exposure map is computed by summing up MOS exposure maps and
those of pn multiplied by the ratio of pn to MOS effective area averaged over
the energy band, which for a power law source with a photon index of 2.0 is
2.76. In the end we combine images and exposure maps of all pointings to
obtain the mosaics of the total field, which allows us to make full use of the
regions located on overlapping pointings.

###  2.2 Source detection

For point source detection we use the XMMSAS tool ewavelet on the 2–7 keV EPIC
mosaic. The wavelet scale is set to 2–4 pixels and the significance threshold
to 5. A total of 860 sources are detected. The much lower number of sources
compared to Brusa et al. (2010) is due to the selection in the hard band only,
while Brusa et al. (2010) also used the 0.5–2 keV band which is intrinsically
more sensitive but highly biased towards unabsorbed objects. We cross-match
our source list with the _Chandra_ COSMOS Legacy survey catalogue that
includes 4016 X-ray sources, 97% of which are optically and near-infrared
(NIR) identified (Civano et al. 2016; Marchesi et al. 2016). We find a match
with the _Chandra_ catalogue with a separation smaller than 10 arcsec for 828
sources. 5 of the remaining unmatched 32 sources are within 20 arcsec of
objects in the cluster catalogue of Gozaliasl et al. (2019). There are also 13
sources outside of the field of the _Chandra_ COSMOS Legacy survey. Therefore
we end up with 14 unmatched sources (1.6%), half of which are close to the
boundary of the _XMM_ -COSMOS survey, so that their positions could be
inaccurate. Among the 828 matched sources, there are 48 sources with 2
counterparts in the _Chandra_ catalogue. To select the real counterpart and
get the correct redshift information, we first look at their fluxes and
relative distances to the XMM-Newton source. If one of the _Chandra_
counterparts is 3 times brighter or 5 arcsec closer than the other one, we
select it as the real counterpart. In this way we identify the counterparts of
21 sources. For the remaining 27 we keep those having only one _Chandra_
counterpart within 3 arcsec, and ignore sources that have both counterparts
either within 3 arcsec or farther away than 3 arcsec, since it is difficult to
distinguish between the two of them. Our final sample contains 819 point
sources matched with a single _Chandra_ counterpart in the 2–7 keV band, 663
(81%) of which have measured spectroscopic redshifts.

Spectra of these sources are extracted directly from the 10-band mosaics using
our C++ aperture photometry tool xphot 222https://github.com/domeckert/xphot
as described in Paper I. The radius of the source region is fixed to 30
arcsec, and the background region is an annulus from 60 arcsec to 90 arcsec
away from the source center, to avoid remaining contamination from the source.
The areas covered by neighboring sources of the catalog are masked when
extracting the background spectrum.

##  3 Simulations

An accurate knowledge of the survey selection function is necessary to
reconstruct the true properties of the parent population. To determine the
selection function, we simulate a large set of sources with varying 2–7 keV
count rates. In addition, as will be explained in detail in Sect. 4, these
simulations will be used to determine the transfer function of our forward-
fitting approach, allowing us to reconstruct the parameter distribution in the
parent population.

Figure 1: Scatter plot of redshift vs intrinsic $L_{\mathrm{X}}$ in 2–10 keV
band; the histograms of these two parameters are plotted on the top and right
sub-panels. The simulated sample is plotted in orange and the _XMM_ -COSMOS
data sample is plotted in black. $L_{\mathrm{X}}$ of the data sample are
calculated using the medians of the PDFs from the results of the spectral
analysis.

###  3.1 Simulation approach

To characterize properly the selection function, we need to reproduce as
precisely as possible the process that leads to the source detection. We first
simulate the source count rate in the 2–7 keV band using the full model
described in Paper I. Source spectral parameters are drawn from their priors
(see Paper I), except that we use here an informative prior for $\Gamma$,
which is a Gaussian distribution with $\mu=1.95$ and $\sigma=0.15$ (Buchner et
al. 2014), to avoid generating completely unrealistic sources. The redshift is
computed from the comoving distance, that is randomly drawn in the comoving
volume from redshift 0 to 4, assuming AGN are uniformly distributed inside,
without evolution of the number density of AGN. At each redshift, the shape of
the distribution of the AGN luminosities is determined following the X-ray
luminosity function (XLF) of Aird et al. (2015) given by the flexible double
power law (FDPL) model.

Next, we randomly select a position in the real 2–7 keV band mosaic, excluding
a circular region of 10 arcsec radius around real sources detected in Sect.
2.2, and a 200-pixel-wide boundary, where the exposure is zero. Although,
ideally, real source regions are supposed to be included in the simulation, it
is quite complex to distinguish the simulated source from the real source when
they are detected as only one source; we therefore exclude sky areas that
would result in such confusion (chosen to be 10–arcsec radius, i.e. twice the
telescope’s FOV-averaged half energy width). Since the source regions defined
in this way only account for less than 1% of the total area, the impact on our
selection function is negligible. Source counts of simulated sources in the
detection band and in all energy bands are computed using the local exposure
time at the position in the exposure map, randomized by Poisson statistics and
further spatially distributed using the point spread function calculated using
the XMMSAS tool psfgen, at 3 keV and 8 arcmin off-axis, which is close to the
FOV-averaged and energy-weighted PSF in our survey. To avoid increasing
noticeably the source density in the image, we only simulate 100 sources at a
time, and repeat the process until enough sources are simulated. The detection
band image is then processed with ewavelet in the same conditions as for the
detection of real sources. The detected source list is then compared with the
simulated list in order to assess the probability of detection.

The spectra of the detected simulated sources are extracted using the same
approach as the one we used for the _XMM_ -COSMOS data by applying the xphot
executable to the twice nine mosaics for the pn and MOS instruments. In this
case, however, the source spectra are simply calculated by adding the Poisson
distributed source counts to the spectral bins extracted from the source
regions, which originally only contain background counts.

###  3.2 Properties of the simulated sample

Figure 2: Histograms of count rates (top) and S/N (bottom) of the simulated
sample (orange) and our _XMM_ -COSMOS (black) sample determined in the 2–7 keV
source detection band. 1-$\sigma$ errors on the histograms of the _XMM_
-COSMOS sample are shown in black. The errors on the simulated sample are not
presented since they are negligible in these plots.

We repeat the simulation process until we get over 5000 detected sources. 47
of them fall into boundaries with exposure from only one detector, so we
exclude them as we do not have such sources in our _XMM_ -COSMOS sample; 4980
sources finally remain in the simulated sample. The $L_{\mathrm{X}}$–$z$
distributions of the simulated sample and of the actual sample are shown in
Fig. 1. The medians of the PDFs of $L_{\mathrm{X}}$ are used in the latter
case (see Sect. 5).

The distributions of count rates and S/N in the source detection band of the
two samples are presented in Fig. 2. S/N is defined as

$\mathrm{S/N}=(N_{\mathrm{s}}-N_{\mathrm{b}}/r)/\sqrt{N_{\mathrm{s}}+N_{\mathrm{b}}/r^{2}},$
(1)

where $N_{\mathrm{s}}$ is the total count in the source region,
$N_{\mathrm{b}}$ is that in the background region and $r$ is the ratio between
background and source areas. Note that the S/N defined here is different from
what is used for source detection, estimated by ewavelet. The count rate
distributions match very well, except for one bin at low count rate around
$5\cdot 10^{-5}$ s-1. This is probably due to statistical fluctuations because
of the much smaller sample size of the _XMM_ -COSMOS data sample, as shown by
the error bars assuming they follow a binomial distribution. The two S/N
distributions also match well, as they both peak around 4.5 and decrease in
the same way below and above the peak. Note that S/N for both samples extend
below zero, because the background counts are estimated from the pixels in the
annulus outside of the source region.

Figure 3: True background count rates vs those estimated from background
regions around the sources of the simulated sample. MOS and pn count rates in
the 2–7 keV source detection band are plotted in blue and pink. The dashed
line is the 1:1 relation, and the dotted lines correspond to 15 percent
deviation from the reference line.

We can actually see how good the background estimation is, since we know the
true background of the simulated sources. In Fig. 3 we show the scatter of the
estimated background count rates versus the true ones for both detectors in
the 2–7 keV detection band. The background estimation is mostly constrained
within 15% uncertainties without noticeable bias for both pn and MOS. However,
there are quite expectedly a few outliers that are the cause of most of the
spectra with negative S/N.

###  3.3 Selection function

The selection function describes the probability of observing a source given
its spectral properties. Since the detection is only a function of the number
of counts in the detection band compared to the number of background counts,
we parameterize the selection function as a function of the count rate $c$.
The effect of varying background counts and exposure times is taken into
account by the randomization of the position of the simulated sources in the
simulated images. The parametrization is made using an error function as
follows:

$P(c|a,b)=0.5\operatorname{erf}a(\log(c)+b)+0.5,$ (2)

where $a$ and $b$ are free parameters. We use a simple least-squares
regression to fit this function to the detected fraction of simulated sources
at different count rates $c$. From the fit we obtain $a=3.0\pm 0.3$,
$b=3.79\pm 0.02$, and the results are shown in Fig. 4. We also compute the
selection function on a single pointing in the central region of the field for
comparison, taking into account the position-dependent PSF, whereas on the
mosaic we use a FOV-averaged PSF. As shown in Fig. 4, if we shift the
selection function of the single pointing to the same depth as that of the
full FOV (same parameter $b$), the shapes of the two match quite well.
Obviously, by combining the data of all the observations we reach much deeper
depths compared to the single pointing, but we see that the effect of using a
FOV-averaged PSF does not significantly affect the shape of the selection
function.

Figure 4: Selection function of the full mosaic (black solid line) fitted to
the detected fraction of our simulated sources (black). The selection function
in the central pointing is shown as an orange solid line. The dashed orange
line shows the central-pointing selection function shifted to the same depth
as that of the full FOV (same parameter $b$ in Eq. (2) for comparison.

##  4 Reconstruction of the parent population’s properties

To infer the distributions of the spectral parameters in the parent population
using the PDFs of each source, the simple approach is to use a point
estimator, for instance the median, or to stack the PDFs. The former can be a
satisfactory approach only in cases where the PDFs are well constrained in a
single peak, which requires good data quality and non-degenerate models. If
the PDFs are widely distributed in the parameter space or are multi-modal,
using point estimators obviously can not extract enough information from the
full PDFs to perform a correct population inference. Stacking PDFs seems
straightforward, but it is actually statistically meaningless, because
probabilities of independent events are not additive, but multiplicative (Malz
2021). So stacking PDFs can only give us a rough estimate of the parent
population, and not a proper inference of the parent population. In addition,
stacked PDFs contain uncertainties due to the measurement process, and provide
too broad parent distributions (Loredo 2004). We explore here two
statistically correct methods to reconstruct the parent distributions, one of
them being an entirely new approach, as far as we know.

###  4.1 Parametric fitting method

The first method is to parameterize the parent distribution and then fit the
model to the PDFs. We follow here J. Speagle’s approach
333https://github.com/joshspeagle/frankenz/blob/master/demos: 5 – Population
Inference with Redshifts.ipynb. Suppose our sample is $\\{s\\}$, and for each
source $s_{i}$ we have an associated parameter estimate $\theta_{i}$ with PDF:
$P(\theta_{i}|s_{i})$. Since our goal is to reconstruct the parent
distribution of parameter $\theta$ given $\\{s\\}$, which can be denoted as
$\phi\equiv P(\theta|\\{s\\})$, we need a likelihood function for the PDFs,
which we can write as $P(\\{p_{\theta}\\}|\phi)$. Now using Bayes theorem, we
can write the posterior as:

$P(\phi|\\{p_{\theta}\\})\propto P(\\{p_{\theta}\\}|\phi)P(\phi)$ (3)

Assuming the PDFs are independent, the likelihood can be rewritten as the
product of all the PDFs marginalized over the parameter space as following:

$P(\\{p_{\theta}\\}|\phi)=\prod_{i}P(p_{\theta_{i}}|\phi)=\prod_{i}\int
P(\theta_{i}|s_{i})P(\theta|\\{s\\})d\theta$ (4)

Namely, the posterior probability of the population distribution is determined
by maximizing its overlap with each PDF $P(\theta_{i}|\theta)$.

Loredo (2004) includes another term in the likelihood function (equation 6),
which asserts that no other sources are detected elsewhere in the parameter
space, to account for the source uncertainties. We choose to simplify it and
ignore this term for the parametric fitting method, because almost all the
parameter space is covered by the PDFs of detected sources (espcially for
$\Gamma$, where this method is applied as presented below), thus the
additional term does not make much difference here.

We infer the parent distribution of the photon index $\Gamma$ as an example of
the application of this method. Since it has been observed that $\Gamma$
follows approximately a Gaussian distribution (Nandra & Pounds 1994; Ueda et
al. 2014; Ricci et al. 2017a), we can parametrize its parent distribution with
a Gaussian with only two free parameters, i.e., mean $\mu$ and standard
deviation $\sigma$. We apply non-informative priors on the two parameters,
including a uniform prior between 1.5 and 2.5 for $\mu$ and a Jeffreys prior
between 0.01 and 1 for $\sigma$. With the likelihood function and prior we
sample the posterior using MultiNest (Feroz & Hobson 2008; Feroz et al. 2009,
2019), which is an implementation of the nested-sampling algorithm (Skilling
2004). MultiNest has the capability to sample multi-modal posteriors, which is
essential to analyse low S/N data with complex models, since more common
sampling approaches like Markov Chain Monte Carlo (MCMC) often have
difficulties in exploring a multi-modal parameter space. We test this approach
using a sub-sample of 660 sources from the total simulated sample introduced
in Sect. 3.2. The true $\Gamma$ distribution follows a Gaussian with
$\mu=1.75$ and $\sigma=0.15$. We select this relatively small mean value of
$\Gamma$, so that it is significantly different from the median of the prior
of $\Gamma$, which is applied for the spectral fitting.

Figure 5: Reconstructed parent distribution of $\Gamma$ of a simulated sample
using the parametric fitting approach (see Sect. 4.1) in blue, with the filled
area showing the 1-$\sigma$ credible interval. The true distribution is shown
as the orange histogram, and the results using medians and stacked PDF are
shown with dashed green (hatched) and magenta lines. The distribution of ML
values is shown with a dash-dotted red line, and the dotted black line is the
prior applied to $\Gamma$ during the spectral fitting step.

The result of the reconstructed Gaussian parent distribution is shown in Fig.
5, as well as the inferred distributions using medians, stacked PDF and best-
fit values from a maximum likelihood method (ML; see Sect. 3.5 in Paper I for
details). The parametric approach in general obtains a good fit with
$\mu=1.82\pm 0.02$, and $\sigma=0.15\pm 0.01$, except that its mean is about
$3\sigma$ too high. The origin of this bias can be two-fold. First, there is a
positive bias in the PDFs at low values of $\Gamma$, which is increasing when
the true value deviates further from the value of 2.0 used as the median of
the prior (see Fig. 10. in Paper I). This induces a bias in the fitted model
as the parametric method cannot deal with biased PDFs. Second, the decrease of
the prior of $\Gamma$ below 1.5 may bias the sources with very hard photon
indices towards softer values, thereby shifting the mean of the fitted
Gaussian. However, the effect is probably insufficient to explain the
difference of more than $3\sigma$ with respect to the true value.

On the other hand, both the point estimator using the medians of the PDFs and
the stacked PDF have trouble recovering the true distribution. The stacked PDF
is too broad, giving too much probability for $\Gamma>2$. As already hinted
above, the main reason is that each PDF is broadened by the uncertainty in the
determination of $\Gamma$ (Loredo 2004), especially for absorbed sources whose
PDFs are hardly constrained (see Fig. 12. in Paper I) and close to the prior
applied on $\Gamma$ during the spectral fitting process. On the other hand,
the medians give a more constrained distribution compared to the stacked PDF.
However, it is still much broader than the true one and, more importantly, it
peaks in the bin $1.9\leq\Gamma<2$, which is significantly higher than the
true mean value (1.75), outside of the 3-$\sigma$ range. The too high location
of the peak is expected, because the median of $Gamma$ for objects with
completely unconstrained posteriors (mostly absorbed) is 2. The result of ML
for spectral fitting shows that it cannot be used to reconstruct the parent
population, as discussed in Paper I, because of the numerous ML fits at the
domain boundaries. In this case it is clear that the parametric fitting method
is much better in recovering the parent distribution when we have loosely
constrained PDFs.

###  4.2 Redistribution matrix population inference (ReMPI)

The parametric method seems to work well for the reconstruction of the parent
distribution of $\Gamma$. However, it assumes a parametric form for the parent
distribution, which may not be adopted to parameters like $N_{\mathrm{H}}$,
for which no obvious parametric form exists. Moreover, the method implicitly
assumes that the PDFs provide the true posterior probabilities of the
parameters, i.e. that they are not biased, which is not the case for
$N_{\mathrm{H}}$. As we discussed in Paper I, the PDF of $N_{\mathrm{H}}$ for
an unabsorbed object is mostly flat below around $10^{22}$ cm-2, and steeply
decreases above it, since the spectrum quickly becomes insensitive to this
parameter in the unabsorbed regime. A similar effect can also be seen for
highly absorbed sources. As a result, the PDFs of $N_{\mathrm{H}}$ of many
sources are highly biased. For instance, the stacked PDF of sources in our
simulated sample that have the true $N_{\mathrm{H}}$ below $10^{20.2}$ cm-2 is
almost uniformly distributed below $10^{21.7}$ cm-2, and more importantly, the
medians of the PDFs are completely off, peaking around $10^{21}$ cm-2, as
shown in Fig. 6. In this case the parametric method is not able to correct
such a bias. On the other hand, both the stacked PDF and the medians of
objects with $N_{\mathrm{H}}$ in the range of $10^{22.7}$ – $10^{23.2}$ cm-2
are much more symmetric, and their peaks are located at the correct position.
Therefore, the parametric method may reconstruct correctly the population
properties for Compton-thin absorbed sources, but not for the less absorbed
ones.

To overcome this difficulty we develop here a new method using a
redistribution matrix to characterize the bias in the PDF. This method is
inspired from the forward-fitting of X-ray spectra through a response, which
is the standard method to analyze X- and gamma-ray spectra with xspec (Arnaud
1996), Sherpa (Freeman et al. 2001), SPEX (Kaastra et al. 1996), etc. X- and
gamma-ray spectral analysis faces indeed a very similar issue, as the recorded
energies of mono-energetic photons present a tail at low energy, biasing the
median recorded energy towards low values. Here, we forward-fit the stacked
PDFs assuming a histogram model for the parent distribution.

Figure 6: Stacked PDFs (hashed) and median distributions (filled) of the
simulated sources with true $N_{\mathrm{H}}$ below $10^{20.2}$ cm-2 (first
bin) shown in orange, and those with $N_{\mathrm{H}}$ from $10^{22.7}$ to
$10^{23.2}$ cm-2 (7th bin) shown in blue. The true $N_{\mathrm{H}}$
distributions of the 2 bins are shown in black (amplitude fixed to 1 for
clarity).

####  4.2.1 Redistribution matrix

Figure 7: Stacked PDFs of sources with true $N_{\mathrm{H}}$ in a single
$N_{\mathrm{H}}$ bin in 4 different redshift sub-samples (top 4 panels), and
in 4 S/N sub-samples (bottom 4 panels). The orange area shows the true
$N_{\mathrm{H}}$ (amplitude fixed to 1 for clarity), and the black line shows
the stacked PDF. Note that sources with negative S/N are formally detected by
ewavelet, since the background is estimated differently.

The classic X-ray spectral fitting technique uses a redistribution matrix
(RMF) that provides the probability of measuring a given energy as a function
of the energy of an incident photon. The RMF can be determined from ground-
based calibration experiments. Similarly, in our case we can determine the
redistribution matrix by computing the PDFs of simulated objects.

Nevertheless, the situation here is more complicated than in the case of the
spectral analysis of X-ray photons, because the redistribution is not only
dependent on the parameter itself, but also on the other spectral parameters.
For example, the PDF of $\Gamma$ for a given object highly depends on
$N_{\mathrm{H}}$. Theoretically, the redistribution matrix should be computed
for every possible set of parameters and provide the redistribution of all of
them simultaneously. This would result however in a high-dimensional matrix
that would be extremely expensive, first to compute and then to use in the
fitting procedure. We can however greatly simplify the problem by
marginalizing over the parameters that are not of immediate interest. For
instance, if one wants to study the distribution of $\Gamma$, we can simulate
a number of sources with fixed $\Gamma$, covering the whole range of interest,
the other source parameters being chosen at random following prior
distributions. Thus this process remains essentially Bayesian, since it
depends on the priors of the marginalized parameters.

The redistribution matrix $M_{\mathrm{red}}$ is defined in the following way.
We represent the true parent population distribution of the parameter of
interest as a 1-dimensional histogram with an n-bin step function expressed by
a vector $A_{\mathrm{parent}}=[p_{1},p_{2},...,p_{n}]^{\mathrm{T}}$, $\sum
p_{i}=1$, where T is the transpose operator and $p_{i}$ represents the
integrated probability in the $i^{\mathrm{th}}$ bin, and another vector for
the stacked PDF $A_{\mathrm{stacked}}=[s_{1},s_{2},...,s_{n}]^{\mathrm{T}}$,
$\sum s_{i}=1$. Then $M_{\mathrm{red}}$ is defined as an $n\times n$ matrix
such that:

$A_{\mathrm{stacked}}=M_{\mathrm{red}}\times A_{\mathrm{parent}}.$ (5)

Since this method uses matrix operations, a histogram model is required to
model the parent distribution. However, a parametric parent distribution can
be used, but it would need to be binned to create $A_{\mathrm{parent}}$. If we
have a joint distribution of multiple parameters, the formalism can be also
extended to higher dimensions, where $A_{\mathrm{stacked}}$ and
$A_{\mathrm{parent}}$ would become $N-$dimensional arrays and
$M_{\mathrm{red}}$ would have dimension $N+1$.

The parent distribution of $N_{\mathrm{H}}$ depends in particular on the
redshift and the data quality, namely S/N. Since our spectral energy bands are
fixed in the 0.6–12 keV range, sources with different redshifts are observed
in different energy bands in their local frame. This impacts the determination
of $N_{\mathrm{H}}$, because photoelectric absorption is more efficient at low
energies. When the spectrum is more and more redshifted, $N_{\mathrm{H}}$ is
more and more difficult to constrain, except for very high absorption. This
effect is shown in the top four panels of Fig. 7, taking as an example the
stacked PDFs of $N_{\mathrm{H}}$ around $10^{22.5}$ cm-2 in four redshift
bins. For sources with redshift smaller than 1, we clearly identify a peak
around the true $N_{\mathrm{H}}$ value in the stacked PDF, but it decreases
quickly as redshift increases, and more and more probability is redistributed
to other bins. The S/N also strongly affects the PDFs, as shown in the bottom
four panels of Fig. 7. The stacked PDF of sources with very high S/N (¿7) is
very well determined, with only a little leakage to lower $N_{\mathrm{H}}$
resulting from a few highly redshifted sources, while those with low S/N show
a much less prominent peak at the true $N_{\mathrm{H}}$ and have clearly
broader distributions. In particular, the stacked PDF of sources with negative
S/N (although formally detected by ewavelet), where the background dominates,
is almost flat across the range.

In order to construct $M_{\mathrm{red}}$, we use the simulated sample
presented in Sect. 3, which has been built so that its basic properties match
as well as possible those of the real _XMM_ -COSMOS data; in particular, Figs
1 and 2 show the excellent agreement between the redshift and S/N
distributions.

####  4.2.2 Forward-fitting

Given $M_{\mathrm{red}}$ and the PDFs of the sources (and hence
$A_{\mathrm{stacked}}$), the next step is to find the parent distribution
$A_{\mathrm{parent}}$ through a forward-fitting procedure. First, we draw the
parameters $p_{i}$ of $A_{\mathrm{parent}}$ from the prior distribution
$P(A_{\mathrm{parent}})$. The prior we choose is an uninformative, flat
histogram, i.e. ${p_{i}}$ represents the coordinates of a simplex in $n$
dimensions ($n$ being the number of bins), whose vertices are the $n$ unit
vectors. This can be achieved using a Dirichlet distribution
$\mathrm{Dir}(\alpha_{i})$, with $\alpha_{i}$ accounting for different sizes
of the binning. In our case the binning is uniform, so $\alpha_{i}=1,\forall
i$. Then we apply the redistribution matrix $M_{\mathrm{red}}$ to achieve the
predicted distribution $A_{\mathrm{red}}=M_{\mathrm{red}}\times
A_{\mathrm{parent}}$.

Secondly, we need to determine the statistical properties of
$A_{\mathrm{stacked}}$, so that we can compute the likelihood, which can be
expressed as $P(A_{\mathrm{stacked}}|M_{\mathrm{red}},A_{\mathrm{parent}})$,
which is equivalent to $P(A_{\mathrm{stacked}}|A_{\mathrm{red}})$. Since the
theoretical distribution of the stacked PDF is unknown, we rely on the data
sample itself to infer their statistics. We perform a bootstrap of the stacked
PDF by drawing individual PDFs with replacement. We compute $10^{6}$
bootstrapped stacked PDFs, allowing us to construct the full statistical
distribution of $A_{\mathrm{stacked}}$. An example of the distributions of
$s_{i}$ of the 10-bin $N_{\mathrm{H}}$ stacked PDFs from a bootstrapped data
sample is shown in Fig. 8. We find that the distributions of $s_{i}$ are
reasonably well approximated with normal distributions. Thus we can determine
Gaussian uncertainties for all $s_{i}$ parameters, as well as their full
covariance matrix using the bootstrapped sample. Finally, we sample the
posterior $P(A_{\mathrm{stacked}}|A_{\mathrm{parent}})$ with the No-U-Turn
Markov Chain Monte Carlo sampler (NUTS) implemented in the Python package
PyMC3 (Salvatier et al. 2016), because of its convenient implementation of the
Dirichlet prior and faster sampling process compared to MultiNest, and we
expect the chance of converging towards local maxima to be small.

Figure 8: Distributions of the parameters $s_{i}$ of $A_{\mathrm{stacked}}$ in
the bootstrapped data sample. The x axis shows the value of $s_{i}$, and the
y-axis shows the number of repetitions. The curves with different colors
correspond to the different bins.

###  4.3 Validation of ReMPI

####  4.3.1 Hydrogen column density $N_{\mathrm{H}}$

We first test this method on $N_{\mathrm{H}}$ using a randomly selected sub-
sample of the simulation with the same number of sources as the _XMM_ -COSMOS
data sample, keeping the same redshift and S/N distributions, so that the same
redistribution matrix can be used. However, we make its true $N_{\mathrm{H}}$
distribution follow a double Gaussian distribution, as shown in the lower
panel of Fig. 9. The reason for choosing this distribution is two-fold. First,
with such complex distribution we can clearly see the performance of the
method on absorbed and unabsorbed sources; and second, it is suggested that
there may be a similar bimodality in the $N_{\mathrm{H}}$ parent distribution
of real sources (Paltani et al. 2008; Ueda et al. 2014; Ricci et al. 2017a).

Figure 9: Lower panel: reconstructed parent distributions of $N_{\mathrm{H}}$
for a sub-sample sharing the same size as the _XMM_ -COSMOS data sample and
being randomly selected from the simulation. The true distribution is plotted
as the filled histogram in orange, and the inferred $N_{\mathrm{H}}$
distribution obtained by applying ReMPI is shown as the black line. The error
bars show the 68.3% credible intervals. For comparison, the stacked PDFs, the
distributions of medians and the ML results are plotted in magenta dashed
line, green dashed line and red dash-dotted line respectively. Upper panel:
corner plot of the parameters of the first four bins ($p_{1}$ – $p_{4}$,
converted to probability density for consistency with the lower panel) in the
unabsorbed regime.

The inferred parent distribution by ReMPI is shown in the same figure,
together with the results obtained using medians of the PDFs, the stacked PDFs
and ML. In the absorbed regime, our approach successfully reconstructs the
parent distribution of $N_{\mathrm{H}}$ with fairly small uncertainties, while
the other methods either have the peak at a wrong place (medians) or severely
underestimate it due to the broadening of the distribution. The bin of
Compton-thick sources (CT; corresponding to sources with $N_{\mathrm{H}}$
above $1.6\cdot 10^{24}$ cm-2) and the valley around $10^{22}$ cm-2 are also
very well constrained, clearly recovering the bimodality. In the unabsorbed
regime, the uncertainties become much larger and strongly correlated, as shown
in the corner plot (Foreman-Mackey 2016) of Fig. 9. The first three bins
($p_{1}$ – $p_{3}$) are basically recovering the Dirichlet prior, while
$p_{4}$ starts to be better constrained and decoupled from $p_{1}$ and
$p_{2}$. However, despite the little information that we obtain for unabsorbed
sources, the result correctly indicates that it is quite difficult to
constrain the distribution of $N_{\mathrm{H}}$ below $10^{22}$ cm-2. By
contrast, the distribution of medians gives us a false peak at
$N_{\mathrm{H}}\sim 10^{21}$ cm-2, while the stacked PDF is not able to
distinguish the bimodality, since it is simply flat until almost $10^{23}$
cm-2 and has a much lower peak. The ML results are even worse, with a very
high peak in the lowest bin due to the fact that almost all unabsorbed sources
are fitted with $N_{\mathrm{H}}$ at the lower boundary.

####  4.3.2 Photon index $\Gamma$

To reconstruct the distribution of the photon index $\Gamma$ with ReMPI, we
need to extend the redistribution matrix to two parameters, because the
measurement of $\Gamma$ strongly depends on $N_{\mathrm{H}}$, as shown in
figure 11 of Paper I, and the flat $N_{\mathrm{H}}$ distribution assumed in
the simulation is probably very different from the reality. While in principle
using more bins could be more accurate, here we simply use an $8\times 2$ two-
dimensional histogram model for the joint distribution of $\Gamma$ and
$N_{\mathrm{H}}$, because our data sample has a limited number of sources and
our main goal for this exercise is to reconstruct the distribution of
$\Gamma$. The two bins in the $N_{\mathrm{H}}$ dimension are separated at the
boundary between absorbed and unabsorbed sources, at $N_{\mathrm{H}}=10^{22}$
cm-2, while for $\Gamma$ the bins are uniform from 1.2 to 2.6.

Figure 10: Lower panel: Same as Fig. 5, but we show in addition the
reconstructed distribution of $\Gamma$ using ReMPI in black, with $\Gamma$
binned to 8 uniform bins in accordance with the model used in ReMPI. Error
bars show the 68.3% credible intervals. Upper panel: corner plot of the
parameters of the second to the fifth bins ($p_{2}$ – $p_{5}$, converted to
probability density for consistency with the lower panel), which include the
vast majority of the sources.

Figure 11: Reconstruction of the absorbed AGN fraction vs redshift (top row)
and $L_{\mathrm{X}}$ (bottom row) for different sub-samples of the simulation.
The black line represents the true absorbed fraction relation in each redshift
or $L_{\mathrm{X}}$ bin, and the reconstructed one obtained by applying our
method is plotted with red dots and vertical error bars showing the 68.3%
credible intervals. The results of using stacked PDF and medians are shown in
yellow and green respectively. The grey histogram represents the
$L_{\mathrm{X}}$ or redshift distributions of the simulated sample.

We test ReMPI on the same sample as introduced in Sect. 4.1, and the result is
illustrated in Fig. 10, together with the one using the parametric method. The
true $\Gamma$ distribution is correctly recovered by ReMPI and the result in
general agrees with that of the parametric one. However, our method corrects
for biased PDFs whereas the parametric method does not. This allows ReMPI to
recover the peak more precisely, as the median of the distribution
reconstructed by ReMPI is around the true value (1.75), whereas the parametric
method biases $\Gamma$ towards larger values, as discussed in Sect. 4.1.

On the other hand, the uncertainties on the histogram distribution of ReMPI
are relatively larger than those on the Gaussian distribution, especially for
the two bins around the peak, but this is expected, because ReMPI has many
more free parameters ($8\times 2$ versus 2), since the parent distribution is
non-parametric and specified with an arbitrary histogram. Also, the assumption
of the Gaussian parent distribution is a much stronger prior than the non-
informative one on the histogram model. It works well when the parent
distribution is indeed a Gaussian, like this test sample, but with ReMPI we
keep the possibility of finding alternative distributions. Actually, in ReMPI
we have the freedom to adjust the informative level of the prior by changing
the concentration parameter $\alpha$ of the Dirichlet distribution. For
example, if we believe that very few objects would have $\Gamma$ close to the
boundaries of the parameter space, we can increase $\alpha$ to give more
weight to the bins around the center and penalize those close to the
boundaries, imposing an approximate Gaussian prior on the parent distribution
of $\Gamma$.

The corner plot (Foreman-Mackey 2016) in Fig. 10 shows the PDFs of $p_{2}$
through $p_{5}$, which contain all the sources, and their correlations. The
adjacent bins are strongly anti-correlated, especially the two bins around the
peak, which is a natural outcome of applying the redistribution matrix. On the
other hand, bins that are not neighbors are clearly uncorrelated. This shows
that the matrix of $\Gamma$ is not far from being diagonal, and mainly
redistributes the probability into the neighboring bins.

####  4.3.3 Absorbed fraction vs $z$ and X-ray luminosity $L_{\mathrm{X}}$

We also test the performance of ReMPI to reconstruct the evolution of the
absorbed fraction with redshift, as well as its dependence on the 2–10 keV
intrinsic X-ray luminosity ($L_{\mathrm{X}}$). The absorbed fraction is
defined as the fraction of sources whose spectral parameter $N_{\mathrm{H}}$
is larger than $10^{22}$ cm-2. In this specific case we focus only on the
classification between absorbed and unabsorbed sources. Therefore, similarly
to $\Gamma$, we only define two bins in the $N_{\mathrm{H}}$ space with
boundary at $10^{22}$ cm-2. The redistribution matrix needs however to be
constructed for several bins in redshift and luminosity, respectively.

To infer the redshift evolution of the absorbed fraction, we divide the total
sample into different redshift bins, and then reconstruct the $N_{\mathrm{H}}$
distribution modelled by the 2-bin histogram model in each redshift bin. The
redistribution matrices of different redshift bins are also derived from the
total simulated sample divided into the same bins. In order to validate this
approach, we randomly draw three sub-samples from the simulation, all of which
share similar redshift and S/N distributions to those of our _XMM_ -COSMOS
data sample, as well as the same sample size. We impose an increasing absorbed
AGN fraction over redshift in the first sample, while the other two have a
decreasing and a flat relation, respectively. We make six bins in the redshift
space, with the first five bins being equally spaced with a bin width of 0.5,
and the last bin having a width of 1.1 to encompass the most distant source
($z=3.6$). This binning is chosen to optimize the resolution in redshift, as
well as the number of sources in each bin, for our _XMM_ -COSMOS sample. The
results, as well as those using medians or stacked PDFs, are shown in the
first row of Fig. 11. For the medians it is straightforward to calculate the
absorbed fraction, while for the stacked PDF, we use the stacked posteriors
integrated above $10^{22}$ cm-2 as an estimate for the absorbed fraction. We
see that, using medians or stacked PDFs, we tend to underestimate the slope of
the relation. This is particularly clear in the case of an increasing absorbed
fraction with redshift, where at high redshift the stacked PDF and medians
substantially underestimate the absorbed fraction (upper left panel of Fig.
11), due to the increasing difficulty of constraining $N_{\mathrm{H}}$ as
redshift increases (see Sect. 5.3 in Paper I). On the other hand, our approach
successfully recovers the true absorbed fraction without measurable bias for
all three samples.

Regarding the absorbed fraction vs $L_{\mathrm{X}}$, we need to reconstruct
the joint $L_{\mathrm{X}}$–$N_{\mathrm{H}}$ distribution, since the intrinsic
$L_{\mathrm{X}}$, which depends on most of the spectral parameters and on the
redshift, is unknown. Although the spectral method recovers relatively well
the true $L_{\mathrm{X}}$ (see Paper I for details), there could still be a
biased relation because of the bias in the classification of absorbed sources.
Therefore, a redistribution matrix that is able to take care not only of the
bias of $N_{\mathrm{H}}$, but also of $L_{\mathrm{X}}$, is required. The same
2-bin model of $N_{\mathrm{H}}$ is applied to derive the absorbed fraction,
and we use a non-uniform 8-bin model for the $L_{\mathrm{X}}$ distribution,
with broad bins at very low and very high $L_{\mathrm{X}}$ (2 bins below
$10^{43}$ erg s-1 and one bin above $10^{45}$ erg s-1), where few sources are
detected, and finer bins in between (5 equal-width bins).

We also test this method using three randomly selected sub-samples from the
simulation, as shown in Fig. 11, which have different absorbed fraction vs
$L_{\mathrm{X}}$ relations, but the same sample size and total
$L_{\mathrm{X}}$ distribution. Specifically, we impose a positive relation in
the first sub-sample, a negative one in the second, and the last one having a
peak of absorbed fraction at luminosity $10^{43.3}$ erg s-1. The results show
that, for the first sub-sample, the inferred fractions using the stacked PDF
or the medians are always biased low, while for the second one, the slope of
the relation is clearly biased towards flatter values. In the third case, the
two approaches also bias the fractions low, but the impact is stronger than in
the previous two cases, as the total relation becomes simply a flat one. In
particular, the recovered value around the peak at $\log L_{X}=43$ is more
than 2-$\sigma$ away from the true value. By contrast, ReMPI is able to
reconstruct the absorbed fraction more accurately with reasonable
uncertainties in the bins with $L_{\mathrm{X}}$ between $10^{42}$ to $10^{45}$
erg s-1 for all three sub-samples, and in the third situation, there is a
clear hint for a peak of the absorbed fraction at the right position, in spite
of the large errors in the two extreme bins due to the scarcity of sources in
the corresponding bins.

##  5 Analysis of the _XMM_ -COSMOS sample

The Bayesian spectral fitting method introduced in Paper I is applied to the
819 sources in our COSMOS data sample (see Sect. 2) to obtain the PDFs of the
spectral parameters of each source. We reconstruct here the parent
distributions of $N_{\mathrm{H}}$ and $\Gamma$, as well as the absorbed
fraction vs 2–10 keV $L_{\mathrm{X}}$ and redshift, using the methods
described in Sect. 4. Unless specified otherwise, the results are corrected
for the selection effects determined based on the simulations (see Sect. 3).

Figure 12: Inferred parent distributions of the photon index $\Gamma$ for
unabsorbed AGN (left) and absorbed AGN (right). Our results are shown in
black, with error bars representing 68.3% credible intervals. The results of
Ricci et al. (2017a) are plotted in red, and the blue dashed lines show those
from the parametric fitting, with the filled area showing also 68.3% credible
intervals. Figure 13: Observed distribution of $N_{\mathrm{H}}$ inferred with
our method (dashed orange line), and corrected for selection effects (solid
black line), with vertical error bars representing 68.3% credible intervals.
The blue and red points are the results from Ueda et al. (2014) and Ricci et
al. (2017a) respectively. The horizontal error bars show the bins. Figure 14:
Absorbed fraction vs 2–10 keV luminosity. Black points are the result of this
work, taking into account selection effects; the vertical error bars show
68.3% credible intervals and the horizontal error bars show the binning. The
blue points are the results of Ueda et al. (2014) for low-redshift sources,
and the orange area shows the results of Burlon et al. (2011). The grey area
is the observed luminosity distribution in our sample. Figure 15: Evolution of
the absorbed fraction with redshift. Our results corrected for selection
effects are shown in black. The results without selection correction are shown
in red with dashed error bars. Vertical error bars show 68.3% credible
intervals; horizontal error bars show the binning. The results of Hasinger et
al. (2007) are shown in blue. The grey area shows the redshift distribution of
our sample, including both spec-z and photo-z.

The inferred distributions of photon indices $\Gamma$ for absorbed and
unabsorbed AGN using ReMPI are shown in Fig. 12, as well as those using the
parametric approach. With ReMPI we find relatively non-Gaussian distributions
for both the absorbed and unabsorbed AGN in our sample, with a dominant peak
in the bin around $\Gamma=2$, containing around 50% of the total probability,
but the rest is distributed differently. The unabsorbed population has an
excess of probability to have a source at $\Gamma$ higher than the peak, but
very low probability below it, while the opposite is true for the absorbed
population. On the other hand, the standard approach obtains a much broader
Gaussian distribution for the unabsorbed population with $\mu=2.05\pm 0.01$
and $\sigma=0.24\pm 0.01$, and a relatively narrower one for the absorbed
population, with $\mu=1.85\pm 0.03$ and $\sigma=0.12\pm 0.02$. It also obtains
in general a harder distribution of $\Gamma$ for absorbed AGN, which is
consistent with the results of ReMPI.

The reconstructed $N_{\mathrm{H}}$ distribution is shown in Fig. 13. We can
see a clear bi-modality with two well separated peaks, one at $10^{20}$ cm-2
and the other at $10^{24}$ cm-2. Most of the unabsorbed sources have
$N_{\mathrm{H}}$ below $10^{20.7}$ cm-2, with a tentative decreasing trend
from $10^{20.2}$ to $10^{21.7}$ cm-2. Between $10^{20.7}$ and $10^{22.2}$
cm-2, there are very few sources and the distribution is compatible with zero,
while it increases steadily after entering the absorbed regime. The
uncertainties are relatively large for bins at both ends of the distribution.
For the unabsorbed ones the uncertainties are mainly the result from
unconstrained PDFs, as discussed in Sect. 4, while for the most absorbed ones
the uncertainties are due to the increasing effect of the selection correction
as absorption gets stronger. The bin with the least probability to have a
source is around $10^{21.5}$ cm-2, which is slightly below the adopted
boundary between absorbed and unabsorbed regimes, conventionally set at
$10^{22}$ cm-2, but agrees with that proposed by Merloni et al. (2014).

The reconstructed relation between the absorbed AGN fraction and the intrinsic
X-ray luminosity in the 2–10 keV band is illustrated in Fig. 14. We find the
well-known anti-correlation between the two in the range of $10^{43}$ –
$10^{45}$ erg s-1. The tentative peak is at $10^{43.2}$ erg s-1, where the
absorbed fraction is around 0.7. The highest luminosity bin ($10^{45}$ –
$10^{46}$ erg s-1) shows a small increase, compatible with no change in the
absorbed fraction considering its large uncertainty. This can be a statistical
fluctuation due to the paucity of very luminous sources in our sample, as
shown by the grey area in Fig. 14. For very faint AGN that are less luminous
than $10^{43.2}$ erg s-1, the behavior is very uncertain because of very few
sources and large selection correction.

The inferred absorbed fraction as a function of redshift is shown in Fig. 15.
The redshift of our data is assumed to be known because our sources have
either spec-z or good photo-z (Salvato et al. 2011; Laigle et al. 2016), whose
errors should mostly be smaller than the bin size, so that we can assign them
unambiguously into redshift bins. As a result, it seems that the absorbed
fraction could be a bit higher at redshift lower than 1, reaching around 0.6.
Beyond redshift 1, the absorbed AGN fraction is relatively stable around 0.5.

##  6 Discussion

###  6.1 Bayesian framework of X-ray spectral analysis and population
inference

We have developed a fully Bayesian framework to analyse a large sample of AGN
X-ray spectra, from the spectral fitting of each source to the inference of
parent distributions of the main spectral parameters, e.g. $N_{\mathrm{H}}$
and $\Gamma$. Since surveys mostly consist of low-S/N data, the key point of
the spectral fitting step is to properly propagate to the parameters of
interest the uncertainties associated with the presence of additional, ill-
constrained components like reflection and soft-excess, that are present in
the X-ray spectra of AGN, but not always detectable. We use a consistent
physical model that includes most of the physical components in AGN to fit the
data, irrespectively of the S/N, and use informative priors derived from deep
surveys to constrain secondary emission components (see Paper I for details).
In the end we marginalize over the nuisance parameters and obtain reliable
PDFs for the main parameters.

However, despite being more reliable, these PDFs can be loosely constrained
and present non-Gaussian shapes, and sometimes even show multi-modality, which
makes it very difficult to make correct population inference based on them.
The standard population inference method presented in Sect. 4.1 is able to
make a proper modelling and take into account the full PDFs, but it assumes
they always provide the true posterior probabilities, i.e., they are unbiased.
We find clear cases this assumption is violated. In particular, the PDFs of
$N_{\mathrm{H}}$ for unabsorbed objects are biased, because the likelihood may
become insensitive to the value of the parameter if it is too low, resulting
in a flat posterior, limited only by the prior. The bias is a priori unknown,
and therefore cannot be easily corrected. In the new method (ReMPI) we propose
here, we model the stacked PDFs by summing the PDFs that are expected from the
parent distribution. These expected PDFs are determined using a transfer
function – called “response” in analogy with the standard way of analyzing
X-ray spectra – that takes into account the biases resulting from the fitting
process. The response is constructed by simulating a large realistic sample of
sources that are added onto the real image data and extracted with the same
source detection tool and the same background as those of real sources, such
that the simulated sources have the same selection function as the real
sources (see Sect. 3 for details).

To reduce the dimensionality of the problem, we construct the response, in the
form of a redistribution matrix, for a specific application by marginalizing
over the other parameters. One may have the concern that the simulation is
based on the assumptions we make on the distributions of these parameters,
thus the redistribution matrix may be not representative of the real data.
However, just like the priors we use in the Bayesian analysis, these
assumptions are also priors that reflect our state of knowledge of this
problem, which are consistent with those we use in spectral fitting. In Sect.
4.3, the relations between the absorbed fraction and redshift and
$L_{\mathrm{X}}$ are successfully recovered for different sub-samples with
different input relations, using the same redistribution matrix computed from
the total simulated sample. This demonstrates that our method is not extremely
sensitive to the simulated sample we use to calibrate the redistribution
matrix. Apart from the priors of the spectral parameters, another prior we put
in the simulation is the X-ray luminosity function, which is essential to
simulate a sample with realistic $L_{\mathrm{X}}$ and redshift distribution.
However, the parameterization we use does not include the dependence on
spectral parameters like $N_{\mathrm{H}}$, hence any difference we find
between absorbed and unabsorbed objects cannot be due to the prior. The
redistribution matrix also takes into account the statistical uncertainty in
the output posterior distribution, so that the method naturally deconvolves
the stacked PDFs from this effect.

The performance of ReMPI in recovering the parent distributions of $\Gamma$,
$N_{\mathrm{H}}$, as well as inferring the redshift and the intrinsic X-ray
Luminosity dependence of the absorbed fraction is tested on a realistic
simulated data set. We see that using point estimators like the median, or
with the simple stacked PDF one cannot recover the true properties of the
parent population. For instance, the medians of the PDFs of $N_{\mathrm{H}}$
indicate a false peak in the unabsorbed regime and another one in the absorbed
regime, which is biased towards stronger absorption, whereas the stacked PDF
hardly recovers the bi-modality and severely underestimates the peak of
absorbed sources (see Fig. 9). However, unlike these two naive approaches, the
parametric fitting approach introduced in Sect. 4.1 is able to recover the
true population, as long as the PDFs are unbiased, e.g., in the case of
$\Gamma$ (see Fig. 5), provided we have a reasonable analytic expression for
the distribution, which is not always available.

###  6.2 Population properties of the AGN X-ray spectral parameters

In this subsection we discuss our results of applying ReMPI to _XMM_ -COSMOS
and compare them with previous studies. We point out however that most other
AGN X-ray population studies make strong assumptions about some parameters,
for instance, no soft excess or no reflection, which is equivalent to using
extremely strong but inappropriate priors, so the comparison with their
results must be taken with some caution.

####  6.2.1 Photon index $\Gamma$

We reconstruct the distributions of $\Gamma$ for absorbed and unabsorbed AGN
separately for two reasons. First, they are found to be different in some
previous works, e.g., Burlon et al. (2011); Ueda et al. (2014), and second, it
is a natural outcome from our method, because we take into account
$N_{\mathrm{H}}$ in the redistribution matrix of $\Gamma$. The two
distributions are simply marginalized over different $N_{\mathrm{H}}$ ranges,
adopting a flat prior on $N_{\mathrm{H}}$.

The results we find using ReMPI from our _XMM_ -COSMOS sample show that the
majority of AGN in both absorbed and unabsorbed populations have $\Gamma$
between 1.9 and 2.1, as shown in Fig. 12, which is about 0.2 higher than the
typical values found in Ricci et al. (2017a) for Swift/BAT AGN, with
$\Gamma=1.80\pm 0.02$ for unabsorbed AGN and $1.76\pm 0.02$ for absorbed AGN.
The reason why we find on average slightly higher $\Gamma$ values could be
that we apply a consistent model, which includes reprocessed emission like
reflection, for all the sources, unlike in Ricci et al. (2017a), where
secondary components are added on a case-by-case basis and many objects have
only an upper limit to the relative reflection. Reflection flattens the slope
of the spectrum ,because it adds photons in the hard X-ray band. Namely, there
is a degeneracy between a single hard power law and a softer one plus
reflection. Quantitatively, if we fit a redshifted cutoff power law model
(zcutoffpl) to a faked XMM-Newton spectrum of a power law plus reflection
(zcutoffpl \+ pexmon) at redshift 0.5 with the relative-reflection parameter
$R=1$ (-1, to use the xspec convention to ignore the direct emission in
pexmon), we will get a fit having a harder photon index by about 0.1. Since we
take into account the possibility of the presence of reflection by applying a
flat prior between 0 and 1.5 on $R$, which is little constrained in most of
our sources, marginalizing over the posterior of $R$ results in larger values
of $\Gamma$. On the other hand, Ueda et al. (2014) fixed $R$ to 1 in all
cases, and found that the unabsorbed AGN have an average photon index
$\Gamma=1.94\pm 0.09$, and the absorbed AGN have $\Gamma=1.84\pm 0.15$, which
are in good agreement with our results, especially those using the parametric
fitting.

Alternatively, the discrepancy between the average photon index estimated here
and that determined by Ricci et al. (2017a) may arise from uncertainties in
the effective area calibration of the various X-ray instruments. Indeed, the
calibration of _XMM-Newton_ and _Chandra_ effective areas is known to be
inconsistent at the 10–15% level. When comparing the fitted temperature of
galaxy clusters, Schellenberger et al. (2015) showed that for high
temperatures the calibration uncertainties can be as large as 20%, with
_Chandra_ returning systematically higher temperatures than _XMM-Newton_. To
translate the uncertainty into a photon index, we simulated EPIC-pn
Bremsstrahlung spectra with temperatures of 8 keV and 10 keV, respectively,
and fitted the simulated spectra with a power law. We found that the
corresponding uncertainty in the photon index is about 0.2, which may explain
the systematic difference observed here. In the majority of cases, Ricci et
al. (2017a) used soft X-ray spectra from _Swift_ /XRT. Using a simultaneous
observation campaign of the quasar 3C 273 with several X-ray observatories,
Madsen et al. (2017) showed that _Swift_ /XRT returns a photon index of 1.5
compared to 1.7 for _XMM-Newton_ /EPIC. The harder photon indices measured by
Ricci et al. (2017a) compared to our work may thus be explained by a
difference in the effective area calibration of _XMM-Newton_ /EPIC with
respect to _Swift_ /XRT.

Another difference between our results and those of Ricci et al. (2017a) is
that we have more constrained distributions for both absorbed and unabsorbed
AGN, while the $\Gamma$ distributions in Ricci et al. (2017a) are broader,
especially for the absorbed AGN. This larger scatter can be explained by the
statistical uncertainties on the photon index measurement of each spectrum,
which are generally of the order of 0.2–0.3, giving rise to the broadening of
the $\Gamma$ distribution. On the other hand, with ReMPI, the redistribution
matrix not only removes the bias of the PDFs, but also the scatter induced by
measurement uncertainties, so that the resulting parent distribution is
minimally affected by the statistical uncertainties of each measurement. The
parametric fitting approach is also able to handle this statistical
uncertainty, but not the bias in the PDFs as discussed above.

Although our results for the $\Gamma$ distributions of absorbed and unabsorbed
objects are comparable to each other to the first order due to the prominent
peak at $\Gamma=2.0$, there is also evidence that the unabsorbed population
contains more soft objects than the absorbed one, while the latter has more
hard objects (see Fig. 12). Using the parametric fitting we obtain two
Gaussian distributions with $\Gamma=2.05\pm 0.24$ for the unabsorbed sample
and $\Gamma=1.85\pm 0.12$ for the absorbed one, showing again that unabsorbed
AGN have on average softer photon indices, which is in agreement with Burlon
et al. (2011) and Ueda et al. (2014). The difference might be a consequence of
a relationship between Eddington rate ($\lambda_{\mathrm{Edd}}$) and
obscuration (Ricci et al. 2017b), combined with the extensively studied
positive correlation between $\Gamma$ and $\lambda_{\mathrm{Edd}}$ (Brandt et
al. 1997; Lu & Yu 1999; Shemmer et al. 2006, 2008; Risaliti et al. 2009;
Brightman et al. 2013; Kawamuro et al. 2016; Ricci et al. 2018), which could
be explained by a more efficient cooling of the X-ray emitting corona at
higher $\lambda_{\mathrm{Edd}}$ (Vasudevan & Fabian 2007; Davis & Laor 2011).

####  6.2.2 Absorbing column density $N_{\mathrm{H}}$

We do not attempt to perform a parametric inference for the $N_{\mathrm{H}}$
distribution, but use only our ReMPI forward-fitting inference approach,
because it is clear from Paper I that the PDFs of $N_{\mathrm{H}}$ are
strongly biased. In addition, it is not clear which functional form should be
used to describe the distribution of $N_{\mathrm{H}}$.

The reconstructed $N_{\mathrm{H}}$ distribution of our _XMM_ -COSMOS sample
shown in Fig. 13 is in good agreement with the results from Ueda et al. (2014)
and Ricci et al. (2017a), both of which use local AGN samples detected by
Swift/BAT in the 14–195 keV energy range, which does not suffer too strongly
from the incompleteness problem at high column densities. Ueda et al. (2014)
utilize the Swift/BAT 9-month sample (Tueller et al. 2008), which contains 137
non-blazar AGN, while Ricci et al. (2017a) adopt the extended 70-month sample
with 733 non-blazar AGN (Baumgartner et al. 2013), and both of them fit the
parameters using complementary observations below 10 keV.

Besides the good match, we can identify very clearly a bimodality in the
$N_{\mathrm{H}}$ distribution that is only suggested by Ueda et al. (2014) and
Ricci et al. (2017a). The two peaks are widely separated, with one peak in the
unabsorbed regime around $10^{20}$ cm-2, and the other in the absorbed regime
around $10^{23}$ – $10^{24}$ cm-2. Even though the uncertainties on the two
peaks are relatively large, the decreasing trend of the $N_{\mathrm{H}}$
distribution of unabsorbed sources from $10^{20}$ to $10^{21.5}$ cm-2 and the
increasing trend from $10^{21.5}$ to $10^{23}$ cm-2 are very clear. Above
$10^{23}$ cm-2, the correction of selection effects starts to be dominant and
even shifts the peak from the original distribution, at $10^{23}$ cm-2, to an
order of magnitude higher. In view of the large uncertainties, the location of
the peak cannot be precisely determined, but it is very likely to be close to
$10^{24}$ cm-2. In addition, our results indicate that there is practically no
source with $N_{\mathrm{H}}$ around $10^{21.5\pm 0.25}$ cm-2, which is in good
agreement with the results of Merloni et al. (2014), where a boundary at
$N_{\mathrm{H}}=10^{21.5}$ cm-2 is suggested to classify absorbed and
unabsorbed AGN. This clear gap is a strong indication that the absorbing
materials of these two populations have different origins (Paltani et al.
2008). The one peaked at low $N_{\mathrm{H}}$ is likely due to the absorption
by the host galaxy itself, while the other one is caused by the denser torus
around the nucleus.

####  6.2.3 Absorbed AGN fraction vs intrinsic X-ray luminosity

The absorbed AGN fraction is an important indicator in the study of AGN
populations, because we can use it to infer the geometry and characteristics
of the absorbing material. For example, the simple unification model of AGN
(Antonucci 1993), in which the amount of obscuration is only an effect of
inclination in a homogeneous population, predicts that the absorbed AGN
fraction should be independent of all other physical parameters. However, we
find a clear decrease of the absorbed fraction with increasing intrinsic 2–10
keV X-ray luminosity from around $10^{43}$ to $10^{45}$ erg s-1, so this
simple model can be ruled out. This dependency was found in previous studies
(Hasinger 2008; Burlon et al. 2011; Ueda et al. 2014; Buchner et al. 2015);
the results of Burlon et al. (2011) and Ueda et al. (2014) are shown in Fig.
14, which are in broad agreement with our results considering the
uncertainties. Similarly, Ricci et al. (2017a) also found that unabsorbed AGN
have typically higher intrinsic luminosities than absorbed ones. Ricci et al.
(2017b) further demonstrated that the dependence disappears when the sources
are separated into bins of $\lambda_{\mathrm{Edd}}$, thus the main driver of
the decreasing absorbed fraction may be $\lambda_{\mathrm{Edd}}$ instead of
$L_{\mathrm{X}}$.

On the high-$L_{\mathrm{X}}$ end, our results are slightly above those of Ueda
et al. (2014) around $10^{44.5}$ to $10^{45}$ erg s-1, and the decreasing
trend may flatten above $10^{45}$ erg s-1. However, this is very uncertain
because of the scarcity of objects. Similarly, on the low-$L_{\mathrm{X}}$
end, because of the large uncertainties of the first two $L_{\mathrm{X}}$
bins, our result is still compatible with all published studies, including
Brightman & Nandra (2011), Burlon et al. (2011) and Buchner et al. (2015),
which found a positive relation when $L_{\mathrm{X}}$ is below $10^{43}$ erg
s-1. As shown in Fig. 14, our result between $10^{42}$ to $10^{43.4}$ erg s-1
match well with those of Burlon et al. (2011), especially regarding the
location of the peak. More data are needed to better constrain the dependence
of the absorbed fraction on $L_{\mathrm{X}}$ below $10^{43}$ erg s-1, which
requires surveys of larger area, e.g., the XXL survey (Pierre et al. 2016),
and/or deeper exposure (Giacconi et al. 2002).

####  6.2.4 Absorbed AGN fraction vs redshift

To study the evolution of the absorbed AGN fraction with redshift, we
reconstruct the $N_{\mathrm{H}}$ distribution in each redshift bin, and then
calculate the respective absorbed fractions. As shown in Fig. 15, to first
order we see a flat distribution of the absorbed fraction as a function of
redshift for $z\geq 1$, with a slight increase at lower redshifts. We compare
our results with those of Hasinger (2008), which were obtained using a
combined sample of 10 independent hard X-ray samples selected in the 2–10 keV
band. Their total sample consists of 1290 sources, 1106 of which have reliable
optical spectroscopic redshifts. Since their results are not corrected for
selection effects, we plot also ours without selection correction in Fig. 15
for comparison. The results indicate that the overall absorbed fraction from
Hasinger (2008) is underestimated because of the selection effects, and even
when selecting in the 2–10 keV range our selection is slightly biased towards
unabsorbed sources. Although the selection functions are different in the two
data samples, both results are in relatively good agreement with each other.
We measure however overall a slightly lower absorbed fraction at $z>1$ and
higher for low redshift ones, but the difference is within $1\sigma$.

As discussed in Hasinger (2008), since in each redshift bin the absorbed AGN
fraction is integrated over luminosity, its correlation with redshift may be
cancelled out by the anti-correlation with $L_{\mathrm{X}}$. It would be
better to reconstruct the $N_{\mathrm{H}}$ distribution in the joint redshift
and luminosity space, in order to break the degeneracies. Following a similar
approach, Hasinger (2008) obtained a positive correlation between the absorbed
fraction and redshift, which was also found in Ueda et al. (2014) and Buchner
et al. (2015). However, in this work the limited size of our COSMOS sample
makes it unfeasible to perform an accurate three-dimensional reconstruction of
the parameter space using ReMPI, since the number of sources per 3-dimensional
bin will be too small to get meaningful measurements. Larger surveys are
required to study the evolution of the absorbed AGN fraction with redshift.

##  7 Summary and conclusions

In this paper we developed ReMPI, a novel non-parametric Bayesian method of
population inference that we applied to the reconstruction of the parent
distributions of AGN X-ray spectral parameters of interest, e.g., $\Gamma$,
$N_{\mathrm{H}}$, and $L_{\mathrm{X}}$, based on the PDFs of the parameters of
each individual spectrum in a large survey where most sources have low to
moderate S/N.

We tested both a parametric fitting approach and ReMPI on simulated data. The
parametric approach works well when a parametric form of the parent
distribution exists, provided that the PDFs are unbiased (e.g., $\Gamma$, see
Sect. 4.1). As some parameters have highly biased PDFs (e.g.,
$N_{\mathrm{H}}$, see Sect. 4.2), we use ReMPI to correct the systematic bias
in the PDFs by applying a redistribution matrix on the stacked PDF of the
total sample that models the bias introduced by the spectral-fitting process,
and by forward-fitting the non-parametric histogram model to the redistributed
stacked PDF. The redistribution matrix $M_{\mathrm{red}}$ is computed from a
large realistic simulation that reproduces the redshift and S/N distributions
of our sample, on which $M_{\mathrm{red}}$ is strongly dependent. The sample
was simulated using the same priors as those of the spectral parameters
applied in the spectral fitting step, so that the other parameters are
marginalized over consistently. The tests on the sub-samples of the
simulation, which have different parent distributions of $N_{\mathrm{H}}$,
$\Gamma$ and so on, show that ReMPI is able to reconstruct the parent
distributions with small biases and reasonable uncertainties, even when using
not well-constrained and sometimes highly biased PDFs (Sect. 4.2).

As a pilot study for our new Bayesian methods, we applied the spectral
analysis method presented in Paper I to the raw XMM-Newton observations of the
COSMOS field, and reconstructed the parent distributions of the main
parameters using ReMPI, as well as the parametric approach for the photon
index $\Gamma$. We found a clear bimodality in the $N_{\mathrm{H}}$
distribution, with widely separated peaks, and an almost total lack of objects
with $N_{\mathrm{H}}$ around $10^{21.5}$ cm-2. The $\Gamma$ distributions of
absorbed and unabsorbed AGN are similar, with the common feature that the
majority of objects have spectral indices close to 2, which is around 0.1–0.2
higher than the typical values found for Swift/BAT AGN Ricci et al. (2017a).
The discrepancy can be possibly explained by calibration issues or by the
difference in the spectral modelling. However, the asymmetries in the
distributions for absorbed and unabsorbed objects indicate that there are more
soft ($\Gamma>2$) objects in the unabsorbed sample than in the absorbed one,
and conversely more hard ($\Gamma<2$) objects in the absorbed sample than in
the unabsorbed one. We also recovered the well-known anti-correlation between
the absorbed AGN fraction and the intrinsic 2–10 keV $L_{\mathrm{X}}$ when
$L_{\mathrm{X}}$ is higher than about $10^{43}$ erg s-1. Finally, we did not
find a clear evolution of the absorbed fraction with redshift, as the absorbed
fraction is mostly flat and close to 0.5, with a small excess below $z=1$.
However, more sources are needed to disentangle the luminosity dependence in
each redshift bin.

With the application of our set of Bayesian methods on the _XMM_ -COSMOS
survey, we demonstrated that it is possible to reconstruct the distributions
of the main AGN X-ray spectral parameters that are essentially consistent with
previous results, albeit in an automatic way and using a single physically
motivated spectral model, and a single sample limited to 819 sources with
mostly low S/N. Our method can be readily applied to larger surveys, e.g., the
XXL (Pierre et al. 2016) and eRASS (Merloni et al. 2012) surveys, which will
allow us to perform a detailed study of the parameters of interest and of
their interplay, which is our next step of the application of our methods.

###### Acknowledgements.

L.G. acknowledges partial support from the Swiss National Science Foundation.

## References

  * Aird et al. (2015) Aird, J., Coil, A. L., Georgakakis, A., et al. 2015, MNRAS, 451, 1892
  * Antonucci (1993) Antonucci, R. 1993, ARA&A, 31, 473
  * Arnaud (1996) Arnaud, K. A. 1996, in Astronomical Society of the Pacific Conference Series, Vol. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes, 17
  * Arnaud et al. (1985) Arnaud, K. A., Branduardi-Raymont, G., Culhane, J. L., et al. 1985, MNRAS, 217, 105
  * Avni & Tananbaum (1986) Avni, Y. & Tananbaum, H. 1986, ApJ, 305, 83
  * Baumgartner et al. (2013) Baumgartner, W. H., Tueller, J., Markwardt, C. B., et al. 2013, ApJS, 207, 19
  * Bianchi et al. (2006) Bianchi, S., Guainazzi, M., & Chiaberge, M. 2006, A&A, 448, 499
  * Blandford et al. (1990) Blandford, R. D., Netzer, H., Woltjer, L., Courvoisier, T. J. L., & Mayor, M. 1990, Active Galactic Nuclei (Springer)
  * Brandt & Alexander (2015) Brandt, W. N. & Alexander, D. M. 2015, A&A Rev., 23, 1
  * Brandt et al. (2000) Brandt, W. N., Laor, A., & Wills, B. J. 2000, ApJ, 528, 637
  * Brandt et al. (1997) Brandt, W. N., Mathur, S., & Elvis, M. 1997, MNRAS, 285, L25
  * Brightman & Nandra (2011) Brightman, M. & Nandra, K. 2011, MNRAS, 414, 3084
  * Brightman et al. (2013) Brightman, M., Silverman, J. D., Mainieri, V., et al. 2013, MNRAS, 433, 2485
  * Brusa et al. (2010) Brusa, M., Civano, F., Comastri, A., et al. 2010, ApJ, 716, 348
  * Buchner et al. (2015) Buchner, J., Georgakakis, A., Nandra, K., et al. 2015, ApJ, 802, 89
  * Buchner et al. (2014) Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125
  * Burlon et al. (2011) Burlon, D., Ajello, M., Greiner, J., et al. 2011, ApJ, 728, 58
  * Cappelluti et al. (2009) Cappelluti, N., Brusa, M., Hasinger, G., et al. 2009, A&A, 497, 635
  * Civano et al. (2016) Civano, F., Marchesi, S., Comastri, A., et al. 2016, ApJ, 819, 62
  * Davis & Laor (2011) Davis, S. W. & Laor, A. 2011, ApJ, 728, 98
  * Donley et al. (2008) Donley, J. L., Rieke, G. H., Pérez-González, P. G., & Barro, G. 2008, ApJ, 687, 111
  * Elvis et al. (2009) Elvis, M., Civano, F., Vignali, C., et al. 2009, ApJS, 184, 158
  * Feroz & Hobson (2008) Feroz, F. & Hobson, M. P. 2008, MNRAS, 384, 449
  * Feroz et al. (2009) Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601
  * Feroz et al. (2019) Feroz, F., Hobson, M. P., Cameron, E., & Pettitt, A. N. 2019, The Open Journal of Astrophysics, 2, 10
  * Ferrarese & Merritt (2000) Ferrarese, L. & Merritt, D. 2000, ApJ, 539, L9
  * Foreman-Mackey (2016) Foreman-Mackey, D. 2016, The Journal of Open Source Software, 1, 24
  * Freeman et al. (2001) Freeman, P., Doe, S., & Siemiginowska, A. 2001, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 4477, Astronomical Data Analysis, ed. J.-L. Starck & F. D. Murtagh, 76–87
  * Ge et al. (2021) Ge, L., Paltani, S., & Eckert, D. 2021, submitted to A&A
  * Ghirardini et al. (2019) Ghirardini, V., Eckert, D., Ettori, S., et al. 2019, A&A, 621, A41
  * Giacconi et al. (2002) Giacconi, R., Zirm, A., Wang, J., et al. 2002, ApJS, 139, 369
  * Gibson et al. (2008) Gibson, R. R., Brandt, W. N., & Schneider, D. P. 2008, ApJ, 685, 773
  * Gozaliasl et al. (2019) Gozaliasl, G., Finoguenov, A., Tanaka, M., et al. 2019, MNRAS, 483, 3545
  * Guainazzi & Bianchi (2007) Guainazzi, M. & Bianchi, S. 2007, MNRAS, 374, 1290
  * Gültekin et al. (2009) Gültekin, K., Richstone, D. O., Gebhardt, K., et al. 2009, ApJ, 698, 198
  * Häring & Rix (2004) Häring, N. & Rix, H.-W. 2004, ApJ, 604, L89
  * Hasinger (2008) Hasinger, G. 2008, A&A, 490, 905
  * Hasinger et al. (2007) Hasinger, G., Cappelluti, N., Brunner, H., et al. 2007, ApJS, 172, 29
  * Hickox & Alexander (2018) Hickox, R. C. & Alexander, D. M. 2018, ARA&A, 56, 625
  * Hinshaw et al. (2013) Hinshaw, G., Larson, D., Komatsu, E., et al. 2013, ApJS, 208, 19
  * Hopkins et al. (2006) Hopkins, P. F., Somerville, R. S., Hernquist, L., et al. 2006, ApJ, 652, 864
  * Iwasawa et al. (2011) Iwasawa, K., Sanders, D. B., Teng, S. H., et al. 2011, A&A, 529, A106
  * Jin et al. (2012) Jin, C., Ward, M., Done, C., & Gelbord, J. 2012, MNRAS, 420, 1825
  * Kaastra et al. (1996) Kaastra, J. S., Mewe, R., & Nieuwenhuijzen, H. 1996, in UV and X-ray Spectroscopy of Astrophysical and Laboratory Plasmas, 411–414
  * Kawamuro et al. (2016) Kawamuro, T., Ueda, Y., Tazaki, F., Ricci, C., & Terashima, Y. 2016, ApJS, 225, 14
  * Kerr (1963) Kerr, R. P. 1963, Phys. Rev. Lett., 11, 237
  * Kormendy & Ho (2013) Kormendy, J. & Ho, L. C. 2013, ARA&A, 51, 511
  * Krolik (1999) Krolik, J. H. 1999, Active galactic nuclei : from the central black hole to the galactic environment (American Journal of Physics)
  * Krolik & Begelman (1988) Krolik, J. H. & Begelman, M. C. 1988, ApJ, 329, 702
  * Laigle et al. (2016) Laigle, C., McCracken, H. J., Ilbert, O., et al. 2016, ApJS, 224, 24
  * Loredo (2004) Loredo, T. J. 2004, in American Institute of Physics Conference Series, Vol. 735, Bayesian Inference and Maximum Entropy Methods in Science and Engineering: 24th International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering, ed. R. Fischer, R. Preuss, & U. V. Toussaint, 195–206
  * Lu & Yu (1999) Lu, Y. & Yu, Q. 1999, ApJ, 526, L5
  * Luo et al. (2017) Luo, B., Brandt, W. N., Xue, Y. Q., et al. 2017, ApJS, 228, 2
  * Maccacaro et al. (1988) Maccacaro, T., Gioia, I. M., Wolter, A., Zamorani, G., & Stocke, J. T. 1988, ApJ, 326, 680
  * Madsen et al. (2017) Madsen, K. K., Beardmore, A. P., Forster, K., et al. 2017, AJ, 153, 2
  * Magorrian et al. (1998) Magorrian, J., Tremaine, S., Richstone, D., et al. 1998, AJ, 115, 2285
  * Malz (2021) Malz, A. I. 2021, Phys. Rev. D, 103, 083502
  * Marchesi et al. (2016) Marchesi, S., Civano, F., Elvis, M., et al. 2016, ApJ, 817, 34
  * Marconi & Hunt (2003) Marconi, A. & Hunt, L. K. 2003, ApJ, 589, L21
  * Marconi et al. (2004) Marconi, A., Risaliti, G., Gilli, R., et al. 2004, MNRAS, 351, 169
  * Merloni et al. (2014) Merloni, A., Bongiorno, A., Brusa, M., et al. 2014, MNRAS, 437, 3550
  * Merloni & Heinz (2008) Merloni, A. & Heinz, S. 2008, MNRAS, 388, 1011
  * Merloni et al. (2012) Merloni, A., Predehl, P., Becker, W., et al. 2012, arXiv e-prints, arXiv:1209.3114
  * Nandra & Pounds (1994) Nandra, K. & Pounds, K. A. 1994, MNRAS, 268, 405
  * Paltani et al. (2008) Paltani, S., Walter, R., McHardy, I. M., et al. 2008, A&A, 485, 707
  * Pierre et al. (2016) Pierre, M., Pacaud, F., Adami, C., et al. 2016, A&A, 592, A1
  * Pounds et al. (1995) Pounds, K. A., Done, C., & Osborne, J. P. 1995, MNRAS, 277, L5
  * Pounds et al. (1990) Pounds, K. A., Nandra, K., Stewart, G. C., George, I. M., & Fabian, A. C. 1990, Nature, 344, 132
  * Ricci et al. (2018) Ricci, C., Ho, L. C., Fabian, A. C., et al. 2018, MNRAS, 480, 1819
  * Ricci et al. (2017a) Ricci, C., Trakhtenbrot, B., Koss, M. J., et al. 2017a, ApJS, 233, 17
  * Ricci et al. (2017b) Ricci, C., Trakhtenbrot, B., Koss, M. J., et al. 2017b, Nature, 549, 488
  * Risaliti et al. (2009) Risaliti, G., Young, M., & Elvis, M. 2009, ApJ, 700, L6
  * Salvatier et al. (2016) Salvatier, J., Wieckiâ, T. V., & Fonnesbeck, C. 2016, PyMC3: Python probabilistic programming framework
  * Salvato et al. (2009) Salvato, M., Hasinger, G., Ilbert, O., et al. 2009, ApJ, 690, 1250
  * Salvato et al. (2011) Salvato, M., Ilbert, O., Hasinger, G., et al. 2011, ApJ, 742, 61
  * Sanders et al. (1989) Sanders, D. B., Phinney, E. S., Neugebauer, G., Soifer, B. T., & Matthews, K. 1989, ApJ, 347, 29
  * Schellenberger et al. (2015) Schellenberger, G., Reiprich, T. H., Lovisari, L., Nevalainen, J., & David, L. 2015, A&A, 575, A30
  * Scoville et al. (2007) Scoville, N., Aussel, H., Brusa, M., et al. 2007, ApJS, 172, 1
  * Shakura & Sunyaev (1973) Shakura, N. I. & Sunyaev, R. A. 1973, A&A, 500, 33
  * Shapiro & Teukolsky (1983) Shapiro, S. L. & Teukolsky, S. A. 1983, Black holes, white dwarfs, and neutron stars : the physics of compact objects (Wiley-VCH)
  * Shemmer et al. (2006) Shemmer, O., Brandt, W. N., Netzer, H., Maiolino, R., & Kaspi, S. 2006, ApJ, 646, L29
  * Shemmer et al. (2008) Shemmer, O., Brandt, W. N., Netzer, H., Maiolino, R., & Kaspi, S. 2008, ApJ, 682, 81
  * Silverman et al. (2008) Silverman, J. D., Green, P. J., Barkhouse, W. A., et al. 2008, ApJ, 679, 118
  * Singh et al. (1985) Singh, K. P., Garmire, G. P., & Nousek, J. 1985, ApJ, 297, 633
  * Skilling (2004) Skilling, J. 2004, in American Institute of Physics Conference Series, Vol. 735, American Institute of Physics Conference Series, ed. R. Fischer, R. Preuss, & U. V. Toussaint, 395–405
  * Tueller et al. (2008) Tueller, J., Mushotzky, R. F., Barthelmy, S., et al. 2008, ApJ, 681, 113
  * Ueda et al. (2014) Ueda, Y., Akiyama, M., Hasinger, G., Miyaji, T., & Watson, M. G. 2014, ApJ, 786, 104
  * Ueda et al. (2003) Ueda, Y., Akiyama, M., Ohta, K., & Miyaji, T. 2003, ApJ, 598, 886
  * Urry & Padovani (1995) Urry, C. M. & Padovani, P. 1995, PASP, 107, 803
  * Vasudevan & Fabian (2007) Vasudevan, R. V. & Fabian, A. C. 2007, MNRAS, 381, 1235
  * Vito et al. (2017) Vito, F., Brandt, W. N., Yang, G., et al. 2017, Monthly Notices of the Royal Astronomical Society, 473, 2378
  * Weaver et al. (2021) Weaver, J. R., Kauffmann, O. B., Ilbert, O., et al. 2021, arXiv e-prints, arXiv:2110.13923
  * Wilms et al. (2000) Wilms, J., Allen, A., & McCray, R. 2000, ApJ, 542, 914
  * Zdziarski et al. (1995) Zdziarski, A. A., Johnson, W. N., Done, C., Smith, D., & McNaron-Brown, K. 1995, ApJ, 438, L63
  * Zdziarski et al. (1996) Zdziarski, A. A., Johnson, W. N., & Magdziarz, P. 1996, MNRAS, 283, 193

Table 1: List of observations. ObsId | COSMOS field # | RA | Dec
---|---|---|---
0203360101 | 1 | 150.61005 | 2.71
0203360201 | 2 | 150.61005 | 2.46
0203360301 | 3 | 150.61005 | 2.21
0203360401 | 4 | 150.61005 | 1.96
0203360501 | 5 | 150.61005 | 1.71
0203360601 | 6 | 150.36000 | 2.71
0203360701 | 7 | 150.36000 | 2.46
0203360801 | 8 | 150.36000 | 2.21
0203360901 | 9 | 150.36000 | 1.96
0203361001 | 10 | 150.36000 | 1.71
0203361101 | 11 | 150.10995 | 2.71
0203361201 | 12 | 150.10995 | 2.46
0203361301 | 13 | 150.10995 | 2.21
0203361401 | 14 | 150.10995 | 1.96
0203361501 | 15 | 150.10995 | 1.71
0203361601 | 16 | 149.85999 | 2.71
0203361701 | 17 | 149.85999 | 2.46
0203361801 | 18 | 149.85999 | 2.21
0203361901 | 19 | 149.85999 | 1.96
0203362001 | 20 | 149.85999 | 1.71
0203362101 | 21 | 149.61000 | 2.71
0203362201 | 22 | 149.61000 | 2.46
0203362301 | 23 | 149.61000 | 2.21
0203362401 | 24 | 149.61000 | 1.96
0203362501 | 25 | 149.61000 | 1.71
0302350101 | B1 | 150.61005 | 2.73
0302350201 | B2 | 150.62670 | 2.46
0302350301 | B3 | 150.61005 | 2.23
0302350401 | B4 | 150.62670 | 1.96
0302350501 | B5 | 150.61005 | 1.73
0302350601 | B6 | 150.34335 | 2.71
0302350701 | B7 | 150.36000 | 2.48
0302350801 | B8 | 150.34335 | 2.21
0302350901 | B9 | 150.36000 | 1.98
0302351001 | B10 | 150.34335 | 1.71
0302351101 | B11 | 150.10995 | 2.73
0302351201 | B12 | 150.12660 | 2.46
0302351301 | B13 | 150.10995 | 2.23
0302351401 | B14 | 150.12660 | 1.96
0302351501 | B15 | 150.10995 | 1.73
0302351601 | B16 | 149.84333 | 2.71
0302351701 | B17 | 149.85999 | 2.48
0302351801 | B18 | 149.84333 | 2.21
0302351901 | B19 | 149.85999 | 1.98
0302352001 | B20 | 149.84333 | 1.71
0302352201 | B22 | 149.62667 | 2.46
0302352301 | B23 | 149.61000 | 2.23
0302352401 | B24 | 149.62667 | 1.96
0302352501 | B25 | 149.61000 | 1.73
0302353001 | B9 | 150.36000 | 1.98
0302353101 | B3 | 150.61005 | 2.23
0302353201 | B24 | 149.62667 | 1.96
0302353301 | B25 | 149.61000 | 1.73
0302353401 | B14 | 150.12660 | 1.96
0501170101 | 20C | 149.84333 | 1.71
0501170201 | 23C | 149.61000 | 2.23

##  Appendix A Observation list