Cosmic metallicity evolution of Active Galactic Nuclei: Implications for optical diagnostic diagrams

We analyze the validity of optical diagnostic diagrams relying on emission-lines ratios and in the context of classifying Active Galactic Nuclei (AGNs) according to the cosmic metallicity evolution in the redshift range 0<z<11.2. In this regard, we fit the results of chemical evolution models (CEMs) to the radial gradients of the N/O abundances ratio derived through direct estimates of electron temperatures (Te-method) in a sample of four local spiral galaxies. This approach allows us to select representative CEMs and extrapolate the radial gradients to the nuclear regions of the galaxies in our sample, inferring in this way the central N/O and O/H abundances. The nuclear abundance predictions for theoretical galaxies from the selected CEMs, at distinct evolutionary stages, are used as input parameters in AGN photoionization models built with the Cloudy code. We found that standard BPT diagnostic diagrams are able to classify AGNs with oxygen abundances 12+logO/H>8.0 [(Z/Zsolar)>0.2) preferably found at redshift z>4. On the other hand, the HeII4685/Hbeta versus [N II]6584/Halpha diagram produces a reliable AGN classification independent of the evolutionary stage of these objects.


INTRODUCTION
Spectral classification of galaxies based on emission line intensities is fundamental for the understanding of the physical processes actuating in the Interstellar Medium (ISM) as well as for the knowledge of the ionizing source, of elemental abundances and of the metallicity of the gas phase of these objects.Baldwin et al. (1981), in their pioneering study, showed that combinations of strong optical emission line intensities can be used to separate objects according to their main excitation mechanism: H ii regions or Star Forming galaxies (hereafter SFs), Planetary Nebulae, Active Galactic Nuclei (AGNs), and objects excited by shock-waves.These authors proposed the [O iii]5007/H versus [N ii]6584/H diagnostic diagram to distinguish SFs from other object classes.Hereafter, this diagram is named as [N ii]-diagram.Veilleux & Osterbrock (1987) revised the method of classification proposed by Baldwin et al. (1981), including the [S ii](6716 + 6731)/H and [O i]6300/H line ratios, commonly known as BPT diagrams.The separation lines between SFs and AGNs in BPT diagrams, called as maximum starburst lines, have been mainly established by using observational data of objects in the local universe ( < 1, e.g.Kauffmann et al. 2003) and relied on photoionization model results (e.g.Kewley et al. 2001).However, optical spectroscopic data of SFs at ★ E-mail: olidors@univap.brhigh redshift ( > 1) have shown an offset of the maximum starburst lines from the local star-forming ones in the [N ii]-diagram, in the sense that SFs at higher redshifts can show emission-line intensity ratios commonly presented by AGNs (e.g.Shapley et al. 2005Shapley et al. , 2015;;Liu et al. 2008;Hainline et al. 2009;Bian et al. 2010;Steidel et al. 2014;Masters et al. 2014;Sanders et al. 2016a;Strom et al. 2017;Kashino et al. 2017;Cameron et al. 2023;Simmonds et al. 2023).
The cause of this offset is under discussion and it has been attributed to higher ionization parameters, harder radiation fields, higher N/O ratios in high redshift SFs and/or due to a bias in sample selection (see Sanders et al. 2023a;Garg et al. 2022).For instance, Kewley et al. (2013) compared a large sample of galaxies in the redshift range 0.5 <  < 2.6 with predictions of photoionization models in the [N ii]-diagram.These authors concluded that the offset is due to the ISM conditions are more extreme at high redshift than those of local galaxies (see also Sugahara et al. 2022).Otherwise, Sanders et al. (2016b), using observations of confirmed SFs (identified by X-Ray or infrared properties) at  ∼ 2.3 taken from the MOSFIRE Deep Evolution Field survey (Kriek et al. 2015), found that the cause of the offset of high- SFs in the [N ii]-diagram is due to the elevated N/O abundance ratio at fixed O/H (see also Shapley et al. 2015;Masters et al. 2014;Pérez-Montero & Contini 2009).Moreover, a combination of metallicity and ionization parameter effects can be the cause of the maximum starburst line displacement (Cowie et al. 2016;Garg et al. 2022;Papovich et al. 2022).
Regarding AGNs, this object class can occupy the SF zone in the [N ii]-diagram even in the local universe, yielding an incorrect classification of galaxy nuclei.Groves et al. (2006), who used photoionization models to examine the effects of metallicity () variations on the narrow emission lines of AGNs, showed that AGNs with (/Z ⊙ ) ≲ 0.5 can occupy similar regions than SFs with high ionization [log([O iii]5007/H) ≳ 0.5].In the same direction, Feltre et al. (2016) compared SF and AGN photoionization model results with observational data of a large sample of objects ( < 0.2) taken from Sloan Digital Sky Survey Data Release 7 (SDSS-DR7; Abazajian et al. 2009) in order to analyze the discrimination by the use of lineratio diagnostics.These authors found that for the low metallicity regime [(/Z ⊙ ) ≲ 0.4], AGN and SF models predict similar values of [O iii]5007/H and [N ii]6584/H, on the side of the [N ii]diagram corresponding to SFs (see also Feltre et al. 2023).Thus, despite there is evidence that only SFs at high- are located above the maximum starburst line, AGNs in the local universe and with low metallicity can occupy the SF zone in BPT diagrams (see also Zhu et al. 2023).This fact produces a bias in the AGN selection and in statistic metallicity studies because active galaxies with low metallicity, misclassified as SFs, are excluded from any selection/analysis.In fact, Osorio-Clavĳo et al. (2023) used optical data taken from the Calar Alto Legacy Integral Field Area survey (CALIFA, Sánchez et al. 2012) and data from the CHANDRA X-ray satellite (Weisskopf et al. 2002) to study the AGN population in the local universe.These authors showed that the AGN population in the CALIFA survey is 4 per cent (Lacerda et al. 2020) when only BPT diagrams are used as a selection criterion.However, this fraction could rise up to ∼ 10 per cent when accounting for both optical and X-ray spectroscopic analyses.Also, Bykov et al. (2023) showed that most part of AGNs hosted in low mass galaxies ( < 10 9.5 M ⊙ ) and selected by using X-ray luminosity are located in the SF zone in the [N ii]-diagram.Moreover, Dors et al. (2020) selected 463 AGNs ( < 0.4) based on BPT diagrams using the SDSS-DR7 (Abazajian et al. 2009) database and classifications from NASA/IPAC Extragalactic Database (NED).These authors estimated the metallicity through several strong-line methods and found that AGNs with (/Z ⊙ ) < 0.6 are only ∼ 10 per cent of the sample.Certainly, if standard BPT diagrams exclude low metallicity AGNs the fraction above is underestimated.
Important insights on the applicability of BPT diagrams at high- are coming from sophisticated analysis relying on hydrodynamic simulations combined with photoionization models.For instance, Garg et al. (2022), by using the simba code (Davé et al. 2019), simulated galaxies contained in a box with a side length of 100 ℎ −1 Mpc.The abundances resulting from simba simulations are assumed as input in SF photoionization models.From this simulation, Garg et al. (2022) showed that SFs at  = 2 either with an increasing of the N/O abundance ratio or with a decreasing ionization parameter at fixed O/H can go through the maximum starburst line in the [N ii]-diagram.Hirschmann et al. (2022) combined cosmological IllustrisTNG 1 simulations (Nelson et al. 2018) with SF and AGN photoionization models, gas ionization by post-AGB stars as well as fast, radiative shocks in order to analyze the validity of diagnostic diagrams.These authors showed that standard BPT diagrams break down for metalpoor ( < ∼ 0.5 Z ⊙ ) AGNs at  > ∼ 1 (see also Hirschmann et al. 2019;Nakajima & Maiolino 2022), hence AGN galaxies move towards the SF zone.
As pointed out by Hirschmann et al. (2022), a source of the inaccuracy of current large-scale cosmological simulations (see also Buck et al. 2021), such as the IllustrisTNG and simba simulations, is that they do not resolve the ISM, i.e. these simulations neither resolve gas properties, such as density in individual ionized regions, nor track the formation, evolution and destruction of dust grains.This problem makes that observed AGN luminosity function, stellar populations and associated scaling relations are not always well reproduced by these models.To circumvent the possible limitation of cosmological hydrodynamic simulations, detailed fitting of chemical evolution models to physical parameter estimations (e.g.star formation rate, abundances, etc) of individual galaxies can produce more reliable constraints to the applicability of the BPT diagrams.For instance, Mollá & Díaz (2005) presented a generalization of the multiphase chemical evolution model (hereafter CEM) applied to a wide set of theoretical galaxies with different masses and evolutionary rates.Over decades, this classical approach of CEMs have been able to predict the radial elemental abundance distributions (with a spatial resolution of few kpc) of galaxies (e.g.Cavichia et al. 2023;Sharda et al. 2021;Mollá et al. 2019;Magrini et al. 2016;Kang et al. 2016;Kubryk et al. 2015;Cavichia et al. 2014;Marcon-Uchida et al. 2010;Schönrich & Binney 2009;Chiappini et al. 2001;Prantzos & Boissier 2000;Diaz & Tosi 1984;Talbot 1980;Alloin et al. 1979).
Another approach of most simulations built to analyze the validity of diagnostic diagrams along the cosmic time is the assumption of the O/H abundance as a metallicity tracer of the ISM (e.g.Garg et al. 2022).However, for AGN photoionization modeling this approach can introduce some uncertainty in the O/H abundance or metallicity values and, consequently, in the predicted emission line intensities used in the BPT diagrams.Radio observations have shown the existence of neutral/molecular gas reservoir in the central parts of galaxies containing AGNs (e.g.Dressel et al. 1982;Hutchings et al. 1987;Bertram et al. 2007;Ho et al. 2008;García-Burillo et al. 2014;Bradford et al. 2018;Ellison et al. 2019;Combes et al. 2019). Recently, do Nascimento et al. (2022), by using MaNGA spectroscopic data of 108 Seyfert nuclei, showed that this neutral/molecular gas reservoir seems to be the cause of the lower (0.16-0.30dex) O/H abundance in AGNs than the expected value from central extrapolation of radial metallicity gradients.Likewise, elemental abundance estimates (e.g.O/H, N/H) can be somewhat uncertain due to several factors.For example, O/H abundances based on direct estimates of electron temperature ( e -method) can be underestimated (up to ∼ 0.2 dex) due to electron temperature fluctuations, possibly present in SFs (e.g.Méndez-Delgado et al. 2023;Oliveira et al. 2008;Peimbert et al. 2007;Hägele et al. 2006;Krabbe & Copetti 2002) and in AGNs (Riffel et al. 2021).On the other hand, it is largely known that photoionization models trend to overestimate elemental abundances (by 0.1-0.4dex) in relation to those via the  e -method (e.g.Kennicutt et al. 2003;Dors & Copetti 2005;Kewley & Ellison 2008).Thus, chemical evolution model predictions and elemental abundance estimations of galaxies involving hydrogen content (e.g.O/H, N/H) can be somewhat uncertain.In this sense, the abundance ratio of heavy elements (e.g.N/O, C/O, Fe/O) trend to be more reliable metallicity indicators in comparison to O/H abundance, since they do not take into account the hydrogen content, as well as the effects of electron temperature fluctuations and the uncertainties in photoionization models trend to be minimized.In fact, Watanabe et al. (2023) compared the values of the abundance ratios (Fe/O, Ar/O, S/O) predicted by chemical evolution models with direct abundance estimates for galaxies in the local universe and in the range of  ∼ 1−10.Important results arose from this comparison as, for instance, winds of rotating Wolf Rayet stars that end up as a direct collapse could explain the high N/O abundance [log(N/O) > −0.5] derived for the galaxy GN-z11 ( ∼ 10.6, Cameron et al. 2023;Senchyna et al. 2023) rather than supernovae ISM enrichment (but see Maiolino et al. 2023 for a different result).
In particular, the N/O abundance ratio has been suggested as a good metallicity tracer and it is correlated to the galaxy mass (e.g.Hayden-Pawson et al. 2022;Florido et al. 2022;Kojima et al. 2017;Douglass & Vogeley 2017;Pérez-Montero et al. 2016;Masters et al. 2016;Andrews & Martini 2013;Pérez-Montero & Contini 2009) as well as it is a useful tool to study the interplay of galactic processes, for instance, star formation efficiency, the timescale of infall, and outflow (e.g.Johnson et al. 2023;Magrini et al. 2018).In addition, recent estimates of radial abundance gradients in local spiral galaxies via direct estimations of the electron temperatures ( e -method 2 ), for example, make available by the CHAOS project 3 (e.g.Rogers et al. 2022Rogers et al. , 2021;;Croxall et al. 2016Croxall et al. , 2015;;Berg et al. 2015) as well as estimates for very high redshift ( > 5) galaxies obtained by the James Webb Space Telescope (e.g.Curti et al. 2023b;Arellano-Córdova et al. 2022;Hsiao et al. 2023) provide reliable constraints to chemical evolution models and, consequently, to the validity of the BPT diagrams in the cosmological context.
Motivated by the large number of direct elemental abundance estimates in galaxies with a wide range of redshift and by the advantage of the use of N/O as a metallicity tracer, in this work, we used predictions of the chemical evolution models proposed by Mollá & Díaz (2005) to reproduce direct estimations of N/O gradients in four local spiral galaxies, whose data were taken from the literature.Extrapolations of the radial gradients to the central parts of the galaxies, predicted by the selected models by Mollá & Díaz (2005), were carried out in order to obtain the central N/O abundances, in a wide range of redshifts ( = 0 − 11).Furthermore, the resulting N/O and the corresponding O/H abundances, for different redshift bins, were assumed as input parameters in AGN photoionization models built with the Cloudy code (Ferland et al. 2013).This methodology makes it possible to investigate the feasibility of optical diagnostic diagrams for the AGN classification in a wide range of redshifts.The structure of this paper is as follows.In Section 2 the methodology employed in the analysis, i.e. observational abundance estimates used as constraints, CEMs and photoionization models descriptions, is presented.In Sect. 3 the results and their discussion are presented.Conclusions are presented in Sect. 4. Throughout this paper we adopt the Planck Collaboration et al. ( 2021) cosmologic parameters: H 0 = 67.4km s −1 Mpc −1 and Ω m = 0.315.

