11institutetext:  Dipartimento di Fisica e Astronomia, Università di Bologna,
Via Gobetti 93/2, I-40129 Bologna, Italy 22institutetext: INAF – Osservatorio
di Astrofisica e Scienza dello Spazio di Bologna, Via Gobetti 93/3, I-40129
Bologna, Italy 33institutetext: Kavli Institute for Cosmology, University of
Cambridge, Madingley Road, Cambridge, CB3 0HA, UK 44institutetext: Cavendish
Laboratory, University of Cambridge, 19 JJ Thomson Avenue, Cambridge, CB3 0HE,
UK 55institutetext: Department of Physics and Astronomy, University College
London, Gower Street, London WC1E 6BT, UK 66institutetext: National
Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan
77institutetext: INAF-Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5,
I-50125, Firenze, Italy 88institutetext: European Southern Observatory, Karl-
Schwarzschild-Strasse 2, 85748 Garching, Germany 99institutetext: Dipartimento
di Fisica e Astronomia, Università di Firenze, Via G. Sansone 1, I-50019,
Sesto F.no (Firenze), Italy

# New AGN diagnostic diagrams based on the [OIII]$\lambda 4363$ auroral line

Giovanni Mazzolari<EMAIL_ADDRESS>112233 Hannah Übler  Roberto
Maiolino 3344 Xihan Ji 334455 Kimihiko Nakajima 3344 Anna Feltre 66 Jan
Scholtz 77 Francesco D’Eugenio 3344 Mirko Curti 3344 Marco Mignoli 88
Alessandro Marconi 227799

The James Webb Space Telescope (JWST) is revolutionizing our understanding of
black hole formation and growth in the early Universe. However, JWST has also
revealed that some of the classical diagnostics, such as the BPT diagrams and
X-ray emission, often fail to identify narrow line TypeII active galactic
nuclei (AGN) at high redshift. Here we present three new rest-frame optical
diagnostic diagrams leveraging the [OIII]$\lambda 4363\,$auroral line, which
has been detected in several JWST spectra. Specifically, we show that high
values of the [OIII]$\lambda 4363\,$$\bigl{/}$H$\gamma\ $ratio provide a
sufficient (but not necessary) condition to identify the presence of an AGN,
both based on empirical calibrations (using local and high-redshift sources)
and a broad range of photoionization models. These diagnostics are able to
separate much of the AGN population from Star Forming Galaxies (SFGs). This is
because the average energy of AGN’s ionizing photons is higher than that of
young stars in SFGs, hence AGN can more efficiently heat the gas, therefore
boosting the [OIII]$\lambda 4363\,$line. We also found independent indications
of AGN activity in some high-redshift sources that were not previously
identified as AGN with the traditional diagnostics diagrams, but that are
placed in the AGN region of the diagnostic presented in this work. We note
that, conversely, low values of [OIII]$\lambda 4363\,$$\bigl{/}$H$\gamma\ $can
be associated either with SFGs or AGN excitation. We note that the fact that
strong auroral lines are often associated with AGN does not imply that they
cannot be used for direct metallicity measurements (provided that proper
ionization corrections are applied), but it does affect the calibration of
strong line metallicity diagnostics.

###### Key Words.:

galaxies: active – galaxies: high-redshift - galaxies: ISM

##  1 Introduction

Thanks to the successful launch of the James Webb Space Telescope (JWST;
Gardner et al., 2023; Rigby et al., 2023), we can now investigate with high
resolution and unprecedented sensitivity both the photometric and
spectroscopic properties of galaxies up to $z\sim 13$ (Curtis-Lake et al.,
2022; Robertson et al., 2022). Within this context, recent studies, exploiting
different kinds of JWST data, have revealed a large population of AGN at high
redshift (Kocevski et al., 2023; Übler et al., 2023a; Matthee et al., 2023;
Maiolino et al., 2024, 2023a; Greene et al., 2023; Bogdán et al., 2023;
Goulding et al., 2023; Kokorev et al., 2023; Furtak et al., 2023; Juodžbalis
et al., 2024; Scholtz et al., 2023a; Chisholm et al., 2024), suggesting a
higher AGN fraction at $z>3$ than previously expected. For instance, Yang et
al. (2023), taking advantage of the JWST-MIRI photometry of the Cosmic
Evolution Early Release Science Survey (CEERS; Finkelstein et al., 2022),
investigated the AGN population using spectral energy distribution (SED)
modelling, and found a black-hole accretion rate density (BHARD) at $z>3$
$\sim 0.5$ dex higher than what was expected from previous X-ray studies (Vito
et al., 2016, 2018). Also Maiolino et al. (2023a) and Harikane et al. (2023),
by selecting broad-line AGN (BLAGN, or TypeI) among NIRSpec spectra, have
found significant AGN excess at $z>4$ with respect to the AGN luminosity
functions derived using X-ray data (Giallongo et al., 2019).

The X-ray emission, produced in AGN by inverse Compton scattering of the UV
photons coming from the SMBH accretion disc, has been traditionally used to
select AGN. Most of the high-z AGN selected using JWST spectroscopy lie in
fields covered by some of the deepest X-ray observations ever performed, like
the 7Ms Chandra Deep Field South (CDFS) (Luo et al., 2017), Chandra Deep Field
North (CDFN) (Xue et al., 2016) and AEGIS-XD (Nandra et al., 2015). However,
most of these newly discovered high-z AGN lack any X-ray emission (e.g.
Maiolino et al. 2024 in prep. Mazzolari et al. 2024 in prep. Scholtz et al.,
2023a; Übler et al., 2023a), with only very few exceptions (Goulding et al.,
2023; Kovacs et al., 2024). This makes the picture even more difficult to
interpret and suggests different scenarios such as Compton thick obscuration
(Gilli et al., 2022) or intrinsic X-ray weakness (Simmonds et al., 2016).

The difficulties in selecting high-z AGN are even more important for TypeII or
narrow-line AGN (NLAGN). Different observational works have already
demonstrated that some of the traditional well-known AGN emission-line
diagnostic diagrams, like the BPT (Baldwin et al., 1981) diagram, are no
longer useful for the high-z Universe (Maiolino et al., 2023a; Scholtz et al.,
2023a; Harikane et al., 2023; Kocevski et al., 2023; Übler et al., 2023a;
Groves et al., 2006). Indeed, at high-z, the environments are systematically
more metal poor and the stellar populations are generally younger, increasing
the ionization parameter ($\rm\log U$) of star-forming galaxies (SFG). The
latter effect, on the one hand, makes SFGs move towards the AGN locus on the
BPT; on the other hand the lower metallicity of the Narrow Line Region (NLR)
of AGN makes its line ratios move towards and overlap with the star formation
locus on these traditional diagnostic diagrams. The detection of high-
ionization emission lines (like N V $\lambda 1242$, [Ne IV]$\lambda 2424$, [Ne
V]$\lambda 3426$) can still be considered as a safe tracer for the presence of
an AGN, but their detection is still difficult even in the deepest JWST
spectroscopic surveys like the JWST Advanced Deep Extragalactic Survey (JADES,
Bunker et al., 2023; Scholtz et al., 2023a).

To have a clear picture of the AGN properties and demographics at high-z, it
is therefore crucial to identify new AGN selection techniques. Thanks to the
large wavelength coverage ($\sim 1-5\mu$m) and sensitivity of JWST it is now
possible to have access to a plethora of rest-UV and optical emission lines in
high-z galaxies that were not accessible before. Among these lines, multiple
JWST spectra have revealed the detection of the [OIII]$\lambda 4363\,$line. In
the metal poor and highly ionized environments of many high-z galaxies,
[OIII]$\lambda 4363\,$is the strongest auroral line, i.e. a collisionally-
excited emission line generated from higher energy levels compared to the
typical nebular lines. Due to its faintness and proximity to the (generally)
stronger H$\gamma\ $line, it has been detected and studied almost only in
galaxies at $z<2$ , before the launch of JWST. Given its sensitivity to the
gas temperature, this line has been extensively used in the local Universe to
determine the metallicity of the interstellar medium (ISM) via the so-called
“direct temperature method”, which exploits the possibility of computing the
metallicity by fixing the gas electron temperature derived from the ratio
between the [OIII]$\lambda 4363\,$and the [OIII]$\lambda 5007\,$nebular lines
(see e.g. review by Maiolino & Mannucci, 2019).

However, Brinchmann (2023) first noticed that the anomalously high
[OIII]$\lambda 4363\,$emission in some of the early NIRSpec spectra at
z$\sim$8 could be indicating the presence of an AGN. More recently, Übler et
al. (2023b) have suggested that a strong ratio between this line and the
H$\gamma\ $line could be a possible indicator for the presence of an AGN. In
their work, by using JWST/NIRSpec Integral Field Spectroscopy (IFS)
observations, they noticed a spatial correspondence between the AGN broad
H$\beta\ $emission and the [OIII]$\lambda 4363\,$line emission. They also
mention the possibility that the enhanced [OIII]$\lambda 4363\,$emission in
correspondence of the AGN could be related to higher ISM temperatures driven
by the AGN activity.

An enhancement of the [OIII]$\lambda 4363\,$emission in the presence of an AGN
is qualitatively expected. Indeed, AGN ionizing photons are more energetic
than ionizing photons produced in star forming regions on average, hence they
can heat the gas more effectively by depositing a larger amount of energy per
unit photon. In this work, we will explore this scenario more quantitatively,
with the ultimate goal of exploring the possibility of selecting AGN using the
[OIII]$\lambda 4363\,$line both through an observational, empirical approach
and via photoionisation models. As a result, we propose three new AGN
diagnostic diagrams that can also be effective at high-z.

