# A High Fraction of Heavily X-ray Obscured Active Galactic Nuclei

Christopher M. Carroll Department of Physics and Astronomy, Washington State
University, 1245 Webster Hall, Pullman, WA 99164, USA  Department of Physics
and Astronomy, Dartmouth College, 6127 Wilder Laboratory, Hanover, NH 03755,
USA  Tonima T. Ananna Department of Physics and Astronomy, Dartmouth College,
6127 Wilder Laboratory, Hanover, NH 03755, USA  Ryan C. Hickox Department of
Physics and Astronomy, Dartmouth College, 6127 Wilder Laboratory, Hanover, NH
03755, USA  Alberto Masini SISSA, Via Bonomea 265, 34151 Trieste, Italy  INAF
– Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, via Gobetti
93/3, I-40129 Bologna, Italy  Roberto J. Assef Núcleo de Astronomía de la
Facultad de Ingeniería y Ciencias, Universidad Diego Portales, Av. Ejército
Libertador 441, Santiago, Chile  Daniel Stern Jet Propulsion Laboratory,
California Institute of Technology, 4800 Oak Grove Drive, Mail Stop 169-221,
Pasadena, CA 91109, USA  Chien-Ting J. Chen Astrophysics Office, NASA Marshall
Space Flight Center, ZP12, Huntsville, AL 35812, USA  Lauranne Lanz Department
of Physics, The College of New Jersey, 2000 Pennington Road, Ewing, NJ 08628,
USA

(Received September 26, 2021; Revised May 6, 2022; March 3, 2023; Accepted
March 9, 2023)

###### Abstract

We present new estimates on the fraction of heavily X-ray obscured, Compton-
thick (CT) active galactic nuclei (AGNs) out to a redshift of $z\leq$ 0.8.
From a sample of 540 AGNs selected by mid-IR (MIR) properties in observed
X-ray survey fields, we forward model the observed-to-intrinsic X-ray
luminosity ratio ($R_{L_{\text{X}}}$) with a Markov chain Monte Carlo (MCMC)
simulation to estimate the total fraction of CT AGNs ($f_{\text{CT}}$), many
of which are missed in typical X-ray observations. We create model
$N_{\text{H}}$ distributions and convert these to $R_{L_{\text{X}}}$ using a
set of X-ray spectral models. We probe the posterior distribution of our
models to infer the population of X-ray non-detected sources. From our
simulation we estimate a CT fraction of $f_{\text{CT}}$ =
$\text{0.555}^{+\text{0.037}}_{-\text{0.032}}$. We perform an X-ray stacking
analysis for sources in Chandra X-ray Observatory fields and find that the
expected soft (0.5–2 keV) and hard (2–7 keV) observed fluxes drawn from our
model to be within 0.48 and 0.12 dex of our stacked fluxes, respectively. Our
results suggests at least 50% of all MIR-selected AGNs, possibly more, are
Compton-thick ($N_{\text{H}}\gtrsim$ 10${}^{\text{24}}$ cm${}^{-\text{2}}$),
which is in excellent agreement with other recent work using independent
methods. This work indicates that the total number of AGNs is higher than can
be identified using X-ray observations alone, highlighting the importance of a
multiwavelength approach. A high $f_{\text{CT}}$ also has implications for
black hole (BH) accretion physics and supports models of BH and galaxy co-
evolution that include periods of heavy obscuration.

Active galactic nuclei (16) — Active galaxies (17) — Luminous infrared
galaxies (946) — Surveys (1671) — X-ray surveys (1824) — X-ray active galactic
nuclei (2035)

††journal: ApJ††facilities: Chandra, GALEX, MMT, NuSTAR, SDSS, UKIRT, WISE,
XMM, 2MASS††software: emcee: The MCMC Hammer (Foreman-Mackey et al., 2013),
xspec (Arnaud, 1996)

##  1 Introduction

Supermassive black holes (SMBHs) are found ubiquitously at the centers of
massive galaxies, growing via intermittent accretion of interstellar material
and occasional SMBH mergers. During intense periods of gas accretion, these
active galactic nuclei (AGNs) are extremely luminous, radiating across the
electromagnetic spectrum. Since the seminal discovery that obscuring dust is
responsible for the observed difference between broad-line and narrow-line
AGNs (Antonucci & Miller, 1985; see Ricci et al., 2017; Ananna et al., 2022a,
b for more recent models of unification), major efforts have been undertaken
to characterize the nature and extent of the obscuring material, with attempts
to determine its composition (e.g., Davies et al., 2015; García-Burillo et
al., 2016) and physical scale (e.g., Hönig et al., 2012; Tristram et al.,
2014; Chen et al., 2015; López-Gonzaga et al., 2016; Panagiotou & Walter,
2019). Obscuration across electromagnetic regimes is also thought to originate
from different sources, namely nuclear gas producing absorption at X-ray
energies and circumnuclear dust grains causing extinction of UV and optical
emission (see Hickox & Alexander, 2018).

Infrared surveys have provided evidence that more than half of the AGN
population exhibit signatures of obscuration (Lacy et al., 2013).
Optical–near-IR (NIR) colors have proven reliable in selecting obscured
sources in large AGN samples (e.g., Hickox et al., 2007; LaMassa et al., 2016)
but lack the power to accurately estimate levels of obscuration, particularly
in heavily obscured sources where the host galaxy dominates at optical and NIR
wavelengths. Mid-IR (MIR) techniques have further improved our ability to
probe heavily obscured AGN that are missed by optical surveys (e.g., Stern et
al., 2005, 2012; Mateos et al., 2013; Assef et al., 2015), as MIR observations
probe emission from dust and are relatively unaffected by dust obscuration.

While MIR observations have contributed reliable techniques for identifying
obscured AGNs, X-ray spectroscopic observations have become the dominant
method of probing nuclear obscuring gas content (see Hickox & Alexander,
2018), the bulk of which need not necessarily trace the distribution of dust.
Though the X-ray regime has become the favored method for measuring obscuring
column densities ($N_{\text{H}}$), recent work has shown that hard X-ray
energies ($\geq$2–10 keV) can remain undetected for heavily obscured AGN in
typical observations (e.g., Lansbury et al., 2015; Yan et al., 2019; Lambrides
et al., 2020). This provides a particular challenge for Compton-thick (CT)
densities (i.e., $N_{\text{H}}\geq$ 1.5$\times$10${}^{\text{24}}$
cm${}^{-\text{2}}$), where scattering optical depths reach order unity (e.g.,
Lansbury et al., 2015). As deep, hard X-ray observations are required to
positively identify and constrain the density of obscuring material
surrounding CT AGNs (e.g., Traina et al., 2021), large-scale searches have
been limited to smaller fields. To date there has been no reliable direct
measure of the complete fraction of CT AGN ($f_{\text{CT}}$), marking a major
uncertainty in the total number of accreting black holes (BHs) in the
universe. New techniques are therefore required to accurately measure the
line-of-sight gas densities for the most heavily obscured AGNs, and enable a
complete view of the population (e.g., Pfeifle et al., 2021).

One compelling approach involves utilizing the empirical relationship between
the X-ray and MIR luminosities ($L_{\text{X}}$ and $L_{\text{MIR}}$
respectively) of AGNs (i.e., Fiore et al., 2009; Stern, 2015; Chen et al.,
2017). The combination of X-ray and MIR luminosities has been shown to provide
a reliable method for identifying extremely obscured AGNs ($N_{\text{H}}>$
10${}^{\text{24}}$ cm${}^{-\text{2}}$) in large samples (e.g., Carroll et al.,
2021, hereafter C21). The $L_{\text{X}}$–$L_{\text{MIR}}$ relation provides a
method of inferring the intrinsic X-ray luminosity of an AGN through its MIR
luminosity ($L_{\text{X}}(L_{\text{MIR}})$), which is reasonably
straightforward to estimate given the quality and depth of all-sky IR data.
The ratio of the observed-to-intrinsic X-ray luminosity ($R_{L_{\text{X}}}$)
then provides insight into the amount of obscuring material present in a given
source and can be calculated directly from observables. X-ray spectral models
provide the final component, allowing us to understand the effects of
increasing $N_{\text{H}}$ on observed X-ray fluxes and ultimately
$R_{L_{\text{X}}}$.

