Star Formation and AGN Activity 500 Myr after the Big Bang: Insights from JWST

We consider the effect of including an active galactic nuclei (AGN) component when fitting spectral energy distributions of 109 spectroscopically confirmed z ≈ 3.5–12.5 galaxies with JWST. Remarkably, we find that the resulting cosmic star formation history is ≈0.4 dex lower at z ≳ 9.5 when an AGN component is included in the fitting. This alleviates previously reported excess star formation at z ≳ 9.5 compared to models based on typical baryon conversion efficiencies inside dark matter halos. We find that the individual stellar masses and star formation rates can be as much as ≈4 dex lower when fitting with an AGN component. These results highlight the importance of considering both stellar mass assembly and supermassive black hole growth when interpreting the light distributions of among the first galaxies to ever exist.


INTRODUCTION
Star formation and the growth of supermassive black holes (SMBHs) are key ingredients in galaxy formation as evidenced by their major contributions to the light distributions of galaxies (Kauffmann & Haehnelt 2000;Schaye et al. 2010).Star formation manifests at many different wavelengths, from the ultraviolet all the way to the radio (Conroy 2013;Davies et al. 2016).While SMBHs are unobservable, their contributions to spectral energy distributions (SEDs) originate from the accretion disk that surrounds the black hole, powering socalled active galactic nuclei (AGN) that radiate at similar wavelength regimes to the star formation (Richards et al. 2006;Shen et al. 2020).Analyses, and interpretations, of star formation in galaxies may therefore be contaminated by AGN because these processes manifest at similar wavelengths and with potentially similar contributions to the total SED of galaxies (Cardoso et al. 2017).
Thus, galaxies with a significant AGN component are often excluded from studies of star formation.Indeed, this was the approach taken by Driver et al. (2018) in calculating the cosmic star formation history (CSFH; Lilly et al. 1996;Madau & Dickinson 2014) where they used a combination of the mid-infrared and radio lumi-nosities to effectively filter out strong AGN from their sample, finding that only a ≈ 0.1 dex uncertainty is induced on the CSFH out to z ≈ 5 as a result of the AGN filtering.While it is encouraging that this selection does not precipitate a greater uncertainty, this approach fails to acknowledge the coevolution of star formation and the growth of SMBHs.
In fact, AGN are thought to inject energy into the host galaxy through the process of accretion-driven feedback.The net effect is to enhance the velocity dispersion of the surrounding gas and/or drive gas outflows (Magorrian et al. 1998;Ferrarese & Merritt 2000) that can affect star formation (Davies et al. 2019;Katsianis et al. 2019;Matthee & Schaye 2019;Davies et al. 2022).In light of these results, if we are to understand how galaxies form and evolve it is crucial to consider the union of these two fundamental processes.This is especially true at the z ≈ 3.5 − 12.5 frontier currently being probed by the James Webb Space Telescope (JWST ), among whose key objectives is to shed light on the formation of the first stars and black holes in the first galaxies (Gardner et al. 2006).Indeed, JWST has already provided valuable insights into the first galaxies and AGN (e.g., Adams et al. 2023;Yan et al. 2023;Curtis-Lake et al. 2023;Arrabal Haro et al. 2 2023;Naidu et al. 2022;Bouwens et al. 2023a;Bunker et al. 2023;Larson et al. 2023;Kocevski et al. 2023;Juodžbalis et al. 2023;Übler et al. 2023).
At this stage, JWST has arguably generated more questions than it has answered about the hitherto unexplored early Universe.Recent results suggest that star formation was more efficient at z ≳ 10 than predicted by models based on typical baryon conversion efficiencies inside dark matter halos (e.g., Harikane et al. 2022Harikane et al. , 2023b)).JWST has also yielded candidates of quenched galaxies at z ≈ 3 whose star formation histories suggest a burst of star formation at z ≈ 11 that appears to be at odds with our current grasp of stellar mass assembly (e.g., Carnall et al. 2023;Glazebrook et al. 2023).
In this work, we investigate these emergent puzzles by exploring the union of star formation and AGN activity at z ≳ 3. We use the SED fitting code ProSpect (Robotham et al. 2020) to determine the physical properties of z ≈ 3.5 − 12.5 galaxies.Critically, ProSpect is able to account for the contribution from both stars and the AGN when fitting the galaxy SED, and thus an invaluable tool to understand the interface of star formation and black hole growth.We build upon the analysis of D 'Silva et al. (2023b) and investigate the effect of AGN on the CSFH.
In Section 2 we describe our JWST data and highredshift galaxy sample.In Section 3 we describe our methods of SED fitting and calculating the CSFH.In Section 4 we present our main results.Our conclusions are summarised in Section 5. We use concordance ΛCDM cosmology throughout, with H 0 = 70 km s −1 Mpc −1 , Ω M = 0.3 and Ω Λ = 0.7.All of our results are derived using a Chabrier (2003) initial mass function.
We used galaxies with confirmed spectroscopic redshifts with the NIRSpec instrument on JWST that were presented in Arrabal Haro et al. (2023) and Nakajima et al. (2023).Our sample includes the z = 12.17 galaxy observed in the GLASS survey whose redshift was spectroscopically confirmed with the Atacama Large Millimetre Array (ALMA; Bakx et al. 2023).Our final sample comprises 109 galaxies with spectroscopic red-shifts from z ≈ 3 to z ≈ 12, with 108 galaxies from the CEERS survey and the z = 12.17 galaxy from GLASS.

