a regime potentially dominated by X-ray binaries (XRBs). In this paper, we aim at constraining the population of faint AGNs, and the relative contribution of the X-ray emission from the AGN and from the XRB population in galaxies of different masses in large-scale cosmological simulations.

XRBs are formed by a compact object, stellar mass BH or neutron star, accreting from a normal star. In a binary system, the mass is transferred from the star to the compact object via Roche lobe overflow or stellar wind mass transfer. The number and luminosity of XRBs are thought to primarily depend on the host galaxy properties, such as mass and specific star formation rate (sSFR; Lehmer et al. 2016Lehmer et al. , 2019;;Fornasini et al. 2018). The galaxy-wide emission from the binary population comes from both high-mass XRBs (HMXBs) and low-mass XRBs (LMXBs). The stellar companion in LMXBs is less massive than the compact object, and LMXB lifetime is longer than 1 Gyr (Lehmer et al. 2010). HMXBs are composed of a neutron star or stellar mass BH and a more massive stellar companion than in LMXBs. As such, their lifetime is shorter than 100 Myr. The galaxy-wide number of LMXBs increases with the galaxy stellar mass.

A scaling of the X-ray emission with the stellar mass of the galaxies and their star formation rate (SFR) was found in several analyses of nearby galaxies (D 50 Mpc; e.g. Grimm, Gilfanov & Sunyaev 2003;Colbert et al. 2004;Gilfanov 2004;Lehmer et al. 2010;Boroson, Kim & Fabbiano 2011;Mineo, Gilfanov & Sunyaev 2012a, b;Zhang, Gilfanov & Bogdán 2012). The most recent scaling relation was derived in Lehmer et al. (2019). In this work, galaxies with X-ray measurements from the Chandra data are decomposed into spatially resolved maps of SFR and stellar mass, and the XRB population statistics is extracted for sSFR bins. The parameters of the scaling relation functional form are fitted in these sSFR bins (i.e. fits over the pixels of all the galaxies together, and galaxy by galaxy). The result scaling relation is then tested by fitting each galaxy of the sample. These studies revealed that the total XRB emission of quiescent galaxies is dominated by LMXBs, for which the emission is proportional to the galaxy mass since LMXBs are long lived. However, the total XRB emission of star-forming galaxies is dominated by young short-lived luminous HMXBs (Fabbiano 2006), and the higher the galaxy SFR the higher the total XRB emission. The X-ray emission from LMXBs and HMXBs is comparable for sSFR ≈ 10 -10 yr -1 (Lehmer et al. 2019).

The X-ray emission of the XRB population could also depend on the metallicity and stellar ages of the galaxies (Fragos et al. 2013;Lehmer et al. 2016;Madau & Fragos 2017). Studies focusing on more distant and diverse galaxies have found deviations from the empirical relations established in the local Universe. This is the case for samples of low-metallicity galaxies or for samples with a wider range of stellar masses (e.g. Kim & Fabbiano 2010;Lehmer et al. 2014;Basu-Zych et al. 2013a, 2016;Douna et al. 2015;Brorby 2016;Tzanavaris et al. 2016). Finally, an evolution with redshift was found in X-ray stacking analyses (e.g. Lehmer et al. 2007Lehmer et al. , 2016;;Basu-Zych et al. 2013b;Kaaret 2014;Aird, Coil & Georgakakis 2017). The emission from LMXBs per unit stellar mass is thought to be higher at higher redshifts. The donor stars of the LMXBs have higher masses at higher redshifts and the LMXBs are more luminous (Lehmer et al. 2016). Moreover, elliptical galaxies with many globular clusters show an excess of LMXBs. Globular clusters can produce LMXBs efficiently through stellar dynamical interactions (Cheng et al. 2018;Lehmer et al. 2019), i.e. by either capture of a stellar object by the second object, or hardening of soft binaries (hard binaries are binaries with bound energy larger than the kinetic energy of the intruding star and stellar encounters involving hard binaries make them harder, see Cheng et al. 2018). The HMXB emission per unit SFR is higher at higher redshifts. This is connected to the lower metallicities of the stellar populations, as observationally shown in Fornasini et al. (2019) and Fornasini, Civano & Suh (2020). The star formation at low metallicities is expected to result in a larger number of massive compact objects and Roche lobe overflow binaries because the mass-loss of massive stars through stellar winds is less effective (Dray 2006;Linden et al. 2010;Fragos et al. 2013;Lehmer et al. 2016).

There is evidence for the presence of X-ray AGN in all types of galaxies, including dwarf galaxies (M 10 9.5 M ; Mezcua et al. 2016;Chilingarian et al. 2018), where the galaxy-wide X-ray emission is higher than expected from the XRB population. In Mezcua & Domínguez Sánchez (2020), the authors identify AGN signatures in 37 local dwarf galaxies in the range log 10 M /M = 8.8 -9.5 using integral field unit (IFU) spectroscopy. These AGNs are faint, sometimes with off-centre X-ray emission, and have bolometric luminosity in the range log 10 L bol /(erg/s) = 38.9 -41.4. BH mass estimates for these AGNs could lie in the range log 10 M BH /M = 5.5 -8. The presence of AGNs in dwarf galaxies could also shed new light on the role of AGN feedback in these galaxies, as studied in simulations (Koudmani et al. 2019;Koudmani, Henden & Sijacki 2021), and observed with the detections of fast AGN-driven outflows (Manzano-King, Canalizo & Sales 2019; Liu et al. 2020).

Signatures of faint AGNs have been identified in more massive galaxies than dwarf galaxies, particularly with stacking analysis of Chandra COSMOS Legacy galaxies (Georgantopoulos et al. 2017;Fornasini et al. 2018). The BHs powering these faint AGNs can have low accretion rates for different reasons: for example, BH accretion rates regulated by BH feedback, BH accretion rates regulated by the feedback in the environment of the BHs, BHs in quenched galaxies depleted in gas, off-centre BHs that do not benefit from the galaxy potential well gas reservoir, and the accretion rate likely depends on BH mass and galaxy properties. Constraining the regime of faint AGNs will help to understand the co-evolution and interplay between BHs and their host galaxies. In this study, we will investigate in which galaxies faint AGNs reside in several cosmological simulations, and whether the low accretion rates on to the simulated faint AGNs can have different causes (e.g. impact of AGN feedback, absence of cold gas), such as outlined above.

The galaxy-wide emission of the XRB population could depend on the SFR of the host galaxies, and in a similar way correlations between AGN activity and the SFR of their galaxies could exist. Even for the regime of the bright AGN (e.g. L AGN 10 42 erg s -1 ), the correlations between AGN activity and the properties of galaxies are still not clear, especially in terms of stellar masses, SFR, and sSFR (e.g. Aird, Coil & Georgakakis 2019). For bright AGN of L AGN 10 44 erg s -1 in massive galaxies, there is evidence for a correlation between AGN luminosity and SFR (Lutz et al. 2008;Bonfield et al. 2011;Mor et al. 2012;Rosario et al. 2012). There are very few observational constraints in the faint AGN regime, and therefore, in order to theoretically investigate possible dependences/correlations of our results with galaxy properties, we will split the simulated galaxies in starburst, star-forming, or quiescent galaxies.

We are also particularly interested in how we can compare the predictions from cosmological simulations to current and future observations to disentangle simulation subgrid models and improve them. In this work, we compare the total X-ray luminosity of the simulated galaxies hosting faint AGNs to observational constraints of stacked galaxies (Fornasini et al. 2018). The high sensitivity of the upcoming X-ray mission Athena (Nandra et al. 2013) and the NASA concept missions Lynx (Gaskin et al. 2018) and AXIS (Mushotzky et al. 2019) will allow us to investigate the properties of AGNs up to high redshift, and will provide us with new constraints on the faint AGN population. Our work paves the way for future investigations on faint AGNs, a goal aligned with the upcoming X-ray missions. As in all studies on the AGN population, obscuration is critical and a major open question. The attenuation of radiation from intervening gas and dust can affect the contribution of both the AGN and the XRBs to the observed galaxy X-ray luminosity. We will build several models for the obscuration of the faint AGN to account for this.

In this paper, we use the four Illustris, TNG100, EAGLE, and SIMBA cosmological hydrodynamical simulations of 100 comoving Mpc (cMpc) box side length. For all these simulations, we model both the AGN luminosity and the XRB population of the simulated galaxies using several empirical scaling relations. We describe the physical models of the simulations, as well as our post-processing modelling of AGN and XRB luminosity in Section 3. In Section 4, we analyse the AGN population of the simulations and in which galaxies they live. In order to derive the total X-ray luminosity of the simulated galaxies, we first analyse in Section 5 the properties of their AGN as a function of their host galaxy's mass and SFR. In Sections 6 and 7, we derive the relative contribution of the AGN and the XRB population to the galaxy total hard X-ray luminosity and we compare these results to recent findings in observations of stacked galaxies from the Chandra COSMOS Legacy survey (Fornasini et al. 2018).

In this section, we describe the Illustris, TNG100, EAGLE, and SIMBA cosmological simulations, and our method to compute AGN luminosity and the X-ray luminosity of the XRB population of the simulated galaxies. We also describe our method to divide the simulated galaxy in three samples with different sSFR.

In the following, we briefly summarize the BH subgrid models, i.e. seeding, accretion, and feedback of the Illustris, TNG100, EAGLE, and SIMBA large-scale cosmological simulations. A more detailed comparison of the modelling of these simulations is presented in Habouzit et al. (2021). Several other aspects differentiate these simulations, for example the presence of magnetic fields in the TNG100 simulation (Pillepich et al. 2018b), the single mode AGN feedback of EAGLE (Schaye et al. 2015), the two mode accretion model of SIMBA (Davé et al. 2019). The detailed models can be found in the references given below.

The Illustris hydrodynamical simulation (Genel et al. 2014;Vogelsberger et al. 2014a, b;Sijacki et al. 2015) simultaneously follows the evolution of dark matter (DM) and baryonic matter, in a volume of (106.5 cMpc) 3 . The Illustris simulation data are publicly available (Nelson et al. 2015). The simulation produces a mix of galaxy morphologies in broad agreement with observations (Vogelsberger et al. 2014a). The equations of gravity and hydrodynamics are evolved with the moving-mesh code AREPO (Springel 2010). The galaxy formation model includes gas cooling, star formation, supernova (SN) feedback, and the physics of BHs. BHs are seeded in haloes with halo mass 7.1 × 10 10 M with seed mass M seed = 1.4 × 10 5 M . BHs can grow in mass with gas accretion and BH mergers. The BH accretion is modelled with the Bondi-Hoyle-Lyttleton formalism as

(1)

where ρ and c s are the density and the sound speed of the surrounding gas, respectively, v BH is the velocity of the BH relative to the gas, and G is the gravitational constant. The boost factor α = 100 accounts for the unresolved gas around the AGN, which tend to underestimate the accretion rates. The boost factor is set to produce a BH population that agrees with the M BH -M diagram at z = 0. A repositioning scheme for BH sink particles is used which ties them to the local minimum gravitational potential. M Edd is the Eddington accretion rate of the BH:

where m p is the proton mass, σ T is the Thomson cross-section, c is the speed of light, and r is the radiative efficiency which is set to 0.2 in the simulation. The Eddington ratio is defined as f Edd = ṀBH / ṀEdd . The feedback from the AGN depends on the BH accretion rate and can operate in high-accretion mode (f Edd > 0.05) or low-accretion mode (f Edd < 0.05). The low-accretion mode model injects highly bursty thermal energy into large 'bubbles' (∼50 kpc) which are displaced away from the central galaxy. In the high-accretion mode, thermal energy is continuously injected into the surrounding gas. Moreover, radiative AGN feedback that modifies the ionization state of the surrounding gas is included as well (Vogelsberger et al. 2013).

