Separating Line Emission from Star Formation, Shocks, and AGN Ionisation in NGC 1068

TL;DR

This paper introduces a new three-dimensional diagnostic diagram to effectively separate and quantify line emission contributions from star formation, shocks, and AGN in NGC 1068, improving understanding of galaxy ionisation sources.

Contribution

The study presents a novel 3D diagnostic method for disentangling emission sources in galaxies, enabling more accurate analysis of ionisation mechanisms.

Findings

Successful separation of emission sources in NGC 1068

Alignment of emission distributions with X-ray, radio, and CO observations

Shocked emission significantly mixes with AGN emission in previous methods

Abstract

In the optical spectra of galaxies, the separation of line emission from gas ionised by star formation and an AGN, or by star formation and shocks, are very well-understood problems. However, separating line emission between AGN and shocks has proven difficult. With the aid of a new three-dimensional diagnostic diagram, we show the simultaneous separation of line emission from star formation, shocks, and AGN in NGC 1068, and quantify the ratio of star formation, shocks, and AGN in each spaxel. The AGN, shock, and star formation luminosity distributions across the galaxy accurately align with X-ray, radio, and CO(3-2) observations, respectively. Comparisons with previous separation methods show that the shocked emission heavily mixes with the AGN emission. We also show that if the H$\alpha$ flux is to be used as a star formation rate indicator, separating line emission from as many sources as possible should be attempted to ensure accurate results.

Tables

Table 1. Table 1 : Metallicity gradient and 100% emission line ratio values for NGC 1068. Solar metallicity ( Z ⊙ subscript 𝑍 direct-product Z_{\odot} ) is set to 12 + log(O/H) = 8.93 (Anders & Grevesse, 1989 ) .

Outer metallicity ($Z_{\odot }$)	1.48
Central metallicity ($Z_{\odot }$)	2.20
100% star formation ELR function value	-0.01
100% shock/AGN ELR function value	0.81

Table 2. Table 2 : Star formation, shock, and AGN fractions for various strong emission lines in NGC 1068 within the S7 field-of-view. Errors on the star formation, shock, and AGN fractions account for differences in final emission line weighting, uncertainties in photoionisation models, and uncertainities in the fits from lzifu . Displayed also are the star formation and AGN fractions using the method from D’Agostino et al. ( 2018 ) , assuming only star formation-AGN mixing. The errors associated with the star formation-AGN mixing fractions are the result of a 0.1 dex uncertainty in the photoionisation models, and the errors from lzifu .

NGC 1068
	New method			Old method
Emission line	Star formation	Shocks	AGN	Star formation	AGN
H$\alpha$	59.7 $\pm$ 8.0%	26.1 $\pm$ 7.4%	14.2 $\pm$ 4.5%	64.4 $\pm$ 7.1%	35.6 $\pm$ 7.1%
H$\beta$	60.0 $\pm$ 8.1%	26.0 $\pm$ 7.4%	14.1 $\pm$ 4.5%	64.6 $\pm$ 7.1%	35.4 $\pm$ 7.1%
[O ii]$\lambda \lambda$3726,3729	35.4 $\pm$ 10.0%	41.2 $\pm$ 12.1%	23.3 $\pm$ 8.8%	38.3 $\pm$ 12.0%	61.7 $\pm$ 12.0%
[O iii]$\lambda$5007	36.1 $\pm$ 10.8%	46.1 $\pm$ 11.5%	17.8 $\pm$ 6.5%	39.6 $\pm$ 12.8%	60.4 $\pm$ 12.8%
[S ii]$\lambda \lambda$6716,6731	39.4 $\pm$ 9.5%	39.1 $\pm$ 11.3%	21.5 $\pm$ 8.0%	42.5 $\pm$ 11.1%	57.5 $\pm$ 11.1%
[N ii]$\lambda$6584	43.9 $\pm$ 9.5%	37.6 $\pm$ 10.3%	18.5 $\pm$ 6.7%	48.0 $\pm$ 10.4%	52.0 $\pm$ 10.4%

Equations

$$\mathrm{ELR}\mathrm{function}=\frac{\log ([\mathrm{N}\text{textsc}\mathrm{ii}]/H\alpha )-{\min }_{\log ([\mathrm{N}\text{textsc}\mathrm{ii}]/H\alpha )}}{{\max }_{\log ([\mathrm{N}\text{textsc}\mathrm{ii}]/H\alpha )}-{\min }_{\log ([\mathrm{N}\text{textsc}\mathrm{ii}]/H\alpha )}}$$

$$\times \frac{\log ([\mathrm{O}\text{textsc}\mathrm{iii}]/H\beta )-{\min }_{\log ([\mathrm{O}\text{textsc}\mathrm{iii}]/H\beta )}}{{\max }_{\log ([\mathrm{O}\text{textsc}\mathrm{iii}]/H\beta )}-{\min }_{\log ([\mathrm{O}\text{textsc}\mathrm{iii}]/H\beta )}}$$

$$L_{\mathrm{Tot}}=\sum _{i=1}^n{L}_i$$

$$L_{\mathrm{source}}=\sum _{i=1}^n{f}_i^{\mathrm{source}}{L}_i$$

$$f_{\mathrm{Tot}}^{\mathrm{source}}=\frac{L_{\mathrm{source}}}{L_{\mathrm{Tot}}}$$

Full text

††pubyear: 2018††pagerange: Separating Line Emission from Star Formation, Shocks, and AGN Ionisation in NGC 1068–Separating Line Emission from Star Formation, Shocks, and AGN Ionisation in NGC 1068

