ALMA Polarimetry Measures Magnetically Aligned Dust Grains in the Torus of NGC 1068

The unresolved (007, 4.2 pc) core is unpolarized at ALMA Band 7 (348.5 GHz, 860 μm) as shown in Figure 1(c). However, we have found a trace of a polarized component at ∼013 (7.8 pc) along a PA of ∼105° east of north. To extract the total flux and polarization signature of the torus, we performed measurements of the total and polarized flux, and the degree and PA of polarization within a 12 × 5 pc region along a PA of 105° centered at the core of NGC 1068 (Figure 2(a)).

Zoom In Zoom Out Reset image size

This region minimizes any potential contamination arising from the polarization of the jet in the northern and southern regions of the core, which maximizes the polarization signature of the torus (Figure 2(a)). The dimensions of this region are comparable with the size of the torus found by previous ALMA observations (i.e., Gallimore et al. 2016; García-Burillo et al. 2016, 2019; Imanishi et al. 2018). Specifically, we coadded the Stokes IQU at each location along the equatorial axis within bins of 25 × 83 mas (1.5 × 4.8 pc). Then, the degree of polarization and polarized flux were estimated and debiased. Figure 2(c) shows the measurements of the total flux, polarized flux, and degree and PA of polarization along a PA of 105° (i.e., extension of the torus).

We measure that the polarized flux and the degree of polarization at the core (within the central beam) are consistent with null polarization. We find that the polarization increases as a function of the distance from the core along the equatorial plane (Figure 2(c)). Specifically, we measure a statistically significant polarized peak of 3.7% ± 0.5% in three consecutive bins from 010 (6 pc) to 018 (10.8 pc) along a PA of 105° east from the AGN. To obtain an independent polarization measurement of the east region of the torus, we estimate the polarization within a 007 × 006 (4.2 × 3.6 pc) aperture centered at 012 and PA = 105° East from the AGN. The degree and PA of polarization were measured to be 2.2% ± 0.3% and 19° ± 2°, respectively. The jet, which has a different PA of polarization, shifts the PA and reduces the degree of polarization of the large aperture measurement when compared to the single bin measurement. The measured PA of polarization of the E-vector is nearly perpendicular to the equatorial plane of the torus, i.e., ${\rm{\Delta }}\mathrm{PA}=| \left(105^\circ \pm 5\right)^\circ -\left(19^\circ \pm 2\right)^\circ | =86^\circ \pm 6^\circ$. This result shows that the measured B-vector is parallel to the equatorial plane of the torus (Figure 2(c)), i.e., the torus of NGC 1068 has a large-scale toroidal B field.

At Band 7, the nonthermal total flux emission is estimated to contribute ∼30%–65% within a 1'' aperture (García-Burillo et al. 2014). In polarized flux, Gallimore et al. (2004) found an unpolarized core at 5 GHz for the nonthermal contribution. Given the high angular resolution of our observations, we are able to disentangle the torus and jet in polarized emission, while we estimate that the total flux of the jet can contribute as much as 15% within a 01.

We identify an emission line in our observations. Figure 3 shows the emission line within the 05 × 05 (30 × 30 pc²) region identified as CN N = 3–2. We resolve the double line fine structure identified as CN N = 3–2 J = 7/2–5/2 and CN N = 3–2 J = 5/2–3/2. The integrated intensity (moment 0) of both lines is also shown in this figure, with the location of the position of the SMBH (black cross) and the polarized region (black circle) shown in Section 3.1. We fit the emission lines with two 1D Gaussian profiles. The J = 7/2–5/2 component has a peak flux of 10.62 ± 0.32 mJy, FWHM of 86 ± 5.5 km s⁻¹, and total flux of 84.13 ± 5.98 mJy centered at 1129 ± 7.1 km s⁻¹. The J = 5/2–3/2 component has a peak flux of 7.84 ± 0.52 mJy, FWHM of 76.8 ± 6.1 km s⁻¹, and total flux of 54.85 ± 5.50 mJy centered at 1334.8 ± 8.3 km s⁻¹. We estimate that ∼15% of the J = 7/2–5/2 component contributes to the J = 5/2–3/2 component and ∼7% of the J = 5/2–3/2 component contributes to the J = 7/2–5/2 component. A detailed analysis to obtain the velocity fields, i.e., moments 1 and 2, of each line is difficult, due to the blend of the fine structure. This emission line is also detected in the CND at scales of hundreds of parsec by Nakajima et al. (2015); however, these authors were not able to resolve the fine structure of this line, due to the lower resolution of their observations.