METHODOLOGY
We obtained optical emission line intensities predicted by photoionization models in order to analyze the reliability of the BPT diagrams for a wide redshift range.For that, we used predicted N/O radial abundances from CEMs proposed by Mollá & Díaz (2005) to reproduce the ones derived through the  e -method in a small sample of local spiral galaxies.After that, the N/O gradient resulting for each theoretical galaxy with distinct evolutionary stages was extrapolated to its central part (galactocentric distance equal to zero) in order to obtain its nuclear abundance.The N/O nuclear abundances and their corresponding O/H abundances were used as input parameters in AGN photoionization models.Finally, the optical emission-line intensities ratios predicted by the photoionization models were employed in order to investigate the galaxy classification for the redshift range 2 For a review of the  e -method see Peimbert et al. (2017) and Pérez-Montero (2017). 3https://www.danielleaberg.com/chaos = 0.0 − 11.0.In the following sections, we present each employed procedure.

Direct radial abundance estimates
We compiled from the literature the N/O radial abundances of the disks of four local ( < 0.1) spiral galaxies derived through direct measurements of the electron temperatures, i.e. by using the  emethod, known as the most reliable method (see e.g.Pilyugin 2003;Hägele et al. 2006Hägele et al. , 2008Hägele et al. , 2012;;Toribio San Cipriano et al. 2017).Since the  e -method requires the presence of faint (about 100 times fainter than H) auroral lines (e.g.[O iii]4363; Díaz et al. 2007), direct estimates of N/O radial gradients are only available for a few spiral galaxies and do not cover the entire optical disks.In fact, only for 20 over the 41 disk H ii regions belonging to the grand design and nearby spiral galaxy M 101 (NGC 5457) observed by Kennicutt & Garnett (1996), it was possible to measure electron temperatures (see Kennicutt et al. 2003).Likewise, more recent spectroscopic data of spiral galaxies obtained by the CHAOS project (e.g.Berg et al. 2015), a project dedicated to measuring direct abundance in disk H ii regions using the Large Binocular Telescope, have shown the difficulty to measure electron temperatures even in nearby galactic disks.Thus, to assume only abundance estimates via the  e -method as a selection criterion penalizes any study, in the sense that few objects are able to be selected (e.g.Dors et al. 2020).
To estimate the radial N/O abundance distributions of spiral galaxies, it is usually assumed a linear regression: where  is the galactocentric distance (in units of kpc), log(N/O) 0 is the extrapolated value of this parameter to the galactic centre ( = 0 kpc) and grad is the slope (in units of dex kpc −1 ).To derive the N/O galactic gradient, H ii regions distributed across the galaxy disk are necessary (e.g.Gusev et al. 2012).Therefore, we selected from the literature only galaxies for which N/O estimates via the  e -method were obtained for more than 10 disk H ii regions with galactocentric distances  reaching up to the isophotal radius  25 (see Pilyugin et al. 2004), defined as the galactocentric distance where the surface brightness is 25 mag arcsec −2 in the -band.The selected objects, the number of disk H ii regions and the authors who calculated the N/O estimates are: NGC 628 (45, Berg et al. 2015), NGC 2403(28, Rogers et al. 2021), M 101 (20, Kennicutt et al. 2003), and M 33 (60, Rogers et al. 2022).
The following methodology was applied by the authors (see Table 1) to estimate the direct abundances: • Two electron temperatures [ e (low) and  e (high)] are derived through the [N ii](6548 + 6584)/5755 and [O iii](4959 + 5007)/4363 line intensities ratios for the zones were the O + -N + and the O +2 ions are originated, respectively (Peimbert 1967).For the cases when it was not possible to apply the previously mentioned procedure to estimate  e (low), it was estimated through theoretical/empirical relations with the  e (high) (see e.g.Hägele et al. 2008;Rogers et al. 2022).
It is worth mentioning that as the N/O estimates for the objects in our sample were taken from different authors, for which distinct observational techniques, reddening law and atomic parameters were assumed, the estimates are obtained adopting not a homogeneous methodology.This could produce artificial scattering or biases in the Figure 1.Comparison between N/O abundance ratio estimates based on the  e -method (black points) and predictions from chemical evolution models by blue points 2005MNRAS.358..521M for the spiral galaxies indicated on each panel.Solid black lines represent linear regressions to the estimates via the  e -method while dashed black lines the uncertainties of these fittings.Blue lines represent linear regressions to the chemical evolution model results.Table 1 lists the main features of the considered galaxies and the references for direct abundance estimates.derived abundances and, consequently, in the CEM interpretations.Effects of the use of distinct methodologies and heterogeneous samples on nebular abundance analysis are discussed, for instance, by Juan de Dios & Rodríguez (2017) and Dors et al. (2023), which could introduce abundance variations in the order of ∼ 0.1 dex, which is similar to the uncertainties introduced by the errors in the optical emission-line measurements (see e.g.Hägele et al. 2008).
In Figure 1, the logarithm of the direct N/O abundances taken from the literature as a function of the galactocentric distances are plotted together with their linear regressions obtained following the Equation 1 for each studied galaxy.The fitting uncertainties are also shown.
In Table 1, the identification of the spiral galaxies considered in our analysis, the logarithm of their masses (in units of the solar mass), the number of disk H ii regions with direct abundance estimates, the galactocentric distance  range in which the N/O values were derived, and their literature references are presented.Zaritsky et al. (1994) showed the existence of correlations between the O/H gradients and galactic macroscopic quantities (e.g.mass, absolute magnitude, Hubble type), being these correlations dependent on the normalization to the galactocentric distance (see also Vila-Costas & Edmunds 1992;Pilyugin et al. 2004Pilyugin et al. , 2019;;Ho et al. 2015).In particular, less massive galaxies trend to have steeper gradients (in units of dex kpc −1 ) in comparison to those in more massive objects (e.g.Ho et al. 2015).Despite this result is still under debate in the literature (see e.g.Belfiore et al. 2017;Poetrodjojo et al. 2018), any study aiming to obtain representative abundance gradients in spiral galaxies must consider a wide galactic mass interval (e.g.Hirschmann et al. 2022).The mass interval of our sample is 9.6 < ∼ log(/M ⊙ ) < ∼ 11, covering practically all the mass range considered by Maiolino et al. (2008), who derived the mass-metallicity relation for galaxies in a wide redshift range.Even though abundance values and metallicity gradients for galaxies with mass values lower than log(/M ⊙ ) < ∼ 9.0 are available in the literature (e.g.Ho et al. 2015;Bresolin 2019;Curti et al. 2023a), N/O gradients via the  e -method in these objects are barely found.
In addition, galaxies located in dense environments (e.g.Dors & Copetti 2006;Mouhcine et al. 2007;Robertson et al. 2012;Lara-López et al. 2022) as well as interacting galaxies (e.g.Krabbe et al. 2008Krabbe et al. , 2011;;Kewley et al. 2010;Rupke et al. 2010;Sánchez et al. 2014;Rosa et al. 2014;Torres-Flores et al. 2014;Zinchenko et al. 2015) trend to have different metallicities and, consequently, abun-dance gradients in comparison to the isolated ones.Moreover, the shape of oxygen abundance profiles, in some cases, can not be represented by a single function (e.g.Martin & Roy 1995;Roy & Walsh 1997;Bresolin et al. 2009Bresolin et al. , 2012;;Marino et al. 2012;Pilyugin et al. 2017;Sánchez-Menguiano et al. 2018), as the one proposed in the Eq. 1.Our sample of spiral galaxies do not present clear signals of interactions with near galaxies, as reported in the works from which the direct estimates were taken.Likewise, as we can see in Fig. 1, for 3/4 objects, the N/O abundance distributions are well represented by a linear fitting, with slopes steeper than −0.02 dex kpc −1 , not indicating any effect of environment/interaction on the abundances.This feature is important because the CEMs proposed by Mollá & Díaz (2005) consider only the cosmic evolution of isolated galaxies.
For 3/4 galaxies in our sample, clear O/H and N/H single linear regressions with negative slopes were derived in the works from which the estimates were obtained (not shown), and thus, the N/O regressions shown in Fig. 1 indicate that the nitrogen has mainly a secondary source of production, i.e. (N) ∼ (O) 2 (e.g.Hamann & Ferland 1993), being  the abundance.It can be also seen in Fig. 1, that for the H ii regions belonging to M 101 and located at galactocentric distances larger than ∼ 20 kpc, a flattened N/O distribution is derived.This result could be due to these H ii regions having low metallicities [12 + log(O/H) < ∼ 8.0], and for this regime the nitrogen has mainly a primary origin (Kennicutt et al. 2003).For the other three galaxies, we are sampling the inner part of the disk and they show similar behavior to M 101 for  < ∼ 20 kpc.For consistency, in M 101, we only consider N/O estimates for the 13 regions with  < ∼ 20 kpc, i.e.only those estimates showing a secondary nitrogen production.Mollá & Díaz (2005) presented multiphase chemical evolution models applied to a wide set of theoretical galaxies with different masses (9 < ∼ log(/M ⊙ ) < ∼ 13), evolutionary rates (age ranging from ∼ 0 to ∼ 13 Gyr) and distinct efficiencies of stellar formation.In these models, the enriched material proceeds from the restitution by dying stars, considering their nucleosynthesis predicted by the stellar products of massive stars by Woosley & Weaver (1995), low and intermediate-mass stars by Gavilán et al. (2005), and an initial mass function (IMF) by Ferrini et al. (1990) that is similar to the one proposed by Scalo (1986).These chemical evolution models predict elemental abundances (e.g.O, N, Fe, etc) along the hypothetical galaxy disks assuming that this is formed by radial regions, being the smallest  value of 1-4 kpc, depending on the galaxy mass.