1 Introduction

One of the biggest mysteries in modern astrophysics is the link between supermassive black hole (SMBH) accretion and the evolution of its host galaxy. Studies have shown SMBH mass correlates with other galaxy properties, such as the velocity dispersion ($M$-$\sigma$ relation; e.g. Ferrarese & Merritt, 2000; Gebhardt et al., 2000; Tremaine et al., 2002; Gültekin et al., 2009; McConnell & Ma, 2013), the stellar mass in the bulge ($M_{BH}$-$M_{*}$ relation; e.g. Magorrian et al., 1998; Marconi & Hunt, 2003; Bennert et al., 2011; McConnell & Ma, 2013), and the luminosity of the bulge ($M_{BH}$-$L$; e.g. Marconi & Hunt, 2003; Gültekin et al., 2009; McConnell & Ma, 2013). However, no such theory linking SMBH accretion and star formation, such as mergers, starburst-driven winds or AGN-driven outflows (e.g. Yuan et al., 2010; Rafferty et al., 2011) have been convincing. As a result, no theoretical model has successfully been able to explain the relationship between star formation and active galactic nuclei (AGN; see the review by Alexander & Hickox, 2012).

Work on uncovering the link between star formation and AGN activity was explored by Kewley et al. (2001) through the use of several emission-line-ratio diagnostic diagrams given by Baldwin et al. (1981) and Veilleux & Osterbrock (1987). Kewley et al. (2001) defined a curve on each of the [O iii]/H$\beta$ vs [N ii]/H$\alpha$ (commonly referred to as the ‘BPT diagram’; Baldwin et al., 1981), [O iii]/H$\beta$ vs [S ii]/H$\alpha$, and [O iii]/H$\beta$ vs [O i]/H$\alpha$ diagrams to represent the theoretical maximum of line emission from star formation, derived from photoionisation modelling using the mappings iii code. Kauffmann et al. (2003) also defined a maximum starburst line on the BPT diagram, albeit empirically from their sample of SDSS galaxies. Contemporary work on star formation-AGN mixing treats line emission in the region below the Kauffmann et al. (2003) line on the BPT diagram as pure star formation, and line emission in the region above the Kewley et al. (2001) line is considered to be dominated by harder sources, such as AGN or shocks. The region on the BPT diagram which occupies the space between the two lines is considered to result from a mixture of star formation with an additional hard component.

The study of star formation-AGN mixing was furthered through the advent of integral field spectroscopy (IFS). Pioneered by Davies et al. (2014a, b), a spatially-resolved spectrum across a galaxy containing an AGN shows a very clear and smooth ‘mixing sequence’ between star formation towards the outer of the galaxy, and AGN activity towards the centre. To calculate the relative contributions to several strong emission lines from star formation and AGN activity across the galaxy, Davies et al. (2014a, b) empirically selected two ‘basis points’ from the data on the BPT diagram to represent 100% star formation and 100% AGN activity. The relative contribution from star formation and AGN in each spaxel was then calculated by finding each spaxel’s ‘star-forming distance’, defined in Kewley et al. (2006). D’Agostino et al. (2018) furthered this technique by selecting five basis points through the use of photoionisation model grids for H ii regions and narrow-line regions (NLRs) using the mappings v photoionisation modelling code.

Despite improving upon the method showcased by Davies et al. (2014a, b), the theoretical calculation of the star formation-AGN fraction shown in D’Agostino et al. (2018) remains incomplete. Considering the line emission from a galaxy to be the result of two sources (star formation and AGN activity) is a vast oversimplification. The BPT diagram, which has been a favoured tool for previous work on star formation-AGN mixing, fails to distinguish between star formation, AGN activity, and other sources of ionisation and excitation. In particular, shocked outflows have been shown to result from starburst galaxies (e.g. Heckman et al., 1987; Rupke et al., 2005), AGN (e.g. Cecil et al., 2002; Rupke & Veilleux, 2011), and galaxy mergers (e.g. Rich et al., 2011, 2014; Rupke & Veilleux, 2013). Furthermore, shocks can result from accretion disks. Simulations and analytical calculations have shown shocks to be prevalent among the accretion disks of black holes (Spruit, 1987; Molteni et al., 1994; Sponholz & Molteni, 1994) and have been observed in the accretion disk of the dwarf nova U Geminorum (Neustroev & Borisov, 1998). Additionally, ionisation and excitation between an AGN and shocks can be further confused through differences in observed phenomena, such as in the AGN spectral energy distribution (SED) or continuum (e.g. Ferland & Mushotzky, 1982; Kraemer et al., 2009), or even gas turbulence (particularly when using [N ii]$\lambda 6584$ as a diagnostic; Gray & Scannapieco 2017). Yet on the BPT diagram, shocked line emission can be found in the same region as AGN-affected spaxels, or within the mixing sequence of the galaxy (for the position of shock features on diagnostic diagrams, see Rich et al., 2010, 2011; Kewley et al., 2013a). Hence, when studying a mixing sequence from a galaxy, simplifying the emission into a ratio between only star formation and AGN activity can lead to misleading if not erroneous results.

A three-component emission decomposition between star formation, shocks, and AGN activity was performed by Davies et al. (2017), by using the BPT diagram as well as the [O iii]/H$\beta$ vs [S ii]/H$\alpha$ diagnostic diagram. Similarly to Davies et al. (2014a, b), the basis points to represent 100% star formation, shock and AGN line emission were selected empirically. Despite utilising another diagnostic diagram, the [O iii]/H$\beta$ vs [S ii]/H$\alpha$ diagram is just as problematic as the BPT diagram when attempting to separate shocked line emission from star formation- and AGN-ionised line emission.

