Hypercubes of AGN Tori (HYPERCAT). I. Models and Image Morphology

Unification of active galactic nuclei (AGNs; e.g., Antonucci 1993; Urry & Padovani 1995) has enjoyed success for three decades. It explained most of the observed dichotomy between type 1 and type 2 AGNs with the simple assumption of an axially symmetric, anisotropic, dusty obscurer—a torus—and its orientation to the observer's line of sight (LOS). AGN infrared (IR) spectral energy distributions (SEDs) were initially modeled as emission from smooth-density dusty tori (e.g., Pier & Krolik 1992, 1993; Granato & Danese 1994; Efstathiou & Rowan-Robinson 1995). However, the dust must be confined in discrete clouds lest the torus be dynamically destroyed very quickly (Krolik & Begelman 1988). In the X-ray regime symmetric variability of the obscuring column density, N_H, suggests a circumnuclear medium composed of discrete gas and dust clouds; individual clouds are witnessed crossing the LOS while orbiting well inside the dusty parts of the torus (Markowitz et al. 2014, and references therein).

Radiative transfer models of clumpy circumnuclear dust are now numerous (e.g., Nenkova et al. 2002, 2008a, 2008b, hereafter N08a, N08b; Hönig et al. 2006; Schartmann et al. 2008; Stalevski et al. 2012; Siebenmorgen et al. 2015). They have all been extensively mined for their resulting SEDs. Studies of the spatially resolved emission morphology, on the other hand, have so far been underappreciated, presumably because resolved observations are only recently possible. Resolved imagery is the only remedy for model degeneracies that are inherent to IR radiative transfer, which is when different dust geometries produce similar (unresolved) SEDs (Vinković et al. 2003). Before the dawn of 30 m class telescopes, only interferometry has the potential to spatially resolve subparsec mid-IR (MIR) emission around the nucleus. The VLTI/MIDI instrument has resolved several nearby bright AGNs at MIR wavelengths (e.g., Jaffe et al. 2004; Poncelet et al. 2006; Tristram et al. 2007; Raban et al. 2009; Burtscher et al.2013; Tristram et al. 2014; Leftley et al. 2018).

MIDI combines the light of only two telescopes; thus, phase closure cannot be achieved, and images cannot be reconstructed without model assumptions. Most authors resorted to modeling directly with a combination of 2D Gaussian components in the plane of the sky the brightness distribution that VLTI has sampled. The semiaxes of the Gaussians were varied, as were their orientations, the effective blackbody temperatures generating each Gaussian, and the optical depth and dust properties along the LOS toward each component. A well-fitting model such as the one applied to VLTI/MIDI observations of Circinus (Tristram et al. 2014) required no fewer than 18 free parameters for three Gaussian components. These empirical models do not afford unambiguous conclusions about the underlying 3D dust distribution.

At longer wavelengths Atacama Large Millimeter/submillimeter Array (ALMA) has recently resolved the cold molecular gas associated with the outer regions of the torus in NGC 1068 and other nearby Seyferts (García-Burillo et al. 2016; Gallimore et al. 2016; Imanishi et al. 2016; Alonso-Herrero et al. 2018, 2019; Combes et al. 2019). In the particular case of NGC 1068 ALMA observations found a torus no larger than ∼11 pc in diameter with tracers of dense molecular gas (e.g., HCO⁺ (4–3)) and 2–3 times larger with lower-density tracers such as CO (2–1) and CO (3–2) (García-Burillo et al. 2019), as well as evidence for nuclear molecular outflow (Gallimore et al. 2016). In addition, ALMA polarimetric observations at 860 μm with an angular resolution of 007 have measured the signature of magnetically aligned dust grains in the torus of NGC 1068 (Lopez-Rodriguez et al. 2020). These results challenge the static nature of the torus in favor of an obscuring structure generated by a hydromagnetic disk outflow, which has been suggested from a theoretical point of view (e.g., Emmering et al. 1992; Elitzur & Shlosman 2006; Wada 2012; Elitzur & Netzer 2016; Dorodnitsyn & Kallman 2017). Recently, more detailed VLTI measurements of two dozen nearby AGNs have complicated the simple unification picture further. In several studied objects most of the MIR emission appears to emanate from polar regions high above the equatorial plane, i.e., not from where the bulk of the dust is thought to reside (Hönig et al. 2012, 2013; Burtscher et al. 2013; López-Gonzaga et al. 2016; Leftley et al. 2018). Tristram et al. (2014) suggested, for the case of Circinus, that the elongation is due to seeing the inner funnel wall of a torus slightly tilted toward us. Others propose that the polar emission is caused by a dusty outflow and have attempted to explain the observed elongations, e.g., by modifying the distribution of dusty clouds in the models (Hönig & Kishimoto 2017), or through interplays between a tilted accretion disk that anisotropically illuminates a hollow, biconical, dusty wind (Stalevski et al. 2017, 2019).

In the near future extremely large single-aperture telescopes (henceforth single-dish for simplicity) such as the Giant Magellan Telescope (GMT), Thirty Meter Telescope (TMT), and Extremely Large Telescope (ELT) will also be able to resolve the nuclear AGN thermal emission in several nearby sources, as we will demonstrate in the companion paper (Nikutta et al. 2021, hereafter Paper II). ALMA continues to improve its capabilities, yielding an ever more detailed view of AGN tori at millimeter wavelengths (Imanishi et al. 2018; García-Burillo et al. 2019; Impellizzeri et al. 2019). Clearly, there is a strong and growing demand for systematic studies and modeling of the resolved AGN emission, i.e., imagery at all relevant wavelengths. Further progress can arrive through brightness maps obtained from radiative transfer in physically motivated dust distributions, but not through models of the directly observed brightness distribution in the sky. The latter have no connection to the astrophysical entities and teach us little about the circumnuclear matter distribution.

In this paper we introduce in Section 2 the n-dimensional hypercubes of brightness and dust maps made with Clumpy, a radiative transfer code for clumpy media (N08a, N08b). In Section 3, using techniques based on image moments, we investigate the morphologies of light and dust maps as a function of model parameters and investigate how such models can generate significant polar elongation in the MIR. We summarize the results in Section 4. A number of appendix sections and a User Manual ¹⁰ distributed with Hypercat shall serve as technical reference and guide users in their work with the Hypercat software and model hypercubes.

In the companion Paper II we will then apply the models to simulate various synthetic observations of NGC 1068, including with 30 m class telescopes.

The need to model integrated-light SEDs from AGNs via radiative transfer was recognized long ago. However, studies of resolved images remain in their infancy. The reasons are twofold. First, until recently no spatially resolved data were available. This has been remedied by VLTI and ALMA, and new facilities such as GMT, TMT, and ELT will add new data in the near future. Second, the computation, storage, and processing requirements for image sets with parameter coverage equivalent to the model sets used in SED studies are prohibitively expensive.

Hypercat is designed to break this second barrier. It is a user-friendly software that abstracts away the underlying complexity of handling very large n-dimensional hypercubes of data—in this case images of emission and projected dust distributions that are complex functions of several independent model parameters. Hypercat enables one to generate an image of the torus emission (or a dust distribution) for any combination of parameters via multilinear interpolation in fractions of a second on off-the-shelf standard computers. It can also simulate observations of these images with single-dish telescopes, performing convolution with any point-spread function (PSF; provided by the user or computed by Hypercat), image transformations (e.g., rotation, scaling in flux and size, resampling), noise addition to specified signal-to-noise ratio, and image deconvolution. Interferometric observations can also be simulated for a set of (u,v) points, generating synthetic visibility measurements, which will be presented in a follow-up manuscript. Furthermore, spatially resolved maps of spectral properties can be simulated, akin to integral field units (IFUs). A module with various image morphology estimators completes the set of high-level tools built into Hypercat. Descriptions of Hypercat software interfaces, and of the (minimal) GUI program, are given in the User Manual.

The software is open-source, written in Python, and designed such that users can easily plug in their hypercubes of model images and perform any analyses facilitated by Hypercat. We document in the User Manual the exact requirements and provide a tool to create such hypercubes from a collection of image files, and we can also assist colleagues upon request.

With this first release of Hypercat we also release our hypercube of model images, generated with Clumpy torus models. They constitute a comprehensive set of thermal emission and dust distribution images, which will allow us to model and analyze data to be obtained by current interferometric facilities and future 30 m class telescopes. We describe the Clumpy/Hypercat models in the remainder of this section.

2.1. Distribution of Dust Clouds around the AGN

Clumpy models the distribution of dust clouds around an AGN with an axisymmetric function of cloud number density (per unit length) that is separable in its radial and polar arguments (r, β),

$$\begin{eqnarray}&&{N}_{C}(r,\beta )=C\,\displaystyle \frac{{{ \mathcal N }}_{T}(\beta )}{{R}_{d}}{\left(\displaystyle \frac{r}{{R}_{d}}\right)}^{-q},\end{eqnarray} \tag{ 1 }$$

where r is the radial distance from the AGN, β is the angle of a radial ray away from the equatorial plane, and ${{ \mathcal N }}_{T}(\beta )$ is the mean number of clouds encountered along such a ray. The free parameter q sets how sharply N_C falls off with distance.