Moreover, C21 showed that the same technique can be used to estimate lower
limits on $N_{\text{H}}$ for sources that are not detected in X-ray
observations, substituting X-ray flux limits in lieu of direct detections.
However, the $R_{L_{\text{X}}}$ distribution for X-ray detected and non-
detected sources differs drastically and cannot be directly compared, as
$R_{L_{\text{X}}}$ for non-detected sources are upper limits. As such,
$N_{\text{H}}$ can only be estimated with any given accuracy for X-ray
detected sources using this methodology. Therefore, a new statistical method
is required to include the X-ray non-detected sources in our estimates of the
distribution of $N_{\text{H}}$.

Previous studies have estimated the space density of AGNs—including CT
sources—predominantly through modeling of the X-ray luminosity function (XLF;
i.e., Ueda et al., 2014; Buchner et al., 2015; Aird et al., 2015; Fotopoulou
et al., 2016; Ananna et al., 2019), which describes the number density of AGNs
as a function of X-ray luminosity and redshift. This approach is dependent on
large samples of X-ray selected AGNs in soft and hard energies. Additionally,
this method is generally reliant on assuming an absorption function to account
for obscuration and/or modeling of obscured and unobscured sources separately.
Multiple observational constraints may be used in the generation of the XLF
through synthesis models, including X-ray number counts, spectral shape of the
integrated cosmic X-ray background (CXB), and the observed fraction of CT
AGNs. The latest synthesis models produce an XLF out to $z$ = 1 that accounts
for all observational constraints (Ananna et al., 2019, hereafter A19),
predicting a CT fraction of 0.50$\pm$0.09 for low- to moderate-luminosity AGNs
(log $L_{\text{X}}/\text{erg s}^{-\text{1}}<$ 43.6), and as high as
0.56$\pm$0.07 at higher luminosities (log $L_{\text{X}}/\text{erg
s}^{-\text{1}}>$ 43.6)—much higher than previous studies.

In this work, we aim to make a direct measurement of $f_{\text{CT}}$ and
determine the underlying $N_{\text{H}}$ distribution of heavily obscured
sources ($N_{\text{H}}\geq$ 10${}^{\text{24}}$ cm${}^{-{\text{2}}}$). Our
approach relies on multiwavelength archival data to account for AGNs which are
either undetected or poorly detected in the X-rays. We use a novel
multiwavelength, non-parametric approach to forward model the $N_{\text{H}}$
distribution by simulating the observed $R_{L_{\text{X}}}$ distribution in a
complete sample of AGNs.

This paper is organized as follows: Section 2 details our photometric sample
and matching X-ray fields; Section 3 describes our modeling procedure; Section
4 details our analysis and results of our modeling; Section 5 summarizes our
results and future work. Throughout the paper, we assume a $\Lambda$CDM
cosmology with parameters $H_{\text{0}}=$ 70 km s${}^{-{\text{1}}}$
Mpc${}^{-{\text{1}}}$, $\Omega_{\text{m}}=$ 0.3, and $\Omega_{\Lambda}=$ 0.7
(Spergel et al. 2007).

##  2 Data

The sample adopted in this paper builds on previous work by C21, having
identified heavily obscured AGNs from multiwavelength archival data. The
details of their data selection and modeling are fully described therein and
briefly summarized here.

Our sample consists of 540 AGNs selected on the basis of MIR colors (Assef et
al., 2018) from the Wide-field Infrared Survey Explorer (WISE; Wright et al.,
2010) and matched to Sloan Digital Sky Survey (SDSS; York et al., 2000)
optical counterparts. Each source in our sample also lies within fields
observed by the Chandra (Evans et al., 2020), XMM-Newton (Rosen et al., 2016),
and/or NuSTAR (Civano et al., 2015; Lansbury et al., 2017; Masini et al.,
2018) X-ray observatories. Although the 540 sources in our final sample all
lie within observed X-ray fields, only 40% have X-ray detections, with the
remaining 60% lacking X-ray counterparts. This working sample was drawn from a
larger parent dataset of over three million sources with redshift measurements
(photometric and spectroscopic), each with matching SDSS–WISE photometry.

Optical photometry for our sample was obtained from a combination of SDSS Data
Release 14 (Abolfathi et al., 2018) and XDQSOz (DiPompeo et al., 2015). In
this work, we chose to replace all MIR photometry from the WISE AllWISE Source
Catalog with unWISE (Lang, 2014). The unWISE Catalog produced WISE forced
photometry at known SDSS source positions, ensuring reliable source
association across catalogs. Use of the unWISE Catalog increases the depth of
WISE observations at wavelengths crucial to disentangling absorbed and non-
absorbed systems, all while maintaining the highest possible signal-to-noise
ratio (S/N). Additional photometry was added in UV and NIR bands where
available from the Galaxy Evolution Explorer (GALEX; Martin et al., 2005), the
UKIRT Infrared Deep Sky Survey (UKIDSS; Lawrence et al., 2007), and the Two
Micron All Sky Survey (2MASS; Skrutskie et al., 2006) (see C21 for additional
details).

To ensure the accuracy of our modeling and minimize systematics, we imposed
several quality selection criteria to our data. First, we required all
available photometry to have a minimum S/N $\geq$ 3.0. We required a minimum
of seven photometric bands per source for adequate coverage to model spectral
energy distributions (SEDs), with the additional requirement of detections in
all four WISE bands to accurately constrain AGN contribution to the MIR. SDSS
photometry was only considered where the PhotoObjAll column clean = 1,
ensuring reliable photometry. Though our data was primarily WISE selected, all
sources are detected in $gri$. Finally, we applied an angular bright stars
mask to remove regions of the sky with possible IR contamination (DiPompeo et
al., 2014). We chose to limit sources to redshifts $z\leq$ 0.8 and only used
SDSS photometric redshifts where the Photoz column photoErrorClass was $-$1,
1, 2, or 3, corresponding to an average root-mean-square error between 0.066
$\leq$ RMSE $\leq$ 0.074. Our final sample consists of sources having 64%
photometric and 36% spectroscopic redshifts.

Figure 1: Distribution of sample properties: redshift ($z$; top left);
observed 2–10 keV X-ray luminosity ($L_{\text{X}}$; top right); MIR AGN 6
$\mu$m luminosity ($L_{\text{MIR}}$; bottom left); intrinsic 2–10 keV X-ray
luminosity estimates ($L_{\text{X}}(L_{\text{MIR}})$; bottom right). In each
panel, the distributions are demarcated as X-ray detected (dashed blue) and
X-ray non-detected (solid orange) sources, along with the comprehensive
distribution (dotted gray).

Broadband SEDs were modeled using a combination of galaxy and AGN templates
(Assef et al., 2010; Kirkpatrick et al., 2015). From our SED modeling, sources
were only considered if they passed a minimum reduced chi-square cut
($\chi^{\text{2}}_{\text{red}}\leq$ 20.0) to reject poor fits. We also
required sources to derive a high fraction of their total MIR SED emission
from AGN contribution ($f_{\text{AGN}}\geq$ 0.7 at 15$\mu$m), and further
limited our final sample to IR-luminous sources ($L_{\text{MIR}}>$
10${}^{\text{42}}$ erg s${}^{-\text{1}}$).

All sources in our final sample contain MIR 6 $\mu$m AGN luminosities
($L_{\text{MIR}}$) drawn directly from our SED modeling. Each object lies
within an X-ray observed field and includes 3$\sigma$ upper-limit X-ray
luminosities ($L_{\text{X-lim}}$) for both X-ray detected and non-detected
sources, drawn from the respective flux limits of each instrument (see Figure
3 of C21). Observed 2–10 keV X-ray luminosities ($L_{\text{X}}$) were
calculated for the X-ray detected sources in our sample. The observed
$L_{\text{X}}$–$L_{\text{MIR}}$ relation of Chen et al. (2017) was used to
estimate intrinsic X-ray luminosities ($L_{\text{X}}$($L_{\text{MIR}}$)) for
all sources, corresponding to their derived MIR AGN luminosities. A breakdown
of the sample properties described in this section are shown in Figure 1 and
are separated for comparison into X-ray detected and non-detected sources.
Lastly, we defined the ratio of observed-to-intrinsic X-ray luminosity as

