11institutetext:  Max Planck Institute for Extraterrestrial Physics,
Giessenbachstrasse, 85741 Garching, Germany 22institutetext: Excellence
Cluster Universe, Boltzmannstr. 2, D-85748, Garching, Germany 33institutetext:
School of Physics and Astronomy, Tel Aviv University, Tel Aviv 69978, Israel
44institutetext: Institute for Astronomy & Astrophysics, National Observatory
of Athens, V. Paulou & I. Metaxa, 11532, Greece 55institutetext: Centre for
Extragalactic Astronomy, Department of Physics, Durham University, United
Kingdom 66institutetext: Institute of Astronomy and Astrophysics, University
of Tübingen, Sand 1, 72076 Tübingen, Germany 77institutetext: Ludwig-
Maximilians-Universität München, Faculty of Physics, Scheinerstrasse 1,
D-81679 München, Germany 88institutetext: Department of Astronomy, University
of Illinois at Urbana-Champaign, Urbana, IL 61801, USA 99institutetext: Kavli
Institute for the Physics and Mathematics of the Universe, The University of
Tokyo, Kashiwa 277-8583 (Kavli IPMU, WPI), Japan

# Genuine Retrieval of the AGN Host Stellar Population (GRAHSP)

Johannes Buchner<EMAIL_ADDRESS>Hattie Starck 11 Mara Salvato 1122
Hagai Netzer 33 Zsofi Igo 11 Brivael Laloux 4455 Antonis Georgakakis 44
Isabelle Gauger 11 Anna Olechowska 1177 Nicolas Lopez 11 Suraj D Shankar 66
Junyao Li 8899 Kirpal Nandra 11 Andrea Merloni 11

(Received xx, Accepted xx)

###### Abstract

Context.  The assembly and co-evolution of super-massive black holes (SMBH)
and their host galaxy stellar population is one of the key open questions in
modern galaxy evolution. Observationally constraining this question is
challenging. Important parameters of galaxies, such as the stellar mass
($M_{\star}$) and star formation rate (SFR), are inferred by modeling the
spectral energy distribution (SED), with templates constructed on the basis of
various assumptions on stellar evolution. In the case of galaxies triggering
SMBH activity, the active galactic nucleus (AGN) contaminates the light of the
host galaxy at all wavelengths, hampering inference of host galaxy parameters.
Underestimating the AGN contribution due to incomplete AGN templates results
in a systematic overestimation of the stellar mass, biasing our understanding
of AGN and galaxy co-evolution. This challenge has gained further impetus with
the advent of sensitive wide-area surveys with millions of newly detected
luminous AGN, including by eROSITA, Euclid and LSST.

Aims.  We aim to estimate robustly the accuracy, bias, scatter and uncertainty
of AGN host galaxy parameters including stellar masses, and improve these
measurements relative to previously used techniques.

Methods.  This work makes two important contributions. Firstly, we present a
new SED fitting code, with an AGN model composed of a flexible power-law
continuum with empirically determined broad and narrow lines and a FeII forest
component, a flexible infrared torus that can reproduce the diverse dust
temperature distributions, and appropriate attenuation on the galaxy and AGN
light components. We verify that this model reproduces published X-ray to
infrared SEDs of AGN to better than 20 per cent accuracy. A fully Bayesian fit
includes uncertainties in the model and the data, making the inference highly
robust. The model is constrained with a fast nested sampling inference
procedure supporting the many free model parameters. Secondly, we created a
benchmark photometric dataset where optically-selected pure quasars are paired
with non-AGN pure galaxies at the same redshift. Their photometry flux is
summed into a hybrid (Chimera) object but with known galaxy and AGN
properties. Based on this data-driven benchmark, true and retrieved stellar
masses, SFR and AGN luminosities can be compared, allowing the evaluation and
quantification of biases and uncertainties inherent in any given SED fitting
methodology.

Results.  The Chimera benchmark, which we release with this paper, shows that
previous codes systematically over-estimate $M_{\star}$ and SFR by 0.5 dex
with a wide scatter of 0.7 dex, at AGN luminosities above
$10^{44}\unit{\per}$. In 20 percent of cases, the estimated error bars lie
completely outside a 1 dex wide band centered around the true value, which we
consider an outlier. In contrast, GRAHSP shows no measurable biased on
$M_{\star}$ and SFR, with an outlier fraction of only about 5 percent. GRAHSP
also estimates more realistic uncertainties.

Conclusions. Unbiased characterization of galaxy hosting AGN enables
characterization of the environmental conditions conducive to black hole
growth, whether star formation is suppressed at high black hole activity, and
identify mechanisms that prevent over-luminous AGN relative to the host galaxy
mass. It can also shed light on the long-standing question whether AGN
obscuration is primarily an orientation effect or related to phases in galaxy
evolution.

##  1 Introduction

The scaling relations of supermassive black holes (SMBH) at the center of
massive, passive galaxies presented in Magorrian et al. (1998), together with
the recognition that these are grown by active galactic nuclei (AGN) episodes
(Soltan, 1982) triggered a revolution in studies of galaxy evolution. It
implies that galaxies and AGN are merely experiencing varying levels of SMBH
accretion activity. Only combined, these phases give a full picture of galaxy
evolution. Further studies extended SMBH scaling relations to all types of
galaxies (e.g., early- and late-type, dwarf and active galaxies Gebhardt et
al., 2000; Ferrarese & Merritt, 2000; Xiao et al., 2011; Gültekin et al.,
2011; Kormendy & Ho, 2013; Schutte et al., 2019), finding correlations between
the black hole mass ($M_{\mathrm{BH}}$) and host physical parameters (e.g.,
stellar mass and velocity dispersion). These facts, combined with the
similarity of the cosmic star formation rate (SFR) and BH accretion history
(e.g. Aird et al., 2012; Madau & Dickinson, 2014) points toward the
interconnection between the SMBHs and their host and require considering the
AGN phase in galaxy evolution studies. Indeed, models have been proposed in
which the energy produced during SMBH growth episodes can modulate the
formation of stars in the host galaxy (e.g., Sanders et al., 1988), thereby
establishing AGN as an important component for our understanding of galaxy
evolution (e.g., Di Matteo et al., 2005; Hopkins et al., 2008). In turn, this
complicates a task that was already difficult for inactive galaxies. Unlike
redshift, estimates of galaxy physical parameters such as stellar mass, star
formation rate, and bolometric luminosity cannot be measured directly. Known
degeneracies like age-metallicity Worthey et al. (1994) or color-redshift
Masters et al. (2015) are hard to break with SED fitting, in part because the
number of parameters in the fit is usually larger than the number of
photometric data points.

For galaxies hosting an active SMBH, the need to account for the nuclear
emission at various wavelengths complicates things further. The relative
host/AGN contribution is unknown, and the AGN contribution at various
wavelengths changes, depending on whether or not the nuclear emission is
obscured by dust. Thus, basic questions such as “Are the physical properties
different for galaxies hosting an inactive/active SMBH?”, “Are the differences
between obscured and unobscured sources dependent on orientation or
evolution?” remain unanswered. Not correctly accounting for the nuclear
component can result in inaccurate and/or bias galaxy parameters and hence
interpretations, and could be the origin of conflicting results reported in
the literature so far. For example, Stemo et al. (2020), using a sample of
about 2500 X-ray selected AGN, found that AGN host galaxies tend to have lower
SFR than normal galaxies of equal mass. This result is in contrast with what
was found in Mountrichas et al. (2022), Pouliasis et al. (2022) Suh et al.
(2019) for the same type of AGN. It has also been reported that the stellar
mass of AGN host galaxies depends on their level of obscuration (suggesting an
evolutionary link; Koutoulidis et al., 2022; Zou et al., 2019) and in other
work that it does not (consistent with the AGN unification model; Masoura et
al., 2021). If there is an evolution in place from obscured to unobscured
objects as suggested by Hopkins et al. (2008), how can it be that the hosts of
obscured AGN are more massive than the host of the unobscured ones
(Koutoulidis et al., 2022)? This uncertainty in the actual stellar mass of AGN
host galaxies indirectly affects other parameters, such as specific accretion-
rate (used to estimate the rate at which the BH accrete material Aird et al.,
2012; Georgakakis et al., 2017) and specific star formation.

The increased availability of wide- and all-sky areas allows the selection of
bright, unobscured AGN. Contrary to typical pencil-beam surveys (e.g., Lockman
Hole (Fotopoulou et al., 2016), Chandra deep-field south (Hsu et al., 2014;
Luo et al., 2017)), where bright, unobscured AGN are virtually absent, they
are abundant in all-sky surveys. To mention already existing surveys, these
include optical spectroscopic by the Sloan Digital Sky Survey SDSS (Eisenstein
et al., 2011), photometry imaging by the Dark Energy Spectroscopic Instrument
Legacy Surveys (Dey et al., 2019), mid-infrared by WISE (Wright et al., 2010;
Mainzer et al., 2011; Meisner et al., 2023), astrometry by Gaia (e.g., Gaia
Collaboration et al., 2023), X-rays by 4XMM (Webb et al., 2020), XMMSlew
(Saxton et al., 2008), and eROSITA (Predehl et al., 2021). Therefore, a major
effort is required to develop algorithms that provide more reliable physical
parameters for unobscured AGN, finally allowing for inclusive galaxy
evolutionary studies. In recent years an increasing number of algorithms have
been developed (and often made public) for this purpose, such as AGNFitter
(Calistro Rivera et al., 2016), MAGPHYS (da Cunha et al., 2008), PROSPECTOR
(Johnson et al., 2021), CIGALE (Boquien et al., 2019), X-CIGALE (extending
CIGALE to X-ray emission; Yang et al., 2020), FAST (Aird et al., 2017), among
others (see review by Pacifici et al., 2023). Figure 1 of Thorne et al. (2022)
nicely summarises the basic ingredients of the most used codes111More
illustrations are available at https://www.astrojess.com/graphics/interactive-
sed-diagram..

This paper demonstrates that the existing approaches are insufficient in the
new parameter space, and present a new code, GRAHSP (Genuine Retrieval of the
AGN Host Stellar Population). GRAHSP differs substantially through a more
complete model description of the AGN component and more realistic uncertainty
estimation. GRAHSP has been developed by the eROSITA team specifically for
measuring the physical parameters of the millions of AGN that eROSITA started
to detect in 2020 (see Merloni et al., 2024, that presents the first all-sky
survey data release). The reliability of measurements obtained with GRAHSP at
different host-to-AGN ratios has been tested thoroughly by measuring the
physical parameters of AGN constructed by combining galaxies with known
stellar mass from COSMOS with a QSO. We have defined this testing sample as
the CHIMERA sample, an ideal benchmark for reference. For non-AGN galaxies we
verify that GRAHSP results are consistent with published results. This ensures
that the use can be generalised to any sample of extragalactic sources. In
addition to stellar mass, GRAHSP estimates key galaxy parameters such as the
star-formation rate, AGN parameters such as the bolometric luminosity, power-
law continuum slope, torus covering factor, and the system’s dust attenuation.

In Section 2 we present GRAHSP. The samples and data sets used in this work
are described in Section 3, which includes AGN-free galaxy samples and
infrared, X-ray, and optically selected AGN samples. Section 4 describe a
novel benchmark data set, where the true decomposition between galaxy and AGN
light is known by construction. In Section 5, the model is vetted with initial
tests. Based on the benchmark data set, section 6 demonstrates the key result:
GRAHSP’s stellar mass estimates are unbiased and outperform literature
methods, even in the difficult case of unobscured, luminous AGN. Throughout
the paper, we assume AB magnitudes unless stated otherwise and consistently
use the cosmology of Planck Collaboration et al. (2020).

##  2 Method

The modeling of multi-wavelength photometry of AGN candidates can be driven by
several goals: (1) the determination of redshifts (photo-z); (2) the
identification of whether an AGN is present; (3) the measurement of AGN
parameters, such as luminosities, and (4) the measurement of galaxy
parameters.

If the goal is to estimate redshifts, the colour-redshift relation can be
mapped with a small set of galaxy, AGN, and hybrid templates (e.g. Salvato et
al., 2009, 2019). If the goal is to identify AGN, light at any wavelength in
excess of that expected by galaxy emission needs to be detected, e.g., by the
statistical preference for an additional blue power-law (associated with the
accretion disk) or the presence of mid-infrared bump components (related
nuclear dust heated by the AGN). Two commonly adopted models for these two
components are those proposed by Richards et al. (2006) and Dale et al.
(2014), respectively; see Thorne et al. (2022) for a recent work. The goal can
also be achieved with diagnostic colour plots in the UV (e.g. Richards et al.,
2009) or infrared (e.g. Lacy et al., 2004).

If the goal is the inference of AGN parameters, such as those of a clumpy
torus (Nenkova et al., 2008a; Hönig & Kishimoto, 2010; Stalevski et al., 2016)
or the estimation of the black hole spin assuming a certain accretion flow
model (Netzer & Trakhtenbrot, 2014), one may assume an AGN model and interpret
the inferred parameters within this framework for the scope of the study. The
most challenging goal is the robust inference of galaxy parameters under the
presence of AGN light. This is the goal tackled in this work.

Suppose, for example, the case of photometry from an extremely luminous quasar
being modeled by an AGN template library222Current codes all try to account
for AGN contribution. Modelling without an AGN, can be considered the extreme
case of the described scenario.. If any residuals remain at this stage due to
an imperfect AGN modelling, and a galaxy template library is added and
optimized to fit the residuals, the stellar population will be extremely
massive (to reach quasar luminosities) and either red (to fix residuals in the
near-infrared) or blue (to fix residuals in the UV-optical). The inferred
parameters will be driven by the AGN mis-modelling and have nothing to do with
the underlying galaxy. In Section 6, we will show that this hypothetical
scenario occurs in real applications leading to biased measurements of key
physical parameters. While the inference would be wrong, this may not be
detectable by residual diagnostics such as reduced $\chi^{2}$, since there may
not be any residuals remaining. Also, the true stellar mass cannot be easily
verified by other, independent methods. The over-estimated galaxy template
normalisation cannot be resolved by refining the AGN template normalisation
further with data from longer and shorter wavelengths, because the AGN
template is already in the right place. This includes incorporating X-ray
information (e.g., Yang et al., 2020). This failure mode has the worrying
potential to induce spurious correlations between host galaxy mass and AGN
accretion rate.

In addition to the challenges above, AGN are ubiquitously variable (Bershady
et al., 1998; Klesman & Sarajedini, 2007; Sesar et al., 2007; Trevese et al.,
2008). While galaxies are, for the most part, static light emitters on human
time scales; and their SEDs can be modelled with a single static template, AGN
break this assumption. Therefore, residuals caused by non-simultaneously
collected photometry should also be considered in the modelling.

Figure 1: Overview of how the individual model components contribute to the
summed emission (black). The AGN power-law continuum (blue; accretion disk, or
BBB) is enhanced by emission lines including an iron forest (red). The disk
model is normalised at the monochromatic luminosity at 5100
$\mathrm{\text{Å}}$ (blue square), here
$L_{\mathrm{AGN}}^{5100$\mathrm{\text{Å}}$}=\qty{e44}{\per}=\qty{e37}{}$. The
torus (yellow dashed curve) typically dominates the disk continuum above
approximately $1$\mathrm{\SIUnitSymbolMicro m}$$. The torus is normalised by
the luminosity ratio 5100 $\mathrm{\text{Å}}$ to 12$\mathrm{\SIUnitSymbolMicro
m}$ (yellow diamond). The galaxy components include a stellar population
(purple), nebular emission lines (gray at the bottom), and galaxy dust
emission (green at the bottom right).  Figure 2: Detailed view of the optical
continuum and emission lines. The model continuum bending power-law (yellow
dashed line) reproduces the polarization measurements of Kishimoto et al.
(2008) (dark blue data points). The full model including a torus and emission
lines is shown in solid lines. Variations of the power-law slope (light grey
to black) reproduce SDSS steep and flat unabsorbed spectra from Richards et
al. (2006) and Selsing et al. (2016) (red, blue and pink lines). Towards the
infrared, the torus component dominates the continuum.