In Section 2 we present the local and high-z observational samples and the set
of photoionization models considered to test the validity of the new
diagnostic diagrams. In Section 3, we describe the new diagnostic diagrams
based on the line ratio [OIII]$\lambda 4363\,$/H$\gamma\ $and the demarcation
lines that can be used to separate the AGN population from SFGs. Lastly, in
Section 4, we provide interpretations on why these diagnostics work and we
show indications on the AGN nature of some newly identified AGN using these
diagnostics. Then, we also discuss the implications of this work on the
metallicity estimates based on the [OIII]$\lambda 4363\,$line and on the
strong-line metallicity diagnostics.

##  2 Methods

We selected several AGN and SFG samples and two different sets of
photoionization models to populate the diagnostic diagrams that we propose to
select AGN. All the observational samples described in the following Sec. 2.1
and Sec. 2.2 are listed in Table 1.

###  2.1 Low Redshift Samples

The sample of low-z SFGs is constituted by both normal SFGs and also by the
so-called local analogues of high-z galaxies, i.e. local galaxies that show
intrinsic properties similar to those observed in high-z galaxies, in terms of
metallicities and emission lines ratios (e.g. Izotov et al., 2018, 2019,
2021). The Sloan Digital Sky Survey (SDSS) DR7 (Abazajian et al., 2009) gives
the most numerous sample, mostly made of $z<0.7$ galaxies. Starting from the
whole catalog provided by MPA/JHU, we considered only sources with a $\rm
SNR>5$ in all the lines involved in the diagnostic diagrams we are going to
present: [OIII]$\lambda 4363\,$, [OIII]$\lambda 5007\,$, [Ne III]$\lambda
3869\ $, [OII]$\lambda 3727\,$, and H$\gamma\ $. Then we distinguished between
AGN and SFGs using the BPT diagram (Baldwin et al., 1981), taking the AGN
demarcation line provided by Kewley et al. (2001) and after considering a
further cut in SNR$>5$ also in the H$\alpha\ $, H$\beta\ $and [N II] lines.
Furthermore, we exclude from the selection AGN classified in the SDSS catalog
as BLAGN, since the catalog does not provide separately the NL and BL
components of the Balmer lines and our diagnostics rely on the NL emission
only. We also correct all the lines for dust attenuation using the Calzetti et
al. (2000) attenuation law. The final SDSS SFG (AGN) sample considered in this
work contains $\sim 2300$ ($\sim 800$) sources.

We also identify several samples of local analogues of high-z galaxies with
reported [OIII]$\lambda 4363\,$line fluxes. From Yang et al. (2017b) we took
the sample of so-called ”blueberries”, which are 40 dwarf starburst galaxies
with small sizes ($<1$kpc), very high ionization ([OIII]$\lambda
5007\,$/[OII]$\lambda 3727\,$$\sim 10-60$) and low metallicities ($\rm
7.1<12+\log(O/H)<7.8$). They were selected from the SDSS DR12 at $z<0.05$ and
followed up spectroscopically using MMT. We also consider 43 ”green pea”
galaxies from Yang et al. (2017a), i.e. nearby SFG with strong [OIII]$\lambda
5007\,$emission line, 2/3 of them also show a strong Ly$\alpha$ emission.

We include in the local analogue sample also a compilation of $\sim 490$ local
SFGs with low metallicities from different works: Izotov et al. (2006), Berg
et al. (2012), Izotov et al. (2018), Izotov et al. (2019), Pustilnik et al.
(2020), Pustilnik et al. (2021), Nakajima et al. (2022). All these samples
were selected considering local sources followed up with different telescopes
to be sensitive to the faint [OIII]$\lambda 4363\,$line, which allowed them to
measure the galaxies’ metallicity via the so-called direct method. The
metallicities of these sources span the range $12+\log(O/H)\sim 6.9-8.9$.

We also include the sample of 165 Extreme Emission Line Galaxies (EELG)
reported in Amorín et al. (2015), selected from the zCOSMOS survey (Lilly et
al., 2007) in the redshift range $0.11<z<0.93$. These are very compact
($r_{50}\sim 1.3$ kpc), low-mass ($\rm M_{*}\sim 10^{7}-10^{10}M_{\odot}$)
galaxies characterized by specific star formation rates (sSFR) above the main
sequence for star-forming galaxies of the same stellar mass and redshift.

The local AGN sample comprises TypeII AGN showing different features and
selected for different aims. The most numerous AGN sample is represented by
the SDSS sample mentioned above. We then include the Dors et al. (2020) TypeII
AGN sample, selected from SDSS DR7 using multiple AGN diagnostic diagrams to
guarantee a high purity and showing a range of metallicities between $\rm
8.0<12+\log(O/H)<9.2$.

We consider the sources described in Perna et al. (2017), selected from SDSS
$z<0.8$ spectra associated with X-ray emission and showing outflows, which can
be spectrally decomposed from the narrow and broad line region (BLR) emission.
We also added both TypeI and TypeII AGN from the S7 survey (Dopita et al.,
2015; Thomas et al., 2017), where 131 nearby AGN were primarily selected
because of their radio-detection and then followed-up using the WiFeS
instrument at high resolution ($R=7000$), allowing to resolve the faint
[OIII]$\lambda 4363\,$line. In Thomas et al. (2017) they provide the narrow
components of of the emission lines, already subtracting the broad component,
if present.

The 35 AGN selected in Armah et al. (2021) are TypeII AGN in the local
Universe ($z<0.06$) with quite high metallicities (the average metallicity is
$\rm 12+\log(O/H)=8.55$), that the authors select to study the neon to oxygen
abundance and its evolution with metallicity.

The last sample of local AGN is taken from the MOSFIRE Deep Evolution Field
(MOSDEF) survey (Kriek et al., 2015), which targeted the rest-frame optical
spectra of $\sim 1500$ H-band-selected Galaxies at $1.37<z<3.8$. We took the
full MOSDEF catalogue and selected only those sources with a [OIII]$\lambda
4363\,$line with a $\rm SNR>3$, and then we distinguished between AGN and SFG
using the BPT diagram. With these SNR cuts, we could identify only 3 sources,
all classified as AGN. This indicates the difficulty in detecting such a faint
line further than in the local Universe in the pre-JWST era.
For all the samples reported above we took dust-corrected fluxes from the
works reported in Tab. 1.

###  2.2 High Redshift Samples

The High-z samples of AGN and SFGs only comprise sources observed
spectroscopically using JWST, the only instrument with enough resolution,
sensitivity and IR-coverage to detect and disentangle the [OIII]$\lambda
4363\,$line at $z>3$. We considered the full publicly released sample of
sources with a medium resolution (MR, $R=1000$) spectrum in the CEERS (PID:
1345; Finkelstein et al., 2022) and JADES (PID: 1210; Bunker et al., 2023)
surveys. For JADES, we considered the emission line fluxes and the AGN
selection described in Scholtz et al. (2023a). Starting from a sample of 110
sources with medium resolution (MR) spectroscopy, a reliable redshift and
sufficient wavelength coverage, the authors performed a NLAGN selection using
multiple AGN rest-frame UV and optical diagnostic diagrams, finally selecting
28 reliable AGN candidates in the redshift range $1.8<z<9.4$. Of all the
initial JADES sources, only 19 show a [OIII]$\lambda 4363\,$detection, 7 of
which AGN. We followed the same approach for the CEERS program: starting from
313 MR spectra published by the CEERS collaboration, we selected 217 sources
with a reliable spectroscopic redshift and we finally selected 45 AGN
candidates (Mazzolari et al. 2024 in prep). For this work, we considered 32
sources with the [OIII]$\lambda 4363\,$line detected among the CEERS sample, 4
of which already classified as AGN. From the same programs, we also collect
the BLAGN identified by Maiolino et al. (2023a) and Harikane et al. (2023) at
$z>4.5$, carefully checking, and in case subtracting, the broad component of
the H$\gamma\ $line.

We further consider the three $z\sim 8$ galaxies in the galaxy cluster SMACS
J0723.3-7327, observed in the Early Release Observations (ERO) JWST programme
(Pontoppidan et al., 2022) and analysed in Curti et al. (2023a). Exploiting
the [OIII]$\lambda 4363\,$line detections in the NIRSpec spectra, the authors
measured metallicities ranging from extremely metal-poor ($\rm
12+\log(O/H)\sim 7$) to about one-third solar. Two of them were later
identified as AGN in Brinchmann (2023) due to the presence of the high-
ionization [Ne IV] emission line in one and based on the high ionization
parameter ($\rm\log U\sim 1$) of the other.

The sample of high-z sources presented in Nakajima et al. (2023) includes
galaxies from three different early JWST programs (CEERS, ERO, GLASS) at
$6<z<9$. We selected only the 10 sources with the [OIII]$\lambda 4363\,$line
detected, and we marked as BLAGN those that were selected based on their broad
H$\alpha\ $line in Harikane et al. (2023). For the latter, we subtract the
broad H$\gamma\ $component, if present.

We also include in the high-z AGN sample the BLAGN reported in Übler et al.
(2023a) at $z\sim 5.55$, whose JWST/NIRSpec Integral Field Spectrograph (IFS)
observation shows very high ionization lines and low metallicity (1/4 solar).
We also include both the components of the dual AGN candidate at $z\sim 7.15$
reported in Übler et al. (2023b). For this source, thanks to IFS spectroscopy,
the authors were able to study the displacement between the position of the
H$\beta\ $broad line region (BLR) and the strong [OIII]$\lambda
5007\,$emission line centroid, which were interpreted as the emissions coming
from two distinct sources 620pc apart from each other. Interestingly, the
authors also found an almost perfect alignment between the [OIII]$\lambda
4363\,$emission line peak and the peak of the H$\beta\ $BLR, suggesting a
correlation between the intensity of the [OIII]$\lambda 4363\,$line and the
presence of an AGN.

