Hypercubes of AGN Tori (HYPERCAT). II. Resolving the Torus with Extremely Large Telescopes

The type 2 Seyfert galaxy NGC 1068 has been extensively studied in the IR using single-dish (e.g., Mason et al. 2006; Alonso-Herrero et al. 2011; Ramos Almeida et al. 2011; Lira et al. 2013; Ichikawa et al. 2015) and interferometric (e.g., Jaffe et al. 2004; Raban et al. 2009; Burtscher et al. 2013; López-Gonzaga et al. 2014) observations. ALMA submillimeter observations have recently resolved the molecular torus in NGC 1068 and constrained its diameter to ∼12 pc using the dust continuum at 432 μm and CO J = 6-5 (García-Burillo et al. 2016), and then found a size 13 × 4 pc in the HCN J = 2-1 transition and 12 × 5 pc for HCO⁺ J = 3-2 (Imanishi et al. 2018). The torus is found to be inclined, and lopsided, at a position angle of its midplane about 100–110° east of north. Including the 432 μm dust continuum observations by ALMA as a constraint in the 1–432 μm nuclear SED fits of NGC 1068, Lopez-Rodriguez et al. (2018, hereafter LR18) derived the Clumpy torus parameters that best reproduced the observations. We adopt these model parameters for our present study; they are listed in Table 1. The position angle of the torus large axis in the plane of the sky, PA ∼138° east of north, was taken from 2 μm polarimetric observations of the nucleus (Packham et al. 1997; Simpson et al. 2002; Gratadour et al. 2015; Lopez-Rodriguez et al. 2015). The 2 μm polarization is thought to arise from the absorption of aligned dust grains in the torus of NGC 1068 where the measured polarization angle indicates the global orientation of the equatorial plane of the torus. This PA has been further confirmed by means of magnetically aligned dust grains across the resolved equatorial axis of the torus using ALMA polarimetric 860 μm observations (Lopez-Rodriguez et al. 2020).

Table 1. Clumpy Parameters and Derived Quantities from the Best-fit SED Model to NGC 1068 by LR18

Model Parameters		Derived & Literature Values
σ (deg)	${43}_{-15}^{+12}$	Rin (pc)	${0.28}_{-0.01}^{+0.01}$
Y	${18}_{-1}^{+1}$	Rout (pc)	${5.1}_{-0.4}^{+0.4}$
${{ \mathcal N }}_{0}$	${4}_{-1}^{+2}$	H (pc)	${3.5}_{-1.3}^{+1.0}$
q	${0.08}_{-0.06}^{+0.19}$	Lbol/1044 (erg s−1)	${5.02}_{-0.19}^{+0.15}$
τV	${70}_{-14}^{+6}$	D a (Mpc)	14.4
i (deg)	${75}_{-4}^{+8}$	PA (deg)	138

Note.

^a At D = 14.4 Mpc, 1'' = 70 pc, adopting H₀ = 73 km s⁻¹ Mpc⁻¹ (Bland-Hawthorn et al. 1997).

Download table as:  ASCIITypeset image

Figure 1 shows the distribution of the thermal emission within the 1–432 μm wavelength range for a Clumpy model with parameters from Table 1. The morphology of the emission changes dramatically across the wavelength regimes as also shown in Figure 5 in LR18. For details on emission morphology please see the companion paper (Paper I).

Zoom In Zoom Out Reset image size

2.2.1. Telescopes, Pupil Images, and PSFs

To produce synthetic observations with a given telescope and instrument combination, the telescope PSF and the instrumental configuration are needed. A telescope PSF can be very complex (e.g., Krist et al. 2011; Perrin et al. 2012) and depends on many factors; for instance, telescope aperture, telescope configuration, seeing conditions, adaptive optics configuration, and wavelength.

Our aim is to provide the best-case scenario of a synthetic observation for a given telescope and wavelength. We assume perfect adaptive optics, and atmospheric conditions without sky or PSF variations. These conditions can be refined by applying specific aberrations relevant to the science cases and telescope/instrument configuration, but this is outside the scope of the present work ¹² . Future releases of Hypercat may incorporate tools to account for these conditions.

Table 2 summarizes the telescopes and instruments available in Hypercat, which include the next generation of 30 m telescopes (GMT, TMT, and ELT) as well as JWST and Keck for comparison. For telescopes without an instrument in a given wavelength range the pixel scale is taken to be equal to the Nyquist sampling (0.5 × λ/D). The FWHM is given by λ/D, except for JWST, whose PSFs have been modeled in detail by the instrument teams.

Table 2. List of Single-dish Telescopes and Instruments Available with Hypercat

