On the Origin of Radio-loudness in Active Galactic Nuclei Using Far-infrared Polarimetric Observations

Active galactic nuclei (AGNs) have been historically divided into radio-loud (RL) and radio-quiet (RQ) based on the relative ratio between their radio and optical luminosities (e.g., Stocke et al. 1992; Visnovsky et al. 1992; Kellermann et al. 1994). The most common definition is given by R = L_5GHz/L _B , where L_5GHz is the luminosity at 5 GHz and L _B is the luminosity in the B band. The typical cutoff is at R = 10. For an optically selected sample, quasars have been found to be mostly RQ, whereas ∼10%–20% are RL (e.g., Visnovsky et al. 1992). These studies typically use the integrated luminosities of the entire galaxy. Although this may be valid for quasars, where the jet outshines the host galaxy, it may not be fully valid for Seyfert galaxies. In addition, beaming effects can also affect both the optical and the radio radiation. Using a morphological classification, RL AGNs are mainly characterized by large-scale (tens of kiloparsecs) radio jets and lobes with their kinetic powers composing a significant fraction of the total bolometric luminosity. Small-scale jets (hundreds of parsecs) with thermal emission by dust and/or accretion disk along the equatorial plane are the dominant structures of RQ AGNs. Aside from the synchrotron component arising from the radio jet, the total flux spectral energy distributions (SEDs) from the infrared (IR) to the X-ray of RQ and RL AGNs are virtually identical (e.g., Sanders et al. 1989; Elvis et al. 1994; Hickox & Alexander 2018, for a review). Thus, the RQ and RL AGN dichotomy is thought to be intrinsically related to the radio jet production and/or stage of activity.

Advected magnetic fields (B-fields) are strongly implicated in all models that explain the launch and collimation of relativistic jets (Blandford et al. 2019, for a review). Within the central parsec, the presence of a radio source at the core and collimated jets imply a manifestation of the strong accretion flow onto the supermassive black hole (SMBH). Accretion flow toward an SMBH can be supported by small-scale (≪pc) B-fields generated by dynamo action and B-field lines stretching in a turbulent shear flow (e.g., Brandenburg et al. 1995; Sano et al. 2004; Beckwith et al. 2011; Simon et al. 2012). The transport of angular momentum is suggested to be operated by magnetorotational instability (MRI; Balbus & Hawley 1991, 1998), which gives rise to a sustained large-scale toroidal B-field and magnetohydrodynamical (MHD) turbulence. The MHD fluctuations within the disk produce a viscous stress responsible for the outward angular momentum transport driving gas accretion.

There are several physical mechanisms sustaining a large-scale B-field in AGNs. Let us first consider the infall of gas to form a disk. The natural outer radius for a disk is given by the radius of influence ${r}_{\inf }={{GM}}_{\mathrm{BH}}/{\sigma }^{2}\sim {10}^{6}$ times the gravitational radius, where the black hole mass is M_BH and the central velocity dispersion is σ. Within this radius, gravity is dominated by the black hole, and a thin disk is supposed to form. The thin disk allows gas to move inward under magnetic and gravitational torques. However, thin accretion disks have been found to be inefficient at depositing magnetic flux close to the SMBH. Thin disks make most of the AGNs RQ. In addition, thin accretion disks do not generate long-term large-scale poloidal B-fields that connect the accretion disk and the SMBH in RL AGNs (Lubow et al. 1994). Alternatively, if the infalling gas has a specific angular momentum >GM_BH/σ, it is likely to form a disk before it reaches ${r}_{\inf }$. The angular velocity will vary more slowly with increasing radius. Under this scenario, a toroidal B-field within a disk beyond the radius of influence could be generated.

Hot accretion in thick disks has been proposed to explain the accretion onto an SMBH (e.g., Ichimaru 1977; Rees et al. 1982; Lubow et al. 1994; Guilet & Ogilvie 2012). Geometrically thick accretion disks (i.e., magnetically elevated disks) are supported vertically against the vertical component of gravity by a strong toroidal B-field generated by the accretion disk dynamo (Begelman & Silk 2017). These accretion disks are thicker and less dense than standard disks, with stratified layers along the vertical axis. The magnetic pressure elevates most of the accretion flow to large heights from the midplane (Begelman & Silk 2017; Zhu & Stone 2018; Mishra et al. 2020). The midplane consists of a thin layer with most of the disk mass but with little contribution to the accretion rate. Under this scheme, the elevated layers can intercept a fraction of the incident radiation from the AGN and regulate the density and motion of the gas. These models may imply that there should be a steady loss of flux through Parker instabilities. These magnetically elevated disks may explain the accretion flow from parsec to subparsec scales, broad-line emission, and the relation between star formation and accretion in AGNs (Begelman & Silk 2017).

Another possibility is that an extended disk drives a large, unipolar MHD wind, which carries off mass, angular momentum, and energy from the disk (Blandford & Globus 2022). The B-field in such a wind will become increasingly toroidal and perhaps unstable. The infall is predominantly happening in the midplane of the disk at large scales (>parsec), while an MHD wind transports outward the energy from the SMBH within the Bondi radius. The MHD wind supports the jet collimation (see example in Figures 1(a) and (b) in Blandford & Globus 2022).

The fundamental difference between RL AGNs and their RQ AGN siblings may indeed be due to the presence or absence of a strong and ordered toroidal B-field. This ordered toroidal B-field enables the collimation of the jet and an effective transfer of angular momentum and energy outward from the SMBH. As the jet power depends on and/or generates the B-fields in the surrounding environment of the AGN, one may expect that the stronger the jet power is, the larger and stronger the ordered toroidal B-field is. This scenario invites the question of the origin and transportation of the B-fields from 1–100 pc scales to subparsec scales, which ultimately connects the host galaxy dynamics with the AGN. Within this framework, this manuscript aims to measure the signature of a magnetized flow at scales of <100 pc in a far-IR (FIR) flux-limited sample of RQ and RL AGNs.

The unified model of AGNs (Antonucci 1993; Urry & Padovani 1995) posits that a dusty and molecular structure (the "torus" ⁹ ) may surround some AGNs. The torus absorbs radiation at all wavelengths and reemits it at IR wavelengths. The postulated torus explains the large variety of AGN types from an orientation point of view. At parsec scales, the gas becomes dusty when the temperature is ≤1300 K (sublimation temperature of graphite), representing the inner wall of a dusty environment surrounding the active nuclei. Seyfert galaxies have been shown to have a dusty and molecular, geometrically thick obscuring disk along the AGN equatorial plane up to a median radius of ∼42 pc with a median molecular gas mass of ∼6 × 10⁵ M_⊙ (García-Burillo et al. 2021). This disk is rotating and resembles the inflowing gas in an accretion disk toward the SMBH with an outflowing component above and below the disk. These results were obtained using continuum emission at 870 μm and CO (3–2) and HCO⁺ (4–3) molecular emission lines with the Atacama Large Millimeter/submillimeter Array (ALMA). Most of the mid-IR (MIR) emission has been found to arise from the thermal emission of dust in parsec-scale outflows above and below the equatorial axis of the torus (e.g., Hönig et al. 2012, 2013; López-Gonzaga et al. 2016; Gámez Rosas et al. 2022). Thus, the dusty and molecular torus comprises a region where outflows and inflows can coexist (e.g., Hönig 2019).