Chemical evolution models
To fit the CEM results to the direct radial abundance for our sample of spiral galaxies the following methodology was adopted.
• Galaxy Mass: The maximum galactic radius ( gal ) is correlated to the galaxy mass ( gal ) by the expression (Lequeux 1983): being  max the maximum rotation velocity.Initially, we selected a set of models (1) whose predicted N/O abundances cover the range [Δ(N/O)] of direct estimates along the radius () of each galaxy of the sample, not taking into account a fitting to the abundance ratio gradient (see also Magrini et al. 2016).Basically, we consider being Δ(N/O) CEM and Δ(N/O) dir. the range of abundance ratios predicted by the CEMs and derived by the  e -method, respectively.Thus, this fitting procedure allows us to find the best set of CEMs  1).In each panel, blue crosses represent predictions from CEM by Mollá & Díaz (2005) for distinct galaxy mass ( gal ), as indicated, for all star formation efficiencies (not discriminated) and assuming an age of 13.2 Gyr.
(1) with distinct masses that describe the direct radial abundances.
For each galaxy of our sample, 1 contains one (unique solution) or at maximum two CEMs with distinct mass.For instance, in Fig. 2, the direct abundance estimates compiled from the literature (see Table 1) for M 101 galaxy is compared with CEM predictions by Mollá & Díaz (2005) assuming four  gal values (which cover the galaxy mass interval for the CEMs) and taking into account all star formation efficiencies (not discriminated).This grand design spiral galaxy has served as the prototype system for studies on gas-phase abundances in extragalactic objects (Kennicutt & Garnett 1996).We can see in Fig. 2 that the CEM with log( gal /M ⊙ ) = 13.07 describes clearly the direct estimates.
• Collapse time-scales (): the  value for each galaxy depends on its total mass through the expression given by Gallagher et al. (1984): being  the age, i.e. galaxies with lower masses have higher  than more massive ones.The age of a given galaxy is difficult to estimate because it requires the chronochemodynamics of stars to reconstruct a timeline of galaxy events (for a review see Soderblom 2010).As there are no age estimations for the galaxies in our sample, we assume only chemical models with 13.2 Gyr.Recent James Webb Space Telescope (JWST) observations have revealed the presence of galactic disks at 3 < ∼  < ∼ 12, corresponding to ∼ 2 and ∼ 0.3 Gyr, respectively, after the Big Bang (e.g.Ferreira et al. 2022a,b;Naidu et al. 2022;Robertson et al. 2023).Thus, it seems to be a good approach to assume for our sample ( ∼ 0.0) ages of 13.2 Gyr.
• Star formation efficiencies (): the CEMs by Mollá & Díaz (2005) assume star formation that takes place in two steps: first, the formation of molecular clouds with efficiency   ; then the cloud-cloud collision with efficiency   .Both   and   are re-lated to the star formation rate (SFR, Mollá & Díaz 2005).The star formation efficiency is defined by: with  constant along the galactic disk.Mollá & Díaz (2005) computed 10 models for each theoretical galaxy, allowing   ranging from 0 to 1 and   fixed according to the ratio ( ln   ln   ) ∼ 0.4.We follow the designation presented in Table 2 of Mollá & Díaz (2005), where CEMs with different   values are defined by the ¨N¨term.Mollá & Díaz (2005) built 440 different models that represent all possible combinations of the parameters described above.
Taking into account the set 1, we establish a new set of CEMs (2), now, comparing the predictions and direct N/O gradient parameters for each galaxy of the sample.To fit the CEM results to the direct radial N/O estimates, for each galaxy of our sample, we assume that a CEM is representative if it yields the minimum value of : being , where the ¨ob.¨and ¨mod.¨indicesindicate the observational and chemical model linear regressions parameters, respectively.The selected CEM results are plotted in Fig. 1 together with their corresponding linear regressions.A unique CEM is found to be representative of each galaxy.
In Fig. 3, direct abundance estimations [log(N/O)] as a function of galactocentric distances (R) are compared with CEM predictions by Mollá & Díaz (2005) for the M 101 galaxy.The CEM predictions are for distinct  values and fixed  gal and  values, as indicated.From this comparison, we can see that the theoretical radial N/O distributions produce a bad fitting to the direct estimates for large  values (see also Magrini et al. 2016;Mollá et al. 2019).The behaviour shown in Fig. 3 for M 101 is also derived for CEM predictions assuming other  gal values (not shown).The discrepancy derived for large  values, probably, is due to the simplistic approach that outer regions evolve at the same rate than inner galactic disk regions (e.g.Magrini et al. 2016), the effect of galactic interactions (more pronounced in outskirt regions, e.g.Toomre & Toomre 1972) and/or wrong SFR suppositions along the hypothetical disk (see Belfiore et al. 2017 and references therein).Recently, Garcia et al. (2023), by using gas-phase metallicity profiles of galaxies predicted by Illustris and IllustrisTNG simulations, showed that breaks in metallicity radial profiles are derived for galaxies at  =0-3 and with a wide mass range, being these due to radial gas mixing.The existence of breaks in the radial N/O linear distributions pointed out the need to assume only linear regressions in the inner parts of the galactic disks in order to infer nuclear abundances, as carried out in the present work.
The comparisons shown in Figs. 2 and 3 clearly indicate the need to impose observational constraints on CEMs.Otherwise, parameters (e.g.masses, star formation efficiencies) no representative to those of real galaxies can be assumed in the CEMs and, consequently, wrong interpretations from these could be obtained.
We need to infer the central intercept N/O abundance [(N/O) 0 ] from the extrapolation of the linear regressions fitted to the theoretical values predicted by the CEMs and along the cosmic time.Due to the possible wrong theoretical predictions for large  values, gradient predictions for the entire disks would produce incorrect (N/O) 0 values.To explore this fact, we study the behaviour of the logarithm of the (N/O) 0 obtained by fitting a linear regression (Eq. 1) to the  1) while solid and dashed black lines represent their fitting and its uncertainties, respectively.Color lines connect CEM results (colored points) for distinct star formation efficiency pairs   ,   according to their N designation in Table 2 by Mollá & Díaz (2005) as indicated.The galactic mass [log(  gal /M ⊙ )] and collapse time-scale () are 13.07 and 13.2 Gyr, respectively.
CEM results (considering the models with N=8 as was selected for M 101 from Fig. 3) in the 0.0-11.2redshift range (see Fig. 4).We fitted these linear regressions considering different galactocentric distances:  < 20 kpc and the entire disk.We also considered, in Fig. 4, the (N/O) values for  = 4 kpc as approximations for the (N/O) 0 values since this is the smallest galactocentric distance with abundance predictions by these CEMs and they are expected to have abundances near to those in nuclear regions.It can be noted in this figure that for  < ∼ 6 the log(N/O) 0 values estimated using the entire disks are ∼ 0.3 dex lower than those estimated using  < 20 kpc.Likewise, predictions assuming the log(N/O) at  = 4 kpc as the log(N/O) 0 values are, as expected, somewhat lower (∼ 0.15 dex) than those estimated via  < 20 kpc.For  > ∼ 6 all the suppositions produce similar (N/O) 0 values due to the flattened of gradients for high redshifts (see below).These results strengthen the assumption of only considering gradient estimates in disk parts where the nitrogen has mainly a secondary origin.Therefore, following this analysis performed for M 101, for the other galaxies in our sample (shown in Fig. 1), we only consider CEM predictions for the inner parts of their galactic disks, i.e. galactocentric distances similar to those for which direct estimates are available.In Table 1 the range of galactocentric distances considered in the CEM fittings are listed.