In the absence of a physical model that accurately reproduces the AGN-plus
galaxy photometry without significant residuals, two broad approaches have
been followed in the literature: (1) The photometry can be degraded (by
increasing the errors) until the model provides a sufficient description. This
approach is implemented in CIGALE v2022.1 (Yang et al., 2022), for example.
(2) A suitable AGN template library can be derived empirically, for instance,
from spectroscopy. Some examples of this approach include the infrared
modelling of Mullaney et al. (2011); Mor & Netzer (2012); Kirkpatrick et al.
(2012), and the recent parametric Seyfert 1 model of Temple et al. (2021).
These templates still have to be constructed to be flexible enough to capture
the full diversity of the AGN phenomenon. While such an empirical AGN model
may not be directly associated with physical properties, it is effective at
“deleting” the AGN light contribution so that the host galaxy light can be
analysed, with realistic uncertainty estimates. The necessary flexibility of
the AGN model leads to many free fitting parameters, in addition to the galaxy
parameters. It is then easily possible that the number of parameters exceeds
the number of data points. Such over-parametrisation leads to degeneracies
which need to be explored computationally. Bayesian inference methods can
address this (see e.g. Calistro Rivera et al., 2016).

We are primarily following the second approach. In Section 2.1, we describe
our flexible parametric AGN model and detail its components. Section 2.2
validates the model’s capability to predict high-quality AGN-galaxy hybrid
spectra. The parameter ranges are also calibrated there. Model parameters are
typeset in monospace, as are code module names in the section headers. For
analysing photometry, GRAHSP builds on some of the open source infrastructure
from CIGALE and accepts similar input and configuration files. However, the
computational methods for achieving rapid model evaluation are substantially
different (Section 2.3). For fitting, GRAHSP is paired with an advanced
inference procedure for uncertainty quantification, which considers model
limitations and AGN variability (Section 2.4).

###  2.1 Model

The spectral model implemented in GRAHSP consists of several components. The
AGN components include a continuum (Section 2.1.1) with emission lines
(Section 2.1.2) and an infrared torus (Section 2.1.3). Bolometric AGN
luminosities are defined in Section 2.1.4. The galaxy components are presented
in Section 2.1.5. Both the AGN and galaxy components can undergo dust
attenuation (Section 2.1.6) and finally redshifting (Section 2.3.1). Figure 1
gives a visual overview of the model components, which are now introduced in
order.

####  2.1.1 The big blue bump continuum (activatepl)

We adopt a flexible model motivated by theoretical and observational
considerations. Thin accretion disk models (Shakura & Sunyaev, 1973) predict a
power-law emission spectrum in the UV to the optical range with a smooth UV
bend-over that depends on the spin, black hole mass, and Eddington ratio (e.g.
Netzer & Trakhtenbrot, 2014). In observed SEDs, the bend towards the far-UV is
wider than predicted (e.g. Blaes et al., 2001), i.e. it is UV-brighter. This
can be partially addressed by assuming strongly inhomogeneous thin disks that
locally fluctuate in temperature (Dexter & Agol, 2011). Because of these
complications, we avoid assuming a specific physical accretion flow model and
instead model the big blue bump (BBB) phenomenologically. From Ly$\alpha$ to
$1$\mathrm{\SIUnitSymbolMicro m}$$, observed stacked spectra (e.g. Zheng et
al., 1997; Selsing et al., 2016) show a power-law spectral density
$F_{\lambda}\propto\lambda^{\alpha}_{\lambda}$ with typical indices
$\alpha_{\lambda}$ between -1.3 and -2.2. There is substantial object-to-
object diversity (Richards et al., 2006) not attributed to dust attenuation.
Additionally, the power-law continuum can change its slope on time scales of
hundreds of days, as demonstrated by optical difference spectra in Ruan et al.
(2014). Above $1\unit{\micro}$, emission by the torus (discussed further in
Section 2.1.3) adds to the SED of the big blue bump. Only thanks to
polarization-spectroscopy by Kishimoto et al. (2008) could this transition be
directly disentangled free of any additional contributions by reprocessing
(emission lines, dusty torus). They showed that the big blue bump emission
continues a power-law decline at least until $2$\mathrm{\SIUnitSymbolMicro
m}$$ with $\alpha_{\lambda}\approx-1.7$ for the observed source.

To encompass this diversity in continua, we adopt the smooth bending power-law
(SBPL) parameterisation of Ryde (1999). This formulation transitions from a
power-law with index $\alpha_{1}$ (uvslope parameter) to one with index
$\alpha_{2}$ (plslope) at a break wavelength $\lambda_{\mathrm{break}}$
(plbendloc). The width of the bend can be controlled with the parameter
$\Lambda$ (plbendwidth). The key L_AGN parameter $\lambda
L_{\mathrm{AGN}}^{5100$\mathrm{\text{Å}}$}$ sets the power-law normalisation
at $\lambda_{0}=5100$\mathrm{\text{Å}}$$. The luminosity spectral density is:

$\displaystyle
L(\lambda)=L_{\mathrm{AGN}}^{5100$\mathrm{\text{Å}}$}\times\left(\frac{\lambda}{\lambda_{0}}\right)^{(\alpha_{1}+\alpha_{2}+2)/2}\times\left(\frac{e^{q}+e^{-q}}{e^{q_{\mathrm{piv}}}+e^{-q_{\mathrm{piv}}}}\right)^{\Lambda\times(\alpha_{2}-\alpha_{1})/2}\times\frac{\lambda_{0}}{\lambda}$
(1)

with $q=\ln\left(\lambda/\lambda_{\mathrm{break}}\right)/\Lambda$ and
$q_{\mathrm{piv}}=\ln\left(\lambda_{0}/\lambda_{\mathrm{break}}\right)/\Lambda$.
We verified that recent state-of-the-art relativistic accretion disk
simulations (e.g., Hagen & Done, 2023) can be approximated with Eq. 1. The
smooth bend-over towards the UV is illustrated by the blue curve in Figure 1.

####  2.1.2 AGN emission lines activatelines

AGN emission lines are superimposed on the continuum. Their contribution to
the observed narrow or broad-band photometry can be substantial (see e.g.
Temple et al., 2021)333Already the same consideration was done in the past for
the contribution of emission lines in the SED fitting of normal galaxies
(e.g., Ilbert et al., 2006; Cardamone et al., 2010; Hsu et al., 2014; Ananna
et al., 2017).. The line luminosity of Hβ is set by default to 2 per cent of
$L_{\mathrm{AGN}}^{5100$\mathrm{\text{Å}}$}$ for the broad line component and
to 0.2 per cent for the narrow line component. This can be modified by the
linestrength_boost_factor parameter. The luminosity ratio of the other lines
relative to Hβ is set based on the broad and narrow line list of Netzer
(1990), listed in Table 1. H$\gamma$ is added based on Rakshit et al. (2020).
With the luminosity ratio defined, the full-width-half-maximum (FWHM;
linewidth parameter) of the lines can be chosen somewhat arbitrarily between
hundreds and tens of thousands $\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$, since
this cannot be distinguished with photometry.

Table 1: List of broad and narrow emission lines, with their centre wavelength, and normalisations relative to H$\beta$. Based on Netzer (1990). These lines are considered in GRAHSP, in addition to FeII. Line |  $\lambda$ [$\mathrm{\text{Å}}$] | broad | narrow
---|---|---|---
Ly$\beta$ | 1026 | 1 | 0.1
OVI | 1035 | 3 | 0.3
Ly$\alpha$ | 1215 | 10 | 35
NV | 1240 | 3 | 3
OI | 1304 | 0.5 | 0
CII | 1336 | 0.3 | 0
SiIV,OIV] | 1400 | 1.5 | 1
NIV] | 1486 | 0.7 | 1
CIV | 1549 | 5 | 8
HeII | 1640 | 0.6 | 0.6
OIII] | 1663 | 0.5 | 0.5
NIII] | 1750 | 0.4 | 0.4
CIII]+SiIII | 1909 | 3 | 3
MgII | 2798 | 3 | 1
[NeV] | 3426 | 0 | 1
[OII] | 3727 | 0 | 3
[NeIII] | 3869 | 0 | 0.5
H$\gamma$ | 4861 | 0.25 | 0.25
HeII | 4686 | 0.1 | 0.1
H$\beta$ | 4861 | 1 | 1
[OIII] | 4959 | 0 | 3
[OIII] | 5007 | 0 | 9
HeI | 5876 | 0.15 | 0.15
[FeVII] | 6087 | 0 | 0.1
[OI] | 6300 | 0 | 0.6
[FeX] | 6374 | 0 | 0.04
H$\alpha$ | 6563 | 4 | 3
[NII] | 6583 | 0 | 1
[SII] | 6716 | 0 | 1
[SII] | 6731 | 0 | 1
CaII | 8498 | 0.1 | 0
CaII | 8622 | 0.1 | 0
[SIII] | 9069 | 0 | 0.2
[SIII] | 9532 | 0 | 0.2
HeI | 10830 | 0.3 | 0.3
P$\alpha$ | 18750 | 0.4 | 0.3

In addition to the individual lines, a FeII forest template is added.
Following Merloni et al. (2010), we select from Bruhweiler & Verner (2008) the
template with density
$n_{\mathrm{H}}=10^{11}\,$\mathrm{/}\mathrm{c}\mathrm{m}^{3}$$,
microturbulence $\xi=20\,$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$$ and ionising
flux $\phi_{\mathrm{H}}=10^{20.5}\,\unit{{}^{-2}{}^{-1}}$, and de-redshift it
from $z=0.004$. The template luminosity is normalised at
$\lambda=4593.4\,$\mathrm{\text{Å}}$$ relative to the Hβ line luminosity (see
above) with ratio $A_{\mathrm{FeII}}$ (AFeII). We also tested inclusion of a
Balmer continuum component, however, the improvements to the presented results
were negligible.

Figure 1 illustrates the importance of these features. The total SED (black)
is elevated in the optical wavelengths above the power-law component (blue) by
the Fe forest and lines (red). Figure 2 compares the AGN model with different
slope (gray shades) to the empirical steep and flat quasar templates by
Richards et al. (2006) (blue and red curves) and the higher redshift stacked
quasar spectrum of Selsing et al. (2016) (pink curve). The key emission
features, shown here with $A_{\mathrm{lines}}=1$, $A_{\mathrm{FeII}}=6$ and
$W_{\mathrm{lines}}=10000\unit{\kilo\per}$, are present in the model and
reproduce the observations. At the shortest wavelengths below $2000\unit{}$,
the current FeII template is potentially missing emission. At wavelengths
above $1\unit{\micro}$, the power law continues downwards (dashed line) as
required by the spectral-polarimetric observations (blue circles), which were
anchored to the dashed line in this figure at the lowest measured wavelengths.
Above $1\unit{\micro}$, the emission begins to be dominated by the emission of
the torus, which is discussed in the next section.

####  2.1.3 AGN torus activategtorus

Figure 3: Overview of the AGN model parameters and how they configure the
spectrum, shown here in $L_{\lambda}$ with arbitrary units. The power law
(blue) is normalised at $5100\unit{}$, where it has a power law slope of
$\beta$. The power law bends over at $\lambda_{\mathrm{bend}}$ towards the UV,
where it has slope $\beta_{UV}$. The width of this transition is set by
$W_{\mathrm{bend}}$. Emission lines with FWHM $W_{\mathrm{lines}}$ and an FeII
template are added and can be further scaled by $A_{\mathrm{lines}}$ and
$A_{\mathrm{FeII}}$, respectively. The torus component (dark yellow) is
normalised by the ratio of $12\mu{}m$ and $5100\unit{}$ $\lambda L_{\lambda}$
luminosities, $f_{\mathrm{cov}}$ (see Eq. 4). It consists of the sum of two
log-quadratic curves, with width $W_{\mathrm{cool}}$ and $W_{\mathrm{hot}}$
and location $\lambda_{\mathrm{cool}}$ and $\lambda_{\mathrm{hot}}$. The peak-
to-peak $\lambda L_{\lambda}$ ratio is set by $f_{\mathrm{hot}}$ (see Eq. 3).
The depth of the Si feature, in emission if positive or in absorption if
negative (here: -1), is set by Si. The flexibility of these 15 parameters is
restricted in Section 2.2.

The AGN infrared emission is associated with the reprocessing of UV emission
by dust (see e.g. the review by Netzer, 2015). The infrared emission is
composed of at least two components, one dominating in the mid-infrared (cold
dust) and one dominating in the near-infrared (hot dust). At long wavelengths,
the infrared emission has recently been imaged in nearby galaxies (e.g.
Garcia-Burillo et al., 2021) and is associated with parsec-scale, cold
molecular dust (the torus). In addition to the cold dust, a near-universal
component is hot dust (e.g. Mor & Netzer, 2012). Through interferometry
observations, the hot dust was associated with polar regions relative to the
torus (Tristram et al., 2014; Asmus et al., 2016). However, Lyu et al. (2017)
found that there is a substantial population (30-40 per cent) of AGN deficient
in warm or hot dust. This suggests that the covering factors, geometries,
and/or dust properties are diverse.

The properties of the dust, chemical composition and grain size distribution
are also not fully understood (e.g. Hoenig, 2013). In empirical studies, the
SED is often approximated by a grey body. Towards short wavelengths, such
parameterisations are flexible, while at the longest wavelength, typically a
$\lambda^{-4}$ power-law decline is assumed (Rayleigh-Jeans tail). The average
shape at the longest wavelength is debated, for example, in Xue et al. (2011)
and Symeonidis (2022). Direct observations with Herschel by Bernhard et al.
(2021) show a wide, log-quadratic bend peaking around $30-50\unit{\micro}$
that is completely dominated by the host galaxy already at $\geq
70\unit{\micro}$.

We use a parameterised empirical template that can emulate this diversity. The
hot dust components can be approximated by black bodies with a temperature
distribution. The peak and low-wavelength distribution is then well described
by a log-quadratic curve. This functional form is sub-optimal at the longest
wavelength, where the Rayleigh-Jeans law motivates a linear rather than
quadratic decline in a log-log plot. However, the tail of the torus hot dust
component is typically dominated by cold galaxy-scale dust emission. As a
result, deviations of the log-quadratic at the long-wavelength regime are
negligible for the modelling (we verify this below in Section 2.2). The same
argument applies to the tail of the hot dust, which is dominated by the cold
dust component. We therefore use the sum of two log-quadratic curves for the
hot and cold dust, which are:

$\begin{split}L_{\mathrm{cool}}&\propto\exp[-(\lambda-\lambda_{\mathrm{COOL}})^{2}/(2W_{\mathrm{COOL}}^{2})],\\\
L_{\mathrm{hot}}&\propto\exp[-(\lambda-\lambda_{\mathrm{HOT}}\,\,\,)^{2}/(2W_{\mathrm{HOT}}^{2})].\end{split}$
(2)

The parameters include the peak $\lambda_{\mathrm{COOL}}$ (COOLlam) and width
$W_{\mathrm{COOL}}$ (in dex; COOLwidth) of the Gaussian-like form for the
spectrum (and analogously for the hot dust component). These parameters are
related to the temperature distribution of the dust. The normalisation of the
hot dust peak relative to the cold dust peak in $\lambda L_{\lambda}$ is a
parameter that we term the hot factor:

$f_{\mathrm{hot}}=\frac{\lambda_{\mathrm{HOT}}L_{\mathrm{HOT}}}{\lambda_{\mathrm{COOL}}}.$
(3)

With the two log-quadratic curves fixed, the normalisation of the torus to the
AGN powerlaw has to be defined. The ratio of the near to mid-infrared
continuum torus template amplitude to that of the UV to optical power-law
template, $f_{\mathrm{cov}}$ (fcov), is commonly referred to as the torus
covering factor:

$\frac{(\lambda L_{\lambda})(12$\mathrm{\SIUnitSymbolMicro m}$)}{(\lambda
L_{\lambda})(5100$\mathrm{\text{Å}}$)}=2.5\times f_{\mathrm{cov}}.$ (4)

However, the translation to a geometric covering factor of a physical
structure is complicated by anisotropy in the emission profile and
distribution of dust clumps (Stalevski et al., 2016). Section 2.2 below
demonstrates that the model can reproduce torus spectral templates derived by
other works, including those of Lyu et al. (2017), Mor & Netzer (2012),
Mullaney et al. (2011) and Kirkpatrick et al. (2012), as well as individual
AGN in the local Universe.

An important feature in the mid-infrared is the absorption and emission by Si
dust at $12$\mathrm{\SIUnitSymbolMicro m}$$. Because the Si feature is not as
pronounced as smooth models suggest, this has been interpreted as evidence for
a clumpy torus (Nenkova et al., 2008b). However, it was later shown that
similar effects can be produced with smooth geometries (Feltre et al., 2012).
Goulding et al. (2012) showed that deep Si absorption features are associated
with edge-on galaxies, suggesting that non-nuclear dust imprints this feature.