$R_{L_{\text{X}}}=\begin{dcases*}\frac{L_{\text{X}}}{L_{\text{X}}(L_{\text{MIR}})}&\quad
for X-ray detected,\\\
\frac{L_{\text{X-lim}}}{L_{\text{X}}(L_{\text{MIR}})}&\quad for X-ray non-
detected,\end{dcases*}$ (1)

which is used as a proxy for nuclear obscuration.

##  3 Modeling

###  3.1 X-Ray Spectral Models

We used a set of X-ray spectral models to estimate the impact of nuclear
obscuration on observed AGN X-ray flux. Our spectral models were constructed
using the X-ray spectral fitting package xspec (Arnaud, 1996) with the
following syntax:

$\textsc{constant}\times\textsc{cutoffpl}\\\
+\,\textsc{zphabs}\times\textsc{cabs}\times\textsc{pexmon}\\\
+\,\textsc{borus02}.$

Each of our models follow a similar prescription, being the sum of
transmitted, reflected, and scattered components to account for the diversity
of AGN X-ray spectra (e.g., Buchner et al., 2014; Ananna et al., 2020).

We began with an initial spectral model cutoffpl, consisting of a 300 keV
high-energy cutoff power law to represent the intrinsic continuum of
UV–optical photons upscattered to X-ray energies through Comptonization by the
AGN corona. To account for Thomson scattering of the X-ray continuum by
photoionized, circumnuclear material, we assumed the variable scattering
fraction of Gupta et al. (2021, hereafter, G21), applied as a multiplicative
constant to the cutoff power law. G21 observe that the scattering fraction
($f_{\text{scatt}}$) decreases with increasing $N_{\text{H}}$ and shows no
dependency on X-ray luminosity or BH mass. To encapsulate any other potential
factors due to AGN diversity in our sample, we chose to incorporate the G21
uncertainties on scattering fraction ($\sigma_{\text{scatt}}$) as a parameter
in our model. Though we technically employ the uncertainties on scattering
fraction ($f_{\text{scatt}}$ \+ $\sigma_{\text{scatt}}$) in our models, to
avoid unnecessary confusion we refer to this parameter as $f_{\text{scatt}}$
except where explicitly necessary.

Next, we considered Compton reflection of X-ray photons by the accretion disk.
A fraction of the high-energy continuum photons in proximity to the colder
material of the accretion disk are thought to be captured and reprocessed,
contributing to an X-ray spectral feature around 30 keV known as the Compton
reflection “hump.” For our models we chose to add the pexmon spectral
component, which assumes a simplistic, slab-like geometry of optically thick,
cold material. The strength, or scaling factor for reflection ($R$), was then
used as an additional parameter in our model. We then multiplied the reflected
component by zphabs and cabs to account for attenuation by photoelectric
absorption and Compton scattering, respectively.

Finally, we handled attenuation of emission by a more distant, circumnuclear
distribution of gas and dust, colloquially referred to as the obscuring torus.
For this, we chose to adopt the radiative transfer code of Baloković et al.
(2018), specifically borus02, representing the AGN torus as a uniform density
sphere with polar cutouts and a cutoff power-law intrinsic continuum. As the
majority of the attenuation in obscured systems is thought to be caused by
nuclear material within the gravitational influence of the central SMBH, the
high upper limit of the torus column density parameter of borus02 allows us to
estimate the effects of extremely high levels of obscuration. The borus02
model also accounts for Compton reflection by the torus, further contributing
to the Compton hump described above.

Starting at $z=$ 0, we increased the line-of-sight column density from
10${}^{\text{21}}$ cm${}^{-\text{2}}$ to 10${}^{\text{25}}$ cm${}^{-\text{2}}$
and calculated the change in the observed-to-intrinsic 2–10 keV flux ratio as
a function of $N_{\text{H}}$. For this calculation, we adopted the flux at
$N_{\text{H}}=$ 10${}^{\text{21}}$ cm${}^{-\text{2}}$ as the intrinsic,
unobscured X-ray flux. The simulated change in the X-ray spectral model flux
ratio provides us with a conversion between $N_{\text{H}}$ and
$R_{L_{\text{X}}}$, which we adopted for the body of our work. Though borus02
obscuration can be set at higher values, the distinction between sources with
column densities higher than the CT limit are not well constrained, thus we
chose $N_{\text{H}}$ = 10${}^{\text{25}}$ cm${}^{-\text{2}}$ as an upper limit
in our models.

We then investigated the effects of different model parameters on
$R_{L_{\text{X}}}$, including scattering fraction $f_{\text{scatt}}$,
reflection strength $R$, photon index ($\Gamma$), and torus opening angle
($\theta_{\text{OA}}$). Following the prescription above, we constructed
spectral models for each set of parameters {$f_{\text{scatt}}$, $R$, $\Gamma$,
$\theta_{\text{OA}}$}. To account for the extremely low scattering fraction
seen in some heavily obscured systems (e.g., Noguchi et al., 2010), we allowed
the scattering fraction to vary uniformly between $-$2.0
$\leq\sigma_{\text{scatt}}\leq$ 0.5 in steps of 0.5. We chose
0.5$\sigma_{\text{scatt}}$ as an upper limit as sources with little or no
obscuration would produce nonphysical results (i.e., total emission scattered
greater than 100%). We allowed the reflection scaling factor to vary between
0.0 $\leq R\leq$ 2.0 in steps of 0.2, where $R$ = 0.0 corresponds to Compton
reflection soley from the back-side of the torus (i.e., no disk reflection).
We allowed photon index to vary in steps of 0.1 between 1.4 $\leq\Gamma\leq$
2.2 to span the full range observed in AGN samples. Finally, we allowed the
opening angle to vary in steps of 10$\arcdeg$, spanning a range from
10$\arcdeg\leq\theta_{\text{OA}}\leq$ 80$\arcdeg$.

Figure 2: Observed-to-intrinsic X-ray luminosity ratio $R_{L_{\text{X}}}$ as a
function of obscuring column density $N_{\text{H}}$. The shaded region of
$R_{L_{\text{X}}}$–$N_{\text{H}}$ shown here represents the continuous
parameter space inferred from our full suite of X-ray spectral models and
parameters {$f_{\text{scatt}}$, $R$, $\Gamma$, $\theta_{\text{OA}}$}. Of these
parameters we find $f_{\text{scatt}}$ to have the most significant effect on
$R_{L_{\text{X}}}$, as seen by the striations, or channels, caused by our
choice of step size in $f_{\text{scatt}}$ parameter space. The remaining
parameters $R$, $\Gamma$, and $\theta_{\text{OA}}$ are responsible for the
narrow scatter within each of the deeper channels attributed to
$f_{\text{scatt}}$. The vertical shaded gray region demarcates Compton-thick
levels of obscuration ($N_{\text{H}}>$ 10${}^{\text{24}}$ cm${}^{-\text{2}}$).

The relation between $R_{L_{\text{X}}}$ and $N_{\text{H}}$ from our models is
shown in Figure 2. Of our choice of parameters, variance on $f_{\text{scatt}}$
had the most significant effect on $R_{L_{\text{X}}}$, causing large
striations, or channels, to develop in $R_{L_{\text{X}}}$–$N_{\text{H}}$.
Aside from $f_{\text{scatt}}$, our other parameters show little to no effect
on $R_{L_{\text{X}}}$. Additionally, Ananna et al. (2020) ruled out regions of
X-ray spectral parameter space where AGNs could not reliably reproduce the
CXB. As these parameters also produce no significant effect on our results, we
chose set values moving forward to save computational costs. For the remainder
of this work, we adopted $R_{L_{\text{X}}}$ values from our X-ray models where
$\Gamma$ = 1.9, $R$ = 0.99, and $\theta_{\text{OA}}$ = 60$\arcdeg$. A subset
of our X-ray spectral models are presented in Figure 3.