Telescope	Diameter	Instrument	Wavelength	Pixel scale	FWHM
	(m)		(μm)	(mas)	(mas)
JWST	6.5	NIRCAM	0.6–5	31–63	31–162 b
		NIRISS	0.6–5	31–63	42–152 c
		MIRI	5–28	110	182–812 d
Keck	10.0	NIRC2	0.9–5.3	10–40	163–241 e
		Nyquist a	8–20	83–206	165–412
GMT	24.5	GMTIFS	0.9–2.5	5	8–21
		Nyquist a	3–5	12–21	25–42
		Nyquist a	8–20	34–84	67–168
TMT	30.0	IRIS	0.85–2.5	4	6–17
		MICHI	3–5	11.9	21–34
		MICHI	8–12	27.5	55–82
ELT	39.1	MICADO	0.8–2.5	4	4–13
		METIS	2.9–5	4	15–26
		METIS	5–14	13	26–74

Notes.

^a For telescopes without an instrument in a given wavelength range, the pixel scale is set to the Nyquist sampling at the shortest wavelength. For FWHM values without a superscript, λ/D was assumed. Otherwise the FWHM values are from: ^b https://jwst-docs.stsci.edu/near-infrared-camera/nircam-predicted-performance/nircam-point-spread-functions ^c https://jwst-docs.stsci.edu/near-infrared-imager-and-slitless-spectrograph/niriss-predicted-performance/niriss-point-spread-functions ^d https://jwst-docs.stsci.edu/mid-infrared-instrument/miri-predicted-performance/miri-point-spread-functions ^e https://www2.keck.hawaii.edu/inst/nirc2/ScamCompare.html

Download table as:  ASCIITypeset image

We compute the PSF as the Fourier transform of the telescope pupil image (see Appendix A.1 for details). The PSF includes the effects of the telescope structures that interact with the science beam. The pupil images are in units of meters, which are then converted to arcsec using the PIXSCALE keyword in the corresponding FITS files. For any given wavelength and telescope, the PSF has a FWHM = λ/D. Figure 2 shows the pupils and PSFs of the telescopes from Table 2.

Zoom In Zoom Out Reset image size

Note that Hypercat also provides two other methods to estimate the PSF: model-PSF, which is a combination of a central Airy disk that contains most of the power of the PSF, and a broad Gaussian halo that encompasses the remaining power of the object; and image-PSF where the user can provide as an image the empirical PSF of an observed standard star for a given telescope/instrument. model-PSF also accounts for Strehl ratio to perform studies of the AO system (Appendix A.1).

To qualitatively explore the PSF for each instrument at a given wavelength we produce monochromatic PSF images (Figure 2) at five representative wavelengths within the 1−12 μm wavelength range for each of the telescopes in Table 2. Note that the throughput of a specific filter may be taken into account for a more realistic PSF. Here, the only diffractive and optical aberration component in our simulations comes from the pupil image.

Figure 3 shows the radial profiles (computed over circular annuli) of the PSFs of GMT, TMT, and ELT at 2.0 μm. As expected, the cleanest (most Airy-like) extended power of the PSF is provided by the GMT, while the highest-resolution core of the PSF is provided by the ELT. As shown by Angel et al. (2003), the best extended PSF will be provided by the telescope configuration with full disks, i.e., GMT, over those with fitted hexagons, i.e., ELT and TMT. This is because the extended regions of the PSF from fitted hexagons have more structures (i.e., spider shape in the wings of the PSF) than those from full circular disks. The instrumental configuration will be thus of great importance when considering specific science goals of the ground-based 30 m class telescopes.

Zoom In Zoom Out Reset image size

2.2.2. Toward a Synthetic Observation

Single-dish observations are degraded by several processes (Appendix A). The fundamental loss of resolution is due to convolution of the incoming wavefront with the PSF of the telescope. The PSF-convolved image is then downsampled to the discrete pixel sizes of the detectors. The PSF-convolved and pixelated image (step E in Figure 4) is then flux calibrated. Specifically, the total flux of the model is equal to the flux density, in units of jansky, from a known observation at the highest spatial resolution at a given wavelength. For our example, we used a flux density of 1 Jy at 3.45 μm by LR18. Finally, noise due to, for instance, fluctuations in detector readout, atmospheric emission, and transmission is taken into account (Appendix A.4). Note that the final synthetic observations have a background set to zero, background subtraction is not required. Figure 4 shows this procedure using the 2D Clumpy torus image of NGC 1068 and the PSF of TMT at 3.45 μm. These synthetic observations can be directly compared with real data. The rightmost column in Figure 4 shows the deconvolved images with and without noise. The improvement in detecting extended emission by deconvolving the noisy image is very conspicuous.