We allow the model to place Si in absorption or emission. An Si template is
created by taking the difference between the average templates of faint (on
average in absorption) and luminous (on average in emission) AGN in the 8 to
18$\mathrm{\SIUnitSymbolMicro m}$ range from Mullaney et al. (2011). The
template is normalised at $12$\mathrm{\SIUnitSymbolMicro m}$$. The Si
parameter then controls the amplitude of this contribution, with zero
corresponding to it not being present. The effect of this parameter is very
localised, as shown by comparing the dotted and solid dark yellow curves in
Fig. 3. Indeed, this parameter is inconsequential if no photometry filter
covers the 8 to $18$\mathrm{\SIUnitSymbolMicro m}$$ rest-frame wavelengths.

####  2.1.4 Bolometric luminosities

The luminosity integrated over the entire wavelength range is a crucial
measure of the radiative energy budget in the circum-nuclear environments of
supermassive black holes. However, since some radiation is absorbed and then
re-emitted (as emission lines or in the torus), some care is needed to define
what to compute to avoid double-counting. Furthermore, because our conversion
between fluxes and luminosities assumes the luminosity distance, we compute
isotropic luminosities, i.e., under the assumption of an isotropic radiation
profile. Converting to more realistic, anisotropic emission luminosities
requires a physical model (see, e.g., Stalevski et al., 2016, for a detailed
discussion).

GRAHSP computes two bolometric AGN luminosities. Both are intrinsic, i.e.,
before applying attenuation, rest-frame, isotropic luminosities. The first,
lumBolBBB, integrates all AGN SED components except for the torus upwards of
91.2nm. The substantial contribution of the ionising far-UV and X-rays are not
included in lumBolBBB. This is because the wavelength range is rarely directly
measured. The user should apply a model-dependent correction factor to obtain
the true bolometric luminosity. The second is the bolometric luminosity of the
torus, lumBolTOR, integrated over the entire wavelength range. While lumBolBBB
is primarily informed by data in UV to optical rest-frame wavelengths,
lumBolTOR is informed by infrared data. Both bolometric luminosities are
reported, as well as the ratio between them, ratioTORBBB=lumBolTOR/lumBolBBB.
How to interpret such a ratio of instantaneous, isotropic luminosities as a
covering factor of a light-year-sized torus that reprocesses the radiation
from an anisotropically emitting accretion disk is discussed extensively in
Stalevski et al. (2016).

####  2.1.5 Galaxy model

For the host stellar population, we adopt standard CIGALE modules Boquien et
al. (2019); Yang et al. (2022). These implement stellar population synthesis
(SPS) by Bruzual & Charlot (2003) (bc03) and Maraston (2005) (m2005) combined
with a parametric star formation history (SFH) model. The stellar population
is scaled to a total stellar mass of $M_{\star}$, a key parameter setting the
normalisation of the galaxy emission spectrum. In the SPS, a metallicity needs
to be assumed. The effects of metallicity on the resulting SED are well known
to be degenerate with stellar age. Instead of exploring this degeneracy, we
follow previous work and choose to measure the effective stellar age, assuming
a fixed metallicity, and interpret the results in this context.

The galaxy evolution field has recently developed non-parametric SFHs (e.g.
Iyer et al., 2019), which can approximate complex SFHs seen in cosmological
simulations. Leja et al. (2019) demonstrated that non-parametric SFHs can more
faithfully reconstruct the stellar mass build-up even in individual galaxies,
given precision photometry. However, both parametric and non-parametric
approaches agree within 0.1 dex on present-time properties, such as the
current stellar mass and the average star formation rate (SFR) in the last 100
Myr (Leja et al., 2019). In this work (see also Ciesla et al., 2015), we
demonstrate that contamination by AGN light creates much larger uncertainties
when inferring SFR. To compare results on stellar mass from the literature, it
is beneficial to stick with a simple parametric SFH. In GRAHSP, any SFH models
implemented in CIGALE (see Boquien et al., 2019) can be enabled. For this
work, we adopt a tau-delayed SFH (sfhdelayed), where the SFH rises linearly
with time $t$ and is then truncated with an exponential cut-off time scale
$\tau$:

$\mathrm{SFR}(t)\propto\frac{t}{\tau^{2}}\exp(-t/\tau)$ (5)

Figure 4 illustrates galaxy spectra of different SFHs (inset) for
completeness. A maximum age for the oldest stars, $t_{0}$, can be set
(main_age), which defines the start of the SFH.

We point out that this parameterisation induces a characteristic SFR prior.
The current SFR (averaged over the last 10 Myr for example) can be read off as
the right-most point of each curve in the inset of Fig. 4. The maximum
(highest $\tau$) values are near SFR(10Myr)=$M_{\star}$/10Myr (and similar for
SFR defined over other time windows), where most stellar mass was built
relatively recently. A uniform or log-uniform grid on $\tau$ then induces an
SFR prior that declines exponentially to arbitrarily low SFR (lowest $\tau$).
The endpoints of the inset of Fig. 4 indeed show a bimodality, with few points
in the middle. This means that even if a random $\tau$ and $M_{\star}$ were
picked, a main sequence can emerge. A red cloud can emerge by imposing a SFR
floor.

As for the AGN component, emission lines may contribute a substantial flux in
galaxy spectra (Ilbert et al., 2006; Hsu et al., 2014). Nebular emission lines
are added with the nebular CIGALE module (Boquien et al., 2013, 2019). Their
effect is clearly seen in recently star-forming templates in Fig. 4.

Figure 4:  Example galaxy SEDs of stellar populations with different star
formation histories. The star formation rate (see inset) rises linearly
towards the present (right), with an exponential cutoff timescale $\tau$. At
low $\tau$, the yellow SFR curve truncates quickly, i.e., is dominated by old
stars. The corresponding yellow SED in the main panel peaks between
$0.3-3\unit{\micro}$. The blue curve in the inset corresponds to continuously
rising star formation. The corresponding blue SED in the main panel is
dominated by luminous young stars, nebular emission lines and infrared dust
emission. Minimal attenuation, E(B-V)=0.01, is applied.

For gauging the relative importance of AGN and host galaxy emission, GRAHSP
also computes an AGN fraction. Following Dale et al. (2014), fracAGNDale is
the AGN luminosity fraction from $5-20\unit{\micro}$. Additionally, the
bolometric AGN fraction fracAGNTOR is the ratio of lumBolTOR to the bolometric
galaxy luminosity. These are also computed before applying attenuation (next
section, Section 2.1.6), and for consistency, the bolometric galaxy luminosity
does not include galaxy dust re-emission. From our key model parameters,
$L_{\mathrm{AGN}}^{5100$\mathrm{\text{Å}}$}$ and $M_{\star}$, we can also
consider a rough contrast ratio of AGN to host galaxy, after converting the
stellar mass with a mass-to-light ratio $\Upsilon$:

$\lambda=L_{\mathrm{AGN}}^{5100$\mathrm{\text{Å}}$}/M_{\star}\times\Upsilon$
(6)

For simplicity, for $\Upsilon$ we adopt the solar bolometric mass-to-light
ratio $\Upsilon_{0}=M_{\odot}/L_{\odot}=M_{\odot}/\qty{3.83e33}{\per}$, which
is close to the Milky Way V-band mass-to-light ratio
$\Upsilon_{\mathrm{V,MW}}=1.5$ (e.g., Flynn et al., 2006). If we assume the
same bolometric corrections for X-ray and $5100\unit{}$ (justified in Section
5.4), from black hole mass scaling relations $\lambda=2.7$ is the approximate
Eddington limit (Aird et al., 2017).

####  2.1.6 Attenuation model (biattenuation)

Figure 5:  Effect of attenuation on the galaxy model. Models are shown from
intrinsic (dark blue) to strongly attenuated (dark red). For illustration, the
extremely attenuated local low-metallicity star-bursting galaxy Haro 11 from
Lyu et al. (2016) is overplotted as a dashed red curve.  Figure 6:  Effect of
attenuation on the AGN model. Models are shown from intrinsic (blue top
curves) to strongly attenuated (red bottom curves). The intrinsic torus model
component (see Section 2.1.3) is shown in dashed black. The dark red curves
illustrate that the most heavily attenuation can also suppress the torus.

The spectra of AGN and galaxies are frequently attenuated by dust along the
line of sight. Dust attenuation can be modeled by a variety of empirical laws.
Salvato et al. (2009) and Hopkins et al. (2004), analysing large photometric
and spectroscopic AGN samples, respectively, find a preference for the dust
attenuation law of Prevot et al. (1984) derived from the Small Magellanic
Cloud (SMC), $A_{\text{SMC}}(\lambda)$. There is still debate whether steeper
and even entirely feature-less (power-law like) attenuation may be preferable
(Fynbo et al., 2013; Zafar et al., 2015). Because the reddening of the
continuum is dependent on the assumed intrinsic AGN continuum model, for which
we chose a flexible parameterisation, we remain with a SMC-like dust
attenuation model. The level of attenuation is parameterised with E(B-V):

$L_{\lambda}^{\prime}(\lambda)=L_{\lambda}(\lambda)\times 10^{E(B-V)\times
A_{\text{SMC}}(\lambda)/-2.5}.$

We approximate the SMC attenuation curve as a broken power law:

$A_{\text{SMC}}(\lambda)=N\times\left(\frac{\lambda}{\lambda_{\mathrm{break}}}\right)^{\gamma}$

where $\gamma=\gamma_{\mathrm{OPT}}$ (default: -1.2) below
$\lambda_{\mathrm{break}}$ and $\gamma=\gamma_{\mathrm{NIR}}$ (default: -3)
above. The normalisation is $N=1.2$ at
$\lambda_{\mathrm{break}}=1100\unit{\nano}$.

In our implementation, the AGN and galaxy components are attenuated
differently. The galaxy attenuation level is parameterized by the E(B-V) color
excess parameter, and the total galaxy luminosity absorbed is recorded.
Enforcing energy balance, the luminosity is then re-emitted in the infrared
following the Dale et al. (2014) model of galactic dust emission
(galdale2014). Figure 5 illustrates this, and how the model approximates a
dusty star-burst galaxy. For the AGN light, the situation is different (see
also Calistro Rivera et al., 2016). The dust emission is already empirically
modeled with the torus component (Section 2.1.3 above). The energy may not be
balanced due to variability and anisotropy, as the line-of-sight absorption
differs from geometrically averaged absorption. Since the AGN is embedded in
the host galaxy dust but also nuclear dust, it may undergo additional
attenuation. Therefore, we attenuate the AGN components not only with E(B-V),
but with E(B-V) \+ E(B-V)-AGN, where E(B-V)-AGN is a parameter giving the
nuclear attenuation color excess.

The E(B-V)-AGN parameter also allows the AGN model to transition from a Sy1 to
a Sy2. Figure 6 illustrates how the blue continuum, its lines and ultimately
also the torus are attenuated by increasing values of E(B-V)-AGN. In our
model, narrow and broad emission lines are both attenuated to the same extent,
while in reality, narrow emission lines should remain visible in Sy2. However,
this approximation is sufficient because in the relevant wavelength range, the
photometry of Sy2 galaxies is dominated by host galaxy light (Hickox et al.,
2017).

###  2.2 Model validation and parameter range calibration

Figure 7:  Atlas of AGN spectra (red) from Brown et al. (2019). The black
curve shows our best-fit GRAHSP model. Dashed curves show model components.