Figure 3: X-ray spectral energy distributions. Variation in the observed-frame
X-ray spectra is shown with changes to the absorbing column density
$N_{\text{H}}$. To assist with visualization, we show three cases:
$N_{\text{H}}$ = 10${}^{\text{21}}$ cm${}^{-\text{21}}$ (top); $N_{\text{H}}$
= 10${}^{\text{23}}$ cm${}^{-\text{2}}$ (middle); and $N_{\text{H}}$ =
10${}^{\text{25}}$ cm${}^{-\text{2}}$ (bottom). The total X-ray SED is shown
(solid blue) as well as the individual components: primary power-law continuum
disk emission (dashed orange), scattered primary emission (dash-dotted green),
and reprocessed torus emission (dotted magenta).

We then evaluated the observed-frame output of the model for a range of
redshifts to cover our sample redshift space (0.0 $\leq z\leq$ 0.8). At each
redshift, we calculated X-ray spectral model fluxes in the observed 0.5–2 and
2–7 keV bands. A comparison between the observed soft and hard fluxes with
rest-frame 2–10 keV flux provides us with a way to estimate 0.5–2 and 2–7 keV
fluxes we would expect to observe for a given model distribution of
$N_{\text{H}}$.

###  3.2 Modeling the Compton-thick Fraction

We used the Bayesian ensemble sampler emcee (Foreman-Mackey et al., 2013) to
forward model the observed $R_{L_{\text{X}}}$ distribution of our sample and
estimate the total fraction of CT AGNs. Using the connection between AGN X-ray
flux and nuclear obscuration explored in Section 3.1, we created model
$N_{\text{H}}$ distributions which we then converted to distributions in
$R_{L_{\text{X}}}$. We then used these model distributions
$M(R_{L_{\text{X}}})$ to infer the likelihood of observing our sample data.

At each step of the MCMC, we used the model parameter $f_{\text{CT}}$ to
create a random sample in $N_{\text{H}}$ covering a range of column densities
(see Figure 4). For column densities up to 10${}^{\text{24}}$
cm${}^{-\text{2}}$, we assumed the intrinsic column density distribution of
Ricci et al. (2017) for non-blazar AGNs. As C21 focused primarily on the
detection of heavily obscured sources, their final sample lacks a significant
fraction of true unobscured, Type-I AGNs ($N_{\text{H}}\approx$
10${}^{\text{20}}$ cm${}^{-\text{2}}$). Considering the lack of dynamic range
in $R_{L_{\text{X}}}$ where $N_{\text{H}}<$ 10${}^{\text{21}}$
cm${}^{-\text{2}}$, we chose to modify the input $N_{\text{H}}$ distribution
by removing all sources with column densities $N_{\text{H}}<$
10${}^{\text{21}}$ cm${}^{-\text{2}}$. We further tested our modeling with a
different input $N_{\text{H}}$ distributions (i.e., the NuSTAR informed
distribution of Lansbury et al., 2015) and found no significant change in our
results, indicating that uncertainties on the fraction of fully unobscured
sources does not play a significant role in our modeling and results. We made
limited assumptions regarding the number of sources and shape of the
underlying distribution for $N_{\text{H}}\geq$ 10${}^{\text{24}}$
cm${}^{-\text{2}}$. Specifically, we set the upper limit of our model to
$N_{\text{H}}$ = 10${}^{\text{26}}$ cm${}^{-\text{2}}$, as densities
approaching this limit have only recently been detected (e.g., Yan et al.,
2019). The fraction of heavily obscured sources was allowed to vary between
0.00 $<f_{\text{CT}}<$ 1.00, and the Compton scattering fraction
$f_{\text{scatt}}$ was allowed to vary within the G21 uncertainties. A flat
prior was assumed on both parameters.

Figure 4: Upper panel: Example $N_{\text{H}}$ model distributions from our
simulation, where each model is generated with varying $f_{\text{CT}}$. This
figure represents a small subset of possible $N_{\text{H}}$ distributions,
presented here in steps of $\Delta f_{\text{CT}}$ = 0.10 with uniformly
distributed number of sources at $N_{\text{H}}>$ 10${}^{\text{24}}$
cm${}^{-\text{2}}$ for clarity. Lower panel: Corresponding $R_{L_{\text{X}}}$
distributions.

Although $f_{\text{CT}}$ and $f_{\text{scatt}}$ are the only free parameters
in our simulation, we also considered the shape of the $N_{\text{H}}$
distribution for sources with high obscuration ($N_{\text{H}}\geq$
10${}^{\text{24}}$ cm${}^{-\text{2}}$). We allowed the fraction of sources
with column densities 10${}^{\text{24}}$–10${}^{\text{25}}$ cm${}^{-\text{2}}$
and 10${}^{\text{25}}$–10${}^{\text{26}}$ cm${}^{-\text{2}}$ ($f_{\text{24--
25}}$ and $f_{\text{25--26}}$, respectively) to vary, where the combined
fractions account for the full number of CT sources in any given model. We
performed multiple runs of our simulation to assess the influence a higher
fraction of extremely obscured sources ($f_{\text{24--25}}<f_{\text{25--26}}$)
would have on our parameter estimates. As shown in Figure 7, we find that
varying $f_{\text{24--25}}$ shows a negligible effect on $f_{\text{CT}}$, but
shares a positive correlation with $f_{\text{scatt}}$. As the number of
heavily obscured sources rises (i.e., $f_{\text{25--26}}$ increases), a slight
increase in $f_{\text{scatt}}$ is required to model the data (i.e., the
uncertainty on scatter fraction $\sigma_{\text{scatt}}$ moves closer to zero).
Variations on $f_{\text{24--25}}$ and $f_{\text{25--26}}$ produce trivial
change in our parameter estimation, therefore we adopt a uniform
$N_{\text{H}}$ distribution where $f_{\text{24--25}}$ = $f_{\text{25--26}}$
for our simulation. Additionally, we ran similar tests for a range of
parameters {$R$, $\Gamma$, $\theta_{\text{OA}}$}, the results of which are
shown in Figures 8 and 9.

Our simulation aimed to estimate $f_{\text{CT}}$ from a sample of 540 IR-
luminous AGNs, many of which lack X-ray counterparts. To accomplish this, we
modeled both the intrinsic number of sources expected in the survey as well as
the expected observed sample size, based on both our model and X-ray flux
sensitivity limits. We assumed a log-likelihood function of the form

$\ln\mathcal{L}=\sum^{N_{\text{det}}}_{i=1}\,\ln\int
M(R_{L_{\text{X}}})\,\mathcal{N}(R_{L_{\text{X}},i},\,0.23)\,\mathop{dR_{L_{\text{X}}}}\\\
-\int
M(R_{L_{\text{X}}})\,\sum^{N_{\text{all}}}_{i=1}\,S_{i}(R_{L_{\text{X-lim}}})\,\mathop{dR_{L_{\text{X}}}}$
(2)

to account for both X-ray detected and non-detected sources. A detailed
description of the two terms in our likelihood function follows. A visual
representation of $M(R_{L_{\text{X}}})$ and both terms of our likelihood are
presented in Figures 10, 11, and 12, respectively.

The first term of Equation 2 sums over all model distributions for the X-ray
detected sources in our sample. Here, $M(R_{L_{\text{X}}})$ represents the
normalized model distribution generated at each step in the simulation given a
choice of parameters {$f_{\text{CT}}$, $f_{\text{scatt}}$}. For each X-ray
detected source, we create a Gaussian distribution
$\mathcal{N}(R_{L_{\text{X}},i},\,\text{0.23})$ centered on the observed
$R_{L_{\text{X}}}$ of the source and an uncertainty of 0.23 dex—the intrinsic
dispersion in the $L_{\text{X}}$–$L_{\text{MIR}}$ relation. The convolution of
these two terms generates a probability density of observing a source given
the current choice of parameters at each step in the simulation. We then
integrated over all possible $R_{L_{\text{X}}}$ and summed over all of the
detected sources in the sample.