We also consider the BLAGN reported in Kokorev et al. (2023) at $z=8.50$ from
JWST UNCOVER Treasury survey (Labbe et al., 2021), showing a robust H$\beta\
$broad component and an unprecedented black hole to host galaxy mass of at
least $\sim 30\%$. An even larger black hole to host galaxy mass ratio was
more recently found in the $z=6.86$ BLAGN analyzed in Juodžbalis et al.
(2024), representing an extreme example of dormant SMBH at high-z. We also
considered the $z=7.04$ BLAGN, triply imaged and lensed by the galaxy cluster
Abell2744-QSO1, reported by Furtak et al. (2023). This source is probably
undergoing a phase of rapid SMBH growth and is also heavily obscured, since
the authors measured an $A_{V}\sim 3$ from the Balmer decrement.

Note that in all these TypeI AGN cases (actually most of them type 1.5-1.9)
the broad component of the Balmer lines is not detected in H$\gamma$, hence
the narrow component of H$\gamma$ can be easily measured.

We added to our high-z sample also the $z=6.1$ extreme SFG reported in Topping
et al. (2024). Using distict UV transitions the authors found an electron
density $\sim 10^{4}-10^{5}\rm cm^{-3}$, a metal poor ionized gas (i.e. $\rm
12+\log(O/H)=7.43$), and a $\rm\log U\sim-1$ based on the ratio [OIII]$\lambda
5007\,$/[OII]$\lambda 3727\,$. The authors ultimately classified this source
as a SFG, based on the ratios between the rest-UV emission lines, but they
also found some indications of possible AGN activity, like broad H$\alpha\
$component accounting for the $\sim 20\%$ of the line flux.

We finally included in the sample also GN-z11, an exceptionally luminous
galaxy at z = 10.6. Its JWST spectrum reveals the presence of an AGN through
the detection of semi-forbidden lines tracing very high densities
(inconsistent with the ISM but typical of the Broad Line Region), other
transitions typical of AGN, fast outflows, ionization cones, and a larger
Ly$\alpha$ halo, consistent with those seen in lower redshift quasars
(Maiolino et al., 2024, 2023b; Scholtz et al., 2023b).

Sample/References | Type | N included | Redshift |
---|---|---|---|---
SDSS DR7 Abazajian et al. (2009) | SFG & AGN | 3188 (2377/811) | $z_{med}\sim 0.1$ |
Yang et al. (2017b) | Blueberries Galaxies | 40 | $z<0.05$ |
Yang et al. (2017a) | Green Pea Galaxies | 43 | $z<0.05$ |
Low-OH compilation | Low-metallicities SFG | 489 | $z<0.1$ |
Amorín et al. (2015) | Extreme Emission Line Galaxies | 165 | $z<1$ |
Dors et al. (2020) | AGN | 26 | $z<0.05$ |
Perna et al. (2017) | AGN | 44 | $z<0.8$ |
S7 Survey Dopita et al. (2015) | AGN | 313 | $z<0.1$ |
Armah et al. (2021) | AGN | 35 | $z<0.06$ |
Kriek et al. (2015) | AGN | 3 | $z\sim 2.5$ |
CEERS (R1000) Finkelstein et al. (2022) | SFG & AGN | 32 (28/4) | $2.3<z<8.8$ |
JADES (R1000) Bunker et al. (2023) | SFG & AGN | 19 (12/7) | $0.66<z<9.43$ |
Curti et al. (2023b) | SFG & AGN | 3 | $z\sim 8$ |
Nakajima et al. (2023) | SFG & AGN | 10 | $6<z<9$ |
Topping et al. (2024) | SFG | 1 | $6.1$ |
Übler et al. (2023a) | AGN | 1 | $5.55$ |
Übler et al. (2023b) | AGN | 2 | $7.15$ |
Juodžbalis et al. (2024) | AGN | 1 | $6.86$ |
Kokorev et al. (2023) | AGN | 1 | $8.50$ |
Furtak et al. (2023) | AGN | 1 | $7.04$ |
Maiolino et al. (2024) | AGN | 1 | $10.6$ |
Table 1: Literature data of the low-z and high-z observational samples
presented in Sec 2.1 and Sec. 2.2 and used to test the effectiveness of the
new diagnostic diagrams.

###  2.3 Photoionization models

The observational data are compared with the results coming from two different
sets of photoionization models. The first set of models was initially
described in Gutkin et al. (2016) (for star formation), and Feltre et al.
(2016) (for narrow-line AGN emission) and updated with more recent stellar
spectra and with a better description of AGN cloud microturbulence (Mignoli et
al., 2019; Hirschmann et al., 2019; Vidal-García et al., 2024). The
photoionisation models are built using the CLOUDY code (version c13.03)
(Ferland et al., 2013) for star formation and AGN narrow-line regions,
assuming a wide range of parameters. In particular they consider gas
metallicities ($Z$, where $[O/H]=\log(Z/Z_{\odot})$) in the range
$10^{-4}<Z<0.07$, ionization parameters of the ionizing source $-4<\log U<-1$,
dust-to-metal ratios ($0.1<\xi<0.5$), hydrogen gas densities
($10^{2}<(n_{H}/cm^{-3})<10^{4}$), different carbon over oxygen abundances
(C/O) and two different initial mass functions (IMF) for star formation
models. For a complete list of the values of the different parameters we refer
to Table 1 of Feltre et al. (2016). These models include predictions for both
the AGN and SFG populations; the main feature differentiating between the AGN
and SFG models is the ionizing spectrum, i.e. the spectral energy distribution
(SED) of the incident ionizing radiation, the first showing a harder radiation
field.

As reported above, the photoionisation models of Gutkin et al. (2016) and
Feltre et al. (2016) cover a wide range of parameters and might actually
include physical conditions rarely found in the general populations of
galaxies and in the local analogs. For example, previous theoretical and
observational studies have suggested a correlation between the metallicity and
ionization parameter in SF regions (e.g., Dopita et al., 2006; Mingozzi et
al., 2020; Ji & Yan, 2022), which implies certain combinations of
metallicities and ionization parameters should be rarely, if not never,
observed in SF regions. Since the demarcation lines in the new diagrams we
propose rely on both the observational data and photoionisation models, we
performed the following check to make sure that the part of the models with
highly unphysical combinations of parameters does not impact our results.

We performed a likelihood analysis comparing the observed line ratios in our
samples of SFG and those predicted by the SFG of Gutkin et al. (2016), in
order to limit the inclusion of highly improbable or unphysical combinations
of the parameters. This method has been already adopted in making model-based
inferences for local galaxies (e.g., Blanc et al., 2015; Mingozzi et al.,
2020; Ji & Yan, 2022).

We considered a total of three sets of line ratios, including [OIII]$\lambda
4363\,$/H$\gamma\ $, [OIII]$\lambda 5007\,$/H$\beta\ $, and
[OII]$\lambda\lambda 3726,3728$/H$\gamma\ $. The inclusion of collisionally
excited lines from different ionization states of oxygen as well as
recombination lines of hydrogen helps break the degeneracy between the
metallicity and ionization parameter in the models. Following the formalism of
Blanc et al. (2015), we calculated the likelihood of each model given each
data point, and combined the likelihood with a flat prior in the log space
spanned by the metallicity and ionization parameter to obtain the posterior
distribution of these two parameters. From the posteriors, we obtained the
weighted average metallicities and ionization parameters for all SFGs. We then
selected a region in the metallicity-ionization parameter space not populated
by the posteriors of our sample of SFGs. We note that while the models of
Gutkin et al. (2016) and Feltre et al. (2016) have a large set of parameters,
the primary drivers of the variations in the predicted line ratios are the
metallicity and the ionization parameter. To verify this point, we repeated
the above calculation for models with different IMFs, densities, dust-to-gas
ratios, and carbon-to-oxygen abundance ratios, which generally resulted in a
similar distribution in the metallicity-ionization parameter space. With these
results, we adopted a conservative cut to select realistic SFG models based on
the $\rm 16^{th}$ percentile of the inferred distribution in the metallicity-
ionization parameter space that produces the highest boundary111The boundary
is determined by the common outer envelopes of SDSS SFGs and local analogs, as
these two samples differ significantly in number and occupy different regions
in the parameter space..

Our final cut removes the region of the SFG parameter space spanned by
$\rm\log(Z/Z_{\odot})$ and $\rm\log U$ according to the following relations:

$\rm\log U>-2.8\ \log(Z/Z_{\odot})-3.4,$ (1)

and

$\rm\log U>0.5\ \log(Z/Z_{\odot})-2.245.$ (2)

In Fig. 5, we plot our fiducial cut together with the inferred distribution of
the metallicities and ionization parameters for our selected SFGs based on
different sets of models. The cut basically removes regions having very high
ionization parameters and high metallicities at the same time. From a physical
point of view, due to the dependence of the ionizing spectra of young stellar
populations and the mechanical feedback from young stellar populations on the
metallicity, it is difficult to maintain a very high ionization parameter at a
high enough metallicity (Dopita et al., 2006; Carton et al., 2017). While
these physically motivated cuts do not impact the location of the AGN
photoionization models in the diagnostic diagrams, the parameter space
occupied by SFG models is partially reduced. We emphasize that 99.9% of the
considered SFGs lie outside of the excluded parameter space, even assuming
different values of $\xi$, $n_{H}$ and $C/O$.

