11institutetext: INAF – Osservatorio Astronomico di Roma, Via Frascati 33,
00040 Monte Porzio Catone (RM), Italy.
11email<EMAIL_ADDRESS>22institutetext: Space Science Data Center,
Agenzia Spaziale Italiana, Via del Politecnico snc, 00133 Roma, Italy.
33institutetext: ASTRON – Netherlands Institute for Radio Astronomy, Oude
Hoogeveensedijk 4, 7991 PD Dwingeloo, The Netherlands. 44institutetext:
Dipartimento di Fisica, Università degli Studi di Roma “La Sapienza”, P.le A.
Moro 5, 00185 Roma, Italy. 55institutetext: INFN, Sezione di Roma I, P.le A.
Moro 2, 00185 Roma, Italy. 66institutetext: Dipartimento di Fisica, Università
di Roma “Tor Vergata”, Via della Ricerca Scientifica 1, 00133 Roma, Italy.
77institutetext: Department of Astronomy, University of Maryland, College
Park, MD 20742, USA. 88institutetext: NASA/Goddard Space Flight Center, Code
662, Greenbelt, MD 20771, USA.

# “Zombie” or active? An alternative explanation to the properties of star-
forming galaxies at high redshift

F. G. Saturni “Zombie” or active? An alternative explanation to the properties
of star-forming galaxies at high redshift“Zombie” or active? An alternative
explanation to the properties of star-forming galaxies at high
redshift“Zombie” or active? An alternative explanation to the properties of
star-forming galaxies at high redshift“Zombie” or active? An alternative
explanation to the properties of star-forming galaxies at high redshift M.
Mancini “Zombie” or active? An alternative explanation to the properties of
star-forming galaxies at high redshift“Zombie” or active? An alternative
explanation to the properties of star-forming galaxies at high redshift E.
Pezzulli “Zombie” or active? An alternative explanation to the properties of
star-forming galaxies at high redshift“Zombie” or active? An alternative
explanation to the properties of star-forming galaxies at high
redshift“Zombie” or active? An alternative explanation to the properties of
star-forming galaxies at high redshift“Zombie” or active? An alternative
explanation to the properties of star-forming galaxies at high
redshift“Zombie” or active? An alternative explanation to the properties of
star-forming galaxies at high redshift“Zombie” or active? An alternative
explanation to the properties of star-forming galaxies at high redshift F.
Tombesi “Zombie” or active? An alternative explanation to the properties of
star-forming galaxies at high redshift“Zombie” or active? An alternative
explanation to the properties of star-forming galaxies at high
redshift“Zombie” or active? An alternative explanation to the properties of
star-forming galaxies at high redshift“Zombie” or active? An alternative
explanation to the properties of star-forming galaxies at high
redshift“Zombie” or active? An alternative explanation to the properties of
star-forming galaxies at high redshift“Zombie” or active? An alternative
explanation to the properties of star-forming galaxies at high
redshift“Zombie” or active? An alternative explanation to the properties of
star-forming galaxies at high redshift“Zombie” or active? An alternative
explanation to the properties of star-forming galaxies at high redshift

(Received 2018 Apr 19; accepted 2018 Jun 18)

###### Abstract

Context. Star-forming galaxies at high redshift show anomalous values of
infrared excess, which can be described only by extremizing the existing
relations between the shape of their ultraviolet continuum emission and their
infrared-to-ultraviolet luminosity ratio, or by constructing ad hoc models of
star formation and dust distribution.

Aims. We present an alternative explanation, based on unveiled AGN activity,
of the existence of such galaxies. The scenario of a weak AGN lends itself
naturally to explain the observed spectral properties of these high-z objects
in terms of a continuum slope distribution and not altered infrared excesses.

Methods. To this end, we directly compare the infrared-to-ultraviolet
properties of high-redshift galaxies to those of known categories of AGN
(quasars and Seyferts). We also infer the characteristics of their possible
X-ray emission.

Results. We find a strong similarity between the spectral shapes and
luminosity ratios of AGN with the corresponding properties of such galaxies.
In addition, we derive expected X-ray fluxes that are compatible with the
energetics from AGN activity.

Conclusions. We conclude that a moderate AGN contribution to the UV emission
of such high-$z$ objects is a valid alternative to explain their spectral
properties. Even the presence of an active nucleus in each source would not
violate the expected quasar statistics. Furthermore, we suggest that the
observed similarities between anomalous star-forming galaxies and quasars may
provide a benchmark for future theoretical and observational studies on the
galaxy population in the early Universe.

###### Key Words.:

galaxies: ISM – infrared: galaxies – ultraviolet: galaxies – quasars: general
– infrared: quasars – ultraviolet: quasars

##  1 Introduction

The estimate of the galactic dust content and properties in objects at high
redshift ($z>5$) is crucial for assessing the details of the star formation
history in galaxies. However, it is not always feasible to allow the
measurement of the far-infrared (FIR) flux emission, due to low instrumental
sensitivity or lack of spectral coverage (e.g., Bouwens et al. 2014; Rogers et
al. 2014; Finkelstein et al. 2015a, b; Bowler et al. 2015; Capak et al. 2015;
Bouwens et al. 2016; Barišić et al. 2017, and references therein). Under the
assumption of radiative equilibrium, it is possible to correlate the FIR
emission to the amount of ultraviolet (UV) flux absorbed by dust and
reprocessed at longer wavelengths since the dust extinction of reasonable
absorption models peaks in the UV band (Weingartner & Draine 2001; Bianchi &
Schneider 2007; Draine 2011). Therefore, it is common practice to quantify the
extinction at such redshifts using the galaxy-emission properties in the rest-
frame UV band (Bouwens et al. 2014, 2015; Finkelstein et al. 2015a).

ID | R.A. (h m s) | dec (d m s) | $z$ |  $L_{\rm FIR}$ $(L_{\odot})$ |  $L_{\rm UV}$ $(L_{\odot})$ | $\beta$ | IRX
---|---|---|---|---|---|---|---
HZ1 | 09 59 53.25 | +02 07 05.43 | 5.69 |  $<$10.32 | $11.210\pm 0.010$ | $-1.92^{+0.14}_{-0.11}$ | $<-0.036$
HZ2 | 10 02 04.10 | +01 55 44.05 | 5.67 |  $<$10.30 | $11.150\pm 0.010$ | $-1.82\pm 0.10$ | $<-0.034$
HZ3 | 10 00 09.43 | +02 20 13.86 | 5.55 |  $<$10.53 | $11.080\pm 0.010$ | $-1.72^{+0.12}_{-0.15}$ | $<-0.022$
HZ4 | 09 58 28.52 | +02 03 06.74 | 5.54 | $11.13\pm 0.54$ | $11.280\pm 0.010$ | $-2.06^{+0.13}_{-0.15}$ | $-0.006\pm 0.292$
HZ5 | 10 00 51.60 | +02 34 57.55 | 5.31 |  $<$10.30 | $11.450\pm 0.010$ | $-1.01^{+0.06}_{-0.12}$ | $<-0.046$
HZ5a | 10 00 51.52 | +02 34 59.29 | 5.31 |  $<$10.30 |  $<$10.370 | — | —
HZ6 | 10 00 21.50 | +02 35 11.08 | 5.29 | $11.13\pm 0.23$ | $11.470\pm 0.100$ | $-1.14^{+0.12}_{-0.14}$ | $-0.013\pm 0.063$
HZ6a | — | — | 5.29 | $10.26\pm 0.23$ | $11.110\pm 0.070$ | $-0.59^{+1.05}_{-1.12}$ | $-0.035\pm 0.058$
HZ6b | — | — | 5.29 | $10.87\pm 0.23$ | $11.000\pm 0.070$ | $-1.50^{+1.05}_{-1.22}$ | $-0.005\pm 0.058$
HZ6c | — | — | 5.29 | $10.79\pm 0.23$ | $10.810\pm 0.070$ | $-1.30^{+0.51}_{-0.37}$ | $-0.001\pm 0.058$
HZ7 | 09 59 30.48 | +02 08 02.81 | 5.25 |  $<$10.35 | $11.050\pm 0.020$ | $-1.39^{+0.15}_{-0.17}$ | $<-0.028$
HZ8 | 10 00 04.06 | +02 37 35.81 | 5.15 |  $<$10.26 | $11.040\pm 0.020$ | $-1.42^{+0.19}_{-0.18}$ | $<-0.032$
HZ8W | 10 00 03.97 | +02 37 36.23 | 5.15 |  $<$10.26 | $10.570\pm 0.040$ | $-0.10\pm 0.29$ | $-0.013\pm 0.002$
HZ9 | 09 59 51.70 | +02 22 42.16 | 5.55 | $11.54\pm 0.19$ | $10.950\pm 0.020$ | $-1.59^{+0.22}_{-0.23}$ | $0.023\pm 0.036$
HZ10 | 10 00 59.30 | +01 33 19.53 | 5.66 | $11.94\pm 0.08$ | $11.140\pm 0.020$ | $-1.92^{+0.24}_{-0.17}$ | $0.030\pm 0.007$
HZ10W | — | — | 5.66 | $11.64\pm 0.08$ | $10.230\pm 0.050$ | $-1.47^{+0.77}_{-0.44}$ | $0.056\pm 0.009$

Table 1: Basic properties of the high-$z$ galaxies discovered by Capak et al.
(2015) and reanalyzed by Barišić et al. (2017). The objects are ordered by ID
number. Values preceded by “$<$” represent an upper limit in the corresponding
quantity, while cells with “—” denote missing data from the original
publications.

The recent detection of FIR emission from high-$z$ galaxies ($z>5$;
Finkelstein et al. 2015a, b; Bowler et al. 2015; Capak et al. 2015; Waters et
al. 2016; Bouwens et al. 2016; Knudsen et al. 2017) provides better
constraints on their dust properties, such as extinction. In particular, the
study of galaxies with non-negligible star formation rates (SFRs) shows that
their intrinsic wavelength spectral index $\beta_{\rm int}$, i.e., the slope
of the power-law spectrum used to describe the galactic UV emission, is nearly
constant with values between $-2.5$ and $-2.3$ (Meurer et al. 1999; Wilkins et
al. 2013; Bouwens et al. 2014; Reddy et al. 2015; Forrest et al. 2016; Mancini
et al. 2016; Narayanan et al. 2018; Reddy et al. 2018; Popping et al. 2017;
Cullen et al. 2017), although Reddy et al. (2015), Casey et al. (2014) and
Forrest et al. (2016) found a value compatible with $\beta>-2.3$ and Talia et
al. (2015) estimated $\beta\sim-3.02$.

Meurer et al. (1999) initially suggested the possibility of correlating the
difference between the observed and the expected $\beta$, which is used to
describe the galactic UV continuum, to the extinction coefficient at 1600 Å
(hereafter $A_{1600}$). This procedure is made possible by the independence of
the intrinsic spectrum of star-forming galaxies from the stellar-population
age and by its weak correlation with the stellar metallicity (see, e.g.,
Meurer et al. 1999; Mancini et al. 2015, and references therein). Furthermore,
it is possible to correlate $A_{1600}$ with the infrared excess (IRX) defined
as

${\rm IRX}=\log{\left(\frac{L_{\rm FIR}}{L_{\rm UV}}\right)},$ (1)

where $L_{\rm FIR}$ is the FIR galactic luminosity integrated from 8 $\mu$m to
1000 $\mu$m and $L_{\rm UV}$ is the monochromatic luminosity at 1600 Å in the
rest frame. The IRX-to-$\beta$ relation obtained in this way only depends on
the dust mass (i.e., the normalization of the dust-extinction law), the dust
properties (i.e., the shape of such a law), and the intrinsic spectral index
$\beta_{\rm int}$ corrected for the dust absorption.