Throughout this work, we use a new three-dimensional diagnostic diagram first described in D’Agostino et al. (2019a), which shows the separation between star formation, shocks, and AGN line emission more clearly than the BPT diagram. Further, we follow D’Agostino et al. (2018) by using theoretical photoionisation modelling to define regions on the 3D diagnostic diagram of 100% star formation, shocks, and AGN. In Sections 2 and 3, we provide detail on the IFU survey data used in this paper, as well as providing information on our test case NGC 1068. The main work of this paper is found in Section 5, which includes demonstrating the 3D diagnostic diagram, as well as calculating the star formation, shock, and AGN fractions in NGC 1068. We compare our results to data in wavelengths outside of the optical range in Section 6, before summarising our findings and describing future applications of this method in Section 7.

2 Siding Spring Southern Seyfert Spectroscopic Snapshot Survey (S7)

The Siding Spring Southern Seyfert Spectroscopic Snapshot Survey (S7) is an optical IFU survey performed between 2013 and 2016. The Wide Field Spectrograph (WiFeS; Dopita et al., 2007; Dopita et al., 2010) was used for the S7 survey, located on the ANU 2.3m telescope at Siding Spring Observatory. WiFeS has a field-of-view of $38\times 25$ arcsec2 composed of a grid of $1\times 1$ arcsec2 spaxels. The S7 is an ideal choice for this work, due to the high spectral resolution and spatial resolution of the survey. In particular, the high spectral resolution ($R\sim 7000$ in the red, 3000 in the blue) allows the study of independent velocity components in the emission lines, providing more information on shocked regions (Ho et al., 2014). Full details of the S7, including wavelength ranges and selection criteria, can be found in Dopita et al. (2015a) and Thomas et al. (2017).

3 NGC 1068

NGC 1068 is the prototypical Seyfert 2 (Osterbrock & Martel, 1993) (R)SA(rs)b (de Vaucouleurs et al., 1991) galaxy, roughly located at a distance of 12.5 Mpc (e.g. Schoniger & Sofue, 1994) in the constellation Cetus. It has an optically-thick torus, which is known to obscure the broad-line region of the galaxy (e.g. Miller & Antonucci, 1983; Antonucci & Miller, 1985; Miller et al., 1991; García-Burillo et al., 2016). However, contemporary work has begun to resolve the circumnuclear disk and nuclear region of NGC 1068 (Marinucci et al., 2016; García-Burillo et al., 2016).

A large-scale biconical outflow with position angle ${\sim}32^{\mathrm{o}}$ is seen from the centre of the galaxy, observed to be moving outward at velocities of ${\sim}3000$ kms*-1* (Cecil et al., 2002). Such an outflow is believed to be the result of radiative acceleration from the central AGN (Pogge, 1988; Cecil et al., 2002; Dopita et al., 2002a; D’Agostino et al., 2018). This conclusion is supported by results at radio wavelengths. Emission from dense molecular gas tracers, in particular CO(3-2), CO(6-5), HCN(4-3), HCO+ (4-3), and CS(7-6), was mapped by García-Burillo et al. (2014) using ALMA. Using these maps, García-Burillo et al. (2014) show a large outflow in all molecular gas tracers in the inner 50 - 400 pc, with a mass of $M_{\mathrm{mol}}=2.7^{+0.9}_{-1.2}\times 10^{7}\;M_{\odot}$. A tight correlation between the motions in the circumnuclear disk (CND), the radio jet, and the ionised gas outflow is also shown by García-Burillo et al. (2014). Such a correlation suggests an AGN-driven outflow. This notion is further supported when considering the outflow rate in the circumnuclear disk. García-Burillo et al. (2014) state that the CND outflow rate of $63^{+21}_{-37}\;M_{\odot}$ yr*-1* is much larger than the star formation rate (SFR) calculated at the same radii in NGC 1068. Furthermore, Tacconi et al. (1994) conclude that gas at the leading edge of the bar in NGC 1068 is shocked, through observations in the near-infrared.

Young stellar clusters towards the centre of NGC 1068 have previously been identified by Storchi-Bergmann et al. (2012), supporting the idea of nuclear star formation. However, the nuclear star formation in NGC 1068 is small, with nuclear SFR estimates of $SFR_{\mathrm{nuclear}}\sim 0.4-1.0\;M_{\odot}$ yr*-1* being calculated between radii of 12 - 140 pc (Davies et al., 2007; Esquej et al., 2014). This is in stark contrast to the SFR calculated away from the nucleus, up to radii of several kpc. Surrounding the CND, Thronson et al. (1989) calculate a SFR of ${\sim}100\;M_{\odot}$ yr*-1*, likely a result of the high molecular mass found in the region (${\sim}2-6\times 10^{9}\;M_{\odot}$). Such a large mass is believed to have been confined to the region surrounding the CND by the bar at the centre of NGC 1068 (Thronson et al., 1989). Hence, the majority of star formation in NGC 1068 is seen in what is aptly named the ‘star-forming ring’ towards the central region of the galaxy. A corollary to such a statement is that extensive star formation towards the outskirts of NGC 1068 is very rarely seen.

4 Photoionisation models

