Torus model properties of an ultra-hard X-ray selected sample of Seyfert galaxies

We characterize for the first time the torus properties of an ultra-hard X-ray (14-195 keV) volume-limited (DL<40 Mpc) sample of 24 Seyfert (Sy) galaxies (BCS40 sample). The sample was selected from the Swift/BAT nine month catalog. We use high angular resolution nuclear infrared (IR) photometry and N-band spectroscopy, the CLUMPY torus models and a Bayesian tool to characterize the properties of the nuclear dust. In the case of the Sy1s we estimate the accretion disk contribution to the subarcsecond resolution nuclear IR SEDs (~0.4'') which is, on average, 46+-28, 23+-13 and 11+-5% in the J-, H- and K-bands, respectively. This indicates that the accretion disk templates that assume a steep fall for longer wavelengths than 1 micron might underestimate its contribution to the near-IR emission. Using both optical (broad vs narrow lines) and X-ray (unabsorbed vs absorbed) classifications, we compare the global posterior distribution of the torus model parameters. We confirm that Sy2s have larger values of the torus covering factor (CT~0.95) than Sy1s (CT~0.65) in our volume-limited Seyfert sample. These findings are independent of whether we use an optical or X-ray classification. We find that the torus covering factor remains essentially constant within the errors in our luminosity range and there is no clear dependence with the Eddington ratio. Finally, we find tentative evidence that even an ultra hard X-ray selection is missing a significant fraction of highly absorbed type 2 sources with very high covering factor tori.


