SOFIA/HAWC+ Traces the Magnetic Fields in NGC 1068

Over the past few decades, astronomers have detected the presence of magnetic fields in galaxies at all spatial scales. These major studies have been performed using optical and radio observations (Kronberg 1994; Zweibel & Heiles 1997; Beck & Gaensler 2004; Beck 2015, for reviews). Radio observations measure the polarization of synchrotron radiation from relativistic electrons; this radiation is sensitive to the cosmic-ray electron population, which does not closely trace the gas mass. Studies of interstellar polarization at optical wavelengths can reveal the magnetic field geometry as a result of magnetically aligned dust grains by radiative alignment torques (Andersson et al. 2015). However, the optical polarization measurements suffer from contamination by highly polarized scattered light. Linear polarization at 6.3 and 21.2 cm using the Effelsberg 100 m telescope and optical polarization (Scarrott et al. 1987) have been compared in the face-on spiral galaxy M51 (Beck et al. 1987; Fletcher et al. 2011). The optical and radio linear polarization shows a similar magnetic field morphology within the eastern and southern quadrants, but with significant variations of the magnetic field direction, up to 60°, in the western quadrant. A study on the face-on spiral galaxy NGC 6946 also shows similarities between the optical and radio polarization (Fendt et al. 1998). Further observations of M51 using near-infrared (NIR) polarization only registered upper-limit polarization across the galaxy, which ruled out the dichroic absorption as the main polarization mechanism (Pavel & Clemens 2012). The scattering cross section of typical interstellar dust declines much faster (∼λ⁻⁴) between 0.55 and 1.65 μm than its absorption (∼λ⁻¹, Jones & Whittet 2015). It is likely that the optical polarization measured in the previous works (i.e., Scarrott et al. 1991) was due to scattering, rather than extinction by dust grains aligned with the interstellar magnetic field. For M51, the expected dichroic absorptive polarization at H-band, based on the measured optical polarization, is ∼0.4%, but an upper limit of 0.05% was measured at H-band (Pavel & Clemens 2012). The dichroic absorptive polarization has shown a higher polarization at H-band than the dust grains measured to be due to magnetically aligned. Observations of the magnetic field geometry more sensitive to the denser gas and dust are needed.

Our understanding of how spiral arms form and their role in galaxy evolution is still incomplete. The most widely accepted theoretical model for spiral arms is the density wave theory (Lindblad 1960; Lin & Shu 1964; Shu 2016). This theory posits that spiral structure can be described as a superposition of waves of enhanced stellar density with constant pitch angles (the angle between the tangent of the wave and circles around the galactic center) and constant pattern speed (Athanassoula 1984). This theory predicts that stars form in the arms as gas moves into the wave and is compressed by its gravitational potential. Under this scheme, the spatial displacement of the spiral arms should be different for different tracers of star formation (e.g., molecular clouds, H ii regions, and newly formed stars) because they appear at different phases of the wave. The spiral arms at optical/NIR wavelengths is expected to trace already born stars, while far-IR (FIR) wavelengths will trace ongoing star formation.

NGC 1068 is the nearest (D_L = 13.5 Mpc, 1'' = 65 pc) grand-design spiral galaxy with both a bright active galactic nucleus (AGN) and a luminous circumnuclear starburst. This galaxy is classified as a Seyfert 2, where the active nucleus is obscured by a dusty structure, and the host is an Sb type. Associated with the AGN is a narrow-line region, i.e., ionization cones, of ∼20'' (∼1.3 kpc) in diameter at a position angle (PA) of ∼40° east of north (i.e., Balick & Heckman 1985) with the northern region protruding toward us out of the plane of the galaxy. Schinnerer et al. (2000) have presented CO observations of the central region of NGC 1068 and discussed the relationship between the stellar and gas dynamics in the galaxy. At radii >15'' (>0.975 kpc), the kinematic axis of the galaxy is approximately east–west, and at very large radii the eccentricities of the galaxy's faintest contours are consistent with the 40° inclination suggested by H i data. However, there is a bright oval structure ∼180'' (∼11.7 kpc) diameter with long axis perpendicular to the kinematic axis that probably corresponds to a large-scale stellar bar (Figure 6 of Schinnerer et al. 2000). The inner Lindblad resonance (ILR) of this bar has a deprojected radius of ∼18'' (∼1.17 kpc), near the location of a compact spiral or ring-like structure in ¹²CO, consistent with theoretical predictions that gas can be transported inward along bars and accumulate at the ILR. NIR observations in the K-band have detected the presence of a 30'' (1.95 kpc) diameter inner-bar (Scoville et al. 1988; Thronson et al. 1989). Telesco & Decher (1988) have suggested that for the case of the NGC 1068 outer-bar the ILR is actually split into two resonances ("inner–inner" or iILR and "outer–inner" or oILR) and that the CO spiral pattern is produced by a spiral density wave between the iILR and oILR that is driven by the inner-bar. We will refer to this annulus as the "starburst ring." In addition to the CO in the starburst ring, there are CO clumps extending down into the northeastern branch of the inner-bar (but not in the southwestern branch) and a compact CO ring at radii ∼1'' (∼65 pc) that may be gas accumulating at the ILR of the inner-bar.