Zoom In Zoom Out Reset image size

2.2.3. Results

Figure 5 shows the synthetic observations for each telescope in the 1.2–11.6 μm wavelength range. They were obtained following the steps shown in Figure 4 through step E, i.e., they are flux calibrated and pixelated, but noise-free and not deconvolved.

Zoom In Zoom Out Reset image size

An outstanding result is that all the upcoming generations of 30 m class telescopes will be able to spatially resolve the dust emission distribution of the torus around NGC 1068 in the 1−12 μm wavelength range. Although the 6.5 m JWST telescope has the most sensitive capabilities, its resolving power, limited by the diameter of the telescope, is smaller than the torus size of NGC 1068 at any given wavelength in the IR. However, JWST can be used to isolate the torus emission from potential contributions of extended diffuse dust emission from the narrow line region, any polar/outflow components, and star-forming regions surrounding the AGN for a large sample of bolometric luminosities and AGN types. We also note that 10 m class telescopes, e.g., Keck, are not able to resolve the torus in NGC 1068 within the 1–12 μm wavelength range, in line with observations to date (e.g., Ramos Almeida et al. 2011; Alonso-Herrero et al. 2011; Asmus et al. 2015; Lopez-Rodriguez et al. 2016). We henceforth focus on the observations resolved by 30 m class telescopes.

In the 1–2.5 μm wavelength range the torus emission arises mostly from the directly irradiated clouds that make up the inner walls of the torus. In all cases, the region of emission is <1 pc (14 mas for NGC 1068), which is resolved by all 30 m class telescopes. Due to the shape of the PSF the extended torus emission is more clearly visible with GMT than with TMT or ELT. However, when deconvolution techniques are applied, TMT and ELT produce resolved images that are more comparable to the 2D original model images (see Figures 6 and 7)

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the 3–5 μm wavelength range the emission morphology becomes more extended, and importantly, the elongation is in the polar direction. If the torus axis is tilted with respect to the observer, the polar component is mainly seen on one side only due to self-obscuration by material in the equatorial region. The emission mainly originates from the inner layers of the extended dusty torus, and has sizes of ∼10 pc (142 mas), which are easily resolved by all 30 m class telescopes. For all these large telescopes the core and the extended emission are resolved. Deconvolution techniques can be applied effectively, and are crucial in providing final images comparable with the 2D models. These results demonstrate the criticality of obtaining not only an excellent but also stable PSF.

In the 8–12 μm wavelength range the torus emission shows an extended component in the polar direction, with a peak of emission in the northern inner wall of the torus. A small amount of emission is found along the equatorial plane of the torus due to the optically thick dust along such lines of sight. The extended emission is marginally resolved by all the 30 m class telescopes.

Apart from telescope size, instrument configuration and sensitivity, the physical/geometrical parameters of the AGN torus determine the resolvability of its emission. See Paper I on how these parameters affect resolvability and the existence of polar elongation of the emission across wavelengths.

Having the resolved images available at any wavelength, we can compute spectral quantities per-pixel, akin to observations with an IFU. We demonstrate this with the 10 μm silicate feature strength S₁₀, defined as

$$\begin{eqnarray}&&{S}_{10}=\mathrm{ln}\displaystyle \frac{F({\lambda }_{10})}{{F}_{\mathrm{cont}}({\lambda }_{10})},\end{eqnarray} \tag{ 1 }$$

i.e., the logarithmic ratio of the flux at the extremum of the 10 μm feature and the continuum flux at that wavelength (interpolated as a cubic spline across the wings of the silicate feature; Spoon et al. 2007).

The silicates that make up the majority of ISM dust have stretching and bending modes that produce emission at 10 and 18 μm. In the context of the AGN torus, the strength of the 10 μm feature is a powerful diagnostic because it is sensitive to several physical and geometrical conditions in the torus. One of these quantities is the inclination, which determines whether we will have a free line of sight (LOS) toward the dust grains that create silicate emission. Another quantity is the optical depth along the LOS, which, if large enough, can lead to the silicate being seen in absorption. It is possible or even likely that there is significant LOS silicate extinction by foreground matter, unrelated to the torus, in many sources (e.g., Alonso-Herrero et al. 2011; Ramos Almeida et al. 2011; Goulding et al. 2012; González-Martín et al. 2013; Prieto et al. 2014; Lopez-Rodriguez et al. 2018).

