X-ray stacking reveals average SMBH accretion properties of star-forming galaxies and their cosmic evolution over 4<~ z<~ 7

With an X-ray stacking analysis of ~ 12, 000 Lyman-break galaxies (LBGs) using the Chandra Legacy Survey image, we investigate average supermassive black hole (SMBH) accretion properties of star-forming galaxies (SFGs) at 4<~ z<~ 7. Although no X-ray signal is detected in any stacked image, we obtain strong 3 sigma upper limits for the average black hole accretion rate (BHAR) as a function of star formation rate (SFR). At z ~ 4 (5) where the stacked image is deeper, the 3 sigma BHAR upper limits per SFR are ~ 1.5 (1.0) dex lower than the local black hole-to-stellar mass ratio, indicating that the SMBHs of SFGs in the inactive (BHAR<~1M_sun yr^{-1}) phase are growing much more slowly than expected from simultaneous evolution. We obtain a similar result for BHAR per dark halo accretion rate. QSOs from the literature are found to have ~ 1 dex higher SFRs and>~ 2 dex higher BHARs than LBGs with the same dark halo mass. We also make a similar comparison for dusty starburst galaxies and quiescent galaxies from the literature. A duty-cycle corrected analysis shows that for a given dark halo, the SMBH mass increase in the QSO phase dominates over that in the much longer inactive phase. Finally, a comparison with the TNG300, TNG100, SIMBA100, and EAGLE100 simulations finds that they overshoot our BHAR upper limits by<~ 1.5 dex, possibly implying that simulated SMBHs are too massive.


