The Emission and Distribution of Dust of the Torus of NGC 1068

NGC 1068 was observed as part of the Guaranteed Time Observations (GTO; PI: T. Herter) on 2016 September 17 using FORCAST (Herter et al. 2012) on the 2.5 m SOFIA telescope. We made observations with the dual-channel mode at the 19.7, 31.5, and 37.1 μm using the two-position chop-nod (C2N) method with symmetric nod-match-chop to remove time-variable sky background and telescope thermal emission and to reduce the effect of 1/f noise from the array. In all observations, we used an instrumental position angle, i.e., long axis of the detector with respect to the north on the sky, of 305°, a chop-throw of 1' with a 30° E of N chop-angle. The on-source times were 427s, 471s, and 343s at 19.7 μm, 31.5 μm, and 37.1 μm, respectively.

SOFIA provided reduced data using the forcast redux pipeline v1.1.3 following the method described by Herter et al. (2013) to correct for bad pixels, "droop" effect, nonlinearity, and cross-talk. The point-spread functions (PSFs) of the observations were estimated using observations of Ceres taken immediately before NGC 1068 observations with the same instrumental configuration and bands. We estimated a full width at half maximum (FWHM) of Ceres of 24, 28, and 29 at 19.7 μm, 31.5 μm, and 37.1 μm, respectively. NGC 1068 was flux-calibrated using the set of standard stars of the observing run, which provides flux uncertainties of 5.0%, 5.2%, and 7.7% at 19.7 μm, 31.5 μm, and 37.1 μm, respectively.

NGC 1068 was observed as part of the GTO (PI: D. Dowell) on 2017 May 06 using HAWC+ (Vaillancourt et al. 2007; D. A. Harper et al. 2018, in preparation) on the 2.5 m SOFIA telescope. We made observations using the Lissajous pattern in the total intensity mode at 53 μm (${\lambda }_{c}=53$ μm, ${\rm{\Delta }}\lambda /{\lambda }_{c}=0.17$ bandwidth). In this new SOFIA observing mode, the telescope is driven to follow a parametric curve at a nonrepeating period whose shape is characterized by the relative phases and frequency of the motion. Figure 1 shows the Lissajous pattern of a single observation at 53 μm of NGC 1068 with a scan rate of 100'' s⁻¹ and a 60'' scan amplitude. We performed a total of five Lissajous scans with relative phases of 5° and 27° with a total on-source time of 455 s.

Zoom In Zoom Out Reset image size

We reduced the data using the Comprehensive Reduction Utility for SHARC II v2.34-3beta (crush, Kovács 2006, 2008)¹⁷ optimized for HAWC+ and the hawc_dpr pipeline v1.1.1. crush estimates and removes the correlated atmospheric and instrumental signals, solves for the relative detector gains, and determines the noise weighting of the time streams in an iterated pipeline scheme. The PSF was estimated using Uranus observations on 2017 May 07 with an FWHM of 49, consistent with diffraction-limited observations at 53 μm. Flux calibrators were not observed during the same flight, thus we cross-calibrated our observations using flux calibrators, i.e., Uranus, from other flights. Although we find a flux calibration accuracy of ∼8% for observations taken within the same flight, the cross-calibration between flights can only ensure a flux accuracy of ∼20%.