Photoionization models
We considered the version 22.00 of the Cloudy code (Ferland et al. 2013(Ferland et al. , 2017) ) in order to build up photoionization model grids representing narrow line regions (NLRs) of AGNs in a wide redshift range ( = 0.0 − 11.2).In this regard, the (N/O) 0 values obtained from the linear fitting to the CEM radial gradient predictions, for dis-  tinct redshifts and for each galaxy of our sample (see Sect. 2.1), and the correspondent O/H values (see below), were assumed as input parameters of the photoionization models.These models were built following the methodology presented by Dors et al. (2019) and the reader is referred to this paper for a complete description.In what follows the input parameters are briefly described.
• Spectral Energy Distribution (SED): the SED was assumed to be the one defined by the Table AGN command in the Cloudy code and it is similar to that deduced by Mathews & Ferland (1987).We compared emission line intensities ratios, considered in BPT diagrams, predicted by our photoionization models with those assuming SEDs derived from observations of the AGNs Mrk 509 and NGC 5548 by Kaastra et al. (2011) and Mehdipour et al. (2015), respectively4 .We found that photoionization models assuming these distinct SEDs produce line intensities ratios that differ by a factor lower than ∼ 0.1 dex (see also Feltre et al. 2016), similar to the uncertainties in observational line ratio intensities (e.g.Hägele et al. 2008).
• Electron density ( e ): in all photoionization models  e was assumed to be constant along the nebular radius and equal to 500 cm −3 , about the mean value found by Dors et al. (2014) for a large number of Seyfert 2 nuclei and derived through the [S ii]6716/6731 line ratio (see also Vaona et al. 2012;Zhang et al. 2013;Flury & Moran 2020;Agostino et al. 2021;XueGuang 2023).Armah et al. (2023) investigated the effect of the presence of a radial density profile (e.g.Congiu et al. 2017;Kakkad et al. 2018;Revalski et al. 2018;Freitas et al. 2018;Mingozzi et al. 2019;Revalski et al. 2021) and high density values ( e > 10 4 cm −3 , see Vaona et al. 2012;Congiu et al. 2017;Cerqueira-Campos et al. 2021) in NRLs on photoionization model predictions.These authors found that density variations and high  e values yield changes in the emission-line intensities ratios used in BPT diagrams lower than 0.1 dex.
• Ionization parameter (): we assume the logarithm of  ranging from −3.5 to −1.0 dex, with a step of 0.5 dex.A similar range of  was derived by Carvalho et al. (2020) from a comparison between photoionization model results and spectroscopic data of ∼ 430 Seyfert 2 (see also Pérez-Montero et al. 2019;Pérez-Díaz et al. 2022).
• Abundance set: the abundance of all metals () was linearly scaled with the solar abundance5 by a factor  , i.e.  =  ×  ⊙ , with the exception of the nitrogen abundance.For each photoionization model representing a galaxy (see Sect. 2.1) at a given ,  is obtained from the (N/O) 0 value derived through the CEM using the N/O radial gradient.
Recently, Johnson et al. (2023), who proposed a multizone galactic chemical evolution model for Milky Way-like galaxies, showed that the (N/O)-(O/H) relation is dependent on the galactocentric distance, being this dependence mainly due to the pAGB stars evolution that enrich the ISM with N but with negligible amounts of O, increasing N/O.Basically, as spiral galaxies begin to form their inner regions before the outer ones in a classical inside-out scheme (e.g.Samland et al. 1997;Portinari & Chiosi 1999;Boissier & Prantzos 2000), the ISM at distinct  values are enriched by stars with distinct evolutionary stages, resulting in a dependence of the (N/O)-(O/H) relation with .In order to study this behavior using the assumed CEMs, we  A1-A4, and used as inputs in our photoionization models.To define the metallicity () through the oxygen abundance x 0 =12+log(O/H) 0 , for each model representing a galaxy at a given , the derived (N/O) 0 was used in Eq. 7 and, thus, its corresponding O/H abundance was obtained.We assume that Eq. 7 is still valid for  = 0 kpc, and not only for the innermost galactocentric distances for which the CEMs were built.A similar, but inverse, procedure was adopted by Garg et al. (2022), where the O/H abundances were retrieved from the simba simulations and N/O from observational estimates obtained by Pilyugin et al. (2012).
The helium abundance in relation to the hydrogen (He/H) was defined, in each model, by the expression w =(0.1215 ± 0.0422) × x 2 − (1.8183 ± 0.6977) × x + (17.6732 ± 2.8798), where w=12+log(He/H) and x=12+log(O/H).This expression was derived by Dors et al. (2022) through the  e -method and by using a large number of AGN and H ii region estimates (see also Dopita et al. 2006).
Recently, Zhu et al. ( 2023) compared dust-free and dusty AGN photoionization model results with observational AGN line intensities ratios in standard BPT diagrams (see also Feltre et al. 2016;Peluso et al. 2023).These authors found some significant discrepancy among the dust-free and dusty models only for the high metallicity regime, i.e. 12 + log(O/H) ≳ 8.7 [(/Z ⊙ ) ≳ 1.0], with discrepancies for the line intensities ratios reaching up to 0.5 dex for 12 + log(O/H) = 9.2 [(/Z ⊙ ) ∼ 3.2].As the abundance of dust in the gas phase of gaseous nebulae is poorly known due to the difficulty of estimating this parameter (e.g.Sofia et al. 1994;Garnett et al. 1995;Peimbert & Peimbert 2010;Brinchmann et al. 2013;Martín-Doménech et al. 2016) and, since our models do not extend to the very high  regime, we followed Nagao et al. (2003Nagao et al. ( , 2006b)), who found a better match between dust-free AGN models to the observed emission line ratios in comparison to dusty models, and we assumed that all AGN photoionization models are dust-free.

