The mass-metallicity relation of local active galaxies

TL;DR

This study measures the gas-phase metallicities and the mass-metallicity relation of local active galaxies using SDSS data, incorporating advanced modeling of emission-line mixing and physical conditions, revealing a systematic offset compared to star-forming galaxies.

Contribution

It introduces a novel method for analyzing active galaxy metallicities by accounting for mixing between emission regions and varying physical parameters, providing new insights into their chemical properties.

Findings

Active galactic nuclei show a metallicity increase with host galaxy mass.

A 0.09 dex offset exists in metallicity between active and star-forming galaxies.

The study highlights potential modeling deficiencies affecting metallicity estimates.

Abstract

We systematically measure the gas-phase metallicities and the mass-metallicity relation of a large sample of local active galaxies for the first time. Observed emission-line fluxes from the Sloan Digital Sky Survey (SDSS) are compared to a four-dimensional grid of photoionization models using the Bayesian parameter estimation code NebulaBayes. For the first time we take into account arbitrary mixing between HII region and narrow-line region (NLR) emission, and the models are also varied with metallicity, ionization parameter in the NLR, and the gas pressure. The active galactic nucleus (AGN) oxygen abundance is found to increase by $\Delta {\rm O/H} \sim 0.1$ dex as a function of host galaxy stellar mass over the range $10.1 < \log M_* / M_\odot < 11.3$. We also measure the metallicity and ionization parameter of 231000 star-forming galaxies for comparison with the sample of 7670 Seyfert 2 galaxies. A systematic offset in oxygen abundance of 0.09 dex is observed between the mass-metallicity relations of the star-forming and active galaxies. We investigate potential causes of the offset, including sample selection and the treatment in the models of diffuse ionized gas, pressure, and ionization parameter. We cannot identify the major cause(s), but suspect contributions due to deficiencies in modeling the ionizing spectra and the treatment of dust physics. Optical diagnostic diagrams are presented with the star-forming and Seyfert data colored by the inferred oxygen abundance, ionization parameter and gas pressure, clearly illustrating the trends in these quantities.

Tables

Table 1. Table 1: The measured mass-metallicity relation parameters γ 𝛾 \gamma , M 0 and θ 0 subscript 𝜃 0 \theta_{0} for different choices in the analysis, along with θ 0 subscript 𝜃 0 \theta_{0} offsets between corresponding star-forming galaxy and AGN mass-metallicity relations.

#	Typea	Lines	$\log U_{\mathrm{HII}}$	$\log P_{\mathrm{HII}}/k^c$	redshift, $z$		EW(H$\beta$)		Low-M⋆ slope, $\gamma$		Turnover, M0		Asymptotic O/H, ${\theta }_0$
		& line		(cm-3 K)			(Å)				($\log$ M${}_{\star }{}^{}/$M⊙)		(12 + $\log$ O/H)
		ratiosb			min	max	min	max	value	error	value	error	value	error	offsetd
1e	SF	A	free	free	-	-	-	-	0.62	0.03	8.76	0.02	8.94	0.009	-
2	AGN	B	-3.50	$P_{\mathrm{NLR}}$	-	-	-	-	”	”	”	”	9.01	0.002	0.07
3e	AGN	B	-3.25	$P_{\mathrm{NLR}}$	-	-	-	-	”	”	”	”	9.03	0.002	0.09
4	AGN	B	-3.00	$P_{\mathrm{NLR}}$	-	-	-	-	”	”	”	”	9.05	0.003	0.11
5	AGN	B	-3.25	5.80	-	-	-	-	”	”	”	”	9.02	0.002	0.08
6	AGN	B	-3.25	6.20	-	-	-	-	”	”	”	”	9.03	0.002	0.09
7	AGN	B	-3.25	6.60	-	-	-	-	”	”	”	”	9.04	0.002	0.10
8f	AGN	Af	-3.25	$P_{\mathrm{NLR}}$	-	-	-	-	”	”	”	”	9.00	0.003	0.06
9	SF	A	free	free	0.000	0.050	-	-	0.63	0.02	8.82	0.02	8.99	0.009	-
10	AGN	B	-3.25	$P_{\mathrm{NLR}}$	0.000	0.050	-	-	”	”	”	”	9.06	0.006	0.07
11	SF	A	free	free	0.050	0.075	-	-	0.54	0.04	8.85	0.04	8.99	0.02	-
12	AGN	B	-3.25	$P_{\mathrm{NLR}}$	0.050	0.075	-	-	”	”	”	”	9.08	0.005	0.08
13	SF	A	free	free	0.075	0.100	-	-	0.68	0.05	8.98	0.03	8.95	0.01	-
14	AGN	B	-3.25	$P_{\mathrm{NLR}}$	0.075	0.100	-	-	”	”	”	”	9.05	0.006	0.10
15	SF	A	free	free	0.100	0.200	-	-	0.76	0.05	9.09	0.02	8.93	0.01	-
16	AGN	B	-3.25	$P_{\mathrm{NLR}}$	0.100	0.200	-	-	”	”	”	”	9.03	0.005	0.10
17	SF	A	free	free	-	-	-30	-10	0.48	0.03	8.81	0.04	8.98	0.02	-
18	AGN	B	-3.25	$P_{\mathrm{NLR}}$	-	-	-30	-10	”	”	”	”	9.08	0.01	0.10
19	SF	A	free	free	-	-	-10	-6	0.66	0.03	8.78	0.02	8.93	0.007	-
20	AGN	B	-3.25	$P_{\mathrm{NLR}}$	-	-	-10	-6	”	”	”	”	9.04	0.005	0.11
21	SF	A	free	free	-	-	-6	-3	0.68	0.02	8.76	0.01	8.93	0.006	-
22	AGN	B	-3.25	$P_{\mathrm{NLR}}$	-	-	-6	-3	”	”	”	”	9.03	0.004	0.10
23	SF	A	free	free	-	-	-3	0	0.88	0.03	8.71	0.02	8.88	0.004	-
24	AGN	B	-3.25	$P_{\mathrm{NLR}}$	-	-	-3	0	”	”	”	”	8.97	0.004	0.09