The magnetic signature in the dusty environment surrounding AGNs can be measured by means of magnetically aligned dust grains. The elongated and paramagnetic dust grains can be aligned by the presence of B-fields described by theories of radiative alignment torques (RATs) and also by intense radiation fields or outflowing media (e.g., Andersson et al. 2015; Hoang & Lazarian 2016). Specifically, dust grains absorb more radiation along one of their axes, making them spin up along their greatest moment of inertia (i.e., short axis). Then, the dust grains acquire a magnetic moment, making them precess along the orientation of the local B-field. The final configuration is one where the short axis of the dust grains is parallel to the local B-field. As radiation propagates through these magnetically aligned dust grains, preferential extinction of radiation along one plane leads to a measurable polarization in the transmission of this radiation, a process called dichroic absorption. For dichroic absorption, the observed position angle (PA) of polarization traces the orientation of the B-field at optical and near-IR wavelengths. The same magnetically aligned dust grains reemit radiation at MIR and FIR (12–500 μm) wavelengths preferentially along their major axis, a process called dichroic emission. For dichroic emission, the observed PA of polarization is perpendicular to the local B-field orientation (e.g., Hildebrand 1988; Aitken et al. 2004).

Recent works have shown that there may be an intrinsic difference in the nuclear thermal polarization between RQ and RL AGNs. The general trends are as follows: (1) some RL AGNs have high intrinsically polarized (∼10%) cores at 2–90 μm dominated by magnetically aligned dust grains in the dusty torus at scales of 1–100 pc (Lopez-Rodriguez et al. 2013, 2015, 2018b; Lopez-Rodriguez 2021), and (2) RQ AGNs have low polarized (≤1%) cores within the 10–90 μm wavelength range (Lopez-Rodriguez et al. 2016, 2018a, 2020b), with their cores showing a variety of physical mechanisms (i.e., dust scattering, self-absorbing dust emission) at scales of approximately tens of parsecs. Thus, the FIR polarimetric results may indicate that the dusty torus represents an interface of a magnetized flow at 1–100 pc scales and a potential clue for the difference between RQ and RL AGNs. Note that some RL AGNs also have highly polarized cores at 2 μm (e.g., Capetti et al. 2000; Tadhunter et al. 2000; Ramírez et al. 2014), where the nuclear polarization arises mainly from scattering off dust located in the surrounding regions of the AGN. These studies analyzed the unresolved nuclear polarization using a single wavelength, which may be affected by several competing polarization mechanisms (i.e., scattering, dichroism, and synchrotron). A multiwavelength polarimetric analysis is required to disentangle the physical structure associated with the dominant polarization mechanism across the SED (Lopez-Rodriguez et al. 2018a; Marin 2018; Marin et al. 2020).

In this paper, we measure the 53–214 μm thermal polarized emission of magnetically aligned dust grains in the dusty environment surrounding RL and RQ AGNs. The results of this work are discussed in the context of a magnetized flow connecting the AGN with the host galaxy and the SMBH. We describe the specifics of our observations in Section 2. Section 3 shows the analysis of the total and polarized SEDs of our AGN sample. The radio-loudness classification and comparison with the observed thermal polarized emission are presented in Section 4. Our discussions are described in Section 5, and our main conclusions are summarized in Section 6. In this paper, we take H₀ = 67.8 km s⁻¹ Mpc⁻¹, Ω_M = 0.308, and Ω_V = 0.692.

To support the analysis of the origin of the polarized thermal emission in our AGN sample, we compile the total nuclear SEDs from 5 μm to 5 GHz. The SEDs have comparable, or better, angular resolution to that from our HAWC+ observations (Table 3). The angular resolution is high for all observations to separate the AGN from extended emission and/or radio lobes associated with the host galaxies. In the 10–30 μm wavelength range, the Spitzer spectroscopic data from the Combined Atlas of Sources with Spitzer IRS Spectra (CASSIS) archive ¹⁰ were used. We took the low-resolution and optimal extraction spectra (Lebouteiller et al. 2011). In the 70–500 μm wavelength range, we took PACS and SPIRE/Herschel observations from the Herschel archive. ¹¹ The FORCAST and HAWC+ observations from SOFIA were used to complement the 30–214 μm wavelength range. The total flux SOFIA observations of these data sets are described in Fuller et al. (2016, 2019) and L. Fuller et al. (2023, in preparation). New total and polarimetric data associated with HAWC+ are described in Section 2. In the 30–500 μm wavelength range, we performed aperture photometry of the SOFIA and Herschel observations using a diameter equal to the beam size at each wavelength as described in Section 2. References for the nuclear fluxes at 5 GHz are shown in Table 2.

Figure 1 shows the compiled 5 μm–5 GHz nuclear total SEDs of the AGN sample. All targets in the sample have similar 10–1000 μm SED behavior. All total SEDs show a bump at ∼100 μm arising from thermal dust emission from the core (Sections 3.1.1 and 3.1.2). This IR bump has been identified as arising from dust emission with a steep decrease toward submillimeter and near-IR (NIR) wavelengths (e.g., Barvainis et al. 1992). The main difference between the continuum total SEDs resides in the radio emission from their jets. This result is in agreement with the classical definition of RL and RQ AGNs and their total SEDs (e.g., Antonucci 1993; Elvis et al. 1994; Urry & Padovani 1995). Thus, the 1–1000 μm total SED range does not necessarily provide further information about the nature of RL and RQ AGNs. We find that the opposite is true in the FIR polarized flux (Section 3.3).

3.1.1. Thermal Emission Component

We characterize the observed IR bump of the total SEDs shown in Figure 1. The dust properties from 20 to 500 μm are estimated using a single-temperature modified blackbody function expressed as

$$\begin{eqnarray}&&{F}_{\lambda }=\frac{{M}_{d}}{{D}_{L}^{2}}{k}_{{\lambda }_{0}}{\left(\frac{\lambda }{{\lambda }_{0}}\right)}^{-\beta }{B}_{\lambda }({T}_{d}),\end{eqnarray} \tag{ 1 }$$

where M _d is the dust mass, D _L is the luminosity distance to the source, ${k}_{{\lambda }_{0}}$ is the dust mass absorption coefficient, which is taken to be 0.29 m² kg⁻¹ at a wavelength of λ₀ = 250 μm (Wiebe et al. 2009), β is the dust emissive index, and B_λ (T _d ) is the blackbody function at a characteristic dust temperature of T_d .

We have three free model parameters, M _d , β, and T _d , that we fit within the 20–500 μm wavelength range. We compute a Markov Chain Monte Carlo (MCMC) 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 range of M _d = [10³, 10¹¹] M_⊙, β = [1, 2], and T _d = [10, 100] K. We run the code using five chains with 5000 steps and a 1000 burn-in per chain, which provides 25000 steps for the full MCMC code useful for data analysis. Final median values and 1σ uncertainties are shown in Table 4, and the best-fit model with the associated 1σ uncertainty is plotted for each source in Figure 1.

Table 4. Best-fit Parameters of the Modified Blackbody Function (Equation (1))

Object	Md	T d	β	${\mathrm{log}}_{10}({L}_{\mathrm{TIR}})$ b	${\mathrm{log}}_{10}({L}_{89\mu {\rm{m}}})$	rd,tk	rd,th	θ89 μm
	(×105 M⊙)	(K)		(erg s−1)	(erg s−1)	(pc)	(pc)	(pc)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)