We also compare the emission line ratios of the observational samples with the
models described in Nakajima & Maiolino (2022). Using CLOUDY, the authors
investigated the emission line ratios of SFG, AGN, Population III stars, and
Direct Collapse black holes (DCBHs) using the BPASS stellar population models
(Eldridge et al., 2017) and a range of other parameters (similar to those
considered in Feltre et al. (2016)) listed in Table 1 of Nakajima & Maiolino
(2022). In this work, we only consider the models for SFG and AGN, as the very
low metallicity values of the other two models ($Z<10^{-4}$) represent extreme
scenarios that are not representative of the sources we want to identify with
this analysis. In particular, for these models, the authors considered values
of $10^{-5}<Z<10^{-3}$ and values of $\log U$ from -0.5 to -3.5. Since SFG
models of Nakajima & Maiolino (2022) always cover regions of the diagnostic
diagrams already covered by the physically-motivated limited Gutkin et al.
(2016) and Feltre et al. (2016) SFG models, confirming the goodness of our
parameters cuts, we did not apply any cut in this case.

##  3 Results

In this section, we present the three new diagnostic diagrams proposed to
select AGN using the [OIII]$\lambda 4363\,$line. In these diagrams we combine
the auroral line, sensitive to temperature, with various other lines which are
sensitive to other additional properties of the ISM, such as ionization
parameter and shape of the ionizing SED, hence can further help to disentangle
the source of excitation. Additionally, these lines are accessible with JWST
at high redshift, and most of them are also close in wavelength, so little
affected by dust reddening.

###  3.1 [OIII]$\lambda 4363\,$$\bigr{/}$H$\gamma\ $vs [OIII]$\lambda
5007\,$$\bigr{/}$[OII]$\lambda 3727\,$

The diagnostic presented in Fig. 1 is based on the ratio between the
[OIII]$\lambda 4363\,$/H$\gamma\ $compared to [OIII]$\lambda
5007\,$/[OII]$\lambda 3727\,$. Hereafter [OII]$\lambda 3727\,$refers to the
sum of [OII]$\lambda\lambda 3726,3728$. In the left panel, where we plot the
observational samples described in Sec. 2, it is noticeable that the
distribution of AGN and SFGs from SDSS are well separated, with AGN occupying
a region with higher [OIII]$\lambda 4363\,$/H$\gamma\ $ratios compared to
normal SFGs. The distribution of the local analogues in this diagram seems to
extend the distribution of the local SDSS galaxies towards higher
[OIII]$\lambda 4363\,$/H$\gamma\ $and higher [OIII]$\lambda
5007\,$/[OII]$\lambda 3727\,$, along the diagonal of the plot. On the
contrary, the local AGN samples cover a wider area, mostly above the
distribution of local analogs and SFGs.

High-z sources are generally distributed in the upper part of the diagnostic
with respect to the local samples. We note that a non-negligible fraction of
high-z sources not classified as AGN lies in the region of the diagnostic
populated by local AGN, while there are also different high-z AGN falling in
the same region covered by local analogues.

On the right panel, we show the distribution of the photoionization models
presented in Sec. 2.3 according to the same line ratios. We see here that SFG
models are predicted to occupy a limited region of this diagram, whose upper
boundary almost perfectly corresponds with the distribution of SDSS SFG and
local analogs, while AGN are expected to cover both the region occupied by
SFGs and also the part of the diagnostic covered by the local AGN samples.
Therefore, SFG models are not able to reach the part of the diagram occupied
by local AGN, as expected also from the distribution of the observational
samples. The opposite is not true, since AGN can occupy the region populated
by SFGs, as observed for some local and high-z AGN. These considerations
remain true even considering the full grid of SFG models computed in Feltre et
al. (2016), as shown in Fig. 7 of the appendix.

As we will discuss, to a first approximation, the ratio [OIII]$\lambda
5007\,$/[OII]$\lambda 3727\,$traces the ionization parameter. The fact that,
at a given [OIII]$\lambda 5007\,$/[OII]$\lambda 3727\,$, AGN photoionization
can reach much higher values of [OIII]$\lambda 4363\,$/H$\gamma\ $than
photoionization from hot stars is likely a consequence, on average, of the
much higher energy of ionizing photons produced by AGN; therefore, at a given
ionization parameter, AGN can heat the gas much more effectively than hot
stars.

It is worth noting the position of some high-z sources in this diagnostic. The
AGN reported in Kokorev et al. (2023) shows an extremely strong [OIII]$\lambda
4363\,$line and it is placed in a region that is not covered by any model,
which probably suggests a very high electron temperature and $\rm\log U$, as
also pointed out in Kokorev et al. (2023). Another significant high-z AGN
falling in the AGN-dominated region of this diagnostic is the source
corresponding to the BLR centroid in the dual source described in Übler et al.
(2023b), where, in the IFS map, the peak of the broad-line emission spatially
corresponds to the peak of the [OIII]$\lambda 4363\,$line emission. The other
component, the [OIII]$\lambda 5007\,$centroid, showing a fainter
[OIII]$\lambda 4363\,$line emission, is close to the border between the local
analogues and the local AGN distribution, together with most of the high-z
BLAGN. The AGN region is also populated by different AGN selected among the
CEERS and JADES MR NIRSPEC spectra, but we also have different sources not
selected as AGN but falling in the AGN-dominated region, as pointed out above.
We demonstrate in Sec. 4 that, from the stack of their spectra, we have an
indication that some of these sources may host a BLAGN too.

An interesting source in this diagnostic is the SFG reported in Topping et al.
(2024). The source shows both strong [OIII]$\lambda 4363\,$/H$\gamma\ $and
[OIII]$\lambda 5007\,$/[OII]$\lambda 3727\,$emission line ratios, following
the extrapolation of the distribution of the local analogues, but at very high
values of the ionization parameters and at low metallicities (see Fig. 6), as
it is indeed found by the analysis performed in Topping et al. (2024). Looking
at the photoionization models, that region is not covered by the models of
Gutkin et al. (2016) and Feltre et al. (2016) (not even taking the whole
parameters grid), but is instead covered by some AGN models according to
Nakajima & Maiolino (2022), potentially suggesting that the main source of
ionising radiation in this extreme galaxy is ambiguous.

Finally, it is worth noting that the effect of dust reddening on the line
ratio [OIII]$\lambda 5007\,$/[OII]$\lambda 3727\,$, whose components can
suffer a different amount of attenuation due to their distance in wavelengths,
would move data points towards the right. This can potentially shift AGNs into
the SFG+AGN region under significant reddening, but not the opposite, i.e.
galaxies dominated by star formation would not contaminate the AGN-only
region.

###  3.2 [OIII]$\lambda 4363\,$$\bigr{/}$H$\gamma\ $vs [Ne III]$\lambda 3869\
$$\bigr{/}$[OII]$\lambda 3727\,$

The diagnostic presented in Fig. 2 is based on the ratio between the
[OIII]$\lambda 4363\,$/H$\gamma\ $compared to [Ne III]$\lambda 3869\
$/[OII]$\lambda 3727\,$. Not surprisingly, the distribution of the local and
high-z sources in this diagnostic is very similar to the one described in the
previous one. Indeed, [Ne III]$\lambda 3869\ $/[OII]$\lambda 3727\,$correlates
closely with [OIII]$\lambda 5007\,$/[OII]$\lambda 3727\,$both in AGN and SFGs,
as shown by Witstok et al. (2021), and as expected given that Neon and Oxygen
are both $\alpha$-elements and given that Ne+ and O+ have similar ionisation
potentials. The disadvantage of [Ne III]$\lambda 3869\ $relative to
[OIII]$\lambda 5007\,$is that it is typically much fainter. The advantage of
[Ne III]$\lambda 3869\ $is that it is at a shorter wavelength, hence
observable at higher redshift (out to z$\sim$12.3 with NIRSpec), and is closer
in wavelength to [OII]$\lambda 3727\,$, hence the [Ne III]$\lambda 3869\
$/[OII]$\lambda 3727\,$ratio is much less affected by dust reddening. As for
the previous diagram, normal SFGs and local analogs are distributed in the
lower-right part of the diagram, in the same region covered by SFG models,
while AGN can be distributed over a wider area, including the upper-left part
that cannot be reached by any SFG models or samples. As for the diagnostic
described in Sec. 3.1, the separation between the two populations holds even
considering the whole grid of the parameter space of the Gutkin et al. (2016)
star forming models, as presented in Fig. 6.

Even if the local analogues are less numerous than in Fig. 1 (simply because
there are fewer reported [Ne III]$\lambda 3869\ $fluxes in the literature),
their distributions seem to be again an extrapolation of the distribution of
local SFGs towards higher $\rm\log U$ and lower metallicities (see Fig. 6).

The sources that were clearly in the AGN region in the previous diagnostic are
still above the SFG distribution in this diagram, and the sources not
identified as AGN but lying in the AGN region of the diagnostic in Fig. 1 are
in the same region also in Fig. 2. In addition to this latter sample, in Fig.
2 there are few other sources from the CEERS and JADES samples that were not
in the diagnostic of Fig. 1, because the [OIII]$\lambda 5007\,$line was not
available, due to the presence of a detector gap. The AGN region is also
populated by three additional sources reported in Nakajima et al. (2023) and
not classified as AGN in the literature. Their location, however, is still
close to the border of the SFG distribution.

###  3.3 [OIII]$\lambda 4363\,$$\bigr{/}$H$\gamma\ $vs [OIII]$\lambda
5007\,$$\bigr{/}$[OIII]$\lambda 4363\,$