Equations

$$\theta (M_{\star })={\theta }_0-\log \left[1+{\left(\frac{M_{\star }}{M_0}\right)}^{-\gamma }\right]$$

Full text

Galaxies: abundances, Galaxies: active, Galaxies: emission lines, Galaxies: ISM, Galaxies: Seyfert

1 Introduction

The gas-phase metallicity of a galaxy is the result of enrichment by stellar processes and the exchange of material with its environment (e.g. Lilly et al., 2013). The assembled stellar mass of a galaxy is known to be strongly positively correlated with its metallicity (e.g. Lequeux et al., 1979; Tremonti et al., 2004).

This mass-metallicity relation (or, similarly, the luminosity-metallicity relation; e.g. Zaritsky et al., 1994) has been measured for star-forming galaxies both locally (e.g. Tremonti et al., 2004) and at a range of redshifts (e.g. Zahid et al., 2014a; Sanders et al., 2015; Wuyts et al., 2016). For galaxies hosting active galactic nuclei (AGN), Matsuoka et al. (2018) recently found a positive correlation between metallicity and host galaxy stellar mass at high redshift. However, to date the mass-metallicity relation has not been cleanly measured for active galaxies in the local universe.

Metallicity measurements in AGN have typically focused on the broad-line regions (BLRs) of Seyfert AGN and quasars (Hamann & Ferland, 1999), although narrow-line region (NLR) abundances have also received some attention. Inferred quasar abundances are ubiquitously high and show negligible evolution with redshift when considering either the BLRs (e.g. Dietrich et al., 2003; Nagao et al., 2006c; Juarez et al., 2009) or the NLRs (Nagao et al., 2006a; Dors et al., 2014). Quasar luminosity is correlated with the metallicity measured in both the BLR (Nagao et al., 2006c) and in the NLR (Nagao et al., 2006a). There is a strong correlation between metallicity-sensitive NLR and BLR line ratios in the same objects (Du et al., 2014).

Optical emission-line metallicity diagnostics for NLRs include those of Storchi-Bergmann et al. (1998) and Castro et al. (2017). Metallicity diagnostics for NLRs are reviewed by Dors et al. (2015). We note that independent estimates of nuclear metallicities in active galaxies may be obtained from the metallicities measured in nearby H II regions and by radial extrapolation of H II region metallicities from the galactic disk (e.g. Storchi-Bergmann et al., 1998; Dopita et al., 2014).

The mixing of H II region and NLR excitation in Seyfert spectra is a key complication in the use of metallicity diagnostics. The analysis in this study compares observations to ‘mixed’ models that account for varying contributions from H II region and NLR emission. Here, using a statistically powerful sample from the Sloan Digital Sky Survey (SDSS; York et al., 2000) and accounting for H II-AGN mixing, we are able to investigate the AGN mass-metallicity relation in detail for the first time.

In Paper 1 (Thomas et al. 2018, accepted to ApJL) we studied the degree of mixing between star-forming and Seyfert narrow-line region (NLR) emission in the SDSS Seyfert spectra. In this study we present further results, including the parameter estimates for the metallicity, pressure and ionization parameter. We additionally expand our analysis to the SDSS star-forming galaxies to contextualize the results of the mixing analysis.