Centaurus A	${\,1.0}_{-0.1}^{+0.1}$	${22.2}_{-1.1}^{+1.0}$	${1.62}_{-0.18}^{+0.18}$	${43.38}_{-0.05}^{+0.06}$	${41.09}_{-0.12}^{+0.14}$	${4.7}_{-0.4}^{+0.4}$	${457}_{-167}^{+217}$	130
Circinus	${11.0}_{-1.4}^{+1.2}$	${37.7}_{-2.2}^{+2.1}$	${1.16}_{-0.14}^{+0.15}$	${44.22}_{-0.09}^{+0.10}$	${42.15}_{-0.11}^{+0.11}$	${5.3}_{-0.5}^{+0.6}$	${185}_{-70}^{+43}$	150
Cygnus A	${557.6}_{-16.4}^{+15.6}$	${47.7}_{-0.9}^{+0.9}$	${1.01}_{-0.01}^{+0.02}$	${45.79}_{-0.03}^{+0.03}$	${43.69}_{-0.03}^{+0.03}$	${20.3}_{-1.8}^{+0.4}$	${417}_{-37}^{+36}$	8880
Mrk 231	${1753.3}_{-331.9}^{+340.3}$	${38.0}_{-3.8}^{+3.8}$	${1.70}_{-0.21}^{+0.23}$	${46.57}_{-0.16}^{+0.18}$	${44.53}_{-0.17}^{+0.17}$	${78.6}_{-13.4}^{+19.7}$	${5802}_{-1699}^{+6773}$	6940
NGC 1068	${52.9}_{-2.0}^{+1.8}$	${28.3}_{-0.5}^{+0.6}$	${1.52}_{-0.05}^{+0.05}$	${45.07}_{-0.03}^{+0.03}$	${42.96}_{-0.04}^{+0.05}$	${23.8}_{-0.8}^{+1.8}$	${1736}_{-158}^{+213}$	560
NGC 1275	${56.6}_{-2.7}^{+2.7}$	${30.7}_{-0.9}^{+1.1}$	${1.02}_{-0.02}^{+0.04}$	${44.99}_{-0.03}^{+0.02}$	${42.84}_{-0.04}^{+0.04}$	${18.3}_{-0.8}^{+0.9}$	${440}_{-7}^{+62}$	2860
NGC 4151 a	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯
NGC 5506	${16.7}_{-1.7}^{+1.4}$	${39.4}_{-2.1}^{+1.9}$	${1.06}_{-0.06}^{+0.10}$	${44.37}_{-0.07}^{+0.05}$	${42.29}_{-0.07}^{+0.06}$	${5.8}_{-0.6}^{+0.6}$	${146}_{-35}^{+19}$	1020
PKS 1345+125	${992.7}_{-232.7}^{+201.1}$	${32.0}_{-4.6}^{+3.9}$	${1.26}_{-0.26}^{+0.38}$	${46.29}_{-0.21}^{+0.16}$	${44.19}_{-0.32}^{+0.21}$	${63.4}_{-12.8}^{+47.4}$	${3346}_{-1784}^{+3967}$	18190

Notes. Column (1): object name. Column (2): dust mass. Column (3): characteristic dust temperature. Column (4): dust emissivity index. Column (5): total IR luminosity. Column (6): 89 μm luminosity. Column (7): the radius of the optically thick dusty environment assuming radiative equilibrium using Equation (2). Column (8): the radius of the optically thin dusty environment assuming radiative equilibrium using Equation (3). Column (9): spatial scale of the beam size of the observations at 89 μm (Table 3).

^a A single modified blackbody within the 30–500 μm wavelength range for NGC 4151 does not reproduce the total flux SED. ^b Thermal IR luminosity in the 8–1000 μm wavelength range.

Download table as:  ASCIITypeset image

We find that the IR bump peaks within the wavelength range of [77, 142] μm, with a median of λ _d = 100 ± 20 μm. The characteristic dust temperature of the IR bump is in the range of [22.2, 47.7] K, with a median of T _d = 34.8 ± 7.8 K. The dust emissive index is in the β = [1.01, 1.70] range, with a median of β = 1.21 ± 0.26. This result is in agreement with the typically assumed value of 1.5 for galaxies in the range of 1–2 (Hildebrand 1983). The associated dust mass of the IR bump has a wide range of 1.0 × 10⁵ to 1.8 × 10⁸ M_⊙. This is mainly due to the volume associated with the measured core, given the fixed angular resolution of our observations. The beam sizes at 89 μm vary from 130 pc in Centaurus A to 18.19 kpc in PKS 1345+125. In addition, the differences may be due to AGN types, for example, Mrk 231 is an ultraluminous IR galaxy with known dust masses up to ∼10⁹ M_⊙ (Lisenfeld et al.2000). We estimate the thermal IR luminosity, L_TIR, in the range of 8–1000 μm to be in the range of ${\mathrm{log}}_{10}({L}_{\mathrm{TIR}}$ [erg s⁻¹]) = [43.38, 46.57].

We compute the effective size of the dust emitting at 89 μm assuming radiative equilibrium at a distance r _d from the central engine as

$$\begin{eqnarray}&&{r}_{d,\mathrm{tk}}=441{\left(\frac{30\,{\rm{K}}}{{T}_{d}}\right)}^{2}{\left(\frac{{L}_{89\mu {\rm{m}}}}{{10}^{12}\,{L}_{\odot }}\right)}^{1/2}\,\mathrm{pc}\end{eqnarray} \tag{ 2 }$$

(Scoville 2013) for a normalized dust temperature of 30 K and a luminosity at 89 μm normalized at 10¹² L_⊙ (Table 4). Note that this relation assumes that the emission arises from optically thick dust. Our best fits suggest that a modified blackbody function best describes the IR bump, which indicates that the dust is optically thin at FIR wavelengths (Section 5.1). The effective source size from optical thin dust can be estimated as

$$\begin{eqnarray}&&{r}_{d,\mathrm{th}}={\left(\frac{1300\,{\rm{K}}}{{T}_{d}}\right)}^{\tfrac{4+\beta }{2}}{\left(\frac{{L}_{89\mu {\rm{m}}}}{{10}^{46}\,\mathrm{erg}\,{{\rm{s}}}^{-1}}\right)}^{1/2}\end{eqnarray} \tag{ 3 }$$

(e.g., Netzer 2008; Baron & Netzer 2019). We show the radii arising from a structure with optically thick and optically thin dust in Table 4. For optically thick dust, r_d,tk is in the range of [5, 79] pc, and for optically thin dust, r_d,th is in the range of [146, 5802] pc. We discuss these results in Section 5.3.

A single-temperature modified blackbody source can explain the IR bump of all galaxies, except for NGC 4151, whose IR total SED is flatter than that of the other galaxies in our sample. NGC 4151 is a type 1.5 AGN with an extended (∼35, ∼255.5 pc) dusty structure cospatial with the narrow-line region. This structure was observed in total intensity at 10 μm (broadband filter N) with Gemini-North/OSCIR (Radomski et al. 2003). Follow-up imaging polarimetric and spectropolarimetric observations with narrowband filters at 8.7, 10.3, and 11.6 μm using CanariCam on the 10.4 m Gran Telescopio CANARIAS (GTC) measured a similar structure in total flux. This extended structure is polarized at a level of 0.5%–1% potentially arising from dust scattering (Lopez-Rodriguez et al. 2018a). Given the complexity of the extended emission within our beam size (8'', 580 pc), several blackbodies may be required to explain the IR total SED of NGC 4151. As we are interested in the AGN, a detailed model of NGC 4151 is outside the scope of our goals.

3.1.2. Nonthermal Emission Component

We follow a similar approach from Lopez-Rodriguez et al. (2018b) to characterize the nonthermal emission from the total SED of AGNs. We assume a self-absorbed synchrotron component with a flat spectrum, ${\nu }^{{\alpha }_{\mathrm{tk}}}$, in the optically thick region and a turnover at a frequency, ν _c , with a flat spectrum, ${\nu }^{{\alpha }_{\mathrm{th}}}$, in the optically thin region. The synchrotron emission is described as

$$\begin{eqnarray}&&{F}_{\mathrm{syn}}=C{\left(\frac{\nu }{{\nu }_{c}}\right)}^{{\alpha }_{\mathrm{tk}}}{e}^{{\left(-\tfrac{\nu }{{\nu }_{c}}\right)}^{{\alpha }_{\mathrm{th}}}},\end{eqnarray} \tag{ 4 }$$

where C is a scale factor.