We aim to obtain the emission from the unresolved core of NGC 1068 at all the observed wavelengths. The SOFIA observations of NGC 1068 show a resolved core (Figure 2) that is thought to arise from an unresolved component and an extended component. To obtain the fractional contribution to the total emission from both unresolved and extended components, we made two different photometric measurements (Table 1). First, the flux in a circular aperture of 10'' (700 pc) diameter was measured, which ensures to enclose the whole flux of an unresolved source at the given wavelength and minimizes the contribution from the diffuse extended emission. Second, the central 20'' × 20'' (1.4 × 1.4 kpc²) emission was fitted with a composite model using the corresponding PSF to each observation and a 2D Gaussian profile. We refer to these methods as "aperture" and "PSF-scaling" photometry, respectively, in the remainder of the paper. The aperture photometry represents the total flux from the observed galaxy at a given wavelength, F_T. In the PSF-scaling method, the total flux from the scaled-PSF, F_PSF, represents the maximum likely contribution from an unresolved nuclear component at the given angular resolution of the observations, while the total flux from the 2D Gaussian profile, F_ext, provides the minimum contribution of the extended component surrounding the central source. We estimated the total flux of the model, ${F}_{T}^{M}$, as the sum of both the PSF and the 2D Gaussian profile. The PSF-scaling method has five free parameters, the amplitudes of both PSF and 2D Gaussian profile, the FWHM of the long, b, and short, a, axes and the position angle (P.A.) of the 2D Gaussian profile. The fitting routine minimizes the residuals (galaxy minus model: scaled-PSF + 2D Gaussian) to a level <5% of the total flux, F_T, within the central 20'' diameter. We also considered a 2D Sérsic profile to fit the extended component. We obtained index profiles ∼0.5 and size parameters of the Sérsic profiles similar to the FWHM of the observations, which is close to the special case of the Sérsic profile tending to a Gaussian profile. Due to this behavior and that the Sérsic profile increases the number of free parameters, we decided to use 2D Gaussian profiles.

Zoom In Zoom Out Reset image size

Table 1.  Measured and Modeled Nuclear Photometry

Photometric Measurementsa							Spectral Decompositionb
${\lambda }_{c}$	FT	F${}_{T}^{M}$	FPSF	Fext	PSF	Extended	Star Formation	Dust at 200 K	Torus
(μm)	(Jy)	(Jy)	(Jy)	(Jy)	%	%	%	%	%
19.7	61.9 ± 3.1	60.9 ± 3.8	22.0 ± 1.4	38.9 ± 2.4	36 ± 4	64 ± 7	10	60	30
31.5	59.4 ± 3.1	58.4 ± 3.7	28.8 ± 1.8	29.6 ± 1.9	49 ± 5	51 ± 6	23	36	41
37.1	59.6 ± 4.6	58.8 ± 5.0	29.7 ± 2.5	29.1 ± 2.5	51 ± 8	49 ± 8	33	26	41
53.0	71.6 ± 14.3	71.5 ± 14.5	23.8 ± 4.8	47.7 ± 9.7	33 ± 15	67 ± 25	64	16	20

Notes. Fractional contribution of emissive components in a 10'' aperture.

^aMeasured and modeled photometry as described in Section 2.3. ^bFractional contribution of the several components used in the spectral decomposition described in Section 4. We estimate a 5% uncertainty for the fractional contribution of each component.

Download table as:  ASCIITypeset image

The uncertainty in the photometry was estimated in the following manner. The aperture photometry uses the flux uncertainties estimated from the flux calibration described in Sections 2.1 and 2.2. For the PSF-scaling photometry, an estimate of the error induced by a variable PSF was obtained by cross-calibrating the standard stars observed on the same or several nights. This error was found to be ∼5%. Another estimate of the error induced by the fitting procedure was estimated to be ∼3%. The total uncertainty for the PSF-scaling photometry was calculated by adding in quadrature these individual contributions.