There is ample evidence for an exceptionally high rate of star formation in the central region of NGC 1068, possibly the result of a recent minor merger (Tanaka et al. 2017). The optical surface brightness of the galactic disk at radii ≤25'' (≤1.63 kpc) (henceforth the "starburst disk") is among the brightest of any galaxy in the local universe (Keel & Weedman 1978; Weedman 1985; Ichikawa et al. 1987). Observations in the FIR have shown that the total luminosity of the starburst disk is >10¹¹L_⊙ (Telesco & Harper 1980; Telesco et al. 1984). Radio observations at 1.465 GHz (Wilson & Ulvestad 1982) show that the same region is extremely rich in the type of nonthermal synchrotron emission produced in supernova explosions of massive young stars. On the basis of 10.8 μm and NIR data, Telesco & Decher (1988) and Thronson et al. (1989), respectively, have argued that almost all the actual star formation takes place within the starburst ring and that the most intense activity occurs near the outer ends of the inner-bar.

These properties have made NGC 1068 a suitable object for spatially resolved polarimetry. Using optical polarimetry of NGC 1068, Scarrott et al. (1991) found a spiral pattern that was interpreted as delineating the magnetic field geometry in the spiral arms of the galaxy. Given the effects of scattering on optical polarimetry measurements and that the host galaxy of NGC 1068 has not been observed using radio polarimetry, this association has to be questioned.

We have performed FIR polarimetric observations to image the central 2' × 2' (7.8 × 7.8 kpc²) of NGC 1068 at wavelengths of 53 and 89 μm. We discuss the 53 μm flux from the AGN in a previous paper (Lopez-Rodriguez et al. 2018b). In this paper, we present our polarimetric results, including an 89 μm image that for the first time reveals galactic spiral structure by means of thermal emission from magnetically aligned dust grains. At FIR wavelengths, scattering and Faraday rotation are not a factor at these scales. The dominant emission is from warm dust, which more closely samples the total gas column density than relativistic electrons producing the synchrotron emission. The paper is organized as follows: Section 2 describes the observations, data reduction, and observational results. In Section 3, we fit our polarimetry data to a spiral galactic magnetic field model, which is then analyzed and discussed in Section 4. In Section 5 we present our conclusions.

NGC 1068 was observed at 53 and 89 μm as part of the guaranteed time observations (GTO, ID: 70_0609) using the High-resolution Airborne Wideband Camera-plus (HAWC+) (Vaillancourt et al. 2007; Dowell et al. 2010; Harper et al. 2018) on the 2.7 m Stratospheric Observatory For Infrared Astronomy (SOFIA) telescope. HAWC+ polarimetric observations simultaneously measure two orthogonal components of linear polarization arranged in two arrays of 32 × 40 pixels each, with pixel scales of 255 and 402 pixel⁻¹, and beam sizes (full width at half maximum) of 485 and 780 at 53 and 89 μm, respectively. For both bands, we performed observations in a four-dither square pattern with a distance of three pixels from the center position in the array coordinate system. In each dither position, four halfwave plate PAs were taken. We used a chop frequency of 10.2 Hz, a nod time of 45 s, a chop-throw of 180'', and a chop-angle of 90° to always be along the short axis of the 32 × 40 pixels array. We used the science instrument reference frame for these observations. At 53 μm, we observed a total of eight dither positions, and a total of 63 dither positions at 89 μm. The total on-source integration times were 0.47^h and 2.03^h at 53 and 89 μm, respectively. The total clock time of the observation's on-source time plus off-source time plus overheads is 1.3^h and 5.2^h, respectively. Based on the morphology of the source using images from the Herschel Space Observatory²⁴ (Pilbratt et al. 2010) and the chop-throw of 180'' used for the HAWC+ observations, the reference beam fluxes should be ≤1/10 times the source flux over a source diameter of 120'' and therefore negligible. Data were reduced using the hawc_dpr pipeline v1.3.0. The pipeline procedure described by Harper et al. (2018) was used to background-subtract and flux-calibrate the data and compute Stokes parameters and their uncertainties. Computation of final degree and PA of polarization accounts for correction of the instrumental polarization (Table 4 by Harper et al. 2018) with typical standard deviations after subtraction of ∼0.3%, and debiasing and PA error estimation (Wardle & Kronberg 1974). Further analyses and high-level displays were done with custom python routines.