The clouds, each with optical thickness τ_V in the V band, exist between the dust sublimation radius

$$\begin{eqnarray}&&{R}_{d}=0.4{\left(\displaystyle \frac{{\text{}}{L}_{\mathrm{bol}}}{{10}^{45}\,\mathrm{erg}\,{{\rm{s}}}^{-1}}\right)}^{1/2}{\left(\displaystyle \frac{{T}_{\mathrm{sub}}}{1500\,{\rm{K}}}\right)}^{-2.6}\,(\mathrm{pc}),\end{eqnarray} \tag{ 2 }$$

given in N08b and Barvainis (1987), and an outer radius R_out, with Y = R_out/R_d a free parameter. An outer radius is not too meaningful a parameter, as it will arbitrarily cut off the distant cloud distribution, but it is a necessary parameter for the sake of numerical simulations.

R_d is set by the AGN bolometric luminosity, L_bol, for a given dust grain sublimation temperature, T_sub (1500 K for a standard ISM silicate/graphite mix). The point source at the center illuminates the dust isotropically with luminosity L_bol. C in Equation (1) is a constant that ensures the correct normalization of ${{ \mathcal N }}_{T}(\beta )=\int {N}_{C}(r,\beta ){\text{}}{dr}$. The default polar-angle density behavior in Clumpy is a soft-edged Gaussian

$$\begin{eqnarray}&&{{ \mathcal N }}_{T}(\beta )={{ \mathcal N }}_{0}\exp (-{\beta }^{2}/{\sigma }^{2}),\end{eqnarray} \tag{ 3 }$$

with ${{ \mathcal N }}_{0}$ the mean number of clouds per radial ray in the equatorial plane and σ the width parameter of the Gaussian, measured in degrees from the equatorial plane. Other prescriptions for ${{ \mathcal N }}_{T}(\beta )$ are possible. For instance, Fritz et al. (2006), Feltre et al. (2012), and Stalevski et al. (2012) use a cosine in the exponential. Whatever the exact nature of the polar edge of the torus, observations require it to be soft and not sharp (N08a).

All free parameters of the model geometry and the external viewing angle of the observer form the vector of model parameters ${\boldsymbol{\theta }}=(\sigma ,i,Y,{{ \mathcal N }}_{0},q,{\tau }_{{\text{}}V})$. Their sampling in Hypercat is listed in Table 1.

Table 1. Parameter Sampling

Parameter	Sampled Values	Units	Dep. a	${N}_{{\theta }_{k}}$ b
σ	15, 30, 45, 60, 75	deg	b, d	5
i	0, 10, ..., 90	deg	b, d	10
Y	5, 6, ..., 20		b, d	16
${{ \mathcal N }}_{0}$	1, 2, ..., 12		b c	12
q	0, 0.5, 1, 1.5, 2		b, d	5
τV	10, 20, 40, 60, 80, 120, 160		b	7
Wave d	1.2, ..., 945	μm	b	25
Total number of combinations without wave:				336 k
Total number of combinations with wave:				8.4 M

Notes.

^a Maps that depend on this parameter: b = brightness map, d = dust map. ^b Number of sampled values for each parameter. ^c ${{ \mathcal N }}_{0}$ is only a multiplicative factor (see text). ^d See Table 2 for wavelength sampling. Other parameter values, for instance, larger Y, can be added in the future based on community needs.

Download table as:  ASCIITypeset image

Note that Equation (1) prescribes a continuous cloud number density at any location (x, y, z). Unlike any other published model of clumpy tori, Clumpy does not place macroscopic clouds at discrete locations within the computational domain, i.e., it does not introduce clumpiness via Monte Carlo simulations. Instead, clumpiness is dealt with analytically by computing the local statistical properties of the stochastic medium. This is computationally very fast but produces images of smooth appearance, despite the clumpiness of the dust. The computed images are simply the limit (mean) of infinitely many random realizations of a discrete model.

We account for two classes of clouds: those irradiated directly by the central source (their source function is S_λ ^dir), and those that do not have a clear LOS to the AGN and are heated only indirectly by an isotropic radiation field provided by other clouds around it (their source function is S_λ ^ind). The directly illuminated clouds have a hot, AGN-facing side and a cooler, averted side. Their source function is therefore highly anisotropic. An observer will see a mix of both the hot and cool sides, depending on the position angle with respect to a cloud, like the phases of the Moon. We precompute the cloud emission source function S_C,λ for a finely sampled range of distances r from the AGN (to account for heating by the incident flux at that distance) and position angles α (to account for the anisotropic illumination). At all other locations the source function is smoothly interpolated. The overall cloud source function (see Equation (10) in N08a) at any location is a mix of the two cloud types

$$\begin{eqnarray}&&{S}_{C,\lambda }(r,\alpha ,\beta )=p(r,\beta ){S}_{\lambda }^{\mathrm{dir}}(r,\alpha )+[1-p(r,\beta )]\,{S}_{\lambda }^{\mathrm{ind}}(r),\end{eqnarray} \tag{ 4 }$$

with the mixing proportion being set by the probability p(r, β) that a cloud at location (r, β) has a clear view of the AGN. This probability is $\exp \{-{ \mathcal N }(r,\beta )\}$, with ${ \mathcal N }(r,\beta )={\int }_{0}^{r}{N}_{C}(r,\beta ){\text{}}{dr}$ the mean number of clouds between the cloud and the AGN. The observed brightness map, integrated flux, and SED are then computed a posteriori and very quickly via ray tracing through a distribution of clouds. We describe this in the next subsection.

2.2. Monochromatic Emission Maps

Figure 1 shows the geometry of ray-tracing through a Clumpy torus model. The plane of the sky is (x, y), with y pointing up and x perpendicular to the page (pointing away from the reader). The AGN-to-observer LOS is z. The torus axis is in the (y, z)-plane, inclined from the z-axis by an angle i. At every coordinate (x, y) in the plane of the sky Clumpy computes the monochromatic intensity I_λ (x, y) along the path s (parallel to the z-axis) as the integral over a product of three functions

$$\begin{eqnarray}&&{I}_{\lambda }(s;x,y)=\int {\text{}}{ds}\,{S}_{C,\lambda }{N}_{C}(s^{\prime} ){P}_{\mathrm{esc}}(s^{\prime} ,s).\end{eqnarray} \tag{ 5 }$$

${S}_{C,\lambda }(s^{\prime} )$ is the photon-generating cloud source function at a given point $s^{\prime}$ along the path. ${N}_{C}(s^{\prime} )$ is the local cloud number density from Equation (1) and describes how many clouds are generating the photons. Finally,

$$\begin{eqnarray}&&{P}_{\mathrm{esc}}(s^{\prime} ,s)=\exp \left\{-{ \mathcal N }(s^{\prime} ,s)(1-\exp (-{\tau }_{\lambda }))\right\}\end{eqnarray} \tag{ 6 }$$

is the probability that a photon generated at $s^{\prime}$ propagates all the way along the rest of the path to the observer without being absorbed again, assuming that the number of clouds is Poisson distributed around the mean (Natta & Panagia 1984). ${ \mathcal N }(s^{\prime} ,s)={\int }_{s^{\prime} }^{s}{N}_{C}{\text{}}{ds}$ is the mean number of clouds between point $s^{\prime}$ and the observer, and τ_λ is the optical depth of a single cloud at wavelength λ. The convergence of Equation (5) is ensured by Clumpy through Romberg integration.

Zoom In Zoom Out Reset image size

Ray-tracing for all sky positions (x, y) in the field of view (FOV) generates a 2D image I_λ (x, y) at wavelength λ. Figure 1 shows two such maps, and their composite, at 9.5 μm (blue) and 860 μm (gold). Correct normalization of these brightness maps, as described in Appendix A of N08b, requires convergence of the area integral ∫dx dy I_λ (x, y); Clumpy handles this via the "Cuhre" integration with globally adaptive subdivision, provided by the Cuba library (Hahn 2005).

2.3. Wavelength Axis and SEDs

2.4. Simplifications and Caveats

2.5. Image Size and Sampling

2.6. Projected Cloud Number Maps

While ray-tracing, Clumpy also computes the integrated number of clouds along the LOS. For every path s parallel to the z-axis in the FOV (i.e., every pixel) this number is ${ \mathcal N }(s)=\int {N}_{C}{\text{}}{ds}$, with N_C from Equation (1). All pixels together make up a map of the projected dust cloud number distribution, C_d (x, y). Because of the same argument about azimuthal symmetry made in Section 2.5, C_d is always left-right symmetric in the image plane, regardless of inclination i.

Unlike the emission maps I_λ (x, y), the cloud maps C_d do not depend on λ or τ_V . Furthermore, from Equations (1) and (3) we see that the cloud distributions at fixed {σ, i, Y, q} are identical in structure for all ${{ \mathcal N }}_{0}$ up to a factor: ${{ \mathcal N }}_{0}$ itself. We can thus store just one cloud map C_d (${{ \mathcal N }}_{0}$ = 1) and later multiply it by the desired ${{ \mathcal N }}_{0}$ when using the cloud maps. This greatly reduces the storage requirements for all cloud maps.