Table 2: Best fit model parameters for the AGN SED Atlas of Brown et al. (2019). The first row indicates the parameter; the second row shows whether the logarithm of the parameter is listed, the third row indicates the units. The last few rows give the minimum, maximum and 10% and 90% quantiles of each column. | $M_{\star}$ | $\tau$ | age | $\alpha$ |  LAGN |  AFeII |  Alines |  Wline | $\beta$ | $\lambda_{\mathrm{bend}}$ |  Wbend | $\beta_{\mathrm{UV}}$ |  EB-V |  E${}_{\mathrm{B-V}}^{\mathrm{AGN}}$ | Si | $f_{\mathrm{cov}}$ | $\lambda_{\mathrm{cool}}$ |  Wcool | $\lambda_{\mathrm{hot}}$ |  Whot | $f_{\mathrm{hot}}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
| log | log | log |  | log | log | log | log |  | log | log |  | log | log |  |  |  | log |  | log |
Galaxy | $M_{\odot}$ | Myr | Myr |  | erg/s |  |  | km/s |  | nm | dex |  | mag | mag |  |  |  $\mu$m | dex |  $\mu$m | dex |
IRAS-11119+3257 | 10.1 | 2.9 | 3.4 | 1.5 | 45.8 | 0.7 | -0.3 | 3.7 | -1.6 | 2.1 | -1.7 | -0.5 | -0.9 | 0.1 | -2.2 | 0.2 | 30.0 | 0.6 | 3.8 | 1.0 | 9.9
IRAS-F16156+0146 | 8.1 | 1.0 | 1.4 | 1.5 | 44.5 | 0.2 | 0.5 | 3.4 | -2.2 | 1.6 | -2.0 | 1.3 | -0.7 | -0.1 | -2.9 | 0.8 | 26.4 | 0.9 | 6.4 | 0.2 | 0.3
F2M1113+1244 | 9.9 | 3.1 | 2.7 | 2.6 | 46.0 | 1.3 | -0.5 | 2.8 | -1.6 | 2.3 | -2.0 | -0.1 | -0.9 | -0.1 | -0.8 | 0.8 | 30.0 | 0.8 | 1.7 | 0.6 | 4.0
PKS-1345+12 | 11.0 | 2.0 | 3.0 | 2.2 | 44.6 | 1.8 | -0.3 | 3.9 | -1.4 | 1.9 | -0.3 | -0.7 | -0.8 | 0.9 | -0.5 | 0.9 | 28.1 | 0.6 | 2.6 | 0.6 | 0.0
MRK231 | 10.6 | 2.5 | 3.3 | 2.7 | 45.1 | 0.6 | -0.4 | 3.9 | -1.7 | 2.2 | -0.7 | 1.7 | -1.1 | -0.1 | -1.6 | 0.7 | 30.0 | 0.2 | 5.1 | 1.0 | 2.5
2MASXJ13000533+1632151 | 7.5 | 1.7 | 1.4 | 1.8 | 45.4 | 1.1 | -0.5 | 4.1 | -2.5 | 2.3 | -0.2 | -1.0 | -0.8 | 0.1 | -0.4 | 0.1 | 18.8 | 0.3 | 4.4 | 0.5 | 2.5
NGC3227 | 7.0 | 3.5 | 3.0 | 1.0 | 42.8 | 0.9 | -0.3 | 3.6 | -1.9 | 2.3 | -1.2 | 1.9 | -0.8 | -0.6 | -0.8 | 0.2 | 19.1 | 0.4 | 4.9 | 0.4 | 0.2
NGC4051 | 8.3 | 1.2 | 3.7 | 1.3 | 41.7 | -0.5 | -0.5 | 3.4 | -1.1 | 2.3 | -2.0 | -1.1 | -1.6 | -1.2 | 0.0 | 0.7 | 12.9 | 0.5 | 2.7 | 0.4 | 0.2
NGC5548 | 7.8 | 1.3 | 3.3 | 2.2 | 43.1 | 0.8 | -0.4 | 3.8 | -1.2 | 1.7 | -1.8 | 0.8 | -2.0 | -1.6 | -0.1 | 0.9 | 12.4 | 0.5 | 1.3 | 0.4 | 0.4
NGC7469 | 8.5 | 3.7 | 3.1 | 1.9 | 43.7 | 1.2 | -0.4 | 4.1 | -1.4 | 1.7 | -1.2 | 1.1 | -2.1 | -1.4 | -0.5 | 1.0 | 24.3 | 0.9 | 1.7 | 0.3 | 0.1
MRK590 | 9.1 | 1.5 | 4.1 | 1.8 | 42.9 | 0.7 | -0.2 | 4.1 | -1.0 | 1.7 | -1.2 | 0.9 | -1.3 | -2.1 | 0.1 | 1.0 | 15.3 | 0.4 | 0.9 | 0.4 | 1.6
MRK1502 | 8.2 | 2.3 | 4.1 | 1.7 | 44.5 | 1.3 | -0.3 | 3.9 | -1.7 | 2.1 | -0.9 | -0.2 | -2.0 | -1.9 | 0.7 | 1.0 | 13.0 | 0.5 | 2.3 | 0.4 | 0.6
MRK817 | 7.7 | 2.1 | 3.9 | 1.9 | 44.0 | 0.9 | -0.4 | 3.8 | -1.5 | 1.9 | -1.7 | -0.0 | -2.0 | -2.0 | 0.3 | 0.6 | 17.9 | 0.5 | 2.0 | 0.4 | 0.3
ARK564 | 9.0 | 1.6 | 3.2 | 1.3 | 43.6 | 1.4 | -0.5 | 3.5 | -1.4 | 1.9 | -1.7 | 0.3 | -1.9 | -1.2 | 0.1 | 0.4 | 21.5 | 0.3 | 3.9 | 0.6 | 0.6
MRK290 | 8.2 | 2.6 | 2.7 | 2.2 | 43.4 | 1.3 | -0.3 | 3.7 | -1.6 | 2.0 | -0.6 | 2.0 | -1.9 | -2.0 | -0.5 | 0.6 | 28.4 | 0.9 | 0.9 | 1.0 | 5.3
NGC3516 | 9.3 | 1.6 | 3.4 | 1.9 | 42.9 | 1.9 | -0.5 | 3.8 | -1.3 | 2.2 | -1.9 | 1.3 | -1.9 | -1.7 | 0.1 | 0.6 | 14.0 | 0.6 | 1.3 | 0.6 | 0.4
PG1415+451 | 8.9 | 2.5 | 4.1 | 1.9 | 44.4 | 0.5 | -0.5 | 3.5 | -1.6 | 2.0 | -1.2 | 0.1 | -2.1 | -1.5 | 0.7 | 0.5 | 14.8 | 0.4 | 1.1 | 0.6 | 0.7
3C390-3 | 7.4 | 3.2 | 1.3 | 1.0 | 44.1 | -0.9 | 0.1 | 3.9 | -2.3 | 2.0 | -0.8 | -0.6 | -1.1 | -1.9 | -0.6 | 0.7 | 23.2 | 0.2 | 1.0 | 0.9 | 3.6
NGC4151 | 7.8 | 1.6 | 2.9 | 3.0 | 42.9 | 1.3 | -0.2 | 3.7 | -2.2 | 1.8 | -0.9 | -0.4 | -1.1 | -2.1 | 0.1 | 0.5 | 15.0 | 0.3 | 2.6 | 0.4 | 0.5
FAIRALL9 | 9.2 | 2.1 | 3.7 | 2.3 | 44.3 | 1.4 | -0.2 | 3.8 | -2.0 | 2.0 | -0.5 | 0.8 | -2.1 | -1.6 | 0.1 | 0.7 | 14.3 | 0.4 | 2.0 | 0.5 | 1.3
MRK493 | 8.3 | 2.1 | 4.0 | 2.1 | 43.5 | 1.6 | -0.5 | 3.6 | -1.8 | 2.0 | -1.5 | 0.5 | -2.1 | -1.8 | 0.3 | 0.5 | 18.4 | 0.6 | 1.1 | 0.6 | 0.8
MRK926 | 9.7 | 3.0 | 3.4 | 2.0 | 44.1 | -0.1 | -0.2 | 3.6 | -1.6 | 2.1 | -0.9 | 0.0 | -2.1 | -1.3 | 0.0 | 0.4 | 10.0 | 0.6 | 2.0 | 0.4 | 2.3
PG1211+143 | 8.3 | 1.3 | 1.9 | 1.2 | 44.8 | -0.4 | -0.3 | 3.5 | -2.2 | 2.0 | -1.3 | 1.5 | -1.9 | -1.1 | -0.1 | 0.4 | 19.3 | 0.2 | 3.4 | 0.6 | 5.5
MRK876 | 8.7 | 2.7 | 4.1 | 1.9 | 45.0 | -0.6 | -0.3 | 3.8 | -1.9 | 2.1 | -1.8 | -0.6 | -2.0 | -1.5 | 0.2 | 0.3 | 14.1 | 0.6 | 2.3 | 0.4 | 1.0
PG2349-014 | 8.0 | 2.8 | 4.1 | 1.9 | 44.9 | 1.4 | -0.3 | 3.9 | -1.8 | 2.1 | -0.7 | -0.3 | -1.9 | -2.0 | -0.2 | 0.4 | 14.5 | 0.4 | 1.9 | 0.5 | 1.3
H1821+643 | 8.5 | 2.7 | 3.2 | 1.8 | 46.0 | 1.2 | -0.1 | 3.8 | -2.7 | 2.2 | -0.2 | 0.2 | -2.0 | -1.9 | 0.1 | 0.5 | 16.6 | 0.3 | 1.4 | 0.6 | 0.7
NGC5728 | 10.4 | 2.0 | 3.7 | 2.6 | 42.6 | 0.0 | -0.5 | 4.1 | -1.3 | 1.7 | -2.0 | -0.1 | -2.1 | -1.2 | -1.6 | 1.0 | 12.2 | 0.8 | 1.2 | 0.2 | 1.5
MRK509 | 7.5 | 1.7 | 1.5 | 2.6 | 44.3 | 1.3 | -0.2 | 3.7 | -1.9 | 2.0 | -0.5 | 1.4 | -2.1 | -1.8 | -0.7 | 0.5 | 21.3 | 0.4 | 1.1 | 0.9 | 3.3
3C120 | 7.9 | 1.8 | 3.3 | 1.0 | 44.1 | 1.2 | -0.1 | 3.6 | -2.5 | 2.0 | -0.3 | -0.6 | -2.1 | -1.9 | 0.1 | 0.6 | 12.1 | 0.4 | 2.0 | 0.4 | 0.8
MRK279 | 8.7 | 1.9 | 3.6 | 2.2 | 44.0 | 0.6 | -0.3 | 3.8 | -1.6 | 2.1 | -0.9 | 0.1 | -2.1 | -1.8 | -0.1 | 0.3 | 13.8 | 0.4 | 3.0 | 0.4 | 0.6
3C351 | 8.1 | 2.8 | 3.8 | 2.2 | 45.9 | 1.0 | -0.5 | 4.1 | -1.8 | 2.0 | -1.9 | 0.3 | -1.9 | -1.5 | 0.6 | 0.3 | 14.1 | 0.4 | 2.4 | 0.3 | 1.0
ARK120 | 8.8 | 2.3 | 4.0 | 2.4 | 44.3 | 1.3 | -0.1 | 3.8 | -2.0 | 1.9 | -0.7 | 0.6 | -1.3 | -2.1 | 0.2 | 0.3 | 11.1 | 0.4 | 1.6 | 0.5 | 2.1
PG1307+085 | 10.3 | 1.7 | 3.3 | 2.3 | 44.8 | 1.5 | -0.3 | 3.6 | -2.2 | 2.0 | -0.1 | 0.3 | -1.8 | -2.0 | -0.1 | 0.4 | 16.1 | 0.2 | 1.6 | 0.7 | 1.3
PG0052+251 | 9.7 | 1.1 | 2.5 | 1.7 | 44.9 | 1.0 | -0.1 | 3.8 | -2.7 | 2.3 | 0.1 | -1.0 | -1.9 | -2.1 | -0.3 | 0.3 | 28.6 | 0.9 | 1.2 | 0.8 | 5.1
MRK110 | 10.1 | 1.6 | 4.1 | 1.4 | 43.8 | 1.3 | -0.1 | 3.3 | -2.4 | 2.1 | -0.8 | -0.1 | -1.9 | -1.4 | 0.2 | 0.3 | 14.0 | 0.3 | 1.9 | 0.5 | 1.3
3C273 | 7.0 | 2.2 | 2.7 | 2.0 | 46.0 | 1.1 | -0.4 | 3.5 | -2.6 | 2.1 | -0.2 | -0.1 | -2.1 | -1.6 | 0.0 | 0.2 | 17.4 | 0.4 | 1.7 | 0.5 | 1.5
TON951 | 9.7 | 1.4 | 4.1 | 2.0 | 44.6 | -0.2 | -0.5 | 3.9 | -1.6 | 2.2 | -0.9 | -1.2 | -2.1 | -1.7 | 0.1 | 0.1 | 12.5 | 1.0 | 7.2 | 0.5 | 6.0
PG0026+129 | 7.8 | 2.3 | 3.5 | 2.8 | 45.2 | 0.5 | -0.5 | 3.8 | -2.0 | 2.0 | -0.4 | -2.0 | -1.5 | -1.7 | 1.1 | 0.1 | 13.4 | 0.3 | 3.0 | 0.3 | 1.3
min | 7.0 | 1.0 | 1.3 | 1.0 | 41.7 | -0.9 | -0.5 | 2.8 | -2.7 | 1.6 | -2.0 | -2.0 | -2.1 | -2.1 | -2.9 | 0.1 | 10.0 | 0.2 | 0.9 | 0.2 | 0.0
10% | 7.5 | 1.3 | 2.0 | 1.2 | 42.9 | -0.2 | -0.5 | 3.5 | -2.6 | 1.7 | -1.9 | -0.9 | -2.1 | -2.0 | -0.8 | 0.2 | 12.4 | 0.3 | 1.1 | 0.3 | 0.3
90% | 10.2 | 3.0 | 4.1 | 2.7 | 45.8 | 1.6 | -0.1 | 4.1 | -1.3 | 2.3 | -0.2 | 1.5 | -0.9 | -0.1 | 0.3 | 0.9 | 28.4 | 0.9 | 4.3 | 0.9 | 5.3
max | 11.0 | 3.7 | 4.1 | 3.0 | 46.0 | 1.9 | 0.5 | 4.1 | -1.0 | 2.3 | 0.2 | 2.0 | -0.7 | 0.9 | 1.1 | 1.0 | 30.0 | 1.0 | 7.2 | 1.0 | 9.9

Table 3: Best-fit model parameters for the AGN SED Atlas of Brown et al. (2019). Only AGN shape parameters are shown. The first row represents an average AGN, while the next two rows represent warm- and hot-dust-deficient sub-populations. These differ in the torus properties listed in the last six columns, which give the covering factors and log-quadratic shape parameters of the cool and hot components. The last row is an AGN torus template starting at 6$\mu{}m$, therefore the first columns are blanked out. | AFeII | Alines | Wline | $\beta$ | $\lambda_{\mathrm{bend}}$ | Wbend | $\beta_{\mathrm{UV}}$ | EB-V | E${}_{\mathrm{B-V}}^{\mathrm{AGN}}$ | Si | $f_{\mathrm{cov}}$ | $\lambda_{\mathrm{cool}}$ | Wcool | $\lambda_{\mathrm{hot}}$ | Whot | $f_{\mathrm{hot}}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
| log | log | log |  | log | log |  | log | log |  |  |  | log |  | log |
Template |  |  | km/s |  | nm | dex |  | mag | mag |  |  |  $\mu$m | dex |  $\mu$m | dex |
Elvis-1994 | 1.8 | -0.5 | 3.6 | -2.7 | 1.7 | 0.2 | 1.9 | -2.0 | -2.0 | 0.1 | 0.6 | 14.0 | 0.4 | 2.5 | 0.5 | 1.1
Lyu-2017-WDD | 0.6 | -0.5 | 3.5 | -2.3 | 2.0 | -1.4 | 0.3 | -2.0 | -1.2 | -0.4 | 0.2 | 18.8 | 0.3 | 2.4 | 0.5 | 3.1
Lyu-2017-HDD | 1.9 | -0.5 | 3.6 | -2.6 | 2.1 | -1.2 | -0.5 | -1.2 | -1.7 | -0.5 | 0.2 | 16.0 | 0.5 | 1.8 | 0.6 | 6.3
Mullaney-2011 | – | – | – | – | – | – | – | – | -1.8 | -0.7 | 0.6 | 14.7 | 0.4 | 2.5 | 0.5 | 0.3
Figure 8:  Relative residuals of the fits of Fig. 7 are shown as curves. Four
individual cases with the largest deviations are highlighted in colours and
labelled. The distribution of residuals is presented on the right-hand side as
a histogram. The vast majority of residuals (96.6%) are concentrated to
smaller than 20 per cent deviations.

This section tests and validates the performance of the GRAHSP model.
Additionally, for each model parameter we identify the range required to
reproduce the observations, to be used as prior knowledge when fitting lower-
quality and noisier photometric observations typical of distant AGN samples.
For this purpose, consistent, high-quality UV to far-IR spectra of galaxies
hosting AGN are needed. Brown et al. (2019) compiled spectroscopic and
photometric data from 0.09 to 30$\mu m$ on a diverse set of 41 local AGN. For
each AGN, they carefully cross-calibrated and aperture-corrected the
observations into one continuous broad-band SED. Figure 7 shows these data-
driven model spectra in red. This AGN SED atlas was created to demonstrate the
diversity of galaxies with AGN, including different host-AGN contrasts and
infrared torus shapes. We use the AGN SED Atlas to identify reasonable ranges
for each fitting parameter to restrict our otherwise very flexible model. We
also include the empirical AGN templates of Elvis et al. (1994), Mullaney et
al. (2011) and the hot-dust deficient and warm-dust deficient templates of Lyu
et al. (2017) derived from 87 Palomar-Green quasars.

We optimize our model parameters to best approximate each spectrum. The best-
fit model is overlaid in black in Figure 7. Individual model components are
presented as dashed curves. Overall, the fit is very good, and the model
captures the diversity of AGN-galaxy SED shapes. The model quality can be
quantified by looking at the fit residuals in Figure 8. These were computed by
averaging the model and data in windows of $\Delta\lambda=0.2\lambda$, as
relevant to broad-band photometry. The residuals remain below 20 per cent
across the entire wavelength range considered, with few (3.4%) exceptions. In
2MASXJ13000533+1632151, the reddened galaxy continuum is not perfectly
approximated near 300nm. IRAS-F16156+0146 shows an extremely deep Si
absorption feature, which is not modelled well in the wings. NGC 7469 and NGC
5728 show complex galaxy polycyclic aromatic hydrocarbon emission features
which are not fully captured by the Dale et al. (2014) templates adopted here.

The best-fit model parameters for each AGN are listed in Table 2. The
diversity of the AGN Atlas can only be captured by allowing all listed
parameters to vary. However, we can use the typical values to restrict the
allowed parameter range. The last rows in Table 2 list the 10 and 90 per cent
quantile for each parameter. For example, for the attenuation of the AGN,
E(B-V) values ranging from 0.01 to 1 mag are found, and for the galaxy, E(B-V)
values ranging from 0.01 to 0.1 mag are found. The best-fit model parameters
for four empirical templates are listed in Table 3. Here we see that the
normal AGN template of Elvis et al. (1994) has a covering factor of 60 per
cent with cool and hot torus components centred at 14 and 2.5$\mu{}m$, of
widths 0.4 and 0.5 dex, respectively. This is very similar to the Mullaney et
al. (2011) template in the last row. In contrast, the warm-dust deficient
template of Lyu et al. (2017) is decomposed with a narrower cold dust energy
distribution, while the hot-dust deficient template contains a hot dust
component centred at shorter wavelengths. The amplitude factor of the hot dust
component (last column) also varies from 0.3 to 6. We cannot rule out the
possibility that parameter degeneracies influence the values listed in Table
3, however, the found ranges from Table 2 and Table 3 motivate the parameter
ranges when fitting the model to data. To achieve efficient fits with the
remaining 21 free model parameters requires an advanced inference engine.

###  2.3 Inference engine

####  2.3.1 Photometric flux prediction

The composed model is redshifted and converted into observable flux
$F_{\lambda}$ using the luminosity distance $D_{L}$:
$F_{\lambda}(\lambda)=L_{\lambda}(\lambda)/(4\pi D_{L}^{2})/(1+z)$. For these
steps, the implementation is taken from CIGALE (Boquien et al., 2013,
redshifting module), which also applies absorption by the inter-galactic
medium (Madau, 1992). Finally, the flux of a photometric filter band $i$ of
interest is predicted with the filter transmission curve $T(\lambda)$ as
$F_{i}=\lambda_{\mathrm{pivot}}^{2}\int{F(\lambda)T(\lambda)\,d\lambda}$ with
the pivot wavelength
$\lambda_{\mathrm{pivot}}=\left(\int{T(\lambda)\times\lambda
d\lambda}\right)/\left(\int{T(\lambda)/\lambda\,d\lambda}\right)$. See
appendix A2 of Bessell & Murphy (2012) for a discussion why the pivot
wavelength, rather than the mean filter wavelength, is relevant here. The
model filter flux $F_{i}$ is a function of all model parameters,
$F_{i}=f(L_{\mathrm{AGN}},M_{\star},\dotsc)$.

