The origin of large emission line widths in massive galaxies at redshifts

We present a sample of 22 massive galaxies with stellar masses > 10 10 𝑀 ⊙ at 3 < 𝑧 < 4 with deep H and K-band high resolution spectra (R=3500-3000) from Keck/MOSFIRE and VLT/KMOS near-infrared spectrographs. We find a large fraction have strong [OIII]5007 and H 𝛽 emission lines with large line widths ( 𝜎 100 – 450 km/s). We measure the sizes of our galaxies from Hubble Space Telescope images and consider the potential kinematic scaling relations of our sample; and rule out an explanation for these broad lines in terms of galaxy-wide kinematics. Based on consideration of the [OIII]5007 / H 𝛽 flux ratios, their location in the Mass–Excitation diagram, and the derived bolometric luminosities, we conclude that Active Galactic Nuclei (AGN) and their Narrow Line Regions most likely give rise to this emission. At redshifts 3 < 𝑧 < 4, we find significantly high AGN fractions in massive galaxies, ranging from 60–70% for the mass range 10 < log ( 𝑀 ★ / 𝑀 ⊙ ) < 11, with a lower limit 30% for all galaxies within that redshift range when we apply our most stringent AGN criteria. We also find a considerably lower AGN fraction in massive quiescent galaxies, ranging from 20-30%. These fractions of AGN point to the period between 3 < 𝑧 < 4 being a time of heightened activity for the development of supermassive black holes in the massive end of the galaxy population and provide evidence for their role in the emergence of the first massive quenched galaxies at this epoch.


INTRODUCTION
Massive galaxies in the early universe appear to assemble extremely rapidly.Spectroscopic detection of massive quiescent galaxies (QGs) of log( ★ / ⊙ ) > 11 at 3 <  < 4 suggests that the stellar content was already in place in the first 1.5 Gyr of the Universe's history (Marsan et al. 2017;Glazebrook et al. 2017;Tanaka et al. 2019;Di Valentino et al. 2021;Forrest et al. 2020;Esdaile et al. 2021;Forrest et al. 2022).The low-redshift stellar mass function requires that these compact galaxies both grow in size (about 3-5 times their effective radius and size-scaling relation) and quench quickly in order to evolve into the early-type galaxies observed in the local universe.Moreover, the number density of massive quiescent galaxies at  > 3 exceeds the predictions from galaxy formation simulations by a factor of up to ten (Schreiber et al. 2018, hereafter S18;Donnari et al. 2021;Lustig et al. 2022).
The physical process that causes the rapid quenching and suppression of star formation in massive QGs is unclear.Both external and internal mechanisms can influence quenching.External phenomena such as ram-pressure stripping, tidal forces, and major or minor mergers can remove or heat the gas in group and cluster environments, ★ E-mail: mmartinezmarin@swin.edu.auaffecting the gas reservoirs of galaxies.Internal processes like active galactic nuclei (AGN) can generate strong winds that can remove gas from the galaxy preventing star formation for a short period or permanently (Moy, E. & Rocca-Volmerange, B. 2002;Feruglio et al. 2010;Fischer, J. et al. 2010;Komugi 2023;Vayner et al. 2024) or may also produce positive feedback, as AGN outflows can trigger star formation by compressing cold, dense gas (Magnelli, B. et al. 2009;Silk & Norman 2009;Zinn et al. 2013;Carniani et al. 2018;Shin et al. 2019).
The efficiency of internal and external quenching mechanisms correlates with the galaxy stellar mass.Low-mass galaxies are more likely to be quenched by environmental effects, while in high-mass galaxies quenching is dominated by internal processes (Peng et al. 2010;Baldry et al. 2006;Somerville & Davé 2015).The star formation activity of AGN hosts can change with redshift.At low redshifts, AGN are predominantly present in galaxies with suppressed star formation activity (Ho 2005;Salim et al. 2007), while this trend weakens at high redshifts, as AGN predominantly reside in star-forming galaxies (Cowley et al. 2016).
AGNs emit radiation across the entire electromagnetic spectrum, as the different emission processes depend on the geometry of the AGN and the line of sight (Antonucci 1993; Urry & Padovani 1995).Therefore they can be detected through many observational tech-niques, such as broad and narrow optical emission lines (Baldwin et al. 1981a;Kewley et al. 2001;Kauffmann et al. 2003), X-ray emission (Brandt & Alexander 2015), mid-infrared emission (Lacy et al. 2004;Stern et al. 2005;Jarrett et al. 2011;Salyk et al. 2019), and radio emission (Heckman & Best 2014;Padovani et al. 2017).
The x-ray and radio properties of typical QGs at  > 3 are barely explored, mainly by individual detected sources that do not trace the entire AGN population and are biased towards high luminosity galaxies.Nevertheless, enhanced AGN activity is found in the QG population with high mass at 0 <  < 5 (Ito et al. 2022).Marsan et al. (2017), from a sample of 6 galaxies, found 5 massive galaxies hosting AGNs at 3 <  < 4, simultaneously fitting the galaxy and possible AGN contribution to the far-ultraviolet to the far-infrared spectral energy distribution and strong [OIII] emission lines.
Recent observations using the NIRSpec spectrograph on JWST (Moy, E. & Rocca-Volmerange, B. 2002) have revealed a massive quiescent galaxy, GS-9209, at  = 4.658 with an AGN identified by broad H emission, implying that AGN feedback plays a significant role in quenching star formation (Carnall et al. 2023;Nanayakkara et al. 2023).The Advanced Deep Extragalactic Survey (JADES) has found 42 type-2 AGN candidates in low mass galaxies up to  = 10, analyzing UV emission lines with NIRSpec.(Scholtz et al. 2023).Nevertheless, due to the limited field of view of JWST, there have been few spectroscopic observations of massive galaxies at these redshifts to date.S18 presented a sample of massive quiescent galaxies observed with the MOSFIRE spectrograph on Keck, ten of which were spectroscopically confirmed.Two of these galaxies had broad emission lines in their spectra, for which they proposed an AGN origin.In this paper we consider a larger sample of 22 massive galaxies at 3 <  < 4 with log ( ★ / ⊙ ) > 10 and similar quality spectroscopy.This extends the analysis to include the massive star-forming population at these redshifts.We will show that broad emission lines are surprisingly common, and we explore how these can arise physically.
This paper is structured as follows.Section 2 describes the ZFOURGE and multi-wavelength data sets, the sample construction, and emission line fitting.In Section 3, we calculate structural properties and consider a scenario where the broad lines arise from galaxy-wide kinematics.In section 4 we present our findings regarding the most plausible explanation for the broad emission lines arising from AGN and discuss the implications in Section 5. Finally, in Section 6, we summarise the implications of the results.