The integral over all pixels in the cloud map yields the total number of clouds in the torus, n_tot = ∫C_d (x, y)dx dy. This number is preserved no matter what viewing angle the torus is inclined at. Note that n_tot changes with q. From Section 2.4 of N08a, n_tot = (N_C /A_C )dV, where the cloud cross section ${A}_{C}\simeq {R}_{C}^{2}$ (and R_C is the cloud size). From Equation (1), N_C ∝ r^−q. With fixed R_d , R_C , and ${{ \mathcal N }}_{0}$ we then obtain n_tot ∝ r^−q dr, i.e., when q increases, n_tot decreases. Figure 2 shows in the first five panels the projected cloud numbers per pixel (per LOS) for a torus with five values of q between 0 and 2, with other parameters fixed at σ = 43°, Y = 20, ${{ \mathcal N }}_{0}$ = 1, and at a viewing angle i = 90°. The radial cloud distribution is visibly concentrated with growing q. At the central pixel, i.e., at zero x and y offsets from the center, ${ \mathcal N }$ = 2 always, because ${{ \mathcal N }}_{0}$ = 1, i.e., there is on average one cloud along a radial ray in the equatorial plane (and one on the far side of the torus). The right panel shows n_tot as a function of q normalized to the maximum value of n_tot, which occurs at q = 0.

Zoom In Zoom Out Reset image size

These physically motivated maps of dust distributions, when compared to the morphology of the thermal emission emanating from the same torus, can help to understand the conditions that cause the polar elongation of the MIR emission observed in some AGN (or equatorial elongations in other sources, for that matter). The behavior of morphology as a function of position angle (PA), wavelength, and all model parameters can be studied. In Section 3 we analyze the dust and emission morphologies in detail.

2.7. Hypercube Model Files

The emergence of instruments powerful enough to resolve the dusty torus in nearby AGNs (e.g., ALMA, 30 m class telescopes) allows us not only to separate the flux contributions of the nucleus and host galaxy but also to quantify the morphology of the observed nuclear brightness distribution. Resolved VLTI observations (e.g., Hönig et al. 2012, 2013; Tristram et al. 2014; López-Gonzaga et al. 2016) show that in several nearby AGNs the MIR emission is elongated in polar directions, with the polar axis usually defined by other means, e.g., the orientation of ionization cones, a nuclear jet, or the perpendicular orientation of a maser disk. In these sources most of the MIR emission appears to be emanating from regions high above the equatorial plane of the "torus," i.e., not from where most of the dust is thought to reside according to AGN unification.

Whatever the mechanism for producing an elongated morphology, it must reconcile the 2D brightness map with the 3D dust distribution consistently. NIR and MIR imaging with extremely large single-dish telescopes will resolve the morphology unambiguously, and free of model assumptions. Phenomenological modeling of the brightness distribution can no longer serve as a surrogate.

In this section we investigate the morphologies of Clumpy brightness maps in several wavelength regimes and of their underlying dust distribution maps. We introduce and measure various morphological quantities, e.g., centroid location, half-light radii (1D), radii of gyration (2D), elongation (aspect ratio) of the morphology, position angle, asymmetry (skewness), and compactness (peakedness). We first investigate the morphological changes to brightness maps as a function of model parameters, with a focus on identifying the part of parameter space that allows for polar elongation of brightness distributions similar to those reported in the literature. We then compare these brightness maps to their corresponding dust maps, inferring as a function of model parameters the localized correspondence (or lack of it) of dust distributions and the emission morphologies generated by them. In Paper II we will then simulate observations with extremely large telescopes and compare to observational results obtained for NGC 1068.

3.1. Image Moments and Moment Invariants

We quantify the morphologies in terms of their image moments and quantities computed from moments, and we compare them to some traditional morphological quantifiers. The so-called central moments are

$$\begin{eqnarray}&&{\mu }_{\mathrm{pq}}=\sum _{x}\sum _{y}I(x,y){\left(x-\bar{x}\right)}^{p}\,{\left(y-\bar{y}\right)}^{q},\quad p,q\in {\mathbb{N}},\end{eqnarray} \tag{ 7 }$$

where $\bar{x}$ and $\bar{y}$ are the image centroid coordinates in pixels. By construction central moments are invariant under translation operations. Combinations of moments can be defined that are invariant under translation, scaling, rotation, or any combination of these operations. We refer to Appendix C for a technical discussion. Functions of image moments can be used to compute various quantities that characterize morphology. Table 4 lists the ones that we make use of in this paper.

Table 4. Useful Morphological Measurements with Moments

Measure	Method
Image mass (sum), or "flux"	M00 or μ00 (identical)
Image centroid location	$\bar{x}={M}_{10}/{M}_{00}$, $\bar{y}={M}_{01}/{M}_{00}$
Radius of gyration about axis	${R}_{g}^{x}=\sqrt{{\mu }_{20}/{\mu }_{00}}$, ${R}_{g}^{y}=\sqrt{{\mu }_{02}/{\mu }_{00}}$
Elongation in the y-direction	$e={R}_{g}^{y}/\,{R}_{g}^{x}$
Image covariance matrix	$\mathrm{Cov}[{\rm{I}}({\text{}}{\rm{x}},{\rm{y}})]=\left[\begin{array}{cc}{\mu }_{20}^{{\prime} } & {\mu }_{11}^{{\prime} }\\ {\mu }_{11}^{{\prime} } & {\mu }_{02}^{{\prime} }\end{array}\right]$
PA of main component	$\displaystyle \frac{1}{2}\arctan \left\{2\,{\mu }_{11}^{{\prime} }/\left({\mu }_{20}^{{\prime} }-{\mu }_{02}^{{\prime} }\right)\right\}\ (\mathrm{radian})$
Skewness	${S}_{x}={\mu }_{03}/\sqrt{{\mu }_{02}^{3}}$, ${S}_{y}={\mu }_{30}/\sqrt{{\mu }_{20}^{3}}$

Note. M_pq raw moment (Equation (C1)), μ_pq central moment (Equation (C3)), ${\mu }_{\mathrm{pq}}^{{\prime} }={\mu }_{\mathrm{pq}}/{\mu }_{00}$ second-order central moment (see Appendix C.2), I(x, y) image.

Download table as:  ASCIITypeset image

3.2. Size, Elongation, and Image Compactness

We wish to quantify the spatial extension of emission morphologies as a function of model parameters. Various measures have been proposed in the literature, e.g., the half-light radius, the FWHM of a fitted 2D Gaussian, or the Petrosian radius (Petrosian 1976; Graham et al. 2005). We show that a quantity based on image moments, the radius of gyration, is ideally suited to measure morphology size. We also introduce the Gini index as a simple way to quantify image compactness. For the remainder of the paper we can assume that there is only one extended source present in the image—the model torus.

3.2.1. Half-light Radius

A commonly used size estimate of an emission region is the half-light radius R_HL of a circular aperture that contains half of the total flux in the image (see Appendix C.5). The aperture is usually anchored on the origin of the coordinate frame. However, the half-light radius is sensitive to where the aperture is centered in the general case that the brightness distribution is not uniform. Figure 3 shows a 10.0 μm image of a Clumpy torus model at several viewing angles i, as indicated. All other model parameters are fixed.

Zoom In Zoom Out Reset image size

Three circular apertures are shown, centered at different loci, and each contains half the light of the image. Because the brightness distribution is not uniform across the image, these half-light radii differ in size.

In pole-on view (i = 0°) all half-light radii are identical, since all apertures are centered on the origin. Note that in pole-on view the "brightest pixel" is ill-defined; it is in fact a ring around the image origin. We thus set the aperture center to the center pixel. In edge-on views the three circles are only identical if the brightest pixel happens to be at the image center. This is the case if the central LOS is optically thin or mildly obscured. If the absorption along the central LOS increases, through either a higher ${{ \mathcal N }}_{0}$ or a higher τ_V , or through a change in wavelength, the emission peak occurs spatially much higher above the equatorial plane, and thus the half-light aperture constructed around the peak pixel is very different from the centroid- and origin-based ones.

Unresolved measurements of the total flux are of course not affected by the choice of the aperture center, since all flux is within the PSF. Interferometers without phase closure, such as VLTI with two telescopes, can resolve the nuclear emission in some nearby AGNs, but the emission cannot be localized within the FOV since absolute astrometry cannot be achieved. However, interferometers with phase closure (e.g., VLTI with MATISSE or GRAVITY, ALMA) and, of course, very large single-aperture telescopes (e.g., GMT, TMT, ELT) do attain absolute astrometry; thus, the locus at which the half-light radius aperture is anchored can be important. We note though that for extragalactic observations of (marginally) resolved nuclear cores, precise image registration is critical (Prieto et al. 2014).