The third diagnostic is presented in Fig. 3 and is based on the ratio between
the [OIII]$\lambda 4363\,$/H$\gamma\ $compared to [OIII]$\lambda
5007\,$/[OIII]$\lambda 4363\,$. This diagnostic was reported also in Übler et
al. (2023b). Looking at the distribution of SDSS SFGs and AGN in the left
panel, we note that there is more overlap between the two populations with
respect to the previous diagnostics, but there is still a region,
characterized by high values of [OIII]$\lambda 4363\,$/H$\gamma\ $and high
[OIII]$\lambda 5007\,$/[OIII]$\lambda 4363\,$, that is populated only by local
AGN. The local analogs distribution seems again to be an extrapolation of the
trend of the SDSS SFG distribution. In this diagram, there are very few high-z
galaxies that are in the AGN-only region of the diagnostic, while most of them
are shifted towards lower values of [OIII]$\lambda 5007\,$/[OIII]$\lambda
4363\,$, even in a region that is not covered by the current photoionization
models, as it is evident from the right panel. This trend is emblematically
represented by the two extreme AGN reported in Kokorev et al. (2023) and in
Furtak et al. (2023). Given the strength of their [OIII]$\lambda 4363\,$line
emission, they are placed in the upper-left part of the plot in a region that
cannot be reproduced by any model. The same trend is shown by most of the
galaxies that were not selected as AGN, but that lie in the AGN region of the
two diagnostics presented above. The distribution of the SFG models matches
the distribution of the local galaxies, in particular in its boundary with the
AGN-only region, while AGN models cover both the region occupied by SFGs and
the one populated only by AGN samples.

In this diagram, assuming the whole parameters’ grid of the Feltre et al.
(2016) models, we would have a complete overlap between the SFG and AGN
models. However, the upper left part of the diagnostic would be covered by SFG
models characterized by high values of $\rm\log U$ and high metallicities, a
condition not observed in any of our galaxy samples (see Appendix A) and in
general in the literature (Kaasinen et al., 2018; Ji & Yan, 2022; Grasha et
al., 2022).

###  3.4 Defining the locus of AGN-only objects

For the three diagnostic diagrams described above, we were able to trace clear
demarcation lines providing a separation between the AGN-dominated regions and
the regions where SFGs and AGN can overlap. To do so, we considered both the
photoionization models and the observational sample distributions. We
ultimately define the demarcation lines in the form:

$Y=a\ X+b,$ (3)

where the values of the parameters $a$ and $b$ for each diagnostic are
reported in Table 2. It is important to note that the AGN selection based on
these demarcation lines is not a necessary condition for a source to be
selected as an AGN, but rather a sufficient condition: objects above the
demarcation line are AGN-dominated with high confidence, while objects below
the demarcation line can either be SFGs or AGN.

| O3Hg vs O32 | O3Hg vs Ne3O2 | O3Hg vs O33
---|---|---|---
a | 0.55 | 0.48 | -1.1
b | -0.95 | -0.42 | 1.47
Table 2: We report the values of the parameters $a$ and $b$ in Eq. 3 defining
the demarcation lines for each of the diagnostics considered.  Figure 1:
Diagnostic diagram showing the [OIII]$\lambda 4363\,$/H$\gamma\ $vs
[OIII]$\lambda 5007\,$/[OII]$\lambda 3727\,$lines ratios. The [OII]$\lambda
3727\,$is the sum of the doublet [OII]$\lambda\lambda 3726,3728$. Left panel:
plot of all the observational samples described in Sec. 2.1 and Sec. 2.2, with
colours and shapes as labelled. Red colours are used for AGN, blue for low-z
SFGs, and grey for high-z sources not classified as AGN. The red and blue
contours show the distribution at SDSS AGN and SFGs, respectively, marking
40%, 70%, and 90% of the populations. Right panel: same plot but showing the
AGN and SFG models computed by Gutkin et al. (2016), Feltre et al. (2016) and
Nakajima & Maiolino (2022), as described in Sec. 2.3. The black dashed line
indicates the demarcation defined in Sec. 3.4. Both the distribution of models
and of observational samples suggest that the dominant ionizing source in
galaxies above the demarcation line is AGN. The error bars reported in the
lower-right corner represent the median errors of the low redshift (left) and
high redshift (right) samples. The effect of dust attenuation on this
diagnostic moves sources towards the right, without contaminating the AGN-only
region with SFGs. Figure 2: Same as in Fig. 1, but for the lines ratios
[OIII]$\lambda 4363\,$/H$\gamma\ $vs [Ne III]$\lambda 3869\ $/[OII]$\lambda
3727\,$. Based on the distribution of observational samples and models, also
this diagnostic diagram identifies a region that can be populated only by AGN.
Figure 3: Same as in Fig. 1, but for the lines ratios [OIII]$\lambda 4363\,$/
H$\gamma\ $vs [OIII]$\lambda 5007\,$$\bigl{/}$[OIII]$\lambda 4363\,$. In the
left panel, the arrow on the AGN reported in Furtak et al. (2023) is for
visualization purposes, since it would be located at O3O3$\sim 0$. The cut in
SFG models described in Sec. 2.3 allowed us to identify an AGN-only region of
the diagnostic to the right of the black dashed line.

##  4 Discussion

By comparing the distribution of the observational samples with the
distribution of the photoionization models, we can explain why these
diagnostic diagrams can separate part of the AGN population from SFGs. In Fig
6 of the appendix, we plot the AGN and SFG models in the same diagnostic
diagrams presented in Sec. 3, but highlighting the variation of $\rm\log U$
and $Z$, i.e. the two main parameters affecting the distribution of the models
in these plots. We know that, since they are very local, SDSS AGN and SFG
samples are populated mostly by solar metallicity and moderate-to-low
ionization parameter sources, as shown also by the models occupying the same
region of these samples in Fig 6. On the contrary, local analogs are
characterized, on average, by lower metallicities and higher ionization
parameters (as high-z galaxies). Indeed, in Fig. 1 and Fig. 2, they fill the
regions that extend from the SDSS SFG towards higher [OIII]$\lambda
4363\,$/H$\gamma\ $and higher [OIII]$\lambda 5007\,$/[OII]$\lambda 3727\,$or
[Ne III]$\lambda 3869\ $/ [OII]$\lambda 3727\,$, in close agreement with the
trend of photoionization models for higher $\rm\log U$ and $Z$ in the same
diagnostics (see Fig. 6).

The reason for the presence, in particular in Fig. 1 and Fig. 2, of a region
that is populated only by AGN, at high [OIII]$\lambda 4363\,$/H$\gamma\
$values, can be explained in the following way. Since the [OIII]$\lambda
4363\,$line is a collisionally excited line generated from high energy levels,
its intensity relative to H$\gamma\ $is directly related to the temperature
and, secondly, to the metallicity of the ISM and ionization parameter. The
main difference between the AGN and SFGs is the SED of the incident ionizing
radiation. For AGN the ionizing radiation is given by the emission of the
accretion disk, while for SFGs it originates from young star clusters. The
former are characterized by a harder spectrum and, therefore, on average
photons are more energetic. Higher energy per ionizing photon, deposited into
the ISM, results into a higher effective heating, hence higher electron
temperature (hence higher [OIII]$\lambda 4363\,$) at a given ionization
parameter.

###  4.1 Stacking of the AGN candidates

To further demonstrate the effectiveness of the diagnostic diagrams proposed
in this work we explored potential tracers of AGN activity in those high-z
galaxies lying in the AGN region of the diagnostics diagrams in Fig. 1 and
Fig. 2 but not previously identified as AGN based on the selection performed
in Scholtz et al. (2023a) and Mazzolari et al. 2024 (in prep). Our stacking
procedure is as follows. We first shift the spectra to the rest-frame and
normalise them to the peak of the H$\alpha\ $line, that is available for all
the sources involved. This ensures that the bright sources are not dominating
the final stacks. We then resample each of the spectra to the common
wavelength grid. We verified that the wavelength grid does not impact our
conclusions by repeating the analysis using wavelength bins ranging from the
best ($\sim 100$km/s) to the worst ($\sim 150$km/s) resolution. Finally,
before stacking the spectra, we fit and subtract the continuum from the rest-
frame rebinned spectra. We performed an inverse variance stacking of 15
sources (11 from CEERS and 4 from JADES).
We first fit simultaneously the H$\alpha\ $and [N II] wavelength region using
only a single Gaussian profile per emission line, with common FWHM and
redshift for all three emission lines. We fixed the [N II] doublet ratio to be
3. We show the narrow emission line fit in the top left panel of Fig. 4. We
detect broad residuals around the H$\alpha\ $emission line. These broad
components can arise from either a broad line region or an outflow.
To model the broad wing in the H$\alpha\ $, we refit the H$\alpha\ $and [N II]
with two additional models - BLR model (by fitting one additional broad
Gaussian component to H$\alpha\ $only) and an outflow model by adding broad
Gaussian components to both H$\alpha\ $and [N II] with fixed FWHM and
centroid. Using the Bayesian Information Criterion (BIC) parameter (Liddle,
2007), we found that the BLR model is strongly preferred to the narrow-only
fit ($\Delta BIC=BIC_{H\alpha\ NL}-BIC_{H\alpha\ BL}=90$). The broad component
has an SNR of 10, making this a solid detection. On the contrary, the fit with
the outflow model did not return any broad component in the [N II] lines.
We then used the same approach with the H$\beta\ $and [O III] $\lambda\lambda
4959,5007\,$complex, to test the possibility that this broad component could
be associated with an outflow on larger scales, determining a broad
[OIII]$\lambda 5007\,$. In the fit shown in the last panel of Fig. 4, we
imposed the broad component of the [O III] $\lambda\lambda 4959,5007\,$to have
the same kinematics as the broad component of the H$\alpha\ $, i.e. same FWHM
and velocity offset. However, we didn’t find a significant broad line, and
therefore outflow, detection in the [O III] $\lambda\lambda 4959,5007\,$. In
particular, by imposing a broad outflow component to the [O III]
$\lambda\lambda 4959,5007\,$and comparing with the NL-only fit, it gives a
$\Delta BIC=BIC_{[OIII]\ NL}-BIC_{[OIII]\ BL}=4$, meaning that the two models
are almost equally preferred, hence a broad [OIII]$\lambda 5007\,$component is
not required.
The absence of a counterpart in [OIII]$\lambda 5007\,$strongly disfavours an
outflow origin of the broad component of H$\alpha\ $. Therefore, the detection
of a broad component in the H$\alpha\ $emission line indicates the presence of
an AGN BLR in the spectra of at least some of these galaxies, which is too
weak to be detected in the individual spectra. This indicates that galaxies
located in the upper region of our diagnostic diagrams, even if originally not
classified as AGN, actually may host faint TypeI AGN. This result provides
confidence that the diagnostic diagrams presented in this work are effective
in identifying AGN above our proposed demarcation lines, although they cannot
distinguish between AGN and SFGs below the demarcation line.