2 Method

Our approach uses the analysis tool NebulaBayes (Thomas et al., 2018), which performs Bayesian parameter estimation by comparing observed emission-line fluxes and errors with grids of theoretical fluxes.

2.1 Observational data

We use a similar sample to the sample used in Paper 1. The parent sample is SDSS DR7 (Abazajian et al., 2009), and we use the emission line fluxes and stellar masses measured by the MPA-JHU group111www.mpa-garching.mpg.de/SDSS/DR7/

http://home.strw.leidenuniv.nl/~jarle/SDSS/ (Tremonti et al., 2004). We apply the empirical scalings to the line flux errors that were derived by that group using duplicate observations. The following cuts were applied to the sample:

We require a signal-to-noise ratio of 5 for the H$\beta$ line flux. This cut is only applied to H$\beta$ because cuts on other lines may introduce biases – in particular, a cut on [O III] may bias against high-metallicity objects with weak [O III]. 2. 2.

We include objects having slightly negative fluxes for other lines, to allow for errors in the continuum subtraction and line fitting – for a flux $F_{i}$ with error $E_{i}$, we require that $F_{i}>-2.5E_{i}$. This cut is applied separately in the star-forming galaxy and Seyfert galaxy analyses on all lines of the chosen set in each case (Section 2.3). 3. 3.

We take only the first of any duplicate observations 4. 4.

Galaxies with an unphysical Balmer decrement are excluded, because this is indicative of errors in calibration and continuum subtraction. We require $F_{\rm H\alpha}/F_{\rm H\beta}>2.7$. 5. 5.

We apply standard cuts on the [N II] and [S II] optical diagnostic diagrams (Kewley et al., 2006). We include only ‘star-forming galaxy’ and ‘Seyfert’ classifications and deliberately exclude LINERs, composite galaxies and ‘ambiguous’ galaxies that have conflicting classifications on the two diagrams. Negative fluxes are set to zero for the purposes of the classification.

We did not apply redshift cuts, but systematically consider the effect of different redshift selections in Section 4. After the above cuts there remained 231429 spectra classified as star-forming and 7669 spectra classified as Seyfert 2. We explore the effects of sample selection in Section 4.

2.2 Model data

We use MAPPINGS V (Sutherland & Dopita, 2017) photoionization models for both the H II regions and NLRs. All models are plane-parallel, one-dimensional, and dusty, with elemental depletions onto dust grains based on $\log({\rm Fe}_{\rm free}/{\rm Fe}_{\rm total})=-1.5$ (Jenkins, 2009, 2014). Total elemental abundances are set using the oxygen-based standard abundance scaling of Nicholls et al. (2017) and the ‘local galactic concordance’ reference abundances, in which the Solar oxygen abundance is $12+\log\,{\rm O/H}=8.76$. The total pressure in a model is equal to the sum of the gas and radiation pressure such that the total pressure increases as radiation is absorbed in each model step.

For our analysis of the Seyfert galaxies we use a ‘mixing’ grid that combines the photoionization models of H II regions and NLRs. The mixing grid has four parameters in total: oxygen abundance $12+\log{\rm O/H}$, NLR ionization parameter $\log U_{\rm NLR}$, gas pressure $\log P/k$, and the fraction of the Balmer flux arising in the NLR as opposed to H II regions, $f_{\rm NLR}$, with the ionizing Seyfert spectrum identical in all models.

Assumptions are employed to reduce the number of free model parameters in the mixing analysis; limiting the number of parameters is necessary because of computational practicalities and the limited number of independent measured line fluxes in each spectrum. The H II and Seyfert nebulae are assumed to have the same metallicity and gas pressure. The ionization parameter in the H II-region component is fixed to $\log U_{\rm HII}=-3.25$, a representative value for high-metallicity galaxies. In Section 4 we investigate the sensitivity of the results to the treatment of gas pressure and to the choice of fixed $\log U_{\rm HII}$.

The Seyfert model components use an ionizing spectrum from Thomas et al. (2016) that parameterizes the energy of the ionizing accretion disk emission by its peak energy, $E_{\rm peak}$. In Paper 1 we concluded that values of $E_{\rm peak}$ in the range $40-50$ eV result in plausible distributions of $f_{\rm NLR}$ measurements. In the present work we assume E_{\rm peak}$$\,=45 eV.