In the edge-on scenario shown in Figure 3, the offset between the image origin and the brightest pixel is about 12.3 R_d . The physical size depends on luminosity; for NGC 1068 it is quite uncertain. Our own fitting (in Paper II) of the K-band image by Gravity Collaboration et al. (2020, hereafter G20) yields L_bol =1.56 × 10⁴⁴ erg s⁻¹ or L_bol = 2.62 × 10⁴⁴ erg s⁻¹, depending on the model; García-Burillo et al. (2014) find L_bol = 4.2 × 10⁴⁴ erg s⁻¹ from SED fitting; Lopez-Rodriguez et al. (2018) report L_bol = 5.02 × 10⁴⁴ erg s⁻¹ when including FIR data (see also Paper II); Alonso-Herrero et al. (2011) adopted L_bol = 1 × 10⁴⁵ erg s⁻¹; Raban et al. (2009) give L_bol = 1.5 × 10⁴⁵ erg s⁻¹; and Burtscher et al. (2013) list L_bol = 1.58 × 10⁴⁵ erg s⁻¹.

This range of luminosities translates to sizes between 1.9 and 6.2 pc for the 12.3 R_d offset and, with a distance of 14.4 Mpc to NGC 1068, between 28 and 89 mas angular extent on the sky. The larger values are well above the diffraction limit of 30–40 m class telescopes. Figure 4 shows the angular size of 12.3 R_d on the sky, as a function of AGN bolometric luminosity (which sets the distance of the dust sublimation radius) and of source distance, together with the diffraction limit of telescopes with 6.5–100 m diameters D (θ = 1.22 λ/D) at 10 μm. The bolometric luminosities for NGC 1068 listed above are marked in the figure. The effect described above is potentially resolvable with 30 and 40 m class telescopes, if good astrometry can be achieved. Smaller telescopes would struggle, unless a source were much brighter than NGC 1068, and/or closer than it (but such sources do not exist). The same is true if the bolometric luminosity of the source is on the lower end of the range. Other nearby AGNs at similar distances have larger tori, which makes the situation easier (e.g., Alonso-Herrero et al. 2018, 2019; Combes et al. 2019).

Zoom In Zoom Out Reset image size

Another drawback of half-light radii is that they cannot measure morphology sizes independently for orthogonal directions on the sky. In the following subsection we remedy this by introducing radii of gyration.

3.2.2. Radii of Gyration

A very robust way to measure image extension is via radii of gyration about the x- and y-axes

$$\begin{eqnarray}&&{R}_{g}^{x}=\sqrt{{\mu }_{20}/{\mu }_{00}},\qquad {R}_{g}^{y}=\sqrt{{\mu }_{02}/{\mu }_{00}}.\end{eqnarray} \tag{ 8 }$$

The gyration radius about an axis is the distance to a line parallel to the axis, at which all the "mass" (i.e., all the pixel values) could be concentrated without changing the second moment about that axis. Like other quantities computed from image moments, the concept of gyration radius is borrowed from mechanics. Figure 3 shows the gyration radii ${R}_{g}^{x}$, ${R}_{g}^{y}$ as the semimajor axes of a red ellipse, plotted over the brightness distributions produced by a Clumpy model. The ellipse is centered on the image centroid $(\bar{x},\bar{y})$.

When estimating morphology sizes, the radii of gyration have several benefits over other measures, such as the half-light radius. One is that they give realistic estimates of both the x and y extensions independently. This allows us then to compute the elongation or aspect ratio of the brightness distribution, which we will explore in the next subsection. Further, in all but the i = 90° case, where the image may show two bright centers of emission offset in the ±y-direction, the ellipse made of gyration radii and anchored on the image centroid follows the bulk of the emission more faithfully than circular apertures pinned to the telescope pointing (usually the suspected locus of the AGN).

Figure 5 shows how the radii of gyration, ${R}_{g}^{x}$ and ${R}_{g}^{y}$, and the half-light radius, ${R}_{\mathrm{HL}}^{\mathrm{origin}}$, depend on (a subset of) the torus model parameters. The analysis was performed on the "infinite" resolution model images (at 10 μm) to illustrate direct physical effects of parameter changes on image morphology. Plotted as a function of viewing angle, the half-light radius varies most extremely among the three shown size measurements. In most cases (${{ \mathcal N }}_{0}$ = 1 being the exception) ${R}_{\mathrm{HL}}^{\mathrm{origin}}$ is the lowest at pole-on viewing angles and the largest at edge-on orientations. That is, the half-light radius likely underestimates/overestimates the characteristic sizes at pole-on/edge-on viewings. The radii of gyration respond more gently to changes in model parameters. Naturally, ${R}_{g}^{y}$ varies more strongly than ${R}_{g}^{x}$ with viewing angle i and torus angular size σ, as these two quantities have the largest influence on the size of the image morphology in the y-direction. ${{ \mathcal N }}_{0}$ (and τ_V , which is not varied in this figure) provides the overall magnitude of the measured size.

Zoom In Zoom Out Reset image size

In the case of very high σ and ${{ \mathcal N }}_{0}$, i.e., an optically and geometrically thick dust "cocoon" scenario (lower right panel in Figure 5), both radii of gyration show very little variation with viewing angle, while the half-light radius still changes strongly between pole-on and edge-on viewing.

3.2.3. Aspect Ratio/Elongation

Because the radii of gyration about the y- and x-axes measure morphological size independently in both directions, their aspect ratio determines the elongation of the morphology

$$\begin{eqnarray}&&e={R}_{g}^{y}/\,{R}_{g}^{x}.\end{eqnarray} \tag{ 9 }$$

In our model images the position angle is always PA = 0°, i.e., ${R}_{g}^{y}$ corresponds to the polar direction and ${R}_{g}^{x}$ to the equatorial direction.

Using this technique, we now want to measure the elongation of torus images produced by Clumpy while varying the model parameters. Figure 6 shows the elongation measured from images for a number of parameter variations, at five wavelengths between NIR and submillimeter, and for two different radial cloud distribution power laws (flat, with q = 0.08 on the left, and steep, with q = 2 on the right). The measurements are considering all pixels in an image unless stated otherwise.

Zoom In Zoom Out Reset image size

Several characteristics are immediately apparent. All curves begin at precisely e = 1.0 for i = 0°. This is so because for pole-on LOSs the torus appears perfectly circular. Conversely, the largest and smallest aspect ratios are realized at or very close to edge-on views (i = 90°).

Another observation is that it is impossible to achieve any level of polar elongation at small values of σ or at the smallest values of ${{ \mathcal N }}_{0}$. A torus whose distribution of dust clouds is disk-like and concentrated in the equatorial plane will always appear elongated in the equatorial direction, i.e., perpendicular to the torus axis. At small q values no net polar elongation can be achieved at FIR and submillimeter wavelengths with any combination of parameters. This is because the total optical depth through the torus is small at these wavelengths, so the emission morphology mostly follows the dust distribution, i.e., it resembles the torus in appearance (see Section 3.4 for more details). For steep distributions (q = 2) no polar elongation is possible even in the N band (MIR).

Significant polar elongation in emission can be easily achieved at intermediate to large viewing angles, moderate to large ${{ \mathcal N }}_{0}$, and wavelengths between the K and N bands (for small q), with the highest aspect ratios observed in the M band, around 4.8 μm. The highest elongations require σ ≈ 45°–60° and i ≈ 90°.

Figure 7 shows the aspect ratio for many parameter variations, but as heat maps. Parameters i = 75°, q = 0.08, and Y = 18 are held fixed, and the other parameters vary as indicated.

Zoom In Zoom Out Reset image size

Regions of significant polar elongation are clearly visible as areas of deep red color in the panels. The highest realized polar elongation for the NGC 1068 case shown in the figure is e = 2.13 and is realized at $(\sigma ,{{ \mathcal N }}_{0},{\tau }_{{\text{}}V},\lambda )=(45^\circ ,12,160,\sim 4.0\,\mu {\rm{m}})$. The smallest ratio in the figure, i.e., the highest in equatorial direction, is e = 0.39. More extreme ratios require more edge-on viewing angles.

MIR observations with long-baseline interferometry on VLTI (Burtscher et al. 2013, and references therein) measured polar elongations of nuclear dust emission for a number of sources: ∼2 (for Circinus), 1.3 for NGC 3783, 1.5 for NGC 424, and 1.3 for NGC 1068, the example object of our study. These elongations are easily achievable with Clumpy torus images. For instance, with i held at 75°, model image elongation ranges up to 1.5 at 10 μm; it can reach well above 2 at higher inclinations. Thus, Clumpy torus models naturally can reproduce the polar elongation of the spatial flux distribution in the N band up to several parsecs away from the AGN.

Interferometric observations that include shorter baselines are sensitive to the brightness distribution on larger spatial scales, i.e., farther away from the nucleus. López-Gonzaga et al. (2016) reported on such observations with VLTI, adding the AT telescopes, and found even larger polar elongations on such scales for several sources: 1.9 for NGC 424 (and for Circinus), 2.3 for NGC 1068, 2.5 for NGC 5506, and 3.1 for NGC 3783. In addition, Leftley et al. (2018) find e = 2.9 for ESO323-G77. Note that for Circinus and NGC 1068 these values are lower limits, as the extended emission is overresolved interferometrically and continues in single-dish observations. NGC 3783 and ESO323-G77 are type 1 AGNs, meaning that our viewing lines to them are most likely not too highly inclined. As we have seen in Figures 6 and 7, high inclinations are necessary to produce significant polar elongation of the MIR morphology. In these type 1 sources a two-component model is therefore almost certainly preferred. For the above measurements, the authors estimated polar elongation by fitting 2D Gaussians, but the FWHM of a Gaussian is equivalent to the radii of gyration we use here (although our moment-based method works for any morphology).