This relation appears to be universal for star-forming galaxies (e.g., Meurer
et al. 1999). However, Capak et al. (2015) observed a population of high-$z$
star-forming galaxies that appear to strongly deviate from the usual
relations, showing lower IRX with respect to that expected on the basis of
their observed spectral slope. These sources were selected to measure [C II]
gas and dust emission in Lyman-break galaxies (LBGs) detected in the 2 deg2
Cosmological Evolution Survey (COSMOS; Scoville et al. 2007) in order to
explore the average properties of the interstellar medium (ISM) at high
redshift. The reanalysis of the Capak et al. (2015) data done by Barišić et
al. (2017) partially relaxed the tension between their IRXs and spectral
slopes, but without eliminating it, especially for the case of upper limits to
the IRX. These objects populate a region of the IRX–$\beta$ plane that is
typical of quiescent galaxies with negligible or no star formation (“dead”)
despite having SFRs from tens to hundreds of $M_{\odot}$ yr-1 as measured by
their [C II] emission (Capak et al. 2015), and could therefore be denoted as
“zombie” galaxies.

Several scenarios have been proposed to explain how star-forming galaxies can
deviate from the usual IRX-to-$\beta$ relations (Bouwens et al. 2016; Ferrara
et al. 2017; see also Ouchi et al. 1999). However, the observed discrepancy
with the predicted IRXs would either imply a strong bias in the determination
of the dust temperature (Bouwens et al. 2016) or extreme physical conditions
in the dust clouds (Ferrara et al. 2017). In this work, we formulate the
alternative hypothesis that the high-$z$ galaxies with anomalous IRX-
to-$\beta$ discovered by Capak et al. (2015) are actually hosting AGN
activity. This scenario allows us to consider objects in which the UV flux is
dominated by accretion-disk processes (e.g., Salpeter 1964; Zel’dovich &
Novikov 1965; Shakura & Sunyaev 1973) and in which the intrinsic UV spectral
slope depends on the properties of the central engine (accretion rate, black
hole mass, black hole spin) and is therefore not limited to steep values
(e.g., Shankar et al. 2016). In this way, the Capak et al. (2015) objects are
no longer constrained to have anomalous dust properties, but are instead
located in the region of the IRX–$\beta$ plane corresponding to their AGN UV
continuum.

Figure 1: IRX–$\beta$ plane for quasars. The data corresponding to the
high-$z$ galaxies by Capak et al. (2015) reanalyzed by Barišić et al. (2017)
are shown (black squares), along with their 1$\sigma$ errors, superimposed on
the HerS/SDSS quasars (filled dots). The empirical IRX-to-$\beta$ relations by
Pettini et al. (1998, dashed lines) and Meurer et al. (1999, solid lines)
computed for values of $\beta_{\rm int}$ from $-2.5$ to $-1$ are also
reported. Left panel: HerS/SDSS quasars segregated in logarithmic $M_{\rm
BH}$. Right panel: HerS/SDSS quasars segregated in logarithmic $M_{\rm gas}$.
Only HerS/SDSS objects with measured $M_{\rm BH}$ and $M_{\rm gas}$ are used.

The paper is organized as follows. In Sect. 2 we exploit the idea developed by
Capak et al. (2015) that AGN activity is responsible for the altered IRX-
to-$\beta$ of the high-$z$ galaxies. We present a tentative statistical
classification of these objects based on their spectral properties in Sect. 3.
We also compute the expected X-ray emission from AGN with properties similar
to these sources in Sect. 4, and compare the fluxes estimated in this way with
the current X-ray upper limits of these sources derived from Chandra
observations of the COSMOS field. Finally, we summarize and discuss our
findings in Sect. 5. Throughout the text, we use “AGN” to indicate both low-
and high-luminosity objects, and adopt a concordance cosmology with $H_{0}=70$
km s-1 Mpc-1, $\Omega_{\rm M}=0.3$, and $\Omega_{\Lambda}=0.7$.

##  2 The IRX-to-$\beta$ relation for AGN

The sample of 16 high-$z$ galaxies discovered by Capak et al. (2015) is
presented in Table 1. It comprises six single sources (HZ1, HZ2, HZ3, HZ4,
HZ7, and HZ9), and four multiple systems (HZ5, HZ6, HZ8, and HZ10) that are
probably composed of interacting objects. Their redshifts were measured
through absorption-line spectroscopy with the instrument DEIMOS at the Keck
Observatory. For all but one of them (HZ5a), the measurements or at the least
upper limits of their UV spectral slope and IRX are given in Barišić et al.
(2017). In the following, we adopt this data set for our analysis of such
objects in search of AGN behavior.

In order to validate the hypothesis that star-forming galaxies at high
redshift with anomalous IRX-to-$\beta$ values are hosts of AGN activity, we
first define a comparison sample of AGN that has measurements of UV to optical
spectral slopes $\beta$, IR luminosities $L_{\rm IR}$, and UV luminosities
$L_{\rm UV}$. We find that Dong & Wu (2016) constructed such a catalog for
quasars, i.e., AGN with $L_{\rm bol}\gtrsim 10^{45}$ erg s-1. A similar
catalog has been constructed by García-González et al. (2016) for nearby
($D_{\rm L}<70$ Mpc) Seyfert galaxies ($L_{\rm bol}\lesssim 10^{45}$ erg s-1).
By comparing both catalogs to the Capak et al. (2015) galaxies, we are able to
span five orders of magnitude in bolometric luminosity from $10^{43}$ to
$10^{48}$ erg s-1, thus covering a wide range of AGN energetics.

We note that the selection criteria used in these AGN catalogs significantly
differ from those adopted by Capak et al. (2015) and Barišić et al. (2017),
due to the diverse purposes for which the samples were constructed (the study
of the link between AGN activity and host-galaxy star formation in the case of
Dong & Wu 2016, and the disentanglement of the FIR emission due to AGN
activity from that due to star formation in García-González et al. 2016).
Nevertheless, a comparison between AGN and the Capak et al. (2015) objects
based on the average properties of such samples is still feasible because we
focus on spectral properties that are nearly unhindered by the nature of the
ISM gas. For instance, the zombie galaxy HZ5, though satisfying the object
selection performed by Capak et al. (2015), is actually a weak Type 1 quasar
(Brightman et al. 2014; Kalfountzou et al. 2014). Therefore, we expect that
the different selection criteria among the three samples do not greatly affect
our findings.

###  2.1 Construction of a quasar comparison sample

The quasar catalog by Dong & Wu (2016) is composed of objects from the Sloan
Digital Sky Survey Data Release 7 (SDSS-DR7; Abazajian et al. 2009; Schneider
et al. 2010; Shen et al. 2011) and Data Release 10 (SDSS-DR10; Ahn et al.
2014; Pâris et al. 2014) that were observed in the FIR domain by the Spectral
and Photometric Imaging Receiver (SPIRE) on board the Herschel Space
Observatory in the framework of the Herschel Stripe 82 Survey (HerS; Viero et
al. 2014). This survey covers $\sim$29% (79 deg2) of the SDSS Stripe 82 (270
deg2 on the celestial equator in the southern Galactic cap; Adelman-McCarthy
et al. 2007) with FIR observations at 250, 350, and 500 $\mu$m.

The joint HerS/SDSS sample includes 207 quasars with the following
characteristics:

  * •

UV to optical ugriz magnitudes from the SDSS;

  * •

near-IR (NIR) $JHK_{\rm s}$ magnitudes from the 2MASS (Skrutskie et al. 2006)
and UKIDSS (Lawrence et al. 2007) surveys;

  * •

mid-IR (MIR) fluxes from the AllWISE catalog (Wright et al. 2010; Mainzer et
al. 2011) ranging from $3.4$ to 22 $\mu$m;

  * •

Herschel FIR flux measurements in the passbands centered at the wavelengths
mentioned above.

The spectral energy distribution (SED) of each object has been decomposed in
several spectral components, each accounting for the emission from a different
source: a power law for the quasar accretion disk, two near- to mid-IR bumps
peaking at around $2-4$ $\mu$m and $10-20$ $\mu$m for the radiation
reprocessed by the dusty torus, two galaxy templates (S$b$ for the young
stellar component and E for the old one) for the host-galaxy contribution, a
Small Magellanic Cloud (SMC) reddening law and a “graybody” radiation with
power-law emissivity for the FIR excess. In this way, Dong & Wu (2016)
measured IR and UV luminosities for each quasar in the sample, along with the
corresponding frequency spectral slope $\alpha$ which is related to the
wavelength slope $\beta$ by

$\beta=-(2+\alpha)$ (2)

Thanks to the SED analysis, the HerS/SDSS catalog also provides information on
several other quasar and host-galaxy parameters, such as the virial black hole
(BH) mass $M_{\rm BH}$ derived from single-epoch relations (Vestergaard &
Peterson 2006; Vestergaard & Osmer 2009; Shen & Kelly 2010) and the mass of
cold gas $M_{\rm gas}$ stored in the AGN surrounding environment. The former
is measured by adopting the quasar continuum luminosities computed from the
power-law emission, while the latter is estimated from the long-wavelength
tail of the dust emission (see Scoville et al. 2014, for details).

Since the SED decomposition in multiple components allows us to measure quasar
intrinsic reddening, this may be used to correct the quantities derived from
the fitting procedure. Therefore, $L_{\rm UV}$ and $\beta$ given in the
HerS/SDSS sample are intrinsic to the quasars, whereas the corresponding
observed quantities (i.e., with no correction for internal reddening applied)
are needed in the construction of the IRX-to-$\beta$ relation. In order to
recover the observed $L_{\rm UV}$ and $\beta$, we apply backward the de-
reddening procedure adopted by Dong & Wu (2016). We individually de-correct
both $L_{\rm UV}$ and $\beta$ through a SMC extinction template (Pei 1992)
scaled to the value of intrinsic reddening $A_{\rm int}$ of each quasar, and
use these reddened quantities in the evaluation of the quasar IRX-to-$\beta$
relation.

The resulting sample of quasar IRX and spectral slopes is shown in Fig. 1. At
first glance, it is evident as all but three high-$z$ galaxies by Capak et al.
(2015) fall inside the region of the IRX–$\beta$ diagram populated by quasars.
For comparison, we compute the IRX-to-$\beta$ relations by Pettini et al.
(1998) and Meurer et al. (1999) for several values of the intrinsic spectral
slope $\beta_{\rm int}$ from $-2.5$ to $-1$, and show them superimposed on the
observed IRX-to-$\beta$ values for high-$z$ galaxies and quasars.

It should be noted that the HerS/SDSS quasars tend to segregate in the
IRX–$\beta$ plane according to both their $M_{\rm BH}$ (left panel of Fig. 1)
and $M_{\rm gas}$ (right panel of Fig. 1). In particular, quasars hosting more
massive black holes have lower IRXs, whereas objects with higher gas masses
show steeper intrinsic slopes coupled to slightly higher IRXs. It is unclear
whether the tendency for quasars with higher $M_{\rm BH}$ to have lower IRXs
is due to a physical mechanism such as quasar feedback (e.g., Cattaneo et al.
2009; Fabian 2012), since the main driver of feedback appears to be the quasar
Eddington ratio $L_{\rm bol}/L_{\rm edd}$ rather than the black hole mass
(Ricci et al. 2017b; see also Ishibashi & Fabian 2016; Ishibashi et al. 2017).
Such a segregation could equally be produced by a selection bias in the
observation of quasars with massive black holes, for example by preferentially
finding objects with bluer spectra at the high end of the $M_{\rm BH}$
distribution due to their lower Eddington ratio.