INTRODUCTION
A tight correlation between the supermassive black hole (SMBH) mass ( BH ) and the bulge mass ( bulge ) seen in local galaxies suggests that SMBHs and galaxies have co-evolved (Kormendy & Ho 2013;McConnell & Ma 2013).When did the correlation emerge?What physical mechanisms made this correlation?To answer these questions and unveil the history of co-evolution, it is needed to investigate the correlation between the two over cosmic time.
Beyond the local universe, however, observations are limited to active galactic nuclei (AGNs) because it is challenging to measure  BH for inactive galaxies owing to the absence of broad lines.It is also hard to measure  bulge for distant AGNs because of the difficulty in imaging their morphology, and thus, the total stellar mass ( ★ ) is often used as a proxy of it.Previous studies of  > 0 AGNs (including QSOs) have obtained mixed results that may depend on redshift.Studies up to  ∼ 2.5 have found no significant evolution of the  BH - ★ relation (e.g., Ding et al. 2020;Suh et al. 2020;Li et al. 2021)  1 .At higher redshift, SMBHs tend to be overmassive with respect to the local relation (e.g., Bogdan et al. 2023;Harikane et al. 2023;Maiolino et al. 2023;Übler et al. 2023), although there ★ E-mail: smatsui@astron.s.u-tokyo.ac.jp 1 Ding et al. (2020) have also measured  bulge and detected positive evolution of  BH / bulge .
An important problem in studies using AGNs is that the obtained results may be biased because AGNs are rare objects with high accretion rates.We cannot rule out the possibility that AGNs, in particular QSOs, are outliers of the correlation scattered toward high  BH (e.g., Lauer et al. 2007).Future observations of fainter sources will provide a less biased view of the correlation.
Given the difficulties in measuring  BH at  > 0, an alternative, although indirect, method for studying co-evolution that can also be applied to inactive galaxies is to use the growth rates of SMBHs and stellar components, i.e., BH accretion rate (BHAR) and star formation rate (SFR).If the local  BH −  ★ relation holds over cosmic time, BHAR/SFR ratios will be close to an average  BH / ★ value.Indeed, a similarity in the evolution of the star formation rate density (SFRD) and the BH accretion rate density (BHAD) seen at  ≲ 2.5 (e.g., Delvecchio et al. 2014;Vito et al. 2018) suggests some connection between star formation and BH accretion.Note, however, that this similarity does not necessarily mean the simultaneous evolution of galaxies and SMBHs on an individual basis because the SFRD and BHAD are integrated quantities over all galaxies.
BHARs at  > 0 can be estimated from X-ray luminosities with less contamination from star formation than at other wavelengths.
Even for galaxies whose X-ray emission is too faint to be detected, we can obtain an average BHAR by stacking their X-ray images and discuss co-evolution in an unbiased way.Stacking analysis has an additional advantage that it essentially provides a time-averaged BHAR of a given galaxy sample by smoothing out possible variability of individual sources on short timescales of 10 2 -10 7 yr (e.g., Novak et al. 2011;Hickox et al. 2014)  2 .
The BHAR -SFR relation is of particular interest at high redshifts where most galaxies, including massive ones, are star-forming galaxies (SFGs) (e.g., Davidzon et al. 2017, their Fig.16).In this study, we examine this relation for SFGs at 4 ≲  ≲ 7 by X-ray stacking.
Several papers have examined the relation between BHAR and host galaxy properties including SFR at redshifts similar to our redshift range by X-ray stacking of a large sample.Fornasini et al. (2018) have used ≃ 75, 000 SFGs at 0.1 <  < 5 in the Chandra COSMOS-Legacy Survey area to find that while BHAR strongly correlates with  ★ , there is no correlation with SFR.Carraro et al. (2020) have conducted a similar analysis on three types of galaxies (SFGs, starburst galaxies [SBGs], and quiescent galaxies [QGs]) at 0.1 <  < 3.5 of a total number of 97,746, finding that BHAR/SFR depends on  ★ and that it does not evolve with redshift.Ito et al. (2022) have examined stacked properties of SFGs and QGs at 0 <  < 5 of a total number of 322,743.They have found that QGs have higher BHAR/SFRs than SFGs, suggesting that AGN activity plays a crucial role in quenching at high redshifts.
In these previous studies, BHAR/SFRs have been found to be lower than expected from the local mass relation, indicating that the growth of  BH in inactive galaxies at the probed redshifts is slower than the prediction of simultaneous co-evolution.However, these studies have primarily focused on the relationship between BHAR and  ★ , without directly comparing and discussing the BHAR and SFR in relation to the local correlation.Their BHAR/SFR values may also be biased because they use stellar-mass limited samples.
In this study, we use a large ( ≃ 12, 000) Lyman-break galaxy (LBG) sample for stacking analysis.This sample is rest-frame FUV continuum magnitude limited, and hence, approximately SFR limited.In addition to SFR, we also examine the relation with the dark halo mass ( h ), or equivalently, the halo accretion rate (HAR), of hosting galaxies using clustering-based dark halo mass estimates available for the sample.Dark halo mass is a primary parameter that determines the evolution of both SMBHs and stellar components (Benson 2010;Vogelsberger et al. 2020); BHAR/ h (or BHAR/HAR) and SFR/ h (SFR/HAR) are measures of their massgrowth efficiencies in dark haloes.By comparing BHAR/SFR and BHAR/HAR ratios with the corresponding local mass relations, we constrain the evolution of SMBHs and host galaxies at earlier cosmic times than probed by the previous studies 3 .
We describe our galaxy sample in Section 2. The X-ray stacking analysis and BHAR calculation methods are described in Section 3. Results are given in Section 4, where we also compare the results with QSOs and with cosmological simulations of co-evolution.We discuss several implications of the results in Section 5. A summary is given in Section 6.Throughout this paper, we use a Chabrier (2003) initial mass function (IMF) and assume a flat cosmology with 2 Short time-scale variability is an uncertain factor when one uses the AGN X-ray luminosity function to constrain cosmological simulations that do not have sufficient time resolution (Habouzit et al. 2022  0 = 70 km s −1 Mpc −1 , Ω Λ = 0.7, and Ω M = 0.3.All magnitude are presented in the AB system (Oke & Gunn 1983).

Star-forming Galaxies Sample
We use the catalog of 12,128 LBGs at 4 ≲  ≲ 7 constructed by Harikane et al. (2022), considering them to be SFGs.These LBGs have been selected by the standard color selection technique (e.g.Steidel et al. 1996;Giavalisco 2002) from , , , , and -band data of the UD-COSMOS field taken in the Subaru/Hyper Suprime-Cam Subaru Strategic Program (HSC SSP) (Aihara et al. 2019).
We derive the SFR for each object as follows (Harikane et al. 2022).First, we calculate the UV absolute magnitude ( UV ) from the observed magnitude "convolvedflux_0_20_mag" (-band for  ∼ 4, -band for  ∼ 5, and -band for  ∼ 6, 7) in Harikane et al. (2022).Then, we estimate the dust attenuation corrected UV luminosity ( UV ) from  UV using the dust attenuation-UV slope relation by Meurer et al. (1999) ( 1600 = 4.43 + 1.99, where  1600 is the dust attenuation in the rest frame 1600 Å and  is the UV slope), where we estimate  from the UV-slope- UV relation shown in Bouwens et al. (2014).Finally, we convert  UV into SFR using the formula given in Madau & Dickinson (2014) (but after rescaling for a Chabrier IMF): SFR Harikane et al. (2022) also provide the average  h as a function of  UV estimated by clustering analysis in their eq.54.

Chandra X-ray Image
For X-ray stacking, we use the data from the Chandra COSMOS Legacy survey (Civano et al. 2016), which is a 4.6 Ms Chandra program that has imaged 2.2 deg 2 of the COSMOS field.The depth of the data reaches 2.2 × 10 −16 erg cm −2 s −1 in the 0.5-2 keV band and 1.5 × 10 −15 erg cm −2 s −1 in the 2-10 keV band.Fig. 1 shows the distribution of our LBGs (parent samples in Table .1) on the 0.5-2 keV image.

X-ray Stacking of Sub-samples
This study aims at constraining the average BH accretion properties of normal SFGs.Even a tiny portion of X-ray detected sources in a sample for stacking can greatly affect the average trend.By crossmatching the LBG catalog with the Chandra COSMOS Legacy Survey source catalog (Civano et al. 2016) with a matching radius of 2 ′′ , we find that 18 out of 9231 LBGs at  ∼ 4 and two out of 2665 LBGs at  ∼ 5 are detected while finding no detected source at  ∼ 6 or 7. We exclude these 20 sources from our X-ray stacking.Their 2-10 keV luminosities after hydrogen absorption correction in a similar manner to stacked samples (see Section 3) are distributed over 1.45 × 10 44 − 9.11 × 10 44 erg s −1 , corresponding to BHARs of 0.37 − 2.31  ⊙ yr −1 .
We also exclude sources brighter than  UV = −23.This is because the UV spectrum of such bright sources may be significantly contaminated by an AGN (Harikane et al. 2022), making  UV -based SFR estimation unreliable.X-ray stacking is performed separately for  ∼ 4, 5, 6, and 7.For  ∼ 4 and 5, we also stack magnitude-divided sub-samples with a 1 mag bin.Table 1 summarises the number of sources used for X-ray stacking as well as those excluded.This paper defines these Chandra X-ray undetected galaxies as 'inactive galaxies'.They constitute the vast majority of the galaxy population at the redshifts in question.

Stacking Procedure
We use the Chandra stacking tool CSTACK v.4.4 4 (Miyaji et al. 2008).CSTACK provides an exposure-time weighted average count rate (with error) at the positions of objects given in an input source catalog after excluding objects that are > 8 ′ from the optical axis or are affected by resolved X-ray sources.An average count rate (CR), which is defined so that 90% of the total counts are included if stacked objects are PSF-like, is calculated for the soft (0.5-2 keV) and hard (2-8 keV) bands separately.The error in the average CR is calculated using bootstrap resampling.CSTACK also generates a stacked 30 ′′ × 30 ′′ image for each band.See, e.g., Ito et al. (2022) for a more detailed explanation of CSTACK and an application to inactive high-redshift galaxies.
We apply CSTACK to the whole sample at each redshift as well as the magnitude-divided sub-samples at  ∼ 4 and 5.No significant (≥ 3) signal is detected in any sample (Fig. 2; only results for the soft band are shown).Therefore, for each sample, we use a value three times the CR error (⟨CR 3 err ⟩) as the 3 upper limit of the average CR of the sources.Because the soft band gives a stronger upper limit to the BHAR than the hard band, we use the results from the soft band in this paper.
(ii) Then, we convert ⟨ 3 X ⟩ to the absorption-corrected average X-ray luminosity in the rest frame 2-10 keV (⟨ 3 2−10keV ⟩) by: where  L is the luminosity distance at the average redshift (⟨⟩) of the sample, and  1 and  2 are the soft-band range (i.e.,  1 = 0.5 keV and  2 = 2 keV).We find that the ⟨ 3 2−10keV ⟩ of each sample is 1.5-2.0times larger than that before internal  H absorption correction.  1) are plotted.
(iii) X-ray binaries (XRBs) can significantly contribute to the ⟨ 3 2−10keV ⟩ of inactive galaxies.We remove the contribution of XRBs from ⟨ 3 2−10keV ⟩ using an empirical relation given in Lehmer et al. (2010): where ⟨ ★ ⟩ and ⟨SFR⟩ are the average  ★ and SFR of the sample, respectively, and obtain the average AGN luminosity, ⟨ 3 X,AGN ⟩, as: We find that the calculated  X,XRB is a factor of 2.6−22 lower than ⟨ 3 X,AGN ⟩, implying that our results described in the next section are robust against the uncertainty in the XRB correction.We also note that using other empirical relations gives similar results.
(iv) We then derive the average bolometric luminosity upper limit (⟨ 3 bol ⟩) by multiplying ⟨ 3 X,AGN ⟩ by the conversion factor,  bol (≡  bol / X ):  (v) Finally, we calculate ⟨BHAR 3 ⟩ as: where  is the speed of light and  is the radiation efficiency of SMBH mass accretion.We adopt  = 0.1 following the previous X-ray stacking studies (e.g.Fornasini et al. 2018;Carraro et al. 2020;Ito et al. 2022).