The second term of Equation 2 acts as a normalizing constant and accounts for
all sources in the sample. Again, $M(R_{L_{\text{X}}})$ represents the current
model at each step in the simulation, but this time before integrating it is
multiplied by the sensitivity function $S_{i}(R_{L_{\text{X-lim}}})$. To
account for observational bias in a flux-limited survey, a complicated
calculation of sample incompleteness is generally required. However, we use
the fact that $R_{L_{\text{X}}}$ is the logarithmic ratio of X-ray
luminosities to our advantage, and instead calculated $R_{L_{\text{X-lim}}}$
using the X-ray flux-limit luminosity of each source (see Equation 1). Using
$R_{L_{\text{X-lim}}}$ in this way allows us to estimate whether an object
would be detected at any given step in our simulation. As such, the
sensitivity function $S_{i}(R_{L_{\text{X-lim}}})$ is a Heaviside function,
shifted by $R_{L_{\text{X-lim}}}$, or

$S_{i}(R_{L_{\text{X-lim}}})=\Theta(R_{L_{\text{X}}}-R_{L_{\text{X-lim}},i}).$
(3)

As before, we convolved $M(R_{L_{\text{X}}})$ with
$S_{i}(R_{L_{\text{X-lim}}})$ and integrated over the full sample. The second
term of Equation 2 therefore serves as the expected observed sample size,
penalizing our likelihood function where $M(R_{L_{\text{X}}})$ does not
produce a comparable number of detected sources.

The ensemble sampler was run with 80 walkers, each with 1000 steps and a burn-
in period of 50 steps. Convergence of the walkers was assessed within emcee
using the Gelman-Rubin statistic. The posterior distributions of our
parameters were estimated by combining the results of all walkers using the
ensemble average (see Figure 5).

Figure 5: Posterior distributions of our modeling parameters, $f_{\text{CT}}$
and $\sigma_{\text{scatt}}$. The mean and covariance of each parameter is
shown alongside its respective distribution.

##  4 Results and Conclusions

The posterior mean $f_{\text{CT}}$ is estimated to be
$\text{0.555}^{+\text{0.037}}_{-\text{0.032}}$, which suggests that at least
half of all MIR-selected AGNs out to $z\leq$ 0.8 are heavily obscured, with
hydrogen column densities on the order of $\sim$10${}^{\text{24}}$
cm${}^{-{\text{2}}}$ or higher. Uncertainties on $f_{\text{CT}}$ were inferred
from our modeling procedure, which accounted for a variety of parameters
(e.g., scattering fraction) and assumptions (e.g., choice of input
$N_{\text{H}}$ distribution). Our results are in remarkable agreement with A19
that estimates a luminosity dependent CT fraction on the order of 0.50–0.56 at
comparable redshifts ($z\leq$ 1 compared to $z\leq$ 0.8 in this work). Most
notable is our arrival at nearly identical results using completely
independent methods—we make no use of AGN XLF, nor adjust our methodology to
account for the integrated CXB.

We further used our results to infer the distribution of AGN obscuration,
including contribution from unidentified CT AGNs. To derive a comprehensive
estimate, we used the intrinsic $N_{\text{H}}$ distribution of Ricci et al.
(2017)—complete with unobscured sources (i.e., $N_{\text{H}}<$
10${}^{\text{21}}$ cm${}^{-\text{2}}$)—as our input distribution. We generated
10,000 $N_{\text{H}}$ distributions following the prescription detailed in
Section 3, but instead fix the CT fraction to $f_{\text{CT}}$ =
$\text{0.555}^{+\text{0.037}}_{-\text{0.032}}$. We then converted
$N_{\text{H}}$ to $M(R_{L_{\text{X}}})$ assuming a Gaussian distribution
$\mathcal{N}(-\text{1},\,\text{0.8})$ for the uncertainties on
$f_{\text{scatt}}$ and the same choice of spectral parameters described in
Section 3.1. Each simulated source was then randomly assigned to an object in
our sample. We compared $M(R_{L_{\text{X}}})$ to $R_{L_{\text{X-lim}}}$ and
flagged each simulated source as “detected” where $M(R_{L_{\text{X}}})\geq
R_{L_{\text{X-lim}}}$ and “non-detected” where
$M(R_{L_{\text{X}}})<R_{L_{\text{X-lim}}}$. Finally, we averaged and
normalized ($N_{\text{H}}<$ 10${}^{\text{24}}$ cm${}^{-\text{2}}$) over all
$N_{\text{H}}$ distributions. Our results modeling the $N_{\text{H}}$
distribution of MIR-selected AGNs is presented in the lower panel of Figure 6.

Figure 6: Modeling results. Top left: Soft (0.5–2 keV) and hard (2–7 keV) band
stacked images for the X-ray non-detected sources in our sample using
StackFast (Ananna, 2023, in prep.). Source extraction regions are represented
by the solid blue circles, while background subtraction regions are depicted
by the dashed green annuli. Top right: Comparison of the X-ray stacked fluxes
(squares) to estimates on observed fluxes from our modeling (stars). Model
fluxes were plotted with horizontal offsets for visual clarity. Bottom panel:
$N{\text{H}}$ distribution with predicted CT fraction of
0.555${}^{+\text{0.037}}_{-\text{0.032}}$ from our simulation (solid line).
The estimated $N_{\text{H}}$ distribution of A19 is shown for reference
(dashed line), where the synthesis model estimates are represented as a
weighted average of our sample sources with low ($\log
L_{\text{X}}(L_{\text{MIR}})/\text{erg s}^{-1}<$ 43.6) and high ($\log
L_{\text{X}}(L_{\text{MIR}})/\text{erg s}^{-\text{1}}\geq$ 43.6) intrinsic
X-ray luminosity. The vertical shaded gray region demarcates Compton-thick
levels of obscuration ($N_{\text{H}}>$ 10${}^{\text{24}}$ cm${}^{-\text{2}}$).
Sources with $N_{\text{H}}<$ 10${}^{\text{20}}$ cm${}^{-{\text{2}}}$ (due to
scatter converting between $R_{L_{\text{X}}}$ and $N_{\text{H}}$) were
reassigned to $N_{\text{H}}$ = 10${}^{\text{20}}$ cm${}^{-{\text{2}}}$ for
clarity.

We then compared our inferred $N_{\text{H}}$ distribution to that of A19. As
their results are luminosity dependent, we averaged over both A19
$N_{\text{H}}$ distributions, weighted by the fraction of our sample sources
in each respective luminosity bin (0.21 and 0.79; split on log
$L_{\text{X}}$/erg s${}^{-{\text{1}}}$ = 43.6). Comparing these results in
Figure 6, we find a striking similarity between our work and that of A19.
Though previous estimates of the CT fraction using the AGN XLF have sometimes
yielded drastically different results, our determination of the CT fraction
and estimated $N_{\text{H}}$ distribution is in close agreement with A19,
using completely different measurement techniques and methodologies.

Table 1: X-ray stacking results and model comparison.

| $t_{\text{exp}}$ | Energy | $N_{\text{src}}$ | $N_{\text{bkg}}$ | $N_{\text{net}}$ | $F_{\text{X}}^{\text{stack}}$ | $F_{\text{X}}^{\text{model}}$
---|---|---|---|---|---|---|---
| [Ms] | [keV] |  |  |  | [erg s-1 cm-2] | [erg s-1 cm-2]
WISE AGN (X-ray non-det.) | 1.49 | 0.5–2 | 740 | 634 | 106 | 5.78$\pm$3.09$\times 10^{-16}$ | 1.74$\pm$1.02$\times 10^{-15}$
|  | 2–7 | 1328 | 1219 | 109 | 1.47$\pm$0.59$\times 10^{-15}$ | 1.94$\pm$1.11$\times 10^{-15}$

Note. — X-ray stacking results for energy ranges of 0.5–2 and 2–7 keV with
StackFast. For each energy range, we present exposure time, photon counts
(source, background, and net), and stacked flux estimates. The final column
depicts mean estimates of observed fluxes from our modeling.