Achieving elongations of 2.5 or even 3 with single-component models (i.e., a torus) may be challenging, or impossible in type 1 orientations. Two-component models, such as a disk and a dusty outflow, provide the necessary morphology by design, but at the cost of higher model complexity. The measured elongation, however, always depends on the pixels considered. Just as in the case of VLTI measurements with UT baselines only, versus UT+AT baselines, the sensitivity to larger spatial scales changes the derived elongations. Similarly, an instrumental flux sensitivity limit will affect which pixels will in fact be considered when estimating morphological parameters.

Figure 8 illustrates this with two Clumpy torus models, one with parameters from SED fitting in LR18, and a version of it with parameters somewhat modified to increase elongation. The elongation measured as a function of wavelength and sensitivity level (i.e., which contour is the limiting isophote) varies strongly. In the LR18 model, at 10 μm, considering the entire image yields e = 1.2, but it reaches 1.43 for the 7% contour. This SED-derived model is elongated enough to explain the value of e = 1.3 measured by VLTI for NGC 1068 (Burtscher et al. 2013) with long baselines alone. It is not elongated enough to explain e = 2.3 derived for VLTI observations López-Gonzaga et al. (2016) performed on larger spatial scales.

Zoom In Zoom Out Reset image size

Modifying some of the model parameters (e.g., higher σ, i, ${{ \mathcal N }}_{0}$, τ_V ) can produce models with larger polar elongations. In Figure 8 the bottom row shows one such model; it produces e = 1.41 at 10 μm if the entire image is considered. However, the same model shows much higher elongations at other sensitivity levels, e.g., e = 3.22 at the 0.5 level of the peak.

To match observations, not only the elongation but also the angular size of the morphology must be considered. In this release of Hypercat our sampling of the torus size (the Y parameter) tops out at Y = 20. With larger values and sufficiently high luminosity it is possible that such elongations could be maintained on larger physical scales.

Somewhat similar to the question of brightness level sensitivity, a potential pitfall in resolved observations is that the measured elongation of a brightness distribution depends on the FOV. This effect may even invert the axis of elongation, depending on the FOV. Figure 9 demonstrates this with a 4.8 μm image of a Y = 20 torus model, looked at with progressively larger FOVs. From small to large FOV, the measured elongation $e={R}_{g}^{y}/{R}_{g}^{x}$ turns from slightly elongated in the equatorial direction to very elongated along the polar axis, as indicated in the panels.

Zoom In Zoom Out Reset image size

3.2.4. Gini Coefficient as Indicator of Image Compactness

The Gini coefficient/index (Gini 1921) was originally developed in the field of economics to characterize the relative inequality of income distributions. It is applicable generally, though, and is a useful measure of image compactness in our case. Applied to the distribution of pixel values, it measures the relative deviation of the distribution from uniformity. In the discrete case, the Gini coefficient of an array I is

$$\begin{eqnarray}&&G=\displaystyle \frac{\sum _{i}(2i-n-1)\cdot {I}_{i}}{n\sum _{i}{I}_{i}},\end{eqnarray} \tag{ 10 }$$

where the n values I_i of the array are sorted in ascending order, and i are the corresponding (one-based) indices of the sorted array. In our case of a 2D image, I(x, y) can be simply flattened prior to computing G. The Gini index is independent of the absolute magnitude of the pixel values; only their relative distribution matters. If all values are identical, G = 0. If the distribution is randomly uniform, G = 1/3. In the most unequal distribution, a delta function (e.g., when exactly 1 pixel has nonzero value), G = 1.

Figure 10 shows G as a function of several parameters that vary as indicated. Y = 18, q = 0, and i = 90° are fixed.

Zoom In Zoom Out Reset image size

Values close to unity (dark-blue colors) are the most compact morphologies, where all flux is concentrated within a small fraction of all pixels. Lower values (red colors) represent broader distributions. For the Y = 18 torus embedded in a (2 Y_max)² = 40 × 40 FOV, the smallest achievable G value is ∼0.36. This would be the case if all the pixels within a Y = 18 circle had the same value while pixels outside that circle were zero (we are ignoring the central dust-free cavity for the sake of simplicity). The smallest G in Figure 10 is 0.40, i.e., there are models whose appearance approaches that of a uniformly bright disk. Not surprisingly, the corresponding model parameters represent the thickest and densest torus of all shown in the figure, namely, $(\sigma ,{{ \mathcal N }}_{0},{\tau }_{{\text{}}V},\lambda )=(75^\circ ,12,160,18\,\mu {\rm{m}})$. The elongation of this morphology is $e={R}_{g}^{y}/{R}_{g}^{x}=1$. On the other hand, the most compact emission morphology is realized with the shallowest torus, smallest number of clouds, lowest optical depth, and shortest wavelength $(\sigma ,{{ \mathcal N }}_{0},{\tau }_{{\text{}}V},\lambda )=(15^\circ ,1,10,2\,\mu {\rm{m}})$. This smallest model has G = 0.97 and elongation e = 0.31. Figure 18 in Appendix C.6 shows both extremes.

In general, compact images are achieved at small ${{ \mathcal N }}_{0}$, small λ, and preferentially small optical depths. For large optical depths, only the very smallest ${{ \mathcal N }}_{0}$ and shortest wavelengths can produce compact emission sizes, i.e., when only the hottest dust close to the center is visible. In Figure 10, at small to moderate τ_V and moderate to high ${{ \mathcal N }}_{0}$, a size enhancement is clearly visible around 10 μm. This is because the dust cross section peaks locally around 10 μm, increasing the optical depth, thus intercepting more photons, also at greater distances, before they can escape the torus. At wavelengths slightly away from the peak of the dust cross-section curve the photon escape probability is significantly higher (see Figure 2 in Nenkova et al. 2008b).

3.3. Image Asymmetry

The interplay between torus angular width σ, viewing angle i, and total optical depth along an LOS can shift the region that emits light (as viewed by an observer) up or down from the image center. This asymmetry in the distribution of emission can be quantified in several ways.

3.3.1. Thermal Emission Distribution

Probably the simplest measure of asymmetry with respect to the x-axis is the fraction of total light contained in the upper image half. This fraction is by design 0.5 in pole-on and edge-in views, but for other inclination angles it will depend on some of the model parameters. Figure 11 shows this for a number of typical model parameter combinations and four wavelengths between 2.2 and 18.5 μm.

Zoom In Zoom Out Reset image size

The asymmetry is largest at short wavelengths (K band), as this is mostly scattered light from the far inner side of the torus. Longer wavelengths can penetrate deeper into the torus in all directions, diminishing the asymmetry in light distribution among the upper and lower image halves. The fraction as a function of viewing angle always has a clear peak, and for all σ > 45° and short wavelengths it is approximately symmetric around i ≈ 45°. At smaller σ or longer wavelengths the peak occurs roughly at 90° − σ, i.e., at the torus half-opening angle. Increasing optical depth via higher ${{ \mathcal N }}_{0}$ increases the overall scale of the emission asymmetry.

The fraction of emission from regions above the nucleus can easily reach 90% or more at 2.2 and 4.8 μm, and over 80% at 10 μm. At 18.5 μm the asymmetry tops out at ∼70%, but it is significantly less in most cases.

In Section 3.4 we will demonstrate that although the distribution of the thermal emission can be highly asymmetric in many configurations, it does not necessarily correlate with the physical distribution of the dust in the torus, whose bulk is still in the equatorial plane.

3.3.2. Centroid Offset

A different method to quantify the asymmetry in emission distribution is by measuring the y position of the brightness centroid, $\bar{y}$. Figure 12 shows this in a similar fashion to Figure 11, but the y-axis in each panel is the vertical offset of the light centroid from y = 0 in units of dust sublimation radius R_d .

Zoom In Zoom Out Reset image size

In K through M bands these curves always peak close to edge-on viewings, but not quite at i = 90° of course, as this is a symmetric morphology. The N band behaves similarly, albeit transitioning to the longer-wavelength case, where the Q band curve shows a much more symmetric shape around i ≈ 45°. The 2.2 and 4.8 μm images show the largest centroid shifts (of very similar magnitude) and can be quite substantial, up to 5 R_d , i.e., almost 30% of the torus overall radius in this Y = 18 R_d case. The centroid offsets at 10.0 μm are a bit smaller, and at 18.5 μm they are the smallest of all. Again, this is due to long-wavelength photons traveling more freely in all directions and into the torus, thus reducing the asymmetry. The dependency of $\bar{y}$ on ${{ \mathcal N }}_{0}$ is stronger than that of the emission fraction in the upper image half discussed in the previous subsection.

3.3.3. Skewness

Using the image moments developed earlier, the asymmetry of emission distribution can be quantified by measuring skewness. The skewness of a distribution of pixel values, projected onto the axes, is a third-order moments function:

$$\begin{eqnarray}&&{S}_{x}={\mu }_{03}/\sqrt{{\mu }_{02}^{3}},\qquad {S}_{y}={\mu }_{30}/\sqrt{{\mu }_{20}^{3}}.\end{eqnarray} \tag{ 11 }$$