Our goal is to show the overall tendency of the synchrotron component rather than a fit to the SEDs. The total flux F_syn is scaled so that it is equal to the total flux at 5 GHz, and the polarized flux PI_syn is equal to the upper limits at radio wavelengths (Section 3.3). We adopt the same spectrum α_tk = + 0.21 and α_th = − 0.9 from the best fit of Cygnus A (Lopez-Rodriguez et al. 2018b). The turnover frequency is adopted to be 543 μm from the best fit of Cygnus A (Lopez-Rodriguez et al. 2018b) and Centaurus A and 2500 μm for PKS 1345+125. These turnover frequencies ensure that the synchrotron emission is always equal to or smaller than the total flux at 500 μm, whose emission may arise from a contribution of both thermal and nonthermal components. For Centaurus A, the break frequency is required to satisfy the upper limit in polarized SED (brown open square). A straight line will produce a higher polarized flux at 850 μm than the measurement. Note that even a single power law toward short wavelengths would produce negligible synchrotron emission at 100 μm.

The total SED model is computed as the linear combination of the modified blackbody function (Equation (1)) and the synchrotron component (Equation (4)), i.e., F _T = F_λ + F_syn. As shown in Figure 1, the synchrotron emission is negligible at ∼100 μm.

Our FIR polarimetric measurements may be diluted by the emission and/or absorption from the dust lane in the host galaxy. For extended sources, we quantify the polarized emission arising from the host galaxy and subtract it from the measured polarization to compute the intrinsic polarization of the nuclear source. Here we show the computation for Centaurus A, Circinus, and NGC 1068. The rest of the galaxies are point-like sources, and the polarized emission from aligned dust grains in the host galaxy cannot be estimated. However, the integrated FIR polarization contribution from the host galaxy is considered to be negligible (Section 5.1). Thus, we take their measured polarization, P _m , shown in Table 3, as the intrinsic polarization.

Centaurus A.—Lopez-Rodriguez (2021) found a 89 μm polarized nuclear source with a measured polarization of P _m = 1.5% ± 0.2% at a PA_E = 61° ± 4° (PA_B = 151° ± 4°). The PA_B is the B-field orientation after a rotation of 90° of the E-vector, i.e., PA_B = PA_E + 90° (Section 5.1). The AGN is embedded in a warped molecular disk with a measured polarization of P_gal ∼ 3.5% and a PA in the range of PA_gal,E = [15°, 57°] (PA_gal,B = [105°, 147°]). Thus, the measured nuclear polarized emission may be the contribution from the AGN and the host galaxy.

The nuclear total flux can be decomposed by emission from the AGN and host galaxy. To obtain the fractional contribution to the total emission from both unresolved (AGN) and extended (host galaxy) components, we follow a similar procedure applied to the HAWC+/SOFIA observations of NGC 1068 by Lopez-Rodriguez et al. (2018c). First, the fluxes in a circular aperture equal to the FWHM at 53 and 89 μm were measured. This flux represents the total flux, F _T , of the nuclear resolved source at the given wavelength and minimizes the contribution from the extended diffuse emission (Table 3). Second, the central 3 × FWHM radius emission was fitted with a composite model using the corresponding HAWC+ PSF to each observation and a 2D Gaussian profile. The total flux from the scaled PSF, ${F}_{\mathrm{PSF}}^{m}$, represents the maximum likely contribution from an unresolved nuclear component at the given angular resolution of the observations. The total flux from the 2D Gaussian profile, ${F}_{\mathrm{gal}}^{m}$, provides the minimum contribution of the extended component surrounding the central source. The total flux model, ${F}_{T}^{m}$, is the sum of both PSF and 2D Gaussian profiles. This method has four free parameters: the amplitude of the PSF, C_PSF; the FWHM of the long (b) and short (a) axes; and the PA, θ, of the 2D Gaussian profile. The amplitude of the 2D Gaussian profile is estimated as the peak of the observed total intensity at each wavelength minus the amplitude of the PSF. This approach ensures that the final model has the same peak flux as the observations.

The fitting routine was performed using an MCMC 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 of ${C}_{\mathrm{PSF}}=[0,\max (\mathrm{Flux})]$, where the upper limit corresponds to the peak flux at 53 and 89 μm; FWHM_a,b = [2, 10] × FWHM_53,89, where FWHM_53,89 corresponds to the FWHM at 53 and 89 μm; and θ = [0°, 180°). Note that the 2D Gaussian profiles are set to be always larger than the PSF at a given wavelength. We run the code using 10 chains with 5000 steps and a 1000 burn-in per chain, which provides 50,000 steps for the full MCMC code useful for data analysis. For the best model, the nuclear fluxes of the PSF, 2D Gaussian, and final model within the beam size at each wavelength were estimated. Appendix A shows the best-inferred parameters and associated fluxes (Table 6), the posterior distributions of the free parameters (Figure 4), and the original and residual images of the PSF-fitting routine (Figure 5).

Our PSF-fitting approach shows that the AGN contributes 47% ± 5% and 32% ± 3% within the central 80 and 130 pc at 53 and 89 μm, respectively. We compute the polarization of the nuclear source corrected by host galaxy dilution as

$$\begin{eqnarray}&&{P}_{\mathrm{agn}}={P}_{m}\left(1+\frac{{F}_{\mathrm{gal}}}{{F}_{\mathrm{PSF}}}\right)\end{eqnarray} \tag{ 5 }$$

(Lopez-Rodriguez et al. 2013, 2015), where P _m is the measured polarization, F_gal is the extended flux of the host galaxy, and F_PSF is the flux of the unresolved core.

After P_agn is estimated, this value and the PA_obs are corrected by the vector quantity arising from the host galaxy as follows. The Stokes parameters q_gal and u_gal from the host galaxy (P_gal and PA_gal) are estimated. Then, they are subtracted from the Stokes parameters q_agn and u_agn, which are estimated using the corrected polarization as shown in Equation (5) and the observed PA as shown in Table 3. The final intrinsic Stokes parameters are estimated such as q_int = q_agn − q_gal and u_int = u_agn − u_gal. We estimate an intrinsic polarization of P_int = 4.0% ± 1.3% and 4.7% ± 0.7% with a PA_int,E = 45° ± 9° and 45° ± 4° at 53 and 89 μm, respectively. The B-field orientation is estimated to be PA_int,B = 135° ± 9° and 135° ± 4° at 53 and 89 μm, respectively.

Circinus.—Within the 53–214 μm wavelength range, the measured nuclear polarization is estimated to be ∼2% with a PA_E ∼ 135° (Lopez-Rodriguez et al. 2022a), where the B-field orientation, PA_B ∼ 45°, is similar to the orientation of the central 5 kpc bar at a PA of ∼45° (Jones et al. 1999). Based on the B-field model by L. Grosset et al. (2023, in preparation), the polarization arises from the spiral arms and starburst ring within the central 1 kpc. These authors found that the AGN at each wavelength is consistent with an unpolarized source. We take that the AGN is intrinsically unpolarized at FIR wavelengths and take an upper limit of 0.4% equal to the absolute minimum polarization measured by HAWC+ (Harper et al. 2018).