For our analysis of star-forming galaxies, the H II-region MAPPINGS grid runs over the three parameters of oxygen abundance, ionization parameter, and gas pressure. The grid is identical to that described by Thomas et al. (2018). A diffuse ionized gas (DIG) component is not included in the models.

2.3 Application of NebulaBayes

In both the H II and mixing analyses we use the Balmer lines H$\alpha$ and H$\beta$ for reddening correction and normalization. We use the ability of NebulaBayes to deredden spectra to match the computed Balmer decrement at every point in the model grid. The dereddening method is described in the appendix of Vogt et al. (2013).

For the calculation of the likelihood in the mixing analysis we use the emission lines [O II]$\,\lambda\lambda\,3726,$3729, [Ne III]$\,\lambda\,3869$, [O III]$\,\lambda\,4363$, H$\beta$, [O III]$\,\lambda\,5007$, He I$\,\lambda\,5876$, [O I]$\,\lambda\,6300$, H$\alpha$, [N II]$\,\lambda\,6583$, and [S II]$\,\lambda\lambda\,6716,$ 6731. These lines include strong optical lines as well as higher-excitation lines that are characteristic of Seyfert nuclei. The fluxes are normalized to H$\beta$ for comparison with the models in the calculation of the likelihood, but we know that other line ratios may effectively constrain parameters of interest. Hence, we use the NebulaBayes ‘line ratio prior’ feature (described in the appendix of Thomas et al., 2018) to apply priors using two diagnostic ratios. These are [N II]/[O II] (a sensitive metallicity diagnostic for metallicites above half Solar; Kewley & Dopita, 2002), which slightly reduces the scatter in the inferred O/H, and [S II] $\lambda 6731$ / [S II] $\lambda 6716$, an electron density diagnostic (e.g. Dopita & Sutherland, 2003; Osterbrock & Ferland, 2006) that improves our constraints on the gas pressure. The priors are calculated by comparing observed and model ratios over the entire parameter space, similarly to the calculation of the likelihood.

To fit the star-forming galaxies, we relied only on diagnostic line ratios by using the NebulaBayes ‘line ratio prior’ as the likelihood (instead of the usual method of likelihood calculation). In particular, we used the [N II]/[O II] ratio (N2O2) to constrain O/H, the [O III]/[O II] ratio to constrain the ionization parameter, and the [S II] $\lambda 6731$ / [S II] $\lambda 6716$ ratio to constrain the pressure, with the N2O2 constraint weighted three times more heavily than each of the other two constraints. This approach recognizes the sensitivity of the N2O2 ratio to metallicity, and was necessary to achieve a smooth metallicity distribution over the entire range in O/H of more than a decade. This analysis used a uniform prior (in the logarithmic space of each parameter).

It is not meaningful to compare a model emission line flux to a negative measured flux, so we set any negative line flux $f_{i}$ with error $e_{i}$ to $f_{i}=+0.05e_{i}$ in the NebulaBayes analysis. This produces effectively the same constraint as $f_{i}=0$ but satisfies the requirement of the code for positive input fluxes.

NebulaBayes allows the inclusion of a relative error on the model grid fluxes, for which we used a value of 0.35 (as a proportion of each model flux; see Section 2.2 of Thomas et al. 2018). This ‘grid error’ is used in both likelihood and prior calculations. An error on predicted fluxes of 35% may often be larger than errors on observed fluxes, but this value is necessary to ‘smooth out’ the constraints and produce a smooth metallicity distribution for the star-forming galaxies. We explore why the large ‘grid error’ is necessary in Section 4.3. For the star-forming galaxy analysis we used cubic interpolation of the model grids.

3 Results

3.1 Parameter estimates

Figure 1 presents optical diagnostic diagrams (Baldwin et al., 1981; Veilleux & Osterbrock, 1987) showing the SDSS data colored by the parameter estimates.

There is systematic variation in measurements of the ionization parameter $U$ in both the AGN and star-forming sequences. The gradient for Seyferts is perpendicular to the mixing line, with lower $U$ inferred for data closer to the LINER region (a region associated with low $U$ in the case of photoionized sources, e.g. Binette et al., 1994). The anticorrelation of $U$ with metallicity in star-forming galaxies has been previously noted (e.g. Nagao et al., 2006b; Dopita et al., 2006).

Figure 2 presents the parameter estimates for the Seyfert spectra as a function of projected distance along the mixing sequence, $d$. Here the mixing sequence (defined in Paper 1) is a straight line on the log-log optical diagnostic diagram designed to follow the upper section of the locus of data that connects the H II and AGN regions. The panel for $f_{\rm NLR}$ is the same as the panel for E_{\rm peak}$$\,=45 eV in Figure 2 of Paper 1.