X-ray Individual Detection Sources
We also calculate the BHAR for the X-ray individual detection sources.There are 18 X-ray detection sources in the  ∼ 4 bin and two sources in the  ∼ 5 bin.Among them, four sources in the  ∼ 4 bin and one source in the  ∼ 5 bin have spectroscopic redshifts within their LBG selection windows, five sources in the  ∼ 4 bin have photometric redshifts within the selection windows, and the remaining sources have redshifts outside the selection windows.We only derive the BHAR for the ten sources within the selection windows. Average stellar mass of the sample calculated from the average SFR using the SFR- ★ relation of main sequence galaxies (Pearson et al. 2018). Average dark halo mass of the sample. Average halo accretion rate derived from ⟨ h ⟩ by eq.6 in Section 4.2. Average X-ray luminosity 3 upper limit. Average AGN X-ray luminosity 3 upper limit.ℎ Average AGN bolometric luminosity 3 upper limit. Average BHAR 3 upper limit.

BHAR versus SFR
In Fig. 3, we plot ⟨BHAR 3 ⟩ against ⟨SFR⟩ for our stacked samples together with QSOs at similar redshifts in the literature.The dotdashed line in each panel corresponds to the local  BH / bulge ratio calculated from the data of early-type galaxies given in Kormendy & Ho (2013) (their Tables 2 and 3), after correction for the mass return fraction ( = 0.41 for a Chabrier IMF), which is the fraction of stellar mass that is ejected back to the interstellar medium (e.g., Ueda et al. (2018)).If SMBHs and their host galaxies are growing in tandem, they are on this line with BHAR/SFR = (1− )( BH / bulge ) local = 4.7 × 10 −3 , where ( BH / bulge ) local = 8.0 × 10 −3 .Here we have assumed that at  = 0, our LBGs have evolved into early-type galaxies with  bulge ≈  ★ and used galaxies with  h ≥ 1 × 10 13  ⊙ for the calculation.This assumption is reasonable because the presentday descendants of our LBGs (11.4 < log(⟨ h ⟩/ ⊙ ) < 13.0) are expected to have  h ≳ 1 × 10 13  ⊙ (e.g.Behroozi et al. 2013), at which masses most of the galaxies hosted are quenched (Behroozi et al. 2019).The ratio is not sensitive to the  h range used for the calculation because the local  BH - bulge relation by Kormendy & Ho (2013) is close to linear ( BH ∝  1.17 bulge ).We find that at  ∼ 4 and 5 the ⟨BHAR 3 ⟩ of our LBGs are ∼ 1.5 dex and ∼ 1 dex lower than the dot-dashed line, respectively, implying that the SMBHs of SFGs at these redshifts are growing much more slowly than expected from simultaneous evolution.At  ∼ 6, 7, the upper limits are only up to ∼ 0.5 dex below the line, possibly because of the small numbers of stacked sources.
QSOs are plotted at  ∼ 5 and 6, 7.Those at  ∼ 5 are bright sources selected from SDSS QSOs while those at  ∼ 6, 7 are discovered by several wide-field optical surveys and include faint sources from the HSC SSP.Their BHARs have been estimated from rest-frame 3000 Å or 1450 Å luminosities in Netzer et al. (2014) and from restframe 1450 Å luminosities in Trakhtenbrot et al. (2017); Nguyen et al. (2020); Venemans et al. (2018) while their SFRs from IR luminosities.We find that the QSOs have ∼ 1 dex higher BHAR/SFR ratios than the  = 0 mass ratio, implying that the SMBHs in QSOs are growing much faster than expected from simultaneous evolution.As a result, the difference in BHAR/SFR between LBGs and QSOs amounts to ∼ 2 dex.
In the left panel of Fig. 3, we also plot the results of three previous stacking studies at their highest-redshift bins: ⟨⟩ = 2.6−2.8(Carraro et al. 2020), ⟨⟩ = 3.2 − 3.6 (Fornasini et al. 2018), and ⟨⟩ = 3.5 − 3.6 (Ito et al. 2022).All three use the same Chandra data as this study.Fornasini et al. (2018)'s and Ito et al. ( 2022)'s estimates and upper limits are lower than our upper limits.This is probably because they can reach fainter luminosity limits at lower redshifts and, in some cases, with larger sub-samples than ours.Carraro et al. (2020)'s results are higher than these two studies, being closer to our upper limits, probably because they have included sources with Xray individual detection.Note that we cannot accurately compare the results of these three studies with ours because they use stellar-mass limited samples and constructed sub-samples by splitting the parent sample by specific SFR and  ★ , not by SFR.