NGC 1068.—Although the 89 μm polarization along the 1 kpc bar is measured to be in the range of 1%–3% (Lopez-Rodriguez et al. 2018b), a minimum in the polarization is found at the nucleus with a PA_E ∼ 141°. The inferred B-field orientation, PA_B ∼ 39°, is similar to the orientation of the bar, PA ∼ 48° (Scoville et al. 1988). At 53 μm, we measure a nuclear polarization ∼1.2% at a PA_E ∼ 33° within the central beam (5'', 70 pc). Note that the PA_E of polarization rotates by ∼119° between 53 and 89 μm. This result can be easily observed in Figure 1 in Lopez-Rodriguez et al. (2022a). This angular rotation may suggest that the thermal polarization at 53 μm arises from absorptive polarization, which indicates that the measured PA of polarization is the orientation of the B-field in the bar. Indeed, the PA of polarization of 33° arising from absorptive polarization at 53 μm is in agreement with the inferred B-field of 39° arising from emissive polarization at 89 μm. Another plausible explanation is that the measured PA of polarization at 53 μm may arise from polarized thermal emission at the origin of the spiral arms coming out from the core. The angular resolution improvement from 560 to 350 pc at 89 and 53 μm respectively, may be able to explain the angular difference. The 8–12 μm imaging polarimetric observations using CanariCam on the 10.4 m GTC within the central ∼100 pc of NGC 1068 showed a highly polarized (∼4%–7%) region with a PA_E ∼ 44° at ∼24 pc north from the core (Lopez-Rodriguez et al. 2016). The polarization was found to arise from dust and gas emission in the ionization cone, heated by the active nucleus and jet, and further extinguished by aligned dust grains in the host galaxy. In addition, these observations found an unpolarized (∼0.4%) core. Using ALMA 860 μm polarimetric observations at a resolution of 007 (4.2 pc), Lopez-Rodriguez et al. (2020a) measured an unpolarized central core and a polarized dusty torus at a level of 3.7% ± 0.5% with a B-field orientation of 108° ± 2° within 5–9 pc east from the core. Given the low angular resolution (560 pc) of the 89 μm HAWC+ observations, we suggest that the nuclear FIR polarization arises from the polarized thermal emission of the aligned dust grains in the bar–rather than the emission from the AGN. We conclude that the AGN of NGC 1068 at 89 μm is intrinsically unpolarized at FIR wavelengths. We take an upper limit of 0.4% as the absolute minimum polarization measured by HAWC+ (Harper et al. 2018). The 53 μm polarimetric observations are inconclusive and will not be used in this work.

In summary, we find that at 89 μm Circinus, Mrk 231, NGC 1068, NGC 4151, and NGC 5506 are consistent with unpolarized sources. We find that Centaurus A, Cygnus A, PKS 1345+125 (4C +12.50), and potentially NGC 1275 (Section 5.4) are consistent with intrinsic polarized sources at a level of 1.9%–11%. Thus, the 89 μm polarization fraction shows differences between AGNs.

We compile the nuclear polarized SEDs from 2 μm to 5 GHz to support the analysis of the origin of the polarized thermal emission from AGNs. The polarized SEDs have angular resolution similar to or better than that from our HAWC+ observations. Appendix B (Table 7) shows the wavelengths, apertures, and polarization fluxes of the polarized measurements shown in Figure 1. Due to the small number of data points in the polarized SEDs, we aim to provide fiducial polarized SEDs to quantify the relative contribution of thermal and nonthermal polarized emission at FIR wavelengths.

For the polarized AGNs, we scaled the polarized thermal emission to be equal to the median polarized flux within the 53–89 μm wavelength range. This scale factor is equal to the median polarization fraction of 4.5% and 10% within 53–89 μm, i.e., 0.045 for Centaurus A and 0.1 for Cygnus A and PKS 1345+125. We scaled the polarized nonthermal emission to be equal to the polarized fluxes at millimeter wavelengths. These scale factors are 0.004, 0.001, and 0.002 for Centaurus A, Cygnus A, and PKS 1345+125, respectively. These AGNs have an intrinsic polarization fraction <0.4% at millimeter wavelengths arising from synchrotron emission. Figure 1 shows the total polarized SEDs (black dashed line) for the polarized AGNs.

For the unpolarized AGNs, we only show the polarized SED of NGC 1068. The reason is that this AGN is the only object that has the dusty torus and surrounded jet resolved within the central 20 pc at a resolution of 4.2 pc (007) at 860 μm (Lopez-Rodriguez et al. 2021b). These observations showed an unresolved core of 4.2 pc in diameter (i.e., FWHM) with an upper-limit polarization fraction of <0.2%. Due to the optically thick dusty torus, the central unpolarized core arises from a self-absorbed synchrotron emission. The dusty torus is resolved and polarized at 3.7% ± 0.5% at 860 μm across the 5–9 pc eastern region of the equational axis. Figure 3 shows the total flux of 12 ± 0.5 mJy (brown filled square), the polarized flux of the dusty torus of 0.44 ± 0.06 mJy, and the unpolarized core with an upper-limit polarized flux of <0.024 mJy. We use the same total SED for the polarized AGNs with a turnover wavelength of 2000 μm. We scaled the polarized flux as follows. The nonthermal emission is equal to the upper-limit polarized flux of the unpolarized core, and the thermal emission is equal to the polarized dusty torus. The scale factors for thermal and nonthermal components are 0.0002 and 0.008, respectively. This result indicates that thermal and nonthermal components are unpolarized at <0.02% at 100 μm.

For the polarized AGNs (i.e., Centaurus A, Cygnus A, and PKS 1345+125) and unpolarized AGNs (NGC 1068), we find that the polarized synchrotron emission is negligible at FIR wavelengths. Figure 1 shows that the nonthermal polarized emission (i.e., synchrotron component) is several orders of magnitude lower than the polarized thermal emission arising from magnetically aligned dust grains at FIR wavelengths. This result holds even if a single power law is used to interpolate the synchrotron emission to FIR wavelengths. As we will discuss in Section 5.4, NGC 1275 is excluded from this analysis because the origin of the measured polarization is inconclusive. For the polarized AGNs with several data points in the FIR (i.e., Centaurus A and Cygnus A), the polarized IR emission has the same shape as the IR bump of the total SED. The polarized thermal component contributes <0.4% of the total SEDs in our unpolarized AGNs (i.e., Circinus, Mrk 231, NGC 1068, NGC 5506).

We conclude that polarized thermal emission by means of magnetically aligned dust grains is the dominant polarization mechanism in the 53–214 μm wavelength range. By combining the total flux and polarized SEDs with the characterization of the IR bump, we can use the FIR polarization measurements to infer the B-field orientation of the dusty environment surrounding polarized AGNs. We show the first observational evidence of an intrinsic difference between AGNs using the 89 μm nuclear polarized SEDs.

In this section, we show that our observations trace the B-field orientation by means of the thermal emission of magnetically aligned dust grains. As a quantitative analysis, we estimate the expected thermal polarized emission based on measurements of the absorptive polarization from magnetically aligned dust grains. Centaurus A has a nuclear polarization of 1.97% ± 0.01% with a PA_E = 13409 ± 016 within an aperture of 8'' at 2.2 μm and a visual extinction of A _v = 16 mag (Packham et al. 1996). The optical depth at K band (2.0 μm) can be estimated such as τ _K = 0.09A _v = 1.44 (Jones 1989). Using the typical extinction curve ¹³ of the Milky Way with an extinction factor of R _v = 3.1 (Weingartner & Draine 2001), we estimate the optical depth at 89 μm to be τ₈₉ = 1.18 × 10⁻² τ _K = 1.69 × 10⁻². This result indicates that the dusty structure associated with the visual extinction to the core is optically thin at 89 μm. Under the optically thin condition, the polarized thermal emission can be estimated as ${P}_{89}^{\mathrm{em}}=-{P}_{89}^{\mathrm{abs}}/{\tau }_{89}$. The negative sign indicates the 90° change from absorptive to emissive polarization (Hildebrand 1988; Aitken et al. 2004). The 90° angular change is satisfied from 2.2 μm, PA_E = 13409 ± 016, to 89 μm, PA_E = 45° ± 4°, where we estimate an angular difference of ∣Δθ∣ = 89° ± 4°. As the extinction curve is representative of the absorptive polarization, the scaled measured 2.2 μm polarization of 1.97% ± 0.01% at 89 μm is estimated to be ${P}_{89}^{\mathrm{abs}}\sim 0.02 \%$. Finally, we estimate an expected emissive polarization of ${P}_{89}^{\mathrm{em}}\sim 1.18 \%$ at 89 μm. This result is in excellent agreement with the measured polarization of 1.5% ± 0.2% at 89 μm (Table 3). Similar results were also found for NGC 1068 (Lopez-Rodriguez et al. 2020b) and NGC 1097 (Lopez-Rodriguez et al. 2021a) and using a detailed total and polarized SED fitting for Cygnus A (Lopez-Rodriguez et al. 2018b). In addition, we show in Section 3.3 that the polarized synchrotron emission is negligible at FIR wavelengths for the polarized AGN. Thus, we conclude that the 89 μm polarization measurements of AGNs arise from the thermal emission of magnetically aligned dust grains. The measured PA_E (E-vector in Table 3) of polarization can be rotated by 90° to show the inferred B-field orientation.