Galaxy mass and SF efficiency
In Fig. 1 the direct ( e -method) N/O abundance estimations as a function of the galactocentric distances for the objects in our sample together with the results of the CEMs that better fit them are shown.In all the cases, it can be seen a good agreement between the linear fittings to the observational data and the CEM results, except for M 33 where the CEM fitting is offset by ∼ 0.05 dex from the one directly derived even though both fittings show the same slope.This difference is in the order of the uncertainties in abundance estimates via the  e -method (∼ 0.1 dex, Kennicutt et al. 2003;Hägele et al. 2008) and lower than those via strong-line methods (∼ 0.2 dex, e.g.Denicoló et al. 2002).Therefore, it has a small influence on the results obtained in the present work.
In Table 1 the CEM fitting parameters are listed.Masses of the theoretical galaxies are higher than those derived from observations by a factor up to ∼ 2.4 dex.Discrepancies between observational and theoretical  gal estimations shown in Table 1 could be mainly due to only the inner parts of the galactic disks (the ones with direct N/O estimates) are considered in the CEM fittings, being any theoretical mass derivation through the CEMs highly uncertain.
We are interested in the nuclear N/O abundance, i.e. (N/O) 0 , in order to provide it (and the corresponding O/H) as an initial parameter for our photoionization models and, thus, to investigate which abundances (and redshift) make AGNs trend to go through the maximumstarburst line in BPT diagrams.The (N/O) 0 values are obtained from radial gradients and they are dependent on the galaxy mass.Mollá & Díaz (2005) showed that O/H radial gradients are steeper for the models with lower maximum rotation velocities ( max ) or lower galaxy masses ( gal , see also Belfiore et al. 2019).Consequently, for disk regions where nitrogen has mainly a secondary production, N/O gradients are also stepper when  gal decreases.Therefore, it is important that the range of (N/O) 0 derived for our sample be representative of spiral galaxies.In order to verify that, in Fig. 7 the distribution of (N/O) 0 values derived through the -method (a method that yields abundances in agreement with those estimated via the  e -method; Pilyugin 2001) for 45 nearby spiral galaxies by Pilyugin et al. (2004) are compared with the (N/O) 0 range predicted by CEMs for our sample of nearby galaxies ( ∼ 0).The log(N/O) 0 for our four spiral galaxies ranges from ∼ −1.0 to ∼ −0.5 dex, being ∼ 70 per cent of the spiral galaxies studied by Pilyugin et al. (2004) have log(N/O) 0 within this range.The extreme (N/O) 0 values estimated by Pilyugin et al. (2004) are not encompassed by our results.More N/O radial estimates via the  e -method for spiral galaxies with a wide mass range could be necessary to extend our results.
Concerning the star formation efficiency , we found that CEMs with  equal to 7 and 8 (see Table 2 by Mollá & Díaz 2005) reproduce the radial N/O gradients of our sample.Magrini et al. (2016) also carried out fittings of the CEMs by Mollá & Díaz (2005) to O/H radial gradients of the spiral galaxies M 31, M 33 and NGC 300.These authors found that CEMs assuming the highest  value (i.e. = 10) are satisfactory to represent the gradients of their galaxies.These three  values (i.e. = 7, 8 and 10) represent star formation with low efficiency in comparison to the  representing the Milky Way (i.e. = 4), according to Mollá & Díaz (2005).The source of this discrepancy is, probably, due to Mollá & Díaz (2005) considering a larger number of constraints in their models (e.g.neutral gas distribution, star formation) in comparison to us and Magrini et al. (2016).