Imaging and photometry
On the basis of our final spectroscopic sample, we made photometric catalogues suitable for SED fitting.Our photometric catalogues are derived from the combination of Hubble Space Telescope (HST ) data from legacy surveys and Near Infrared Camera (NIRCam; Rieke et al. 2005;Rieke et al. 2023) JWST imaging.
For the NIRCam data, all available UNCAL files from the two programs were downloaded from the Mikulski Archive for Space Telescopes (MAST) and processed through to the CAL stage using the JWST Calibration Pipeline 1 .
Our CEERS reductions used pipeline version 1.10.2(Bushouse et al. 2023a) and CRDS pmap=1084.Our GLASS reductions used pipeline version 1.11.2 (Bushouse et al. 2023b) and pmap=1100.We do not expect significant differences between these data as a result of these slightly different calibration versions, especially since JWST's NIRCam performance has significantly converged since its first light (e.g., Adams et al. 2023).As a test, we remade a subset of the CEERS data with the updated pipeline version and pmap, finding negligible differences in the imaging data for all NIR-Cam modules and filters between versions.The most significant update between pipeline versions is improved snowball and outlier detection.For reference, CRDS version 1084 was implemented on 2023 May 5 and CRDS version 1100 was implemented on 2023 July 14.
We performed a custom pipeline step that corrects for 1/f noise and wisp artefacts and subtracts the sky background in the images.These custom steps were built on the source extraction tool ProFound (Robotham et al. 2018).Further details on the 1/f removal can be found in Windhorst et al. (2023) and details on the wisp removal step can be found in Robotham et al. (2023a).
These artefact-free and sky-subtracted CAL files were then used to make mosaics via inverse-variance weighted stacking with the novel code ProPane2 (for further details, please see Robotham et al. 2023b).Hot pixels were patched with the equivalent pixels from a median stack.We draw attention to this stacking step as it is an entirely novel stacking code, unlike many of the custom JWST pipelines (e.g., GRIZLI, PENCIL) that rely on DrizzlePac3 .
HST observations were folded in to improve our coverage of the SEDs.For CEERS we downloaded the 3 HST mosaics provided by the CEERS team, while for GLASS we downloaded Hubble Multivisit Mosaics that overlapped with our JWST NIRCam observations.The CEERS mosaics are built from CANDELS observations in the Extended Groth Strip (Koekemoer et al. 2011).The GLASS HST mosaics are built from many legacy surveys (e.g., Frontier Fields Lotz et al. 2017) where all observations that overlap in cells projected on the sky are combined.These HST observations are aligned to the JWST observations using ProPane and the pixel scale is 0.06 arcsec/pixel for all HST and JWST imaging.
Source detection was performed on an all-filter inverse-variance weighted stack of the NIRCam images, spanning a wavelength range of ≈ 1.0µm → 4.5µm.ProFound was used to create the initial source detection map and extract multiband measurements from each individual filter.ProFound has been used to make the multiband catalogues for the Galaxy and Mass Assembly (GAMA) Survey (Bellstedt et al. 2020a) and the Deep Extragalactic VIsible Legacy Survey (DEV-ILS; Davies et al. 2021).
Briefly, ProFound calculates a grid of sky statistics to model the sky in our detection image and flux peaks are then first identified if they are above some user-defined multiple of the sky RMS and number of connecting pixels (in our case we use 1.5σ skyRMS and seven connecting pixels).ProFound then searches for connecting pixels to grow the segmentation map for the sources.Overblending is mitigated through the use of 'watershed' deblending in ProFound where flux peaks of neighbouring pixels need to be sufficiently distinct to be considered multiple sources (in our case blending only occurs if the saddle point in flux between two neighbouring sources is < 3σ skyRMS of the height of the highest flux peak).To capture all of the flux, ProFound then iteratively dilates the detected segments until the change in added flux with each dilation approaches zero.This differs from the commonly adopted methodology of SourceExtractor (e.g., Bertin & Arnouts 1996) as the segments are not confined to a fixed shape (they conform to the natural morphology of the galaxy) and they capture most, if not all, of the flux, eschewing any need for aperture corrections.
Finally, to obtain realistic photometric errors, circular apertures of varying radii were placed in blank parts of each of our multiband images and the flux errors on similarly sized segments were scaled up by the value in the equivalent aperture.This step is important as spatial covariance introduced in the mosaicking stage is not captured by the Poisson error calculated inside of the galaxy segment.This sky sampling also allows us to ac-count for cross terms in the error budget induced from detector and sky fluctuations.Additionally, we quadrature add a 10 per cent uncertainty on top of the scaled ProFound errors to further account for potential systematic errors.Ultimately, we aim to establish credible errors on our photometry to obtain robust constraints on the final physical parameters derived from SED fitting.All of our SEDs span from ≈ 0.6µm → 4.5µm.
The full suite of processing, source detection and multiband photometry codes, tailored for JWST, is referred to as Jumprope and is publicly available via Github (D 'Silva et al. 2023) 4 .

ProSpect fitting
With our combined, multiwavelength catalogue we conducted SED fitting with ProSpect to determine physical quantities of the galaxies.A key advantage of ProSpect is its ability to distinguish the relative contributions to the SED from star formation and AGN.
Stellar emission is modelled by the Bruzual & Charlot (2003) templates with the Chabrier ( 2003) initial mass function and the metallicity-dependent line energy tables from Levesque et al. (2010) for nebular emission lines.Dust attenuation and reemission is modelled by the Charlot & Fall (2000) attenuation law and the Dale et al. (2014) far-infrared templates.Star formation histories (SFHs) are modelled parametrically as skewed normal distributions (eqs.1-5 in D 'Silva et al. 2023b).While parametric SFHs cannot capture random bursts of star formation, skewed normal distributions have been shown to follow the average shape of SFHs when applied to galaxy formation simulations (e.g., Robotham et al. 2020).A key advantage of ProSpect is its treatment of evolving metallicity, where metallicity is linearly mapped to the stellar mass growth and encapsulates chemical enrichment of the interstellar medium.The implementation of skewed normal parametric SFHs and physically motivated metallicity evolution in ProSpect is crucial for reproducing the correct CSFH when forensically constructed from z ≈ 0 SFHs, suggesting that ProSpect accurately recovers the epochs at which stars formed (Bellstedt et al. 2020b).AGN are characterised by the models of Fritz et al. (2006).Emission of the primary source is modelled as a combination of power laws.The dusty torus is modelled as a flared disk with a range of dust grain sizes.The opening angle and the angle of observation are also incorporated in the Fritz et al. (2006) models.ProSpect has previously been shown to accurately recover the star-forming main sequence (SFMS), stellar mass function (SMF) and AGN bolometric luminosity function from z ≈ 0 − 5 by e.g., Thorne et al. (2021Thorne et al. ( , 2022)); D 'Silva et al. (2023b).
ProSpect was specifically designed to fit broadband SEDs of large samples; as such, only the HST and JWST photometry, and not the NIRSpec spectra, were fitted.The redshift was fixed to be the spectroscopic redshift of the sample to effectively eliminate potential redshift uncertainty on the fitted astrophysical properties.
We used two flavours of ProSpect runs that are highlighted in Tab. 1. Essentially, Stellar+AGN accounts for a stellar and AGN component in fitting the SED, while Stellar only accounts for a stellar component.The idea with using both of these SED fits is to bracket the extent to which an AGN may contribute to the SED.
Fig. 1 shows both Stellar+AGN and Stellar fits to the confirmed broad-line AGN in CEERS first presented by Larson et al. (2023).We see that a comparison of the maximum log-probability (LP) values shows that there is a factor of 10 1.3 ≈ 20 times better fit probability when accounting for an AGN component.From the comparison of the observed photometry to the fit photometry in the bottom panel, there is overall good agreement highlighting that we are adequately fitting the SED of this galaxy.Larson et al. (2023) also provide estimates of the bolometric luminosity of the AGN, L AGN ≈ 10 45.1±0.2 erg s −1 , that is calculated from the broad component Hα flux.Comparing this to our measurement from the multiband photometry, L AGN ≈ 10 45.02±2.35erg s −1 , the agreement is remarkable.Although, we note that the uncertainty on the AGN bolometric luminosity from Stellar+AGN is much greater than that derived from the broad component of the Hα line in Larson et al. (2023).Ultimately, this consistency reassures us that ProSpect produces reliable SED fits when simultaneously accounting for AGN and stars.These ProSpect SED fits are then used to calculate the CSFH.

Calculating CSFH
Usually, estimates of the CSFH rely on integration of the SFR distribution function or, as often is the case at high redshift, integration of the UV luminosity function as a proxy for SFR (e.g., Adams et al. 2023;Harikane et al. 2023b;Bouwens et al. 2015).This approach requires that the distribution function be well known with careful consideration of the sample selection and effective volumes (e.g., Weigel et al. 2016).Essentially, the idea is to account for all sources of star formation in the Universe and not underestimate the star formation by missing the faint sources, e.g., Malmquist bias.
Unfortunately, the gain in using spectroscopically confirmed z ≳ 3.5 galaxies, as opposed to photometric candidates, is somewhat offset by the loss in knowledge of the selection criteria used to observe those galaxies.It is difficult to directly compute the effective volumes of this sample and, as such, hard to compute the distribution function.This difficulty is highlighted in Harikane et al. (2023a) who use a similar spectroscopic sample to ours to estimate the UV luminosity function and CSFH.In this work, we instead pursued an alternative approach to calculate the SFR density distribution and CSFH.The logic is where ϕ(M ⋆ ) is the SMF and f (M ⋆ ) is a function that describes the correlation between SFR and stellar mass.The SMF in Eq. 1 is used only to obtain the effective volumes of the galaxies.Thus, in this work, f (M ⋆ ) is the independent variable observed between Stellar+AGN and Stellar ProSpect runs, and the effect on the CSFH when AGN is taken into account can be explored.

Fitting the SFMS
Therefore, this alternative approach to derive the CSFH hinges on the existence of the SFMS, the tight sequence of SFR and stellar mass on which the bulk of galaxies tend to lie (Noeske et al. 2007;Salim et al. 2007;Speagle et al. 2014).
The left panels of Fig. 2 show the SFR-M ⋆ of our galaxies, separated into three equal bins of redshift, from z ≈ 3.5 − 12.5.In this figure we show the results computed from both Stellar and Stel-lar+AGN.In both runs we clearly observe the existence of the SFMS, and so we computed the function f (M ⋆ ) as a log-linear function of stellar mass, i.e.,   2023) report LAGN ≈ 10 45.1±0.2 erg s −1 , while we find LAGN ≈ 10 45.02±2.35erg s −1 .The bottom panel shows the error weighted difference between the observed and fitted photometry, with the same colour scheme used for the fitted SED curves above.
The motivation for this functional form is twofold.First, this function has the minimum degrees of freedom (dof = 2) to fit the SFMS.Second, previous works have used double power laws, which amount to two linear functions separated by a break mass in log-space (dof = 4), to fit the SFMS, but have shown that deviation from a single log-linear component only becomes significant at z < 2.6 (e.g., Thorne et al. 2021).We used the fully Bayesian fitting tool hyper-fit (Robotham & Obreschkow 2015) to fit the SFMS, where Markov Chain Monte Carlo is employed to sample the posterior distributions of the slope and intercept.Because of the limited number of data points in each bin, we first used linear smoothed-splines to first roughly estimate the slope and intercept, and then assumed Gaussian priors about each of those rough estimates of widths 0.5 and 1.0, respectively.The solid lines in Fig. 2 are from the most likely fit parameters, and the shaded regions are the 16th-84th percentile ranges of the posterior distributions, indicating that we are adequately fitting the SFMS.
While not the focus of this Letter, simply the existence of the SFMS at z ≈ 3.5 − 12.5 is striking.The SFMS is thought to be evidence for self-regulation in galaxies where episodes of gas accretion and condensation are balanced against stellar and/or AGN feedback (Tacchella et al. 2016 Thorne et al. 2021).The grey shaded region shows the exclusion range up to M⋆ = 10 8 M⊙ where the stellar masses may not be robust.The inset panel in the highest-redshift bin shows a zoom-in of the data points and the SFMS.Right: SFR density distributions per unit stellar mass at z ≈ 3.5 → 12.5.The colour scheme is the same as in the left panels.
Thus, we interpret the existence of the SFMS to be a result of regulation being established as early as z ≈ 10 (the median redshift of our highest-redshift bin), when the Universe was only ≈ 500 Myr old.The orange and green dashed lines in Fig. 2 are predictions of the SFMS at z ≳ 5 from the FLARES and SHARK simulations (e.g., D 'Silva et al. 2023a;Lovell et al. 2021;Vijayan et al. 2021;Lagos et al. 2018).The black line shows a log-linear fit to the z = 7 − 10 SFMS from Heintz et al. (2023), who also use JWST observations.It should also be mentioned that Nakajima et al. (2023) also show the existence of the SFMS at z ≳ 4 with JWST.The key point is that the existence of the SFMS within the first few hundred Myr since the Big Bang is supported by galaxy formation simulations and independent analyses of JWST observations, further encouraging us of our methodology.
We compare our results to the z ≈ 5 SFR-M ⋆ from the DEVILSD10 survey (Thorne et al. 2021), shown with grey points, where we see that our points exhibit similar trends, giving us confidence in the robustness of the SED fits.We do see that some of the Stellar+AGN points in Fig. 2 are hitting a limit that causes their final stellar mass to cluster at M ⋆ ≲ 10 7 M ⊙ .In these cases the SED fits of these galaxies are completely dominated by the AGN component.To ensure that our results are unbiased against potentially erroneous SED fits, we only fitted the SFMS with galaxies of M ⋆ ≥ 10 8 M ⊙ as a conservative threshold.

SFR distribution function
We drew on previous determinations of the SMF at z ≳ 3.5 to estimate the effective volumes of our spectroscopic sample.Using Spitzer /IRAC observations in the CANDELS and Hubble Ultra Deep fields, Stefanon et al. ( 2021) compute z = 6 − 10 SMFs as single-component Schechter functions (e.g., Schechter 1976;Weigel et al. 2016), and, as of yet, no other z ≈ 10 observations of the SMF exist besides those.Stefanon et al. (2021) only compute the SMF down to z ≈ 6 and so we require a similar SMF to use in our lowest-redshift bin, z = 3.5 − 6.5, to be consistent with the results in the two higher-redshift bins.
We first investigated using the z ≈ 5 − 8 SMFs from Song et al. (2016)  While exploring the differences between the underlying SMFs is not the focus of this work, we remark that the adopted SMF will affect the final CSFH.The origin of the difference likely resides in the assumed selection function and resulting effective volumes for the galaxies.These SMFs also used photometric redshifts, and so we are potentially introducing implicit photometric redshift uncertainty into our final CSFH.While it is not ideal to use a mixture of SMFs in different redshift bins, we stress that the aim of this work is to explore the relative difference between the Stellar+AGN and Stellar CSFH.As we used the same SMFs for each of these samples, the relative difference does not depend on the exact choice of SMF.Nevertheless, in a follow-up work, we intend to use our own photometric redshifts and selection functions based on joint HST +JWST photometry to compute the direct SFR density distribution (along with the SMF).In which case, we hope to establish convergence on the true CSFH in the z ≳ 3.5 Universe by means of many independent approaches.
The right panels of Fig. 2 show the SFR density distributions per unit stellar mass.As we wish to highlight the relative effect of the AGN component on the CSFH, we used the same SMFs without an assumed uncertainty for both the Stellar and Stellar+AGN results.Thus, the quantiles show the 16th-84th uncertainty intervals only obtained from the Monte Carlo sampling of the 1000 spline fits of the SFMS.
There are systematic differences on the final SFR density distribution between the Stellar and Stel-lar+AGN results, with the Stellar showing higher normalisation at M ⋆ ≈ 10 8 − 10 10 M ⊙ and a shallower turnover to lower stellar masses at all redshifts shown.We see that the 9.5 ≤ z < 12.5 bin exhibits diverging SFR density at M ⋆ ≲ 10 9 , meaning that we can only calculate an upper limit on the final CSFH.We note that the stellar mass density functions from Stefanon et al. (2021) are themselves diverging at low stellar masses.Nevertheless, this final SFR density satisfies the integrand of Eq. 1 needed to calculate the CSFH.

Cosmic star formation history
While we only fitted the SFMS for galaxies with M ⋆ ≥ 10 8 M ⊙ , we extrapolated the log-linear SFMS for the integration of the SFR density function.We integrated the SFR density curves in the domain 10 5 M ⊙ ≤ M ⋆ < 10 12 M ⊙ .The bounds of the integration were chosen to encompass a range of galaxies that contribute to the CSFH.The lower limit of 10 5 M ⊙ is motivated by stellar mass estimates of Segue 2 (Belokurov et al. 2009;Kirby et al. 2013) that is the least massive galaxy discovered, and is likely a tidally stripped dwarf galaxy.The upper limit is motivated by theoretical considerations of the most massive galaxy that can exist in ΛCDM (e.g., Lovell et al. 2023a).Note, that we find little enhancement in the CSFH when we instead integrate up to M ⋆ = 10 13 M ⊙ .
Fig. 3 shows the resulting CSFH obtained from our sample.As a validation of our method to calculate the CSFH, we used the z ≈ 0 − 4 SMFs and SFMSs from Thorne et al. (2021) and compare the final CSFH with the results of D 'Silva et al. (2023b)  range.The results from D 'Silva et al. (2023b) are direct detections of the CSFH where the full SFR distribution function is integrated.Both of these data sets were built from the same sample of GAMA+DEVILS galaxies and the same ProSpect fits meaning it is correct to compare them.Overall, the two results are consistent, and so we build confidence that this same method, using the SMF and SFMS, can recover the CSFH for our JWST sample.
Moving on, our results at z ≳ 3.5 agree with previous estimates from Bouwens et al. (2015Bouwens et al. ( , 2023a)); Harikane et al. (2023b); Adams et al. (2023).There is a tentative systematic offset between the Stellar and the Stel-lar+AGN results, ranging from ≈ 0.1 dex at z ≈ 5 to ≈ 0.4 dex at z ≈ 10.5, indicating that the AGN component can potentially significantly contaminate the resulting star formation rates in z > 3.5 galaxies.Estimates of the CSFH at z > 3.5 generally rely on the UV luminosity function.Indeed, this is the case for the z ≳ 10 estimates of the CSFH from Harikane et al. (2023b); Bouwens et al. (2023b); Adams et al. (2023).But the UV luminosity function is known to be derived from a combination of stellar and AGN UV photons (Finkelstein & Bagley 2022), and so using it to compute the CSFH hinges on a careful delineation between the two processes.

AGN growth at z ≳ 10
In light of this, we remark that there is a greater offset between the Stellar and Stellar+AGN estimates in the z ≈ 9.5 − 12.5 bin than at the lower-redshift bins.Though, we see a higher uncertainty in this highestredshift bin, on account of our SFR density distributions basically diverging at M ⋆ ≲ 10 8 M ⊙ .Nevertheless, this increasing offset at z ≳ 9.5 is reminiscent of the offset reported by Harikane et al. (2023b) between their CSFH and the constant efficiency models of the CSFH presented in Harikane et al. (2022), which we show by the solid line in Fig. 3.As such, intuiting these reported offsets to be of common origins, the higher star formation efficiency observed at z ≳ 10 may instead be AGN masquerading as star formation.
Indeed, JWST has confirmed that SMBHs are common at these redshifts (e.g., Larson et al. 2023;Kocevski et al. 2023;Übler et al. 2023), can accrete at super-Eddington rates (Schneider et al. 2023) and outpace the growth of the host galaxy's stellar mass (Maiolino et al. 2023).Furthermore, Harikane et al. (2023c) found 10 Hα broad-line (FWHM ≈ 1000 − 6000 km s −1 ) type-1 AGN at z = 3.8 − 8.9 from CEERS and GLASS NIR-Spec spectra.They suggest that the fraction of galaxies hosting type-1 AGN is higher at z = 4 − 7 than it is at z ≈ 0, and that AGN are abundant in the early Universe.We performed a coordinate match between our sample and theirs and found all possible matches.As we were able to find all type-1 AGN identified by Harikane et al. (2023c) that overlap with our imaging, this is further evidence of the high AGN contribution at these redshifts.
On the contrary, Harikane et al. (2023b) concluded that 9/10 of their z ≈ 12 − 16 galaxies cannot possess significant AGN, and thus cannot explain the offset with the constant efficiency model as a result of AGN contamination, on account of their extended morphologies.Indeed, recent results with JWST have shown that galaxy disks were already in place at z ≈ 5 − 8 (Ferreira et al. 2023;Ormerod et al. 2023).This implies that the relative weight to the entire morphology of a central point source (AGN) compared to an extended structure (SFR) would be low in these galaxies.This aligns with the low AGN fraction at z ≈ 5 reported in D 'Silva et al. (2023b), where AGN fraction is the ratio of AGN contributed flux between rest-frame 5−20µm to the total bolometric flux of the galaxy.However, this is only a relative assessment of star formation/AGN contribution: a galaxy at high redshift can have a low AGN fraction but a bolometrically luminous AGN.Thorne et al. (2022) suggest that galaxies with a significant AGN contribution have an AGN fraction of 10 per cent, meaning that significant AGN can even inhabit galaxies where 90 per cent of the light is constrained to an extended morphological structure.In other words, while extended morphology is a clear discriminant of AGN in cases where the AGN either contributes to all or none of the SED (e.g., AGN fraction = 1 or 0), its significance when we consider the union of star formation and AGN is less so.Thus, the spatial light distributions need not necessarily be centrally concentrated point sources for z ≳ 9.5 galaxies hosting powerful AGN because extended star formation can be just as active during the lead up to the peak of CSFH.

AGN effect on host galaxies
On the basis that our galaxies can in fact host powerful AGN, we further investigate the effect of the AGN on the properties of the host galaxies.
Fig. 4 shows the difference in stellar mass between the Stellar and Stellar+AGN SED fits to show what effect accounting for an AGN component has on the fitting.Galaxies with M ⋆ ≳ 10 8 M ⊙ computed from both ProSpect runs have similar derived stellar masses.To further explore the trends, we fitted the relations with third-order smoothing splines, where the three degrees of freedom allow us to see any inflection beyond a loglinear correlation.For the same reasons that we expressed in fitting the SFMS in Section 3.2.1,we only performed fitting for galaxies with M ⋆ ≥ 10 8 M ⊙ (from both ProSpect runs).Overall, we see that in the M ⋆ ≳ 10 9 M ⊙ quadrant at all redshifts, the results agree within uncertainties.There is a tentative inclination for Stellar+AGN to produce more massive galaxies, but this is observed only well into the extrapolated fit and so we cannot make robust conclusions here.
In the M ⋆ ≲ 10 9 M ⊙ quadrant we see that that relation turns up showing that Stellar+AGN produces systematically lower stellar masses than Stellar for these galaxies.This is consistent with the study of the bright (M UV = −21.6AB mag), z = 10.6 galaxy, GN-z11, whose nuclear morphological component contributes 2/3 of the total UV luminosity but only 1/5 of the stellar mass (Tacchella et al. 2023).Thus, for galaxies that host powerful AGN, which are convincingly abundant at z ≳ 5, the UV luminosity predominately originates from accelerating charged particles in the accretion disk around the SMBH and not from star formation, where only the latter process contributes to stellar mass.
Turning our attention to the exact difference in stellar mass between runs, accounting for an AGN, incredibly, can shift the final stellar mass by up to ≈ 4 dex.For these galaxies, the SED is basically fit entirely by the AGN component as evidenced by the greatest offset for galaxies containing the most luminous AGN.Though, we repeat that this is highly suggestive that the Stellar+AGN fitting is hitting a limit in stellar mass and SFR.The limit is due to the mass-to-light ratio approaching the minimum allowable value from ProSpect's set up of Bruzual & Charlot (2003) stellar population models when the SED is modelled overwhelmingly by an AGN component.As a consequence of this limit these fits may not be robust.Nevertheless, such a significant difference between ProSpect runs indicates to us that AGN is an important factor in interpreting the SEDs of z ≈ 3.5 − 12.5 galaxies.Tantalisingly, these results may offer an avenue to explain massive, quenched galaxies at z ≳ 3 whose formation times are challenging to reconcile with the stateof-the-art galaxy formation models (e.g., Lovell et al. 2023a;Carnall et al. 2023;Glazebrook et al. 2023).Sure enough, recent theoretical works have shown that quiescent galaxies can exist at these redshifts in quantities that agree with observations and that among the main drivers of star formation suppression is feedback from the AGN (Lovell et al. 2023b;Lagos et al. 2023).In essence, alternative prescriptions that account for the AGN, both from the perspective of SED modelling in the observations and subgrid physics implementation in the simulations, can explain the existence of these massive galaxies.
Offsets in stellar mass between the two ProSpect runs is the cumulative effect of the difference in the underlying SFR.In Fig. 5 we show the difference in SFR between the two runs as a function of AGN bolometric luminosity.A clear trend is observed where galaxies that Stellar+AGN has determined to host powerful, L AGN ≈ 10 43 erg s −1 , AGN experience as much as a ≈ 4 dex difference in SFR between runs.We also see that these galaxies have tiny stellar masses, M ⋆ ≲ 10 6.5 M ⊙ , when determined from the Stellar+AGN fits.This difference in SFR toward higher AGN bolometric luminosity is the genesis of the difference in our CSFH, seen in Fig. 3.
What does this mean in terms of the physics of star formation and AGN activity at z ≳ 3.5?Well, we interpret the offset in CSFH between the two ProSpect runs to mean that star formation does not overwhelm AGN activity (with respect to the relative contributions to the bolometric luminosity of the galaxy), especially at z ≳ 9.5.Indeed, as the gas supply is abundant at high redshift we can intuit that the same gas reservoir fuelling star formation is also feeding the AGN.This agrees with the picture outlined in D 'Silva et al. (2023b) who showed that the common depletion of the gas supply from z ≈ 2 → 0 produces a coeval decline in both the star formation and AGN activity.
Whether this coevolution persists at z ≳ 2 remains elusive.Recent results with JWST tend to indicate that SMBHs are overmassive compared to their host galaxy's stellar mass (e.g., Maiolino et al. 2023) on account of super-Eddington accretion onto the SMBH (e.g., Schneider et al. 2023), pointing to an unequal growth of stellar mass compared to black hole mass.Interestingly, these objects may follow a formation mechanism where massive, metal-poor gas clouds directly collapse into SMBHs whose accretion disk dominates the contribution to the SED compared to the stellar component (e.g., Bromm & Loeb 2003;Agarwal et al. 2013;Kroupa et al. 2020).
To further unveil the intricacies of star formation and AGN activity at z ≳ 2 − 5, in a future work, D 'Silva et. al (in preparation) 2023).The work presented in this Letter is thus a stepping stone and anchor point to address the coevolution of star formation and AGN within the first galaxies.

CONCLUSION
Central to the science aspirations of JWST is shedding light on the nature of stellar mass assembly and the growth of SMBHs within the first few billion years after the Big Bang.JWST has so far generated more questions than it has answered about the z ≈ 3.5 − 12.5 frontier.A key discussion freshly sparked by JWST is whether star formation is more efficient than we previously expected.In this work, we performed SED fitting on 109 spectroscopically confirmed z ≈ 3.5 − 12.5 to see if AGN masquerading as star formation can alleviate this tension.Our main findings are: • The CSFH is tentatively as much as ≈ 0.4 dex lower when computed using SFRs derived from the Stellar+AGN fits that simultaneously model the stars and the AGN.Thus, a hidden AGN component can potentially explain the excess of CSFH compared to constant star formation efficiency models.
• Stellar+AGN fits produce as much as ≈ 4 dex lower stellar masses and SFRs compared to Stellar.This is essentially the origin of the offset between the CSFH.Thus, a hidden AGN component may be a promising avenue to explain massive quenched galaxies at z ≳ 3 that are seemingly at odds with current theories of stellar mass assembly.However, the large difference between ProSpect runs is suggestive that we are hitting a limit in SFR and stellar mass in the SED fitting.Nevertheless, this highlights that AGN may be an important factor in interpreting the SEDs of z ≈ 3.5 − 12.5 galaxies.We hope that this work highlights the importance of considering the union of star formation and AGN activity within the first galaxies to ever exist.
6. DATA AVAILABILITY Analysis scripts used in this work are available at https://github.com/JordanDSilva/project2jwst csfh and here 10.5281/zenodo.10148456.SED fits, photometry and imaging used in this work will be available upon reasonable request to the corresponding author.JWST data used in this work can be found at 10.17909/t317-yp54 and 10.17909/akh9-5057.

Figure 1 .
Figure1.ProSpect fit to CEERS 1019, a confirmed broad-line AGN at z = 8.679 first presented inLarson et al. (2023).The black points with error bars show the ProFound photometry.The error bars show the 16th-84th percentile ranges of the photometry.The Stellar fit is shown in blue, while the Stellar+AGN fit is shown in red.The log-probability (LP), being the maximum of the log-posterior distribution, is also included.The blue and red squares show the fit photometry for Stellar and Stellar+AGN.We show the relevant filter transmission curves from HST +JWST with rainbow colours.Larson et al. (2023) report LAGN ≈ 10 45.1±0.2 erg s −1 , while we find LAGN ≈ 10 45.02±2.35erg s −1 .The bottom panel shows the error weighted difference between the observed and fitted photometry, with the same colour scheme used for the fitted SED curves above.

Figure 2 .
Figure2.Left: SFMS at z ≈ 3.5 → 12.5.Blue points and lines show the results from Stellar, while red shows the results from Stellar+AGN.Grey dots are the z ≥ 5 SFR − M⋆ from the DEVILS survey (e.g.,Thorne et al. 2021).The grey shaded region shows the exclusion range up to M⋆ = 10 8 M⊙ where the stellar masses may not be robust.The inset panel in the highest-redshift bin shows a zoom-in of the data points and the SFMS.Right: SFR density distributions per unit stellar mass at z ≈ 3.5 → 12.5.The colour scheme is the same as in the left panels.
that are computed from similar Spitzer /IRAC observations in the CANDELS and Hubble Ultra Deep fields asStefanon et al. (2021).We also investigated the z ≈ 5 − 7 SMFs fromWeaver et al. (2023) that use the COSMOS2020 catalogue.We compared these two SMFs with theStefanon et al. (2021) SMFs at z = 6 − 7, which is the redshift range where the three studies overlap.The normalisation and the pivot mass of the Weaver et al. (2023) SMFs, instead of theSong et al. (2016) SMFs, best agree withStefanon et al. (2021), as the systematic offset of these fitted parameters between studies is lower and the parameters themselves are better constrained byWeaver et al. (2023).The low-mass powerlaw slope fromSong et al. (2016), however, agrees better withStefanon et al. (2021).The low-mass power law slope is fixed byWeaver et al. (2023), being ≈ 30 per cent shallower than the results fromStefanon et al. (2021), and the resulting CSFH is ≈ 0.4 dex less than that computed with the Song et al. (2016) SMF.Because of the fixed slope inWeaver et al. (2023), we decided to simply use theSong et al. (2016) SMF to compute the CSFH in our lowest-redshift bin.

Figure 3 .
Figure 3. CSFH at z ≈ 0 − 15.Blue triangles show the z ≈ 3.5 − 12.5 Stellar results, while red inverted triangles show the results for Stellar+AGN.The error bars show the 16th-84th percentile from the SFMS fits in Fig. 2. The arrows indicate where the results are upper limits.The blue and red points are offset by 0.05 in z for clarity.The brown dashed line shows the Madau & Dickinson (2014) fitted CSFH.The brown solid line shows the constant star formation efficiency model of the CSFH from Harikane et al. (2022).The brown shaded region shows the CSFH calculated at z ≈ 0 − 4 using the SFMSs and SMFs from Thorne et al. (2021).The grey points show previous literature results of the CSFH from D'Silva et al. (2023b); Bouwens et al. (2015); Harikane et al. (2023b); Bouwens et al. (2023b); Adams et al. (2023).

Figure 4 .
Figure 4. Left: Stellar M⋆ against Stellar+AGN M⋆ at z ≈ 3.5 − 12.5.The colour of the points corresponds to the AGN bolometric luminosity determined by ProSpect, where blue is weaker AGN and red is stronger.The purple line and shading show the third-order spline fit.The grey shaded region shows exclusion range up to M⋆ = 10 8 M⊙ where the stellar masses may not be robust.The dashed line shows the one-to-one relation.The inset panel in the highest-redshift bin shows a zoom-in of the data points and fitted relation.Right: ∆M⋆ = Stellar+AGN − Stellar.The colour scheme is the same as the left panels.The horizontal dashed line shows ∆M⋆ = 0.

Figure 5 .
Figure 5. ∆SFR = Stellar+AGN − Stellar as a function of AGN bolometric luminosity.We show the results binned in 2 dex intervals of Stellar+AGN stellar mass from M⋆ ≈ 10 4.5 M⊙ → 10 10.5 M⊙.The colour of points indicates the redshift, where blue corresponds to low redshift and red to high.The horizontal dashed line shows ∆SFR = 0.

Table 1 .
Description of our two ProSpect runs.We refer to these descriptions throughout the text.