Figure 2: IRX–$\beta$ plane for Seyfert galaxies. The data by Barišić et al.
(2017) are shown (black squares), along with their 1$\sigma$ errors,
superimposed on the sample of Seyfert galaxies by García-González et al.
(2016, magenta crosses). For comparison, the theoretical IRX-to-$\beta$
relation by Meurer et al. (1999) with $\beta_{\rm int}=-2.23$ (solid line) and
the relation by Pettini et al. (1998) with $\beta_{\rm int}=-2$ (dashed line)
are also reported.

Conversely, a high $M_{\rm gas}$ may be associated with gas-rich quasar
surrounding environments, resulting in high-luminosity objects which tend to
have bluer spectra and hence steeper $\beta_{\rm int}$ (e.g., Shankar et al.
2016). Since high gas masses also enhance the star formation and dust content
of the host galaxy (e.g., Kennicutt & Evans 2012), the tendency of quasars
with higher $M_{\rm gas}$ to segregate towards steeper UV spectral slopes and
higher IRXs is explained by the simultaneous presence of a higher AGN and
stellar luminosity, and a higher dust mass. We note that future observations
of high-$z$ galaxies could benefit from using these tendencies as diagnostics
to infer the presence of a centrally located AGN powering part of the UV
emission.

###  2.2 Construction of a Seyfert comparison sample

The Seyfert comparison sample has been constructed in García-González et al.
(2016) by selecting objects from the Revised Shapley – Ames (RSA) catalog
(Sandage & Tammann 1987) that were identified as Seyfert galaxies by Maiolino
& Rieke (1995). These objects were selected in order to have
Herschel/Photometric Array Camera (PACS) imaging in at least two bands and
SPIRE photometric observations. The final catalog by García-González et al.
(2016) comprises 33 galaxies (15 Sy 1 and 18 Sy 2), of which we only use those
with near- and far-UV (NUV and FUV) magnitudes available in the literature in
order to compute UV spectral slopes and luminosities. We therefore select
Seyfert galaxies from this catalog by performing a search on the NASA
Extragalactic Database111Available at ned.ipac.caltech.edu. (NED) for objects
with both NUV and FUV measurements.

Our final sample is composed of 19 Seyferts with

  * •

NUV and FUV AB magnitudes measured by the Galaxy Evolution Explorer (GALEX)
satellite for UV surveys (Martin et al. 2005; Gil de Paz et al. 2007);

  * •

MIR spectroscopy with high angular resolution (see García-González et al.
2016, and refs. therein);

  * •

Herschel FIR photometric measurements from 70 $\mu$m to 500 $\mu$m.

All these objects have IR SEDs fitted by García-González et al. (2016) with
graybody components which allow us to compute $L_{\rm FIR}$. Coupling these
luminosities to the UV magnitudes, we are able to compute the IRX for this
subsample, which we show in Fig. 2 as a function of the corresponding $\beta$
compared to the high-$z$ galaxies by Capak et al. (2015). A visual inspection
reveals that these Seyferts have shallower UV spectral slopes and higher IRXs
with respect to quasars, thus populating the region at
$-1\lesssim\beta\lesssim 0$ and ${\rm IRX}\gtrsim 0$ which contains neither
high-$z$ galaxies nor HerS/SDSS quasars (see Fig. 1). This suggests that a
possible AGN activity in the Capak et al. (2015) objects must have properties
closer to efficient high-luminosity AGN than to weak ones.

Figure 3: Left panel: theoretical IRX-to-$\beta$ relations for mixed AGN and
host-galaxy UV contributions. A $\beta_{\rm int}^{\rm(q)}=-2$ AGN dominance of
30% (solid line) to 70% (dashed line) has been adopted. For comparison, the
theoretical IRX-to-$\beta$ relation by Meurer et al. (1999) with
$\beta^{\rm(g)}_{\rm int}=-2.23$ used to evaluate the UV host-galaxy
contribution (dot-dashed line) and the relation by Pettini et al. (1998) for
the UV AGN contribution with $\beta^{\rm(q)}_{\rm int}=-2$ (dotted line) are
also shown superimposed to the data sets by Barišić et al. (2017, black
squares), Dong & Wu (2016, cyan dots) and García-González et al. (2016,
magenta crosses). Right panel: least-squares fit of Eq. 8 (solid line) to the
high-$z$ galaxy data by Barišić et al. (2017) and the AGN samples by Dong & Wu
(2016) and García-González et al. (2016), along with its boundaries at
$1\sigma$ (dashed lines) and $2\sigma$ level (dotted lines).

###  2.3 Estimate of the galactic UV contribution

The scenario in which the high-$z$ galaxies by Capak et al. (2015) host AGN
activity with average spectral properties (i.e., $\beta_{\rm int}\sim-2$ once
a correction for the dust reddening is taken into account in evaluating the
spectral slope; Dong & Wu 2016) holds well in describing the observed IRX-
to-$\beta$ distribution of such objects. However, a residual galaxy
contribution to the total UV flux cannot be neglected, since models of AGN and
host-galaxy coupling, derived from numerical simulations (e.g., Schneider et
al. 2015), show that the AGN emission accounts on average for 30% to 70% of
the system’s UV energetics. Therefore, we tried to reproduce the observed
values of IRX and $\beta$ for the high-$z$ galaxies with a hybrid relation
that takes into account both the AGN and host-galaxy UV contribution.

Equation 1 can be rewritten in a form that only depends on the fraction
$D(\tau_{\lambda},\beta)$ of reprocessed UV light $L_{1600}$ into IR radiation
and on the UV dust reddening $S(\tau_{\lambda},\beta)$. In doing this, we
obtain $L_{\rm UV}=L_{1600}\cdot S(\tau_{\lambda},\beta)$, $L_{\rm
FIR}=L_{1600}\cdot D(\tau_{\lambda},\beta)$ and

${\rm IRX}=\log{\left[\frac{L_{\rm 1600}\cdot D(\tau_{\lambda},\beta)}{L_{\rm
1600}\cdot
S(\tau_{\lambda},\beta)}\right]}=\log{\left[\frac{D(\tau_{\lambda},\beta)}{S(\tau_{\lambda},\beta)}\right]}$
(3)

The reprocessed UV fraction $D(\tau_{\lambda},\beta)$ and the reddening
$S(\tau_{\lambda},\beta)$ are respectively given by

$D(\tau_{\lambda},\beta)=\int_{1500}^{3000}\left(\frac{\lambda}{1600\mbox{
\AA}}\right)^{\beta}\left(1-e^{-\tau_{\lambda}}\right)d\lambda$ (4)

and

$S(\tau_{\lambda},\beta)=\int_{1500}^{3000}\left(\frac{\lambda}{1600\mbox{
\AA}}\right)^{\beta}e^{-\tau_{\lambda}}d\lambda,$ (5)

where $\beta$ is the intrinsic slope of the UV spectrum and $\tau_{\lambda}$
is the wavelength-dependent optical depth of the dust distribution.

Since the AGN (through accretion-disk emission) and the host galaxy (through
star formation) both heat the dust that reprocesses the UV into FIR emission,
we can separate the contributions to the IRX by writing

${\rm IRX}=\log{\left\\{\frac{L^{\rm(q)}_{\rm 1600}\cdot
D\left[\tau^{\rm(q)}_{\lambda},\beta_{\rm q}\right]+L^{\rm(g)}_{\rm 1600}\cdot
D\left[\tau^{\rm(g)}_{\lambda},\beta_{\rm g}\right]}{L^{\rm(q)}_{\rm
1600}\cdot S\left[\tau^{\rm(q)}_{\lambda},\beta_{\rm q}\right]+L^{\rm(g)}_{\rm
1600}\cdot S\left[\tau^{\rm(g)}_{\lambda},\beta_{\rm g}\right]}\right\\}},$
(6)

where $\beta_{\rm q}$ and $\beta_{\rm g}$ are the AGN and host-galaxy UV
spectral slopes respectively. Defining $f_{\rm q}$ as the fraction of energy
injected in the dust by the AGN, we get

$f_{\rm q}=\frac{L^{\rm(q)}_{\rm 1600}\cdot
S\left[\tau^{\rm(q)}_{\lambda},\beta_{\rm q}\right]}{L^{\rm(q)}_{\rm
1600}\cdot S\left[\tau^{\rm(q)}_{\lambda},\beta_{\rm q}\right]+L^{\rm(g)}_{\rm
1600}\cdot S\left[\tau^{\rm(g)}_{\lambda},\beta_{\rm g}\right]}$ (7)

and finally, combining Eqs. 3, 6, and 7,

${\rm IRX}=\log{\left[f_{\rm q}\cdot 10^{{\rm IRX}_{\rm q}}+\left(1-f_{\rm
q}\right)\cdot 10^{{\rm IRX}_{\rm g}}\right]}$ (8)

We adopt the functional form by Pettini et al. (1998) with $\beta_{\rm
int}=-2$ in order to describe the ${\rm IRX}_{\rm q}$. This model is in fact
valid if the dust distribution is foreground to the UV source (e.g., Gordon et
al. 1997; Inoue 2005; Gallerani et al. 2010), which is suggested for AGN by
unified scenarios where a clumpy distribution of obscuring dust surrounds the
central engine (e.g., Urry & Padovani 1995; Risaliti et al. 2011; Krumpe et
al. 2014). For ${\rm IRX}_{\rm g}$, we adopt the Meurer et al. (1999) IRX-
to-$\beta$ relation with $\beta_{\rm int}=-2.23$, since in galaxies the dust
is expected to be mixed with stars rather than being (almost) completely
foreground (Calzetti 1997; Pettini et al. 1998; Meurer et al. 1999). In order
to cover the whole average range of AGN contributions to the total UV
emission, we choose for $f_{\rm q}$ the extreme values of 0.3 and 0.7
(Schneider et al. 2015).

Figure 3 shows such combined IRX-to-$\beta$ relations, along with the pure
relations by Pettini et al. (1998) and Meurer et al. (1999). The relations
with 30% to 70% UV contribution from AGN emission do not fully describe any of
the data by Barišić et al. (2017), while the Meurer et al. (1999) relation
partially intercepts the observations at IRX $\gtrsim 0$ and the Pettini et
al. (1998) relation partially succeeds in describing the bulk of observations
at $-1<$ IRX $<0$ and $-2<\beta<-1$. Overall, an IRX-to-$\beta$ relation with
both an AGN and galactic contribution to the UV emission appears to be the
most suitable to simultaneously describe the spectral properties of Seyferts,
quasars, and star-forming galaxies at high redshift. The posterior
distribution of least-squares fits of Eq. 8 performed over bootstrap
realizations of the total data set by Barišić et al. (2017, high-$z$ galaxies
with anomalous IRX), Dong & Wu (2016, HerS/SDSS quasars), and García-González
et al. (2016, local Seyferts) with $f_{\rm q}$ as the only free parameter
gives $f_{\rm q}=0.73\pm 0.11$, compared with a value of $f_{\rm q}=0.70\pm
0.12$ obtained by fitting only the AGN-dominated objects.

###  2.4 Intrinsic reddening estimates of the high-$z$ galaxies