Zoom In Zoom Out Reset image size

The integrated intensity profile is concentrated along the west–east direction with an elongation of ∼03 × 015 (18 × 9 pc) at a PA of 100°–110°. The peak of the integrated intensity is slightly shifted to the west, and there is more elongation toward the eastern than the western part. This result shows similar physical conditions to the HCN J = 3–2 and HCO⁺ J = 3–2 lines (Imanishi et al. 2018). These lines show a highly inhomogeneous torus along the equatorial axis with low velocity dispersion in the eastern part of the torus, while a large velocity dispersion and strong emission are found in the western region of the torus. García-Burillo et al. (2016) found that the less dense molecular structure is lopsided. We conclude that our observed polarization region at 8 pc east from the core is cospatial with (1) a low velocity dispersion region and (2) a less dense molecular structure in the outermost eastern region of the molecular torus.

In polarized flux, we find resolved polarized emission structures along the N–S direction (Figure 1(b)) within the central 1 × 1.5 arcsec² (60 × 90 pc²). We measure polarization levels from ∼2% to ∼11% with the E-vectors of the PA of polarization mostly in the N–S direction. We overlaid (Figure 4) the 860 μm polarized flux and polarization vectors with the radio jet at 5 GHz observed with MERLIN (Gallimore et al. 2004). We find that the 860 μm polarized extended emission regions are spatially anticorrelated with the location of the radio knots, i.e., NE, C, and S2, of the radio jet at 5 GHz.

Zoom In Zoom Out Reset image size

We also overlaid the 860 μm polarized flux and polarization vectors with the polarization degree at 8.7 μm observed with CanariCam on the 10.4 m Gran Telescopio Canarias (Lopez-Rodriguez et al. 2016). The polarization degree and polarized flux at 8.7 μm show a similar morphology (Figures 2 and 3 of Lopez-Rodriguez et al. 2016). We find that the 860 μm polarized flux is located at the boundaries of the eastern region of the MIR polarized extended emission identified as the "North Knot." The North Knot is the region at ∼05 (30 pc) north from the core of NGC 1068, where the jet changes direction from N to S to ∼45° east of north due to the interaction with a giant molecular cloud. At the center of the North Knot, we find a highly polarized, ∼7%, region, which is at a similar level to the measured polarization, ∼6.5%, at 8.7 μm. However, the PA of polarization at 860 μm differs from the PA of polarization at 8.7 μm, i.e., ${\rm{\Delta }}\mathrm{PA}=| {\mathrm{PA}}_{860\mu {\rm{m}}}-{\mathrm{PA}}_{8.7\mu {\rm{m}}}| =| 175^\circ \pm 7^\circ -44^\circ \pm 3^\circ | \,=131^\circ \pm 8^\circ$. Section 5.5 discusses the origin of the submillimeter polarization of these regions.

A lot of effort has been put into successfully applied approaches (e.g., Lee & Draine 1985; Fiege & Pudritz 2000; Planck Collaboration XX 2015b; Chen et al. 2016; King et al. 2018) that have computed the synthetic observations of the polarization emission in magnetohydrodynamical (MHD) simulations. We here follow these standard approaches to estimate the expected polarization of the torus at 860 μm by dust emission of magnetically aligned dust grains. We can express the Stokes parameters in terms of the local magnetic field B = (B_x, B_y, B_z) for optically thin, τ_λ ≪ 1, dust by the form

$$\begin{eqnarray}&&I=\int n\left(1-{p}_{0}\left(\displaystyle \frac{{B}_{x}^{2}+{B}_{y}^{2}}{{B}^{2}}-\displaystyle \frac{2}{3}\right)\right){ds}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&Q={p}_{0}\int n\left(\displaystyle \frac{{B}_{y}^{2}-{B}_{x}^{2}}{{B}^{2}}\right){ds}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&U={p}_{0}\int n\left(\displaystyle \frac{2{B}_{x}{B}_{y}}{{B}^{2}}\right){ds}\end{eqnarray} \tag{ 3 }$$