INTRODUCTION
Active galactic nuclei (AGNs) are powered by accretion of material onto supermassive black holes (SMBHs), which release energy in the form of radiation and/or mechanical outflows to the host galaxy interstellar medium. Although they comprise just a small fraction of the galaxies in the local universe (∼10%), AGNs are now considered to be a short but recurrent phase in the overall lifetime of galaxies. Accordingly, galaxies are observed as AGN during an active phase when their SMBHs are accreting material at a relatively high rate (e.g. Bennert et al. 2011). Several studies found a correlation between the SMBH and host galaxy bulge mass (e.g. Kormendy & Ho 2013 and references therein) which is interpreted as a sign of co-evolution of AGNs and their host galaxies. However, the study of the AGN-host galaxy connection is difficult due to the very different spatial and temporal scales involved. Therefore it is of great importance to investigate the innermost regions of AGN to better understand this connection (see Ramos Almeida & Ricci 2017 and references therein).
The key piece of the AGN unified model (Antonucci 1993) is a dusty molecular torus that obscures the central engines of type 2, and allows a direct view in the case of type 1 sources. This dusty torus absorbs part of the AGN radiation and reprocesses it to emerge in the infrared (IR).
To correctly separate the nuclear emission from the foreground galaxy emission and be able to characterize the properties of the nuclear obscurer the highest possible spatial resolution is required. Since Seyfert (Sy) galaxies are intermediate-luminosity AGNs, and, in general, are relatively nearby, they are one of the best astrophysical laboratories to study the inner regions of active galaxies.
The torus radius has been constrained to be compact (∼0.1-10 pc) in the mid-IR (MIR; ∼5-30 µm). For example, using MIR direct imaging,  and Radomski et al. (2008) found for Circinus and Centaurus A that the MIR size of the torus is less than ∼4 pc (diameter). The modelling of MIR interferometric data shows a relatively compact torus of r< 10 pc (e.g. Jaffe et al. 2004;Tristram et al. 2007;Tristram et al. 2009;Burtscher et al. 2009;Raban et al. 2009;Burtscher et al. 2013;López-Gonzaga et al. 2016). Atacama Large Millimeter/submillimeter Array (ALMA) observations of the archetypal Seyfert 2 galaxy NGC 1068 have spatially resolved for the first time the submillimeter (sub-mm) counterpart of the putative torus (García-Burillo et al. 2016;Gallimore et al. 2016;Imanishi et al. 2018). This is a disk of ∼7-10 pc diameter. More recently, Alonso-Herrero et al. (2018) and Combes et al. (2018) have found even larger nuclear molecular disks for other Seyfert galaxies and lowluminosity AGNs. Thus, as theoretically predicted (e.g. Schartmann et al. 2008;Stalevski et al. 2012), the radii measured in the sub-mm for the dusty and molecular torus are found to be larger than those inferred from IR observations. Therefore, to constrain the properties of the warm dust, we still need to compare torus models to the observed SEDs.
Torus models can be broadly grouped in two categories: physical (e.g. Wada & Norman 2002;Schartmann et al. 2008;Wada 2012) and geometrical (ad-hoc; e.g. Pier & Krolik 1992;Efstathiou & Rowan-Robinson 1995;Nenkova et al. 2008a,b;Hönig et al. 2010;Stalevski et al. 2012;Siebenmorgen, Heymann, & Efstathiou 2015;Hönig & Kishimoto 2017). Some of the geometrical models also include a polar component in the MIR range (e.g. Hönig & Kishimoto 2017). However, this polar emission has been detected so far in six Seyfert galaxies of the 23 observed using IR interferometry (López-Gonzaga et al. 2016;Leftley et al. 2018) and, therefore, more observations are needed in order to study whether this is a common feature in AGNs. The physical models are more realistic since they include important processes, such as supernovae and AGN feedback. However, they require large computational times and therefore it is more difficult to compare with observations. On the other hand, geometrical torus models are more degenerate, but they can be easily compared with the observations, assuming various geometries and compositions of the dust (see Ramos Almeida & Ricci 2017 for a review).
In this work we use the Nenkova et al. (2008a,b) clumpy torus models, known as CLUMPY, and the Bayesian tool BayesClumpy (Asensio  to fit the nuclear IR emission of an ultra-hard X-ray selected sample of Seyfert galaxies. Our aim is to study the torus properties that are driving the Seyfert type classification, the difference in the dusty torus of the various Seyfert types and how they vary with the central engine properties. The paper is organized as follows. Section 2 and 3 describe the sample selection, the observations and data compilation, respectively. The nuclear IR SED construction and modelling are presented in Section 4. In Section 5, we compare the torus properties for the different Seyfert subgroups. Finally, in Section 6 we summarize the main conclusions of this work. Throughout this paper, we assumed a cosmology with H0=73 km s −1 Mpc −1 , Ωm=0.27, and ΩΛ=0.73, and a velocity-field corrected using the Mould et al. (2000) model, which includes the influence of the Virgo cluster, the Great Attractor, and the Shapley supercluster.
We selected all the Seyfert galaxies in the nine month catalog with luminosity distances DL <40 Mpc. We used this distance limit to ensure a resolution element of 50 pc in the MIR, considering the average angular resolution of 8-10 mclass ground-based telescopes (∼0.3 ′′ at 10 µm). This constraint provides us with a sample of 24 local Seyfert galaxies (hereafter BCS40 sample; GB16) containing 8 Sy1 (Sy1, Sy1.2 and Sy1.5), 6 Sy1.8/1.9 and 10 Sy2 galaxies. This sample covers an AGN luminosity range log(L 2−10 keV int )∼40.5-43.4 erg s −1 . See GB16 for further details on the sample selection. The properties of the BCS40 sample are shown in Table 1.

OBSERVATIONS
Our aim is to construct high angular resolution IR SEDs for the whole sample. In the following we describe the new and archival MIR and near-IR (NIR; ∼1-5 µm) observations used in this work.

Gran Telescopio CANARIAS/CanariCam
We obtained subarcsecond resolution N-band spectra (7.5-13 µm) of two Seyfert galaxies (NGC 4138 and UGC 6728) using the low spectral resolution (R∼175) grating available in the instrument CanariCam (CC; Telesco et al. 2003;), on the 10.4 m Gran Telescopio CA-NARIAS (GTC). CC is a MIR (7.5-25 µm) imager with spectroscopic, coronagraphic and polarimetric capabilities. It uses a Si:As detector, which covers a field of view (FOV) of 26×19 arcsec 2 on the sky and it has a pixel scale of 0.08 arcsec. NGC 4138 and UGC 6728 were observed in 2016 March and the slit, of width 0.52 arcsec, was oriented at PA= 145 and 150 degrees, respectively. The total on-source integration times were 1061 and 1415 s, respectively. In both cases, the standard MIR chopping-nodding technique was used with chop and nod throws of 15 arcsec (see Table 2). The data were taken on 2016 March 14 and 15 as part of a Director's Discretionary Time program (GTC04-15B DDT; PI: I. García-Bernete). Using the acquisition images of the standard stars used for NGC 4138 (HD 95121) and UGC 6728 (HD 105943), we measured for the standard stars full width at half-maximum (FWHM) values of 0.28 ′′ (at λ=10.3 µm) and 0.34 ′′ (at λ=8.7 µm), respectively.

Gran Telescopio CANARIAS/CIRCE
We obtained NIR imaging data (J, H & K bands) with the Canarias InfraRed Camera Experiment (CIRCE; Garner et al. 2014) on the 10.4-m GTC. The instrument was equipped with an engineering grade Hawaii2RG detector with a total FOV of 3.4×3.4 arcmin 2 and a plate scale of 0.1 arcsec pixel −1 . Note that all the observations were taken using a 5 dither pattern. See Table 3 for observation details.
We performed the data reduction by using the IDL (Interactive Data Language) routines employed in D' Ammando et al. (2017). The first step in the data processing includes the subtraction of dark current frames. From twilight sky exposures, we obtained an illumination correction to compensate a decrease of about 40 per cent from the centre to the border of the FOV. At this point, we introduced a correction to remove a pattern of inclined stripes related to reading amplifiers. Once this pattern was removed, the images corresponding to each dither cycle were median combined to form a sky frame, which was subtracted for each frame of the cycle. We then combined all sky-subtracted images with the commonly used shift-and-add technique. During the combination of these frames, we applied a bad-pixel mask, which includes the two vertical bands corresponding to non-functional amplifiers. Finally, we obtained the photometric calibrations relative to photometric standard PSF stars using their Two-Micron All-Sky Survey (2MASS) photometry.
To estimate the NIR nuclear fluxes in the J, H & K bands we used the PSF subtraction method (see GB16 and references therein), which consists of subtracting the PSF star from the galaxy profiles. This method has been widely used in ground-based IR images (e.g. Soifer et al. 2000;Radomski et al. 2002Radomski et al. , 2003García-Bernete et al. 2015).

Archival Data
We downloaded the fully reduced NIR imaging data of NGC 4138 (unpublished, to our knowledge) from the ESA Hubble Legacy Archive 1 . This Seyfert 1.9 galaxy was observed in February 2008 with the Near Infrared Camera and Multi-Object Spectrometer (NICMOS) and the narrow F190N filter (λc=1.9 µm). This observation was taken using the NIC3 camera, which has a FOV 51.2×51.2 arcsec 2 on the sky and a pixel scale of 0.2 ′′ . This image was taken 1 http://archives.esac.esa.int/hst/ as part of the Hubble programs GO11080 (cycle:15, PI: D. Calzetti) and the exposure time was 13474 s. In order to accurately subtract the unresolved AGN component, first, we generated a theoretical Tiny Tim 2 PSF (Krist 1995;Krist, Hook, & Stoehr 2011) for the NIC3 camera F190N filter and, then, we used the PSF subtraction method.

Literature High Angular Resolution IR Observations
We compiled the highest angular resolution IR (∼1-30 µm) nuclear fluxes available from the literature for our sample. The compiled nuclear NIR fluxes are from both groundand space-based (i.e. Hubble Space Telescope; HST) data (see Table 2). In the case of the MIR nuclear fluxes, we used the measurements of the unresolved MIR emission (angular resolutions ranging from 0.2 ′′ to 0.6 ′′ ) calculated in GB16, where we employed the PSF subtraction method on high angular resolution MIR images from 8-10 m-class ground-based telescopes (GTC/CanariCam, VLT/VISIR, Gemini/T-ReCS and MICHELLE; see Table 2 of GB16). We retrieved 31.5 µm high angular resolution (FWHM∼3.1 ′′ ) nuclear fluxes of six Seyfert galaxies (see Table 4), which were observed with the long-wavelength camera (LWC; λ > 25 µm) within the Faint Object Infrared Camera for the SOFIA Telescope (FORCAST; Herter et al. 2012) on the 2.5 m SOFIA telescope. These observations were obtained using the 31.5 µm filter (∆λ= 5.7 µm). See Fuller et al. (2016) for further details on the observations, data reduction and obtention of unresolved nuclear fluxes.
Finally, we compiled N-band spectra (7.5-13 µm) for the majority of the sample (17/24 sources), which were obtained with different instruments (GTC/CC, VLT/VISIR, Gemini/T-ReCS and MICHELLE). Details on these observations are given in Table 2, and we used the fully reduced and flux calibrated spectra noted.

Name
Flux Density (mJy)   Table 4. Summary of the nuclear MIR emission derived from the AGN contribution based on spectral decomposition of Spitzer/IRS spectra and the SOFIA/FORCAST 31.5 µm fluxes. Column 1 corresponds to the galaxy name. Columns from 2 to 6 list the 4.5, 5.5, 18.0, 25.0 and 30.0 fluxes, respectively. The final column 7, corresponds to the SOFIA/FORCAST 31.5 µm fluxes, reported by Fuller et al. (2016). Note that † corresponds to nuclear fluxes calculated at 28 µm instead of 30 µm. See Section 4.1 for further details.

SED construction
To construct the entire nuclear IR SEDs sampling similar physical scales, we use NIR nuclear fluxes from our own GTC/CIRCE observations, HST archival data, or the highest angular resolution nuclear IR fluxes available in the literature. For those cases in which the angular resolution available is greater than 1 ′′ or there is evidence of a possible extra contribution from the host galaxy we used them as upper limits (see Table 2). When available, we used the subarcsecond nuclear spectra extracted as a point source, resampling them to 50 points, following the same methodology as in previous works using N-band nuclear spectra and clumpy torus models (e.g. AH11; Ramos Almeida et al. 2014;García-Bernete et al. 2015). In general, there is a good agreement between the flux calibration of the nuclear spectra and the N-band nuclear fluxes. However, for consistency, we systematically scaled the spectra to the N-band nuclear fluxes, unless there is any evidence to discard them due to the possible contribution of either emission lines or PAH features in the specific spectral window of the filters (e.g. NGC 7582). We estimated a ∼15% total uncertainty for the nuclear spectra by quadratically adding the errors in the flux calibration and point source extraction.
In addition, we estimated the AGN contribution at 5.5, 25 and 30 µm for all the galaxies based on spectral decomposition of Spitzer/IRS galaxies (see Table 4) 3 . To do so, we first scaled the AGN component to the N-band fluxes and then calculated homogeneous nuclear fluxes at 5.5, 25 and 30 µm using a 1 µm window in the scaled AGN component, using the same method as in GB16. We remark that when a specific rest-frame AGN template extends down to ∼4 µm, which occurs for roughly half of our sample (11/24 sources), we also derived the 4.5 µm nuclear fluxes (see Table 4). Finally, for those sources without Q-band (17-25 µm) photometry (e.g. NGC 4388), we calculated the 18 µm fluxes using the same methodology.
Five sources lack high angular resolution nuclear spectra (NGC 4395, NGC 6300, NGC 6814, NGC 7314 and ESO 005-G004). Nevertheless, we have high angular resolution photometry in the N-and Q-bands, and we then used the scaled AGN components derived from the IRS spectra to obtain N-band "pseudo-nuclear" spectra (e.g. Hernán-Caballero et al. 2015). For consistency with the other 19 nuclear IR SEDs, we restricted the scaled AGN component to have the same wavelength range (7.5-13 µm) as the ground-based spectra and resampled to 50 points. Note that we also use the "pseudo-nuclear" spectra for NGC 4138, NGC 4945, NGC 7172 and UGC 6728. In the case of NGC 4945 and NGC 7172, their nuclear spectra show a strong contribution from the host galaxy, while those of NGC 4138 and UGC 6728 are very noisy and practically identical in spectral shape to the AGN component.

SED observational properties
In Fig. 1 we present the nuclear IR SEDs (∼1-30 µm) of our sample of Seyfert galaxies. In these plots we compare the spectral shapes and the average nuclear IR SEDs for the different Seyfert types considered in this study (Sy1, Sy1.8/1.9 and Sy2 galaxies). The average nuclear IR Sy1, Sy1.8/1.9 and Sy2 templates were constructed using the nuclear IR SEDs described in Section 4.1, but excluding the lowest angular resolution data (i.e. upper limits). For consistency, we used the same wavelength grid for all the photometry (1.6, 2.2, 5.5, 8.8, 18.0, 25.0, 30.0 µm). To do so, we performed a quadratic interpolation of nearby measurements for each galaxy. In this process, we avoid using L-and M-bands due to the large number of upper limits at these wavelengths. Note that we computed the interpolated fluxes for the sole purpose of deriving the average Seyfert templates. In addition, we used N-band spectra, either the subarcsecond angular resolution or the "pseudo-nuclear" spectra (see Section 4.1).
We measured the NIR (1.6-8 µm), MIR (8-18 µm) and total IR (1.6-25 µm) spectral indices (fν α ν α ), the H/N and N/Q flux ratios, and the strength of the silicate feature (9.7 µm) for each galaxy in the sample. We also repeated these measurements in the derived Sy1, Sy1.8/1.9 and Sy2 templates, which are representative of each group of SEDs (see Table 5). We find steeper IR slopes for Sy2 than for Sy1, and the Sy1.8/1.9 and Sy1 slopes are very similar. Steeper IR slopes for type-2 AGN have been previously reported in the literature for Seyfert galaxies (e.g. Alonso-Herrero et al. 2003, RA11 and references therein) and more luminous AGNs (e.g. Mateos et al. 2016). In addition, we measured practically the same MIR slopes for the three groups (αMIR ∼-2) within the errors, in good agreement with the results reported by RA11. Following the same methodology as in RA09 and RA11, we also compare the spectral shapes of the different Seyfert types using the H/N and N/Q flux ratios. In agreement with the values reported by the latter authors, we found similar N/Q flux ratios (∼0.3-0.2). On the other hand, we found that Sy1 (0.11±0.05) and Sy1.8/1.9 (0.15±0.12) galaxies have slightly larger values of the H/N flux ratio than those of Sy2 (0.04±0.05), but the values are consistent within the errors.
Taking advantage of the spectroscopy data we compare the strength of the silicate feature (9.7 µm) for the different Seyfert types (see Table 5). The latter is computed as S Sil = ln (fcont/f9.7), where fcont and f9.7 are the flux densities of the continuum and the feature, which we measured at 9.7 µm. As can be seen from the top-left panel of Fig.  1, the majority of the Sy1 galaxies show weak or moderate emission (S Sil >0; the only exception is NGC 3227, which has a value of -0.2 and it could be related to the emission of PAHs), whereas Sy1.8/1.9 and Sy2 galaxies have relatively deep silicate features (S Sil =-0.3 and -1.0, respectively). This feature is normally observed in weak emission or absent in Sy1 and in shallow absorption in type 2 Seyfert galaxies when observed in subarcsecond resolution data (e.g. Alonso-Herrero et al. 2016, García-González et al. 2017 and references therein).

Accretion disk fitting
The NIR emission of AGN is mainly produced by the emission of very hot dust and the direct emission from the AGN (i.e. accretion disk) in the case of type 1s, although another important contribution can be stellar emission from the host galaxy. The contribution from the accretion disk declines with increasing wavelength. According to both theoretical models (e.g. Hubeny et al. 2001) and polarized light observations (e.g. Kishimoto et al. 2001) the NIR emission of the accretion disk can be explained by a power-law extension of the optical/UV spectrum to the NIR range. This power-law extrapolation is commonly used to fit the AGN direct emission in Seyfert 1 galaxies (e.g. Stalevski et al. 2012). However, the clumpy torus models of Nenkova et al. (2008a,b) assume a steep fall of the disk spectrum for wavelengths longer than 1 µm. We note that the CLUMPY models cannot reproduce the NIR bumps observed in the SEDs of some Sy1s (e.g. Mor, Netzer, & Elitzur 2009;RA11;AH11). For example, Mateos et al. (2016) successfully reproduced the IR SEDs of a sample of X-ray selected quasars using a nontruncated disk component and the CLUMPY torus models.
In order to quantify the contribution from the accretion disk to the nuclear NIR emission, we follow the same procedure as described in Hernán-Caballero et al. (2016) using optical, NIR and MIR photometry (see Tables 2, 4 and 6) to fit the accretion disk emission for all Sy1 galaxies 4 in our sample. This method used a semi-empirical model consisting of a single template for the accretion disk and two blackbodies for the dust emission.
In Fig. 2, we present the fitting results and in Table  6 we list the fractional contribution of this component to the nuclear NIR emission. Using only the fits with subarcsecond resolution data, we find that the average contribution of the accretion disk to the J-, H-and K-band emission are 46±28, 23±13 and 11±5% in ∼0.4 ′′ apertures, which are in good agreement with the values reported by Hernán-Caballero et al. (2016) for the rest-frame J-, Hand K-band (48±16, 27±14 and 17±1%) using a sample of luminous quasars. We note that the largest contribution from the accretion disk to the NIR emission is found for NGC 4151. This is in agreement with previous works on this galaxy (e.g. Swain et al. 2003;Kishimoto et al. 2007;Riffel, Storchi-Bergmann, & McGregor 2009).
Since we find a significant contribution of the accretion disk emission in the NIR range of Sy1, we subtracted this component in all Sy1 galaxies prior to fitting the nuclear IR SEDs with torus models.

SED modelling with the CLUMPY torus models
Using the CLUMPY models and BayesClumpy, we fit all the nuclear NIR-to-MIR SEDs in our sample (See Appendix A). A detailed description of the CLUMPY model parameters (see Table 8) can be found in Nenkova et al. (2008a,b). For approximately half of our sample (13/24 sources; see Figure 1. Observed nuclear IR SEDs for the Sy1, Sy1.8/1.9 and Sy2 galaxies in the BCS 40 sample. Note that different colours and symbols correspond to the galaxies labelled in each panel.  Table 5. Spectral shape information of the nuclear IR SEDs. The strength of the 9.7 µm silicate feature is computed as S Sil = ln (fcont/f 9.7 ), where fcont and f 9.7 are the flux densities of the continuum and the feature, which we measured at 9.7 µm.  Table 7 for further details). In addition, we used the IR extinction curve of Chiar & Tielens (2006) of the local ISM to account for any possible foreground extinction from the host galaxy. This curve covers the range ∼1-35 µm and accounts for the two silicate features at 9.7 and 18 µm. We used different priors for the foreground extinction from the host galaxy (A f or V ) for the various Seyfert types, taking into account the values avail-able in the literature (see Table 1). We used A f or V =[0,2] mag for Sy1 and [0,8] mag for Sy1.8/1.9 and Sy2. Finally, we used uniform priors for the rest of the parameters. When the observed data introduce sufficient information into the fit, the resulting posteriors will clearly differ from the input uniform distributions, either showing trends or being centered at certain values within the considered intervals.
We note that for this study we used the updated version (October 2014) of the Nenkova et al. (2008a,b) clumpy Figure 2. Accretion disk emission fits of Sy1s. Blue circles, red triangles, and black diamonds represent broadband photometry in the rest-frame optical, NIR, and MIR, respectively. We show the dust components (red dashed lines), accretion disk component (blue dot-dashed line) and best fits (black solid lines). Note that the dust component is modelled as a linear combination of two black-bodies at adjustable temperatures (green and orange dotted lines). The vertical dotted lines mark the rest-frame wavelength of Hα and Paα.      (2000); i) Bryant & Hunstead (1999); j) Wilson, Baldwin, & Ulvestad (1985). † Derived from NLR kinematics modelling.
torus models 5 . Older versions of these models used the optical depth along the slab normal for the synthetic clouds. However, in a recent comparison with spherical clouds (3D radiative transfer), the calculations showed that the effective optical depth through a cloud was two times higher than in the former approach (Heymann, Nikutta, and Elitzur, in 5 https://www.clumpy.org/  Table 8. Clumpy torus model parameters. i=0 • is face-on and i=90 • is edge-on. We note that the foreground extinction is unrelated to the torus. preparation). Although the absorption caused by clouds is not affected by this, the cloud emission does change since its source function is wavelength-dependent. As a consequence, a moderate change in the spectral shape has been reported on the CLUMPY webpage (less than 20% at any given wavelength).
In Appendix A, we present the results of the nuclear IR SED fitting process with the CLUMPY models (see Section 4.4), which are the marginal posterior distributions of the six parameters that define these models plus the foreground extinction and vertical shift. This shift scales with the AGN bolometric luminosity. We can also translate the posterior distributions of the parameters into a best-fitting model described by the combination of parameters that maximizes the posterior (maximum-a-posteriori; MAP) and a median model, computed with the median value of each posterior (see Appendix A). We found different average models of each subgroup from the median fitted nuclear IR SEDs. The Sy1 average model including the accretion disk emission component (black dotted line of Fig. 3) shows a flat NIR slope and the shape for the Sy2 average model (red solid line of Fig.  3) is very steep. The Sy1.8/1.9 average model lies between those of the Sy1 and Sy2 models.  Table 9. AGN and torus model properties derived from the fits and X-ray properties. The torus model properties were derived from the median values of the marginal posterior distributions. Columns 2, 3, 4, 5, 6 and 7 list the AGN bolometric luminosity (L model bol ), torus bolometric luminosity (L torus bol , obtained by integrating the corresponding model torus emission), torus gas mass (M T orus ), torus outer radius (r T orus ), escape probability (Pesc) and torus covering factor (C T ) derived from the torus model. Note that median values are listed with their corresponding ±1σ values around the median. Columns 8, 9 and 10 correspond to the hydrogen column density, intrinsic 2-10 keV and absoption-corrected 14-195 keV X-ray luminosities taken from Ricci et al. (2017). Column 11 and 12 list the BH masses with their references and the derived Eddington Ratio following the same methodology as in Ricci et al. (2017c). References for M BH /M ⊙ : a) This work; b) Koss et al. (2017); c) Vasudevan et al. (2010). See Appendix C for further information about the estimation of the BH masses.
In general, the CLUMPY models provide good fits (χ 2 /dof (degrees of freedom) <2.0) to the majority (19/24) of the nuclear IR SEDs (see Appendix A). While the MIR emission is well fitted for practically all the SEDs, we found that 5/24 galaxies (i.e. NGC 3783, NGC 4395, NGC 5506, NGC 7172 & NGC 7314) show a clear excess of emission in the NIR. This likely indicates an extra hot dust component is needed to reproduce their IR SEDs (see also Mor, Netzer, & Elitzur 2009). We note that the main goal of this work is to obtain a global statistical analysis of the clumpy torus model parameters of the various Seyfert galaxy types, rather than focussing on the individual fits (see Appendix A). As a sanity check, we repeated the analysis using only those galaxies with the best (χ 2 /dof <1.0; ∼63% of the sample) and good (χ 2 /dof <2.0; ∼79% of the sample) fits and we find the same results within 1σ.