We use photoionisation models produced with mappings v.1 (Sutherland et al., 2018) to aid in the separation of star formation, shock and AGN line emission. Models of H ii regions and NLRs are used to represent line emission from star formation, and shocks/AGN respectively. We use identical H ii region and NLR models to those described in D’Agostino et al. (2018), with the exception of using a NLR model pressure of $P/k=2\times 10^{7}\;\mathrm{cm}^{-3}\;\mathrm{K}$, corresponding to assumptions of an initial NLR temperature of 20,000 K and initial electron density of $n=1000\;\mathrm{cm}^{-3}$. An initial electron density of $n_{e}=1000\;\mathrm{cm}^{-3}$ for the NLR of NGC 1068 is justified through the electron density distribution shown in D’Agostino et al. (2018), calculated using the [S ii] doublet ratio. An initial temperature of 20,000 K is considered to be within the typical range for NLRs (e.g. Woltjer, 1959; Camenzind & Courvoisier, 1983; MacAlpine, 1986). We note however that an electron temperature of ${\sim}20,000$ K is calculated using the [O iii] ratio (e.g. Vaona et al., 2012), and the use of different emission line ratios (e.g. the [O ii] or [S ii] ratios) may lead to a broader range of possible calculated temperatures for the NLR (Taylor et al., 2003; Vaona et al., 2012). For a detailed explanation of the parameters associated with the H ii region and NLR models, see D’Agostino et al. (2018).

5 Separating line emission from star formation, shocks, and AGN

5.1 The need for a new diagnostic

Previous studies which considered separating line emission from star formation and an additional high energy component, such as shocks (e.g. Rich et al., 2010, 2011, 2014), or AGN (e.g. Kewley et al., 2001, 2006; Kewley et al., 2013a; Kewley et al., 2013b; Davies et al., 2014a, b; D’Agostino et al., 2018) utilised the BPT diagram most favourably, as well as other diagnostic diagrams such as the [O iii]/H$\beta$ vs [S ii]/H$\alpha$ and the [O iii]/H$\beta$ vs [O i]/H$\alpha$ diagnostic diagrams introduced by Veilleux & Osterbrock (1987). The radiation field produced by sources such as shocks or AGN is harder than that from star formation, leading to increased excitation of collisionally-excited emission lines, such as [N ii]$\lambda$6584, and [O iii]$\lambda$5007. Hence, separation of star formation and line emission from harder sources is straightforward on an emission line diagnostic diagram.

These diagnostic diagrams fail however, when attempting to separate line emission between shocks and AGN. Shocks and AGN can produce very similar values of the line ratios on the aformentioned diagnostic diagrams, leading to the two sources being indistinguishable if both shocks and AGN are present in a single galaxy. Furthermore, shocked spaxels on the BPT diagram in particular may coincide with spaxels along the galaxy’s mixing sequence (see Davies et al., 2014a, b; D’Agostino et al., 2018), leading to further confusion. Hence, in order to simulataneously separate line emission from star formation, shocks, and AGN activity in a single galaxy, a new diagnostic diagram must be explored.

5.2 The 3D diagnostic diagram

To simultaneously separate line emission from star formation, shocks, and AGN, we utilise the 3D diagnostic diagram first shown by D’Agostino et al. (2019a). We consider the 3D diagnostic diagram to be an extension of the BPT diagram, combining the BPT line ratios with distance and velocity dispersion. The information from the BPT diagram is retained by using a function of the two line ratios, defined as the emission-line-ratio (ELR) function. The ELR functional form is given in Equation 1. The ELR function is defined in such a way that any two hardness-sensitive line ratios may be used instead of those from the BPT diagram. We have chosen the BPT line ratios in particular, as the BPT diagram has been historically the most favoured diagnostic diagram for studying the separation of star formation, and AGN or shocks. Other diagnostic line ratios such as [O i]/H$\alpha$ or [S ii]/H$\alpha$ may be easily substituted into the ELR function. The purpose of the ELR function is to simply order the spaxels on the 3D diagram in terms of their combined [O iii]/H$\beta$ and [N ii]/H$\alpha$ ratios. We emphasise that the ELR function is data-dependent, and the range of values from 0 to 1 is arbitrary. The endpoints of 0 and 1 do not represent any physical phenomena, such as 100% star formation emission and 100% AGN emission respectively.

Prior to the data being displayed on the 3D diagnostic diagram, we first reduce and flux-calibrate the emission-line fluxes in the raw cubes using lzifu (Ho et al., 2016). lzifu is capable of fitting multiple Gaussian components to each emission line, thereby resolving multiple velocity components. Once each emission line is fit with a maximum of three Gaussian components, the neural network LZComp (Hampton et al., 2017) is used to determine the recommended number of Gaussian components for each emission line in each spaxel. Maps showing the distribution across the galaxy of each resolved velocity component in each spaxel are shown in Figure 1. Full details of the S7 data reduction are provided in Thomas et al. (2017).

An example of the 3D diagnostic diagram showing data from NGC 1068 is shown in Figure 2. The data points shown on the 3D diagram represent each individual velocity component in the data. For every component present in each spaxel, the velocity dispersion of the individual component (hereafter referred to as the ‘single-component velocity dispersion’) is combined with the total flux (‘zeroth’ component) and the radial value for the spaxel to define the position on the figure. Full details of the 3D diagram can be found in D’Agostino et al. (2019a). Data points which contain the first-component, second-component, and third-component fits to the velocity dispersion are shown on the 3D diagram in Figure 2 in red, blue, and green respectively. The first, second, and third components are numbered in order of the narrowest to broadest velocity dispersion. Hence, higher-order components will have greater values of the single-component velocity dispersion. We include the higher-order component fits to the velocity dispersion on the 3D diagram, as shocks are typically diagnosed and categorised by high velocity dispersions (Rich et al., 2010, 2011, 2014; Ho et al., 2014). All relevant emission lines are corrected for extinction, in accordance with Cardelli et al. (1989). The extinction calculation from Cardelli et al. (1989) uses the colour excess $E(B-V)$, which itself uses the Balmer decrement H$\alpha$/H$\beta$.