In the context of the clumpy torus, S₁₀ is also sensitive to the number of discrete clouds, and to their distribution in radial and angular directions, as clumpiness fundamentally permits edge-on configurations that afford a clear LOS toward silicate-emitting hot dust at the far inner side of the torus. S₁₀ has been measured in moderate absorption in type 2 Seyfert galaxies (Roche et al. 2006; Hao et al. 2007; Wu et al. 2009; Hönig et al. 2010; Alonso-Herrero et al. 2016). Deep absorption features are a hallmark of smooth dust distributions, and are observed in ultra-luminous infrared galaxies, which are shrouded in copious amounts of dust on large scales, but not in Seyferts (Levenson et al. 2007). Silicates have been observed in emission in type 1 Seyferts (Hao et al. 2005; Siebenmorgen et al. 2005; Sturm et al. 2005; Hönig et al. 2010; Alonso-Herrero et al. 2016), even in galaxies that were classified as type 2 for their lack of broad emission lines (Teplitz et al. 2006; Mason et al. 2009; Nikutta et al. 2009).

For the NGC 1068 model we adopted here (with parameters σ = 45°, i = 75°, Y = 18, ${{ \mathcal N }}_{0}=4$, q = 0.08, τ_V = 70) Figure 8 shows spatially resolved maps of S₁₀, i.e., the silicate feature strength is computed per pixel, using model images generated between 9 and 12 μm. Panel (b1) shows the S₁₀ map for a torus-only model (no intervening dust screen). In panels (c1) and (d1) the same torus model is viewed through a screen of ISM dust (same composition as our dusty torus clouds). We note that inhomogeneous dust screens could alter the results we derive below due to light leakage (Krügel 2009), but we do not consider them here. Instead, we consider homogeneous screens with A_V = 9 mag of extinction in panel (c1) (similar to results in Lopez-Rodriguez et al. 2018), and A_V = 15 mag in panel (d1). The mean profiles of S₁₀ in panel (e1), computed along the y direction over a ±50 mas wide band (indicated only in panel (d1) as a gray stripe) show how a foreground screen effectively reduces the strength of silicate emission, deepens any silicate absorption features, and increases the area that is observed in absorption.

Zoom In Zoom Out Reset image size

How much of the spatial S₁₀ map could be actually detected depends on the level of flux captured by each spaxel. Figure 8 therefore also shows in panels (b2), (c2), and (d2) flux-weighted versions of the S₁₀ maps. For the no-screen model only the central 1–2 pc generate moderate silicate absorption (down to S₁₀ ≈ −0.5). Polar regions above and below the torus equatorial plane produce mild silicate emission (S₁₀ ≲ 0.1). A foreground dust screen deepens the central silicate absorption, and dampens further the polar silicate emission. Our modeling shows that by A_V = 33 mag no single spaxel exhibits silicate emission anymore.

The behavior described above is in line with observations. For example, Mason et al. (2006) measure a nuclear silicate feature with S₁₀ = −0.4. Further away from the nucleus the feature becomes more shallow (see their Figure 6). One qualification is that our model FOV spans just about 04, i.e., comprises only the very center of the measurements in Mason et al (2006). Dilution of the central absorption feature with larger beams could lead to the observed reduction of the weak polar emission areas. A second difference is that in Mason et al. (2006) the feature strength observed along profiles just 04 away from the nucleus is S₁₀ ≲ −0.2, i.e., the feature is still in absorption, but weaker. This can be achieved if in addition to the nuclear absorption caused by the torus dust there is an extrinsic absorption screen along the LOS, generating an additional Δ S₁₀ ≈ −0.2. Such a screen has been invoked to explain observed data (e.g., Raban et al. 2009; Alonso-Herrero et al. 2011; López-Gonzaga et al. 2014), and can be physically associated with dust in the inner bar of NGC 1068.

For the fixed-origin model (I) all but two fitted parameters converge robustly. The results suggest a geometrically thin structure (σ = 246 for our Gaussian soft-edge torus), viewed at i = 852 inclination, and with a fitted PA = 519 of the torus axis (E of N, CCW), or PA = 1419 of the torus equatorial plane. The steepness of the radial cloud number distribution is ∝ r^−1.35. The mean number of clouds along radial equatorial rays, ${{ \mathcal N }}_{0}$, and the optical depth per cloud at visual, τ_V, do drift toward the smallest values present in our model hypercube, but do not appear to stabilize before. While we can thus not constrain them further, they do occupy the optically thin regime preferred by the G20 ring model, with the 2.2 μm emission emanating also from our near side of the dusty disk/torus. In fact, in our best-fit model with σ = 246, and taking ${{ \mathcal N }}_{0}=1$, there are on average only 0.96 clouds along our LOS to the central source (at i = 852), according to Equation (3) from Paper I. If their optical depth at visual (0.55 μm) is 10, then the overall LOS optical depth in the K band is τ_2.2 = 0.89 when assuming our standard silicate-graphite ISM dust, where τ_2.2/τ_0.55 = 0.093.