(King et al. 2018), where x and y define the plane of the sky, z is parallel to our line of sight (LOS), $N=\int {nds}$ is the column density, and n is the volume density. We here define the quantity N₂ as

$$\begin{eqnarray}&&{N}_{2}=\int n\left(\displaystyle \frac{{B}_{x}^{2}+{B}_{y}^{2}}{{B}^{2}}-\displaystyle \frac{2}{3}\right){ds},\end{eqnarray} \tag{ 4 }$$

thus Stokes I can be defined as

$$\begin{eqnarray}&&I=N-{p}_{0}{N}_{2},\end{eqnarray} \tag{ 5 }$$

where the term N₂ is a corrective factor that accounts for the decrease of emission by dust grains at a given inclination with respect to the plane of the sky (Fiege & Pudritz 2000). This parameter provides a quantification of the level of dust grain alignment. For all grains aligned in the line of sight, ${N}_{2}=-\tfrac{2}{3}N$, while for all grains aligned in the plane of the sky, ${N}_{2}=\tfrac{1}{3}N$. In general, N₂ represents a small correction to the Stokes I and is of similar order to the dimensionless factor p₀. Thus, Stokes I (Equation (5)) is approximately the column density of the object. p₀ relates the grain cross sections and grain alignment properties, and it can be approximated as the maximum polarization by dust emission at the given wavelength. At submillimeter wavelengths, the polarization by dichroic emission is not greater than 10%, thus we adopt p₀ = 0.1 (Fiege & Pudritz 2000).

The polarization fraction, p, and the polarization angle in the plane of the sky, χ, are given by

$$\begin{eqnarray}&&p=\displaystyle \frac{\sqrt{{Q}^{2}+{U}^{2}}}{I}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&\chi =0.5\arctan \left(\displaystyle \frac{U}{Q}\right).\end{eqnarray} \tag{ 7 }$$

The overall magnetic field of the torus is currently unknown. We rely on magnetohydrodynamical (MHD) simulations to assume the overall structure of the magnetic field in the torus. Recent MHD simulations by Dorodnitsyn & Kallman (2017) assume an initial poloidal field in the torus. This poloidal field quickly and efficiently generates a toroidal field to maintain a long-lived torus due to differential rotation. More complex magnetic field configurations can also be assumed (e.g., Aitken et al. 2002), but here we will minimize the number of variables by assuming a toroidal structure to estimate the general behavior of the polarization as a function of the radius and optical depth. The toroidal configuration also presents a straightforward interpretation in terms of emission and absorption across the torus. For the magnetic field configuration of the torus, we assume a toroidal magnetic field in Cartesian coordinates of the form

$$\begin{eqnarray}&&{B}_{x}={B}_{0}(r)\cos (\phi )\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{B}_{y}={B}_{0}(r)\sin (\phi )\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{B}_{z}={B}_{0}(r),\end{eqnarray} \tag{ 10 }$$

where B₀(r) is the magnetic field as a function of the radial distance to the center, and ϕ is the azimuthal angle. When this magnetic field configuration is viewed at an inclination, i, and tilt angle, θ, projected on the plane of the sky, the magnetic field at the observer's frame is given by

$$\begin{eqnarray}&&{B}_{{x}^{s}}={B}_{x}\cos \theta +({B}_{y}\cos i-{B}_{z}\sin i)\sin \theta \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{B}_{{y}^{s}}=-{B}_{x}\sin \theta +({B}_{y}\cos i-{B}_{z}\sin i)\cos \theta \end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{B}_{{z}^{s}}={B}_{y}\sin i+{B}_{z}\cos i\end{eqnarray} \tag{ 13 }$$

where (x^s, y^s, z^s) are the major axis, minor axis, and the LOS, respectively. The face-on view corresponds to i = 0° and edge-on view to i = 90°. The tilt axis has a reference point, θ = 0°, along the north–south direction and positively increases east from north. We assume a magnetic field distribution with a constant strength as a function of the radial distance to the center, B₀(r) = const. (Section 5.3).