Figure 1 shows the total flux image with overlaid polarization vectors in a 2' × 2' (7.8 × 7.8 kpc²) field-of-view (FOV) at 89 μm. Polarization vectors have been rotated by 90° to show the morphology of the magnetic field. We note that the polarization orientations are ambiguous by 180°, thus the displayed polarization vectors for all figures should be interpreted as magnetic field lines. Each polarization vector shows a statistically independent polarization measurement with P/σ_P ≥ 3. Polarization vectors are shown every 402 (Nyqvist sampling at 89 μm). Throughout the paper we assume that dust grain alignment is perpendicular to the direction of the magnetic field. We observe extended dust emission along the spiral arms of the galaxy up to a diameter of ∼100'' (∼6.5 kpc) centered at the nucleus while the central ∼1 kpc radius emission is extended along PA of ∼45° and cospatial with the inner-bar. The polarization map shows a spiral shape with a diameter of ∼50'' (∼3.25 kpc). This spiral shape can be clearly identified in Figure 2 (left) in the image made using the line integral convolution contour (LIC) technique (Cabral & Leedom 1993). This figure uses all polarization vectors regardless of signal-to-noise ratio measured in our observations. Figure 2 (right) shows combined HST observations at V-band, I-band, and H_α with overlaid magnetic field vectors (blue vectors) from our 89 μm HAWC+ observations. The polarized flux is shown in Figure 1 (right). The polarized flux shows extended emission along the E–W direction with a peak shifted by 4'' (Section 4.4) NE from the peak of the total flux intensity.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Given the short integration time at 53 μm, only polarization at the nucleus (one independent vector is measured) of NGC 1068 was detected, and the polarization map is not shown. For the total intensity image at 53 μm, we refer to the HAWC+ observations using the Lissajous observing mode presented by Lopez-Rodriguez et al. (2018b). We here report the measured nuclear degree of polarization of 1.3% ± 0.3% and PA of polarization of 139° ± 6° (E-vector) within a 5'' (325 pc) diameter at 53 μm. For 89 μm, we measured a nuclear polarization of 0.6% ± 0.3% and PA of polarization of 156° ± 13° (E-vector) within a 8'' (520 pc) diameter. For both bands, the measured nuclear polarization is below the instrumental polarization and above the residual polarization after instrumental polarization correction. Thus, the measured polarization is consistent with a low polarized core. As the fractional contribution of the AGN increases at short wavelengths and the dust emission from the host galaxy increases at long wavelengths (Figure 6 by Lopez-Rodriguez et al. 2018b), we expect to measure a decrease of the nuclear degree of polarization as wavelength increases. Our measurements are compatible with this behavior.

The observational results shown in Section 2 suggest the presence of a large-scale magnetic field along two spiral arms. We here produce a model of the magnetic field to characterize this structure. We use an axisymmetric spiral structure (i.e., Braun et al. 2010; Ruiz-Granados et al. 2010) defined in Cartesian coordinates by the form

$$\begin{eqnarray}&&{B}_{x}=-{B}_{0}(r)\sin (\phi +{\rm{\Psi }})\cos \chi (z)\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{B}_{y}={B}_{0}(r)\cos (\phi +{\rm{\Psi }})\cos \chi (z)\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{B}_{z}={B}_{0}(r)\sin \chi (z)\end{eqnarray} \tag{ 3 }$$

where B₀(r) is the magnetic field as a function of the radial distance from the core of the galaxy, and Ψ is the pitch angle defined as the angle between the azimuthal (ϕ) direction and the magnetic field direction. χ(z) is the out-of-plane magnetic field taken as $\chi (z)={\chi }_{0}\tanh (z/{z}_{0})$ (Ruiz-Granados et al. 2010).

When this magnetic field configuration is viewed at an inclination, i, and tilt angle,²⁵ θ, of the galaxy 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{ 4 }$$

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

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