####  2.3.2 Computational optimisations

Figure 9: GRAHSP computation pipeline. One sub-pipeline (yellow) produces the
galaxy SED (large yellow box), one sub-pipeline (red large box) produces the
AGN SED. Each sub-pipeline has a separate cache. Subsequent modules are
applied on-the-fly, including the AGN power-law, attenuating the AGN and the
galaxy components, summing and transforming into observed-frame fluxes.

GRAHSP includes several optimisations that allow for the fitting of a single
galaxy in 2-3 minutes on a typical single computing core. This includes the
generation of visualisations of the SED fit and the inferred parameter
uncertainties. This speed enables processing large samples, or experimentation
and variation of models and data processing of individual objects. For rapid
evaluation of a model with a substantial number of components, and fitting
that fully explores the degeneracies between the many model parameters, GRAHSP
departs substantially from existing approaches. We review some of these first.

In CIGALE, the model is implemented as a pipeline of modules that successively
build up the SED. Each module has parameters. Each parameter can take a set of
possible values, which is specified by the user as a list of numbers. CIGALE
then creates an n-dimensional grid of all possible parameter combinations and
instantiates each model spectrum. Because of the curse of dimensionality and
limited computing resources (in memory and computation budget), this requires
the user to sparsely choose the grid points for each parameter. This problem
is exacerbated by the fact that an AGN+galaxy model is a linear combination of
two models. However, CIGALE is unaware of this and thus builds the full model
grid inefficiently, with the frac_agn specifying the relative normalisation of
the AGN and galaxy components. The AGNFitter (Calistro Rivera et al., 2016)
and FortesFit (Rosario, 2019) instead first build two grids and linearly
combine them during fitting, with stellar mass and AGN luminosity as
normalisations. The parameterisation can be important because a log-uniform
prior on stellar mass and AGN luminosity, each, prefers a much wider specific
accretion rate distribution than setting a uniform prior on frac_agn. The
latter can induce a bias to a narrow range of specific accretion rates, which
may affect galaxy-AGN co-evolution studies. To test a new model variant (e.g.,
for a new filter set and/or redshift) with the two-grid precomputation
approach, the user needs to write custom code to create the grids, and to
balance the possible science goals with the grid size and computational cost.
Another approach is to generate models on-the-fly, including instantiating the
stellar population, attenuation, redshifting, and filter flux computation.
Crucially, this enables full exploration of degenerate parameters, such as the
E(B-V) and power-law slopes, with arbitrarily fine resolution.

GRAHSP takes a hybrid approach. It uses the pipeline-of-modules approach of
CIGALE, but computes models on-the-fly. Models are cached using a hash map,
which retrieves the output of the last module where all parameters matched.
Therefore, while sampling, GRAHSP progressively builds up an implied grid.
However, the grid is not necessarily fully constructed. The cache keeps only a
user-defined maximum number of models (the CACHE_MAX environment variable),
dismissing from memory models that have become uninteresting as the fit
proceeds. The caching speeds up computation when the available machine memory
is large. More specifically, the main run script of GRAHSP, dualsampler.py,
entertains two pipelines, as illustrated in Figure 9. One pipeline creates the
AGN components with mock galaxy components and a second one the galaxy
components with mock AGN components. The mocked components only permit a
single value for each module parameter. This means that the implied grid for
each of the two pipelines is, as in AGNfitter or FortesFit, only concerned
with either the galaxy or AGN parameter space. In the final model evaluation,
the AGN components from the AGN pipelines and the galaxy component from the
galaxy pipeline are scaled and added together. These pipelines include the
module list that can be freely chosen by the user at run time. After the
combination of the galaxy and AGN components, several modules can be applied
relatively cheaply on-the-fly afterwards (see Fig. 9). This includes flux
computation, redshifting, biattenuation and also activepl. The expensive
computation of attenuation for each filter band is disabled during fitting,
and re-enabled during the analysis of the results. Since these modules are
evaluated on-the-fly, their parameter grids can be arbitrarily fine. This
enables a detailed exploration of the degeneracies between continuum shape
(activatepl) and its attenuation (biattenuation). It also enables the
incorporation of redshift uncertainties, either as a free fitting parameter or
with an informative prior from photometric redshift estimation.

GRAHSP is also parallelised. SED analysis can be trivially split across
computing cores without communication when analysing more objects than
computing cores (known as an embarrassingly parallel problem). The user can
set the number of cores to use (\--cores command-line option) and choose among
several parallelisation backends provided by the joblib and multiprocessing
python libraries. Here, a further modification to CIGALE is necessary. CIGALE
stores its database of model templates in a sqlite file on disk. However,
mutual locking of multiple parallel runs accessing the sqlite file can cause
substantial delays on shared file systems, as typically happens in computing
centres. To avoid this, if the environment variable DB_IN_MEMORY is set to 1,
the database ($2$\mathrm{GB}$$) is copied to memory on start-up, which brings
substantial speed-ups on large machines.

####  2.3.3 Module Pipeline

The order of execution of the modules is: sfhdelayed, bc03 (or m2005),
nebular, activate, activatelines, activatetorus, activatepl, biattenuation,
galdale2014, redshifting. The module galdale2014 has to be placed after
biattenuation to receive the attenuated luminosity.

####  2.3.4 Parameters and priors

Table 4 summarizes the modules (in square brackets) and their parameters
introduced in the previous sections. Figure 1 gives an impression of the
shapes of the model components with arbitrarily chosen galaxy and AGN
parameters. For stellar mass, AGN luminosity and redshift, a continuous
parameter space is explored, and the user can specify priors through the data
file and command line options. The redshift can be constrained with a fixed
(spectroscopic) value, or with uncertainties444Supporting arbitrary
photometric redshift probability distributions is planned for future releases.
(photometric), or left unconstrained. The AGN luminosity can be constrained by
providing a log-normal flux prior, for example from an X-ray detection. To
test fits without a host galaxy, the upper limit on the stellar mass can be
lowered through a command-line option. The other parameters are specified
through a configuration file (pcigale.ini). The meaning of the AGN parameters
is illustrated in Fig. 3. The effect of the two E(B-V) parameters is shown in
Figs. 5 and 6. As in CIGALE, the specified grid parameter values also define a
prior. numpy functions can be used to generate e.g., log-uniform grids.

The values chosen for the tests in this work are listed in the right-most
column of Table 4. For the galaxy, these follow previous CIGALE-based fits by
Ciesla et al. (2015) and Yang et al. (2020). For the AGN component, we
generally follow the ranges found in Section 2.2.

Table 4: Modules and their parameters. GRAHSP allows the user to specify the parameter ranges or values as they see fit. The last column lists the choices for this work. In cases of discrete parameters, grid values are specified. For continuous parameters, probability density functions are adopted. Parameter | Unit |  Prior
---|---|---
Galaxy components: |  |
stellar_mass | $M_{\odot}$ |  log-uniform between $10^{5}$ and $10^{\mathtt{mass\\_max}}M_{\odot}$
mass_max | $\log M_{\odot}$ |  15 (command-line parameter)
[sfhdelayed] |  |
tau_main | Myr |  100, 200, 500, 1000, 3000, 5000, 7000, 10000
age | Myr |  158 to 10000 in 18 log-uniform steps
sfr_A |  |  1
[m2005] |  |
imf |  |  0 (Salpeter)
metallicity |  |  0.02 (solar)
separation_age | Myr |  10
[nebular] |  |
logU |  |  -2.0
f_esc |  |  0.0
f_dust |  |  0.0
lines_width | km/s |  300
[galdale] |  |
alpha |  |  uniform between 0.75 and 2.75 in steps of 0.0625
AGN components |  |
[activate] |  |
L_AGN | erg/s |  log-uniform between $10^{38}$ and $10^{50}$, or
|  |  log-normal with parameters taken from the data file columns ‘FAGN’ (mean) and ‘FAGN_err’ (standard deviation), shifted with the luminosity distance at redshift.
[activatelines] |  |
AFeII |  |  0.6, 1, …, 32 (10 logarithmic steps)
Alines |  |  0.3, 0.5, 0.7, 1, 1.5, 2, 4, 10, 20
linewidth | km/s |  10000
[activategtorus] |  |
fcov |  |  uniform between 0.05 and 0.95 in steps of 0.05
Si |  |  uniform between -4 and +4 in steps of 0.2
COOLlam |  |  uniform between 10 and 30 in steps of 0.01
COOLwidth |  |  uniform between 0.2 and 0.65 in steps of 0.05
HOTfcov |  |  0.04, 0.1, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 2.0, 2.5, 3, 5, 10
HOTlam |  |  uniform between 1 and 5.5 in steps of 0.01
HOTwidth |  |  uniform between 0.2 and 0.65 in steps of 0.05
[activatepl] |  |
uvslope |  |  0
plslope |  |  uniform between -2.7 and -1 in steps of 0.01
plbendloc | nm |  50\. 80. 90. 100. 120. 150.
plbendwidth | dex |  log-uniform between 0.1 and 10 with 10 steps
attenuation & redshifting |  |
[biattenuation] |  |
Law |  |  Prevot
E(B-V) |  |  log-uniform between 0.01 and 10 with 80 steps
E(B-V)-AGN |  |  log-uniform between 0.01 and 0.1 with 80 steps
[redshifting] |  |
redshift |  |  spec-z mode: fixed value, from ‘redshift’ column of input file (used in this work)
|  |  photo-z mode: if ‘redshift_err’ is also given: Weibull distribution matched in mean and standard deviation to ‘redshift’ and ‘redshift_err’, except
|  |  no-z mode: if ‘redshift_err’ is negative: uniform distribution between 0.001 and 6.

###  2.4 Fitting method

This section describes the procedure for constraining model parameters from
observational data. In many fields of astronomy, low signal-to-noise ratio
data can be interpreted by fitting parametric physical models using Bayesian
fitting algorithms. In contrast, in SED fitting and today’s precision
photometry surveys, it is common that the data quality exceeds that of the
model. This situation is known as inference under model mis-specification.
There is little benefit from advanced algorithms when the results are
primarily limited by the quality of the model. Instead, it is more important
to achieve fast fits, test robustness by varying assumptions, and understand
the systematics in the results. In this work, we include improvements in the
model (previous section), the likelihood (section 2.4.1), model uncertainties
(section 2.4.2) and the Bayesian sampling algorithm exploring the parameter
space (section 2.4.3). Finally, section 2.4.4 illustrates the outputs of the
fit.

####  2.4.1 Likelihood

First, we review the $\chi^{2}$ likelihood implemented in CIGALE and LePhare
(Arnouts et al., 1999). Assuming photometric fluxes $f_{i}$ were measured with
associated uncertainties $\sigma_{i}$ for a set of bands $i\in\mathcal{B}$,
the usual Gaussian likelihood $\mathcal{L}$ is adopted:

$\displaystyle{\mathcal{L}}=\prod_{i}\frac{1}{\sqrt{2\pi\sigma^{2}_{i}}}\times\exp\left[{-\frac{1}{2}\left(\frac{F_{i}-f_{i}}{\sigma_{i}}\right)^{2}}\right].$
(7)

If the $\sigma_{i}$ are fixed during the analysis, one can rewrite this as

$\displaystyle\chi^{2}=-2\times\log{\mathcal{L}}=\sum_{i}{\left(\frac{F_{i}-f_{i}}{\sigma_{i}}\right)^{2}}+\mathrm{const}.$
(8)

To maximise the likelihood (and minimise the $\chi^{2}$), the optimal model
normalisation factor can be found analytically. The profile likelihood variant
$\chi_{\mathrm{prof}}^{2}$ is then:

$\chi_{\mathrm{prof}}^{2}=\sum_{i}{\left(\frac{s\times
F_{i}-f_{i}}{\sigma_{i}}\right)}$

with

$s=\sum_{i}\frac{F_{i}\times
f_{i}}{\sigma_{i}^{2}}/\sum_{i}\frac{F_{i}^{2}}{\sigma_{i}^{2}}.$

This is how CIGALE and LePhare optimise the normalisation analytically without
exploring this additional parameter. In particular, for each model, CIGALE
computes the optimal normalisation, and scales proportional parameters,
including the stellar mass and AGN luminosity. Thus, the posterior
distribution of stellar masses is determined primarily by the model template
shapes, but does not fully consider the flux uncertainties. This can be
demonstrated by fitting a single photometry data point with 20% flux error in
CIGALE with a single template model (all parameters fixed). The uncertainties
on the stellar mass should then trivially be 20%. However, CIGALE finds an
uncertainty of 0%, which is then heuristically raised to 5%555see line 141 in
pcigale/managers/results.py, commit 805612292. This shows that the CIGALE
uncertainties on normalisation-related parameters can be severely under-
estimated. In GRAHSP, we instead use a Bayesian approach and resolve this
issue with the full likelihood (eq.7) instead of a profile likelihood. The
normalisation(s) are free model parameter(s), as e.g., in FortesFit or
AGNFitter.

GRAHSP supports flux upper limits. Upper limits are a property of the
detection process (see Kashyap et al., 2010), and describe the probability
distribution that a source of a flux $f$ escaped detection by chance. We
assume that the cumulative probability distribution up to flux $f$ can be
described by the cumulative probability density of a Gaussian with mean
$f_{j}$ and width $\sigma_{j}$, for each band $j\in\mathcal{B}$ with a non-
detection. These are consistently incorporated in the likelihood as additional
multiplication terms:

$\displaystyle{\mathcal{L}}^{\prime}={\mathcal{L}}\times\prod_{j}\int_{0}^{f_{j}}{\frac{1}{\sqrt{2\pi\sigma_{j}^{2}}}\times\exp{\left[-\frac{1}{2}\left(\frac{F_{j}-f}{\sigma_{j}}\right)^{2}\right]}\,df}.$
(9)

This is the same treatment as implemented in CIGALE. Users provide upper
limits by setting in the data file the flux uncertainty to $-\sigma_{j}$
(negative values) and the flux to $f_{j}$. Forced photometry at a position
known from another band, may indicate a non-significant flux detection for
that band. Following Kashyap et al. (2010), such upper bounds are treated as
normal flux measurements (see the beginning of this section). In addition to
measurement uncertainties, additional systematic uncertainties have to be
considered.

####  2.4.2 Variability and model uncertainties

This section discusses the three types of uncertainties that GRAHSP considers:
measurement uncertainties, model uncertainties, and stochasticity introduced
by AGN variability. In the implementation, the total $\sigma$ used in the
Gaussian likelihood (eq. 7 and 9) is a combination with Bienaymé’s identity of
three contributions:

$\displaystyle\sigma_{\mathrm{tot}}=\sqrt{\sigma_{\mathrm{obs}}^{2}+\sigma_{\mathrm{sys}}^{2}+\sigma_{\mathrm{var}}^{2}}.$
(10)

The flux measurement uncertainty $\sigma_{obs}$ is provided by the input
photometric catalogue. The model uncertainty $\sigma_{sys}$ is computed for
each filter as:

$\sigma_{i,\mathrm{sys}}=f_{\mathrm{systematic\\_deviation}}\times
F_{\mathrm{obs}}+f_{\mathrm{sys}}\times F_{i,\mathrm{AGN}}.$ (11)

The first term is already used in CIGALE to account for systematic uncertainty
in the galaxy model, with $f_{\mathrm{systematic\\_deviation}}$ set to 10 per
cent by default and also here. In addition, in GRAHSP the systematic
fractional AGN model uncertainty, $f_{\mathrm{sys}}$, is a free-fitting
parameter (labeled systematics in the output). It scales the uncertainty with
the photometric band’s model flux, $F_{i,\mathrm{AGN}}$, considering only AGN
components. Large values of $f_{\mathrm{sys}}$ provide small $\chi^{2}$
values, however, due to the normalisation terms in eq.7, these are disfavored
by their low probability density. We assign $f_{\mathrm{sys}}$ an exponential
prior distribution centered at zero with a user-definable scale, set to
systematics_width=0.2% in this work. The heavy tails allow poorly fitted data
to increase the parameter to values much larger than systematics_width (see
e.g. discussion in Gelman, 2006).

