The Star-Formation Properties of the Observed and Simulated AGN Universe: BAT vs EAGLE

AGN can be identified in many different bands of the electromagnetic (EM) spectrum such as the optical, Infra-red (IR) or radio (Kauffmann et al., 2003; Smolčić, 2009; Lacy et al., 2013). Selection of AGN in the ultra hard X-ray band ($E\gtrsim$ 10 keV), however, can be advantageous. Firstly, there are relatively few sources other than AGN which produce Ultra hard X-ray emission, reducing the chance of contamination in both sample selection and luminosity measurements (Hopkins et al., 2009). Secondly, absorption of Ultra hard X-rays only occurs through high column densities of gas along the line of sight (N${}_{\mathrm{H}}\gtrsim 10^{24}$ cm^-2, e.g. Wilms et al., 2000). Other parts of the EM spectrum such as the Ultraviolet (UV) or optical may be blocked at much lower column densities when gas and dust are present (see Hickox & Alexander, 2018, and references therein).

Selection in the Ultra hard X-ray regime therefore provides a relatively unbiased sample of AGN compared to IR or radio selection (Brandt et al., 2000; Mushotzky, 2004). It can also lead to the discovery of new, highly obscured AGN (Koss et al., 2010; Baumgartner et al., 2013; Lansbury et al., 2017).

We use the 58 month catalogue¹¹1https://swift.gsfc.nasa.gov/results/bs58mon/ of sources from the NASA Swift satellite’s Burst Alert Telescope (hereafter Swift-BAT)²²2https://swift.gsfc.nasa.gov/, which operates in the Ultra-hard X-ray band (14 - 195 keV), to select our sources.

In order to reduce the contamination that Galactic sources may have on any targets in both the the optical and IR photometry described in Section 2.2, objects within $\pm$ 10^∘ of the Galactic plane were removed. In order to reduce statistical biases that can emerge due to volume driven effects (the larger the volume the more statistically likely it is to have objects with extreme luminosities, stellar masses, SFRs etc.) we matched the observational volume to the EAGLE simulation volume (in this study we use the 100 cMpc box size, described in Section 3.1) by applying a cut-off in luminosity distance of 67.8 Mpc to the parent sample in the Swift-BAT 58 month catalogue.

The final observational sample contains 72 objects, which are hereafter referred to as our BAT AGN sample. These differences in sources can be seen in the left hand panel of Figure 1, where the black points represent our BAT AGN sample and the red points represent sources cut from the 58 month catalogue (after removing those sources within the Galactic plane).