Finally, we considered whether $f_{\text{CT}}$ and the $N_{\text{H}}$
distribution predicted by our modeling would produce similar soft and hard
X-ray fluxes comparable to measurements obtained from X-ray stacking. We
performed an X-ray stacking analysis of sources in Chandra fields using
stackfast (Ananna, 2023, in prep.). To ensure accuracy, we restricted our
comparison to sources in fields observed by Chandra, rather than XMM or
NuSTAR, as the high angular resolution of Chandra allows for more robust X-ray
stacking measurements. Of the 540 sources in our dataset, 151 sources lie
within observed Chandra fields (60 X-ray detected; 91 X-ray non-detected). For
each of our X-ray non-detected sources we extracted photon counts and exposure
times in both soft (0.5–2 keV) and hard (2–7 keV) energies. We computed fluxes
using count-rate-to-flux ratios characteristic of our X-ray spectral models
(see Section 3.1), and assumed a power law spectrum with photon index $\Gamma$
= 1.8 and typical Galactic absorption ($N_{\text{H}}$ = 10${}^{\text{20}}$
cm${}^{-\text{2}}$) as is typical of observed AGN X-ray spectra. We adopted
count-rate-to-flux conversion factor of 3.61$\times$10${}^{-\text{12}}$ erg
cm${}^{-\text{2}}$ count${}^{-\text{1}}$ and 7.44$\times$10${}^{-\text{12}}$
erg cm${}^{-\text{2}}$ count${}^{-\text{1}}$ for 0.5–2 and 2–7 keV,
respectively. Fluxes calculated from our count-rate-to-flux conversion are
representative of Chandra ACIS-I Cycle 12 responses—the current cycle for the
majority of the observations in our stacking analysis. Uncertainties on the
stacked fluxes were estimated by bootstrap resampling of input sources. The
calculated uncertainties therefore reflect the distribution of input fluxes
and are generally larger than Poisson uncertainties. Although our choice of
X-ray spectral models in Section 3.1 includes additional components, the
spectral shape chosen for our count-rate-to-flux conversion is consistent in
the energy range considered in this work (2–10 keV) and was deemed sufficient
for our analysis.

To estimate the observed X-ray fluxes from our model, we randomly assigned
each simulated sources previously flagged as non-detected to one of the 91
X-ray non-detected Chandra sources in our sample. We extracted rest-frame 2–10
keV luminosities for each model source using the combination of
$R^{\text{model}}_{L_{\text{X}}}$ and $L_{\text{X}}(L_{\text{MIR}})$. We then
converted the rest-frame 2–10 keV luminosities to observed 0.5–2 and 2–7 keV
fluxes using the X-ray spectral conversions discussed in Section 3.1. To
account for possible contamination in the data, we randomly set the X-ray flux
for 10% of all model sources to zero. We chose to remove 10% of model sources
to match the WISE AGN R90 catalog (Assef et al., 2018), which assesses 90%
reliability in their AGN selection methods. Though recent work has suggested
that the level of contamination may be higher (e.g., Hainline et al., 2016;
LaMassa et al., 2019), we found that varying the fraction of contamination
between 5%–20% had a negligible effect on our results. We repeated this
process 100 times and calculated the mean flux over all model sources.

The results of both our X-ray stacking and model estimates are presented in
Table 1, while a side-by-side visual comparison is also shown in the top panel
of Figure 6. We find our model fluxes to be in good agreement with the X-ray
stacking results to within 0.48 dex at soft energies (0.5–2 keV) and 0.12 dex
at hard energies (2–7 keV). While our hard X-ray fluxes are in excellent
agreement with our stacking results, our soft X-ray fluxes estimates are three
times higher than our stacking results. An increase in X-ray photons at soft
energies could possibly be due to overestimation of the reflection strength
parameter $R$ used in our spectral models. In fact, one postulated explanation
of the soft X-ray excess found in nearly half of all nearby AGNs is
upscattered X-ray photons reflecting back off of the accretion disk. Though
this model may be losing favor (e.g., Done & Nayakshin, 2007; Ursini et al.,
2020), it stands that reflection does contribute to an increase in soft
X-rays. The discrepancy between the soft X-ray fluxes of our model and
stacking may be caused by a relatively high amount of Compton scattering,
generally attributed to Type-I AGNs.

We also considered whether the X-ray selection of targeted sources may bias
our results. To test this we removed targets close to field centers (at both
0.5′ and 1′) and reran our analysis to track any changes. In the end, we found
no appreciable differences in our results. We find that the fraction of
detected sources actually increases linearly with distance from field center.
Because of this, we do not believe there is any significant bias introduced by
combining X-ray targeted sources with serendipitous sources. It is possible
that detection of point sources in the vicinity of targeted bright X-ray
objects could be limited by diffuse emission, but this is not something we can
currently account for.

##  5 Discussion

Our approach differs from previous work in distinct ways: modeling of the
observed $R_{L_{\text{X}}}$ distribution (i.e., a measurable proxy for nuclear
obscuration) in lieu of assumptions on $N_{\text{H}}$ distribution; the choice
of generalized X-ray spectral models, accounting for fewer parameters, degrees
of freedom, and potential degeneracies; constraints from X-ray stacking of
individual sources rather than an integrated spectrum, such as the CXB.
Whereas synthesis models may fit obscured and unobscured AGNs separately to
estimate the overall XLF, our technique relies solely on direct observables
(i.e., X-ray and infrared observations) to infer the overall distribution of
CT AGNs from simulated $N_{\text{H}}$ distributions. Our approach of forward
modeling the observed $R_{L_{\text{X}}}$ distribution allows us to deduce the
fraction of CT AGNs without assumptions and uncertain completeness corrections
on the number of sources or the shape of the $N_{\text{H}}$ distribution for
heavily obscured sources.

There are a number of broad implications for a high fraction of CT AGNs, but
one direct consequence involves BH accretion physics. A significant increase
in the number of heavily obscured AGNs, missing in typical X-ray observations,
may signify an increase in the average cosmic accretion efficiency driving BH
growth (i.e., Comastri et al., 2015). The BH mass function, which provides a
record of the mass of SMBHs in the local universe, in turn provides
information on the accretion history of SMBHs. An increase in the integrated
energy density due to AGNs would necessitate an increase in the efficiency of
accretion, under the assumption that the empirically derived BH mass function
is a statistically complete relation (i.e., Blecha et al., 2016). As the AGN
community converges on the true underlying population of actively accreting
SMBHs in the Universe, we can take further steps to painting a more complete
picture of galaxy and AGN co-evolution and the role AGNs play in other
observables, such as the CXB.

Furthermore, significant efforts have been made to understand the importance
of major mergers on the activity and duty cycle of AGNs. Though the
significance of mergers as a primary driver of SMBH growth remains a matter of
active investigation, clear evidence is seen for an increase in the fraction
of disturbed host morphologies in heavily obscured and CT AGNs (e.g., Kocevski
et al., 2015). While mergers are predicted to provide gas which may obscure
AGNs at host galaxy scales (e.g., Blecha et al., 2018), the results presented
in this work deal primarily with nuclear attenuation. A follow-up study on the
host morphologies of our sample, which include many CT AGNs yet undetected at
X-ray energies, could provide valuable insight on potential biases due to AGN
selection methods in such studies. If the “missing” CT AGNs in our sample are
found to have particular galaxy morphologies, this would provide further
evidence for a scenario in which CT AGN can be obscured by material on scales
beyond the central torus and may represent a distinct phase of SMBH growth
(Sanders et al., 1988).

Future work may attempt to extrapolate our estimate of the Compton-thick
fraction from a MIR-selected sample to the full AGN population. A full
accounting of observational biases must be considered to ensure a complete
working sample of AGNs collected from different selection techniques at
various wavelengths.