Figure 2 shows the NGC 1068 observations, scaled PSF, Model (PSF+2D Gaussian), and residuals (NGC 1068–Model) of the central 50'' × 50'' (3.5 × 3.5 kpc²) observations at 19.7, 31.5, 37.1, and 53 μm. Table 1 shows the measured and modeled nuclear photometry for each photometric method and model component in the central 10'' aperture. The fractional contribution of the PSF and 2D Gaussian profile in the central 10'' diameter is also shown. Our flux density at 19.7 μm estimated by using PSF-scaling of 22.0 ± 1.4 Jy is in excellent agreement with the flux density of 20.2 ± 3.4 Jy in a 04 aperture by Tomono et al. (2001). We took PACS/Herschel spectroscopic data of NGC 1068 from the Herschel Archive and we obtained a nuclear flux density of ∼70 Jy at 60 μm per spaxel, where a spaxel is $9\buildrel{\prime\prime}\over{.} 4\,\times \,9\buildrel{\prime\prime}\over{.} 4$ (658 × 658 pc²). Despite the difference in wavelength, this result is in good agreement with our total flux of 72 ± 14 Jy at 53 μm using HAWC+. At all wavelengths, the total flux model, ${F}_{T}^{M}$, is <2% of the total flux of the observations, F_T, and the P.A. of the extended emission is $44\buildrel{\circ}\over{.} 7\pm 1\buildrel{\circ}\over{.} 3$ with a decrease of the ratio of the short and long axis, a/b, from 0.85 at 19.7 μm to 0.61 at 53 μm. Our extended diffuse emission is spatially coincident with the large scale, 32'' (1.92 kpc) inner bar at a P.A. of $48^\circ \pm 3^\circ$, the so-called NIR bar (Scoville et al. 1988; Schinnerer et al. 2000; Emsellem et al. 2006). Despite any contribution of diffuse extended emission within the PSF of SOFIA in the 20–53 μm wavelength range (Section 4), we found a turnover of the unresolved emission, F_PSF in Table 1, within the 31.5–53 μm wavelength range. Specifically, the PSF fractional contribution to the total flux decreases from ∼50% in the 30–40 μm wavelength range to <40% at shorter and longer wavelengths. This result is in agreement with the observational constraint by Fuller et al. (2016), who suggested that the turnover of the torus emission does not occur until wavelengths >31.5 μm for a sample of nearby AGNs.

We here describe the details of the fitting to the nuclear SED of NGC 1068. The nuclear SED is composed of our PSF-scaling photometry (F_PSF) in conjunction with the 04 aperture photometry from 2 to 20 μm photometry by Tomono et al. (2001), the 8–13 μm nuclear spectrum in a 04 aperture using Michelle on the 8.1 m Gemini-north Telescope by Mason et al. (2006), and the 432 μm ALMA observation by García-Burillo et al. (2016). We fitted the nuclear SED using the clumpy torus models of Nenkova et al. (2008) and the BayesClumpy approach developed by Asensio Ramos & Ramos Almeida (2009). This approach has been successfully applied to this and other Seyfert galaxies (e.g., Alonso-Herrero et al. 2011; Ichikawa et al. 2015). The free parameters of the model were set with a flat prior distribution, with the exception of the foreground visual extinction to the core, A_V, was set to be in the 0–10 mag range. The best inferred model with the 1σ uncertainty region is shown in Figure 3, and the posterior distributions of each model parameter are shown in the Appendix. Table 2 shows the output values of each torus model parameter.

Zoom In Zoom Out Reset image size

Table 2.  Clumpy and Smooth Torus Model Parameters

clumpy Torus			Smooth Torus
Parameter	Symbol	Value	Parameter	Symbol	Value
Angular width	σ	${43}_{-15}^{+12}$°	Opening angle	${\theta }_{{OA}}$	${37}_{-8}^{+23}$°
Radial thickness	Y	${18}_{-1}^{+1}$	Radial thickness	Ys	${20}_{-4}^{+4}$
Number clouds along the equatorial plane	N0	${4}_{-1}^{+2}$	⋯	⋯	⋯
Index of the radial density profile	q	${0.08}_{-0.06}^{+0.19}$	Index of the radial density profile	qs	1 (fixed)
Optical depth of each cloud	${\tau }_{v}$	${70}_{-14}^{+6}$	Optical depth of the torus, LOS	${\tau }_{v,s}$	${250}_{-10}^{+20}$
Viewing angle	i	${75}_{-4}^{+8}$°	Viewing angle	is	${79}_{-10}^{+7}$°
Inner radius	rin	${0.28}_{-0.01}^{+0.01}$ pc		${r}_{\mathrm{in},s}$	${0.41}_{-0.02}^{+0.05}$ pc
Outer radius	rout	${5.1}_{-0.4}^{+0.4}$ pc		${r}_{\mathrm{out},s}$	${8.5}_{-0.7}^{+7.9}$ pc
Height	H	${3.5}_{-1.3}^{+1.0}$ pc		Hs	${4.2}_{-0.2}^{+0.5}$ pc
Bolometric luminosity (erg s−1)	Lbol	${5.02}_{-0.19}^{+0.15}\times {10}^{44}$		${L}_{\mathrm{bol},s}$	${1.11}_{-1.23}^{+0.28}\times {10}^{44}$