The TNG100 simulation builds directly on the Illustris framework and has a volume of (110.7 cMpc) 3 (Marinacci et al. 2018;Naiman et al. 2018;Nelson et al. 2018Nelson et al. , 2019b;;Springel et al. 2018;Pillepich et al. 2018a). The simulation includes magnetic fields. The BH seed mass is M seed = 1.1 × 10 6 M . The BH seed mass was increased in comparison with the Illustris seed mass by nearly one order of magnitude (Pillepich et al. 2018b). The BH accretion rate is described by the Bondi-Hoyle-Lyttleton accretion rate limited by the Eddington accretion rate (Weinberger et al. 2017). However, in the TNG model, the additional boost factor α from equation ( 1) is removed. The effective sound speed c 2 s = c 2 s,therm + (B 2 /4πρ) includes the thermal and the magnetic signal propagation. The transition between the AGN feedback modes in the TNG model depends on the BH accretion rate and BH mass. A BH is assumed to be in the high accretion state as long as (Weinberger et al. 2017)

The low accretion feedback is modelled as kinetic outflows from the BHs (Weinberger et al. 2017), while the high accretion mode is modelled as injection of thermal energy in the BH surroundings (Weinberger et al. 2017).

The EAGLE simulation (Crain et al. 2015;Schaye et al. 2015) has a volume of (100 cMpc) 3 , and was run with the code ANARCHY. This modified version of GADGET3 (Springel 2005) is based on the smoothed particle hydrodynamics (SPH) method. The simulation includes subgrid models for radiative cooling, star formation, stellar mass-loss, metal enrichment, and energy feedback from star formation (Schaye et al. 2015). The BH seeding mass is M seed = 1.48 × 10 5 M and BHs are seeded in haloes with a mass of at least M h = 1.48 × 10 10 M . The accretion rate is determined by the Bondi-Hoyle-Lyttleton model (Rosas-Guevara et al. 2015, 2016):

with ṀBondi and ṀEdd defined as in equations ( 2) and ( 3) and ṀBH = (1r ) Ṁacc . The radiative efficiency is r = 0.1. Here, c s represents the sound speed of the surrounding gas, V is the rotation speed of the gas around the BH, and C visc is a free parameter. It is related to the viscosity of the accretion disc. The correction factor min ((c s /V ) 3 /C visc , 1) multiplied with the Bondi rate accounts for a lower accretion rate for gas with angular momentum. In that case, the accretion is not spherically symmetric and proceeds through an accretion disc (Schaye et al. 2015). The simulation uses a singlemode AGN feedback model. A fraction of the accreted gas on to the BHs is stochastically injected in the surroundings as thermal energy (Schaye et al. 2015).

SIMBA (Davé et al. 2019) has a volume of (147 cMpc) 3 and builds on its predecessor Mufasa (Davé, Thompson & Hopkins 2016), using the code GIZMO (Hopkins 2015(Hopkins , 2017) ) in its 'meshless finite mass' hydrodynamics mode. The BH seeding mass is M seed = 1.4 × 10 4 M and BHs are seeded in galaxies with M ∼ 10 9.5 M . This is motivated by FIRE simulations showing that stellar feedback strongly suppresses BH growth in lower mass galaxies (Anglés-Alcázar et al. 2017b;C ¸atmabacak et al. 2020), in agreement with other models (e.g. Habouzit, Volonteri & Dubois 2017). Because of the seeding, we will not draw any conclusion on BH properties in galaxies with M < 10 9.5 M for the SIMBA simulation. The BH accretion rate is given by ṀBH = (1r ) min( ṀBondi , ṀEdd ) + min( ṀTorque , 3 ṀEdd ) , (6) with ṀBondi and ṀEdd defined as in equations ( 2) and (3). The radiative efficiency is r = 0.1. ṀTorque describes the gas inflow rate driven by gravitational instabilities from the scale of the galaxy to the accretion disc of the BH (Hopkins & Quataert 2011;Anglés-Alcázar et al. 2017a). This accretion is only evaluated for the cold gas (T < 10 5 K). For the hot gas (T > 10 5 K), the Bondi-Hoyle-Lyttleton model is applied (Davé et al. 2019). AGN feedback in SIMBA is modelled as an injection of kinetic energy following a two-mode approach, with high-mass loading outflows in the radiative 'quasar' mode and lower mass loading but faster outflows at low Eddington ratios in the jet mode (Davé et al. 2019). BHs begin to transition into jet mode for f Edd < 0.2, reaching full speed jets at f Edd ࣠ 0.02. Also, X-ray feedback is included for galaxies with M gas < 0.2 M and full speed jets, following the implementation of Choi et al. (2012). See Thomas et al. (2019Thomas et al. ( , 2021) ) and Habouzit et al. (2021) for previous studies of BHs in SIMBA.

We compute in post-processing the luminosity of the BHs, following the model of Churazov et al. (2005), i.e. explicitly distinguishing radiatively efficient and radiatively inefficient AGNs. The bolometric luminosity of radiatively efficient BHs (f Edd > 0.1) is given by

BHs with small Eddington ratios (f Edd ≤ 0.1) are considered radiatively inefficient and their bolometric luminosities are given by (Habouzit et al. 2021)

We compute the hard (2-10 keV) X-ray luminosities of the AGN with the bolometric correction of Hopkins

In the following, we present the luminosity of the AGN in erg s -1 .

AGNs can be obscured by gas and/or dust, either physically close to the AGN and at larger scales within the host galaxies. We generally refer to two types of obscured AGNs, either Compton-thin AGN with column densities of N H 10 22 cm -2 , or more heavily obscured Compton-thick AGN with N H 10 24 cm -2 . While the hard X-ray band suffers less from obscuration than softer bands, the observed luminosity of the AGN could still be attenuated. The obscuration might depend on the intrinsic AGN luminosity and the redshift: there is more gas present at higher redshifts and therefore, the obscuration could be stronger (Gilli et al. 2014;Merloni et al. 2014;Vito et al. 2014Vito et al. , 2016)). We build several models for the fraction of obscured AGNs and apply them to the simulations. Our models depend on redshift and hard X-ray (2-10 keV) AGN luminosity for AGNs with L AGN 10 41 erg s -1 (Habouzit et al. 2019, for more details), and only on AGN luminosity for fainter AGNs. Obscuration could in principle also depend on galaxy SFR, a parameter that we do not consider here.

For the AGN with L x 10 41 erg s -1 , we build the models based on the observational constraints of Merloni et al. (2014) and Ueda et al. (2014): the coloured lines of Fig. 1 are a fit to the observations. These constraints assume that the obscured AGNs mostly include Compton-thin AGNs, but that the presence of Compton-thick AGNs cannot be ruled out. Here, we assume that our models based on these observations include all types of obscured AGNs, and we apply Table 1. Parametrizations from the literature of the hard X-ray luminosity of the XRB population in galaxies, as a function of their total stellar mass M (M ), SFR, and redshift z. In this paper, we mostly use the scaling relation from Lehmer et al. (2019). The M , SFR, and z ranges describe the sample of galaxies used to derive the empirical relations.

Scaling relations M ranges SFR ranges (M yr -1 ) Redshift ranges L10 Lehmer et al. (2010) L XRB = 10 28.96 M + 10 39.21 SFR M = 10 10-11 M SFR = 10 1 -2 z ∼ 0 L16 Lehmer et al. (2016) L XRB = 10 29.37 (1 + z) 2.0 M + 10 39.28 (1 + z) 1.3 SFR M = 10 9-12 M SFR = 10 -2 -3 z = 0-4 L19 Lehmer et al. (2019) L XRB = 10 29.15 (1 + z) 2.0 M + 10 39.73 (1 + z) 1.3 SFR M = 10 9-11.5 M SFR = 10 -2 -1 z ∼ 0 A17 Aird et al. (2017) L XRB = 10 28.81 (1 + z) 3.9 M + 10 39.50 (1 + z) 0.7 SFR 0.86 M = 10 8.5-11.5 M SFR = 10 -1 -3 z = 0.1-4 F18 Fornasini et al. (2018) L XRB = 10 29.98 (1 + z) 0.62 M + 10 39.78 (1 + z) 0.2 SFR 0.84 M = 10 9.5-11.5 M SFR = 10 -1 -3 z = 0.1-5 them on the full AGN simulated samples. Since we are lacking observational constraints on the population of faint AGNs, and even more so on the fraction of those that are obscured, we build four models that cover a broad range of possible fractions of obscured AGNs (Fig. 1). In our model M1, a large fraction (80 per cent) of the faint AGN with L AGN 10 41 erg s -1 are obscured, and this fraction does not depend on AGN luminosity. For M2, M3, and M4, the fraction does depend on L AGN and can progressively go down to 60 per cent (M2), 40 per cent (M3), or 20 per cent (M4) of AGNs being obscured. Our choice of low fractions of obscured AGNs among the faint AGN with L AGN 10 40 erg s -1 is motivated by the fact that in simulations if these BHs do not have high accretion rates, their close surrounding is depleted of gas. Consequently, if the gas reservoir is small there may be little room for the gas/dust to obscure the emission from the AGN. Since this is highly uncertain, we also have our model M1 with a high fraction of obscured AGNs to bracket the range of possibilities. In the following, we apply our obscuration models only when specified in the text and figures.

We parametrize the X-ray emission (2-10 keV band) from the XRB population of the simulations (which is not modelled in the simulations), including both LMXBs and HMXBs. To do so, we employ XRB emission scaling relations that have been derived from observational samples. These relations provide the luminosity of the XRB population for a given galaxy as a function of its SFR, its stellar mass, and for some of them its redshift:

For example, the parameters log 10 α LMXB = 29.37 ± 0.15 erg/s/M , log 10 β HMXB = 39.28 ± 0.05 erg/s/(M /yr), γ = 2.03 ± 0.60, δ = 1.31 ± 0.13 have been described, e.g. in Lehmer et al. (2016Lehmer et al. ( , 2019)). The relations described in Aird et al. (2017) and Fornasini et al. (2018) use a power law of SFR 0.86 and SFR 0.84 , respectively. While it has been shown that the metallicity of the galaxies can play a role in the amplitude of the XRB emission, especially in some specific regimes of stellar mass and redshift, we do not include any metallicity dependence in our analysis below.