[TABLE]

5.2.1 Effects of beam smearing on the velocity dispersion

The finite beam size of a telescope inevitably leads to a phenomenon known as beam smearing (e.g. Bosma, 1978; Begeman, 1987). The smearing of line emission across adjacent spaxels produces broader line profiles, which ultimately can be erroneously interpreted as a higher velocity dispersion in a given spaxel. This could prove problematic, because the 3D diagram from D’Agostino et al. (2019a) uses the velocity dispersion in a spaxel as a major shock diagnostic. Large velocity dispersions towards the centre of an AGN are predominantly the result of outflows from the NLR (e.g. Freitas et al., 2018). Thus, the impact of beam smearing may be small, yet potentially significant nonetheless.

To verify that the large velocity dispersions seen in NGC 1068 are not solely the result of beam smearing, we use the code 3Dbarolo (Di Teodoro & Fraternali, 2015) to recover the rotational kinematics of the galaxy. The results of using 3Dbarolo are shown in Figure 3. The velocity field is fit using the H$\alpha$ line in the 3Dbarolo model. The model panels in Figure 3 (right-hand-side panels) show maps of the line-of-sight velocity field and velocity dispersion of the derived rotation curve of the galaxy. We are only concerned with the central region of the galaxy, because beam smearing is associated with strong velocity gradients. Strong velocity gradients are more likely found towards the centre of a galaxy due to disk rotation, where the velocity field will transition between positive and negative velocities rapidly. Figure 3 shows that in the ideal model case, beam smearing will only increase the velocity dispersion at the centre of NGC 1068 slightly (seen as the small region of ${\sim}100$ kms*-1* spaxels in the model velocity dispersion panel). The data however, shows a much larger region with very high velocity dispersions. This provides confidence that the high velocity dispersions at the centre of NGC 1068 are physical, and not the result of beam smearing. Furthermore, high spatial resolution images from Hubble Space Telescope (HST) (e.g. Cecil et al., 2002) show high-velocity regions (up to ${\sim}3000$ kms*-1*) within the central arcsecond of NGC 1068, supporting our claim of the physical nature of these high velocity dispersions.

5.2.2 Calculating the distance measurement in each spaxel

In order to assign each spaxel a distance value, we first check for local maxima in the [O iii]/H$\beta$ distribution across the galaxy. Local maxima in the [O iii]/H$\beta$ distribution are identified, and shown as red crosses in Figure 4a. Each maximum is treated as a centre, and a deprojected radius is calculated from each maximum, assuming the same values for the position angle and axis ratio for the galaxy. In the event where multiple maxima are detected (as in the case of NGC 1068), the deprojected radii calculated from each maximum are multiplied together, and then raised to the power of $1/n_{\mathrm{maxima}}$, where $n_{\mathrm{maxima}}$ is the number of maxima detected (the ‘geometric mean’). If only a single peak is detected in the [O iii]/H$\beta$ map, the peak is treated as the centre of the galaxy, and the distance measurement defaults to the deprojected galactocentric radius. A map showing the distribution of distance values across NGC 1068 is shown in Figure 4b. Hereafter, the distance measurement may simply be referred to as ‘distance’.

5.3 Defining the star formation, shock, and AGN extrema

The 100% regions in the 3D diagram indicate where emission is expected to be solely from one ionising source. Each of the 100% star formation, shock, and AGN regions of the 3D diagram are defined by a value of the ELR function, similarly to the spaxels on the BPT diagram. We calculate the 100% star formation, shock, and AGN ELR function values through photoionisation modelling using mappings v.1.

The 100% star formation and AGN emission line ratios for NGC 1068 are obtained by using the BPT diagram directly. We first construct new grid lines of constant metallicity on the H ii region and NLR grids, after calculating the metallicity gradient of the galaxy using the Kewley & Dopita (2002) (hereafter KD02) metallicity diagnostic. The [N ii]/[O ii] diagnostic described by KD02 is a good choice when attempting to calculate the metallicity gradient of galaxies which contain an AGN. Kewley & Ellison (2008) show that the relative contribution from an AGN has a very minor effect ($\leq 0.04$ dex at ${\sim}15$% AGN contribution) on the [N ii]/[O ii] ratio, in spaxels below the Kewley et al. (2001) theoretical maximum starburst line. The central and outer metallicities for NGC 1068 are shown in Table 1. The newly-constructed line of constant metallicity in the H ii region grid uses the outer metallicity of the galaxy, whilst the line constructed on the NLR grid uses the central metallicity. This is because the mixing sequence of the galaxy traces increasing AGN contribution from H ii regions outside of the AGN extended narrow-line region (ENLR) towards the AGN at the centre of the galaxy (see Kewley et al., 2001; Kewley et al., 2013a; Davies et al., 2014a; D’Agostino et al., 2018).