The inclusion of this systematic uncertainty has important implications. In
standard practice, when the model produces a poor fit (large $\chi^{2}$
values), the model parameter uncertainties are typically too small. The fit
then either has to be discarded (like in LePhare) or the resulting parameter
distributions cannot be trusted. Here, instead, when the model fits the data
poorly, $f_{\mathrm{sys}}$ becomes large, which in turn produces larger
uncertainties on the fit parameters. Therefore, even in situations where the
model poorly describes the data, the estimates and uncertainties by GRAHSP can
be meaningful. This enables systematic analyses of samples and the direct use
of the fitting results. Additionally, $f_{\mathrm{sys}}$ can indicate when the
model is a poor fit.

The third term in Eq. 10 is the variance introduced by AGN variability. It is
active when variability_uncertainty is set to True in the configuration. This
term is also proportional to the AGN model fluxes,
$\sigma_{var}=F_{i,\mathrm{AGN}}\times F_{\mathrm{var}}(L_{\mathrm{AGN}})$.
From studying the year-to-year multi-band photometry variability of a sample
of X-ray selected AGN in Pan-STARRS1, Simm et al. (2016) found a relation for
the fractional variability independent of wavelength and black hole mass,
which we adopt as:

$\displaystyle F^{2}_{\mathrm{var}}=\mathrm{NEV}=\min(0.1,10^{-1.43-0.74\times
l_{45}}).$ (12)

Here, $l_{45}$ is the logarithm of the bolometric AGN BBB luminosity
$L_{\mathrm{bol}}$ in units of $10^{45}\,\mathrm{erg/s}$. To summarise, the
AGN variability uncertainty term depends on the AGN luminosity, with low-
luminosity AGN having the largest variations. The practical implication is
that photometry in similar wavelengths collected over several years from
multiple surveys can be consistently included. Indeed, if the photometry shows
statistically significant variability, this will make the model fit prefer an
AGN-dominant model, because only the AGN components can induce additional
variance.

####  2.4.3 Bayesian sampling algorithm

The continuous and discrete parameters listed in Table 4, plus the systematic
model uncertainty parameter, define a challenging 21-dimensional parameter
space. The model parameter space is likely to involve non-trivial degeneracies
that should be fully explored. This is evident from the flexibility of the
model illustrated in Fig. 3 pinned down by perhaps only a small number of
photometric observations.

The nested sampling (Skilling, 2004; Ashton et al., 2022) Monte Carlo
algorithm can perform well in this setting. Nested sampling estimates the
parameter posterior probability distribution and the Bayesian evidence (useful
for model comparison), by maintaining a population of live points sampled from
the prior. At each nested sampling iteration, the lowest likelihood live point
is discarded, and the sampling continues, however, with the constraint that
any new prior samples must exceed the likelihood of the discarded point. This
has two effects: at each step, the likelihood threshold is raised such that
the prior volume shrinks by a factor of K/(K+1), where K is the number of live
points. Secondly, the population of live points initially globally samples the
parameter space, but gradually focuses more towards the likelihood peak. The
posterior probability distribution is estimated by posterior samples. These
are the discarded points, weighted by the volume discarded at the
corresponding iteration i, $V_{i}=1/K\times(K/(K+1))^{i}$, and the likelihood
of the point, $w_{i}=\mathcal{L}_{i}\times V_{i}$. The evidence is estimated
as $Z=\int{\mathcal{L}\times\pi(\theta)\,d\theta}\approx\sum_{i}w_{i}$. In
practice, after many iteration the likelihood peak is reached and
$\mathcal{L}$ plateaus while $V_{i}$ declines exponentially; thus the
iterations can be stopped. The remainder of the live points is also added to
the posterior and evidence integrals. For rapid progress towards the posterior
peak, we choose K=50 live points.

Iterative progress is made in nested sampling by adding a new live point
sampled from the prior, with the requirement that its likelihood is higher
than that of the discarded live point. In the first few nested sampling
iterations, the likelihood-restricted prior is similar to the entire prior.
Therefore, it is efficient to sample directly from the prior using rejection
sampling. We adopt the ellipsoidal rejection sampling technique of Mukherjee
et al. (2006) for the first 10000 model evaluations, and then switch to slice
sampling. Slice sampling (e.g. Jasa & Xiang, 2012; Handley et al., 2015)
efficiently samples from a likelihood-restricted prior in high-dimensional
parameter spaces. We start from a randomly selected live point and perform 20
slice sampling steps before the point is accepted as a new live point. A slice
sampling step proposes a direction vector, and along this direction, proposal
points are drawn until the proposal exceeds the likelihood threshold. This can
be made efficient (see Kiatsupaibul et al., 2011) by a stepping-out procedure
identifying slice end points. For a literature review of techniques, see
Buchner (2021a). For proposing a direction vector, Buchner (2022) compared
various proposals with a diverse set of test problems. The most rapidly
converging proposal alternates randomly between differential vectors of random
pairs of live points and randomly chosen principal axes of the live point
distribution. We adopt this slice sampling technique until the posterior
weight of the remaining live points is negligible (¡1%).

The nested sampling procedure described above is implemented in the open-
source Python library
UltraNest666https://johannesbuchner.github.io/UltraNest/. We connect GRAHSP to
UltraNest to obtain posterior samples, Bayesian evidence estimates,
diagnostics and visualisations of the posterior. UltraNest allows GRAHSP to
resume from previously started runs, which enables modifying post-processing
analysis outputs (see the next section) without refitting from scratch.

####  2.4.4 Uncertainty quantification

For each analysed source, an output folder is created. UltraNest places the
posterior samples and other nested sampling outputs777see
https://johannesbuchner.github.io/UltraNest/performance.html#output-files
there. From the posterior samples the corresponding SED model is calculated.
For computational efficiency, by default GRAHSP process only 50 posterior
samples (\--num-posterior-samples command line option). Posterior predictions
of the individual model components are made and visualised. For each
photometric band selected in the configuration file, predicted model fluxes
for the entire model, and the AGN and galaxy components alone, are computed.
This allows comparisons to other datasets not used in the fit or predictions
for future observations. The fit parameters and marginal posterior
distributions of derived quantities, such as the star formation rate within
the last 100 Myr, the $12\,$\mathrm{\SIUnitSymbolMicro m}$$ AGN luminosity or
the NEV, are also created, and corner plots can be optionally created. For all
visualisations, the information is stored in text files, so that the user can
replot the data. Finally, an output file listing all the analysed sources with
the inferred parameter values and derived quantities, is created. The
posterior probability distribution is summarised with the median (suffix _med
to the analysed parameter name), standard deviation (_std), arithmetic mean
(_mean), geometric mean (_lmean), and 2-sigma equivalent lower and upper
quantiles (_lo, _hi). The geometric mean can be beneficial for log-quantities
such as luminosities and masses.

##  3 Data

We aim to thoroughly test our code and model with diverse data sets. These
include pure galaxies (without detected AGN activity) from deep pencil-beam
surveys, and large samples of AGN selected in X-ray, optical and infrared
wavelengths. Table 6 in Appendix C presents an overview of the sample sizes,
photometric bands and depths.

###  3.1 Dataset from COSMOS

The 2 deg2 COSMOS survey field (Scoville et al., 2007) is one of the prime
extragalactic survey fields. Deep multiwavelength observations from UV/optical
(down to 26 mag) and infrared enable the study of faint and distant galaxies.
In X-rays, it has been fully observed with XMM first (Hasinger, 2008) and
Chandra later (Civano et al., 2015) down to a flux limit of
$\qty{2.2d16}{\per\centi^{2}\per second}$ in the 0.5-2 keV band, making it
easy to detect also faint AGN (Brusa et al., 2010; Marchesi et al., 2016).
Furthermore, extensive and deep spectroscopy campaigns (COSMOS Deep, zCOSMOS,
MOSFIRE, DEIMOS, LEGA-C, just to name a few) were conducted over the years,
making COSMOS the ideal benchmark for testing SED fitting techniques. For the
spectroscopy, we used the compilation that was available within the COSMOS
collaboration in 2017 and considered only extragalactic sources with secure
redshift (reliability 100%), following the ranking introduced by Lilly et al.
(2009). We take UV to mid-IR photometry from the COSMOS2015 catalog (Laigle et
al., 2016), namely near-UV from Galex, optical from Subaru broad (BVrizy),
intermediate and narrow-band filters, near-infrared from the wide-field camera
(WFCAM, J band), Canada-France Hawaii Telescope (CFHT, H and Ks band), and
mid-infrared bands 1-4 of Spitzer’s infrared camera (IRAC). We further cleaned
the sample by discarding all sources with a flag indicating any flux
extraction problems such as saturation, and those with neighbors within 2 arc
seconds in HST i-band imaging to avoid deblending problems. Fluxes from 3
arcsec apertures are preferred, which is sufficiently large to encompass most
galaxies. AB magnitudes are converted into fluxes. Magnitude errors are
conservatively converted into flux errors by comparing the difference between
the flux of the nominal magnitude and the flux obtained for a 1$\sigma$ lower
magnitude. Flux measurements are corrected for Milky Way extinction with the
values listed in Table 3 of Laigle et al. (2016) for each band, multiplied by
the catalogued galactic E(B-V) values.

After this basic cleaning, we have created two sub-samples: $\bullet$ COSMOS-
pure-galaxies: in this sub-sample, all AGN candidates have been removed based
on either (a) broad emission lines in the optical VLT (Lilly et al., 2009) and
DEIMOS (Hasinger et al., 2018) spectra, (b) detection by XMM-Newton (Brusa et
al., 2010) or Chandra (Civano et al., 2012; Marchesi et al., 2016), (c) an
entry in Simbad marking it as an AGN or AGN candidate, (d) the IRAC-based
infrared AGN selection criterion of Donley et al. (2012). While some AGN may
remain in this sample, they are likely limited to low luminosities. The final
pure galaxy sample contains 3367 galaxies.
$\bullet$ X-ray detected AGN: The second sub-sample instead consists of X-ray
selected AGN taken from the homogeneous catalog of Marchesi et al. (2016)
based on the Chandra (Civano et al., 2012) and Chandra-Legacy (Civano et al.,
2015) surveys of the COSMOS area. Out of the 4016 sources listed there, we
retain 225 AGN.

The sample of X-ray-detected AGN is used to quantify the stellar mass
estimation as a function of: a) redshift accuracy: Often spectroscopic
redshifts are not available and photometric redshift are challenging to obtain
for AGN (e.g., Salvato et al., 2019). For the spectroscopic redshift sample,
we compute masses adopting the spectroscopic redshift. Additionally, to test
the impact of uncertainty induced by photometric redshifts, we also estimated
the masses adopting the photometric redshift value from Marchesi et al. (2016)
computed as in Salvato et al. (2011). Their accuracy and fraction of redshift
outliers are at the few per cent level. b) depth of the multiwavelength data:
Unlike deep pencil-beam surveys, wide-area surveys are characterized by a
limited number of shallow photometric points. To emulate this situation, the
photometry of these 225 AGN is modified in two steps: first, when the observed
flux falls below the flux limits of current all-sky surveys, the measurement
was discarded. The limiting depth in AB magnitudes was assumed to be:
Galex/near-UV 20.8, B (23.15), r (22.70), i (22.20), z (20.70), yHSC
(20.5+0.634888Vega to AB conversion), J (20+0.938), Hw (18.8+1.379), Ksw
(18.4+1.9), IRAC1 (16.6+2.699), IRAC2 (15.6+3.339), IRAC3 (11.3+5.174) and
IRAC4 (8+6.620). Secondly, the photometric errors were increased by setting
the flux errors to a fraction of the measured flux. That fraction is chosen
from the mean fractional flux error in COSMOS of the respective band. This
simple approach follows the assumption that both COSMOS and the emulated
surveys are signal-to-noise limited, and the vast majority of sources are
found at the faint end.

Figure 10: r magnitude as a function of redshift for the samples considered in
this work. The COSMOS X-ray selected AGN (big blue circles) are, on average,
fainter and at higher redshifts than all-sky WISE or X-ray selected quasars
(red). The Chimera objects (crosses) combine COSMOS galaxies (pink) with
different levels (colour bar) of quasar light (from optically selected
quasars). Figure 11: WISE diagnostic plot for the same samples as shown in
Fig. 10. Sources dominated by galaxy light tend to lie towards the bottom-
right. With increasing AGN luminosities (colour bar), sources move to the
upper-left. Magnitudes are in AB.

###  3.2 All-sky selected AGN

The AGN detected in X-ray pencil-beam surveys are not capturing the entire
range of the AGN population, especially those with the most vigorously
accreting black holes. For this reason, we test GRAHSP on AGN samples from
wide- or all-sky surveys, selected in infrared, X-rays and optical. The
following Section 3.2.1 and Section 3.2.2 provide the details for the
construction of the samples and the extraction of photometry.

####  3.2.1 DR16QWX & eFEDS: all-sky infrared and X-ray AGN

In this section, we identify a complementary X-ray and mid-infrared sample of
spectroscopic type 1 quasars. For our DR16QWX sample, the starting point is
the QSO data release 16 (DR16Q) from SDSS, presented in Lyke et al. (2020).
The DR16Q lists quasars targeted because of infrared and X-ray selection, but
also includes previously targeted quasars that have X-ray counterparts. In
terms of infrared selection, the DR16Q includes WISE (Wright et al., 2010) AGN
selected in the Stripe82X survey field (LaMassa et al., 2019). In terms of
X-ray selection, the DR16Q includes quasars that Salvato et al. (2019)
identified as WISE counterparts to the X-ray all-sky ROSAT/2RXS (Boller et
al., 2016) and XMMSLEW2 (Saxton et al., 2008; XMM-SSC, 2018). For combining
XMMSLEW2 with DR16Q, we find the Salvato et al. (2019) counterpart by a 1”
position match of the respective WISE counterparts. For 2RXS, the Salvato et
al. (2019) counterparts are listed directly in DR16Q. In addition, the DR16Q
also includes the optical counterparts to XMM sources from XMMPZCAT (Ruiz et
al., 2018) and the 3XMM-DR8, both based on the third version of the XMM-Newton
Serendipitous Source Catalog (Rosen et al., 2016). Therefore, from the DR16Q,
we select quasars selected by WISE, ROSAT or XMM. Only quasars (classified as
normal quasar or broad-absorption-line quasar in Table 2 of Lyke et al., 2020)
were kept. Furthermore, the b-band galactic extinction is required to be less
than 0.1 magnitudes, which limits the sample to high galactic latitudes. Of
the remaining sources, all show high redshift confidence (above 0.35). The
vast majority of these (80%) were not targeted by the optical selection of
SDSS DR7 (Schneider et al., 2010), making the selection complementary to
bright, blue quasars. We restrict the sample to high quality redshifts
(ZWARNING=0 and SN_MEDIAN_ALL¿1.6, see e.g. Menzel et al. (2016)), and to not
be in DR7 (IS_QSO_DR7Q¡1). This leaves 2349 sources. Finally, to focus on non-
jetted AGN, sources within one arc second of source positions listed in the
CRATES flat-spectrum radio source catalog (Healey et al., 2007) or the BZCAT
blazar catalog (Massaro et al., 2015) are removed. Our final DR16QWX sample
includes 2296 AGN.

For an eROSITA sample, we adopt the extragalactic sources from the eROSITA
Final Equatorial-Depth Survey (eFEDS Brunner et al., 2021; Salvato et al.,
2022). We focus on sources with spectroscopic redshifts above 0.002, and
remove jetted AGN as described above.

Next, the UV to infrared photometry is collated with a one arcsecond search
radius around the optical position. We carefully construct photometry from
near-UV to mid-IR with matched apertures. Our implementation is released with
this paper as RainbowLasso999https://github.com/JohannesBuchner/RainbowLasso.
Near-UV magnitudes from Galex (Bianchi et al., 2017) are converted into fluxes
(see Section 3.1). Measurements below 20.8 AB mag, the nominal Galex survey
depth, are discarded. Co-extracted optical (g, r, i, z) and WISE (W1-W4)
photometry is obtained from the Legacy Survey data release 10 (LS10 Dey et
al., 2019). Sources were discarded if the photometry extraction fit was
forced, hit limits, or the parameters had to be held fixed101010see
https://www.legacysurvey.org/dr10/bitmasks/#fitbits. Heavily blended sources
are also discarded. If a multi-source fit assigned less than 10 per cent of
flux in the g, r or W2 bands to other sources (fracflux_* columns), the source
is considered isolated. Final sample sizes are listed in the first row of
Table 6. Because of the large point-spread function in the W3 and W4 bands, we
discard W3 and W4 photometry if other sources dominate in either band or the
contribution in W1 or W2 exceeds 10%. We consider point sources those where a
point-spread function model provided the best fit (TYPE=“PSF”), and the
remainder extended sources. For extended sources, we use fluxes extracted from
5 arcsecond apertures (apflux column 7 for the optical bands and column 5 for
WISE bands), while for point sources, we use the model flux. Fluxes are
corrected for galactic attenuation by dividing by the respective
MW_TRANSMISSION columns. Near-infrared photometry is added from the Visible
and Infrared Survey Telescope for Astronomy (VISTA) Hemisphere survey (VHS)
data release 5 (McMahon et al., 2013, 2021), the United Kingdom Infrared
Telescope (UKIRT) and InfraRed Deep Sky Survey (UKIDSS) (Lawrence et al.,
2007; Warren et al., 2007). For point sources, we use aperture-corrected PSF
model fluxes extracted from 2.8 arcsecond diameter apertures (apermag4), while
for extended sources, we use not-aperture corrected fluxes from 5.6” diameters
(aper*corr6). We also retain the X-ray flux from the source X-ray catalogs,
where available, and convert it to a $2-10\,$\mathrm{keV}$$ luminosity
assuming a power law index of 1.8, as most sources are only mildly absorbed.

