SOFIA/HAWC+ traces the magnetic fields in NGC 1068

TL;DR

This study presents the first detection of galactic spiral structure via thermal dust emission polarization in NGC 1068, revealing magnetic field alignment with spiral arms and detailed magnetic topology near the galaxy's core and AGN.

Contribution

It provides new high-resolution polarimetric imaging of NGC 1068 at 89 μm, demonstrating magnetic field structures aligned with spiral arms and modeling the magnetic field geometry in an active galaxy.

Findings

Magnetic fields closely follow spiral arms in NGC 1068.

Polarization pattern consistent with a logarithmic spiral model.

Magnetic field deviations near the starburst ring and AGN.

Abstract

We report the first detection of galactic spiral structure by means of thermal emission from magnetically aligned dust grains. Our 89 $\mu$m polarimetric imaging of NGC 1068 with the High-resolution Airborne Wideband Camera/Polarimeter (HAWC+) on NASA's Stratospheric Observatory for Infrared Astronomy (SOFIA) also sheds light on magnetic field structure in the vicinity of the galaxy's inner-bar and active galactic nucleus (AGN). We find correlations between the 89 $\mu$m magnetic field vectors and other tracers of spiral arms, and 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. The observations can be fit with a logarithmic spiral model with pitch angle of $16.9^{+2.7}_{-2.8}$$^{\circ}$ and a disk inclination of $48\pm2^{\circ}$. We infer that the bulk of the interstellar medium from which the polarized dust emission originates is threaded by a magnetic field that closely follows the spiral arms. Inside the central starburst disk ($<1.6$ kpc), the degree of polarization is found to be lower than for far-infrared sources in the Milky Way, and has minima at the locations of most intense star formation near the outer ends of the inner-bar. Inside the starburst ring, the field direction deviates from the model, becoming more radial along the leading edges of the inner-bar. The polarized flux and dust temperature peak $\sim 3-6$" NE of the AGN at the location of a bow shock between the AGN outflow and the surrounding interstellar medium, but the AGN itself is weakly polarized ($< 1$%) at both 53 and 89 \um.

Equations

$$B_x$$

$$B_y$$

$$B_z$$

$$B_{x^s}$$

$$B_{y^s}$$

$$B_{z^s}$$

$$\log (\frac{M_{BH}}{M_{\odot }})=(8.21\pm 0.16)-(0.062\pm 0.009)\Psi$$

Full text

infrared: galaxies - techniques: polarimetric - galaxies: individual (NGC 1068) - galaxies: magnetic fields

††facilities: SOFIA (HAWC+), Herschel (PACS, SPIRE), HST (WFPC2, ACS)††software: astropy (Astropy Collaboration et al., 2013), PyMC3 (Salvatier, Wiecki & Fonnesbeck, 2016), aplpy (Robitaille & Bressert, 2012), matplotlib (Hunter, 2007).