COMPARISON OF THE TORUS PROPERTIES
In this section, we investigate the main differences between the clumpy torus model parameters for the BCS40 sample. Table 9 reports the main derived properties of the torus from the model parameters and the X-ray measurements for comparison.

Optical Classification
In this section we discuss the global posterior distributions and their mean values for Sy1, Sy1.8/1.9, and Sy2. To this end we apply the hierarchical Bayesian approach also used in Ichikawa et al. (2015). To be more precise, we assume that the global properties of the objects are extracted from common prior distributions and we infer the hyperparameters of these priors. Because of its flexibility, we decided to use beta distributions as these prior distributions. To take advantage of the already computed sampling of the posterior for each individual object we leverage the importance sampling trick developed by Brewer & Elliott (2014). Although one should ideally sample from the full hierarchical probabilistic model, we consider this approximate technique as sufficient for our purposes.
As can be seen from Figure 4, the majority of the global distributions are clearly different. To quantify these differences we use the Kullback-Leibler divergence (KLD; Kullback & Leibler 1951) as in RA11. This approach takes into account the overall shape of the posterior distribution and it always has a positive value. In the case of two identical distributions it is equal to zero and the larger the values the more different the distributions. RA11 suggested that for values larger than one (bold face in Table 10), two posterior distributions may be considered to be significantly different. Following this, we find that the differences in σ, N0 and τV (see Fig. 4) between Sy1 and Sy2 are significant according to the KLD. The same applies to σ and τV between Sy1.8/1.9 and Sy2 galaxies. We note that RA11 found essentially the same differences between Sy1 and Sy2 galaxies. All these results are in good agreement with previous works (e.g. RA11,   AH11 and Ichikawa et al. 2015). However, we find smaller values of the cloud optical depth for Sy2 (τV ∼56) than for Sy1 and Sy1.8/1.9 galaxies (τV ∼94-114) and smaller values of i for Sy1 (i=19±16 • ).