Download table as:  ASCIITypeset image

In general, the full family of clumpy torus model solutions when using the SOFIA observations from 20 to 53 μm provides a tighter 1σ dispersion than previous studies (i.e., Alonso-Herrero et al. 2011, Figure 5; García-Burillo et al. 2016, Figure 4, our Figure 4). This result is due to the better sampling of the turnover of the torus emission in the 30–40 μm as in comparison with previous studies. The median value for the foreground visual extinction was found to be ${A}_{V}={9}_{-1}^{+1}$ mag. We notice that if the extinction to the core was set to be negligible, ${A}_{V}\lt 5$ mag, then the fitting tends to obtain viewing angles, i, of the torus in the $30^\circ \mbox{--}40^\circ$ range, and does not fit the SED in the 1–5 μm wavelength range. Packham et al. (1997) found that a visual extinction of 36 mag to the core of NGC 1068 can explain the absorptive polarization at 2.0 μm, compatible with the expected null polarization observed at 10 μm by the emissive polarization of the torus found by Lopez-Rodriguez et al. (2016). Both works also found a visual extinction by the central dust lane to be ∼8 mag. Thus, we expect that the torus emission is extinguished by a column of dust into our LOS with a visual extinction in the range of 8–36 mag, our computed visual extinction of ${A}_{V}={9}_{-1}^{+1}$ mag is in agreement within that range. We obtained a viewing angle, $i={75}_{-4}^{{+8}\circ}$ compatible with the H₂0 maser observations in the central parsec of NGC 1068, which suggests a torus with an almost edge-on view, ∼90°.

Zoom In Zoom Out Reset image size

Based on the best inferred clumpy torus model, we can estimate torus morphological parameters as the outer radius, ${r}_{\mathrm{out}}={r}_{\mathrm{in}}Y$ pc, where r_in is the inner radius of the torus defined by the distance of the sublimation temperature of dust grains, T, as a function of the bolometric luminosity, L_bol, as ${r}_{\mathrm{in}}\,=0.4{({L}_{\mathrm{bol}}/{10}^{45}{\text{erg s}}^{-1})}^{0.5}(T/1500\,{\rm{K}})$ pc (Barvainis 1987), and the torus scale height as $H={r}_{\mathrm{out}}\sin \sigma$ pc. The estimated bolometric luminosity from our clumpy torus model, ${L}_{\mathrm{bol}}\,={5.02}_{-0.15}^{+0.19}\times {10}^{44}$ erg s⁻¹, yields an inner torus radius of ${r}_{\mathrm{in}}\,={0.28}_{-0.01}^{+0.01}$ pc for dust grains at a temperature of 1500 K. Using $Y={18}_{-1}^{+1}$ and $\sigma ={43}_{-15}^{+12}$°, the torus radius and scale height are estimated to be ${r}_{\mathrm{out}}={5.1}_{-0.4}^{+0.4}$ pc and $H={3.5}_{-1.3}^{+1.0}$ pc, respectively. Table 2 summarizes these results. These results are in agreement with the resolved 7–10 pc torus extension (García-Burillo et al. 2016), 12 × 7 pc (Gallimore et al. 2016), and 12 × 5 pc (Imanishi et al. 2018) torus diameter using ALMA observations. We can estimate the visual extinction into our LOS, ${A}_{{\rm{v}}}^{\mathrm{LOS}}$, as ${A}_{{\rm{v}}}^{\mathrm{LOS}}=1.086{N}_{0}{\tau }_{{\rm{v}}}\exp (-{(i-90)}^{2}/{\sigma }^{2})$ mag. From our best clumpy torus model, we estimate ${A}_{{\rm{v}}}^{\mathrm{LOS}}\,={248}_{-142}^{+201}$ mag. Using the standard Galactic ratio ${A}_{{\rm{v}}}/{N}_{{\rm{H}}}=5.23\,\times {10}^{-22}$ mag cm² (Bohlin et al. 1978), we estimate a column density of ${N}_{{\rm{H}}}={4.7}_{-2.7}^{+3.9}\times {10}^{23}$ cm⁻².