Figure 4: Different fits of the spectrum derived from the stack of the sources
not identified as AGN in the literature but lying in the AGN region of the
diagnostic diagrams presented in Sec. 3. The stacked spectrum in these plots
is resampled at the best resolution among those of the single spectra involved
(i.e. $\sim 100$km/s), but the results do not change considering the worst
resolution (i.e. $\sim 150$km/s). Top panels: fit of the H$\alpha\ $and [N II]
complex with only narrow component (left) and adding a broad component to the
H$\alpha\ $line (right). The global fit is presented in red, while the narrow
and broad H$\alpha\ $components are in green and yellow, respectively. In the
upper panel are shown the residuals of the fit compared with the distribution
of the $1\sigma$ (blue) and $3\sigma$ (red) errors on the fluxes. On the upper
right part of the plots are reported the FWHM and the velocity offset of the
different components considered in the fits. Lower panels: same as above but
for the H$\beta\ $and [O III] $\lambda\lambda 4959,5007\,$doublet complex. In
particular, on the right, we added to the fit of the [OIII]$\lambda 5007\,$a
broad component with the same FWHM and velocity offset as the broad component
of the H$\alpha\ $, but here it is not required by the fit.

###  4.2 Ruling out [FeII]$\lambda 4360$ contamination

We investigate the possibility that the [OIII]$\lambda 4363\,$emission line in
the high-z sources located in the AGN-only region of the diagnostics could be
contaminated by the [FeII]$\lambda 4360$ emission line (Curti et al., 2017;
Arellano-Córdova et al., 2022). Most of the high-z [OIII]$\lambda
4363\,$detections come from JWST MR spectra, whose resolution is not enough to
disentangle the two lines. Curti et al. (2017) showed that this contamination
is possible, but most likely in high-metallicity galaxies, while our high-z
sample is generally metal poor.
To investigate the possible iron contamination of the [OIII]$\lambda
4363\,$line, we consider the close [FeII]$\lambda 4288$ emission line, which
is isolated and originates from the same energy level as the [FeII]$\lambda
4360$. We ran a grid of Cloudy photoionization models over a wide range of
metallicities (with $\rm-2\leq log(Z/Z_{\odot})\leq 0.5$) and ionization
parameters (with $\rm-4\leq log(U)\leq-1$), and found that the flux ratio of
[FeII]$\lambda 4288$ / [FeII]$\lambda 4360$222Our calculation takes into
account 6 lines from $\rm Fe^{+}$ blended at 4288 Å and 8 lines from $\rm
Fe^{+}$ blended at 4360 Å, although the majority of the fluxes is contributed
by [FeII]$\lambda 4287.39$ and [FeII]$\lambda 4359.33$. is roughly a constant
of $\sim 1.25$. Therefore, we can use the [FeII]$\lambda 4288$ emission line
to reliably constrain the intensity of the [FeII]$\lambda 4360$. We did not
detect the [FeII]$\lambda 4288$ line in any of the high-z sources in the AGN-
only regions of the diagnostics. Also performing a weighted spectral stacking
of all the spectra in the same region, the [FeII]$\lambda 4288$ emission line
remains undetected. Therefore we conclude that any possible contribution to
the [OIII]$\lambda 4363\,$detections given by the [FeII]$\lambda 4360$
emission line is negligible.

###  4.3 Impact on metallicity estimates and strong-line metallicity
diagnostics

The [OIII]$\lambda 4363\,$emission line has been widely used in the literature
to derive metallicities via the so-called $T_{e}-$method (see Maiolino &
Mannucci, 2019, for a review), although Marconi et al. (2024) have shown that
these metallicities may be biased low. The metallicities inferred from the
$T_{e}-$method have then been used to calibrate the so-called strong-line
metallicity diagnostics, i.e. ratios of (optical) emission lines typically
much stronger than the auroral lines, which can be used to estimate the
metallicities on larger samples of galaxies. These calibrations have been
inferred both by using local samples of galaxies (e.g. Maiolino et al., 2008;
Bian et al., 2018; Curti et al., 2020) and, especially with JWST, at high
redshift (Sanders et al., 2023; Laseter et al., 2024).

In most of these past studies care was taken in using only SFGs, or assuming
that the gas ionization was dominated by young hot stars. The finding of
consistent populations of AGN at high redshift by recent JWST studies, and
especially our own finding in this paper that a significant number of strong
[OIII]$\lambda 4363\,$detections at high redshift are due to AGN heating, may
appear as concerning when using this transition to infer the gas metallicity.
However, the [OIII]$\lambda 4363\,$line (together with [OIII]$\lambda 5007\,$,
or [Ne III]$\lambda 3869\ $as a proxy when the latter is not available)
provides a measurement of the temperature in the O+2 zone, regardless of the
nature of the ionizing source. So its reliability remains unaffected in the
case of photoionization by AGN (modulo the potential biases discussed in
Marconi et al., 2024, which anyway are a concern also for SFGs). What is
obviously important is to properly use this information to infer the
metallicity, and in particular by applying the proper ionization correction
factors, which are different in the case of AGN and SFGs. In particular, once
the abundance of $O^{+2}/H^{+}$ is derived, one has to estimate the
contribution to the abundances from $O^{+3}/H^{+}$ and $O^{+}/H^{+}$. The
former is very difficult to assess, even in AGN, as there are no [OIV] strong
transitions in the wavelength ranges typically accessible. Yet, Dors et al.
(2020) estimate that, even in AGN, the contribution of $O^{+3}/H^{+}$ is
typically negligible. The contribution from $O^{+}/H^{+}$ is estimated from
the [OII]$\lambda 3727\,$doublet; however, the key issue is that the
temperature of the $O^{+}$ region is different from the temperature of the
$O^{+2}$ region probed by the [OIII]$\lambda 4363\,$line. If [OII] auroral
lines are not available (as in the vast majority of the cases at high
redshift), then one has to assume a relation between $T(O^{+})$ and
$T(O^{+2})$. This relation has been extensively explored in local SFGs,
however it cannot be applied to AGN. Indeed, assuming in AGN the same
temperature relation as for SFGs results into a large underestimation of the
metallicities (Dors et al., 2015). This issue is greatly mitigated once AGN-
bespoke relations between $T(O^{+})$ and $T(O^{+2})$ are adopted (Dors et al.,
2020).

In summary, the fact that at high-redshift [OIII]$\lambda 4363\,$is boosted by
AGN heating does not prevent it from being used for measuring the metallicity
in those galaxies, provided that the proper ionization correction factors are
adopted and, in particular, the adequate $T(O^{+})$ and $T(O^{+2})$ is adopted
(in absence of [OII] auroral lines).

However, the presence of AGN ionization, excitation and heating, certainly
affects the empirical calibrations of the strong line diagnostics. Indeed, at
a given gas temperature, the AGN radiation (with different ionizing shape and,
typically, higher ionization parameter) results in a different ionization
structure of the gas clouds, and also different collisional excitation rates
of the various transitions typically used in the strong line diagnostics
([OIII], [OII], [NeIII], [NII], [SII]). As a consequence, separate and
different empirical metallicity calibrations should be inferred for SFGs and
AGN host galaxies. This has been attempted in the local universe (Dors, 2021).
However, the same effort should be undertaken at high redshift, given the
large abundance of AGN. Inferring empirical strong-line metallicity
calibrations differentiating AGN and SFGs will be the focus of a separate
dedicated paper.

##  5 Conclusions

While JWST has revealed that some of the classical AGN diagnostics break down
at high redshift, it has also opened the opportunity to explore new
diagnostics. In this work, we studied the possibility of selecting AGN using
the [OIII]$\lambda 4363\,$auroral line, whose detection in large numbers of
high-z galaxies has become possible with JWST. In particular, we proposed
three new diagnostic diagrams, and three corresponding demarcation lines, that
allow to identify a large population of AGN from SFGs, by providing a
sufficient (but not necessary) condition to claim the presence of an AGN.

To demonstrate the effectiveness of these diagnostics, we used multiple
observational samples of local and high-z SFGs and AGN as well as
photoionization models from Feltre et al. (2016) and Nakajima & Maiolino
(2022).

We specifically proposed the following three diagnostic diagrams:

  * •

[OIII]$\lambda 4363\,$/H$\gamma$ vs [OIII]5007/[OII]3727. This is the most
thoroughly explored diagram from the empirical perspective (it has the largest
observational test sample), and can be used with NIRSpec out to z$<$9.4.

  * •