This optically thin configuration results in a K-band morphology that very closely resembles the underlying dust distribution, apparent in the two right panels in the top row of Figure 9, including the brightened edges of the central dust-free cavity. It seems challenging to reconcile such an optically thin disk with the amount of emission observed, and with the existence of coplanar and colocated maser spots (e.g., Gallimore et al. 2004) which require high densities.

The fitted zoom factor 1.037 of the images, needed in comparison to their original resolution in the model hypercube, is small. With the native Clumpy image resolution of 6 pix/R_d and a pixel scale of the GRAVITY image of 0.3634 mas/pix, a zoom factor of 1.037 means that the dust sublimation radius R_d subtends 2.26 mas on the sky. If we adopt the canonical distance to NGC 1068, D = 14.4 Mpc, then the physical size of R_d is just 0.158 pc. This is in line with some previously derived values (e.g., Burtscher et al. 2013), but on the lower end of measurements reported in the literature. Equation (2) then yields a sublimation-setting luminosity of L = 1.56 × 10⁴⁴ erg s⁻¹. This value is even smaller than the lower limit from modeling by G20. This apparently small luminosity may have several causes, as follows.

(i) During fitting we convolve a Clumpy image with the effective beam of the telescope/instrument configuration at the time of observation, which G20 estimate as a 3.1 × 1.1 mas FWHM elliptical beam, with the short axis inclined 486 E from N. This inclination is unfortunately almost identical to the PA of the 2.2 μm emitting dust structure; the convolution of the model image with the beam smears out the emission along its long axis, effectively increasing its apparent angular size. All the while the underlying model is smaller. In our tests this effect can increase the apparent angular size by up to 15%–20%.

(ii) The dust sublimation temperature also enters Equation (2), and it is often assumed at 1500 K. Increasing T by just a few hundred Kelvin, e.g., through a slight overabundance of graphites, and/or by having a slightly larger average grain size, can easily produce smaller R_d at standard luminosities. For instance, L = 0.4 × 10⁴⁵ erg s⁻¹ and T = 1800 K (the sublimation temperature of pure graphite grains, which are likely to be the ones to survive closest to the AGN) produce R_d = 0.16 pc, in line with our modeling. See, e.g., Mor et al. (2009), Heymann & Siebenmorgen (2012), Xie et al. (2015) for further discussion on the influence of dust chemistry on the dust sublimation temperature and sublimation radius in AGNs.

(iii) In our models the central illuminating source is isotropic. If it were instead emitting like an accretion disk with a $\cos \theta$ intensity falloff with polar angle θ (see, e.g., Netzer 1987), then the dust sublimation radius is a function of θ. Kawaguchi & Mori (2010) have shown that in this case the average R_d is only about a third of the isotropic case, and close to edge-on views it can be as small as 10% of the isotropic radius. The dusty disk in NGC 1068 is seen at an angle closer to edge-on, and therefore Equation (2) may not be applicable as-is, but requires an angle-dependent luminosity correction.

For model (II) all sampled parameter ranges were kept the same, and we allowed the 2.2 μm emission pattern to move up to ±10 pixels in the x and y directions within the 69 pixel wide field of view (FOV) used for fitting. A majority of the parameters have again converged robustly. ${{ \mathcal N }}_{0}$ and τ_V do not converge to a fixed value, but for later iterations of the loss function minimization their long-term mean is stable. For these two parameters we therefore adopt the mean of the final third of all iterations (over 12,000 iterations), i.e., ${{ \mathcal N }}_{0}=10.9$ clouds and τ_V = 143.5.

The freedom to change a bit the locus of the AGN allowed for optically thick models, where both the number of clouds along the LOS, and the optical depth of individual clouds are higher. Such models, when viewed at an angle somewhat above the equatorial plane (i = 693 in our case), produce thermal emission morphologies that are vertically offset from the midplane. For our best-fit model the shifts amount to 1.37 and 0.74 mas to the west and south, respectively.

The resulting morphology, after convolution (second-from-left panel in the lower row in Figure 9), is strikingly similar to that of model (I); the RMSE is ∼6% lower than in model (I). However, the rightmost panel in Figure 9 reveals just how different the underlying dust distribution is in this case, and by how much it had to be offset to cause this effect. In the unconvolved image (third-from-left panel) it is evident how the emission emanates just above the inner, thin-disk like part of the dust cloud distribution.