Because the beam sizes of our 89 μm observations are >130 pc (Table 4) and we mostly have a single wavelength observation per object (Table 3), it is very challenging to analyze the dust grain properties. We use the resolved, ∼5 pc, observations of 30 Doradus, a star-forming region in the Large Magellanic Cloud, as an example of thermal polarization emission properties in an external galaxy. Using 89 μm polarimetric observations with HAWC+, Tram et al. (2021) measured polarization fractions up to 15% in 30 Doradus. The highest polarization arises from dust grains located at the lowest column density with dust temperatures ≤37 K around the star-forming regions. These authors argued that this polarization could be explained under the RAT (radiative alignment torques) theory (Andersson et al. 2015; Hoang & Lazarian 2016), where the polarization fraction depends on the dust grain alignment size. Specifically, the dust grain alignment size depends on the gas volume density, gas temperature, and radiation field in the interstellar medium (ISM). In 30 Doradus, typical grain sizes of 0.01–0.82 μm could explain the 15% polarization fraction. Despite the physical differences of 30 Doradus and at ≤130 pc around AGNs, these results show that high levels of polarized thermal emission can be observed under extreme conditions. In addition, these levels of polarization can be explained with nominal dust grain sizes and compositions. Note that for higher polarization fractions, such as those up to 25% measured by Planck (Planck Collaboration et al. 2015, 2020), other dust grain compositions (e.g., "astrodust") may be required (Draine & Hensley 2021; Hensley & Draine 2021).

Based on the results of emissive polarization in galaxies using FIR wavelengths, we conclude that typical dust grain sizes and composition in the ISM can explain our measured high polarization fractions of AGNs. Note that "core" refers to spatial scales smaller than the beam size of our observations, i.e., from <130 pc to <18.19 kpc. Our results suggest that the FIR polarimetric observation may trace an ordered B-field within a dusty structure with low dust temperature and low surface brightness around the active nuclei.

We estimated the effective radius of the dusty structure emitting at 89 μm assuming radiative equilibrium to be r_d,tk = [5, 79] pc and r_d,th = [146, 5802] pc for optically thick and optically thin dust, respectively (Section 3.1.1). We interpret r_d,tk as a lower-limit radius of the dusty structure surrounding the AGN. The r_d,th is substantially larger than r_d,tk because the AGN radiation can travel to larger distances through the dusty medium without being totally absorbed. Note that r_d,th is larger than or similar to the physical scale of the beam size of our observations at 89 μm for most of our AGNs. However, the full range of distances up to the maximum of r_d,th cannot produce the polarization levels measured in RL AGNs. We discuss the origin of the highly polarized dusty regions in RL AGNs and the estimations of the effective radius.

Based on FIR polarimetric observations of resolved galaxies, spiral galaxies have a median 154 μm polarization fraction of 3.3% ± 0.9% (Lopez-Rodriguez et al. 2022b). The lowest polarization arises from regions with high turbulence and/or tangled B-fields in the ISM owing to galaxy interaction and star formation activity. These physical conditions generate an isotropic turbulent B-field and/or tangled B-fields along the LOS, resulting in a net null polarization within the beam (typically ∼300 pc diameter) of the observations. The highest polarization mainly arises from low-density, low surface brightness, low column density, low dust temperature, and low gas turbulence regions of the ISM. These conditions increase the polarization efficiency of dust grains with the local B-field in the host galaxies without invoking other physical conditions than typical dust parameters in the ISM. For both face-on and edge-on spiral galaxies without an (or a weak) AGN, the cores (beam size ≤1 kpc) are measured to be unpolarized at FIR wavelengths (Borlaff et al. 2021; Lopez-Rodriguez et al. 2022b). The median 89 μm polarization fraction across the resolved ∼3 kpc dust lane of Centaurus A is ∼4% (Lopez-Rodriguez 2021). These results show individual polarization measurements with spatial resolutions from 90 pc to 1 kpc across the host galaxies. Note that RL AGNs are typically hosted by early-type galaxies (Kauffmann et al. 2003), which may have different extended properties from those of spiral galaxies compared above. In addition, inactive galaxies have unpolarized cores at FIR wavelengths. These results suggest that, most likely, the high polarization in RL AGNs may not arise from physical structures associated with the host galaxy below the beam size of our observations.

The integrated FIR polarization of spiral galaxies is <1% (Table 4 and Figure 3 by Lopez-Rodriguez et al. 2022b). This estimation was performed for the entire host galaxy covering spatial scales of 2–10 kpc. This result rules out that the high polarization in unresolved RL AGNs (i.e., Cygnus A, PKS 1345+12) may be attributed to the integrated polarization arising from the entire host galaxy.

Knowing that the entire galaxy and physical structures within the host galaxy cannot produce high polarization, and using Centaurus A as the example owing to its proximity, we suggest that the polarized source in the core of RL AGNs should arise from a structure r_d,th > 5 pc and <130 pc in radius directly radiated by the AGN.

The estimations of r_d,tk and r_d,th assume that the RL AGN solely heats the polarized dust. In reality, the total flux at 89 μm may arise from a combination of dust heated by the AGN, dust heated by star formation around the AGN (i.e., starburst ring), and/or dust in the ISM of the host galaxy. However, not all these physical components can produce the high levels of FIR polarization measured in RL AGNs. As the RQ AGNs are unpolarized, we cannot constrain the physical sizes of the polarized structure, and the assumption that the RQ AGN solely heats the polarized dust may not apply to these objects. We explore these scenarios to explain the 5%–11% polarization at 89 μm in RL AGNs.

Tadhunter et al. (2007) measured that the AGN radiation, traced by [O iii] λ5007 emission, is the primary heating mechanism at 24 and 70 μm. The scatter measured in [O iii] and 70 μm may arise from starburst activity heating the cool dust at the lower end of luminosity ranges, ${\mathrm{log}}_{10}({L}_{70\mu {\rm{m}}}[\mathrm{erg}\,{{\rm{s}}}^{-1}])\lt 42.6$, but only for 20%–30% of objects. Star formation activity decreases the FIR polarization owing to tangled B-fields and/or increased turbulence in the ISM. Galactic shocks, as measured in the intersection of the galactic bars and starburst rings of NGC 1097 (Lopez-Rodriguez et al. 2021a), can compress the B-fields in the dense and cold ISM. Although a constant B-field of polarization parallel to the starburst ring would be measured for an edge-on view of the starburst ring, the 89 μm polarization is measured to be ≤3%.