BHAR versus HAR
We then compare the mass growth rate of SMBHs with that of hosting dark haloes (HAR: halo accretion rage) in Fig. 4. We estimate the ⟨HAR⟩ of each LBG sample from its ⟨ h ⟩ using the "mean" version of Fakhouri et al. (2010)'s eq.2: which has been obtained from Millennium and Millennium-II simulation results (Springel 2005;Boylan-Kolchin et al. 2009).The ⟨ h ⟩ of each sample is calculated from the  UV −  h relation in Harikane et al. (2022).The dot-dashed line in each panel corresponds to the local  BH / h ratio (BHAR/HAR = (1−)( BH / h ) local = 1.2×10 −4 ).We obtain this ratio using the data of local galaxies given in Kormendy & Ho (2013) (their Tables 2 and 3), limiting to those with  h ≥ 1 × 10 13  ⊙ for the same reason as for  BH / bulge .
QSOs  2018)'s QSOs are taken from Shimasaku & Izumi (2019).Shaded areas correspond to QSO samples with clustering-based  h estimates (Shen et al. 2007;He et al. 2018;Timlin et al. 2018).We find that the QSOs are mostly above the  = 0 mass ratio except for Arita et al. (2023)'s sample at  ∼ 6, with individual objects being distributed ∼ 1 − 2 dex above the line, implying that the SMBHs in QSOs are growing faster than hosting dark matter haloes.This result appears to be qualitatively in accord with observations that highredshift QSOs have overmassive SMBHs with respect to the local  BH - h relation (Trainor & Steidel 2012;Shimasaku & Izumi 2019).The difference in BHAR/HAR between LBGs and QSOs reaches ∼ 2 dex or more.

SFR versus HAR
We also examine SFR versus HAR for LBGs and QSOs in Fig. 5.We convert the  = 0  bulge / h ratio (Kormendy & Ho 2013) to SFR/HAR = (1 − )( bulge / h ) local = 2.5 × 10 −2 in a similar manner to BHAR/SFR, where ( bulge / h ) local = 1.5 × 10 −2 .We find that the ⟨SFR⟩/⟨HAR⟩ of LBGs agree with this ratio within ∼ 0.3 dex except for the most massive bin at  ∼ 4.This means that in SFGs at these redshifts the stellar component and the dark halo are growing roughly in tandem ('co-evolving') on average, in contrast to the stellar component and the SMBH.On the other hand, most QSOs have higher SFR/HAR than the local mass ratio.The median SFR/HAR of  ∼ 5 ( ∼ 6, 7) QSOs is found to be ∼ 1 dex (∼ 0.5 dex) higher than the local mass ratio.

Comparison with Large-scale Cosmological Simulations
We compare our results with the TNG300, TNG100 (Marinacci et al. 2018;Naiman et al. 2018;Nelson et al. 2018;Pillepich et al. 2018;Springel et al. 2018), SIMBA100 (Davé et al. 2019), and EAGLE100 (Crain et al. 2015;Schaye et al. 2015;McAlpine et al. 2016) simulations, all of which are large-scale cosmological hydrodynamical simulations and are open to the public.The TNG300 and TNG100 assume  = 0.2, whereas the SIMBA100 and EAGLE100 adopt    = 0.1.Because our ⟨BHAR 3 ⟩ values are for  = 0.1, we recalculate them with  = 0.2 when comparing with TNG300 and TNG100.For the accretion model, all but SIMBA adopt the Bondi-Hoyle-Lyttleton formalism or a modified version of it, while SIMBA uses a torque-limited accretion model for cold gas and the Bondi-Hoyle-Lyttleton model for hot gas.For a brief summary of these simulations (e.g., box size, spatial and mass resolutions, BH seed mass, model of accretion and feedback), see, e.g., Habouzit et al. (2021).

Comparing BHAR versus SFR with simulations
First, we compare our ⟨BHAR 3 ⟩ versus ⟨SFR⟩ with the simulations in Fig. 6.As mentioned earlier regarding the difference in , we divide the plot into two parts: TNG300 and TNG100 in one, and SIMBA100 and EAGLE100 in the other.
We first look at the results for  ∼ 4 where our upper limits are most stringent.We find that all four simulations predict much higher BHARs than the data.The largest discrepancy is seen in TNG100, with ∼ 1−1.5 dex overprediction at SFR ≃ 10−100 ⊙ yr −1 , but the other three simulations also overshoot the data by ∼ 1 dex.Note that our data are upper limits, meaning that true differences will be even larger.The discrepancy is reduced at  ∼ 5, but TNG100 still predicts ∼ 1 dex higher values.At  ∼ 6, 7, all simulations are consistent with our 3 upper limits.

Comparing 𝑆𝐹 𝑅 versus 𝑀 h with Simulations
To constrain the cause of the overprediction found above, we compare the ⟨SFR⟩ versus ⟨ h ⟩ of our LBGs with the simulations in Fig. 7.We find that the simulations correctly predict ⟨SFR⟩ as a function of ⟨ h ⟩, with differences from the observations within ∼ 0.3 dex.Note also that the simulations are consistent with Behroozi et al. (2019)'s empirical relations.Therefore, the discrepancy seen in BHAR versus SFR is not due to the star formation physics but the modeling of BH accretion.In the inactive phase, the simulations feed the central SMBH in a given dark halo too efficiently while producing stars in accord with observations.

DISCUSSION
In the previous section, we find that at 4 ≲  ≲ 7, our LBGs' ⟨BHAR 3 ⟩/⟨SFR⟩ and ⟨BHAR 3 ⟩/⟨HAR⟩ are ∼ 1−1.5 and ∼ 1−2 dex lower than the corresponding  = 0 mass ratios.This result leads us to conclude that at these redshifts, SFGs cannot obtain as high  BH / bulge and  BH / h as the local mass ratios without experiencing a rapid accretion phase (e.g., QSO phase).In Section 5.1, we compare the contribution of SFGs and QSOs to the growth of SMBHs.In Section 5.2, we discuss the star formation and BH accretion activities of SFGs, QSOs, and two other galaxy populations among which an evolutionary link has been suggested.We also find that the simulations overpredict the inactive galaxies' BHAR in Section 4.4.1, so we discuss implications of this result in Section 5.3.MeanBHAR as a function of SFR predicted by TNG300 (brown dotted lines) and TNG100 (brown solid lines) with  = 0.2 for  ∼ 4 (left panel), 5 (middle), and 6.5 (right).Bottom panels.Same as the top panels but for SIMBA100 (magenta lines) and EAGLE100 (limegreen lines at  ∼ 4, 5, 6 and a green line at  ∼ 7) with  = 0.1.Since EAGLE100 does not provide a  = 6.5 snapshot, we plot  = 6 and 7 results.When calculating these functions, we only use SFR > 1 ⊙ yr −1 sources for TNG300, TNG100, and EAGLE, SFR > 10 ⊙ yr −1 sources for SIMBA100 to ensure that each simulation is complete in terms of SFR.We examine the completeness using the SFR versus  ★ plot.The meaning of diamonds with down arrows and dashed and dot-dashed lines is the same as in Fig. 3.