Spectral Sample
Our sample comprises massive galaxies originally identified using photometric redshifts from the medium band near-infrared FourStar Galaxy Evolution Survey (ZFOURGE; Straatman et al. 2016).This survey utilizes the near-infrared FOURSTAR instrument on the Magellan telescope across the Chandra Deep Field South (CDFS; Giacconi et al. 2002), COSMOS (Scoville et al. 2007), and UDS (Grogin et al. 2011;Koekemoer et al. 2011) legacy fields.FOURSTAR has medium-bandwidth filters from 1-1.8 µm to calculate robust photometric redshifts, the imaging reaches 5 point-source limiting depths of ∼ 25 AB mag in Ks (Straatman et al. 2016;Spitler et al. 2014).ZFOURGE is supplemented with existing data from CAN-DELS Hubble Space Telescope (HST) WFC3/F160W imaging and Spitzer/Infrared Array Camera (IRAC), as well as other ground-based imaging, generating multiwavelength catalogs spanning 0.3-8 µm.
We selected galaxies from ZFOURGE with log( ★ / ⊙ ) > 10, and photometric redshift 3 <  < 4, and  < 24.5.We only used the galaxies with a quality flag use=1, obtaining a sample of 130 massive galaxies within the 11'x11' ( total of 450 sq.arcmin) coverage of the CDFS, COSMOS and UDS fields.We search the KOA and ESO archives gathering spectra for 54 of the 130 massive galaxies at 3 <  < 4 from MOSFIRE and KMOS spectrographs.
The spectral data was taken from 4 different programs (PI: Glazebrook, Oesch, Illingworth, Kewley), including all the spectra from S18, where the observations, integration time, and data reduction are described.The observations were made with the MOS-FIRE (McLean et al. 2012) spectrograph installed on the Keck I telescope at the summit of Mauna Kea in Hawaii, and the KMOS spectrograph installed at the VLT at Cerro Paranal in the Atacama Desert in Chile (Davies et al. 2013).MOSFIRE is a multi-object infrared spectrograph with a field of view of 6' × 3' that can be used to observe up to 46 slits per mask, with a resolving power of R ∼ 3500 in a single band from Y to K bands.KMOS is also a K-band Multi-Object Spectrograph capable of performing spectroscopy in the near-infrared band of up to 24 targets simultaneously, with a spectral resolution of R∼ 3000 − 4000.We note our sample includes all galaxy spectra from S18.
All galaxies were observed in the H and K-band on several masks across multiple nights, all observed with standard ABBA exposures, nodding along the slit with a 0.7 arcsec slit for mask configurations.Each observing program has different integration times in the K band, described in detail in S18.MOSFIRE and KMOS have a wavelength coverage of 1.93-2.45µm and 2.038-2.290µm in the K band.Therefore, galaxies at z=3.8 reach the detection limit for [OIII]5007 emission lines.The noise also increases at longer wavelengths, and there is a loss due to atmospheric transmission at 2.35 µm.The minimal detection limit for the [OIII]5007 emission line is about 4.7×10 −18 erg/s/c 2 .
Out of the 54 galaxies with spectra, we have secured spectroscopic redshifts for 33 galaxies.In Section 2.2, we fit a sample of 22 objects with strong [OIII] and H emission lines, which we then analyze further.Of the remaining 11, 3 lie outside 3 <  < 4, and we discard these.A further 5 have redshifts from S18 with weak or no emission lines, and another one has redshifts from weak emission lines or absorption lines, which we determined using slinefit (Schreiber 2017).One further galaxy does not display the [OIII]5007 emission line necessary for our analysis, showing only H emission line.Lastly, the slit of one of the galaxies is contaminated with a second object and it has been rejected.21 galaxies have no spectroscopic redshift that can be determined from the spectra.

Spectra and Emission line fluxes
We measure the flux density of the [OIII]5007 and H emission lines with the non-linear least-squares minimization and curve fitting python library LMfit-py (Newville et al. 2016).LMfit allows us to fit all three lines ([OIII]5007, [OIII]4959 and H) together, constraining the center and velocity dispersions consistently with redshift, utilizing a Gaussian model centered on the emission line's wavelength.A visual inspection does not find any cases where H emission lines have a different width than [OIII] emission lines.We used the Specutils python package to calculate the standard equiva-lent width (EW) measurement.We employ a continuum-normalized spectrum to obtain the absorption line width of the adjacent continuum that has the same area as that of the absorption line.
We fit the emission lines with LMFit-py and remove galaxies from the sample with [OIII]5007 emission line signal-to-noise (S/N) <3.We remove from the sample galaxies with a standard deviation larger than 30% in the full width at half maximum (FWHM) measurements, reducing the sample to 26 galaxies.We focus only on galaxies with detected [OIII]5007 and H emission lines.Therefore we remove four other galaxies (COSMOS-4417, COSMOS-15056, COSMOS-17779, and COSMOS-12000) with H4862 absorption lines or which lack detectable H emission lines, further reducing the sample to 22 galaxies.
Four galaxies in the MOSFIRE archives were covered by more than one spectrum (COSMOS-4214, COSMOS-12173, COSMOS-20169, and UDS-6023) with different slit positions and orientations, which can lead to different velocity profile measurements.Therefore we modeled the spectra jointly, so the H and [OIII] doublet could be identified.
We show an example 2D and 1D spectrum with the best fit for all the object in each field -COSMOS, UDS, and CDFS -in Apendix Figures A1-3, with the best fit for the [OIII]4959, [OIII]5007 and H4862 lines obtained with LMfit.