Predicted radial gradients
In Fig. 8, bottom panel, the slopes of the N/O radial gradients [grad(N/O)] versus the redshift () predicted by the CEMs by Mollá & Díaz (2005) are shown.Points for  = 0 are obtained from the fittings shown in Fig. 1.Points for  > 0 are obtained by carrying out the fitting of Eq. 1 to the CEM predictions (not shown) following the same procedure than for  = 0 (see Sect. 2.2).The error bars indicate the uncertainties in the linear fittings to the CEM N/O gradients for a given  value.We can see that grad(N/O) flattens out as the  increases showing negative values from  = 0 to  ∼ 6 ( ∼ 0.9 Gyr), depending on the galaxy mass.For  > ∼ 6, all the gradients become almost flat, until an inversion (positive slopes) for  > ∼ 8.
Unfortunately, abundance gradients of galaxies at high redshift have been rarely derived and only for O/H, mainly due to instrumental constraints to observe the red optical spectral range and/or separate [N ii] lines from H, making almost unknown the grad(N/O)- relation at  > 1.In fact, Simons et al. (2021) derived the O/H (or ) gradients, relying on strong-line methods, for a sample of 238 star-forming galaxies at 0.6 <  < 2.6.These authors found that a large fraction of galaxies have flat and positive O/H gradients (see also Sharda et al. 2021;Wang et al. 2022).Moreover, even in the local universe ( < 1), flat and positive O/H gradients are derived in some objects (e.g.Krabbe et al. 2008Krabbe et al. , 2011;;Carton et al. 2018;do Nascimento et al. 2022;Boardman et al. 2023), being this result due to, for instance, interactions of galaxies (e.g.Rosa et al. 2014), metal-poor gas accretion (e.g.Kereš et al. 2005) and metal-enriched of the circumgalactic medium by stellar feedback and later gas reaccretion (e.g.Tumlinson et al. 2017).In any case, N/H gradients trend to be steeper than O/H gradients (e.g.Pilyugin et al. 2004), resulting in negative N/O gradients in practically all spiral galaxies in the local universe.Belfiore et al. (2017), by using strong-line methods, derived the O/H and N/O radial gradients for a sample of 550 nearby ( < 0.15) galaxies.These authors found that N/O gradients do not flatten in the innermost regions of galaxies, where a flattening of the oxygen abundance gradient is observed in some objects.Thus, probably, future direct estimates of N/O gradients in high- galaxies could confirm our findings, i.e.N/O gradients in inner disk of spiral galaxies only flatten at very high redshift values ( > ∼ 5).In Fig. 8, top panel, the logarithm of (N/O) 0 versus the redshift predicted by the CEM fittings are shown.The hatched area represents the log(N/O) value for which the nitrogen has mainly a primary origin, i.e. the N/O plateau [log(N/O) = −1.41± 0.09] as derived by Berg et al. (2019) through the  e -method and by using observational optical spectroscopic data of local dwarf galaxies.It can be seen that for  > ∼ 5, the N/O nuclear abundance variations with the redshift become smaller trending the (N/O) 0 values to be in the shadowed area.Therefore, the nitrogen ISM enrichment in nuclear regions of spiral galaxies at very high redshift is predominately of primary origin.This result has a deep consequence on  estimates of AGNs through strongemission line ratios involving nitrogen lines (e.g.[N ii]6584/H, Carvalho et al. 2020), in the sense that, in this case, there is a weak (or nonexistent) dependence of such ratios with .Thus, for  estimates of AGNs with  > ∼ 5 it is advisable to use calibrations based on, for instance, the  23 =([O ii]3727+[O iii]4959 + 5007)/H index, as suggested by Dors (2021).
Up to now, due to the lack of N/O estimates for AGNs at high- in the literature, it is only possible to compare our N/O nuclear abundance predictions with those derived from integrated spectra of SFs.In this regard, in Fig. 9 our (N/O) 0 predictions are compared with abundances for high redshift SFs derived through the  e -method and compiled from the literature.Moreover, it is worth mentioning that we are comparing abundance estimates in distinct object classes, which have different star formation histories (see Riffel et al. 2023).Taking into account the high observational uncertainties in the N/O direct estimations (∼ 0.2 dex) of the SFs and the different nature of AGNs and SFs, we are able to assume that most part (9/12) of the N/O direct estimates are reproduced by our CEM predictions.
The (N/O) 0 - relation obtained from our nuclear abundance predictions is represented by: where a = 0.007 ± 0.001, b = −0.13± 0.01 and c = −0.83± 0.02.This relation is represented in Fig. 9 by a black curve.In Fig. 10, bottom panel, our predicted CEM O/H values for the  and it is plotted as a black line.With the goal of comparing our results, we overlapped on this figure other estimates taken from the literature.
In the bottom panel of Fig. 9, our O/H abundance predictions as a function of  are compared with: • Estimates of O/H abundances obtained by using the  e -method for 48 SFs at the redshift range 1.4 < ∼  < ∼ 9.5.These estimates are represented by pink points while the pink line is the liner fitting to them given by 12 + log(O/H) SF = −0.03± 0.02 ×  + 7.89 ± 0.11. (12) • O/H abundances for SF galaxies at the redshift range 0.07 < ∼  < ∼ 3.5, with stellar masses ranging from 9.0 < ∼ (log( * /M ⊙ ) < ∼ 11.0 derived through strong-line methods by Maiolino et al. (2008) are included in the figure as a region delimited by orange lines.
• Estimates (blue points) for narrow line regions of AGNs (1.2 < ∼  < ∼ 3.8) based on strong emission lines observed in the ultraviolet wavelength range and obtained from comparison with photoionization model predictions by Dors et al. (2019).
In Fig. 9, top panel, our O/H abundance predictions as a function of  are compared with: • Estimates (pink points) obtained through absorption lines for 151 Damped L and sub-Damped L galaxies (DLA) at a redshift range of 0.2 < ∼  < ∼ 5.1.A linear fit to these estimates yields: 12 + log(O/H) DLA = −(0.37 ± 0.03) ×  + (8.33 ± 0.09).( 13)  et al. 2001).It is worth noting that, for  < ∼ 5, our predicted nuclear O/H abundances are in consonance with those for SFs and AGNs.Only a few active nuclei show O/H values higher, by ∼ 0.2 dex, than our CEM estimates.However, for  > ∼ 5, our predictions yield O/H values systematically lower than those via the  e -method for SFs, showing the latter a small O/H abundance decrease with .The origin of this discrepancy is likely due to an observational bias, in the sense that the galaxies with higher luminosity (and metallicity) are easier to observe (e.g.Nagao et al. 2006c).The DLA estimates show lower O/H values and a steeper decrease with  than our predictions.This behavior was also found by Dors et al. (2014), who compared  estimates of AGNs and DLAs, and it is probably due to the use of distinct methods (i.e. estimates via emission and absorption lines) rather than an observational bias.
In summary, from the comparisons presented in Figs. 9 and 10, we are able to conclude that our nuclear N/O and O/H abundances are in consonance with those for the most part of SFs located at  < ∼ 5.For  > ∼ 5, we found a discrepancy that could be due to the existence of a bias in galaxy observations.In the case of DLAs the disagreement could be attributed to the distinct methods applied to estimate .
Finally, in Fig. 11, our CEM oxygen abundance predictions [(O/H) 0 ] are compared with the extrapolated O/H abundances to the nuclear regions ( = 0 kpc) of 14 high- galaxies observed by Simons et al. (2021).These galaxies have masses in the 8.6 < ∼ log( * /M ⊙ ) < ∼ 10 range, redshifts in the 1 < ∼  < ∼ 2 range and O/H radial gradients derived through the izi Bayesian photoion-  A1-A4.The blue curve represents the fitting to our estimates (Eq.10).ization fitting code (Blanc et al. 2015).It can be seen in Fig. 11 that our CEM estimates for nuclear regions of galaxies are in most part in agreement with those from Simons et al. (2021).Again, more estimations of nuclear metallicities or radial abundance gradients for high- ( > 2) galaxies are necessary to produce a better validation of our estimates.
It is worth mentioning that, in general, supersolar metallicities are derived for Broad Line Regions of AGNs (e.g.Dietrich et al. 2003;Nagao et al. 2006a;Juarez et al. 2009;Sameshima et al. 2017), reaching up ∼ 18  ⊙ in N-loud quasi-stellar objects (QSOs, e.g.Batra & Baldwin 2014).Recent metallicity estimates for high-redshift (5.8 <  < 7.5) quasars have also found high metallicities (∼ 2−4 ⊙ , Lai et al. 2022), indicating no evolution with the redshift (e.g.Dors et al. 2014;Onoue et al. 2020).As shown by Dors et al. (2019),  values in BLRs are a factor of about 2-3 higher than those in NLRs and are inconsistent with predictions of chemical evolution models (see Fig. 11).The origin of this discrepancy has been attributed to star formation in accretion disks (e.g.Fan & Wu 2023;Huang et al. 2023;Cheng & Loeb 2023;Wang et al. 2023;Dittmann et al. 2023), uncertainties in BLR metallicity estimates from some UV lines (e.g.Matsuoka et al. 2017) and the fact that broad lines originate in a small region with a radius lower than 1 pc (e.g.Kaspi et al. 2000;Suganuma et al. 2006), which may evolve more rapidly than the NLRs (e.g.Matsuoka et al. 2018).It is beyond of the scope of this study to investigate the BLR/NLR metallicity discrepancy.We only pointed out that, apparently, NLR estimates seem to be better suited as a proxy for the properties of the host galaxy since the spatial extent of the NLR region (10 1−4 pc) is larger than the BLR ( < 1 pc, Armah et al. 2023).2023).Blue points: AGN metallicity estimates via photoionization models by Dors et al. (2019).Orange lines delimit the region occupied by SF galaxies [9.0 < ∼ log(  * /M ⊙ ) < ∼ 11.0] whose O/H abundances were derived through strong-line methods by Maiolino et al. (2008).Top panel: Pink points represent metallicity estimations for Damped L and sub-Damped L galaxies via absorption lines by Rafelski et al. (2013), Fox et al. (2007) and Kulkarni et al. (2005).Black points are as in bottom panel.