where (${x}^{s},{y}^{s},{z}^{s}$) are the major axis, minor axis, and the line of sight (LOS), respectively, in the sky coordinate system.

The PA of polarization in the plane of the sky, which is parallel to the direction of the magnetic field in the (${x}^{s},{y}^{s},{z}^{s}$) coordinate system, is computed as ${\mathrm{PA}}_{B}=\arctan ({B}_{{y}^{s}}/{B}_{{x}^{s}})$. The degree of polarization as a function of the inclination and magnetic field direction in the plane of the sky is computed as $P={p}_{o}(1-{({B}_{{z}^{s}}/{B}_{s})}^{2})$, where p_o is a constant factor, and ${B}_{s}^{2}={B}_{{x}^{s}}^{2}+{B}_{{y}^{s}}^{2}+{B}_{{z}^{s}}^{2}$ is the total magnetic field. We assume a negligible contribution of the out-of-plane magnetic field, χ₀ = 0, thus our final magnetic field is coplanar with the disk of the galaxy. This model is purely a geometrical description and does not fulfill the divergence-free vector field as B_z is neglected.

Since the FIR polarimetry is not directly sensitive to the magnetic field strength, there are three free model parameters: pitch angle, Ψ, inclination angle, i, and tilt axis PA, θ. The face-on view corresponds to i = 0° and edge-on to i = 90°. Tilt axis has a reference point, θ = 0, along the north–south direction and positively increases east from north. To fit for model parameters, we computed a Markov Chain Monte Carlo approach using the differential evolution metropolis sampling step in the python code PyMC3 (Salvatier et al. 2016). The prior distributions are set to flat within the ranges Ψ = (0°, 50°), i = (0°, 90°), and θ = (0°, 90°). We run the code using 15 chains with 6000 steps and a 1000 step burn-in per chain. Because the central 83 × 83 (0.54 × 0.54 kpc²) region of NGC 1068 is dominated by the inner-bar, we have excluded this region in this analysis. One could perhaps question the size of the exclusion zone, but the overall effect on the model fit (presented above) of the small number of vectors within the innermost annulus will be small.

Figure 3 shows the marginalized posterior distributions of the free parameters and their maximum a posteriori (MAP) with their 1σ estimations. This figure also shows, as a function of azimuthal angle, the measured degree of polarization, the measured PA of polarization Δϕ (black dots) in radial coordinates centered on the galaxy nucleus, the best-fit model (orange solid lines), and 1σ uncertainties (blue shadow region). If Δϕ = 0°, the angle is perfectly radial, and if Δϕ = ±90°, the angle is perfectly azimuthal. To generate the data points in the figure, the Stokes I, Q, and U images were binned by azimuthal angle in bins of 10°. Then, debiased degree and PA of polarization were estimated for each bin.

Zoom In Zoom Out Reset image size

We find a symmetric pattern of the polarization as a function of the azimuthal angle, reflecting the relationship among degree of polarization, polarization PA, and model parameters shown above. An axisymmetric spiral structure with a pitch angle of ${\rm{\Psi }}={16.9}_{-2.8}^{+2.7\circ}$, inclination $i={48.1}_{-1.9}^{+1.8\circ}$, and tilt angle of θ = 378 ± 17 best describes our data. The host galaxy is inclined by 40° (de Vaucouleurs et al. 1991), and using the H i velocity field, Brinks et al. (1997) estimated an inclination angle of 40° ± 3°, which both are fairly close to our estimated values. The main difference may be that de Vaucouleurs et al. (1991) estimated the galaxy inclinations using isophotes on the optical images, and that Brinks et al. (1997) used the rotation curve of the H i emission observations while our estimation is based on a smaller region of the galaxy dominated by the FIR emission. The PA of the kinematic major axis of H i is estimated to vary from 270° ± 3° in the inner 30''–70'' disk to 286° ± 5° beyond 100'', which is almost east–west (Brinks et al. 1997). We further discuss the pitch angle in Section 4.3.

Figure 4 (left) shows the magnetic field configuration (blue vectors) of the best inferred parameters (Figure 3) with the overlaid measured PA of polarization rotated by 90° (black vectors) of the central 3 × 3 kpc² of NGC 1068. Figure 4 (right) shows the total flux (color scale) image of NGC 1068 with overlaid magnetic field model (streamlines) using LIC techniques and the measured polarization vectors (black lines). All polarization vectors have the same length, and their orientation shows the magnetic field morphology. The large-scale magnetic field morphology agrees with an axisymmetric spiral structure. The largest difference between our model and the observations is located within ∼1 kpc of the nucleus (Section 4.4). As our axisymmetric magnetic field model does not fulfill the divergence-free requirement, we cannot reproduce the field in the inner areas of the galaxy. The deviations within ∼1 kpc around the nucleus may indicate a limitation of the model and/or a different magnetic field morphology as described in Section 4.4.