From Figure 2 we note the following:

Metallicity: The high $12+\log{\rm O/H}$ measurements (medians of $1.4-1.7$ solar) are consistent with the ‘mixing’ with the high-metallicity end of the H II-region sequence (Figure 1). The ${\sim}0.1$ dex increase in median metallicities with $d$ is presumably related to a systematic offset between inferred H II-region and NLR metallicities (Section 4.1). 2. 2.

Mixing fraction: The $f_{\rm NLR}$ estimates show that the Seyfert-classified SDSS spectra are highly contaminated by H II-region emission, and $f_{\rm NLR}$ increases smoothly along the mixing sequence (the $f_{\rm NLR}$ results are the subject of Paper 1). 3. 3.

Ionization parameter: The measured $U_{\rm NLR}$ is typically in the tight range $-3.3<\log U_{\rm NLR}<-2.7$, and shows little variation as a function of $d$. The gradient in $\log U_{\rm NLR}$ in Figure 1 is not evident in Figure 2 because the gradient is approximately orthogonal to the mixing line. The lower inferred $U$ at low $f_{\rm NLR}$ (low $d$) may be due to LINER or DIG contamination. 4. 4.

Pressure: The median inferred $\log P/k$ values increase by approximately 1.2 dex along the mixing sequence.

The gas pressure is measured to be higher in galaxies with emission dominated by Seyfert excitation than in galaxies of similar masses that are dominated by star-forming excitation. This is evident in the differences between the top and bottom ends of the mixing sequence in Figure 1 and Figure 2. The higher pressures in NLRs are presumably caused by the effects of radiation pressure (e.g. Dopita et al., 2002; Groves et al., 2004) and potentially also the ram pressure of AGN-induced outflows. The single pressure parameter in our mixing analysis is effectively a ‘mean pressure’, which naturally increases as $f_{\rm NLR}$ increases and higher-pressure NLRs dominate the emission.

3.2 Mass-metallicity relations

Figure 3 presents the mass-metallicity relation resulting from our oxygen abundance measurements for both the H II- and Seyfert-classified spectra. Our derived AGN mass-metallicity relation is similar to that of star-forming galaxies, albeit covering only the high-mass end of the relation.

We fit the following functional form (Moustakas et al., 2011; Zahid et al., 2014b) to the mass-metallicity relations:

[TABLE]

where $\theta=12+\log{\rm O/H}$ is the oxygen abundance and $M_{\star}$ is the galaxy stellar mass. The three parameters are $\gamma$, the power-law slope of the relation at low stellar masses; $M_{0}$, a characteristic stellar mass at which the mass-metallicity relation begins to flatten; and $\theta_{0}$, the asymptotic oxygen abundance at high stellar masses. Equation 1 increases monotonically with $M_{\star}$.

Our approach is to firstly fit the mass-metallicity sequence of the star-forming galaxies, and then apply the resulting values of $\gamma$ and $M_{0}$ in the fit to the AGN sequence, such that only $\theta_{0}$ may vary. This approach is necessary because the AGN sequence does not extend to sufficiently low stellar masses to effectively constrain $\gamma$ and $M_{0}$. The fits were performed over the ranges $8.5<\log M_{*}/M_{\odot}<11.3$ and $10.1<\log M_{*}/M_{\odot}<11.3$ for the star-forming and active galaxies, respectively.

Our derived AGN mass-metallicity relation is similar in shape to the relation for star-forming galaxies, albeit covering only the high-mass end of the relation. We detect an increase in AGN metallicity with galaxy mass, with the median oxygen abundance increasing by ${\sim}0.1$ dex over the fitted range of 1.3 dex in stellar mass. This is consistent with the observation of Groves et al. (2006) that [N II] emission generally increases with the stellar mass of SDSS active galaxies, interpreted by the authors as an increase in metallicity with mass.

The scatter about the derived mass-metallicity relations, characterized by the median absolute residual, is 0.07 dex in O/H for both mass-metallicity relations. The NebulaBayes errors on individual metallicity measurements are generally more than three times larger than the residuals; these errors effectively include both statistical and systematic unceratinties because of the 35% error used on the grid fluxes (Section 2.3; Section 4.3).

3.3 Comparison to the literature

Figure 4 compares our results to some other calculated mass-metallicity relations. We apply the recent N2O2 metallicity calibration of Castro et al. (2017) to our Seyfert sample and include the resulting relation in the figure. We also show a selection of relations calculated by Kewley & Ellison (2008) using various diagnostics and a similar star-forming galaxy sample to our own. The well-known systematic offsets between different metallicity calibrations (e.g. Kewley & Ellison, 2008) are apparent. We compare our results to those of Tremonti et al. (2004) in more detail in Section 4.2.