Contribution of SFGs and QSOs to SMBH Growth
We compare the contribution of SFGs and QSOs to SMBH growth using our LBG samples and literature QSO samples with duty-cycle (  duty ) estimates, where  duty is defined as the fraction of QSOhosting dark haloes among all haloes over a certain halo mass range.Only He et al. (2018)'s andShen et al. (2007)'s samples at  ∼ 4 have simultaneous estimates of BHAR,  h (hence HAR), and  duty .The former sample consists of faint sources with a median of  1450 = −23.49while the latter consisting of bright sources with a median of  1450 = −26.35.The mean  h (with uncertainty) of these samples, log( h / ⊙ ) = 12.2 − 12.5 for He et al. (2018) and log( h / ⊙ ) = 12.8 − 12.9 for Shen et al. (2007), are roughly comparable to those of our  ∼ 4 LBGs.
In Fig. 8, we shift the two QSO samples' original BHAR ranges downwards by multiplying them by  duty = 0.03 (He et al. 2018) and 0.17 (Shen et al. 2007), respectively, to show effective, or time-  2007)'s with a duration indicated by each  duty estimate.We find that the effective BHARs are still higher than the ⟨BHAR 3 ⟩ of LBGs with similar HAR by ∼ 0.7 dex (He et al. 2018) and ∼ 2 dex (Shen et al. 2007).Although the  duty estimates have large uncertainties (∼ 1 − 1.5 dex [see Fig. 11 of He et al. (2018)]), the  BH increase in the QSO phase is likely to dominate over that in the much longer, inactive (BHAR ≲ 1 ⊙ yr −1 ) period.We also find that the effective BHAR/HAR of Shen et al. (2007)'s sample is comparable to the local mass ratio while that of He et al. (2018)'s is much lower.This implies that if a galaxy experiences a bright QSO phase like Shen et al. (2007)'s QSOs, the mass increase in its SMBH is roughly consistent with simultaneous evolution.

Star Formation and BH Accretion Activities of Various Galaxy Populations
Hopkins et al. ( 2008) have proposed a co-evolution model that links SFGs, starburst galaxies (SBGs), QSOs, and elliptical galaxies on the assumption that major, gas-rich mergers cause QSO activity.Regardless of the correctness of this model, it is interesting to quantitatively compare star formation and BH accretion activities between these populations at the redshifts covered by our LBG sample.To do so, we use SFR/ h and BHAR/ h because both stars and SMBHs grow in their common hosting dark haloes.We consider LBGs to be normal SFGs in the star-formation main sequence.For SBGs, we use the catalog of IR-detected galaxies in COSMOS provided by Jin et al. (2018), in which SFRs and AGN luminosities have been estimated by SED fitting to   to IR (or radio) photometry.They are essentially dusty SBGs with SFR ≃ 80 − 2000 ⊙ yr −1 , typically being located ∼ 1 dex above the star-formation main sequence.For ellipticals, we use X-ray stacked BHAR estimates for QGs at  ∼ 2 − 5 given by Ito et al. ( 2022) although their average redshift,  ∼ 3.2 − 3.4, is slightly lower than ours.We estimate the  h of SBGs and QGs from their  ★ using the conversion formula given by Behroozi et al. (2019).
The bottom panels of Fig. 9 plot BHAR/ h against SFR/ h for the four populations (LBGs [SFGs], SBGs, QSOs, and QGs), obtained from the top (BHAR versus  h ) and middle panels (SFR versus  h ).The QSOs plotted here have −25.86 >  1450 > −27.87 at  ∼ 5 (Trakhtenbrot et al. 2017;Netzer et al. 2014;Nguyen et al. 2020) and −22.83 >  1450 > −29.3 at  ∼ 6, 7 (Shimasaku & Izumi 2019, and reference therein).The red solid line in the left panel represents the average relation for all  ∼ 4 galaxies with  h = 10 11.5 − 10 13.3  ⊙ calculated from the empirically obtained SFR versus  h relation (Behroozi et al. 2019) and BHAR versus  h relation (Yang et al. 2018).Three dot-dashed lines in each panel correspond to the local  BH / bulge ,  BH / h , and  bulge / h ratios from Kormendy & Ho (2013), respectively.If, for example, a source is located above the  BH / h line, then the SMBH of this source is growing faster than the dark halo relative to simultaneous evolution.
We find that the QSOs and SBGs are located in a similar region, with ∼ 1 dex higher SFR/ h and ≳ 2 dex higher BHAR/ h than the LBGs, indicating that in QSOs and SBGs, BH accretion activity is more enhanced than star formation activity with respect to SFGs (For the QSOs, we are focusing on the median values of all objects at each redshift).The ∼ 1 dex enhancement of SFR/ h in QSOs and SBGs compared with LBGs implies that if every SFG experiences a QSO phase and a starburst phase like plotted here, the stellar mass acquired in these phases will dominate over that by normal star formation if the  duty of these phases is higher than ∼ 10%.
A detailed comparison between the QSOs and SBGs shows that, on average, the SBGs have slightly higher SFR/ h at a fixed BHAR/ h .This result is consistent with the finding by Andonie et al. (2022) for  < 3 infrared QSOs that obscured ( H > 10 22 cm −2 ) systems have ≈ 3 times higher SFRs than unobscured ones with similar  ★ .The BHAR/SFR of our SBGs is comparable to the local  BH / bulge ratio, meaning that their SMBH and stellar component are growing at paces roughly consistent with simultaneous evolution.
We also find that the QGs have comparable or higher BHAR/ h to the LBGs despite having ∼ 1 dex lower SFR/ h .If these QGs have evolved from SBGs and QSOs like plotted here, the increase in  BH in QGs appears to be insignificant unless the  duty of SBGs and QSOs is less than 1%.
We then compare  H among the four populations.The  H of our LBGs are predicted to be in the range log( H /cm −2 ) = 23.25−23.96(see Section 3.1).The  H of optically-selected QSOs (broad-line QSOs) are log( H /cm −2 ) ≲ 22 (e.g., Hasinger 2008;Martocchia et al. 2017), while the QGs have log( H /cm −2 ) ∼ 23.5 (Ito et al. 2022).To estimate the  H of the SBGs, we perform an X-ray stacking analysis in a similar manner to the LBGs.No significant signal is detected at any redshift bin.By comparing the average 3 flux upper limits thus obtained (log⟨ 3 X,stacked ⟩) with the average fluxes expected from super-deblended SED-fitting based AGN bolometric luminosities (log⟨ X,expected ⟩), we obtain lower limits of  H as shown in Table 3.The SBGs have extremely high column densities of log( H /cm −2 ) ≳ 24.6 − 24.9, being an order of magnitude higher than expected from Gilli et al. ( 2022)'s relation (log( H /cm −2 ) = 23.25 − 23.96) and well into the Compton-thick regime.Thus, the  H of the four galaxy populations plotted here are different by at least ∼ 2.5 dex although this large difference may be partly due to selection biases of the QSOs and SBGs.
The distribution of the four populations in Fig. 9 and their  H values are qualitatively consistent with Hopkins et al. (2008)'s model in which major mergers of isolated SFGs end up with ellipticals after experiencing a dusty starburst phase with obscured AGN and a subsequent QSO phase.

