Infrared emission of z \sim 6 galaxies: AGN imprints
F. Di Mascia, S. Gallerani, C. Behrens, A. Pallottini, S. Carniani, A. Ferrara, P. Barai, F. Vito, T. Zana
MMNRAS , 1–22 (2020) Preprint February 19, 2021 Compiled using MNRAS L A TEX style file v3.0
Infrared emission of z ∼ galaxies: AGN imprints F. Di Mascia (cid:63) , S. Gallerani , C. Behrens , A. Pallottini , S. Carniani ,A. Ferrara , P. Barai , , F. Vito , T. Zana Scuola Normale Superiore, Piazza dei Cavalieri 7, I-56126 Pisa, Italy Institut f¨ur Astrophysik, Georg-August Universit¨at G¨ottingen, Friedrich-Hundt-Platz 1, 37077, G¨ottingen, Germany N´ucleo de Astrof´ısica - Universidade Cidade de S˜ao Paulo Universidade Cruzeiro do Sul, Rua Galv˜ao Bueno 868, S˜ao Paulo, 01506-000, Brasil
Accepted XXX. Received XXX; in original form XXX
ABSTRACT
We investigate the infrared (IR) emission of high-redshift ( z ∼ >
100 M (cid:12) yr − ) galaxies, with/without Active Galactic Nuclei (AGN), using asuite of cosmological simulations featuring dust radiative transfer. Synthetic Spec-tral Energy Distributions (SEDs) are used to quantify the relative contribution ofstars/AGN to dust heating. In dusty ( M d (cid:38) × M (cid:12) ) galaxies, (cid:38) (cid:38)
100 pc. In runs withAGN, a clumpy, warm ( ≈
250 K) dust component co-exists with a colder ( ≈
60 K)and more diffuse one, heated by stars. Warm dust provides up to 50% of the total IRluminosity, but only (cid:46) .
1% of the total mass content. The AGN boosts the MIR fluxby 10 − × with respect to star forming galaxies, without significantly affecting theFIR. Our simulations successfully reproduce the observed SED of bright ( M UV ∼ − z ∼ z ∼
6. Fi-nally, the MIR-to-FIR flux ratio of faint ( M UV ∼ −
24) AGN is > × higher than fornormal star forming galaxies. This implies that combined JWST/ORIGINS/ALMAobservations will be crucial to identify faint and/or dust-obscured AGN in the distantUniverse. Key words: methods: numerical - dust - galaxies: evolution - galaxies: high-redshift- galaxies: ISM - quasars: supermassive black holes - infrared: general
Gas accretion onto super massive black holes (SMBH, M BH ∼ − M (cid:12) ) residing in the center of most massivegalaxies ( M (cid:63) ∼ − M (cid:12) ; e.g. Kormendy & Richstone1995; Magorrian et al. 1998; Marconi et al. 2004; Kormendy& Ho 2013) turns them into active galactic nuclei (AGN).A large fraction ( ∼ − (cid:63) [email protected] Over the last decade, thanks to several optical/nearinfra-red (NIR) surveys, such as the Sloan Digital Sky Sur-vey (SDSS; Fan et al. 2006; Jiang et al. 2009), the UKIDSSLarge Area Survey (Venemans et al. 2007), the Canada-France High-z Quasar Survey (CFHQS; Willott et al. 2007),the VISTA Kilo-Degree Infrared Galaxy Survey (VIKING;Venemans et al. 2013, 2015), Pan-STARRS1 (Ba˜nados et al.2014), the Very Large Telescope Survey Telescope ATLASsurvey (Carnall et al. 2015), the Dark Energy Survey (DES;Reed et al. 2015, and the Subaru High-z Exploration of Low-luminosity Quasars (SHELLQs; Kashikawa et al. 2015; Mat-suoka et al. 2016), more than 200 quasars have been discov-ered at the most distant redshifts probed so far ( z ∼ − II and C IV ) produced by Broad LineRegion clouds have confirmed that these sources are pow- © a r X i v : . [ a s t r o - ph . GA ] F e b F. Di Mascia et al. ered by ∼ − M (cid:12) BHs (Fan et al. 2000; Willott et al.2003; Kurk et al. 2007; Jiang et al. 2007; Wu et al. 2015).The challenge is to understand how SMBHs have formed in < z ∼
6. Theoret-ical models of black hole accretion are in fact facing seriousdifficulties in explaining such a rapid growth (e.g. Volonteriet al. 2003; Tanaka & Haiman 2009; Haiman 2013; Pacucciet al. 2015; Lupi et al. 2016), also including the rather un-certain formation mechanism of SMBH seeds (Shang et al.2010; Schleicher et al. 2013; Latif et al. 2013; Ferrara et al.2014; Latif & Ferrara 2016).The problem is exacerbated by the unsuccessful searchfor high- z AGN powered by ∼ − M (cid:12) BHs (e.g. Xue et al.2011; Cowie et al. 2020). Whether these sources are too rare(Pezzulli et al. 2017), and/or too faint to be detected bycurrent optical/NIR survey (Willott et al. 2010; Jiang et al.2016; Pacucci et al. 2016; McGreer et al. 2018; Matsuokaet al. 2018; Wang et al. 2019; Kulkarni et al. 2019), and/ortheir optical/UV emission is obscured by dust, remains un-clear. This latter hypothesis is supported by at least twoobservational results: (i) multi-wavelength studies of ∼ α ab-sorption profiles of distant quasars (e.g. Davies et al. 2019).Both these facts resonate with the expectation that earlygrowth of SMBHs, typically characterized by low accretionrates, is buried in a thick cocoon of dust and gas (e.g. Hickox& Alexander 2018, for a review on this subject).In this scenario a certain fraction of UV photons areabsorbed and/or scattered by dust grains in gas clouds inthe host galaxy. By transferring energy and momentum tothe surrounding dusty environment, AGN radiation can sub-stantially affect the conditions of the interstellar (ISM) andcircumgalactic (CGM) medium of the host galaxy in sev-eral ways. UV radiation heats the dust, leading grains tore-emit in the far-infrared. Moreover, radiation pressure ondust grains may drive powerful outflows (e.g. Fabian 1999;Murray et al. 2005; Wada et al. 2016; Venanzi et al. 2020)that push away the gas surrounding the black hole, clean upthe line of sight, and prevent further accretion onto the BH(Di Matteo et al. 2005; Sijacki et al. 2007; Barai et al. 2018).In addition to that, it is unclear whether star formation inthe host galaxy might be quenched (Schawinski et al. 2006;Dubois et al. 2010, 2013; Teyssier et al. 2011; Schaye et al.2015; Weinberger et al. 2018) or triggered (De Young 1989;Silk 2005; Zubovas et al. 2013; Zinn et al. 2013; Cresci et al.2015a,b; Carniani et al. 2016) by AGN-driven outflows.Signatures of such a complex interplay betweenAGN/stellar radiation and dust grains remain imprintedin the rest-frame UV-to-FIR spectral energy distribution(SED) of galaxies. Therefore, multi-wavelength SED anal-ysis of galaxies and AGN can be used to infer informationon their dust properties (mass, temperature, grain size dis-tribution, composition), to shed light on their star formationand nuclear activities, and to quantify the relative contribu-tion of stars and AGN radiation to dust heating (Bongiornoet al. 2012; Pozzi et al. 2012; Berta et al. 2013; Gruppi- oni et al. 2016). Telescopes sensitive to Mid-Infrared (MIR,5 (cid:46) λ RF (cid:46) µ m), like Spitzer (Werner et al. 2004) andHerschel (Pilbratt et al. 2010), and to Far-Infrared (FIR,45 (cid:46) λ RF (cid:46) µ m) wavelengths (e.g. ALMA, NOEMA)have made possible to study the panchromatic SED of bright( M UV (cid:46) −
26) quasars at z ∼ < λ RF < µ m).The result of this study is that star formation may con-tribute 25 −
60% to the bolometric TIR luminosity, withstrong variations from source to source. In particular, Leip-ski et al. (2014) performed a multi-component SED analysison a sample of 69 z > ∼ ∼
50 K)dust component for the FIR emission. This work shows that,in addition to the standard AGN-heated component, a largevariety of dust conditions is required to reproduce the ob-served SED. Yet these kinds of studies are limited to a smallsample of bright sources. Future facilities in the rest-frameMIR, such as the proposed Origins Space Telescope (OST;Wiedner et al. 2020) with a sensitivity ∼ z ∼ M dyn ∼ − M (cid:12) )are characterized by high star formation rates ( SF R ∼ − (cid:12) yr − ), and large amount of molecular gas( ∼ M (cid:12) ) and dust ( ∼ M (cid:12) ), that are typically dis-tributed on galactic scales ( (cid:46) ∼ −
30 kpc), possibly driven by fast outflowinggas ( v out (cid:38) − ) with extreme mass outflow rate of˙ M out ∼ M (cid:12) yr − .These FIR data, combined with X-ray and UV obser-vations, have shown the presence of galaxy/AGN compan-ions in the field of z ∼ ∼ −
50% out to z ∼ . z ∼
6, there are only tentative X-ray detec-tions of double systems (e.g. Vito et al. 2019a; Connor et al.2019, but see also Connor et al. 2020.)The occurrence of UV detected and sub-mm galaxy(SMG) companions is instead more frequent: Marshallet al. (2020) detected up to nine companions with − (cid:46) M UV (cid:46) −
20 in the field of view of six quasars at z ∼ MNRAS000
20 in the field of view of six quasars at z ∼ MNRAS000 , 1–22 (2020)
R emission of early galaxies: AGN imprints ported the [CII] line and 1 mm continuum ( F cont ) detectionof SMGs close to 4 (out of ∼
20) quasars at z ∼
6, with0 . (cid:46) F cont (cid:46) . ∼ ∼
60 kpc.ALMA data of z ∼ T ∼ − K ,Behrens et al. 2018; Arata et al. 2019; Sommovigo et al.2020). The origin of warm dust in early galaxies can betraced back to their (i) large SFR surface densities thatfavour an efficient heating of dust grains (Behrens et al.2018) and (ii) more compact structure of molecular clouds(MC) that delays their dispersal by stellar feedback, imply-ing that a large fraction ( ∼ ∼ M (cid:12) ) BHs, whose UV luminosity is absorbed by dustlocated in the ISM of the host ( (cid:46) ∼ − < M UV < − (i) If high- z galaxies contain an obscured AGN, doesthis imply warmer dust temperatures? (ii) Is there a relationbetween the dust temperature and the BH accretion rate?(iii)What are the most promising spectral ranges and observa-tional strategies to detect obscured AGN? To answer thesequestions it is necessary to build up a model that follows theco-evolution of BHs with their host galaxy from their birthup to the formation of SMBHs powering z ∼ z ∼ z ∼ z ∼ z ∼ − M (cid:12) at z = 6) andthe impact of different AGN feedback prescriptions on theirhost galaxies, residing in a ∼ M (cid:12) dark matter halo.The paper is organised as follows: in Section 2 we il-lustrate both the hydrodynamical simulations (Section 2.1)and the model adopted for the radiative transfer calcula-tions (Section 2.2). We present our results in Section 3 andwe compare them with observations in Section 4. We thenmake predictions for the proposed mission ORIGINS in Sec-tion 5. Finally we summarise our results in Section 6 alongwith our conclusions. We describe the main characteristics of the hydrodyami-cal simulations adopted in this work in Section 2.1 and wepresent the Radiative Transfer (RT) post-processing analy-sis runs performed in Section 2.2, where we also discuss thedetails of the numerical setup and the assumptions made forthe dust properties and emitting sources.
The hydrodynamical cosmological zoom-in simulations usedin this work are described in details in B18 and we sum-marise the main points in the following.B18 use a modified version of the Smooth Particle Hy-drodynamics (SPH) N-body code gadget-3 (Springel 2005)to follow the evolution of a comoving volume of (500 Mpc) ,starting from cosmological initial condition (IC) gener-ated with music (Hahn & Abel 2011) at z = 100 andzooming-in on the most massive dark matter (DM) haloinside the box down to z = 6 . The mass resolution is m DM = 7 . × M (cid:12) and m gas = 1 . × M (cid:12) for DMand gas particles, respectively. For these high-resolutionDM and gas particles the gravitational softening length is1 h − kpc comoving. For the gas, the smoothing length isdetermined at each time step according to the local densityand typically ranges from 300 pc in the ISM ( n ≈
100 cm − )to 6.5 kpc in the CGM ( n ≈ − cm − ).The code accounts for radiative heating and cool-ing according to the tables computed by Wiersma et al.(2009), which also include metal-line cooling. Star forma-tion in the ISM is implemented following the multiphase A flat ΛCDM model is assumed with the following cosmologicalparameters (Planck Collaboration et al. 2016): Ω M , = 0 . Λ , = 0 . B , = 0 . H = 67 .
74 km s − Mpc − . In the low-resolution DM-only simulation, the most massivehalo at z = 6 has a mass of M halo = 4 . × M (cid:12) (virial radius R = 511 kpc comoving), massive enough to host luminousAGN, as suggested by clustering studies (e.g. Allevato et al. 2016).MNRAS , 1–22 (2020) F. Di Mascia et al. model by Springel & Hernquist (2003), adopting a den-sity threshold for star formation of n SF = 0 .
13 cm − anda Chabrier (2003) initial mass function (IMF) in the massrange 0 . −
100 M (cid:12) . Stellar evolution and chemical enrich-ment are computed for the eleven element species (H, He, C,Ca, O, N, Ne, Mg, S, Si, Fe) tracked in the simulation, fol-lowing Tornatore et al. (2007). Kinetic feedback from super-novae (SN) is included by relating the wind mass-loss rate( ˙ M SN ) with the star formation rate ( ˙ M (cid:63) ) as ˙ M SN = η ˙ M (cid:63) and assuming a mass-loading factor η = 2. The wind kineticenergy is set to a fixed fraction χ of the the SN energy: ˙ M SN v = χ(cid:15) SN ˙ M (cid:63) , where v SN = 350 km s − is the windvelocity and (cid:15) SN = 1 . × erg M − (cid:12) is the average energyreleased by a SN for each M (cid:12) of stars formed .In the simulation each BH is treated as a collisionlesssink particle and the following seeding prescription is used.When a DM halo – that is not already hosting a BH –reaches a total mass of M h = 10 M (cid:12) , a M BH = 10 M (cid:12) BHis seeded at its gravitational potential minimum location.BHs are allowed to grow by accretion of the surroundinggas or by mergers with other BHs. Gas accretion onto theBH is modelled via the classical Bondi-Hoyle-Littleton ac-cretion rate ˙ M Bondi (Hoyle & Lyttleton 1939; Bondi & Hoyle1944; Bondi 1952) and it is capped at the Eddington rate˙ M Edd . The final BH accretion rate ˙ M BH reads as follows:˙ M BH = min( ˙ M Bondi , ˙ M Edd ) . (1)To avoid BHs moving from the centre of the halo in whichthey reside because of numerical spurious effects, we im-plement BH repositioning or pinning (see also e.g. Springelet al. 2005; Sijacki et al. 2007; Booth & Schaye 2009; Schayeet al. 2015): at each time-step BHs are shifted towards theposition of minimum gravitational potential within theirsoftening length. During its growth a BH radiates away afraction of the accreted rest-mass energy, with a bolometricluminosity L bol = (cid:15) r ˙ M BH c , (2)where c is the speed of light and (cid:15) r is the radiative efficiency.B18 set (cid:15) r = 0 .
1, a fiducial value for radiatively efficient,geometrically thin, optically thick accretion disks around aSchwarzschild BH (Shakura & Sunyaev 1973). A fraction (cid:15) f = 0 .
05 of this energy is distributed to the surroundinggas in a kinetic form .In this work we consider the following three runs per-formed by B18, starting from the same ICs: • noAGN : control simulation without BHs. • AGNsphere : simulation accounting for BH accretionand AGN feedback. The kinetic feedback is distributed ac-cording to a spherical geometry. • AGNcone : same as the
AGNsphere run, but with kineticfeedback distributed inside a bi-cone with an half-openingangle of 45 ◦ . In the ISM multiphase model adopted here (Springel & Hern-quist 2003), kicked particles mimicking stellar winds are tem-porarily hydrodynamically decoupled. This procedure may affectboth the properties of the resulting outflows and the structure ofthe surrounding ISM (e.g. Dalla Vecchia & Schaye 2008). We refer to B18 for details about the choice of the value for (cid:15) f and the numerical implementation of the kinetic feedback. In Table 1 we report the main physical properties of thezoomed-in halo at z = 6 . ≈
60 kpc) centred on the halo’s centre of mass. This choiceallows to have an overview of all the relevant dynamicalstructures around the central galaxy, i.e. satellites, clumps,filaments, star forming regions.In Fig. 1 we show the hydrogen column density (toprow) and the star formation rate (middle row) for thezoomed-in halo in the three simulations for a line of sightaligned with the angular momentum of the particles insidethe selected region. In the following, this is our reference lineof sight. From the top row, it can be seen that the centralregion, corresponding to the main galaxy, is characterisedby the highest column density in all the runs. It reachesvalues of N H ∼ × cm − in the noAGN run, whereasit is an order of magnitude lower when AGN feedback isincluded. This is because kinetic feedback kicks gas awayfrom the accreting BHs. In turn, the decreased gas densityquenches the overall SFR density. In fact, star formationrate densities Σ SFR as high as Σ
SFR ≈
600 M (cid:12) yr − kpc − are found in the noAGN run, in sharp contrast with thosein the AGNsphere ( ≈
130 M (cid:12) yr − kpc − , characterized bya total BH accretion rate ˙ M BH = 3 . (cid:12) yr − ), and AGN-cone ( ≈
50 M (cid:12) yr − kpc − , ˙ M BH = 89 M (cid:12) yr − ) cases. Thesame trend is observed also for the total SFR, as reportedin Table 1. We post-process the snapshots at z = 6 . skirt (Baes et al. 2003; Baes & Camps2015; Camps & Baes 2015; Camps et al. 2016). skirt solvesthe continuum radiative transfer problem in a dusty mediumwith a Monte-Carlo approach, by sampling the SED of thesources with a finite number of photon packets (in the fol-lowing simply referred to as photons ). Photons are scat-tered and/or absorbed by dust grains in the simulation vol-ume according to their properties. Dust grains, after be-ing heated up, thermally re-emit the absorbed energy atIR wavelengths. One of the main advantages of the skirt code is its flexibility: it allows the user to handle input datafrom different numerical codes (e.g. Adaptive Mesh Refine-ment and Smooth Particle Hydrodynamic codes), to accountfor different dust properties (i.e. grain size distribution andcomposition), to implement different SEDs for the radiat-ing sources (e.g. stars and accreting BHs), to include manyphysical mechanisms (e.g. dust stochastic heating and self-absorption).To relate the energy absorbed by dust with itswavelength-dependent emissivity we adopt the dust mod-els described in Section 2.2.1. We describe the SED adoptedin different RT runs for stars and accreting BHs in Section2.2.3. Version 8, .MNRAS000
50 M (cid:12) yr − kpc − , ˙ M BH = 89 M (cid:12) yr − ) cases. Thesame trend is observed also for the total SFR, as reportedin Table 1. We post-process the snapshots at z = 6 . skirt (Baes et al. 2003; Baes & Camps2015; Camps & Baes 2015; Camps et al. 2016). skirt solvesthe continuum radiative transfer problem in a dusty mediumwith a Monte-Carlo approach, by sampling the SED of thesources with a finite number of photon packets (in the fol-lowing simply referred to as photons ). Photons are scat-tered and/or absorbed by dust grains in the simulation vol-ume according to their properties. Dust grains, after be-ing heated up, thermally re-emit the absorbed energy atIR wavelengths. One of the main advantages of the skirt code is its flexibility: it allows the user to handle input datafrom different numerical codes (e.g. Adaptive Mesh Refine-ment and Smooth Particle Hydrodynamic codes), to accountfor different dust properties (i.e. grain size distribution andcomposition), to implement different SEDs for the radiat-ing sources (e.g. stars and accreting BHs), to include manyphysical mechanisms (e.g. dust stochastic heating and self-absorption).To relate the energy absorbed by dust with itswavelength-dependent emissivity we adopt the dust mod-els described in Section 2.2.1. We describe the SED adoptedin different RT runs for stars and accreting BHs in Section2.2.3. Version 8, .MNRAS000 , 1–22 (2020)
R emission of early galaxies: AGN imprints simulation run AGN feedback M gas [M (cid:12) ] M (cid:63) [M (cid:12) ] SFR [M (cid:12) yr − ] ˙ M BH [M (cid:12) yr − ] M UV [mag] noAGN no 2 . × . ×
600 -
AGNsphere spherical 2 . × . ×
312 3 . AGNcone bi-conical 1 . × . ×
189 89 -27.97
Table 1.
Summary of the hydrodynamic runs of B18 used in this work. For each run, we indicate the feedback model used in thesimulation and the main physical properties of the zoomed-in halo at z = 6 . ∼
50% of the virial radius): gas mass ( M gas ), stellar mass ( M (cid:63) ), star formation rate (SFR, averaged over the last 10 Myr), and the sumof the accretion rate of all the black holes (BHs) in the selected region ( ˙ M BH ). We further associate to ˙ M BH an intrinsic UV magnitude M UV (see Appendix C). Figure 1.
Morphology of the most massive halo at z = 6 . noAGN (left column), AGNsphere (middle column) and
AGNcone (right column). The top, middle and bottom panel show thehydrogen column density, the star formation rate and the dust surface density (assuming a dust-to-metal ratio f d = 0 .
08, see Section2.2.1), respectively. White empty circles show the location of BHs accreting at ˙ M BH > (cid:12) yr − .MNRAS , 1–22 (2020) F. Di Mascia et al.
Dust formation, growth and destruction processes are nottracked in the hydrodynamic simulations considered here.Similarly to other RT works (Behrens et al. 2018; Arataet al. 2019; Liang et al. 2019), we derive the dust mass dis-tribution by assuming a linear scaling with the gas metal-licity (Draine et al. 2007), parametrizing the mass fractionof metals locked into dust as: f d = M d /M Z , (3)where M d is the dust mass and M Z is the total mass of allthe metals in each gas particle in the hydrodynamical sim-ulation (see Section 2.1). The choice of f d directly affectsthe total dust content. The RT calculation is sensitive tothe f d value, which is poorly constrained by high-redshiftgalaxies observations (see Wiseman et al. 2017 and refer-ences therein) and theoretical models (Nozawa et al. 2015).In particular, recent theoretical works (Asano et al. 2013a;Aoyama et al. 2017) suggest that f d is constant in the earlystages of galaxy evolution and then it grows with metallicityup to the Milky-Way (MW) value of f d = 0 . f d , and focus our attention on how the dust content ofgalaxies affects their panchromatic SED.We adopt two different f d values for the normalization:i) a MW like value ( f d = 0 . f d = 0 . z ∼ f d = 0 .
08 case is shown in the bottom rowof Figure 1. High dust surface density regions correspondto active star forming regions where gas metal enrichmentis more pronounced. Therefore, gas and dust density, andSFR are generally correlated in our simulations, as can beseen in Fig. 1.The properties of dust as chemical composition andgrain size distribution are not known in early (AGN-host)galaxies. The nature and origin of dust at high redshift is infact a widely debated topic (e.g. Valiante et al. 2009; Strattaet al. 2011; Asano et al. 2013b; Hirashita et al. 2015; Hi-rashita & Aoyama 2019). Some works (Maiolino et al. 2004;Gallerani et al. 2010) have suggested that z (cid:38) of Weingartner & Draine(2001). We defer the inclusion of a SN-type extinction curveto a future work. Throughout this paper the gas metallicity is expressed in solarunits, using Z (cid:12) = 0 .
013 as a reference value (Asplund et al. 2009). We consider the revised optical properties evaluated in Draine(2003a,b,c). skirt
Dust is distributed in the computational domain in an octreegrid with a maximum of 8 levels of refinement for high dustdensity regions, achieving a spatial resolution of ≈
230 pcin the most refined cells, comparable with the softeninglength in the hydrodynamic simulation ( ≈
200 physical pcat z = 6 . . We verify in App. A that the number of re-finement levels adopted in our fiducial setup is sufficient toachieve converge of the results. Adopting an SMC-like dust,the grain size distribution of graphite and silicates is sam-pled with 5 bins for each component. Gas particles hotterthan 10 K, are considered dust-free as at these tempera-tures thermal sputtering is very effective at destroying dust(Draine & Salpeter 1979; Tielens et al. 1994; Hirashita et al.2015). This assumption does not affect the main results ofour work, as discussed in App. B.Grain temperature and emissivity are evaluated by im-posing energy balance between the local radiation field anddust re-emission. By default, when dust emission photonspropagate, skirt accounts for the self-absorption by dust,but it does not take this absorption into account when com-puting the dust temperature, unless the self-absorption flagis turned on. As this effect may be relevant if dust is IR-optically thick, we have enabled a self-consistent evaluationof the dust temperature, iterating the RT calculation fordust absorption and re-emission until the dust IR luminos-ity converges within 3%.We also include non-local thermal equilibrium (NLTE)corrections to dust emission, which include the contributionfrom small grains that are transiently heated by individualphotons. In this case grains of different sizes are no longer ata single equilibrium temperature, but follow a temperaturedistribution. Behrens et al. (2018) found in their calculationsthat stochastic heating affects mostly the MIR portion of theSED (rest-frame wavelength (cid:46) µ m) but it has a minorimpact on the FIR and (sub)mm emission.We do not include heating from CMB radiation. As dis-cussed in Section 3.2, only a small fraction of dust grains isat a temperature comparable to T CMB . We expect this effectto be negligible, as seen a posteriori from the RT results.We do not include any subgrid model for dust clumpi-ness. Recent works (e.g. Camps et al. 2016; Trayford et al.2017; Liang et al. 2021) that account for subresolution struc-tures of birth clouds harboring young stars (Jonsson et al.2010), whose typical scales are not resolved by the hydrody-namical simulations, are based on SED templates (Groveset al. 2008) not consistent with our fiducial set up. The stel-lar emission in the Groves et al. (2008) template is, in fact,calculated from Starbust99 models (Leitherer et al. 1999), by When distributing the dust content derived from the hydro-dynamical simulation into an octree grid, a kernel-based inter-polation is required in order to convert the dust content froma particles-based distribution into an octree geometry. This pro-cedure leads to a discrepancy between the total amount of dustcarried by the SPH gas particles imported from the hydrodynami-cal simulation and the effective dust content in the computationaldomain used for the RT calculation. Therefore, it is important tocheck that the structure of the dust grid adopted achieves suffi-cient convergence relative to the overall dust content. We find thatthe relative difference in the overall dust content is within 0 . . .
4% for noAGN , AGNsphere and
AGNcone , respectively.MNRAS000
AGNcone , respectively.MNRAS000 , 1–22 (2020)
R emission of early galaxies: AGN imprints assuming a Kroupa (2002) initial mass function, whereas wemodel stellar emission using the Bruzual & Charlot (2003)model (see Section 2.2.3), based on the Chabrier (2003) IMF.Moreover, they include PAH molecules in the dust composi-tion that are instead not considered in our work. We noticethat Liang et al. (2021) find that the Groves et al. (2008)template mainly affects the IR emission from PAHs whichis enhanced up to 50% (see Fig. 23 of their paper). Giventhat in our model we adopt an SMC dust composition (i.e.no PAHs), we do not expect that the inclusion of a subgridmodel that accounts for dust clumpiness would significantlyaffect the main results of our work. The ultraviolet (UV) radiation field mainly responsible fordust heating is provided by stellar sources and black holes.We describe in the following how the two components areimplemented in our model.Stellar particles in the simulation represent a SingleStellar Population (SSP), i.e. a cluster of stars formed at thesame time and with a single metallicity. Given the mass, ageand metallicity of the stellar particle imported, skirt buildsthe individual SEDs according to the Bruzual & Charlot(2003) family of stellar synthesis models, placing the sourcesat the locations of the stellar particles.Black holes are treated as point source emitters as thetypical sizes of the accretion disk and the dusty torus aremuch smaller ( (cid:46)
10 pc) than the width of the most refinedgrid cells ( ≈
230 pc, see Section 2.2.1). We implement theiremission in skirt adopting a SED as described in Section2.2.4.The radiation field is sampled using a grid covering the rest-frame wavelength range [0 . − ] µ m. The choice ofthe lower limit is quite common for RT simulations in dustygalaxies (Schneider et al. 2015; Behrens et al. 2018) and itis motivated by the fact that codes like skirt typically donot account for the hydrogen absorption of ionising photons( λ <
912 A ◦ ). The choice of 10 µ m as the upper limit of thewavelength grid is motivated by the fact that the intrinsicemission from stars and BHs is negligible above this limit.The base wavelength grid is composed of 200 logarithmicallyspaced bins.A total of 10 photon packets per wavelength bin islaunched from each source, i.e. stellar particles and BHs .We collect the radiation escaping our computational domainfor the 6 lines-of-sight perpendicular to the faces of the cubiccomputational domain. The SED of an AGN is shaped by the numerous physicalmechanisms involved in the process of gas accretion ontothe BH (see Netzer 2015 for a comprehensive review onthis topic). AGN SED templates are typically based bothon theoretical arguments and observations (e.g. Shakura &Sunyaev 1973; Vanden Berk et al. 2001; Sazonov et al. 2004;Manti et al. 2016; Shen et al. 2020), possibly including thedusty torus modelling (Schartmann et al. 2005; Nenkovaet al. 2008; Stalevski et al. 2012, 2016). For this work, weadopt a composite power-law for the AGN emission writtenas: L λ = c i (cid:18) λµ m (cid:19) α i (cid:18) L bol L (cid:12) (cid:19) L (cid:12) µ m − , (4)where i labels the bands in which we decompose the spec-tra and the coefficients c i are determined by imposing thecontinuity of the function based on the slopes α i . The coeffi-cients c i and α i adopted and the relative bands are reportedin Table 2 and they are chosen as described in the following.For the X-ray band, based on the results by Pi-concelli et al. (2005) and Fiore et al. (1994) inthe hard (2 −
10 keV, α X , hard = − . ± .
1) and soft(0 . − − . < α X , soft < .
3) band, respectively, weconsider α X , hard = − . α X , soft = − .
7. Consistentlywith Shen et al. (2020), in the wavelength range 50 −
600 A ◦ we use α = 0 . ◦ for continuity). For the ExtremeUV band (EUV, 600 < λ <
912 A ◦ ) we use α EUV = − . − . < α EUV < − .
3) based on the analysis of near-zonesobserved around high redshift quasars.The analysis of a large sample (4576) of z (cid:46) . (cid:46) λ (cid:46) ◦ has shown that the spectral slopes aredistributed in the range ( − . < α < − .
2) and peak around α = − .
6. In the 912 < λ < ◦ band, Lusso et al.(2015) have constructed a stacked spectrum of 53 quasarsat z ∼ . α = − . ± .
01. Moreover, Galleraniet al. (2010) have analysed 33 quasars in the redshift range3 . (cid:46) z (cid:46) . α = − . ± .
5, whereas reddened quasars prefersteeper slopes ( α < − . F ν ∝ ν / , which trans-lates into α = − .
3. Given the uncertain value of the slopefor wavelengths longer than 912 A ◦ , we consider two possi-ble values for the slope in the range from the UV to NIR: α UV = − .
5, which is representative of unreddened quasars, The total AGN bolometric luminosity is distributed from theX-ray to the IR according to the SED adopted (see Section 2.2.4).The choice of the wavelength range adopted in our simulationsaffects the fraction of the AGN bolometric luminosity that ef-fectively enters in the calculation (see Fig. 2). For the fiducial
SED introduced in Sec. 2.2.4, this fraction is ≈ ≈
40% for the
UV-steep
SED. We verified that the number of packets used is sufficient toachieve numerical convergence by comparing the results with con-trol simulations with 5 × photon packets per wavelength bin.MNRAS , 1–22 (2020) F. Di Mascia et al.
Figure 2.
AGN SED for a bolometric luminosity L bol = 10 L (cid:12) : the fiducial SED ( α UV = − .
5) is shown with a blue thick line; the
UV-steep
SED ( α UV = − .
3) is shown with a red thick line. We plot the SED template derived in Shen et al. (2020) for comparisonwith a thick green line, re-scaling it in order to have the same L ◦ of the fiducial SED. The SEDs differ mainly at wavelength longerthan the UV band, with the UV-steep SED dropping faster than the other two. The fiducial SED is in very good agreement with theShen et al. (2020) SED up to ≈ µ m, from where the contribution by dust in the torus and in the galaxy included in their IR templatebegins to dominate the emission in their SED. As a reference for the SED plots in the following, we also plot our two SEDs as the F ν (in µJy ) vs λ with dotted lines, keeping the same colour legend. and α UV = − .
3. We will refer to these two models as the fiducial and
UV-steep model, respectively.At longer wavelengths, the intrinsic AGN emission isexpected to follow the Rayleigh-Jeans tail regime F ν ∝ ν ,which corresponds to α IR = −
4. The transition betweenthe UV slope and the IR one increases with the blackhole mass (Shakura & Sunyaev 1973; Pringle 1981; Sazonovet al. 2004). In this work, we adopt a transition wavelength λ trans = 5 µ m. This component represents the IR emissionfrom the accretion disk only. We did not include the emis-sion from the hot dust component from the torus becausewe cannot resolve the scales (1 −
10 pc) of the torus itself.We discuss how this affects our results in Section 5.1.The fiducial and UV-steep SEDs adopted in this workare shown in Fig. 2 with blue and red lines, respectively. Wealso report with a green line the SED derived by Shen et al.(2020). Our bolometric corrections (reported in Table 3)are consistent with the ones by Shen et al. (2020) (reportedin the top panel of their Fig. 2), for L bol ≈ erg s − . Wefurther calculate the α OX = 0 .
384 log L ν (2keV) /L ν (2500A ◦ )index for our SEDs and find that it is in agreement with ob- Consistently with Shen et al. (2020), we express the UV bandluminosity as ν ◦ L ν ◦ , the B band as ν ◦ L ν ◦ ,whereas the soft [hard] X-ray luminosity is the integrated lumi-nosity in the 0.5-2 [2-10] keV band. servations of z ∼ We perform RT calculations on the three hydrodynamic sim-ulations presented in section 2.1. For each hydro-simulationwe vary the dust to metal ratio from f d = 0 .
08 to f d = 0 . AGNcone run we consider both the AGN SEDs de-scribed in Section 2.2.4. We end up with a total of 8 post-processed runs, as reported in Table 4.In this section we present the results obtained throughour RT calculations. We first present in Section 3.1 the mor-phology of the ultraviolet (UV, 1000 − ◦ ) and totalinfrared (TIR, 8 − µ m) emission and discuss how it isaffected by the presence of the AGN and total dust con-tent. Then, in Section 3.2 we derive the dust temperaturein the different runs. Finally, we discuss in Section 3.3 thesynthetic SEDs resulting from our calculations. In Fig. 3, we show the UV (top row) and TIR (middle row)emission maps derived for the runs noAGN (left column),
AGNsphere (middle column),
AGNcone (right column) for f d = 0 .
08. In Fig. 4 we show the same maps but for f d = 0 . MNRAS000
08. In Fig. 4 we show the same maps but for f d = 0 . MNRAS000 , 1–22 (2020)
R emission of early galaxies: AGN imprints hard X soft X X to EUV EUV UV to NIR NIR to FIR[2 −
10] keV [6 . −
60] A ◦ [50 − ◦ [600 − ◦ [0 . − µ m [5 − ] µ m c (fiducial) 2 0.042 14.133 1.972 0.111 6.225 α (fiducial) -1.1 -0.7 0.4 -0.3 -1.5 -4.0 c (UV-steep) 0.003 0.066 22.499 3.140 0.026 0.402 α (UV-steep) -1.1 -0.7 0.4 -0.3 -2.3 -4.0 Table 2.
Coefficients of our AGN SEDs models as expressed in eq. 4. The slopes α i and the ranges of the piece-wise decomposition werechosen as explained in Section 2.2.4. Imposing the continuity of the function determines the coefficients c i . The SED built in this way isby construction normalised to the bolometric luminosity of the source expressed in L (cid:12) according to eq. 4.SED model L bol L X , hard L bol L X , soft L bol L UV L bol L B α OX fiducial
130 130 3.4 6.0 -1.65
UV-steep
80 81 3.1 13.6 -1.51
Table 3.
Bolometric corrections ( L bol /L band ) and α OX for the fiducial ( α UV = 1 .
5) and
UV-steep ( α UV = 2 .
3) AGN SED mod-els adopted in this work. The bands used to compute the bolo-metric corrections are defined as: hard X-ray [2 −
10] keV, softX-ray [0 . −
2] keV, UV [0 . − . µ m. L B is defined as λL λ at λ = 4400 A ◦ . The two models mostly differ for the luminosity inthe B band. For a bolometric luminosity L bol = 10 erg s − , ourbolometric corrections are consistent with the observational con-straints reported in the top panel of Fig. 2 by Shen et al. (2020). By comparing the TIR maps with the dust surface den-sity (Fig. 1, bottom row) we see that the morphology of theTIR emission matches the dust distribution, as expected.Moreover, the brightest TIR spots in the noAGN (AGNruns) correspond to the locations of the most highly starforming regions (accreting BHs), responsible for the dustgrains heating. We discuss in more details the dust temper-ature in Section 3.2.For what concerns UV emission, in the noAGN case,its distribution correlates with the star formation surfacedensity (see middle row in Fig.1); in the AGN runs, thebrightest spots are located in correspondence of the AGNpositions, identified by white circles. Noticeably, whereas inthe
AGNcone simulation with f d = 0 .
08 three peaks appearin the UV emission map (labelled as A, B, and C in Fig. 3),corresponding to the AGN positions , in the case f d = 0 . noAGN run 84 −
94% of the total UV emissionis extincted by dust if f d = 0 . − .
3. For comparison, inthe
AGNcone and
AGNsphere runs, the same fraction is77 −
99% and 54 − M d (cid:38) × M (cid:12) ) a largefraction ( (cid:38) The accretion rate quoted in Table 2.1 for the case
AGNcone is in fact the sum of the accretion rates of the most active blackholes in the simulations ( ˙ M BH ≈
32, 7 and 50 M (cid:12) yr − , for thesources A, B, and C, respectively). The range reported for the UV reprocessed luminosityin Tab. 5 refers to the variation occurring along differentlines of sights: the minimum and maximum values differ bya factor that can be as high as ∼ (cid:46)
100 pc) and complex internal structure ( ∼ −
10 pc,Padoan & Nordlund 2011; Padoan et al. 2014; Vallini et al.2017) are not resolved by of our simulations.According to the Unified Model (Urry & Padovani1995), the classification between TypeI (unobscured) andTypeII (obscured) AGN is based on the presence ofa dusty, donut-like shaped structure that is responsiblefor anisotropic obscuration in the circum-nuclear region( <
10 pc). Our results show that large UV luminosity vari-ations with viewing angle, in addition to the ones due to thetorus, arise from the inhomogeneous distribution of dustygas surrounding the accreting BH, on ISM scales ( (cid:38)
200 pc;see also Gilli et al. 2014).
One of the key physical quantities derived from RT cal-culations is the mass-weighted dust temperature ( (cid:104) T d (cid:105) M ).In what follows, we first describe how we compute theluminosity-weighted dust temperature ( (cid:104) T d (cid:105) L , see Behrenset al. 2018; Sommovigo et al. 2020) and compare this valuewith (cid:104) T d (cid:105) M ; then we discuss how the dust temperature isaffected by the total amount of dust, and different types ofUV sources (stars vs AGN). T d To compute (cid:104) T d (cid:105) L , we assume that each dust cell emits as agrey body L TIR ∝ M d T β d d , where β d is the dust emis-sivity index . (cid:104) T d (cid:105) L , finally depends on the total amount ofdust M d in the simulation, determined by our choice of f d .Runs with f d = 0 . ∼
10% lower with respect to the corresponding This approximation holds only for dust cells that are opticallythin to IR radiation, although we caveat that a small number ofcells in the simulation is actually optically thick. The actual value of β d depends on the RT calculation. Forexample, Behrens et al. (2018) found a value of β d = 1 .
7. Forcomputing the luminosity-weighted temperature, we assume β d =2. This choice does not significantly affect the final results: theestimate of (cid:104) T d (cid:105) L varying 1 . < β d < . , 1–22 (2020) F. Di Mascia et al.
RT run name Hydro run name Radiation field AGN SED f d noAGN
008 noAGN stars 0 . noAGN
03 noAGN stars 0 . AGNsphere
008 AGNsphere stars + BHs fiducial 0 . AGNsphere
03 AGNsphere stars + BHs fiducial 0 . AGNcone
008 AGNcone stars + BHs fiducial 0 . AGNcone
03 AGNcone stars + BHs fiducial 0 . AGNcone
UVsteep
AGNcone stars + BHs UV-steep 0 . AGNcone UVsteep
AGNcone stars + BHs UV-steep 0 . Table 4. skirt post-processing runs performed. The first column labels the RT simulation, the second column indicates the correspondinghydrodynamical run, the third column specifies the radiation field included (e.g. stars with or without black holes), the fourth columnspecifies the AGN SED used (if black holes are present), and the fifth column contains the dust to metal ratio f d adopted. Figure 3.
UV (top row), TIR (middle row), luminosity-weighted dust grain temperature (bottom row) maps for the runs with f d = 0 . AGNcone runs, which will be discussed in more details in Section 4.1. MNRAS000
UV (top row), TIR (middle row), luminosity-weighted dust grain temperature (bottom row) maps for the runs with f d = 0 . AGNcone runs, which will be discussed in more details in Section 4.1. MNRAS000 , 1–22 (2020)
R emission of early galaxies: AGN imprints Figure 4.
Same as in figure 3 but for for f d = 0 . runs with f d = 0 .
08. This is because the same UV energy isdistributed over a larger amount of dust mass.In Fig. 5 we show the (cid:104) T d (cid:105) L PDF (blue histograms),compared with the mass-weighted (cid:104) T d (cid:105) M one (red his-tograms) for the noAGN , AGNsphere and
AGNcone sim-ulations with f d = 0 .
08, as a reference case. In each run, thePDF of (cid:104) T d (cid:105) M peaks at lower dust temperatures with re-spect to (cid:104) T d (cid:105) L . The difference between the mass-weightedand luminosity-weighted temperatures is particularly evi-dent in the runs in which AGN radiation is included. Inparticular, the spikes of the luminosity-weighted histogramscorrespond to dust cells in the immediate proximity of ac-creting BHs. This dust component constitutes only a smallfraction of the total mass, but it provides a significant con-tribution to the overall luminosity, as further discussed inthe next section. The brightest TIR spots in Fig. 3 and 4 in the noAGN (AGNruns) correspond to the locations of the most highly starforming regions (accreting BHs). Whereas in the noAGN run the maximum T d value is about 60 K, in the AGNruns, dust grains reach luminosity-weighted temperatures T d (cid:38)
200 K close to BHs, and T d ≈
60 K in the diffuse gas.We underline that in the noAGN run (cid:104) T d (cid:105) L is up to4 times lower with respect to the AGN runs despite hav-ing a star formation rate 3 times higher. These results in-dicate the dominant role played by AGN radiation in thedust heating. This is particularly evident if we compare inmore details the run noAGN and AGNsphere . In the noAGN case, L intrUV = L UV , stars = 4 . × L (cid:12) ; in the AGNsphere case, L intrUV = L UV , stars + L UV , BH = (2 . . × L (cid:12) =3 . × L (cid:12) .Thus, although the UV budget in the AGNsphere run
MNRAS , 1–22 (2020) F. Di Mascia et al.
RT run L UV L TIR M d (cid:104) T d (cid:105) L T min / maxd L intr UV τ UV [10 L (cid:12) ] [10 L (cid:12) ] [10 M (cid:12) ] [K] [K] [10 L (cid:12) ] noAGN
008 5 . − . . − . . ± −
64 4 . . − . noAGN
03 3 . − . . − . ± −
57 4 . . − . AGNsphere
008 7 . −
17 3 . − . . ±
27 17 −
179 3 . . − . AGNsphere
03 2 . − . . − . ±
25 15 −
178 3 . . − . AGNcone
008 27 −
90 43 −
50 3 . ±
78 22 −
282 39 1 . − . AGNcone
03 5 . −
32 54 −
71 13 182 ±
69 20 −
272 39 2 . − . Table 5.
Overview of the main physical properties of the galaxies for the RT runs performed (see Table 1). The table contains: (firstcolumn) the name of the run, (second column) the processed UV (integrated in the band 1000 − ◦ ) luminosity L UV , (third column)the processed total infrared (integrated in the band 8 − µ m) luminosity L TIR , (fourth column) the total dust mass contained in thesimulated region M d , (fifth column) the luminosity-weighted temperature of the dust grains (cid:104) T d (cid:105) L , reported as the mean of the PDF withinone standard deviation, (sixth column) the minimum and maximum value the dust grains temperature, (seventh column) the intrinsic(i.e, not dust-processed) UV luminosity L intrUV , (eighth column) the effective UV optical depth τ UV , estimated as as e − τ UV = L UV /L intrUV .For the dust-processed UV, TIR luminosities and UV optical depth we report the range bracketed by the six line of sights considered foreach simulation. is mostly provided by stars, and the total UV intrinsic lumi-nosity is comparable to the noAGN case, T d peaks at highertemperature values if BH accretion is present. In Fig. 6, wecompare the fraction of mass (left panel) and TIR luminos-ity (right panel) from dust with a temperature above acertain threshold for the three runs. In the noAGN run theTIR luminosity is arising from dust with T d (cid:46)
50 K. In the
AGNsphere ( AGNcone ) run >
50% of the TIR luminosity isarising from dust with T d (cid:38)
70 K ( T d (cid:38)
150 K); this warmdust only constitutes 0.1% of the total dust mass. This con-firms that a small mass fraction of warm dust dominates theIR emission, as expected from the scaling L d ∝ M d T β d d . Fig. 7 shows the fraction of dust mass and infrared lumi-nosity as a function of the distance from the regions withthe highest star formation for the noAGN case and from theBHs with the highest accretion rate for the run AGNsphere and
AGNcone .In the noAGN case, the dust mass within r (cid:46)
300 pcrepresents ∼ .
3% of the total dust, and it provides ∼
3% ofthe total IR luminosity. In the
AGNsphere ( AGNcone ) case,only ∼ .
1% ( ∼ .
06 %) of the total dust mass is foundat r (cid:46)
300 pc from an accreting BH but it contributes 20%( ∼ The luminosity is computed assuming β d = 2 for consistencywith the temperature PDF. The resulting luminosity varying1 . < β d < . ≈
10% from the quoted values. Given that there are multiple accreting BHs, we selected the 2(3) most active ones in the
AGNsphere AGNcone run and 2 mostaccreting star forming regions (the main galaxy and its largestsatellite) in the noAGN run. For each cell containing dust in theoctree grid we evaluate the distance from each reference sourceand then we consider the minimum one for this calculation.
Fig. 8 shows the intrinsic flux density from stars (dashedline) and AGN (dotted line) for the first six runs reportedin Table 4. The higher value of the flux density from starsin the noAGN run with respect to both AGN runs is dueto the negative AGN feedback that in the AGN simulationsquenches the star formation rate in the host galaxy (seeSection 3.7 of B18 for an extensive discussion on this topic).This effect is more pronounced in the
AGNcone run sinceit is characterised by a black hole accretion rate that is afactor of ∼
30 higher than in
AGNsphere (see Table 1).The total intrinsic flux (dotted-dashed line) is comparablebetween noAGN and
AGNsphere (see also Table 5).We now analyse the differences between the reprocessedflux density (observed, solid line) resulting from our calcu-lations, focusing on the rest-frame NIR (1 (cid:46) λ RF (cid:46) µ m),MIR (5 (cid:46) λ RF (cid:46) µ m), and FIR (40 (cid:46) λ RF (cid:46) µ m)wavelength ranges.The intrinsic NIR flux is suppressed by ≈
10 times in allruns; the highest rest-frame UV attenuation is seen in the
AGNcone run, with some (all) lines of sight showing a fluxreduced by ≈
100 times for f d = 0 .
08 ( f d = 0 . AGNcone run still providesthe highest rest-frame UV flux. In this wavelength range,the SED is nearly constant in the noAGN run whereas itincreases toward larger wavelengths in the runs with AGN,as a consequence of the contribution from accretion. Theobserved optical-NIR flux depends both on the radiationfield and dust content.For what concerns the MIR, at short wavelengths,( λ RF ∼ − µ m), the SED is dominated by the almostunattenuated emission from stars and/or AGN; the AGN-cone
SED is ∼
30 times brighter than the
AGNsphere oneas a consequence of its higher BH activity. At longer wave-lengths ( λ RF > µ m), the observed flux arises from heateddust IR emission. The flux density in this wavelength rangeis the result of the sum of multiple greybodies, each emit-ting at different temperatures, according to the luminosity- MNRAS000
AGNsphere oneas a consequence of its higher BH activity. At longer wave-lengths ( λ RF > µ m), the observed flux arises from heateddust IR emission. The flux density in this wavelength rangeis the result of the sum of multiple greybodies, each emit-ting at different temperatures, according to the luminosity- MNRAS000 , 1–22 (2020)
R emission of early galaxies: AGN imprints Figure 5.
Mass-weighted (red histograms) and luminosity-weighted (blue histograms) dust grains temperature PDF. Thepanels refer to: (top) noAGN , (middle)
AGNsphere and (bottom)
AGNcone . As a reference case, we show the results for f d = 0 . weighted dust temperature PDF discussed in Section 3.2.1.The warm dust in AGN runs produces a MIR excess withrespect to the noAGN run, and shifts the peak of theemission toward shorter wavelengths: λ peak noAGN = 59 . µ m, λ peak AGNsphere = 54 . µ m and λ peak AGNcone = 27 . µ m.Finally, the Rayleigh-Jeans tail of the FIR emission ismostly sensitive to the total dust content. In fact, by com-paring the f d = 0 .
08 and f d = 0 . z ∼ QUASAR DATA
To test the results of our model (SPH simulation post-processed with RT calculations), we compare in Fig. 9 ourpredictions from the
AGNcone run ( M UV = − .
97) withmulti-wavelength (NIR to FIR) observations of z ∼ − (cid:46) M UV (cid:46) −
26) quasars (see Table 6).In the NIR, our predicted SEDs are underluminous withrespect to the flux of TNG/GEMINI spectra (grey linesin Fig. 9). This mismatch cannot be solved by decreasingthe dust content, since by assuming f d < .
08 the syntheticSEDs would become underluminous in the FIR with respectto ALMA data. We instead suggest that a better agreementwith observations can be obtained by assuming an extinc-tion curve flatter than the SMC (Gallerani et al. 2010, DiMascia et al in preparation).For what concerns the comparison in the MIR, modelswith α fidUV are in good agreement both with Spitzer/Herschelphotometric data and with the slope/shape of templatesby Hern´an-Caballero & Hatziminaoglou (2011) resemblingSpitzer/IRS spectra. Vice-versa, models with α steepUV are bothunder-luminous with respect to Spitzer/IRAC observationsat λ obs = 24 µ m (namely λ RF ∼ µ m at z ∼
6) and show aslope in the MIR that does not agree with observed spectra.We underline that the model with α steepUV can be possi-bly reconciled with observations if the torus is included. Infact, a dust component with temperature close to sublima-tion ( ∼ . We give a first estimate ofthe impact of the torus emission on our predicted SEDs inSection 5.1 (Fig. 11) and we defer the inclusion of the torusinto our model to a future study.By comparing our predicted SEDs with FIR observa-tions, we note that both models ( f d = 0 . − .
3) providea reasonable match with FIR data, independently on theassumed UV slope (fiducial vs UV-steep). We find that themodels with a larger dust-to-metal ratio f d = 0 . f d = 0 .
08 case we can only explain theless luminous FIR sources.Hereafter, we consider as fiducial the model with α fidUV and f d = 0 . The most massive halo in the
AGNcone run at z = 6 . ∼
70% of the total UV flux. The emission of a greybody at temperature T d and with β d = 2peaks at λ peak = (2 . × ) /T d µ mMNRAS , 1–22 (2020) F. Di Mascia et al.
Figure 6.
Mass fraction (left) and TIR luminosity fraction (right) of dust with a temperature T d > T as a function of the temperature T for the runs noAGN (blue line), AGNsphere (green line),
AGNcone (red line). Results for f d = 0 .
08 are shown.
Figure 7.
Cumulative mass fraction (left) and TIR luminosity fraction (right) of the dust at a distance r from AGN location or moststar forming regions. The lines show the results for noAGN (blue line), AGNsphere (green line) and
AGNcone (red line) with f d = 0 . However, it does not correspond to the most accreting BH,which is instead powering source C, distant ∼
10 kpc from A.Despite having the highest intrinsic UV budget, this sourceis fainter than A in the UV because it is enshrouded by dust:source C is in fact the most luminous IR source of the sys-tem and provides ∼ −
80% ( ∼ (Marshall et al. 2020; Decarli et al. 2017), we foundthat sources A, B and D would be detectable and resolvedwith HST; for what concerns the FIR, given the angularresolution of current ALMA data (i.e. 1” that corresponds We do not consider constraints from MIR observations sinceindividual sources cannot be resolved at these wavelengths as aconsequence of the poor angular resolution. We further refer toVito et al. in preparation for a detailed comparison with X-rayobservations. to ∼ z = 6 . M UV (cid:54) − z ∼ (cid:38) (cid:12) yr − ) and starforming galaxies (e.g. source D). Deeper and higher resolu-tion ALMA data and JWST observations are required tobetter characterize the properties of galaxy companions inthe field of view of z ∼ Given the good agreement between our results and cur-rently available z ∼ MNRAS000
80% ( ∼ (Marshall et al. 2020; Decarli et al. 2017), we foundthat sources A, B and D would be detectable and resolvedwith HST; for what concerns the FIR, given the angularresolution of current ALMA data (i.e. 1” that corresponds We do not consider constraints from MIR observations sinceindividual sources cannot be resolved at these wavelengths as aconsequence of the poor angular resolution. We further refer toVito et al. in preparation for a detailed comparison with X-rayobservations. to ∼ z = 6 . M UV (cid:54) − z ∼ (cid:38) (cid:12) yr − ) and starforming galaxies (e.g. source D). Deeper and higher resolu-tion ALMA data and JWST observations are required tobetter characterize the properties of galaxy companions inthe field of view of z ∼ Given the good agreement between our results and cur-rently available z ∼ MNRAS000 , 1–22 (2020)
R emission of early galaxies: AGN imprints Figure 8.
Intrinsic and processed (observed) SED for the first six runs in Table 4. The first column refers to noAGN , the second to
AGNsphere and the third to
AGNcone , whereas the first row to f d = 0 .
08, and the second to f d = 0 .
3. The solid line shows the observedflux for the reference line of sight, whereas the shaded area brackets the scatter in the observed SED between the six lines of sights usedfor the computation. The intrinsic flux is also shown with a a dot-dashed line. In the runs with AGN, the individual components arealso shown: radiation from stars is denoted with dashed lines, radiation from AGN with dotted lines. In noAGN runs, thin solid linesindicate the flux that would be observed if the galaxy is magnified by a factor µ = 5. Sensitivity bands of JWST and ALMA are shownas yellow and cyan shaded regions respectively. The grey lines indicate the sensitivity reached by the ORIGINS telescope at 5 σ in 1 hr ofobserving time. The colored rectangles and horizontal lines indicate the sensitivity of the two instruments of the SPICA telescope: SMI( λ obs = 27 µ m, red rectangle), and SAFARI, in photometric mapping mode at short (SW, λ obs = 45 µ m, blue rectangle), mid (MW, λ obs = 72 µ m, green rectangle), long (LW, λ obs = 115 µ m, orange line) and very long wavelengths (LLW, λ obs = 185 µ m, violet line).The upper sides of rectangles represent the sensitivity that will be reached by SPICA at 5 σ in 1 hr of observing time. The bottom sideof rectangles represents the maximum sensitivity reachable with SPICA, and it is obtained by considering the confusion limit flux at 3 σ (such a high sensitivity can be reached in the case of follow-up observations). If the confusion limit is reached in less than 1 hr, it isshown as a single line. Figure 9.
Left panel : Comparison between synthetic fiducial SEDs and observations of z ∼ f d = 0 .
08 and f d = 0 .
3, respectively. Grey lines represent dust-reddened z (cid:38) A > .
8, see Table 1 in Gallerani et al. 2010). The spectra are calibrated by usingthe measures of λL λ at 1450 A ◦ provided in Table 1 by Juarez et al. (2009). We further show with orange lines the template obtainedfrom the analysis of 125 quasar spectra at 0 . < z < .
355 taken with Spitzer/IRS (the tree lines correspond to the median SED, the25 th and 75 th percentiles, in the case of the luminous, log( λL ) > .
55, sub-sample by Hern´an-Caballero & Hatziminaoglou (2011)),and with red lines the IR AGN SED derived by Xu et al. (2020) from 42 quasars at z < . z ∼ Right panel : same as left panel, but for the UV-steep AGN SED models.MNRAS , 1–22 (2020) F. Di Mascia et al. source z M UV ReferenceJ1030+0524 6.31 -27.12 [1-2,8]J1048+4637 6.23 -27.60 [1-2,8]J1148+5251 6.43 -27.85 [1-8]J1306+0356 6.03 -26.76 [1-2,8-9]J1602+4228 6.07 -26.85 [1-2,8]J1623+3112 6.25 -26.71 [1-2,8]J1630+4012 6.07 -26.16 [1-2,8]J0353+0104 6.07 -26.56 [8]J0818+1722 6.00 -27.44 [8]J0842+1218 6.08 -26.85 [8,9]J1137+3549 6.01 -27.15 [8]J1250+3130 6.13 -27.18 [8]J1427+3312 6.12 -26.48 [8]J2054-0005 6.04 -26.15 [8]P007+04 6.00 -26.58 [9]P009-10 6.00 -26.50 [9]J0142-3327 6.34 -27.76 [9]P065-26 6.19 -27.21 [9]P065-19 6.12 -26.57 [9]J0454-4448 6.06 -26.41 [9]P159-02 6.38 -26.74 [9]J1048-0109 6.68 -25.96 [9]J1148+0702 6.34 -26.43 [9]J1207+0630 6.04 -26.57 [9]P183+05 6.44 -26.99 [9]P217-16 6.15 -26.89 [9]J1509-1749 6.12 -27.09 [9]P231-20 6.59 -27.14 [9]P308-21 6.23 -26.30 [9]J2211-3206 6.34 -26.65 [9]J2318-3113 6.44 -26.06 [9]J2318-3029 6.15 -26.16 [9]P359-06 6.17 -26.74 [9]J0100+2802 6.33 -29.30 [10]P338+29 6.66 -26.01 [14]J0305-3150 6.61 -26.13 [15]
Table 6.
Quasars used for the comparison with the predictionby our model. Columns indicate: (first) source name, (second)redshift, (third) M UV and (fourth) references for the photometricdata used in the comparison, according to the legend. [1] Galleraniet al. (2010); [2] Juarez et al. (2009); [3] Walter et al. (2003); [4]Bertoldi et al. (2003b); [5] Riechers et al. (2009); [6] Galleraniet al. (2014); [7] Stefan et al. (2015); [8] Leipski et al. (2014); [9]Venemans et al. (2018); [10] (Wang et al. 2016); [11] (Venemanset al. 2012); [12] (Venemans et al. 2017b); [13] (Willott et al.2017); [14] (Mazzucchelli et al. 2017); [15] (Venemans et al. 2016). simulations to make predictions for the proposed OriginsSpace Telescope (OST ; Wiedner et al. 2020). OST cov- OST is a concept study for a 5 . . Figure 10.
Comparison between the SEDs extracted fromsources A, B, C, D in the field of view of the
AGNcone runwith f d = 0 .
3, keeping the same color legend as in Fig. 4. Thetotal SED is instead plotted with a blue solid line. Black pointsindicate rest-frame UV limits from deep HST observations byMarshall et al. (2020), whereas the grey point FIR fluxes for star-forming companion galaxies around quasars from Decarli et al.(2017). ers the wavelength range 2 . − µ m, and is designedto make broad-band imaging (Far-IR Imager Polarimeter,FIP), low resolution ( R ∼ R ∼ − z ∼ . − µ m, and is designed to make high-resolution ( R ∼ − µ m)and mid-infrared (30–37 µ m) broad band mapping, andsmall field spectroscopic and polarimetric imaging at 100,200 and 350 µ m. These are the characteristics of the SpaceInfrared Telescope for Cosmology and Astrophysics (SPICA;e.g. Spinoglio et al. 2017; Gruppioni et al. 2017; Egami et al.2018; Roelfsema et al. 2018), an infrared space mission, ini-tially considered as a candidate for the M5 mission, but can-celled in October 2020 (Clements et al. 2020).The noAGN case is detectable by ORIGINS in fivebands, corresponding to ≈ − µ m rest-frame. ORIGINSwould be able to probe the SED of highly star forming galax-ies (SFR ∼
600 M (cid:12) yr − ) at wavelengths shorter than thepeak wavelength, which is crucial in order to have a soliddetermination of the dust temperature (Behrens et al. 2018;Sommovigo et al. 2020). The noAGN case falls just belowthe SPICA sensitivity threshold. We thus consider the possi-bility of observing lensed galaxies with SPICA; the thin solidSED in the noAGN panels in Fig. 8 accounts for a magnifica-tion factor µ ∼
5. Our results show that highly star forminggalaxies (SFR ∼
600 M (cid:12) yr − ) without an active AGN willbe at the SPICA reach if lensed by a factor µ (cid:38) AGNsphere case, the simulatedrun corresponds to a faint AGN ( M UV = − .
4; X-ray lumi-nosity L X ∼ erg s − ). This kind of sources is not easilydetectable through UV and X-ray observations: (i) less than20 z ∼ M UV = − .
75 have been dis-
MNRAS000
MNRAS000 , 1–22 (2020)
R emission of early galaxies: AGN imprints covered so far (Matsuoka et al. 2018); (ii) none z ∼ L X < × erg s − has been detected so far withChandra (Vito et al. 2019b). Our predictions show that theSED of a faint AGN is instead well above ORIGINS’ sen-sitivities at all wavelengths and also above the sensitivitiesof two SPICA’s bands for all the simulations we performed.This result emphasises the important role that future MIRfacilities would have in studying the faint-end of the UV andX-ray luminosity function in z ∼ AGNcone runs show that quasars with M UV < − ∼ z (cid:38) > M UV < − M UV = −
26 have been detected so far at mmwavelengths at > σ only in two z (cid:62) Combined ALMA data with follow-up JWST and/or ORIGINS observations will be crucial to discoverfaint/obscured AGN and to distinguish them from galaxieswithout an active nuclei. In fact, by comparing the predictedfluxes in ORIGINS band 1 and/or MIRI band at 29 µ m( F µ m ) with the ones in ALMA band 7 F band7 , we find: F µ m ( AGNsphere ) /F band7 ( AGNsphere ) F µ m ( noAGN ) /F band7 ( noAGN ) ≈ − , meaning that we expect a a MIR-to-FIR excess of one or-der of magnitude in the case of a faint AGN host galaxy( AGNsphere ) with respect to a star forming galaxy with-out AGN ( noAGN ). This result shows that by following upwith JWST and/or ORIGINS sources already detected withALMA it will be possible to discriminate between star form-ing galaxies and faint/obscured AGN.We note that, given the limited resolution ( ∼
200 pc)of the hydrodynamical simulations adopted in this work, wecannot resolve the torus ( ∼ . −
10 pc) that is, therefore,not included in our modelling. The presence of a dusty torussurrounding accreting BHs provides an additional source ofMIR emission boosting the MIR excess expected in AGN.This can increase both the detectability of faint quasars witha SPICA-like telescope and the possibility of exploiting thesynergy between ALMA and MIR facilities to unveil dust-obscured AGN. For example, we qualitatively show in Fig.11 how our predicted SEDs would change with the inclusionof the emission from the dusty torus. For this comparison weconsider the
AGNsphere case ( M UV = − . The James Webb Space Telescope is planned to fly on October31, 2021, with 10 years of operation goal. The proposed ORIGINSmission is planned for launch in the early 2030s, so it will ideallycontinue the work of JWST. emission as a single-temperature T dust greybody, with dustmass M dust , and β dust = 2. We consider different models tocover the range in masses (10 − M (cid:12) ) and temperatures(200 − (cid:46)
200 pc scales di-rectly irradiated by the AGN and possibly the IR emissioncoming from high temperature ( T d ∼ −
300 K) regions.We plan to include the torus emission in a consistent way inour model in a future work and to further examine its impacton our results and on the potential of future facilities.
In this work, we have considered a suite of zoom-in cos-mological hydrodynamic simulations of a massive halo( ∼ M (cid:12) ) at z ∼ ∼
600 M (cid:12) yr − ) without BHs (called noAGN run), and two simulations with accreting BHs that accountfor AGN kinetic feedback distributed according to a spheri-cal ( AGNsphere run) and bi-conical (
AGNcone run) geome-try. These two different feedback prescriptions result in dif-ferent SFRs of the host galaxy ( ∼
300 and ∼
200 M (cid:12) yr − in AGNsphere and
AGNcone , respectively), and AGN activity( ∼ ∼
90 M (cid:12) yr − in AGNsphere and
AGNcone , re-spectively).We performed dusty radiative transfer calculations ofthe three runs in post-process by exploiting the code skirt (Baes et al. 2003; Baes & Camps 2015; Camps & Baes 2015;Camps et al. 2016) with the aim of understanding the im-pact of radiative feedback on the observed spectral energydistributions (SEDs) of z ∼ F λ ∝ λ α constrained through observational and theoreticalarguments; (ii) SMC dust properties (grain size distribu-tion and composition); (iii) different total dust mass content(parametrized in terms of the dust-to-metal ratios f d = 0 . f d = 0 . • In dusty galaxies ( M d (cid:38) × M (cid:12) ) a large fraction( (cid:38) • Large UV luminosity variations with viewing angle, canbe at least partially due to the inhomogeneous distributionof dusty gas on scales (cid:38)
100 pc. • Simulations including AGN radiation show the presence
MNRAS , 1–22 (2020) F. Di Mascia et al.
Figure 11.
Predicted SEDs with the inclusion of the dusty torus emission to our models for the
AGNsphere case with f d = 0 .
3. Thegreen dashed line refers to the original model, the red dashed line to the torus and the black solid line to the sum of the two. The torusemission is modelled as a greybody with M dust and T dust as specified in the panel, and β dust = 2. of a clumpy, warm ( ≈ −
300 K) dust component, in ad-dition to a colder ( ≈ −
70 K) and more diffuse dustymedium, heated by stars; warm dust provides up to 50%of the total infrared luminosity, though constituting only asmall fraction ( (cid:46) . z ∼ M UV (cid:46) −
26) quasars,the only class of sources for which multi-wavelength ob-servations, ranging from the optical-NIR to the mm, areavailable so far. For what concerns the intrinsic SEDs, wehave considered two variations for the rest-frame UV band:a fiducial value α fidUV = − . α steepUV = − . • We find a good agreement between simulations andboth MIR (Spitzer/Herschel) and millimetric (ALMA) data,in the case of α fidUV = − .
7. In the rest-frame UV, our pre-dicted SEDs are underluminous with respect to data, sug-gesting peculiar extinction properties (Gallerani et al. 2010,see also Di Mascia et al. in preparation). • The case α steepUV = − . λ obs = 24 µ m and show a slopein the MIR that does not agree with Spitzer/IRS spectra.This discrepancy can be possibly alleviated by adding toour model the emission arising from a dusty torus with T d ∼ • Quasars powered by SMBHs are part of complex, dust-rich merging systems, containing both multiple accretingBHs and star forming galaxies that, because of strong dustabsorption, are below the detection limit of current deepoptical-NIR observations (Mechtley et al. 2012), but ap-pear as SMG companions, consistently with recent shallowALMA data Decarli et al. (2017). Deeper ALMA and futureJWST observations are required to study the environmentin which z ∼ • Highly star forming galaxies (SFR ∼
600 M (cid:12) yr − )without an active AGN will be easily detected by ORIGINS.It will also be able to probe the peak of the dust emission,allowing a solid estimate of the dust temperature in starforming galaxies at high redshift. These galaxies would alsobe detected by a SPICA-like telescope, if lensed by a factor µ (cid:38) • Bright high- z quasars ( M UV < −
26) are detectablewith ORIGINS/SPICA at a signal-to-noise ratio highenough to get high quality spectra even in these very distantsources. • The FIR/MIR flux ratio in star forming galaxies is oneorder of magnitude higher with respect to AGN hosts, evenin the case of low accretion rates ( ˙ M BH ∼ (cid:12) yr − ). Byfollowing up with ORIGINS/SPICA galaxies already de-tected with ALMA it will be possible to unveil faint and/ordust-obscured AGN, whose fraction is expected to be large( > ACKNOWLEDGEMENTS
FD thanks Alessandro Lupi, Laura Sommovigo and MilenaValentini for helpful discussions and Peter Camps for codesupport. We acknowledge fruitful discussions with Paola An-dreani, Sarah Bosman, Eiichi Egami, Carlotta Gruppioni,Francesca Pozzi, Luigi Spinoglio, Christian Vignali. SG ac-knowledges support from the ASI-INAF n. 2018-31-HH.0grant and PRIN-MIUR 2017 (PI Fabrizio Fiore). AF andSC acknowledge support from the ERC Advanced GrantINTERSTELLAR H2020/740120. Any dissemination of re-sults must indicate that it reflects only the author’s view
MNRAS , 1–22 (2020)
R emission of early galaxies: AGN imprints and that the Commission is not responsible for any use thatmay be made of the information it contains. Generous sup-port from the Carl Friedrich von Siemens-Forschungspreisder Alexander von Humboldt-Stiftung Research Award iskindly acknowledged (AF). We acknowledge usage of thePython programming language (Van Rossum & de Boer1991; Van Rossum & Drake 2009), Astropy (Astropy Collab-oration et al. 2013), Cython (Behnel et al. 2011), Matplotlib(Hunter 2007), NumPy (van der Walt et al. 2011), pynbody (Pontzen et al. 2013), and SciPy (Virtanen et al. 2020). DATA AVAILABILITY
Part of the data underlying this article were accessed fromthe computational resources available to the CosmologyGroup at Scuola Normale Superiore, Pisa (IT). The deriveddata generated in this research will be shared on reasonablerequest to the corresponding author.
APPENDIX A: CONVERGENCE OF THE DUSTGRID
The dust content derived from the hydrodynamical simula-tions is distributed in an octree grid with a maximum of 8level of refinement, achieving a maximum resolution of ∼ ∼
200 pc at z = 6. In this Section, we check ifthe number of refinement levels adopted in our fiducial setupis sufficient to achieve converge of the results. We performthree control simulations, in which the maximum refinementlevels are 6, 7 and 9, corresponding to a spatial resolution of937 pc, 469 pc and 117 pc respectively. In Fig. A1, we showthe SED plot for the AGNcone run, adopting f d = 0 .
08 andthe fiducial AGN SED, for the aforementioned values of themaximum refinement levels. The four SEDs mainly differ inthe MIR range (6 − µ m rest-frame). The MIR emissionincreases when increasing the number of refinement levels,because dust around AGN, which is heated to the highesttemperatures, is better resolved. However, the variation be-tween our fiducial model and the model at the highest reso-lution is less than 30% in the MIR band, thus we concludethat the spatial resolution of the dust grid adopted in ourcalculations is sufficient to achieve reasonable numerical con-vergence. APPENDIX B: DUST THERMALSPUTTERING
We have assumed that dust grains with temperature abovea given threshold (
T > K) are destroyed by thermalsputtering (Draine & Salpeter 1979; Tielens et al. 1994; Hi-rashita et al. 2015), as commonly done in simulations (Lianget al. 2019; Ma et al. 2019). However, this dust destructionprocess may be inefficient in the proximity of AGN, becauseof grain charging (Tazaki & Ichikawa 2020; Tazaki et al.2020). To quantify how this assumption affects our resultswe re-run the
AGNcone model with the lower dust content,i.e. f d = 0 .
08 (fifth row in Table 4), after removing thethreshold on the dust temperature. In this case, the mass
Figure A1.
Spectral Energy Distribution of the
AGNcone run( f d = 0 .
08, fiducial AGN SED) for different numbers of the max-imum refinement levels: 6, 7, 8 (fiducial) and 9, corresponding to937 pc, 469 pc, 234 pc (fiducial) and 117 pc, respectively.
Figure B1.
Comparison of the SEDs of the
AGNcone run, as-suming f d = 0 .
08 (red), f d = 0 . f d = 0 .
08 withoutdust sputtering (grey). of emitting dust is a factor ∼ M d = 6 × M (cid:12) ). In Fig. B1 we compare theSEDs obtained with f d = 0 .
08 and f d = 0 . f d = 0 .
08 and f d = 0 . APPENDIX C: ˙ M BH − M UV RELATION
For a radiation efficiency (cid:15) r = 0 .
1, the bolometric luminosity L bol can be related to the BH accretion rate as follows: L bol ≈ . × (cid:18) ˙ M BH M (cid:12) yr − (cid:19) L (cid:12) . (C1) MNRAS , 1–22 (2020) F. Di Mascia et al.
Using the bolometric corrections reported in Table 3, wecan convert the accretion rate into an UV luminosityby multiplying the bolometric luminosity by a factor f UV = L UV /L bol . Then, we adopt the definition of the ABmagnitude m AB = − . F ν − . , where F ν is in cgs units, and we express M UV in terms ofthe product λL λ : M UV = 89 . − . (cid:18) λL λ erg s − (cid:19) , (C2)where ν and λL λ are evaluated at λ = 1450 A ◦ . By com-bining the previous equations we obtain: M UV = − . − . (cid:18) ˙ M BH M (cid:12) yr − (cid:19) . (C3) References
Allevato V., et al., 2016, ApJ, 832, 70Aoyama S., Hou K.-C., Shimizu I., Hirashita H., Todoroki K.,Choi J.-H., Nagamine K., 2017, MNRAS, 466, 105Arata S., Yajima H., Nagamine K., Li Y., Khochfar S., 2019,MNRAS, 488, 2629Asano R. S., Takeuchi T. T., Hirashita H., Inoue A. K., 2013a,Earth, Planets, and Space, 65, 213Asano R. S., Takeuchi T. T., Hirashita H., Nozawa T., 2013b,MNRAS, 432, 637Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A,47, 481Astropy Collaboration et al., 2013, A&A, 558, A33Ba˜nados E., et al., 2014, AJ, 148, 14Ba˜nados E., et al., 2018, Nature, 553, 473Baes M., Camps P., 2015, Astronomy and Computing, 12, 33Baes M., et al., 2003, MNRAS, 343, 1081Bakx T. J. L. C., et al., 2020, MNRAS, 493, 4294Barai P., Gallerani S., Pallottini A., Ferrara A., Marconi A., Ci-cone C., Maiolino R., Carniani S., 2018, MNRAS, 473, 4003Behnel S., Bradshaw R., Citro C., Dalcin L., Seljebotn D., SmithK., 2011, Computing in Science Engineering, 13, 31Behrens C., Pallottini A., Ferrara A., Gallerani S., Vallini L.,2018, MNRAS, 477, 552Berta S., et al., 2013, A&A, 551, A100Bertoldi F., Carilli C. L., Cox P., Fan X., Strauss M. A., BeelenA., Omont A., Zylka R., 2003a, A&A, 406, L55Bertoldi F., et al., 2003b, A&A, 409, L47Bianchi S., Schneider R., 2007, MNRAS, 378, 973Blecha L., Snyder G. F., Satyapal S., Ellison S. L., 2018, MNRAS,478, 3056Bondi H., 1952, MNRAS, 112, 195Bondi H., Hoyle F., 1944, MNRAS, 104, 273Bongiorno A., et al., 2012, MNRAS, 427, 3103Booth C. M., Schaye J., 2009, MNRAS, 398, 53Bruzual G., Charlot S., 2003, MNRAS, 344, 1000Camps P., Baes M., 2015, Astronomy and Computing, 9, 20Camps P., Trayford J. W., Baes M., Theuns T., Schaller M.,Schaye J., 2016, MNRAS, 462, 1057Carilli C. L., Walter F., 2013, ARA&A, 51, 105 In the case of the fiducial SED, f UV ≈ .
29. The results arehowever very similar in the case of the UV-steep SED. In this expression we do not include k-corrections for the dis-tance modulus ( µ ) calculations. At z = 6 .
3, the difference between µ and the k-corrected one µ k is µ = µ k + 2 .
1. Carnall A. C., et al., 2015, MNRAS, 451, L16Carniani S., et al., 2016, A&A, 591, A28Carniani S., et al., 2019, MNRAS, 489, 3939Chabrier G., 2003, Publ. Astr. Soc. Pac., 115, 763Chakrabarti S., Whitney B. A., 2009, ApJ, 690, 1432Chakrabarti S., Cox T. J., Hernquist L., Hopkins P. F., RobertsonB., Di Matteo T., 2007, ApJ, 658, 840Cicone C., et al., 2015, A&A, 574, A14Clements D. L., Serjeant S., Jin S., 2020, Nature, 587, 548Connor T., et al., 2019, ApJ, 887, 171Connor T., et al., 2020, ApJ, 900, 189Cowie L. L., Barger A. J., Bauer F. E., Gonz´alez-L´opez J., 2020,ApJ, 891, 69Cresci G., et al., 2015a, A&A, 582, A63Cresci G., et al., 2015b, ApJ, 799, 82Dalla Vecchia C., Schaye J., 2008, MNRAS, 387, 1431Davies F. B., Hennawi J. F., Eilers A.-C., 2019, ApJL, 884, L19De Young D. S., 1989, ApJL, 342, L59Decarli R., et al., 2017, Nature, 545, 457Di Matteo T., Springel V., Hernquist L., 2005, Nature, 433, 604Draine B. T., 2003a, ARA&A, 41, 241Draine B. T., 2003b, ApJ, 598, 1017Draine B. T., 2003c, ApJ, 598, 1026Draine B. T., Salpeter E. E., 1979, ApJ, 231, 77Draine B. T., et al., 2007, ApJ, 663, 866Dubois Y., Devriendt J., Slyz A., Teyssier R., 2010, MNRAS, 409,985Dubois Y., Gavazzi R., Peirani S., Silk J., 2013, MNRAS, 433,3297Egami E., et al., 2018, Publ. Astr. Soc. Australia, 35, 48Fabian A. C., 1999, MNRAS, 308, L39Fan X., et al., 2000, AJ, 120, 1167Fan X., et al., 2006, AJ, 131, 1203Ferrara A., Salvadori S., Yue B., Schleicher D., 2014, MNRAS,443, 2410Ferrara A., Viti S., Ceccarelli C., 2016, MNRAS, 463, L112Fiore F., Elvis M., McDowell J. C., Siemiginowska A., WilkesB. J., 1994, ApJ, 431, 515Gallerani S., et al., 2010, A&A, 523, A85Gallerani S., Ferrara A., Neri R., Maiolino R., 2014, MNRAS,445, 2848Gallerani S., Fan X., Maiolino R., Pacucci F., 2017a, Publ. Astr.Soc. Australia, 34, e022Gallerani S., et al., 2017b, MNRAS, 467, 3590Garc´ıa-Burillo S., et al., 2016, ApJL, 823, L12Garc´ıa-Burillo S., et al., 2019, A&A, 632, A61Gilli R., et al., 2014, A&A, 562, A67Groves B., Dopita M. A., Sutherland R. S., Kewley L. J., FischeraJ., Leitherer C., Brandl B., van Breugel W., 2008, ApJS, 176,438Gruppioni C., et al., 2016, MNRAS, 458, 4297Gruppioni C., et al., 2017, Publ. Astr. Soc. Australia, 34, e055Hahn O., Abel T., 2011, MNRAS, 415, 2101Haiman Z., 2013, The Formation of the First Massive Black Holes.Astrophysics and Space Science Library, p. 293Hern´an-Caballero A., Hatziminaoglou E., 2011, MNRAS, 414, 500Hickox R. C., Alexander D. M., 2018, ARA&A, 56, 625Hirashita H., Aoyama S., 2019, MNRAS, 482, 2555Hirashita H., Nozawa T., Villaume A., Srinivasan S., 2015, MN-RAS, 454, 1620Hjorth J., Vreeswijk P. M., Gall C., Watson D., 2013, ApJ, 768,173Hopkins P. F., Richards G. T., Hernquist L., 2007, ApJ, 654, 731Hoyle F., Lyttleton R. A., 1939, Proceedings of the CambridgePhilosophical Society, 35, 405Hunter J. D., 2007, Computing in Science Engineering, 9, 90Jiang L., Fan X., Vestergaard M., Kurk J. D., Walter F., KellyB. C., Strauss M. A., 2007, AJ, 134, 1150MNRAS000
1. Carnall A. C., et al., 2015, MNRAS, 451, L16Carniani S., et al., 2016, A&A, 591, A28Carniani S., et al., 2019, MNRAS, 489, 3939Chabrier G., 2003, Publ. Astr. Soc. Pac., 115, 763Chakrabarti S., Whitney B. A., 2009, ApJ, 690, 1432Chakrabarti S., Cox T. J., Hernquist L., Hopkins P. F., RobertsonB., Di Matteo T., 2007, ApJ, 658, 840Cicone C., et al., 2015, A&A, 574, A14Clements D. L., Serjeant S., Jin S., 2020, Nature, 587, 548Connor T., et al., 2019, ApJ, 887, 171Connor T., et al., 2020, ApJ, 900, 189Cowie L. L., Barger A. J., Bauer F. E., Gonz´alez-L´opez J., 2020,ApJ, 891, 69Cresci G., et al., 2015a, A&A, 582, A63Cresci G., et al., 2015b, ApJ, 799, 82Dalla Vecchia C., Schaye J., 2008, MNRAS, 387, 1431Davies F. B., Hennawi J. F., Eilers A.-C., 2019, ApJL, 884, L19De Young D. S., 1989, ApJL, 342, L59Decarli R., et al., 2017, Nature, 545, 457Di Matteo T., Springel V., Hernquist L., 2005, Nature, 433, 604Draine B. T., 2003a, ARA&A, 41, 241Draine B. T., 2003b, ApJ, 598, 1017Draine B. T., 2003c, ApJ, 598, 1026Draine B. T., Salpeter E. E., 1979, ApJ, 231, 77Draine B. T., et al., 2007, ApJ, 663, 866Dubois Y., Devriendt J., Slyz A., Teyssier R., 2010, MNRAS, 409,985Dubois Y., Gavazzi R., Peirani S., Silk J., 2013, MNRAS, 433,3297Egami E., et al., 2018, Publ. Astr. Soc. Australia, 35, 48Fabian A. C., 1999, MNRAS, 308, L39Fan X., et al., 2000, AJ, 120, 1167Fan X., et al., 2006, AJ, 131, 1203Ferrara A., Salvadori S., Yue B., Schleicher D., 2014, MNRAS,443, 2410Ferrara A., Viti S., Ceccarelli C., 2016, MNRAS, 463, L112Fiore F., Elvis M., McDowell J. C., Siemiginowska A., WilkesB. J., 1994, ApJ, 431, 515Gallerani S., et al., 2010, A&A, 523, A85Gallerani S., Ferrara A., Neri R., Maiolino R., 2014, MNRAS,445, 2848Gallerani S., Fan X., Maiolino R., Pacucci F., 2017a, Publ. Astr.Soc. Australia, 34, e022Gallerani S., et al., 2017b, MNRAS, 467, 3590Garc´ıa-Burillo S., et al., 2016, ApJL, 823, L12Garc´ıa-Burillo S., et al., 2019, A&A, 632, A61Gilli R., et al., 2014, A&A, 562, A67Groves B., Dopita M. A., Sutherland R. S., Kewley L. J., FischeraJ., Leitherer C., Brandl B., van Breugel W., 2008, ApJS, 176,438Gruppioni C., et al., 2016, MNRAS, 458, 4297Gruppioni C., et al., 2017, Publ. Astr. Soc. Australia, 34, e055Hahn O., Abel T., 2011, MNRAS, 415, 2101Haiman Z., 2013, The Formation of the First Massive Black Holes.Astrophysics and Space Science Library, p. 293Hern´an-Caballero A., Hatziminaoglou E., 2011, MNRAS, 414, 500Hickox R. C., Alexander D. M., 2018, ARA&A, 56, 625Hirashita H., Aoyama S., 2019, MNRAS, 482, 2555Hirashita H., Nozawa T., Villaume A., Srinivasan S., 2015, MN-RAS, 454, 1620Hjorth J., Vreeswijk P. M., Gall C., Watson D., 2013, ApJ, 768,173Hopkins P. F., Richards G. T., Hernquist L., 2007, ApJ, 654, 731Hoyle F., Lyttleton R. A., 1939, Proceedings of the CambridgePhilosophical Society, 35, 405Hunter J. D., 2007, Computing in Science Engineering, 9, 90Jiang L., Fan X., Vestergaard M., Kurk J. D., Walter F., KellyB. C., Strauss M. A., 2007, AJ, 134, 1150MNRAS000 , 1–22 (2020)
R emission of early galaxies: AGN imprints Jiang L., et al., 2009, AJ, 138, 305Jiang L., et al., 2016, ApJ, 833, 222Jonsson P., Groves B. A., Cox T. J., 2010, MNRAS, 403, 17Juarez Y., Maiolino R., Mujica R., Pedani M., Marinoni S., NagaoT., Marconi A., Oliva E., 2009, A&A, 494, L25Kashikawa N., et al., 2015, ApJ, 798, 28Kormendy J., Ho L. C., 2013, ARA&A, 51, 511Kormendy J., Richstone D., 1995, ARA&A, 33, 581Koss M., Mushotzky R., Treister E., Veilleux S., Vasudevan R.,Trippe M., 2012, ApJL, 746, L22Kroupa P., 2002, Science, 295, 82Kulkarni G., Worseck G., Hennawi J. F., 2019, MNRAS, 488,1035Kurk J. D., et al., 2007, ApJ, 669, 32Laporte N., et al., 2017, ApJL, 837, L21Latif M. A., Ferrara A., 2016, Publ. Astr. Soc. Australia, 33, e051Latif M. A., Schleicher D. R. G., Schmidt W., Niemeyer J., 2013,MNRAS, 433, 1607Leipski C., et al., 2013, ApJ, 772, 103Leipski C., et al., 2014, ApJ, 785, 154Leitherer C., et al., 1999, ApJS, 123, 3Li Y., et al., 2008, ApJ, 678, 41Li J., et al., 2020, ApJ, 889, 162Liang L., et al., 2019, MNRAS, 489, 1397Liang L., Feldmann R., Hayward C. C., Narayanan D.,C¸ atmabacak O., Keres D., Faucher-Gigu´ere C.-A., HopkinsP. F., 2021, MNRAS,Lupi A., Haardt F., Dotti M., Fiacconi D., Mayer L., Madau P.,2016, MNRAS, 456, 2993Lusso E., Worseck G., Hennawi J. F., Prochaska J. X., VignaliC., Stern J., O’Meara J. M., 2015, MNRAS, 449, 4204Lyu J., Rieke G. H., Alberts S., 2016, ApJ, 816, 85Ma X., et al., 2019, MNRAS, 487, 1844Magorrian J., et al., 1998, AJ, 115, 2285Maiolino R., Schneider R., Oliva E., Bianchi S., Ferrara A., Man-nucci F., Pedani M., Roca Sogorb M., 2004, Nature, 431, 533Maiolino R., et al., 2005, A&A, 440, L51Maiolino R., et al., 2012, MNRAS, 425, L66Manti S., Gallerani S., Ferrara A., Feruglio C., Graziani L.,Bernardi G., 2016, MNRAS, 456, 98Marconi A., Risaliti G., Gilli R., Hunt L. K., Maiolino R., SalvatiM., 2004, MNRAS, 351, 169Marshall M. A., et al., 2020, ApJ, 900, 21Matsuoka Y., et al., 2016, ApJ, 828, 26Matsuoka Y., et al., 2018, ApJ, 869, 150Mazzucchelli C., et al., 2017, ApJ, 849, 91McGreer I. D., Fan X., Jiang L., Cai Z., 2018, AJ, 155, 131Mechtley M., et al., 2012, ApJL, 756, L38Mortlock D. J., et al., 2011, Nature, 474, 616Murray N., Quataert E., Thompson T. A., 2005, ApJ, 618, 569Nanni R., Vignali C., Gilli R., Moretti A., Brand t W. N., 2017,A&A, 603, A128Nenkova M., Sirocky M. M., Ivezi´c Z., Elitzur M., 2008, ApJ, 685,147Netzer H., 2015, Annual Review of Astronomy and Astrophysics,53, 365Novak M., et al., 2019, ApJ, 881, 63Nozawa T., Asano R. S., Hirashita H., Takeuchi T. T., 2015,MNRAS, 447, L16Ono Y., et al., 2018, Pub. Astron. Soc. Japan, 70, S10Pacucci F., Volonteri M., Ferrara A., 2015, MNRAS, 452, 1922Pacucci F., Ferrara A., Grazian A., Fiore F., Giallongo E., Puc-cetti S., 2016, MNRAS, 459, 1432Padoan P., Nordlund a., 2011, ApJ, 730, 40Padoan P., Federrath C., Chabrier G., Evans N. J. I., JohnstoneD., Jørgensen J. K., McKee C. F., Nordlund a., 2014, in Pro-tostars and Planets VI. p. 77 Pallottini A., Ferrara A., Bovino S., Vallini L., Gallerani S.,Maiolino R., Salvadori S., 2017, MNRAS, 471, 4128Pezzulli E., Valiante R., Orofino M. C., Schneider R., GalleraniS., Sbarrato T., 2017, MNRAS, 466, 2131Piconcelli E., Jimenez-Bail´on E., Guainazzi M., Schartel N.,Rodr´ıguez-Pascual P. M., Santos-Lle´o M., 2005, A&A, 432,15Pilbratt G. L., et al., 2010, A&A, 518, L1Planck Collaboration et al., 2016, A&A, 594, A13Pontzen A., Rovskar R., Stinson G. S., Woods R., Reed D. M.,Coles J., Quinn T. R., 2013, pynbody: Astrophysics Simula-tion Analysis for PythonPozzi F., et al., 2012, MNRAS, 423, 1909Pringle J. E., 1981, ARA&A, 19, 137Reed S. L., et al., 2015, MNRAS, 454, 3952Ricci C., et al., 2017, Nature, 549, 488Richards G. T., et al., 2003, AJ, 126, 1131Riechers D. A., et al., 2009, ApJ, 703, 1338Roebuck E., Sajina A., Hayward C. C., Pope A., Kirkpatrick A.,Hernquist L., Yan L., 2016, ApJ, 833, 60Roelfsema P. R., et al., 2018, Publ. Astr. Soc. Australia, 35, e030Sazonov S. Y., Ostriker J. P., Sunyaev R. A., 2004, MNRAS, 347,144Schartmann M., Meisenheimer K., Camenzind M., Wolf S., Hen-ning T., 2005, A&A, 437, 861Schawinski K., et al., 2006, Nature, 442, 888Schaye J., et al., 2015, MNRAS, 446, 521Schleicher D. R. G., Palla F., Ferrara A., Galli D., Latif M., 2013,A&A, 558, A59Schneider R., Bianchi S., Valiante R., Risaliti G., Salvadori S.,2015, Astronomy and Astrophysics, 579, A60Shakura N. I., Sunyaev R. A., 1973, A&A, 500, 33Shang C., Bryan G. L., Haiman Z., 2010, MNRAS, 402, 1249Shen X., Hopkins P. F., Faucher-Gigu`ere C.-A., Alexander D. M.,Richards G. T., Ross N. P., Hickox R. C., 2020, MNRAS, 495,3252Sijacki D., Springel V., Di Matteo T., Hernquist L., 2007, MN-RAS, 380, 877Silk J., 2005, MNRAS, 364, 1337Silverman J. D., et al., 2020, ApJ, 899, 154Snyder G. F., Hayward C. C., Sajina A., Jonsson P., Cox T. J.,Hernquist L., Hopkins P. F., Yan L., 2013, ApJ, 768, 168Sommovigo L., Ferrara A., Pallottini A., Carniani S., GalleraniS., Decataldo D., 2020, MNRAS, 497, 956Spinoglio L., et al., 2017, Publ. Astr. Soc. Australia, 34, e057Springel V., 2005, MNRAS, 364, 1105Springel V., Hernquist L., 2003, MNRAS, 339, 289Springel V., Di Matteo T., Hernquist L., 2005, MNRAS, 361, 776Stalevski M., Fritz J., Baes M., Nakos T., Popovi´c L. C., 2012,MNRAS, 420, 2756Stalevski M., Ricci C., Ueda Y., Lira P., Fritz J., Baes M., 2016,MNRAS, 458, 2288Stefan I. I., et al., 2015, MNRAS, 451, 1713Stratta G., Gallerani S., Maiolino R., 2011, A&A, 532, A45Tanaka T., Haiman Z., 2009, ApJ, 696, 1798Tazaki R., Ichikawa K., 2020, ApJ, 892, 149Tazaki R., Ichikawa K., Kokubo M., 2020, ApJ, 892, 84Teyssier R., Moore B., Martizzi D., Dubois Y., Mayer L., 2011,MNRAS, 414, 195Tielens A. G. G. M., McKee C. F., Seab C. G., Hollenbach D. J.,1994, ApJ, 431, 321Todini P., Ferrara A., 2001, MNRAS, 325, 726Tornatore L., Borgani S., Dolag K., Matteucci F., 2007, MNRAS,382, 1050Trayford J. W., et al., 2017, MNRAS, 470, 771Urry C. M., Padovani P., 1995, Publ. Astr. Soc. Pac., 107, 803Valiante R., Schneider R., Bianchi S., Andersen A. C., 2009, MN-RAS, 397, 1661MNRAS , 1–22 (2020) F. Di Mascia et al.
Vallini L., Ferrara A., Pallottini A., Gallerani S., 2017, MNRAS,467, 1300Van Rossum G., Drake F. L., 2009, Python 3 Reference Manual.CreateSpace, Scotts Valley, CAVan Rossum G., de Boer J., 1991, CWI Quarterly, 4, 283Vanden Berk D. E., et al., 2001, The Astronomical Journal, 122,549–564Venanzi M., H¨onig S., Williamson D., 2020, ApJ, 900, 174Venemans B. P., McMahon R. G., Warren S. J., Gonzalez-SolaresE. A., Hewett P. C., Mortlock D. J., Dye S., Sharp R. G., 2007,MNRAS, 376, L76Venemans B. P., et al., 2012, ApJL, 751, L25Venemans B. P., et al., 2013, ApJ, 779, 24Venemans B. P., et al., 2015, MNRAS, 453, 2259Venemans B. P., Walter F., Zschaechner L., Decarli R., De RosaG., Findlay J. R., McMahon R. G., Sutherland W. J., 2016,ApJ, 816, 37Venemans B. P., et al., 2017a, ApJ, 845, 154Venemans B. P., et al., 2017b, ApJL, 851, L8Venemans B. P., et al., 2018, ApJ, 866, 159Vignali C., et al., 2018, MNRAS, 477, 780Virtanen P., et al., 2020, Nature Methods, 17, 261Vito F., Gilli R., Vignali C., Comastri A., Brusa M., CappellutiN., Iwasawa K., 2014, MNRAS, 445, 3557Vito F., et al., 2018, MNRAS, 473, 2378Vito F., et al., 2019a, A&A, 628, L6Vito F., et al., 2019b, A&A, 630, A118Volonteri M., Haardt F., Madau P., 2003, ApJ, 582, 559Wada K., Schartmann M., Meijerink R., 2016, ApJL, 828, L19Walter F., et al., 2003, Nature, 424, 406Walter F., Riechers D., Cox P., Neri R., Carilli C., Bertoldi F.,Weiss A., Maiolino R., 2009, Nature, 457, 699Wang R., et al., 2013, ApJ, 773, 44Wang R., et al., 2016, ApJ, 830, 53Wang F., et al., 2018, ApJL, 869, L9Wang F., et al., 2019, ApJ, 884, 30Wang F., et al., 2021, ApJL, 907, L1Weinberger R., et al., 2018, MNRAS, 479, 4056Weingartner J. C., Draine B. T., 2001, ApJ, 548, 296Werner M. W., et al., 2004, ApJS, 154, 1Wiedner M. C., et al., 2020, arXiv e-prints, p. arXiv:2012.02731Wiersma R. P. C., Schaye J., Smith B. D., 2009, MNRAS, 393,99Willott C. J., McLure R. J., Jarvis M. J., 2003, ApJL, 587, L15Willott C. J., et al., 2007, AJ, 134, 2435Willott C. J., et al., 2010, AJ, 139, 906Willott C. J., Bergeron J., Omont A., 2017, ApJ, 850, 108Wiseman P., Schady P., Bolmer J., Kr¨uhler T., Yates R. M.,Greiner J., Fynbo J. P. U., 2017, A&A, 599, A24Wu X.-B., et al., 2015, Nature, 518, 512Wyithe J. S. B., Bolton J. S., 2011, MNRAS, 412, 1926Xu J., Sun M., Xue Y., 2020, ApJ, 894, 21Xue Y. Q., et al., 2011, ApJS, 195, 10Younger J. D., Hayward C. C., Narayanan D., Cox T. J., Hern-quist L., Jonsson P., 2009, MNRAS, 396, L66Zafar T., Watson D., Fynbo J. P. U., Malesani D., Jakobsson P.,de Ugarte Postigo A., 2011, A&A, 532, A143Zafar T., et al., 2018, MNRAS, 480, 108Zinn P. C., Middelberg E., Norris R. P., Dettmar R. J., 2013,ApJ, 774, 66Zubovas K., Nayakshin S., King A., Wilkinson M., 2013, MNRAS,433, 3079van der Walt S., Colbert S. C., Varoquaux G., 2011, Computingin Science Engineering, 13, 22This paper has been typeset from a TEX/L A TEX file prepared bythe author. MNRAS000