The mass-metallicity relation we derive using the Castro et al. (2017) diagnostic is very similar to that produced by our NebulaBayes analysis in its shape, scatter, and absolute oxygen abundance values (after correcting for the 0.07 dex difference in assumed solar oxygen abundance between the Castro et al. (2017) models and our own). The Castro et al. (2017) photoionization models neglect mixing between H II-region and NLR emission, as well as variations in pressure. Despite these simplifications, the Castro et al. (2017) diagnostic has produced an almost identical Seyfert mass-metallicity relation to our more sophisticated analysis, which suggests that N2O2 is a very robust metallicity diagnostic, both for Seyfert galaxies and in general.

4 Further analysis and discussion

We confirm the results of Groves et al. (2006) and others that sub-solar AGN are rare in the SDSS sample. Low-metallicity AGN should lie between the ‘wings’ of the BPT ‘butterfly diagram’ (e.g. Kawasaki et al., 2017); in our sample objects in this region appear to have a combination of high ionization parameter and relatively low metallicity (Figure 1).

4.1 What causes the systematic offset between star-forming and active galaxies?

There is a systematic offset in oxygen abundance of $0.09$ dex between the star-forming and active galaxies in Figure 3. We do not expect that the gas-phase metallicities of star-forming galaxies and AGN should be different for samples of galaxies at the same mass, although there is at least one mechanism that speculatively may produce abundance enhancements in active nuclei – star formation with a top-heavy IMF in the accretion disk (e.g. Nayakshin & Sunyaev, 2005; Collin & Zahn, 2008; Bartko et al., 2010). In this section we explore potential methodological causes of the offset by reanalyzing the data in a series of experiments. The H II region and NLR models used the same photoionization code (MAPPINGS V) and abundances, so these could not have contributed to the offset. The reanalyses involved varying the lines and mixing grids used by NebulaBayes, as well as varying the selection of the sample. The experiments are summarized in Table 1.

4.1.1 HII region ionization parameter

The ionization parameter of the H II region component of the mixed models was set to $\log U_{\rm HII}=-3.25$, a value chosen to be representative of high-metallicity galaxies. The choice of fixed $U_{\rm HII}$ should affect the inferred AGN metallicities because measurements of ionization parameter and metallicity are interdependent, because $U$ and O/H are observed to be anticorrelated (Figure 1), and because the H II-region emission contribution is significant for most SDSS AGN.

Rows 2 and 4 of Table 1 present the results of varying $U_{\rm HII}$ by 0.25 dex below and above the assumed value in row 3. The inferred AGN metallicities (and hence $\theta_{0}$ estimates) are found to increase with $U_{\rm HII}$. It appears that the $\theta_{0}$ offset could possibly be eliminated by setting an unphysically low value of $U_{\rm HII}$. Nevertheless, we have determined that the choice of the fixed value of $U_{\rm HII}$ within the plausible range of values is not a significant contributor to the $\theta_{0}$ offset.

4.1.2 Treatment of pressure

In the mixing analysis only one parameter was used to account for varying gas pressure between galaxies. This approach assumes that the H II region and NLR pressures are equal and the parameter therefore measures an ‘effective’ pressure.

To test the sensitivity to this particular choice of pressure treatment, we performed the experiments presented in rows 5, 6 and 7 of Table 1. In these experiments we fixed the pressure of the HII-region contribution to the mixed models, such that the free pressure parameter was for the gas pressure in the NLR contribution only. Inspection of the pressure results for these experiments showed that the distributions of inferred NLR pressures were implausible (not smooth and with an implausibly large range). The mass-metallicity relation offsets remain of a similar size with this different pressure treatment. We therefore conclude that the treatment of pressure in the models is unlikely to make a significant contribution to the offset.

4.1.3 Choices of lines, line ratios and weightings

The most important choices for an analysis with NebulaBayes are in the selection of emission lines, emission line ratios and corresponding weightings for the comparison between models and observations. The star-forming galaxy and mixing analyses were very different in these choices (Section 2.3). In particular, the star-forming galaxy analysis relied purely on diagnostic line ratios and the mixing analysis took advantage of high-excitation lines that are only observed in AGN.

Row 8 of Table 1 presents the results of an experiment in which our preferred analysis (row 3) was repeated, but using the lines, line ratios and weightings that were applied in the star-forming galaxy analysis. This approach produced a reduction in the $\theta_{0}$ offset from 0.09 to 0.06, indicating that the treatment of lines, line ratios and weightings does substantially contribute to the $\theta_{0}$ offset.