Spectral classification of the sample
The ( −,  −) rest frame color-color diagram (hereafter UVJ) has become a standard technique to characterize high-redshift galaxies (Franx et al. 2003;van Dokkum et al. 2003;Williams et al. 2009) and in particular to classify galaxies into quiescent, low obscuration star-forming and high-obscuration star-forming categories based on the overall shape of the spectral energy distribution (SED).Use of this technique requires extensive multiwavelength photometry and accurate redshifts, which we can find in ZFOURGE.
For our sample, we find galaxies with strong emission lines (see Table 4); the contribution of these lines to the photometry in the K band, can result in incorrect classifications with UVJ which assumes the broad band photometry traces the stellar continuum.To correct for this we follow the procedure from S18 to correct the photometry for emission line contamination.Equation 1 below, assuming a constant continuum flux density within the filter and a constant filter response over the spectral extent of the lines: Where   is the original flux density, () is the response curve of the corresponding filter,   obs the observer-frame equivalent width, and   the central wavelength of the line .
In practice, the emission lines only impact the K band, and after correcting the photometry, all our galaxies dropped in flux, with a median offset of 32 % and CDFS-9833 being the most affected, with 74% of its flux removed.In Figure 1 we show the UVJ diagram of the final sample after the emission lines from the broadband photometry were subtracted.It can be seen that our galaxies mostly fall in the star-forming region of the diagram (both obscured and unobscured); we can also observe the two galaxies that overlap with S18.
In Figure 1, we show the distribution of the specific star formation rate (sSFR) of our final galaxy sample and compare it to the massive quiescent galaxies from S18. Star-forming galaxies show an empirical relation between specific star formation rates and stellar mass, widely known as star-forming main sequence (MS) (Elbaz et al. 2007;Noeske et al. 2007).S18 uses a definition for 'quiescence' of sSFR=sSFR  = 0.15Gyr −1 , valid for galaxies  > 3, and correspond to a 'relative' quiescence, as galaxies below this threshold at  > 3 are forming stars with rates an order of magnitude lower than the bulk of the star-forming population.We can observe that most of the galaxies in our sample reside close to the MS in comparison with the MS from ZFOURGE and above the sSFR  except for 3 galaxies.One of them is UDS-8197 from S18, that analysis proposed as an AGN candidate from the emission lines in its spectra.One of the other two galaxies UDS-8092 was detected as an AGN in X-rays by Cowley et al. (2016) and the final one CDFS-6691, even though it appears quiescent in the UVJ plot and lower than sSFR  , shows [OIII]5007 emission lines, with low H emission.We will keep displaying these galaxies in red to be able to trace them during the development of the paper.

Structural Properties
We measure the half-light radius along the semi-major axis   and Sersic index  of our galaxy sample from H-band WFC3/F160W HST images, utilizing the profile modeling software GALFIT (Peng et al. 2002).We use single components Sersic profiles, providing an initial   , Sersic index, and axis ratio  = /.Point Spread Function (PSFs) were taken from (Skelton et al. 2014;Grogin et al. 2011;Koekemoer et al. 2011) for COSMOS, GOODS, and UDS fields.
We compare our measurements obtained with GALFIT to the galaxies with a good fit flag from CANDELS/3D-HST catalogs (van der Wel et al. 2014), finding a good agreement (2-3%).We find four galaxies with bad fit flags in our sample (COSMOS-20133, UDS-6023, CDFS-7535 and CDFS-9332), where inputs for axis ratio (b/a) and position angle (PA) from (van der Wel et al. 2014) were inaccurate, creating poor models and therefore unphysical   and Sersic index for those galaxies.We corrected those galaxies that resulted in a bad fit and summarized the results in Table 4.We also noticed no sign of multiple sources or merging systems in the GALFIT analysis.Purple crosses are massive quiescent galaxies at 3 <  < 4, from S18. Grey dots are the first selection of galaxies with spectra from MOSFIRE and KMOS archives within ZFOURGE.The diamonds are the emission line galaxies with flux correction (see section 2.3).Red, orange, and blue correspond to quiescent, dusty, and non-dusty star-forming galaxies correspondingly.We have two galaxies in common with S18, identified as AGN candidates.Right: Logarithm of the specific star formation rate (sSFR) versus stellar mass.The black line indicates the main sequence (MS) for galaxies in the ZFOURGE calculated by S18.The entire ZFOURGE galaxy sample is shown as dark contours.The dashed grey line indicates the quenching limit (sSFR  ) defined by S18.

Stellar Masses
The CANDELS/3D-HST catalogs also provide stellar masses for our galaxy sample.These masses were calculated with photometric redshifts, and do not take for account the presence of AGN in the galaxies, which could originate the emission lines that we observe in the spectra.Therefore, we re-calculated them with the new spectroscopic redshifts (see Table 4 and ZFOURGE photometry).
We use the code CIGALE (Boquien et al. 2019) to re-compute the stellar masses, utilizing the (Bruzual & Charlot 2003) stellar population model.We assume a (Chabrier 2003) initial mass function (IMF) with a solar metallicity ( = 0.02).For the dust attenuation, we adopt dustatt_modified_starburst in CIGALE (Calzetti et al. 2000) with   up to 6 mag, and for the dust emission models we use (Casey 2012).For star formation history (SFH), we use a delayed SFH with optional exponential burst (except for galaxy UDS-3417, which uses a delayed SFH with an optional constant burst/quench in order to find a suitable fit with a χ2 < 2), we allow the e-folding time and stellar age varying from 0.3-1 Gyr to 0.5--2 Gyr, respectively.For the AGN component, we adopt the skirtor2016 module based on a clumpy torus model SKITOR AGN from (Stalevski et al. 2012(Stalevski et al. , 2016)).The relative strength between the AGN and galaxy components are set by the fracAGN parameter, we allow fracAGN to vary from 0.1 to 0.9, and the 9.7m optical depth to be 3, 7 and 11.The stellar masses calculated with CIGALE show a median offset of 0.242 and a DEX scatter of 0.14 compared to the CANDELS/3D-HST catalogue.(See Table 4).