Zoom In Zoom Out Reset image size

4.1. Temperature and Column Density Maps

To support the analysis of the following sections, we have constructed temperature and column density maps of NGC 1068. We have combined our 89 μm HAWC+ observations with archival Herschel Space Observatory (Pilbratt et al. 2010) observations at 70, 160, and 250 μm taken with the PACS (Poglitsch et al. 2010) and SPIRE (Griffin et al. 2010) instruments. We binned each observation to the pixel scale, 6'', of the 160 μm Herschel data. Then we extracted the intensity values of each pixel associated with the same part of the sky at each wavelength. Finally, we fit a modified blackbody function assuming a dust emissivity of ${\epsilon }_{\lambda }\propto {\lambda }^{1.6}$ (e.g., Boselli et al. 2012) and with the temperature and optical depth, τ, as free parameters. We estimate the column density using the relation ${N}_{{\rm{H}}{\rm{I}}+{{\rm{H}}}_{2}}=\tau /(k\mu {m}_{{\rm{H}}})$, where τ is the optical depth at 250 μm from the best-fit modified blackbody function at a given pixel, k = 0.1 cm² g⁻¹ (Hildebrand 1983), μ = 2.8 is the mean molecular mass per H molecule (Kauffmann et al. 2008), and m_H is the hydrogen mass. Figure 5 shows the estimated temperature and column density distributions within the same FOV (2' × 2', 7.8 × 7.8 kpc) as Figure 1.

Zoom In Zoom Out Reset image size

The estimated dust temperatures lie in the range of 20–34 K and peak in a ridge along the inner-bar. There is a strong peak temperature (34 ± 2 K) displaced ∼7'' (∼0.45 kpc) NE from the position of the AGN (black cross). The dust temperature has another peak ∼45'' (∼2.93 kpc) SW spiral arm at the location of a large complex of H ii regions and giant molecular clouds (see Figure 2 (right), Kaneko et al. 1989; Tosaki et al. 2017). Figure 6 shows the location of the H ii regions (bright compact sources) using the HST/F555W (λ_c = 0.54 μm, Δλ = 0.11 μm) filter to compare with our inferred magnetic field morphology.

Zoom In Zoom Out Reset image size

The column densities in the mapped area lie in the range ${10}^{20}\mbox{--}{10}^{22.4}$ cm⁻², corresponding to optical extinctions of A_V ∼ 0.05–14 mag. The distribution is extended along the inner-bar but is not as strongly peaked along the ridge as the temperatures. There is a weak maximum SW of the nucleus, and at the lowest contour levels there is an extension along the SW spiral arm. The morphology is consistent with H i absorption maps observed with the Very Large Array with a beam size of 21 × 10 (Gallimore et al. 1994). They found column densities in the range [1–4] × 10²¹ cm⁻² in the region SW of the nucleus. Integrating over the column density map, we derive a mass of $[6.1\,\pm \,1.4]\,\times \,{10}^{8}$ M_⊙ in the central 2' × 2' (7.8 × 7.8 kpc).

4.2. On the Role of Magnetic Fields in the Spiral Arms

The main result of these observations is the detection of coherent magnetic fields over a 3 kpc diameter region. Our column density map shows that all points have optical depths <0.2 at 89 μm. We converted the estimated hydrogen density, N_H, from Figure 5 to the visual extinction, A_V, using the standard ratio, ${A}_{V}/{N}_{{\rm{H}}}=5.35\times {10}^{-22}$ mag cm⁻² (Bohlin et al. 1978). Then, we used the typical extinction curve²⁶ of the Milky Way for R_V = 3.1 (Weingartner & Draine 2001) to estimate the conversion ${\tau }_{0.55}/{\tau }_{89}=4\times {10}^{-3}$ between optical depths at 0.55 and 89 μm. Hence, at each LOS it is likely that we are sampling the integrated flux from aligned dust grains all the way through the galactic disk. In the outer regions of the map, the optical extinctions are also low to moderate, and there is a close alignment between visual tracers of spiral arms and our magnetic field vectors (see Figure 1 right, and Figure 6 left). Within a diameter of 10 kpc, CO is also aligned with the optical spiral structure, and H_α residual velocities reveal streaming motions along the arms (Dehnen et al. 1997).