We select two basis points each on both the new H ii region and NLR grid lines, to account for the spread in ionisation parameter in the mixing sequence on the BPT diagram. The basis points are intended to signify the regions of 100% star formation and 100% AGN on the H ii region and NLR grids respectively. These basis points are located along the new lines of constant metallicity, and are shown in light blue on Figure 5. The mean value of the pair of basis points then determines the 100% ELR function value for the star-forming and AGN regions. The point representing the mean value for each pair is shown in light green on Figure 5. We use the 100% AGN ELR function value for the 100% shock ELR function value also. This is because shocks may be prevalent in the NLR of galaxies, produced from sources such as the radiation pressure from accretion disk emission, and the pressure from the radio jets (e.g. Dopita, 1995; Wilson & Raymond, 1999; Dopita et al., 2015b, and references therein). The values of the ELR function used to represent 100% star formation, AGN, and shock emission are given in Table 1.

5.4 Calculating the star formation-shock-AGN fraction

The positions of the spaxels on the 3D diagram are important in quantifying the contribution to emission from star formation, shocks, and AGN in each spaxel. Seen in Figure 2, the spread of spaxels on the 3D diagram shows two distinct and clear sequences. The first sequence contains spaxels at low velocity dispersions, and shows very little change in velocity dispersion with increasing values of the ELR function. The second sequence however displays a clear increase in velocity dispersion as the value of the ELR function increases in the spaxels. D’Agostino et al. (2019a) show that the first and second sequences represent the star formation-AGN and star formation-shock mixing respectively (see D’Agostino et al. 2019a for complete details on the diagram and its interpretation).

Each data point on the 3D diagnostic diagram is assigned a star formation-shock-AGN fraction. The star formation-shock-AGN fraction for each data point is calculated by considering its distance along the two sequences. Points containing high ELR function values, found towards the top of the first and second sequences, are assigned high AGN and shock fractions respectively. Conversely, points with low ELR function values will be assigned higher star formation fractions. The proximity of the data point to each sequence influences its final AGN and shock fractions. Depending on the position of the data points on the 3D diagnostic diagram, it is possible that any given data point will be assigned a fraction of 0% for one or more of the ionising sources.

5.4.1 Results

The results of calculating AGN, shock, and star formation fractions on the 3D diagnostic diagram for NGC 1068 are shown in Figures 6, 7, and 8, respectively. We also include 2D projections of the fractions in the distance-$\sigma$ and ELR function-$\sigma$ planes in Figures 6, 7, and 8. Unlike in Figure 2, we do not include the ELR function-distance plane projection, due to the obscuration of many spaxels preventing its usefulness. A line of best fit is applied to both the star formation-AGN and star formation-shock sequences of spaxels to show the direction of star formation-AGN and star formation-shock mixing respectively. The fractions are calculated for all spaxels on the 3D diagnostic diagram; some spaxels may at this point have multiple values for their star formation-shock-AGN fractions, depending on the number of velocity components present in the spaxel. As a result, for a spaxel with multiple velocity components, the final fractions are calculated as the weighted average of the fractions in the individual components, weighted by the contribution of each component to the total [O iii]$\lambda 5007$ flux. Maps of the final star formation, shock, and AGN fractions across the galaxy in the S7 field-of-view are shown in Figure 9. BPT diagrams of NGC 1068, with each spaxel colour-coded to the final star formation, shock, and AGN fractions, are shown in Figure 10. We have chosen [O iii]$\lambda 5007$ as the line by which we calculate the relative weight of each component, because it is a strong emission line associated with the presence of an AGN and shocks. Hence, it is a reliable emission line to use for the analysis in this work. However, other strong lines are present in the data, such as H$\alpha$ and [N ii]$\lambda 6584$. We calculate that the choice of which emission line is used to calculate the weighting can impact the final fractions by up to ${\sim}8$%. It is important to note that if two galaxies are to be compared equally, the emission line by which the final fractions are weighted must be identical. When calculating the relative contribution of each ionising source to the flux of strong emission lines in Table 2 (described below), the errors on the fractions account for the possible differences in value as a result of weighting by a different emission line.

The star-forming ring in NGC 1068 can be seen surrounding the nucleus in the star formation ratio map from Figure 9. The star formation ratio map also shows that ionisation in the nucleus is dominated by non-stellar sources. In particular, the biconical outflow most prominently seen in the [O iii]$\lambda$5007 is largely visible in the AGN ratio map. This agrees with previous work, claiming that the bicone is the result of material radiatively accelerated by the AGN (e.g. Pogge, 1988; Cecil et al., 2002; Dopita et al., 2002a; D’Agostino et al., 2018). The shock ratio map shows that shocks are seen to dominate the central few arcseconds of NGC 1068, however, due to the location of the high shock fractions, these shocks are likely from the AGN. An accreting supermassive black hole will produce thermal X-rays which provide a radiation pressure on nearby gas, causing expulsion of the gas in the form of a wind. This high-velocity wind will cause shocked material upon interaction with the ISM (Zubovas & King, 2012). Hence, the shocked structure seen in the nucleus of NGC 1068 is likely the result of the AGN.

The star formation, shock, and AGN contributions to the flux of several strong lines are shown in Table 2. For comparison, also in Table 2, we show the star formation and AGN contributions to the same strong emission lines, calculated using the star formation-AGN mixing method from D’Agostino et al. (2018). The method used to calculate the star formation, shock, and AGN contributions to each emission line listed in Table 2 is identical to the method described in Davies et al. (2014a) and D’Agostino et al. (2018):

The total luminosity of any emission line in an IFU field-of-view is given by:

[TABLE]