It measures the deviation from a symmetric distribution in the given direction. Since Clumpy images are always symmetric about the y-axis, the skewness S_y is always zero. We thus only consider S_x , the image skewness about the x-axis (i.e., in the y-direction). Negative values of skewness indicate asymmetry toward the positive direction of an axis. Note that a zero-valued skewness does not guarantee symmetry, but a symmetric distribution will always have skewness zero.

Figure 13 plots S_x for Clumpy model images at three wavelengths and as a function of several torus model parameters that vary as indicated. For all but the very smallest value of ${{ \mathcal N }}_{0}$, the skewness curves appear very similar in shape for a given σ and wavelength. The overall amplitude of the curve depends very strongly on ${{ \mathcal N }}_{0}$. The wavelength determines the degree of oscillations in the curves, with the shortest wavelength showing the most variation with viewing angle (we did not plot the curve for the K band since its range is much larger than that of longer wavelengths). At short wavelengths (4.8 μm) the skewness curve increases first with growing viewing angle i (i.e., the skewness is toward the negative y-axis). This is due to the lateral "horns" of the dust-free cavity in the brightness map protruding below the image x-axis. As the viewing angle continues to increase, the curve then quickly plunges to negative values, indicating skewness of the image toward the positive y-axis. The viewing angle i at which this turnover occurs depends on the torus angular width σ. Around 10 μm very little evidence remains of positive skewness values, and the negative excursions are of smaller amplitude than at 4.8 μm. At long wavelengths (18.5 μm and above) almost no skewness is perceivable, except at the smallest σ; the images have almost no skew, regardless of viewing angle.

Zoom In Zoom Out Reset image size

3.4. Morphologies of Projected Dust Maps

The Hypercat model cubes also contain 2D maps of the projected dust cloud distribution in the plane of the sky, C_d (x, y). These are maps of the number of clouds along each LOS (pixel). Whenever the LOSs through the torus are optically thick, the morphology of the torus emission cannot be expected to follow the dust cloud distribution. Figure 14 shows this clearly for one model, at three viewing angles and a number of wavelengths. The leftmost column shows the underlying dust cloud distribution, with labeled contours of constant number of clouds along an LOS overplotted. Only at FIR wavelengths above ≈18 μm does the torus become globally optically thin, and the dust emission morphologies densely fill out the cloud number contours.

Zoom In Zoom Out Reset image size

The methods developed earlier for brightness maps can be used to quantify the dust morphologies. We remind the reader that (i) the dust maps do not depend on τ_V ; (ii) they depend on ${{ \mathcal N }}_{0}$ only as a multiplicative factor, i.e., ${C}_{d}({{ \mathcal N }}_{0})={{ \mathcal N }}_{0}\times {C}_{d}({{ \mathcal N }}_{0}=1)$; and (iii) moment-based morphological measurements on these dust maps, such as radii of gyration and elongations, are independent of ${{ \mathcal N }}_{0}$ as well. They only depend on σ, i, Y, and q, i.e., on the truly geometrical parameters. We can thus compare the extension of the thermal emission maps with the extension of the dust distribution map. Figure 15 shows as orange lines the morphological elongation of the 2D dust maps for a number of model parameter combinations. σ (columns) and λ (rows) vary from panel to panel. In every column the orange curve is the same, since the dust distribution does not depend on wavelength. In all cases the elongation of the dust distribution is in the equatorial direction (i.e., ratios smaller than 1:1), which is expected given the range of σ values sampled.

Zoom In Zoom Out Reset image size

The other lines show the ratio of this elongation of the emission to the elongation of the underlying dust morphology, at four different values of ${{ \mathcal N }}_{0}$. In almost all cases the elongation of the emission distribution is more polar than the elongation of the associated dust distribution. The only marginal exception is at the smallest σ and smallest ${{ \mathcal N }}_{0}$. The emission-to-dust elongation ratios are higher at shorter wavelengths. At large σ and wavelengths longer than ∼12 μm the elongations of the emission and dust morphologies are quite close to each other.

Figure 16 demonstrates spatial 2D maps of the ratio of local emission to dust content. The dust map is the same for all panels and corresponds to the best-fit NGC 1068 model as before, with σ = 43°, i = 75°, Y = 18, q = 0.08, and τ_V = 70. ${{ \mathcal N }}_{0}$ and λ vary between panels as indicated. Each panel shows the ratio (computed per pixel) of the light emission map and the generating dust map, which both were first normalized to unit sum. That is, both are fractional maps of the total emission and total dust cloud number. The scaling of the ratio maps is symmetrically logarithmic. A single black contour line divides regions where the fractional emission dominates over the fractional dust content (blue colors) from those regions where there is fractionally more dust present than light produced (red colors).

Zoom In Zoom Out Reset image size

These maps give quick insight of the interplay between the global morphology of photons generated by the distribution of dust. For our model with viewing angle i = 75°, i.e., tilted by 15° toward the observer, it is immediately clear that at MIR wavelength the emission has the upper hand over dust contribution in regions high above the equatorial plane of the torus, along the system axis. When ${{ \mathcal N }}_{0}$ is very small, at nearly all wavelengths the emission is concentrated to a small circular region around the image center. The same is true at any value of ${{ \mathcal N }}_{0}$ when the wavelength is very long (submillimeter).

To date only VLTI and ALMA, both interferometers, have been able to spatially resolve the IR emission from nuclear dusty environments in some nearby AGNs. In the future a number of extremely large single-dish telescopes will be pointing their 24–39 m diameter (GMT, TMT, ELT) primary mirrors at AGNs as well. They will resolve the IR emission in these sources, providing answers about the morphologies of light emission and dust distribution that are free from model assumptions. Resolved and model-free imagery is the only technique capable of breaking degeneracies intrinsic to IR radiative transfer that can produce similar SEDs from different dust geometries.

A major observational finding in recent years was that in several nearby AGNs the MIR emission resolved with VLTI appears elongated along the polar direction of the system. That is, the MIR light emanates from regions above the equatorial plane of the AGN torus. Several attempts at explaining these observations were proposed, some of which significantly deviate from simple AGN unification models.

We developed techniques based on image moments to quantify image morphology, and we show that thermal emission maps produced by simple "classical" toroidal models such as Clumpy can generate significant polar elongation in the MIR. Typical polar-to-equatorial size ratios of up to 2:1 found by interferometric observations on physical scales of several parsecs are easily achievable by Clumpy torus models and can be even larger depending on sensitivity levels. We constrain the ranges of parameters that produce such elongations and find that they are rather standard. A slight tilt of the torus system axis toward the observer is enough to produce elongations compatible with many observations. We further find that it is impossible to obtain any polar elongation at small values of σ and i, at smallest values of ${{ \mathcal N }}_{0}$, or at wavelengths longer than ∼14–18 μm. If the radial cloud distribution is compact (higher q values), polar elongations are only achievable up to ∼5 μm.

While the clumpy distributions that we model here fundamentally provide the optically and geometrically thick material that AGN unification calls for, and alone can account for much of the observed polar elongation, additional components may still be present, especially on larger physical scales (tens of parsecs). For example, extended winds (Hönig & Kishimoto 2017; Izumi et al. 2018; Asmus 2019) and photoionization of existing material may also contribute, including at MIR wavelengths. In some sources, such as NGC 1068, MIR interferometric data suggest that there are additional components, e.g., from a putative dusty outflow, which are extended on larger scales, or even disassociated (López-Gonzaga et al. 2014); outflows on even larger scales, and crossing to the galactic domain, are seen in submillimeter, molecular gas observations (e.g., Aalto et al. 2020, and references therein). These very extended components may not be explainable with a torus model alone, and their physical connection with the torus, or lack thereof, needs to be studied.

Evidence from radiation-hydrodynamical simulations also strongly suggests that two-component models may be preferable, especially on larger scales, where outflowing winds and "failed fountains" may provide a good fraction of the observed MIR emission on scales of tens of parsecs and be elongated in polar directions (e.g., Schartmann et al. 2014; Wada et al. 2016; Williamson et al. 2020). Note that only large dust grains could survive in such polar regions (e.g., Tazaki & Ichikawa 2020). The full picture likely requires magnetohydrodynamics as well (Blandford & Payne 1982; Emmering et al. 1992) given the recent direct detection of B-fields along the equatorial direction of the dusty torus in NGC 1068 (Lopez-Rodriguez et al. 2020). These measurements of magnetically aligned dust grains using ALMA polarization at 860 μm support the picture that magnetic fields up to a few parsecs contribute to the accretion flow onto active nuclei.

We also compared the emission morphologies to their generative projected (2D) dust maps. A key result is that the emission morphology does not simply trace the dust distribution in the system. Optical depth effects are important at most wavelengths, including in the MIR. In the FIR and submillimeter regimes the torus is mostly optically thin, even for a high number of clouds along the LOS. There, the observed emission does trace the underlying dust distribution.