X-ray Classification
So far, we have compared the torus properties for Seyfert galaxies with different optical classifications. In this section we obtained the global posterior distributions of the sample divided into unabsorbed (NH <10 22 cm −2 ) and absorbed (NH >10 22 cm −2 ) Seyfert galaxies in X-rays.
In Fig. 5 we show these distributions. We find essentially the same trends as when we divide the sample into Sy1 and Sy2 using an optical classification, but with less significance (see Table 10). This is due to the fact that the majority of optically classified Sy1 and Sy2 correspond to the unabsorbed and absorbed subgroups, respectively. Half of the Sy1.8/1.9 galaxies are classified as unabsorbed and the other half as absorbed, while only one Sy1 galaxy is absorbed (i.e. NGC 4151).

Torus Size, Angular Width and Mass
In this section we discuss the main torus model properties: torus size, angular width and mass, which can be compared to those derived from high angular resolution observations from the Atacama Large Millimeter/submillimeter Array (ALMA; e.g. García-Burillo et al. 2016 andAlonso-Herrero et al. 2018). We can derive the physical radius of the clumpy torus (Ro) by using the radial extent of the torus (Y=Ro/R d ), the bolometric luminosity and the dust sublimation radius (R d ) definition (equation 1). Note that we use bolometric luminosities derived by using the 14-195 keV band and a fixed bolometric correction (L bol /L 14−195 keV =7.42). We obtained this factor from the commonly used bolometric correction of 20 (e.g. Vasudevan & Fabian 2009) at 2-10 keV and assuming a power-law slope of 1.8 (see Appendix B). Finally, we can also estimate the torus gas mass for each AGN, using the Galactic dust-to-gas ratio from Bohlin, Savage, & Drake (1978) and σ, N0, τV , R d and Y (equation 3), where Iq = 1, Y/(2 ln Y) and 1/3 for q =2, 1 and 0, respectively. In agreement with previous works using the CLUMPY torus models (e.g. RA09; RA11; AH11;Lira et al. 2013;Ichikawa et al. 2015;Fuller et al. 2016), we find relatively compact torus sizes for all the Seyfert galaxies in our sample (Ro <15 pc), with median values of 2.8±1.2, 1.9±1.2 and 3.5±3.9 pc for torus radius of Sy1, Sy1.8/1.9, and Sy2 galaxies. ×10 5 M⊙ for the torus gas masses of Sy1, Sy1.8/1.9, and Sy2 galaxies. Lira et al. (2013) found values of the torus masses ranging from 10 4 -10 6 M⊙ using a sample of 48 Sy2 galaxies from the extended 12µm Galaxy Sample (Rush, Malkan, & Spinoglio 1993). For 5/8 Sy1 galaxies, we can compare the torus gas masses with measurements corresponding to the inner 30 pc (radius) reported by Hicks et al. (2009). They derived masses of M H 2 gas =3-20×10 6 M⊙ for NCG 3227, NCG 3783, NCG 4051, NCG 4151 and NCG 6814. These masses were obtained from the H2 1-0S(1) emission line at 2.12 µm. We do not have any Sy2 galaxy in common with Hicks et al. (2009), but we com- pare with the two Sy2 galaxies in their sample (Circinus and NGC 1068). For Circinus using a smaller radius of 9 pc, they found a M H 2 gas =1.9 ×10 6 M⊙. For NGC 1068 the latter authors reported a mass of M H 2 gas =2.3 ×10 7 M⊙ in the inner 30 pc. As expected, we found that the gas masses inferred from the fit of the nuclear IR SEDs in a smaller radius (∼0.5-15 pc) are smaller than those measured in the inner ∼30 pc. Using the CO(6-5) line observed with ALMA, García-Burillo et al. (2016) reported a smaller gas mass (1.2×10 5 M⊙) for the inner ∼7-10 pc of NGC 1068. Finally, using ALMA/CO(2-1) data, Alonso-Herrero et al. (2018) found that the nuclear disk of the Sy2 NGC 5643 is a factor of ∼10 more massive and larger (∼26 pc of diameter) than that of NGC 1068. Therefore, we obtain comparable values of the torus mass with those derived from the highest angular resolution data.