While we discussed above some of the shortcomings of mapping magnetic fields with optical polarimetry, the V-band vectors of Figure 2 in Scarrott et al. (1991) or the data shown by Ward-Thompson (1987) line up well with our mathematical model. In the outer regions (r > 15'', outside the starburst ring) that is also true for our FIR polarimetry. Inside the starburst ring, where the optical depths are high (Figure 5-top), the V-band degree of polarization is very small, but in general the directions follow our model. The FIR emission, on the other hand, is more highly polarized and there are significant deviations from the model (see Section 4.4). From the above, we infer (1) that our FIR polarimetry results are, indeed, tracing spiral structure, (2) that, as expected, they seem to offer significant advantages over optical polarimetry for measuring the bulk properties of the ISM, particularly along lines of sight with high optical depth, and (3) that the field lines appear to follow the directions of gas flows along spiral arms, which may indicate that the arms have strong shocks otherwise the field would cross the arms at a given angle.

The magnetic field in the ISM is often modeled using a combination of constant and turbulent components. The trend of decreasing fractional polarization with increasing column density (Hildebrand et al. 1999; Jones 2015; Jones & Whittet 2015) provides an indirect measurement of the effect of the turbulent component. For maximally aligned dust grains along an LOS with a constant magnetic field direction, the degree of polarization, P, in emission will be constant with optical depth τ. If there is a region along the LOS with some level of dust grain misalignment, it will result in a reduced degree of polarization. If the magnetic field direction in a region varies completely randomly along the LOS with a well-defined scale length in τ, we have $P\propto {\tau }^{-0.5}$ (Jones 2015). For any combination of constant and random components, the slope will be between these two limits (Jones et al. 2015). However, if (1) grain alignment decreases in denser regions, or (2) there is coherent cancellation of polarization (e.g., crossed field directions), or (3) there are regions of very high turbulence on very small scale lengths, the polarization can decline faster with optical depth than a slope of −1/2.

Figure 7 shows log/log plots of polarization versus column density and surface brightness for all polarization vectors with σ_P < 5%. The column density and dust temperature for each polarization vector is shown in color scale. It is clear that the polarization is higher when the temperature, column density, and fluxes are low, although the spatial anticorrelation of temperature and column density (Figure 5) are the more fundamental relationships. This behavior has been previously observed in molecular clouds complexes in our Galaxy and Galactic center (e.g., Hildebrand et al. 1999; Chuss et al. 2003), which shows a depolarization effect as flux increases. The upper envelope of polarization measurements provides the least depolarized LOSs on the source, which allows us to find the overall trend of the optical depth index. We use a linear fit to the log–log plot along the upper envelope of the measurements, which gives an optical depth index of −0.5 (black dashed line). This result indicates that most of the polarization measurements suffer from depolarization with an optical depth index steeper than −0.5, which may be due to high turbulence at small scales, effects of dust grain alignment toward dense regions and/or cross fields on the LOS.

Zoom In Zoom Out Reset image size

The degree of polarization in our map is typically <3%, with a peak of ∼7% (with P/σ_P ∼ 3). The polarization is systematically higher in the outer regions than in the starburst disk at radius ≤25'' (1.6 kpc). FIR polarimetric observations of Galactic clouds typically measure degree of polarization of 3%–6% in the 60–100 μm wavelength range (Dotson et al. 2000). For our observations of NGC 1068, although we have an ordered polarization spiral pattern, we identify a potential disorder in the magnetic field on scales smaller than the resolution element, 8'' (520 pc), of our observations. In a sense, this should not be surprising. The forces creating large-scale structures like spiral arms result from gravity and rotation. Where there are additional local sources of energy such as those arising from the assembly and disruption of star-forming complexes, there will be multiple additional mechanisms for distorting preexisting magnetic fields or generating new ones with orientations that differ from the large-scale field. The larger the fraction of the ISM involved in these processes, the greater the average misalignment with the general field and the greater the depolarization of the emergent FIR radiation. This is consistent with studies that show that Galactic star-forming complexes with scale lengths of tens of parsecs are misaligned with the Galaxy's field on kiloparsec scales (Stephens et al. 2011). One of the most significant aspects of our observations may be that in spite of the extremely high star formation rate in the starburst disk of NGC 1068, there is still a clearly discernible spiral field over much of its extent.