Under the assumption that the Capak et al. (2015) galaxies are actual
observations of AGN hosts that exactly follow an IRX-to-$\beta$ relation with
AGN contribution, it is now possible to derive the unattenuated AGN spectral
slopes $\beta_{\rm int}^{\rm(q)}$ using the measured IRX and the observed UV
slopes $\beta_{\rm obs}$. In order to obtain these values, we evaluate Eq. 8
for each high-$z$ galaxy with $f_{\rm q}$ between 0.3 and 1, fixing the
intrinsic slope for the galactic contribution $\beta_{\rm int}^{\rm(g)}$ at
$-2.23$ (Meurer et al. 1999) and adopting as our fiducial estimate the value
of $\beta_{\rm int}^{\rm(q)}$ for which the composite IRX-to-$\beta$ relation
intercepts each data point. In particular, we find that the lowest AGN
contribution (i.e., $f_{\rm q}=0.3$) is the most reasonable choice to
reproduce the observed IRX-to-$\beta$ values for the four leftmost objects in
the IRX-to-$\beta$ diagram, namely HZ4, HZ9, HZ10, and HZ10W. This choice of
$f_{\rm q}$ predicts values of $\beta_{\rm int}^{\rm(q)}$ between $-5.7$ and
$-2.8$ for such objects; higher values of the AGN contribution would predict
even steeper intrinsic AGN UV slopes or would fail to reproduce the observed
data.

The remaining objects are instead well reproduced if we set a complete AGN
domination in the UV flux (i.e., $f_{\rm q}=1$). Lower values of $f_{\rm q}$
lead the composite IRX-to-$\beta$ relation to intercept the region between the
IRX-to-$\beta$ fiducial values below $\beta_{\rm int}^{\rm(q)}\sim-2$. An AGN-
dominated UV spectrum allows us instead to reproduce the data of HZ1, HZ2,
HZ3, HZ6b, and HZ6c with a single relation at $\beta_{\rm int}^{\rm(q)}=-2$;
we recall that this is the average value of the intrinsic UV slopes for AGN
(Dong & Wu 2016). Furthermore, such a relation also describes the remaining
objects assuming $-2\lesssim\beta_{\rm int}^{\rm(q)}\lesssim 0$.

We then use the $\beta_{\rm int}$ derived in this way along with the HST/WFC3
data by Barišić et al. (2017) in order to obtain an estimate of the UV
extinction $A_{1600}$ affecting the AGN emission. To do this, we adopt the
simplified assumption that the bolometric luminosities of these objects are
approximated using the optical/UV bolometric corrections by Runnoe et al.
(2012b). Therefore, we estimate the optical monochromatic luminosities at 3000
Å from the F160W magnitudes by Barišić et al. (2017), and apply the
corrections listed in Table 1 of Runnoe et al. (2012a) to derive the de-
reddened UV luminosities at 1600 Å. The choice of indirectly computing such
luminosities through the correction of the optical magnitudes is motivated by
the fact that, at least for low to moderate values of $A_{1600}$, the
luminosity at optical wavelengths is less affected by extinction than the UV
flux.

Figure 4: Extinction coefficient $A_{1600}$ as a function of the intrinsic
spectral slope $\beta_{\rm int}^{\rm(q)}$ for the high-$z$ galaxies by Capak
et al. (2015). The average extinction (dot-dashed line) is indicated along
with its 1$\sigma$ confidence interval (gray band).

Under these hypotheses, we define the values of $A_{1600}$ by assuming that
the intrinsic UV slopes of the objects by Capak et al. (2015) are reddened by
an SMC-like extinction. Such extinction coefficients are shown in Fig. 4 as a
function of $\beta_{\rm int}^{\rm(q)}$, suggesting that the possible AGN
emission of these high-$z$ galaxies is affected on average by low to moderate
reddening along the line of sight ($A_{1600}\lesssim 1.2$ mag).

###  2.5 Possible issues in the comparison between different samples of
extragalactic sources

Although the presented analysis focuses on the comparison among the average
emission properties of the sample of high-$z$ star-forming galaxies by Capak
et al. (2015) with those of two AGN classes, the dependence of the latter on
quantities derived from SED fitting can bias the results due to (i) the
availability of different models for the paramount spectral components , (ii)
the possible parameter degeneracies arising in the SED construction, and (iii)
the adoption of peculiar extinction curves. In the following, we briefly
describe the expected impact of these biases on the main findings reported in
this section.

#### Model dependencies.

Model dependencies are particularly important in fitting the main sources of
UV (e.g., the details of the AGN process and its ionized surrounding
environment) and IR emission (e.g., the reprocessing of radiation by dust in
the AGN’s surrounding environment). For instance, the AGN emission may be
reproduced by either using a UV template to take into account the additional
light from resolved emission lines and Fe II blending (García-González et al.
2016) or adopting a simple power-law functional form (Dong & Wu 2016).
However, the difference in such approaches generates discrepancies of
$\lesssim$10% only in the determination of the UV flux from AGN, i.e., the
average contribution of emission features to the total emission with respect
to the continuum (e.g., Vanden Berk et al. 2001; Glikman et al. 2006), and
hence is not expected to have a great impact on the derivation of the source’s
properties unless tailoring the selection for objects with extreme line
emission. On broader wavelength ranges, the differences introduced by SEDs
related to different emission models can directly affect the determination of
IRX, UV slope, and dust extinction. Duras et al. (2017) show that there are
some differences among quasar SEDs, but these are in general compatible with
each other within at most 0.5 dex on a wavelength range of 100 to $10^{6}$ Å
(see, e.g., their figures 2 and 3), apart from some peculiar choices (e.g.,
Stalevski et al. 2016). Additionally, the largest discrepancies in the IR
emission are due to the opening angle of a dusty AGN torus (Stalevski et al.
2016; Duras et al. 2017), but they are again limited to $\sim$0.5 dex. A shift
of $\sim\pm 1$ dex in the IRX of AGN introduced by SED differences is still
compatible with the observed IRXs of the zombie galaxies by Capak et al.
(2015), thus not greatly biasing our main results.

#### Degeneracies.

Another problem is the degeneracies among the free parameters involved in the
SED-fitting procedure, which can be explored in detail only with the
availability of rich data sets in order to give statistical significance to
the fit results (see, e.g., Mullaney et al. 2011). For instance, the UV SED of
a quasar affected by intrinsic reddening could be reproduced by either
assuming a shallow UV slope (thus neglecting the dust extinction) or an
appropriate extinction law applied to a steep spectrum ($\beta\sim-2$; e.g.,
Dong & Wu 2016). In the FIR, both the dust emission from the host galaxy and
the emission from a possible old stellar population peak at $\sim$100 $\mu$m,
thus making it difficult to disentangle quasar and galaxy light. Dong & Wu
(2016) overcome these issues by adopting quasar-galaxy mixing diagrams (Hao et
al. 2013) that allow us to distinguish AGN-dominated sources from galaxy-
dominated and reddened SEDs. Additionally, they also fix some SED parameters
to typical literature values, such as the AGN torus temperature at 1300 K and
the slope of the graybody emissivity for galactic dust at 1.6 in order to
reduce possible secondary degeneracies. On the other hand, García-González et
al. (2016) rely on spectro-photometric indicators such as the
$f_{\nu}(\mbox{70 }\mu{\rm m})/f_{\nu}(\mbox{160 }\mu{\rm m})$ flux ratio and
the equivalent width of polycyclic aromatic hydrocarbon (PAH) features to
discriminate between galaxy emission and AGN activity. Therefore, we can
assume that parameter degeneracies are already taken into account in the
adopted AGN samples, thus not affecting their comparison with the properties
of the zombie galaxies.

#### Dust extinction.

Regarding the choice of the extinction curve, we have already outlined that
taking into account different extinction curves and dust distribution
geometries do not help reconcile the expected behavior of the zombie galaxies
in the IRX–$\beta$ plane, as widely described in the literature (see, e.g.,
Mancini et al. 2016; Narayanan et al. 2018; Ferrara et al. 2017; Behrens et
al. 2018, and references therein). In particular, the work by Narayanan et al.
(2018) suggests that in order to have a galaxy that covers the same locus of
the zombie galaxies these object would need a spectral emission dominated by
an old stellar population. Alternatively, Ferrara et al. (2017) suggest that
it is possible to reduce the discrepancy between the observed and expected
IRXs by considering a stellar population completely embedded in a molecular
cloud phase. However, even in this strong condition, it is not possible to
completely explain the FIR upper limits measured by Capak et al. (2015) and
Barišić et al. (2017). Another possibility is to consider a highly
inhomogeneous dust distribution produced by dust accretion of surface metals
(e.g., Mancini et al. 2016), which would produce higher values of $\beta$ at
the same IRXs. Nevertheless, this effect is not strong enough to reconcile the
theoretical expectations of UV slopes and FIR emission with the observed
behavior of the zombie galaxies, particularly the upper limits in the IRX at
$\beta\gtrsim-1.5$ (see Mancini et al. 2016, for further details). Therefore,
we remark that the existence of a weak AGN activity at the center of such
high-$z$ objects appears to be a more reasonable explanation for their
peculiar IRX-to-$\beta$ values.

##  3 Color-magnitude relation for AGN

The difficulty to reproduce the IRX-to-$\beta$ relation for the sources by
Capak et al. (2015) within the scenario of galaxy-dominated objects can be
avoided by assuming that these sources are actually hosting AGN activity with
UV spectral emission depending on the details of the accretion process.
Figures 1, 2, and 3 show how a UV spectrum with $\beta_{\rm int}$ not fixed to
galactic values may explain the bulk of the data after the reanalysis by
Barišić et al. (2017). If this alternative scenario holds, these objects
should also follow typical relations that are known for AGN, such as showing
steeper spectral slopes at higher luminosities (e.g., Shankar et al. 2016).

Figure 5 shows $\beta$ as a function of $M_{\rm UV}$ for the HerS/SDSS
quasars, the Seyfert galaxies by García-González et al. (2016), and the
high-$z$ galaxies by Capak et al. (2015). In order to quantify the
$\beta$-to-$M_{\rm UV}$ relation for AGN, we perform a Monte Carlo linear
regression by running the script PyMC3 (Salvatier et al. 2016) on the joint
sample of quasar and Seyfert spectral slopes $\beta$ as a function of their UV
absolute magnitude $M_{\rm UV}$. We choose a standard likelihood defined as

$\mathcal{L}\propto\exp{\left(-\frac{\chi^{2}}{2}\right)}=\exp{\left\\{-\frac{\left[\beta-(a\cdot
M_{\rm UV}+b)\right]^{2}}{2\sigma^{2}}\right\\}}$ (9)

with very broad priors for the parameters $a$, $b$, and $\sigma$:

$\begin{array}[]{lr}a\sim\mathcal{N}_{(0,10)},\\\
b\sim\mathcal{N}_{(0,10)},\\\ \sigma\sim\mathcal{N^{+}}_{(0,1)},\\\
\end{array}$ (10)

where $\mathcal{N}_{(\mu,k)}$ is the Gaussian distribution of UV magnitudes
with a given average $\mu$ and variance $k^{2}$, and $N^{+}_{(\mu,k)}$ is its
positive side. In this way, we find

$\beta_{\rm AGN}=\left(0.24\pm 0.02\right)\cdot\left(M_{\rm
UV}+19.5\right)-\left(0.47\pm 0.03\pm 0.14\right),$ (11)

where we adopt the average uncertainty in the measurement of $\beta$ for
HerS/SDSS objects as an estimate of the systematic error on the intercept of
the color-magnitude relation.