####  3.2.2 SDSS-DR7Q: pure quasar sample

We construct a pure quasar sample that complements the COSMOS pure-galaxy
sample in Section 3.1. By “pure” we mean here that the emission is dominated
by AGN processes at essentially all wavelengths, and galaxy contribution is
negligible. Ultraviolet and optically selected quasars were targeted by SDSS
up to data release 7 catalogue (DR7 Schneider et al., 2010). Shen et al.
(2011, 2012) selected sources where the i-band absolute magnitude exceeds -22
and the SDSS spectra show at least one broad emission line. This resulted in
105,783 extremely luminous, unobscured quasars. We assume that these are
completely dominated by the AGN process and can be treated as point sources.
We use the Milky Way-corrected AB magnitudes in the catalogue of Shen et al.
(2011). Besides the optical bands (ugriz), the catalogue also provides near-
infrared Vega magnitudes (JHKs) from the Two Micron All Sky Survey (2MASS)
(Skrutskie et al., 2006) and infrared Vega magnitudes from WISE. The
magnitudes and errors (see Section 3.1) were converted into fluxes.

##  4 The Chimera benchmark data set

To date, no benchmark exists for characterising the bias and variance of AGN
host galaxy properties. These include fundamental physical quantities such as
stellar masses, star formation rates, and stellar ages. The fundamental reason
is that, unlike in the problem of estimating redshifts from photometry, true
labels are generally not available even with more costly observations
(spectroscopy). Comparisons of outputs from various codes and different models
indicate that the results vary with model assumptions (Mobasher et al., 2015;
Santini et al., 2012). Works on AGN hosts usually exclude from the analysis
the unobscured AGN (e.g., Aird et al., 2015), because they present large
uncertainties (e.g., Bongiorno et al., 2012). Ciesla et al. (2015)
demonstrated with simulated unobscured AGN how stellar mass and SFR are
respectively over and under-estimated, already with a 40 per cent light
contribution by the AGN. Reliable retrieval becomes even more challenging the
fewer photometric bands are available. However, because the model assumed for
simulating data was the same as the model for inference, such tests provide
limited insight into the true bias and scatter of estimated host parameters.

To account for these issues, we aim to design a benchmark dataset in which
fiducial truth physical galaxy properties are known. The Chimera111111In
biology, a chimera is an individual consisting of tissues of diverse genetic
constitution. In Greek mythology, a chimera is a fire-breathing she-monster
having a lion’s head, a goat’s body, and a serpent’s tail. benchmark for AGN-
galaxy decomposition is derived from real-world data without assuming a true
model for AGN. We rely on the inference of physical properties in pure
galaxies (non-AGN), as validated in Mobasher et al. (2015). We verified that
the stellar masses determined by GRAHSP are consistent with those presented in
the COSMOS2015 catalog (see the next section) and other codes (Section 5.5),
so that any of these methods can be used as a reference value. This
establishes fiducial truth physical galaxy properties, in the sense that we
know what would have been measured if the galactic nucleus were inactive. In
the second step, pure quasars from the DR7Q dataset (Section 3.2.2) are
randomly paired with pure galaxies (Section 3.1), under the constraint that
the spectroscopic redshifts are within 0.01 of each other. The fluxes in each
photometry band are then summed. For this step, the optical and near-infrared
bands (urizJHK) can be summed directly, as the filter curves are in close
agreement. For IRAC and WISE bands 1 and 2, tiny conversions from Stern et al.
(2012) are applied. The higher IRAC bands substantially differ from WISE 3 and
4, and we do not consider them. The fluxes are summed, with a weighting factor
applied to the quasar light, which controls how luminous a quasar is added.
Uncertainties are propagated assuming Gaussian noise. Fluxes are only provided
for bands where both the galaxy and the quasar had a measurement.

The Chimera benchmark dataset is enriched with reference galaxy and AGN
properties. These are the galaxy properties (stellar mass, SFR, stellar age,
attenuation, etc.) from COSMOS2015 and AGN properties (bolometric and
$5100$\mathrm{\text{Å}}$$ quasar luminosities, black hole mass) from Shen et
al. (2011). The Chimera benchmark is
released121212https://doi.org/10.5281/zenodo.8431646 with this paper, so that
it can be used as a reference for comparing physical parameter estimates for
AGN in codes currently in use or under development.

Figure 12: Colour diagnostic plot. The GRAHSP templates (coloured dots as a
function of redshift) cover the entire plane. The observations are illustrated
by coloured contours containing 99% of each sample. At the centre-left,
quasars dominate (red outline, DR16QWX), while the eFEDS sample includes a
wider distribution. The dashed diagonal line demarks the separation of Salvato
et al. (2022) between galaxies (above) and stars (below the line).

##  5 Preliminary Tests on GRAHSP

###  5.1 Model plausibility

As a first check, model colours are verified against colours of observed
samples of AGN and galaxies. Figure 12 shows the g-r-z-W1 colour-colour
diagnostic plot presented in Salvato et al. (2022). We draw random samples
from the GRAHSP model prior, construct the template and present the model
colours as dots in Fig. 12. These cover the entire parameter space and beyond.
The model dots envelop the distribution of quasars from DR16QWX and eFEDS AGN,
which are illustrated as contours in Fig. 12. The locus of stars (see Salvato
et al., 2022) is below the dashed line, and as expected, few galaxy model
templates are in this region. Appendix D extensively tests the model colours
against four observational samples (DR16QWX, eFEDS, COSMOS AGN and COSMOS pure
galaxies). There it is shown that emission lines and the FeII template are an
important contributor, without which the observed colours cannot be
reproduced. This has also been noted in Temple et al. (2021), which studied
the mean quasar colour over a wide redshift range. Here, we demonstrate that
the GRAHSP model is capable of reproducing not only the mean colour but also
the wide diversity of colours in the galaxy and AGN population.

###  5.2 Retrieving stellar masses for COSMOS-pure-galaxies sample

Figure 13: Stellar masses from COSMOS2015 are compared with those from GRAHSP,
with the AGN component turned OFF (top panel) or ON (middle panel), for
various redshift bins. In both cases, the agreement is extremely good, with
very little scatter around the 1:1 relation. For the case with the AGN
component turned ON, the GRAHSP AGN luminosity is plotted as a function of
stellar mass. As expected, in all sources the AGN component is well below 10%
in all cases but two.

We first test GRAHSP on the COSMOS-pure-galaxy sample, using only galaxy
models. The top panel of Fig. 13 compares the resulting stellar mass estimates
to those published by Laigle et al. (2016). To ensure consistent assumptions,
we show only sources where their adopted redshift agrees with the
spectroscopic redshift within 1 per cent. The stellar mass estimates agree
well at all redshifts and at all stellar masses.

Secondly, we test the impact of introducing our AGN model components, even if
it is not required. On the same COSMOS-pure-galaxy sample, the middle panel of
Fig. 13 compares two runs: once without AGN from the previous experiment and
once with the full GRAHSP setup described in section 2. The estimates show
close agreement, indicating that the added model complexity does not bias the
measurement. In most cases, the upper error bar of the AGN luminosity is well
below ${10}^{42}\,\mathrm{erg/s}$, allowing for low-luminosity AGN, at most.
In the SED, the AGN is always sub-dominant compared to the host. This can be
quantified with the AGN-to-galaxy light contrast ratio,
$L_{\mathrm{AGN}}^{5100$\mathrm{\text{Å}}$}/M_{\star}/\Psi$, where $\Psi$ is
the solar light-to-mass ratio (see Section 2.1.4). The bottom panel of Fig. 13
shows that the contrast lies well below 10% in the vast majority of cases.
From the entire sample, only few have a blue SED that can be fit with a star-
burst or a luminous AGN, giving rise to large (but correct) uncertainties in
the contrast ratio. These can also be seen in the middle panel.

###  5.3 Template fit test

Figure 14: Example of an AGN-dominated SED fit from the DR16QWX sample. The
total model (black) and the individual components (colours, see legend) are
presented with the posterior mean (solid line) and 2 sigma equivalent
uncertainties (shaded areas). Observed fluxes are shown as blue data points
with blue squares and wide blue capped $3\sigma$ error bars. The overlaid blue
error bars with thin caps show the enlarged total error budget after the fit.
The filter curves of observed photometric bands are shown at the top, grouped
to the same colours by instrument. Predicted model fluxes are shown as red
points. These can deviate from the black curve because of the averaging over
the filter curve. The bottom panel presents the relative residuals of model
and observations. The title lists the source ID, redshift, $\chi^{2}$ (Eq. 8)
per number of data points and Bayesian model evidence. On the right, key
values are shown, including the stellar mass (in $M_{\odot}$), star formation
rate (in $M_{\odot}$/yr), the galactic and AGN extinction, and the
$5100\unit{}$ AGN luminosity (in erg/s). Lower case variables are logarithmic.
Below $2\unit{\micro}$, the model is dominated by the AGN disk, and above by
the torus. The third data point (z-band) is raised compared to the neighboring
data points, which can be understood by the $H\alpha$ emission line at that
redshift. The galaxy component has an extremely wide luminosity range
(purple), i.e., the stellar mass is unconstrained (see the first number on the
right) in this luminous AGN (see the fourth number on the right).  Figure 15:
Example of a host-galaxy-dominated SED, with elements as in Fig. 14. Here, the
SED is dominated by galaxy components (purple and green curves), including
nebular emission lines (grey). The stellar mass and star formation rate are
well constrained (numbers on the right), while the AGN luminosity has more
than an order of magnitude uncertainty around a low value of
$10^{41}\unit{\per}$. Here, model uncertainties (thin cap error bars) were
only slightly increased without requiring an AGN component.  Figure 16:
Parameter posterior probability distributions for two different objects: One
AGN-dominated source (blue) of Fig. 14 and one host-galaxy-dominated source
(red) of Fig. 15. Panels with thick histograms are model parameters. Panels
with thin histograms are derived parameters. The large bottom-right panel
shows the joint posterior distribution of AGN luminosity and stellar mass as
posterior samples.

Next, the ability to describe SEDs from a diverse sample of AGN by the models
in GRAHSP, is tested. We use here the DR16QWX sample of infrared and X-ray
selected AGN (section 3.2.1) and show three extreme cases. Firstly, Fig. 14
shows an AGN-dominated case. The model follows the data points closely. The
H$\alpha$ emission line at 0.653$\mu$m substantially contributes to the
measured flux in the relevant filter (blue square), which is also reflected in
the model flux (red point). This demonstrates the importance of emission
features in the templates. Figure 14 shows the total model and its components.
Because the fit is fully Bayesian with continuous amplitude parameters, the
posterior (shaded areas) is smooth. The corresponding parameter constraints
are shown in Figure 16 as blue histograms. In particular, the stellar mass
(top left panel) posterior is a flat distribution reaching up to approximately
$10^{10}M_{\odot}$, indicating that for this AGN-dominated source only an
upper limit can be estimated. A host-galaxy-dominated source is presented in
Fig. 15 as a second example. The galaxy template is constrained well.
Comparing the thin cap error bars with the wide cap blue error bars, we see
that the data uncertainties were slightly enlarged by the fit (see Section
2.4.4). The parameter constraints are shown in red in Fig. 16. The panels
marked ‘log_stellar_mass‘ and ‘sfh.sfr‘ show that stellar mass and star
formation rate could be constrained, respectively. The AGN luminosity (second
panel in the top row) is constrained to less than $\qty{e44}{\per}$.
Degeneracies between all model parameters can also be explored. The bottom
right panel of Fig. 16 only shows the conditional distribution of stellar mass
and AGN luminosity, with the obtained posterior samples.

Figure 17: Distribution of posterior medians of the systematic model error
parameter $f_{\mathrm{sys}}$ (see eq.10).

We investigate this a bit closer by looking at the inferred systematic
uncertainty parameter values. Figure 17 presents the distribution of posterior
medians on the systematic model error parameter fsys for the entire DR16QWX
sample. The median of the prior is 0.2 (vertical black dashed line). This
indicates that, for most sources, no additional model uncertainty had to be
introduced during fitting. For reference, the cases from Fig. 14 and Fig. 15
have a median fsys of approximately 0.08.

###  5.4 Extension to X-ray data

Figure 18:  Top panel: The AGN luminosity at $5100$\mathrm{\text{Å}}$$
inferred from GRAHSP is compared to the $2-10$\mathrm{keV}$$ luminosity from
X-ray observations in the DR16QWX sample. Points are colour-coded by redshift,
with yellow indicating the highest redshifts. Bottom panel: The ratio is shown
as a black histogram. For comparison, the grey dashed curve represents the
mean ratio and scatter observed by Koss et al. (2018).

Yang et al. (2022) investigated the use of X-ray constraints to improve the
model fit. As illustrated with the galaxy-dominated case (Fig. 15), in
particular low-luminosity and obscured AGN sources can benefit from additional
constraints on the AGN component. GRAHSP supports placing priors on the AGN
luminosity. Here we investigate whether X-ray observations are appropriate to
constrain the $5100$\mathrm{\text{Å}}$$ luminosity parameter of our model.

The upper panel of Fig. 18 shows the correlation between SED-inferred
$5100\unit{}$ luminosities and X-ray-based luminosities. The X-ray
luminosities were not used in the fit and are taken by conversion from X-ray
survey catalog fluxes (see section 3.2.1). They show a 1:1 correlation with
some scatter, particularly due to sources with large uncertainties. We
investigate the intrinsic scatter by considering the luminosity ratio (blue
histogram in the bottom panel of Fig. 18). The distribution is compared to
results from the local Universe in the literature, where Koss et al. (2018)
reported a $0.43\pm 0.47$ log-ratio distribution between the
$5100$\mathrm{\text{Å}}$$ and $14-195$\mathrm{keV}$$ luminosities, which the
authors relate to the $2-10$\mathrm{keV}$$ luminosity with a shift of 0.42
dex. Fortunately, the two shifts almost cancel out, making the
$5100$\mathrm{\text{Å}}$$ luminosity approximately equal the
$2-10$\mathrm{keV}$$ luminosity, with a scatter of 0.5 dex (gray distribution
in the bottom panel of Fig. 18). This agrees with our distribution of mean
ratios. This demonstrates that measured X-ray fluxes, blurred with a 0.5 dex
Gaussian, can be used as AGN luminosity priors in GRAHSP fits. Practically,
this is achieved by providing the necessary columns in the input file (see
Table 4).

###  5.5 Literature model setups

In the subsequent sections we test GRAHSP. We compare the results from GRAHSP
to those derived from three alternative methodologies that have been described
in the literature. All of these are based on CIGALE, and for our analysis we
adopt the latest version of that code (2022.1). We adopt broadly the same
setups that are described in the literature except that in all cases we use a
Planck cosmology Planck Collaboration et al. (2020) and the Maraston (2005)
SPS model, noting however that adopting the Bruzual & Charlot (2003) models
produces similar results.

The first model, “D14”, adds the AGN model of Dale et al. (2014) to the
stellar population. The slope is allowed to vary (values: 0.0625, 0.5, 1.0,
1.5, 2.0, 2.5, 3.0, 3.5, 4.0), and so is the AGN fraction (values: 0.01, 0.1,
0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8,
0.85, 0.9, 0.99). Since the AGN model has only two parameters, the galaxy
parameters can be evaluated on a fine grid while retaining an acceptable
computational cost, with $\tau$ taking values of 1, 50, 500, 1000, 2000, 2500,
3000, 3500, 4000, 4500, 5000, 5500, 6000, 7000, 8000 Myr and the age ranges
from 100 to 13 Gyr in 15 logarithmic steps. A modified starburst attenuation
law is applied.