To compute the column density of the torus per LOS, we use the surface brightness and cloud distribution computed by the radiative transfer code clumpy torus (Nenkova et al. 2008). Specifically, we use the HyperCubes of AGN Tori (hypercat;¹¹ R. Nikutta et al. 2020, in preparation). hypercat uses the clumpy torus models with any combination of torus model 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 and a torus tilt angle on the plane of the sky of θ = 110° east of north based on the molecular torus morphology of our integrated flux of the CN line, the total flux continuum at Band 6 by García-Burillo et al. (2016), and the molecular torus morphology observed at HCN J = 3–2, HCO⁺=3–2 by Imanishi et al. (2018). To compute the 2D images of the dust continuum emission, we took the inferred torus parameters using the SED of NGC 1068 by Lopez-Rodriguez et al. (2018): angular width σ = 43°, radial thickness Y = 18, number of clouds along the equatorial plane N₀ = 4, index of the radial density profile q = 0.08, optical depth of each cloud at the V band τ_V = 70, and inclination angle i = 75°. These authors estimated a torus diameter of 10.2 ± 0.4 pc (∼02) with a height of ${3.5}_{-1.3}^{+1.0}$ pc (∼006), which is compatible with the sizes of the molecular torus shown in Figure 3 and by García-Burillo et al. (2016) and Imanishi et al. (2018). We use the cloud distribution, i.e., the number of clouds per LOS, multiplied by the optical depth at 860 μm, as a proxy of the column density, N. We convert the optical depth per cloud from the V band, τ_V, to 860 μm by using the synthetic extinction curve¹² for R_V = 5.5 computed by Weingartner & Draine (2001).

Using Equations (1)–(3), a toroidal magnetic field configuration inclined at i = 75° and tilted at θ = 110°, and the 2D images of the continuum dust emission using the inferred torus parameters of Lopez-Rodriguez et al. (2018), we compute the Stokes IQU at a full resolution of 1.2 mas (Figure 5). The polarized flux is estimated as the Stokes I multiplied by the degree of polarization, p. To obtain the synthetic observations at the spatial resolution and instrumental configuration of the ALMA observations, we use the Stokes IQU at full resolution and convolved with the observed beam, 77 × 58 mas (4.6 × 3.5 pc) and PA = 7513, and pixelated to 25 mas pixel⁻¹. We obtain the synthetic Stokes IQU observations using CASA v5.4,¹³ which computes the synthetic beam for the antenna configuration at the day of observations, takes into account the total integration time, and convolves it with the Stokes IQU at full resolution. Then, the final images are pixelated to the pixel scale, 25 mas pixel⁻¹, of the ALMA observations to obtain the synthetic Stokes IQU (Figure 5). Finally, we compute the degree and PA of polarization and polarized flux. The polarization vectors (E-vectors) and magnetic field pattern (B-vectors) inferred from the synthetic Stokes IQU images are shown in Figure 6. The length of E-vectors shows the degree of polarization, and their orientation shows the PA of polarization. In polarized dust emission, the measured E-vector is perpendicular to the direction of the magnetic field. We derive the B-vector map by rotating the E-vectors by 90°. Following the approach with our ALMA observations (Section 3.1), we compute the profiles of the Stokes I, degree of polarization, and polarized flux as a function of the radius at a PA of 105° using the synthetic observations (Figure 6).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We find that a toroidal magnetic field with an inclination of 75°, tilted at 110° E of N, and the cloud distribution of the torus inferred from the nuclear SED (Lopez-Rodriguez et al. 2018) are able to explain the general polarization behavior of the 860 μm polarimetric observations across the equatorial plane of the torus within the central $0\buildrel{\prime\prime}\over{.} 3\times 0\buildrel{\prime\prime}\over{.} 3$ (18 × 18 pc) of NGC 1068. By default our polarimetric model is axisymmetric, thus the measured asymmetry across the equatorial plane is not explained by our model. This asymmetry is explained due to variations of a turbulent magnetic field in Section 5.2. Figure 6 shows that the measured trend of the polarization across the equatorial plane of the torus can be explained by magnetically aligned dust grains arising from an optically thin layer at the outer edge of the torus.