Implications of Too High BHARs in Simulations
In Fig. 6, we find that all four simulations predict too high BHARs.The sub-grid models of BH accretion of these simulations use the Bondi-Hoyle-Lyttleton model (see Table 1 of Habouzit et al. (2021)): or a modified version of it, where  is the gravitational constant and  and   are, respectively, the density and sound speed of the gas around the BH.SIMBA uses a torque-limited accretion model for cold gas as well.Because all simulations similarly overpredict the observed BHARs despite using different detailed formulations of accretion including the calculation of  and  s , a possible cause of the overprediction would be that simulated SMBHs are too massive.
On a simple assumption of BHAR simul /BHAR obs = ( simul BH / true BH ) 2 (including SIMBA, in which the torque model for cold gas dominates at high redshifts (Habouzit et al. 2021)), we calculate the average 3 upper limit of the true  BH (⟨ 3 BH ⟩) value at each of our ⟨BHAR 3 ⟩ data points.Fig. 10 shows the results for  ∼ 4 where the discrepancies of BHAR are largest.The true ⟨ 3 BH ⟩ of our LBGs at  ★ ≳ 10 10  ⊙ are expected to be > 3−5 times lower than simulated.Because the simulated  BH −  ★ relations are below the local relation, this result suggests that at  ∼ 4, the SMBHs of SFGs are greatly undermassive compared with simultaneous evolution.On the other hand, overmassive SMBHs have been found in high-redshift QSOs as mentioned in Section 1.The distribution of  BH at a fixed  ★ may be much wider than expected.The discussion here is sketchy since simulated galaxies will not exactly follow BHAR ∝  2 BH even at fixed gas properties especially in SIMBA, and since decreasing  BH will change gas properties accordingly.In addition, we cannot rule out other causes.For example, the accretion rate decreases if the SMBH has a velocity, , relative to the ambient gas, meaning that the denominator of eq.7 becomes (c 2 s +  2 ) 3/2 .Although  is considered only in EAGLE (Crain et al. 2015), it might be non-negligible or underestimated.A magnetic field in the gas may also reduce the accretion rate (e.g., Kaaz et al. 2023).Investigating these causes is, however, beyond the scope of this paper.In any case, the simulations appear to need some improvement to reproduce our low ⟨BHAR 3 ⟩.
In the above comparison, we have adopted a constant radiative efficiency of  = 0.1 or 0.2.However,  may be lower for low activity (i.e., low Eddington ratio,  edd =  bol,AGN / edd , with  edd = 1.26× 10 38 ( BH / ⊙ ) erg s −1 ) systems.To examine the influence of such dependence on the comparison with simulations, we calculate  bol for simulated galaxies considering  edd dependence of  and compare them with observed ⟨ 3 bol,AGN ⟩ values in Fig. 11.We have used the formula of  ( edd ) given by Churazov et al. (2005)  √︁ Ω M (1 + ⟨⟩ ) 3 + Ω Λ , so that HAR is exactly proportional to  h .This approximation is accurate within 0.1 dex over  h = 10 11 − 10 13  ⊙ in which almost all sources plotted here are included.Purple dot-dashed lines indicate the local  BH / bulge ratio (same as in Fig. 3), navy dot-dashed lines the local  BH / h ratio (same as in Fig. 4), and black dot-dashed lines the local  bulge / h ratio (same as in Fig. 5).A solid red line in the left panel shows the average relation of all galaxies over log(  h / ⊙ ) = 11. Average AGN X-ray luminosity calculated as ⟨ AGN,SED ⟩  / bol,SED  . Expected average X-ray flux calculated from ⟨ X,SED ⟩  and ⟨⟩.
Average 3 upper limit of  X obtained by X-ray stacking analysis, after correction of the Galactic hydrogen absorption.
SBGs'  H needed to achieve ⟨ x,expected ⟩  .open symbols are for total samples while filled symbols for magnitude divided samples) calculated on a simple assumption of ⟨BHAR⟩ simul /⟨BHAR⟩ obs = (  simul BH / true BH ) 2 for the four simulations (TNG300, TNG100, SIMBA100, and EAGLE100 from left).We adopt the same color choice as in Fig. 6 to represent each simulation.Dot-dashed lines indicate the local  BH - bulge relation (Kormendy & Ho 2013).In SIMBA100, the all, 25-26 mag, and 26-27 mag samples are not plotted because SIMBA100 is incomplete at SFRs corresponding to these faint apparent magnitudes.Similarly, in EAGLE100, the 23-24 mag sample is not plotted because EAGLE100 does not contain such bright (i.e., high-SFR) galaxies.
(i)  edd > 0.1: (ii)  edd ≤ 0.1: We find that the simulations again overshoot the data, although with slightly smaller discrepancies than seen in BHAR.In fact, at  ∼ 4−7, even inactive (BHAR ≲ 1 ⊙ yr −1 ) simulated galaxies like our LBGs have relatively high  edd .Thus, the conclusions obtained from the comparison of BHARs with simulations are not sensitive to  edd dependence of .Schirra et al. (2021) have found that all four simulations they have examined (Illustris, TNG100, EAGLE, and SIMBA) overestimate the  tot (=  X,AGN +  X,XRB ) of SFGs at high redshift ( ∼ 3, 4) by comparing these simulations with Fornasini et al. (2018)'s SFG results in the  tot - ★ plane.Habouzit et al. (2022) have compared six cosmological simulations (Illustris, TNG100, TNG300, Horizon-AGN, EAGLE, and SIMBA) with observed AGN X-ray luminosity functions, finding that all but EAGLE produce too many AGNs with low X-ray luminosities of  X,2−10keV ∼ 10 43−44 erg s −1 .These findings appear to be qualitatively in line with what we find.These 'overaccretion' problems may have the same origin.