4.1.4 Redshift cuts and aperture effects

It is possible that aperture effects and changes in sample properties with redshift may have systematically affected our analysis. To explore these effects, we considered subsamples with different ranges of redshift. Rows 9 to 16 of Table 1 present the results of fitting the star-forming galaxy and AGN mass-metallicity relations for galaxies in four different redshift bins. The subsamples have $\theta_{0}$ offsets of 0.07-0.10 dex, with the offset increasing by ${\sim}0.03$ as a function of redshift over the four bins. The variation in the offset suggests that aperture and selection effects could make a small contribution to the systematic offset in our preferred analysis, but cannot explain the majority of the offset.

4.1.5 Omission of DIG in the models

Emission from DIG is ubiquitous and contributes ${\sim}10-70\%$ of the H$\alpha$ luminosity in star-forming galaxies, with later (earlier) types showing contributions at the lower (higher) end of this range (e.g. Lacerda et al., 2018). Emission from DIG is associated with lower pressure, lower ionization parameter, and harder ionising spectra than H II-region emission, and DIG contamination may cause systematic offsets in metallicity measurements (Zhang et al., 2017; Pilyugin et al., 2018). We expect the mixing analysis to be less sensitive to DIG contamination because NLRs in Seyfert-classified objects are likely far more luminous than DIG, and the NLR models include DIG-like emission from the partially ionized zone.

To explore potential systematic effects of DIG contamination on the analysis, we considered subsamples defined to have different ranges of the equivalent width (EW) of the H$\beta$ emission line. The Balmer line EW is used as an indicator of DIG emission, with relatively small (large) EW associated with higher (lower) DIG contributions to the line emission (e.g. Cid Fernandes et al., 2011; Belfiore et al., 2016; Lacerda et al., 2018). Rows 17 to 24 of Table 1 present the results of reanalysing the data in each H$\beta$ EW bin. The four subsamples have $\theta_{0}$ offsets of $0.09-0.11$ dex, and the offset does not vary suggestively as a function of the H$\beta$ EW of the bin. We are therefore confident that the treatment of DIG in the models and DIG contamination of the SDSS spectra do not significantly contribute to the offset.

4.1.6 Treatment of dust

Dust destruction is not included in our models but does occur in AGN NLRs, as evidenced by observations of the strength of the coronal lines in Seyfert spectra (e.g. Dopita et al., 2015). Dust destruction in NLRs would lead to more free species contributing to line cooling than in models at the same metallicity, resulting in overestimation of metallicities. It is possible that dust destruction in Seyferts may lead to systematically different reddening compared to star-forming galaxies, but this does not appear to be the case: Seyferts have uncorrelated $f_{\rm NLR}$ and $A_{v}$ (Pearson $r=-0.07$); $f_{\rm NLR}$ and O/H are also uncorrelated ($r=0.08$). However, issues of geometry make it difficult to use extinction data to draw conclusions about varying dust destruction between NLRs and H II regions. Firstly, dust screens may be significantly physically separated from the relevant nebulae. Secondly, the observed extinction depends on the relative geometry of sources and dust (e.g. Calzetti et al., 1994), with NLRs having strongly differing geometries and extents to H II regions.

4.1.7 The ionizing spectra

Ionizing stellar spectra are uncertain (e.g. Morisset et al., 2004) because of model sensitivity to properties such as stellar rotation, helium dredge-up, mass-loss rates, and binarity, and additional issues with stellar track and stellar atmosphere models such as poor parameter-space coverage, interpolation, and inconsistent abundance sets. The ionizing spectra in NLRs are also uncertain, although our $f_{\rm NLR}$ parameter may partially compensate for variations in the ionizing Seyfert spectrum. If the offset primarily arises due to issues with the assumed ionizing spectra, it is unclear whether the problem is due to the stellar spectra, AGN spectra, or a combination of the two.

4.1.8 Summary

The major cause of the $0.09$ dex systematic offset in measured oxygen abundance between the star-forming and active galaxies has not been identified. Our analyses suggest that the offset is not primarily due to the treatment of pressure or the H II-region ionization parameter in the mixing analysis, to a sample selection effect, or to contamination of the observed spectra by DIG. We determined that a non-negligible contribution to the offset (${\sim}0.3$ dex) was made by the differences between the choices of emission lines in the two analyses. It remains unclear what causes the remainder of the offset.

The unexplained ${\sim}0.06$ dex offset may be due to real physical differences between the observed H II regions and NLRs that we do not account for in the models, including issues of geometry and the treatment of dust. We also suspect that the highly uncertain ionizing spectra make significant contributions to the offset.