We estimated that the polarized synchrotron emission is negligible at 89 μm (Sections 3.1.2 and 3.3) and that the B-field orientation is highly offset from the radio jet axis in RL AGNs (Section 5.4). The radio jets in our RL AGN sample are mostly perpendicular to the galaxy disks. The radio jet in NGC 1068 interacts with a giant molecular cloud at ∼20 pc north from the AGN producing high levels of polarization: ∼7% at 10 μm (Lopez-Rodriguez et al. 2016), 5% at 860 μm (Lopez-Rodriguez et al. 2020a), and up to 11% at 860 μm in several knots within the central 100 pc. However, the polarization mechanisms differ from those at 89 μm. At 10 μm, the polarization arises from dichroic absorption of a polarized region through the dust lane, while synchrotron polarization dominates at 860 μm. In addition, the 89 μm polarization of the core and the bar of NGC 1068 is <2% (Lopez-Rodriguez et al. 2020b). All these results rule out the possibility that the polarization in RL AGNs arises from dust heated by jet interaction in the ISM and star formation activity.

For the polarized objects of our sample, we compute the angular difference between the measured B-field orientation at 89 μm and the nuclear jet axis, i.e., jet–disk (Table 5). Figure 3 shows our measured B-field orientation at 89 μm and the 5 GHz radio jet orientation over the WFPC2/HST images. We find that 75% (three out of four objects) of the polarized cores have a jet–disk angular difference in the range of Δθ = 64°–84°, with a median of 〈Δθ〉 = 65° ± 26°. The exception is NGC 1275 with an angular difference of 19° ± 2°. Despite the small sample, we find several interesting results.

Zoom In Zoom Out Reset image size

Table 5. Jet Axis and Magnetic Field Orientation at 89 μm

Object	${\mathrm{PA}}_{\mathrm{jet}}^{1}$	PAB	Δθ	References
	(deg)	(deg)	(deg)
Centaurus A	51	135 ± 4	84 ± 4	Israel (1998)
Cygnus A	105	39 ± 7	66 ± 7	Sorathia et al. (1996)
NGC 1275	160	141 ± 2	19 ± 2	Pedlar et al. (1990)
PKS 1345+125	157	41 ± 8	64 ± 5	Lister et al. (2003)

Download table as:  ASCIITypeset image

The inferred B-field orientation is highly offset with respect to the jet axis, but not exactly at 90°. Centaurus A is known to have a warped disk from 2 to 6500 pc (Quillen et al. 2010). Using our estimated beam size of 8'' (radius of 4'' = 64 pc) at 89 μm, we find that the PA of the velocity component is 1248 at a radius of 399 from the core (Table 1 in Quillen et al. 2010). The associated ring of matter at a radius of 64 pc around the center of Centaurus A is almost parallel, ∣135°–1248∣ = 102, to the measured B-field orientation. This result indicates that the measured B-field may be spatially associated with the gas flows at 64 pc from the core of Centaurus A. Using 2 μm polarimetric observation with NICMOS/HST, Tadhunter et al. (2000) found that the ionization cone of Cygnus A is only polarized along one of the edges up to 1.3 kpc from the core. The 2 μm polarization is thought to arise from dust scattering in the ionization cone. These authors argued that the observed polarization is due to intrinsic anisotropic radiation caused by a warped disk. A disk of 80–100 pc in radius roughly perpendicular to the radio jet axis was observed using the Very Long Baseline Array (VLBA) in HI absorption (Struve & Conway 2010). These results indicate that the measured B-field may be spatially associated with the warped disk at 24.7 pc from the core of Cygnus A. The radio jet of PKS 1345+125 is thought to be precessing with a cone half-angle of 23° (Lister et al. 2003). Indeed, misaligned disks have been measured in the nuclear structures of AGNs (e.g., Lawrence & Elvis 2010).

These results show that the most plausible interpretation of our measurements is that the measured B-field orientation at 89 μm may be associated with a toroidal B-field along the equatorial axis of the AGN. The angular difference between the equational and jet axes, Δθ, can be explained by intrinsic anisotropies of the magnetized structure and jet axis (i.e., warped disk, precessing jets, MHD wind).

The exception is NGC 1275, where the inferred B-field orientation seems almost parallel, Δθ = 19° ± 2°, to the jet axis. Our 89 μm dust continuum emission observations show extended emission at an 11σ detection with a PA ∼125° out to a radius of ∼12 kpc (Appendix C, Figure 6). NGC 1275 shows thermal emission features extended at PAs of ∼120° and ∼75° using 70–160 μm Herschel observations with beam sizes of 5''–12'' (Mittal et al. 2012). The thermal emission at 450 μm with a beam size of 91 × 88 using the James Clerk Maxwell Telescope (JCMT) also shows an extended emission at ∼115° up to a radius of ∼20 kpc (Irwin et al. 2001). These observations indicate that a large amount of dust is present in this galaxy and the surrounded Perseus cluster (∼6 × 10⁷ M_⊙ with T _d ∼ 30 K and β = 1.3; Irwin et al. 2001). Given the low angular resolution of our observations, it is unclear how to spatially identify a single structure with the measured B-field orientation. We speculate that a physical structure associated with the extended dusty emission along the NW region at a PA of ∼125° may contain a dominant ordered B-field. This magnetized structure may be associated with one of the Hα (Caulet et al. 1992; Conselice et al. 2001) and CO filaments (Nagai et al. 2019). In addition, the contribution of the synchrotron emission from these filaments at FIR wavelengths is unknown. Further high angular resolution (<1'') IR total and polarimetric observations are required to identify the associated physical structure with the inferred B-field orientation at 89 μm. The 89 μm polarimetric observations of NGC 1275 are inconclusive and will not be used for the physical interpretation discussed in Section 5.6.

At 10–100 pc scales, we have shown that the nuclear thermal dust polarization at 89 μm arises from magnetically aligned dust grains (Sections 3.3 and 5.1) and increases with the nuclear radio-loudness (Section 4). In addition, we show the presence of a 5–130 pc scale ordered toroidal B-field in RL AGNs but not in RQ AGNs (Section 5.4).

At 1–10 pc scales, the inner wall of the dusty torus at scales of 1–2 pc is dominated by an ordered toroidal B-field in the RL AGN IC 5063 (Lopez-Rodriguez et al. 2013) and RQ AGN NGC 1068 (Lopez-Rodriguez et al. 2015). The B-field orientation was measured to be parallel to the equatorial axis of the dusty torus. The B-field strength in the RL AGN IC 5063 (∼41 mG; Lopez-Rodriguez et al. 2013) was measured to be a factor 10 larger than for the RQ AGN NGC 1068 (∼3 mG; Lopez-Rodriguez et al. 2020a) at 1 pc radius from the active nuclei. These results are based on the characterization of hot (∼1300 K) dust at scales of 1 pc from the central engine. ALMA polarimetric observations at 860 μm found a magnetized dusty torus in NGC 1068 (Lopez-Rodriguez et al. 2020a). The polarization arises from magnetically aligned dust grains, where the inferred B-field orientation was measured to be parallel to the equatorial axis within the 3–8 pc eastern side of the dusty torus, i.e., toroidal B-field. In the same region of the torus, the B-field strength was estimated to be ${0.67}_{-0.31}^{+0.94}$ mG. These results provide the B-field configuration and physical mechanisms at the inner boundaries (∼parsec scales) of our findings.

Following these results and assuming conservation of the magnetic energy (i.e., "flux-freezing") and perfect conductivity (i.e., without diffusion), the toroidal B-field strength decreases with distance from the core as B ∝ r⁻¹ (for a toroidal B-field B_ϕ ∼ ρ r; Begelman et al. 1984). For RL AGNs, we take B = 41 mG at 1 pc (Lopez-Rodriguez et al. 2013) and estimate a B-field strength of 8–0.3 mG at 5–130 pc. For RQ AGNs, we take the measured B-field strength of ∼0.67 mG at 5 pc in NGC 1068 (Lopez-Rodriguez et al. 2020a) and estimate a B-field strength of 28 μG at 130 pc from the core. Thus, a stronger B-field associated with the gas flow in RL AGNs may be present at larger distances from the core than in RQ AGNs.