To obtain the most up-to-date and uniformly assessed redshifts and X-ray luminosities, the data were matched to the BAT AGN Spectroscopic Survey (BASS³³3https://www.bass-survey.com/, Koss et al., 2017), an optical and IR spectroscopic survey of Swift-BAT sources, estimating redshifts from the optical spectra and collecting the X-ray luminosities from the combined processing of Swift/BAT and softer X-ray datasets from XMM-Newton, Chandra, etc. (Ricci et al., 2017). The distribution in X-ray luminosity and redshift of the BAT AGN sample can be seen in the right hand panel of Figure 1. A table of the sample and its estimated properties can be found in Appendix A.

Host galaxy properties such as stellar mass and SFR can be estimated via Spectral Energy Distribution (SED) fitting. This procedure is described in Section 2.3. In order to fit an SED to an AGN and its host galaxy to estimate these properties, multiple photometric data points spanning the optical and IR are needed. Due to the all-sky nature of the Swift-BAT sample, data from various surveys were collected.

To obtain uniform counterpart photometry in the Near and Mid Infra-red (NIR and MIR) for as many of our observed AGN as possible, we used the all-sky datasets of the Wide-field Infra-red Survey Explorer (WISE)⁴⁴4https://www.nasa.gov/mission_pages/WISE/main/index.html (Wright et al., 2010) and the 2 Micron All Sky Survey (2MASS)⁵⁵5https://www.ipac.caltech.edu/2mass/ (Skrutskie et al., 2006). We compiled MIR photometry from the online database of the WISE mission, served on the NASA/IPAC website, spanning four bands at 3.5, 4.6, 12 & 22 $\mu$m. Since our targets are well-resolved nearby galaxies, we used the GMAG photometric measurements in all cases. These are photometry using apertures scaled to the 2MASS Extended Source Catalog (XSC) (Jarrett et al., 2000) size measurements, with a correction for the Point Spread Function of the WISE images. We converted WISE magnitudes to fluxes using the conversions from the online WISE Explanatory Supplement. The all-sky 2MASS XSC provides photometry in 3 NIR bands (J, H, K_s) for all the BAT AGN. We used the isophotal curve-of-growth magnitudes in this work, and employ standard conversions from 2MASS magnitudes to fluxes from the online 2MASS Explanatory Supplement.

To provide better constraints on the star-forming component of each AGN (as explained in Section 2.3), Far Infra-red (FIR) photometric data were added. We obtained Herschel photometry for the BAT AGN from the data release described in Shimizu et al. (2017), based on the analysis presented in Meléndez et al. (2014) (PACS) and Shimizu et al. (2016) (SPIRE). We used their measurements for the integrated photometry of the sources, covering 5 bands: PACS 70 & 160 $\mu$m and SPIRE 250, 350 & 500 $\mu$m.

The collection of uniform and accurate optical photometry for the BAT AGN sample is a more complex process than for the IR. The NASA Extragalactic Database (NED) supplies a heterogeneous dataset of photometry from the literature. We obtained B band photographic photometry of a subset of galaxies from the NED version of the RC3 catalog (de Vaucouleurs, 1991; Corwin et al., 1994). The SDSS (York et al., 2000) and PanSTARRS (Chambers et al., 2016) optical imaging surveys yielded magnitudes in 6 bands from $u$–$y$. From the photometric catalogs of these surveys, we used the CMODEL magnitudes, which combine two different light profile models of extended sources to produce best-effort integrated photometric measurements. Where possible, we replace the SDSS database magnitudes with those from the NASA Sloan Atlas (NSA), a value-added resource for nearby galaxies with more accurate photometry for extended sources. We also used, where available, HyperLEDA measurements in five bands (UBVRI) which are homogenised curve-of-growth photometry as presented in Prugniel & Heraudeau (1998). Finally, for a small number of galaxies, curve-of-growth galaxy photometry in the BVRI bands is available from the Carnegie-Irvine Galaxy survey (CGS; Ho et al. (2011)).

To ensure the highest quality integrated photometry that is estimated uniformly, we restricted our optical photometry to only one of these datasets for each galaxy, taking the following priority in decreasing order: The Carnegie-Irvine galaxy survey⁶⁶6https://cgs.obs.carnegiescience.edu/CGS/Home.html, the HyperLEDA extragalactic survey⁷⁷7http://leda.univ-lyon1.fr/, the NSA⁸⁸8http://www.nsatlas.org/, SDSS⁹⁹9https://www.sdss.org/, Pan-STARRs¹⁰¹⁰10https://panstarrs.stsci.edu/ and finally RC3 from NED¹¹¹¹11https://ned.ipac.caltech.edu/.

SFRs can be estimated using various techniques; for example conversions from H$\alpha$, total UV light and radio emission from supernova to name a few (see Kennicutt & Evans, 2012, for a summary and references therein). Most of these methods can be significantly contaminated by emission from the AGN itself. The UV can be dominated by emission from the accretion disk, an effect especially prominent for unobscured AGN, estimations of the SFR in the radio can be contaminated by non-thermal radiative processes such as synchrotron radiation emitted from relativistic jets and lobes and/or the central regions and H$\alpha$ can be contaminated due to photoionisation from the AGN.

Star formation contributes strongly at FIR and sub-mm wavelengths. This is due to the re-processing of UV or optical light from young stars via the dust rich clouds that typically surround star forming regions. Similarly, the AGN dusty torus also emits in the IR, however due to the proximity of the torus to the accretion disk the temperature of the dust is higher than that of the dust around star forming regions (see Netzer et al., 2007, and references herein). Therefore the emission peaks in the MIR, dropping off sharply as we reach the FIR, where re-emitted light from star formation dominates (Mullaney et al., 2011; Rosario et al., 2012).

In order to estimate SFRs and stellar masses, we performed template fitting to the photometry ranging from the UV to FIR that comprises the overall SEDs of our sources. We used the FortesFit¹²¹²12https://github.com/vikalibrate/FortesFit SED template fitting routine for this task. FortesFit is a python-based fully Bayesian fitting routine. It contains a library of parameterised model SEDs that work additively in order to try and fit the full SED to the data points provided. Priors of certain parameters can be provided, such as the AGN type classification, in order to better constrain the expected contribution of various components, thereby minimising the chance of mis-fitting. The fitting routine then provides posterior probability distributions for various values such as the AGN IR luminosity component, stellar mass and the integrated luminosity in the IR due to star formation. For further details of the SED modelling, refer to Scholtz et al. (2020)

To convert the integrated IR luminosity component (in erg s^-1, integrated between 8-1000 $\mu$m) due to star formation into the SFR, we use Equation 1 below (equation 6 from Rosario et al., 2016, see the paper for further details of the full assumptions and models). We also assume an intrinsic uncertainty of 0.25 dex in this calibration (Rosario et al., 2016), which we take into account in further analysis.

These estimates and other data for our BAT AGN can be found in Appendix A.

The Evolution and Assembly of GaLaxies and their Environments (EAGLE) is a suite of hydrodynamical simulations ¹³¹³13http://icc.dur.ac.uk/Eagle/ (Crain et al., 2015; Schaye et al., 2015). The simulations use an improved version of the GADGET-3 Smoothed Particle Hydrodynamics (SPH) code (Springel, 2005), with the cosmological parameters from the Planck mission (Planck Collaboration et al., 2015) implemented, whereby $\Omega_{\Lambda}=0.693$, $\Omega_{m}=0.307$, $\Omega_{b}=0.04825$ and $H_{0}=67.77$ km s^-1 which we also adopt for any calculations used in this paper. To account for all processes which operate below the numerical resolution limit, the simulations employ sub-grid physics prescriptions relying only on local hydrodynamic properties. This accounts for radiative cooling and photo-ionisation heating (Wiersma et al., 2009a), star formation (Schaye & Dalla Vecchia, 2008), stellar mass loss (Wiersma et al., 2009b) and stellar feedback (Dalla Vecchia & Schaye, 2012) and black hole growth and feedback. These parameters are then calibrated to reproduce the stellar mass function, the black hole mass - bulge mass relation and the galaxy stellar mass - size relation at $z=0.1$. A full description of the suite of simulations and calibrations can be found in Schaye et al. (2015) and Crain et al. (2015).

For the purposes of this work, star formation and black hole accretion are the most relevant sub-grid physical processes in EAGLE. Star formation is modelled according to Schaye & Dalla Vecchia (2008) as a stochastic process based on the pressure dependent Kennicutt-Schmidt law with a Chabrier (2003) IMF. Black holes are seeded in the densest gas particle in every dark matter halo greater than 1.48 $\times$ 10¹⁰ M_⊙. The black holes then grow through accretion of other gas particles. To prevent over-cooling, the energy released from accretion, assuming a radiative efficiency $\epsilon=0.1$, is not released immediately into the gas particles surrounding the black hole, but stored as an energy reservoir for a period of time given by an injection probability distribution.

EAGLE has been shown to also reproduce many properties of galaxies in the local and higher redshift universe without requiring specific tuning of the simulation parameters, such as galaxy colours, the Tully-Fisher relation, the evolution of the cosmic average SFR, the passive galaxy fraction, rotation curves and metallicities (Schaller et al., 2015; Lagos et al., 2015; Trayford et al., 2015; Furlong et al., 2015; Schaye et al., 2015). Similar to these studies, our comparison serves as a somewhat independent test of the EAGLE model, as the simulations have not been calibrated specifically to reproduce the star formation properties of AGN host galaxies.

The simulation also has various simulation box sizes, resolutions and physical models (Rosas-Guevara et al., 2016) such as higher AGN heating temperature and increased black hole viscosity. Further details of these variations can be found on the EAGLE database (McAlpine et al., 2016). The standard model referred to hereafter, incorporates AGN feedback and standard values for multiple parameters which can be found in Crain et al. (2015) as RefL0100N1504.

To obtain the properties of galaxies within EAGLE we queried the public database using the SQL interface (McAlpine et al., 2016). The key properties in low mass haloes are increasingly inaccurate due to the limitations of numerical resolution, so we only considered galaxies with stellar mass M${}_{\mathrm{Stellar}}>10^{8}$M_⊙ (i.e the number of particles $N\gtrsim 100$). To ensure uniform measurements of properties such as SFR or BHAR across the simulated galaxy sample we also applied an aperture of 30 kpc. This is the same aperture as used in previous studies and defines the maximal distance of particles directly associated with a galaxy (McAlpine et al., 2016).

Stellar mass, black hole mass and halo mass are directly available from the database. EAGLE also reproduces the shape of the cosmic evolution of the average SFR well, however Furlong et al. (2015) found that EAGLE predicts a $\simeq$ 0.2 dex lower average SFR over all epochs, possibly the result of calibration uncertainties. Therefore, following McAlpine et al. (2017), we also applied a + 0.2 dex offset to all the SFRs from the EAGLE database (see also Scholtz et al., 2018, for more discussion). We make it clear at this point that by applying an offset in the SFR measurements, we are not representing the true star formation history (SFH) of a galaxy and thereby not treating the stellar mass measurements self-consistently. We leave a fuller discussion of this, however, to Section 5.

Most of the galaxies in EAGLE, however, are inactive galaxies or AGN that would be too faint to be detected by Swift-BAT. These objects are usually lower in SFR and/or stellar mass and hence introduce a bias into our results. Stellar mass and/or X-ray luminosity matched samples may provide good comparison samples, however in order to investigate to what degree EAGLE reproduces the local universe we needed to apply the same selection to EAGLE as that of the Swift-BAT sample, namely the sensitivity or hard X-ray flux limit of the instrument.

To calculate the hard X-ray luminosity, we followed McAlpine et al. (2017), converting the BHAR ($\dot{m}$) given by the database into bolometric luminosity using the following conversion: L${}_{\mathrm{bol}}=\epsilon\dot{m}c^{2}$, assuming a radiative efficiency of $\epsilon=0.1$ and where L_bol is the bolometric luminosity of the AGN. We then converted the bolometric luminosity to hard X-ray luminosity with the prescription given by Trakhtenbrot et al. (2017), whereby L${}_{\mathrm{bol}}=$ 8.5 $\times$ L_14-195keV. We also accounted for extinction caused by obscuration of the AGN as EAGLE does not contain a prescription for the high levels of extinction that may occur near an AGN.

Ricci et al. (2015) estimated the fraction of AGN with various column densities, ranging from N_H = 10²⁰ cm^-2 - 10²⁵ cm^-2 (figure 4 of Ricci et al., 2015) and the fraction of intrinsic X-ray flux which is absorbed along the line of sight as a function of column density (figure 1 of Ricci et al., 2015). We used these distributions and randomly assigned a column density to each AGN according to the probability fraction. We then adjusted the X-ray luminosity value according to the extinction expected from this assigned column density. This process reduces the simulated AGN counts from an average of 159 sources to 132 sources ($\sim$ 15 per cent reduction), once the flux matching process is applied as described below.

To convert the X-ray luminosity to X-ray flux we calculated the distances from each galaxy to the centre of the simulation box. We then used these distances in conjunction with the converted and corrected X-ray luminosities in order to find the expected X-ray flux. Swift-BAT does not have homogeneous coverage of the sky (i.e the flux sensitivity varies for certain areas of the sky). 50% of the sky is covered down to a flux limit of 1.1 $\times$ 10^-11 erg s^-1 cm^-2, 40% of the remaining sky down to 1.5 $\times$ 10^-11 erg s^-1 cm^-2 and the remaining 10% to 1.7 $\times$ 10^-11 erg s^-1 cm^-2. We therefore used a random number generator to assign one of these flux limits to each galaxy according to the relative probability. If the estimated X-ray flux was above the randomly assigned flux limit, then the galaxy was retained. If the estimated X-ray flux was below the assigned limit, the galaxy was discarded. This process yields what we hereafter call a single realisation of the EAGLE AGN. We then repeat this process 10,000 times to fully utilise the statistical power of EAGLE (results seen in Section 4.1), hereafter the EAGLE AGN sample. Overall we find, on average, 132 EAGLE AGN per realisation compared to our 72 BAT AGN. We discuss this factor $\sim 2$ discrepancy further, as well as other bolometric to X-ray conversions, in Section 5.

The results of this process can be seen in Figure 2. The black points represent our BAT AGN sample, the blue points represent the AGN retained for a single realisation of the comparison sample (hereafter a single realisation of the EAGLE AGN) and the red points represent the galaxies in EAGLE that were discarded from the flux matching process. The solid red line shows the approximate flux limit of Swift-BAT. To investigate how much variance our process of flux correction and matching introduces, and to fully utilise the full statistical power of EAGLE, we also run the entire flux matching process 10,000 times.

One advantage of simulations is that each galaxy and its properties can be followed from the present day up to high redshift in various snapshots. This allows the exploration of the evolution of a population of galaxies over time and the properties of their progenitors. We explored the host galaxy evolution of a single realisation of our EAGLE AGN compared to the AGN selected at high redshift from EAGLE in the work of Scholtz et al. (2018). This allows us to investigate if AGN selected at the peak of cosmic AGN and star forming activity ($z\sim 2$), are inherently different, i.e. if key growth processes, as modelled by the simulation, differ at these two epochs.

The individual histories of galaxies in EAGLE can be traced over 28 snapshots from $z=0$ to $z=20$ using the merger tree system. To evaluate the entire assembly history of a galaxy at $z=0$, the paths between the present day galaxy and all of its progenitors can be traced, creating an entire tree. Alternatively, to simplify the evolutionary history, the main branch of the galaxy can be traced from the present day galaxy to the galaxy with the TopLeafID, i.e. the path of the ’main’ or most massive progenitor at each snapshot. This is a more meaningful comparison as the main branch is where the majority of the overall growth of both stellar material and SMBH mass, capturing the gross statistical trends.

A major merger may greatly affect the properties of a particular simulated galaxy between any two snapshots, i.e. a sudden and significant increase in stellar mass and/or star formation rate. This could be a potential source for biases when comparing the average galaxy properties of two different snapshots or epochs. The simulated major merger rate is low however, with half of all galaxies in the simulations with stellar mass $10^{9.5}<$ M_⊙ $<10^{10.5}$ ($10^{10.5}<$ M_⊙ $<10^{11}$, $10^{11}<$ M_⊙ $<10^{12}$) having not undergone a major merger since $z\sim 4$ ($z\sim 3$, $z\sim 1$, Qu et al., 2017). This means that when we compare the evolution of characteristic host galaxy properties of an entire AGN sample, the biases that any one major merger may contribute is likely to be washed out.

The EAGLE AGN sample from Scholtz et al. (2018) was selected to have a stellar mass M${}_{\mathrm{Stellar}}>$ 10¹⁰ M_⊙, an X-ray luminosity L_2-10keV > 10⁴³ - 10⁴⁵ and a redshift range of $1.4<z<3.6$, yielding 292 AGN (hereafter S18-EAGLE-AGN). We compared the evolution of the host galaxy properties of the S18-EAGLE-AGN to a single realisation of the EAGLE AGN; we note that although we only selected one realisation, the overall trends are not significantly dependent on different realisations.

To allow the S18-EAGLE-AGN sample to be tracked back far enough in cosmic time before the epoch of selection, we chose a redshift cut-off of $z=8$. This redshift cut-off also has the advantage that the majority of AGN host galaxies in both samples still have a sufficient number of stellar particles to avoid the statistical inaccuracies introduced by low counts. We queried the SFRs, stellar masses, black hole masses, and halo masses at each snapshot between $0<z<8$. In order to quantify the characteristic host galaxy properties at each snapshot, the median values of each quantity and in each redshift snapshot were calculated and the 16^th and 84^th percentiles of the distributions calculated. Finally, to investigate under what conditions AGN feedback plays a role in quenching star formation within EAGLE (i.e. as indicated by a decrease in star formation relative to average star forming galaxies) we also compared our EAGLE AGN sample and the S18-EAGLE-AGN to their respective star forming main sequences (hereafter SF-MS).

For the SF-MS we used the prescription from Schreiber et al. (2015). This study utilised a mixture of photometry from previous surveys and stacked galaxy imaging from $Herschel$, spanning multiple stellar mass bins (between 10^9.5 < M${}_{\mathrm{Stellar}}<10^{11.5}$ M_⊙) and redshifts (0.3 $<z<$ 5). They then selected star forming galaxies using $UVJ$ colour selection and fitted SEDs from the UV to the FIR in order to estimate galaxy properties. These estimates were then utilised to construct an average SF-MS as a function of stellar mass and redshift.

We also found that 31 galaxies within the EAGLE simulations are contained in both samples (i.e. 23% of our EAGLE AGN galaxies selected at $z=0$ are selected as AGN at $1.4<z<3.6$, according to the selection criteria of Scholtz et al., 2018). In order to further investigate this sub-population of galaxies which display significant re-current AGN episodes at both high and low redshift, we also tracked their host galaxy property evolution throughout cosmic time in the same manner as the EAGLE AGN and the S18-EAGLE-AGN.

In order to investigate how well EAGLE reproduces the star forming properties of local AGN galaxies we directly compare our BAT AGN and a single realisation of the EAGLE AGN sample. The first key tests are EAGLEs predictions of the trends between SFR and X-ray luminosity and the respective distributions: any major discrepancies indicates possible problems with the treatment of AGN in the simulation. The trends between these two variables can also be used to assess whether previously found trends are reproduced in both the observational and simulated data sets (Rosario et al., 2012; Stanley et al., 2015; Shimizu et al., 2017).

Figure 3 shows the results of X-ray luminosity versus SFR for a single realisation of our EAGLE sample. We find good overall agreement between the two samples, with no distinct outliers.

In order to reveal trends between the SFR and X-ray luminosity for the BAT AGN sample, we use Linmix, a Bayesian linear regression code originally written in IDL (Kelly, 2007). It has been shown to outperform other fitting codes, accounts for errors and also allows for upper limits and non-detections in the estimations. A python version was created and is publicly available on GitHub¹⁴¹⁴14https://github.com/jmeyers314/linmix.

The trend yielded by the fitting process for this specific single realisation is represented in Figure 3 by the solid blue line with 1$\sigma$ uncertainties given by the shaded areas. We find a slight trend of increasing SFR with X-ray Luminosity, as seen in Shimizu et al. (2017). When we run Linmix on our EAGLE AGN sample we find that EAGLE (dashed blue line) recovers the trend well and within the uncertainties (shaded regions).

The median, maximum and minimum X-ray luminosities of the EAGLE AGN lie $\approx$ 0.1 dex lower than the BAT AGN, within the observational uncertainties. The SFR distributions (Figure 3, right hand panel) are in good agreement for this single choice of realisation, as can be seen in the right hand panel, with the differences in the median within 0.1 dex and with similar ranges. Most realisations of the EAGLE AGN yield similar results; we perform a further quantitative statistical analysis in Section 4.2.

We then compare the overall distributions in SFR, stellar mass, sSFR and X-ray luminosity of the 10,000 realisations of the EAGLE AGN with the respective BAT AGN distributions, shown in Figure 4. The solid black line represents the host galaxy property distributions of the BAT AGN, the solid blue line gives the host galaxy distributions of the EAGLE AGN sample after 10,000 realisations and the shaded blue shows the distribution of a single realisation of the EAGLE AGN sample, in order to provide a comparison of the variance of a single realisation to the average. The median of each host galaxy property distribution for all three samples as well as the 16^th and 84^th percentiles are given in Table 1.

In order to statistically compare the similarity of the distributions of both a single realisation and the EAGLE AGN distribution to the BAT AGN we performed a Kolomogorov-Smirnov (K-S) test using the SciPy implementation. From the calculated supremum, a p-value is derived, providing a statistic of the probability that the two distributions are drawn from the same parent distribution. We choose a threshold of 5% (a p-value of 0.05 or below) that the two distributions of host galaxy properties could be drawn from the same parent distribution. The results of this test for both one realisation and the median, 16^th and 84^th percentiles for the full 10,000 realisations, with the associated suprema values, can also be found in Table 1.

Although the probabilities given by a K-S test for most distributions are relatively low, the probabilities in the SFR distributions is above the 5% threshold we set. The median, 16^th and 84^th percentiles of the distributions in SFR and X-ray luminosity are in agreement, with differences of $\sim$ 0.1 in the medians and similar values in the percentiles. These values are within expected observational uncertainties. We see a difference in the medians of the stellar mass of 0.2 dex. This combined with the differences in the SFR estimations, causes $\sim$ 0.4 dex differences in the sSFR distributions. This leads to very low probability yielded by a K-S test that the two distributions are drawn from the same parent distribution.

By investigating the position of AGN host galaxies relative to the SF-MS, we may reveal clues about their evolutionary stage, and thereby how an AGN may be affecting a host galaxy. By also investigating the width of the sSFR distributions, we can also investigate possible signs of AGN feedback as described in Scholtz et al. (2018), whereby implementing AGN feedback in the EAGLE simulations drives broader sSFR distributions in the galaxy population than when AGN feedback is not implemented.

Following Scholtz et al. (2018), we split the BAT AGN sample into low and high stellar mass bins using the median stellar mass and then took the median values of the sSFR distributions. For the EAGLE AGN sample, we split the samples into 16 stellar mass bins, each of 0.2 dex ranging in Log₁₀ M_Stellar $=$ 9 - 12 M_⊙ and plotted the median. These values were then compared to the SF-MS according to the prescription of Schreiber et al. (2015). We also calculate the width of the sSFR distribution of each bin, excluding those bins with only 1 galaxy. These results can be seen Figure 5. We also compare this to the results found in figure 6 of Scholtz et al. (2018), to see if the two studies yield similar results despite utilising differently selected samples of AGN at different redshifts.

In the top panel, the BAT AGN samples (black points) lie below the $z\sim 0$ SF-MS (lower red dotted line) by 0.4 and 0.3 dex respectively at constant stellar mass, outside of observational uncertainties and the expected scatter of $\sim$ 0.3 dex in the SF-MS (Schreiber et al., 2015). The EAGLE AGN sample (blue line) lies $\approx$ 0.5 dex below the SF-MS at most stellar masses with a difference of 0.6 dex and 0.4 dex lower at the same stellar masses as the two BAT AGN points. We compare this to Scholtz et al. (2018), where the observed AGN (hereafter S18 X-ray AGN, green squares) lie 0.3 dex and 0.6 dex below the $z\sim 2$ main sequence (upper red dotted line). The S18-EAGLE-AGN (orange line) are 0.4 dex below at a stellar mass of 10^10.6 M_⊙, the median stellar mass of the low mass S18 X-ray AGN sample. This shows the behaviour in the observations and simulations is similar across different reshifts and selections, however that the differences in the sSFR distributions found in Section 4.2 translate to the sSFR vs mass trends, whereby the EAGLE AGN show lower sSFRs than expected for the observations.

In the bottom panel we see the width of the sSFR distribution as a function of stellar mass. We see similar results between the simulations at different redshifts, which both display similar widths ranging between $\sim$ 0.4 dex and 0.6 dex at all stellar masses, with no obvious trend with stellar mass. We see that the widths from the observational points are within 0.1 dex of the simulated data, in good agreement.

Many modestly luminous simulated AGN in EAGLE are excluded from our analysis by our application of the BAT survey limits calculated from the central point of the simulation box (see Section 3.2 for details). Due to the relatively small size of the simulation box (100 cMpc per side), this could introduce biases due to cosmic variance. In order to check that the centre of the simulation box (or indeed any point in the simulation box) is not biased, we applied Monte-Carlo (MC) methods to draw multiple samples of AGN from the EAGLE box.

By assuming the distribution of galaxies in the total EAGLE volume is approximately homogeneous and isotropic we can model the number of galaxies per bin in proper distance as a power law. This is not unreasonable as we expect only one massive cluster within the observational volume we are using and therefore that our sample is dominated by field AGN and galaxies. To do this, we use the random number generator routine from NumPy based on a power law to generate a random distribution of distances with the amount of galaxies in the simulation box above a stellar mass of 10⁹ M_⊙ and a maximum distance matched to our BAT AGN sample, 67.8 Mpc. We then assign one of these randomised distances to each galaxy. Once the distances are assigned, assuming the same Planck cosmology applied to EAGLE, we can calculate the expected AGN X-ray flux from each galaxy. This allows us to repeat the same flux correction due to X-ray obscuration described in Section 3.2 on the new X-ray flux values and select a new sample of AGN from the simulation as if observing from a different area of the simulation box, while saving significantly on computational power.

This technique allows us to randomly draw a different sample of lower luminosity AGN while retaining a large number of higher luminosity AGN which we would expect to see from most positions in the simulation box, thereby utilising the full statistical power of EAGLE. We then compare the medians and scatter of the distributions and run a K-S test on the BAT AGN SFR (sSFR, stellar mass, X-ray luminosity) distribution and the respective newly drawn EAGLE distributions in order to recover the probability of the two distributions being drawn from the same parent sample. We run this process of random assignment, calculation, application of the BAT survey limits and comparison 10,000 times, collecting the host galaxy property distributions and the probabilities (and associated suprema values) between the observational sample and each run of the MC analyis. The median of each host galaxy property distribution, with the 16^th and 84^th percentiles given in the first 3 columns of Table 2. We also show the K-S results from a comparison of the single realisation of the EAGLE AGN via MC analyses and the median probabilities with 16^th and 84^th percentiles from the full 10,000 realisations.

Table 2 shows that differences in the median and scatter of the SFR and X-ray luminosity distributions of the two samples are unchanged and remain of the order of 0.1 dex, within the observational uncertainties. The median probabilities of the two samples being drawn from the same parent distribution are also higher at 45% and 10% respectively. The differences in the stellar mass distributions and therefore the sSFR distributions remain, however, with differences of $\sim$ 0.2 dex in the stellar mass distributions and 0.4 dex in the sSFR distributions, with low percentages of agreement from the K-S tests. These results show that comparing the BAT AGN to the EAGLE AGN in the centre of the box or drawing random samples is relatively unbiased as far as the SFRs, sSFRs and stellar masses are concerned. The X-ray luminosities, however, go from below our level of acceptance provided by a K-S test to above the threshold; we discuss this further in Section 5.1.

Figure 6 shows the results of our investigation into the evolution of various characteristic host galaxy properties between $z=0$ and $z=8$ of our EAGLE AGN sample and the S18-EAGLE-AGN sample as well as those AGN host galaxies contained in both samples; i.e., those systems that display significant AGN activity above the X-ray flux limit at both $z=0$ and high redshift. In the top left hand panel we can see that our EAGLE AGN sample (blue solid line) follows the respective SF-MS (blue dashed line) very tightly until $z\simeq 2.5$. We calculate the value of the SF-MS at all points using the redshift and median stellar mass, according to the prescription in Schreiber et al. (2015). At this point the average SFR decreases slowly relative to the main sequence until we reach the present day. Our average sample lies significantly (0.5 dex) below the SF-MS at $z=0$. Similarly with the S18-EAGLE-AGN (orange solid line), the galaxies are on the respective SF-MS (orange dashed line) until significant quenching starts to have an effect at $z\simeq 3.5$, by $z=0$ the average offset of the S18-EAGLE-AGN from the SF-MS is $\sim$ 0.6 dex.

Having established that the star formation is indeed quenching at different epochs for the two simulated AGN samples, we now try to establish the host galaxy properties influencing or following this behaviour. In the top right hand panel we have the black hole mass as a function of redshift. We see that black holes are seeded at 10⁵ M_⊙ and undergo very little initial growth. The black holes then reach a point at which they grow very rapidly and non-linearly. This phase is followed by a weakening of this rate of growth towards low redshifts. This behaviour is in agreement with previous findings in Bower et al. (2017) and McAlpine et al. (2018), whereby they propose a three stage growth scenario of black holes within EAGLE which is halo mass and redshift dependent. First, an initial stellar feedback dominated phase, where the black hole experiences little growth due to a mixture of stellar and supernovae feedback inhibiting gas accretion onto the black hole. Second, a non-linear rapid black hole growth phase, where the gravitational potential well of the galaxy overcomes the stellar and supernovae feedback and causes significant amounts of gas to migrate to the centre of the galaxy, and finally an AGN feedback regulated growth phase where the rate of growth slows.

In the bottom panels of Figure 6 we see the growth of the stellar (left hand panel) and halo mass (right hand panel). The stellar mass grows non-linearly, with a slight correlation with the black hole growth, i.e. the change in growth gradient occurs at the same points as for the black hole, although this behaviour is much weaker. The rate of growth is significantly less at low redshifts, especially for the S18-EAGLE-AGN host galaxies. The halo masses grow fairly linearly throughout redshift, however there are hints of a deceleration of this growth at very low redshifts.

The black line represents those galaxies contained in both samples as described in Section 3.3, i.e. galaxies which display significant AGN activity at both high and low redshift within the simulation (23% of all AGN galaxies selected at $z=0$). We notice that their behaviour is similar to that of the S18-EAGLE-AGN sample at high redshifts. At low redshifts, however, these objects appear to continue growing in black hole mass, stellar mass and halo mass. The median of the SFR of this sample is 0.7 dex above the median value of the S18-EAGLE-AGN, with a 0.3 dex difference in the black hole mass, 0.3 dex difference in the stellar mass and a 0.4 dex difference in the halo mass. With an average black hole mass of 10^8.6 M_⊙, stellar mass of 10^11.3 M_⊙ and halo mass of 10^13.1 M_⊙, these objects are some of the most massive in the simulation. These differences and potential drivers are discussed further in the following section.

We find reasonable agreement between the SFR distributions of our volume and flux matched BAT AGN and EAGLE AGN comparison samples, with the caveat that we applied an offset of 0.2 dex to the SFRs, following previous studies (e.g. McAlpine et al., 2016; McAlpine et al., 2018; Scholtz et al., 2018). The median and 16^th and 84^th percentiles of the SFR distributions are similar and the statistics recovered from a K-S test provide a moderate probability that the two SFR distributions are drawn from the same parent sample.

We find $\sim$ 0.1 dex differences in the median and scatter of the X-ray luminosity distributions, and a weak probability that the two distributions are drawn from the same sample, based on a K-S test. We note that although the probability yielded by a K-S test is low, the differences of 0.1 dex in the median and scatter of the distributions are within the observational uncertainties, which should be the main indicator of the similarity of the distributions.

The sSFR distributions are less well reproduced (differences of $\sim$ 0.4 dex in the median and scatter of the distributions). This is caused by EAGLE, on average, overpredicting the stellar masses of AGN host galaxies by $\gtrsim$ 0.2 dex, as seen from the distributions in Figure 4. We note, however, that observational uncertainties of the stellar masses are of the order of 0.2 dex, which could contribute to some of the differences between the two distributions. K-S tests yield very low probabilities that the observed and simulated distributions are drawn from the same parent sample for both the stellar mass and sSFR. We also note that these differences would be increased had we not applied the 0.2 dex offset in the SFR.

The results from this comparison of the AGN host galaxy property distributions are reinforced by our MC methods, which were carried out to verify that the central point of observation in the simulation box is not biased. The differences in the median and scatter of the X-ray luminosity and SFR distributions remain at $\sim$ 0.1 dex and within the observation uncertainties. The probability that the SFRs are drawn from the same parent sample are higher. The differences in the median and scatter of the sSFR and stellar mass distributions remain when we consider the results from the MC methods compared to selecting our sample from the centre of the simulation box. The probabilities of the sSFR and stellar mass distributions being drawn from the same parent sample remain low, justifying our conclusion that these results hold regardless of the point of observation in the simulation box.

The X-ray luminosity distributions, however, go from a low probability of being drawn from the sample parent sample in the centre of the simulation box ($\sim 1$ per cent) to a higher one from our MC methods ($\sim 10$ per cent, which is above our 5 per cent threshold). To test if the centre of the simulation box is biased in the X-ray luminosities we selected AGN using the same flux matching techniques as earlier, however from the eight corners of the simulation box, with a distance cut of 70 Mpc. This re-creates the Swift-BAT volume, but from a different set of independent vantage points in the simulation and thereby check if the the centre of the simulation box is biased. When we run the same techniques and tests, we find that the probabilities of the SFR, sSFR and stellar mass distributions being drawn from the same parent sample remain similar to before, however the probabilities for the luminosity distributions are 5.2 per cent. We postulate that the change in probabilities could be due to two scenarios. First that there is a slight bias in the centre of the simulation box compared to other locations within the simulation box. Second that the Milky Way may be in a biased point within the nearby universe compared to the simulation box. A thorough investigation, however, is beyond the scope of this study.

One scenario which could explain the 0.2 - 0.4 dex difference in the average stellar masses (and thereby tensions of 0.4 dex in the average sSFRs) found in our study has been postulated in a recent study carried out by Leja et al. (2020). They compare stellar mass estimates using parametric SFHs to those estimated using non-parametric SFHS. They find that parametric SFHs may underestimate empirically derived stellar masses in observations on the order of 0.1 - 1 dex, which could account for the $\gtrsim$ 0.2 dex differences in median and scatter of the observed and simulated stellar masses distributions. This result also affects the observed stellar mass function, so future work in this field is also likely to have implications to the modelling of simulations, namely that if stellar masses are systematically underestimated then the stellar mass function needs revising. As this is one of the most commonly used calibrators for simulations, this would require changes in the sub-grid physics.

An alternative solution to the tensions in the stellar masses can be provided by McAlpine et al. (2018). This study found that the rapid growth phase of SMBHs in EAGLE was triggered by increasing fractions of major mergers towards low redshifts. Major mergers can deliver significant amounts of stellar mass to the central galaxy, thereby biasing the comparison sample towards higher stellar masses. If there are indeed differences in the fueling mechanisms between the observations and simulations (such as if AGN fuelling is triggered less frequently by mergers in the observed universe), this could be a source of the biases in the stellar masses in simulated AGN host galaxies. An in depth comparison of AGN fueling mechanisms between simulations and observations, however, is beyond the scope of this work.

As stated in Section 3.2, once we apply our X-ray luminosity correction, we find that EAGLE predicts 132 AGN per sample on average over 10,000 realisations, compared to the 72 in our BAT AGN sample, an overestimation in the number counts. We note that there are assumptions we have made that may impact these results. Our conversion of the fraction of luminosity from an AGN (L_bol) radiated in the Ultra hard X-ray band is taken from Trakhtenbrot et al. (2017) (a value of 8.5). When we assume a bolometric conversion such as that applied by Koss et al. (2017) (with a value of 8), we find an average of 142 AGN over 10,000 realisations. Conversely, when we assume a bolometric correction factor of 10 such as that in Scholtz et al. (2018), albeit for a different X-ray energy band, an average of 120 AGN is yielded over 10,000 realisations. This shows that the AGN counts are extremely sensitive to the bolometric to ultra hard X-ray conversion factor used.

We also note that the X-ray luminosity correction from Ricci et al. (2015) accounts for fractions of AGN at various opening angles, which is not accounted for in the simulations. Finally, EAGLE itself could be overpredicting the number count of AGN, similar to overpredicting the stellar mass of the host galaxy. The source of the differences would need to be investigated in further work.

Although the different bolometric conversions yield varying counts of AGN on average in the sample selection, when implemented and the simulated samples compared to the BAT AGN sample, we find very little difference in the results. The difference in the percentages yielded by a K-S test for the probability of the two samples being drawn from the same parent distribution change by less than 3% in the SFR distributions, 0.5% for the X-ray luminosity, 0.1% for the sSFR distributions and 0.3% for the stellar mass, in both the original method and in the MC methods. Similarly the medians of SFR, sSFR and stellar mass distributions change by $<$ 0.1 dex, not significantly different to the results from our original comparison.

From our use of the merger trees for AGN samples selected at different epochs, our EAGLE AGN sample and the S18-EAGLE-AGN sample, we see how the selection criteria can affect the expected host galaxy properties of AGN and how similarly their host galaxy properties can evolve. Downsizing, whereby less massive objects undergo physical processes later in the history of the universe (Brandt & Hasinger, 2005), is clearly displayed in the SFRs, black hole masses, stellar masses and halo masses, as shown in Figure 6.

The EAGLE AGN (which have on average smaller stellar or halo masses across all epochs) leave the SF-MS later than the S18-EAGLE-AGN ($z=3.5$ compared to $z=2.5$), in other words their star formation is quenched at later times. The black hole masses in the EAGLE AGN sample also grow at later epochs than the S18-EAGLE-AGN, especially their phase of most rapid growth. The point at which they start this rapid growth phase also approximately coincides with the point at which both populations leave the SF-MS. This shows the connection between the time at which the AGN injects proportionally the most energy and the quenching of star formation, as expected. We also see this growth slow down at later epochs, confirming the overall behaviour of black hole growth seen by McAlpine et al. (2018).

We also see that the star formation starts deviating from the main sequence for both our EAGLE AGN sample and the S18-EAGLE-AGN when the halo masses reach a critical point of M_Halo $\approx$ 10¹² M_⊙. Studies by Bower et al. (2017) and McAlpine et al. (2017) explain that the halo is massive enough that the stellar and supernovae feedback is no longer efficient enough to eject gas from the galaxy. At this point gas starts flowing into the centre of the galaxy, accreting around the black hole, thereby growing it. The black hole at this point then starts injecting energy back into the host galaxy, thereby quenching star formation.

AGN galaxies which are contained in both samples and therefore display significant AGN activity at both high and low redshift are of interest, as this phenomena cannot be observed. We find that 23% of all AGN galaxies selected at $z=0$ displayed significant activity at high redshift. These objects tend to follow similar behaviour to the S18-EAGLE-AGN sample, however at low redshift they continue to grow in stellar, halo and black hole mass.

One problem with some hydrodynamic simulations is the inability to fully reproduce the quenched fraction of massive galaxies, driven by higher gas fractions than in the observations (see e.g. Nelson et al., 2015). A possible scenario this may lead to is a sub-sample of massive galaxies at high redshift able to undergo an AGN fuelling phase and then a low redshift continuing to be fuelled by gas, leading to higher than expected levels of star formation and accretion onto the SMBH. This then leads to higher BH, stellar and halo masses and an AGN phase at $z\sim 0$. We cannot fully rule out, however, that an AGN in the observed universe would display this behaviour due to the duty cycle or triggering from mergers.