POSSIBLE ORIGINS OF BROAD LINES
After correction for emission line contributions, most of our galaxies lie in the star-forming regions, both obscured and unobscured, of the UVJ diagram.Broad emission lines are common with 13/22 having velocity dispersions > 200 km −1 (with several > 400 km −1 and rest-frame equivalent widths up to 2700 Å; see Table 4).Such a high frequency of strong, broad emission lines are not typically seen in lower-redshift samples of massive galaxies (Pagotto, Ilaria et al. 2021;Freeman et al. 2019).
Broad lines in galaxies can arise due to two principle mechanisms.The first is the kinematics of ionized gas distributed throughout a galaxy and moving in its gravitational potential.For example, the integrated line width of emission lines was used to study the Tully-Fisher scaling relations in high-redshift galaxies (Glazebrook 2013;Christensen & Hjorth 2017) before the advent of integral field spectroscopy.In this dynamical scenario, we would expect to be able to estimate an implied dynamical mass from the galaxy size and velocity width, and this should show a reasonable correlation with stellar mass (e.g.Straatman et al. 2022;Kriek et al. 2009).There is a basic underlying assumption here of virial equilibrium conditions; either pressure or rotationally supported.Given the sizes of our galaxies and that the spectroscopy is with ground natural seeing it is not possible to determine the spatial extent of the emission line flux; so we use the HST continuum sizes measured in subsection 2.4 instead.We note that our line widths are much larger than typically seen in lowredshift kinematic surveys (Franx et al. 2008).Nevertheless, some of our galaxies are very compact, and that can give rise to comparable large dispersions in the emission lines (van Dokkum et al. 2009;Nelson et al. 2014).
The second common mechanism we consider is that the emission comes from a small region at the galaxy center due to the accretion of matter by a supermassive black hole, i.e. an AGN.Our line widths are comparable with those seen in Type 2 AGN1 .We can test this scenario by analyzing the ionization mechanism using line ratios.This is commonly done using the 'Baldwin, Phillips & Terlevich' (BPT) diagram (Baldwin et al. 1981b), which can classify galaxies into star-forming or AGN, from their emission-line intensity ratios.This is extensively used to diagnose ionization mechanisms in galaxies; in particular, AGN are strongly enhanced in [OIII]/H and mildly enhanced in [NII]/H which serves to distinguish strong AGN from star-forming galaxies (Kauffmann et al. 2003;Heckman et al. 2004;Yan et al. 2006).In our case the redshifted [NII]/H lies outside the red-end of our ground-based spectra so instead we utilize the massexcitation (MEx) diagram (Juneau et al. 2011).The MEx diagram leverages the relationship between stellar mass and metallicity in place of the [NII]/H ratio.By measuring the H and [OIII] ratios, we can identify AGNs through the MEx-diagram.We note that we may still find a weak correlation between a dynamical mass (calculated under the galaxy-wide assumption) and stellar mass, as the velocity dispersion from the [OIII] line will trace the size of the narrow line region (NLR), which extends up to a few parsecs, as its radius is proportional to the square root of the [OIII] luminosity (Bennert et al. 2002;Schmitt et al. 2003) and hence be correlated with the mass of the SMBH which also correlates with stellar mass (Kormendy & Ho 2013).

Exploring a dynamical origin
In a dynamical scenario, the emission lines in our galaxy spectra would trace the kinematics of ionized gas distributed throughout the galaxy.In this case, we can expect a strong correlation between the estimated dynamical masses and the stellar masses (Schreiber et al. 2018;Esdaile et al. 2021).Thus, star formation would be enough to explain the emission lines in our galaxies.
We use the   and the F160W band   measurements to calculate 'dynamical masses' with the following equation (2): Where () is the virial coefficient used to approximately account for structural and orbital non-homology (Cappellari et al. 2006): (3)   is the half-light radius,   is the line of sight velocity dispersion (LOSVD) within the   and G is the gravitational constant.This equation essentially assumes virial equilibrium and has been commonly used in high-redshift kinematic papers (Cappellari et al. 2013;Schreiber et al. 2018;Esdaile et al. 2021).
The derived dynamical masses are in the range 9 < log 10 (/ ⊙ ) < 11.6.We note that some of these values are rather high, and in the light of our conclusions below favoring the AGN scenario they are not likely to represent true dynamical masses.
After obtaining spectral classification, stellar mass, effective radius, velocity dispersion, and dynamical mass for our galaxies, we analyze our sample and compare it with those of massive galaxies at various redshifts.
We examine the size-velocity-stellar mass scaling relations; and then compare the implied dynamical mass derived from size and velocity to stellar mass.Our principle comparison is with other samples of massive galaxies at redshifts 3 <  < 4 with kinematic analysis that have a good correlation between stellar and dynamical mass.In particular in Figure 2 we compare our galaxy sample with galaxies at  ∼ 0 from SDSS (Maraston 2005;Thomas et al. 2005), galaxies from Mendel et al. (2020) 2014), who observed that star-forming galaxies are 3.2±1.3times larger than quiescent galaxies (r  =0.63±0.18kpc) at z∼4.Notably, we also detected a wide range in size, and poor correlations of size vs mass and size vs velocity, with some galaxies exhibiting similar effective radii as those at lower redshifts and others as large as those at z∼0.This suggests that the diversity of the galaxy population was established early in the universe.Overall our galaxy sample extends to lower stellar masses than the QGs at the same redshift range, we interpret this as a selection effect as absorption line kinematic studies can only be done for the brightest and most massive objects.
Comparing the velocity dispersion and the stellar mass of the galaxies, we find a weak correlation with a Pearson correlation factor of 0.4 ± 0.01 after 1000 Monte Carlo sampling.We note that the comparison samples measured velocity dispersions using absorption lines, whereas we used emission lines.However, these should be suitable for testing a galaxy-wide kinematic origin of such lines.In low redshift star-forming galaxies, emission and absorption largely trace the same kinematics (Cortese et al. 2016).
In the dynamical origin hyopthesis, the velocity dispersion of the emission lines detected will trace the kinematics of ionized gas produced by star formation distributed throughout a galaxy due to its gravitational potential.In this scenario, we should expect an agreement between the stellar and dynamical mass of the galaxy.Comparing the dynamical mass with the stellar mass (see Fig 2-d), we found a weak Pearson correlation factor of 0.29 ± 0.01 after 1000 Monte Carlo sampling.
This weak correlation could arise from the similar weak correlation of velocity dispersion and stellar mass which could also arise in the AGN scenario (see Section 4.1).To understand this better, we need to see if the combination of size and velocity dispersion which goes into the dynamical mass calculation, provides more information than from velocity alone.To test this, we ran 1000 Monte Carlo simulations where we replaced the size of each galaxy with a size value picked randomly from the full sample and recomputed the dynamical mass.Doing this, we found a correlation coefficient of 0.32 ± 0.02.This is not significantly different from the values using the true size, so we conclude that the weak correlation we see is simply due to the velocity-stellar mass correlation, which can also arise under the AGN scenario.
Putting this all together we conclude that the velocity dispersion of the [OIII] 5007 emission line does not appear to be tracing the gravitational potential of the galaxy, thus is unlikely to be due (solely) to star formation.2022), at 3 <  < 4. The blue and pink triangles are star-forming and quiescent galaxies from (Straatman et al. 2014) respectively.Green diamonds are star-forming galaxies from our sample at 3 <  < 4.

THE PRESENCE OF AGN IN OUR SAMPLE
AGNs emit radiation across the entire electromagnetic spectrum, so in order to further test the AGN hypothesis, we first look for counterpart detections of our galaxies at other wavelengths.We cross-match our galaxy sample with the Cowley et al. (2016) AGN catalog (using a tolerance of ±1 arcsec), which considers infrared, X-ray, and radio AGN detections in the ZFOURGE sample.Only one of the galaxies in our sample is detected as an infrared AGN in this catalog (COSMOS-12173), with a counterpart detected in X-rays.We find one galaxy detected in both X-ray and radio (CDFS-13954) and three detected only in X-ray (CDFS-9332, CDFS-8428 and UDS-8092).Two QGs were earlier identified as AGN candidates by S18 (COSMOS-20133, UDS-8197), due to their significant [OIII] emission.In summary, we only have counterpart observations that indicate the presence of AGNs for 5 of our 22 galaxies.