We measured a median offset of 65° ± 26° between the inferred B-field orientation and the radio jet for RL AGNs (Section 5.4). Recent MHD simulations have shown that the infalling gas is not rotationally supported and moves radially inward, which may generate a highly asymmetric disk at scales of 10–100 pc (Takasao et al. 2022). The accreting material from large scales (∼100 pc) to the accretion flow at scales of <10 pc may produce a warped disk (e.g., Lawrence & Elvis 2010). Furthermore, Centaurus A (Graham 1979; Struve et al. 2010), Cygnus A (Struve & Conway 2010), and PKS 1345+125 (Gilmore & Shaw 1986; Lister et al. 2003) show signs of mergers, which may cause a warped disk with a warped B-field orientation (Lopez-Rodriguez 2021) and may also trigger episodes of hot accretion (Sikora & Begelman 2013). Any of these mechanisms (i.e., disk instabilities, mergers) may reproduce the observed offsets between the inferred B-field orientation and the radio jet axis in the core of RL AGNs (Section 5.4).

Our results may be associated with a physical structure in the interface between the sphere of influence of the SMBH and jet production at subparsec scales and the host galaxy at kiloparsec scales. We need to connect our measurements with the physical mechanisms and structures associated with generating B-fields at subparsec and kiloparsec scales.

At kiloparsec scales, the B-fields in the dense, N_HI+H2 = [10²⁰, 10²³] cm⁻², and cold, T _d = [19, 48] K, ISM of nearby galaxies hosting AGNs have been recently traced using FIR polarimetric observations (e.g., Lopez-Rodriguez et al. 2022b). For spiral galaxies, a kiloparsec-scale B-field cospatial with the optical spiral arms has been measured for NGC 1068 (Lopez-Rodriguez et al. 2020b) and NGC 1097 (Beck et al. 2005). NGC 1068 and NGC 1097 have an inner bar of a few kiloparsecs in length, where the B-field is compressed with an orientation parallel to the long axis of the bar. Within the central 1 kpc of NGC 1097, the B-field in the dense ISM (traced in FIR wavelengths) is compressed by a shock driven by the galactic bar dynamics around the starburst ring. The B-field in the diffuse ISM (traced in radio wavelengths) is thought to dominate the infalling material from kiloparsec scale to 100 pc scale toward the AGN (Beck et al. 2005; Lopez-Rodriguez et al. 2021a). Specifically, the mixing of the radio emission from relativistic electrons and the diffuse ionized ISM responds to shear flows, which drag the diffuse gas via an MHD dynamo toward the active nucleus. The B-field strength is estimated to be 60 μG at ∼500 pc from the AGN in NGC 1097 (Beck et al. 2005). These observations show that an MHD dynamo with a kiloparsec-scale spiral B-field allows the flow of gas at kiloparsec scales from the galaxy hosting AGNs to move inward toward the core up to scales of ∼100 pc. These results provide the B-field configuration and physical mechanism at the outer boundaries of our findings, i.e., an ordered spiral B-field on an MHD dynamo feeding gas from kiloparsec to ∼100 pc scales toward the SMBH.

At subparsec scales, we describe two main scenarios. On the one hand, magnetically elevated accretion disks provide the physical structure for the transfer of angular momentum and the production of jets and outflows by the combination of MRI and accretion disk dynamo (Lubow et al. 1994; Guilet & Ogilvie 2012; Begelman & Silk 2017). A weak magnetic field can grow in time to significantly influence the dynamics of the accretion disk. Specifically, the jet power depends on the location of the emitting region along the equatorial plane and the strength and geometry of the B-field at that location. For a sample of 229 blazars, Ghisellini et al. (2014) found that the dominant emitting region along the equatorial plane of the AGN is located at smaller scales than the radius of the BLRs for 85% of the sources and between the BLR and the inner radius of the dusty torus for 15% of the sources. These authors stated that if the emitting region is located at larger distances (≥parsec scales), where the emitting structure is unaffected by external radiation, the jet power will increase. Thus, a strong and ordered toroidal B-field at scales of 1–10 pc will enhance the accretion flow and therefore the jet power.

Under this first scenario, an MHD accretion flow at large radii (up to ∼100 pc) may be the dominant mechanism for the RQ and RL dichotomy. Our results indicate that the stronger the jet power, the stronger and larger toroidal B-field is present affecting the accretion flow at larger distances from the SMBHs. We showed that a 5–130 pc scale ordered toroidal B-field is present in RL AGNs, which we suggest may be the outer layers of the accretion flow of the AGN and the cause of the radio-loudness. Whether the polarized thermal emission structure spatially corresponds with the outer edge of the dusty torus (<50 pc) or a larger (50–100 pc) circumnuclear dusty environment may be an issue of semantics.

This first scenario has several caveats. The dynamics of the gas flow and how the MHD dynamo from large scales interpolates to <130 pc, making the dynamo more efficient, are unknown. In addition, for the magnetically elevated disk to be efficient at these distances, the energy from organized vertical B-fields has to become comparable with the turbulent and thermal energy, and the MRI has to operate in the molecular and dusty disk at <130 pc. Sikora & Begelman (2013) proposed a scenario where hot accretion events can trap magnetic flux in the disk that then is released owing to the loss of angular momentum from the SMBH. Further model development of the effect of MRI in a rapidly rotating flow is required to test this scenario.

On the other hand, the radio jet generates a strong and large-scale toroidal B-field with a current pointing outward (e.g., Begelman et al. 1984). It has also been suggested that, under certain conditions, some of the rotational power of the spinning black hole can be transmitted to the disk to form an MHD wind that drives away most of the mass (Blandford & Globus 2022). This scenario implies that there is an infall of matter at large distances (i.e., Bondi radius) from the SMBH, sufficient to retain the magnetic flux to extract the energy from the SMBH. The median of gas mass in the tori using a sample of 19 Seyfert galaxies was estimated to be ∼6 × 10⁵ M_⊙ (∼10³⁹ g; García-Burillo et al. 2021). We estimated a B-field strength of 8–0.3 mG at a radius of 5–130 pc from the SMBH in RL AGNs. These results are larger than the minimum requirements described in Blandford & Globus (2022), which may satisfy this condition. If the infalling gas has a specific angular momentum >GM_BH/σ, a toroidal B-field could be generated within a disk over the region observed by HAWC+/SOFIA. This scenario also requires a large-scale and toroidal B-field in the presence of a powerful MHD wind, but this will likely arise within the Bondi radius. Under this second scenario, the fact that RL AGNs already have powerful jets may imply that the radio-loudness causes the measured B-field from our observations.

This second scenario has several caveats. The B-field morphology and strength at parsec scales depend on the current being unscreened. It is unknown how far away the outward current, generated by the jet power, closes when the jet propagates into the ISM. In addition, the accretion disk may be highly turbulent, and the B-field outside the jet may diffuse rapidly into a free–free configuration (e.g., Begelman et al. 1984). As an example, the resolved B-fields in the dusty torus of NGC 1068 using 860 μm polarimetric observations showed highly turbulent gas motions on the western side of the equatorial torus. Interestingly, the ordered and toroidal B-field in the dusty torus of NGC 1068 was only measured in the eastern region, which is associated with low turbulent gas motions in the disk (Lopez-Rodriguez et al. 2021b). Finally, this second scenario has been applied to M87, and SOFIA/HAWC+ observations of M87 are not available to further test this model.

The fact that there is an ordered toroidal B-field at scales of 5–130 pc in RL AGNs implies that the MHD models of accretion in SMBHs and jet production may need to be connected with the MHD dynamo in the host galaxy. Our results show that the AGN dusty torus needs to be studied within an MHD framework.