In full resolution, Stokes I, a proxy of the column density N, shows that the equatorial plane and the core contain a thin layer of optically thick dust at 860 μm, where only continuum dust emission is observed in the outer parts of the torus. The synthetic observations show an unresolved torus with a slightly larger FWHM than the beam of the observations due to the contribution of emission in the equatorial direction. The east–west extended emission is clearly observed in the polarized flux. The core is unpolarized, while the outermost part of the torus along the equatorial plane shows measurable polarized flux. The total flux profile (blue open circles in Figure 6) is flatter than the observed (blue closed circles) (Figure 2(b)). We attribute this difference to the contribution of the subparsec-scale jet in the N–S direction. This component is not assumed in our model where both (1) the excess of total flux emission and (2) the drop of polarization within the beam can be due to the mixed contribution of the polarization from the torus and from the jet. For (2), Gallimore et al. (2004) measured an unpolarized core from the radio jet at 5 GHz, and we estimated a total flux contribution of <15% in the core of NGC 1068 at 860 μm. Although a multiwavelength analysis using higher angular resolution observations can potentially disentangle both the jet and torus components within the 77 × 58 mas (4.6 × 3.5 pc) region, we estimate that the measured polarized flux at ∼8 pc east from the core at PA ∼105° is dominated by dust emission of the torus.

In this section, we aim to put further constraints on the torus model geometry by exploring the synthetic polarimetric observations of NGC 1068 generated with different parameters of the clumpy torus models of Nenkova et al. (2008). This study provides an independent approach to modeling the nuclear IR SED and MIR spectroscopy using the same models (e.g., Alonso-Herrero et al. 2011; Ichikawa et al. 2015; Lopez-Rodriguez et al. 2018). Using a toroidal magnetic field in all cases, the nominal torus parameters used in Section 4.2 were explored, and the figures in the Appendix show the most characteristic synthetic observations of each parameter range.

We find that the inclination i, angular width σ, and the index radial density profile q produce significant changes in the polarization pattern of the synthetic observations. The optical depth of each cloud, τ_V, and the number of clouds, N, have no implications for the synthetic observations because the torus is mostly optically thin at 860 μm.

For the inclination i, the toroidal magnetic field was also inclined for each model as shown in Equations (12)–(14). We find that for a torus and a toroidal magnetic field at an inclination i > 65°, the observed polarization is reproduced with an uncertainty in the PA of polarization <10°. At lower inclinations, the E-vectors deviate by ∼90° from the observed PA of polarization (Figure A1 and A2). The degree of polarization for low inclinations, i.e., almost face-on views, produces a double dip in the degree of polarization due to the difference in polarization angles as a function of the azimuthal angle. The torus inclination must be close to edge on, i ∼ 90°, to reproduce the observations.

For angular widths σ < 30° (Figure A3), high (>5%) polarization is measured across the torus and the PA of polarization differs by >10° from our measurements. Given our 7σ polarization detection, the uncertainty in our measured PA of polarization is ±4°. The height of the torus is compact, and most of the dust is in a thin layer across the equatorial plane. Our observations have enough sensitivity to resolve this configuration; however, a torus with σ < 30° does not reproduce the observations. For angular widths σ > 60° (Figure A4), the height of the torus is wide, which produces low (<1%) polarization across the torus. Although this configuration produces a comparable degree of polarization, the PA of polarization at the core differs by >10° from our measurements and does not reproduce the observations. Based on the measured polarization, we find that the torus angular width must be in the range of 33°–50° to reproduce the observations.

For index radial profiles q > 1, the torus became very compact, with most of the dust concentrated within the central few parsecs (Figure A5). This configuration produces a highly polarized and unresolved core, which does not reproduce the observations. The torus index radial profiles must be q ≪ 1 to reproduce the observations.

For radial thickness Y < 10, the torus size will be <6 pc (01), which will be unresolved by our observations as well as incompatible with previous ALMA observations (i.e., García-Burillo et al. 2016; Imanishi et al. 2018). Our observations implied a torus with a radial thickness of Y ≥ 17, which implies a torus radius ≥9 pc.

To sum up, we find that to reproduce the 860 μm polarimetric observations, the torus (1) is highly inclined, i > 65°, (2) has intermediate angular widths σ = 33°–50°, (3) has a flatter radial distribution of clouds, q ≪ 1, and (4) has a radius ≥9 pc (Y ≥ 17).