Figure 5: Color-magnitude relation for the HerS/SDSS quasars by Dong & Wu
(2016, cyan dots), the Seyfert galaxies by García-González et al. (2016,
magenta crosses), and the high-$z$ galaxies by Barišić et al. (2017, black
squares). The $\beta$-to-$M_{\rm UV}$ relation derived by fitting the
HerS/SDSS quasars (blue solid line) is shown superimposed to the data. For
comparison, the relation valid for galaxies at $z\sim 5$ (Bouwens et al. 2014,
red dashed line) is also shown.

ID | $\chi_{\rm AGN}^{2}$ | $\chi_{\rm gal}^{2}$ | $P_{\rm AGN}$
---|---|---|---
HZ1 | 9.663 | 4.481 | 0.07
HZ2 | 9.175 | 2.649 | 0.04
HZ3 | 6.967 | 0.597 | 0.04
HZ4 | 11.25 | 7.939 | 0.16
HZ5 | 1.001 | 6.270 | 0.93
HZ5a | – | – | –
HZ6 | 0.255 | 2.319 | 0.74
HZ6a | 0.181 | 0.796 | 0.58
HZ6b | 0.190 | 0.008 | 0.48
HZ6c | 0.732 | 0.649 | 0.49
HZ7 | 1.842 | 0.927 | 0.39
HZ8 | 1.994 | 0.577 | 0.33
HZ8W | 3.348 | 28.95 | 1.00
HZ9 | 4.105 | 0.015 | 0.11
HZ10 | 7.697 | 2.208 | 0.06
HZ10W | 2.073 | 0.406 | 0.30

Table 2: Values of $\chi^{2}$ obtained by comparing the Barišić et al. (2017)
data with Eq. 11 and 12, respectively, along with the corresponding Bayesian
probabilities of AGN domination (see text for details).

Conversely, the dependence of $\beta$ on UV emission in star-forming galaxies
is mainly determined by the amount of extinction. Therefore, the expected
trend is to find shallower spectral slopes in bigger, more luminous, more
chemically enriched objects with respect to fainter ones (e.g., Meurer et al.
1999; Bouwens et al. 2014; Finkelstein et al. 2015a, b). Here we adopt the
relation inferred by Bouwens et al. (2014) at $z\sim 5$ to describe the
relation between $\beta$ and $M_{\rm UV}$ in galaxy-dominated objects,

$\beta_{\rm gal}=-(0.14\pm 0.02)\cdot\left(M_{\rm UV}+19.5\right)-(1.91\pm
0.02\pm 0.06),$ (12)

which exhibits an opposite trend with respect to the relation for AGN.

The sample of objects by Capak et al. (2015) is numerically small if compared
to the HerS/SDSS ensemble (207 objects) or the Bouwens et al. (2014) catalog
($>$4000 objects). Furthermore, a visual inspection of Fig. 5 reveals that
this sample tends to cluster around the $\beta$-to-$M_{\rm UV}$ relation for
high-$z$ galaxies. Nevertheless, we aim to test statistically the probability
that these objects are AGN- or galaxy-dominated in terms of relative proximity
to one relation with respect to the other. In order to do this, we define a
$\chi^{2}$ variable as

$\chi^{2}_{\rm x}=\frac{\left(\beta-\beta_{\rm
x}\right)^{2}}{\sigma_{\beta}^{2}+\sigma_{\rm x}^{2}},$ (13)

where $\beta$ and $\sigma_{\beta}$ refer to the observed slopes with their
average errors by Barišić et al. (2017), while $\beta_{\rm x}$ are the values
of UV spectral slopes expected for the Capak et al. (2015) objects using
either Eq. 11 ($\beta_{\rm x}=\beta_{\rm AGN}$) or Eq. 12 ($\beta_{\rm
x}=\beta_{\rm gal}$) with associated standard deviations $\sigma_{\rm x}$. The
Bayesian probability $P_{\rm AGN}$ for a data point by Barišić et al. (2017)
to come from an AGN-dominated object is then

$P_{\rm AGN}=\frac{1}{1+\exp{\left(\Delta\chi^{2}/2\right)}},$ (14)

with $\Delta\chi^{2}=\chi^{2}_{\rm AGN}-\chi^{2}_{\rm gal}$ for an object
having $\chi^{2}_{\rm AGN}$ with respect to Eq. 11 and $\chi^{2}_{\rm gal}$
with respect to Eq. 12. The complementary probability $P_{\rm gal}=1-P_{\rm
AGN}$ is thus the probability that the data point represents a galaxy-
dominated object. We note that we assume a flat prior on both relations (i.e.,
the objects by Capak et al. 2015 are not forced to follow a priori one between
Eq. 11 or 12).

We report the probabilities in Table 2 for each of the high-$z$ objects by
Barišić et al. (2017). Only HZ5 and HZ8W can be classified as AGN hosts at a
90% confidence level, with HZ5 already recognized as an X-ray selected quasar
(Capak et al. 2015), whereas HZ1 to HZ3 and HZ10 are likely to be normal
galaxies. The other targets show intermediate values of $P_{\rm AGN}$, and
hence it is not possible to strictly classify them in one category, though all
of the objects are enclosed within the scatter of the joint AGN color-
magnitude relation. Therefore, no firm conclusion may be drawn about the
relative spectral domination of AGN or galaxy emission in the sample of
objects by Capak et al. (2015). In order to better distinguish AGN-dominated
objects from galaxies, we define as $\Delta\beta$ the distance of the observed
$\beta$ from the $\beta$-to-$M_{\rm UV}$ relation for AGN described by Eq. 11,

$\Delta\beta=\beta^{i}-0.24\cdot\left(M^{i}_{\rm UV}+19.5\right)+0.47,$ (15)

where $\beta^{i}$ and $M^{i}_{\rm UV}$ are the UV spectral slope and absolute
magnitude for the $i$-th object.

We show this quantity in Fig. 6 as a function of the IRX for the Dong & Wu
(2016), the García-González et al. (2016) and the Barišić et al. (2017) data.
It can be seen at first glance that the AGN are clustered by construction
around a value of $\Delta\beta\sim 0$, and scattered over a range of IRX
$\in[-2,3]$ with the bulk of the data comprised between IRX $\sim-1$ and $2$.
The objects from Capak et al. (2015) are superimposed on the region covered by
actual AGN, with the tendency to populate the lower part of the diagram.
Interestingly, segregating these objects by their classification based on
Bayesian probability shows that the galaxy-dominated objects (i.e., those
objects with $P_{\rm gal}>0.90$) cluster around $\Delta\beta\sim-0.8$, while
the AGN hosts exhibit positive values of $\Delta\beta$ up to $\sim$1\. This
diagram synthesizes all the available information about the emission
properties of the Capak et al. (2015) galaxies as reanalyzed by Barišić et al.
(2017) since it correlates $M_{\rm UV}$ with the $\rm IRX$ and $\beta_{\rm
obs}$. We note that this kind of analysis could be used for a preliminary
discrimination between AGN hosts and normal star-forming galaxies where a
measurement of the FIR emission is feasible, although a more thorough analysis
of their UV and FIR properties is needed. Furthermore, the IRX–$\Delta\beta$
plane might also be useful for analyzing the evolution of the local galaxies.

Paper | $m$ | $q$ | Sample
---|---|---|---
Pettini et al. (1998) | $0.39$ | $1.11$ |  $z\sim 3$ Ly-break galaxies
Meurer et al. (1999) | $0.74$ | $2.11$ | Local starbursts
Penner et al. (2012) | $0.87$ | $2.71$ |  $z\sim 2$ DOGs
Casey et al. (2014) | $0.75$ | $0.54$ | Local low-IRX galaxies
Talia et al. (2015) | $0.41$ | $1.81$ |  $1<z<3$ SFG spectra
Forrest et al. (2016) | $0.88$ | $1.75$ |  $1<z<3$ composite SEDs

Table 3: Coefficients of the IRX-to-$\Delta\beta$ relations specialized from
literature results on different samples of galaxies.

To directly compare the effect of removing the AGN dependence on $M_{\rm UV}$
from samples of different objects (quasars, Seyferts, local starbursts,
medium-$z$ star-forming galaxies), we rewrite the IRX-to-$\beta$ relations
available in the literature as a function of $\Delta\beta=\beta_{\rm
gal}-\beta_{\rm AGN}$. Substituting Eq. 12 into Eq. 11 to eliminate the
dependence on $M_{\rm UV}$ and then inverting the result to write $\beta_{\rm
gal}$ as a function of $\Delta\beta$, we find

$\Delta\beta=2.71\beta_{\rm gal}+3.74$ (16)

Therefore, the IRX-to-$\Delta\beta$ relations assume the form

${\rm
IRX}\left(\Delta\beta\right)=\log{\left[10^{0.4\left(m\cdot\Delta\beta+q\right)}-1\right]}+0.076,$
(17)

with the coefficients $m$ and $q$ depending on the analyzed sample of objects.
Unlike in the previous sections of this paper, here we do not restrict
ourselves to the relations by Pettini et al. (1998) and Meurer et al. (1999).
Instead, we focus on several relations available in the literature that are
valid for samples of different classes of galaxies. In Table 3, we list the
values of the coefficients $m$ and $q$ for the IRX-to-$\Delta\beta$ relations,
along with the source paper of the original IRX-to-$\beta$ relation and a
synthetic description of the corresponding type of astrophysical objects.

Figure 6 shows these IRX-to-$\Delta\beta$ relations superimposed on the AGN
samples and on the high-$z$ galaxies by Capak et al. (2015). Apart from the
curve corresponding to the IRX-to-$\beta$ relation derived by Casey et al.
(2014) for star-forming galaxies at $z<0.085$, it is evident that none of the
considered relations is able to intercept the Capak et al. (2015) galaxies
with $\rm IRX<0$. This further strengthens the AGN-host hypothesis for the
high-$z$ galaxies with anomalous IRXs, and again suggests that deepening the
investigation of the IRX–$\Delta\beta$ plane should be a future tool for
object discrimination.

Figure 6:  IRX as a function of the slope difference $\Delta\beta$ for the
HerS/SDSS quasars by Dong & Wu (2016, blue dots), the Seyfert galaxies by
García-González et al. (2016, blue crosses), and the high-$z$ galaxies by
Barišić et al. (2017). The high-$z$ galaxies are classified according to the
statistical analysis presented in Sect. 3 (red triangles: actual galaxies,
yellow squares: unclassifiable objects, green stars: actual quasars). The
curves available in the literature that describe the IRX-to-$\Delta\beta$
relation for various classes of galaxies (Pettini et al. 1998; Meurer et al.
1999; Penner et al. 2012; Casey et al. 2014; Talia et al. 2015; Forrest et al.
2016) are also shown.

ID | $\log{\left(L_{\rm UV}/L_{\odot}\right)}$ | SFR ($M_{\odot}$/yr) | $\log{\left(L_{\rm X}/L_{\odot}\right)}$ | $\log{\left(L^{\star}_{\rm X}/L_{\odot}\right)}$
---|---|---|---|---
HZ1 | 11.21 | 24 | 10.31 | 8.19
HZ2 | 11.15 | 25 | 10.26 | 8.20
HZ3 | 11.08 | 18 | 10.21 | 8.10
HZ4 | 11.28 | 51 | 10.36 | 8.38
HZ5∗ | 11.45 |  $<$3 | 10.82∗∗ |  $<$7.59
HZ5a |  $<$10.37 |  $<$3 |  $<$9.66 |  $<$7.59
HZ6 | 11.47 | 49 | 10.50 | 8.35
HZ6a | 11.11 | — | 10.23 | —
HZ6b | 11.0 | — | 10.15 | —
HZ6c | 11.81 | — | 10.00 | —
HZ7 | 11.05 | 21 | 10.19 | 8.11
HZ8 | 11.04 | 18 | 10.18 | 8.06
HZ8W | 10.57 | 6 | 9.82 | 7.76
HZ9 | 10.95 | 67 | 10.11 | 8.46
HZ10+ | 11.14 | 169 | 10.25 | 8.72
HZ10W+ | 10.23 | — | 9.54 | —
∗No [C II] redshift measurement available (Capak et al. 2015), photometric
redshift is used.
∗∗Luminosity in the [0.5–2] keV energy range observed in the C-COSMOS field.
+These objects fall outside the C-COSMOS field.