To quantify morphologies, we advocate the use of moment-based techniques (e.g., radii of gyration) over more traditional methods (e.g., half-light radius). Moment-based morphological quantities are independent of the image brightness. Thus, in our case, maps of the dust cloud distribution are independent of ${{ \mathcal N }}_{0}$ and τ_V , since these are just multiplicative factors, and the brightness maps are independent of the total flux. Furthermore, morphological size measurements using radii of gyration provide size estimates along both axes independently, making it trivial to measure elongations. Moment-based quantities can estimate higher-order morphological quantities such as image asymmetry, skewness, and kurtosis.

In anticipation of the extremely large telescopes currently planned or under construction (GMT, TMT, ELT), and to facilitate studies of image morphology, we have introduced with Hypercat a hypercube of clumpy torus images spanning a wide range of wavelengths and model parameters. The sampled wavelengths include important NIR, MIR, FIR, and submillimeter bands. Hypercat also comprises a software suite whose key functionality is to interpolate the images continuously at any set of model parameters. Another capability of Hypercat is to simulate several observing modes, e.g., single-dish observations with PSF convolution, ¹³ detector pixelization, noise modeling, and image deconvolution to recover some of the original signal in the images. The processing software may also be applied to analyze different emission models, which need not be limited to AGN tori.

In the companion Paper II we apply the capabilities of the Hypercat software to the Clumpy image hypercube and simulate observations of NGC 1068 with several single-dish telescopes (JWST, Keck, GMT, TMT, ELT) to investigate the spatial resolvability of the NIR to MIR emission. We also simulate IFU-like observations to analyze the 10 μm silicate emission feature. Finally, we fit directly with the Clumpy images the K-band image of NGC 1068 recently published by G20, deriving likely model parameters.

The hypercubes of Clumpy images and cloud distribution maps (as a function of all model parameters and wavelengths) are provided to the community through FTP download. See Appendix A, and the User Manual, for download instructions.

We also provide the Hypercat software, which comprises user-friendly tools to operate on the hypercubes, simulate observations, and investigate morphology. All codes, scripts, and telescope pupil images can be found in the GitHub repository at https://github.com/rnikutta/hypercat/.

The repository also contains the User Manual, several instructional Jupyter notebooks, and API documentation. The User Manual contains high-level usage examples and also showcases low-level operations that can be performed with Hypercat.

We kindly ask for citation of this paper if the reader decides to make use of any of the Hypercat data cubes, software functions, scripts, or telescope pupil images.