4.2 Comparison of method with Tremonti et al. (2004)

We compare our results to those of Tremonti et al. (2004), whose mass-metallicity relation is a classic result in the literature and who followed a methodology similar to our own (Brinchmann et al., 2004). The shape of our star-forming galaxy mass-metallicity relation resembles that derived by Tremonti et al. (2004) (Figure 4; Tremonti et al. (2004) fitted a cubic polynomial). The scatter is also similar, with Tremonti et al. (2004) reporting a 1$\sigma$ spread in O/H of approximately $\pm 0.1$ dex ($\pm$0.07 dex at high galaxy masses), which compares to a median absolute residual of $\pm$0.07 dex for both our star-forming and active galaxy relations. The similarity in the scatter is consistent with the similarities between the methods of Tremonti et al. (2004) and our methods.

There is a systematic O/H offset between the star-forming galaxy mass-metallicity relation of Tremonti et al. (2004) and the relation derived in this work (Figure 4), which we suspect is due predominantly to differences in the modeling. Tremonti et al. (2004) used theoretical models by Charlot & Longhetti (2001) that had been calibrated to observed line ratios in nearby galaxies, whereas our models are purely theoretical. Another difference is that we fitted diagnostic line ratios in our star-forming galaxy analysis, whereas Tremonti et al. (2004) used line fluxes instead of flux ratios (similarly to our AGN analysis). However, the choice of lines, ratios and priors did not contribute most of the offset between our star-forming and active galaxy relations (Section 4.1.3), so we suspect the larger offset considered here is also mostly independent of these choices.

4.3 Evaluation of H II-region models

We have compared the observed and best-fit model fluxes for the star-forming galaxies in bins of measured O/H. The comparison gives insight into why the star-forming galaxy analysis required a different approach to the Seyfert analysis (Section 2.3).

Considering measurements of a given line for all galaxies in an O/H bin, we find typical offsets of ${\sim}20-30\%$ between the median dereddened observed flux and median predicted flux. The offsets are similar to the $0.1-0.15$ dex quoted by Kewley & Ellison (2008) as a common estimate for model uncertainties. However, an exception is [O III]$\lambda 5007$, which is systematically underpredicted by a factor of ${\sim}3-5$ at supersolar metallicities. It is likely that the difficulties fitting the H II-region models were primarily attributable to this [O III] offset. The offsets may have many causes (such as those considered in Section 4.1), but the underprediction of [O III] in particular suggests that the model ionising stellar spectra may be inaccurate (Section 4.1.7).

In general, the systematic offsets explain why it was necessary to use the NebulaBayes ‘grid error’ input to obtain a smooth star-forming galaxy metallicity distribution (Section 2.3). A consequence of using this parameter is that the derived uncertainties on individual parameter estimates are effectively a combination of both statistical uncertainties and model-related systematic uncertainties, rather than being solely statistical uncertainties. We note that the grid error also allows for systematic errors in the observations, e.g. due to difficulties with continuum subtraction.

5 Conclusions

We have used a Bayesian method to analyze the line emission in SDSS AGN and star-forming galaxies, in particular making systematic measurements of the gas-phase metallicities. The method accounts for mixing between H II-region and NLR emission in AGN for the first time.

We conclude that the inferred AGN metallicities increase with the galaxy stellar mass, as is the case for star-forming galaxies. The SDSS AGN are located in high-metallicity and high-mass galaxies, and are rare or difficult to detect in lower-mass galaxies.

The shape of the AGN mass-metallicity relation is consistent with the shape of the high-mass end of the star-forming galaxy relation. Both our star-forming and active galaxy relations have a median absolute residual (scatter) of 0.07 dex in oxygen abundance. The zero points of the mass-metallicity relations differ, however, with the inferred oxygen abundances for Seyfert galaxes being systematically 0.09 dex higher than for star-forming galaxies. We have investigated some potential causes of this offset, and although the difference remains largely unexplained we suspect that the model ionizing spectra make a significant contribution. Nevertheless, we conclude that our method produces robust relative AGN metallicities.

We look forward to applying similar analyses across large samples of AGN with spectroscopic data that is spatially-resolved, beginning with data from the S7 survey (Dopita et al., 2015; Thomas et al., 2017).

This research was conducted by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. B.G. acknowledges the support of the Australian Research Council as the recipient of a Future Fellowship (FT140101202). M.D. and R.S. acknowledge support from ARC discovery project #DP160103631. This research is supported by an Australian Government Research Training Program (RTP) Scholarship.