Table 4: Luminosity properties of the objects described by Capak et al.
(2015). The expected X-ray luminosities in the observer-frame soft band have
been estimated by adopting a bolometric correction to the UV luminosity of 5.2
(Runnoe et al. 2012b).

##  4 Expected X-ray luminosities

We now present the analysis of the X-ray emission expected for these objects
within the AGN-host scenario. In fact, it has been clearly established that
the accretion process around massive black holes can power the emission of
X-rays from structures (jets, coronae) collateral to the ionized disk (e.g.,
Halpern 1982; Begelman et al. 1983; Begelman & McKee 1983; Zdziarski 1985).
Therefore, it is possible to select and characterize AGN by studying the
details of their high-energy spectra (e.g., Giustini et al. 2008; Lusso et al.
2010; Lusso & Risaliti 2016; Page et al. 2017). In order to give a complete
view of the X-ray energetics of the Capak et al. (2015) galaxies, we both
evaluate the expected X-ray fluxes of such objects from their measured UV-to-
IR properties and perform a stacking analysis on the Chandra X-ray imaging of
the COSMOS field at the relevant sky coordinates.

###  4.1 Estimated upper limits on the X-ray fluxes

The field of the Capak et al. (2015) galaxies has been observed by the Chandra
X-ray observatory in the framework of the Chandra COSMOS Survey222An extension
of the C-COSMOS program, the 4.6 Ms Chandra COSMOS Legacy Survey (e.g.,
Marchesi et al. 2016), is present in the literature, but its imaging data are
not publicly available for stacking analysis. (C-COSMOS; Elvis et al. 2009;
Puccetti et al. 2009; Civano et al. 2012; Lanzuisi et al. 2013). This survey
has imaged for 1.8 Ms the central region of the COSMOS field around 10h +02∘
with an average flux limit of $7.3\cdot 10^{-16}$ erg s-1 cm-2 in the full
band ($0.5-10$ keV). Of the objects listed in Table 1, HZ10 and HZ10W fall
outside the area of 0.9 deg2 covered by the C-COSMOS program, and the weak
quasar HZ5 is detected as a resolved source. None of the other galaxies in the
Capak et al. (2015) sample have been detected.

Assuming that the emission of the non-detected sources in the X-ray comes from
AGN-dominated objects, we follow Runnoe et al. (2012b) in computing the
bolometric luminosities $L_{\rm bol}$ of such sources from the values of
$L_{\rm UV}$ reported in Table 1. Therefore, we adopt $L_{\rm
bol}=f{\alpha_{\rm UV}}L_{\rm UV}$, where the UV-to-bolometric correction is
$\alpha_{\rm UV}=5.2$, and $f=0.75$ is a factor taking into account the
anisotropic emission from an accretion disk. The luminosities in the hard
($2-10$ keV) rest-frame X-ray band $L_{\rm X}$ can then be inferred from the
relation obtained by Marconi et al. (2004),

$\log\left(\frac{L_{\rm bol}}{L_{\rm
X}}\right)=1.54+0.24\mathfrak{L}+0.012\mathfrak{L}^{2}-0.0015\mathfrak{L}^{3},$
(18)

where $\mathfrak{L}=\log(L_{\rm bol})-12$. These X-ray luminosity evaluations
must be interpreted as upper limits on the intrinsic energetic output of these
objects since they would overestimate the real values in the case of non-
negligible contribution of the host-galaxy emission in the UV-to-optical band.

Nevertheless, a significant fraction of the observed X-ray emission,
especially at low X-ray luminosities ($L_{\rm X}\sim 10^{43}$ erg s-1), may be
powered by the cumulative radiation produced by high-mass X-ray binaries
(HMXBs) in the host galaxy (Ranalli et al. 2003, 2005; Lehmer et al. 2010,
2016). Therefore, we estimate this contribution $L^{\star}_{\rm X}$ in the
hard X-ray band following Volonteri et al. (2017), and adopting the redshift-
dependent relation found by Lehmer et al. (2016) for galaxies in the redshift
range $z\in[0,7]$

$\log\left(L^{\star}_{\rm X}\right)=39.82+0.63\log({\rm SFR})+1.31\log(1+z),$
(19)

where the selected values for $z$ and SFR are shown in Table 1 and Table 4,
respectively.

We report the predicted luminosities $L_{\rm X}$ and $L^{\star}_{\rm X}$ for
the analyzed sample in Table 4. It is clearly visible that, on average, the
foreseen luminosities assuming an AGN-driven X-ray emission are $\sim$2 orders
of magnitude higher than the values obtained under the hypothesis of X-ray
flux produced by an HMXB population. Even accounting for a host-galaxy
contribution of $\sim$70%, the two predicted emissions still remain very
different.

We further note that, at high redshifts, X-ray fluxes can be absorbed by
highly obscured environments, such as large amounts of gas and dust
surrounding accreting BHs (Ueda et al. 2003; Gilli et al. 2007; Treister et
al. 2009). The two main attenuation processes in the X-rays are the
photoelectric absorption, accounted for by the cross section $\sigma_{\rm
ph}$, and the Compton scattering of photons against free electrons,
$\sigma_{\rm T}$. The photoelectric cross section decreases with increasing
energy, and the contribution of the dust to $\sigma_{\rm ph}$ is relatively
modest for $E\gtrsim 1$ keV (Morrison & McCammon 1983). Therefore, at higher
energies ($E\gtrsim 10$ keV for gas metallicities $Z\sim Z_{\odot}$), the
energy-independent $\sigma_{\rm T}$ becomes dominant. Softer X-ray photons are
thus expected to be more absorbed than harder ones. Furthermore, X-ray
observations of AGN at $z\gtrsim 4$ typically sample the rest-frame hard X-ray
emission, i.e., in the $3-12$ keV range on average for the sample by Capak et
al. (2015). Even assuming a strong obscuration produced by gas column
densities $N_{\rm H}\sim 10^{23}$ cm2, the unabsorbed AGN SED for $E\gtrsim 3$
keV and $Z\lesssim 0.1Z_{\odot}$ is at most 1.2 times higher than the absorbed
one (see, e.g., fig. 3 in Pezzulli et al. 2017). The theoretical simulations
(Mancini et al. 2016) and observations (Mannucci et al. 2010; Ma et al. 2016)
both suggest that the expected metallicity for this kind of objects is
$\lesssim 0.1Z_{\odot}$. This suggests that neither absorption nor a dominant
contribution of the stellar emission to the observed optical photometry is
able to reduce the predicted AGN-driven X-ray emission enough to make $L_{\rm
X}$ and $L^{\star}_{\rm X}$ comparable. Deeper X-ray observations may greatly
help to discriminate between zombie and active galaxies.

Figure 7: Results of the stacking analysis performed on the C-COSMOS field at
the coordinates of the Capak et al. (2015) high-$z$ galaxies given by Barišić
et al. (2017, see Table 1). Upper panels: $100^{\prime\prime}\times
100^{\prime\prime}$ stacked map of the $0.5-2$ keV Chandra passband (left) and
corresponding count-rate histogram (right). Lower panels: the same for the
$2-8$ keV Chandra passband. In both histograms, the count-rate frequencies
(blue histogram) are shown along their cumulative distributions (red
histogram) and their nominal mean (green dot) with 1$\sigma$ uncertainty.

###  4.2 Stacking analysis of the Chandra X-ray imaging

We perform the stacking analysis of the C-COSMOS data at the positions of the
remaining sources in order to reveal possible statistical signals, or at least
get a flux upper limit that may be useful in order to plan future deep
observations with X-ray facilities. For this purpose, we use the software
CSTACK, a web-based tool for the stacking analysis of Chandra data developed
by T. Miyaji333Available at http://cstack.ucsd.edu/cstack/ or
http://lambic.astrosen.unam.mx/cstack/. (see Appendix A for further details).
We only stack sources with an upper limit in $L_{\rm FIR}$, i.e., HZ1, HZ2,
HZ3, HZ5a, HZ7, HZ8, and HZ8W, because the objects with FIR detection could
host heavily obscured AGN in which the X-ray emission is almost totally
suppressed by the presence of high column densities of gas and IR-emitting
dust (e.g., Ricci et al. 2017a). In addition, such objects are the most
distant from the traditional IRX-to-$\beta$ relations adopted for galaxies,
further suggesting that their spectral properties can only be explained by AGN
domination.

We show the CSTACK output in Fig. 7. While no stacked signal is detected in
the hard band, the soft band exhibits an average count rate of $(1.95\pm
1.86)\cdot 10^{-5}$ cts s-1 in the source region which is not compatible with
zero at the $1\sigma$ level. However, due to the poor statistics of such a
measurement, we are forced to treat it as an upper limit to the X-ray flux
rather than a detection. We recall that at $z\sim 5.4$, the observer-frame
soft band at $0.5-2$ keV actually observes hard $3-13$ keV rest-frame X-rays
(for comparison, the observer-frame hard band at $2-8$ keV would detect
photons produced in the rest-frame range $13-51$ keV). Adopting a mean
redshift of $5.4$, an AGN-typical X-ray photon index $\Gamma\sim 1.8$ (e.g.,
Trakhtenbrot et al. 2017) and an average Galactic $N_{\rm H}\sim 1.9\cdot
10^{20}$ cm-2 (Dickey & Lockman 1990; Kalberla et al. 2005; McClure-Griffiths
et al. 2009; Willingale et al. 2013; HI4PI Collaboration et al. 2016) in the
direction of the COSMOS field, the soft-band count rate upper limit of $\sim
3.8\cdot 10^{-5}$ cts s-1 derived from the stacking analysis corresponds to an
observed flux $\lesssim 5.7\cdot 10^{-16}$ erg s-1 cm-2 in the full Chandra
band. Interpreting the soft-band excess as a detection, it corresponds to a
hard X-ray luminosity of $\sim 4.6\cdot 10^{10}$ $L_{\odot}$ at $z\sim 5.4$,
in excellent agreement with the luminosities reported in Table 4 for the case
of AGN-dominated sources.

##  5 Discussion

In this paper, we have presented and validated the hypothesis of star-forming
galaxies with anomalous IRXs at high redshift (the zombie galaxies) hosting
unveiled AGN activity. This alternative scenario to the dust temperature bias
(Bouwens et al. 2016) and molecule-rich dust clouds (Ferrara et al. 2017)
provides a natural explanation to the IRX-to-$\beta$ values of such objects in
terms of variable intrinsic spectral shape and dust distribution with respect
to inactive galaxies.

