A&A593,A109(2016) Astronomy DOI:10.1051/0004-6361/201526880 & (cid:13)c ESO2016 Astrophysics Disentangling star formation and AGN activity in powerful infrared luminous radio galaxies at 1 < z < 4 G.Drouart1,2,3,B.Rocca-Volmerange3,C.DeBreuck4,M.Fioc3,M.Lehnert3,N.Seymour2,D.Stern5,andJ.Vernet4 1 DepartmentofEarthandSpaceScience,ChalmersUniversityofTechnology,OnsalaSpaceObservatory,43992Onsala,Sweden 2 InternationalCentreforRadioAstronomyResearch,CurtinUniversity,Perth,Australia e-mail:[email protected] 3 Institutd’AstrophysiquedeParis,98bisboulevardArago,75014Paris,France 4 EuropeanSouthernObservatory,KarlSchwarzschildStraße2,85748GarchingbeiMünchen,Germany 5 JetPropulsionLaboratory,CaliforniaInstituteofTechnology,MailStop169-221,Pasadena,CA91109,USA Received3July2015/Accepted31May2016 ABSTRACT High-redshiftradiogalaxiespresentsignsofbothstarformationandAGNactivity,makingthemidealcandidatestoinvestigatethe connectionandcoevolutionofAGNandstarformationintheprogenitorsofpresent-daymassivegalaxies.Wemakeuseofasample of11powerfulradiogalaxiesspanning1<z<4whichhavecompletecoverageoftheirspectralenergydistribution(SED)fromUV toFIRwavelengths.UsingHerscheldata,wedisentangletherelativecontributionoftheAGNandstarformationbycombiningthe galaxyevolutioncodePÉGASE.3withanAGNtorusmodel.Wefindthatthreecomponentsarenecessarytoreproducetheobserved SEDs:anevolvedandmassivestellarcomponent,asubmmbrightyoungstarburst,andanAGNtorus.Wefindthatpowerfulradio galaxiesformatveryhigh-redshift,butexperienceepisodicandimportantgrowthat1<z<4asthemassoftheassociatedstarburst variesfrom5to50%ofthetotalmassofthesystem.Thepropertiesofstarformationdifferfromsourcetosource,indicatingno generaltrendofthestarformationpropertiesinthemostinfraredluminoushigh-redshiftradiogalaxiesandnocorrelationwiththe AGNbolometricluminosity.Moreover,wefindthatAGNscatteredlighthaveaverylimitedimpactonbroad-bandSEDfittingon oursample.Finally,ouranalysisalsosuggestsawiderangeinoriginsfortheobservedstarformation,whichwepartiallyconstrain forsomesources. Keywords. galaxies:active–galaxies:evolution–galaxies:high-redshift–galaxies:starformation–quasars:general– galaxies:starburst 1. Introduction mid-infrared (MIR) emission, which showed hot dust close to the supermassive black hole (Haasetal. 2008; DeBreucketal. High-redshift radio galaxies (HzRG) are among the most lu- 2010). In the orientation-based Unification scheme (e.g. minous objects in the Universe and are the target of a consid- Antonucci 1993), the hot dust near central active galactic nu- erable number of investigations aiming to understand galaxy clei (AGN) is completely obscured by material along the line evolution(Miley&DeBreuck2008,forareview).Whilepow- of sight in type 2 AGN, and becomes gradually less obscured erful radio emission betrays the presence of a supermassive when transitioning to type 1 AGN (e.g., Leipskietal. 2010; black hole, near-infrared (NIR) images reveal that their host Drouartetal. 2012; Rawlingsetal. 2013). Spitzer observations galaxies are amongst the brightest galaxies in the Universe also revealed that high-redshift radio galaxies are preferen- at any redshift (the K-z relation, e.g., Lilly&Longair 1984; tiallyfoundindenserenvironments,consistentwithgalaxyclus- Inskipetal. 2002). Modelling the K-band emission suggests ters in formation (e.g., Mayoetal. 2012; Galametzetal. 2012; that radio galaxies are extremely massive galaxies, 1012M Wylezaleketal.2013a,b). (cid:12) (Rocca-Volmerangeetal.2004).Opticalimagesrevealthatmost The availability of sensitive submm bolometer arrays, such radio galaxies present a surprising alignment between the ul- as SCUBA or LABOCA (Hollandetal. 1999; Siringoetal. traviolet(UV)/optical emission and the radio jet axis (e.g., 2009)allowedthefirstdetectionofdustystarformationathigh- Chambersetal. 1987; McCarthyetal. 1987; McCarthy 1993; redshift through the negative K-correction (Blainetal. 1999). Pentericcietal.1999,2001).Moregenerally,thediscoverythat Extensive surveys at z > 3 showed that radio galaxies exhibit the mass of the central supermassive black hole is related to high submm flux, indicating vigorous, on-going star formation its host galaxy and dark matter halo (Magorrianetal. 1998; (e.g., Archibaldetal. 2001; Stevensetal. 2003; Reulandetal. Gebhardtetal. 2000; Ferrarese&Merritt 2000; Häring&Rix 2004), similar to submm galaxies (SMGs, e.g., Borysetal. 2004) implicates radio galaxies as important objects for under- 2003; Coppinetal. 2006; Weißetal. 2009) at similar redshift. standingtheformationandevolutionofmassivegalaxies. However, single-dish observations are limited in spatial reso- The Spitzer Space Telescope (Werneretal. 2004) provided lution (typical 20(cid:48)(cid:48)), and only a handful of sources have been host galaxy masses of these objects, that place them among observed with interferometric facilities (DeBreucketal. 2005; the most massive galaxies in the Universe (Seymouretal. Ivisonetal.2012),whichmakesgeneralconclusionsontheen- 2007). Spitzer also allowed the first characterisation of their tiresampleofhigh-redshiftradiogalaxiesdifficult. ArticlepublishedbyEDPSciences A109,page1of26 A&A593,A109(2016) The Herschel Space Observatory (Pilbrattetal. 2010) cov- It is difficult to simultaneously characterise the UV-to- ersthe60−700µmrange,allowingthefirstsensitiveexploration submmSEDwithAGNandgalaxyevolutioncodes:(i)themul- oftheentireinfrared(IR)spectralenergydistribution(SED)of tiplecomponentsmeananincreasednumberoffreeparameters; galaxies in the distant Universe. This part of the SED is of ut- (ii)anefficientfittingprocedureischallenging;(iii)asampleof most importance because most of the energy of star-forming sourcesneedstobeconstrainedthatincludestheAGNandhost galaxies is emitted at IR wavelengths. The IR comes from the contribution at similar levels; and (iv) homogeneous coverage emission of dust that has absorbed a significant fraction of the inwavelengthforthesesamplesisnecessary.Forthesereasons, UVandopticalemissionfromyoungstars(e.g.,Chary&Elbaz only a few attempts exist in the literature (although the num- 2001; Sanders&Mirabel 1996). AGN emission by dust also berisincreasingthankstotheavailabilityofnewtechniquesand contributes mainly in the IR, which can make it difficult to de- data).Mostoftheseattemptssufferfromsimplifyingonecom- termine the origin of the IR emission. While the submm emis- ponenttostudytheSED,eitherthehostgalaxyortheAGN,de- sion is generally associated with cold dust and therefore with pendingonthesampleorthescientificgoalsofthestudy(e.g., star formation, the complete IR SED allows for a simultane- Trichasetal. 2012; Feltreetal. 2013; Rocca-Volmerangeetal. ous characterisation of the AGN and star formation emission 2013; Banerjietal. 2014; Karouzosetal. 2014; Leipskietal. (e.g., Tadhunteretal. 2007; Dickenetal. 2009; Mullaneyetal. 2014;Cieslaetal.2015). 2011;Feltreetal.2012;Rosarioetal.2012;Drouartetal.2014; This paper is structured as follows. We present the sample Leipskietal. 2014). The transition regime between dust heated andthedatainSect.2.Section3describesthefittingprocedure by the AGN and dust heated by star formation is at 40-60µm and the model used in this analysis. We present the results in andisthereforeimportanttoobservebecauseitaffectsthetotal Sect.4anddiscusstheirimplicationsinSect.5.Throughoutthis infraredemissionandcharacterisationofeachcomponent. paper,weassumethestandardconcordancecosmologicalmodel As the wavelength coverage of galaxies increased over (H0 = 70kms−1Mpc−1, ΩΛ = 0.7, ΩM = 0.3, Spergeletal. the past decade, more advanced models were developed to 2003). predict the emission of galaxies and AGN for given evo- lutionary scenarios and/or morphologies. The literature now offers a large number of these models; we only mention here the models using stellar libraries and scenarios of 2. Sampleanddata evolution to predict the emission of the integrated emis- sion of the galaxy over a wide wavelength range (e.g., To reliably separate the different emission components of ra- Leitherer&Heckman 1995; Fioc&Rocca-Volmerange 1997; dio galaxies, a wide and homogeneous coverage of the SED Maraston 1998; Bruzual&Charlot 2003). As the high-redshift is of prime importance. The HeRGÉ sample (Herschel Radio infrared observations are becoming increasingly common, the Galaxy Evolution sample; Drouartetal. 2014) is ideally suited models have extended their coverage to take the dust emis- to estimate the relative contributions of AGN and their host sion into account (e.g. Leithereretal. 1999; Burgarellaetal. galaxies thanks to the existing Spitzer, Herschel and submm 2005;Maraston2005;daCunhaetal.2008;Grovesetal.2008; data. We briefly summarise here the criteria selected to build Nolletal.2009, forarecentreview,seeConroy2013)Themain thissample.Theradiogalaxieswereselectedtocoverhomoge- characteristicofPÉGASE.3 weuseherearetocoherentlypre- neouslytheradioluminosity-redshiftplane,applyingthecriteria dictattenuatedmetal-dependentSEDsextendedtothedustemis- L3GHz > 1026WHz−1, where L3GHz is the total luminosity at a sion in coherence with metal abundances. The code has been rest-frame frequency of 3GHz (Table 1; Seymouretal. 2007). usedinRV13andwillbedetailedinFiocetal.(inprep.).Simul- To complete our SED coverage from UV to submm, we there- taneously, numerous AGN torus models have been developed foredefinedasubsamplewiththefollowingcriteria.Eachsource to characterise nearby AGN. These models predict the UV to (i) must have at least four Herschel detections; (ii) must have far-infrared (FIR) SED as a function of torus geometry, solv- atleasttwobroad-bandphotometricobservationsbracketingthe ing the radiative transfer equations. Different approaches have discontinuity at 4000 Å restframe; and (iii) the filter response beenusedtodescribethetorusproperties,particularlytorepro- curves should avoid strong emission lines as much as possible. duce the 9.7µm silicate feature. The two model types assume We aim with this selection to obtain more than ten data points either a continuous dust distribution or a clumpy distribution from UV to submm wavelengths as homogeneously as possi- (e.g., Pier&Krolik 1992; Hönigetal. 2006; Fritzetal. 2006; ble to allow for a reliable fit. In total, 11 sources remain out of Nenkovaetal. 2008; Stalevskietal. 2012, for a recent review, 70 from the HeRGÉ sample (see Fig. 1). In addition, a polari- seeAntonucci2012). sation measurement at restframe 1500 Å is beneficial, if avail- In this paper, we present the first in-depth analysis of the able,tomeasurethescatteredAGNlight.Ofourelevensources, different emission components present for a sample of power- seven have polarisation measurements (from broad-band pho- ful radio galaxies, handling data continuously covering UV to tometry or spectroscopy, see Table 1 and the detailed discus- submmwavelengths.Disentanglingthedifferentcomponentsof sion in Sect. 4.4). The sources span 1 < z < 4, with a median a galaxy is particularly challenging in the UV-NIR domain be- redshift (cid:104)z(cid:105) = 2.5. The requirement of Herschel detections bi- cause several components can contribute simultaneously to the ases this sample towards the brightest IR emitters, with LIR > emissionandtohighlyvariableproportions.TheFIRbringsan 4×1012L (Drouartetal.2014).WhenLIRistakenintoaccount, (cid:12) essential piece of information as, thanks to the energy balance these radio galaxies cannot be distinguished from ULIRGs de- between absorption and emission, constraints on each compo- tected by other samples (see Fig. 1). The two highest redshift nent are added. We aim in this paper to take advantage of the sources,4C41.17andTNJ2007-1316areextensivelyanalysed broadwavelengthcoverageofboththedataandmodelstoreli- in Rocca-Volmerangeetal. (2013, hereafter RV13). They had ably disentangle the different components of a sample of high- been selected on the basis of their weak AGN contribution in redshiftradiogalaxies.Toinvestigatethedifferentcomponents, theUVandopticaltotestthemethod,whichisfurtherimproved we use two models, one to model the host galaxy and one to in this paper. We now extend the sample to 11 sources with a modeltheAGN. widerrangeofAGNcontributions. A109,page2of26 G.Drouartetal.:DisentanglingSFandAGNinluminousHzRGs Table1.Summarisedcharacteristicsofoursample. Name z Ref. RA Dec Morph.peculiarities Radiojetaxis logL Polarisation Ref. 3GHz spec. [J2000] [J2000] [WHz−1] [%] 3C368 1.132 a 18:05:06.37 +11:01:33.1 twocomponents aligned 27.63 7.6±0.9(V) g 3C470 1.653 b 23:58:36.00 +44:04:45.0 twocomponents misaligned 28.24 ... ... MRC0324-228 1.894 c 03:27:04.54 −22:39:42.1 oneclosecompanion ... 28.49 6.5+3.8(R) h −3.3 PKS1138-262 2.156 c 11:40:48.38 −26:29:08.8 multiplecompanions aligned 28.14 ... ... MRC0406-244 2.427 d 04:08:51.46 −24:18:16.4 multiplecomponents aligned 28.11 2.5+2.3(I) h −2.5 MRC2104-242 2.491 e 21:06:58.28 −24:05:09.1 twocomponents misaligned 27.88 ... ... USS0828+193 2.572 c 08:30:53.42 +19:13:15.7 multiplecomponents aligned 27.47 10.0±2.0(spec) i 4C28.58 2.891 c 23:51:59.20 +29:10:29.0 multiplecomponents aligned 27.89 ... ... USS0943-242 2.923 c 09:45:32.73 −24:28:49.7 oneclosecompanion aligned 27.95 6.6±0.9(spec) i 4C41.17a 3.792 c 06:50:52.23 +41:30:30.1 multiplecomponents aligned 28.17 <2.4(spec) j TNJ2007-1316a 3.840 f 20:07:53.26 −13:16:43.6 oneclosecompanion ... 27.79 ∼3(spec) k Notes. Columns2and3indicatethespectroscopicallyconfirmedredshiftsandtheirreferences.Column6referstomorphologicalpeculiarities at optical/NIR wavelengths: companion: source close to the object identified as the radio galaxy (few arcsecond) and components: potential substructureswhendetectedintheimages.Column7indicatesthealignmentoftheradiojetaxiscomparedtotheoptical/NIRextensionwhen detectable in the images. Column 8 lists the restframe 3GHz luminosity from DeBreucketal. (2010). The penultimate column indicates the polarisationandtheband. References. (a)Meisenheimer&Hippelein(1992);(b)Hewitt&Burbidge(1991);(c)Roettgeringetal.(1997);(d)McCarthyetal.(1996);(e) McCarthyetal.(1990);(f)Bornancinietal.(2007);(g)diSeregoAlighierietal.(1989);(h)Buchard(2008);(i)Vernetetal.(2001);(j)Deyetal. (1995);(k)Rocca-Volmerangeetal.(2013). 2.1. BuildingtheSEDfrombroad-bandphotometry Weanalysedgalaxieswithmulti-wavelengthdata,thereforethe consistency of the photometry throughout the entire spectral rangeisessential.Specifically,thedatafromtheoptical/NIRand FIR domain have very different spatial resolutions. The optical datahavesub-arcsecresolution(HST),whiletheHerschelbeam can reach ∼35(cid:48)(cid:48) (FWHM). For this reason, we report the total fluxwhenwecompilephotometryfromtheliterature.Thetotal fluxwascalculatedbyapplyingcorrectionstoaperturephotom- etrybasedontheusermanualorprimaryreferenceofeachcor- responding instrument, or by selecting a large aperture, that is, 0 >30kpc(>4(cid:48)(cid:48))or64kpc(∼8(cid:48)(cid:48)). Basedonthefilterresponseandredshiftofthesource,strong emission lines (i.e., Lyα, [OII], [OIII], Hα...) often fall in a band, increasing the measured flux relative to the pure contin- uumcontribution(byupto∼40%).Ifspectroscopicobservations Fig. 1. Total infrared luminosity versus redshift for several samples areavailableforthecontaminatedbands,wesubtractedtheesti- availableintheliterature.ThelightgreendotsshowtheCOSMOSsam- mated flux of the line based on the spectroscopy. All reported plefromKartaltepeetal.(2010),makinguseofSpitzerdata.Theblue fluxes in the tables are corrected and therefore correspond to squares denote the GOODS samples from Elbazetal. (2011) and the pure continuum emission for the remainder of this paper (see darkgreendotsrepresenttheselectionfromSymeonidisetal.(2013), Tables D.1−D.11). Only MRC 0406-244 does not have the in- both making use of Herschel data. The red circles indicate the sub- formationonlinecontamination,andfluxesareusedasgiven. sample considered in this paper. The red diamonds show the sample ofHzRGspresentedinDrouartetal.(2014). 2.2. Notesonindividualsources The data set for each source is mainly rich and of high quality. Especially,morphologicalpeculiaritiescanbeisolatedandhelp to understand the overall picture in term of evolution in each emission from specific components of the radio galaxy. The fi- source. We report the notes on each source in the Appendix C. nalgoalistodisentangletheAGNandstellaremissionthrough Wedescribethecombinedhigh-resolutionimageandradiocon- SED fitting. We made use of two models, the AGN model of tours (when electronically available) in Sect. 5 along with the Fritzetal.2006andthestellarpopulationmodelPÉGASE.3for discussionontheevolutionarystatusofthegalaxy. thehostgalaxy(Fiocetal.,inprep.).Wepresentthefittingpro- cedure,consideredmodelsandtemplatelibrariesforeachcom- ponent (AGN and galaxy) in the following sections. We also 3. Models,fittingprocedure,andtemplatelibraries explain in detail the creation of the new templates and the as- The broad wavelength coverage from UV to FIR allows us to sumptions made on both models to limit our total parameter probemultipleregimesofenergyandthereforetobetterisolate space. A109,page3of26 A&A593,A109(2016) Table2.Fixedandfittedparametersforourlibraryofevolved,SB,andAGNtemplates. Evolved SB AGN Hubbletype[Sa,S0,E,E2] SSPscenario Openingangleθ=[40◦] Fixedparameters IMF[Kroupa] IMF[Kroupa] Radialdensityprofile[α=0] Azimuthaldensityprofile[γ=2] Normalisation Normalisation log10fA2G0Nµm=[−3,−2,−1.7,−1,−0.7,0,0.7,1,2,3] Fittedparameters Age,tevolved[<tUniverse,z] Age,tSB[<tevolved] Inclinationi=[40,50,60,70,80,90] ColumndensityfactorK=[1,10] SizeY=[10,60,150] Initialmet.Zinit=[0,0.0005,0.001,0.005,0.01] Opacityτ9.7µm=[0.1,1,10] Notes.WecombinedtheAGNandSBlibrariestocreatehybridtemplates.Consideredrangesandvaluesinthefittingarereported.Inparticular, theageofthestarburst,t ,mustbeyoungerthantheageoftheevolvedcomponent,t ,whichinturnmustbeyoungerthantheageofthe SB evolved Universeattheconsideredredshift.Thesinglestellarpopulation(SSP)iswithoutinfall,hasnogalacticwinds,andisinstantaneous. 3.1. Fittingprocedure stochasticheatingofgrains.ThecontinuoussyntheticSEDsthat extend from the far-UV to submm wavelengths, were built by ThefittingprocedureispresentedinRV13.Thecodeallowsfit- minimising the input parameter numbers to avoid large degen- tingofuptotwostellarcomponentssimultaneouslybasedonχ2 eracies in the solutions. We characterised each star formation minimisationbyautomaticallyexploringalibraryoftemplates. scenariowithonlyfourfreeparameters:starformationrate,ini- ToadaptthefittingprocedurefromPÉGASE.2tothenewcode tialmassfunction(IMF),inflowsandoutflows.WerefertoFioc PÉGASE.3,wemadethefollowingmodifications:(i)weadded et al. (in prep.) for a more detailed presentation of the code anestimateoftheuncertaintiesonthefittedquantities;(ii)modi- PÉGASE.3 and its detailed documentation, which is available fiedthegraphicaldisplayingproceduretoincludetheAGNcom- onthePÉGASEwebsite1.Inthepresentversion,usedhereafter ponent; (iii) extended the filter database; and (iv) updated the asinRV13,thestarformationparametersbyspectraltypesare compatibilitywiththenewoutputsofPÉGASE.3. taken from LeBorgne&Rocca-Volmerange (2002), fitting the Adding the AGN component was challenging in the con- HubblesequencegalaxypropertiesinthelocalUniverse.These textofthetwo-componentfittingprocedurewithoutmodifyinga parameters are set to reproduce observed local galaxies from largepartofthecoreprocedure.Wesolvedthisdifficultybycre- SDSS (Tsalmantzaetal. 2009, 2012), including evolution with atinghybridtemplatescomposedofthesumoftheevolvingstar- aformationatzform =10.Wefocusedourlibraryonfourgalaxy burst (SB) and AGN components (see Sect. 3.2). The principal types, E, E2, S0, and Sa to test our sensitivity to the galaxy problem in adding new components stems from the degenera- typewhenmodellingtheradiogalaxyhost.Despitetheirnames, cies induced by the increase in the number of free parameters. thesetemplatesrefertodifferentstarformationhistoriesandnot We addressed this specific problem by combining a broad and to specific morphologies, see LeBorgne&Rocca-Volmerange homogeneousSEDcoverage(with12−19datapoints)andcare- (2002). fullychoosingmodelsandparameters:theself-consistenttreat- Thestarburst(SB)templatesaredefineddifferentlythanthe mentoftheopticalandIRemissionoftheAGNandPÉGASE.3 evolved, main component. They are single stellar populations models(calculatedwithradiativetransfercodes),enabledusto (SSPs), assuming a short formation in one step of time, that leave a significant number of parameters free with only lim- is, 1 Myr (δ function). SSPs do not take into account any in- ited parameter degeneracies (see Sect. 5 for a discussion of the fall or galactic winds. We discuss the effect of this 1 Myr for- uncertainties). mation assumption in detail in Sect. 4.2. We adopted a Kroupa IMF(Kroupa2001)andarangeofinitialmetallicitiesofthegas Our approach uses grid-model fitting and requires brute- (Z ). RV13 showed that the column density plays a large role forcecalculationthroughtheparameterspace.Wethereforeset init intheobservedSED. some parameters to physically justified values and took repre- Weexploredthesamequantitiesfortheevolvedpopulation sentativemodelstoensurethatthecoderemainedexecutableon with the SB component, that is, the age and mass of the stel- a desktop machine in a reasonable computing time. More ex- lar populations, but we also added the initial metallicity and tensiveexplorationoftheparameterspacewouldrequireanew the column density factor K (see RV13, Sect. 3). Table 2 re- fitting algorithm using a Bayesian approach and Monte Carlo portstherangeofvaluesconsideredinourlibrary.Wehereex- samplingoftheparameterspace,butthisisbeyondthescopeof plore the age and mass of the starburst to determine whether this paper because it requires a thourough redesign of the core the ten times higher column density preferentially found in ra- procedure.Althoughnotcomplete,thecoverageoftheparame- ter space is still satisfying, with >107 templates tested for each diogalaxiesbyRV13isconfirmedforalargersample.Wealso note that the fitting procedure requires the evolved component source. to be younger than the age of the Universe at a given redshift (t <t )andtheyoungcomponenttobeyoungerthan evolved Universe,z 3.1.1. PÉGASEmodelandlibrary theevolvedcomponent(tSB <tevolved;seealsoTable2). We used the PÉGASE.3 code to generate a library of galaxies 3.1.2. AGNmodelandlibrary and starbursts.The newcode PÉGASE.3is acoherent spectro- chemical evolution model predicting simultaneously the metal- The AGN component was incorporated using the AGN model licityenrichment,thecorrespondingSED,andtheattenuationor fromFritzetal.(2006).Tobrieflydescribethismodel:acentral, emission by dust grain models calculated with radiative trans- ferMonteCarlosimulations,takingintoaccountscatteringand 1 www2.iap.fr/pegase A109,page4of26 G.Drouartetal.:DisentanglingSFandAGNinluminousHzRGs point-like source emits an SED composed of a sum of power Table 3. Best χ2 for one (evolved only), two (evolved and starburst), laws.Thetorusisdefinedinsphericalcoordinatesbyitsgeom- andthree(evolved,starburst,andAGN)components. etry (size and opening angle), dust profile density (radial and azimuthal),opacityat9.7µmandinclinationwithrespecttothe Name 1comp. 2comp. 3comp. observer. The torus is divided into cells assuming an axisym- #freeparameters 2 6 10 metric geometry and the radiative transfer equations are solved 3C368 3615[17] 101[13] 73[9] usingtheΛ-iterationtechnique(seeFritzetal.2006,foracom- 3C470 4000[12] 71[8] 32[4] pletedescriptionofthemodel).Feltreetal.(2012)showedthat MRC0324-228 2152[11] 57[7] 25[3] theFritzetal.(2006)modelisdirectlycomparabletotheclumpy PKS1138-262 14896[14] 834[10] 79[6] modelfromNenkovaetal.(2008)foracertainsetofparameters MRC0406-244 1208[14] 65[10] 19[6] (onlypartiallycoveredinourcurrentlibrary). MRC2104-242 167[10] 20[6] 9[2] We adopted the Fritzetal. (2006) model for the present USS0828+193 5235[12] 619[8] 36[4] study, leaving the comparison with clumpy models for further 4C28.58 247[13] 87[9] 29[5] publications. Moreover, a precise characterisation of the torus USS0943-242 250[13] 16[9] 5[5] propertiesisbeyondthescopeofthispaper,becauseitrequires 4C41.17 491[16] 39[12] 11[8] spectroscopyaroundthesilicatefeatureat9.7µmrestframe.We TNJ2007-1316 148[13] 73[9] 9[5] focusonthelargestimpactparameters,thatis,thesize,theopac- ity, and the inclination of the torus. The size, Y, is the ratio Notes. We also indicate the number of degrees of freedom (d.o.f.) in Rmax/Rmin, where Rmin is set as the sublimation radius of the brackets,d.o.f.=N−p−1whereNisthenumberofdatapointsandp dust (which, in turn depends on the luminosity of the central isthenumberoffreeparameters. source).Theopacity,τ ,correspondstotheintegratedopac- 9.7µm ityat9.7 µmalongtheequatorialaxis(i = 90◦).Thedustden- brute-force calculation through our parameter space, therefore sity(radialandazimuthal)ofthetorusiscalculatedwithrespect each fit possesses a χ2 value, that indicates the goodness of fit. tothesize,opacityandprofilechosenbytheuser.Thesevalues We can therefore consider the χ2 distribution to gain insights were set to their default values here, radially constant and az- into the confidence interval of the adopted free parameters. We imuthallydecreasingasthesquareofthealtitude(theyalsohave a limited effect on the SED). The inclination, i, is the orienta- used the description from Numerical recipes (Pressetal. 1992, Chap. 15) to estimate the uncertainties on the parameters, us- tion of the torus with respect to the observer. An edge-on view corresponds to i = 90◦, and i = 0◦, corresponds to a face-on ing the complementary incomplete gamma function to calcu- view. The opening angle, θ, was set to θ = 40◦, as suggested late the corresponding ∆χ2 corresponding to a 68.3% confi- dence interval. For our ten free parameters, this corresponds to by previous studies (e.g., Barthel 1989; Drouartetal. 2012). In practice, all inclinations i < 50◦present a line of sight free of ∆χ2 = 11.5. We defined the confidence interval as χ2min +∆χ2 (seeFigs.D.1−D.10andTable4). dustandthereforesimilarSEDs.Tosavecomputingtimebutin- cludethepossibilityofatype1AGN,weonlyusedi = 40◦ as Theuncertaintiesonstellarmassesandluminositieswerees- timatedwitha differenttechniquebecausethey aredirectlyex- atype1AGNtemplate.Table2presentsalltheassumedvalues pressedfromthenormalisationofthedifferentcomponents.We for the fixed and fitted parameters. The objective of this selec- considered the extreme values of each mass within the allowed tion is to explore typical torus configurations and evaluate how 68%confidenceintervalofallotherparametersanddefinedthis ourmethodhandlestheadditionofthisthird(AGN)component. asthelowerandupperuncertainties.FortheAGNluminositywe addedtheuncertaintyfromtheAGNfraction(f20µm)inquadra- 3.2. HybridAGN-SBcomponent tureandthenormalisationofthehybridcompoAnGenNt.Finally,we add a special caution of interpreting these uncertainties. They Wecreatedhybridtemplates,correspondingtothesumoftheSB onlyrepresentthestatisticaluncertaintiesontheparametersin- templatesproducedbyPÉGASE.3andanAGNtorus.Foreach duced from the noise in the observations and do not take into evolvingSBtemplate,wenormalisedanAGNtorus(foragiven accountpotentialsystematicsassociatedwiththemodelchoices inclinationi,sizeY,opacityτ )toafractionofthetotalflux 9.7µm and assumptions (see Appendix A for a more extensive discus- at20µmfollowingtheequation: sionofuncertaintiesandparametercalculation). F (t)= F (t)+F (i,Y,τ )× f20µm, (1) hybrid SB AGN 9.7µm AGN where F isthesumoftheAGNandtheSBtemplates, F hybrid SB 4. Results istheSBtemplate, F istheAGNtemplateand f20µm isthe AGN AGN WeperformedtheSEDfittingdescribedintheprevioussections relativefractionoftheAGNfluxcomparedtotheSBtemplateat 20µm2(seeTable2).Weperformedastep-wiseχ2-minimisation on the 11 sources in our sample, taking one (evolved compo- nentonly),two(evolvedandstarburst),andthree(evolved,star- thatpermittingdeterminingofintermediatesolutionsandcheck- ingfordegeneracies.Intotal,>107templateswerefittedforeach burst, and AGN) components into account. Table 3 shows the best χ2 depending on the number of components. The addition source. of extra components clearly dramatically increases the quality of the fit. This improvement is particularly strong from one to 3.3. Uncertaintiesofparameters,68%confidenceintervals two components and is still significant from two to three com- ponents(ameliorationbyafactorsuperiortotwo).Wetherefore Calculating uncertainties on the parameters can be challenging only consider the three-component fit in the remainder of this for a multi-dimensional analysis. However, as we performed paper. Table 4 summarises the results of the fitting, including 2 Weevaluated f20µmwithanidealisedfiltercentredon20µmof4µm theirassociated(68%)uncertainties.Figure2showsanexample bandwidthtominiAmGNiseeffectsfromsharpgradientsand/oremissionand of the best-fit model for 4C 41.17, while the results for the re- absorptionlines(i.e.PAHs). mainingthesamplearepresentedinFigs.D.1−D.10.Wepresent A109,page5of26 A&A593,A109(2016) 4C 41.17, z=3.792 Fig. 2. Best fit for 4C 41.17 (z = 3.792). The orange, blue and green dashed lines symbolise the evolved, starburst, and AGN components, respectively.Thesumofthecomponentsisthedarkline,theblackdiamondscorrespondtobroad-bandphotometry,theverticallinestothe1σ uncertainties and horizontal lines to the FWHM of the filters. The downward triangles represent the 3σ upper limits. The insets show the χ2 distributionforsevenofthefreeparametersofthefit.Fromtoplefttobottomright:ageoftheevolvedcomponent(orangelineintheSED),age ofthestarburst(bluelineintheSED),inclinationoftheAGNtorus(greenlineintheSED),fractionoftheAGNat20µm,equatorialopacityof thetorus,sizeofthetorusandinitialmetallicityofthestarburst.Black,blueandredindicatethethreedifferentapproachesweusedtoestimate the polarisation effect, without, with the lowest and the highest contamination, respectively. We note that in this case the lowest polarisation contribution is zero, hence the black and blue lines are superposed. The horizontal coloured lines correspond to the 68% confidence interval describedinSect.3.3. here the overall trends on the different components for the en- mass distribution seen in large field galaxy samples (e.g., tire sample. Overall, the fittings are satisfying for all sources Marchesinietal. 2009; Ilbertetal. 2013). Interestingly, these even if some discrepancies are still observed (see Sect. 5.1). It masses are also comparable to SMGs at similar redshift is important to note that three components are necessary to re- (e.g., Borysetal. 2003). This suggests that powerful ra- producetheSEDofthepowerfulradiogalaxiesinthissample, dio galaxies are drawn from the same sample of massive and each appears to be dominating at a particular wavelength galaxies, formed at high-redshift (z > 6) on a relatively range. short timescale. Similarly as in RV13, the evolved compo- nent clearly dominates the Spitzer photometry, consistent with the shift of the 1µm peak of the most evolved early-type 4.1. Evolvedcomponent galaxies. We call the oldest stellar component the evolved component TherightpanelofFig.3confirmsthisresultbecauseallred (>0.5Gyr), which can be defined as the radio galaxy host. points are clustered along the black lines (plain and dashed), Figure 3 shows the mass of the evolved stellar population as consistentwithanearlyformationepoch.Wealsonotethatthe a function of redshift (red points). The host galaxy appears fourHubbletypesusedinouranalysisprovidesimilarχ2values massive at z = 4, and this evolved massive component is within the 68% confidence ranges. We therefore conclude that present in every source throughout our redshift range (M > low-resolution broad-band photometry alone cannot constrain stel 1011M , except for MRC 0324-228 with 1010.7M ). Previ- theexactstarformationhistoryofthegalaxies,butstillfavoursa (cid:12) (cid:12) ous studies have shown that the most powerful radio galax- formationatearlytimesintheUniverse.Weemphasisethatthe ies present exceptionally high masses (e.g. DeBreucketal. ages determined by our method are well constrained, as illus- 2002;Rocca-Volmerangeetal.2004;Seymouretal.2007).This tratedinFigs.D.1−D.10wherethebestχ2 solutionsarelocated places powerful radio galaxies at the top of the galaxy indeepminima. A109,page6of26 G.Drouartetal.:DisentanglingSFandAGNinluminousHzRGs Fig.3. Twostellarpopulations,youngandevolved,shownasbluecirclesandreddiamonds,respectively.Theverticallinesindicatethe68% confidenceintervalsfromthefit.Whenpolarisationdataareavailable,wereportthefitwithlowestandhighestcontributionfromscatteredAGN light as the downward and upward empty triangles respectively. For clarity, we do not add the 68% confidence intervals for the polarisation- subtractedfits.Left:massesofthetwostellarpopulationsversusredshift.Right:agesofthetwostellarpopulationswheretheblacklineshows theoldestuniverseageasafunctionofredshiftassumingthestandardcosmologicalconcordancemodel(Sect.1).Thedashedlinerepresentsthe ageofagalaxyforz =10.Wenotethatthe68%confidenceintervalcanbesmallerthanpointsizeandpolarisation-subtractedfitsaremainly form containedinthe68%confidenceintervalswhenscatteredAGNlightisignored. 4.2. Intenseandmassivestarburst andlocalisedminima.Moderatelyhigh-resolutionspectroscopic observationsarenecessarytofurtherexplorethispartofthepa- Figure 3 and Table 4 show two remarkable features of the rameter space because the emission lines are the result of the youngcomponentrequiredtoreproducethestrongHerscheland ISMphotoionisationfromboththeyoungstarsandtheAGN. submmemissioninoursample.Firstly,thestarburstisextremely We make a final note about the star formation rate related massive,M >1010M .Ineachcase,itrepresentsasignificant stel (cid:12) tothestarburstcomponent.Bydefinition(inthemodel),allthe fractionofthemassofthesystem(5%to50%).Thesehigh-mass star formation is taking place in a single time step (we discuss fractions indicate that even if the bulk of the stars has already thedurationofSBinSect.5).Whenthefitconvergestoanage been formed at high-redshift (Sect. 4.1), the radio galaxies in >1Myr, the related star formation is therefore null. Moreover, our sample still experience vigorous star formation3. Secondly, theabsenceofinfraredimagingatarcsecresolutionandtheun- Fig. 3 also shows that this starburst spans a wide range in age, certainties about the duration of the burst make deriving a star 5Myr < tSB < 100Myr, indicating that star formation is on- formationratefromourfittingdifficult.Wethereforeusedsimple goingorintheaftermathofaviolentevent. physicalassumptionsandempiricalrelationstoestimatestarfor- It is also interesting to note that all but one (MRC 2104- mationratesfromthestarburstcomponent(seeSect.5.2). 242)ofthegalaxieshaveacolumndensityfactor K = 10.This indicatesthatstarburstsassociatedwithradiogalaxieshaveanon average ten times higher column density than the Milky Way4. 4.3. BolometricallyluminousAGNwithanopaquetorus This indicates that high-redshift radio galaxies, or at least the most IR-luminous radio galaxies, favour higher covering dust OurSEDfittingallowsustoderivesomeoftheAGNandtorus fractions,indicatingahighobscuration,similartoSMGsinthe properties of HzRGs, such as the size, opacity and inclination. sameredshiftrange(e.g.,Chenetal.2015). WithEq.(1)wecalculatedthecorrespondingintrinsicluminos- ity of the AGN (Table 4, last column). All our AGN are lumi- About half of the solutions favour a high initial metallic- ity (see Table 4). This could be expected for a starburst form- nous,withLAinGtN >∼1012 L(cid:12),mostwithLAinGtN ≈1013−14 L(cid:12).Un- fortunately,thesmallspanofpropertiesofoursample(interms ing from IGM gas, which is expected to be metal-rich at high- ofluminosityormass),usuallyonlycoveringoneorderofmag- redshift(e.g.,Pierietal.2014).Interestingly,theotherhalftends nitude,doesnotallowustoidentifyanycorrelationbetweenthe tohavealow(ornull)initialmetallicity,indicatingamorepris- AGN and host properties. Interestingly, these values are sim- tine gas source. The evolution of the starburst implies ISM en- ilar to optically luminous quasars observed at similar redshift richmentwhenSNexplodeandthereforetheobservedmetallic- (e.g.,Kollmeieretal.2006).Thisisconsistentwithradiogalax- ity will be higher than the one reported here (initial metallicity ies experiencing a strong radiation-dominated accretion phase ofthegas).Giventhatweusedonlybroad-bandphotometry,the withcopiousaccretionontothecentralsupermassiveblackhole. derivedinitialmetallicitiesshouldbetreatedwithagrainofcau- InFigs.2andD.1−D.10,weshowthattherelativecontribution tion, even though they do provide some clues about the overall at20µmandthereforetheAGNluminosityiswellconstrained. metallicity of the initial gas of the starburst (null, sub-solar, or super-solar).ThebroadtrendiswellshowninFigs.D.1−D.10, Except for 3C368, the tori in our sample appear opaque, wherethebest-fitsolutionsareshownnottobelocatedindeep with τ9.7µm ≥ 1.0. This opacity corresponds to AV > 18 using theDraine(2003)conversion.Theobscurationisalsorelatedto 3 We note that our sample selection requires bright IR emission (≥4 theinclination.Theinclinationrepresentstheamountofduston Herschel detections), implying the sample is likely biased towards thelineofsight.Mostofourbestfitsconvergetosolutionswith strongstar-formingsystems. i>40◦,indicatingType2AGN(edge-onview),asexpectedfor 4 N ∼6.8×1021atomscm−2isthestandardgalacticcolumndensity radiogalaxies.Theoneexceptionis3C368,thesourcewiththe HI inPÉGASE.3. lowest redshift in our sample, with i = 40◦. The best solutions A109,page7of26 A&A593,A109(2016) TN3aCbalme36e48.SummHartuyySbpaobeflethelofi1gt1Mt1i.n03a[g+−sEM00srve..10(cid:12)so]ullvt4se.5d[0MAc0ogy+−me52r00]p00o0ne0n.0Zt1∗12 lo1g1L1.60to[+−tL01..(cid:12)23] d light [%]10800 M4RCC T 40N 14J.021607-0274-4131M6RUCS S0 3302C9 443-362-822842USS 0828+193 3MCR4C700324 ES2a 1110..77+−+000...113 31500000+−+012000000 00..00029624 1121..08+−+000...370 arise 60 PKS1138 S0 12.3−+10..21 2000−+450000 0.0142 12.7−+00..02 pol −0.0 −200 −1.0 n MMRRCC02410064 EE22 1111..10+−+000...143 11220000++−14030000 00..00225599 1111..87+−+000...330 n of u 40 USS0828 E2 11.50+.00.0 140−02+000 0.0257 12.20+.00.3 tio 4C2858 S0 11.60+.00.1 2000+−00 0.0142 12.0−+00..02 rac 20 USS0943 E 11.3−+00..01 500+−120000 0.0114 12.1−+10..12 F θ=45° θ=90° −0.4 −100 −0.2 4TCNJ41210707 SE0 1111..59+−+000...211 1700000+−+312000000 00..00015478 1122..34+−+000...262 00 5 10 15 20 25 30 −0.1 −100 −0.7 Observed Polarisation [%] SBcomponent∗ Fig. 4. Highest and lowest polarisation diagram at 1500 Å. The grey shaded area represents the allowed fraction of unpolarised light de- Name Mass Age Zinit K Ltot pendingontheinclination(seetext).Polarisationmeasurementsofour log10[M(cid:12)] [Myr] log10[L(cid:12)] sourcesarereportedasverticallines withouttheuncertaintiesforvis- 3C368 10.1+0.1 7+1 0.01 10 12.7+0.3 ibility(i.e.theuncertaintiesincreasetheallowedrangeofunpolarised −0.2 −2 −0.4 lightaccordingly). 3C470 10.0+0.4 5+3 0.01 10 12.9+0.1 −0.2 −1 −0.5 MRC0324 10.8+−00..32 18+−74 0.01 10 12.9+−00..15 vary from Y = 10 to Y = 150 for the size. As we calculated PKS1138 10.7+0.3 20+15 0.0005 10 13.2+0.1 theintrinsicluminosity,wecanestimatethephysicalsizeofthe −0.3 −6 −0.4 MRC0406 11.3+0.4 45+45 0.0 10 13.1+0.1 torus (see Eq. (1) of Fritzetal. 2006). Taking the median of −0.2 −15 −0.4 MUSRSC02812084 1100..63+−+000...873 906+−+163650 00..001 110 1123..92+−+000...711 oanudrsRammapxle∼, L4Ai1nGtpNc,∼∼22.456×pc1,01a3nLd(cid:12)∼t6ra1n6splactefsorinYto=Rm1in0,∼604,.1anpdc 4C2858 10.8−+00..13 30−+130 0.005 10 13.3−+00..60 150, respectively. We note the difficulties in constraining the −0.5 −12 −0.1 toruspropertiesareexpectedbecauseevenwithdedicatedsam- USS0943 11.4+0.5 90+110 0.01 10 13.0+0.3 −0.8 −60 −0.1 ples in the nearby Universe, some parameters are only loosely 4C4117 11.4+−00..02 30+−55 0.001 10 13.4+−00..35 constrained (e.g., AsensioRamos&RamosAlmeida 2009). In- TNJ2007 11.0+0.4 45+95 0.001 10 13.0+0.2 terestingly, even after chosing the most relevant parameters for −1.0 −27 −0.1 the AGN torus, the properties of the torus vary from source to AGNcomponent source, with no clear trend on the entire sample (see Figs. 2 Name Y τ i Ltot andD.1−D.10). 9.7µm [degrees] log [L ] 10 (cid:12) 4.4. ModerateeffectofAGN-scatteredlightinbroad-band 3C368 150+0 0.1+0.0 40+0 11.7+0.3 SEDfitting −0 −0.0 −0 −0.4 3C470 60+90 10.0+0.0 50+10 13.2+0.5 −0 −0.0 −0 −0.4 The AGN emission is partly absorbed by the dusty torus or es- MRC0324 60−+900 1.0+−00..09 90+−010 12.3+−00..42 capes along the open part of the torus. Escaping photons can PKS1138 10+0 1.0+0.0 50+10 13.5+0.5 bescatteredbydustpresentintheionisationconetowardtheob- −0 −0.0 −0 −0.3 MRC0406 10+0 10.0+0.0 50+20 13.3+0.5 server,andtherebycontaminatetheoverallUVandopticallight. MRC2104 150−+00 10.0−+00..00 80−+010 13.7−+00..49 Thiscontributionespeciallyaffectstype2AGN(seenedge-on), USS0828 10+−090 1.0+−00.0.0 70−+3100 13.5−+00..15 andcanbestrongintheUVandopticalrest-frame(upto100%– 4C2858 60−+00 10.0−+00.0.0 90−+00 14.3−+00..67 diSeregoAlighierietal. 1996; Cimattietal. 1998; Vernetetal. −0 −0.0 −0 −0.1 2001). It is important to quantify this effect in our SED fitting. USS0943 150+0 10.0+0.0 70+10 13.4+0.6 −140 −0.0 −20 −0.0 WeusedthetechniquefromVernetetal.(2001)toestimatethe 4C4117 10−+1040 10.0+−00..00 80+−1100 13.9+−00..35 range of possible contributions of the optical AGN scattered TNJ2007 10+0 1.0+0.0 50+20 13.0+0.8 light, as illustrated in Fig. 4. We stress that this range relies −0 −0.0 −0 −0.1 on several assumptions regarding the geometry and inclination Notes. Theχ2 isgiveninthelastcolumnofTable3.Thetableisdi- oftheobscuringdustdistribution,whichareconsistentwiththe vided in three parts corresponding to each fitted component. For the adoptedtorusmodels(θ=40◦,inclinationi). evolvedcomponent:Z referstotheaveragemetallicitylockinstarsfor ∗ One fundamental property assumed here is the grey be- thegivenage.FortheSBcomponent:Z referstotheinitialmetallicity init haviour of polarisation, that is, the fact that polarisation from ofthegas,K isthefactorrelatingtothecolumndensity.FortheAGN dust scattering is approximately independent of wavelength component:Yreferstothesizeratiobetweentheouterandinnerradius (Vernetetal. 2001). This allows us to estimate the range of ofthetorus,τ theopacityat9.7µm,anditotheinclinationwith 9.7µm respecttotheobserver.Asymmetricuncertaintiesofeachparameterare AGN contribution in the UV and optical domain from a sin- providedasexplainedinSect.3.3.Wenotethatzerovaluesoftheun- gle polarisation measurement at 1500Å. When a polarisation certaintiesareduetothelackofcoverageintheparameterspaceand measurementwasavailable,wethereforeconsideredthefollow- arecontainedinthe68%uncertainties,refertoFigs.D.1−D.10. ingthreecases:(i)withoutpolarisation;(ii)withlowestpossible A109,page8of26 G.Drouartetal.:DisentanglingSFandAGNinluminousHzRGs AGN contamination; and (iii) with the highest possible AGN contamination. We applied the following correction to the UV USS 0828+193, z=2.572 andopticalphotometrybeforethefitting:wescaledthe1500Å contributionfromFig.4andsubtractedthecorrespondingfrac- tion from our broad-band photometry5 with a typical UV and opticalcontinuumquasarspectrumfromCristiani&Vio(1990). The latter represents a typical UV and optical SED of a type 1 AGN (i.e. without obscuration). We recall that these are the boundary cases and that the true contamination is located be- tweenthesetwoextremes. The most polarised source of our sample, USS 0828+193, with P = 10.0±2.0% is predicted to have a contribution from unpolarised light of between 0% and ∼65%. The remaining AGNcontributionisthusbetween∼35and100%(illustratedin Fig.5).Subtractingthelowest(35%)AGNcontributionhasal- mostnoeffectonthefit(χ2 = 35.7orχ2 = 37.1,respectively). This suggests either that (i) the addition of this component is notparticularlyjustified,orthat(ii)thetruecontaminationisnot USS 0828+193, z=2.572 close to the lowest 35% contribution of the AGN. As the high- estpossibleAGNcontributionforthissourceis100%intheUV, thiswouldleavenoroomforanystellarcontribution.Suchanex- ceptionalcase(onlyoneinoursample)cannotbeproperlytaken intoaccountinourfittingprocedure.Wethereforeconsiderthe lowest possible AGN contribution as an upper limit on the ef- fect on USS 0828+193. For this lowest contribution, the effect onthefittingislimitedtoafactorof∼2forthemassand∼2for the age of the stellar components (see Fig. 3). For the remain- ing six sources with polarisation measurements, the difference inducedbypolarisationistypicallycontainedinour68%confi- denceintervals(seeemptysymbolsinFig.3).Theeffectonthe otherparametersispresentedintheinsetsintheSEDsinFigs.2 and D.1−D.10 (blue and red points for the lowest and highest Fig.5.SEDsofUSS0828+193showingtheimpactofthepolarisation contamination,respectively). onourfitting.ThecolourcodingisthesamethanFig.2.Thetoppanelis We conclude that the effect of polarisation for our broad- withoutpolarisation,whilethebottompanelistheminimalsubtraction bandSEDfittingisweakerthanourglobaluncertaintiesoneach case.Themaximalcaseofa100%subtractionisnotreported(seetext). parameter. Polarisation therefore has a negligible effect on our main conclusion, and we ignore this in the remainder of this paper. themorphologyofoursourcesmakinguseofacombinationof multi-wavelengthdata. 5. Discussion 5.1. Modellinguncertainties,limitations, Our fitting procedure enables us to carefully disentangle the andmulti-wavelengthdataset threemainspectralcomponentsofoursources:anevolvedstel- lar component, an SB component and an AGN. We focus on The self-consistency in our models is a key point of our mod- the new information provided by the Herschel data in the far- elling:forinstance,thecomputedFIRemissionisdependenton infrared(FIR),allowingustodisentangletheSBfromtheAGN thedustpropertiesdefinedbytheevolutionaryscenarioandpro- contribution. We find that the two stellar components are mas- duced from the stellar population evolution (ISM enrichment, sive(>1010M )andthatboththeAGNandtheSBareverylu- (cid:12) see Sect. 3.1.1). We discuss here the limitations and degenera- minous (>1010L ). The properties of the evolved stellar pop- (cid:12) cieswithinourlibrariesandfromourobservationsonthefitting ulation confirm a very high formation redshift of the galaxy procedure. We refer to Sect. 3 and Appendix A for the discus- host, leading to local massive early-type galaxies as previously sion of the model assumptions and systematics introduced by proposedinliterature(e.g.Lilly&Longair1984).Interestingly, ourchoiceofcodeandlibraries. our sources have similar properties as quasars and SMGs in ThePÉGASE.3modelonlyexploresalimitednumberofstar thesameredshiftrange(e.g.,Alexanderetal.2005;Wangetal. formationhistories.Fortheevolvedcomponent,thestarforma- 2011; Leipskietal. 2014, see Sect. 4). We discuss the limita- tionhistoryissetbytheHubbletype,althoughthedatapresented tionsofourmodellinginafirststepandthenexploretheprop- hereareunabletopreferasingletypewithintheconfidencein- erties of star formation making use of physical arguments and tervals (Sect. 3.3). However, the effect on the age and mass of empiricallawsfromSEDfitting.Finally,weclassifyanddiscuss stellar populations is limited, within a factor of 2 (see RV13 for better details). For the starburst, the star formation is con- 5 Weareinterestedonlyinthecontinuumemission,thereforewere- sidered instantaneous (in one time step). Forming ∼1010M in moved the strong emission lines from the templates by linearly inter- (cid:12) polatingthecontinuumbelowtheline.Inaddition,linepolarisationis suchashorttimeappearshighlyunlikely.Nevertheless,wenote typicallyclosetonullbecausethephysicalprocessesinvolvedinline thatonlyashortburstcanreproducethestrongsubmmemission emissionaredifferentfromthecontinuum. in our sample. We stress that star formation is not necessarily A109,page9of26 A&A593,A109(2016) spatially concentrated, but is most probably temporally con- V = σ/0.6(e.g.,Rixetal.1997),wecalculatethedynamical circ nected. Therefore, the SB component may be interpreted as a timeast =2r /σ . dyn SB SB collectionofsmallerclumps,formedinashortperiod.Inaddi- The projected star formation per surface unit is a power- tion,RV13havediscussedtheeffectofthisδfunction,showing ful proxy to estimate the physical conditions of the starburst thattheagefromthefittingismeasuredwithinafactorof∼2for (e.g., Lehnert&Heckman 1996). Based on the physical sizes constantstarformationoveralongerperiodoftime. and timescales derived above, we calculate the star formation TheFritzetal.(2006)AGNmodelisoneofseveralmodels density, Σ , with the following formula Σ = SFR/(2πr2 ) SFR SFR SB availableintheliterature.Directcomparisonsbetweendifferent withSFR=M /t . SB dyn modelsaredifficultbecauseofthevariationinthenumberoffree Itisinterestingtonotethattheoverallpropertiesfromthese parametersandbecauseoftheassumptionsonthetorusproper- empirical calculations indicate (i) a high efficiency of forming ties itself (smooth, clumpy, bi-phased). More investigations are starswithansSFR>20Gyr−1inallcases,consistentwithmea- reserved for further papers, taking different AGN models into surements from Drouartetal. (2014); (ii) a dense star-forming account(e.g.Podigachoskietal.,inprep.).Wethereforedonot region Σ > 5M yr−1 kpc−2, similar to the local starburst SFR (cid:12) furtherdiscusstheimplicationsofthesmoothtorusmodelonour (e.g.Bigieletal.2008)andlargelyexceedingthelimitforstar- analysis, but we warn that the results provided here are depen- burst driven winds (e.g. Lehnert&Heckman 1996); and (iii) dentonthischoice.Wenotethatpartofourbestfithasdifficul- half (5/11 sources) of our sample presents an SB mass of the ties to reproduce the 30µm restframe data (e.g. PKS 1138-262 same order as the host component (Table 5), strongly suggest- andUSS0828+193).Themostlikelyexplanationisthelackof ing merging activity. These estimates have important implica- exploration in the parameter space, especially on the size pa- tions in terms of star formation properties because it suggests rameter.Alargertorus,havingdustatlargerradii,wouldbetter that the star formation is very efficient in very confined re- reproducethiscoldemission. gions. However, the large projected radius (>1kpc), suggests Finally,observationsthemselvescanbeaffectedbysystem- thatthestarformationarisesinseverallocationswithinthesame atics. This is particularly true for our set of multi-wavelength galaxy. observationsthatwereobtainedwitharangeofinstrumentsand telescopes, each presenting their own calibration uncertainties, 5.3. Localisationofthestarformationfromtheimages especiallyatmoderatesignal-to-noiseratio.TheFIRandsubmm domain are particularly challenging observations and the pho- ThelimitedspatialresolutionofHerschelmakesitdifficulttoan- tometry estimation method can have a significant effect on the swerthekeyquestionsofthelocationofstarformationandofthe final measurement, especially in the <5σ range (by a factor of process that leads to the observed star formation. Complemen- ∼20%; Popessoetal. 2012). A good illustration is the 160µm tarytothespectralinformation,spatialinformationcanbeused photometryin4C28.58(Fig.D.8),wherethe160µmfluxmight to understand the localisation of star formation and the overall beexpectedtobehigher.Similarly,sky-subtractionstrategy(on- evolutionarystatusoftheradiogalaxiesinoursample.Wemade offversusraster)mayaffectthe850µmphotometry,possiblyex- use of the high-resolution HST and radio imaging, along with plainingthedifferencesinthecolddustpartoftheSED.These the moderate resolution of the NIR data, to investigate the rel- uncertainties are mitigated by the good sampling of our SEDs, ative location of the emission at different wavelengths6. We vi- however. suallyclassifiedoursourcesintothreenon-exclusivecategories as follows (and report their likelihood in Table 5 as none, low, medium,andhigh): 5.2. StarformationpropertiesfromtheSEDfitting Intheprevioussectionsweestimatedseveralcrucialproperties 1. “exsitustarformation”,wherestarformationoccursoutside relatedtothestarformationinoursample,suchastheage,aver- thegalaxy,at>5kpcoftheradiogalaxy(late-stagemergers agecolumndensity,massand,tosomeextent,initialmetallicity. andstar-formingcompanionsfallintothiscategory); Evenifthestarformationmostlikelyoccurssimultaneouslyin 2. “in situ star formation”, which refers to star formation that different clumps, we treat this component as a single idealised occurswithintheradiogalaxy,at<5kpc(rotatingstructure starburst (completely localised and isolated). This provided us andnuclearstarformationfallintothiscategory); withinformationusefulforcomparingoursampletoothergalax- 3. “AGN-driven star formation” corresponding to star forma- ies.Assumingthestarbursthasaprojectedcirculargeometry(at tion linked to AGN activity, independently of spatial scale constantdensity),weestimatethesizeofthestarburstregionas (bothjet-inducedandshock-inducedstarformationfallinto follows: thiscategory). (cid:32) (cid:33)1 rSB ∝ πMKNSB 2 , (2) TtichaellaimnditNoIfR57kdpocmisaidnrsivaenndbtyhethteypreicsaollustiizoenoofbataminaesdsivinetghaelaoxpy- HI atz>1(Fanetal.2010). where r is the projected radius, M is the mass of the star- SB SB burst,KisthecolumndensityfactorandN thecolumndensity HI (here N = 6.8×1021 atomscm−2). These values are reported 5.3.1. Exsitustarformation HI in the third column of Table 5. With an estimated size, we can Ex situ star formation can be divided into two subcategories, a defineatypicalvelocitydispersionassumingequilibrium: massive star-forming companion or a multiple system of less (cid:114) GM massive sub-units, which can be interpreted as a signature of σ ∼ SB, (3) SB r SB 6 The HST data are downloaded from the HLA. The Herschel data where σ is the velocity dispersion of the gas, and G is the SB are unresolved when compared to optical or NIR images and are not gravitational constant. We recall that when the system is col- reportedforsimplicity(butavailablein Drouartetal.2014). lapsing,thisvalueisanupperlimit,whileforanexpandingsys- 7 Werecallthatatz > 1,1(cid:48)(cid:48) correspondsroughlyto8kpc,assuming tem, this value becomes a lower limit. Assuming the relation thestandardcosmologicalmodel. A109,page10of26
Description: