BAT AGN Spectroscopic Survey-XXIII. A New Mid-Infrared Diagnostic for Absorption in Active Galactic Nuclei

In this study, we use the SWIFT/BAT AGN sample, which has received extensive multiwavelength follow-up analysis as a result of the BAT AGN Spectroscopic Survey (BASS), to develop a diagnostic for nuclear obscuration by examining the relationship between the line-of-sight column densities ($N_{\rm{H}}$), the 2-10 keV-to-$12\,\rm{\mu m}$ luminosity ratio, and WISE mid-infrared colors. We demonstrate that heavily obscured AGNs tend to exhibit both preferentially ''redder'' mid-infrared colors and lower values of $L_{\rm{X,\,Obs.}}$/$L_{12\,\rm{\mu m}}$ than less obscured AGNs, and we derive expressions relating $N_{\rm{H}}$ to the $L_{\rm{X,\,Obs.}}$/$L_{12\,\rm{\mu m}}$ and $L_{22\,\rm{\mu m}}$/$L_{4.6\,\rm{\mu m}}$ luminosity ratios as well as develop diagnostic criteria using these ratios. Our diagnostic regions yield samples that are $\gtrsim80$% complete and $\gtrsim60$% pure for AGNs with log($N_{\rm{H}})\geq24$, as well as $\gtrsim85$% pure for AGNs with $\rm{log}(N_{\rm{H}})\gtrsim23.5$. We find that these diagnostics cannot be used to differentiate between optically star forming galaxies and active galaxies. Further, mid-IR contributions from host galaxies that dominate the observed $12~\rm{\mu m}$ emission can lead to larger apparent X-ray deficits and redder mid-IR colors than the AGNs would intrinsically exhibit, though this effect helps to better separate less obscured and more obscured AGNs. Finally, we test our diagnostics on two catalogs of AGNs and infrared galaxies, including the XMM-Newton XXL-N field, and we identify several known Compton-thick AGNs as well as a handful of candidate heavily obscured AGNs based upon our proposed obscuration diagnostics.