SUMMARY
With an X-ray stacking analysis on a large ( ≃ 12, 000) LBG sample, we have constrained the average BHAR of inactive SFGs at 4 ≲  ≲ 7 as functions of SFR and HAR.We have also compared the results with observations of QSOs, SBGs, and QGs at similar redshifts and with cosmological simulations of co-evolution.Our main results and discussion are summarized as follows: (i) At  ∼ 4 and 5, the 3 BHAR upper limits of LBGs are ∼ 1.5 dex and ∼ 1 dex lower than expected from the local  BH / ★ ratio, respectively, implying that the SMBHs of SFGs at these redshifts are growing much more slowly than expected from simultaneous evolution.At  ∼ 6 − 7, the upper limits are closer to the local mass ratio, possibly because of the small numbers of stacked sources.
(ii) Similarly, at  ∼ 4, LBGs' ⟨BHAR 3 ⟩/⟨HAR⟩ is ∼ 1 − 2 dex lower than the  = 0  BH / h ratio.The difference is smaller at higher redshifts, but the ⟨BHAR 3 ⟩/⟨HAR⟩ is still ∼ 1 dex lower at  ∼ 5 and ∼ 0.5 dex lower at  ∼ 6 − 7. Mean  bol,AGN versus SFR predicted by the four simulations (lines: TNG300, TNG100, SIMBA100, and EAGLE100), compared with the observed 3 upper limits of  bol,AGN of LBGs (diamonds with down arrows: open symbols for total samples and filled ones for magnitude-divided samples).We calculate the  bol,AGN of simulated galaxies from their BHAR taking into account the value of  edd (see Section 5.3).We adopt the same color choice as in Fig. 6 to represent each simulation.
(iii) QSOs have ∼ 1 dex higher BHAR/SFR ratios than the  = 0  BH / ★ ratio, and hence ∼ 2 dex higher ratios than the LBGs, implying that the SMBHs in QSOs are growing much faster than expected from simultaneous evolution.We find a similar trend in BHAR versus HAR.
(iv) Effective (i.e., duty-cycle corrected) BHARs of  ∼ 4 QSOs suggest that the  BH increase in the QSO phase dominates over that in the much longer, inactive (BHAR ≲ 1 ⊙ yr −1 ) period.
(v) QSOs and SBGs have ∼ 1 dex higher SFR/ h and ≳ 2 dex higher BHAR/ h than LBGs, suggesting that if the duty-cycle of QSOs' and SBGs' phases is ≳ 10%, the stellar mass acquired in these phases dominates over that by normal star formation.
(vi) SBGs have more enhanced star-forming activity than QSOs at a fixed BHAR/ h .The BHAR/SFR of SBGs is comparable to the local  BH / ★ ratio, meaning that their SMBH and stellar components are growing at paces roughly consistent with simultaneous evolution.
(vii) QGs have comparable or higher BHAR/ h to LBGs despite having ∼ 1 dex lower SFR/ h .If these QGs have evolved from SBGs and QSOs, the increase in  BH in QGs appears to be insignificant unless the duty-cycle of SBGs and QSOs is ≲ 1%.
(ix) A comparison of our BHAR upper limits with four cosmological simulations (TNG300, TNG100, SIMBA100, and EAGLE100) finds that they overpredict the BHAR at a fixed SFR up to ∼1.5 dex.This result may imply that simulated SMBHs are too massive.Similar overpredictions are also found for  bol calculated on an assumption that low- edd objects effectively have low .
The large differences in BHAR/SFR and BHAR/HAR between LBGs and QSOs, as well as LBGs' much lower BHARs than predictions of cosmological simulations, suggest diverse co-evolution paths at  ≳ 4, the era when the BHAD no longer follows the SFRD.We plan to utilize deep observational data of our LBGs by large telescopes such as the JWST (e.g., COSMOS-Web), along with advanced simulations, to impose stronger constraints on co-evolution.

)
Yang et al. (2018) provide  bol ( X ) for log( X /erg s −1 ) > 42.4 sources, which is a modified version of Lusso et al. (2012)'s and is

Figure 2 .
Figure 2. Stacked soft-band X-ray images of our LBG samples.No signal is detected in any sample.The size of each image is 30 ′′ × 30 ′′ .
from the literature are also plotted in each panel.Small open symbols indicate individual QSOs; the  h of Trakhtenbrot et al. (2017); Netzer et al. (2014); Nguyen et al. (2020)'s QSOs are calculated in the same manner as Shimasaku & Izumi (2019) and those of Decarli et al. (2018)'s, Izumi et al. (2018)'s, and Venemans et al. (