We conclude that the measured 860 μm polarization arises from magnetically aligned dust grains in the torus. Our magnetic field model suggests that a constant toroidal magnetic field distribution inclined at 75° with a size up to 18 pc diameter produces the alignment of the elongated dust grains, and thus the measured polarization. This configuration produces a coherent PA of polarization (E-vectors) along the equatorial plane, which is consistent with our measured PA of polarization in the eastern part of the torus. Our polarization model suggests that the degree of polarization varies as a function of the optical depth, a known behavior of polarization emission in molecular clouds from far-IR to submillimeter (Hildebrand et al. 1999; Hildebrand & Kirby 2004). The polarization reaches a minimum at the location of the core where the AGN is highly extinguished, i.e., high optical depths. However, the degree of polarization increases at large radii and heights, i.e., low optical depths. The E-vector polarization map (Figure 6) shows this behavior. Although the expected theoretical polarization at the core is above the minimum detectable polarization by ALMA at Band 7, the effects of (1) optical depth, (2) angular resolution, and (3) contamination of the subparsec jet make the core unpolarized. Only the outermost regions of the torus are less influenced by these effects because the torus is spatially resolved and optically thin.

The magnetic field has both a constant and turbulent component (e.g., Planck Collaboration XIX 2015a). The effects of the turbulent component can be studied through the variations of the PA of polarization using a structure function (Hildebrand et al. 2009), as well as trends in the degree of polarization with the optical depth (Hildebrand et al. 1999). Our observations only detect a significant polarization in the eastern part of the torus, which is not fully explained by only using our magnetic field model.

As mentioned in the introduction and Section 3.2, the CN J = 5–2, HCN J = 3–2, and HCO⁺ J = 3–2 lines indicate that the torus is highly inhomogeneous along the equatorial axis. Specifically, low velocity dispersion is found in the eastern part of the torus, while large velocity dispersion and strong emission are found in the western region of the torus (Imanishi et al. 2018). Both observations and theoretical models (Meijerink et al. 2007) suggest that CN is found in abundance at early stages of molecular cloud evolution, or when the region is irradiated by intense UV and X-ray photons. Under these conditions, CN can be found in several scenarios: (1) in high-temperature conditions, the reaction CN + ${{\rm{H}}}_{2}\,\to$ HCN + H efficiently converts CN into HCN (Harada et al. 2010), and (2) in regions with high turbulence, it can efficiently mix ionized/atomic gas and increase the amount of CN. In general, any perturbance in the velocity fields introduces turbulence in the magnetic field. The detection of CN by our observations and the fact that the CN is associated with regions with high-turbulence, high-density clouds, and intense UV and X-ray radiation, may explain the behavior of the degree of polarization across the equatorial plane. Specifically, the measured gradient of the velocity dispersion along the equatorial plane makes the turbulent component of the magnetic field dominate in those areas with high velocity dispersion. For regions where the velocity dispersion is low, the coherent magnetic field may dominate, which can enhance the grain alignment and produce a higher degree of polarization. We conclude that the asymmetry in the degree of polarization along the equatorial axis can be explained as an increase in the level of turbulence, i.e., increase of the magnetic field turbulence at subparsecc scales, from the eastern to the western side of the torus.

The torus of NGC 1068 has been observed to have a gradient in velocity dispersion along the equatorial plane that may produce a variation of the turbulent magnetic field. To infer the relative contribution of both components, Hildebrand et al. (2009) show that the ratio of turbulent to large-scale components of the magnetic field strength, $\langle {B}_{t}^{2}{\rangle }^{1/2}/{B}_{o}$, is a function of the dispersion of polarization angles on the plane of the sky, α, given by

$$\begin{eqnarray}&&\displaystyle \frac{\langle {B}_{t}^{2}{\rangle }^{1/2}}{{B}_{0}}=\displaystyle \frac{\alpha }{\sqrt{2-{\alpha }^{2}}}.\end{eqnarray} \tag{ 14 }$$

Based on our 860 μm polarimetric observations, we estimate a dispersion of the PA of polarization of α = 5° ± 2° (0.0873 ± 0.0349 rad). We estimated the dispersion and its associated error using the variations of the measured PA of polarization per bin in the eastern side of the torus from 005 (3 pc) to 020 (12 pc) shown in Figure 6. The turbulent to constant ratio of the magnetic field strength is estimated to be $\langle {B}_{t}^{2}{\rangle }^{1/2}/{B}_{o}=0.062\pm 0.025$ within a 005–015 (3–8 pc) region of the eastern part of the torus.