The alignment efficiency of dust grains may be an issue in very dense molecular cloud cores where the extinction exceeds A_V ∼ 20 and no embedded young stellar objects (YSOs) illuminate the dust (Jones et al. 2015). To determine the effect of loss of grain alignment on our observations, we need to calculate what fraction of the flux in our 53 and 89 μm beams could arise from starless or pre-stellar molecular cores. Using ¹³CO and 1.1 mm dust continuum observations of the Milky Way, Battisti & Heyer (2014) find that only about 10% or less of the mass in Giant Molecular Clouds (GMC) is contained in compact regions with ${A}_{v}\gt 10$ and temperatures of T ∼ 10 K. Given that some of these regions may contain a YSO that could produce enough radiation to align grains (e.g., Jones et al. 2016), this must be considered an upper limit to the mass fraction in starless cores. Given this upper limit on the mass fraction of GMCs in the form of starless cores and the 20–34 K range of temperatures we derive for the dust emission in our beam (Figure 5) for NGC 1068, the fraction of the flux in our beam from cloud cores with unaligned grains must be less than 1%. Thus, in our 53 and 89 μm beams (5''–8''), loss of grain alignment cannot be a significant contributing factor to the trend seen in Figure 7. This leaves the coherent cancellation of polarization and/or regions of very high turbulence on small scale lengths as potential explanations for the slope in Figure 7.

4.3. The Pitch Angle and Masses

As the density wave passes through the host galaxy, it enhances the star formation. The spiral density wave theory predicts that the pitch angle may vary with wavelength, as described in Section 1. NGC 1068 has been studied across the electromagnetic spectrum where the pitch angle of the arms has been estimated to be 206 ± 45 at 0.46 μm (Berrier et al. 2013), 173 ± 22 using R band images with the 2.5 m Las Cumbres Observatory (LCO; Seigar et al. 2008), 15° using the H_α velocity field (Emsellem et al. 2006), and 7°–10° using the CO (J = 1 − 0) emission (Planesas et al. 1991). We estimated a pitch angle of ${16.9}_{-2.8}^{+2.7\circ}$ in the magnetic field inferred from 89 μm observations, which is between the pitch angles estimated in the spiral arms at 0.46 μm and the CO (J = 1 − 0) emission. The material in the galaxy which is sampled by our FIR measurements is spatially coincident with the H_α velocity field (see Figure 7 by Emsellem et al. 2006). The analysis by Davies et al. (1998) computed that the H_α velocity field pattern is spatially correlated with the starbursting H ii regions which was enhanced by a burst in the past 30 Myr. Later, Emsellem et al. (2006) found that the H_α and CO intensity distributions are not spatially coincident due to star formation and extinction. We conclude that the matter sampled by our estimated pitch angle follows the starbursting regions along the spiral arms.

The spiral density wave theory has also shown that the concentration of host galaxy mass and nuclear mass determines the pitch angle (Roberts et al. 1975). As a consequence, recent studies have gone further and empirically shown that the spiral arm pitch angle, Ψ, and the central black hole mass, M_BH, may be correlated in the modal density wave (e.g., Seigar et al. 2008; Berrier et al. 2013) by the form

$$\begin{eqnarray}&&\mathrm{log}\left(\displaystyle \frac{{M}_{\mathrm{BH}}}{{M}_{\odot }}\right)=(8.21\pm 0.16)-(0.062\pm 0.009){\rm{\Psi }}\end{eqnarray} \tag{ 7 }$$

(Berrier et al. 2013). Using our estimated pitch angle of ${16.9}_{-2.8}^{+2.7\circ}$, we estimate a nuclear mass of $\mathrm{log}(M/{M}_{\odot })={7.16}_{-0.51}^{+0.46}$. Our result is in agreement, within the uncertainties, with the black hole mass of $\mathrm{log}({M}_{\mathrm{BH}}/{M}_{\odot })=6.95\pm 0.02$ estimated using maser modeling on the central pc of NGC 1068 (Lodato & Bertin 2003).

4.4. The Central Region of NGC 1068

We have presented the first FIR polarimetry of NGC 1068. We found a large-scale (∼3 kpc) spiral pattern that we attribute to thermal emission from magnetically aligned grains. There is a spatial and morphological correspondence between the 89 μm magnetic field vectors and other tracers, i.e., O iii, H_α, optical, of spiral arms. We found a symmetric polarization pattern as a function of the azimuthal angle arising from the projection and inclination of the disk field component in the plane of the sky. We are able to explain this behavior with an axisymmetric spiral polarization pattern with pitch angle ${16.9}_{-2.8}^{+2.7\circ}$, inclination ${48.1}_{-1.9}^{+1.8\circ}$, and tilt angle 378 ± 17. The matter sampled by our estimated pitch angle follows the starbursting regions along the spiral arms, and, using the pitch angle-mass relationship, predicts a nuclear mass of $\mathrm{log}(M/{M}_{\odot })={7.16}_{-0.51}^{+0.46}$.