The Covering Factor
Clumpy torus models imply that the differences between type-1 and type-2 AGN depend of whether there is a direct view of the broad line region (BLR; e.g. Nenkova et al. 2008a,b). Therefore, the observed classification is the result of the probability for an AGN-produced photon to escape through the torus along a viewing angle without being absorbed. As this probability is always non-zero, it is always possible to have a direct view of the BLR, regardless of the torus orientation. Therefore, the larger the covering factor (CT ) the larger the probability of classifying an AGN as type 2. In fact, the geometrical covering factor gives the type 2/total fraction (e.g. Mateos et al. 2017).
In order to compare the covering factors for the three subgroups we derived the combined probability distributions. To do so, we concatenated together the individual arrays of the CT values returned Bayesian modelling for all objects in subgroups and we computed the combined probability distributions since it is a nonlinear function of the torus model parameters (σ and N0; see left panel of Fig. 6). We note that we do not use the hierarchical Bayesian approach since if we use the generalized beta distribution as the prior we would introduce an extra prior in the CT derived quantities. Using our ultra-hard X-ray volume-limited sample of Seyfert galaxies, we confirm the results first reported by RA11 that the covering factors of Sy2 are larger than those of Sy1 galaxies. Indeed, using the optical classification, we find that Sy2 have larger median values of the covering factor combined probability distributions (CT =0.95± 0.04 0.18 ) than Sy1 (CT =0.66± 0.16 0.52 ) and Sy1.8/1.9 (CT =0.53± 0.21 0.37 ) which is the same as for Sy1s, within the errors. These results are in good agreement with previous works (e.g. RA11, AH11, Ichikawa et al. 2015 andMateos et al. 2016). We also repeat the global posterior distribution for the covering factor using the X-ray classification (see Section 5.1.2) and we find the same trend (see right panel of Fig. 6 and Table 10).