The radial cloud distribution is much flatter than in model (I), with q = 0.47. This is favorable for producing emission patterns preferentially elongated in polar directions, but in this case, together with the small value of σ ≈ 15°, is to first order not enough to generate the MIR elongation of 2.3:1 measured for NGC 1068 (López-Gonzaga et al. 2016); in fact, the corresponding 10 μm image of this model only produces e ≈ 0.4–0.5 at the 0.1 contour level (of the peak). We discuss the issue of polar elongation in detail in Section 3 of Paper I.

The zoom factor found with this shifting-origin model (II) is 1.344, and translates to R_d = 2.93 mas ≡ 0.2 pc, and L_bol =2.62 × 10⁴⁴ erg s⁻¹ (all for distance D = 14.4 Mpc). The same caveats (i)–(iii) for small luminosities apply as for model (I), but to a significantly lesser degree.

Overall the GRAVITY image is striking both in its unprecedented spatial resolution, but also in the fact that the circumnuclear NIR emission appears to be concentrated in only a few point-like peaks of sufficient intensity to be detected. This calls for a cautious interpretation. We emphasize that the instrumental beam smears out the observations and causes them to appear very similar for different underlying geometries. These could be either an elongated structure suggestive of a highly inclined emission ring, or a geometrically thin but optically thick flared disk where the emission arises from a narrow strip of hot cloud surface layers on the far inner side of the torus funnel. As noted in G20, extremely precise NIR astrometry of future instruments could help distinguish between the two scenarios, especially with quantitative comparison of model images such as Clumpy.

Hypercat provides three methods to estimate the PSF for a given telescope and wavelength.

Model-PSF (Airy function + Gaussian): The simplest PSF model of a circular telescope aperture is an Airy disk, i.e., the Fourier transform of the aperture. It is sufficient for seeing-limited observations. Focusing on 30 m class and space-based telescopes, our modeled observations will be nearly diffraction limited in the NIR to MIR. A better PSF model is then a combination of a central Airy disk, which contains most of the power of the PSF, and a broad halo that encompasses the remaining power of the object (e.g., Hardy & Thompson 2000). This is a semi-empirical function determined by fitting the PSF of adaptive optics observations at optical and NIR wavelengths. The free parameters of the model are the primary mirror diameter D, the Strehl ratio S, and the wavelength of the observation λ.

The FWHM of the Airy core FWHM_A = λ/D is defined as the location of the first minimum such that

$$\begin{eqnarray}&&{\mathrm{PSF}}_{{\rm{A}}}={I}_{0}^{A}{J}_{1}({\mathrm{FWHM}}_{{\rm{A}}}/2)\end{eqnarray} \tag{ A1 }$$

where I^A ₀ is the peak of the core of the PSF and J₁ is the normalized first order Bessel function used by the 2D Airy function from astropy.

The normalized intensity of the peak "p", for a Strehl ratio S_p , is defined as ${I}_{0}^{A}={e}^{-{\sigma }_{p}^{2}}$ (see Hardy & Thompson 2000), where σ_p is the mean-square wavefront error. For typical diffraction-limited observations ${\sigma }_{{p}_{0}}=0.0745$ at a Strehl of ${S}_{{p}_{0}}=0.8$. Using these values as normalization factors of the mean-square wavefront error, ${\sigma }_{{p}_{0}}/{S}_{{p}_{0}}=0.0745/0.8=0.093125$, and we can then define the peak of the core as a function of the Strehl ratio, S_p , as ${I}_{0}^{A}={e}^{-{(0.093125/{S}_{p})}^{2}}$. Note that the peak of PSF_A is always <1 because the missing flux is distributed within the halo (Gaussian profile), which is taken into account in Equation (A4).

In the general 2D case the halo is given by a double Gaussian profile

$$\begin{eqnarray}&&\begin{array}{l}{\mathrm{PSF}}_{{\rm{G}}}={I}_{0}^{G}\exp \left[-4\mathrm{ln}2\left\{{\left(\displaystyle \frac{x-{x}_{0}}{{\mathrm{FWHM}}_{{\rm{G}}}^{{\rm{x}}}}\right)}^{2}\right.\right.\\ \ \,\left.\left.+\,{\left(\displaystyle \frac{y-{y}_{0}}{{\mathrm{FWHM}}_{{\rm{G}}}^{{\rm{y}}}}\right)}^{2}\right\}\right],\end{array}\end{eqnarray} \tag{ A2 }$$

which in the circular case simplifies to

$$\begin{eqnarray}&&{\mathrm{PSF}}_{{\rm{G}}}={I}_{0}^{G}\exp \left[-\displaystyle \frac{4\mathrm{ln}2}{{\mathrm{FWHM}}_{{\rm{G}}}^{\,2}}\left\{{(x-{x}_{0})}^{2}+{(y-{y}_{0})}^{2}\right\}\right].\end{eqnarray} \tag{ A3 }$$