[OIII]$\lambda 4363\,$/H$\gamma$ vs [NeIII]3869/[OII]3727. This diagram is the
least sensitive to dust extinction given the wavelength proximity of both line
ratios. Additionally, it can be used with NIRSpec out to z$<$10.9.

  * •

[OIII]$\lambda 4363\,$/H$\gamma$ vs [OIII]5007/[OIII]$\lambda 4363\,$. Among
the three, this is the one with the smallest region where AGN can be safely
identified, but it is also relatively dust insensitive and requires the
detection of only three lines instead of four, making it applicable to larger
samples.

In each of these cases, we have provided a demarcation line above which (i.e.
with higher [OIII]$\lambda 4363\,$/H$\gamma$ values) objects can be safely
identified as hosting an AGN that dominates the nebular emission lines
excitation.

We have illustrated that at least some of few objects falling into the AGN-
only region and not previously identified as AGN, do show AGN signatures when
stacked. This further supports the effectiveness of the diagnostics. At the
same time, we did not find any indication of [FeII]$\lambda 4360$
contamination of the [OIII]$\lambda 4363\,$emission line in these sources,
which could in principle artificially increase the [OIII]$\lambda 4363\,$line
due to blending effects.

The physics behind these diagnostics is tightly linked to the primary property
of AGN, i.e. the hardness of their ionizing spectrum. The average energy of
AGN’s ionizing photons is much higher than the energy of hot, young stars in
SFGs. Therefore, at a given ionizing radiation field intensity, AGN photons
are more effective in heating the ionized gas than in SF regions, hence
yielding higher temperature, hence higher [OIII]$\lambda 4363\,$relative to
H$\gamma$ (which normalizes by the overall photoionizing radiation field).

We stress that being above the demarcation line in the proposed diagrams is a
sufficient but not necessary condition for an object to be identified as an
AGN. Galaxies located below such demarcation lines can be powered by star
formation and/or AGN.

Finally, we note that the fact that strong auroral lines are often associated
with AGN does not imply that they cannot be used for direct metallicity
measurements (provided that proper ionization corrections are applied), but it
does affect the calibration of strong line metallicity diagnostics, calling
for new, AGN-specific, calibrations of these diagrams.

## Acknowledgements

GM acknowledges useful conversations with Roberto Gilli and Marcella Brusa.
FDE, JS and RM acknowledge support by the Science and Technology Facilities
Council (STFC), from the ERC Advanced Grant 695671 “QUENCH”, and by the UKRI
Frontier Research grant RISEandFALL. RM also acknowledges funding from a
research professorship from the Royal Society. HÜ gratefully acknowledges
support by the Isaac Newton Trust and by the Kavli Foundation through a
Newton-Kavli Junior Fellowship.

## References

  * Abazajian et al. (2009) Abazajian, K. N., Adelman-McCarthy, J. K., Agüeros, M. A., et al. 2009, ApJS, 182, 543
  * Amorín et al. (2015) Amorín, R., Pérez-Montero, E., Contini, T., et al. 2015, A&A, 578, A105
  * Arellano-Córdova et al. (2022) Arellano-Córdova, K. Z., Berg, D. A., Chisholm, J., et al. 2022, ApJ, 940, L23
  * Armah et al. (2021) Armah, M., Dors, O. L., Aydar, C. P., et al. 2021, MNRAS, 508, 371
  * Baldwin et al. (1981) Baldwin, J. A., Phillips, M. M., & Terlevich, R. 1981, PASP, 93, 5
  * Berg et al. (2012) Berg, D. A., Skillman, E. D., Marble, A. R., et al. 2012, ApJ, 754, 98
  * Bian et al. (2018) Bian, F., Kewley, L. J., & Dopita, M. A. 2018, ApJ, 859, 175
  * Blanc et al. (2015) Blanc, G. A., Kewley, L., Vogt, F. P. A., & Dopita, M. A. 2015, ApJ, 798, 99
  * Bogdán et al. (2023) Bogdán, Á., Goulding, A. D., Natarajan, P., et al. 2023, Nature Astronomy [arXiv:2305.15458]
  * Brinchmann (2023) Brinchmann, J. 2023, MNRAS, 525, 2087
  * Bunker et al. (2023) Bunker, A. J., Cameron, A. J., Curtis-Lake, E., et al. 2023, arXiv e-prints, arXiv:2306.02467
  * Calzetti et al. (2000) Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682
  * Carton et al. (2017) Carton, D., Brinchmann, J., Shirazi, M., et al. 2017, MNRAS, 468, 2140
  * Chisholm et al. (2024) Chisholm, J., Berg, D. A., Endsley, R., et al. 2024, arXiv e-prints, arXiv:2402.18643
  * Curti et al. (2017) Curti, M., Cresci, G., Mannucci, F., et al. 2017, MNRAS, 465, 1384
  * Curti et al. (2023a) Curti, M., D’Eugenio, F., Carniani, S., et al. 2023a, MNRAS, 518, 425
  * Curti et al. (2020) Curti, M., Maiolino, R., Cirasuolo, M., et al. 2020, MNRAS, 492, 821
  * Curti et al. (2023b) Curti, M., Maiolino, R., Curtis-Lake, E., et al. 2023b, arXiv e-prints, arXiv:2304.08516
  * Curtis-Lake et al. (2022) Curtis-Lake, E., Carniani, S., Cameron, A., et al. 2022, arXiv e-prints, arXiv:2212.04568
  * Dopita et al. (2006) Dopita, M. A., Fischera, J., Sutherland, R. S., et al. 2006, ApJS, 167, 177
  * Dopita et al. (2015) Dopita, M. A., Shastri, P., Davies, R., et al. 2015, ApJS, 217, 12
  * Dors (2021) Dors, O. L. 2021, MNRAS, 507, 466
  * Dors et al. (2015) Dors, O. L., Cardaci, M. V., Hägele, G. F., et al. 2015, MNRAS, 453, 4102
  * Dors et al. (2020) Dors, O. L., Maiolino, R., Cardaci, M. V., et al. 2020, MNRAS, 496, 3209
  * Eldridge et al. (2017) Eldridge, J. J., Stanway, E. R., Xiao, L., et al. 2017, PASA, 34, e058
  * Feltre et al. (2016) Feltre, A., Charlot, S., & Gutkin, J. 2016, MNRAS, 456, 3354
  * Ferland et al. (2013) Ferland, G. J., Porter, R. L., van Hoof, P. A. M., et al. 2013, Rev. Mexicana Astron. Astrofis., 49, 137
  * Finkelstein et al. (2022) Finkelstein, S. L., Bagley, M. B., Arrabal Haro, P., et al. 2022, ApJ, 940, L55
  * Furtak et al. (2023) Furtak, L. J., Labbé, I., Zitrin, A., et al. 2023, arXiv e-prints, arXiv:2308.05735
  * Gardner et al. (2023) Gardner, J. P., Mather, J. C., Abbott, R., et al. 2023, PASP, 135, 068001
  * Giallongo et al. (2019) Giallongo, E., Grazian, A., Fiore, F., et al. 2019, ApJ, 884, 19
  * Gilli et al. (2022) Gilli, R., Norman, C., Calura, F., et al. 2022, A&A, 666, A17
  * Goulding et al. (2023) Goulding, A. D., Greene, J. E., Setton, D. J., et al. 2023, ApJ, 955, L24
  * Grasha et al. (2022) Grasha, K., Chen, Q. H., Battisti, A. J., et al. 2022, ApJ, 929, 118
  * Greene et al. (2023) Greene, J. E., Labbe, I., Goulding, A. D., et al. 2023, arXiv e-prints, arXiv:2309.05714
  * Groves et al. (2006) Groves, B. A., Heckman, T. M., & Kauffmann, G. 2006, MNRAS, 371, 1559
  * Gutkin et al. (2016) Gutkin, J., Charlot, S., & Bruzual, G. 2016, MNRAS, 462, 1757
  * Harikane et al. (2023) Harikane, Y., Zhang, Y., Nakajima, K., et al. 2023, arXiv e-prints, arXiv:2303.11946
  * Hirschmann et al. (2019) Hirschmann, M., Charlot, S., Feltre, A., et al. 2019, MNRAS, 487, 333
  * Izotov et al. (2019) Izotov, Y. I., Guseva, N. G., Fricke, K. J., & Henkel, C. 2019, A&A, 623, A40
  * Izotov et al. (2021) Izotov, Y. I., Guseva, N. G., Fricke, K. J., et al. 2021, A&A, 646, A138
  * Izotov et al. (2006) Izotov, Y. I., Stasińska, G., Meynet, G., Guseva, N. G., & Thuan, T. X. 2006, A&A, 448, 955
  * Izotov et al. (2018) Izotov, Y. I., Thuan, T. X., Guseva, N. G., & Liss, S. E. 2018, MNRAS, 473, 1956
  * Ji & Yan (2022) Ji, X. & Yan, R. 2022, A&A, 659, A112
  * Juodžbalis et al. (2024) Juodžbalis, I., Maiolino, R., Baker, W. M., et al. 2024, arXiv e-prints, arXiv:2403.03872
  * Kaasinen et al. (2018) Kaasinen, M., Kewley, L., Bian, F., et al. 2018, MNRAS, 477, 5568
  * Kewley et al. (2001) Kewley, L. J., Dopita, M. A., Sutherland, R. S., Heisler, C. A., & Trevena, J. 2001, ApJ, 556, 121
  * Kocevski et al. (2023) Kocevski, D. D., Onoue, M., Inayoshi, K., et al. 2023, arXiv e-prints, arXiv:2302.00012
  * Kokorev et al. (2023) Kokorev, V., Fujimoto, S., Labbe, I., et al. 2023, ApJ, 957, L7
  * Kovacs et al. (2024) Kovacs, O. E., Bogdan, A., Natarajan, P., et al. 2024, arXiv e-prints, arXiv:2403.14745
  * Kriek et al. (2015) Kriek, M., Shapley, A. E., Reddy, N. A., et al. 2015, ApJS, 218, 15
  * Labbe et al. (2021) Labbe, I., Bezanson, R., Atek, H., et al. 2021, UNCOVER: Ultra-deep NIRCam and NIRSpec Observations Before the Epoch of Reionization, JWST Proposal. Cycle 1, ID. #2561
  * Laseter et al. (2024) Laseter, I. H., Maseda, M. V., Curti, M., et al. 2024, A&A, 681, A70
  * Liddle (2007) Liddle, A. R. 2007, MNRAS, 377, L74
  * Lilly et al. (2007) Lilly, S. J., Le Fèvre, O., Renzini, A., et al. 2007, ApJS, 172, 70
  * Luo et al. (2017) Luo, B., Brandt, W. N., Xue, Y. Q., et al. 2017, ApJS, 228, 2
  * Maiolino & Mannucci (2019) Maiolino, R. & Mannucci, F. 2019, A&A Rev., 27, 3
  * Maiolino et al. (2008) Maiolino, R., Nagao, T., Grazian, A., et al. 2008, A&A, 488, 463
  * Maiolino et al. (2023a) Maiolino, R., Scholtz, J., Curtis-Lake, E., et al. 2023a, arXiv e-prints, arXiv:2308.01230
  * Maiolino et al. (2024) Maiolino, R., Scholtz, J., Witstok, J., et al. 2024, Nature, 627, 59
  * Maiolino et al. (2023b) Maiolino, R., Uebler, H., Perna, M., et al. 2023b, arXiv e-prints, arXiv:2306.00953
  * Marconi et al. (2024) Marconi, A., Amiri, A., Feltre, A., et al. 2024, arXiv e-prints, arXiv:2401.13028
  * Matthee et al. (2023) Matthee, J., Naidu, R. P., Brammer, G., et al. 2023, arXiv e-prints, arXiv:2306.05448
  * Mignoli et al. (2019) Mignoli, M., Feltre, A., Bongiorno, A., et al. 2019, A&A, 626, A9
  * Mingozzi et al. (2020) Mingozzi, M., Belfiore, F., Cresci, G., et al. 2020, A&A, 636, A42
  * Nakajima & Maiolino (2022) Nakajima, K. & Maiolino, R. 2022, MNRAS, 513, 5134
  * Nakajima et al. (2023) Nakajima, K., Ouchi, M., Isobe, Y., et al. 2023, arXiv e-prints, arXiv:2301.12825
  * Nakajima et al. (2022) Nakajima, K., Ouchi, M., Xu, Y., et al. 2022, ApJS, 262, 3
  * Nandra et al. (2015) Nandra, K., Laird, E. S., Aird, J. A., et al. 2015, ApJS, 220, 10
  * Perna et al. (2017) Perna, M., Lanzuisi, G., Brusa, M., Cresci, G., & Mignoli, M. 2017, A&A, 606, A96
  * Pontoppidan et al. (2022) Pontoppidan, K. M., Barrientes, J., Blome, C., et al. 2022, ApJ, 936, L14
  * Pustilnik et al. (2021) Pustilnik, S. A., Egorova, E. S., Kniazev, A. Y., et al. 2021, MNRAS, 507, 944
  * Pustilnik et al. (2020) Pustilnik, S. A., Kniazev, A. Y., Perepelitsyna, Y. A., & Egorova, E. S. 2020, MNRAS, 493, 830
  * Rigby et al. (2023) Rigby, J., Perrin, M., McElwain, M., et al. 2023, PASP, 135, 048001
  * Robertson et al. (2022) Robertson, B. E., Tacchella, S., Johnson, B. D., et al. 2022, arXiv e-prints, arXiv:2212.04480
  * Sanders et al. (2023) Sanders, R. L., Shapley, A. E., Topping, M. W., Reddy, N. A., & Brammer, G. B. 2023, arXiv e-prints, arXiv:2301.06696
  * Scholtz et al. (2023a) Scholtz, J., Maiolino, R., D’Eugenio, F., et al. 2023a, arXiv e-prints, arXiv:2311.18731
  * Scholtz et al. (2023b) Scholtz, J., Witten, C., Laporte, N., et al. 2023b, arXiv e-prints, arXiv:2306.09142
  * Simmonds et al. (2016) Simmonds, C., Bauer, F. E., Thuan, T. X., et al. 2016, A&A, 596, A64
  * Thomas et al. (2017) Thomas, A. D., Dopita, M. A., Shastri, P., et al. 2017, ApJS, 232, 11
  * Topping et al. (2024) Topping, M. W., Stark, D. P., Senchyna, P., et al. 2024, MNRAS, 529, 3301
  * Übler et al. (2023a) Übler, H., Maiolino, R., Curtis-Lake, E., et al. 2023a, A&A, 677, A145
  * Übler et al. (2023b) Übler, H., Maiolino, R., Pérez-González, P. G., et al. 2023b, arXiv e-prints, arXiv:2312.03589
  * Vidal-García et al. (2024) Vidal-García, A., Plat, A., Curtis-Lake, E., et al. 2024, MNRAS, 527, 7217
  * Vito et al. (2018) Vito, F., Brandt, W. N., Yang, G., et al. 2018, MNRAS, 473, 2378
  * Vito et al. (2016) Vito, F., Gilli, R., Vignali, C., et al. 2016, MNRAS, 463, 348
  * Witstok et al. (2021) Witstok, J., Smit, R., Maiolino, R., et al. 2021, MNRAS, 508, 1686
  * Xue et al. (2016) Xue, Y. Q., Luo, B., Brandt, W. N., et al. 2016, ApJS, 224, 15
  * Yang et al. (2023) Yang, G., Caputi, K. I., Papovich, C., et al. 2023, ApJ, 950, L5
  * Yang et al. (2017a) Yang, H., Malhotra, S., Gronke, M., et al. 2017a, ApJ, 844, 171
  * Yang et al. (2017b) Yang, H., Malhotra, S., Rhoads, J. E., & Wang, J. 2017b, ApJ, 847, 38