The Davis–Chandrasekhar–Fermi method (DCF method hereafter; Davis 1951; Chandrasekhar & Fermi 1953) provides an empirical estimation of the plane-of-the-sky magnetic field strength. This method relates the magnetic field strength, B, with the dispersion of the polarization angles, α, of the constant component of the magnetic field, and the velocity dispersion of the gas, σ_V. This method assumes equipartition and a constant component of the magnetic field. A modified version of the DCF method takes into account the characteristic turbulent-to-ordered ratio, which allows us to estimate the strength of the large-scale magnetic field as

$$\begin{eqnarray}&&{B}_{o}\simeq \sqrt{8\pi \rho }\displaystyle \frac{{\sigma }_{v}}{\alpha }\end{eqnarray} \tag{ 15 }$$

(Hildebrand et al. 2009), which is valid when B_t ≪ B_o, and where ρ is the volume mass density in g cm⁻³, and σ_V is the velocity dispersion in cm s⁻¹.

The velocity dispersion of the molecular gas at the location of the polarized region in our observations is estimated to be ${\sigma }_{v,\mathrm{HCN}3-2}=40$ km s⁻¹ and ${\sigma }_{v,{\mathrm{HCO}}^{+}3-2}=10$ km s⁻¹ using ALMA observations (Imanishi et al. 2018). As mentioned in the previous section, the measured polarization spatially corresponds with regions of low velocity dispersion. We take a velocity dispersion of σ_V = 10 km s⁻¹. García-Burillo et al. (2016) estimated the dust mass of the torus to be (1.6 ± 0.2) × 10³ M_⊙, which assumes that other mechanisms rather than thermal emission contribute up ∼18%. Lopez-Rodriguez et al. (2018) estimated that a torus diameter of 10.2 ± 0.4 pc and a height of ${3.5}_{-1.3}^{+1.0}$ pc yield a median volume mass density of $({1.4}_{-0.62}^{+1.43})\,\times {10}^{-22}$ g cm⁻³. Based on our 860 μm polarimetric observations, we estimate a dispersion of the PA of polarization of α = 5° ± 2° (0.0873 ± 0.0349 rad). We estimate a magnetic field strength of $B={0.67}_{-0.31}^{+0.94}$ mG within the 3–8 pc region of the eastern side along the equatorial plane of the torus. Note that we have used the median volume mass density of the whole torus, while the location of the polarized region of the torus corresponds to the less dense areas of the torus. Thus, we estimate an upper limit of the large-scale magnetic field strength in the eastern size of the torus.

Our estimate of $\langle {B}_{t}^{2}{\rangle }^{1/2}/{B}_{o}$ shows that within the eastern part of the torus, at 3–8 pc from the core, there seems to be an ordered region with small turbulent components (∼0.06 times B_o). The magnetic field strength at a distance of 0.4 pc of NGC 1068 was estimated to be ${52}_{-8}^{+4}$ mG (Lopez-Rodriguez et al. 2015) with $\langle {B}_{t}^{2}{\rangle }^{1/2}/{B}_{o}$ estimated to be 0.1. To conserve the magnetic field flux, assuming perfect conductivity, the magnetic field strength varies as a function of the radial distance as B ∝ r⁻². Taken the values estimated by Lopez-Rodriguez et al. (2015), we estimate a magnetic field strength at 8 pc of ∼0.14 mG, which is in agreement within the uncertainties of our measurements. We note that a magnetic field distribution with a variable strength as a function of the radial distance to the center (B₀(r) in Section 4.1) will be required to study future high-angular-resolution polarimetric observations as well as MHD modeling of AGN.

The results presented in this section indicate that at subparsec scales, the turbulent component of the magnetic field is higher than at parsec scales. The radiative pressure dominates at subparsec scales, while at larger distances (few parsecs) the magnetic field is more coherent. We conclude that a large-scale magnetic field seems to be present in the torus from 0.4 to 8 pc, with a decrease in the turbulence component as the radius increases along the equatorial plane of the torus. At parsec scale, the magnetic fields can dominate the dynamics of the torus.