I^G ₀ is the peak of the Gaussian profile

$$\begin{eqnarray}&&{I}_{0}^{G}=\displaystyle \frac{1-{e}^{-{(0.093125/{S}_{p})}^{2}}}{1+{\left(\tfrac{D}{{\rho }_{0}}\right)}^{2}}.\end{eqnarray} \tag{ A4 }$$

The peak of the Gaussian profile accounts for the power missed by the central PSF, where ${I}_{0}^{A}+{I}_{0}^{G}=1$. FWHM_G is defined as

$$\begin{eqnarray}&&{\mathrm{FWHM}}_{{\rm{G}}}=1.22\,\displaystyle \frac{\lambda }{D}{\left[1+{\left(\displaystyle \frac{D}{{\rho }_{0}}\right)}^{2}\right]}^{1/2}.\end{eqnarray} \tag{ A5 }$$

The coherent length of turbulence for short exposures is

$$\begin{eqnarray}&&\rho ={\rho }_{0}\left[1+0.37{\left(\displaystyle \frac{{\rho }_{0}}{D}\right)}^{1/3}\right],\end{eqnarray} \tag{ A6 }$$

where ρ₀ is ∝ λ^6/5, with typical values of 0.1−0.2 m at 0.5 μm. We use the normalization of 0.15 m at 0.5 μm. The final model-PSF is computed as PSF = PSF_A + PSF_G with a normalized total peak flux. Figure 10 shows an example of a PSF modeled this way for a 30 m telescope at 2.2 μm, FWHM_A of 15 mas and Strehl ratios of 0.8 and 0.2. As expected, for larger Strehl ratios more power is accumulated in the core.

Zoom In Zoom Out Reset image size

Pupil-PSF (Fourier transform of the pupil): Several telescope structures such as the central obscuration of the secondary mirror, the thick spiders holding it, edges of the aperture, as well as the segmentation of the primary mirror for large telescopes, create a PSF with features far from the ideal case described above. To model the effects of these geometrical features on the PSF we can use the pupil image of the telescope. The pupil image has the shape of the primary mirror and the rest of telescope structures in front of it. The pupil images of all telescopes in Table 4 are available as FITS files in Hypercat. ¹⁵ Where they are available, Hypercat uses these pupil images by default. Figure 2 shows the pupil images and the PSFs of the telescopes from Table 2. Figures 2–7 show in their upper rows the pupil images of several telescopes. The PSF of a real telescope can then be estimated as the Fourier transform of the pupil image

$$\begin{eqnarray}&&{\rm{PSF}}({x}^{{\rm{{\prime} }}},{y}^{{\rm{{\prime} }}})=\int {\rm{Pupil}}(x,y){e}^{-j2\pi ({x}^{{\rm{{\prime} }}}x+{y}^{{\rm{{\prime} }}}y)}\,{dx}\,{dy},\end{eqnarray} \tag{ A7 }$$

where (x,y) and $(x^{\prime} ,y^{\prime} )$ are the coordinates of the object and PSF plane, respectively, and the integrals are over the image domain. The PSF plane, $(x^{\prime} ,y^{\prime} )$, is then scaled to the pixel scale in the plane of the sky at a given wavelength as λ/D, in units of arcsecond. Finally, Hypercat also generates a pixelated PSF with the same pixel scale of the source.

Image-PSF (PSF image provided externally): The previous two models might lack additional components of the PSF. Some telescopes provide forward-modeled PSF images at specific wavelengths, e.g., PSFs for HST computed with the Tiny Tim ¹⁶ tool (Krist et al. 2011), or for JWST with the WebbPSF ¹⁷ tool (Perrin et al. 2012). Telescopes also take into account observations of point sources to characterize their PSFs at several instrument configurations. Hypercat can take such PSF images, together with a pixel scale, and apply them to the 2D Clumpy model images. A drawback of this method is that a PSF image must be provided at the exact wavelength of interest (the wavelength of the model image). Interpolation of PSF images between sampled wavelengths is not advisable as the PSF behavior can be very complex, and interpolation artifacts may be of larger magnitude than the differences between neighboring wavelengths.

This method accepts the PSF taken by a given instrumental configuration can be used to produce synthetic observations. The PSF should be provided in FITS format, where the header must have the pixel scale in units of arcsecond per pixel. To compute the convolution with the model at full resolution, the PSF must also have the same FOV as the model image. Table 4 shows the available PSF modes for each telescope.