Figure 4 .
Figure 4. BHAR plotted against HAR for  ∼ 4 (left panel),  ∼ 5 (middle), and  ∼ 6, 7 (right).The meaning of diamonds with down allows, small open symbols, black dots with error bars, and dashed lines is the same as in Fig 3.The  h of individual QSOs, which are needed to calculate HAR, are taken from (or calculated in the same manner as) Shimasaku & Izumi (2019).Shaded areas on the left and right panels correspond to QSO samples with clustering-based  h estimates (Timlin et al. 2018; He et al. 2018; Shen et al. 2007; Arita et al. 2023); for a given QSO sample, the range of the HAR corresponds to the uncertainty in  h estimate while the range of the BHAR corresponds to the central 68 percentile of the  1450 distribution, where we have calculated  bol from  1450 using eq.1 of Venemans et al. (2016).The median  1450 of these QSO samples are  1450 = −24.58(all) and  1450 = −24.39(faint) (Timlin et al. 2018),  1450 = −23.49(He et al. 2018),  1450 = −26.35(Shen et al. 2007), and  1450 = −23.92(Arita et al. 2023).Dot-dashed lines indicate the local  BH / h ratio (BHAR/HAR = 1.2 × 10 −4 ) (Kormendy & Ho 2013).

Figure 6 .
Figure 6.Top panels.MeanBHAR as a function of SFR predicted by TNG300 (brown dotted lines) and TNG100 (brown solid lines) with  = 0.2 for  ∼ 4 (left panel), 5 (middle), and 6.5 (right).Bottom panels.Same as the top panels but for SIMBA100 (magenta lines) and EAGLE100 (limegreen lines at  ∼ 4, 5, 6 and a green line at  ∼ 7) with  = 0.1.Since EAGLE100 does not provide a  = 6.5 snapshot, we plot  = 6 and 7 results.When calculating these functions, we only use SFR > 1 ⊙ yr −1 sources for TNG300, TNG100, and EAGLE, SFR > 10 ⊙ yr −1 sources for SIMBA100 to ensure that each simulation is complete in terms of SFR.We examine the completeness using the SFR versus  ★ plot.The meaning of diamonds with down arrows and dashed and dot-dashed lines is the same as in Fig.3.

Figure 7 .
Figure 7. Mean SFR as a function of  h predicted by the four simulations (TNG300, TNG100, SIMBA100, and EAGLE100) at  ∼ 4 (left panel),  ∼ 5 (middle), and  ∼ 6, 7 (right).Our LBGs are shown as open diamonds (total samples) and filled diamonds (magnitude divided samples).Red dashed lines represent the average relation of all galaxies at each redshift given in Behroozi et al. (2019).

Figure 8 .
Figure 8. BHAR versus HAR for He et al. (2018)'s (yellow shades) and Shen et al. (2007)'s (blue shades) QSO samples at  ∼ 4 before (light color) and after (dark color)  duty correction.The correction is made by multiplying the observed BHAR by  duty (0.03 for He et al. (2018) and 0.17 for Shen et al. (2007)).Also plotted are our LBGs (diamonds with down arrows) and the local  BH / h ratio (BHAR/HAR = 1.2 × 10 −4 ) of Kormendy & Ho (2013) (dot-dashed line).The meaning of a dashed line is the same as in Fig. 3.

Figure 9 .
Figure 9. Top panels.BHAR versus  h for LBGs (SFGs), QSOs (only those with SFR measurements), SBGs, and QGs at three redshift bins.The meaning of diamonds with down arrows, black dots with error bars, small open symbols, and dashed lines is the same as in Fig. 3.The newly added orange solid line represents the BHAR versus  h relation for all galaxies at  = 4 obtained by Yang et al. (2018).Open squares represent QGs at  ∼ 3.2 − 3.4 from Ito et al. (2022) and filled stars are for SBGs from Jin et al. (2018).We calculate the  h of QGs (SBGs) from their  ★ using the average  ★ −  h relation for QGs (SFGs) at  ∼ 4 given in Behroozi et al. (2019).Middle panels.SFR versus  h for LBGs, QSOs, SBGs, and QGs at three redshift bins.The meaning of diamonds and red dashed lines is the same as in Fig. 7. Bottom panels.Distributions of LBGs, QSOs, SBGs, and QGs in the BHAR/ h -SFR/ h plane at three redshift bins calculated from the top and middle panels.In each panel, the upper  axis shows SFR/HAR while the right  axis showing BHAR/HAR; we calculate HAR from  h for these two axes by approximating the term  h 10 12 ⊙ 1.1 X-ray to bolometric luminosity conversion factor calculated from Yang et al. (2018)'s Fig.8.

Figure 10 .
Figure 10.Comparison at  ∼ 4 of simulated mean  BH −  ★ relations (lines) with the 3 upper limits of  BH of our LBGs (diamonds with down arrows:open symbols are for total samples while filled symbols for magnitude divided samples) calculated on a simple assumption of ⟨BHAR⟩ simul /⟨BHAR⟩ obs = (  simul BH / true BH ) 2 for the four simulations (TNG300, TNG100, SIMBA100, and EAGLE100 from left).We adopt the same color choice as in Fig.6to represent each simulation.Dot-dashed lines indicate the local  BH - bulge relation(Kormendy & Ho 2013).In SIMBA100, the all, 25-26 mag, and 26-27 mag samples are not plotted because SIMBA100 is incomplete at SFRs corresponding to these faint apparent magnitudes.Similarly, in EAGLE100, the 23-24 mag sample is not plotted because EAGLE100 does not contain such bright (i.e., high-SFR) galaxies.

Figure
Figure11.Mean  bol,AGN versus SFR predicted by the four simulations (lines: TNG300, TNG100, SIMBA100, and EAGLE100), compared with the observed 3 upper limits of  bol,AGN of LBGs (diamonds with down arrows: open symbols for total samples and filled ones for magnitude-divided samples).We calculate the  bol,AGN of simulated galaxies from their BHAR taking into account the value of  edd (see Section 5.3).We adopt the same color choice as in Fig.6to represent each simulation.

Table 1 .
LBG samples and sub-samples used in this study.Number of LBGs outside the Chandra image.Number of LBGs within the Chandra image.Number of LBGs brighter than  UV = −23.Number of LBGs in each sub-sample.Number of LBGs in each sub-sample used for X-ray stacking.,,, The number in parentheses indicates the number of individual X-ray detected LBGs.