Revisiting the complex kinematics of ionized gas at the central region of NGC 1068: evidence of an additional active galactic nucleus?
DD RAFT VERSION F EBRUARY
11, 2021
Preprint typeset using L A TEX style emulateapj v. 12/16/11
REVISITING THE COMPLEX KINEMATICS OF IONIZED GAS AT THE CENTRAL REGION OF NGC 1068:EVIDENCE OF AN ADDITIONAL ACTIVE GALACTIC NUCLEUS? J AEJIN S HIN , , J ONG -H AK W OO , M INJIN K IM , AND J UNFENG W ANG Department of Astronomy and Atmospheric Sciences, Kyungpook National University, Daegu 41566, Republic of Korea Astronomy Program, Department of Physics and Astronomy, Seoul National University, Seoul, 08826, Republic of Korea Department of Astronomy, Xiamen University, Xiamen, 361005, China
Draft version February 11, 2021
ABSTRACTWe present a spatially resolved analysis of ionized gas at the nuclear region of the nearby galaxy NGC1068. While NGC 1068 has been known to have gas outflows driven by its active galactic nucleus (AGN),more complex kinematical signatures were recently reported, which were inconsistent with a rotation or simplebiconical outflows. To account for the nature of gas kinematics, we performed a spatially resolved kinematicalstudy, finding a morphologically symmetric pair of approaching and receding gas blobs in the northeast region.The midpoint of the two blobs is located at a distance of 180 pc from the nucleus in the projected plane. Theionized gas at the midpoint shows zero velocity and high velocity dispersion, which are characteristics of anoutflow-launching position, as the two sides of a bicone, i.e., approaching and receding outflows are superposedon the line of sight, leading to no velocity shift but high velocity dispersion. We investigate the potentialscenario of an additional AGN based on a multiwavelength data set. While there are other possibilities, i.e.,X-ray binary or supernova shock, the results from optical spectropolarimetry analysis are consistent with thepresence of an additional AGN, which likely originates from a minor merger.
Subject headings: galaxies: active — galaxies: individual (NGC 1068) — ISM: jets and outflows — tech-niques: imaging spectroscopy INTRODUCTIONNGC 1068 at a distance of 14.4 Mpc (Bland-Hawthorn et al.1997) is a prototype Seyfert 2 galaxy with a well-studied ac-tive galactic nucleus (AGN). The first detection of the hid-den broad-line region (BLR) was reported based on spec-tropolarimetric observations (e.g., Antonucci & Miller 1985).Since then, it has been widely accepted that the dichotomy ofSeyfert 1 and Seyfert 2 galaxies is due to an orientation ef-fect caused by a dusty molecular torus surrounding the BLR.The presence of molecular tori has been observationally con-firmed by high spatial resolution observations (e.g., García-Burillo et al. 2016; Imanishi et al. 2016, 2018; Combes et al.2019). The mass of the central black hole was reliably mea-sured as M BH = 0.8–1.7 × M (cid:12) based on spatially resolvedmaser kinematics (Greenhill et al. 1996; Huré 2002; Lodato& Bertin 2003; Impellizzeri et al. 2019), while the bolometricluminosity was determined as L Bol = 0.4–4.7 × ergs − invarious ways using multiwavelength data (Gravity Collabora-tion et al. 2020, and references therein). The correspondingEddington ratio ranges from 0.19 to 4.7, indicating a high ac-cretion rate.One intriguing feature of NGC 1068 is that ionized gaskinematics shows complexity in the central region (i.e., < 500pc) as reported by Walker (1968) and subsequent studies. Us-ing various spectroscopic data, including HST observations,detailed analysis has been conducted to constrain the origin ofthe gas kinematics. Biconical gas outflows have mainly beeninterpreted as being driven by the central AGN and manifestas an approaching gas blob in the northeast (NE) region and areceding gas blob in the southwest (SW) region from the nu-cleus (Cecil et al. 1990; Arribas et al. 1996; Axon et al. 1998;Crenshaw & Kraemer 2000; Cecil et al. 2002; Das et al. 2006; Author to whom any correspondence should be addressed
Gerssen et al. 2006).In contrast, there were reports of an additional receding(redshifted) gas blob, which was detected at 2.5–4.5 (cid:48)(cid:48) (i.e.,180–320 pc) NE of the nucleus (e.g., Axon et al. 1998; Cecilet al. 2002). Because this gas blob disagrees with the trend ofthe well-known biconical outflows (i.e., approaching in NE),more complex mechanisms are likely to be responsible for thechange in the velocity sign in the NE region. While the reced-ing gas blob in the NE region was claimed to be the result of1) the lateral expansion of radio jets (Axon et al. 1998; Ce-cil et al. 2002), or 2) escaped radiation through a patchy dusttorus or scattered radiation (Das et al. 2006), the origin of thecomplex gas kinematics is yet to be clearly understood.To constrain the origin of the complex gas kinematicsat the nuclear region, we performed spatially resolvedspectroscopic analysis by utilizing the VLT/MUSE data. Wedescribe the data and analysis in Section 2. The main resultsand discussion are followed in Section 3 and 4. Conclusionsare presented in Section 5. We adopt a cosmology with H = 70 km s − Mpc − , Ω Λ = 0 . Ω m = 0 . DATA AND ANALYSIS2.1.
Data
NGC 1068 was observed with the VLT/MUSE as a partof the MAGNUM survey (094.B-0321(A) PI: A. Marconi),which covered a 1 (cid:48) × (cid:48) (i.e., 4.3 kpc × a r X i v : . [ a s t r o - ph . GA ] F e b four exposures (observed on 2014 Oct 6) and 100 s for eightexposures (observed on 2014 Dec 1). For this work, we onlyutilized the eight 100 s exposures since [O III ] and H α are sat-urated at a few central pixels in the 500 s exposures. The datawere retrieved from the ESO archive and reduced using thestandard reduction pipeline of the VLT instruments, ESORE-FLEX (Freudling et al. 2013). The seeing of the individualexposures ranged from 0.81 to 0.98 (cid:48)(cid:48) , corresponding to ∼ Analysis
We performed a spectral fitting analysis by focusing on thecentral 14 (cid:48)(cid:48) × (cid:48)(cid:48) (or 1 kpc × ∼ (cid:48)(cid:48) (or ∼
350 pc) region, where ionized gas showed high veloc-ity as well as high velocity dispersion (i.e., V [O III] > 1000 kms − and σ [O III] >1000 km s − , e.g., Cecil et al. 2002; Das et al.2006). Outside the central ∼ (cid:48)(cid:48) region, gas velocity and ve-locity dispersion are typically less than 100 km s − , which arecomparable to those of the stellar component of NGC 1068(e.g., Emsellem et al. 2006). We checked the whole FOVof the MUSE data (1 (cid:48) × (cid:48) ) and confirmed that gas kinemat-ics in large scales generally follows the rotation of the galaxywithout showing a strong outflow signature (i.e., high veloc-ity). We analyzed optical emission lines, i.e., [O III ] λ α , and [Fe VII ] λ λλ III ] & H α ), withthe MPFIT package (Markwardt 2009). To reproduce the ob-served emission lines, we adopted multiple Gaussians (up tothree), whose amplitude-to-noise ratios are all larger than 3.The number of Gaussians was determined with the reducedchi-square of the fitting results. In the fitting, we tied the ve-locity and velocity dispersion for the (1) [O III ] doublet and(2) H α and [N II ] doublet, respectively. We also fixed fluxratios of [N II ] λλ III ] λλ ∼ (cid:48)(cid:48) × (cid:48)(cid:48) region, wherea strong nonstellar continuum is present, leading to unreliablefitting results of the stellar component (see Figure 6 and 7 ofGerssen et al. 2006). The MUSE data also showed a strongnonstellar continuum in the same central region. Therefore,we used another wavelength window, 8400-8800Å, where theCaII triplet is clearly detected representing stellar kinematics.Even in this wavelength range, we also detected a non-stellar continuum and the features of emission lines, whichare associated with e.g., [Cl II] λ λ λ λ RESULT3.1.
Stellar component
In Figure 1, we present the maps of the flux, velocity,and velocity dispersion of the stellar component analyzed inthe two wavelength ranges, (1) 4800-6800Å, and (2) 8400-8800Å. Throughout the paper, we consider the highest fluxpoint in the stellar continuum flux map as the center (0, 0)of NGC 1068. Without any astrometry correction, the cen-ter is offset by only ∼ (cid:48)(cid:48) (i.e., the half of the pixel size ofthe MUSE data) from the position of the peak of X-ray emis-sion at the nucleus of NGC 1068 (RA= 2:42.40.71, Dec= –00:00:47.7, Young et al. 2001). We find that both of themcommonly show a radially decreasing flux and a rotation pat-tern. We also find that the mean difference between the twovelocity maps is 12.5 ±
20 km s − . This means that both veloc-ity maps can be used as a reference for the systemic velocity.In the case of the velocity dispersion maps, however, thereis a significant discrepancy (the right panels of Figure 1).Using the 4800-6800Å range, we find high velocity dis-persion along the NE-SW direction. In contrast, when the8400-8800Å region was used, we detect a trend of σ dropat the central circular region within a radius of ∼ (cid:48)(cid:48) . Thecentral σ drop was detected in a number of galaxies includingNGC 1068 (Emsellem et al. 2004, 2006; Falcón-Barrosoet al. 2006; Gerssen et al. 2006; Comerón et al. 2008; Kanget al. 2013) and its origin has been discussed as recent starformation at the nuclear region (e.g., Emsellem et al. 2001;Comerón et al. 2008; Storchi-Bergmann et al. 2012). We notethat the reason for the high velocity dispersion along the NE–SW in the former one (the top-right panel of Figure 1) wouldbe mainly due to the unreliable fitting result of the stellarcontinuum at the central region as described in Section 2. As F IG . 1.— Flux, velocity, and velocity dispersion maps of the stellar component analyzed in two wavelength ranges, 1) 4800-6800Å (upper) and 2) 8400-8800Å(lower). The center of each map is defined with the the highest flux point (X) in the stellar continuum flux map. Other marks (e.g., R1) will be described inSection 3.2. also presented in Gerssen et al. (2006), the stellar continuumfitting is unreliable at the region with high velocity dispersion.3.2. Ionized gas
We present the flux, velocity, and velocity dispersion mapsof the two major optical emission lines, [O
III ] and H α , as wellas a coronal emission line, [Fe VII ] in Figure 2. [O
III ] andH α represent the ionized gas kinematics in the NLR, while[Fe VII ] with a high ionization potential (99.1 ev) traces thecoronal line region (CLR; e.g., Rodríguez-Ardila et al. 2006;Rodríguez-Ardila & Fonseca-Faria 2020). In the velocitymaps, we use relative velocity (i.e., V gas – V ∗ ) at each spaxelin order to show the relative motion of gas with respect tostellar component, for which velocity is measured from stel-lar absorption line in the 8400-8800Å range (bottom-middlepanel of Figure 1). Because stellar velocity is much smallerthan gas velocity, we find no significant change in the velocitymap with/without subtracting stellar velocity.First, the flux maps show an elongated distribution alongthe NE-to-SW direction, which represents an ionizing cone asdiscussed in the previous studies (e.g., Capetti et al. 1997).Second, the velocity maps present a complex structure, witha wide range of line-of-sight velocity from –700 km s − to+400 km s − ). The much higher gas velocity than that stel-lar velocity indicates that gas kinematics is mainly governed by non gravitational effect, i.e., AGN-driven outflows, insteadof the gravitational potential of the host galaxy. All velocitymaps show a common kinematic structure, by presenting fourhigh-velocity gas blobs (i.e., | V | > 200 km s − ): one pair ofblueshifted (B1) and redshifted (R1) blobs in the NE direc-tion, and another pair of blueshifted (B2) and redshifted blobs(R2) in the SW direction as marked in the [O III ] velocity map.Interestingly, R1 and B1 show morphologically similar tail-like structures, which are elongated roughly in the horizontaldirection, implying that R1 and B1 are physically connectedas a pair. The four gas blobs are also detected in the [Fe
VII ]velocity map, suggesting that the high-velocity gas blobs aredriven by AGN rather than star formation because [Fe
VII ]has a high ionization potential, i.e., 99.1 eV. The positions ofthese gas blobs are slightly different depending on the tracer.For example, there is a ∼ (cid:48)(cid:48) offset of the center of R1 be-tween the [O III ] and H α velocity maps. The spatial offset ispresumably due to the difference in ionizing sources. WhileH α (and also H β ) is ionized by AGN as well as star formation,[O III ] and [Fe
VII ] are mainly ionized by AGN. In addition,the possible uncertainty of the [Fe
VII ] velocity, which can bemuch larger than that of H α and [O III ] because of the lowerflux of [Fe
VII ] (see Figure 2), can be partly responsible forthe offset. Because the exact position of these gas blobs is notthe main interest for our analysis, we determined the locationof gas blobs using the [O
III ] velocity map. F IG . 2.— Flux, velocity, and velocity dispersion maps of [O III ] (top), H α (middle), and [Fe VII ] (bottom). Four high-velocity gas blobs (i.e., R1, R2, B1,and B2) and a region with zero velocity and high velocity dispersion (X1) are marked based on the velocity and velocity dispersion maps of [O
III ]. The nucleusof NGC 1068 is indicated as X. The center of each map is defined with the the highest flux point in the stellar continuum flux map as same as Figure 1.
The high-velocity gas blobs are spatially resolved in theHST narrowband image (Capetti et al. 1997). For example,the gas blob B1 is resolved into a couple of smaller blobs.Without velocity information, however, it is difficult to con-strain how the substructure of B1 detected in the narrowbandimage is related to the outflows. Previous studies of gas kine-matics based on the HST/STIS data showed consistent fea-tures compared to those presented in this work (Cecil et al.2002; Das et al. 2006). For example, a cloud located at thecenter of B1 also showed negative velocities, ∼
600 km s − (see Figure 3 of Cecil et al. 2002).The velocity structure presented in this work is very dif-ferent compared to the biconical outflows in Seyfert galax-ies, which typically show a pair of blueshifted and redshiftedgas blobs (e.g., NGC 5728; Durré & Mould 2019; Shin et al. 2019). While the previous studies treated B1 and R2 as ap-proaching and receding gas outflows, respectively, driven bythe AGN at the nucleus of NGC 1068 (e.g., Das et al. 2006),there are apparently additional gas blobs (R1 and B2), whosegas kinematics is inconsistent with a single biconical out-flows.To better understand the intriguing kinematical structure,we present one-dimensional velocity distributions by locatingtwo pseudo-slits in the [O III ] velocity map: 1) one throughthe nucleus of NGC 1068 and R1 (slit 1, solid line, PA= 44 ◦ )and 2) another along the orientation of the radio jet (slit 2,dotted line, PA= 30 ◦ , Wilson & Ulvestad 1983). As shown inFigure 3, the position angle of Slit 1 is tilted by ∼ ◦ from theorientation of a radio jet, which has PA=30 ◦ (Wilson & Ulves-tad 1983). We find a dramatic variation of gas velocity along F IG . 3.— Top: velocity and flux map of [O III ] with two pseudo-slits, which (1) through the nucleus of NGC 1068 and R1 (slit 1, solid line, PA= 44 ◦ ) and (2)along the orientation of the radio jet (slit 2, dotted line, PA= 30 ◦ , Wilson & Ulvestad 1983). One-dimensional velocity and flux distributions along slit 1 (middle)and 2 (bottom) are shown with the positions of high-velocity gas blobs (i.e., R1). Note that the center of R2 is not exactly on slit 2, while we simply mark R2 atthe position of the foot of perpendicular from the center of R2 to the slit 2. Slit 1 (i.e, B1 to R1), which indicates the change of outflowdirection at the radius of ∼ –190 pc (and also ∼
200 pc). Thechange of gas velocity along Slit 2 (i.e., jet direction throughR2) is slightly less, but clearly shows a change in velocitysign. Moreover, we detect a significant flux bump around R1(middle right panel of Figure 3). While it is difficult to quan-tify the additional flux, the [O
III ] flux around R1 is somewhatincreased ( < ∼ (cid:48)(cid:48) (i.e., ∼
350 pc) from thenucleus, the velocity dispersion (>300 km s − ) is significantlyhigher than stellar velocity dispersion ( ∼
150 km s − ), indi- cating the strong nongravitational kinematics of ionized gas.We find that the velocity dispersion is even higher (up to 800km s − ) at the midpoint between B1 and R1 (marked as X1in Figure 2). In contrast, the gas velocity is almost zero atX1. Because high velocity dispersion and zero velocity aretypically detected at the launching point of gas outflows innearby AGNs (e.g., NGC 5728; Durré & Mould 2019; Shinet al. 2019), these features suggest that the R1 and B1 gasblobs represent outflows, which are launched at the positionof X1. 3.3. Emission-line flux ratios
In order to understand the ionization mechanism, we in-vestigate emission-line flux ratios and identify major ioniz-ing source (i.e., AGN, star-forming, low-ionization nuclearemission-line region, (LINER), or composite region), using F IG . 4.— BPT morphology maps with three diagnostics ([N II ], [S II ], and [O I ]). Red, yellow, cyan, and blue indicate AGN, LINER, composite region, andstar-forming region, respectively. Marks are the same as in the top-middle panel of Figure 2. the BPT diagrams with three diagnoses ([N II ], [S II ], and[O I ], e.g., Baldwin et al. 1981; Kewley et al. 2006). Weconfirm that AGN and LINER are the major ionizing sourceswithin 5 (cid:48)(cid:48) (or ∼
350 pc), while star formation ionizes gas inthe outskirts of the FOV as previously reported (see Figure4; D’Agostino et al. 2019; Mingozzi et al. 2019). The high-velocity gas blobs are mainly ionized by AGN, suggesting thatoutflows in these blobs are likely driven by AGN.In Figure 5, we present the [O
III ]/H β ratio map to traceionization parameter. As expected from the BPT maps, theratio is high (> ∼
5) in the ionized region, where AGN isthe main ionizing source. However, the [O
III ]/H β ratioranges from ∼ ∼
10, and the R1 and B1 blobs showrelatively high [O
III ]/H β ratios. We find two interestingfeatures in the ratio map. First, the highest ratio is detectedin the gas blob B1 (i.e., ∼ (cid:48)(cid:48) NE from the position X),instead of the nucleus. By reporting this feature, previousstudies discussed that shocks from the interaction betweenthe ratio jet and ISM play as an additional ionizing source inelevating the [O
III ]/H β ratio in the NE region (e.g., Kraemeret al. 1998; Kraemer & Crenshaw 2000; Cecil et al. 2002).Second, the region with a high [O III ]/H β ratio (i.e. > ∼ DISCUSSION4.1.
The origin of intriguing kinematics
Based on the spatially resolved analysis of the central re-gion of NGC 1068, we present the complex kinematical sig-nature of ionized gas, which is inconsistent with either thehost galaxy gravitational potential or a simple biconical out-flow scenario. We find a pair of two distinct gas blobs (i.e.,R1 and B1), which show a kinematically similar morphologyin the gas velocity map and share a consistent radial veloc-ity structure with a similar increasing and decreasing trend.Based on the kinemorphological symmetry, we interpret thatthe two blobs represent a pair of receding and approachinggas blobs, which are centered at a launching position (X1) ata projected distance of 180 pc from the galaxy center.While the gas blob B1 is interpreted to be the approachingside of outflows, which are launched from the nucleus (e.g.,Cecil et al. 2002; Das et al. 2006), the gas blob R1 is mov-ing in the opposite direction compared to B1. To explain the F IG . 5.— [O III ]/H β ratio map. Red dashed lines represent a bicone shapecentered at X1. Marks are the same as in the top-middle panel of Figure 2. origin of R1, two scenarios have been suggested: (1) the lat-eral expansion of radio jets (e.g., Axon et al. 1998; Cecil et al.2002) and (2) escaped radiation through a patch torus or scat-tered radiation (Das et al. 2006). These scenarios may explainthe opposite orientation of R1 compared to B1 (see, e.g., Fig-ure 6 of Cecil et al. 2002).However, neither scenario is acceptable because such achange of velocity sign in outflows has not been reported inany other AGNs observed with gas outflows and radio jets.For example, in the cases of well-studied NGC 4151 andNGC 5728, the velocity structure of gas outflows is consis-tent with biconical outflows, manifesting as approaching gasin one side and receding gas in the other side, without chang-ing the velocity sign in one direction (e.g., Das et al. 2005;Durré & Mould 2019; Shin et al. 2019).Moreover, we find that the orientation from the nucleus (X)to the position X1 (i.e., PA= 44 ◦ ) is tilted by ∼ ◦ from thePA of the radio jet, which has PA= 30 ◦ . If the lateral expan-sion of radio jets and/or escaped radiation is the origin of R1,it is natural to expect a similar gas blob, which is laterally ex-panded from the jet and located in the opposite direction (i.e.,PA= 16 ◦ ; see the slit position in Figure 3). In other words,a symmetric gas blob would be located on the other side ofthe jet. In contrast, we find no structure in that position (at ∆ RA= ∼ (cid:48)(cid:48) , ∆ Dec= ∼ (cid:48)(cid:48) ) while it is possible that the lack ofionized gas in that location may explain why the lateral ex-pansion or escaped radiation is morphologically asymmetric(top-right panel of Figure 3). In addition, it was pointed outthat the radio jets are not energetic enough to drive strong gasoutflows (Das et al. 2006). Therefore, the two aforementionedscenarios cannot fully explain the observed flux and velocitystructures at around R1, implying that the AGN at the nucleusof NGC 1068 may not be responsible for the kinematics of thegas blob R1.We propose a new scenario where there is an additional(second) mass-accreting black hole at around the position X1,and the outflows driven by the 2nd AGN at X1 are manifestedby a pair of receding (R1) and approaching (B1) gas blobs.We present several supporting evidence. First, as presented inSection 3.2, ionized gas at X1 shows zero velocity and highvelocity dispersion, which are typical observational character-istics of an outflow-launching point in nearby Seyfert galax-ies (e.g., NGC 5728; Durré & Mould 2019; Shin et al. 2019).This trend is due to the fact that the two sides (i.e., recedingand approaching) of the gas outflows are superposed alongthe line of sight at their launching point because of the beamsmearing effect, hence, the cancellation between positive andnegative velocities leads to net zero velocity, while the broadvelocity distribution results in high velocity dispersion (e.g.,Durré & Mould 2019; Shin et al. 2019). Therefore, our sce-nario of a 2nd AGN can naturally explain the gas kinematicsat X1. Second, the morphological symmetry (i.e., the elon-gated tail-like structures) between R1 and B1 in the velocitymap suggests that they are physically related, sharing a sameorigin (i.e., the 2nd AGN). Third, we detect the radial veloc-ity pattern of acceleration followed by deceleration for bothof the gas blobs R1 and B1, which are centered at the po-tential launching position X1 (see Figure 3). This increasingand decreasing velocity pattern is a signature of gas outflowsdetected in a number of nearby AGNs (Crenshaw & Krae-mer 2000; Crenshaw et al. 2000; Müller-Sánchez et al. 2011;Durré & Mould 2019). While the radial velocity pattern, aswell as the relatively high gas velocity in NGC 1068 werealready noticed in previous studies (Crenshaw & Kraemer2000; Müller-Sánchez et al. 2011), gas outflows are assumedto be launched from the nucleus. Fourth, we find a biconicalshape centered at X1 in the [O III ]/H β ratio map, which is pre-sumably due to an ionizing source at X1 (i.e., the 2nd AGN).Note that the approaching gas blob B1 was considered to bedriven by the 1st AGN at the nucleus in the previous studies(e.g., Cecil et al. 2002; Das et al. 2006; Gerssen et al. 2006).However, the gas blob B1 is more extended than R1, coveringthe two slit positions shown in Figure 3, and it is likely thatB1 is the combination of the approaching gas blobs, which arelaunched form the nucleus and the position X1, respectively(see the schematic view of two pairs of biconical outflows inFigure 6).We investigate whether the kinematics of gas blobs R1 andB1 can be explained by a rotating disk, which is centered atX1. As shown in the lower panel of Figure 3, the peak velocityof B1 and R1 reaches ∼
500 km s − at the projected distance ofonly ∼
50 pc from X1. To produce a rotational velocity of 500km s − at 50 pc, the required enclosed mass needs to be a min-imum of 10 . M (cid:12) when the disk is assumed to be edge on. For the case of NGC 1068, however, the enclosed mass withinthe central 50 pc is only 10 . M (cid:12) , as estimated based on thestellar rotational velocity (i.e., ∼
30 km s − at 50 pc from thenucleus) and the inclination (i.e., 40 ◦ , Bland-Hawthorn et al.1997). Because the expected mass of the gas disk at X1 is afactor of ∼
200 higher than that based on stellar kinematics,a rotating gas disk is not likely responsible for the kinematicsof gas blobs B1 and R1. Previous studies of nearby AGNs, in-cluding NGC 1068, pointed out the same conclusion that suchhigh velocity (> 300km s − ) of ionized gas is due to gas out-flows, not a galactic disk rotation (e.g., Müller-Sánchez et al.2011; Durré & Mould 2019).Alternatively, supernova (SN) may be the origin of R1 asthere was a report of high-velocity molecular gas outflows( V ∼ − ) driven by a SN in NGC 6240 (Treister et al.2020). However, SN is not likely responsible for R1, becauseof the following two reasons. First, as shown in Figure 4, themajor ionizing source for R1 is AGN (see also D’Agostinoet al. 2019; Mingozzi et al. 2019). While D’Agostino et al.(2019) found significant shock contribution in the central re-gion ( ∼ × (cid:48)(cid:48) NE ofthe nucleus) is mainly photoionized rather than collisionallyionized, implying that AGN (not SN) is the major ionizingsource of the region. In fact, our BPT analysis clearly showedthat the photoionization of R1 is dominated by AGN (see Fig-ures 4 and 5). Note that we cannot rule out the SN origin ofR1 if the SN-driven photoionization is relatively weak and un-detectable as the BPT map requires. In this case, it is not clearhow the SN provides a strong impact on the gas kinematics.One may expect a light excess around X1 in the flux map, ifthere is indeed a mass-accreting black hole. While we find nostrong evidence, we detect an excess of [O
III ] flux around X1in the flux map of ionized gas (see Figure 3). In contrast, thereis no signature of dynamical disturbance in the stellar flux andkinematical maps, which may suggest that the galaxy-galaxymass ratio is very large and the contribution of the mergedcomponent is negligible and not detected in the flux-weightedstellar absorption lines. The possibility of minor mergers inNGC 1068 has been discussed based on the detections of (1)ultradiffuse objects around NGC 1068 (Tanaka et al. 2017),(2) an off-centered circumnuclear disk (CND, García-Burilloet al. 2014), and (3) the non-Kerplerian motion of moleculargas in the molecular torus (e.g., Imanishi et al. 2018). We notethat García-Burillo et al. (2014) also discussed the presenceof another AGN in the SW of the nucleus ( ∆ RA= –1.50 (cid:48)(cid:48) , ∆ Dec= –1.00 (cid:48)(cid:48) ) to explain the off-centered CND even thoughit is yet to be observationally confirmed.Recently, Wang et al. (2020) claimed a close binary super-massive black hole (SMBH) with a separation of ∼ X X1 X X1 X X1 F IG . 6.— Schematic view of gas outflows in NGC 1068 overlaid on the velocity map of [O III ]. Gas outflows launched from the nucleus, X (left, PA= 30 ◦ ),X1 (middle, PA= 44 ◦ ), and both of them (X and X1, right) are presented. Solid and dotted lines are the same as in Figure 3. this work.4.2. Counterpart in multiwavelength data
To investigate whether the 2nd AGN scenario is consistentwith other observational signatures, we use X-ray and radiodata and search for evidence of multiple cores at the cen-tral region. Starting with X-ray data, we check the Chan-dra high-resolution camera (HRC) data of NGC 1068 (Ob-sID: 12705; PI: Fabbiano), which were previously presentedby Wang et al. (2012). As shown in the top-right panel of Fig-ure 1 of Wang et al. (2012), the Chandra image with the PSFdeconvolution resolves another putative X-ray point source at3.6 (cid:48)(cid:48)
NE of the nucleus ( ∆ RA= +2.00 (cid:48)(cid:48) , ∆ Dec= +2.75 (cid:48)(cid:48) . Thedetection of the point source is significant at the >7 σ levelwhen we take into account the Poisson noise of the adjacentpixels. The position of the X-ray point source is close to thecenter of the pair of gas blobs R1 and B1 (i.e., X1 position),albeit with an offset of ∼ (cid:48)(cid:48) to the NE of X1, and it is likelyto represent the 2nd AGN around X1. The spatial offset canbe explained by various effects (i.e., dust obscuration and ISMimpact), which can cause the positional shift of the kinematiccenter of the gas outflows (i.e., X1) from the position of an X-ray point source, where an AGN is located (see e.g., Durré &Mould 2018). For example, a ∼ (cid:48)(cid:48) separation was found be-tween the kinematic center of gas outflows and the position ofthe central X-ray point source in NGC 5728 (Durré & Mould2019; see also Shin et al. 2019).To investigate the possibility of the X-ray point-source asa potential AGN, we measure the X-ray luminosity in 2–10keV by converting the X-ray photon counts (133) at the pixelof the X-ray point source position, using the PSF deconvolvedChandra HRC image. For this practice, we calculate a conver-sion factor from WEBPIMMS using the effective area curve of Chandra
Cycle 12, a power law index ( Γ = 1.8), and a Galac-tic column density ( N H = 2 . × cm − , Murphy et al.1996). As a result, we obtain the 2-10 keV X-ray luminos-ity as ∼ . ± . × erg s − . Note that the uncertainty isfrom the Poisson statistics. The X-ray luminosity is not suf-ficiently high to confirm the presence of an additional AGNbecause X-ray binaries can be the origin of the X-ray emis-sion. Note that the Eddington luminosity of a 10 M (cid:12) stellarblack hole is ∼ erg s − , while the luminosity of high-massX-ray binaries is typically less than 10 - 10 erg s − (e.g., http://cxc.harvard.edu/toolkit/pimms.jsp Grimm et al. 2003; Fabbiano 2006; Mineo et al. 2012).We investigate the possibility that the luminosity of the X-ray point source is underestimated as the potential 2nd AGNis likely to be type 2 without presenting broad emission linesin the optical spectrum. In this case, we do not directly mea-sure the intrinsic luminosity because of the obscuration by adust torus. It has been known for the 1st AGN in NGC 1068that the transmitted X-ray emission from the central sourceis completely obscured by a dust torus as the column den-sity of this Compton-thick source is N H ∼ cm − . Thus,it is difficult to measure the intrinsic X-ray luminosity as weonly observe the X-ray emission reflected by the inner wallof the dust torus and ionized gas (Colbert et al. 2002; Ogleet al. 2003; Matt et al. 2004; Bauer et al. 2015, but see alsoMarinucci et al. 2016). To determine the intrinsic X-ray lumi-nosity, for example, Colbert et al. (2002) used the [O III ] lu-minosity as a proxy for the bolometric luminosity of the cen-tral AGN, estimating the intrinsic X-ray luminosity of ∼ × erg s − , which is a factor of 250 higher than the re-flected X-ray luminosity. More reliably, Bauer et al. (2015)performed X-ray spectral analysis, reporting the intrinsic X-ray luminosity of the central AGN as 2.2 × erg s − (seeMarinucci et al. 2016).In the case of the potential 2nd AGN, we cannot use thesame correction factor obtained for the 1st AGN because the2nd AGN is not likely a Compton-thick source. Instead, aproper spectral analysis is required to determine the intrinsicX-ray luminosity. Note that the AGN X-ray spectrum is com-plex as it is composed of various components, including softX-ray access and a power-law component. Due to the lack ofspectral information on the X-ray source around X1, we esti-mate the intrinsic X-ray luminosity from the observed X-rayluminosity by assuming a typical column density for type 2AGNs and calculate the correction factor with WEBPIMMS .We obtain a range of the obscuration correction factor from1.06 and 7.14, depending on the column density N H = 10 toN H = 10 . As an example, if we multiply the correction fac-tor of 1.56, corresponding to N H = 10 , to the observed X-rayluminosity of ∼ . × erg s − , the intrinsic X-ray lumi-nosity is 1 . × erg s − . While we cannot measure theexact correction factor for the X-ray source due to the lackof X-ray spectral information, it is likely that the intrinsic X-ray luminosity is not sufficiently high to rule out a high-massX-ray binary as an origin. https://heasarc.gsfc.nasa.gov/cgi-bin/Tools/w3pimms/w3pimms.pl In addition, we analyze Advanced CCD Imaging Spectrom-eter (ACIS) data (ObsID= 370, 0.4 s frame, exposure time=11.5 ksec), which were presented by Young et al. (2001), inorder to check whether the X-ray spectrum of the X-ray pointsource shows a power-law continuum, as a power-law contin-uum is another observational signature of AGNs. Note thatwe do not investigate deeper data (ObsID= 344, 3.2 s frame,exposure time= 47.4 ksec), which were heavily affected bypile-up at around X1 (and also the nucleus) because of thehigh count rate (see Young et al. 2001). We extract the ACISspectrum from the position of the X-ray point source as de-tected in the Chandra HRC image with a radius of 0.5 (cid:48)(cid:48) . Then,we fit the spectrum with a various combinations of 1) a single-temperature plasma model (MEKAL), 2) bremsstrahlung, 3)power law, and 4) a thermal plasma component (APEC). Allof the models include the absorption by the Galactic columndensity ( N H (Gal) = 2.99 × cm , Murphy et al. 1996) assimilarly done by Young et al. (2001). However, we find thatthe inclusion of a power law does not improve the fitting re-sults compared to those without a power law.To further investigate the origin of the X-ray point source,we measure the [O III ]/soft X-ray ratio, which is a tracer ofionization states (Bianchi et al. 2006). By investigating the ra-tio for several X-ray knots in the central region of NGC 1068,Wang et al. (2012) reported that a few of them were originatedfrom shocks. Using the same data set and method of Wanget al. (2012), we find the ratio of ∼ . ± .
86 at the pixelof the X-ray point-source position (i.e., within the extractionarea of 0 . (cid:48)(cid:48) × . (cid:48)(cid:48) ). This result indicates photoionizationrather than shock. Note that the flux ratios, [O III ]/H β and[N II ]/H α , [S II ]/H α , or [O I ]/H α obtained around the X1position are also consistent with AGN photoionization in theBPT diagrams. Thus, we conclude that X-ray emission is un-likely to originate from shocks unless [O III ] flux is stronglydominated by the photoionization by the 1st AGN.Secondly, we check the ALMA data of CO (3-2) (García-Burillo et al. 2014), and dense molecular gas (i.e., HCN (3-2)and HCO + (3-2)), presented by Imanishi et al. (2018), as wellas the VLA radio data by used by Wilson & Ulvestad (1983);Gallimore et al. (1996), while we do not find any strong sig-nature (i.e., radio core) of the 2nd AGN around X1. This canbe explained with the low luminosity and/or the low mass ofthe host galaxy of the 2nd AGN, as discussed in the previoussection.Overall, we find no clear evidence of a second AGN basedon the multiwavelength data because the X-ray luminosity isrelatively low and comparable to high-mass X-ray binaries.Nevertheless, we find no evidence to rule out the possibilityof the 2nd AGN.4.3. Hidden broad emission line
One of the characteristics of an AGN is the presence ofbroad emission lines (i.e., FWHM > 2000 km s − ), which areemitted from the sub-pc scale broad line region. For Seyfert2 galaxies, the broad line region is obscured by a dust torus,but it can be detected through spectropolarimetric observa-tions. As described in Section 1, NGC 1068 is the first type2 Seyfert galaxy whose hidden broad Balmer lines were de-tected (e.g., Antonucci & Miller 1985). Since then, a numberof spectropolarimetric results have been presented for NGC1068. Among them, spatially resolved spectropolarimetricstudies (Inglis et al. 1995, see also Miller et al. 1991) found F IG . 7.— One-dimensional distribution of the FWHM (upper) and lineflux (lower) of the hidden broad component of H α . A dotted line denotes theposition of X1. All measurements are adopted from Inglis et al. (1995). the interesting result that the line-width of the hidden broadBalmer lines detected in the NE (e.g., ∼ ±
400 km s − at 5 (cid:48)(cid:48) , or 360 pc NE) of the nucleus is significantly smallerthan that detected at the nucleus ( ∼ ±
300 km s − ). In-glis et al. (1995) presented the line width measurements of thehidden broad H α line, which were detected at multiple pointsfrom the nucleus to the NE direction (i.e., 2.5, 5, and 7.5 (cid:48)(cid:48) offset, corresponding to 180, 360 and 540 pc distance fromthe nucleus), show a consistency ( ∼ − ) within 1- σ (see the upper panel of Figure 7). This result may suggest thatthe hidden broad H α detected in the NE has a different ori-gin from that at the nucleus. One possible explanation is thatthere is an additional AGN in between 2.5–7.5 (cid:48)(cid:48) (or 180–540pc) NE of the nucleus, whose hidden broad lines are intrinsi-cally narrower than those of the 1st AGN at the center becauseof the smaller black hole mass of the 2nd AGN.As shown in the lower panel of Figure 7, the line flux ofthe hidden broad H α line rather increases at 360 pc NE of thenucleus, which is in disagreement with the trend of radiallydecreasing flux if it originated from the nucleus. This trend isconsistent with the scenario of a 2nd AGN at a distance of 360 ±
90 pc NE from the nucleus. Note that the kinematic centerof R1 and B1 is X1, which is positioned at 180 pc from the nu-cleus; however, the true location of an AGN can be somewhatoffset from the kinematical center as discussed in Section 4.2.It is worth noting that imaging polarimetric observations havedetected a highly polarized knot at 4.7 (cid:48)(cid:48)
NE of the nucleus(Miller et al. 1991; Scarrott et al. 1991; Simpson et al. 2002).It could be also associated with the second AGN while it is 2 (cid:48)(cid:48) offset from X1.In previous works, the origin of the broad Balmer emis-sions and high polarized intensity in the NE region wasinterpreted as the scattered light from the nucleus (e.g., Milleret al. 1991; Inglis et al. 1995) and molecular gas, whichprevents the expansion of the radio jets (Simpson et al. 2002),respectively. However, these scenarios cannot clearly explainthe difference of the line widths in the NE region as well asthe bump in the flux distribution at 360 pc. In contrast, thescenario with the 2nd AGN positioned at around X1 mayprovide a better explanation.0 4.4.
Mass outflow rate and Energetics
We discuss whether the energetics of the potential 2nd AGNis eligible to drive gas outflows manifested by gas blobs B1and R1, by investigating the mass outflow rate and kinetic en-ergy output. We use two different methods, namely, geometricand luminosity approaches (see e.g., Revalski et al. 2018).4.4.1.
Geometric approach
For the geometric approach, the mass outflow rate and en-ergy budget can be calculated as follows (see, e.g., the equa-tion 1 of Müller-Sánchez et al. 2011), ˙ M out = 2 m p N e v max A f (1) ˙ E out = 1 / ˙ M out ( v + σ ) (2)This approach is based on a geometric assumption by deter-mining five parameters: (1) the electron density, N e , (2) thedeprojected maximum outflow velocity, v max , (3) the lateralsurface area of outflows at the maximum velocity position, A ,(4) the averaged velocity dispersion in gas outflows, σ , and (5)the filling factor, f . Because we interpret B1 as a superposedoutflows, respectively launched from X and X1, we only useR1 to calculate outflow parameters. First, we estimate theelectron density based on the [S II ] λ II ] λ (cid:48)(cid:48) radius from the maximum velocity position in R1.By assuming the electron temperature to be 10000 K, we ob-tained an electron density of 907 cm − , which is consistentwithin ∼
10% with the previously reported value Mingozziet al. (2019, i.e., ∼
800 cm − ). Second, we adopt the (pro-jected) maximum velocity of [O III ] in R1 (i.e., 291 km s − )as the lower limit from the velocity map (see Figures 2 and3) because the inclination of the outflows is unknown. Third,we estimate the lateral surface by determining the size of themajor axis of R1. If outflows have a biconical shape, the lat-eral surface would be round, which is shown as an ellipse inthe plane of the sky. Thus, the size of the major axis can beused to calculate the area of a circle. Based on this idea, wedetermine the size of the major axis to be 2 (cid:48)(cid:48) (i.e., 142 pc)by drawing a line through the maximum velocity position inR1, which is perpendicular to the gas outflow direction (i.e.,Slit 1; see Figure 3). To account for the hollow shape of out-flows (e.g., Müller-Sánchez et al. 2011; Bae & Woo 2016;Shin et al. 2019), we adopt the hollow cone geometry of the1st gas outflows, which were constrained as outer angle= 27 ◦ and inner angle= 14 ◦ ) based on biconical gas outflow model-ing by Müller-Sánchez et al. (2011). Using these estimates,we calculate the lateral surface area as 1.3 × pc . Fourth,the average velocity dispersion in R1 is measured as 485 kms − using the same spectrum used for electron density estima-tion. Finally, we adopt the filling factor of 0.11 from Storchi-Bergmann et al. (2010), who calculated the filling factor forgas outflows in NGC 4151. Because we adopt the projectedmaximum velocity, these parameters provide the lower limitof the mass outflow rate and energy outflow rate as ˙ M out =18.1 M (cid:12) yr − and ˙ E out = 1.9 × erg s − , respectively.For a consistency check, we calculate the energetics of thegas outflows (R2) driven by the 1st AGN. Applying the samemethod to R2, we measure the averaged electron densityin R2 (479 cm − ), the maximum velocity (321 km s − ), thelateral surface area (0.8 × pc ), and the averaged velocitydispersion in R2 (544 km s − ). By adopting the same filling factor of 0.11, we obtain the lower limit of the mass outflowrate as 6.5 M (cid:12) yr − and the energy outflow rate as 8.2 × erg s − . If we adopt the inclination of the 1st gas outflows(i.e., 9 ◦ ) presented by Müller-Sánchez et al. (2011), thedeprojected (intrinsic) maximum velocity becomes 2057 kms − , which is consistent within 10% with that constrained byMüller-Sánchez et al. (2011). We find that the mass outflowrates and kinetic energy are ˙ M out = 38.8 M (cid:12) yr − and ˙ E out =4.8 × erg s − , respectively, which are consistent within anorder of magnitude with those reported by Müller-Sánchezet al. (2011). 4.4.2. Luminosity approach
The luminosity approach calculates mass outflow rate andenergy budget as follows, ˙ M out = 3 M gas v out R out (3) ˙ E out = 12 ˙ M out v , (4)where M gas is the ionized gas mass, v out is the averaged out-flow velocity (i.e., v out = (cid:113) v + σ ) within the outflows, and R out is the outflow size (e.g., Karouzos et al. 2016b; Baeet al. 2017; Kang & Woo 2018). First, the ionized gasmass is converted from the [O III ] luminosity with an esti-mated electron density using the equation of M gas = 0 . × M (cid:12) × ( L [OIII] , )(100 cm − / n e ) (see, e.g., Carniani et al.2015; Karouzos et al. 2016b). In the previous works by ourgroup (e.g., Karouzos et al. 2016b; Bae et al. 2017), the ion-ized gas mass was estimated using the integrated [O III ] lu-minosity within the outflow size. If NGC 1068 hosts anadditional AGN along with its corresponding gas outflows, itis not trivial to measure the [O
III ] emission solely due to thephotoionization by the 2nd AGN. For simplicity, we constrainthe upper limit of the ionized gas mass by combining [O
III ]flux around X1. We calculate twice the integrated [O
III ] lu-minosity in R1 (4.6 × erg s − ) within a 1 (cid:48)(cid:48) radius as the[O III ] luminosity and calculate the ionized gas mass to beM gas = 2 × M (cid:12) using the average electron density in thisregion (907 cm − ). Note that assuming a bicone, we multi-ply a factor of 2 to the [O III ] luminosity of R1. Becausea large fraction of the warm gas is photoionized by the 1stAGN, the derived gas mass should be considered as an upperlimit. While the [O
III ] map in Figure 2 shows a dominanceof photoionization by the 1st AGN even at the R1 area, wedetect additional flux at X1 and R1 areas as shown in Fig-ure 3, which is presumably due to the photoionization by the2nd AGN. While it is difficult to estimate the fraction of theionized gas flux due to the 2nd AGN with respect to the total[O
III ] flux in R1 area, the fraction should not be negligiblebecause we clearly detect distinct kinematical signatures inthe velocity and velocity dispersion maps based on the flux-weight [O
III ] line profile. If we assume the [O
III ] flux dueto the 2nd AGN is ∼ − ) and veloc-ity dispersion (485 km s − ) in R1. Third, we determine theoutflow size as 140 pc (or ∼ (cid:48)(cid:48) ), which is estimated at the po-sition where the [O III ] velocity approaches 100 km s − (see ∼ –330 pc in the middle panel of Figure 3). Note that the size1of outflows is uncertain, and we simply use this value as ourbest approximation. As a result, we obtain the upper limit ofthe gas mass and energy outflow rate as ˙ M out = 0.24 M (cid:12) yr − .and ˙ E out = 2.2 × erg s − . The actual gas mass and energyoutflow rate should be corrected by the fraction of the [O III ]flux due to the 2nd AGN over the total [O
III ] flux.4.4.3.
Comparison of outflow energetics with bolometricluminosity
Based on the two aforementioned methods, we estimate themass outflow rate and energy outflow rate. However, thereis a large discrepancy between them by ∼ ˙ E out /L bol ). To estimate the bolomet-ric luminosity of the 2nd AGN, we utilize two different indi-cators, X-ray luminosity and [O III ] luminosity. If we use theluminosity of the X-ray source around X1 with an obscurationcorrection factor of 1.56 for N H = 10 cm − and a bolometriccorrection factor of 10.85, we obtain the bolometric luminos-ity of the 2nd AGN as L bol = 1.5 × erg s − . Here the bolo-metric correction factor is calculated based on the Equation 2of Duras et al. (2020) with the coefficients for type 2 AGNs.Note that L bol could be lower by a factor of ∼ ∼
7, if N H is 10 cm − or 10 cm − , respectively.Also, note that the bolometric correction for the hard X-ray isrelatively uncertain. On the other hand, if we use the upperlimit of [O III ] luminosity around R1, along with a bolometriccorrection of 3500 (Heckman et al. 2004), we derive the up-per limit of L bol = 1.6 × erg s − . If we assume 1% of thetotal [O III ] in R1 is photoionized by the 2nd AGN, we obtainL bol = 1.6 × erg s − . Note that the large difference be-tween X-ray based and [O III ] based bolometric luminositiesis mainly due to various uncertain factors, including the ob-scuration correction factor and bolometric correction of hardX-ray luminosity as well as the uncertain fraction of the [O
III ]flux due to the 2nd AGN around R1 area.Comparing with these estimates of bolometric luminosity,the energy outflow rate determined by the geometric approachis challenging because the kinetic energy carried by the out-flow is much larger than the X-ray based bolometric lumi-nosity, with the kinetic coupling efficiency of ∼ ∼ III ] flux in the R1 area is due to the 2nd AGN, we obtain ˙ E out = 2.2 × ,which is ∼ − of the X-ray based and [O III ]-basedbolometric luminosities, respectively. These results are withinthe range of the kinetic coupling efficiencies of warm ionizedgas reported in the literature (Harrison et al. 2018). Albeitwith large systematic uncertainties of outflow energetics andAGN bolometric luminosity due to a number of assumptions,we find no strong evidence to rule out the scenario of outflowsdriven by a 2nd AGN. CONCLUSIONIn this work, we investigate the complex kinematics ofionized gas at the central region of NGC 1068. The mainresults are summarized below. • We detect a pair of blueshifted and redshifted gas blobs,which are located in the NE region at ∼
110 and ∼