The comparison of the IR and UV spectral properties of the high-$z$ galaxies
discovered by Capak et al. (2015) and reanalyzed by Barišić et al. (2017) with
those of quasars (Dong & Wu 2016) and Seyferts (García-González et al. 2016)
shows that high-$z$ star-forming galaxies and AGN with bolometric luminosities
from $10^{43}$ to $10^{48}$ erg s-1 populate the same region of the IRX-
to-$\beta$ diagram. Furthermore, we find that the bulk of the data by Barišić
et al. (2017) is described by an IRX-to-$\beta$ relation suitable for AGN-
dominated objects with intrinsic spectral slopes around $\beta_{\rm
int}\sim-2$ and low galactic contribution to the UV emission (i.e., $f_{\rm
q}\sim 0.7$), and that this AGN UV flux is moderately reddened
($A_{1600}\lesssim 1.2$ mag). Model dependencies and degeneracies are already
taken into account in the SED fitting of the AGN comparison samples, and are
therefore not expected to greatly bias these conclusions (see Sect. 2.5).

We further investigated the possibility of ultimately classifying the
anomalous star-forming galaxies from Capak et al. (2015) as actual AGN hosts
based on the correlation between their spectral properties, and on their
possible X-ray emission. Though superimposing within the scatter to the AGN
color-magnitude correlation derived by the Dong & Wu (2016) data, the analyzed
objects also lie at the boundary of the calibration region of the color-
magnitude relation for inactive galaxies by Bouwens et al. (2014). Therefore,
a clear distinction between galaxy-dominated objects and AGN hosts is hard to
obtain on the basis of their statistical proximity to one relation with
respect to the other. We are in fact only able to classify four objects (HZ1
to HZ3 and HZ10) as normal galaxies and two (HZ5 and HZ8W) as AGN.

Finally, since an analysis of the X-ray emission is crucial in order to
definitely assess the nature of such sources, we attempted to derive their
expected X-ray luminosities in the opposite cases of AGN domination and galaxy
domination, and to detect an X-ray signal from the FIR-undetected objects
through stacking analysis. Interestingly, a signal not compatible with zero at
the 1$\sigma$ level arises from stacking the positions of the Capak et al.
(2015) objects in the Chandra soft band at $0.5-2$ keV, hinting at a possible
rest-frame hard X-ray emission from at least some of the sources. According to
the luminosities derived in Table 4, such an emission is compatible with the
flux expected from AGN-dominated sources at $z\sim 5.4$. Hence, future X-ray
follow-ups of these objects or surveys through Athena are key to disentangling
the two populations of objects.

Our work is based on the hypothesis that all of the zombie galaxies host an
AGN at their center. Is this a overly optimistic scenario that would
ultimately lead to overestimating the actual number of AGN at $z>5$? To test
this scenario, we computed the expected AGN amount in the redshift interval
probed by the analyzed objects. We adopted the quasar luminosity function by
Assef et al. (2011) for the interval $5<z<5.85$ (see, e.g., their figure 7).
At those redshifts, the value of the quasar space density is between $\sim
3.4\cdot 10^{-6}$ and $\sim 1.1\cdot 10^{-8}$ Mpc-3 mag-1 at the 1$\sigma$
level for sources with $-25<M_{J}<-22.7$ (i.e., the $J$-band absolute
magnitude of the zombie galaxies under the hypothesis of AGN domination),
which corresponds to a number of up to 55 high-$z$ quasars in the COSMOS field
of view. Restricting the area to the $C$-COSMOS value of 0.9 deg2 implies a
decrease of this amount to $\sim$25 objects. It should be noted that a NED
search for actual quasars at $z\sim 5.4$ in the COSMOS field resulted in two
sources (SDSS J095728.55+022037.8 and SDSS J100100.04+021300.6; Richards et
al. 2009). Therefore, the sample of (at most) 16 high-$z$ objects discovered
by Capak et al. (2015) fits inside the currently expected statistics of
quasars. This topic will be further exploitable taking advantage of the
upcoming releases of COSMOS spectroscopic redshifts (Salvato et al., in
preparation)444We refer to http://cosmos.astro.caltech.edu/page/
specz for updates..

In addition, performing high-precision measurements with the Very Large
Baseline Interferometer (VLBI) is an option to probe the radio emission of
those objects, aiming at detecting possible jets associated with the accreting
AGN (see, e.g., Herrera Ruiz et al. 2017, and references therein). In
conclusion, a multiwavelength analysis of such high-redshift objects is not
only desirable, but required in order to correctly estimate the star formation
history and the evolution of galaxies. To this end, the James Webb Space
Telescope (JWST; Gardner et al. 2006; Kalirai 2018), the Wide Field Infrared
Survey Telescope (WFIRST; Content et al. 2013), the Atacama Large
Millimeter/submillimeter Array (ALMA), and many other facilities play a key
role in acquiring new data at $z>5$, thus opening new windows to better
understand the vast, complex panorama of (mostly unknown) astronomical sources
in the high-redshift Universe.

###### Acknowledgements.

We are grateful to our anonymous referee for the helpful comments. We thank M.
Ginolfi and R. Schneider (University of Rome La Sapienza) for useful
discussions. FT acknowledges support from the Programma per Giovani
Ricercatori – Anno 2014 “Rita Levi Montalcini”. This research has made use of
the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet
Propulsion Laboratory, California Institute of Technology, under contract with
the National Aeronautics and Space Administration. The development of CSTACK
has been supported by CONACyT Grants No. 83564/179662, UNAM-DGAPA Grant PAPIIT
IN110209/IN104113/IN104216 and Chandra Guest Observer Support Grant
GO1-12178X. The Chandra X-Ray Center (CXC) is operated for NASA by the
Smithsonian Astrophysical Observatory. Reproduced with permission from
Astronomy & Astrophysics, © ESO.