Dependence with AGN luminosity and Eddington ratio
To investigate the relation between the bolometric luminosity derived from the X-rays (2-10 and 14-195 keV; see Appendix B) and the covering factor (see top panels of Fig. 7), we divided our sample into several luminosity bins. In the first bin we included the three sources with log(L AGN bol )<43, while the rest of the sample was divided in two bins of equal logarithmic width (1 dex). We find that the σ parameter remains essentially constant within the errors, throughout our luminosity range (log(L AGN bol )∼41-45 erg s −1 ; see top panels of Fig. 7). The same applies for the covering factor (see top panels of Fig. 8). We find slightly larger values of CT ( 0.5) in the log(L AGN bol )∼44-45 erg s −1 luminosity range because the majority of the sources in that bin are Sy2s. Thus we do not find a statistically significant trend in the covering factor with AGN luminosity, which is in good agreement with recent studies (e.g. Mateos et al. 2016Mateos et al. , 2017Netzer et al. 2016;Stalevski et al. 2016;Lani, Netzer, & Lutz 2017;Ichikawa et al. 2018).
More recently it has been suggested that the Eddington ratio is the key parameter determining the covering factor, instead of the bolometric luminosity (e.g. Buchner & Bauer 2017;Ricci et al. 2017c). Ricci et al. (2017c) found that the covering factor rapidly decreases at higher Eddington ratios (see the orange solid line of Fig. 9). We derived Eddington ratios using the 2-10 keV bolometric X-ray luminosities and the black hole mass estimates from Ricci et al. (2017) and Koss et al. (2017), respectively. For the remaining sources we estimate the black hole masses (see Appendix C) as in Koss et al. (2017). The only exceptions are NGC 7213 and ESO 005-G004, for which we take their black hole masses from Vasudevan et al. (2010). The Eddington ratios of the sample are listed in Table 9.
Although we find higher values of the CT for lower Eddington ratios (see Fig. 9), we do not find a statistically significant dependence of the torus covering factor with the Eddington ratio. This result suggests, albeit for a small luminosity range and a limited number of galaxies, that the Eddington ratio would not be driving the geometrical covering factor.