1 Introduction

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, Lazarian & Vaillancourt, 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 M 51 (Beck, Klein & Wielebinski, 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 show similarities between the optical and radio polarization (Fendt, Beck & Neininger, 1998). Further observations of M 51 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 $(\sim\lambda^{-4})$ between $0.55$ and $1.65~{}\micron$ than its absorption ($\sim\lambda^{-1}$, 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 $\sim 0.4$%, but an upper limit of $0.05\%$ was measured at H-band (Pavel & Clemens, 2012). The dichroic absorptive polarization should has shown a higher polarization at H-band than the measured to be due to magnetically aligned dust grains. 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, HII 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 $\sim 20\arcsec$ ($\sim 1.3$ kpc) in diameter at a position angle of $\sim 40^{\circ}$ 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\arcsec$ ($>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 is consistent with the $40^{\circ}$ inclination suggested by HI data. However, there is a bright oval structure $\sim 180\arcsec$ ($\sim 11.7$ kpc) diameter with long axis perpendicular to the kinematic axis that probably corresponds to a large-scale stellar bar (fig. 6 of Schinnerer et al., 2000). The inner Lindblad resonance (ILR) of this bar has a de-projected radius of $\sim 18\arcsec$ ($\sim 1.17$ kpc), near the location of a compact spiral or ring-like structure in 12CO, 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\arcsec$ ($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 $\sim 1\arcsec$ ($\sim 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, Yagi & Taniguchi, 2017). The optical surface brightness of the galactic disk at radii $\leq 25\arcsec$ ($\leq 1.63$ kpc)  (henceforth the “starburst disk”) is among the brightest of any galaxy in the local universe (Keel & Weedman, 1978; Weedman, 1985; Ichikawa, Okamura & Kaneko, 1987). Observations in the FIR have shown that the total luminosity of the starburst disk is $>10^{11}L_{\odot}$ (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 non-thermal synchrotron emission produced in supernova explosions of massive young stars. On the basis of 10.8 $\mu$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\arcmin\times 2\arcmin$ ($7.8\times 7.8$ kpc2) of NGC 1068 at wavelengths of 53 and 89 $\mu$m. We discuss the 53 $\mu$m flux from the AGN in a previous paper (Lopez-Rodriguez et al., 2018a). In this paper, we present our polarimetric results, including an 89 $\mu$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.

2 Observations and Data Reduction

NGC 1068 was observed at $53$ and $89$ $\mu$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\times 40$ pixels each, with pixel scales of 2$$\farcs$$55 and 4$$\farcs$$02 pixel*-1*, and beam sizes (full width at half maximum, FWHM) of 4$$\farcs$$85 and 7$$\farcs$$80 at $53$ and $89$ $\mu$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 position angles (PA) 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^{\circ}$ to always be along the short axis of the $32\times 40$ pixels array. We used the science instrument reference frame (SIRF) for these observations. At $53$ $\mu$m, we observed a total of $8$ dither positions, and a total of $63$ dither positions at $89$ $\mu$m. The total on-source integration times were $0.47$h and $2.03$h at $53$ and $89$ $\mu$m, respectively. The total clock time of the observations—on-source time plus off-source time plus overheads—is 1.3h and 5.2h, respectively. Based on the morphology of the source using images from the Herschel Space Observatory111Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA (Pilbratt et al., 2010) and the chop-throw of $180\arcsec$ used for the HAWC+ observations, the reference beam fluxes should be $\leq 1/10$ times the source flux over a source diameter of $120\arcsec$ 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 $\sim 0.3\%$, and de-biasing 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\arcmin\times 2\arcmin$ ($7.8\times 7.8$ kpc2) field-of-view (FOV) at $89$ $\mu$m. Polarization vectors have been rotated by $90^{\circ}$ to show the morphology of the magnetic field. We note that the polarization orientations are ambiguous by $180^{\circ}$, 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/\sigma_{P}\geq 3$. Polarization vectors are shown every $4.02\arcsec$ (Nyqvist sampling at 89 $\mu$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 $\sim$100″ ($\sim 6.5$ kpc) centered at the nucleus while the central $\sim 1$ kpc radius emission is extended along PA of $\sim 45^{\circ}$ and co-spatial with the inner-bar. The polarization map shows a spiral shape with a diameter of $\sim 50$″ ($\sim 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 $\mu$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\arcsec$ (Section 4.4) NE from the peak of the total flux intensity.

Given the short integration time at $53$ $\mu$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 $\mu$m, we refer to the HAWC+ observations using the Lissajous observing mode presented by Lopez-Rodriguez et al. (2018a). We here report the measured nuclear degree of polarization of $1.3\pm 0.3$% and PA of polarization of $139\pm 6^{\circ}$ (E-vector) within a $5$″ ($325$ pc) diameter at $53$ $\mu$m. For 89 $\mu$m, we measured a nuclear polarization of $0.6\pm 0.3$% and PA of polarization of $156\pm 13^{\circ}$ (E-vector) within a $8\arcsec$ ($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 (fig. 6 by Lopez-Rodriguez et al., 2018a), we expect to measure a decrease of the nuclear degree of polarization as wavelength increases. Our measurements are compatible with this behavior.

3 Model of the large-scale galactic magnetic field

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. Ruiz-Granados et al., 2010; Braun, Heald & Beck, 2010) defined in cartesian coordinates by the form

[TABLE]

where $B_{0}(r)$ is the magnetic field as a function of the radial distance from the core of the galaxy, and $\Psi$ is the pitch angle defined as the angle between the azimuthal ($\phi$) direction and the magnetic field direction. $\chi(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 angle222Also described as the PA of the major axis of the projected galaxy plane., $\theta$, of the galaxy projected on the plane of the sky, the magnetic field at the observer’s frame is given by

[TABLE]

where ($x^{s},y^{s},z^{s}$) are the major axis, minor axis, and the line of sight, 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 $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^{2}_{s}=B^{2}_{x^{s}}+B^{2}_{y^{s}}+B^{2}_{z^{s}}$ is the total magnetic field. We assume a negligible contribution of the out-of-plane magnetic field, $\chi_{0}=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, $\Psi$, inclination angle, $i$, and tilt axis position angle, $\theta$. Face-on view corresponds to $i=0^{\circ}$ and edge-on to $i=90^{\circ}$. Tilt axis has a reference point, $\theta=0$, along the north-south direction and positively increases east from north. To fit for model parameters, we computed a Markov Chain Monte Carlo (MCMC) approach using the differential evolution metropolis sampling step in the python code PyMC3 (Salvatier, Wiecki & Fonnesbeck, 2016). The prior distributions are set to flat within the ranges $\Psi$ = ($0^{\circ}$,$50^{\circ}$), $i$ = ($0^{\circ}$,$90^{\circ}$), and $\theta$ = ($0^{\circ}$,$90^{\circ}$). We run the code using 15 chains with 6000 steps and a 1000 steps burn-in per chain. Because the central $8\farcs 3\times 8\farcs 3$ ($0.54\times 0.54$ kpc2) 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$\sigma$ estimations. This figure also shows, as a function of azimuthal angle, the measured degree of polarization, the measured PA of polarization $\Delta\phi$ (black dots) in radial coordinates centered on the galaxy nucleus, the best fit model (orange solid lines), and 1$\sigma$ uncertainties (blue shadow region). If $\Delta\phi=0^{\circ}$, the angle is perfectly radial, and if $\Delta\phi=\pm 90^{\circ}$, 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.

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 $\Psi=16.9^{+2.7}_{-2.8}$∘, inclination $i=48.1^{+1.8}_{-1.9}$∘, and tilt angle of $\theta=37.8\pm 1.7$∘ best describes our data. The host galaxy is inclined by 40*∘* (de Vaucouleurs et al., 1991), and using the HI velocity field, Brinks et al. (1997) estimated an inclination angle of $40\pm 3^{\circ}$, 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 HI emission observations while our estimation is based on a smaller region of the galaxy dominated by the FIR emission. The position angle of the kinematic major axis of HI is estimated to vary from $270\pm 3^{\circ}$ in the inner $30-70\arcsec$ disk to $286\pm 5^{\circ}$ beyond $100\arcsec$, 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 (Fig. 3) with the overlaid measured PA of polarization rotated by $90^{\circ}$ (black vectors) of the central $3\times 3$ kpc2 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 $\sim 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 $\sim 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.

4 Discussion

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$ $\mu$m HAWC+ observations with archival Herschel Space Observatory (Pilbratt et al., 2010) observations at $70$, $160$, and $250$ $\mu$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$ $\mu$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, $\tau$, as free parameters. We estimate the column density using the relation $N_{HI+H_{2}}=\tau/(k\mu m_{H})$, where $\tau$ is the optical depth at 250 $\mu$m from the best fit modified blackbody function at a given pixel, $k=0.1$ cm2 g*-1* (Hildebrand, 1983), $\mu=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\arcmin~{}\times~{}2\arcmin$, $7.8~{}\times~{}7.8$ kpc) as Fig. 1.

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\pm 2$ K) displaced $\sim 7\arcsec$ ($\sim 0.45$ kpc) NE from the position of the AGN (black cross). The dust temperature has another peak $\sim 45\arcsec$ ($\sim 2.93$ kpc) SW spiral arm at the location of a large complex of HII regions and giant molecular clouds (see Fig. 2-right, Kaneko et al., 1989; Tomoka et al., 2017). Fig. 6 shows the location of the HII regions (bright compact sources) using the HST/F555W ($\lambda_{c}=0.54$ $\mu$m, $\Delta\lambda=0.11$ $\mu$m) filter to compare with our inferred magnetic field morphology.

The column densities in the mapped area lie in the range $10^{20}-10^{22.4}$ cm*-2*, corresponding to optical extinctions of A${}_{V}\sim 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 HI absorption maps observed with the VLA with a beam size of $2\farcs 1~{}\times~{}1\farcs 0$ (Gallimore et al., 1994). They found column densities in the range $[1-4]~{}\times~{}10^{21}$ cm*-2* 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\arcmin~{}\times~{}2\arcmin$ ($7.8~{}\times~{}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$ $\mu$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_{H}=5.35\times 10^{-22}$ mag cm*-2* (Bohlin, Savage & Drake, 1978). Then, we used the typical extinction curve333The extinction curve used for the optical depth conversion can be found at https://www.astro.princeton.edu/~draine/dust/dustmix.html 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 $\mu$m and 89 $\mu$m. Hence, at each line of sight  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 Fig. 1, right and Fig. 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\arcsec$, outside the starburst ring) that is also true for our FIR polarimetry. Inside the starburst ring, where the optical depths are high (Fig. 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 & Whittet, 2015; Jones, 2015) provides an indirect measurement of the effect of the turbulent component. For maximally aligned dust grains along a line of sight with a constant magnetic field direction, the degree of polarization, $P$, in emission will be constant with optical depth $\tau$. If there is a region along the line of sight 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 line of sight with a well-defined scale length in $\tau$, 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 $\sigma_{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 anti-correlation 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 line-of-sights on the source, which allows 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 suffers of 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 towards dense regions and/or cross fields on the LOS.

The degree of polarization in our map is typically $<3$%, with a peak of $\sim 7$% (with $P/\sigma_{P}\sim 3$). The polarization is systematically higher in the outer regions than in the starburst disk at radius $\leq 25\arcsec$ ($1.6$ kpc). FIR polarimetric observations of Galactic clouds typically measure degree of polarization of $3-6$% in the $60-100$ $\mu$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\arcsec$ ($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 at those arising from the assembly and disruption of star-forming complexes, there will be multiple additional mechanisms for distorting pre-existing 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}\sim 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 $\mu$m beams could arise from starless or pre-stellar molecular cores (Caselli, 2013). Using 13CO 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}>10$ and temperatures of T$\sim 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 $\mu$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 cancelation 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 $20.6\pm 4.5^{\circ}$ at 0.46 $\mu$m (Berrier et al., 2013), $17.3\pm 2.2^{\circ}$ using R band images with the 2.5-m Las Cumbres Observatory (LCO) (Seigar, 2008), $15^{\circ}$ using the $H_{\alpha}$ velocity field (Emsellem et al., 2006), and $7-10^{\circ}$ using the CO (J=1-0) emission (Planesas, Scoville & Myers, 1991). We estimated a pitch angle of $16.9^{+2.7}_{-2.8}$∘ in the magnetic field inferred from 89 $\mu$m observations, which is between the pitch angles estimated in the spiral arms at 0.46 $\mu$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_{\alpha}$ velocity field (see fig. 7 by Emsellem et al., 2006). The analysis by Davies, Sugai & Ward (1998) computed that the $H_{\alpha}$ velocity field pattern is spatially correlated with the starbursting HII regions which was enhanced by a burst in the past 30 Myr. Later, Emsellem et al. (2006) found that the $H_{\alpha}$ 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, Roberts & Shu, 1975). As a consequence, recent studies have gone further and empirically shown that the spiral arm pitch angle, $\Psi$, and the central black hole mass, $M_{BH}$, may be correlated in the modal density wave (e.g. Seigar, 2008; Berrier et al., 2013) by the form

[TABLE]

(Berrier et al., 2013). Using our estimated pitch angle of $16.9^{+2.7}_{-2.8}$∘, we estimate a nuclear mass of $\log(M/M_{\odot})=7.16^{+0.46}_{-0.51}$. Our result is in agreement, within the uncertainties, with the black hole mass of $\log(M_{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 found several features between the measured PA of polarization and the magnetic field model within $\sim 1$ kpc of the nucleus. In the map of polarized flux (right-hand panel of Figure 1), we identify:

• •

a low polarization, $0.6\pm 0.3$%, within an $8\arcsec$ ($0.52$ kpc) diameter beam centered on the AGN,

• •

a strong peak in the polarized flux at PA $\sim 22^{\circ}$ at a radius of $\sim 4\arcsec$ ($\sim 0.26$ kpc) from the position of the AGN,

• •

three distinct minima at (a) at PA $\sim 22^{\circ}$, $8\arcsec~{}\leq~{}r~{}\leq~{}16\arcsec$ ($0.52-1$ kpc); (b) at PA $\sim 212^{\circ}$, $3\arcsec\leq r\leq 12\arcsec$ ($0.20-0.78$ kpc); and (c) at PA $\sim 240^{\circ}$, $14\arcsec\leq r\leq 18\arcsec$ ($0.91-1.17$ kpc), and

• •

significant discrepancies between our observed polarization directions and our spiral model.

The polarization of the nuclear region has been extensively studied from ultraviolet (UV) to MIR wavelengths. The UV-Optical polarization is very high, $~{}15$% (intrinsic polarization), and it is attributed to dust scattering from the NE regions of the NLR within the central kpc (e.g. Antonucci & Miller, 1985; Antonucci, Hurt, & Miller, 1994; Kishimoto, 1999). The NIR polarization, $\sim 7$% (intrinsic polarization), is attributed to dichroic absorption of aligned dust grains in the pc-scale dusty obscuring torus around the active nucleus (e.g. Bailey et al., 1988; Young et al., 1995; Lumsden et al., 1999; Packham et al., 1997; Lopez-Rodriguez et al., 2015; Gratadour et al., 2015). The MIR polarization, $<1$%, is attributed to self-absorbed dust emission from aligned dust grains in the pc-scale optically thick dusty torus (Packham et al., 2007; Lopez-Rodriguez et al., 2016). We attribute our measured FIR polarization of $0.6\pm 0.3$% at the core due to depolarization from 1) beam, $\sim 8$″, averaging of a more complex underlying field, 2) the high column density toward the core (Section 4.1), and/or 3) turbulent environment caused by the jet.

It is interesting to note the difference between the nuclear polarization of the radio-quiet AGN NGC 1068 from the recently reported highly polarized radio-loud AGN, Cygnus A (Lopez-Rodriguez et al., 2018b). At 89 $\mu$m, Cygnus A has a nuclear polarization of $\sim 10$% with a PA of polarization (B-vector) along the equatorial axis of the torus. Although both galaxies are at different distances and their angular diameters in the plane of the sky are different, we can estimate the polarization of NGC 1068 at the same spatial scale of Cygnus A. We measure the polarization of NGC 1068 at the same physical scale ($8$ kpc diameter) of Cygnus A, and we obtain a polarization level of $\sim 0.7$% with a PA of polarization (E-vector) of $\sim 144^{\circ}$, i.e. perpendicular to the overall extension of the polarized inner-bar (Fig. 1). The core of NGC 1068 at 89 $\mu$m seems to be affected by depolarization arising from the inner-bar and/or an intrinsically weak magnetic field, while Cygnus A may have an intrinsically polarized core. Further magnetohydrodynamical models are required to test the influence of magnetic fields in AGN tori to compare between radio-loud and radio-quiet objects.

The peak of the polarized flux is located $\sim 4$″ ($\sim 0.26$ kpc) NE from the AGN near the position of maximum dust temperature. This feature coincides with the northeastern radio lobe at 4.9 GHz (Wilson & Ulvestad, 1982), which is associated with high-velocity OIII-emitting clouds within a bow-shock that separates the AGN outflow from the local ISM. Young, Wilson & Shopbell (2001) observed a bright X-ray emission at the same position. The X-ray spectrum of this feature is a smooth continuum bremsstrahlung spectrum plus emission lines. Observations at $19.7-53$ $\mu$m show excess emission at $3\arcsec-6\arcsec$ ($0.20-0.39$ kpc) NE of the AGN after removal of an unresolved source at the position of the AGN (Lopez-Rodriguez et al., 2018a). Scarrott et al. (1991) also found the region to be polarized at a level of $3$%, which they attributed to the interaction of the jet with the ionization cone. At the angular resolution of the current FIR data, it is not possible to differentiate between the case in which the observed magnetic field at the position of peak polarized flux is generated by the bow shock or is intrinsic to the interstellar material with which the outflow is interacting.

The minima in polarized flux at $22^{\circ}$ (a) and $240^{\circ}$ (c) coincide with peaks in the $53$ $\mu$m map of Lopez-Rodriguez et al. (2018a), the locations of the most intense star formation in the starburst ring. The most likely explanation for the low polarization is that at these positions, the coherent large-scale field is completely dominated by randomly oriented fields generated during the collapse and destruction of a large number of star-forming cores.

The minimum in polarized flux at $212^{\circ}$ (b), in the direction of the counter-jet, may stem from a deficit in ISM density and CO intensity in that direction. The CO cloud density along the NE branch of the NIR bar and near the NE ionization cone appears to be greater than along the SW branch (see, e.g. Schinnerer et al., 2000; Tosaki et al., 2017).

The largest discrepancies between the measured polarization and the model of Section 3 (see Figure 4) lie at $5\arcsec\leq r\leq 15\arcsec$ at PA $\sim 70^{\circ}$ and PA $\sim 270^{\circ}$, in the directions of gaps in molecular cloud density along the starburst ring (Schinnerer et al., 2000). The gaps occur between the regions of most active star formation (see the $53$ $\mu$m map of Lopez-Rodriguez et al. (2018a)) and the points at which the outer spiral arms diverge from the starburst ring (PA $\sim 112^{\circ}$ and PA $\sim 292^{\circ}$). The deviations of the observed field directions from the model are in the sense of a larger radial component in the observed vectors. In the sector toward PA $\sim 70^{\circ}$, the polarized flux is particularly strong (see Fig. 1-right) and the vectors line up well with both the NE branch of the inner-bar (PA $\sim 45^{\circ}$) and visible dust lanes (see Fig. 2-right). If the magnetic field traces the gas flow, this would be consistent with inward transport of gas at the leading edges of the inner-bar, with the highest present-day flow rates occurring along the NE branch of the bar. Magnetic fields tracing inward gas transport have been observed in the inner-bar of NGC 1097 in radio synchrotron polarization (Beck et al., 2005).

5 Conclusions

We have presented the first FIR polarimetry of NGC 1068. We found a large-scale ($\sim 3$ kpc) spiral pattern that we attribute to thermal emission from magnetically-aligned grains. There is a spatial and morphological correspondence between the 89 $\mu$m magnetic field vectors and other tracers, i.e. OIII, 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.7}_{-2.8}$∘, inclination $48.1^{+1.8}_{-1.9}$∘, and tilt angle $37.8\pm 1.7^{\circ}$. 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 $\log(M/M_{\odot})=7.16^{+0.46}_{-0.51}$.

Outside the starburst disk (radius $\leq 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 position angles centered at $\sim 70^{\circ}$, 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 $\sim 4\arcsec$ ($\sim 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\pm 0.3$% within $5\arcsec$ ($0.33$ kpc) diameter at 53 $\mu$m and $0.6\pm 0.3$% within $8\arcsec$ ($0.52$ kpc) diameter at 89 $\mu$m. This is much lower than the $\sim 10$% FIR polarization of the radio-loud AGN in Cygnus A (Lopez-Rodriguez et al., 2018b). 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 M 82 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 that 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).