We wish to thank Leonard Burtscher, Konrad Tristram, Marko Stalevski, Moshe Elitzur, Ric Davies, and Tanio Diaz Santos for illuminating discussions on the subjects of this paper. We are thankful to the referee, whose comments helped improve the manuscript. R.N. acknowledges early support by FONDECYT grant No. 3140436. E.L.-R. acknowledges support from the Japanese Society for the Promotion of Science (JSPS) through award PE17783, the National Observatory of Japan (NAOJ) at Mitaka, and the Thirty Meter Telescope (TMT) Office at NAOJ-Mitaka for providing a space to work and great collaborations during the short stay in Japan. K.I. acknowledges support by the Program for Establishing a Consortium for the Development of Human Resources in Science and Technology, Japan Science and Technology Agency (JST), and partial support by the Japan Society for the Promotion of Science (JSPS) KAKENHI (20H01939; K. Ichikawa). C.P. acknowledges support from NSF grant No. 1616828. S.F.H. acknowledges support by the EU Horizon 2020 framework program via the ERC Starting grant DUST-IN-THE-WIND (ERC-2015-StG-677117). A.A.-H. acknowledges support through grant PGC2018-094671-B-I00 (MCIU/AEI/FEDER,UE). A.A.-H. work was done under project No. MDM-2017-0737 Unidad de Excelencia "María de Maeztu" - Centro de Astrobiología (INTA-CSIC). R.N., E.L.-R., K.I. are very grateful to NOAO (now part of NSF's NOIRLab), SOFIA Science Center, and the Program for Establishing a Consortium for the Development of Human Resources in Science and Technology, Japan Science and Technology Agency (JST), for providing travel grants that made three on-site project meetings possible. We are grateful to colleagues at several telescope collaborations and consortia who were willing and able to provide us with the latest versions of their pupil images. These are, in order of increasing telescope diameter: Marshall Perrin (JWST), Andrew Skemer (Keck), Warren Skidmore and Christophe Dumas (TMT), Rebecca Bernstein (GMT), and Suzanne Ramsey (ELT). With their permission we are here publishing these pupil images as FITS files (see Section 5, and the project repository https://github.com/rnikutta/hypercat/).

Facilities: JWST - James Webb Space Telescope, Keck - , VLTI (GRAVITY). -

Software: Hypercat (this work), astropy (Astropy Collaboration et al. 2013; Price-Whelan et al. 2018), h5py (Collette 2019), matplotlib (Hunter 2007), numpy (Harris et al. 2020), r2py (Gautier 2014), scipy (Virtanen et al. 2020), urwid (Ward et al. 2014).

Readers interested in using the Clumpy hypercubes will need to download one or more hdf5 files from https://www.clumpy.org/images/, or from the current FTP location: ftp://ftp.tuc.noirlab.edu/pub/nikutta/hypercat/. Table 3 lists the available options. The ∗_all.hdf5.gz file contains the image hypercube at all sampled wavelengths. This is the maximally compressed version of the hdf5 file, which must be uncompressed on the user's computer system. To reduce the peak storage required on the target computer, both steps can be executed in one go (all commands in a single line):

lftp -e "set net:timeout 10; cat /pub/nikutta/hypercat/hypercat_20200830_all.hdf5.gz; bye" ftp.tuc.noirlab.edu ∣ gunzip > hypercat_20200830_all.hdf5

The program lftp must be installed on the target system, and 914 GB of space must be available on it (but only 271 GB of compressed data will be downloaded). The reader should also download the checksums file ftp://ftp.tuc.noirlab.edu/pub/nikutta/hypercat/hypercat_20200830.md5 and verify the hypercube file:

# this can take 30 minutes even on a modern computer

md5sum ----ignore-missing -c hypercat_20200830.md5

hypercat_20200830_all.hdf5: OK

# or on MacOS and BSD variants

md5 hypercat_20200830_all.hdf5

#... and compare the printed hash with the one in the .md5 file

The download was successful if the computed checksum is identical to the one in the MD5 sums file. The same procedure applies to the other ∗hdf5.gz files analogously (which hold subsets of the wavelength sampling).

The Hypercat model images are embedded within the same FOV, corresponding to the largest sampled Y value, ${Y}_{\max }$. The total number of pixels in all stored brightness maps is therefore ${N}_{\mathrm{pix}}^{\mathrm{img}}={N}_{x}\,{N}_{y}\times {\prod }_{k}{N}_{{\theta }_{k}}$, where ${N}_{{\theta }_{k}}$ is the number of sampled values per model parameter θ_k , $k\in \{\sigma ,i,Y,{{ \mathcal N }}_{0},q,{\tau }_{{\text{}}V},\lambda \}$, as compiled in Table 1. Thus, ${\prod }_{k}{N}_{{\theta }_{k}}$ is the total number of realized parameter combinations. If full images are stored, the total number of pixels in the hypercube is ${N}_{\mathrm{pix}}^{\mathrm{img}}={\left(2\,{Y}_{\max }\,\eta +1\right)}^{2}\times {\prod }_{k}{N}_{{\theta }_{k}}$. If half-images are stored, as is the case in Hypercat because of image symmetry, then ${N}_{\mathrm{pix}}^{\mathrm{img}}=\left({Y}_{\max }\,\eta +1\right)\left(2\,{Y}_{\max }\,\eta +1\right)\,\times {\prod }_{k}{N}_{{\theta }_{k}}$. That is, the number of computed pixels grows quadratically with ${Y}_{\max }$, quadratically with the image resolution per unit linear size η, and linearly with the sampling of any model parameter including wavelength. Clumpy also computes integrated dust maps during ray-tracing, and with the parameter sampling of our grid, 4000 such maps were produced (see Sections 2.6 and 2.7). With half-image ray-tracing and η = 6, the final number of stored pixels (images + cloud maps) is then N_pix ^tot = N_pix ^img + N_pix ^cld ≈ 2.45 × 10¹¹ pixels, which at single precision (32 bits = 4 bytes) requires ≈913 GB of storage. We deliver this hypercube of half-images as a compressed HDF5 file (.hdf5.gz) of 270 GB size, ¹⁴ and we have also produced four smaller HDF5 files with partial wavelength coverage. Table 3 lists their properties.

Such file sizes can be handled relatively easily with regard to network download and possibly even storage on a researcher's personal workstation, but they are 1–3 orders of magnitude too large to load the entire set into memory (for instance, for the purpose of n-dimensional interpolation). We must therefore devise routines for slicing the hypercube along any parameter axes of interest, and for on-the-fly access to any image map within the envelope of the sampled parameter space (see User Manual for details). Together with the hypercubes, we provide functions to transparently extract full-sized (square) images via n-dimensional interpolation on the hypercube and (automatic) left-right mirroring about the y-axis.

The Hierarchical Data Format version 5 (HDF5) is well established for packaging structured and heterogeneous data sets and metadata, and it is agnostic about the operating system. We use the Python package h5py to construct our hypercubes, and the Hypercat software uses h5py to access said data. HDF implementations in other languages are readily available. ¹⁵ The structure of every Hypercat HDF5 file is simple. Listing 1 shows the hierarchy tree of the main HDF5 file released with this paper.

Listing 1. 

There are Nhypercube = 2 hypercubes in the file. Their names are imgdata and clddata, each stored in a separate HDF group. Group imgdata contains the thermal emission maps, stored in a nine-dimensional ($\sigma ,i,Y,{{ \mathcal N }}_{0},q,{\tau }_{{\text{}}V},\lambda ,x,y$) array imgdata/hypercube. Several other small data sets of auxiliary information are also stored inside the group, e.g., the mapping of parameter values to the hypercube grid vertices (theta). The hypercube shape is given in hypercubeshape, and every array axis corresponds to one list in theta and to one parameter name in paramnames. theta in the listing above shows the sampling of the parameters sig and i, which are the torus angular thickness σ and the viewing angle i. Other parameter samplings are truncated here for brevity. The clddata group does the same for the 2D maps of cloud number along the LOS. These maps do not depend on wavelength and cloud optical depth τ_V , and ${{ \mathcal N }}_{0}$ is just a multiplicative factor (see Section 2.6); thus, their hypercube is six-dimensional (σ, i, Y, q, x, y).

Various methods exist to quantify image morphology. One very powerful approach, due to its generality, is the framework of image moments. For digitized images I(x, y) with x and y pixel indices, the (p + q)-order raw or geometric moment is defined as

$$\begin{eqnarray}&&{M}_{\mathrm{pq}}=\sum _{x}\sum _{y}I(x,y){x}^{p}\,{y}^{q},\end{eqnarray} \tag{ C1 }$$

with a standard power basis. Other bases, e.g., orthogonal functions, Zernike polynomials, etc., are possible and are being used in numerous applications. For binary images the raw moment M₀₀ is the area (sum) of the one-valued pixels. For grayscale images, in the astronomical context, it is the total light contained in all pixels.

Figure 17 demonstrates the effects of various image operations (using 2D Gaussians as images) on the measured moments, from left to right and top to bottom: (a) translation, (b) stretching, (c) scaling, (d) and (e) skewing, and (f) rotation.

Zoom In Zoom Out Reset image size

Panel (a) shows a 2D Gaussian with σ_x = 5, σ_y = 10 pixels (in a 101 × 101 pixel frame) being shifted from y₀ = −20 through +20 pixels, and the corresponding centroid locations $\bar{x}$ (not changing) and $\bar{y}$ (changing).

Panel (b), left axis, shows how ${R}_{g}^{y}$ grows when the 2D Gaussian is being stretched vertically, and panel (c), left axis, shows how ${R}_{g}^{x}$ and ${R}_{g}^{y}$ both increase linearly with the increasing scale of a 2D Gaussian. This linear behavior makes radii of gyration very well suited for measuring the characteristic sizes of the pixel distribution.

Panel (b), right axis, also shows how elongation is measured, using a 2D Gaussian as an example whose x-size σ_x is 9 pixels always, while σ_y is varied from 3 to 15 pixels. That is, the Gaussian spans a range of aspect ratios from 1:3 to 5:3. The measured elongation (solid line) grows linearly with σ_y .

Panels (d) and (e) show how skewness can be measured using Equations (11). For this demonstration we use the 2D version (i.e., k = 2) of the multivariate skew-normal distribution as defined by Azzalini & Dalla Valle (1996), with density function

$$\begin{eqnarray}&&f(z)=2\,{\phi }_{k}(z;{\rm{\Omega }}){\rm{\Phi }}({\alpha }^{T}z),\quad (z\in {{\mathbb{R}}}^{k}),\end{eqnarray} \tag{ C2 }$$

where ϕ( · ) is the k-dimensional normal density with zero mean and correlation matrix Ω, Φ( · ) is the N(0,1) distribution function, and α is a k-dimensional vector of skewness coefficients. In our case, α = (α_x = 0, α_y ). ¹⁶ In panel (d) a symmetric 2D Gaussian is modified by changing the skewness parameter α_y . The skewness S_x about the x-axis, i.e., in the y-direction, changes as a result, while S_y remains constant. In panel (e) the same operations are applied to an initially elongated 2D Gaussian; the resulting S_x curve is identical, i.e., the skewness about an axis is independent of the distribution of pixels with respect to the other axis.

C.1. More on Image Moments

Versions of central moments that are invariant to both translation and scale operations are possible by properly normalizing the central moments

$$\begin{eqnarray}&&{\eta }_{{pq}}=\displaystyle \frac{{\mu }_{{pq}}}{{\mu }_{00}^{\left(1+\tfrac{p+q}{2}\right)}}\qquad \forall p+q\geqslant 2.\end{eqnarray} \tag{ C3 }$$

Scale invariance in this context means similitude, which is the resizing of a 2D shape while preserving its proportions. Invariants simultaneously under translation, scaling, and rotation can be defined as combinations of moments (Hu 1962; Reiss 1991; Flusser & Suk 1993). For our applications we do not need to concern ourselves with rotation, since we are always able to set the positional angle PA = 0 deg in our images. Some useful observations on moments:

• (a)  
  Zero-order central and raw moments are identical: μ₀₀ = M₀₀.
• (b)  
  First-order central moments are zero: μ₁₀ = μ₀₁ = 0.
• (c)  
  Second-order central moments μ₂₀ and μ₀₂ measure the distribution of mass w.r.t. the image axes, i.e., can be used to measure source elongation and orientation (see Section 3.2.3).
• (d)  
  Third- and fourth-order central moments can measure higher-order image properties such as skew/asymmetry (see Section 3.3) and kurtosis/peakedness (see Section C.3).

Note that some of the moment-based functionality is available in the photutils package for Python. Unfortunately, at the time of writing, there exist several incompatibilities between Astropy's models package, the scikit-image module, and photutils. For that reason, we rely on our own moment-based implementation. It is straightforward to use, is performant, and does not have many dependencies. This morphology module is part of the Hypercat source code distribution.

C.2. Image Elongation via Covariance Matrix

An equivalent method exists to measure the elongation e of a morphology using the eigenvalues of the image covariance matrix. The eigenvalues and eigenvectors of the covariance matrix Σ of an image I encode the magnitude and direction of the principal components of I, i.e., the 2D extension of the pixel distribution. In our case we can always set the position angle to PA = 0 deg, i.e., we are only interested in elongations along the orthogonal e_x and e_y directions of the image. Thus, the eigenvalues measure the magnitude of the pixel distribution along the x- and y-directions.

The covariance matrix Σ can be computed in several ways, e.g., from second-order central moments

$$\begin{eqnarray}{\rm{\Sigma }}\equiv \mathrm{Cov}[I(x,y)]=\left[\begin{array}{cc}{\mu }_{20}^{{\prime} } & {\mu }_{11}^{{\prime} }\\ {\mu }_{11}^{{\prime} } & {\mu }_{02}^{{\prime} }\end{array}\right],\end{eqnarray} \tag{ C4 }$$

where ${\mu }_{\mathrm{pq}}^{{\prime} }={\mu }_{\mathrm{pq}}/{\mu }_{00}$. The eigenvalues of Σ can be obtained analytically, or conveniently by computer packages (e.g., linalg.eigvals in numpy). The elongation of the intensity distribution in I(x, y) is then the square root of their ratio, i.e., $e=\sqrt{{\lambda }_{y}/{\lambda }_{x}}$. The same elongation can also be computed using the gyration radii as $e={R}_{g}^{y}/\,{R}_{g}^{x}$. Thus, ${\lambda }_{x}={({R}_{g}^{x})}^{2}$ and ${\lambda }_{y}={({R}_{g}^{y})}^{2}$.

C.3. Kurtosis

Kurtosis is a fourth-order measure that expresses the "peakedness" of a distribution

$$\begin{eqnarray}&&{K}_{x}=\displaystyle \frac{{\mu }_{40}}{{\mu }_{20}^{2}}-3,\qquad {K}_{y}=\displaystyle \frac{{\mu }_{04}}{{\mu }_{02}^{2}}-3.\end{eqnarray} \tag{ C5 }$$

In the context of image morphology, it can be used to quantify the compactness of the emission pattern. K_x,y are zero for a Gaussian distribution, greater than zero for distributions more peaked than Gaussians, and less than zero for flatter distributions.

C.4. Image Orientation

C.5. Half-light Radius

A circular aperture with half-light radius R_HL contains half of the flux in an image I, such that (see, e.g., Burtscher et al. 2013)

$$\begin{eqnarray}&&\displaystyle \frac{1}{{F}_{\mathrm{tot}}}{\int }_{0}^{{R}_{\mathrm{HL}}}{\text{}}{dr}\,I\,2\pi r=\displaystyle \frac{1}{2},\end{eqnarray} \tag{ C6 }$$

where R_HL is the root of this equation. The solution is not unique if the circle need not be centered in the image.

C.6. Gini Coefficient

Figure 18 shows the smallest and largest morphologies within the parameter space from Figure 10.

Zoom In Zoom Out Reset image size