Missing Obscured Seyferts?
In Fig. 10 we show the Sy2 fraction in our sample as estimated by using two covering factor bins (0.5-0.8 & 0.8-1.0) 6 . To estimate the uncertainties, we used the bootstrap error estimation generating 10 6 mock samples of Sy1, Sy1.8/1.9 and Sy2s by randomly selecting sources using replacements, with their corresponding covering factor distributions from the original samples. Note that the number of Sy1, Sy1.8/1.9 and Sy2s in each mock sample keep constant the observed number of Seyferts (i.e. nSy1+n Sy1.8/1.9 +nSy2). Finally, for each source, we calculate the obscured fraction in each bin by integrating its probability distribution (see e.g. Mateos et al. 2017).
Our data points should follow the 1:1 blue line shown in Fig. 10 if our sample did not miss any high covering factor source (covering factor values ∼1). However, the Sy2/total fraction is always below the 1:1 line. In general, we found that the most highly absorbed sources are the ones with higher torus covering factors (see also RA11, AH11 and Mateos et al. 2016). All this suggests that even an ultra hard X-ray (14-195 keV) Swift/BAT selection is missing a significant fraction of highly absorbed type 2 sources with very high covering factor tori. This is expected since at column densities NH >10 23.5 cm −2 , even high energy photons (14-195keV X-ray band) are absorbed. An example of these missing sources could be NGC 4418, which is a very highly obscured AGN at ∼30 Mpc. It has a compact IR bright core with the deepest known silicate absorption but it is not detected in the Swift/BAT hard X-ray band (e.g. Roche, Alonso-Herrero, & Gonzalez-Martin 2015 and references therein). The result presented here agrees with those reported in Ricci et al. (2015) and Koss et al. (2016) at energies >10 keV and Mateos et al. (2017) at energies >4.5 keV. The latter authors inferred the existence of a population of X-ray undetected objects with high torus covering factor, especially at high bolometric luminosities (>10 44 erg s −1 ).

CONCLUSIONS
We present for the first time a detailed modelling of the nuclear dust emission of an ultra-hard X-ray (14-195 keV) volume-limited (DL <40 Mpc) sample of 24 Seyfert galaxies. We selected our targets from the Swift/BAT nine month catalog. Our sample covers an AGN luminosity range log(L 2−10 keV int )∼40.5-43.4 erg s −1 . We fitted the nuclear IR SEDs obtained with high angular resolution data (∼1-30 µm) with the CLUMPY models using a Bayesian approach. From these fits, we derived torus model parameters for the individual 24 galaxies. In the case of Seyfert 1s we took special care to subtract the accretion disk contribution from the observed nuclear SEDs using the type 1 QSO accretion disk template from Hernán-Caballero et al. (2016). The main goal of this work was to obtain a global statistical analysis of the clumpy torus model parameters of type 1 and 2 Seyfert galaxies. We used both optical (broad vs narrow lines) and X-ray (unabsorbed vs absorbed) classifications for our analysis. Using these classifications, we compared the global posterior distribution of the torus model parameters, rather than focusing on the individual fits.
We verified our previous results that type 2 Seyferts have tori with larger widths and more clouds than type 1/1.8/1.9s. These findings are independent of whether we use an optical or X-ray classification. We found that the covering factor is likely the main parameter driving the classification of Seyfert galaxies. We derived compact torus sizes (radius <15 pc), and gas masses in the 10 4 -10 6 M⊙ range for both types.
We derived geometrical covering factors for the individual galaxies and globally for Sy1s and Sy2s. In clumpy torus models the geometrical covering factor is a function of the angular size and the number of clouds. Using these distributions, we confirmed that Seyfert 2 galaxies have larger values of the covering factor (CT =0.95± 0.04 0.18 ) than type 1s (CT =0.66± 0.16 0.52 ) using, for the first time, an ultra-hard Xray selected sample of Seyferts. We found that the torus covering factor remains constant within the errors in our luminosity range and no clear dependence with the Edding-ton ratio. Finally, we compared the derived covering factor with the observed type 2 fraction for our sample. From this comparison, we found tentative evidence that even an ultra hard X-ray selection is missing a significant fraction of highly absorbed type 2 sources with very high covering factor tori, as also concluded by Mateos et al. (2017) at lower X-ray energies using a more distant and luminous sample of AGN.
We note that detailed studies such as this, carried out not only using larger samples of galaxies but covering wider luminosity and redshift ranges are needed to improve the statistics of the results we report here. In the future, this methodology may be applied to AGN samples using high angular resolution and sensitive MIR data, observed with the combined spectral coverage of NIRSpec and MIRI aboard the James Webb Space Telescope (JWST).