We explore other polarization mechanisms to explain the measured 860 μm polarization at the core of NGC 1068. The polarization signature of synchrotron polarization is dominated by a high degree of polarization and a PA of polarization of the E-vector perpendicular to the magnetic field direction. Although the measured PA of polarization, 19° ± 2°, is tentatively parallel to the north–south direction of the jet, the location of the measured polarization, ∼8 pc east from the core, and the measured degree of polarization, 2.2% ± 0.3%, rule out this mechanism. Self-scattering by dust has been proposed as a dominant polarization mechanism in protoplanetary disks observed with ALMA (e.g., Yang et al. 2016; Stephens et al. 2017). In Band 7 and for HL Tau (a protoplanetary disk inclined at ∼45°), a degree of polarization with a dip and PA of polarization almost perpendicular to the long axis of the disk are observed. These studies suggest that self-scattering by large dust grains of sizes >50 μm are responsible for the observed polarization. Most of the studies characterizing the dust grains sizes in AGN point to submicron to micron dust grain sizes (i.e., Li et al. 2008; Shao et al. 2017; Xie et al. 2017). Although the measured polarization and its pattern are compatible to our polarimetric observations, the requisite physical conditions of the dust grain sizes are not present, and we can rule out self-scattering mechanisms as the dominant polarization.

Although the goal of this paper is a detailed analysis of the polarization signature of the dusty torus, we find highly polarized, up to 11%, emission in the extended regions of NGC 1068. A thorough analysis of this extended polarized emission is outside the scope of this paper, but here we identify several polarization regions and provide an interpretation of the potential dominant polarization mechanisms.

We found an almost uniform PA of polarization along the N–S direction with degree of polarization in the range of 2% to 11% from 01 to 03 (6 − 18 pc) North from the core. García-Burillo et al. (2019) derived the power-law index (α, S_ν ∝ ν^α) in the range of −0.79 to −0.26, which is consistent with the α = −0.67 at 1–12 GHz by Gallimore et al. (1996). These estimations were performed within the 344.5–694 GHz continuum emission using ALMA observations over an area of 011 × 011 (6.6 pc²) centered at knot C. Using these power-law indexes, we estimate an intrinsic synchrotron polarization to be close to the theoretical 65% polarization. However, the measured polarization from our observations range between 2% and 11%. This area is cospatial with intense UV emission (e.g., Kishimoto 1999) and dust emission at MIR wavelengths (e.g., Tomono et al. 2001; Lopez-Rodriguez et al. 2016). There may be selective extinction across this area that decreases the intrinsic polarization by synchrotron emission. The inferred magnetic field morphology of the synchrotron emission is perpendicular to the measured PA of polarization. A helical magnetic field can likely explain the morphology right above the core, while a more complex structure, i.e., bow shock, at the arc-shape area of this knot is required. Another alternative can be radiative dust grain alignment (RAT; Hoang & Lazarian 2009) due to the intense UV emission from the AGN. Further modeling is required to test these hypothesis.

We measure a highly polarized, ∼7%, region spatially coincident with the North Knot detected at 8.7 μm. The North Knot polarization decreases from ∼7% to ∼4% in the 8.7–11.6 μm wavelength range (Lopez-Rodriguez et al. 2016), which was interpreted as polarization arising from dust and gas emission heated by the post-shock gas after the jet–molecular cloud interaction. The dust emission is then further extinguished by aligned dust grains in the inner bar at ∼45° E of N of the host galaxy. At submillimeter wavelengths, the dust in the inner bar is optically thin. Therefore, the measured 860 μm polarization of the North Knot is the intrinsic polarization by thermal emission from magnetically aligned dust grains.

We found several highly polarized, up to 11%, regions at ∼07 (∼42 pc) NE from the core. The NE knot has a spectral index of −1.04 at 1–12 GHz (Gallimore et al. 1996), and −1.6 at 694–344.5 GHz (García-Burillo et al. 2019), which is consistent with synchrotron emission and comparable to the spectral index at the large-scale of the jet (Wilson & Ulvestad 1987). The measured 860 μm polarization is in the range of 5%–11% with a PA of polarization almost N–S, with a change in the degree of polarization along the western part of the NE knot (where the peak polarization is located). Using a spectral index of −1.6, we estimated an optically thin synchrotron emission polarization of ∼60%, and ∼9% for optically thick synchrotron emission. The difference in polarization under the optically thin regime and the comparable polarization of the optically thick synchrotron emission indicates a strong depolarization due to Faraday rotation on the NE knot. The difference in the degree and PA of polarization in the western part of the knot may indicate the presence of a shock by the jet with the ISM and/or a disruption of the knot in the radio jet.