To study the impact of the modelling of the XRB luminosity, we use five different empirical relations (Lehmer et al. 2010(Lehmer et al. , 2016(Lehmer et al. , 2019;;Aird et al. 2017;Fornasini et al. 2018). We report these relations in Table 1. We compare the scaling relations in Fig. 2, as a function of galaxy stellar mass (top panels), SFR (middle panels), and sSFR (bottom panels). The left-hand panels show redshift z = 0, and the right-hand panels z = 2. The difference between the different models can be up to one order of magnitude in luminosity at z = 0, and more than an order of magnitude at higher redshift. At fixed SFR, the luminosity of the XRB population increases with the stellar mass of galaxies. We note that for galaxies with high SFR of log 10 SFR/(M /yr) ∼ 2 the luminosity is almost constant (Lehmer et al. 2010(Lehmer et al. , 2016(Lehmer et al. , 2019;;Aird et al. 2017;Fornasini et al. 2018), at z = 0 (lefthand panels) and z = 2 (right-hand panels). The difference between the models can be up to one order of magnitude in hard X-ray (2-10 keV) luminosity at z = 0, and more at higher redshift. The XRB luminosity increases with M (at fixed SFR), and with SFR (at fixed M ), and with redshift for most of the models. LMXBs dominate the hard X-ray (2-10 keV) luminosity of the galaxies with log 10 sSFR/yr -10.5. For higher sSFR, the X-ray emission is dominated by HMXBs.

with stellar mass, for all the scaling relations studied here. The XRB population luminosity also increases with SFR, at fixed galaxy stellar mass. The model of Lehmer et al. (2010) does not evolve with redshift. For the other models, the normalization of the scaling relations increases with increasing redshift. In general, the models of Figure 3. SFR -M plane with hexabins colour coded by the number of galaxies in the bins. Galaxies with SFR 10 -4 M yr -foot_1 are shown with this value. We define three samples of galaxies: the high-sSFR sample with galaxies 0.5 dex above the star-forming main sequence, the intermediate-sSFR sample with galaxies on the star-forming sequence, and the low-sSFR sample with galaxies below 0.5 dex of the main sequence. The high-sSFR, intermediate-sSFR, and low-sSFR samples are, respectively, composed of starburst, main-sequence, and quiescent galaxies. Lehmer et al. (2019) and Fornasini et al. (2018) provide the highest normalizations of the L XRB relation at z = 0. At higher redshift, the model of Aird et al. (2017) also provides high luminosities. For low sSFR galaxies with log 10 sSFR/yr < -10.5, the XRB population luminosity is dominated by LMXBs whose luminosity depends on the stellar mass. But for galaxies with higher sSFR of log 10 sSFR/yr > -10.5, the XRB luminosity is dominated by HMXBs whose luminosity scales with the SFR (Lehmer et al. 2019).

In this paper, we will present our results as a function of galaxy properties, i.e. their stellar mass and SFR. In the following, we define the simulation star-forming main sequence and build three different galaxy samples from the simulations.

We show the SFR -M relation for galaxies from the Illustris, TNG100, EAGLE, and SIMBA simulations in Fig. 3. We define the star-forming main sequence of the simulations as the mean SFR in fixed width bins (same width for all the simulations) of galaxy stellar mass for 10 9 M 10 10 M 1 (see also Genel et al. 2014;Vogelsberger et al. 2014b;Furlong et al. 2015;Sparre et al. 2015;Bluck et al. 2016;Terrazas et al. 2017;Davé et al. 2019;Donnari et al. 2019;Matthee & Schaye 2019). We purposely exclude more massive galaxies to compute the main sequence as many of these galaxies are quenched. The main sequence can be defined as a power law:

with SFR MS the SFR of the galaxies on the main sequence, in M yr -1 . The parameters α and β are computed for each simulation and redshift and are given in Table A1. The simulations present some differences in the galaxy population in Fig. 3, such as the exact normalization and slope of their star-forming main sequence, or the formation of the population of quenched galaxies, i.e. galaxies with low or null SFR. We note that the SIMBA simulation has the steepest main sequence, but overall we find a good agreement between the simulations for the normalization of the star-forming main sequence.

In observations there is, at the moment, no real consensus on the exact slope and normalization of the star-forming main sequence, and the parameters depend on the observed samples. The slope of the star-forming main sequence of observed galaxies can be slightly shallower (Habouzit et al. 2019;Hahn et al. 2019). When looking at the population of quenched galaxies, we can identify some features of the AGN feedback modelling such as the sharp decrease of the SFR in TNG100 galaxies of M a few 10 10 M (but still in broad agreement with observational data, see Donnari et al. 2021b). This shows the transition at M BH ∼ a few 10 8 M between the high accretion state and the low accretion state of the AGN feedback modelling (see equation 4). The low accretion mode feedback is efficient, and both regulates the BHs and quenches their host galaxies (Weinberger et al. 2018;Habouzit et al. 2019). Towards z = 0, we also note a strong decrease of the SFR of a large fraction of the galaxies with M 10 9.5 M in all the simulations, which reflects the quenching of satellite galaxies (Donnari et al. 2021a).

For the purpose of our analysis, we divide the galaxies into three subsets, which are simulation-and redshift-dependent. The high-sSFR subset (blue background in Fig. 3) consists of starburst galaxies with SFR higher than half a dex above the main sequence. The intermediate-sSFR subset (green background in Fig. 3) includes starforming galaxies on the main sequence, i.e. within 1 dex. Finally, the low-sSFR subset (red background in Fig. 3) includes all galaxies with SFR below half a dex of the main sequence, i.e. quiescent galaxies or galaxies on their way to quiescence.

We find that a low fraction of the galaxies of M 10 9 M are in the high-sSFR sample, typically less than 5 per cent in all the simulations and at all redshifts. The intermediate-sSFR represents the largest fraction of galaxies ( 90 per cent, z = 4) at high redshift in all the simulations. With time, the number of galaxies in the intermediate-sSFR samples decreases, for all the simulations. The percentage of galaxies in the low-sSFR samples increases with time for all the simulations, as more and more galaxies quench. Most of the galaxies in the low-sSFR samples are massive with M 10 10 M and have reduced SFR because of AGN feedback. For example, in TNG there is a sharp decrease of the SFR in galaxies of M 10 10.5 M due to the kinetic low-accretion mode of the AGN feedback model. Still, some of the galaxies present in the low-sSFR samples have lower masses (M < 10 10 M ), and can have lower SFR due to gas starvation, SN feedback, and environmental quenching. The final percentage of galaxies in each sample at z = 0 varies from simulation to simulation. We present these numbers in Table A2. We use these three samples in Section 5.

From now on, we only consider galaxies of total stellar mass M 10 9 M which host a BH. The BH can be an AGN, i.e. accreting mass, or a non-accreting BH, and in that case L AGN = 0. The mass of the BHs is the mass of individual BHs in all the simulations, except in TNG and Illustris for which the BH mass is the sum of the mass of all the BHs within a galaxy. In practice, only a couple of galaxies per output host several BHs at the same time. The total stellar mass that we use here for all the simulations is not exactly computed in the same way in all the simulations, but we prefer to use each simulation definition. The total stellar mass of the galaxies is computed, e.g. as twice the stellar mass in the half mass radius for the Illustris and TNG100 simulation, within an aperture of 30 kpc in EAGLE. In SIMBA, galaxies are identified using a friends-of-friends galaxy finder, assuming a spatial linking length of 0.0056 times the mean interparticle spacing (equivalent to twice the minimum softening length). Discussion on the impact of these different definitions can be found in Pillepich et al. (2018a) for TNG. Finally, we do not distinguish between central and satellite galaxies, and all our results group these two types of galaxies.

In this section, we first investigate the properties of galaxies hosting faint AGNs. We divide the galaxies in subsamples depending on their AGN luminosity: L AGN = 10 37.5-38.5 erg s -1 , L AGN = 10 39.5-40.5 erg s -1 , and L AGN = 10 41.5-42.5 erg s -1 . For each subsample as well as for the whole galaxy sample, we show the distributions of stellar mass, BH mass, and sSFR in Fig. 4. Galaxies with log 10 sSFR/yr -1 < -13 are shown with this value.

To quantify the findings described in the following, we perform several Kolmogorov-Smirnov (KS) tests for the distributions of the faint AGN with log 10 L AGN /(erg s -1 ) ∼ 38 and log 10 L AGN /(erg s -1 ) ∼ 42, for all the simulations and for M , M BH , and SFR (Table B2). We cannot reject the null hypothesis that two distributions were drawn from the same distribution for high values of p, i.e. p ≥ 0.01. However, if p < 0.01 the two distributions can be considered as statistically different.

In the Illustris simulation, faint AGNs with log 10 L AGN /(erg s -1 ) ∼ 38 are mostly located in low-mass galaxies of M 10 10 M at any redshift. These AGNs also tend to be low-mass BHs of M BH 10 8 M . More luminous AGNs with L AGN ∼ 10 42 erg s -1 are also located in low-mass galaxies at high redshift, but with time they start deviating from the L AGN ∼ 10 38 erg s -1 AGNs and are largely located in massive galaxies of M ∼ 10 11 M at z = 0. These AGNs are also mostly massive BHs of M BH ∼ 10 9 M at z = 0.

For Illustris, we find with our KS test that the distributions of M , M BH , and sSFR for log 10 L AGN /(erg s -1 ) ∼ 38 and log 10 L AGN /(erg s -1 ) ∼ 42 are always statistically different, except the sSFR distribution at z = 3.

In the TNG100 simulation, the faint AGNs are located in galaxies with generally different properties. From high redshift to low redshift, the AGNs with log 10 L AGN /(erg s -1 ) ∼ 38 are mostly located in massive galaxies of M 10 10.5 M . There is also a population of AGN with log 10 L AGN /(erg s -1 ) ∼ 38 in lower mass galaxies of M ∼ 10 9 M , similar to the population found in Illustris. Brighter AGNs with log 10 L AGN /(erg s -1 ) ∼ 42 are mostly hosted in lower mass galaxies with M 10 10 M . Since the M BH -M relationship is tight in TNG100 (Habouzit et al. 2021, and Fig. 5 below), the log 10 L AGN /(erg s -1 ) ∼ 38 AGNs are also mostly massive BHs of M BH ∼ 10 9 M . More luminous AGNs are BHs of M BH 10 7 M at high redshift (z ≥ 2) and BHs of M BH 10 8 M at z ≤ 2. The distributions of the log 10 L AGN /(erg s -1 ) ∼ 38 and the log 10 L AGN /(erg s -1 ) ∼ 42 AGN are statistically different (KS test). At z = 3, there are only few AGNs with log 10 L AGN /(erg s -1 ) ∼ 42 and the KS test cannot be performed.

The fact that faint AGNs (log 10 L AGN /(erg s -1 ) < 40) are found in massive galaxies of a few 10 10 M and that they are powered by BHs of a few 10 8 M is due to the efficient low-accretion mode of the TNG feedback model. These AGNs are found in quenched galaxies with reduced SFR. The population of quenched galaxies can also be observed in the sSFR distribution. Normalized distributions of galaxy stellar mass, BH mass, and sSFR, for z = 0, 1, and 3 (top, middle, and bottom panels). Red distributions only include galaxies with L AGN ∼ 10 38 erg s -1 (bin 10 37.5 -10 38.5 erg s -1 ), L AGN ∼ 10 40 erg s -1 (bin 10 39.5 -10 40.5 erg s -1 ) for the green distributions, and L AGN =∼ 10 42 erg s -1 (bin 10 41.5 -10 42.5 erg s -1 ) for the blue ones. Distributions are normalized by the total number of galaxies in the specific luminosity bin. The grey shaded histograms show the distributions for the full galaxy population. Faint AGNs of different luminosities are powered by different BHs and reside in different stellar mass galaxies in the different simulations. For example, the AGNs with L AGN ∼ 10 38 erg s -1 reside in low-mass galaxies in Illustris and in galaxies with much higher masses in TNG100 and SIMBA. 2004) (lower normalization of the relation). In these simulations, the median BH mass of low-sSFR galaxies is higher than for higher sSFR galaxies, at fixed stellar mass.

In the EAGLE simulation, the AGNs with log 10 L AGN /(erg s -1 ) ∼ 42 are mostly located in galaxies of M ∼ 10 10 M at z = 3. The distribution of the stellar mass of their host galaxies becomes broader with time, with AGNs located in the range M = 10 9.5 -10 11 M . At high redshift z ≥ 2, the distribution of the stellar mass of the galaxies hosting AGNs with log 10 L AGN /(erg s -1 ) ∼ 38 is different from the log 10 L AGN /(erg s -1 ) ∼ 42 AGNs. The log 10 L AGN /(erg s -1 ) ∼ 38 AGNs are located in less massive galaxies. However, at lower redshift (z < 1) the M distribution of log 10 L AGN /(erg s -1 ) ∼ 38 extends to more massive galaxies and become very similar to the distribution of more luminous AGNs with log 10 L AGN /(erg s -1 ) ∼ 42. For the mass of the BHs powering these AGNs, we note some differences at high redshift (z ∼ 3) between AGNs of log 10 L AGN /(erg s -1 ) ∼ 38 and log 10 L AGN /(erg s -1 ) ∼ 42, with fainter AGNs being lower mass BHs of M BH ∼ 10 6 M . The M BH distributions have broad humps for the two types of AGNs at high redshift. At z = 0 when the distributions are broader, the peak of the distributions is still different but the distributions are less distinguishable.

In EAGLE, the faintest AGNs are found in two BH populations: predominantly in low-mass BHs of M BH 10 6 M , but also in more massive BHs of M BH ∼ 10 7 M . The first BH population does not power luminous AGNs because they reside in low-mass galaxies that are regulated by SN feedback. The second BH population consists of more massive BHs located in more massive galaxies, and are regulated by AGN feedback.

For the EAGLE simulation, by applying a KS test we find that all the log 10 L AGN /(erg s -1 ) ∼ 38 and log 10 L AGN /(erg s -1 ) ∼ 42 distributions are significantly different, except the sSFR distributions at z = 0 and z = 3 and the M distribution at z = 0.

The distributions of the faint AGN with log 10 L AGN /(erg/s) ∼ 38 and log 10 L AGN /(erg/s) ∼ 42 of the SIMBA simulation are qualitatively similar to the distribution of TNG100. However, the M distributions of TNG100 extend to lower values than the M distributions of SIMBA. Our KS test indicates that the distributions of log 10 L AGN /(erg/s) ∼ 38 and log 10 L AGN /(erg/s) ∼ 42 in SIMBA are statistically different. At z = 3, there are only few AGNs with log 10 L AGN /(erg/s) ∼ 42 and we cannot perform the KS test.

AGNs with log 10 L AGN /(erg/s) ∼ 38 are located in more massive galaxies than the AGNs with log 10 L AGN /(erg/s) ∼ 42. The peak of the two M corresponding distributions are separated by less than an order of magnitude. The M distributions are less separated at higher redshift. For the BHs, we note that at z = 1, the peak of the BH mass is very similar, but separate at later time. At z = 0, simulated AGNs with log 10 L AGN /(erg/s) ∼ 42 are found mostly as BHs of M BH ∼ 10 7 M . In contrast, the distribution of log 10 L AGN /(erg s -1 ) ∼ 38 AGN peaks at M BH ∼ 10 8 M .

The sSFR distributions are very different for the different AGN luminosities in SIMBA. The AGNs with log 10 L AGN /(erg s -1 ) ∼ 38 have a peak in their sSFR distribution which is nearly two orders of magnitude smaller than the AGNs with log 10 L AGN /(erg s -1 ) ∼ 42. In SIMBA, the faintest AGNs are in galaxies with a low sSFR, they are quenched, and the AGNs with log 10 L AGN /(erg s -1 ) ∼ 42 are in star-forming galaxies. There is a large population of galaxies that quench with time, and with them a large population of fainter AGNs at z = 0.

In this first section, we demonstrated that the BHs powering the faint AGN in the Illustris, TNG100, EAGLE, and SIMBA simulations, as well as the properties of their host galaxies, are different among the different simulations. Indeed, faint AGNs of L AGN ∼ 10 38 erg s -1 can be powered by relatively massive BHs and be located in massive galaxies (M 10 10 M ) with reduced SFR (TNG100, SIMBA), or be powered by lower mass BHs in less massive galaxies (M 10 10 M ) still forming stars (Illustris, EAGLE). Both galaxy and BH mass appear to be fundamental quantities to understand the faint AGN populations in simulations. Moreover, in some simulations we already see that the sSFR of the host galaxies can play an important role. In the next section, we go further and divide the full simulated galaxy samples into subsamples as a function of the galaxies distance from the star-forming main sequence.

We show the median of the mass of the BHs located in the three galaxy samples in Fig. 5, as a function of the total stellar mass of the galaxies and for several redshifts from z = 0 (top row) to z = 4 (bottom row).

On average we find that galaxies with lower SFR (low-sSFR samples, red lines) host more massive BHs at fixed stellar mass; the effect is mild and in most cases within the 15th-85th percentiles of the sample distributions. We find this for the Illustris and TNG100 simulations at all stellar masses, and for galaxies of M 10 10 M in the EAGLE and SIMBA simulations. At low redshift, the EAGLE high-sSFR galaxies with M < 10 10 M (blue lines) host BHs slightly more massive than galaxies with lower SFR. In SIMBA, we cannot really assess the behaviour of BH mass of galaxies with M 10 10 M as the seeding takes place in galaxies of M ∼ 10 9.5 M . Most simulations have a tight M BH -M relation (not shown here, but see Habouzit et al. 2021), TNG100 having the tightest relation, which explains the very mild difference in the median of BH mass for the three samples at all redshifts. We do not show it here, but the M BH -M diagram of the simulations at z = 0 is in broad agreement with the observational samples of e.g. Reines & Volonteri (2015). To guide the eye, we show in grey shaded area in the z = 0 panels a region enclosing several M BH -M bulge empirical scaling relations (e.g. Häring & Rix 2004;Kormendy & Ho 2013;McConnell & Ma 2013) that have been used to calibrate the simulations. It should be noted that the calibrations are done on the entire galaxy population of the simulations and not on one given sSFR subsample, whereas in this paper we are testing the simulation outcomes to deeper and more constraining details. It may hence be not surprising that some discrepancies as the ones revealed above are in place between e.g. the high-sSFR M BH -M relation and the empirical scaling relations.

In this section, we showed that, on average, lower sSFR galaxies (quenched or on their way to quiescence) tend to host more massive BHs in cosmological simulations. Since the M BH -M relation is tight in simulations, the effect is small. This is in agreement with the results of Thomas et al. (2019) for SIMBA. Moreover, this is consistent with what is found in observations (Terrazas et al. 2017), but such observational samples with estimates of both dynamical BH mass and SFR are still restrictively small. Our previous paper (Habouzit et al. 2021) provides more information and comparisons between the simulated population of BHs in these simulations and observations.

In this paper, we aim at deriving the galaxy total X-ray luminosity of galaxies hosting faint AGNs and comparing it to observational constraints, for different types of galaxies. As a first step we derive the X-ray emission of the AGN, a crucial component of the total luminosity of galaxies.

We present in Fig. 6 the median of the AGN luminosity as a function of the stellar mass of galaxies. In order to focus on analysis on the faint AGN, we only consider AGNs with hard X-ray luminosity below L AGN ∼ 10 42 erg s -1 . In practise, we use the luminosity limit of the COSMOS survey to detect AGNs as individual sources. This allows us to evolve the luminosity limit with redshift in a way that we will be able to compare our results with observations in Section 7. In Section 7, we will also apply this limit and describe in detail how we compute it. For our analysis we remove all AGNs that are brighter than this limit. In this section, we also do not consider AGN obscuration. The shaded regions indicate the 15th-85th percentiles of the distribution. A median value is only calculated for stellar mass bins with more than 10 galaxies, otherwise individual galaxies are shown as points.

In general, the median AGN luminosity increases with redshift. Furthermore, the galaxies with a higher sSFR (e.g. the blue lines in Fig. 6) have higher luminosities in the simulations than lower sSFR galaxies (red lines), at fixed stellar masses. There is an exception in Illustris at redshifts z = 0, 1. Here, the AGN luminosity of the low-sSFR sample (in red) is higher for stellar masses between M = 10 10 M and a few 10 10 M . The median AGN luminosity drops for more massive galaxies in Illustris.

In the panels of the TNG100 simulation, the low-accretion mode of the AGN feedback (kinetic mode) leads to a sharp decrease of the AGN luminosity (Weinberger et al. 2018) for galaxies with stellar mass log 10 M /M = 10.25 at redshifts z = 0, 1, 2, 3 in the low- Individual galaxies are shown as points in bins with less than 10 galaxies. The shaded regions indicate the 15-85th percentiles. We do not consider AGN obscuration here. All AGNs that could be detected individually by the COSMOS survey are excluded (black dashed line). AGNs are on average brighter in galaxies with high sSFR. AGN feedback is responsible for the decrease of luminosity in massive galaxies, in some simulations.

sSFR sample. We note here that the sharp decrease is present at z = 0, but masked by the low L AGN in lower mass galaxies due to the rarefaction of cold gas (lower accretion rates on to the BHs) at low redshift.

A large shaded region indicates a broad distribution of luminosities in a given stellar mass bin. The size of the shaded regions decreases with increasing redshift. Moreover, the shaded regions are large for low-sSFR galaxies, smaller for intermediate-sSFR galaxies, and even smaller for high-sSFR galaxies. Many different luminosities are possible for low-sSFR galaxies depending on how much gas there is still available for the central BH. The amount of gas and its properties (e.g. temperature) can be affected by both the SN feedback in low-mass galaxies and by AGN feedback primarily in more massive galaxies (e.g. Habouzit et al. 2017;Nelson et al. 2019a;Zinger et al. 2020). In general, this leads to broad distributions.

In the previous section, we derived the AGN luminosity of the host galaxies of faint AGNs. In this section, we derive the luminosity of the XRB population of the simulated galaxies, and evaluate whether the XRB luminosity can be higher than the AGN luminosity, and for which galaxies.

In Fig. 7, we show as an example the X-ray luminosity of the XRB population in the simulated galaxies of Illustris, for different redshifts. We use the relation of Lehmer et al. (2019), which is the most recent among the models discussed in this paper. The error bars indicate the uncertainties on the parameters in equation ( 11). At high redshift (z ∼ 4), most of the simulated galaxies form stars efficiently and have high sSFR. For the Illustris simulation and this given empirical XRB relation, this corresponds to hard X-ray luminosity of the XRB populations in the range L XRB = 10 41 -10 43 erg/s. With time, there are more and more galaxies with lower SFR. As a consequence the diagram extends towards low sSFR values, and also towards lower XRB population luminosity. At z = 0, the Xray luminosity of the XRB population in the Illustris galaxies is L XRB 10 41 erg/s.

Fig. 8 compares the luminosity contributions from the AGN to the luminosity contribution from XRBs. Here, we use the model of Lehmer et al. (2019). We show the contributions for different redshift with different colours. There are 50 bins (per dimension) in the 2D histogram and a contour line is drawn if at least five galaxies are present. We draw a grey solid line to show equal contributions from the XRB population and the AGN. Since we only consider galaxies with M ≥ 10 9 M , the minimum XRB luminosity is L XRB = 10 38.15 erg s -1 (minimum possible value from the empirical XRB scaling relations). At high redshift (z = 4), we find that in the Illustris, TNG100, and SIMBA simulations the X-ray luminosity of the AGN is higher than the luminosity of the XRB population for most of the galaxy population. At z = 4, almost all the galaxies are on the main sequence, independently of their stellar mass. The emission of the XRB population is more important for galaxies with higher sSFR, at fixed stellar mass (see Fig. 2). The XRB populations have high luminosities of L XRB = 10 41 -10 43 erg/s. In these highredshift star-forming galaxies, there is also gas available to feed the central BHs, and therefore, we also find bright AGNs in these galaxies with L AGN = 10 41 -10 45 erg/s, on average. At z = 4, > 90 per cent of the galaxies in Illustris, TNG100, and SIMBA have a higher relative AGN luminosity L AGN > L XRB . We note that the picture is different in the EAGLE simulation. The XRB population in the EAGLE simulation covers the same luminosity range. However, there is a larger diversity of luminosity for the AGN, and they have (much) lower luminosity than the other simulations, on average. As a result, at z = 4 a large fraction of the galaxies have L XRB > L AGN (below the grey line). Only 13 per cent of the galaxies have high AGN luminosities compared to their XRB populations.

With time, the luminosity of the XRB population generally decreases, as well as the luminosity of the AGN. In all the simulations, more and more galaxies reach a regime with L XRB > L AGN . At z = 0, most of the galaxies have L XRB > L AGN in EAGLE (only 5 per cent of galaxies have L AGN > L XRB ). This is also the case for a large fraction of the Illustris simulation with 10 per cent of galaxies with L AGN > L XRB . Finally, we note that in TNG100 and SIMBA a large number of galaxies still host bright AGNs (∼ 63 per cent for both simulations), and therefore are still in the L AGN > L XRB regime (above the grey line). Now if we only consider faint AGNs in Fig. 8, i.e. by only considering the galaxies below an horizontal line at L AGN 10 40 erg/s or L AGN 10 42 erg/s (shown as grey dashed lines in Fig. 8), we find that a significant fraction of the galaxies are dominated by the XRB populations. We quantify this in Table 2. The percentage of XRB dominated galaxies with AGN luminosities L AGN ≤ 10 42 erg s -1 is 92 per cent in Illustris at redshift z = 0 (compared to 90 per cent in the full sample). This increases further if we only include galaxies with AGN luminosities L AGN ≤ 10 40 erg/s (98 per cent). We find similar trends for TNG100, EAGLE, and SIMBA.

In the previous section, we have considered all galaxies independently of their properties (M , SFR). In the following, we Figure 8. Hard X-ray (2-10 keV) luminosity of the AGN and XRB population of the simulated galaxies. Galaxies with AGN luminosity smaller than 10 30 erg s -1 are plotted at that luminosity. Different colours indicate different redshifts. We only consider galaxies with M ≥ 10 9 M . We draw a grey solid line to show equal luminosities from the XRB population and the AGN. The fraction of galaxies with L XRB > L AGN increases at lower redshifts. In the faint AGN regime (illustrated by the regions below the horizontal dashed lines) more and more galaxies have higher XRB luminosity than AGN luminosity with time.

Table 2. Percentage (per cent) of all the simulated galaxies with M 10 9 M with L XRB > L AGN , at redshift z = 4, 2, 0. In the second and third columns, we only consider galaxies with a hard X-ray (2-10 keV) AGN luminosity ≤10 42 erg s -1 or ≤10 40 erg s -1 and show the percentages of galaxies with L XRB > L AGN again. We neglect AGN obscuration here. investigate the relative luminosities of the AGN compared to the XRB populations as a function of galaxy properties.

For illustration we show in Fig. 9 the L AGN /L XRB ratio between AGN and XRB luminosity as function of stellar mass and SFR for all simulations at z = 0. The median L AGN /L XRB vary from one simulation to another. In Illustris and EAGLE, there is a significant population of galaxies for which the XRB population outshine the AGN, and even more so for lower mass galaxies and for galaxies with lower SFR. In TNG and SIMBA, instead, AGNs dominate the galaxy X-ray luminosity for M 10 10.5 M and also in galaxies with log 10 SFR/(M /yr) -1.

We quantify our findings and the redshift evolution below, and in Table 3. In Illustris, almost all low-mass galaxies have a higher XRB luminosity than AGN luminosity, for all the sSFR galaxy groups. Indeed, we find that 81 per cent of the low-mass high-sSFR galaxies with M = 10 9 -10 9.5 M have a higher XRB luminosity than AGN luminosity. We find 96 per cent in the intermediate-sSFR sample and > 99 per cent in the low-sSFR sample. For more massive galaxies in the range M = 10 9.5 -10 10.5 M , the XRB contribution dominates in 75 per cent in the high-sSFR sample, 91 per cent in the intermediate-sSFR sample, and 93 per cent in the low-sSFR sample. The corresponding values for the other simulations are given in Table 3. At high redshift, the galaxy total luminosity is dominated by the AGN emission in general for all the simulations, except EAGLE which produces fainter AGNs. More galaxies are dominated by the XRB emission with time, for all the Illustris, TNG100, EAGLE, and SIMBA simulations. In more detail, the number of Illustris and EAGLE galaxies dominated by XRB at low redshift decreases at higher galaxy stellar mass, for all the sSFR galaxy subsets. However, in TNG100 the number of galaxies with L XRB > L AGN increases for more massive galaxies. This is because at low redshift a significant fraction of massive galaxies experiences quenching. In general for all the simulations, we find that the XRB emission is more likely to dominate the total galaxy emission in low-sSFR galaxies than in higher sSFR galaxies.

In this section, we demonstrated that while the simulations can have similar trends (more low-sSFR galaxies with L XRB > L AGN ), the XRB emission can dominate already at high redshift (EAGLE), or not (Illustris, TNG100, SIMBA), can dominate more in low-redshift low-mass galaxies (Illustris, EAGLE) or more in massive galaxies (TNG100).

Our findings here have a large implication for the detection of AGNs in dwarf galaxies. As indicated in Table 3, we find that all the simulations presented here have a very large fraction of low-mass galaxies (M 10 9.5 M ) with L XRB > L AGN at z = 0. This makes the confirmation of the presence of an AGN in these galaxies challenging, and is in agreement with results from the current search for these AGNs (Reines, Greene & Geha 2013;Baldassare et al. 2015;Mezcua et al. 2016Mezcua et al. , 2018;;Mezcua, Suh & Civano 2019;Greene et al. 2020;Reines et al. 2020). Interestingly, some simulations predict a large fraction of L XRB > L AGN galaxies independently of the SFR of these galaxies (Illustris, EAGLE), some other simulations predict that AGN emission dominates in star-forming main-sequence and starburst galaxies (TNG100, SIMBA). Faint AGNs in low-mass galaxies should be more detectable in hard X-rays in galaxies forming stars more efficiently. If the trend identified in TNG100 is correct, an enhancement of X-ray emission due to the AGN would be observed more Table 3. The percentages of all galaxies dominated by XRB emission (L XRB > L AGN ) are given for different sSFR and stellar mass samples. Percentages of galaxies hosting faint AGNs with L AGN ≤ 10 42 erg s -1 with L XRB > L AGN are also given. Percentages for subsamples with less than 10 galaxies are written in parentheses. We do not consider AGN obscuration here.

Percentage of all galaxies and faint AGN hosts with L XRB > L AGN ( often in low-redshift star-forming galaxies than in more quiescent galaxies.

In the previous sections, we have analysed the luminosity of the faint AGN, but also the luminosity of the XRB population in their host galaxies, as a function of galaxy properties. In this section, we predict the average total hard X-ray (2-10 keV) luminosity of the faint AGN host galaxies, from the Illustris, TNG100, EAGLE, and SIMBA simulations. We compare our predictions to recent observations of stacked galaxies (Fornasini et al. 2018).

We compare the total X-ray luminosity of the galaxies, defined as L total = L AGN + L XRB , in the different simulations. We do not include the emission from the hot gas. Here, we purposely work in the hard X-ray band (2-10 keV), in which the emission from the hot gas is thought to be the smallest contribution and much lower than the XRB contribution (Lehmer et al. 2016). The hot gas contributes to the X-ray emission with a diffuse, soft thermal component (Fornasini et al. 2018). We address the contribution of the hot gas further in Section 7. In Fig. 10, we show the median of the total X-ray luminosity of the simulated galaxies for different stellar mass bins for the four simulations in coloured lines. We show the median AGN luminosity (without obscuration) as black solid lines. Similarly as in Fig. 6, we only consider the galaxies hosting faint AGNs, below the individual detection limit of the X-ray COSMOS survey. In the regime of faint AGN host galaxies, we have demonstrated in the previous section that the luminosity of the XRB populations could overshine those of the AGN. This, of course, also depends on the X-ray scaling relations that we use to compute the XRB population luminosity. To understand the impact of these different relations, we show the median values of the galaxy total luminosity after applying the different scaling relations (indicated by different linestyles for the coloured lines). The shaded regions and individual galaxies (plotted as points) are only given for the XRB scaling relation of Lehmer et al. (2019).

In general, the median galaxy total hard X-ray luminosity increases with stellar mass. The shape of the total luminosity as a function of the galaxy stellar mass follows the shape of the median AGN luminosity, when the latter is brighter on average than the XRB population. In Illustris and EAGLE, the median total luminosity increases with stellar mass, up to the most massive galaxies. In TNG100, we find a strong evidence for the efficient low-accretion mode AGN feedback effect in massive galaxies, especially in the low-SFR sample. We find the same feature in the SIMBA populations. Depending on the scaling relations for the XRB that we employ, the decrease due to AGN feedback of the galaxy total hard X-ray luminosity in massive galaxies can be washed out by high luminosities of the XRB or can also still be there in case of lower XRB luminosities.

Different XRB scaling relations have a strong influence on the total luminosity, changing it by more than one order of magnitude depending on the relations employed. In Illustris, TNG100, and SIMBA, the choice of the XRB scaling relation particularly affects the low redshifts z ≤ 1. We note that it affects TNG100 and SIMBA even more in the low-sSFR galaxies (red lines): particularly in massive galaxies of M 10 10.2 M in TNG100 at z ≥ 1 and all 1. Black solid lines indicate the median AGN luminosity (without obscuration) as shown in Fig. 6, and entering in the computation of L total . We also show the effect of AGN obscuration with the other black lines:

(dashed for model M1, dotted for M2, dash-dotted for M3, and loosely dotted for M4). A larger fraction of obscured faint AGN leads to a more linear shape of the L total -M relation, because the total luminosity is dominated by the XRB contribution.

galaxies at z = 0, and both in low-mass galaxies of M 10 9.5 M and massive galaxies of M 10 10.5 M in SIMBA at all redshifts. In EAGLE, the choice of the XRB scaling relation impacts all galaxy types (high-sSFR, intermediate-sSFR, and low-sSFR samples), at all redshifts. The differences that we find here in the impact of the scaling relation are inversely related to the median of the AGN luminosity (black lines in Fig. 10). If the AGNs are luminous and brighter than the XRB population the choice of the scaling relation does not make a huge difference, because the XRBs do not contribute significantly to the total luminosity of the galaxies. This is the case in most of the simulations at high redshift. Instead, if the median AGN luminosity is small, as it is the case in EAGLE at all redshifts, the contribution of the XRBs becomes more important, and the choice of the XRB scaling relation as well.

The obscuration of a given fraction of the faint AGN could reinforce the contribution of the XRB in the total X-ray luminosity of the galaxies. In Fig. 10, we show the median hard X-ray luminosity of the AGN populations with different models for their obscuration in dashed (model M1), dotted (M2), dash-dotted (M3), and loosely dotted (M4) lines. We only plot the models for which we decrease the luminosity of the obscured faint AGN by one order of magnitude. For the Illustris, TNG100, and SIMBA simulations, the largest impact of our obscuration models is found at the highest redshift (z = 4). At this redshift the obscuration can lead to almost one order of magnitude decrease in the AGN median luminosity. In general, we find very little difference between our models M1, M2, M3, and M4 in the median of the AGN luminosity. At lower redshifts, the effect of obscuration slightly diminishes but is still significant. In some simulations, the effect of obscuration decreases with the host galaxy stellar mass. In practice, if there is a large fraction of obscured AGNs among the faint AGNs, the relative contribution of the XRB population to the observable galaxy X-ray luminosity would be larger. Consequently the shape of the L total -M relation would be almost completely driven by the XRB in these simulations, for all galaxy sSFR groups, and at any stellar mass M (except e.g. the low-mass end of TNG100 and SIMBA).

Here, we compare the population of simulated faint AGNs to observations of stacked galaxies presented in Fornasini et al. (2018). The observational study presents the total galaxy X-ray luminosity extracted from ∼75 000 stacked star-forming galaxies in the redshift range 0.1 < z < 5, using the Chandra COSMOS Legacy survey.

We follow Fornasini et al. (2018), and we divide the simulated galaxy populations into three new subsets. The high-sSFR subset now consists of galaxies with an sSFR higher than 10 -8.5 yr -1 . The intermediate-sSFR subset includes galaxies whose sSFR range from 10 -9.5 to 10 -8.5 yr -1 . Finally, the low-sSFR subset includes all galaxies with a sSFR lower than 10 -9.5 yr -1 . To mimic the observation selection, we also now exclude galaxies with sSFR 10 -11 M /yr from our low-sSFR sample. For the simulations, this means that we remove galaxies with very low sSFR, which in some simulations corresponds to a non-negligible number of quenched galaxies. The new samples are shown in Fig. A1 for all the simulations. Compared to our previous samples, the most important difference is that the SFR range covered by galaxies in a given sample will not change with redshift. For example, galaxies in the intermediate-sSFR sample are mostly on the star-forming main sequence of the simulations at z = 4 or z = 2, while they are clearly among the most star-forming galaxies at z = 0. We provide a summary of the sample definitions employed in this section and in the previous sections in Table 4.

As before, we only include galaxies whose AGN could not be detected as a point source by the X-ray instrument Chandra (COS-MOS survey), i.e. galaxies with AGN luminosities smaller than the sensitivity limit of the COSMOS survey.

The COSMOS survey luminosity sensitivity to detect objects depends on a set of parameters, including their distances to the instrument (i.e. their redshift), and the sensitivity of the instrument. Here, we use the flux-area sensitivity curve from Civano et al. (2016) (their fig. 16, for the 0.5 -2 keV band). The sensitivity limit of the COSMOS survey depends on the area of sky covered by the instrument, with a lower sensitivity at the edges of the pointings that are combined to build the survey. We compute the 0.5 -2 keV flux sensitivity of the COSMOS survey for the Illustris, TNG100, EAGLE, and SIMBA simulations at redshifts z = 0, 1, 2, 3, 4 (using the simulations cosmology). We convert these sensitivity limits from the 0.5 -2 keV band to the 2-10 keV band with the following kcorrection:

with γ = 1.4 (De Luca & Molendi 2004). For redshift z = 0, we directly use the same flux sensitivity limit as Fornasini et al. (2018). The flux received from a given object with an intrinsic luminosity L depends on the redshift of this object: F = L / 4πd 2 L , with d L the luminosity distance. The luminosity distance d L is computed from the angular diameter distance with d L = (1 + z) 2 d A . We invert the flux/luminosity equation and convolve it with the instrument flux sensitivity limit to compute the luminosity limit for individual AGN detection in the COSMOS survey. We use this luminosity upper limit in the following: we only keep galaxies with AGNs fainter than this limit. At z = 0, the hard Xray COSMOS limit is L AGN ∼ 2 × 10 41 erg/s for all the simulations. The limit is L AGN ∼ 6.2-6.5 × 10 42 , 1.9-2.8 × 10 43 , 3.8-5.9 × 10 43 , 6.2-8.0 × 10 43 erg/s for z = 1, 2, 3, 4, respectively. The exact value for each simulation depends on the cosmology and volume of the given simulation, but the variations do not affect our findings.

The effect of the different cuts that we use is shown in Fig. 11, at z = 2 for illustration. Solid lines indicate the median AGN luminosity in the three galaxy samples high-sSFR, intermediate-sSFR, and low-sSFR, without any cuts. We first apply the upper limit for individual detection in the COSMOS survey (dashed lines). In Illustris, this cut leads to a decrease of the median luminosity for the intermediate-sSFR galaxies (green lines) with log 10 M /M 10 up to one order of magnitude at log 10 M /M ≈ 10.8. In TNG100, the median luminosity decreases for log 10 M /M 9.5. The effect of this cut increases with the galaxy stellar mass. This cut does not affect significantly the low-sSFR galaxies. For the high-sSFR and intermediate-sSFR galaxies, the cut has a similar effect in the SIMBA simulation. We note a stronger effect for the low-sSFR galaxies than in the TNG100. The effect of the upper luminosity cut is the smallest in the EAGLE simulation.

The lower sSFR 10 -11 yr -1 cut only affects the median values for low-sSFR galaxies (red lines) and its effect can be seen as the difference between the dashed line (when we apply the COSMOS cut) and the dotted line (COSMOS cut + sSFR cut). In Illustris, the sSFR cut does not have any effect on the median luminosity. In TNG100, the sSFR cut leads to an increase of two orders of magnitude of the median AGN luminosity in the stellar mass bin 10.25 < log 10 M /M < 10.5 and to an increase of one order of magnitude in the stellar mass bin 10.5 < log 10 M /M < 10.75. This corresponds to the stellar mass ranges in which a significant fraction of these massive galaxies have reduced SFR due to the efficient quenching of the kinetic AGN feedback mode in the TNG model (see Fig. 3). Many low-sSFR galaxies with small luminosities ,4,5,6.1 and Figs 3,5,6,10 Galaxy samples defined as redshift-and simulation-dependent. High-sSFR sample 0.5 dex above the star-forming main sequence. Intermediate-sSFR sample 1 dex around the star-forming main sequence. Low-sSFR sample 0.5 dex below the star-forming main sequence.

Comparison to Fornasini et al. (2018) Used in Section 6.2 and Fig. 11, 12, 13, B1, C1, C2

Galaxy samples identical to Fornasini et al. (2018).

High-sSFR sample sSFR/yr > 10 -8.5 Intermediate-sSFR sample 10 -9.5 < sSFR/yr < 10 -8.5 Low-sSFR sample 10 -11 < sSFR/yr < 10 -9.5

Cuts to mimic the study of Fornasini et al. (2018) sSFR cut sSFR 10 -11 yr -1 Upper luminosity cut Detection of individual AGN with the COSMOS survey in the 2-10 keV band. Depends on redshift and volume of the simulation. AGN brighter than this upper limit are removed from the samples.

Figure 11. Impact of the COSMOS sensitivity cut for individual AGN detection (dashed lines) and the sSFR cut (dotted lines) on the median AGN hard X-ray luminosities (2-10 keV, solid lines) for the three galaxy samples of all the simulations (blue for high-sSFR galaxies, green for intermediate-sSFR ones, and red for low-sSFR galaxies). Only galaxies with sSFR 10 -11 yr -1 , and hosted AGN fainter than individual detections in the COSMOS surveys are considered here. Median luminosity values lower than 10 38 erg s -1 are shown at this value. We only show the results for z = 2, for which we have statistics in the three galaxy samples, but we find similar results for other redshifts. The COSMOS sensitivity cut leads to a lower median AGN luminosity and the sSFR cut leads to a higher median AGN luminosity in the low-sSFR sample.

are excluded which leads to an increase in the median luminosity. We find a similar behaviour in the SIMBA simulation. There is almost no impact of the sSFR cut on the median AGN luminosity in the EAGLE, only a small increase of the luminosity for the low-sSFR galaxies.

We show the median of the total X-ray luminosity of the galaxies L total = L AGN + L XRB (excluding the hot gas component) in Fig. 12 for the high-sSFR sample (blue, left-hand panels), the intermediate-sSFR galaxies (green, middle panels), and the low-sSFR galaxies (red, right-hand panels). Here, we decided to work with the median of the galaxy luminosity but the mean luminosity provides almost identical values. All the limits/cuts explained above are applied.foot_3

Figure 12. Median galaxy total hard X-ray (2-10 keV) luminosities (L total = L + L XRB ). Different line styles are used for the Illustris, TNG100, EAGLE, and SIMBA simulations. The galaxies are divided into three sSFR subsets: red for sSFR < 10 -9.5 yr -1 , green for 10 -9.5 yr -1 <sSFR < 10 -8.5 -1 , and blue for sSFR > 10 -8.5 yr -1 . Furthermore, galaxies with sSFR < 10 -11 yr -1 are excluded as well as galaxies with AGN fainter than the individual COSMOS detection limit which is indicated by the black dashed line (we only show the limit computed for the Illustris simulation for simplicity). The lines indicate the median values for stellar mass bins with at least 10 galaxies, otherwise we show the individual galaxies (coloured circles for Illustris, stars for TNG100, downward triangles for EAGLE, and upward triangles for SIMBA). The shaded regions indicate the 15-85th percentile range of each subset. Here, the XRB correction from Lehmer et al. (2019) is used. The results are compared to the observed total luminosities of stacked galaxies from Fornasini et al. (2018;black dots). The triangle symbols indicate upper limits of the observations. The uncertainty of the observational data is smaller than the width of the symbols used here.

The results are compared to the observed total luminosities of stacked galaxies from Fornasini et al. (2018;black dots). The contribution from the hot gas is included in the total X-ray luminosity in Fornasini et al. (2018), and not in the simulations. The triangle symbols indicate upper limits of the observations. The observed X-ray luminosities are corrected for obscuration based on the measured X-ray hardness ratios which provide a rough measure, and mostly impact the low-redshift data (Fornasini et al. 2018). We discuss this further in the following section.

In the observations, the total X-ray luminosity of the galaxies increases with their stellar mass for both the galaxies in the high-sSFR and intermediate-sSFR samples. For the low-sSFR sample, the correlation is less obvious (except at z = 0, for which an increase with M is favoured), and the total X-ray luminosity is more or less constant for massive galaxies with M 10 10.5 M . In the following, we describe our general findings, first for these models ignoring the impact of obscuration, on the agreement between the simulations and the observations for each of the galaxy sSFR subsets.

High-sSFR samples (sSFR > 10 -8.5 yr -1 ):

We find that all the simulations have a higher median of the X-ray total luminosity of their galaxy population for the high-sSFR samples, compared to the observational constraints of Fornasini et al. (2018). The luminosity for the simulations does not include the hot gas component, and therefore, the difference with the observations could be potentially larger. All the simulations have an increasing galaxy total luminosity with the stellar mass of the galaxies, in agreement with the trend in the observations.

Intermediate-sSFR samples (10 -9.5 yr -1 < sSFR < 10 -8.5 yr -1 ):

The trend obtained for the intermediate-sSFR sample of simulated galaxies is also similar to the observed one. We find a good agreement between the EAGLE simulation and the observations, for all redshift except z = 0. The other simulations overpredict the total X-ray luminosity of their intermediate-sSFR galaxies, on average. The overprediction is stronger towards less massive galaxies, and can be more than one order of magnitude.

Low-sSFR samples (10 -11 yr -1 < sSFR < 10 -9.5 yr -1 ): Finally, for the low-sSFR sample galaxies, we find that on average, the TNG100 and SIMBA simulations overpredict the total luminosity for the faint AGN low-mass host galaxies. A better general agreement is found for the EAGLE and SIMBA simulations, which show an increasing total luminosity with increasing stellar mass. In more detail, we find that the best agreement at z = 3, 2 is obtained for the SIMBA simulation. At z = 1, the Illustris, EAGLE, and SIMBA simulations provide a good agreement with the observations. While the Illustris and EAGLE have increasing galaxy total luminosity with stellar mass, the X-ray median luminosity is decreasing in SIMBA up to M ∼ 10 10 M , and then increasing slightly to M ∼ 10 11 M . For more massive galaxies, the median luminosity of SIMBA stays constant, as in the observations. At z = 0, the best agreement is found for the Illustris simulation. This is particularly true for the low-mass galaxies M 10 10 M where the median luminosity of the Illustris simulation is lower than the SIMBA simulation.

For different redshifts and galaxy types, we find that the highest/lowest galaxy total X-ray luminosity is not always produced by the same simulations. We develop this in the following. For the three galaxy samples and all the redshifts, we find that the EAGLE simulation produces the lowest median total X-ray luminosity (dotted lines). This is because EAGLE is the simulation with the largest population of faint AGN (Rosas-Guevara et al. 2016). The AGN Xray luminosity function of EAGLE is below all the other simulations (e.g. Sijacki et al. 2015;Habouzit et al. 2019;Thomas et al. 2019). The luminosity function is in agreement with Aird et al. (2015) and Buchner et al. (2015) in the range log 10 L x,2-10 keV 10 44 erg/s at high redshift for z ≥ 2, in the range log 10 L x,2-10 keV 10 43 erg/s at z = 2, and is underestimated otherwise. At z = 0, the EAGLE simulation is below the constraints of Aird et al. (2015) and Buchner et al. (2015). These constraints do not cover AGN with fainter luminosity than log 10 L x,2-10 keV = 10 42 erg/s. We find that for the high-sSFR samples the highest median is in Illustris for M 10 10 M , and TNG100 for M ∼ 10 10 M , and in SIMBA for more massive galaxies. For the intermediate-sSFR samples, TNG100 has the highest median for galaxies with 10 9.5 M M 10 10.5 M at z = 2, 3, 4, while SIMBA has the highest median at z = 1 and for galaxies with M 10 10.5 M at z = 1, 2, 3, 4. TNG100 produces the highest total luminosity median for low-sSFR sample galaxies of M a few 10 10 M , and SIMBA and Illustris for more massive galaxies.

Interestingly, we know from Fig. 11 that the fact that observations here exclude galaxies with sSFR < 10 -11 yr -1 can artificially boost the galaxy total X-ray luminosity for the low-sSFR galaxies. This would particularly affect the median in the bins containing massive galaxies. In our analysis we use the same sSFR cut as in the observations of Fornasini et al. (2018), however, this does not ensure that the sSFR distributions of the simulated and observed samples are the same. Depending on the number of galaxies with sSFR close to the limit sSFR < 10 -11 yr -1 in the observations, we may still be counting too many simulated galaxies with very low sSFR. This could explain the lower total luminosity in massive galaxies, particularly for the TNG100 and SIMBA simulations (as they were already affected by the sSFR < 10 -11 yr -1 cut in Fig. 11).

Observing faint systems is challenging. We cannot rule out that the faintest systems that we see in simulations are actually missed by the observations. Observations are always limited when trying to observe fainter and fainter systems. The completeness of the observational samples is particularly difficult to address at the faint end of the galaxy distribution. This can potentially create significant discrepancies with the samples of simulated galaxies, since in the latter we have access to the faintest galaxies and AGNs. Therefore, in a separate test, we applied different lower luminosity cuts to our samples of simulated galaxies. In practice, to mimic what could take place in the observations we assumed that (i) galaxies with total X-ray luminosity below 10 40 or 10 41 erg/s (for z = 0, 1 and z = 2, 3, 4, respectively, following Fornasini et al. 2018) are not detected, or (ii) that galaxies with total luminosity two or three orders of magnitude below the COSMOS Legacy limit are not detected. We found that missing the faintest galaxies when observing the simulations would mostly impact the total luminosity of lowmass galaxies with M 10 10.5 M for the strongest luminosity cut. The median/mean luminosity of the detected galaxies would be increased compared to the median/mean luminosity of the full galaxy population with the same physical characteristics, but it would not strongly impact our main conclusions regarding the agreement between simulations and the observations of Fornasini et al. (2018).

The observations of Fornasini et al. (2018) already include some rough estimate for obscuration, assuming that the XRB and the AGN are subject to the same obscuration (which may or may not be a valid assumption, see the discussion in Fornasini et al. 2018). Their correction mostly impacts the total luminosity at low redshift. An upper limit on the impact of the correction is ∼0.4 dex difference for z < 0.6, 0.2 dex difference in the range 0.6 < z < 2.3, and negligible differences at higher redshift. Given the uncertainties on this average correction applied to the XRB and AGN, and in order to quantify the role of AGN obscuration on the total luminosity, we apply our correction to the simulations.

In Fig. 13, we show the same figure as Fig. 12, i.e. the median of the galaxy hard X-ray luminosity, but this time we account for obscured AGN. We note that there are no strong differences between the model variation M1, M2, M3, M4 in our results, which means that the fractions of obscured AGN with L AGN 10 41 erg/s does not impact the total galaxy luminosity. In practice, this is because the median of the galaxy X-ray luminosity is often higher than 10 41 erg/s. We apply the M1 obscuration model to the simulated AGN samples. More precisely, we show two sets of obscuration models, for which we either decrease by one order of magnitude the AGN luminosity (fainter AGN models), or either assume L AGN = 0 (missed AGN). We describe below the impact of our two sets of obscuration models on the total galaxy luminosity, and re-examine our previous conclusions. Typically, the median (or mean) L total can be affected by about half an order of magnitude when we account for AGN obscuration (fainter AGN model), and even more if we assume that the obscured AGN would be completely missed by the observations (missed AGN model). With strong AGN obscuration models, the AGN features that were visible in the galaxy total luminosity (i.e. non-linearity of the L tot -M relation, peak due to large fraction of bright AGN, decrease of L tot because of AGN feedback) can be completely erased. Here, we used the XRB scaling relation of Lehmer et al. (2019) which predicts the highest XRB luminosity among all XRB models. Using other XRB models or having more obscured AGN leads to smaller total galaxy luminosity.

In general, we find that obscuration can decrease the total galaxy luminosity by up to one order of magnitude for the fainter AGN models. When we assume that we would completely miss the obscured AGN (missed AGN models), the impact can be higher and reach two orders of magnitudes at stellar masses for which simulations produce a large fraction of AGN. This is the case in the TNG100 and SIMBA simulations, which have a large fraction of AGN for galaxies with stellar masses of M 10 10.5 M . When we set L AGN = 0 for the obscured AGN, the total galaxy luminosity is fully driven by the XRB population.

For the high-sSFR galaxies, adding the obscuration models does not impact our previous conclusions: all the simulations overpredict the total galaxy luminosity. For all the simulations except EAGLE, the galaxy luminosity of the high-sSFR subsets is dominated by the AGN luminosity, which the decrease of the luminosity when we obscure some of the AGN. We note that even when we fully remove the obscured AGN, i.e. when the galaxy luminosity is dominated by the XRB population, the luminosity is higher than for the observations. This could mean that the empirical scaling relations for the luminosity of the XRB population overpredict the median/mean hard X-ray galaxy luminosity.

For galaxies in the intermediate-sSFR group and for most of the simulations, our fiducial model with the XRB scaling relation of Lehmer et al. (2019) and the fainter-AGN obscuration models still produce higher galaxy luminosity (see Fig. 13). A better agreement is found if we assume lower XRB scaling relations (Lehmer et al. 2010(Lehmer et al. , 2016;;Aird et al. 2017;Fornasini et al. 2018) or that the obscured faint AGNs are completely missed by the observations.

For galaxies in the low-sSFR subsets in Illustris and EAGLE, we find that the total galaxy luminosity is in general smaller than in the observations. For these two simulations, the discrepancy would be larger for other XRB scaling relations (as shown in Fig. 10) or obscuration models accounting for more obscured AGN. In SIMBA, the agreement between the total galaxy luminosity from simulations (with the L19 XRB relation and without AGN obscuration) and observations is good for z > 1. We note some differences at z = 0, with a higher galaxy luminosity median for simulated galaxies with M 10 10 M , and lower median for M 10 10.5 M than in observations. Obscuration models assuming more obscured AGN would in general lower the good agreement of SIMBA with observations. TNG100 has a particular shape of the total galaxy luminosity median, which is not clearly visible in observations, with a higher luminosity median in low-mass galaxies, and lower luminosity median in massive galaxies, compared to observations. Assuming more obscured AGN in the simulations reduces the differences with observations.

From our analysis we find that for some simulations such as Illustris or EAGLE, whose total luminosity at low redshift is dominated by XRB, are not impacted by our modelling of AGN obscuration. For the other simulations with a higher contribution of the AGN luminosity to the total galaxy luminosity, without considering our obscuration correction we often find an excess of luminosity in the simulations compared to the observations. We have demonstrated that this excess can be reduced by obscuration at low redshift for all intermediate-sSFR galaxies and high-sSFR galaxies with M 10 10.5 M -10 11 M (depending on the simulation). Similarly, the excess could be due to an overestimate of the obscuration in the observation data (Fornasini et al. 2018) at low redshift.

In this section, we discuss different aspects that could impact our results, and our comparison with current observations of galaxy total hard X-ray luminosity.

The BH subgrid models of these simulations have been broadly calibrated against one of the empirical M BH -M bulge relations available in the literature. While the Illustris, TNG100, and EAGLE simulations adjust the efficiency parameter of the AGN feedback model, the simulation SIMBA instead calibrates the accretion efficiency. None of the simulations studied here were calibrated against the AGN luminosity function.

At z = 0, these simulations have X-ray luminosity functions in good overall agreement with the constraints of, e.g. Buchner et al. (2015). TNG100 has a higher normalization of the luminosity function (Habouzit et al. 2019), and EAGLE a lower normalization, for log 10 L AGN /(erg/s) 43 (Rosas-Guevara et al. 2016). Illustris and SIMBA lie within the constraints of Buchner et al. (2015) for this luminosity regime (see Sijacki et al. 2015;Thomas et al. 2019, for the analysis of AGN properties with observational constraints). At higher redshift, all the simulations except EAGLE overpredict the observational constraints in the range log 10 L AGN /(erg/s) 42.5 -43 (Habouzit et al. 2021). The regime that we investigate here (log 10 L AGN /(erg/s) 42) is below the range which is usually constrained by observations of the X-ray luminosity function. The differences between the simulations and the observational constraintsfoot_4 on the X-ray luminosity function could affect the galaxy total X-ray luminosity in a non-trivial way, for the different galaxy sSFR samples.

To make the analysis of the simulations consistent with the observations from Fornasini et al. (2018) (stacked galaxies from the Chandra COSMOS Legacy survey), we have adopted the same sSFR cut (excluding simulated galaxies with sSFR 10 -11 yr -1 ), and the same upper luminosity limit to exclude galaxies with AGN detectable by the Chandra COSMOS Legacy survey. These same cuts do not ensure that the sSFR distributions of the observation and simulation samples are similar. For example, if a high-sSFR simulation sample has a distribution peaking at larger sSFR values than the corresponding observational sample, this will lead to higher total galaxy luminosities in the simulations than in observations. Therefore, we cannot exclude that some discrepancies between the observations and the simulations are due to different distributions of their respectives sSFR samples.

Two additional caveats for our comparison in this paper are the X-ray emission from hot gas and the obscuration.

In our analysis, we have neglected the contribution to the X-ray emission from the hot interstellar medium (ISM). This contribution was, however, studied in the observations considered here. The hot gas is expected to contribute to the X-ray emission with a diffuse, soft thermal component (Fornasini et al. 2018). The hot gas component is thought to be significant and to dominate the X-ray emission over the XRB emission for rest-frame energy of <1.5 keV (Lehmer et al. 2016, their fig. 4). For high energies of >1.5 keV, as it is the case in our analysis, the XRB emission is expected to dominate (Lehmer et al. 2016). There is currently no definitive quantification (from observations or theory) across galaxy masses, galaxy apertures, types, and redshifts. The contribution of the hot gas could not only depend on energy band, but could also decrease with increasing redshift, and could be important mostly for z < 1 (Lehmer et al. (2016), results obtained assuming that the spectral energy distribution (SED) of star-forming galaxies do not change strongly with redshift. According to Mineo et al. (2012b), the hot gas X-ray emission could increase with the SFR of a galaxy. In that case, the influence of the hot gas emission would be higher in the high-sSFR galaxy sample at a fixed stellar mass and at higher stellar masses at fixed sSFR. From a numerical and theoretical perspective, definitive assessments of the X-ray emission from the hot gas is still needed, but some preliminary estimates of the gas X-ray emission within galaxies in the soft bands have been derived from the Illustris and TNG simulations (Truong et al. 2020). We postpone the task of including the contribution of the X-ray emitting hot gas within galaxies to future work, as this requires a careful assessment of the dependence on X-ray wavelength and on aperture within which the mock or real observations are taken.

Gas and dust obscuration can play a crucial role when comparing observations to simulations. It is not clear yet if the obscuration originates from material close to the galaxies nuclear region. In that case, the luminosity of the XRB population distributed within the galaxies would not be strongly obscured, but the AGN would be. The obscuration could also originate from material in the whole galaxy, and in that case it would lead to the obscuration of both the XRB population and the AGN. Buchner & Bauer (2017) find that the gas on the galaxy scale is only responsible for a part of the Compton-thin AGN, and does not provide Compton-thick lines of sight. The heavily obscured (Compton-thick) AGNs would therefore mostly result from obscuration in the nuclear region. If this is the case, the hard X-ray emission from the galaxy-wide XRB population would be less impacted than the emission from the AGN.

Simulations do not consistently capture the obscuration of the AGN, as obscuration can arise from regions close to the AGN on spatial scales below the simulation resolution. Therefore, we have tested the role of AGN obscuration by applying four different models to the simulated AGN samples. Our models depend on redshift and hard X-ray (2-10 keV) AGN luminosity for AGN with L AGN 10 41 erg/s (Habouzit et al. 2019, for more details), and only on hard X-ray AGN luminosity for fainter AGN. Obscuration could also depend on galaxy SFR, a parameter that we do not consider here. We have not applied any further obscuration model to the XRB populations.

In practice, the galaxy total hard X-ray luminosity can be impacted significantly depending on the fraction of obscured AGN that we assume. We find that the more the faint AGNs are obscured the more the shape of the total galaxy luminosity as a function of galaxy stellar mass is driven by the XRB luminosity (since we do not apply any obscuration model to the XRB emission). If the linearity of the XRB empirical L XRB -M scaling relations is a good estimate (as found in e.g. Lehmer et al. 2019) and in the presence of a large population of obscured faint AGN (or just a population of very faint AGN, such as in EAGLE), we should observe a linear L total -M relation in the observations, independently of the galaxy sSFR. In that case, we find that any deviation from a linear L total -M relation would be due to features of the AGN populations.

It is now demonstrated that AGNs can exist in dwarf galaxies (Reines et al. 2013(Reines et al. , 2020;;Baldassare et al. 2015;Mezcua et al. 2016Mezcua et al. , 2018Mezcua et al. , 2019;;Greene et al. 2020). Quantifying the fraction of galaxies hosting BHs and the BH mass distribution in dwarf galaxies can constrain the theoretical models of BH formation (Greene 2012). Since these galaxies have not evolved much over cosmic times compared to their massive counterparts, they could have retained the initial properties of BH formation: initial BH mass and initial BH formation efficiency. While the BH occupation fraction in low-mass galaxies in such large-scale simulations may not be relevant/accurate because of the simple seeding of BHs in massive galaxies or haloes, the AGN occupation fraction is fundamental to address the connection between the AGN and their host galaxies (e.g. the correlations between AGN activity and the SFR of their galaxies).

The AGN found in dwarf galaxies in observations can generally be qualified as faint AGN (Mezcua et al. 2016;Chilingarian et al. 2018;Mezcua & Domínguez Sánchez 2020, and references therein). In this work, we have shown that in the Illustris, TNG100, and EAGLE simulations the XRB population in galaxies of M 10 9.5 M can outshine the AGN emission in hard X-rays. This is a significant issue when trying to detect an AGN in X-rays; detection in X-ray is one of the most common method to detect low-mass AGN in low-mass galaxies today. What is interesting is that the simulations do not predict the same trend with SFR. In Illustris and EAGLE, the XRB population outshine the AGN in > 90 per cent of the galaxies, whether these galaxies form stars efficiently (starburst) or not (below the star-forming main sequence). However, in TNG100, the XRB population outshine the AGN only in galaxies below the main sequence, but not in main-sequence galaxies or starburst galaxies. There, AGN activity is enhanced when SFR activity is enhanced. Confronting current and future observations of AGN in dwarf galaxies to our results on the AGN population predicted by cosmological simulations will help us to understand the observations and at the same time to constrain our modelling of BH and galaxy physics in simulations.

We have analysed the properties of the faint AGN (in the hard X-ray band 2-10 keV, L AGN 10 42 erg/s) and their host galaxies of the four large-scale cosmological simulations Illustris, TNG100, EAGLE, and SIMBA. We have modelled the contribution from the XRB population and from the AGN (including their possible obscuration) to the total galaxy hard X-ray luminosity. We summarize below our main findings.

(i) The properties of the faint AGN host galaxies vary from simulation to simulation (Fig. 4). Faint AGN of L AGN ∼ 10 38 erg/s can be powered by relatively massive BHs and be located in massive galaxies (M 10 10 M ) with reduced SFR (TNG100, SIMBA), or be powered by lower mass BHs in less massive galaxies (M 10 10 M ) still forming stars (Illustris, EAGLE). (ii) We find that the two possible behaviours described above depend on the effectiveness of AGN feedback in massive galaxies. In TNG100 and SIMBA, the efficient feedback taking place in massive galaxies reduces both their sSFR (Fig. 4), but also the ability of their BHs to accrete efficiently. AGN feedback drives the build up of the faint AGN population in these simulated galaxies.

(iii) In all the simulations, except EAGLE, most of the galaxies have brighter AGN than the XRB population, at high redshift (z > 2, Fig. 8). With time, the AGN number density decreases and consequently more and more galaxies have a brighter XRB population than their AGN. The general fainter population of AGN in EAGLE (at all redshifts) compared to the other simulations leads to a significant number of galaxies with a brighter XRB population than AGN.

(iv) The relative contribution of the AGN and the XRB population to the XRB + AGN total galaxy hard X-ray luminosity depends on the stellar mass and the SFR of the galaxies (Fig. 10). Starburst galaxies host brighter AGN in all simulations and across redshift: the AGN luminosity dominates over the XRB luminosity in most of these galaxies. At low redshifts (z ≤ 3), the efficient AGN feedback in TNG100 and SIMBA leads to a strong decrease of the median AGN luminosity in massive galaxies (M 10 10 M ) with reduced sSFR, and more galaxies are dominated by XRB emission.

(v) In low-mass galaxies of M 10 9.5 M at z = 0, we find that the XRB emission always outshines the AGN emission in low-sSFR galaxies in all the simulations (Table 3). The XRB still dominates in main-sequence and startburst galaxies in Illustris and EAGLE, but does not outshine the AGN in TNG100. This has important implications for the search of AGN in dwarf galaxies.

(vi) The total AGN + XRB hard X-ray luminosity of faint AGN host galaxies (i.e. neglecting the hot ISM X-ray emission) increases with increasing M , for all redshifts and all the simulations (Fig. 10). We note a turnover for the massive TNG100 and SIMBA galaxies for which the lower AGN median luminosity (due to AGN feedback) propagates to the total galaxy hard X-ray luminosity.

(vii) We find that a non-linear L total -M relation in faint AGN galaxies (Fig. 10 and Fig. 13) is explained by a non-linear L XRB -M scaling relation (in that case XRB luminosity models need to be updated), or by peaks of AGN activity at some stellar masses. We find that the obscuration of faint AGN can completely erase these AGN signatures in the L total -M relation (see Fig. 13). In that case, the shape of the L total -M relation would be fully driven by the XRB emission.

(viii) The simulations, with our modelling of AGN and XRB luminosity, tend to overestimate the total AGN + XRB galaxy X-ray luminosity in the high-sSFR sample and for most of the simulations in the intermediate-sSFR sample (neglecting the hot gas ISM X-ray emission) compared to the observations of the COSMOS Legacy stacked galaxies (Fornasini et al. 2018). Simulated galaxies with sSFR > 10 -9.5 yr -1 are too bright. Galaxies with 10 -9.5 yr -1 < sSFR < 10 -8.5 yr -1 are also too bright, except a good agreement for EAGLE. For low-sSFR galaxies of 10 -11 yr -1 < sSFR < 10 -9.5 yr -1 , we find that some simulations underestimate or overestimate the median galaxy luminosity (Fig. 12).

(ix) In both simulations and observations (Fornasini et al. 2018), high-sSFR galaxies have higher total galaxy X-ray luminosity than low-sSFR galaxies at fixed stellar mass, in general (Fig. 12).

(x) The empirically driven XRB scaling relations used in this work span 0.5 dex in luminosity (at fixed M ), which is about the same order of magnitude as some of our obscuration models (Fig. 13). These two aspects are highly degenerate and further observational constraints will be needed to disentangle them.

Our work and predictions pave the way for upcoming and concept space missions such as Athena, AXIS, and Lynx. These missions will increase by several orders of magnitude the sensitivity of the current X-ray instruments, and will allow us to make promising progress on our understanding of faint AGN, a luminosity regime that as we have demonstrated can be dominated by XRBs for specific sSFR and M regimes.

In Table A1, we present the parameters α and β of the main sequence of each simulation for different redshifts as discussed in Section 2.4.1. A2. Percentage of galaxies (per cent) and number density of galaxies n (10 -5 cMpc -3 ) in the three high-sSFR, intermediate-sSFR, and low-sSFR samples, for redshift z = 0, 1, 2, 3, 4. Only galaxies with M 10 9 M are included. This table refers to the division shown in Fig. 3 and Fig. A1.

In Table A2, we present the percentages and number densities for the three different subsets that we used in Sections 3, 5, 6, and 7. We only consider galaxies that are well resolved in all the simulations, i.e. M 10 9 M . The high-sSFR sample groups starburst galaxies 0.5 dex above the main sequence. The intermediate-sSFR galaxies represent the galaxies on the main sequence (within 1 dex). Finally, low-sSFR galaxies are galaxies below 0.5 dex of the main sequence.

To compare the simulations to the results from observations of Fornasini et al. (2018), we have changed the definition of our galaxy sSFR samples. High-sSFR galaxies are defined by sSFR/yr > 10 -8.5 , intermediate-sSFR galaxies by 10 -9.5 < sSFR/yr < 10 -8.5 , and low-sSFR galaxies by 10 -11 < sSFR/yr < 10 -9.5 . Fig. A1 shows the SFR as a function of the stellar mass with these new definitions. The background colours show the three sSFR samples. Contrary to Fig. 3, the division does not depend on redshift or simulation and galaxies on the star-forming main sequence of the simulations are not always included in the intermediate-sSFR sample. In Table A2, we add the percentages and number densities of the different subsets for our second set of definitions.

. SFR as a function of galaxy stellar mass for the Illustris, TNG100, EAGLE, and SIMBA simulations. Hexabins are colour coded by the number of galaxies in bins. Finally, we show the three samples that we use with background colours: blue for galaxies with sSFR > 10 -8.5 yr -1 , green for galaxies with 10 -9.5 < sSFR < 10 -8.5 yr -1 , and red for 10 -11 < sSFR < 10 -9.5 yr -1 . Galaxies with sSF R < 10 -11 yr -1 are shown at log 10 SFR/(M /yr) = -4, but are not considered in the comparison with observations. These definitions follow Fornasini et al. (2018), and are different from the ones used in the first sections of the paper and showed in Fig. 3, especially for z < 2.