## References

  * Abazajian et al. (2009) Abazajian, K. N., Adelman-McCarthy, J. K., Agüeros, M. A., et al. 2009, ApJS, 182, 543
  * Adelman-McCarthy et al. (2007) Adelman-McCarthy, J. K., Agüeros, M. A., Allam, S. S., et al. 2007, ApJS, 172, 634
  * Ahn et al. (2014) Ahn, C. P., Alexandroff, R., Allende Prieto, C., et al. 2014, ApJS, 211, 17
  * Assef et al. (2011) Assef, R. J., Kochanek, C. S., Ashby, M. L. N., et al. 2011, ApJ, 728, 56
  * Barišić et al. (2017) Barišić, I., Faisst, A. L., Capak, P. L., et al. 2017, ApJ, 845, 41
  * Begelman & McKee (1983) Begelman, M. C. & McKee, C. F. 1983, ApJ, 271, 89
  * Begelman et al. (1983) Begelman, M. C., McKee, C. F., & Shields, G. A. 1983, ApJ, 271, 70
  * Behrens et al. (2018) Behrens, C., Pallottini, A., Ferrara, A., Gallerani, S., & Vallini, L. 2018, MNRAS, 477, 552
  * Bianchi & Schneider (2007) Bianchi, S. & Schneider, R. 2007, MNRAS, 378, 973
  * Bouwens et al. (2016) Bouwens, R. J., Aravena, M., Decarli, R., et al. 2016, ApJ, 833, 72
  * Bouwens et al. (2014) Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2014, ApJ, 793, 115
  * Bouwens et al. (2015) Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2015, ApJ, 803, 34
  * Bowler et al. (2015) Bowler, R. A. A., Dunlop, J. S., McLure, R. J., et al. 2015, MNRAS, 452, 1817
  * Brightman et al. (2014) Brightman, M., Nandra, K., Salvato, M., et al. 2014, MNRAS, 443, 1999
  * Calzetti (1997) Calzetti, D. 1997, AJ, 113, 162
  * Capak et al. (2015) Capak, P. L., Carilli, C., Jones, G., et al. 2015, Nature, 522, 455
  * Casey et al. (2014) Casey, C. M., Scoville, N. Z., Sanders, D. B., et al. 2014, ApJ, 796, 95
  * Cattaneo et al. (2009) Cattaneo, A., Faber, S. M., Binney, J., et al. 2009, Nature, 460, 213
  * Civano et al. (2012) Civano, F., Elvis, M., Brusa, M., et al. 2012, ApJS, 201, 30
  * Content et al. (2013) Content, D., Aaron, K., Abplanalp, L., et al. 2013, in Proc. SPIE, Vol. 8860, UV/Optical/IR Space Telescopes and Instruments: Innovative Technologies and Concepts VI, 88600E
  * Cullen et al. (2017) Cullen, F., McLure, R. J., Khochfar, S., Dunlop, J. S., & Dalla Vecchia, C. 2017, MNRAS, 470, 3006
  * Dickey & Lockman (1990) Dickey, J. M. & Lockman, F. J. 1990, ARA&A, 28, 215
  * Dong & Wu (2016) Dong, X. Y. & Wu, X.-B. 2016, ApJ, 824, 70
  * Draine (2011) Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium (Princeton Univ. Press)
  * Duras et al. (2017) Duras, F., Bongiorno, A., Piconcelli, E., et al. 2017, A&A, 604, A67
  * Elvis et al. (2009) Elvis, M., Civano, F., Vignali, C., et al. 2009, ApJS, 184, 158
  * Fabian (2012) Fabian, A. C. 2012, ARA&A, 50, 455
  * Ferrara et al. (2017) Ferrara, A., Hirashita, H., Ouchi, M., & Fujimoto, S. 2017, MNRAS, 471, 5018
  * Finkelstein et al. (2015a) Finkelstein, S. L., Ryan, Jr., R. E., Papovich, C., et al. 2015a, ApJ, 810, 71
  * Finkelstein et al. (2015b) Finkelstein, S. L., Song, M., Behroozi, P., et al. 2015b, ApJ, 814, 95
  * Forrest et al. (2016) Forrest, B., Tran, K.-V. H., Tomczak, A. R., et al. 2016, ApJ, 818, L26
  * Gallerani et al. (2010) Gallerani, S., Maiolino, R., Juarez, Y., et al. 2010, A&A, 523, A85
  * García-González et al. (2016) García-González, J., Alonso-Herrero, A., Hernán-Caballero, A., et al. 2016, MNRAS, 458, 4512
  * Gardner et al. (2006) Gardner, J. P., Mather, J. C., Clampin, M., et al. 2006, Space Sci. Rev., 123, 485
  * Gil de Paz et al. (2007) Gil de Paz, A., Boissier, S., Madore, B. F., et al. 2007, ApJS, 173, 185
  * Gilli et al. (2007) Gilli, R., Comastri, A., & Hasinger, G. 2007, A&A, 463, 79
  * Giustini et al. (2008) Giustini, M., Cappi, M., & Vignali, C. 2008, A&A, 491, 425
  * Glikman et al. (2006) Glikman, E., Helfand, D. J., & White, R. L. 2006, ApJ, 640, 579
  * Gordon et al. (1997) Gordon, K. D., Calzetti, D., & Witt, A. N. 1997, ApJ, 487, 625
  * Halpern (1982) Halpern, J. P. 1982, PhD thesis, Harvard University, Cambridge, MA.
  * Hao et al. (2013) Hao, H., Elvis, M., Bongiorno, A., et al. 2013, MNRAS, 434, 3104
  * Herrera Ruiz et al. (2017) Herrera Ruiz, N., Middelberg, E., Deller, A., et al. 2017, A&A, 607, A132
  * HI4PI Collaboration et al. (2016) HI4PI Collaboration, Ben Bekhti, N., Flöer, L., et al. 2016, A&A, 594, A116
  * Inoue (2005) Inoue, A. K. 2005, MNRAS, 359, 171
  * Ishibashi et al. (2017) Ishibashi, W., Banerji, M., & Fabian, A. C. 2017, MNRAS, 469, 1496
  * Ishibashi & Fabian (2016) Ishibashi, W. & Fabian, A. C. 2016, MNRAS, 457, 2864
  * Kalberla et al. (2005) Kalberla, P. M. W., Burton, W. B., Hartmann, D., et al. 2005, A&A, 440, 775
  * Kalfountzou et al. (2014) Kalfountzou, E., Civano, F., Elvis, M., Trichas, M., & Green, P. 2014, MNRAS, 445, 1430
  * Kalirai (2018) Kalirai, J. 2018, ArXiv e-prints [arXiv:1805.06941]
  * Kennicutt & Evans (2012) Kennicutt, R. C. & Evans, N. J. 2012, ARA&A, 50, 531
  * Knudsen et al. (2017) Knudsen, K. K., Watson, D., Frayer, D., et al. 2017, MNRAS, 466, 138
  * Krumpe et al. (2014) Krumpe, M., Markowitz, A., & Nikutta, R. 2014, in IAU Symposium, Vol. 304, Multiwavelength AGN Surveys and Studies, ed. A. M. Mickaelian & D. B. Sanders, 265–265
  * Lanzuisi et al. (2013) Lanzuisi, G., Civano, F., Elvis, M., et al. 2013, MNRAS, 431, 978
  * Lawrence et al. (2007) Lawrence, A., Warren, S. J., Almaini, O., et al. 2007, MNRAS, 379, 1599
  * Lehmer et al. (2010) Lehmer, B. D., Alexander, D. M., Bauer, F. E., et al. 2010, ApJ, 724, 559
  * Lehmer et al. (2016) Lehmer, B. D., Basu-Zych, A. R., Mineo, S., et al. 2016, ApJ, 825, 7
  * Lusso et al. (2010) Lusso, E., Comastri, A., Vignali, C., et al. 2010, A&A, 512, A34
  * Lusso & Risaliti (2016) Lusso, E. & Risaliti, G. 2016, ApJ, 819, 154
  * Ma et al. (2016) Ma, X., Hopkins, P. F., Faucher-Giguère, C.-A., et al. 2016, MNRAS, 456, 2140
  * Mainzer et al. (2011) Mainzer, A., Bauer, J., Grav, T., et al. 2011, ApJ, 731, 53
  * Maiolino & Rieke (1995) Maiolino, R. & Rieke, G. H. 1995, ApJ, 454, 95
  * Mancini et al. (2016) Mancini, M., Schneider, R., Graziani, L., et al. 2016, MNRAS, 462, 3130
  * Mancini et al. (2015) Mancini, M., Schneider, R., Graziani, L., et al. 2015, MNRAS, 451, L70
  * Mannucci et al. (2010) Mannucci, F., Cresci, G., Maiolino, R., Marconi, A., & Gnerucci, A. 2010, MNRAS, 408, 2115
  * Marchesi et al. (2016) Marchesi, S., Civano, F., Elvis, M., et al. 2016, ApJ, 817, 34
  * Marconi et al. (2004) Marconi, A., Risaliti, G., Gilli, R., et al. 2004, MNRAS, 351, 169
  * Martin et al. (2005) Martin, D. C., Fanson, J., Schiminovich, D., et al. 2005, ApJ, 619, L1
  * McClure-Griffiths et al. (2009) McClure-Griffiths, N. M., Pisano, D. J., Calabretta, M. R., et al. 2009, ApJS, 181, 398
  * Meurer et al. (1999) Meurer, G. R., Heckman, T. M., & Calzetti, D. 1999, ApJ, 521, 64
  * Morrison & McCammon (1983) Morrison, R. & McCammon, D. 1983, ApJ, 270, 119
  * Mullaney et al. (2011) Mullaney, J. R., Alexander, D. M., Goulding, A. D., & Hickox, R. C. 2011, MNRAS, 414, 1082
  * Narayanan et al. (2018) Narayanan, D., Davé, R., Johnson, B. D., et al. 2018, MNRAS, 474, 1718
  * Ouchi et al. (1999) Ouchi, M., Yamada, T., Kawai, H., & Ohta, K. 1999, ApJ, 517, L19
  * Page et al. (2017) Page, M. J., Carrera, F. J., Ceballos, M., et al. 2017, MNRAS, 464, 4586
  * Pâris et al. (2014) Pâris, I., Petitjean, P., Aubourg, É., et al. 2014, A&A, 563, A54
  * Pei (1992) Pei, Y. C. 1992, ApJ, 395, 130
  * Penner et al. (2012) Penner, K., Dickinson, M., Pope, A., et al. 2012, ApJ, 759, 28
  * Pettini et al. (1998) Pettini, M., Kellogg, M., Steidel, C. C., et al. 1998, ApJ, 508, 539
  * Pezzulli et al. (2017) Pezzulli, E., Valiante, R., Orofino, M. C., et al. 2017, MNRAS, 466, 2131
  * Popping et al. (2017) Popping, G., Puglisi, A., & Norman, C. A. 2017, MNRAS, 472, 2315
  * Puccetti et al. (2009) Puccetti, S., Vignali, C., Cappelluti, N., et al. 2009, ApJS, 185, 586
  * Ranalli et al. (2003) Ranalli, P., Comastri, A., & Setti, G. 2003, A&A, 399, 39
  * Ranalli et al. (2005) Ranalli, P., Comastri, A., & Setti, G. 2005, A&A, 440, 23
  * Reddy et al. (2015) Reddy, N. A., Kriek, M., Shapley, A. E., et al. 2015, ApJ, 806, 259
  * Reddy et al. (2018) Reddy, N. A., Oesch, P. A., Bouwens, R. J., et al. 2018, ApJ, 853, 56
  * Ricci et al. (2017a) Ricci, C., Bauer, F. E., Treister, E., et al. 2017a, MNRAS, 468, 1273
  * Ricci et al. (2017b) Ricci, C., Trakhtenbrot, B., Koss, M. J., et al. 2017b, Nature, 549, 488
  * Richards et al. (2009) Richards, G. T., Deo, R. P., Lacy, M., et al. 2009, AJ, 137, 3884
  * Risaliti et al. (2011) Risaliti, G., Salvati, M., & Marconi, A. 2011, MNRAS, 411, 2223
  * Rogers et al. (2014) Rogers, A. B., McLure, R. J., Dunlop, J. S., et al. 2014, MNRAS, 440, 3714
  * Runnoe et al. (2012a) Runnoe, J. C., Brotherton, M. S., & Shang, Z. 2012a, MNRAS, 427, 1800
  * Runnoe et al. (2012b) Runnoe, J. C., Brotherton, M. S., & Shang, Z. 2012b, MNRAS, 422, 478
  * Salpeter (1964) Salpeter, E. E. 1964, ApJ, 140, 796
  * Salvatier et al. (2016) Salvatier, J., Wiecki, T. V., & Fonnesbeck, C. 2016, PeerJ Computer Science, 2, e55
  * Sandage & Tammann (1987) Sandage, A. & Tammann, G. A. 1987, A Revised Shapley-Ames Catalog of Bright Galaxies
  * Schneider et al. (2010) Schneider, D. P., Richards, G. T., Hall, P. B., et al. 2010, AJ, 139, 2360
  * Schneider et al. (2015) Schneider, R., Bianchi, S., Valiante, R., Risaliti, G., & Salvadori, S. 2015, A&A, 579, A60
  * Scoville et al. (2007) Scoville, N., Aussel, H., Brusa, M., et al. 2007, ApJS, 172, 1
  * Scoville et al. (2014) Scoville, N., Aussel, H., Sheth, K., et al. 2014, ApJ, 783, 84
  * Shakura & Sunyaev (1973) Shakura, N. I. & Sunyaev, R. A. 1973, A&A, 24, 337
  * Shankar et al. (2016) Shankar, F., Calderone, G., Knigge, C., et al. 2016, ApJ, 818, L1
  * Shen & Kelly (2010) Shen, Y. & Kelly, B. C. 2010, ApJ, 713, 41
  * Shen et al. (2011) Shen, Y., Richards, G. T., Strauss, M. A., et al. 2011, ApJS, 194, 45
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163
  * Stalevski et al. (2016) Stalevski, M., Ricci, C., Ueda, Y., et al. 2016, MNRAS, 458, 2288
  * Talia et al. (2015) Talia, M., Cimatti, A., Pozzetti, L., et al. 2015, A&A, 582, A80
  * Trakhtenbrot et al. (2017) Trakhtenbrot, B., Ricci, C., Koss, M. J., et al. 2017, MNRAS, 470, 800
  * Treister et al. (2009) Treister, E., Urry, C. M., & Virani, S. 2009, ApJ, 696, 110
  * Ueda et al. (2003) Ueda, Y., Akiyama, M., Ohta, K., & Miyaji, T. 2003, ApJ, 598, 886
  * Urry & Padovani (1995) Urry, C. M. & Padovani, P. 1995, PASP, 107, 803
  * Vanden Berk et al. (2001) Vanden Berk, D. E., Richards, G. T., Bauer, A., et al. 2001, AJ, 122, 549
  * Vestergaard & Osmer (2009) Vestergaard, M. & Osmer, P. S. 2009, ApJ, 699, 800
  * Vestergaard & Peterson (2006) Vestergaard, M. & Peterson, B. M. 2006, ApJ, 641, 689
  * Viero et al. (2014) Viero, M. P., Asboth, V., Roseboom, I. G., et al. 2014, ApJS, 210, 22
  * Volonteri et al. (2017) Volonteri, M., Reines, A. E., Atek, H., Stark, D. P., & Trebitsch, M. 2017, ApJ, 849, 155
  * Waters et al. (2016) Waters, D., Wilkins, S. M., Di Matteo, T., et al. 2016, MNRAS, 461, L51
  * Weingartner & Draine (2001) Weingartner, J. C. & Draine, B. T. 2001, ApJ, 548, 296
  * Wilkins et al. (2013) Wilkins, S. M., Bunker, A., Coulton, W., et al. 2013, MNRAS, 430, 2885
  * Willingale et al. (2013) Willingale, R., Starling, R. L. C., Beardmore, A. P., Tanvir, N. R., & O’Brien, P. T. 2013, MNRAS, 431, 394
  * Wright et al. (2010) Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868
  * Zdziarski (1985) Zdziarski, A. A. 1985, ApJ, 289, 514
  * Zel’dovich & Novikov (1965) Zel’dovich, Y. B. & Novikov, I. D. 1965, Soviet Physics Doklady, 9, 834

##  Appendix A CSTACK setup

To estimate the amount of the residual X-ray emission from the zombie galaxies
by Capak et al. (2015) described in Sect. 4.2, we adopt the following CSTACK
setup:

  1. 1.

use of the C-COSMOS public data set;

  2. 2.

standard Chandra $0.5-2$ and $2-8$ keV energy bands;

  3. 3.

maximum off-axis angle of $15^{\prime}$;

  4. 4.

size of $100^{\prime\prime}\times 100^{\prime\prime}$ for the stacked image;

  5. 5.

radius of $1^{\prime\prime}.5$ for the source region (corresponding to three
times the Chandra angular resolution);

  6. 6.

inner background radius of $2^{\prime\prime}$;

  7. 7.

union of all energy bands for the exclusion of resolved sources;

  8. 8.

resolved sources excluded from both the source and the background region.

Points 1 and 2 are CSTACK defaults. Point 3 ensures that the whole $C$-COSMOS
field is explored, avoiding for example the exclusion of regions lying at the
field-of-view border. Point 4 only generates images large enough to visually
identify possible signal excesses concentrated at the image center. Points 5
to 8 allow CSTACK to automatically run on source positions and background
regions that do not contain any resolved object in the C-COSMOS field: for
instance, the weak quasar HZ5 is excluded from the stacking, but the nearby
source HZ5a is included due to an angular separation of $1^{\prime\prime}.8$.