We thank the referee for exceptionally thoughtful and constructive comments
which led to a substantial improvement of this work. This research has made
use of X-ray data from the following: the NuSTAR mission, a project led by the
California Institute of Technology, managed by the Jet Propulsion Laboratory,
and funded by the National Aeronautics and Space Administration (NASA); the
Chandra Source Catalog 2.0, provided by the Chandra X-ray Center (CXC) as part
of the Chandra Data Archive; XMM-Newton, a European Space Agency (ESA) science
mission with instruments and contributions directly funded by ESA Member
States and NASA. This research also made use of data products from the Wide-
field Infrared Survey Explorer, which is a joint project of the University of
California, Los Angeles, and the Jet Propulsion Laboratory/California
Institute of Technology, and funded by NASA. This work is also based in part
on NIR data obtained as part of the UKIRT Infrared Deep Sky Survey, as well as
data products from the Two Micron All Sky Survey, which is a joint project of
the University of Massachusetts and the Infrared Processing and Analysis
Center/California Institute of Technology, funded by NASA and the National
Science Foundation. This research also made use of software provided by the
High Energy Astrophysics Science Archive Research Center (HEASARC), which is a
service of the Astrophysics Science Division at NASA/GSFC and the High Energy
Astrophysics Division of the Smithsonian Astrophysical Observatory. C.M.C.
acknowledges support from the John Templeton Foundation through LSSTC Catalyst
Fellowship grant number 62192. C.J.C., R.C.H., and T.T.A. acknowledge support
from the NSF through CAREER award number 1554584. R.C.H. and T.T.A.
acknowledge support from NASA through ADAP grant number 80NSSC19K0580. T.T.A.
also acknowledges support from NASA ADAP grant number 80NSSC23K0557. R.J.A.
was supported by FONDECYT grant number 1191124 and 1231718, and by ANID BASAL
project FB210003.

##  Appendix A Additional Parameter Tests

Figure 7: Effects on parameter estimates for a variable distribution of CT
sources. Left: Estimates on $f_{\text{CT}}$ for increasing fraction of heavily
obscured sources $f_{\text{25--26}}$. Right: Similar estimates on
$\sigma_{\text{scatt}}$.

Figure 8: Additional parameter tests—$f_{\text{CT}}$. Left: Compton-thick
fraction $f_{\text{CT}}$ as a function of $f_{\text{25--26}}$ {0.3, 0.5 ,0.7}
for an opening angle of $\theta_{\text{OA}}$ = 30$\arcdeg$. Variations on
reflection strength $R$ {0.5, 1.0} (open and closed symbols) and photon index
$\Gamma$ {1.8, 1.9, 2.0} (squares, circles, and stars) show negligible impact
on $f_{\text{CT}}$. Symbols are plotted with horizontal offset for visual
clarity. Right: Similar plot, but for an opening angle of $\theta_{\text{OA}}$
= 60$\arcdeg$.

Figure 9: Additional parameter tests—$f_{\text{scatt}}$. Left: Uncertainties
on scattering fraction $\sigma_{\text{scatt}}$ as a function of $f_{\text{25--
26}}$ {0.3, 0.5 ,0.7} for an opening angle of $\theta_{\text{OA}}$ =
30$\arcdeg$. Variations on reflection strength $R$ {0.5, 1.0} (open and closed
symbols) and photon index $\Gamma$ {1.8, 1.9, 2.0} (squares, circles, and
stars) show negligible impact on $\sigma_{\text{scatt}}$. Symbols are plotted
with horizontal offsets for visual clarity. Right: Similar plot, but for an
opening angle of $\theta_{\text{OA}}$ = 60$\arcdeg$.

##  Appendix B Log-Likelihood Function

Figure 10: Example $M(R_{L_{\text{X}}})$ given a set of parameters
{$f_{\text{CT}}$, $f_{\text{scatt}}$}.

Figure 11: Integrand of the likelihood function—first term. Left: Observed
$(R_{L_{\text{X}},i},\,\text{0.23})$ distribution. Right: Observed
$(R_{L_{\text{X}},i},\,\text{0.23})$ distribution convolved with the model
$M(R_{L_{\text{X}}})$. Note that each curve represents an X-ray detected
source from our sample. The first term of the likelihood function (Equation 2)
is the sum of the logarithm of the integral of these functions.

Figure 12: Integrand of the likelihood function—second term. Left: Sensitivity
function $S_{i}(R_{L_{\text{X-lim}}})$. Right: Sensitivity function
$S_{i}(R_{L_{\text{X-lim}}})$ convolved with the model $M(R_{L_{\text{X}}})$.
Note that each curve represents a different object in the sample (both X-ray
detected and non-detected). The second term of the likelihood function
(Equation 2) is the sum of these functions, integrated over
$R_{L_{\text{X}}}$.