The second model, “C15”, follows the setup of (Ciesla et al., 2015). The AGN
is modelled with the Fritz et al. (2006) template library. The full list of
parameters is given in Appendix A. Due to the high number of AGN parameters,
the galaxy parameters need to be more coarsely gridded in CIGALE, with $\tau$
taking values of 500, 1000, 3000, 5000, 10000 Myr and age taking values of
500, 1000, 3000, 5000, 7000, 10000. A Calzetti-Leiterer attenuation law is
applied (Calzetti et al., 2000; Leitherer et al., 2002).

The third model, “Y19”, follows the setup of Yang et al. (2019). The AGN model
is SKIRTOR (Stalevski et al., 2016), which also has many free parameters
listed in Table 5 in the appendix. The galaxy parameters are coarsely gridded,
with $\tau$ taking values of 100, 500, 1000, 5000 Myr and age taking values of
500, 1000, 3000, 5000, 7000.

The comparison of GRAHSP and the three CIGALE model setups is shown in Section
6. For a fair comparison, we compute the Chimera fiducial truth from COSMOS
with the galaxy model setup exactly matching the respective model setup.
However, adopting the COSMOS2015 results as fiducial truth does not change the
results, as there is good agreement on pure galaxy inference (see, e.g., the
top panel of Fig. 13).

GRAHSP has many more free parameters than the above setups do. For example,
SKIRTOR is the model with the largest number of free parameters, but for
computationally feasible fits with CIGALE, many need to be fixed. Appendix E
investigates whether having all parameters free can resolve any differences
between GRAHSP and CIGALE. In particular, similar to Section 2.2, we fit
SKIRTOR against the Brown et al. (2019) Atlas of AGN spectra in Appendix E,
and find substantial residuals.

###  5.6 Validation of the Chimera benchmark

Figure 19: Comparison of stellar masses measured from COSMOS-pure-galaxies
sample to that measured by GRAHSP for the Chimeras built using the same
galaxies and with low AGN contribution. The measurements are consistent; for
some cases, the uncertainties are very large.

(a) original COSMOS galaxy$M_{\star}=\qty{1.4e9}{M_{\odot}}$

(b) C15 model fit$M_{\star}={8\pm 3}\times\qty{e9}{M_{\odot}}$

(c) GRAHSP fit

Figure 20: Example of an outlier at the lowest Chimera quasar weight (0.0001),
which should be galaxy-dominated with very low AGN luminosity. The top panel
shows the original COSMOS galaxy. The observed fluxes (purple) are reproduced
with a high reduced $\chi^{2}$ (see title) by the galaxy model (black curve),
with substantial residuals in the near-infrared and optical wavelengths. The
middle panel presents the constructed Chimera fluxes (purple), fitted with the
C15 model’s galaxy (yellow curve) and AGN (orange curve) mixture model, using
CIGALE. The AGN component models the data in the infrared. The stellar mass
has become eight times larger than in the top panel. The bottom panel presents
the analysis of the same Chimera data with GRAHSP. Here the fit is dominated
by a low-luminosity AGN, and the stellar mass is poorly constrained. The
x-axis range is wider in this panel.

Before applying GRAHSP to the Chimera benchmark and evaluating the results, we
perform a sanity check on the created low-luminosity AGN objects. At
luminosities of $L(\mathrm{AGN})<{10}^{42}\,\mathrm{erg/s}$ (i.e., typically
corresponding to quasar flux weights of 0.1 per cent), the Chimera objects
have essentially the same fluxes as those from the original COSMOS-pure-galaxy
sample, but with two differences. Firstly, the number of filters is much more
limited compared to the COSMOS filter set used to establish the fiducial truth
stellar mass. Secondly, the filter curves are taken from SDSS, which do not
exactly correspond to the COSMOS filters where the galaxy fluxes are from
(Subaru and CFHT). To establish that this is not a concern, and that the
Chimeras have usable fluxes, Fig. 19 compares the stellar mass computed by
GRAHSP on these Chimeras to the original stellar mass from Laigle et al.
(2016). The stellar masses from the two analyses are consistent, with a mean
of 0.0 dex and a scatter of 0.2 dex. Similar results are found for the
comparison of CIGALE with the D14, C15 and Y19 models. There are a few (0.2%)
outliers, however.

One outlier is shown in Fig. 20. The top panel shows the analysis of the
original COSMOS data with CIGALE. The subsequent panel shows the Chimera
object with CIGALE (middle) and with GRAHSP (bottom panel). The COSMOS data
are fitted with a strongly star-forming galaxy template, which is required
because of the longest-wavelength infrared data points. Heavy attenuation then
suppresses the UV (compare the blue curve to the black curve). The model is
not a good fit, with a reduced $\chi^{2}$ of 3.9. The infrared colours are not
indicative of an AGN, following Assef et al. (2013). Instead, this indicates
that the galaxy model may be incomplete. The stellar mass in units of
$10^{10}\,M_{\odot}$, when analysing the COSMOS photometry, is 5 in the
COSMOS2015 catalogue, but 0.14 in CIGALE and 0.18 in GRAHSP when run without
an AGN component. When an AGN component is used to fit the Chimera object, the
infrared excess is modelled with an AGN component (Fig. 20, CIGALE: middle
panel; GRAHSP: bottom panel). This differs slightly in shape between the D14,
C15 and Y19 models, but the results are similar. The middle panel of Fig. 20
demonstrates this point: galaxy light is assigned to the AGN component. The
resulting stellar masses are now an order of magnitude higher, between 1 and
$\qty{2.2e10}{M_{\odot}}$, as much less attenuation is needed.

We include such cases where even the galaxy is not perfectly modelled. These
are realistic, and would be difficult to recognise if the uncertainties were
larger or fewer bands were used. They are relevant as the inference is
unstable when measuring a stellar mass of galaxies for comparison to AGN
samples. Nevertheless, such cases are rare. An order of magnitude or more
difference to the COSMOS reference masses are present in four per cent of the
Chimeras in the case of CIGALE, and six per cent with GRAHSP, as shown in
(Fig. 20). Their rarity implies that the COSMOS data analysis is valid as a
fiducial truth for Chimera objects. In the subsequent sections, larger AGN
luminosities are also considered, and the impact on the stellar mass estimate
reliability is determined.

##  6 Results

Figure 21: Retrieval of stellar masses on the Chimera benchmark data set. The
y-axis always shows the ratio of the estimated stellar mass to the true
stellar mass (x-axis), in log units of $M_{\odot}$. The first three rows are
AGN models in CIGALE, the last row is GRAHSP. The sample is split into four
bins (panel columns) of true AGN bolometric luminosity (in erg/s, log).
Results with large error bars ($>1$ dex) are in grey. If the measured stellar
mass deviates from the expectation (dotted grey horizontal lines) by more than
0.5 dex, it is considered an outlier and indicated with a circle. In each bin,
the mean and standard deviation are quoted, together with the fraction of
outliers in brackets. CIGALE models show a 0.2-0.5 dex bias at L$>$43\. For
GRAHSP, the bias is below 0.1 dex. The GRAHSP fraction of outliers is also
lower. The downward-pointing triangles at the very top of the figure indicate
where the AGN and galaxy luminosities are approximately equal ($\lambda=1$ in
Eq. 6). The CIGALE models start to deviate to the left of these markers,
corresponding to $\lambda>1$.

(a) original COSMOS galaxy$M_{\star}=10^{8}\,M_{\odot}$

(b) D14 model fit, $\chi^{2}_{\mathrm{red}}=1.11$

(c) C15 model fit, $\chi^{2}_{\mathrm{red}}=1.68$

(d) Y19 model Cigale fit, $\chi^{2}_{\mathrm{red}}=4.26$

(e) GRAHSP fit

Figure 22: Example of outlier in CIGALE-based analyses from the
$10^{43-44}\,\mathrm{erg/s}$ AGN luminosity bin. Panel (a) shows the original
COSMOS data (CosmosId: 420968 at z=0.33), where a
$M_{\star}=10^{8}\,M_{\odot}$ galaxy-only model (black) reproduces the
observed fluxes (purple) with small residuals ($\chi^{2}_{\mathrm{red}}=0.4$).
Panels (b), (c) and (d) present CIGALE fit to the Chimera fluxes (purple) with
D14, C15 and Y19 settings. In all three models, the galaxy template is
incorrectly placed at ${10}^{36}\,\mathrm{W}$ instead of
$<{10}^{35}\,\mathrm{W}$ as in panel (a). In the C15 and Y19 models, the best-
fit model is dominated by the AGN component (orange). The bottom panel (e)
shows the fit by GRAHSP for the same object. The galaxy component has very
large uncertainties (purple shade at the bottom), but includes the true
stellar contribution (compare the purple curve here to the black curve in the
top left panel). The error budget has slightly increased during the analysis
from the blue error bars with wider caps to those with thin caps. The power-
law AGN model (blue) dominates the fit at all wavelengths.

The Chimera benchmark described in Section 4 provides galaxy properties true
values that can be used for testing the capabilities of SED fitting methods.
In this section we test GRAHSP. We also compare it to CIGALE with the three
model setups (D14, C15, Y19; described in Section 7.3).

###  6.1 Stellar mass retrieval

First, we investigate the retrieval of stellar masses. We split the Chimera
sample into bins of true AGN bolometric luminosity. In Fig. 21, we compare the
true values against results from runs of GRAHSP and three CIGALE setups. At
low luminosities (first column), all estimates have relatively small
uncertainties, and centre around the true value (near zero deviation
$\Delta=\log M_{\mathrm{estimated}}/M_{\mathrm{true}}$). The content of the
bottom left panel was already presented in Fig. 19 as a correlation. In the
following, the mean of the deviations is referred to as the bias, and the
standard deviation as the scatter of the estimator. We weight the mean and
standard deviation by the error bar size. Finally, we mark outliers with
circles. These are defined as cases where the error bar excludes the 1 dex
wide band centred around the true value, illustrated with grey dotted
horizontal lines. Section 5.6 already discussed one example outlier in the
lowest luminosity bin from the D14, C15 and Y19 models.

The bias, scatter and outlier fraction are noted in the upper-right corner of
each panel. At higher luminosities (last column of Fig. 21), the settings most
commonly used in literature show a positive bias of +0.5 dex. This is most
pronounced at lower stellar masses, where the mass is over-estimated by an
order of magnitude or more. The outlier fraction is 20 per cent. In contrast,
GRAHSP (bottom row), has outlier fractions below 5 per cent and near-zero bias
($\leq 0.2\,\mathrm{dex}$). The scatter is comparable across all methods.

An example fit result for such an outlier is illustrated in Fig. 22. The top-
left panel shows the true COSMOS galaxy, modelled with the many photometric
bands in COSMOS. In the subsequent panels, the constructed Chimera object is
shown, with all luminosities about one order of magnitude higher
($\qty{1e36}{}$ instead of $\qty{1e35}{}$), due to the AGN. Even though the
galaxy should therefore be 10x fainter than the AGN, the Chimera object is
fitted by the D14, C15 and Y19 models (see panels) with approximately equal
luminosity between AGN (orange) and galaxy (yellow) contribution. The bottom
panel shows the result from GRAHSP. The more flexible AGN model of GRAHSP,
together with its fine parameter space sampling and consideration of
systematic uncertainties, can describe the data with an AGN-dominated
combination (bottom panel of Fig. 22), where the host galaxy stellar mass is
consistent with that of the COSMOS fiducial truth.

We caution that such failures could not be recognised from the fit alone.
Firstly, Fig. 22 shows $\chi^{2}$ values comparable to unity for some of the
models. Secondly, we implemented the cleaning criteria proposed by Mountrichas
et al. (2021), which require the maximum-likelihood stellar mass to agree with
the posterior mean stellar mass within a factor of five. This occurs in 100
per cent of the sample. Therefore, the criterion does not improve the outlier
fractions, which range from 5 to 20 per cent. Thirdly, agreement between
different model assumptions also does not resolve incorrect solutions, as
illustrated by the example of Fig. 22 and Fig. 20, which give across three
different CIGALE AGN models a consistent but incorrect high stellar mass.

Figure 23: Star formation rate retrieval on the Chimera benchmark data set. As
in Fig. 21. Each panel shows the deviation of the estimated SFR from the true
SFR (in log units of $M_{\odot}/\mathrm{yr}$, x-axis). The first three rows’
right panels show the performance of current models at moderate AGN
luminosities. Over-estimates by 1-2 dex are common at low star formation
rates, with outlier fractions between 10-30 per cent. The bias of GRAHSP (last
row) is lower and the outlier fraction consistently below 10 per cent.

###  6.2 Star formation rate

The next host galaxy property considered is the star formation rate (SFR). SFR
is known to be difficult to estimate in individual galaxies (e.g., Ciesla et
al., 2015). Nevertheless, works such as Aird et al. (2019) rely on sample-
averaged stellar masses, so quantifying the bias and scatter of estimators is
important. Here, the SFR in the last Myr is expressed in
$M_{\odot}/\mathrm{yr}$, and again compared to the true SFR from the Chimera
benchmark sample in Fig. 23. We only consider those objects where in the
reference COSMOS galaxy, the fiducial truth SFR could be constrained to better
than an order of magnitude by all models.

At the lowest true SFRs, the inferred SFR is overpredicted by one to three
orders of magnitude. At the highest SFR ($>10\,M_{\odot}$, the retrieval is
unbiased at least to the right order of magnitude. At the highest AGN
luminosities, GRAHSP shows the lowest bias. The D14 model (top panels) shows
an erratic behaviour at the highest SFR for AGN luminosities of $\log L>43$.
The other two CIGALE models, have a larger positive bias than GRAHSP. In the
CIGALE models, the outlier fraction lies between 20 and 30 per cent, while the
outlier fractions of GRAHSP remain below 10 per cent even at the highest
luminosity bin.

###  6.3 Luminosity-dependence of estimators

Figure 24: Comparison of codes and models on the Chimera benchmark data sets.
The panel rows show the outlier fraction (see Fig. 21), the bias and scatter
in dex. The first column shows stellar mass estimation (from Fig. 21), the
second column star-formation rate estimation (from Fig. 23). Each curve as a
function of AGN bolometric luminosity (log, in erg/s) represents one method.
The outlier fraction of GRAHSP is lower than those of other methods. The bias
of GRAHSP is near-zero. Other methods systematically overestimate the masses
in luminous AGN (middle left panel), and star-formation rates (middle right
panel).  Figure 25: AGN luminosity retrieval on the Chimera benchmark data
set. The top panel shows the deviation from the bolometric luminosity (in log
units of erg/s). Current methods (top panel) show an upward bias at
$L_{\mathrm{bol}}<\qty{e43}{\per}$, and an outlier fraction of 15 per cent.
The bolometric luminosity (see Section 2.1.4) from the BBB and the torus are
presented in the second and third panels, respectively. The bottom panel
compares the retrieval of $5100\AA$ luminosity. In all cases, the average bias
of GRAHSP is close to zero, with an outlier fraction lower than 5 per cent.

The trends with AGN luminosity are summarised in Fig. 24. The left panels show
the outlier fraction, bias and scatter of the stellar mass estimator, from top
to bottom panel. The outlier fraction is much lower in GRAHSP than the other
methods, and does not strongly increase with luminosity. Indeed, the outlier
fraction (top left panel in Fig. 24) approximately corresponds to the
probability expected outside the error bars. This is not the case for the Y19,
C15 and D14 CIGALE models, for which the outlier rates rise up to 20 per cent
at high luminosities. While the bias (centre left panel in Fig. 24) of the
D14, C15 and Y19 methods starts near zero in low-luminosity AGN, it rises
steeply with luminosities, up to 0.4 dex. This means that any sample-averaged
stellar masses are systematically raised by this amount.

The right column of Fig. 24 presents trends of the SFR estimator. The range
(right axis) of the bias panel had to be increased compared with the stellar
masses (left column), indicating how much more uncertain this measurement is.
While the scatter is comparable between methods, the bias and outlier fraction
are lowest in GRAHSP. The latter indicates that CIGALE-based reported SFR
error bars are too small, whereas GRAHSP produces realistic error bars.

###  6.4 AGN luminosity

The last comparison from the Chimera benchmark compares the retrieval of
luminosities. The top panel of Fig. 25 shows the retrieval of bolometric