INTRODUCTION
Known to reside at the centers of most galaxies (e.g. Ferrarese & Merritt 2000;Kormendy & Ho 2013), super-massive black holes (SMBHs) grow and evolve through periods of activity characterized by the accretion of large quantities of gas. Classically, these active galactic nuclei (AGNs) are categorized based upon the characteristics of their optical spectroscopic emission lines, where the apparent differences between AGNs may be reconciled through a unification scheme involving a dusty obscuring torus (e.g. Antonucci 1993;Urry & Padovani 1995;Netzer 2015;Ramos Almeida & Ricci 2017), for which different inclination angles of the torus correspond to the observation of different AGN classes.
One ubiquitous observational signature of accretion onto SMBHs is X-ray emission, produced very close to the accretion disk (Fabian et al. 2009) due to inverse Compton scattering of optical and ultraviolet (UV) photons from the accretion disk by hot electrons in the corona (Haardt & Maraschi 1991, 1993. In the X-ray band, the line-of-sight gas column density, N H , is largely transparent to the 2-10 keV X-ray flux even up to column densities of a few times 10 23 cm −2 , however, significant attenuation and reprocessing of the 2-10 keV Xray emission does occur for Compton-thick (CT) AGNs, which possess gas column densities of 10 24 cm −2 (e.g. Lansbury et al. 2014;Ricci et al. 2015;Bauer et al. 2015;Puccetti et al. 2016;LaMassa et al. 2019;Toba et al. 2020). CT AGNs even at low redshift have thus proven to be very difficult to find and characterize (Alexander & Hickox 2012) using lower energy X-ray observatories such as Chandra and XMM-Newton, since the X-ray flux below 10 keV suffers significant photoelectric absorption and Compton scattering. This prevents the detection of some sources, and even those detected have fewer observed photons, which reduces the accuracy of spectral modeling. Important spectral signatures used to characterize the nature of AGNs can be missed without higher energy X-ray observations (Lansbury et al. 2015;Koss et al. 2016;Ricci et al. 2017a).
CT AGNs are particularly important among the general AGN population, as large fractions of CT AGNs are required to reproduce the observed Cosmic X-ray Background (CXB, Gilli et al. 2007;Buchner et al. 2015;Ananna et al. 2019). CT AGNs likely represent a significant fraction of the total intrinsic AGN population in the local Universe (∼20%, Burlon et al. 2011;∼27%, Ricci et al. 2015), with the most recent SMBH synthesis model developed by Ananna et al. (2019) suggesting that CT AGNs represent 50 ± 9% and 56 ± 7% of the total intrinsic AGN population up to both z ∼ 0.1 and z ∼ 1.0, respectively, though many previous synthesis models predicted lower intrinsic fractions of CT AGNs than this, including but not limited to Gilli et al. (2007); Treister et al. (2009); Akylas et al. (2012); Ueda et al. (2014). Furthermore, questions remain regarding the exact nature of the obscuring structure and how it relates to, for example, the host environment: some CT AGNs could represent the evolutionary phase predicted to occur as a result of galaxy mergers (e.g. Hopkins et al. 2008;Kocevski et al. 2015;Ricci et al. 2017b;Blecha et al. 2018). Identifying further cases of CT AGNs is crucial for providing a full census of accreting SMBHs, placing constraints on the CXB, as well as placing constraints on evolutionary and unification models.
UV radiation from the accretion disk is also reprocessed by the obscuring dusty torus, wherein the radiation is scattered and absorbed by the dust grains and reemitted thermally with a peak usually at mid-infrared (mid-IR) wavelengths. While the classically accepted origin of the mid-IR emission is the dusty torus itself, recent high angular resolution infrared observations of AGNs suggest that the mid-IR emission is in fact dominated by a dusty polar outflow rather than the torus itself (e.g. Hönig et al. 2012Hönig et al. , 2013Tristram et al. 2014;Asmus et al. 2016;Hönig & Kishimoto 2017;Stalevski et al. 2018;Hönig 2019). Other recent studies (Baron & Netzer 2019a,b) have actually attributed mid-IR emission to dusty outflows located on the order of tens to hundreds of parsecs from the centers of the galaxies.
A correlation between the intrinsic (unabsorbed) hard X-ray 2-10 keV luminosity (L X ) and mid-IR luminosity (L MIR ) of AGNs was reported as early as Elvis et al. (1978). Universally, studies find no difference (< 0.3 dex) in the ratios of L X to L MIR between Type 1 and Type 2 AGNs (Lutz et al. 2004;Levenson et al. 2009;Gandhi et al. 2009;Ichikawa et al. 2012;Mateos et al. 2015;Asmus et al. 2015) and this ratio is also insensitive to the neutral gas column density along the lineof-sight (Gandhi et al. 2009;Mateos et al. 2015;Asmus et al. 2015), even in the case of CT AGNs after correcting the X-rays for absorption (Gandhi et al. 2009). Many previous studies have also pointed out that this relation can serve as a useful tool to select obscured, particularly CT, AGNs, because CT AGNs tend to exhibit severe deficits in their absorbed X-ray emission when compared to their mid-IR emission, and thus CT AGNs fall significantly off of the L X -to-L MIR relation (Alexander et al. 2008;Goulding et al. 2011;Georgantopoulos et al. 2011;Rovilos et al. 2014;Asmus et al. 2015). Moreover, Asmus et al. (2015) also demonstrated that the ratio of L X /L MIR can be used to predict column densities for significantly obscured objects, deriving an equation that relates L X /L MIR to log(N H /cm −2 ) for log(N H /cm −2 ) > 22.8.
It is important to gather large samples of powerful AGNs across a range of column densities to test the utility of the X-ray to mid-IR relation as a tracer of nuclear obscuration. Hard X-ray selection provides one of the least biased methods of identifying powerful AGNs, as hard X-ray emission is largely unaffected by the lineof-sight obscuration for column densities < 10 24 cm −2 (see Fig. 1 from Ricci et al. 2015). The Swift/BAT ultrahard X-ray (14-195 keV) all-sky survey has dramatically increased the number of known hard X-ray extragalactic sources (Baumgartner et al. 2013;Oh et al. 2018), and has therefore been the focus of a large multiwavelength follow-up campaign (the BAT AGN Spectroscopic Survey, or BASS 1 ) designed to characterize the most powerful AGNs in the local Universe Ricci et al. 2017a;Lamperti et al. 2017). A second release of optical spectroscopy (BASS DR2) will also soon be publicly available (Koss et al., in prep; Oh et al., in prep). Ricci et al. (2017a) presented a detailed X-ray analysis of 838 ultra hard X-ray detected Swift/BAT AGNs, providing constraints on column densities (N H ) as well as the absorbed (observed) and unabsorbed 2-10 keV luminosities, while Ichikawa et al. (2017) provided mid-IR to far-infrared photometric data and corresponding IR luminosities for the 604 mid-IR-detected Swift/BAT AGNs. These two catalogs provide precisely the information needed to test the relationship between the absorbed hard X-ray and mid-IR emission with respect to the line-of-sight obscuration in AGNs.
In this paper, we present a new diagnostic for absorption in AGNs which combines the power of the known L X /L 12 µm correlation with WISE mid-IR colors. Using multiwavelength catalogs available for the Swift/BAT hard X-ray-selected sample of AGNs, we show this diagnostic reliably identifies the most obscured AGNs, at least for nearby X-ray bright AGN. Our proposed diagnostics could prove valuable in the search for obscured AGNs in the ongoing eROSITA survey (Predehl et al. 2010;Merloni et al. 2012). In Section 2 we describe our sample. In Section 3 we describe the analysis of the sample and propose our new absorption diagnostics as well as develop expressions that constrain column densities. In Section 4 we explore the emission ratios of opticallyselected star forming galaxies, we compare our diagnostics for obscuration to other recent studies, and we apply our diagnostics to the XMM-Newton XXL North field (Pierre et al. 2016(Pierre et al. , 2017. In Section 5 we detail our conclusions. Throughout this manuscript we assume the following cosmological values: H 0 = 70 km s −1 Mpc −1 , Ω M = 0.3, Ω Λ = 0.7. All luminosities quoted in this work are given in units of erg s −1 .

SAMPLE CONSTRUCTION
We selected our sample from the 70-month Swift/BAT X-ray properties catalog (Ricci et al. 2017a) which de-1 https://www.bass-survey.com/ tails the broadband 0.3-150 keV X-ray spectral properties of the 838 AGNs detected in the ultra hard X-ray 14-195 keV band by Swift/BAT and reported in the 70month source catalog (Baumgartner et al. 2013). We matched this catalog to the 70-month Swift/BAT infrared catalog of Ichikawa et al. (2017), which provides the complete near-infrared to far-infrared photometry for 604 Swift/BAT non-beamed AGNs at high Galactic latitudes (|b|>10 • ). We refer the reader to Ricci et al. (2017a) and Ichikawa et al. (2017) for further details on the construction of these catalogs. These catalogs yielded a parent sample of 604 non-beamed AGNs; any systems flagged as beamed in the Ricci et al. (2017a) catalog were removed during the matching process.
In order to conduct our analysis, we required X-ray and mid-IR detections in all four WISE bands for the AGNs in our sample, which excluded another 78 AGNs from the final sample 2 . We adopted the hard X-ray 2-10 keV luminosities (observed, uncorrected for intrinsic absorption) from Ricci et al. (2017a) and the infrared luminosities in all four WISE bands (3.6 µm, 4.6 µm, 12 µm, and 22 µm) from Ichikawa et al. (2017), which are not corrected for any host galaxy contributions. We further limited the sample to z < 0.1 AGNs to avoid redshift effects; this redshift cut removed another 70 AGNs from the sample. The Ricci et al. (2017a) catalog includes independent estimates of N H from a torus model in the event that the column density found with the phenomenological model was ≥ 10 24 cm −2 . In these cases, we instead use the column density inferred from the torus model, because torus models more accurately account for the 2-10 keV emission of CT AGNs.
Following the matching of the catalogs and the application of the above requirements, the final sample was composed of 456 nearby, non-beamed AGNs with a median redshift of z 0.032 and mid-IR luminosities in the range 1.6 × 10 39 ≤ L 12 µm /erg s −1 ≤ 8.5 × 10 44 .

DATA ANALYSIS
Comparing the observed X-ray 2-10 keV luminosities (L X, Obs. ) to the 12 µm luminosities (L 12 µm ), we see in Figure 1 that unobscured Swift/BAT AGNs tend to follow the relation between the intrinsic (unabsorbed) 2-10 X-ray and nuclear 12 µm luminosities (dashed black line) established by Asmus et al. (2015), whereas obscured AGNs appear X-ray suppressed when compared  Figure 1. The observed 2-10 keV X-ray vs. 12 µm luminosities. We color code the data points using the derived column density from the Ricci et al. (2017a) catalog, where the color map is denoted on the auxiliary axis. AGNs with NH upper limits of log(NH/cm −2 ) ≤ 20.0 are denoted with gray triangles. When comparing the 2-10 keV X-ray luminosity to the 12 µm luminosity for our sample of 456 Swift/BAT AGNs, we see a general decrease in the X-ray-to-mid-IR ratio with increasing column density, as expected. The relation between the AGN intrinsic 2-10 keV luminosity (corrected for absorption) and the nuclear 12 µm luminosity derived by Asmus et al. (2015) is represented by a dashed black line, whereas we plot the logarithmic ratio of log(L X, Obs. /L12 µm) = −1.3 with a black dotted line. Most heavily obscured (CT) AGNs exhibit luminosity ratios log(L X, Obs. /L12 µm) < −1.3.
to L 12 µm . This decrease in L X, Obs. compared to L 12 µm with increasing column density is expected, of course, since the X-ray emission will suffer greater attenuation than the mid-IR emission, and we color code the data in Figure 1 according to the column densities adopted in Section 2. As a population, CT AGNs are generally the furthest offset from the Asmus et al. (2015) relation, exhibiting luminosity ratios of log(L X, Obs. /L 12 µm ) < −1.3 (this was previously discussed in, e.g. Alexander et al. 2008); we plot this ratio between log(L X, Obs. ) and log(L 12 µm ) as a dotted black line in Figure 1. As has been done in previous works (e.g. Alexander et al. 2008;Goulding et al. 2011;Asmus et al. 2015), we can use this ratio between L X, Obs. and L 12 µm as a diagnostic tool to differentiate between less obscured and heavily obscured or CT AGNs. Ratios of two mid-IR luminosities, sufficiently separated in wavelength, could exhibit the same trend as is observed for the L X, Obs. /L 12 µm ratio, in that the shorter wavelength mid-IR emission could appear suppressed compared to the longer wavelength emission due to the obscuring material surrounding the AGN. We test a new diagnostic tool based on the ratio of the WISE 22 µm to 4.6 µm luminosity (L 22 µm /L 4.6 µm ) in Figure 2, and we find that the logarithmic ratio increases with column density, with the most significant increase in WISE ratio corresponding to the highest column densities (although with significant scatter). For example, we find a mean luminosity ratio and 1σ uncertainty of log(L 22 µm /L 4.6 µm ) = 0.50 ± 0.32 for AGNs with log(N H /cm −2 ) ≥ 24, whereas we find a mean ratio and 1-σ uncertainty of log(L 22 µm /L 4.6 µm ) = 0.06 ± 0.25 for unobscured AGNs with column densities log(N H /cm −2 ) < 22. This suggests that WISE colors can also be used as a diagnostic tool for identifying heavily absorbed AGNs. We combine the L X, Obs. /L 12 µm and L 22 µm /L 4.6 µm ratio diagnostics in Figure 3, plotting the two diagnostic ratios against one another in Panel A and color coding the data points by column density on the auxiliary axis. As the column density increases, the AGNs tend to exhibit lower values of L X, Obs. /L 12 µm and higher values of L 22 µm /L 4.6 µm , and thus we find that the most heavily obscured AGNs predominantly occupy the lower right portion of the parameter space (i.e. the largest X-ray deficits and highest L 22 µm /L 4.6 µm ratios). In Panel B of  The logarithmic L X, Obs. /L12 µm and L22 µm/L4.6 µm luminosity ratios (Panel A) with each point color coded to indicate the column density (as indicated by the auxiliary axis). Sources with NH upper limits of log(NH/cm −2 ) < 20.0 are denoted with gray triangles. Binning by the L22 µm/L4.6 µm luminosity ratio (Panel B), we see a general trend of decreasing X-ray-to-mid-IR ratio with increasing values of the L22 µm/L4.6 µm ratio. Binning by the L22 µm/L4.6 µm luminosity ratio and comparing it to the column density (Panel C), we see a general trend of increasing WISE luminosity ratios with increasing column density. Solid error bars in Panels B and C represent the standard error of the mean while dashed error bars represent the standard deviation computed for the respective bin.  . In an identical fashion to that shown in Figure 3, we show the logarithmic L X, Obs. /L12 µm and L12 µm/L4.6 µm luminosity ratios (Panel A) with the column density given on the auxiliary axis. We find very similar results to those shown in Figure 3 when binning and comparing the L X, Obs. /L12 µm and L12 µm/L4.6 µm luminosity ratio (Panel B) and when binning and comparing the L12 µm/L4.6 µm ratio and NH (Panel C).
and examine the scatter in the column density by bin; we notice a positive correlation between the column density and the mid-IR luminosity ratio. All solid error bars in Panels B and C represent the standard error of the mean while dashed error bars represent the standard deviation computed for the respective bin. Considering these three panels together, we may define a parameter space using the L X, Obs. /L 12 µm and the L 22 µm /L 4.6 µm luminosity ratio, which could be used to identify the most heavily absorbed AGNs. We repeated this analysis for an alternative WISE luminosity ratio of L 12 µm /L 4.6 µm and found very similar results (see Figure 4). We explored several AGN selection methods to potentially mitigate or at least account for the scatter observed, which we discuss in the Appendix. Ultimately, we did not apply any selection criteria to our sample of Swift/BAT AGNs during the analysis described below. Of important note, however, is the interesting result that the correlation between the mid-IR color and N H holds true for both WISE (selected via W 1 − W 2 > 0.8; Stern et al. 2012) as well as non-WISE AGNs.
3.1. A Relation for N H as a Function of L X, Obs. /L 12 µm Since the L X, Obs. /L 12 µm ratio is known to correlate with obscuring column (e.g. Ichikawa et al. 2012;Asmus et al. 2015;Yan et al. 2019), we derived an expression to describe this relationship using the Swift/BAT sample studied here. Note that, in contrast to Asmus et al. (2015) (see Section 4.2 for further details), we do not exclude sources with log(N H /cm −2 ) < 22.8 from the fitting process.
For our fitting process, we incorporate the luminosity ratios and the associated uncertainties. We pull the uncertainties in the mid-IR fluxes from Ichikawa et al. (2017) and Ichikawa et al. (2019); while uncertainties in the observed X-ray luminosities were not available in the Ricci et al. (2017a) catalog, we conservatively adopt uncertainties of log(L X, Obs. ) = ± 0.1 for AGNs with log(N H /cm −2 ) < 24 and log(L X, Obs. ) = ± 0.3 for AGNs with log(N H /cm −2 ) ≥ 24. In order to take into account the asymmetric uncertainties associated with log(N H /cm −2 ) and log(L X, Obs. /L 12 µm ), we employed the following Monte Carlo fitting routine: 1. Bootstrap: To fully account for the effect of outliers, we generated a new data set as a sample of the original, allowing repeats.
2. Monte Carlo: For each point in the bootstrapped data set, we generate a new point given the uncertainties in log(N H /cm −2 ) and log(L X, Obs. /L 12 µm ). We incorporate asymmetric error bars by randomly drawing from a Gaussian distribution separately for the negative and positive error bars. For upper limits in log(N H ) 3 , we generate a new point using a uniform distribution from the limit to an arbitrarily small log(N H /cm −2 ) value below it. After experimenting a number of different values, each giving similar results, we settled for a lower bound of log(N H /cm −2 ) = 19, which is only marginally lower than the smallest log(N H /cm −2 ) value through the Milky Way according to the maps by (Kalberla et al. 2005).

Orthogonal Distance Regression: We fit each
Monte Carlo-bootstrapped data set with a function of the form y = a · 10 b·(x−20) + c using orthogonal distance regression (from the Python scipy package odr (Brown et al. 1990, pp. 186;Virtanen et al. 2020). This ensures we account for minimisation in both variables during the fitting procedure.

Parameter Estimation:
Steps 1 -3 were performed many times, giving a distribution for the parameters a, b and c. For each parameter we report the 50 th percentile and uncertainties derived from the the 84 th and 16 th percentiles.
The relationship between L X, Obs. /L 12 µm and N H can be expressed as: We plot the Swift/BAT sample and Equation 1 (orange line) in Figure 5. To visualize the uncertainty in this relation, we show the 1000 realizations of the best fit as grey lines. We attempted to fit directly for log(N H ) as a function of log(L X, Obs. /L 12 µm ) but the fitting routine could not successfully converge on a reasonable fit to the data. Instead, we chose to invert Equation 1 to recover the obscuring column density as a function of the X-ray to mid-IR luminosity ratio: To better quantify the uncertainty in log(N H /cm −2 ) associated with this relation, we tabulated values of Relation for the column density vs. log(L X, Obs. /L12 µm) for the Swift/BAT sample (white circles). New data samples simulated during the MC fitting routine are displayed as blue dots. The best-fit trendline (orange line, given by Equation 1), represents the median of 1000 realizations (gray lines) obtained through orthogonal distance regression during the MC fitting process (see Section 3.1). The best-fit relation from Asmus et al. (2015) is shown as a red dashed line. log(N H ) and the uncertainty for specific values of log(L X, Obs. /L 12 µm ), as shown in Table 1. Using the parameter distributions calculated during the MC fitting routine, we derived a distribution of log(N H /cm −2 ) values using an expression of the form shown in Equation 2 (inverted from Equation 1) and use the 50 th , 84 th , and 16 th percentiles of the distribution to list the median column density and associated uncertainty for each chosen value of log(L X, Obs. /L 12 µm ). While strictly empirical, Equations 1 and 2 reproduce the observed trend of decreasing values of log(L X, Obs. /L 12 µm ) with increasing column density. We do note that this relation is largely insensitive to ratios of log(L X, Obs. /L 12 µm ) > −0.3, with the relation appearing nearly flat for column densities of log(N H /cm −2 ) < 22.5. Equations 1 and 2 are most effective for obscured AGNs with column densities of log(N H /cm −2 ) ≥ 22.5, and readers should keep this in mind when using these expressions to derive estimates for column densities.
Column 1: logarithmic ratio of the 2-10 keV to 12 µm luminosities; Equation 2 is insensitive to ratios of log(L X, Obs. /L12 µm) > −0.3, which are largely exhibited by AGNs with log(NH/cm −2 ) < 23. Column 2: column density and associated error, derived by Equation 2; the parameter distributions found during the MC fitting process described in Section 3.1 were read into an inverted expression of the form in Equation 2, from which we retrieved the median column density and upper and lower bounds using the 50 th , 84 th , and 16 th percentiles of the resulting log(NH/cm −2 ) distribution. Column 3: column density derived using Equation 6 from Asmus et al. (2015).

Relation between N H and the Mid-Infrared Colors
Given the correlation between L 22 µm /L 4.6 µm and N H , as shown in Figure 3, we followed the same fitting procedure as discussed above to develop an expression relating the mid-IR luminosity ratio and N H . The equation relating these properties can be expressed as: In Figure 6 we plot the Swift/BAT sample and the bestfitting trendline (orange line) along with 1000 realizations of the best fit found during the fitting procedure (shown as gray lines). As mentioned in Section 3.1, we attempted to fit directly for N H as a function of L X, Obs. /L 12 µm , but the fitting routine could not successfully converge on a reasonable line of best fit. We therefore simply invert Equation 3 like before to recover the column density as a Relation for column density vs. log(L22 µm/L4.6 µm) for the Swift/BAT sample (white circles). The fitting process for this relation is identical to that in Figure 5. The best-fit trendline (orange line, given by Equation 1) represents the median of 1000 realizations (gray lines) obtained during the MC fitting process. The black dashed line represents a theoretical prediction for a simple model consisting of monochromatic radiation passing through a homogeneous screen of dust (see Section 4.1.) function of the mid-IR luminosity ratio: In an identical fashion to the calculations in Section 3.1, we tabulated values of log(N H /cm −2 ) and the associated uncertainties for specific values of log(L 22 µm /L 4.6 µm ) in Table 2. Equation 4 is not sensitive to values of log(L 22 µm /L 4.6 µm ) < 0 and, as in the case of Equation 2, readers should use this relation cautiously and bear in mind that it is really only effective for obscured AGNs with log(N H /cm −2 ) ≥ 22.5.

Diagnostic Regions for Heavily Absorbed AGNs
The correlations between L X, Obs. /L 12 µm , L 22 µm /L 4.6 µm , and N H discussed above suggest that these relationships may be combined to help differentiate between Swift/BAT AGNs according to the levels of obscuration. We probed this potential diagnostic for heavily absorbed AGNs by plotting L 22 µm /L 4.6 µm vs. L X, Obs. /L 12 µm and binning the full sample by absorbing column as Column 1: logarithmic ratio of the 22 µm to 4.6 µm luminosities; Equation 2 is insensitive to ratios of log(L22 µm/L4.6 µm) > 0. Column 2: column density and associated error, derived by Equation 4; the median value and associated uncertainties were derived in an identical fashion to that described in Table 1 .
shown in Figure 7. The sample was divided into four bins in obscuration corresponding to: We also split the two Compton-thin bins into two subbins each in increments of log(N H /cm −2 ) = 0.5, as summarized in Table 3. Contours were computed for each bin and/or sub-bin and in each case are designed to encompass ∼ 68% of the population of each respective bin. The binned subplots shown in Figure 7 demonstrate much more clearly that, in general, the heavily obscured sources (Panel D) tend to occupy a separate region of space than the unobscured (Panel A) or Compton-thin 'lightly obscured' (Panel B) sources. While there is some overlap between the CT and Compton-thin 'moderately obscured' populations, this is predominantly due to Compton-thin 'moderately obscured' sources with 23.5 ≤ log(N H /cm −2 ) < 24. In light of this, appropriate selection criteria can be used to construct diagnostic regions in this parameter space which separate heavily obscured and less obscured AGN populations. Here we define the completeness of selection criteria as 'the fraction of true heavily obscured AGNs selected,' while we define purity as 'the fraction of selected AGNs which are heavily obscured,' i.e. the fractional contribution of heavily obscured AGNs to a diagnostic region. These definitions can be extended in an analogous fashion to the other column density bins.
Defining a horizontal cut in this parameter space: which is shown as a black dashed line in Figure 7, provides a simple yet robust method of differentiating between the most heavily obscured AGNs and less obscured AGNs in the Swift/BAT sample. We report the sample statistics for this cut in Table 3. To derive the population statistics for Table 3 (and all percentages quoted hereafter), we calculated the median (50 th percentile) value for each population in question, while the uncertainties on the fractions are the 16 th and 84 th quantiles of a binomial distribution, all computed following Cameron (2011). The criterion in Equation 5 yields 88.1 +4.5 −5.7 % completeness for the heavily obscured AGNs, a 60.5 +6.3 −6.5 % pure sample, and a mean column density of log(N H /cm −2 ) = 24.0 ± 0.1 for the diagnostic region. It is important to note that the majority of impurities selected with Equation 5 arise from AGNs with column densities 23.5 ≤ log(N H /cm −2 ) < 24.0; the diagnostic region is in fact ∼ 88% pure for AGNs with log(N H /cm −2 ) ≥ 23.5 and suffers minimal impurities from AGNs with lower column densities.
In a similar fashion, we could define a vertical cut based on the WISE colors in this space: and while this criterion also yields a highly complete sample (88 +4.5 −5.7 %) of heavily obscured AGN, the selected sample is only 13.1 +2.1 −2.0 % pure for heavily obscured AGNs and is greatly contaminated by moderately obscured, lightly obscured, and even unobscured AGNs. This is not at all surprising, given the scatter in the mid-IR ratio as demonstrated in Figures 3 and 5. Mid-IR selection alone is therefore not sufficient when attempting to select both a highly complete and fairly pure sample of heavily obscured sources.
Next we defined a slightly more stringent box region (gray dash-dotted line and black dotted line in Panels A-D of Figure 7), which encompasses the majority of the most heavily absorbed sources with minimal overlap with the unobscured and Compton-thin bins, using the following relations: Note-Statistics derived from the log(L Obs. X /L12µm) < −1.3 threshold, defined in Section 3 (Equation 5), for various NH bins and sub-bins. Column 1: NH bin. Column 2: Completeness, or the fraction of AGNs selected (per column density bin). Column 3: Purity of the sample, or the percentage contribution to the diagnostic box. Note-Statistics derived from the mid-IR log(L22 µm/L4.6 µm) ≥ 0.1 criteria (without invoking any cut in log[L Obs. X /L12µm]) defined in Section 3 (Equation 6), for various NH bins and sub-bins. Columns 1-3: The same as Table 3.
We report the population statistics for this diagnostic box in Table 5. This box offers a completeness of 83.0 +5.4 −6.4 % for the heavily obscured AGN population, a 62.4 +6.5 −6.8 % pure sample, and a mean column density of log(N H /cm −2 ) = 24.1 ± 0.1. As with Equation 5, the largest source of impurities within this region are AGNs with 23.5 ≤ log(N H /cm −2 ) < 24; the region is ∼ 90% Note-Statistics derived from the diagnostic box developed for the L22 µm/L4.6 µm luminosity ratio (defined by Equation 7 in Section 3) for various NH bins and sub-bins. Columns 1-3: The same as Table 3.
pure for AGNs with log(N H /cm −2 ) ≥ 23.5 and suffers few impurities from AGNs of lower column densities.
We repeated this analysis for the alternative L 12 µm /L 4.6 µm luminosity ratio, and we show the contoured populations binned by column density along with an alternative diagnostic box for heavily absorbed sources in Panels E-H of Figure 7. We construct this box with the following relations: and we find that this box yields a completeness of 80.5 +5.7 −6.7 % for heavily absorbed AGNs, a purity of 59.4 +6.5 −6.8 %, and a median column density of log(N H /cm −2 ) = 24.0 ± 0.1.
Again, AGNs with column densities 23.5 ≤ log(N H /cm −2 ) < 24 contribute the most to the impurity of the sample, while AGNs with lower column densities do not contribute as significantly.
The diagnostic metrics defined above -developed using the well-constrained X-ray and mid-IR properties previously found for Swift/BAT AGNs (e.g. Ricci et al. 2017a;Ichikawa et al. 2017) -carve out parameter spaces that yield fairly complete and fairly pure samples of heavily obscured AGNs, offering an efficient and effective method for identifying heavily obscured or CT AGN candidates, particularly in large samples of AGNs. The observed trend of the 22 µm/4.6 µm and 12 µm/4.6 µm WISE ratios increasing with column density can be readily understood from dust absorption and emission properties and basics of the radiation transfer. For media optically thin to the mid-IR radiation, the shape of the resulting SED will be determined predominantly by the dust temperature and its gradient throughout the dusty structure. If the dusty medium is optically thick to its own radiation, the outgoing emission will be reshaped due to a number of reasons: (i) Warm dust emission at shorter wavelengths will be absorbed and re-emitted at longer wavelengths. (ii) Dust emission at shorter wavelengths will suffer more extinction than the long-wavelength emission, owing to the wavelengthdependent extinction for typical AGN dust (Laor & Draine 1993). (iii) Warm dust emission originates closer to the inner rim of the torus, while colder emission originates farther out. As a consequence, longer wavelength mid-IR radiation has to travel a shorter path through the dust before reaching us and thus, suffers even less extinction than the emission of shorter wavelengths. These three effects result in an increased ratio of longerto-shorter wavelength mid-IR emission and scale with the column density of the medium through which the Xray radiation is traversing. (iv) Additionally, a disk-like molecular structure in hydrostatic equilibrium is expected to have a vertical gradient (Hönig 2019). In this case, the observed trend of increasing WISE ratios with increasing N H can be explained simply as an inclination effect: the closer our viewing angle is to the equator, the higher is the column density along the line-of-sight and, at the same time, the dust emission becomes "redder". All these effects contribute to the observed trend of increasing WISE ratios with N H and also explain why the effect is more pronounced in the 22 µm/4.6 µm than at 12 µm/4.6 µm luminosity ratio (for illustration, see torus model SEDs in, e.g., Hönig & Kishimoto 2010;Stalevski et al. 2012).
There are a few caveats: The dust emission is often degenerate, as SEDs of similar shape can be produced by different combinations of the geometrical and physical parameters of the torus, some of which can conspire to work against or hide the trend in luminosity ratios. However, the above reasoning should hold in general since it relies on universal radiative transfer effects. Another deviation can be introduced by the presence of silicate dust grains which exhibit a strong increase of absorption efficiency around 10 and 18 µm. The apparent strength of these features appear in an SED depends on several factors, including the amount of silicates, grain size distribution, and radiative transfer effects.
We illustrate these effects in Figure 6 with a black dashed line, which represents a theoretical expectation for a very simple model: monochromatic radiation pass- . The Compton-thin bins are broken into two sub-bins each: 22 ≤ log(NH/cm −2 ) < 22.5 (black points), 22.5 ≤ log(NH/cm −2 ) < 23 (blue points), 23 ≤ log(NH/cm −2 ) < 23.5 (black points), and 23.5 ≤ log(NH/cm −2 ) < 24 (orange points). Contours were computed to encompass ∼ 68% of the population for each respective bin or sub-bin. Note in general that the most obscured sources tend to populate a different region of the parameter space than do the unobscured sources and Compton-thin 'lightly obscured' sources, while there is some overlap with Compton-thin 'moderately obscured' sources.
ing through a homogeneous screen of dust. For this example, we assumed a typical Galactic interstellar dust mixture of silicates and graphite (e.g., Stalevski et al. 2016). The grain size distribution are from Mathis et al. (1977) and optical properties from Laor & Draine (1993) and Li & Draine (2001). The conversion between the optical depth and N H assumes Galactic relation between extinction and column density found by Predehl & Schmitt (1995). We see that the theoretical curve is following the trend of the data at lower column densities, but is reaching the breaking point sooner. This is because the simple dust screen model does not account for a number of radiative transfer effects (self-consistent absorption and re-emission of the thermal IR radiation), which together with geometry of the dusty medium and orientation shape the resulting SED, and thus, the observed trend of luminosity ratios with column density. Asmus et al. (2015) Using sub-arcsecond resolution mid-IR observations -which enabled the isolation of the nuclear mid-IR emission (F nuc 12 µm ) -Asmus et al. (2015) found a sig-nificant correlation between log(F nuc 12 µm / F obs 2−10 keV ) and log(N H /cm −2 ) for 53 AGNs with reliable X-ray observations and column densities log(N H /cm −2 ) > 22.8, which was expressed as (see Equation 6 in Asmus et al. 2015):

Comparison to
log N H 22.8 cm −2 = (0.14 ± 0.11) This expression is plotted in Figure 5 (dashed red line) along with our Equation 2 (solid orange line) and the Swift/BAT sample. Along with values of log(N H /cm −2 ) derived using Equation 2 in Section 3.1, we use the relation from Asmus et al. (2015) to derive values of log(N H /cm −2 ) and the uncertainties for specific values of log(L X, Obs. /L 12 µm ) in order to compare to our own results. Despite the fact that Asmus et al. (2015) removed AGNs with log(N H /cm −2 ) < 22.8 and utilized subarcsecond-resolution mid-IR emission (whereas in this work we utilized lower angular resolution mid-IR photometry for the Swift/BAT sample), it does appear that the two relations generally agree (within the uncertainties) for ratios of −0.75 log(L X, Obs. /L 12 µm ) −3.0. The two relations differ more severely for higher ratios of log(L X, Obs. /L 12 µm ), though this is expected due to (1) the wide range of ratios that unobscured AGNs exhibit and (2) the fact that our relation turns over to account for less obscured sources while the Asmus et al. (2015) relation does not take into account less obscured sources.

Diagnosing Column Densities with Uncertain Dust Heating Sources
While the diagnostic boxes defined in Section 3 provide a reliable way to identify the most heavily obscured AGNs, star formation activity can contribute non-negligibly to the mid-IR colors of an AGN host. The mid-IR colors assumed to originate from the AGN itself could therefore be overestimated without performing detailed spectral energy decomposition (SED) fitting to differentiate between the AGN and star formation contributions to the mid-IR continuum. Furthermore, Satyapal et al. (2018) demonstrated, using Cloudy (Ferland et al. 2013 radiative transfer models, that heavily obscured star formation activity can actually mimic the mid-IR colors of AGNs. These two points suggest that our diagnostic boxes defined in Section 3 may (1) misdiagnose the column density of an AGN if significant star formation is present, as the contaminating stellar emission could lead to much redder colors than the AGN intrinsically exhibits, or (2) mislead us to think an AGN is present in cases where the dominant dust heating sources are actually stellar-related rather than AGN-related . To investigate this potential contamination of the diagnostic boxes, we constructed a catalog of optically-selected galaxies whose optical spectroscopic line ratios suggest star formation dominates the observed emission, and we examined methods -for example mid-IR or X-ray selection criteria -through which this contamination could be mitigated.
Beginning with the MPA-JHU catalog of galaxy properties (from the SDSS data release 8, Aihara et al. 2011), we first selected systems with redshifts z < 0.1 and included only systems which are classified as star forming systems ("BPTClass" = 1) based upon their Baldwin, Phillips, Telervich (BPT; Baldwin et al. 1981) optical spectroscopic emission line ratios. We also removed any systems with QSO and AGN flags within the "TARGETTYPE," "SPECTROTYPE," and "SUB-CLASS" columns, and then narrowed the sample to only systems with WISE counterparts and X-ray counterparts from the 4XMM point source catalog (Webb et al.,  8) for the L22 µm/L4.6 µm (L12 µm/L4.6 µm) ratio with a population of star forming galaxies (see Section 4.3) overlaid as triangles. The horizontal dashed black line is given by Equation 5. The observed X-ray luminosity is denoted on the auxiliary axis. The diagnostics presented here cannot unambiguously differentiate between AGNs and star forming systems since a significant fraction of the star forming galaxy population falls within the absorption diagnostic box. This contamination can be mitigated with a mid-IR WISE cut of W 1 − W 2 > 0.8 (Stern et al. 2012); only six optically normal galaxies satisfy this mid-IR criteria (red open circles), and these systems fall outside of the diagnostic box defined by Equation 7 (Equation 8).
submitted). These criteria yielded a full parent sample of 448 galaxies which we assume are 'purely' star forming systems based upon optical spectroscopic measurements. We make no distinction between morphological classes of galaxies.
We plot our population of optically-selected star forming galaxies (color coded according to the observed X-ray luminosity) along with the L 22 µm /L 4.6 µm and L 12 µm /L 4.6 µm diagnostic boxes in Figure 8. While star formation dominated galaxies tend to exhibit lower ratios of log(L X, Obs. /L 12 µm ) than the majority of the Swift/BAT sample, they do tend to exhibit similar Xray deficits as well as L 22 µm /L 4.6 µm and L 12 µm /L 4.6 µm mid-IR colors as those exhibited by the heavily obscured Swift/BAT AGNs; in fact, 44.4±2.3 % of this star forming population overlaps the L 22 µm /L 4.6 µm diagnostic region, while 47.1 +2.3 −2.3 % of the population overlaps the L 12 µm /L 4.6 µm region. We therefore caution that this diagnostic is emphatically not designed to differentiate between star forming and AGN-dominated systems and should not be used as a diagnosis of the dominant dust heating source. Nevertheless, the contamination from optically-selected star forming systems can be mitigated through the use of reliable mid-IR and X-ray selection criteria traditionally used for identifying AGNs.
We applied the two band WISE AGN selection cut (W 1[3.4 µm]-W 2[4.6 µm] > 0.8) from Stern et al. (2012) to the sample of star forming galaxies, which removed all but six systems (red empty circles, shown in Figure 8). As expected, requiring traditional mid-IR AGN selection criteria eliminates virtually all contamination by optically-selected star forming systems within the diagnostic boxes. While six star forming galaxies (1.5 +0.6 −0.5 % of the total population) succeed in meeting the Stern et al. (2012) cut, these do still fall outside of our diagnostic boxes. 4 Thus, use of mid-IR AGN selection tools could be used to avoid misdiagnosing the dominant photoionization process of the sources within the diagnostic regions. However, it is important to bear in mind that the relation between L X, Obs. /L 12 µm , L 22 µm /L 4.6 µm (and L 12 µm /L 4.6 µm ), and N H holds true for both WISE and non-WISE selected AGNs (see Appendix), and therefore requiring a WISE cut could in general remove true AGNs as well as star forming systems. For example, imposing the W 1-W 2 > 0.8 cut on the Swift/BAT sample examined in this work would remove 240 systems (or 52.6 +2.3 −2.3 % of the parent sample of 456 AGNs); in the parent sample of 456, there are 71 AGNs which possess column densities in excess of 5 × 10 23 cm −2 , and 40 of these would be removed with this mid-IR cut.
We also found that requiring an observed X-ray luminosity of L X, Obs. > 10 42 erg s −1 removes nearly all of the optically star forming population from the diagnostic region (5 galaxies, or 1.3 +0.6 −0.5 %, remain within the L 22 µm /L 4.6 µm diagnostic box), though this method must also be used judiciously to avoid removing heavily obscured AGNs, which could exhibit lower X-ray luminosities.
Ideally, the usage of this diagnostic should be limited to systems whose dominant photoionization processes are unambiguous or for which detailed SED fitting can be performed to differentiate between AGN and host emission. Otherwise, we recommend proceeding cautiously, taking into account the various caveats outlined above to avoid inaccurate estimations of the obscuration along the line-of-sight.

Mid-IR Emission Contributions from Galaxies Hosting Obscured AGNs
The realization in Section 4.3 that optically star forming galaxies, which presumably do not host AGNs, can exhibit luminosity ratios similar to those exhibited by more heavily obscured or CT Swift/BAT AGNs raises the intriguing point of how much host galaxies may contribute to the observed luminosity ratios derived for the Swift/BAT AGNs. Figure 9 shows the log(L X, Obs. /L 12 µm ) and log(L 22 µm /L 4.6 µm ) ratios of the Swift/BAT AGNs and is color-coded according to the fractional contribution by the AGN to the observed 12 µm emission (f 12 µm AGN ) derived through detailed SED fitting in Ichikawa et al. (2019). While host-dominated systems at 12 µm (f 12 µm AGN < 0.5) can be found across this parameter space, a significant fraction (43.3 +6.9 −6.7 %, 22 out of 51) of AGNs within the diagnostic region (Equation 7) reside in host-dominated systems. This suggests that the host galaxies could contribute significantly to the observed mid-IR colors of the heavily obscured AGN population in particular, presumably via dust emission heated through star formation. Ichikawa et al. (2019) provided decomposed logarithmic AGN 12 µm luminosities for the Swift/BAT AGNs, which we can use here to examine how the log(L X, Obs. /L 12 µm ) ratios may change if we use the AGN 12 µm luminosity (L 12 µm, AGN ) instead of the total observed 12 µm luminosity (L 12 µm ). Figure 10 (top panel) shows that host-dominated systems are predominantly occupied by AGNs with log(N H /cm −2 ) 22.5; half of the CT Swift/BAT AGNs reside in host dominated systems. After recalculating the luminosity ratio using L 12 µm, AGN instead (bottom panel), hostdominated systems exhibit a shift toward higher luminosity ratios, with an average difference in ratio of ∆log(L X, Obs. /L 12 µm ) ≈ 0.5 for AGNs with f 12 µm AGN < 0.5, although we note that these shifts are not limited The logarithmic L X, Obs. /L12 µm vs. L22 µm/L4.6 µm ratios of the Swift/BAT AGN sample, where the data are color-coded according to the fractional contribution of the AGN 12 µm emission to the total observed 12 µm emission (f 12 µm AGN ). While there is no clear offset in this parameter space between AGN-dominated and host-dominated systems, several of the most heavily obscured AGNs (i.e. sources with significant X-ray deficits and log(L22 µm/L4.6 µm) > 0.1) do reside in systems where the host dominates the 12 µm emission.
only to heavily obscured AGNs. Here we have assumed that the X-ray emission is AGN-dominated, rather than host-dominated; in reality, if some portion of the X-ray emission is due to the host as well, the observed ratio shifts will not be as large. Figure 9 demonstrates that host galaxies do indeed contribute significantly to the diagnostic ratios probed in this work, at least for systems in which the host dominates the mid-IR emission at 12 µm. In these cases, the host contribution to the mid-IR leads to a perceived larger X-ray deficit for the AGN at a given column density. It does appear, though, that generally this effect actually works in our favor when attempting to identify CT AGNs, as these more severe X-ray deficits and presumably 'redder' mid-IR colors aid in separating this population from less obscured populations in color space. Decomposed AGN 22 µm and 4.6 µm luminosities were not included in Table 1 of Ichikawa et al. (2019) and therefore could not be examined in a similar fashion here. While it is beyond the scope of this paper, an analysis of the interplay between the mid-IR colors, host galaxy emission, and AGN emission with regard to the selection diagnostics presented in this work should be performed more rigorously in a future study. . The auxiliary axes represent the fractional contribution by the AGN to the total observed 12 µm emission. (Top) As is already known, log(L X, Obs. /L12 µm) decreases with increasing column density, however a significant number of the obscured and CT AGNs contribute less than 50% of the total observed 12 µm (f 12 µm AGN ). (Bottom) The relationship between log(L X, Obs. /L12 µm) and log(NH/cm −2 ) is still present when recalculating the ratio using the decomposed AGN 12 µm luminosity (L12 µm, AGN), although obscured and CT AGNs exhibit smaller deficits than when using the total 12 µm luminosity. Host galaxies can therefore contribute significantly to the observed X-ray to mid-IR ratios of AGNs, especially CT AGNs.
CT AGNs using mid-IR and far-IR photometry. They report, as we do here in Section 3 of this study, a shift in infrared colors (specifically mid-IR and far-IR) toward 'redder' colors with increasing column density, and they define a physically motivated color-color diagram (see Figure 11,  To further test their diagnostic, they applied this color selection criteria to the AKARI infrared galaxy catalog developed in Kilerci Eser & Goto (2018), which contains over 17,000 galaxies, and recover one known CT AGN (NGC 4418, e.g. Sakamoto et al. 2013). 5 The remainder of the infrared galaxy sample of Kilerci Eser & Goto (2018) is represented with blue contours in F 11 that partially overlap a significant number of Swift/BAT CT, obscured, and unobscured AGNs, suggesting that some portion of these infrared galaxies may in fact host heavily obscured AGNs. Despite finding a few cases of CT AGNs between the Swift/BAT and Kilerci Eser & Goto (2018) samples, the diagnostic criteria set forth in Kilerci-Eser et al. (2020) does not appear to provide a complete or reliable (see Section 5.4 of Kilerci-Eser et al. 2020) method of selecting CT AGNs.
As an additional comparison between our selection method and that proposed by Kilerci-Eser et al. (2020), we turned our attention to the infrared galaxy catalog from Kilerci Eser & Goto (2018). We matched this sample to the AllWISE catalog and the 4XMM DR9 XMM-Newton Serendipitous Source Catalog (Webb et al. submitted), using a match radius of 10 for each, which yielded a sample of 401 local (z < 0.1) infrared galaxies with XMM-Newton and WISE detections. In Figure 11, we show the resulting sample of infrared galaxies (red stars), along with our diagnostic criteria from Equations 5 and 7. As in Section 4.3, it is impossible to discern whether or not any of these galaxies host AGNs without a reliable method for removing star formation dominated systems. We tried four different mid-IR AGN selection criteria, defined in Jarrett et al. (2011), Stern et al. (2012, Assef et al. (2018), and Satyapal et al. 5 The selection criteria also recovers two other sources: NGC 7714, an unobscured AGN (Gonzalez-Delgado et al. 1995;Smith et al. 2005), and NGC 1614, which has no clear evidence of an AGN (e.g. Xu et al. 2015;Pereira-Santaella et al. 2011;Herrero-Illana et al. 2014).  Table 7 for more details on these AGNs.
(2018) (with the understanding that some heavily obscured AGNs will be missed with this simple approach), and in all four cases we recover a significant number of candidate heavily obscured or CT AGNs. We show mid-IR AGNs selected as a result of the Stern et al. (2012) cut in Figure 11 (blue squares); these candidate CT AGNs likely inhabited the blue contoured prominence that overlapped the Swift/BAT AGNs in F 11, but were missed due to fact that they did not satisfy the criteria proposed by Kilerci-Eser et al. (2020). Table 7 in the Appendix lists the 36 candidate CT AGNs selected from Kilerci Eser & Goto (2018) using the L 22 µm /L 4.6 µm diagnostic region and at least one of the mid-IR selection cuts listed above. We include in the table the source coordinates, redshifts, luminosity ratios, and alternative identifiers, and we also categorize the column densities of the sources (using the column density bins defined in Section 3.3) based upon any available measurements in the literature. Of the candidates that have inferred or directly measured column densities in the literature (24/36), we find eight CT AGNs, six moderately obscured AGNs, two lightly obscured AGNs, and three unobscured AGN, while a remaining five AGNs have conflicting measurements of N H in the literature (all five of which have been reported as CT at least once in the past). 6 Therefore, we conclude that our diagnostic criteria proposed in Equations 5, 7, and 8 offer a more reliable method for identifying candidate CT AGNs than the mid-IR to far-IR color criteria proposed by Kilerci-Eser et al. (2020).
In light of Section 4.4, we caution that some portion of these AGNs may not dominate the observed 12 µm emission and, therefore, may exhibit larger X-ray deficits and redder mid-IR colors than might be expected for the AGN alone due to additional mid-IR contributions from the host galaxy.

Diagnosis of N H in the XXM-XXL Field
To test the power of our absorption diagnostic, we turn our attention to its application in the XMM XXL North field (Pierre et al. 2016(Pierre et al. , 2017. Menzel et al. (2016) presented a rigorous multiwavelength analysis of 8445 X-ray sources detected by XMM-Newton in an 18 deg 2 area of the XMM XXL North field (hereafter XXL-N), with a limiting flux of F 0.5−10 keV > 10 −15 erg cm −2 s −1 , providing optical Sloan Digital Sky Survey (SDSS) and mid-IR WISE counterparts to the XMM-Newton sources. In a complementary investigation, Liu et al. (2016) presented a thorough X-ray spectral analysis of the 2512 XXL-N AGNs, deriving the spectral properties (e.g. photon index Γ, N H , L X, Obs. ) for those AGNs using a Bayesian statistical approach contained within the Bayesian X-ray Astronomy (BXA) software package (Buchner et al. 2014).
We combined the catalogs from Menzel et al. (2016) and Liu et al. (2016) to obtain the observed 2-10 keV luminosities and mid-IR WISE magnitudes, from which we derived the relevant luminosity ratios examined in Section 3. Initially we limited the XXL-N sample to only local (z < 0.1) AGNs, and as with the Swift/BAT AGNs, we did not employ any mid-IR or X-ray selection criteria. We plot the resulting luminosity ratios of the low redshift AGNs from the XXL-N field in the left panel of Figure 12 along with our L 22 µm /L 4.6 µm diagnostic box (Equation 7) and horizontal cut in L X, Obs. /L 12 µm (Equation 5). The data and auxiliary axis in Figure 12 are color coded to represent the derived 50 th -percentile N H values from the Liu et al. (2016) catalog and the markers denote the obscuration bin for each AGN. In the right panel of Figure 12 we plot the L X, Obs. /L 12 µm 6 There are other moderately and heavily obscured WISE AGNs in this sample which did exhibit log(L X /L 12 µm )< −1.3 but fall outside of our more stringent diagnostic region. ratio against the derived log(N H /cm −2 ) values from Liu et al. (2016), where the error bars represent the 16 th and 84 th percentiles and the data points use the same marker and color scheme as the left panel. A dearth of low redshift AGNs in the XXL-N field is immediately apparent, and while it appears that most of the more obscured AGNs do exhibit "redder" colors, there is some overlap in L X, Obs. /L 12 µm ratios exhibited by AGNs of starkly different obscuration bins, e.g. heavily obscured and unobscured AGNs. Unfortunately, we cannot draw definitive conclusions about the reliability of our diagnostic for the XXL-N field with such poor AGN statistics.
We repeated this analysis for the entire sample of XXL-N AGNs, breaking the sample into redshift bins of ∆z = 0.5 each. While a small number of obscured AGNs overlap with the diagnostic region defined by Equation 7, the majority of the heavily absorbed AGNs still occupy the same parameter space as the Compton-thin and unobscured AGNs, even in the local (z < 0.5) redshift bin, in stark contrast to the results found with the Swift/BAT AGNs. We find a very similar result when examining the L 12 µm /L 4.6 µm diagnostic ratio. Due to the fact that at higher redshift the L 22 µm /L 4.6 µm and L 12 µm /L 4.6 µm diagnostics do not correspond to the same wavelength ranges as they do at local redshifts, we then turned to the luminosity ratio of L 12 µm /L 3.4 µm . Yet again the heavily absorbed AGNs generally occupy the same parameter space as unobscured AGNs. The spectral curvature method (Koss et al. 2016) may provide a more effective means of selecting heavily obscured AGNs at higher redshift in the XXL-N field, and indeed Baronchelli et al. (2017) demonstrated its effectiveness in selecting high redshift (z > 2) CT AGNs in both the Chandra Deep Field South and the Chandra COSMOS legacy survey.
There are a few explanations for why the luminosity ratios of the XXL-N AGNs do not as clearly differentiate between obscuration levels like that seen for the Swift/BAT AGNs. For one, there may simply not be enough local redshift XXL-N AGNs for a proper statistical comparison to the results found for the Swift/BAT AGNs. Secondly, it is possible that the chosen redshift bins are not fine enough and we are including AGNs across too large of redshift ranges. After splitting the z < 0.5 bin into five sub-bins, however, we find the same result as before: the unobscured and heavily obscured AGNs coexist within the parameter space.
Another explanation lies in the reliability of the results of the X-ray spectral fitting performed in Liu et al. (2016): while the Bayesian statistical framework employed in that work is a powerful method for constrain-  Liu et al. (2016), and the data points are color-coded according to the auxiliary axis to the left. To aid the reader, we use different markers to denote the obscuration bin for each source. Unlike the Swift/BAT AGNs, there is not as clear of relationship between the luminosity ratios and column density for the XXL-N AGNs, though this comparison should be viewed with caution because of the large uncertainties on log(NH/cm −2 ) and because these AGNs were selected with a softer energy band (0.3-10 keV).
ing the spectral properties for AGNs with low counts, our results suggest that the column densities derived for the AGNs may still be inaccurate due to the low counts acquired. For example, for all AGNs detected in the z < 0.5 bin, the median number of counts detected in the epic pn and in the two epic mos detectors is only 16.8 and 15.1 counts, respectively. Furthermore, identification of Compton-thick AGNs via XMM-Newton spectroscopy is quite difficult due to the softer X-ray passband (0.3-10 keV) probed with the XMM-Newton imaging. It can be very difficult to distinguish between the scenario in which the source is heavily obscured, and the X-ray emission is dominated by reprocessed radiation from the circumnuclear material, and that in which the source is unobscured but the emission is dominated by relativistic reflection from the accretion disk when using low signal-to-noise X-ray spectra alone due to the strong model-dependent degeneracies involved (see e.g., Gandhi et al. 2009, Treister et al. 2009). Our mid-IR ratio predictions for the highest Xray column densities can hence be very complementary to help classify sources in which unobscured and obscured reflection models can fit an observed X-ray spectrum equally well.
Future, deeper XMM-Newton or NuSTAR follow-up observations could improve upon the current photon statistics and could provide robust constraints on the column densities for at least the local XXL-N AGNs which lie within the diagnostic boxes. For now, we identify these systems as candidate heavily obscured or CT AGNs, and we list their spectral properties in Table 8.

CONCLUSIONS
Using the well-studied Swift/BAT sample of ultra hard X-ray selected AGNs, we have presented an analysis of the AGN L X, Obs. /L 12 µm luminosity ratio and two different mid-IR luminosity ratios: L 22 µm /L 4.6 µm and L 12 µm /L 4.6 µm . Using the well constrained X-ray (Ricci et al. 2017a) and mid-IR (Ichikawa et al. 2017) properties of the Swift/BAT AGNs, we probed the utility of these luminosity ratios as tools for inferring line-ofsight column densities and identifying the most heavily obscured AGNs in the local Universe (z < 0.1). We summarize the results of our analysis as follows: • We have derived expressions relating the column density N H to both the L X, Obs. /L 12 µm and L 22 µm /L 4.6 µm luminosity ratios, which are defined by Equations 2 and 4. These expressions can be inverted to give N H as a function of the luminosity ratios. We provide these expressions again here: • We have demonstrated that unobscured and heavily obscured AGNs tend to exhibit different L X, Obs. /L 12 µm , L 22 µm /L 4.6 µm , and L 12 µm /L 4.6 µm luminosity ratios. All three of our diagnostic regions (Equations 5, 7, and 8) identify (in general) the most heavily absorbed AGNs, with average column densities of log(N H /cm −2 ) ≥ 24.0 for each defined parameter space. These regions are all 80% complete and 60% pure for AGNs with log(N H /cm −2 ) ≥ 24. The greatest impurities arise due to AGNs with 23.5 log(N H /cm −2 ) < 24; these regions are 85% pure for AGNs with log(N H /cm −2 ) 23.5.
• While optically star forming systems can fall within our diagnostic regions, this contamination can be virtually eliminated via mid-IR or X-ray selection criteria. Such selection criteria should be used judiciously to avoid removing non-mid-IR AGNs. These diagnostic regions should not be used to differentiate between AGNs and galaxies dominated by star formation.
• Swift/BAT AGNs which do not dominate the total observed 12 µm emission tend to exhibit redder colors and larger X-ray deficits with increasing column density, suggesting that host galaxy contributions to at least the mid-IR emission can be a significant factor in the luminosity ratios examined here, particularly in the case of mildly obscured and CT AGNs selected with our diagnostics. However, it appears that this effect actually aids in the identification of CT AGNs, as the host contributions result in a larger separation in color space between less obscured and more obscured AGNs.
• We find that the selection criteria proposed here are more reliable at identifying obscured and CT AGNs than the mid-IR and far-IR selection criteria proposed by Kilerci-Eser et al. (2020). We identify several known obscured and CT AGNs, as well as several candidate CT AGNs, within the IR galaxy catalog of Kilerci Eser & Goto (2018) (see Table 7).
• We applied our diagnostics to the XMM-Newton XXL-N field and found, in contrast to Swift/BAT AGNs, that obscured and unobscured XXL-N AGNs do not appear to exhibit distinctly different luminosity ratios. This disparity could be due to poor photon statistics or the softer X-ray energies probed for the XXL-N AGNs, which could lead to inaccurate column density values. Although, given the small number of z < 0.1 XXL-N AGNs and the large errors associated with several N H values, this comparison should be viewed with caution.
Identifying heavily obscured AGNs remains an important yet difficult task, though the study of such AGNs is an important step in the development of our understanding of the evolution of AGNs. In a future study, we could expand our analysis to include diagnostic regions appropriately modified to differentiate between unobscured and heavily obscured sources at higher redshift, although it is difficult to speculate at the moment how the emission ratios may change with redshift, as both star formation activity and AGN activity are expected to increase with redshift.
The selection criteria presented here offers a complementary approach to the spectral curvature method developed in Koss et al. (2016), which is very effective at selecting heavily obscured AGNs at local-z, with the caveats that one must already have hard X-ray (> 10 keV) measurements with NuSTAR or Swift/BAT and that it is most effective for brighter AGNs. Softer X-ray missions can only be utilized for higher redshift sources (z ∼ 3) for which the hard X-ray emission has shifted into the rest frame 10-30 keV passband. The diagnostics presented here, on the other hand, do not require higher energy passbands in order to select local-z sources and can take advantage of softer X-ray missions such as Chandra and XMM-Newton. The synergy between these two approaches is best summed up by the fact that they select many of the same sources using different passbands, and that they select heavily obscured AGNs missed by one another (see Figure 15 in the Appendix), yielding a more complete census of heavily obscured Swift/BAT AGNs overall.
The diagnostic regions proposed in this study, as well as the expressions derived relating N H to the L X, Obs. /L 12 µm and L 22 µm /L 4.6 µm luminosity ratios, could be used to differentiate between unobscured and heavily obscured AGNs in future, large samples of AGNs, such as those now being detected by the eROSITA all-sky survey (Predehl et al. 2010;Merloni et al. 2012). In particular, the eROSITA survey will provide the first all-sky X-ray imaging survey at energies up to 10 keV, yielding a highly complementary catalog to those of other all-sky missions, such as WISE . Future works could cross-match the WISE and eROSITA catalogs and use the diagnostics presented here to identify many more cases of CT AGN candidates, select targets for deeper follow-up multiwavelength observations, and to compute the CT fraction for the future sample, all of which will be crucial in the quest to construct a Software: odr (Brown et al. 1990;Virtanen et al. 2020), pandas (McKinney 2010), scipy (Virtanen et al. 2020), numpy (Oliphant 2006;van der Walt et al. 2011;Oliphant 2015), matplotlib Hunter (2007), BXA (Buchner et al. 2014) .
This publication makes 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, funded by the National Aeronautics and Space Administration. Funding for SDSS-III has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, and the U.S. Department of Energy Office of Science. The SDSS-III web site is http://www.sdss3.org/.
This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. In the process of our analysis, we explored whether any quality cuts or selection cuts could be applied to the data to reduce the scatter observed in Panel A of both Figures 3 and 4. Our parent sample was divided into four sub-samples as shown in Figure 13: • AGN-dominated systems with W 1-W 2 > 0.8 (Stern et al. 2012, blue squares).
• AGNs with > 300 spectral counts in the X-ray spectra, which provides a statistically significant number of counts to constrain N H in the X-ray spectral fitting analysis (Ricci et al. 2017a, inverted cyan triangles).
In Figures 13 and 14 we compare these sub-samples for the L 22 µm /L 4.6 µm and L 12 µm /L 4.6 µm luminosity ratios, respectively, and how they correlate with N H . The resulting mean values per bin for each different sub-sample are consistent with the values (the error bars represent the standard deviation of each subsample) originally found for the parent sample; interestingly, the observed correlation between WISE color and column density holds for AGNs that satisfy the Stern et al. (2012) criterion as well as AGNs which do not satisfy that criterion. We observe a larger difference in the L 22 µm /L 4.6 µm ratios when moving to higher column densities than for the L 12 µm /L 4.6 µm ratios. Additionally, as shown previously in Figure 4, while there is a large amount of scatter in the lowest-N H bin of the bottom panel of Figure 14, this scatter is likely due to the low number of sources within that bin, and this still does not overlap the bin probing the highest obscuring columns. Due to the consistency between the results for the sub-samples and that found for the parent sample, we do not implement any of these cuts during our analysis. From Figure 13 and 14, it becomes clear that there a number of WISE AGNs (W1-W2 > 0.8) within the Swift/BAT sample which exhibit extremely red mid-IR colors, with log(L 22 µm /L 4.6 µm ) > 0.5, suggesting significant obscuration. We tabulate these WISE AGNs in Table 9. Indeed, the majority of these AGNs (18/23) are moderately to heavily obscured with log(N H /cm −2 ) > 23, though there are a few exceptions, notably: • HS0328+0528, an unobscured Seyfert 1.
WISE selection based on the cut defined in Stern et al. (2012) is not, however, a necessarily good method for selecting obscured over unobscured AGNs. As we discussed in Section 4.3, 40/71 of AGNs with log(N H /cm −2 ) > 5.0 × 10 23 cm −2 in the Swift/BAT sample studied here do not meet a color cut of W1-W2 > 0.8, reinforcing our choice to not invoke such a cut on the parent sample.  Figure 13. We applied four different cuts to our full Swift/BAT sample to explore the origin of the scatter observed in Figure 3. Here we present two different comparisons of the derived values for NH from Ricci et al. (2017a) and the L22 µm/L4.6 µm mid-IR ratios, where we have (top panel) binned by log(NH) and (bottom panel) binned by the L22 µm/L4.6 µm ratio. We observe more scatter when binning by the mid-IR ratio than we do when binning instead by NH, but nonetheless the general trend is the same: we observe increasing mid-IR ratios of L22 µm/L4.6 µm with increasing column density regardless of the sub-sample.  Figure 13 except here we examine the alternative mid-IR diagnostic ratio which depends upon L12 µm and L4.6 µm. We observe an increase in the mid-IR ratio of L12 µm/L4.6 µm with column density as was observed with L22 µm/L4.6 µm, although with a large amount of scatter in the NH values for the lowest mid-IR ratio bin. However, given that our focus is on the most heavily obscured sources (> a few times 10 23 cm −2 ) this scatter is not a concern.   Note-All eight low-redshift (z < 0.1) candidate heavily obscured or CT AGNs within the XMM -XXL-N field pulled from the Menzel et al. (2016) and Liu et al. (2016) catalogs selected using the diagnostic box defined by Equation 7. Column 1: unique X-ray identification string used by both Menzel et al. (2016) and Liu et al. (2016). Column 2-3: X-ray source right ascension (RA) and declination (Dec) values (uncorrected for systematic offsets) given in degrees. Column 3: redshift. Column 4: Hard X-ray 2-10 keV fluxes. Column 5: 50 th -percentile logarithmic line-of-sight column density derived by Liu et al. (2016). Error bounds are calculated using the 84 th and 16 th -percentile column density values. Column 6-7: values for the log(L X, Obs. /L12 µm) and log(L22 µm/L4.6 µm) luminosity ratios. Note-WISE -selected Swift/BAT AGN (W 1−W 2 > 0.8, Stern et al. 2012) which exhibit mid-IR ratios of log(L22 µm/L4.6 µm) > 0.5. Column 1: SWIFT I.D. number. Column 2-4: Right ascension, declination, and redshift. Column 5-6: log(L X, Obs. /L12 µm) and log(L22 µm/L4.6 µm) luminosity ratios. Column 7: Alternative I.D. drawn from Ricci et al. (2017a). Column 8: Column density derived in Ricci et al. (2017a). CT (Koss et al. 2016) AGNs selected via Eq. 4 Figure 15. A comparison between the X-ray and mid-IR selection criteria introduced in this work and the spectral curvature method from Koss et al. (2016). The spectral curvature of each point was calculated using observations from NuSTAR, whereas the column densities come from Ricci et al. (2017a). Only AGNs with spectral curvature errors of < 0.2 are included in this plot. The dashed grey line represents a spectral curvature value of 0.4, above which an AGN is considered to be CT. Red points are AGNs selected using the diagnostic box defined in Equation 7, whereas blue points are not selected via Equation 7. The spectral curvature method and the diagnostic defined in this work find several of the same sources, and in fact these two approaches also find CT AGNs missed by one another: a few CT AGNs fall below the 0.4 spectral curvature cutoff, yet the diagnostic presented here selects them, meanwhile the spectral curvature method recovered a CT AGN missed by our new diagnostic. There is one outlier which is selected by Equation 7 and exhibits a relatively high spectral curvature, yet it possesses a column density of only ∼ 10 22 cm −2 . This source, NGC 1365, is a well known variable absorber and has gone through massive absorption transitions. Recent high signalto-noise observations have shown that the column density remains substantial, above 10 23 cm −2 (e.g. Risaliti et al. 2009a,b,c;Maiolino et al. 2010;Walton et al. 2010;Brenneman et al. 2013) and occasionally increases to the extent of becoming CT (Risaliti et al. 2005). It is possible that the mid-IR emission is tracing a higher absorption period seen in past observations.