##  Appendix A Cut in photoionization models

In Fig. 5 we show the distribution of the SDSS SFG and local analogues samples
in the $\rm\log U$ vs $\rm\log OH$ plane as derived following the likelihood
procedure with the Gutkin et al. (2016) models described in Sec. 2.3. Even
considering different values of the C/O, $\xi$, and $n_{H}$ parameters in the
models, the chosen excluded region contains less than $0.1$% of the sources
(located very close to the borders).

Figure 5: Distribution of the SDSS SFG sample (circles) and local analogues
(triangles) in the $\rm\log U$ vs $\rm\log OH$ plane, according to the
likelihood procedure with the Feltre et al. (2016) models described in Sec.
2.3. In each panel, we also show, with a red dashed line, the region of the
parameter space that we decided to exclude, a posteriori, from the same
photoionization models. The upper panels show the variation of the
distribution of the sources by considering models with different electron
densities $\rm n_{H}$, in the middle panels different dust-to-metal ratios
$\xi$, and in the bottom panels different $C/O$ abundance. The excluded region
is never significantly populated by any source in all the panels.

##  Appendix B Photoionization models properties

In Fig. 6 we show the distribution of the SFG and AGN photoionization models
from Gutkin et al. (2016), Feltre et al. (2016), and Nakajima & Maiolino
(2022) in the three diagnostic diagrams presented in Sec. 3. In particular,
the plots highlight the variation of $\rm\log U$ and $Z$ of the models across
the diagnostics.

Figure 6: Distribution of the photoionization models computed in Feltre et al.
(2016) (circles) and Nakajima & Maiolino (2022) (stars) with respect to the
line ratios reported in Fig. 1 (top panel),Fig. 2 (central panel), Fig. 3
(bottom panel). The photoionization models for SFG are reported on the right,
while the AGN models are on the left. The points are colour-coded according to
their metallicity (the stronger the colour, the higher the metallicity), and
the marker size depends on the ionization parameter (the larger the marker,
the higher the ionization parameter). The maximum and minimum values of $Z$
and $\rm\log U$ for the two classes of models are reported in the top part of
the top panels.

##  Appendix C Diagnostic diagrams including all photoionization models

The demarcation lines presented in Sec. 3.4 were defined considering both the
distribution of the observational samples and of photoionization models with
the cut described in Sec. 2.3 and shown in Fig. 5. If we consider the entire
grid of parameters for the models of Gutkin et al. (2016) and Feltre et al.
(2016), we obtain the distributions shown in Fig. 7. In this case, the
demarcation lines for the O3Hg-O32 and O3Hg-Ne3O2 diagnostic plots still hold,
except for a very small fraction of SFG models populating the AGN-only part of
the two diagnostics. On the contrary, in the O3Hg-O33 diagnostic plot, there
is an almost complete superposition between the SFG and AGN models. Even if
this superposition is determined by highly unlikely combinations of parameters
in the models, we cannot completely rule out the possibility that the AGN-only
region of the O3Hg-O33 diagnostic diagram could actually be partially
populated also by SFG.

Figure 7: From left to right, same plots reported in the right panels of Fig.
1, Fig. 2 and Fig. 3, but showing all the models computed in Gutkin et al.
(2016), Feltre et al. (2016) and Nakajima & Maiolino (2022), i.e. without the
cut described in Sec. 2.3. We note that the demarcation lines for the O3Hg-O32
and O3Hg-Ne3O2 still hold, while there is an almost complete superposition
between SFG and AGN models in the O3Hg-O33 diagnostic diagram.