Emission line diagnostic diagrams
By analyzing certain emission line intensity ratios, it is possible to determine whether the ionized regions belong to a galaxy that is actively forming stars or if the lines are being powered by a more intense ionizing spectrum linked to an AGN (Kewley et al. 2001;Kauffmann et al. 2003;Brinchmann et al. 2004 The MEx plot replaces the ratio of the redder emission lines ([NII]6583/H) with the stellar mass to fill gaps where these lines fall in especially noisy regions or outside of the atmospheric cutoff.This diagram utilizes a probabilistic approach to split galaxies into sub-categories and finds intrinsically weak and heavily absorbed AGNs, taking advantage of the known stellar mass-metallicity (MZ) relation.One might expect the MEx diagram to be more sensitive to changes in the stellar-mass-metallicity (MZ) relation compared to the traditional BPT.Juneau et al. (2011) presented the MEx as an AGN diagnostic tool applicable to galaxies out to  = 1.However they explored possible evolutionary effects, and they found some evidence that the intermediate-redshift galaxies may be offset in the MEx diagram relative to the lower-redshift sample (Juneau et al. 2014, J14 hereafter;Coil et al. 2015;Strom et al. 2017).These differences can be attributed to variations in the underlying stellar populations, as well as differences in the interstellar medium conditions (ISM).J14 shifted the original criterion by Δlog(M ★ /M ⊙ )=0.25 expanding the capabilities of the MEx diagram up to  = 2.
In Figure 3, we show the log([OIII]/H) ratios versus the stellar mass of our galaxy sample.We plot with a black line the J14 limits.Based on this limit, all galaxies in our sample are AGN candidates, with only one galaxy (UDS-8844) likely to be in the MEx intermediate region, indicating both star formation and an AGN.
The boundary between AGN and star-forming galaxies in emission line diagnostics diagrams is fuzzy; (Kewley et al. 2013) (hereafter K13) used a theoretical approach to make predictions about the locations of galaxies and AGNs in the BPT diagram up to  = 3.They considered four different evolutionary scenarios, which involved changing the ISM conditions (ranging from normal to extreme) and the metallicity (from metal-poor to metal-rich) of the AGN narrow line region at high redshift.An "extreme" ISM condition could result from a larger ionization parameter, a denser ISM, or a harder ionizing radiation field.In a metal-rich scenario, the gas near the AGN is enriched to a higher metallicity than is typically found in the host galaxy, while in the metal-poor scenario, the gas near the AGN has the same metallicity as the gas in the host galaxy.Regardless of the scenario, all predictions for galaxies at  = 3 indicate that sources with log 10 ( [OIII]/H) > 1 are ionized by AGN, but there is a lack of observational data to test these predictions at this redshift.(Coil et al. 2015) (hereafter C15) tested the MEx diagram with a sample of galaxies identified as AGNs either through X-ray or IR emission at  ∼ 2.3 from the MOSFIRE Deep Evolution Field survey (MOSDEF) sample (Kriek et al. 2015).They propose a substantially higher shift of Δlog(M ★ /M ⊙ )=0.75 for the J14 limits.We have included a dashed line on the graph to represent the stricter classification criteria used in C15.According to this approach, three galaxies (CDFS-9833, UDS-14075, CDFS-7752) would be considered star-forming, and the other three galaxies (UDS-8844, COSMOS-4214, UDS-6639) would fall into the MEx mixing sequence.Galaxies along the mixing sequence are called composites or transition objects.They may have a mix of star formation, shock excitation, and/or AGN activity.However, we also plot with grey dots the AGN sample from C15, and two galaxies with a log(M ★ /M ⊙ ) < 10.5, confirmed as AGNs with X-rays, which are under the C15 limits (Fig. 3).We also show with a dot-dashed line the theoretical limit from K13, noticing that the five X-ray detected galaxies in our sample (plotted with a star symbol) are above this limit, located in the AGN region of the MEx diagram.The line ratios diagnostic methods can be sensitive to metallicity and can move the position of the AGN to lower metallicities (lower masses in the MEx diagram; Kewley et al. 2013).Although low-metallicity AGNs are extremely rare among local galaxies, we are dealing with a sample of galaxies at high redshift, which can be an important factor in explaining the presence of our galaxies in this mixing sequence (Harikane et al. 2023a).
This galaxy sample could be a missed population of heavily obscured AGN, or the X-ray emission is too faint to be detected at this redshift (Straatman et al. 2014;Cowley et al. 2016;Forrest et al. 2020).

[OIII] luminosities and bolometric correction
We can calculate the bolometric luminosity of our galaxy sample to compare them with galaxies that are known to host an AGN at lower redshifts and see if they have similar luminosity ranges.To obtain the bolometric luminosity of our galaxy sample, we used the [OIII]5007 luminosity and applied the luminosity-dependent [OIII] bolometric correction factor,  OIII from Lamastra et al. (2009),   ∼ 454 ×  OIII , assuming that the observed [OIII] line luminosities are a lower limit to the intrinsic luminosity of [OIII].
From Table 4 we can see that the luminosities range from 5 × 10 44 erg s −1 to 2 × 10 46 erg s −1 , which is consistent with them hosting powerful hidden AGNs.In Figure 3 we compare the bolometric luminosity ( bol ) of our galaxy sample with the  bol from the catalog of optical spectral properties for all X-ray selected SPIDERS AGNs at redshifts 0 <  < 2.5 (Coffey, D. et al. 2019) and star forming galaxies from SDSS-III DR9 spectra between 0 <  < 0.75 (Brinchmann et al. 2004;Kauffmann et al. 2003;Tremonti et al. 2004).The SPIDERS AGN and the SDSS-III DR9 catalogues, estimated the bolometric luminosities using the bolometric corrections from (Richards et al. 2006) and (Shen et al. 2011) for  3000Å and  5100Å .We utilize the [  ] 5007 luminosities from the value-added catalogs from the MPA-JHU (Brinchmann et al. 2004;Kauffmann et al. 2003;Tremonti et al. 2004) to calculate the bolometric luminosities of the star-forming galaxies of the SDSS-III DR9 catalogs.From the plot, it is evident that the luminosity distribution of our galaxy sample is greater than that of the star-forming galaxy region.Moreover, our sample lies at the upper end of  bol of the AGN region.As shown in other works, this can provide further evidence that our sample galaxies likely host AGN (Marsan et al. 2017;Harrison et al. 2016).
The radiative transfer code SKITOR (Stalevski et al. 2016), implemented in CIGALE, can model the distribution of dust in the torus as a two-phase medium.This allows us to estimate the infrared radiation emitted by the dust to obtain the ratio of the torus to the AGN luminosity, as well as the fraction of AGN IR luminosity to total IR luminosity (      ).We estimated the      from the SED fitting of our galaxy sample (see Table 4) and obtained a      > 0.2 for the entire sample.A      < 0.2 is considered a non-dominant AGN, while a      > 0.85 is a highly-dominant AGN (Ciesla et al. 2015;Ramos-Padilla et al. 2021).Overall, the galaxies in our sample show an important      , supporting the presence of AGN in our galaxies.

AGN Fraction
The fraction of AGNs (  AGN ) within a galaxy population is an indicator of AGN activity and supermassive black hole growth at galaxy /H from K13.The stars represent the galaxies detected either in X-rays or radio.On the right side, we plot the velocity dispersion against the   ∼ 454 ×     of the galaxies.The orange contours are star-forming galaxies from the SDSS-III DR9 at 0 <  < 0.75 and the purple contours are AGN galaxies from the SDSS IV SPIDERS at 0 <  < 2.5.
centers.The energy output of AGNs is a powerful feedback mechanism that can effectively regulate gas inflow and star formation in large galaxies.As such, assessing the level of AGN activity in massive galaxies with respect to redshift is a valuable means of evaluating the impact of AGN feedback at any cosmic epoch.The  AGN for massive galaxies with log(M ★ /M ⊙ )>11 seems to evolve across cosmic time, as different surveys show that  AGN changes with redshift, showing a dramatic transition in the  AGN happening at  ∼ 2.5 when the universe was ∼3Gyr old (Kauffmann et al. 2003;Kriek et al. 2009;Marsan et al. 2017).
We calculate the fraction of AGNs in our sample assuming all of our galaxies are AGNs, utilazing the J14 limits of the MEx diagram as follows: Where the numerator is the number of AGN detected by broad [OIII] and high [OIII]/H within 3 <  < 4, and the denominator, the number of galaxies with spectra with spectroscopic redshift detected within 3 <  < 4 (See Figure 4).We use the beta distribution (Gupta 2011) to derive errors on the ratios.We note that our sample may be biased towards the strongest emission lines as these are easier to detect and that we have 21 galaxies with unconfirmed spectroscopic redshifts between 3 <  < 4. We summarise the  AGN in Table 2.
In order to gain a better understanding of the  AGN , it is important to compare our results to similar stellar mass ranges and also see the evolution of the  AGN with cosmic time.So, we compare our  AGN in Figure 5 to different samples, in three stellar mass ranges (log(M ★ /M ⊙ )>10, 10<log(M ★ /M ⊙ )<11 and log(M ★ /M ⊙ )>11), and to different  AGN between 0 <  < 4. Our sample comes originally from ZFOURGE, therefore we also calculate the  AGN for galaxies with stellar mass 10<log(M ★ /M ⊙ )<11 and log(M ★ /M ⊙ )>11 of the ZFOURGE cat-alog at different redshift ranges to trace the cosmic evolution of the  AGN up to  = 4. Nevertheless, as the ZFOURGE catalog provides AGN from multi-wavelength detections, this time the numerator is the number of Xray+FIR+Radio AGN detections, and the denominator is the number of objects in the ZFOURGE catalog with flag use=1.We summarize the AGN fraction in table 2.
Our sample is primarily made up of MOSFIRE spectra, which allows us to compare our AGN fraction to that of the MOSDEF survey.Aird et al. (2017) utilize a multiwavelength approach of optical, IR, and Xray data to find the AGNs in the MOSDEF sample.They identified 55 AGNs out of their 482 galaxy sample, where 24 AGNs were detected by the line rations of the BPT diagram obtained from MOSFIRE spectra.The MOSDEF sample covers galaxies from 1.4 <  < 3.8, nevertheless, most of the optical AGNs are at 2 <  < 3. So, we just compare a rough  AGN = 0.134 ± 0.009 for galaxies between 1.4 <  < 3.8 and log(M ★ /M ⊙ )>10 (in blue colors) in Figure 5.
Finally, for galaxies with log(M ★ /M ⊙ )>11, we compare (in purple colours) our  AGN = 0.714 ± 0.   We have considered the evident from line widths, luminosities and [  ]/  ratio lying in the MEX-AGN region as pointing to all the galaxies having an origin of line emission from AGN. However not all galaxies meet all these criteria simultaneously.If we conservatively consider much stricter criteria considering AGNs with line ratios of [OIII]/H  > 10 and a large velocity dispersion of  5007 > 300 km/s, we find an  AGN = 0.29 ± 0.09 We can calculate a robust lower limit on  AGN if we make the extreme assumption that all galaxies with photometric redshifts between 3 <  < 4 (see Fig 4) actually have spectroscopic redshifts in this range but lack emission lines and do not harbor AGN.This gives us an  AGN = 0.361 ± 0.088 and  AGN = 0.14 ± 0.05 with the stricter criteria.
We can also calculate the f   for quiescent galaxies in our sample and consider quiescent galaxies from S18.Of the 22 galaxies in our sample, 3 were classified as quiescent galaxies in the UVJ plot (See Fig 1b).From S18 we have 12 quiescent galaxies, 2 of them in common with our sample, and one confirmed (UDS-8197) as a quiescent galaxy.Therefore we have 13 quiescent galaxies spectroscopically confirmed, where 3 of them host an AGN.This gives us an f   = 0.23 ± 0.11 for quiescent galaxies in ZFOURGE for galaxies at 3 <  < 4.This is significantly less than the f   for star-forming galaxies.We notice that the three AGNs in quiescent galaxies are within the stellar mass range of 10 10 < ⊙ <10 11 , so the AGN fraction for these stellar mass range is  AGN = 0.3 ± 0.13.

Shocks
Finally, we briefly consider shocks as an explanation for the line ratios.In 4.1, we applied the MEx line diagnosis method to determine the excitation sources of galaxies using ratios of strong optical emission lines of [OIII] and H.AGNs lie at large [OIII]/H ratios on this diagram, where we found 64% of our galaxy sample according to C15 criteria, or 95% according to J14 (See fig 3. Nevertheless, the dependence on metallicity of the line ratios diagnostic methods can drag the position of the AGN to lower masses in the MEx diagram Kewley et al. (2013).Separating between star formation, AGNs, and shocks is always difficult.Integral Field Spectroscopy is the most successful technique for separating these components among the various diagnostic diagrams proposed (Allen et al. 1998;Best et al. 2000;Moy, E. & Rocca-Volmerange, B. 2002;Groves et al. 2004;Feltre et al. 2016;Jaskot & Ravindranath 2016).Nevertheless, we only measure [OIII] and H emission lines, so we cannot apply these diagnostics.Nevertheless, we calculate the bolometric luminosity of the galaxy sample (See fig 3 ) utilizing the [0III]5007 line, which provides a robust tracer of AGN power and is less contaminated by emission due to star formation activity.The high bolometric luminosity obtained for our sample makes it unlikely that the broad emission lines from our galaxies are powered only by shocks, as bolometric luminosities found for shocks in detailed studies with IFU data for the galaxy NGC 1068 are about L  =7-24 ×10 41 erg/s (D'Agostino et al. 2019).

DISCUSSION
The high-AGN fraction -70% of our star-forming galaxies and 30% of our quiescent galaxies -shows that the period spanning 3 <  < 4 seems to be a time of great activity for the growth of supermassive black holes for the massive galaxy population.The observed AGN fractions exhibit an upward trend regardless of the stellar mass range.However, this trend is particularly remarkable for galaxies with log(M ★ /M ⊙ )>11.These findings align with similar studies with smaller samples of massive galaxies, where a f   =0.8 (5/6) is found at 3 <  < 4 by the strength of the line rations of their spectra, in combination with X-ray and radio detections (Marsan et al. 2017), highlighting the evolution of  AGN in cosmic time.
The fraction of time that an AGN is active (duty cycle) in massive quiescent galaxies (log(M ★ /M ⊙ )>10) decreases at high stellar masses, while it increases with stellar mass for star-forming galaxies.At the same time, the fraction of galaxies with black holes accreting mass strongly evolves with redshift being comparable with starforming galaxies at  = 2 (Aird et al. 2017).This could explain the different AGN fractions found in our sample regarding star-forming and quiescent galaxies, as short duty cycles make it harder to find AGNs in action at higher redshifts.
Hydrodynamical simulations can also give us constraints on the AGN population, although their calibration, sub-grid models of black holes and galaxy formation vary.(Habouzit et al. 2021) compared the six Illustris, TNG100, TNG300, Horizon-AGN, EAGLE and SIMBA large-scale cosmological simulations (Nelson et al. 2017;Marinacci et al. 2018;Springel et al. 2017;Naiman et al. 2018;Pillepich et al. 2017;Kaviraj et al. 2017;Schaye et al. 2014;Davé et al. 2019).Finding a population of AGN that agrees with observational constraints at all redshifts and luminosities is challenging for simulations.Nevertheless, they found that the fraction of AGN is always higher at high redshifts for all the galaxy stellar mass bins.TNG and SIMBA simulations present similar trends to observations, although their AGN number densities peak at different redshifts.For simulated AGN with L  > 10 44 erg s −1 in massive galaxies of log(M ★ /M ⊙ )>11 between 3 <  < 4, TNG300 and SIMBA show an AGN fraction of about 80% (Habouzit et al. 2021).
Although the triggering mechanism for the AGN activity still seems stochastic, AGNs are known to alter galaxies' star formation and be a driver of quenching for galaxies at low redshifts.The high AGN fraction found suggests it might be the dominant quenching mechanism at high redshifts (Marchesini et al. 2009;Aird et al. 2017).
There are two primary processes through which AGN feedback is thought to occur.One is called the "quasar-mode feedback" (Silk & Rees 1998).In this process, wind ejects gases from galaxies within this mode, thereby impeding star formation.This process is widely believed to be prevalent in high-luminosity AGNs, such as those featuring QSOs, which are close to the Eddington limit.The other mode of feedback is known as "radio-mode feedback" or kinetic-mode feedback, as described by (Fabian 1996).In this mode, low-luminosity AGNs, which make up less than one percent of the Eddington luminosity, heat the circumgalactic and halo gas through their radio jets, thereby preventing the gas from cooling.Compared to the quasar-mode feedback, radio-mode feedback is expected to keep the quiescence rather than reduce the star formation.
Low luminosity AGN have been found in X-ray and radio studies supporting the critical role of AGNs in star formation and quenching of massive galaxies (Ito et al. 2022).The low range of  bol = 10 44 − 10 46 erg s −1 in our sample of galaxies, would support the kinetic-mode feedback scenario as a quenching mechanism.Cos-mological models and simulations require a feedback mechanism to replicate the properties of the massive galaxy population (Springel et al. 2005;Croton et al. 2006;Somerville et al. 2008;Choi et al. 2015;Weinberger et al. 2018;Davé et al. 2019), as it is challenging to achieve the necessary energy levels through stellar feedback alone (Bower et al. 2006).
AGN feedback has been incorporated in various ways (Somerville & Davé 2015), a common prescription in semianalytic models involves major mergers that trigger AGN-driven winds, expelling gas and eventually truncating star formation (Hopkins et al. 2008).Cosmological simulations such as IllustrisTNG suggest that AGN feedback is necessary as early as  ∼ 6 to produce massive quiescent galaxies at  = 4 (Hartley et al. 2023;Qin et al. 2017).
The first census of type-1 AGNs at  > 4 (Harikane et al. 2023a) identified by JWST/NIRSpec deep spectroscopy indicates that about 5 % of the galaxies at  = 4 − 7 harbor faint type-1 AGNs, while studies of local AGNs imply that only 1-2 % of galaxies with similar bolometric luminosities are type 1 AGNs.A high fraction of the broad-line AGNs indicates that the number density of such faint AGNs is higher than an extrapolation of the quasar luminosity function, implying a large population of AGNs, including type 1 and type 2, in the early universe.

CONCLUSIONS
We analyze a sample of 22 massive galaxies with masses 10 < log( ★ / ⊙ ) < 11 at redshift 3 <  < 4 from the ZFOURGE catalog.The galaxies' spectra show high velocity dispersions and large equivalent widths in their [OIII] and H emission lines.We conducted tests to identify the sources of broad emission lines in our galaxy sample.We evaluated two scenarios of dynamical and AGN origin: • We find that a dynamical scenario is improbable based on the weak correlation between implied dynamical and stellar masses, i.e. that we are not seeing star-forming gas distributed around a galaxywide potential.
• Five of our 22 galaxies were detected either in X-ray or radio in the ZFOURGE catalog, unequivocally indicating the presence of an AGN.
• We apply line diagnostic MEx diagrams to further distinguish between star-forming galaxies and AGNs.These show that in the majority of cases (19/22) the sources are likely to be AGN.This is primarily driven by the high [OIII]/H values.The other three objects lie in intermediate regions where the classification depends on the choice of MEx boundary but are still highly likely to be AGN.
• We derive bolometric luminosities from the [OIII] luminosities and these are in the range typical of AGN at lower redshifts.
• From the CIGALE SED fitting, we obtained a fraction of AGN IR luminosity to total IR luminosity higher than 0.2 for the entire sample, indicating AGN contribution in all the galaxies of the sample.
• The AGN fraction is high; overall we find fractions of 65.3±8.8% in the low mass bin increasing to 71.4 ±14.9% in the high mass bin.This is higher than previous work and higher than at lower redshifts.
• Applying a stricter criteria to the AGNs considering line ratios of [OIII]/H  > 10 and a large velocity dispersion of  5007 > 300 km/s, we find a lower limit of  AGN = 29 ± 9%.
• The lower limit for the AGN fraction is 36.1 ±8%, assuming that all galaxies with photometric redshift between 3 <  < 4 have spectroscopic redshifts in this range but lack emission lines and do not harbour AGN and  AGN = 14 ± 5% with stricter criteria.
Overall our results point to a ubiquity of AGN in massive galaxies at 3 <  < 4 and is consistent with a picture in which the massive quenched galaxies at this epoch would have been produced by AGN quenching.Future work could improve this.For our sample it is now possible to access H, [NII] and [SII] lines using the NIRSPEC spectrograph on the James Webb Space Telescope thus accessing the full BPT diagram.Multiple emission lines can be measured from rest-frame UV to near-infrared which will allow improved diagnostics of the ionizing mechanisms, the density and pressure of gas in the narrow line region, and its metallicity.We predict that the emission lines in our sample would arise from point sources at the center of each galaxy; this could be tested using the NIRSPEC Integral Field Spectrograph or with ground based IFS.The latter could also measure galaxy-wide kinematics and hence derive scaling relations between true dynamical mass and central black hole mass.

Figure 1 .
Figure1.Left: UVJ colour-colour plot.Purple crosses are massive quiescent galaxies at 3 <  < 4, from S18. Grey dots are the first selection of galaxies with spectra from MOSFIRE and KMOS archives within ZFOURGE.The diamonds are the emission line galaxies with flux correction (see section 2.3).Red, orange, and blue correspond to quiescent, dusty, and non-dusty star-forming galaxies correspondingly.We have two galaxies in common with S18, identified as AGN candidates.Right: Logarithm of the specific star formation rate (sSFR) versus stellar mass.The black line indicates the main sequence (MS) for galaxies in the ZFOURGE calculated by S18.The entire ZFOURGE galaxy sample is shown as dark contours.The dashed grey line indicates the quenching limit (sSFR  ) defined by S18.
at 1.4 <  < 2.1, and quiescent massive galaxies from Esdaile et al. (2021), Forrest et al. (2022) and Tanaka et al. (2019) at 3 <  < 4. Considering size-mass relation, our analysis indicates that the galaxy sample has an average effective radius of r  =2.9±0.15kpc, as shown in Fig 2-a and 2-b.Interestingly, this finding aligns with Straatman et al. (

Figure 3 .
Figure3.Light blue diamonds are non-dusty star-forming galaxies, orange diamonds are dusty star-forming galaxies and red diamonds are quiescent galaxies from our sample.On the left side, we show the Mass excitation (MEx) plot.The grey lines are the limits of the J14 MEx plot, the black dashed line is the C15 criteria for the MEx diagram.The dot-dashed line corresponds to the theoretical BPT ratio for [OIII]/H from K13.The stars represent the galaxies detected either in X-rays or radio.On the right side, we plot the velocity dispersion against the   ∼ 454 ×     of the galaxies.The orange contours are star-forming galaxies from the SDSS-III DR9 at 0 <  < 0.75 and the purple contours are AGN galaxies from the SDSS IV SPIDERS at 0 <  < 2.5.
149 to (Marsan et al. 2017) with  AGN = 0.83 ± 0.231 in the same redshift range 3 <  < 4. (Marsanet al. 2017) analyze a sample of 5 galaxies, where 4 of them were identified as AGNs from strong [O III] emission lines in the IR SED properties and X-ray and radio detections.We also compare our AGN fraction to the ZFOURGE catalog  AGN = 0.391 ± 0.0698.We also add the  AGN from SDSS(Kauffmann et al. 2003) at  = 0 and the  AGN from(Kriek et al. 2009) at 1.7 <  < 2.7, to trace the  AGN through redshift.Regarding the lower stellar mass range (10<log(M ★ /M ⊙ )<11), we can observe the  AGN remains similar through redshift, except for our recent AGN detections that boosted the  AGN = 0.0561 to  AGN = 0.653 ± 0.088 between 3 <  < 4. On the other hand, in

Figure 4 .
Figure 4. Summary Schema of the  AGN .We use a beta distribution to calculate the error bars.

Figure 5 .
Figure 5. AGN fraction evolution between 0 <  < 4. Blue, orange and purple colors, differentiate the different stellar mass ranges.The dots show the    for ZFOURGE Cowley et al. (2016), the diamond show the    for SDSS Kauffmann et al. (2003), the triangle show the    for Kriek et al. (2009), the hexagon show the    for MOSDEF, the cross show the    Marsan et al. (2017) and the stars show the    of this work.

Figure A1 .
Figure A1.Best fit of the H, [OIII]5007 and [OIII]4959 emission lines for each galaxy of COSMOS field.For each example, we show on the top left panel the 2D spectrum from keck/MOSFIRE, on the middle panel the spectrum with the best fit in red, on the bottom panel the residuals, and in the right corner a WFC3/F160W HST image stamp.

Figure A2 .
Figure A2.Best fit of the H, [OIII]5007 and [OIII]4959 emission lines for each galaxy of CDFS field.For each example, we show on the top left panel the 2D spectrum from keck/MOSFIRE, KMOS for CDFS-6691 and CDFS-7535 on the middle panel the spectrum with the best fit in red, and on the bottom panel the residuals.In the right corner, we show the WFC3/F160W HST image stamp.

Figure A3 .
Figure A3.Best fit of the H, [OIII]5007 and [OIII]4959 emission lines for each galaxy of UDS field.For each example, we show on the top left panel the 2D spectrum from keck/MOSFIRE, on the middle panel the spectrum with the best fit in red, and on the bottom panel the residuals.In the right corner, we show the WFC3/F160W HST image stamp.

Table 3 .
ZFOURGE and SDSS AGNs were selected by X-ray+FIR+Radio, (Marsan et al. 2017) AGNs were selected by X-ray, [OIII] and H emission lines, Kriek+07 and our work AGNs were selected by [OIII] and H emission lines.

Table 4 .
Final properties of the galaxy sample.