Optical diagnostic diagrams
The main goal of this work is to analyze which are the physical parameters (i.e.O/H and ), and the correspondent  values, that let AGNs fall below the maximum starburst lines in diagnostic diagrams.In this regard, in Figs. 12 and The O/H gradients were calculated by these authors through the izi Bayesian photoionization fitting code (Blanc et al. 2015).
Sky Survey (SDSS) DR176 (Abdurro'uf et al. 2022) database are shown.The procedures applied to these data are the same as in Dors et al. (2020).
Firstly, we can see in Fig. 12 that AGN models with 12 + log(O/H) > ∼ 9.0 [(/Z ⊙ ) > ∼ 2.0)], independently of the assumed  value, are located near or just below the maximum starburst lines of the classical BPT diagrams (see also Feltre et al. 2016;Nakajima & Maiolino 2022).However, similar very high oxygen abundance (or metallicity) is obtained in some few AGNs in the local universe ( < 1).In fact, O/H estimates for a large sample (463 objects) of local ( < 0.4) AGNs obtained through strong-line methods and inferred extrapolating radial abundance gradients by Dors et al. (2020) showed that objects with 12 + log(O/H) > ∼ 9.0 comprehends only ∼ 10 per cent of their sample.In any case, strong-line methods based on [N ii]/[O ii] and [N ii]/H lines ratios proposed by Castro et al. (2017) and Carvalho et al. (2020), respectively, can be used to distinguish between metal-rich and metal-poor AGNs.
Secondly, in Fig. 12, it can be seen the known result (see e.g.Feltre et al. 2016) that low metallicity AGNs occupy a large area in the BPT diagrams where SFs are located.In each diagram, we indicated the average O/H abundance for which AGNs pass through the maximum starburst line, i.e.AGNs with 12 + log(O/H) < ∼ 8.0 or (/Z ⊙ ) < ∼ 0.2 are located in BPT regions occupied by SFs.In each diagram of Fig. 13 we indicated the average redshift, according to our CEM results, for which AGNs reach the oxygen abundance of 12+log(O/H)≃8.0.
• AGN galaxies with these abundances are preferably found at  > ∼ 1.
The abundance discrepancy between our results and those obtained by Hirschmann et al. (2022) is, probably, due to distinct photoionization parameters [e.g.SED, (N/O)-(O/H) relation, dust content, atomic parameters] assumed in both works.In any case, assuming the abundance uncertainty of ∼ 0.2 dex for strong-line methods (e.g.Denicoló et al. 2002), the O/H values derived by us and by Hirschmann et al. (2022) are in agreement to each other.
Regarding the redshift value discrepancy, it could be due to the different  −  relations derived in both works.Unfortunately, this relation is not presented by Hirschmann et al. (2022).However, CEM predictions by Hirschmann et al. (2017), relied on the SPHGal code (Hu et al. 2014) for nuclear regions of three hypothetical galaxies with log( gal / * ) ∼ 11, indicate that galaxies with redshift around  = 4 reach abundances that result in emission-line ratios in the SF-zone, in good agreement with our present estimation.
We compare our  −  relation with others derived by considering some CEMs available in the literature.For that, we consider: • Chemical Evolution models by Pei et al. (1999) calibrated by using optical imaging surveys of quasars and the COBE DIRBE and FIRAS extragalactic infrared background measurements.These models are valid for 0.0 ≤  ≤ 5.0.
• Semi-analytic models of galaxy formation set within the cold dark matter merging hierarchy built by Somerville et al. (2001).These models are able to reproduce the metallicity of high-redshift Lyman-break galaxies in the redshift range  = 0.0 − 4.0.
• Metallicity predictions from IllustrisTNG simulations by Torrey et al. (2019) for the redshift range  = 0.0−10.0 and for hypothetical galaxies with masses in the range 8.0 ≤ log(/ * ) ≤ 10.5.We consider the average O/H value predicted by this mass interval.Although Hirschmann et al. (2022) used the predictions of these models, these authors assumed as input in their photoionization models metallicity values obtained in a co-moving sphere of 1 kpc radius around the hypothetical nuclear galaxy region, about the same approach that the one assumed in the present work.
• Recent multi-zone galactic chemical evolution models proposed by Johnson et al. (2023) which reproduces the (N/O)-(O/H) relation found for the Milky Way stars and for the gas phase of external galaxies.We selected CEMs by Johnson et al. (2023) considering radial stellar migration proceeds and the O/H abundance predicted by the innermost galactocentric distance , i.e  = 4 kpc.This procedure was adopted in order to obtain the nearest nuclear abundances for galaxies with distinct evolutionary stages.These models are appropriate for Milky Way mass galaxies with 10.5 ≤ log( gal / * ) ≤ 11.
Except for Johnson et al. (2023), the CEMs listed above predict an average metallicity for the ISM for the whole galaxy at a given redshift, i.e. they do not predict nuclear metallicity estimates in the hypothetical galaxies.In Fig. 14, the O/H abundance as a function of the redshift predicted by the CEMs listed above is shown.Also in this figure, the redshift at which AGNs reach the abundance 12+log(O/H)=8.0 [(/Z ⊙ ) = 0.2] are indicated for the distinct CEMs (green points).Hereafter, we define this redshift value, i.e. the redshift in which AGNs reach 12+log(O/H)=8.0, as the  AGN -limit.We can see that the  AGN -limit is strongly dependent on the CEMs considered, ranging from  ∼ 1.4 to  ∼ 10 for all the models except the model by Johnson et al. (2023).Since the nuclear galaxy abundance trends to be higher than the one of the whole galaxy (e.g.Hirschmann et al. 2017) or than that of the ISM, the  AGN -limit indicated in Fig. 14 must be interpreted as a lower limit.In this scenario, our result that AGNs with 12 + log(O/H) < ∼ 8.0 are preferable found at  > ∼ 3.7 is in consonance with the estimated range for three of the CEMs from the literature.The CEM results from Johnson et al. (2023) predict higher O/H abundances along the redshift interval considered.Possibly, it is due to the models proposed by Johnson et al. (2023) being designed for massive galaxies (see above).
Finally, we discuss our results for the He ii/H versus [N ii]/H diagram shown in the lower right panels of Figs. 12 and 13.As was mentioned before, this diagram was suggested by Shirazi & Brinchmann (2012) to separate AGNs from SFs.These authors defined the maximum starburst line by using observational data of galaxies ( < ∼ 0.2) from the Sloan Digital Sky Survey (SDSS) DR7 (York et al. 2000;Stoughton et al. 2002) previously classified through classical BPT diagrams.Since the nebular He ii emission has a high ionization potential (54.4 eV), it is preferably produced by a hard ionizing spectrum, which may indicate AGN activity.Recently, Tozzi et al. (2023), who used observational data from Mapping Nearby Galaxies at APO survey (MaNGA DR15, Bundy et al. 2015), found an overall increased fraction (2 per cent) of AGNs in MaNGA galaxies when the He ii diagram is used (11 per cent) in comparison to the BPT-only census (9 per cent).It is worth to mention that the He ii/H versus [N ii]/H diagram suffers a limitation due to the difficulty to measure the He ii 4685 emission line (about 4-10 times weaker than H, see Fig. 12).In fact, Shirazi & Brinchmann (2012), based on SDSS spectroscopic data ( < ∼ 0.4), showed that the He ii 4685 is most frequently measured in SF galaxies at low redshift ( < ∼ 0.1, see also Sartori et al. 2015;Bär et al. 2017;Koss et al. 2017;Wang & Kron 2020;Mayya et al. 2020;Oh et al. 2022;Tozzi et al. 2023).However, recent observations with the JWST (e.g.Übler et al. 2023) have shown the reliability of measuring this line in AGNs at very high redshift ( > ∼ 5).Thus, the new generation of extreme large and space telescopes will increase the feasibility of diagnostic diagrams using He lines, such as the He ii/H versus [N ii]/H (Shirazi & Brinchmann 2012) and [O iii]/[O ii] versus He ii/He i (Dors et al. 2022).
We can see in Figs. 12 and 13 that our cosmological predictions indicate that, independently of the metallicity and of the redshift, AGNs are always above the maximum starburst line in the He ii/H versus [N ii]/H diagram, indicating that this diagram is the most appropriated to classify objects at any redshift, independently of the AGN nebular parameters (i.e.O/H and ).It is worth to mention that, as can be seen in Figs. 12 and 13, our theoretical results indicate a similar AGN-SF separation line in agreement with that proposed by Molina et al. (2021), who presented a sample of 81 dwarf galaxies ( * ≤ 3 × 10 9  ⊙ ) with detectable [Fe x]6374 coronal emission line, indicating accretion onto massive black holes.The criterion proposed by Molina et al. (2021) suggests that objects with log(He ii4686/H) ≳ −1.0, are classified as AGNs, otherwise, SFs.In Figs. 12 and 13 one can see that the photoionization model results do not reproduce some SDSS observational data.This is due to the limited sample of galaxies considered by us (which define the space of model parameters) in comparison to the data from SDSS sample.It is worth mentioning that, for a better match, is necessary a larger sample of spiral galaxies with direct radial N/O gradients.