## References

  * Abolfathi et al. (2018) Abolfathi, B., Aguado, D. S., Aguilar, G., et al. 2018, ApJS, 235, 42, doi: 10.3847/1538-4365/aa9e8a
  * Aird et al. (2015) Aird, J., Coil, A. L., Georgakakis, A., et al. 2015, MNRAS, 451, 1892, doi: 10.1093/mnras/stv1062
  * Ananna (2023, in prep.) Ananna, T. T. 2023, in prep., ApJ
  * Ananna et al. (2020) Ananna, T. T., Treister, E., Urry, C. M., et al. 2020, ApJ, 889, 17, doi: 10.3847/1538-4357/ab5aef
  * Ananna et al. (2019) —. 2019, ApJ, 871, 240, doi: 10.3847/1538-4357/aafb77
  * Ananna et al. (2022a) Ananna, T. T., Weigel, A. K., Trakhtenbrot, B., et al. 2022a, ApJS, 261, 9, doi: 10.3847/1538-4365/ac5b64
  * Ananna et al. (2022b) Ananna, T. T., Urry, C. M., Ricci, C., et al. 2022b, ApJ, 939, L13, doi: 10.3847/2041-8213/ac9979
  * Antonucci & Miller (1985) Antonucci, R. R. J., & Miller, J. S. 1985, ApJ, 297, 621, doi: 10.1086/163559
  * 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
  * Assef et al. (2018) Assef, R. J., Stern, D., Noirot, G., et al. 2018, ApJS, 234, 23, doi: 10.3847/1538-4365/aaa00a
  * Assef et al. (2010) Assef, R. J., Kochanek, C. S., Brodwin, M., et al. 2010, ApJ, 713, 970, doi: 10.1088/0004-637X/713/2/970
  * Assef et al. (2015) Assef, R. J., Eisenhardt, P. R. M., Stern, D., et al. 2015, ApJ, 804, 27, doi: 10.1088/0004-637X/804/1/27
  * Baloković et al. (2018) Baloković, M., Brightman, M., Harrison, F. A., et al. 2018, ApJ, 854, 42, doi: 10.3847/1538-4357/aaa7eb
  * Blecha et al. (2018) Blecha, L., Snyder, G. F., Satyapal, S., & Ellison, S. L. 2018, MNRAS, 478, 3056, doi: 10.1093/mnras/sty1274
  * Blecha et al. (2016) Blecha, L., Sijacki, D., Kelley, L. Z., et al. 2016, MNRAS, 456, 961, doi: 10.1093/mnras/stv2646
  * Buchner et al. (2014) Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125, doi: 10.1051/0004-6361/20132297110.48550/arXiv.1402.0004
  * Buchner et al. (2015) —. 2015, ApJ, 802, 89, doi: 10.1088/0004-637X/802/2/89
  * Carroll et al. (2021) Carroll, C. M., Hickox, R. C., Masini, A., et al. 2021, ApJ, 908, 185, doi: 10.3847/1538-4357/abd185
  * Chen et al. (2015) Chen, C.-T. J., Hickox, R. C., Alberts, S., et al. 2015, ApJ, 802, 50, doi: 10.1088/0004-637X/802/1/50
  * Chen et al. (2017) Chen, C.-T. J., Hickox, R. C., Goulding, A. D., et al. 2017, ApJ, 837, 145, doi: 10.3847/1538-4357/837/2/145
  * Civano et al. (2015) Civano, F., Hickox, R. C., Puccetti, S., et al. 2015, ApJ, 808, 185, doi: 10.1088/0004-637X/808/2/185
  * Comastri et al. (2015) Comastri, A., Gilli, R., Marconi, A., Risaliti, G., & Salvati, M. 2015, A&A, 574, L10, doi: 10.1051/0004-6361/201425496
  * Davies et al. (2015) Davies, R. I., Burtscher, L., Rosario, D., et al. 2015, ApJ, 806, 127, doi: 10.1088/0004-637X/806/1/127
  * DiPompeo et al. (2014) DiPompeo, M. A., Myers, A. D., Hickox, R. C., Geach, J. E., & Hainline, K. N. 2014, MNRAS, 442, 3443, doi: 10.1093/mnras/stu1115
  * DiPompeo et al. (2015) DiPompeo, M. A., Myers, A. D., Hickox, R. C., et al. 2015, MNRAS, 446, 3492, doi: 10.1093/mnras/stu2341
  * Done & Nayakshin (2007) Done, C., & Nayakshin, S. 2007, MNRAS, 377, L59, doi: 10.1111/j.1745-3933.2007.00303.x
  * Evans et al. (2020) Evans, I. N., Primini, F. A., Miller, J. B., et al. 2020, in American Astronomical Society Meeting Abstracts, Vol. 235, American Astronomical Society Meeting Abstracts #235, 154.05
  * Fiore et al. (2009) Fiore, F., Puccetti, S., Brusa, M., et al. 2009, ApJ, 693, 447, doi: 10.1088/0004-637X/693/1/447
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Fotopoulou et al. (2016) Fotopoulou, S., Buchner, J., Georgantopoulos, I., et al. 2016, A&A, 587, A142, doi: 10.1051/0004-6361/201424763
  * García-Burillo et al. (2016) García-Burillo, S., Combes, F., Ramos Almeida, C., et al. 2016, ApJ, 823, L12, doi: 10.3847/2041-8205/823/1/L12
  * Gupta et al. (2021) Gupta, K. K., Ricci, C., Tortosa, A., et al. 2021, MNRAS, 504, 428, doi: 10.1093/mnras/stab839
  * Hainline et al. (2016) Hainline, K. N., Reines, A. E., Greene, J. E., & Stern, D. 2016, ApJ, 832, 119, doi: 10.3847/0004-637X/832/2/119
  * Hickox & Alexander (2018) Hickox, R. C., & Alexander, D. M. 2018, ARA&A, 56, 625, doi: 10.1146/annurev-astro-081817-051803
  * Hickox et al. (2007) Hickox, R. C., Jones, C., Forman, W. R., et al. 2007, ApJ, 671, 1365, doi: 10.1086/523082
  * Hönig et al. (2012) Hönig, S. F., Kishimoto, M., Antonucci, R., et al. 2012, ApJ, 755, 149, doi: 10.1088/0004-637X/755/2/149
  * Kirkpatrick et al. (2015) Kirkpatrick, A., Pope, A., Sajina, A., et al. 2015, ApJ, 814, 9, doi: 10.1088/0004-637X/814/1/9
  * Kocevski et al. (2015) Kocevski, D. D., Brightman, M., Nandra, K., et al. 2015, ApJ, 814, 104, doi: 10.1088/0004-637X/814/2/104
  * Lacy et al. (2013) Lacy, M., Ridgway, S. E., Gates, E. L., et al. 2013, ApJS, 208, 24, doi: 10.1088/0067-0049/208/2/24
  * LaMassa et al. (2019) LaMassa, S. M., Georgakakis, A., Vivek, M., et al. 2019, ApJ, 876, 50, doi: 10.3847/1538-4357/ab108b
  * LaMassa et al. (2016) LaMassa, S. M., Civano, F., Brusa, M., et al. 2016, ApJ, 818, 88, doi: 10.3847/0004-637X/818/1/88
  * Lambrides et al. (2020) Lambrides, E. L., Chiaberge, M., Heckman, T., et al. 2020, ApJ, 897, 160, doi: 10.3847/1538-4357/ab919c
  * Lang (2014) Lang, D. 2014, AJ, 147, 108, doi: 10.1088/0004-6256/147/5/108
  * Lansbury et al. (2015) Lansbury, G. B., Gandhi, P., Alexander, D. M., et al. 2015, ApJ, 809, 115, doi: 10.1088/0004-637X/809/2/115
  * Lansbury et al. (2017) Lansbury, G. B., Alexander, D. M., Aird, J., et al. 2017, ApJ, 846, 20, doi: 10.3847/1538-4357/aa8176
  * Lawrence et al. (2007) Lawrence, A., Warren, S. J., Almaini, O., et al. 2007, MNRAS, 379, 1599, doi: 10.1111/j.1365-2966.2007.12040.x
  * López-Gonzaga et al. (2016) López-Gonzaga, N., Burtscher, L., Tristram, K. R. W., Meisenheimer, K., & Schartmann, M. 2016, A&A, 591, A47, doi: 10.1051/0004-6361/201527590
  * Martin et al. (2005) Martin, D. C., Fanson, J., Schiminovich, D., et al. 2005, ApJ, 619, L1, doi: 10.1086/426387
  * Masini et al. (2018) Masini, A., Civano, F., Comastri, A., et al. 2018, ApJS, 235, 17, doi: 10.3847/1538-4365/aaa83d
  * Mateos et al. (2013) Mateos, S., Alonso-Herrero, A., Carrera, F. J., et al. 2013, MNRAS, 434, 941, doi: 10.1093/mnras/stt953
  * Noguchi et al. (2010) Noguchi, K., Terashima, Y., Ishino, Y., et al. 2010, ApJ, 711, 144, doi: 10.1088/0004-637X/711/1/144
  * Panagiotou & Walter (2019) Panagiotou, C., & Walter, R. 2019, A&A, 626, A40, doi: 10.1051/0004-6361/201935052
  * Pfeifle et al. (2021) Pfeifle, R. W., Ricci, C., Boorman, P. G., et al. 2021, arXiv e-prints, arXiv:2102.04412.  https://arxiv.org/abs/2102.04412
  * Ricci et al. (2017) Ricci, C., Trakhtenbrot, B., Koss, M. J., et al. 2017, Nature, 549, 488, doi: 10.1038/nature23906
  * Rosen et al. (2016) Rosen, S. R., Webb, N. A., Watson, M. G., et al. 2016, A&A, 590, A1, doi: 10.1051/0004-6361/201526416
  * Sanders et al. (1988) Sanders, D. B., Soifer, B. T., Elias, J. H., et al. 1988, ApJ, 325, 74, doi: 10.1086/165983
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163, doi: 10.1086/498708
  * Spergel et al. (2007) Spergel, D. N., Bean, R., Doré, O., et al. 2007, ApJS, 170, 377, doi: 10.1086/513700
  * Stern (2015) Stern, D. 2015, ApJ, 807, 129, doi: 10.1088/0004-637X/807/2/129
  * Stern et al. (2005) Stern, D., Eisenhardt, P., Gorjian, V., et al. 2005, ApJ, 631, 163, doi: 10.1086/432523
  * Stern et al. (2012) Stern, D., Assef, R. J., Benford, D. J., et al. 2012, ApJ, 753, 30, doi: 10.1088/0004-637X/753/1/30
  * Traina et al. (2021) Traina, A., Marchesi, S., Vignali, C., et al. 2021, ApJ, 922, 159, doi: 10.3847/1538-4357/ac1fee
  * Tristram et al. (2014) Tristram, K. R. W., Burtscher, L., Jaffe, W., et al. 2014, A&A, 563, A82, doi: 10.1051/0004-6361/201322698
  * Ueda et al. (2014) Ueda, Y., Akiyama, M., Hasinger, G., Miyaji, T., & Watson, M. G. 2014, ApJ, 786, 104, doi: 10.1088/0004-637X/786/2/104
  * Ursini et al. (2020) Ursini, F., Petrucci, P. O., Bianchi, S., et al. 2020, A&A, 634, A92, doi: 10.1051/0004-6361/201936486
  * Wright et al. (2010) Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868, doi: 10.1088/0004-6256/140/6/1868
  * Yan et al. (2019) Yan, W., Hickox, R. C., Hainline, K. N., et al. 2019, ApJ, 870, 33, doi: 10.3847/1538-4357/aaeed4
  * York et al. (2000) York, D. G., Adelman, J., Anderson, Jr., J. E., et al. 2000, AJ, 120, 1579, doi: 10.1086/301513