We have also used smooth torus models (Efstathiou & Rowan-Robinson 1995) to fit the nuclear SED of NGC 1068. The best-fit model and the output parameters are shown in Figure 3 and Table 2, respectively. In general, the smooth torus models reproduce well the nuclear SED of NGC 1068, except for the FIR (20–60 μm) wavelength range. In this spectral range, the smooth torus models underestimate the measured nuclear fluxes. We note that if we force these models to fit the FIR range, then the smooth torus models overpredict the submillimeter fluxes by a factor of 10 or more, and the torus size increases to a few tens of parsecs. Thus, we use the MIR spectroscopic observations, i.e., the 10 μm silicate feature, and the ALMA observations to find the best fit of the smooth torus model.

We find a smooth torus with similar physical characteristics as the clumpy torus (Table 2), except for the outer radius, which is larger in the case of the smooth torus. This difference is mainly due to the sublimation temperature used by both models, the smooth torus models use a maximum temperature of dust grains of 1000 K, in comparison with the 1500 K used by the clumpy torus models. Although we exclusively used the smooth torus models of Efstathiou & Rowan-Robinson (1995), we speculate that other smooth torus models (i.e., Schartmann et al. 2005; Fritz et al. 2006) may similarly fail to account for the complete SED of the torus. This may be due to the fact that smooth models generally have much less flexibility in the specification of the distribution of the dust in the torus and especially its outer part to which the far-infrared and submillimeter observations are more sensitive. This certainly merits further study and we plan to pursue this in future work.

We here investigate the emission and distribution of dust in the torus of NGC 1068 using several SED coverages. Direct comparison between previous studies are not straightforward due to the development of the clumpy torus models through the past several years.¹⁸ Our nuclear SED is constructed using the mentioned previous studies in Section 3.1. Thus, to avoid any potential misinterpretation of the physics of the torus that a slightly different version of the models could introduce, we here reanalyze the nuclear SED as a function of the SED coverage using the most updated version of the Clumpy torus models. Table 3 shows the output parameters of the best inferred Clumpy torus model for an SED using 1–20 μm imaging and spectroscopic observations (labeled as NIR+MIR), and for an SED using 1–20 μm imaging and spectroscopic observations and ALMA observations (labeled as NIR+MIR+ALMA). The posterior distributions for each parameter with their median value and 1σ error are shown in the Appendix. For each SED coverage, Figure 4 shows the best inferred clumpy torus model and their 1σ uncertainty.

Table 3.  clumpy Torus Model Parameters As a Function of SED Coverage

σ (°)	Y	N0	q	${\tau }_{v}$	i (°)	rout (pc)	SED
20${}_{-3}^{+5}$	${13}_{-3}^{+4}$	${11}_{-3}^{+2}$	${0.22}_{-0.13}^{+0.20}$	${28}_{-6}^{+10}$	${75}_{-6}^{+4}$	${3.5}_{-0.9}^{+1.3}$	1–20 μm SED+MIR Spectroscopy
31${}_{-8}^{+20}$	${19}_{-1}^{+1}$	${5}_{-2}^{+3}$	${0.06}_{-0.04}^{+0.08}$	${59}_{-13}^{+16}$	${71}_{-3}^{+5}$	${5.5}_{-0.4}^{+0.4}$	1–20 μm SED+MIR Spectroscopy+ALMA

Download table as:  ASCIITypeset image