CONCLUSIONS
We investigated the reliability of optical diagnostic diagrams, based on emission-line ratios [4000 < (Å) < 7000], in the framework of classifying Active Galactic Nuclei (AGNs) according to the metallicity evolution in the redshift range 0 ≤  ≤ 11.

Figure 2 .
Figure2.Logarithm of the N/O abundance ratio versus the galactocentric distance  (in kpc) for the M 101 galaxy.Black points represent direct abundance estimates compiled from the literature (see Table1).In each panel, blue crosses represent predictions from CEM byMollá & Díaz (2005) for distinct galaxy mass ( gal ), as indicated, for all star formation efficiencies (not discriminated) and assuming an age of 13.2 Gyr.

Figure 3 .
Figure3.Logarithm of the N/O abundance ratio versus the galactocentric distance  (in kpc) for the M 101 galaxy.Black points represent direct abundance estimates compiled from the literature (see Table1) while solid and dashed black lines represent their fitting and its uncertainties, respectively.Color lines connect CEM results (colored points) for distinct star formation efficiency pairs   ,   according to their N designation in Table2byMollá & Díaz (2005) as indicated.The galactic mass [log(  gal /M ⊙ )] and collapse time-scale () are 13.07 and 13.2 Gyr, respectively.

Figure 4 .
Figure 4. Logarithm of (N/O) 0 (extrapolated abundance to the galactic center) at different evolution times considering predictions from the CEMs by Mollá & Díaz (2005) for M 101 and assuming  =8 as selected through Fig. 3. Points represent (N/O) 0 assuming radial gradients estimated by fitting the Eq. 1 to the CEM predictions for the galactocentric distance  < 20 kpc, entire disk and assuming (N/O) 0 equal to (N/O) for  = 4 kpc (see text), as indicated.Error bars represent the uncertainty in (N/O) 0 .

Figure 7 .
Figure 7. Distribution (in black) of the central intercept N/O abundance [log(N/O) 0 ] values derived by Pilyugin et al. (2004) using the -method(Pilyugin 2001) for 45 local ( < 0.02) galaxies.Red hatched area represents the log(N/O) 0 range derived from our CEM fittings to direct radial gradients of four local spirals (see Table1).

Figure 8 .
Figure 8.Chemical evolution model predictions of the radial gradient slope (in units of kpc dex −1 ) for N/O (bottom panel) and for the logarithm of (N/O) 0 (top panel).Distinct colors indicate the predictions of CEMs by Mollá & Díaz (2005) for the spiral galaxies as indicated.The dashed line in bottom panel represents the zero value of the N/O gradient.Hatched area in top panel represents the N/O plateau [log(N/O) = −1.41± 0.09] derived by Berg et al. (2019).
13, the [O iii]/H versus [O i]/H, [N ii]/H and [S ii]/H BPT diagrams containing our photoionization model results are shown.Also in these figures, the He ii/H versus [N ii]/H diagram, proposed by Shirazi & Brinchmann (2012) to separate AGNs from SFs, are shown.In Fig. 12, the colour bars represent the O/H abundance of the photoionization models, while in Fig. 13 the corresponding redshifts for which the theoretical AGNs reach the O/H abundance.The model results assuming different log  values are represented by distinct symbols, as indicated.In each diagram the maximum starburst lines, taken from the literature, are also plotted.Likewise, lines to separate AGNs from low-ionization nuclear emission-line regions (LINERs, represented by, simplicity, by LIN) are plotted in the [O iii]/H versus [N ii]/H and [S ii]/H diagrams.Also in Fig. 12 and Fig. 13, observational emission line ratio intensities, corrected by reddening and after apply stellar population continuum subtraction, of galaxies taken from Sloan Digital

Figure 12 .
Figure 12.Diagnostic diagrams involving optical emission-line ratios as indicated.In the [O iii]/H versus [O i]/H, [N ii]/H and [S ii]/H the black curves represent the maximum starburst line proposed by Kewley et al. (2001).The dashed line separates AGNs from LINERs (refereed as LIN) as suggested by Kewley et al. (2006).The [S ii]6725 corresponds to the sum [S ii]6716 + 6731.In the He ii/H versus [N ii]/H diagram, the line represents the maximum starburst line proposed by Molina et al. (2021).In all diagrams, points represent results of our photoionization models (see Sect.2.3) by using as input parameters the abundances predicted by CEMs (see Sect. 2.2) of Mollá & Díaz (2005).Symbols represent photoionization model with iso-, as indicated.Color bars indicate the oxygen abundance assumed in the photoionization models.AGNs with oxygen abundance values lower than the ones indicated in each of the three BPT diagrams [average 12+log(O/H) values] are located below the maximum starburst lines.Contours are galaxies whose data were taken from SDSS DR-17 (Abdurro'uf et al. 2022).
Figure 12.Diagnostic diagrams involving optical emission-line ratios as indicated.In the [O iii]/H versus [O i]/H, [N ii]/H and [S ii]/H the black curves represent the maximum starburst line proposed by Kewley et al. (2001).The dashed line separates AGNs from LINERs (refereed as LIN) as suggested by Kewley et al. (2006).The [S ii]6725 corresponds to the sum [S ii]6716 + 6731.In the He ii/H versus [N ii]/H diagram, the line represents the maximum starburst line proposed by Molina et al. (2021).In all diagrams, points represent results of our photoionization models (see Sect.2.3) by using as input parameters the abundances predicted by CEMs (see Sect. 2.2) of Mollá & Díaz (2005).Symbols represent photoionization model with iso-, as indicated.Color bars indicate the oxygen abundance assumed in the photoionization models.AGNs with oxygen abundance values lower than the ones indicated in each of the three BPT diagrams [average 12+log(O/H) values] are located below the maximum starburst lines.Contours are galaxies whose data were taken from SDSS DR-17 (Abdurro'uf et al. 2022).

Figure 13 .Figure 14 .
Figure 13.As the Fig. 10 but with the color bars indicating the redshift value of each photoionization model representing an AGN.AGNs with redshift values ( AGN -limit) lower than the ones indicated in each BPT diagram (average z values) are located below the maximum starburst lines.

Table 1
).nuclear galaxy regions are plotted against the redshift.The expected tendency of the decrement of O/H (or ) with the redshift is easily noticeable.This relation can be represented by a linear regression given by: 12 + log(O/H) CEM = −(0.24±0.01 × ) + 8.64 ± 0.06(11) 2. With this aim, we fit results of Chemical Evolution Models (CEMs) to the radial abundance gradients derived through direct estimates of electron temperatures ( e -method) in a sample of four local spiral galaxies.Unlike the majority of previous studies, we consider as metallicity tracer the N/O abundance ratio to select the representative CEMs for each object belonging to our galaxy sample.The N/O abundance radial predictions provided by the CEMs were extrapolated to the nuclear regions (galactocentric distance  = 0 kpc) in order to infer N/O and O/H abundances in theoretical galaxies at distinct evolutionary stages.These resulting abundance ratios were used as input parameters in AGN photoionization models built with the Cloudy code.Extensive comparison between CEM predictions and direct abundance estimates, in a wide range of redshift, was performed in order to validate our cosmic abundance values.From our analysis we conclude that: Our CEM predictions show that AGNs reach the oxygen abundance of 12 + log(O/H) = 8.0 preferably at redshift  ∼ 4, indicating that BPT diagrams break down for higher redshift values.•We found that the He ii4685/H versus [N ii]6584/H diagram is able to separate AGNs from star-forming regions in the redshift range 0 ≲  ≲ 11.2.

Table A1 .
Mollá & Díaz (2005)ntersect values of the N/O regressions estimated using the predictions of the CEMs built byMollá & Díaz (2005)for different redshift.The term grad represents the slope of the gradients in units of dex kpc −1 .The term log(N/O) 0 represents the extrapolated abundance ratio to the galactocentric distance () equal to zero.The central oxygen abundance 12+log(O/H) 0 is estimated using log(N/O) 0 and Eq.7.