ACKNOWLEDGMENTS
IGB acknowledges financial support from the Instituto de Astrofísica de Canarias through Fundación La Caixa. IGB also acknowledges Oxford University and Durham University for their hospitality during his stays in 2017 August when this project was started and 2018 May, respectively. IGB also acknowledges Cardiff University for their hospitality from 2018 May to August. CRA and IGB acknowledge financial support from the Spanish Ministry of Science and Innovation (MICINN) through project AYA2016-76682-C3-2-P. IGB, AAH and FJC also acknowledge financial support through grant PN AYA2015-64346-C2-1-P (MINECO/FEDER). Funded by the Agencia Estatal de Investigación, Unidad de Excelencia María de Maeztu. CRA also acknowledges the Ramón y Cajal Program of the Spanish Ministry of Economy and Competitiveness. MPS acknowledges support from STFC through grant ST/N000919/1. OGM thanks to the PAPIIT UNAM project IA103118. S.M. acknowledges financial support through grant AYA2016-76730-P (MINECO/FEDER). IM and JM acknowledge financial support from the research project AYA2016-76682-C3-1-P (AEI/FEDER, UE). L.F. and C.P. acknowledge support from the NSF-grant number 1616828.
Finally, we are extremely grateful to the GTC staff for their constant and enthusiastic support, and to the anonymous referee for useful comments. each parameter around the median. The galaxies with subarcsecond angular resolution N-band spectra are labelled as 'nuclear spectrum', while those labelled as 'AGN template' are sources that either do not have high angular resolution nuclear spectra or are noisy/include a strong contribution from the host galaxy (see Section 4.1). In these cases we used the N-band "pseudo-nuclear" spectra.
We note that while the MIR photometry and spectroscopy are well fitted in the majority of the cases, we found for 5/24 galaxies (i.e. NGC 3783, NGC 4395, NGC 5506, NGC 7172 & NGC 7314; see Figs. A1-A4) a NIR excess that the CLUMPY models cannot reproduce. This suggests that an extra component of very hot dust is needs to reproduce their IR SEDs (see also Mor, Netzer, & Elitzur 2009). Therefore, we repeat the global posterior distribution of each subgroup considered here for only the best (χ 2 /dof <1.0; ∼63% of the sample) and good (χ 2 /dof <2.0; ∼79% of the sample) fits and we find the same results within 1σ (see Fig.  A5 and Fig. A6).

APPENDIX B: AGN LUMINOSITY
We can compare the AGN bolometric luminosities derived from the torus model fits with those derived from X-ray measurements. With this aim we compiled X-ray luminosities (2-10 and 14-195 keV; see Table 9) from Ricci et al. (2017) and used the fixed bolometric correction factor of 20 (Vasudevan & Fabian 2009) and 7.42 for the 2-10 keV and 14-195 keV bands, respectively. The latter was obtained from the fixed bolometric factor at 2-10 keV by assuming a power-law of 1.8 as in Trakhtenbrot et al. (2017). We found that the relationship between bolometric luminosites derived from both X-ray bands and those derived from the clumpy torus models show an offset (see Fig. B1) which had not been found in previous works (e.g. AH11). As a sanity check, we compare the bolometric luminosities derived from the older and new version of CLUMPY models and we find that all the sources are in the 1:1 line. This new finding should be mainly related to the CLUMPY scaling factors. Therefore, in Section 5.2, we used the bolometric luminosities derived from the harder X-ray band (14-195 keV) to obtaining tori sizes and masses.

APPENDIX C: BLACK HOLE MASSES
In this appendix, we estimate the black hole masses for the remaining sources in our sample not included in the work of Ricci et al. (2017c). The only exceptions are NGC 7213 and ESO 005-G004, for which we take their black hole masses from Vasudevan et al. (2010). To compare with the results reported by Ricci et al. (2017c), we follow the same methodology as in Koss et al. (2017) to estimate the black hole masses. To do so, we use broad lines (e.g. Pa β , Hα, and H β ) and the stellar velocity dispersion (σ * ), from the 0.85 µm calcium triplet (hereafter, CaT) and CO (H-and K-bands) absorption features, for Sy1s and Sy1.8/1.9/2s, respectively. See Table C1, C2 and C3.
To calculate black hole masses from the broad Pa β line, we used the formula from La Franca et al. (2015) (equation B1), which is based on the FWHM and luminosity of Pa β . In the case of broad Hα and H β lines, we used the formulas from Greene & Ho (2005) (equation B2) and Trakhtenbrot & Netzer (2012) (equation B3), respectively. For Sy1.8/1.9 and Sy2 galaxies, we used the formula from Kormendy & Ho (2013) (equation B4) that is based on the stellar velocity dispersion. We note that the relation from La Franca et al. (2015) was calibrated using a virial factor f = 4.31. We prefer to use the stellar velocity dispersion measurements derived from the CaT band, rather than those from the CO bands, when possible. Although the extinction at 0.85 µm (CaT) should be larger (∼5 times) than at 2.3 µm (CO K−band ), the CaT feature trace an old stellar population that is more representative of the galaxy dynamical mass (e.g. Riffel et al. 2015). CO absorption features may have a stronger contribution from young stars than the CaT feature. Indeed, the effect of the young stellar population on the stellar velocity dispersion measured from the CaT feature is practically insignificant (see e.g. Riffel et al. 2015).  Figure A4. Nuclear IR SED of the Sy2 galaxies in the sample normalized at 11.2 µm. Solid red and dashed blue lines correspond to the MAP and median models respectively. Grey curves are the clumpy models sampled from the posterior and compatible with the data at 1σ level. A 4.88 · · · · · · NGC 5506 · · · 8.29 · · · Sy2 galaxies MCG-05-23-016 · · · 7.98 · · · NGC 4945 · · · 7.96 7.47 NGC 5128 · · · · · · 7.94 NGC 6300 7.01 · · · · · · NGC 7582 7.52 8.02 8.01  Sy1.8/1.9 Sy2 Figure A6. Comparison between the torus covering factor (C T ) combined probability distributions for each Seyfert galaxy subgroup using only good fits (χ 2 /dof <2.0). Blue dot-dashed, green dashed and red solid lines represent the parameter distributions of Sy1, Sy1.8/1.9, and Sy2, respectively. Figure B1. Bolometric luminosities derived from the torus models versus Swift/BAT bolometric X-ray luminosities derived from the 2-10 keV and 14-195 keV X-ray band. We plot the 1:1 line for comparison. Blue diamonds, green squares and red circles are Sy1, Sy1.8/1.9 and Sy2 galaxies, respectively.