Outside the starburst disk (radius ≤1.6 kpc) in NGC 1068, the degree of polarization is similar to that seen in Milky Way sources. At smaller radii, it decreases with flux, dust temperature, and column density, and it has two minima near the locations of the ends of the NIR inner-bar in the starburst ring. This trend is consistent with dilution of the net spiral field by random fields created by injection of kinetic energy into the ISM by active star formation.

Inside the starburst ring (<1.6 kpc), we found evidence for large-scale coherent magnetic fields that align with visual tracers but not with our model field. The discrepancies are largest along the leading edges of the inner-bar. If the magnetic field traces gas flows, this is consistent with inward transport of gas induced by the bar. The intensity of polarized flux is stronger along PAs centered at ∼70°, suggesting that at the present time the strongest flows are along the NE branch of the inner-bar.

A peak in polarized flux intensity and dust temperature occurs at ∼4'' (∼0.26 kpc) NNE of the AGN near the location of the bow shock separating the AGN outflow from the surrounding ISM. This is consistent with the hypothesis that the magnetic field has been amplified at the shock interfaces along the edges of the outflow cavity, but our current angular resolution is insufficient to rule out the possibility that the field is intrinsic to the interstellar clouds with which the outflow is interacting. There is a minimum in polarized flux to the SSW of the AGN at the location of the counter-jet, possibly because of a current deficit of ISM in that direction.

The degree of FIR polarization at the position of the AGN is low, 1.3% ± 0.3% within 5'' (0.33 kpc) diameter at 53 μm and 0.6% ± 0.3% within 8'' (0.52 kpc) diameter at 89 μm. This is much lower than the ∼10% FIR polarization of the radio-loud AGN in Cygnus A (Lopez-Rodriguez et al. 2018a). The degree of polarization of NGC 1068 integrated over an 8 kpc diameter region comparable to the physical scale of the region encompassed by the Cygnus A measurement is only 0.7%.

The results presented here, along with our prior studies of M82 and NGC 253 (Jones et al. 2019), provide evidence that FIR polarimetry can be a valuable tool for studying magnetic field structure in external galaxies, particularly in regions of high optical depth.

We are grateful to the anonymous referee who helped to clarify and improve the text. Based on observations made with the NASA/DLR Stratospheric Observatory for Infrared Astronomy (SOFIA) under the GTO Program. SOFIA is jointly operated by the Universities Space Research Association, Inc. (USRA), under NASA contract NAS2-97001, and the Deutsches SOFIA Institut (DSI) under DLR contract 50 OK 0901 to the University of Stuttgart. Portions of this work were carried out at the Jet Propulsion Laboratory, operated by the California Institute of Technology under a contract with NASA. PACS has been developed by a consortium of institutes led by MPE (Germany) and including UVIE (Austria); KU Leuven, CSL, IMEC (Belgium); CEA, LAM (France); MPIA (Germany); INAF-IFSI/OAA/OAP/OAT, LENS, SISSA (Italy); IAC (Spain). This development has been supported by the funding agencies BMVIT (Austria), ESA-PRODEX (Belgium), CEA/CNES (France), DLR (Germany), ASI/INAF (Italy), and CICYT/MCYT (Spain). SPIRE has been developed by a consortium of institutes led by Cardiff University (UK) and including Univ. Lethbridge (Canada); NAOC (China); CEA, LAM (France); IFSI, Univ. Padua (Italy); IAC (Spain); Stockholm Observatory (Sweden); Imperial College London, RAL, UCL-MSSL, UKATC, Univ. Sussex (UK); and Caltech, JPL, NHSC, Univ. Colorado (USA). This development has been supported by national funding agencies: CSA (Canada); NAOC (China); CEA, CNES, CNRS (France); ASI (Italy); MCINN (Spain); SNSB (Sweden); STFC, UKSA (UK); and NASA (USA).

Facilities: SOFIA (HAWC+) - , Herschel (PACS - , SPIRE) - , HST (WFPC2 - , ACS). -

Software: astropy (Astropy Collaboration et al. 2013), PyMC3 (Salvatier et al. 2016), aplpy (Robitaille & Bressert 2012), matplotlib (Hunter 2007).