Let ${I}_{\lambda }^{\mathrm{mod}}(x^{\prime} ,y^{\prime} )$ be the 2D Clumpy model torus image at a given wavelength λ, at the full spatial sampling provided by Hypercat (as described in Section 2.5 of Paper I). ${\mathrm{PSF}}_{\lambda }(x^{\prime} ,y^{\prime} )$ be the telescope PSF (Section A.1) at the same wavelength. Then the synthetically observed image is given by the convolution of the full-resolution model and the PSF

$$\begin{eqnarray}\begin{array}{rcl}{I}_{\lambda }^{\mathrm{obs}}(x,y) & = & {I}_{\lambda }^{\mathrm{mod}}(x^{\prime} ,y^{\prime} )\,* \,{\mathrm{PSF}}_{\lambda }(x^{\prime} ,y^{\prime} )\\ & = & \displaystyle \int {I}_{\lambda }^{\mathrm{mod}}(x^{\prime} ,y^{\prime} ){\mathrm{PSF}}_{\lambda }(x-x^{\prime} ,y-y^{\prime} ){dx}^{\prime} \,{dy}^{\prime} ,\end{array}\end{eqnarray} \tag{ A8 }$$

where the PSF has been normalized to unity

$$\begin{eqnarray}&&{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\mathrm{PSF}}_{\lambda }(x,y){dx}\,{dy}=1\end{eqnarray} \tag{ A9 }$$

and has the same spatial sampling as the 2D Clumpy torus image, ${I}_{\lambda }^{\mathrm{mod}}(x^{\prime} ,y^{\prime} )$. The User Manual shows in detail how the pupil images and PSFs can be incorporated into a workflow.

While the PSF-convolved image is at the same (high) angular resolution as the model image, the detector can only generate discrete elements of resolution given by its pixel scale. We produce a pixelated image given the instrumental configuration shown in Table 2. Specifically, we downsample the PSF-convolved image to the desired detector pixel scale using a spline interpolation of order 3. The PSF-convolved image with a N_x × N_y = 241 pix × 241 pix and pixel scale Δpx_full is then pixelated to achieve the detector pixel scale ${\rm{\Delta }}{\mathrm{px}}_{\det }$, by using the pixel conversion ${\rm{\Delta }}{\mathrm{px}}_{\mathrm{obs}}={\rm{\Delta }}{\mathrm{px}}_{\det }/{\rm{\Delta }}{\mathrm{px}}_{\mathrm{full}}$. The final pixelated image has ${N}_{x}^{{\prime} }\times {N}_{y}^{{\prime} }={N}_{x}/{\rm{\Delta }}{\mathrm{px}}_{\mathrm{obs}}\times {N}_{y}/{\rm{\Delta }}{\mathrm{px}}_{\mathrm{obs}}$ pixels. Figure 4 (middle row) shows an example of this procedure using the 2D PSF-convolved Clumpy model image of NGC 1068 obtained in Appendix A.2, and pixelated to the Nyquist sampling of TMT at 3.45 μm.

In a real measurement noise must be considered due to, for instance, detector fluctuations, atmospheric emission and transmission, etc. In Hypercat we approximate the (optional) noise contribution as background-limited. The added noise pattern is a Gaussian with mean μ = 0 and standard deviation $\sigma =f\times max({\rm{image}})/{\rm{SNR}}$, that is, the noise level is specified by the desired SNR, for a signal given as fraction f of the peak pixel value. The default fraction is f = 1, i.e., the actual peak pixel value, but it can be set within $\left]0,1\right]$. This is useful for instance if the observer is interested in the SNR in the extended emission of the source, rather than at its peak. Figure 4 panel G shows an example of this procedure using the convolved and pixelated synthetic observations of NGC 1068 with TMT at 3.45 μm (see Section A.3), with an SNR of 20 at the peak pixel (and noise level fraction f = 1).

Deconvolution techniques make it possible to remove the PSF contribution from (semi-)resolved observations to obtain finer details of the astrophysical object. The Richardson–Lucy algorithm (Richardson 1972; Lucy 1974) is one of the most commonly used methods for image deconvolution in astronomy. We use it on the synthetic observations generated by Hypercat, both for the pixelated (Section A.3) and the noisy (Section A.4) images. The user can introduce the number of iterations required to produce a final deconvolved image based on their own criteria. Figures 6 and 7 show examples of this procedure using the pixelated synthetic and noisy observations of NGC 1068 with 10 iterations. For larger iterations the algorithm produces artifacts at the lower surface brightness, while for small iterations the algorithm did not have time to converge. Our criterion was to stop the iterations once the final image did recover the overall shape of the contour at the level of 0.01× the peak flux of the full-resolution image shown in Figure 6.