When clumpy torus model fitting is used with data only in the 1–20 μm wavelength range, (1) the torus is smaller and more compact (large q values) than the current resolved torus of ∼10 pc diameter of NGC 1068 by the ALMA observations, and (2) the turnover of the torus emission peaks at shorter wavelengths than when the SED coverage includes observations at longer wavelengths. This result implies that 1–20 μm observations are not able to probe the full extent of the torus. Despite the angular resolution, $2\buildrel{\prime\prime}\over{.} 4\mbox{--}4\buildrel{\prime\prime}\over{.} 9$, Figure 4 shows that the turnover of the torus emission occurs in the range of 30–40 μm, which corresponds to a characteristic temperature of 70–100 K. This result indicates that (1) the amount of cold dust and/or (2) the radiation from indirectly radiated clouds is substantial to shift the peak emission of the torus toward longer wavelengths.

We use the radiative transfer code clumpy torus (Nenkova et al. 2008) to compute the surface brightness and cloud distributions of the dusty torus as a function of wavelength for each set of parameters shown in Table 3. Specifically, we use the HyperCubes of AGN Tori (HyperCAT¹⁹ ; R. Nikutta et al. 2018, in preparation). HyperCAT uses the clumpy torus models with any combination of parameters to generate physically scaled and flux-calibrated 2D images of the dust emission and distribution for a given AGN. We use a distance of 14.4 Mpc, a torus orientation on the plane of the sky of ∼138° east of north based on the IR polarimetric signature of the nucleus (e.g., Packham et al. 1997; Simpson et al. 2002; Gratadour et al. 2015; Lopez-Rodriguez et al. 2015), and the bolometric luminosities and model parameters by each inferred model from Table 3. Figure 5 shows the dust emission distribution from 2.2 to 432 μm for the 1–20 μm SED+MIR Spectroscopy (labeled as NIR+MIR), 1–20 μm SED+MIR Spectroscopy+ALMA (labeled as NIR+MIR+ALMA), and this work, whose torus parameters are summarized in Table 2. The last column shows the cloud distribution for each clumpy torus.

Zoom In Zoom Out Reset image size

For all cases, the clouds are distributed in the equatorial plane with major differences in their torus sizes, angular widths, and radial density profile (Figure 4, last column). These differences affect the morphology of the dust emission as a function of wavelength. We find that the 2.2 μm dust emission is concentrated on the inner edge of the torus where dust is directly radiated by the central engine, while the 8–12 μm dust emission is along the polar direction as the high opacity in the equatorial direction is absorbing most of the radiation from the central engine. At longer wavelengths, >30 μm, the dust emission is along the equatorial plane, where the 432 μm truly describes the bulk of dust distribution in the torus.

We point out the tight morphological and size similarities between our 2D clumpy torus image at 432 μm, using the well sampled SED from 1 to 432 μm, with the observed torus emission by the ALMA observations (García-Burillo et al. 2016; Imanishi et al. 2018). However, when ALMA observations are compared with the 2D images produced by the 1–20 μm SED, the inferred torus model is smaller, more compact, and thinner, which supports the discussion above regarding the 1–20 μm observations, which underestimate the true size of the torus.

Using our 2D clumpy torus image at 12.0 μm, we find that the dust emission along the polar direction is ∼5 times that of the central dust emission. We estimate the dust emission along the polar direction using an ellipse of semimajor axis of 139 mas (9.7 pc) and an eccentricity of 0.91, while the central dust emission was estimated in a circular aperture of 50 mas (3.5 pc) diameter. At this wavelength, the north emission is more prominent than that coming from the south, due to the inferred torus inclination of 75°. We estimate that the fractional contribution along the polar direction increases at larger inclinations, although further extinction by the host galaxy and/or dust in the NLR will play an important role in the observed emission as a function of the location in the galaxy. López-Gonzaga et al. (2014), using IR interferometric observations with MIDI/VLTI, found that dust emission in the polar regions at scales of 5–10 pc contributes four times more at 12 μm than dust located in the torus. They also found a dependency on flux density as a function of the location along the north–south direction attributed to extension by the host galaxy. Our results using a solely 2D clumpy torus model and those by IR interferometric observations are of comparable order of magnitude, which tentatively indicates that the IR interferometric observations of NGC 1068 may be observing the dust emission from the optically thin dust of the torus.