where $L_{i}$ is the luminosity of the emission line in spaxel $i$. The total luminosity of the emission line attributable to any source (star formation, shocks, or AGN) can be calculated by

[TABLE]

where $f_{i}^{\mathrm{source}}$ is the source fraction in spaxel $i$. It follows that the relative fraction of emission attributable to each source for the given emission line is provided by

[TABLE]

In addition to accounting for the choice of emission line used for the final weightings, the errors associated with each fraction given in Table 2 also consider the error in the photoionisation models described in Section 4. Assuming an error of ${\sim}0.1$ dex in both the [O iii]/H$\beta$ and [N ii]/H$\alpha$ ratios (e.g. Kewley et al., 2001; D’Agostino et al., 2019b), the resulting change in the position of the basis points on the 3D diagram can alter the final fractions by up to ${\sim}10$%. Finally, the errors on the fractions in Table 2 also consider the uncertainties in the emission line fits from lzifu. The error analysis performed on the star formation-AGN mixing fractions also in Table 2 is updated from that in D’Agostino et al. (2018), considering the 0.1 dex uncertainty in the photoionisation models.

The method from D’Agostino et al. (2018) assumes the emission in each spaxel can be expressed as a combination of star formation and AGN activity, without accounting for shock emission. The fractions in Table 2 show that if shocks are not considered as a possible source of ionisation, then the contribution to photoionisation from an AGN may be severely overestimated. The star formation fractions in each of the strong lines listed are similar between the two methods, implying that the majority of shock emission mixes with the emission from AGN photoionisation. This is further illustrated in Figure 11. We show a BPT diagram and map of NGC 1068 coloured by AGN fraction, assuming only a star formation-AGN mixing regime, against a BPT diagram and map coloured by the shock $+$ AGN fraction, calculated by adding the fractions seen in Figures 9 and 10. The results are very similar. This is unsurprising, given that shocks can produce emission line ratios similar to that of AGN. This is particularly true on the BPT diagram, where emission from shocks and AGN can be found within the same region of the diagram. Hence, using the BPT diagram to calculate star formation and AGN fractions, and only assuming a star formation-AGN mixing regime as in Davies et al. (2014a, b), Davies et al. (2016), and D’Agostino et al. (2018), will likely overestimate the true contribution from the AGN to emission line fluxes.

Davies et al. (2014a, b), Davies et al. (2016), and D’Agostino et al. (2018) all show that using H$\alpha$ as a SFR indicator in an AGN galaxy may lead to an overestimate of the true SFR, as the AGN is responsible for a fraction of the total H$\alpha$ luminosity. This is also seen in Table 2, as the AGN contribution to the H$\alpha$ flux in NGC 1068 is non-zero, irrespective of the inclusion of shocks to the decomposition. However, once shocks are accounted for in the decomposition, the star formation fractions of each of the strong lines, while similar, are systematically lower. Therefore, it is evident that emission from star formation does not solely mix with AGN, but rather multiple sources of ionisation and excitation. In order to calculate a SFR as accurately as possible using the H$\alpha$ flux, as many sources of ionisation and excitation as possible must be simultaneously separated. Failure to do so will result in the SFR consistently being overestimated. Using the SFR(H$\alpha$) relation given by Kennicutt et al. (1994), we calculate a SFR of 3.1 $M_{\odot}$ yr*-1* in the S7 field-of-view (radius ${\sim}1$ kpc), after calculating the H$\alpha$ luminosity from star formation to be $3.9\times 10^{41}$ erg s*-1*.

6 Comparison with other data

To verify our star formation, shock, and AGN fractions in NGC 1068, we use the structure seen in other wavelengths. X-ray emission is known to be ubiquitous in AGN (e.g. Gandhi, 2005), and nuclear shocks may be caused by the relativistic jets from the SMBH, which are very visible in the radio part of the spectrum (Blandford & Königl, 1979). We also study the CO(3-2) distribution in the nucleus of NGC 1068 to verify our star formation results. García-Burillo et al. (2014) show that the majority of CO(3-2) flux measured by ALMA (${\sim}63$% of the total) is detected in the star-forming ring surrounding the nucleus.

Shown in Figure 12 are maps of the [O iii]$\lambda 5007$ luminosity from the AGN, and the 0.25-7.5 keV Chandra X-ray distribution from Young et al. (2001). Contours of each map are shown on the other. The X-ray contours accurately trace the regions of high luminosity in the AGN map. Similarly, we find large agreement with the distribution of the [O iii]$\lambda 5007$ shock luminosity, and the 3mm radio continuum contours from García-Burillo et al. (2017) in Figure 13. The 3mm radio continuum data from García-Burillo et al. (2017) clearly shows the distribution and position of the radio jets. Finally, in Figure 14 we show the CO(3-2) contours from García-Burillo et al. (2014) over the H$\alpha$ luminosity distribution from star formation, and the star formation ratio map from Figure 9. The CO(3-2) contours accurately trace regions of high star formation surrounding the nucleus, in agreement with the findings of García-Burillo et al. (2014).

We also use of the He ii$\lambda 4686$/H$\beta$ distribution map to verify our results. He ii is a high-ionisation line (ionisation potential of ${\sim}54$ eV), thus requiring a hard radiation field (the kind associated with AGN) for ionisation. The line ratio He ii/H$\beta$ is not contained in our emission line ratio function shown in Equation 1, and so provides an independent test of our results. The He ii/H$\beta$ map from the S7 data is shown in Figure 16, normalised to the maximum value. Regions of high He ii/H$\beta$ can be seen to align with regions of high AGN fraction in Figure 9, as expected.

