MNRAS 000, 1–24 (2018) Preprint 30 November 2023 Compiled using MNRAS LATEX style file v3.0
beagle-agn I: Simultaneous constraints on the properties of
gas in star-forming and AGN narrow-line regions in galaxies
A. Vidal-García1,2?, A. Plat3†, E. Curtis-Lake4‡, A. Feltre5, M. Hirschmann6,7, J. Chevallard8
and S. Charlot9
1 Observatorio Astronómico Nacional, C/ Alfonso XII 3, 28014 Madrid, Spain
2 École Normale Supérieure, CNRS, UMR 8023, Université PSL, Sorbonne Université, Université de Paris, F-75005 Paris, France
3 Steward Observatory, 933 N. Cherry Avenue, University of Arizona, Tucson, AZ 85721, USA
4 Centre for Astrophysics Research, Department of Physics, Astronomy and Mathematics, University of Hertfordshire, Hatfield AL10
9AB, UK
5 INAF - Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Via P. Gobetti 93/3, 40129 Bologna, Italy
6 Institute of Physics, GalSpec Laboratory, Ecole Polytechnique Federale de Lausanne, Observatoire de Sauverny, Chemin Pegasi 51,
1290 Versoix, Switzerland
7 INAF, Osservatorio Astronomico di Trieste, Via G. B. Tiepolo 11, 34131 Trieste, Italy
8 Sub-department of Astrophysics, Department of Physics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford OX1
3RH, UK
9 Sorbonne Université, CNRS, UMR7095, Institut d’Astrophysique de Paris, F-75014, Paris, France
ABSTRACT
We present the addition of nebular emission from the narrow-line regions (NLR) surrounding active galactic nuclei
(AGN) to beagle (BayEsian Analysis of GaLaxy sEds). Using a set of idealised spectra, we fit to a set of observables
(emission-line ratios and fluxes) and test the retrieval of different physical parameters. We find that fitting to stan-
dard diagnostic-line ratios plus [O ii]λ3726, λ3729/[O iii]λ5007, Hβ/Hα, [O i]λ6300/[O ii]λ3726, λ3729 and Hα flux,
degeneracies remain between dust-to-metal mass ratio (ξnlrd ) and ionisation parameter (U
nlr
s ) in the NLR gas, and
between slope of the ionizing radiation (αpl, characterising the emission from the accretion disc around the central
black hole) and total accretion-disc luminosity (Lacc). Since these degeneracies bias the retrieval of other parameters
even at maximal signal-to-noise ratio (S/N), without additional observables, we suggest fixing αpl and dust-to-metal
mass ratios in both NLR and H ii regions. We explore the S/N in Hβ required for un-biased estimates of physical
parameters, finding that S/N(Hβ) ∼ 10 is sufficient to identify a NLR contribution, but that higher S/N is required
for un-biased parameter retrieval (∼ 20 for NLR-dominated systems, ∼ 30 for objects with approximately-equal Hβ
contributions from NLR and H ii regions). We also compare the predictions of our models for different line ratios to
previously-published models and data. By adding [He ii]λ4686-line measurements to a set of published line fluxes for
a sample of 463 AGN NLR, we show that our models with −4 < logUnlrs < −1.5 can account for the full range of
observed AGN properties in the local Universe.
Key words: galaxies:ISM, galaxies:active, HII regions.
1 INTRODUCTION
In searching for obscured (termed type-2) Active Galac-
tic Nuclei (AGNs), astronomers often exploit ‘BPT’ di-
agrams, which are named after the founding paper by
Baldwin et al. (1981). The scheme proposed in Bald-
win et al. (1981), which employs ratios between strong
emission line fluxes, was later updated by Veilleux
& Osterbrock (1987), who settled on three diagnos-
tic diagrams: [O iii]λ5007/Hβ versus [S ii]λ6717, λ6731/Hα,
[O iii]λ5007/Hβ versus [N ii]λ6584/Hα and [O iii]λ5007/Hβ
? Email: a.vidal@oan.es
† Email: plat@iap.fr
‡ Email: e.curtis-lake@herts.ac.uk
versus [O i]λ6300/Hα. These diagrams provide clear separa-
tion of the different ionizing sources powering the line emis-
sion. However, at low metallicities, the regions within the
diagnostic diagrams occupied by star-formation and AGN-
driven line emission increasingly overlap, muddying the clas-
sification (Groves et al. 2006; Feltre et al. 2016; Hirschmann
et al. 2019). Therefore, at low masses, or higher redshifts,
the standard BPT diagrams will be less helpful to distinguish
AGN from star-forming activity. It is therefore important to
develop new tools to probe the impact of AGNs on galaxy
evolution with the new era of large spectroscopic surveys at
high redshifts (e.g. with the James Webb Space Telescope).
Beyond simple classification of galaxy-wide star-formation
and AGN-driven line emission, there has been significant
© 2018 The Authors
2 Vidal-García et al.
progress in deriving physical properties of the star-forming
regions themselves (e.g. Kewley et al. 2019; Maiolino & Man-
nucci 2019, and references therein). However, progress on the
side of AGN gas properties is more limited. The first emission-
line calibrations to derive oxygen abundances of type-2 AGNs
were supplied by Storchi-Bergmann et al. (1998). The sub-
sequent efforts have proceeded along three avenues: further
research of observables to derive physical properties (e.g. from
rest-frame UV lines, Dors et al. 2014; or the rest-frame op-
tical, Castro et al. 2017); comparison of emission-line ratios
with photoionization models, often searching for a minimum
χ2 solution (e.g. for rest-frame UV lines, Nagao et al. 2006b;
Matsuoka et al. 2009, 2018); and full Bayesian parameter es-
timation (Pérez-Montero et al. 2019; Thomas et al. 2018b;
Mignoli et al. 2019).
Dors et al. (2020a) used a large sample of type-2 AGNs
selected from the Sloan Digital Sky Survey (SDSS, York et al.
2000) to provide the first comparison of methods for deriving
oxygen abundances of type-2 AGNs. These methods include
the calibrations of Storchi-Bergmann et al. (1998), Castro
et al. (2017) and the Bayesian-style method employed by the
Hii-Chi-mistry code (Pérez-Montero et al. 2019). They find
fairly poor agreement between the different methods as well
as no significant trend between a galaxy’s stellar mass and
the metallicity of the narrow-line region (NLR) surrounding
the AGN.
Dors et al. (2020a) also test a ‘direct’ method, that
estimates the electron temperature (Te) within the high-
ionization zone using the RO3 = [Oiii]λ4959, λ5007/λ4363
ratio. They find that their implemented ‘direct’ method sig-
nificantly under-estimates the derived metallicities compared
to the other methods. In Dors et al. (2015), they inferred
that the ‘direct’ metallicities were significantly lower than
extrapolated radial metallicity gradients derived from gas-
phase abundances. Dors et al. (2020b) subsequently identify
the cause of the discrepancy between various calibrations and
the ‘direct’ method. They update the ‘direct’ method and
demonstrate that the previous discrepancies with other cali-
brations are much reduced. Further to this, Dors (2021) pro-
vide the first strong-line calibration against ‘direct’ metallic-
ity estimates for AGNs.
Despite large discrepancies in different oxygen abundance
estimates, there have been efforts to characterise the type-2
AGN population by studying line emission. The first studies
investigating the statistics of gas metallicity in the narrow-
line region of AGNs mostly indicated a lack of evolution with
redshift (Matsuoka et al. 2009; Dors et al. 2014), although
they did reveal a luminosity-metallicity relation, where less
luminous AGNs (characterised by the [He ii]λ1640 luminosity
in both these studies) display lower metallicities than their
bright counterparts. In a later study, Mignoli et al. (2019)
investigated the properties of the NLR gas of a sample of
Type-2 AGNs selected in a homogeneous way and found sig-
nificant evolution with redshift. This study employed a wider
selection of lines, as well as updated photoionization models,
which better reproduce the [Nv]λ1240 line without resort-
ing to the very high metallicities evoked in previous stud-
ies of quasars (e.g., Hamann & Ferland 1993; Dietrich et al.
2003; Nagao et al. 2006a) and narrow-line Seyfert-I galax-
ies (Shemmer & Netzer 2002). Recently, do Nascimento et al.
(2022) used strong-line calibrations to study the NLR in type-
2 AGNs within the MaNGA survey (Mapping Nearby Galax-
ies at Apache Point Observatory; Bundy et al. 2015), finding
that the central NLR region typically has lower metallicity
than the surrounding HII regions. They posit that this may
be due to accretion of metal-poor gas at the centre of these
galaxies, which is feeding the central black hole.
One limitation of all the methods mentioned thus far is
the lack of any accounting for H ii region contribution to the
line emission. Thomas et al. (2018b) explicitly investigate
the degree of H ii and NLR mixing in the line emission in
a sample of SDSS galaxies using their Bayesian code, Neb-
ularBayes. They demonstrate that even in regions of the
BPT diagram where AGNs are determined to be cleanly se-
lected, the Balmer lines can have a significant contribution
from H ii regions ionized by young stars (see also Agostino
et al. 2021). Further to this, Thomas et al. (2019) measure the
mass-metallicity relation of type-2 AGNs finding a moderate
increase in oxygen abundance with increasing stellar mass.
With this paper, we present the incorporation of the Feltre
et al. (2016) NLR models into beagle (BayEsian Analysis of
GaLaxy sEds), a tool to model and interpret galaxy spectral
energy distributions (Chevallard & Charlot 2016) (Section 2).
This addition allows the mixing of line emission from young
stellar birth clouds (Section 2.1) with that from the NLR of
type-2 AGNs (Section 2.2). beagle also self-consistently in-
cludes stellar emission and attenuation by interstellar dust
(Section 2.3), which are not explicitly modelled in Nebu-
larBayes. In Section 3, we take a pedagogical approach to
defining what parameters of our model can be constrained
by fitting a given set of observables (Section 3.2) in idealised
spectra (Section 3.1). Using such spectra, we quantify the
S/N required to constrain H ii-region and NLR-gas parame-
ters for different NLR contributions to the total Hβ flux of
a galaxy (Sections 3.3 and 3.4). With this work, we focus
on rest-frame optical observables, though beagle can also
be used to study emission lines from the rest-frame UV. Fi-
nally, in Section 4, we compare the results of our model with
those obtained using previously published methods. This last
section of our work focuses, therefore, on comparing the pre-
dictions of the model grid included in beagle to different
ones in the literature. We expand on the comparison work of
Dors et al. (2020a) by explaining how different methods will
derive different oxygen abundances (Section 4.1). In particu-
lar, in Section 4.2 we compare several emission-line ratios and
several free parameters in the models and in a set of type-
2 AGN observations. In order to explain different trends for
the ionization parameter, in Section 4.3 we present new mea-
surements of the [He ii]λ4686 fluxes in a few hundred type-2
AGNs from DR7 SDSS and we compare the measured data
to the fluxes predicted by the F16 models. We also show the
predictions of the fluxes for sulfur and nitrogen lines in Sec-
tion 4.4, and end this discussion (Section 4.5) by comparing
the Dors (2021) empirical Te-based 12 + log(O/H) calibra-
tion to the different model grids. The paper ends in Section 5
with a summary of the main findings of this work. This is
the first in a series of three papers, to be followed by a paper
on fitting of a sample of type-2 AGNs with beagle, and a
study of the extent to which line emission from shocks and
post-AGB stars may affect our inferences.
MNRAS 000, 1–24 (2018)
beagle-agn 3
2 MODELLING THE EMISSION FROM STARS,
H ii REGIONS AND AGN NLRS IN beagle
beagle (Chevallard & Charlot 2016, hereafter CC16) is a
Bayesian SED-fitting code, which allows efficient exploration
of a wide grid of physical parameters affecting the light emit-
ted by a galaxy. The code employs MultiNest (a Bayesian
analysis tool based on the Nested Sampling algorithm of
Skilling et al. 2006; see Feroz et al. 2009) to sample from the
posterior probability distributions of physical parameters. We
refer the reader to CC16 for more details.
A main feature of beagle is the incorporation of physically
consistent stellar continuum and nebular emission models,
which trace the production of starlight and its transmission
through the interstellar medium self-consistently. This paper
is concerned with the extension of beagle to allow for the
interpretation of mixed emission-line signatures of AGNs and
stars. In Section 2.1 below, we briefly review the modelling
of the emission from stars and H ii regions in beagle, while
in Section 2.2, we describe our incorporation of the emission
from AGN narrow-line regions.
2.1 Emission from stars and H ii regions
beagle employs the latest version of the Bruzual & Char-
lot (2003) stellar population synthesis models to compute the
emission from stars. This version of the code computes the
integrated spectrum between 5.6Å and 3.6 cm of a stellar
population with an age between 0.01Myr and 13.8Gyr and a
metallicity in the range 0.0001 6 Z 6 0.04 using evolutionary
tracks which include stars with masses up to 350 M popu-
lating the Chabrier Initial Mass Function (see Vidal-García
et al. 2017, for more modeling details).
The line and continuum emission of H ii regions ionized
by stars younger than 10Myr is computed following the pre-
scription of Gutkin et al. (2016, see also Charlot & Longhetti
2001) for a grid of gas parameters: the metallicity of the H ii-
region gas, Zhiigas (noted Zism in Gutkin et al. 2016); the ion-
ization parameter, logUs1; and the dust-to-metal mass ratio
within the H ii regions, ξd. Nitrogen abundances are scaled
to oxygen abundances according to the formula:
N/H ' 0.41 O/H
[
10−1.6 + 10(2.33+logO/H)
]
, (1)
which well matches the relation between N/O and
12 + log(O/H) in observed H ii regions (see equation 11 and
figure 1 of Gutkin et al. 2016).
For simplicity, we fix here the carbon-to-oxygen ratio of
the H ii regions to solar, (C/O) = 0.44, and the hydrogen
density to nh = 100 cm−3. In the stellar population mod-
els, we fix the upper limit of the stellar initial mass function
(IMF) to 100 M. An important feature of this grid is the
coupling of stellar light to its transmission through the inter-
stellar medium (ISM), via not only the nebular emission, but
also the attenuation by dust in the H ii regions themselves
(see Section 2.3).
1 Note that logUs is defined as the ionization parameter at the
Strömgrem radius, which differs from the volume-averaged ioniza-
tion parameter, < U > according to < U >= 9/4Us
12.5 13.0 13.5 14.0 14.5 15.0 15.5 16.0 16.5 17.0
log(  / Hz)
22
21
20
19
18
17
16
15
lo
g(
L
 / 
er
g 
s
1  H
z
1 )
, a
rb
itr
ar
y 
sc
al
in
g
13.6eV
F16, PL = 1.7
F16, full range of PL
Mathews & Ferland 87
Perez-Montero+19 ox = 0.8
Perez-Montero+19 ox = 1.2
Thomas+2018, Epeak: 17.8eV
Figure 1. The incident radiation used to compute different NLR
emission models. The models of F16 are used in this work, and
allow a range of slope (αpl) of the ionizing portion of the spec-
trum. The F16 models use a simplified parametric form which
describes the thermal emission from the accretion disc surround-
ing the central black hole. We use a model grid that includes a
range of spectral indices, αpl, in the rest-frame UV to soft X-ray
[log(ν/Hz) & 15], as indicated. The shape is described in Schart-
mann et al. (2005) and Feltre et al. (2012). We also plot the inci-
dent radiation used in a number of different works for comparison,
as specified in the legend. We compare the emission line ratios of
these different models to observations in Section 4.
2.2 Emission from AGN narrow-line regions
We appeal to the Feltre et al. (2016, hereafter F16) models of
AGN-NLR emission to introduce simultaneous fitting of the
physical properties of NLR and H ii regions with beagle.
The F16 NLR models were produced with cloudy c13.03
(Ferland et al. 2013, the same version as used by Gutkin
et al. 2016 to produce the H ii-region models). Three impor-
tant aspects must be considered to incorporate the treatment
of NLR emission into beagle: the shape of the AGN ioniz-
ing continuum; the properties of the NLR gas itself; and the
integration of the F16 model grid with the H ii-region model
grid within beagle. In the next paragraphs, we summarise
the salient features of our approach to achieve these goals.
In the F16 model, the radiation illuminating the NLR gas is
the thermal emission from the accretion disc surrounding the
central black hole. For simplicity, this emission is described by
a broken power law, as shown in Fig. 1 (see also equation 5 of
F16). The strengths of emission lines emerging from the NLR,
as well as the ratios between them, are primarily sensitive to
the slope, αpl, of the power law at high frequencies (short
wavelengths). F16 allow this parameter to vary in the range
−1.2 < αpl < −2.0, which are generally used as “standard”
values for the power-law slope (e.g. Groves et al. 2004).
This simplified description of the purely thermal emission
from the accretion disc does not attempt to model the hard
X-ray component thought to arise from the Comptonization
of photons scattered in the hot corona surrounding the disc
(Haardt & Maraschi 1991). A soft X-ray component is also
often observed in AGNs, the exact nature of which is still de-
bated (e.g., disc reflection, Crummy et al. 2006; or optically-
thick Comptonized disc emission, which can explain the soft
MNRAS 000, 1–24 (2018)
4 Vidal-García et al.
X-ray excess in narrow-line Seyfert-Is, Done et al. 2012). In
practice, the flexibility in αpl allows the exploration of dif-
ferent contributions of soft X-rays to the incident ionizing
radiation.
Once the spectrum of the accretion-disc radiation reaching
the NLR is defined, cloudy can be used to model the trans-
fer of this radiation through the NLR gas. The F16 NLR grid
employs the option of ‘open geometry’, suitable for gas with
low covering fraction. The models are computed for a typical
accretion luminosity (i.e., the integral of the broken power law
in Fig. 1) of Lacc = 1045 erg s−1 cm−2 (denoted LAGN in F16).
The equation 1 of this reference also shows that the number
of ionizing photons, Q(H), is proportional to Lacc, for a given
αpl. 2 The version of the F16 model incorporated into beagle
includes further updates by Mignoli et al. (2019) which better
account for the observed [Nv]λ1240/[He ii]λ1640 emission-
line ratios in a sample of high-redshift type-2 AGNs. Specifi-
cally, the inner radius of the gas is set to rin ≈ 90 pc, and the
internal micro-turbulence velocity to vmicr = 100 km s−1.
To incorporate nebular emission from AGN NLR into bea-
gle, we use a grid of F16 models covering full ranges in
αpl and other NLR-gas parameters: the metallicity, noted
Znlrgas ; the ionization parameter, logUnlrs ; and the dust-to-
metal mass ratio, ξnlrd . In this work, we fix the hydrogen
density in the NLR to nnlrh = 1000 cm−3, as typically mea-
sured from optical line doublets (Osterbrock & Ferland 2006).
Furthermore, as for the H ii regions (Section 2.1), we fix the
C/O abundance ratio in the NLR gas to solar and use equa-
tion (1) to describe the dependence of N/H on O/H. Table 1
lists the full set of physical parameters pertaining to the H ii
and NLR regions in our models.
In combining the H ii-region and NLR model grids in bea-
gle, we allow for the normalisation of the NLR models to
change. In practice, we achieve this by changing the param-
eter Lacc, which controls the absolute luminosities of NLR
emission lines (Table 1). We note that this is compatible with
the adoption of a fixed Lacc = 1045 erg s−1 cm−2 in the pho-
toionization calculations of F16. Indeed, we find that, at fixed
other parameters, the emission-line luminosities computed
with cloudy for values of Lacc in the range from 1042 to
1048 erg s−1 cm−2 differ by only 0.1 to 6 per cent (depending
on the considered line) from those obtained by simply scal-
ing the luminosities of a model with Lacc = 1045 erg s−1 cm−2.
This is because we assume a fixed scaling relation between
the inner radius of the ionized gas (rin) and Lacc, given by
Lacc/(4pir
2
in) ' 102 erg s−1 cm−2 (Feltre et al. 2016). 3 We
must also account for the inhomogeneous distribution of the
narrow-line emitting gas around the central AGN. This can
be achieved by multiplying Lacc by a gas covering fraction,
whose recovery from observables will however be degenerate
with Lacc. To avoid introducing this extra degeneracy when
fitting to data with beagle, we fix the covering fraction of
NLR gas to 10 per cent, which is within the range of values
2 For completeness, log Q(H) = 54.7787, 54.7429, 54.6364, 54.4972
for αpl = −1.2,−1.4,−1.7,−2.0, respectively and Lacc =
1045 erg s−1 cm−2.
3 This negligible influence of Lacc on model predictions at fixed
other parameters (including logUnlrs ) arises from a degeneracy be-
tween Lacc and the volume-filling factor of the gas entering the
definition of the ionization parameter (see equation 4 of F16).
(2-20%) obtained for the covering fraction of the NLR gas of
the Palomar-Green quasar sample (Baskin & Laor 2005).
Since we are modelling type-2 AGNs, we are focused on
modelling the narrow line-emitting region. Continuum emis-
sion from the accretion disc itself is assumed to be completely
obscured by the dust surrounding the AGN centre, often re-
ferred to as the AGN torus. We use an open geometry in
cloudy, which means that any non-ionizing continuum light
(of the form in Fig. 1) that is reflected from the NLR is in-
cluded in the models (see Sections 2.2.4 and 2.3 of Hazy c17
manual), however, its contribution to most SEDs is negligible
relative to the stellar continuum in practice. For simplicity, we
assume that the radiation by the accretion disc and surround-
ing NLR gas are spherically symmetric. This is a reasonable
assumption when teamed with the covering fraction, but does
not take account of scattering or absorption of the accretion
disc radiation by the dust in the torus, nor of the anisotropic
nature of the emission from the accretion disc itself.
2.3 Attenuation by dust
In this work, we adopt the Charlot & Fall (2000, hereafter
CF00) two-component dust model. This model accounts for
the enhanced dust content in stellar birth clouds compared
to the diffuse ISM. The dust within birth clouds is split be-
tween the inner H ii and outer H i regions, with the dust
in H ii regions being self-consistently accounted for in the
nebular emission grid of Gutkin et al. (2016). As of beagle
v0.27.1, this dust component is accounted for self-consistently
within the implementation of the CF00 model, as described
in Curtis-Lake et al. (2021).
When combining NLR emission with the emission from
stars we account for the dust both within the NLR and the
diffuse ISM. As for the H ii region models, dust is included in
the F16 cloudy models of the NLR. Hence, within beagle,
in the framework of the CF00 dust model, light emitted from
the NLR is subject only to further attenuation by dust in the
diffuse ISM.
3 RETRIEVAL OF PARAMETERS FROM
IDEALISED MODELS
3.1 Idealised models of active galaxies
To investigate how well we can identify type-2 AGNs and con-
strain their NLR properties at different redshifts, we choose
to model two ‘typical’ star-forming galaxies, one at z = 0 and
another one at z = 2, to which we add various levels of AGN
contribution. We select the input parameters of these models
by appealing to the IllustrisTNG-100 cosmological, magneto-
hydrodynamic simulation of galaxy formation (Springel et al.
2018; Pillepich et al. 2018; Naiman et al. 2018; Nelson et al.
2018; Marinacci et al. 2018). See Hirschmann et al. in prep.
for the details of how these parameters have been retrieved
(see also Hirschmann et al. 2017, 2019). Fig. 2 shows his-
tograms of selected physical properties of galaxies with stellar
masses at least 3× 109 M4 extracted from this simulation.
4 The stellar masses that we study are significantly higher than
the mass-particle resolution of the IllustrisTNG-100.
MNRAS 000, 1–24 (2018)
beagle-agn 5
Table 1. List of physical parameters used in this work to simulate and fit to (active and inactive) galaxy spectra.
Parameter Description
Mtot/M Integrated SFH.
M/M Stellar mass, including stellar remnants.
ψ/Myr−1 Current star formation rate.
Zhiigas/Z, Znlrgas/Z Metallicity of gas in H ii regions and the NLR, respectively.
τˆv Total V -band attenuation optical depth in the ISM.
µ Fraction of τˆv arising from dust in the diffuse ISM, fixed to 0.4.
logUs, logUnlrs Effective gas ionization parameter in H ii regions and the NLR, respectively.
ξd, ξnlrd Dust-to-metal mass ratio in H ii regions and the NLR, respectively.
nh/cm−3, nnlrh /cm−3 Hydrogen gas density in H ii regions (fixed to 100) and the NLR (fixed to 1000), respectively.
(C/O)/(C/O) Carbon-to-oxygen abundance ratio in units of (C/O) = 0.44, fixed to unity.
mup/M Upper mass cutoff of the IMF, fixed to 100.
t/Gyr Age of the oldest stars.
αpl Slope of the ionizing radiation from thermal emission of the accretion disc around the central black hole.
Lacc Integrated thermal emission from the accretion disc.
−2 0 2
10−3
10−2
10−1
100
N
/N
m
a
x
−1 0 1 −4 −3 −2 −1 −1 0 1 −5 −4 −3 −2 −1
−2 0 2
log(ψ / M¯ yr−1)
10−3
10−2
10−1
100
N
/N
m
a
x
−1 0 1
log(Zhiigas/Z¯)
−4 −3 −2 −1
logUs
−1 0 1
log(Znlrgas /Z¯)
−5 −4 −3 −2 −1
logUnlrs
z ∼ 0
z ∼ 2
Figure 2. Histograms of the physical parameter space covered by galaxies in IllustrisTNG-100 at z=0 (top row) and z=2 (bottom row).
Dashed vertical lines indicate the mean value of each quantity, while the solid black lines indicate the values used.
Table 2. Properties of the idealised models of active galaxies at
z = 0 and z = 2 used to test parameter retrieval with beagle.
Parameter z = 0 z = 2
log(M/M) 10.5 10.5
ψ/Myr−1 0.4 1.3
Zhiigas 0.014 0.007
logUs −3.6 −3.0
nh/cm−3 100 100
ξd 0.3 0.3
t/Gyr 10 1
Znlrgas 0.030 0.012
logUnlrs −3.5 −2.75
nnlrh /cm−3 1000 1000
ξnlrd 0.3 0.3
τˆv 0.3 0.3
µ 0.4 0.4
αpl −1.7 −1.7
We fix the stellar masses of both z = 0 and z = 2 galax-
ies to log(M/M) = 10.5, corresponding to just below the
turnover of the stellar mass function of star-forming galaxies
measured over this redshift range (e.g. Tomczak et al. 2014).
Setting the ages of the oldest stars to 10 Gyr at z = 0 and
1 Gyr at z = 2 ensures that the objects are younger than the
Universe. We assume a constant star formation history to
which we add a 10-Myr old burst of star formation. We take
the current star formation rates of the z = 0 and 2 galaxies to
correspond to the mean values from the Illustris simulation
at these redshifts, i.e., ψ ≈ 0.4 and 1.3 M yr−1, respectively.
The distributions of Zhiigas and logUs for star-forming galaxies
at z ∼ 0 and 2 in the Illustris simulation are fairly broad and
allow for some freedom in choosing typical values. We adopt
values of Zhiigas and logUs such that the purely star-forming
models lie on the sequence of SDSS star-forming galaxies in
the classical [N ii]λ6584/Hα BPT diagram (see Fig. 3). For
the AGN component, we adopt the mean Znlrgas from Illustris
at each redshift and tweak logUnlrs to ensure that the models
sample the AGN portion of the [O i]λ6300/Hα BPT (rather
than the LINER/shock region at lower [O iii]λ5007/Hβ and
higher [O i]λ6300/Hα).
Among the other parameters remaining to be specified, the
C/O abundance ratio is fixed to solar, the hydrogen densities
in the H ii-region and NLR gas to log(nh/cm−3) = 2 and 3,
respectively, and the upper mass cutoff of the IMF to 100M
MNRAS 000, 1–24 (2018)
6 Vidal-García et al.
Figure 3. The [N ii]λ6584, [S ii]λ6717, λ6731 and [O i]λ6300 BPT diagrams showing the coverage of galaxies within SDSS as grey filled
contours. Violet symbols show pure NLR models chosen to represent typical AGNs at z = 0 (circles) and z = 2 (stars), while orange
symbols show the models chosen to represent typical star-forming galaxies at the same redshifts. The parameters chosen to represent
these galaxies are motivated by the IllustrisTNG-100 simulation (see text for details). The points spanning between the two extremes
show different contributions of NLR model to the corresponding star-forming model, parametrized by the fractional contribution of NLR
to the Hβ flux (as described in the legend).
(see Sections 2.1 and 2.2). We further fix the dust-to-metal
mass ratio to ξd = ξnlrd = 0.3, the V-band attenuation optical
depth and its fraction arising from the dust in the ambient
ISM to τˆv = 0.3 and µ = 0.4 respectively. Finally, we set
the power-law slope of the accretion-disc emission model to
αpl = −1.7, in the middle of the range probed by the AGN-
NLR model (Section 2.2). Table 2 summarises the parameters
of these idealised models.
When combining the H ii-region and NLR models, we con-
sider three types of galaxies at each redshift: one for which
line emission is dominated by star formation, one for which it
is dominated by the NLR, and one with equal contribution.
We set the ratio of Hβ luminosities from H ii regions and the
NLR be HβHII/HβNLR = 0.1, 1 and 10. This ensures that
our models bridge the space between the regions dominated
by star-forming galaxies and AGNs in the BPT. We ensure
the highest value of Lacc (obtained for HβHII/HβNLR = 0.1)
is lower than the Eddington luminosity,
LEdd =
4piGcmp
σe
Mbh = 1.26× 1038(Mbh/M) erg s−1 , (2)
whereMbh is the black-hole mass, c the speed of light,mp the
proton mass and σe the Thomson-scattering cross-section for
electrons. In the local Universe, the mass of the central black
hole in active galaxies tends to scale with total stellar mass
as Mbh ∼ 0.025M? (e.g. Reines & Volonteri 2015). Hence,
the condition Lacc < LEdd requires the galaxy stellar mass to
be at least ∼ 40 times the corresponding black-hole mass in
equation (2), which is the case in our models at both redshifts.
From these simulated spectra, we make a mock catalogue
using beagle, which provides simulated line fluxes. We add
noise to these line fluxes, incorporating a flat noise model
with the noise level for each line set according to the desired
S/N in Hβ. We fit to noisy line fluxes, rather than to the
full spectrum, since we are primarily interested in the star-
forming and narrow-line region contributions to the emission
lines themselves.
Although we only fit the noisy model line fluxes, we show
in Fig. 4 the window of the full optical spectrum of the z =
0 galaxy where the transitions used in the BPT diagrams
arise. In particular, we show spectra with two different noise
levels, S/N(Hβ)=3 and 40, for the three galaxies modeled as
a function of the contribution of the different regions to the
line emission. In all the cases, all the emission lines are visible
for the S/N(Hβ)=40 spectra, but only the strongest lines are
still detectable for the noisier ones with S/N(Hβ)=3.
At low S/N(Hβ), the perturbed line fluxes for weak lines
can fluctuate significantly. To mitigate this effect, we produce
10 spectra for each S/N(Hβ) value and average the results
after having fitted the line fluxes with beagle.
3.2 Fitting with beagle
Our aims when fitting with beagle are two-fold. First we
hope to distinguish objects with an obscured-AGN compo-
nent contributing to the emission-line fluxes. Additionally we
wish to constrain the gas properties of the NLR and star-
forming regions. Our base set of observables are the lines
used to make the three emission-line diagnostic diagrams
identified by Baldwin et al. (1981) for distinguishing the
dominant ionizing source, namely [O iii]λ5007, [O i]λ6300,
Hα, Hβ, [N ii]λ6584 and [S ii]λ6717, λ6731. We add the
[O ii]λ3726, λ3729 emission line doublet to gain constraints
on the ionization state of the gas (most directly together with
[O iii]λ5007 for constraints on logUs and logUnlrs ).
From a first analysis we found that fitting to indi-
vidual line fluxes produces biased estimates of AGN
parameters, and so we construct the likelihood using line
ratios particularly sensitive to the estimated physical
MNRAS 000, 1–24 (2018)
beagle-agn 7
(a) SF-dominated (b) AGN-dominated
(c) Equal contributions
Figure 4. Window of the optical spectrum of the z = 0 galaxy with different levels of noise, parametrized in terms of the observed Hβ
signal-to-noise ratio, S/N(Hβ), for the three types of galaxies considered.
parameters. In particular, we use: the BPT line ratios
[O i]λ6300/Hα, [N ii]λ6584/Hα, [S ii]λ6717, λ6731/Hα
and [O iii]λ5007/Hβ; Hβ/Hα to provide constraints
on dust attenuation; [O i]λ6300/[O ii]λ3726, λ3729 and
[O ii]λ3726, λ3729/[O iii]λ5007 to provide constraints on the
ionization parameters. To constrain the absolute normali-
sation of the line fluxes, we also include the comparison of
Hα flux in the likelihood calculations. This, in effect, means
that we are fitting directly to these emission line ratios (plus
Hα line flux), not to the individual line fluxes. The chosen
line ratios produce improved estimates of various parameters
thanks to the higher weighting of lines in the partially
ionized zone ([N ii]λ6584, [S ii]λ6717, λ6731, [O i]λ6300) to
the likelihood evaluation. When these lines are faint and not
fitted as line-ratios but rather individually, the information
in the likelihood is dominated by stronger lines which
contain less information about the NLR properties (Hα, Hβ,
[O ii]λ3726, λ3729). We make the simplifying assumption
that the line ratios are independent. The method of using
line ratios in this way is similar to that adopted by other
works, such as Pérez-Montero et al. (2019).
Within the likelihood determination, we assume that the
measurement of the ratio is Gaussian-distributed, with µ the
mean and σ the standard deviation of the distribution. The
ratio between two Gaussian-distributed variables, Z = X/Y
[where X ∼ N (µx, σ2x) and Y ∼ N (µy, σ2y)], is strictly de-
scribed by a Cauchy distribution when the mean of each dis-
tribution is zero. However, the distribution of the ratio of two
independent normal random variables can be approximated
as a Gaussian with variance:
σ2z =
µ2x
µ2y
(
σ2x
µ2x
+
σ2y
µ2y
)
(3)
when 0 < σx/µx < λ 6 1 and 0 < σy/µy 6
√
λ2 − σ2x/µ2x <
λ, where λ is a constant between zero and one (Díaz-Francés
& Rubio 2013), and σ/µ is equal to the inverse of the S/N on
the given parameter. We always choose ratios with the lower
S/N line as the numerator, which safely respects the condition
on σx/µx5. Re-writing the central part of the second condi-
tion using σy/µy = 1N σx/µx (which is equivalent to saying
the signal-to-noise on y is a multiple of N times higher than
the signal-to-noise on x), we have
√
(N2 + 1)σx/µx < Nλ,
which is trivially true if N > 1 and σx/µx < λ. Therefore, we
are generally within the regime where we can safely model
the distribution of the line ratios as a Gaussian. In practice,
the lines on the numerator do not always satisfy σx/µx < 1.
However, when they have very low S/N (high σx/µx), their
weighting within the likelihood is much lower and so this
should not significantly affect the results.
5 We note here that this condition should be preserved also work-
ing with real objects, where the heterogeneous line intensities from
object to object might imply switching the numerator and the de-
nominator, when needed.
MNRAS 000, 1–24 (2018)
8 Vidal-García et al.
Table 3 presents the prior limits and fixed values of dif-
ferent parameters. Note that the mass and age of the stel-
lar population are not displayed in the table because they
act as nuisance parameters of the gas properties, which are
the main focus of the following discussion. The table shows
the three different parameter configurations explored: vary-
ing αpl while fixing ξnlrd , varying ξ
nlr
d while fixing αpl, and
fixing both of these parameters.
3.3 High-quality spectra with S/N(Hβ)=100
3.3.1 AGN versus star-formation dominated spectra
We start by investigating the idealised scenario of high-
S/N spectra with full coverage of the lines of interest. From
our grid of models, we test the retrieval of physical pa-
rameters for the AGN and star-formation dominated cases
(Hβsf/Hβnlr = 0.1 and 0.9, respectively) with high S/N spec-
tra [S/N(Hβ)=100]. At this S/N in Hβ we reach S/N > 10
in all lines used in the analysis. We fit all components (H ii,
stellar continuum and NLR) to each spectrum to ensure that
the parameters of the dominant component are not biased
by the expanded parameter space including the rest of the
components.
Fig. 5 summarises the joint posterior probability distribu-
tions derived from beagle fits to the z = 0 and z = 2 AGN-
dominated spectra. We display only parameters that affect
the NLR emission in the model, even though H ii and stellar
parameters were fitted. The plot is an effective average of the
posteriors derived for each of the 10 realisations at each grid
point (where we take the same input spectrum and produce
10 realisations with random noise). For each Monte Carlo
realisation, we take samples from the posterior probability
(derived using beagle) for each fitted parameter before com-
puting the logarithm of the ratio to the corresponding input
value. We then fit a bi-variate Gaussian to the joint distribu-
tion of logarithmic ratios for each parameter pair. The mean
of the 10 bi-variate Gaussian centres are plotted as crosses.
The ovals show the average of the bi-variate Gaussian distri-
butions of the fits to the 10 noisy realisations (specifically,
the covariance matrices were constructed using the mean of
the corresponding entries in the individual bivariate Gaussian
covariance matrices). We show an illustration of the creation
of this type of average triangle plot in Fig. 6. We also show
an example triangle and marginal plot for one of the noisy
realisations to the z = 0 AGN dominated spectrum in Fig 7.
In Appendix A, we show the triangle plots of one of the 10
noisy realisations of the same spectrum for different param-
eter configurations in objects at different redshifts and with
different contributions of the star-forming and NL regions to
the Hβ line flux.
The results are displayed for three different parameter con-
figurations; varying αpl while fixing ξnlrd (yellow), varying ξ
nlr
d
while fixing αpl(blue), and fixing both of these parameters
(red). We find that varying ξnlrd leads to logU
nlr
s and Znlrgas
being under-estimated for both the z = 0 and z = 2 objects.
ξnlrd is clearly poorly constrained and there are degeneracies
between ξnlrd and logU
nlr
s . The depletion onto dust grains
can, in principle, be constrained using observables from ele-
ments with different refractory properties. For example, Oxy-
gen is depleted onto dust grains, but Nitrogen is not. A stan-
dard line ratio used for determining gas-phase N/O ratio is
[N ii]λ6584/[O ii]λ3726, λ3729, because both ions have simi-
lar ionization energies. We tried explicitly including this ratio
in the fits with little effect on the results. This is because, in
the presence of ionizing radiation from an accretion disc, the
[N ii]λ6584/[O ii]λ3726, λ3729 ratio is heavily affected by the
hardness of the incident radiation, and the residual depen-
dence on gas-phase N/O abundance is insufficient to con-
strain ξnlrd .
Allowing αpl to vary leads to over-estimated Lacc while
still under-estimating Znlrgas . There are no clear degeneracies
that would lead to such a biased Znlrgas estimate. There is a
clear degeneracy between αpl and Lacc, resulting from the
definition of Lacc, the integration of the thermal accretion
disc model (Fig. 1). Increasing αpl while keeping everything
else constant will increase Lacc. We find the best recovery of
Lacc, Znlrgas and logUnlrs when fixing both ξnlrd and αpl. Feltre
et al. (2016) show that increasing αpl pushes the models to
higher values of [O iii]λ5007/Hβ in the [N ii]λ6584/Hα BPT
diagram (see their figure 2), so it may be appropriate to vary
αpl for objects with more extreme measured line ratios. See
paper II for a further analysis of how the variation of αpl
affects line ratios in the models and the fitting to real data.
The fraction of the V-band attenuation optical depth aris-
ing from the dust in the ambient ISM, µ, is not well con-
strained by this set of observables for the z = 0 galaxy nor
for the z = 2 one.
Fig. 8 displays the stellar and H ii region parameters for
the star-formation dominated z = 0 and z = 2 galaxies. We
show results when nebular ξd is fixed and fitted. We find that
ξd is poorly constrained, leading to over-estimated logUs and
Zhiigas, as well as tighter constraints on the biased estimates.
The results are somewhat less biased for the z = 2 galaxy.
Finally, we note that for the star-forming dominated case,
both in the z = 0 and the z = 2 cases, µ is well constrained.
3.3.2 Equal contributions by AGN and star formation
We proceed with the most challenging case of equal contribu-
tions of NLR and H ii regions to the Hβ flux. Figs. 9 and 10
show the average joint posteriors for all fitted parameters for
the z = 0 and z = 2 galaxy, respectively. We display the
constraints for different configurations of fitted parameters,
as indicated in the legend. Notably, the constraints are quite
different for the two different objects. For instance, αpl is
well recovered for the z = 0 object, but is significantly bi-
ased for the z = 2 object. The biased estimates of αpl for the
z = 2 object consequently lead to significantly biased esti-
mates of logUnlrs , Znlrgas and Lacc. We therefore learn that the
accuracy of the αpl estimates depends heavily on the region
of the observable parameter space that is being probed. In
other words, αpl can be well retrieved for objects in certain
regions of the [N ii]λ6584 BPT diagram, probably dependent
on the sampling of the parameter space6 by the models and
the degeneracies present in that particular region. We there-
fore suggest that αpl should be fixed as standard when using
this set of observables, unless reaching an acceptable fit re-
quires this parameter to be varied.
A parameter which was found to be problematic in the
6 An irregular sampling of the parameter space is a natural out-
come from physical motivated model grids.
MNRAS 000, 1–24 (2018)
beagle-agn 9
Table 3. Uniform prior limits and fixed values applied to the three different fits for different parameters in beagle.
Parameter Fit αpl Fit ξnlrd Fixed αpland ξ
nlr
d
Ψ = log(ψ/Myr−1) Uniform ∈ [−4., 4.] Uniform ∈ [−4., 4.] Uniform ∈ [−4., 4.]
log(Zhiigas/Z) Uniform ∈ [−2.2, 0.24] Uniform ∈ [−2.2, 0.24] Uniform ∈ [−2.2, 0.24]
logUs Uniform ∈ [−4,−1] Uniform ∈ [−4,−1] Uniform ∈ [−4,−1]
ξd Fixed to 0.3 Fixed to 0.3 Fixed to 0.3
τˆv Uniform ∈ [0, 2] Uniform ∈ [0, 2] Uniform ∈ [0, 2]
µ Uniform ∈ [0, 1] Uniform ∈ [0, 1] Uniform ∈ [0, 1]
log(Lacc L) Uniform ∈ [40, 48] Uniform ∈ [40, 48] Uniform ∈ [40, 48]
ξnlrd Fixed to 0.3 Uniform ∈ [0.1, 0.5] Fixed to 0.3
logUnlrs Uniform ∈ [−4,−1.5] Uniform ∈ [−4,−1.5] Uniform ∈ [−4,−1.5]
log(Znlrgas/Z) Uniform ∈ [−2.2, 0.24] Uniform ∈ [−2.2, 0.24] Uniform ∈ [−2.2, 0.24]
αpl Uniform ∈ [−2.0,−1.2] Fixed to -1.7 Fixed to -1.7
Lacc−0
.05
0.0
0
0.0
5
Z
N
L
R
ga
s
Lacc−0
.1
0.0
U
N
L
R
S
ZNLRgas
Lacc−0
.1
0.0
ξ
N
L
R
d
ZNLRgas U
NLR
S
Lacc
−1
0
1
τˆ V
ZNLRgas U
NLR
S ξ
NLR
d
Lacc
−2
0
µ
ZNLRgas U
NLR
S ξ
NLR
d τˆV
−0
.1 0.0 0.1
Lacc
0.0
00
0.0
25
0.0
50
α
P
L
−0
.050.0
0
0.0
5
ZNLRgas
−0
.1 0.0
UNLRS
−0
.1 0.0
ξNLRd
−1 0 1
τˆV
−2
.5 0.0
µ
axes: log(x/xin)
x =
x
=
ξNLRd fitted, αPL fixed
ξNLRd fixed, αPL fitted
ξNLRd fixed, αPL fixed
(a) z = 0
Lacc−0
.05
0.0
0
0.0
5
Z
N
L
R
ga
s
Lacc−0
.1
0.0
U
N
L
R
S
ZNLRgas
Lacc−0
.1
0.0
ξ
N
L
R
d
ZNLRgas U
NLR
S
Lacc
−1
0
1
τˆ V
ZNLRgas U
NLR
S ξ
NLR
d
Lacc
−2
0
µ
ZNLRgas U
NLR
S ξ
NLR
d τˆV
−0
.1 0.0 0.1
Lacc
0.0
00
0.0
25
0.0
50
α
P
L
−0
.050.0
0
0.0
5
ZNLRgas
−0
.1 0.0
UNLRS
−0
.1 0.0
ξNLRd
−1 0 1
τˆV
−2
.5 0.0
µ
axes: log(x/xin)
x =
x
=
ξNLRd fitted, αPL fixed
ξNLRd fixed, αPL fitted
ξNLRd fixed, αPL fixed
(b) z = 2
Figure 5. Triangle plots of the parameter retrieval for the z = 0 (left) and z = 2 (right) NLR-dominated galaxies with S/N∼ 100 in
Hβ. Only NLR parameters are displayed even though SF parameters are also included in the fits. The triangle plots display the average
parameter constraints relative to the input values from fits to 10 realisations of noisy spectrum. The crosses show the average bias for each
parameter pair, while the ovals show the 1σ contour of the bi-variate Gaussian fitted to the joint posteriors (see text for details about how
this plot is made). These show an approximate representation of the biases and degeneracies between different parameters. We display
the results for different combinations of fixed or fitted ξnlrd and αpl, as specified in the legend.
AGN-dominated fits is ξnlrd (Section 3.3.1), which again is
poorly constrained in case of equal contributions by NLR and
H ii regions. ξnlrd is degenerate with Ψ, Lacc and logU
nlr
s ,
meaning that when ξnlrd is un-constrained, the recovery of
these parameters will also be biased. Fixing ξnlrd somewhat
improves the retrieval of the NLR parameters, but for the
z = 0 galaxy, ξd is still poorly constrained. This parameter
is degenerate with logUs which is, in turn, degenerate with
logUnlrs . Therefore fixing both ξnlrd and ξd provide less bi-
ased estimates of the ionization parameters of the NLR and
H ii region emission. We therefore obtain best constraints for
both objects when ξd and ξnlrd are both fixed.
Finally, µ is well constrained for the z = 0 object for the
different configurations of the fitted and fixed parameters but
it is not so well constrained for the same configurations for
the z = 2 object. In order not to be biased by the incertitude
on the retrieval of this parameter, we fix the value of µ in the
following sections.
3.4 Influence of S/N on parameter retrieval
We investigate the S/N in emission lines required to derive
un-biased parameter estimates. Based on the results for the
S/N(Hβ) ∼ 100 fits, we fix ξd, ξnlrd , αpl and µ. For the
AGN-dominated galaxies, the z = 0 object constraints de-
grade with S/N more noticeably than for the z = 2 object,
so to be conservative, we only show the z = 0 results in
Fig. 11. These results suggest that the AGN parameters are
well-constrained at S/N(Hβ) & 20. With decreasing S/N, de-
generacies between Lacc and τˆv, as well as Lacc and logUnlrs
MNRAS 000, 1–24 (2018)
10 Vidal-García et al.
−0.05 0.00
−0.09
−0.08
−0.07
−0.06
−0.05
−0.04
−0.03
lo
g(
U
N
L
R
S
/i
np
ut
U
N
L
R
S
)
−0.05 0.00 −0.05 0.00 −0.05 0.00 −0.05 0.00 −0.05 0.00 −0.05 0.00 −0.05 0.00 −0.05 0.00 −0.05 0.00 −0.05 0.00
log(ZNLRgas /input ZNLRgas )
Figure 6. Demonstration of the creation of the average triangle plots, showing the 2D posterior probability distribution between logUnlrs
and Znlrgas for each Monte Carlo realisation (first 10 panels) for the NLR-dominated z ∼ 2 galaxy. The beagle fits correspond to those with
ξnlrd fitted and αpl fixed (orange ellipses in Fig.5). In each of the first ten panels, the blue points show a sampling from the 2D posterior
resulting from the beagle fits, while the orange contours show the best-fit Gaussian fit to the 2σ credible region. The final (eleventh)
panel shows each of the individual Gaussian fits to the credible intervals in orange (repeated from each of the 10 previous panels) as
well as the average in black. The covariance matrix of the black ellipse is calculated from the mean of the corresponding values in the 10
individual orange ellipses.
become more problematic, causing increasingly biased esti-
mates of the NLR properties. The specific behaviour with
low S/N will be different for each object as our results in
Section 3 show how dependent parameter retrieval is on the
point in parameter space occupied by the object, because of
the very non-linear and irregular coverage of the models in
observable space.
The H ii region properties for the star-formation domi-
nated objects are well constrained and un-biased down to
S/N(Hβ) ∼ 10, as displayed in Fig. 12. Disentangling the con-
tributions from star-forming and NLR components requires
higher S/N, however, as illustrated by Fig. 13, which shows
the parameter retrieval for the case with equal contributions
by SF and NLR to Hβ. Parameter retrieval starts to become
biased and poorly constrained for S/N(Hβ) . 30.
3.5 NLR contribution to the spectrum
beagle can be used to identify the potential contribution
by an AGN to the emission-line luminosities of a galaxy. To
test how well beagle can detect the presence of an AGN, we
show the retrieved fractional NLR contribution to Hβ flux
as a function of the S/N on Hβ in Fig. 14. If we wish only
to identify objects with a significant contribution from the
NLR to Hβ, a S/N of 3 in Hβ is sufficient for our adopted
properties of a typical z = 0 galaxy. However, this S/N is
insufficient to reliably estimate the fractional contribution to
Hβ, which is over-estimated if stars and the AGN contribute
in equal amounts to the line flux.
For our z = 2 galaxy, the fractional contribution of the
NLR to Hβ is completely unconstrained at S/N(Hβ) = 3.
For S/N(Hβ) ∼ 10, the NLR contribution is over-estimated
for the mixed case. However, only a small probability of a
high NLR contribution to the star-formation dominated spec-
trum (seen from the low fractional area of the yellow curve
above 0.5) means that even at this reasonably low S/N, a
fairly clean sample of objects with significant contribution
of NLR to Hβ flux can selected. Being more conservative,
S/N∼ 20− 30 is required to firmly identify objects with pos-
sible NLR contribution, as well as constraining that contribu-
tion. It is harder to distinguish NLR contribution in the z = 2
case because the metallicity is lower, and so key line fluxes
such as [O iii]λ5007, [O ii]λ3726, λ3729 or [O i]λ6300 are also
much fainter with respect to Hβ. This leads to much more
overlap between the star-forming and NLR models in this
region of the BPT diagram (see Fig. 3), making separation
of the NLR and star-forming contributions more challenging.
The addition of other line ratios, especially those of transi-
tion in the UV of highly ionized species, such as C2+, C3+,
N3+, N4+, Ne3+ and Si2+, may improve the ability of finding
NLR contribution at low metallicity (Feltre et al. 2016).
4 DISCUSSION
Other works have investigated the derivation of NLR physical
parameters from emission lines. Throughout the discussion,
we compare the modelling assumptions we make to those of
other works, with the aim of outlining when and why dif-
ferent models will provide different (or similar) properties for
the same data sets. The works we compare to are: the strong-
line calibrations of Storchi-Bergmann et al. (1998) and Dors
(2021), and the Bayesian approaches of Pérez-Montero et al.
(2019) and Thomas et al. (2018b). Before making compar-
isons, we summarise the different approaches in more detail
here.
4.1 Different approaches to deriving NLR physical
properties
4.1.1 Storchi-Bergmann et al. 1998
The first calibration of oxygen abundances based on strong
optical line emission from NLRs was provided by Storchi-
Bergmann et al. (1998). They used cloudy to produce a
grid of emission line fluxes for an incident radiation field de-
scribed by the segmented power-law of Mathews & Ferland
(1987) (which is shown in comparison to that used in the
F16 models in Fig. 1). The grid covers a range of oxygen
abundances and ionization parameters. Although not explic-
itly stated, the ionization parameter is likely that at the inner
MNRAS 000, 1–24 (2018)
beagle-agn 11
Figure 7. An example triangle plot for one of the noisy realisations used to produce the average triangle plot in Fig. 5(a). The diagonal
panels show the 1D marginalised posterior probability distribution for each individual parameter. The grey shaded regions show the
marginalised 68% credible region. The off-diagonal panels show 2D posterior probability distributions with the dark, medium and light
blue contours showing the 68%, 95% and 99% credible intervals, respectively. The upper right inset shows the marginal plot for the
observables, displayed as the ratio between the fitted quantity (whether line ratio or line flux) divided by the observed value (y-axis) for
each line ratio or line flux used in the fit (x-axis labels). The blue shaded regions show the coverage of the models, while the coral points
and error bars show the observed values and associated uncertainties.
MNRAS 000, 1–24 (2018)
12 Vidal-García et al.
ψ
0.0
0.2
U
S
ψ
0.0
0.2
ξ d
US
ψ−
0.1
0.0
0.1
Z
H
II
ga
s
US ξd
ψ−1
0
τˆ V
US ξd ZHIIgas
−0
.25 0.0
0
0.2
5
ψ
−1
0µ
0.0
0
0.2
5
US
0.0
0
0.2
5
ξd
−0
.1 0.0 0.1
ZHIIgas
−1 0
τˆV
axes: log(x/xin)
x =
x
=
ξd fitted
ξd fixed
(a) z = 0
ψ
0.0
0.2
U
S
ψ
0.0
0.2
ξ d
US
ψ−
0.1
0.0
0.1
Z
H
II
ga
s
US ξd
ψ−1
0
τˆ V
US ξd ZHIIgas
−0
.25 0.0
0
0.2
5
ψ
−1
0µ
0.0
0
0.2
5
US
0.0
0
0.2
5
ξd
−0
.1 0.0 0.1
ZHIIgas
−1 0
τˆV
axes: log(x/xin)
x =
x
=
ξd fitted
ξd fixed
(b) z = 2
Figure 8. As in fig. 5, but showing triangle plots of the parameter retrieval for the z = 0 (left) and z = 2 (right) star-formation dominated
galaxies. We display only SF parameters, though NLR parameters were also included in the fits. Results when ξd is both fixed and fitted
are displayed according to the legend.
edge of the cloud (rather than at the radius of the Strömgrem
sphere as for the F16 model grid). They account for both pri-
mary and secondary nitrogen, i.e. dependent on the oxygen
abundance, following log(N/O) = 0.96[12 + log(O/H)]−9.29.
They also account for depletion of refractory elements from
the gas phase onto dust grains following the values of the
observed abundance of the interstellar medium from Cowie
& Songaila (1986). The following equations describe the two
calibrations for the gas-phase oxygen abundance derived by
Storchi-Bergmann et al. (1998):
(O/H)SB98,1=8.34 + (0.212x)− (0.012x2)− (0.002 y)
+(0.007xy)− (0.002x2y) + (6.52× 10−4 y2)
+(2.27× 10−4 xy2) + (8.87× 10−5 x2y2),
(4)
where x = [N ii]λ6548, λ6584/Hα and y =
[O iii]λ4959, λ5007/Hβ and
(O/H)SB98,2 = 8.643− (0.275u) + (0.164u2)
+ (0.655 v)− (0.154uv)− (0.021u2v)
+ (0.288v2) + (0.162uv2) + (0.0353u2v2),
(5)
where u = log([O ii]λ3726, λ3729/[O iii]λ4959, λ5007) and v
= log([N ii]λ6548, λ6584/Hα). All these line ratios make use
of de-reddened emission-line fluxes in their work. The term
(O/H) above corresponds to 12+log(O/H). Both calibrations
are valid for 8.4 6 12 + log(O/H) 6 9.4.
The calibrations are appropriate for nh = 300 cm−3, but
they also provide a correction factor for different densities:
(O/H)final = (O/H)− 0.1 log(nh/300) (6)
4.1.2 Dors 2021
The first strong-line calibration based on gas-phase oxygen
abundance estimates derived from the Te method (or elec-
tron temperature method) for NLRs was presented in Dors
(2021). From the Te-derived gas-phase oxygen abundances of
a sample of AGNs7, Dors (2021) provide a 2D surface fitted
to the objects within a 3D space defined by 12 + log(O/H),
R23 and P (defined below):
R23 =
[O ii]λ3726, λ3729 + [O iii]λ4959, λ5007
Hβ
(7)
P =
[O iii]λ4959, λ5007/Hβ
R23
(8)
The R23 ratio is often used as an oxygen abundance indicator
but is known to be dependent on ionization parameter, as
well as the hardness of the ionizing radiation (Pilyugin 2000),
while P is sensitive to radiation hardness (Pilyugin 2001).
The line fluxes measurements are also reddening corrected in
this case. The resulting calibration is defined as:
Znlrgas = (−1.00± 0.09)P + (0.036± 0.003)R23 + (8.80± 0.06)
(9)
where Znlrgas = 12 + log(O/H).
The oxygen abundances used to derive the above calibra-
tion do still rely on photoionization models, though indirectly.
In particular, when only one auroral line is detected (in the
case of Dors 2021, the [Oiii]λ4363 line), the temperature for
7 Selected from the location of the galaxies in the standard BPT
diagrams, see Dors et al. (2020a) for the details of the target se-
lection.
MNRAS 000, 1–24 (2018)
beagle-agn 13
ψ
−0
.5
0.0
0.5
Z
H
II
ga
s
ψ
0.0
0.5U
S
ZHIIgas
ψ
0.0
0.5
ξ d
ZHIIgas US
ψ−0
.5
0.0
0.5
L a
cc
ZHIIgas US ξd
ψ−0
.1
0.0Z
N
L
R
ga
s
ZHIIgas US ξd Lacc
ψ−
0.2
5
0.0
0
U
N
L
R
S
ZHIIgas US ξd Lacc ZNLRgas
ψ−0
.25
0.0
0
ξ
N
L
R
d
ZHIIgas US ξd Lacc ZNLRgas U
NLR
S
ψ
−1
0
τˆ V
ZHIIgas US ξd Lacc ZNLRgas U
NLR
S ξ
NLR
d
ψ−2
0µ
ZHIIgas US ξd Lacc ZNLRgas U
NLR
S ξ
NLR
d τˆV
−0
.5 0.0 0.5
ψ
−0
.5
0.0
0.5
α
P
L
−1
.0
−0
.50.0 0.5
ZHIIgas
0.0 0.5
US
0.0 0.5
ξd
−0
.5 0.0 0.5
Lacc
−0
.10
−0
.050.0
0
0.0
5
ZNLRgas
−0
.25 0.0
0
UNLRS
−0
.25 0.0
0
ξNLRd
−1 0
τˆV
−1 0 1
µ
axes: log(x/xin)
x =
x
=
ξNLRd fitted, αPL fixed, ξd fitted
ξNLRd fixed, αPL fitted, ξd fitted
ξNLRd fixed, αPL fixed, ξd fitted
ξNLRd fixed, αPL fixed, ξd fixed
Figure 9. As in fig. 5 but now showing the triangle plots of the parameter retrieval for the z = 0 galaxy with equal contribution from SF
and NLR to Hβ flux. We display the results for all fitted parameters, with a range of different configurations for fixing or fitting to ξnlrd ,
αpl.
part of the nebula is measured, and abundance estimates of
the ionization species found within that region is well de-
termined. Yet, to derive the total oxygen abundance, one
needs to account for all ionization species that exist in re-
gions with different effective temperatures. The temperature
of these other regions must, therefore, be inferred using other
methods. Dors et al. (2020b) created a grid of photoionization
models to derive the relationship between the temperature of
the high ionization region, t3 (where oxygen is doubly ion-
ized, traced by [O iii]λ4959, λ5007 and [Oiii]λ4363), and the
low ionization region, t2 (where oxygen is singly ionized). To
derive a relation between t2 and t3 suitable for the NLR, Dors
et al. (2020b) created a grid of models using ionizing spectra
with power-law slopes of αpl = −0.8,−1.1,−1.4, ionization
parameters over the range −3.5 to −0.5 (again, defined at
the inner radius of the ionized cloud), and metallicities over
MNRAS 000, 1–24 (2018)
14 Vidal-García et al.
ψ−2
.5
0.0
Z
H
II
ga
s
ψ−2
.5
0.0
2.5
U
S
ZHIIgas
ψ−
0.5
0.0
0.5
ξ d
ZHIIgas US
ψ
0.0
0
0.2
5
L a
cc
ZHIIgas US ξd
ψ
−0
.1
0.0
0.1
Z
N
L
R
ga
s
ZHIIgas US ξd Lacc
ψ−0
.5
0.0
U
N
L
R
S
ZHIIgas US ξd Lacc ZNLRgas
ψ−
0.5
0
−0
.25
ξ
N
L
R
d
ZHIIgas US ξd Lacc ZNLRgas U
NLR
S
ψ−2
0
2
τˆ V
ZHIIgas US ξd Lacc ZNLRgas U
NLR
S ξ
NLR
d
ψ−2
.5
0.0µ
ZHIIgas US ξd Lacc ZNLRgas U
NLR
S ξ
NLR
d τˆV
−2
.5 0.0
ψ
0.0
5
0.1
0
α
P
L
−2
.5 0.0
ZHIIgas
−2
.5 0.0 2.5
US
−0
.5 0.0 0.5
ξd
0.0
0
0.2
5
Lacc
−0
.1 0.0 0.1
ZNLRgas
−0
.25 0.0
0
UNLRS
−0
.50
−0
.25
ξNLRd
−1 0 1
τˆV
−2
.5 0.0
µ
axes: log(x/xin)
x =
x
=
ξNLRd fitted, αPL fixed, ξd fitted
ξNLRd fixed, αPL fitted, ξd fitted
ξNLRd fixed, αPL fixed, ξd fitted
ξNLRd fixed, αPL fixed, ξd fixed
Figure 10. As in fig. 5 but now showing the triangle plots of the parameter retrieval for the z = 2 galaxy with equal contribution from
SF and NLR to Hβ flux.
the range 0.2 < Znlrgas/Z < 2.0, and for a range of different
electron densities from 100 to 3000 cm−3. In this work, they
follow the relation log(N/O) = 1.29[12 + log(O/H)] − 11.84
to assign a given nitrogen abundance. They also fit a relation
between t2 and t3 for the whole grid, as well as for different
electron densities. This relation may therefore be reasonable
for ‘typical’ NLRs (if the grid presented in Dors et al. (2020b)
is appropriate for ‘typical’ NLRs), but individual galaxies will
have very different t2/t3 values, as demonstrated by the wide
spread in t2/t3 for their fiducial grid (their figure 5). Treat-
ment of depletion onto dust grains is not referred to in their
publications.
4.1.3 Pérez-Montero et al. 2019
Pérez-Montero et al. (2019) extended the Bayesian-like Hii-
MNRAS 000, 1–24 (2018)
beagle-agn 15
Lacc
−0
.2
0.0
Z
N
L
R
ga
s
Lacc
−0
.2
0.0
0.2
U
N
L
R
S
ZNLRgas
0.0 0.5
Lacc
−1
0
1
τˆ V
−0
.25 0.0
0
ZNLRgas
−0
.25 0.0
0
UNLRS
axes: log(x/xin)
x =
x
=
S/N(Hβ)=40
S/N(Hβ)=30
S/N(Hβ)=20
S/N(Hβ)=10
S/N(Hβ)=3
Figure 11. Triangle plots of the NLR parameter retrieval for the
z = 0 AGN-dominated objects. Different colors represent different
S/N in emission lines. For a given parameter, the axes show the
logarithmic fraction of the output by the input value.
ψ
−1
0
1
Z
H
II
ga
s
ψ
−1
0
1
U
S
ZHIIgas
−2
.5 0.0 2.5
ψ
−1
0
1
τˆ V
−1 0 1
ZHIIgas
−1 0 1
US
axes: log(x/xin)
x =
x
=
S/N(Hβ)=40
S/N(Hβ)=30
S/N(Hβ)=20
S/N(Hβ)=10
S/N(Hβ)=3
Figure 12. Same as Fig. 11 but for the parameter retrieval of the
H ii region for SF-dominated objects at z = 0.
ψ
−2
.5
0.0
Z
H
II
ga
s
ψ
0.0
2.5
U
S
ZHIIgas
ψ
0
1
L a
cc
ZHIIgas US
ψ
−0
.5
0.0
Z
N
L
R
ga
s
ZHIIgas US Lacc
ψ
−0
.5
0.0
U
N
L
R
S
ZHIIgas US Lacc ZNLRgas
−5 0
ψ
−1
0
1
τˆ V
−2
.5 0.0
ZHIIgas
0.0 2.5
US
0 1
Lacc
−0
.5 0.0
ZNLRgas
−0
.5 0.0
UNLRS
axes: log(x/xin)
x =
x
=
S/N(Hβ)=40
S/N(Hβ)=30
S/N(Hβ)=20
S/N(Hβ)=10
S/N(Hβ)=3
Figure 13. Same as Fig. 11 but for the parameter retrieval of
the H ii region and the NLR for objects at z = 0 with an equal
contribution of both regions to Hβ.
0.0
0.5
1.0
fNLR(Hβ)=90 fNLR(Hβ)=50 fNLR(Hβ)=10
0 10 20 30 40
S/N Hβ
0.0
0.5
1.0
m
ea
su
re
d
f N
L
R
(H
β
)
Figure 14. Retrieved fractional contribution of the NLR to the
total Hβ flux for each object in our simulated grid. The upper
(lower) panel shows the retrieval for the typical z = 0 (z = 2)
objects with different input fractional NLR contribution to Hβ as
displayed in the legend.
Chi-mistry code8 (Pérez-Montero 2014) to include charac-
terization of the NLR gas using a grid of photoionization
8 This code establishes a Bayesian-like comparison between the
predictions of the lines emitted in the ionized gas from a grid of
photoionization models covering a large range of input parameters.
In this comparison, the code does not assume any fixed relation
between secondary and primary elements.
MNRAS 000, 1–24 (2018)
16 Vidal-García et al.
models produced with cloudy v.17.01. The incident ion-
izing radiation consists of two components: that character-
izing the ‘big blue bump’ [the curved component peaking
at log(ν/Hz) ∼ 13.5 in Fig. 1], and a power-law slope be-
tween 2keV (6Å) and 5eV (2500Å) of αox = −0.8 (they also
compare to a model with αox = −1.2); power-law emission
with αx = −1 representing the non-thermal X-ray emission.
We show the incident radiation in Fig. 1. Pérez-Montero
et al. (2019) add dust within the NLR and model the gas
as homogeneously-distributed with filling factor 0.1 and hy-
drogen density, nh = 500 cm−3. All element abundances are
assumed to scale to solar, except for nitrogen, which is an
additional free parameter.
Hii-Chi-mistry uses standard line ratios to derive
physical parameters, and works by first constraining
N/O using the de-reddened ratios between [N ii]λ6584
and [O ii]λ3726, λ3729 and between [N ii]λ6584 and
[S ii]λ6717, λ6731, before constraining metallicity and
ionization parameter. The chemical composition of the
gas is characterized by the total oxygen abundance,
which covers the range 6.9 < 12 + log(O/H) < 9.1,
−2 < log(N/O) < 0 and the released version of the code
covers −2.5 < logUnlrs < −0.5. The ionization parameter
in these models is defined the same way as for beagle
(following description in section 3.1 of Pérez-Montero 2014).
The models consider default grain properties and relative
abundances, and all elements but Nitrogen are scaled to the
solar values given by Asplund et al. (2009) considering the
cloudy default depletion factors.
4.1.4 Thomas et al. 2018
NebularBayes, presented in Thomas et al. (2018a), is the
only code, other than beagle, which simultaneously fits the
star-forming and NLR contributions to nebular emission. The
NLR models are based on physically motivated prescriptions
for the incident radiation produced by the accretion disc
(Thomas et al. 2016). The accretion disc model has three
parameters: the energy of the peak of the accretion disc emis-
sion, Epeak; the photon index of the inverse Compton scat-
tered power-law tail, Γ (non-thermal X-ray radiation); and
the proportion of the total flux that goes into the non-thermal
tail, pNT . The last two parameters are fixed to the fiducial
values of Γ = 2.0 and pNT = 0.15.
The NLR grid included in NebularBayes is calculated for
constant pressure models using mappings v5.1 (Sutherland
& Dopita 2017). This is different from the constant density
models produced for the other works. Groves et al. (2004)
demonstrated the large difference in coverage of the BPT
plane that dusty isobaric models introduce. In particular, for
the optical line ratios explored in this paper, their isobaric
dusty model curves also better sample the space covered by
the data than their dust-free, constant density models. This
is due to the absorption of ionizing photons by dust, which
increases the radiation pressure (and hence density), espe-
cially at high values of the ionization parameter. The four
parameters of the NLR that can be derived using Nebular-
Bayes are: oxygen abundance of the gas; ionization param-
eter, logUnlrin , defined at the inner edge of the ionized gas;
peak of the ionizing radiation, Epeak; and pressure.
The grid of models includes 12 oxygen abundances in the
range −1.7 6 log(O/H) 6 0.54, using the oxygen abundance
scaling from Nicholls et al. (2017). Depletion onto dust grains
is considered based on iron being 97.8% depleted. The ion-
ization parameter at the inner edge of the nebula can take 11
uniformly spaced values in the range −4.2 6 logUnlrin 6 −0.2.
The initial gas pressure samples 12 values in the range
4.2 6 logP/k(cm−3K) 6 8.6 and values of Epeak sample six
values in the range −2.0 6 logEpeak(keV) 6 −0.75.
Thomas et al. (2018b) demonstrates Epeak cannot be
constrained when fitting to de-reddened [O ii]λ3726, λ3729,
[N iii]λ3869, Hβ, [Oiii]λ4363, [O iii]λ4959, λ5007, [O i]λ6300,
He iλ5876, Hα, [N ii]λ6584 and [S ii]λ6717, λ6731 (those avail-
able in the SDSS DR7 MPA JHU measured line fluxes), but
find that values in the range 40–50 eV give reasonable frac-
tions9 of NLR contribution, and so, in the following, we fix
the value of this parameter to Epeak = 45 eV.
4.2 Comparison to other works and data
We start by comparing our NLR-emission models (F16)
to those of Thomas et al. (2018a) and Pérez-Montero
et al. (2019), as well as observed line ratios of 44 ob-
served Seyfert II galaxies collated by Dors et al. (2017).
Fig. 15 shows gas-phase 12 + log(O/H) against two differ-
ent emission line ratios: R23 (defined in equation 7) and
log([O ii]λ3726, λ3729/[O iii]λ4959, λ5007). The lines and
matching coloured symbols show the coverage of different
models in these plots, while the black and grey symbols show
measured line ratios and estimated abundances for the ob-
served galaxies. For each of the observed galaxies, we show
three different abundance estimates: those from Dors et al.
(2017) (filled black circles) which were estimated by compar-
ing directly to cloudy photoionization models; and those
derived with each of the Storchi-Bergmann et al. (1998) cal-
ibrations, with filled grey (open black symbols) showing the
measurements obtained with the calibration in equation 4
(5). As noted by Storchi-Bergmann et al. (1998), the second
calibration (equation 5) leads to systematically higher oxygen
abundance estimates than the other two methods.
The models in Fig. 15 span a range of ionization parameters
ranging from high to low values with dark-to-faint symbols.
For the F16 models, we set αpl = −1.7. Thomas et al. (2018b)
fix Epeak to 45 eV, which is somewhat higher than the peak
incident radiation within our models. We investigated values
in the range 8 < Epeak/eV < 56 but in order not to crowd the
plot too much, we show the results for a central value within
this range of ∼ 17.8 eV. We also explored the full pressure
range 4.2 < log(P/k) < 8.6, but plot those with the central
value of log(P/k) = 7. Then, since the Pérez-Montero et al.
(2019) model grid encompasses a wide range of N/O abun-
dances, we fix to a value of log(N/O)=−0.5 (the same consid-
ered as the ‘standard’ one by Pérez-Montero et al. (2019)),
although the diagrams in Fig. 15 are not very sensitive to the
value chosen.
Fig. 15 (a) shows the oxygen abundance versus R23, which
is defined in equation 7. The F16 models cover the full range
of the observations, with the exception of a few objects with
the lowest values of R23. The Pérez-Montero et al. (2019)
9 Higher Epeak values give a very low fractional contribution of
NLR to Hβ line flux, whereas the values they choose span 0.1-1.0
along the AGN branch in the BPT diagrams.
MNRAS 000, 1–24 (2018)
beagle-agn 17
−1.5 −1.0 −0.5 0.0 0.5 1.0 1.5 2.0
log R23
7.6
7.8
8.0
8.2
8.4
8.6
8.8
9.0
9.2
12
+l
og
(O
/H
) Dors+2017
Dors+2017, SB1
Dors+2017, SB2
Feltre + 16 :
log UNLRs =−1.5,−2.5,−3.5,−4.0
Thomas + 18a(log p/k = 7.0, log Epeak =−1.75) :
log UNLRin =−0.6,−1.4,−2.6,−4.2
Perez−Montero + 19(N/O =−0.5) :
log UNLRs =−0.5,−1.5,−2.5
(a)
−2.5 −2.0 −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5 2.0
log([OII]λ3726,λ3729/[OIII]λ4959,λ5007)
7.6
7.8
8.0
8.2
8.4
8.6
8.8
9.0
9.2
12
+l
og
(O
/H
)
(b)
Figure 15. Gas-phase oxygen abundance (12 + log(O/H)) versus different emission line ratios, log(R23) (defined in equation 7) and
log([O ii]λ3726, λ3729/[O iii]λ4959, λ5007). In these figures we show the coverage of different NLR emission line models via lines and
symbols with a range of metallicity and ionization parameters. The symbols and lines range from fainter to darker colours from lower
to higher logUnlrs or logUnlrin respectively, as indicated in the legend. The black, grey and open circles show the measured emission
line ratios and inferred NLR 12 + log(O/H) values for a set of type-2 AGNs collated in Dors et al. (2017). Black circles display the
12 + log(O/H) estimates from Dors et al. (2017), while the grey and open circles display the 12 + log(O/H) estimates calculated with the
two Storchi-Bergmann et al. (1998) NLR calibrations, defined in equations 4 and 5.
models provide good coverage of the observations with low
values of R23, but provide only marginal coverage of the high-
est R23 values. We display the Thomas et al. (2018b) models
with a single value of Epeak, although have tested a range
from 8 to 56 eV, finding that the models cover the full range
of observed values, except for those few observations with
the lowest R23. It is important to note that the abundance
estimates of the observed objects plotted are not objectively
more secure than other abundance estimates presented in this
work. They are subject to their own modelling assumptions,
so we focus mostly on the coverage of the plotted models over
the emission line ratio space.
Fig. 15 (b) shows the gas-phase oxygen abundance
for the same models and data as in Fig. 15 (a), but
plotted against [O ii]λ3726, λ3729/[O iii]λ4959, λ5007. The
Thomas et al. (2018a) and Pérez-Montero et al. (2019)
models do not include [Oiii]λ4959, which we assume to
have 1/3 the flux of [O iii]λ5007. We see that for the
range of ionization parameters plotted, our models cover
the parameter space spanned by the observations for all
but the highest [O ii]λ3726, λ3729/[O iii]λ4959, λ5007 ratios.
The Pérez-Montero et al. (2019) models cover the high-
est [O ii]λ3726, λ3729/[O iii]λ4959, λ5007 values, but obser-
vations with [O ii]λ3726, λ3729/[O iii]λ4959, λ5007 . −0.5
are not covered. The Thomas et al. (2018a) models cover
the full range of [O ii]λ3726, λ3729/[O iii]λ4959, λ5007 (.-
0.5) from the observations.
Both panels of Fig. 15 show that the Pérez-Montero
et al. (2019) models give the opposite dependence of
R23 and [O ii]λ3726, λ3729/[O iii]λ4959, λ5007 on logUnlrs
to those of F16 and Thomas et al. (2018a). To
understand this, we plot logUnlrs as a function of
[O ii]λ3726, λ3729/[O iii]λ4959, λ5007 in Fig. 16. Pérez-
Montero et al. (2019) chose to publish only the highest
logUnlrs models, as the behaviour between logUnlrs and
[O ii]λ3726, λ3729/[O iii]λ5007 is double-valued in their mod-
els. We do not see the same behaviour in the F16 models
with αpl = −1.7, but when we plot our models with shal-
lower αpl = −1.2, we start to see a similar behaviour at
the highest logUnlrs and 12 + log(O/H) values (red symbols).
The strong increase in [O ii]λ3726, λ3729/[O iii]λ5007 with
logUnlrs of the Pérez-Montero et al. (2019) models is there-
fore presumably due to the choice of a hard incident ionizing
spectrum (Fig. 1). This is also likely the reason for the mod-
erate offset to higher [O ii]λ3726, λ3729/[O iii]λ5007 values at
given logUnlrs compared to the F16 models (given that the
partially ionized zone will be more extensive).
Including only the higher part of the logUnlrs fork, the
Pérez-Montero et al. (2019) models will measure an oppo-
site 12 + log(O/H) versus logUnlrs to that measured with
the Thomas et al. (2018b) or F16 models, which will in
turn provide systematically different oxygen abundances.
The data used in the Pérez-Montero et al. (2019) fit-
ting could not distinguish which branch of logUnlrs ver-
sus [O ii]λ3726, λ3729/[O iii]λ4959, λ5007 the data inhabit.
It would clearly be beneficial to better constrain the ion-
ization state of NLRs in observed type-2 AGN to moti-
vate the range of logUnlrs the models must span. Fig. 17
displays the ion fraction for a number of different species
(hydrogen, helium and oxygen) from the F16 models with
αpl = −1.2. The ion fraction are shown as a function of the
total hydrogen column density, NH(R) = nnlrh R. The high
energy photons from hard ionising radiation of the accretion
disc penetrates further through the gas than lower energy
photons, creating a partially ionized region. We might ex-
pect [O ii]λ3726, λ3729/[O iii]λ5007 to decrease with increas-
ing logUnlrs . However, as explained in Pérez-Montero et al.
(2019), and shown clearly in Fig. 17, for high ionization pa-
MNRAS 000, 1–24 (2018)
18 Vidal-García et al.
−2.5 −2.0 −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5
log([OII]λ3726,λ3729/[OIII]λ4959,λ5007)
−4.0
−3.5
−3.0
−2.5
−2.0
−1.5
−1.0
−0.5
lo
g
U
N
L
R
s
FeltreαPL =−1.7 :
12+log(O/H)=8.2,8.52,8.66,9.0
FeltreαPL =−1.2 :
12+log(O/H)=8.3,8.6,8.7,9.2
Perez-Montero (log N/O=-0.5):
12+log(O/H)=8.2,8.6,9.0
Figure 16. logUnlrs or logUnlrin versus
[O ii]λ3726, λ3729/[O iii]λ4959, λ5007 for the F16 and Pérez-
Montero et al. (2019) models for a range of different 12 + log(O/H)
values. We show the F16 models with two values of αpl: -1.2 (red)
and -1.7 (violet), as indicated in the legend.
rameters and hard ionizing radiation, this is not the case as
at small radii within the Strömgrem radius, [O iii]λ5007 is
lost to higher ionization states of oxygen such as O3+ and
O4+, while O+ is still present in the partially ionized zone.
There are no visible emission lines from these higher ioniza-
tion states in the wavelength range of SDSS. However, He ii
has a similar ionization energy to O3+ (54.4eV versus 54.9eV),
and the SDSS spectra provide coverage of the [He ii]λ4686 re-
combination line.
4.3 Nebular [He ii]λ4686 fluxes to constrain the
ionization state of the NLR in type-2 AGNs
within SDSS
We measure the nebular [He ii]λ4686 fluxes for 463 confirmed
type-2 AGNs presented in Dors et al. (2020a) from DR7 SDSS
spectra (Abazajian et al. 2009). [He ii]λ4686 is not delivered
in the MPA-JHU group10 distribution of measured line fluxes
used in Dors et al. (2020a), since measurements of this line
can be complicated by stellar wind features and broaden-
ing. We measured [He ii]λ4686 using the python package em-
cee (Foreman-Mackey et al. 2013). Since we wish to gain
constraints on the nebular (narrow-line) contribution to the
[He ii]λ4686 line, we tie the width of the line11 to the width
of [O iii]λ5007 in the same spectra, by simultaneously fitting
to the two lines. This effectively uses information from the
higher S/N [O iii]λ5007 line to provide a strong prior on the
shape of [He ii]λ4686 if it is to arise in the same physical re-
gion. We define continuum regions on either side of the two
lines and subtract a linear fit to the continuum before per-
forming the line fits. In the case where no obvious nebular
[He ii]λ4686 emission is present, the tied fitting will provide
firm constraints on the limits on the [He ii]λ4686 flux with the
10 https://wwwmpa.mpa-garching.mpg.de/SDSS/DR7/
11 We do not tie any other quantity, so, for instance, the velocity
offset is different for both lines.
Table 4. [He ii]λ4686 line fluxes measured for the sample of type-2
AGNs presented by Dors et al. (2020a). Full sample will be avail-
able online.
MJD-PLATE-FIBER [He ii]λ4686 flux width
(/10−16 erg s−1 cm−2) (/km s−1)
51882-0442-061 0.826± 0.270 151± 3
52056-0603-545 0.209± 0.145 130± 15
51691-0340-153 0.258± 0.192 268± 55
52079-0631-387 0.231± 0.152 125± 3
52017-0516-232 1.01± 0.26 194± 8
51915-0453-002 1.17± 0.30 231± 5
52049-0618-535 0.261± 0.184 121± 15
52000-0335-428 0.135± 0.116 117± 21
delivered uncertainties. We report the measured line fluxes
and line widths in Table 4.
We show [O ii]λ3726, λ3729/[O iii]λ4959, λ5007
against [He ii]λ4686/[O ii]λ3726, λ3729 (left) and
[He ii]λ4686/[O iii]λ4959, λ5007 (right) in Fig. 18, with
the F16 models for comparison. For this figure, we
used the MPA-JHU line fluxes for [O ii]λ3726, λ3729 and
[O iii]λ4959, λ5007. The sample of type-2 AGNs clearly follow
the general trend of the F16 models. Unfortunately the mod-
els of Pérez-Montero et al. (2019) do not include [He ii]λ4686,
so we cannot compare to their predicted ratios directly. If
[O ii]λ3726, λ3729/[O iii]λ4959, λ5007 were increasing with
logUnlrs (as for the Pérez-Montero et al. 2019 models in
Fig. 16), we would expect [He ii]λ4686/[O ii]λ3726, λ3729
to increase with [O ii]λ3726, λ3729/[O iii]λ4959, λ5007, as
well as much higher values of [He ii]λ4686/[O ii]λ3726, λ3729
for the observed galaxies. The observations show low
values of [He ii]λ4686/[O ii]λ3726, λ3729 and a negative
trend between [O ii]λ3726, λ3729/[O iii]λ4959, λ5007 and
[He ii]λ4686/[O ii]λ3726, λ3729 (Fig. 18, a), suggesting that
the set of observed type-2 AGNs have significantly lower
ionization parameter than covered by the Pérez-Montero
et al. (2019) models.
4.4 Comparison of the nitrogen and sulfur emission
lines
Because the [S ii]λ6717, λ6731 doublet is routinely used to
derive the gas density in photoionized regions, we have also
compared the predictions for this doublet from the different
NLR models and data presented in the previous sections. In
particular, we show in Fig. 19 the gas-phase oxygen abun-
dance as a function of the ratio [S ii]λ6717, λ6731/Hα. In
this case, we only show the data with the predictions for
12 + log(O/H) from Dors et al. 2017 in order to not crowd
the figure too much. The Pérez-Montero et al. 2019 models
reproduce the data except for the points with the highest
oxygen abundances and lowest [S ii]λ6717, λ6731/Hα ratios.
The coverage of the F16 and Thomas et al. 2018a is very
similar and almost complete except for three objects with
high [S ii]λ6717, λ6731/Hα ratios. We have noted (not shown
here) that for lower values of αpl F16 models cover better
this region of the parameter space.
Abundance studies of galactic and extragalactic H ii re-
gions find that nitrogen has a primary origin for the low-
metallicity regime (12 + log(O/H). 8.2) and a secondary
one for the high-metallicity regime (e.g. Mouhcine & Con-
tini 2002; Pérez-Montero & Contini 2009). Groves et al. 2004
MNRAS 000, 1–24 (2018)
beagle-agn 19
Figure 17. The ion fractions of H+ (black dashed), He2+ (black), O+ (red), O2+ (orange), O3+ (green) and O4+ (blue) as a function
of the total hydrogen column density. The NLR models have Znlrgas = 0.030, αpl = −1.2 and logUnlrs = −1.5 (left) and -2.5 (right).
−2.5 −2.0 −1.5 −1.0 −0.5 0.0
log(HeIIλ4686/[OII]λλ3726, 3729)
−2.0
−1.5
−1.0
−0.5
0.0
0.5
1.0
1.5
lo
g(
[O
II
]λ
λ
37
26
,3
72
9/
[O
II
I]
λ
λ
49
59
,5
00
7)
O/H=8.2,8.52,8.66,9.0
−2.5 −2.0 −1.5 −1.0 −0.5 0.0
log(HeIIλ4686/[OIII]λλ4959, 5007)
−2.0
−1.5
−1.0
−0.5
0.0
0.5
1.0
1.5
lo
g(
[O
II
]λ
λ
37
26
,3
72
9/
[O
II
I]
λ
λ
49
59
,5
00
7)
O/H=8.2,8.52,8.66,9.0
Figure 18. Assessing the ionization state of the gas within the type-2 AGN sample from Dors et al. (2020a). With our measurements of
the [He ii]λ4686 line from the SDSS spectra of these objects (see text for details), we compare log([O ii]λ3726, λ3729/[O iii]λ4959, λ5007)
versus log([He ii]λ4686/[O ii]λ3726, λ3729) [log([He ii]λ4686/[O iii]λ4959, λ5007)] to the F16 models with αpl=-1.7 in the left [right]
panels. The highest logUnlrs models have low [O ii]λ3726, λ3729/[O iii]λ4959, λ5007 ratios. The data shows a clear trend of decreas-
ing [O ii]λ3726, λ3729/[O iii]λ4959, λ5007 with increasing [He ii]λ4686/[O ii]λ3726, λ3729 or decreasing [He ii]λ4686/[O iii]λ4959, λ5007, as
shown by the F16 models.
find that the nitrogen abundance from the primary and sec-
ondary origin can be related to that of oxygen using an equa-
tion of the form of Equation 1. To capture the richness of
the variation of the N abundance as a function of O, we have
expanded the grid of NLR models. The fiducial grid has a
N/O dependence on the oxygen abundance following:
log(N/O)tot ≈ log{10−1.732 + 10[log(O/H)tot+2.19]}+ C (10)
with C = −0.25. Now we are parametrizing N/O in terms
of the total abundances, (N/O)tot (i.e. gas + dust) and its
dependence on the total oxygen abundance, (O/H)tot (gas +
dust). In expanding the grid we add two new values of the
normalisation, C, with C = 0 (the Nicholls et al. 2017 rela-
tion) and C = 0.5 (the upper limit). See Plat et al. (in prep.)
for the details of how these models have been computed.
We test the nitrogen abundance in the different grid of
models using the N2 and O3N2 line ratios, defined as follows:
N2 = log
(
[N ii]λ6584
Hα
)
(11)
O3N2 = log
(
[O iii]λ5007
Hβ
Hα
[N ii]λ6584
)
(12)
We show in Fig. 20 the gas phase oxygen abundance as a
function of these line ratios. For completeness, we show F16
and Pérez-Montero et al. 2019 models with different N/O val-
ues. In particular, for the Pérez-Montero et al. 2019 models,
we add the ones with the lowest and highest value of N/O
(cyan and dark blue points) and for the F16 models, the ones
with the highest (C = 0.5, pink points with dashed lines)
and lowest (C = −0.25 that of equation 1) scaling factor of
MNRAS 000, 1–24 (2018)
20 Vidal-García et al.
−3.0 −2.5 −2.0 −1.5 −1.0 −0.5 0.0 0.5
log ([S II]λ6717,λ6731/Hα)
7.6
7.8
8.0
8.2
8.4
8.6
8.8
9.0
9.2
12
+l
og
(O
/H
)
Feltre + 16 :
log UNLRs =−1.5,−2.5,−3.5,−4.0
Thomas + 18a(logp/k = 7.0, logEpeak =−1.75) :
log UNLRin =−0.6,−1.4,−2.6,−4.2
Perez−Montero + 19 :
log UNLRs =−0.5,−1.5,−2.5
Figure 19. Gas-phase oxygen abundance (12 + log(O/H)) as a
function of [S ii]λ6717, λ6731/Hα for the same models as in Fig. 15.
the N/O relations explored. We recall that nitrogen abun-
dance is not a free parameter in the Thomas et al. 2018a
grid. Most likely due to the fact that nitrogen abundances is
fixed, we see in Fig. 20a how Thomas et al. 2018a models are
not able to sample the data points with the highest nitrogen
abundances. However, increasing the N/O value allows for a
full and almost full coverage of the data for Pérez-Montero
et al. 2019 and F16 models respectively. In Fig. 20b the three
grids of models sample well the space of the diagram where
the data lie. We note however, that in the case of Thomas
et al. 2018a models, for data points with low O3N2 values, the
logUnlrs parameter might be underestimated. This might be
due to the lack of higher Nitrogen abundances that will tend
to be compensated in the O3N2 ratio with weaker [O iii]λ5007
lines, and therefore, lower logUnlrs values.
beagle does not yet allow for sampling over N/O, though
this comparison to observations shows that it is vital to ade-
quately explain measured line ratios including Nitrogen. It is
beyond the scope of this paper to implement and test retrieval
of parameters including the new N/O grid with our pedagog-
ical approach, but this functionality will be implemented for
future work.
4.5 Comparing to the Dors 2021 Te-based empirical
emission-line calibration
We plot the different theoretical models against the Dors
(2021) empirical emission-line calibration in Fig. 21. The cal-
ibration describes a surface in 12 + log(O/H)–P–R23 space
(as described in Eqs. 7 and 8 in Section 4.1.2). This provides
a quick comparison of the models to Te-based NLR oxygen-
abundance estimates without running full fits with beagle.
We will present full beagle fits to SDSS type-2 AGNs in
a future paper in this series. The black symbols display the
Te-based gas-phase oxygen abundance measurements for the
set of AGNs used to produce this empirical calibration, while
the grey plane shows the empirical calibration.
Fig. 21 (a) shows the F16 models plotted for a range of
logUnlrs , 12 + log(O/H), and two values of αpl (−1.2 and
−1.7). For a given αpl and 12 + log(O/H), different colours
represent different logUnlrs values, with dark-to-light colours
spanning high-to-low values. We see that at very low metal-
licities our models sit below the plane describing the Dors
(2021) empirical relation. Within the region occupied by the
observed galaxies, the models agree well with the empirical
relation until high metallicities, where the models curve up
from the empirical plane. The observations also preferentially
scatter above the relation. The behaviour of the models sug-
gests this scatter to higher oxygen abundances is real, not
simply intrinsic scatter about the empirical plane. There-
fore, objects above the relation have higher oxygen abun-
dances than can be explained by their position in the P-R23
plane, while the models can capture this behaviour. This be-
haviour is also observed in the Thomas et al. (2018a) constant
pressure models. The curvature from the plane is reminis-
cent to the effect of the addition of dust in NLR in Groves
et al. (2004) as the change from constant density to constant-
pressure models. However, our constant density models show
a similar behaviour, indicating that the observed curvature
might be due to the presence of dust, self consistently mod-
elled in cloudy.
The αpl = −1.7 models do not cover all the observed data
points in Fig. 21 (this might be best appreciated by view-
ing the interactive plot online, where different models can be
removed from the plot for clarity https://chart-studio.
plotly.com/~AlbaVidalGarcia/15/#/). A range of αpl ap-
pears to be required to fully cover the data. αpl = −1.2
covers the data well. With our idealized grid exploration in
Section 3, we found that the lines chosen for that analysis
could not fully constrain αpl, and that we recommend the
user to set it to αpl = −1.7 unless there is a significant im-
provement in the goodness of fit when αpl is allowed to vary.
We defer the search of other lines that can constrain better
this parameter to a later paper.
Panel (b) of Fig. 21 shows the same models, observations
and empirical relation as shown in panel (a), but we in-
clude also the Thomas et al. (2018a) and Pérez-Montero
et al. (2019) models. In particular, for the Thomas et al.
(2018a), we plot a range of −4.2 < logUnlrin < −0.2, with
dark-to-light colours representing high-to-low values, and for
three values of Epeak, as indicated in the legend. In particu-
lar, we show the one taken as standard throughout this pa-
per, -1.75 (green points) and a lower and a higher one, -2.0
(green diamonds) and -1.5 (green squares) respectively. For
the Pérez-Montero et al. (2019) models, we plot a range of
−2.5 < logUnlrs < −0.5 but we do not vary the N/O values
because this variation has a negligible impact on the P and
R23 values. This figure shows that all models curve off the
empirical relation to high metallicities.
This comparison to the Dors (2021) empirical metallicity
calibration shows the limitations of both the empirical diag-
nostic itself and our models. The empirical diagnostic cannot
capture the variation with metallicity above the plane. Our
models require a range of αpl to account for the data, a pa-
rameter we cannot constrain from the set of line ratios we
tested. Our tests in Section 3 further demonstrate how biased
other NLR gas parameters become when allowing αpl to vary.
However, below the empirical relation (lower metallicities),
a region of the parameter space we expect to find at higher
redshifts with missions such as JWST, the difference between
the two values of αpl is less significant for the F16 models.
It is beneficial to use Te-based O abundance estimates, em-
pirical diagnostics, and these types of sophisticated radiative
MNRAS 000, 1–24 (2018)
beagle-agn 21
−5 −4 −3 −2 −1 0 1
N2
7.6
7.8
8.0
8.2
8.4
8.6
8.8
9.0
9.2
12
+l
og
(O
/H
)
Feltre+16 (C=-0.25):
log UNLRs =−1.5,−4.0
Thomas+18 (log p/k=7.0, log Epeak=-1.75):
log UNLRin =−0.6,−4.2
Perez-Montero+19 (N/O=-0.5):
log UNLRs =−0.5,−2.5
Feltre+16 (C=0.5)
log UNLRs =−1.5,−4.0
Perez-Montero+19 (N/O=-2.0):
log UNLRs =−0.5,−2.5
Perez-Montero+19 (N/O=0.0):
log UNLRs =−0.5,−2.5
(a)
−3 −2 −1 0 1 2 3 4 5
O3N2
7.6
7.8
8.0
8.2
8.4
8.6
8.8
9.0
9.2
12
+l
og
(O
/H
)
(b)
Figure 20. Gas-phase oxygen abundance (12 + log(O/H)) as a function of N2 (a, see equation 11) and O3N2 (b, see equation 12) for the
same models as in Fig. 15 and different variations of the N/O relation: the Gutkin et al. 2016 relation (C = −0.25) and the upper limit
from eq. 10 (C = 0.5) for the F16 models and the highest and lowest values for nitrogen abundance of Pérez-Montero et al. 2019 grid.
The data points are from Dors et al. 2017.
(a) (b)
Figure 21. Dors (2021) empirical Te-based 12 + log(O/H) calibration, shown as the grey surface in 12 + log(O/H)–P–R23 space (where
P and R23 are defined in Eqs. 8 and 7, respectively). Both panels show a sample of observed AGNs with 12 + log(O/H) estimated
from the Te method by Dors (2021) (black points). The F16 models are shown in panel (a) for two values of αpl, and for a range of
ionization parameters spanning −4 < logUnlrs < −1.5 with dark-to-light colours spanning high-to-low logUnlrs values. Panel (b) includes
the Thomas et al. (2018a) and Pérez-Montero et al. (2019) models (see text for details). The 3D version of these figures can be found at
https://chart-studio.plotly.com/~AlbaVidalGarcia/15/#/.
transfer models together to tease out the limitations of each
approach. Until Te measurements are possible for large sam-
ples of type-2 AGNs at high redshift, the radiative transfer
models can capture the evolution in the population more ef-
fectively than an empirical calibration based on objects in a
small range of metallicities.
5 SUMMARY
In this paper we present the incorporation of the Feltre et al.
(2016) NLR models into beagle. This addition allows the
mixing of emission from H ii regions with that from the NLR
of type-2 AGNs. Dust attenuation is applied self-consistently
taking into account the dust within the NLR, H ii regions and
the diffuse ISM.
We take a pedagogical approach to defining which pa-
rameters of our model (from both H ii and NLR regions)
MNRAS 000, 1–24 (2018)
22 Vidal-García et al.
can be constrained by fitting a given set of observables in
the optical range. To determine this we work with idealised
z ∼ 0 and z ∼ 2 galaxy spectra with varying contribu-
tion of the NLR to Hβ flux. We fit to a set of line ra-
tios ([O i]λ6300/Hα, [N ii]λ6584/Hα, [S ii]λ6717, λ6731/Hα,
[O iii]λ5007/Hβ, Hβ/Hα, [O i]λ6300/[O ii]λ3726, λ3729 and
[O ii]λ3726, λ3729/[O iii]λ5007) as well as the Hα line flux.
We then quantify the S/N required to constrain the param-
eters of our model for different NLR contributions. Finally,
we compare the results of our model with those obtained us-
ing previously published methods. The main results are the
following:
• With the set of observables used, the retrieval of physical
parameters for the AGN dominated spectra with high S/N
shows degeneracies between ξnlrd and logU
nlr
s that lead to
an underestimation of logUnlrs and Znlrgas if ξnlrd is left as a
free parameter. We obtain similar results when adding to the
observables a ratio sensitive to the N/O abundance. Allowing
αpl to vary leads to an over-estimated retrieval of Lacc while
still under-estimating Znlrgas . We find the best recovery of Lacc,
Znlrgasand logUnlrs when fixing both ξnlrd and αpl.
• The retrieval of physical parameters for the star-
formation dominated spectra is also biased by the nebular
ξd. We find that ξd is poorly constrained and when left free
leads to an over-estimation of the logUs and Zhiigas, as well as
tighter constraints on the biased estimates.
• We also test the case of equal-contributions of NLR and
H ii regions to the Hβ flux. In this case, the accuracy of the
αpl parameter depends heavily on the region of the observ-
able parameter space that is being probed, so we suggest that
it should be fixed, unless varying it is required to reach an
acceptable fit. Both ξd and ξnlrd are degenerate with several
parameters when left free, biasing the recovery of these pa-
rameters. We obtain the best constraints for the objects when
both ξd and ξnlrd are fixed.
• We have investigated the S/N in Hβ required to derive
un-biased parameter estimates (while fixing ξd, ξnlrd and αpl).
AGN parameters are well constrained with a S/N(Hβ) ∼ 20
for the AGN-dominated galaxies. H ii region parameters are
well constrained with S/N(Hβ) ∼ 10 for star-formation dom-
inated galaxies. However, we need a higher S/N(Hβ) ∼ 30
to disentangle the contributions from star-forming and NLR
components when they have similar contributions to Hβ flux.
• We test how well beagle can detect the presence of an
AGN by retrieving the fractional contribution of the NLR to
the total Hβ flux of a galaxy. For the z = 0 object, we find
that a S/N∼10 already allows for identification of objects that
have a significant contribution to the flux, even for the most
challenging case of equal contributions to the flux from NLR
and H ii regions. For the z = 2 object, a NLR contribution
is detected but overestimated at S/N∼10, though we suggest
this is still a reasonable limit for which to find objects with
a significant NLR contribution.
• We have compared the predictions of our model to
those from different sets of models and data compiled in
the literature. These include the models from Thomas et al.
(2018a) and Pérez-Montero et al. 2019 and the calibra-
tions to derive gas-phase oxygen abundance from Storchi-
Bergmann et al. 1998 and Dors 2021 for the data com-
piled by these last authors. We find that in general, all
the models have a good coverage of the parameter space
where the data lie. However, we find an opposite depen-
dence of the Pérez-Montero et al. 2019 models of R23 and
[O ii]λ3726, λ3729/[O iii]λ4959, λ5007 on logUnlrs with re-
spect to those of F16 and Thomas et al. 2018a. This is pre-
sumably due to the choice made in Pérez-Montero et al. 2019
to publish only the highest logUs models, as the behaviour
between logUs and [O ii]λ3726, λ3729/[O iii]λ4959, λ5007 is
double-valued in their models.
• To constrain the ionization state of NLRs in observed
type-2 AGN we measure the [He ii]λ4686 recombination line
that has a similar ionization energy to Oiv, in 463 confirmed
type-2 AGNs from DR7 SDSS spectra. From the comparison
to the data we see how the data follow the trend of the F16
models and they suggest that these data have a siginificantly
lower ionization parameter values than those covered by the
Pérez-Montero et al. 2019 models.
• To better capture the dispersion of the N/H-O/H rela-
tion in the data, we create a grid of different scaling factors of
this relation that we compare to the data and the same set of
models from the literature. Whereas F16 and Pérez-Montero
et al. 2019 models sample well the region of the parameter
space where the data lies, the fact that Thomas et al. 2018a
models fix the Nitrogen abundance might lead to an under-
estimation of the logUnlrs parameter when fitting to nitrogen
emission lines, because models with extremely low logUnlrs
values are needed to sample the data.
• We also compare the models to the empirical emission-
line calibration from Dors 2021. Within the region sampled
by observations, the models agree well with the empirical
calibration. However, at high metallicities, the models curve
up from the plane that the empirical calibration draws in
the 12 + log(O/H)–P–R23 space. The observations are also
scattered preferentially above the plane in the region cov-
ered by the models. Therefore, it is likely the models can
capture the behaviour of the observations that also scatter
above the plane. At lower metallicities, the models lie be-
low the Dors 2021 calibration. Additionally, the difference
between models with different αpl values is less significant,
which will help constraining the gas parameters in galaxies
at high redshifts, such as those expected to find with mis-
sions such as the JWST. In general, the radiative transfer
models can capture the behaviour of objects with different
oxygen abundances more effectively than an empirical cali-
bration based on objects with a limited range of metallicities.
This is the first in a series of three papers, that will be
followed by a paper on fitting of a sample of type-2 AGNs
with beagle, and a study of the extent to which emission
from shocks and post-AGB stars may affect our inferences.
STATEMENT OF CONTRIBUTIONS BY
AUTHORS
This is a highly collaborative project led by ECL where ev-
ery author has had an important contribution. In particular
for this paper, AF and ECL incorporated the models into
BEAGLE. AVG and AP defined and fitted to the idealised
models in Section 3 from the parameters that MH provided
from the simulations. AVG worked on the discussion com-
paring to other models which led to ECL measuring HeII
line fluxes in SDSS data and also the need of creating a new
N/O grid done by AP and AF. The writing has been done by
MNRAS 000, 1–24 (2018)
beagle-agn 23
ECL with the help of AVG, and editing by SC. SC, MH and
JC have contributed to the scientific discussions and followed
the project from the beginning. If the field of astronomy had
the mechanism for assigning multiple first authors (as is done
in other fields), AVG, AP and ECL would have been assigned
as lead authors. We reflect this with the fact that they are
assigned as corresponding authors.
ACKNOWLEDGEMENTS
We thank the anonymous referee for the thorough reading
of the manuscript and their timely suggestions that have
improved this work, in particular, making it more under-
standable. The authors acknowledge support from the Euro-
pean Research Council (ERC) via an Advanced Grant un-
der grant agreement no. 321323-NEOGAL. AVG also ac-
knowledges support from the European Research Council
through the Advanced Grant MIST (FP7/2017-2020, No
742719). ECL acknowledges support of an STFC Webb Fel-
lowship (ST/W001438/1). MH acknowledges financial sup-
port from the Swiss National Science Foundation via the
PRIMA grant “From cosmic dawn to high noon: the role of
BHs for young galaxies” (PROOP2_193577). AF acknowl-
edges support from grant PRIN MIUR 2017 20173ML3WW.
JC acknowledges funding from the ERC Advanced Grant
789056 “FirstGalaxies” (under the European Union’s Hori-
zon 2020 research and innovation programme).PySpecLines
uses the python package PySpecKit https://pyspeckit.
readthedocs.io/en/latest/index.html
DATA AVAILABILITY
The astrophysical data used in this publication is supplied
within the paper and the online supplementary materials.
The idealised spectra analysed in Section 3 will be shared
with reasonable request to the corresponding authors.
REFERENCES
Abazajian K. N., et al., 2009, ApJS, 182, 543
Agostino C. J., et al., 2021, ApJ, 922, 156
Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A,
47, 481
Baldwin J. A., Phillips M. M., Terlevich R., 1981, PASP, 93, 5
Baskin A., Laor A., 2005, MNRAS, 358, 1043
Bruzual G., Charlot S., 2003, MNRAS, 344, 1000
Bundy K., et al., 2015, ApJ, 798, 7
Castro C. S., Dors O. L., Cardaci M. V., Hägele G. F., 2017, MN-
RAS, 467, 1507
Charlot S., Fall S. M., 2000, ApJ, 539, 718
Charlot S., Longhetti M., 2001, MNRAS, 323, 887
Chevallard J., Charlot S., 2016, MNRAS, 462, 1415
Cowie L. L., Songaila A., 1986, ARA&A, 24, 499
Crummy J., Fabian A. C., Gallo L., Ross R. R., 2006, MNRAS,
365, 1067
Curtis-Lake E., Chevallard J., Charlot S., Sandles L., 2021, MN-
RAS, 503, 4855
Díaz-Francés E., Rubio F. J., 2013, Statistical Papers, 54, 309
Dietrich M., Hamann F., Shields J. C., Constantin A., Heidt J.,
Jäger K., Vestergaard M., Wagner S. J., 2003, ApJ, 589, 722
Done C., Davis S. W., Jin C., Blaes O., Ward M., 2012, MNRAS,
420, 1848
Dors O. L., 2021, MNRAS, 507, 466
Dors O. L., Cardaci M. V., Hägele G. F., Krabbe Â. C., 2014,
MNRAS, 443, 1291
Dors O. L., Cardaci M. V., Hägele G. F., Rodrigues I., Grebel
E. K., Pilyugin L. S., Freitas-Lemes P., Krabbe A. C., 2015,
MNRAS, 453, 4102
Dors O. L. J., Arellano-Córdova K. Z., Cardaci M. V., Hägele G. F.,
2017, MNRAS, 468, L113
Dors O. L., et al., 2020a, MNRAS, 492, 468
Dors O. L., Maiolino R., Cardaci M. V., Hägele G. F., Krabbe
A. C., Pérez-Montero E., Armah M., 2020b, MNRAS, 496, 3209
Feltre A., Hatziminaoglou E., Fritz J., Franceschini A., 2012, MN-
RAS, 426, 120
Feltre A., Charlot S., Gutkin J., 2016, MNRAS, 456, 3354
Ferland G. J., et al., 2013, Rev. Mex. Astron. Astrofis., 49, 137
Feroz F., Hobson M. P., Bridges M., 2009, MNRAS, 398, 1601
Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013,
PASP, 125, 306
Groves B. A., Dopita M. A., Sutherland R. S., 2004, ApJS, 153,
75
Groves B., Heckman T., Kauffmann G., 2006, MNRAS, 371, 1559
Gutkin J., Charlot S., Bruzual G., 2016, MNRAS, 462, 1757
Haardt F., Maraschi L., 1991, ApJ, 380, L51
Hamann F., Ferland G., 1993, ApJ, 418, 11
Hirschmann M., Charlot S., Feltre A., Naab T., Choi E., Ostriker
J. P., Somerville R. S., 2017, MNRAS, 472, 2468
Hirschmann M., Charlot S., Feltre A., Naab T., Somerville R. S.,
Choi E., 2019, MNRAS, 487, 333
Kewley L. J., Nicholls D. C., Sutherland R. S., 2019, ARA&A, 57,
511
Maiolino R., Mannucci F., 2019, A&ARv, 27, 3
Marinacci F., et al., 2018, MNRAS, 480, 5113
Mathews W. G., Ferland G. J., 1987, ApJ, 323, 456
Matsuoka K., Nagao T., Maiolino R., Marconi A., Taniguchi Y.,
2009, A&A, 503, 721
Matsuoka K., Nagao T., Marconi A., Maiolino R., Mannucci F.,
Cresci G., Terao K., Ikeda H., 2018, A&A, 616, L4
Mignoli M., et al., 2019, A&A, 626, A9
Mouhcine M., Contini T., 2002, A&A, 389, 106
Nagao T., Marconi A., Maiolino R., 2006a, A&A, 447, 157
Nagao T., Maiolino R., Marconi A., 2006b, A&A, 447, 863
Naiman J. P., et al., 2018, MNRAS, 477, 1206
Nelson D., et al., 2018, MNRAS, 475, 624
Nicholls D. C., Sutherland R. S., Dopita M. A., Kewley L. J.,
Groves B. A., 2017, MNRAS, 466, 4403
Osterbrock D. E., Ferland G. J., 2006, Astrophysics of gaseous
nebulae and active galactic nuclei
Pérez-Montero E., 2014, MNRAS, 441, 2663
Pérez-Montero E., Contini T., 2009, MNRAS, 398, 949
Pérez-Montero E., Dors O. L., Vílchez J. M., García-Benito R.,
Cardaci M. V., Hägele G. F., 2019, MNRAS, 489, 2652
Pillepich A., et al., 2018, MNRAS, 475, 648
Pilyugin L. S., 2000, A&A, 362, 325
Pilyugin L. S., 2001, A&A, 369, 594
Reines A. E., Volonteri M., 2015, ApJ, 813, 82
Schartmann M., Meisenheimer K., Camenzind M., Wolf S., Hen-
ning T., 2005, A&A, 437, 861
Shemmer O., Netzer H., 2002, ApJ, 567, L19
Skilling J., et al., 2006, Bayesian Analysis, 1, 833
Springel V., et al., 2018, MNRAS, 475, 676
Storchi-Bergmann T., Schmitt H. R., Calzetti D., Kinney A. L.,
1998, AJ, 115, 909
Sutherland R. S., Dopita M. A., 2017, ApJS, 229, 34
Thomas A. D., Groves B. A., Sutherland R. S., Dopita M. A.,
Kewley L. J., Jin C., 2016, ApJ, 833, 266
Thomas A. D., Dopita M. A., Kewley L. J., Groves B. A., Suther-
land R. S., Hopkins A. M., Blanc G. A., 2018a, ApJ, 856, 89
MNRAS 000, 1–24 (2018)
24 Vidal-García et al.
Thomas A. D., Kewley L. J., Dopita M. A., Groves B. A., Hopkins
A. M., Sutherland R. S., 2018b, ApJ, 861, L2
Thomas A. D., Kewley L. J., Dopita M. A., Groves B. A., Hopkins
A. M., Sutherland R. S., 2019, ApJ, 874, 100
Tomczak A. R., et al., 2014, ApJ, 783, 85
Veilleux S., Osterbrock D. E., 1987, ApJS, 63, 295
Vidal-García A., Charlot S., Bruzual G., Hubeny I., 2017, MNRAS,
470, 3532
York D. G., et al., 2000, AJ, 120, 1579
do Nascimento J. C., et al., 2022, arXiv e-prints, p.
arXiv:2203.08602
APPENDIX A: TRIANGLE PLOTS OF THE
DIFFERENT RETRIEVED PARAMETERS
In this Section we show the 1D marginalised and 2D posterior
probability distributions of one of the noisy realisations of
the 10 fits populating the ellipses of Figs.5-10. Fig. A1 shows
the PDF of the fits with ξnlrd and αpl fixed for the NLR-
dominated galaxies at z = 0 and z = 2. Fig. A2 shows the
fit with ξd fixed for the SF-dominated galaxies at z = 0 and
z = 2. Finally, we show the triangle plot of the parameter
retrieval for the galaxy with equal contribution from SF and
NLR to Hβ flux at z = 0 in Fig. 9 and at z = 2 in Fig. 10.
MNRAS 000, 1–24 (2018)
beagle-agn 25
45.325 45.375 45.425
0.28 0.29
3.995 3.990
0.2 0.4 0.6 0.8
0.2 0.4 0.6 0.8
0.28
0.29
lo
g(
Z
N
LR
ga
s
/Z
)
3.995
3.990
lo
gU
N
LR
s
0.2
0.4
0.6
0.8
lo
g
V
45.325 45.375 45.425
Lacc
0.2
0.4
0.6
0.8
0.28 0.29
log(ZNLRgas /Z )
3.995 3.990
logUNLRs
0.2 0.4 0.6 0.8
log V
(a) z = 0
46.50 46.55
0.15 0.12
2.55 2.52
0.4 0.8 1.2 1.6
0.2 0.4 0.6 0.8
0.15
0.12
lo
g(
Z
N
LR
ga
s
/Z
)
2.55
2.52
lo
gU
N
LR
s
0.4
0.8
1.2
1.6
lo
g
V
46.50 46.55
Lacc
0.2
0.4
0.6
0.8
0.15 0.12
log(ZNLRgas /Z )
2.55 2.52
logUNLRs
0.4 0.8 1.2 1.6
log V
(b) z = 2
Figure A1. Triangle plots of the parameter retrieval for one of the noisy realisations for the z = 0 (left) and z = 2 (right) NLR-dominated
galaxies with S/N∼ 100 in Hβ. Only NLR parameters are displayed even though SF parameters are also included in the fits. We display
the results for the fit with ξnlrd and αpl fixed, which corresponds to the purple ellipses of Fig. 5.
0.400.440.480.520.56
3.992 3.976
0.02 0.02 0.06
0.1 0.2 0.3 0.4
0.2 0.4 0.6 0.8
3.992
3.976
lo
gU
s
0.02
0.02
0.06
lo
g(
Z
H
II
ga
s/Z
)
0.1
0.2
0.3
0.4
lo
g
V
0.400.440.480.520.56
log( /M yr 1)
0.2
0.4
0.6
0.8
3.992 3.976
logUs
0.02 0.02 0.06
log(ZHIIgas/Z )
0.1 0.2 0.3 0.4
log V
(a) z = 0
1.44 1.48 1.52 1.56
2.58 2.46 2.34
0.5 0.4 0.3
0.12 0.18 0.24 0.30
0.2 0.4 0.6 0.8
2.58
2.46
2.34
lo
gU
s
0.5
0.4
0.3
lo
g(
Z
H
II
ga
s/Z
)
0.12
0.18
0.24
0.30
lo
g
V
1.44 1.48 1.52 1.56
log( /M yr 1)
0.2
0.4
0.6
0.8
2.58 2.46 2.34
logUs
0.5 0.4 0.3
log(ZHIIgas/Z )
0.12 0.18 0.24 0.30
log V
(b) z = 2
Figure A2. As in Fig. A1, but showing triangle plots of the parameter retrieval for the z = 0 (left) and z = 2 (right) star-formation
dominated galaxies. We display only SF parameters, though NLR parameters were also included in the fits. These are the results of the
fits when ξd is fixed, which correspond to the blue ellipses of Fig. 8.
MNRAS 000, 1–24 (2018)
26 Vidal-García et al.
44.32 44.40
0.276 0.288
3.996 3.988
0.12 0.18 0.24 0.30
0.2 0.4 0.6 0.8
0.42 0.45 0.48 0.51
3.992 3.976
0.00 0.05
log(ZHIIgas/Z )
0.276
0.288
lo
g(
Z
N
LR
ga
s
/Z
)
3.996
3.988
lo
gU
N
LR
s
0.12
0.18
0.24
0.30
lo
g
V
0.2
0.4
0.6
0.8
0.42
0.45
0.48
0.51
lo
g(
/M
yr
1 )
3.992
3.976
lo
gU
s
44.32 44.40
Lacc
0.00
0.05
lo
g(
Z
H
II
ga
s/Z
)
0.276 0.288
log(ZNLRgas /Z )
3.996 3.988
logUNLRs
0.12 0.18 0.24 0.30
log V
0.2 0.4 0.6 0.8 0.42 0.45 0.48 0.51
log( /M yr 1)
3.992 3.976
logUs
Figure A3. As in Fig. A1 but now showing the triangle plots of the parameter retrieval for the z = 0 galaxy with equal contribution
from SF and NLR to Hβ flux. We display the results for all fixed parameters, corresponding to the black ellipses of Fig. 9.
MNRAS 000, 1–24 (2018)
beagle-agn 27
45.52 45.60
0.16 0.08
2.52 2.40
0.080.160.240.320.40
0.2 0.4 0.6 0.8
1.321.381.441.501.56
2.7 2.6 2.5 2.4
0.6 0.3
log(ZHIIgas/Z )
0.16
0.08
lo
g(
Z
N
LR
ga
s
/Z
)
2.52
2.40
lo
gU
N
LR
s
0.08
0.16
0.24
0.32
0.40
lo
g
V
0.2
0.4
0.6
0.8
1.32
1.38
1.44
1.50
1.56
lo
g(
/M
yr
1 )
2.7
2.6
2.5
2.4
lo
gU
s
45.52 45.60
Lacc
0.6
0.3
lo
g(
Z
H
II
ga
s/Z
)
0.16 0.08
log(ZNLRgas /Z )
2.52 2.40
logUNLRs
0.080.160.240.320.40
log V
0.2 0.4 0.6 0.8 1.321.381.441.501.56
log( /M yr 1)
2.7 2.6 2.5 2.4
logUs
Figure A4. As in fig. A1 but now showing the triangle plots of the parameter retrieval for the z = 2 galaxy with equal contribution from
SF and NLR to Hβ flux.
MNRAS 000, 1–24 (2018)