To further verify our shock results, we consider the distribution of the [O iii]$\lambda 4363$/[O iii]$\lambda 5007$ ratio across the galaxy. The [O iii] ratio is a temperature-sensitive line ratio, and increases in the line ratio indicate an increase in the electron temperature. Shocked gas can reach temperatures on the order of $10^{6}$ K (for a fully-ionised plasma; e.g. Dopita & Sutherland, 2003), compared to the NLR [O iii] temperatures of ${\sim}20,000$ K discussed in Section 4. Hence, an increase in the [O iii] line ratio is expected to coincide with regions of high shock fractions in Figure 9.

The [O iii] line ratio map is shown in Figure 16. Only spaxels which contain significant detection of both the [O iii]$\lambda 4363$ and [O iii]$\lambda 5007$ emission lines are shown in Figure 16, corresponding to a signal-to-noise greater than 3. Such spaxels are concentrated towards the centre, where the [O iii] luminosity is high, seen for example in Figures 12 and 13. The variation in the [O iii] ratio amongst included spaxels is low, and hence the values in Figure 16 have been enhanced $\times 100$ to better highlight the variation. Clearly, an increase in the [O iii] ratio is seen towards the top-left of the included spaxels. Similarly, the shock fraction map in Figure 9 also shows a slight increase in the shock fractions in spaxels towards the top-left of centre. Although, it should be noted that in this region, the AGN fractions of the spaxels also seen in Figure 9 are considerably high (${\sim}40$%). Furthermore, photoionisation modelling with mappings shows that AGN photoionisation can lead to an increase in the [O iii] ratio, in the range from that expected from fast shocks of ${\sim}0.1-0.01$ (e.g. Dopita et al., 2002b; Allen et al., 2008). The variation seen in the [O iii] ratio map is better aligned with the features in the shock + AGN map in Figure 11. Thus, we conclude that the distribution of the [O iii] ratio does not show sufficient variation to make confident assumptions about the shock-to-AGN ratio.

The shock-dominated regions are primarily located close to the nucleus of NGC 1068. Thus, possible mechanisms for their creation include NLR outflows, or the associated radio structures. We suggest that the shocks are the result of NLR outflows. The radio jet structure seen in HST imaging from Cecil et al. (2002) shows a radio structure at the nucleus of NGC 1068 with an extent of ${\sim}1-2$ arcseconds. This is far smaller than the extent of the shock-dominated spaxels, which can extend beyond 10 arcseconds, seen in Figure 9. In addition, our Figure 1, and Figure 4 from D’Agostino et al. (2019a) show the extent of the high-$\sigma$ spaxels to be far beyond ${\sim}1-2$ arcseconds. Hence, the small-scale radio jet must not be the dominant shock-production mechanism.

7 Conclusions and future work

We have demonstrated a new method to simultaneously separate line emission from gas ionised by star formation, shocks, and AGN in IFU data. Our 3D diagnostic diagram shows two clear mixing sequences, indicating star formation-AGN and star formation-shock mixing within a single galaxy. Using the information in these two clear mixing sequences, we have quantified the ratio of star formation-, shock-, and AGN-ionised line emission in each spaxel of the Seyfert galaxy NGC 1068.

We have shown that if shocks are not considered as an ionising source, then the relative contribution from an AGN to emission line fluxes may be greatly overestimated. The shock + AGN fractions for each emission line are very similar to the AGN fractions calculated assuming only a star formation-AGN mixing regime, using the method from D’Agostino et al. (2018). This indicates that shocked line emission mixes heavily with the AGN-ionised line emission. Including shocks as a possible ionising source also systematically lowers the star-forming fractions calculated in each emission line. If the flux of H$\alpha$ is to be used as a SFR indicator, it is recommended to separate as many ionsing sources as possible.

We have compared our results to data in various wavelengths. These comparisons are summarised:

(i) The luminosity distribution of AGN-ionised line emission closely resembles the 0.25-7.5 keV X-ray photon map from Young et al. (2001).

(ii) The luminosity distribution of shock-ionised line emission is accurately aligned with the radio jets in NGC 1068, seen in the 3mm radio continuum map from García-Burillo et al. (2017).

(iii) The luminosity distribution of star formation-ionised line emission closely traces the CO(3-2) molecular line flux from García-Burillo et al. (2014). The star-forming ring in NGC 1068 is host to the majority (${\sim}63$%) of the CO(3-2) flux in the ALMA field-of-view (García-Burillo et al., 2014).

(iv) Regions of AGN-dominated line emission correlate very accurately with regions of high He ii/H$\beta$ ratios. He ii$\lambda 4686$ is a high-ionisation line, and thus requires a hard radiation field (such as that from an AGN) to produce large emission line fluxes.

(v) Regions of shock-dominated line emission show a slight correlation with regions of increased [O iii]$\lambda 4363$/[O iii]$\lambda 5007$, although the AGN contribution in these regions is also considerable. This indicates that the shock and AGN separation may need further advancements.

In a future publication, we aim to use our star formation-shock-AGN separation method to quantify the star formation, shock, and AGN fractions of merging galaxies, at various merger stages.

Acknowledgements

Parts of this research were conducted by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013.

Support for AMM is provided by NASA through Hubble Fellowship grant #HST-HF2-51377 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555.

The authors wish to also acknowledge the contribution from the reviewer, who facilitated great discussion amongst the co-authors and thus greatly improved the paper.