BlueTides simulation: establishing black hole-galaxy relations at high-redshift

The scaling relations between the mass of supermassive black holes ($M_{\bullet}$) and host galaxy properties (stellar mass, $M_{\star}$, and velocity dispersion, $\sigma$), provide a link between the growth of black holes (BHs) and that of their hosts. Here we investigate if and how the BH-galaxy relations are established in the high-$z$ universe using \textsc{BlueTides}, a high-resolution large volume cosmological hydrodynamic simulation. We find the $M_{\bullet}-M_{\star}$ and $M_{\bullet}-\sigma$ relations at $z=8$: $\log_{10}(M_{\bullet}) = 8.25 + 1.10 \ \log_{10}(M_{\star}/10^{11}M_{\odot})$ and $\log_{10}(M_{\bullet}) = 8.35 + 5.31 \ \log_{10}(\sigma/200kms^{-1})$ at $z=8$, both fully consistent with the local measurements. The slope of the $M_{\bullet}-\sigma$ relation is slightly steeper for high star formation rate and $M_{\star}$ galaxies while it remains unchanged as a function of Eddington accretion rate onto the BH. The intrinsic scatter in $M_{\bullet}-\sigma$ relation in all cases ($\epsilon \sim 0.4$) is larger at these redshifts than inferred from observations and larger than in $M_{\bullet}-M_{\star}$ relation ($\epsilon \sim 0.14$). We find the gas-to-stellar ratio $f=M_{\rm gas}/M_{\star}$ in the host (which can be very high at these redshifts) to have the most significant impact setting the intrinsic scatter of $M_{\bullet}-\sigma$. The scatter is significantly reduced when galaxies with high gas fractions ($\epsilon = 0.28$ as $f<10$) are excluded (making the sample more comparable to low-$z$ galaxies); these systems have the largest star formation rates and black hole accretion rates, indicating that these fast-growing systems are still moving toward the relation at these high redshifts. Examining the evolution (from $z=10$ to 8) of high mass black holes in $M_{\bullet}-\sigma$ plane confirms this trend.


INTRODUCTION
Over the last three decades, scaling relations between mass of supermassive black holes (SMBHs) and several stellar properties of their host galaxies such as bulge stellar mass and bulge velocity dispersion (Magorrian et al. 1998;Häring & Rix 2004;Gebhardt et al. 2000;Tremaine et al. 2002;Gültekin et al. 2009;Kormendy & Ho 2013;McConnell & Ma 2013;Reines & Volonteri 2015) have been discovered and measured for galaxies with black holes (BHs) and active galactic nuclei (AGN) from z = 0 and up to z ∼ 2 (using different techniques).
Many theoretical models have been developed to understand the origin of these relations. Several cosmological simulations that follow the formation, growth of BHs and their host galaxies have successfully reproduced the ⋆ E-mail: kuanweih@andrew.cmu.edu scaling relations at low-z; these include recent simulations such as the Illustris simulation (Vogelsberger et al. 2014;Sijacki et al. 2015), the Magneticum Pathfinder SPH simulation (Steinborn et al. 2015), the Evolution and Assembly of GaLaxies and their Environment (EAGLE) suite of SPH simulation (Schaye et al. 2015), and the MassiveBlackII (MBII) simulation (Khandai et al. 2015;DeGraf et al. 2015). Thus, the scaling relations from observations and simulations agree with each other at low-z, linking the growth of SMBHs to the growth of their hosts via AGN feedback.
A popular way to interpret the scaling relations is by invoking AGN feedback. Many models (and simulations) show that the SMBHs regulate their own growth with their hosts by coupling a fraction of their released energy back to the surrounding gas (Di Matteo et al. 2005). The BHs grow only until sufficient energy is released to unbind the gas from the local galaxy potential (Silk & Rees 1998;King 2003;Springel 2005;Bower et al. 2006;Croton et al. 2006;Di Matteo et al. 2008;Ciotti et al. 2009;Fanidakis et al. 2011). However, there are also models which have been proposed to explain the scaling relations without invoking the foregoing coupled feedback mechanism. For instance, it has also been shown that dry mergers can potentially drive BHs and their hosts towards a mean relation (Peng 2007;Hirschmann et al. 2010;Jahnke & Macciò 2011). Regardless, studying the scaling relations in both observation and simulation is essential for understanding the coupled growth of galaxies and BHs across cosmic history.
An important related question is when the scaling relations are established, and if they still persist at higher redshifts when the first massive BHs form (z > 6). To understand this, galaxies with AGN play a key role in observations (Bennert et al. 2010;Merloni et al. 2010;Kormendy & Ho 2013). A strong direct constraint on the high-redshift evolution of SMBHs comes from the luminous quasars at z ∼ 6 in SDSS (Fan et al. 2006;Jiang et al. 2009;Mortlock et al. 2011). Most recently, the earliest quasar is discovered at z = 7.5 (Bañados et al. 2017) in ALLWISE, UKIDSS, and DE-CaLS. However, it is still not established as to whether these objects follow the local BH-galaxy relations and whether there is a redshift evolution, because of the systematic uncertainties (Woo et al. 2006) and selection effects (Lauer et al. 2007;Treu et al. 2007;Schulze & Wisotzki 2011. At high-z, AGN is our only proxy for studying the BH mass assembly. The luminosity functions (LFs) of AGN however remain uncertain. For example, BH mass function at z = 6 has been inferred from optical AGN LFs in Willott et al. (2010). On the other hand, several works (Wang et al. 2010;Volonteri & Stark 2011;Fiore et al. 2012;Volonteri & Reines 2016) have argued that there are large populations of obscured, accreting BHs at high-z. For instance, the BH mass density at z = 6 from X-ray observations of AGN has been shown to be greater than that inferred from optical quasars by an order of magnitude or more (Treister et al. 2011;Willott 2011). A luminosity dependent correction for the obscured fraction is proposed in Ueda et al. (2014) and the obscured fraction tends to increase with redshift up to z ∼ 4 (Merloni et al. 2014;Vito et al. 2014;Buchner et al. 2015;Vito et al. 2018). Regardless of the exact amount of the obscured AGNs, it is certain that a fraction of AGNs is obscured, and therefore missed by observations. As a result, quantities such as the BH mass function, BH mass density, and BH accretion rate density are still uncertain at high-z.
Here, we use the BlueTides simulation  to make predictions for both the global BH mass properties and the scaling relations (M • − M ⋆ and M • − σ relations) from z = 8 to z = 10. BlueTides is a largescale and high-resolution cosmological hydrodynamic simulation with 2 × 7040 3 particles in a box of 400h −1 M pc on a side, which includes improved prescriptions for star formation, BH accretion, and associated feedback processes. With such high resolution and large volume, we are able to study the scaling relations and the global properties of BH mass at high-z for the first time. So far, various quantities measured in BlueTides have been shown to be in good agreement with all current observational constraints in the high-z universe such as UV luminosity functions Waters et al. 2016a,b;Wilkins et al. 2017), the first  Di Matteo et al. 2017;Tenneti et al. 2018), the Lyman continuum photon production efficiency , galaxy stellar mass functions (Wilkins et al. 2018), and angular clustering amplitude (Bhowmick et al. 2017). The paper is organized as follows. In Section 2, we briefly describe the BlueTides simulation and several physics implementations. In Section 3, we report BH mass properties: mass function, mass density, and accretion rate. In Section 4, we demonstrate the scaling relations between M • and M ⋆ , and σ. In Section 5, we study the selection effects for several galaxy properties on the scaling relations. In Section 6, we investigate the assembly history of how BHs evolove on M • − M ⋆ and M • − σ planes. In Section 7, we summarize the conclusion of the paper.

BlueTides hydrodynamic simulation
The BlueTides simulation has been carried out using the Smoothed Particle Hydrodynamics code MP-Gadget on the Blue Waters system at the National Center for Supercomputing Applications. The hydrodynamics solver in MP-Gadget adopts the new pressure-entropy formulation of smoothed particle hydrodynamics (Hopkins 2013). This formulation avoids non-physical surface tensions across density discontinuities. BlueTides contains 2 × 7040 3 particles in a cube of 400h −1 M pc on a side with a gravitational smoothing length ǫ = 1.5h −1 k pc. The dark matter and gas particles masses are M DM = 1.2 × 10 7 h −1 M ⊙ and M gas = 2.36 × 10 6 h −1 M ⊙ , respectively. The cosmological parameters used were based on the Wilkinson Microwave Anisotropy Probe nine years data (Hinshaw et al. 2013) (see Table 1 for a brief summary of the parameters). With an unprecedented volume and resolution, BlueTides runs from z = 99 to z = 8. BlueTides contains approximately 200 million star-forming galaxies, 160000 of which have stellar mass > 10 8 M ⊙ , and 50 thouthand BHs, 14000 of which have BH mass > 10 6 M ⊙ (the most massive BH's mass ∼ 4 × 10 8 M ⊙ ). A full description of BlueTides simulation can be found in Feng et al. (2016).

Sub-grid physics and BH model
A number of physical processes are modeled via sub-grid prescriptions for galaxy formation in BlueTides. Below we list the main features of the sub-grid models: • Star formation based on a multiphase star formation model (Springel & Hernquist 2003)   • Gas cooling through radiative processes (Katz et al. 1996) and metal cooling (Vogelsberger et al. 2014).
• Formation of molecular hydrogen and its effects on star formation (Krumholz & Gnedin 2011).
• Type II supernovae wind feedback (the model used in Illustris ).
• A model of 'patchy' reionization (Battaglia et al. 2013) yielding a mean reionization redshift z ∼ 10 (Hinshaw et al. 2013), and incorporating the UV background estimated by Faucher-Giguère et al. (2009); • Black growth and AGN feedback. BHs grow in mass by gas accretion and by merging with other BHs.
We model BH growth and AGN feedback in the same way as in the MassiveBlack I & II simulations, using the SMBH model developed in Di Matteo et al. (2005) with modifications consistent with Illustris. BHs are seeded with an initial seed mass of M •,seed = 5 × 10 5 h −1 M ⊙ (commensurate with the resolution of the simulation) in halos more massive than 5 × 10 10 h −1 M ⊙ while their feedback energy is deposited in a sphere of twice the radius of the SPH smoothing kernel of the black hole. Gas accretion proceeds −3/2 according to Hoyle & Lyttleton (1939); Bondi & Hoyle (1944); Bondi (1952), where ρ and c s are the density and sound speed of the gas respectively, α is a dimensionless parameter, and v is the velocity of the BH relative to the gas. We allow for super-Eddington accretion but limit the accretion rate to three times of the Eddington rate: where m p is the proton mass, σ T is the Thomson cross-section, and η is the radiative efficiency. Fixing η to 0.1 according to Shakura & Sunyaev (1973) throughout the simulation, the BH is assumed to radiate with a bolometric luminosity (L B ) proportional to the accretion rate (

Kinematic decomposition
As σ or M ⋆ in observational studies of M • − σ or M • − M ⋆ relations are often measured from the bulge component of galaxies, we perform a kinematic decomposition for the stellar particles of the galaxies in BlueTides as in Feng et al. (2015). This allows us to determine which stars are on planar circular orbits and which are associated with a bulge, in each galaxy (Vogelsberger et al. 2014;Tenneti et al. 2016), providing kinematically classified disks and bulges, and a disk to total (D/T) ratio for our galaxies. We perform this analysis following Abadi et al. (2003): a circularity parameter is defined for every star particle as κ = j z / j (E), where j z is the specific angular momentum around a selected z-axis and j (E) is the possible maximum specific angular momentum of the star with the specific binding energy E. The star particle with κ > 0.7 is identified as a disk component according to Vogelsberger et al. (2014) and Tenneti et al. (2016). Thus, the D/T ratio for the stellar component of each galaxy is obtained, allowing us to calculate the bulge stellar mass and the bulge velocity dispersion for our galaxies. In Figure 1, we show the relation between B/T = 1 − D/T (bulge to total ratio) and total stellar mass (M ⋆ ) color coded according to number of star particles for each galaxy. According to the standard assumption D/T < 0.3 is considered a bulge dominated galaxies Tenneti et al. 2016). For galaxies with M ⋆ > 10 9 M ⊙ , the number of star particles is higher than 1000. We require this minimum number of star particles to have a reliable kinematic decomposion. For this reason, for the rest of the analysis we will only consider objects with M ⋆ > 10 9 M ⊙ .

THE GLOBAL PROPERTY OF BH MASS
We begin by investigating the global properties of BH mass (M • ) from z = 8 to z = 12 in BlueTides. We choose a BH population with M • > 1.5 × 10 6 M ⊙ which is roughly twice the BH seed mass ( M •,seed = 7.2 × 10 5 M ⊙ ), in order to minimize any possible influence of the seeding prescription on our analysis.

BH mass function and bolometric luminosity
We first look at the bolometric luminosity (L B ) of BH population in BlueTides at z = 8 in the left panel in Figure 2. The brown dashed and dotted lines are X-ray luminosity L X = 10 42.5 and 10 43 erg/s, which are calculated by the bolometric correction in Marconi et al. (2004). These two values will be used as thresholds when studying other global properties of BH mass in this section. Statistically, there are more than 76 and 16 percent of our BHs with L X > 10 42.5 and 10 43 erg/s respectively. This indicates that the global quantities of BH mass is sensitive to L X when measured from X-ray survey in observation. In addition, L B with the mean Eddington ratio in our BH population (λ Edd = 0.29) is shown (the green solid line).
The right panel in Figure 2 shows BH mass functions (BHMFs) in BlueTides from z = 8 to z = 12 (the solid curves), as well as the ones with thresholds L X > 10 42.5 and 10 43 erg/s (the dashed and dotted curves respectively). We also show the BHMFs inferred from optical quasars at z = 6 in Willott et al. (2010) (W10 hereafter; the purple dotted curve) and the theoretical prediction combined with observed Lyman break galaxy population (Stark et al. 2009;Volonteri & Stark 2011) (the red squares). The slope of BHMFs in BlueTides are generally steeper than the one in W10 but similar to theoretical predictions (see Volonteri & Stark (2011) for more detail and comparison), particularly at the low mass end (which is currently unconstrained at these redshifts). In addition, the normalization of the BHMFs in BlueTides suggests that there is a larger BH population than those which have been observed from optical quasars, consistent with the claim that there is a large population of obscured accreting BHs at high-z.

BH mass density and stellar mass density
The left panel in Figure 3 shows BH and galaxy stellar mass density in BlueTides from z = 8 to z = 12 with observations from z = 4 to z = 7. For stellar mass density (SMD), BlueTides (the solid blue curve) agrees with the trend from Grazian et al. (2015) (the orange stars). Galaxies with M ⋆ > 10 8 M ⊙ are selected for both cases. For BH mass density (BHMD), we report the results with two different M • thresholds: M • > 1.5 × 10 6 M ⊙ ( double of M •,seed in BlueTides) and M • > 10 7 M ⊙ ) and with L X > 10 43er g/s (the red solid, olive dash-dotted, and red dash curves respectively) to show the influence from M • and L X thresholds. Current upper limits from X-ray observation are also presented: Salvaterra et al. (2012) (cosmic X-ray background (XRB) ) and Treister et al. (2013) (the brown diamonds and gray triangles respectively). BHMD in BlueTides complies with those upper limits (at least at z = 8), pointing that our BH mass function is just steeper than the one in W10 but still within the upper limits from X-ray observation. In addition, results in Volonteri & Reines (2016) (the blue shade) are also included to support our BHMF, arguing that the integrated BH density depends on the M • − M ⋆ relation. It has been discussed that the stellar mass density exceeds the BH mass density roughly by a factor of 10 3 for low-z. To understand the ratio of these two quantities at higher redshifts, we show the ratio by normalizing SMD to BHMD in the left panel in Figure 3 (the green solid curve). Overall, SMD grows more rapidly than the BHMD at early times. Parameterizing the ratio by an evolutionary factor (1 + z) α , we find that SM D/BH M D = 1.3(1 + z) 3.1 .

BH accretion rate density and SFR density
After BHMD and SMD, we investigate their assembly rate: the BH accretion rate density and the SFR density. The right panel in Figure 3 shows the BH accretion rate density (BHAD) again with M • and L X thresholds M • > 1.5×10 6 M ⊙ , M • > 10 7 M ⊙ , and L X > 10 43 erg/s (the red solid, olive dash-dotted, and red dash curves respectively) and the SFR density (SFRD; the blue solid curve) in BlueTides from z = 8 to z = 12. For SFRD, we show observational result in Vito et al. (2016) (the orange stars), and on the other hand for BHAD, we not only show results from observation (Vito et al. (2018); the green dots) but also from simulation (Sijacki et al. (2015); the dotted purple curve) because it has been noticed that the prediction from simulation tends to be higher than current observation by more than an order of magnitude. It is possibly due to the difficulty of observing AGNs from deep X-ray surveys. Similar to Section 3.2, we compare these two quantities by their ratio via normalizing SFRD to BHAD (the green curve). The ratio of the SFRD and the BHAD increases as z increases and the order of which (ranging from 10 3 to 10 4 ) is close to the order of the ratio of the BHMD and the SMD in our simulation. Again, we fit the ratio as a function of (1 + z) α : SFRD/BH ARD = 3.4(1 + z) 2.7 .

Measuring M ⋆ and σ
The scaling relations between BH mass and their host galaxy properties have been measured by total stellar mass (M ⋆,total )  In simulations, total or half stellar mass is available as a proxy for M ⋆,bulge . A proxy for σ often used is the velocity dispersion within half-light (or mass) radius (as for example in Sijacki et al. (2015) and DeGraf et al. (2015)). With the dynamical disk-bulge decomposition (see Section 2.3) for the stellar components of galaxies in BlueTides, we directly have M ⋆,bulge and σ bulge in our galaxies. The top panel in Figure 4 shows the comparison between the M ⋆,bulge and M ⋆,total color coded according to the number of galaxies. We find that about 82% of objects are bulge-dominated with D/T < 0.3. The bottom panel in Figure 4 shows the comparison between σ bulge and σ hm color coded according to the number of galaxies. Again, there is certainly a strong correlation between the two but an increased scatter above the one-to-one relation for a small number of objects, for which we have larger values of σ bulge for a given σ hm . In particular, we find that over 94% and 97% of our galaxies have the difference between σ bulge and σ hm less than 10% and 20% respectively. We shall see in the next section how the detailed dynamical decomposition and the resulting σ bulge and M ⋆,bulge impact the scaling relations.    (1), the total number of data points N , and the standard deviation of residuals ǫ of the scaling relations at each redshift.
where M • is in units of M ⊙ , and X is M ⋆ /10 11 M ⊙ or σ/200kms −1 . The fitting coefficients (normalization α and slope β) are summarized in Table 2, including the total num-ber of data points N and the standard deviation of the residuals ǫ. The left panels in Figure 5 show the M • − M ⋆ relations. The red and blue lines show the best-fitting relation with M ⋆,total and M ⋆,bulge respectively while the gray lines show the observations: Häring & Rix (2004), McConnell & Ma (2013) and Kormendy & Ho (2013) with bulge stellar mass and elliptical samples while Volonteri & Reines (2016) with total stellar mass. Our simulation provides the M • − M ⋆ relation in the form of log 10 (M • ) = 8.25+ 1.10 log 10 (M ⋆ /10 11 M ⊙ ) with M ⋆,total at z = 8, suggesting that the slopes are consistent with the observations but the normalizations are lower than most observations except for the one in Häring & Rix (2004). Both α and β with M ⋆,total (the red lines) are lower than the ones with M ⋆,bulge (the blue lines) across all three redshifts, and both get steeper as z is higher. The standard deviation of the residuals (ǫ) is shown as the shaded area.
The right three panels in Figure 5 show the M • − σ relations. The red and green lines show the best-fitting relation with σ hm and σ bulge respectively while the gray lines show the observations in McConnell & Ma (2013). Our simulation provides the M • − σ relation with σ hm as log 10 (M • ) = 8.35 + 5.31 log 10 (σ/200kms −1 ) at z = 8, which is consistent with the results of McConnell & Ma (2013). We note that both α and β using σ hm (the red lines) are lower by ∼ 3% and ∼ 10%, respectively than the ones with σ bulge (the blue lines) across all three redshifts, and both get steeper with increasing z . Moreover, M • − σ relations with σ bulge are higher than local measurements. ǫ is shown as the shaded area and, more importantly, M • −σ relation shows a larger scatter than the M • − M ⋆ relation (ǫ ∼ 0.4 and ǫ ∼ 0.1 respectively) in our simulations. We will examine this in Section 5.
For most observational results, the scaling relations are established with bulge-dominated galaxies. Here, we report how the relations change in our simulation with bulgedominated galaxies (D/T < 0.3). In the top two panels in Figure 6, we show both relations color coded according to the D/T ratio. In the bottom panels, we show α of M • −M ⋆ relations and β and ǫ of M • −σ relations as functions of limiting D/T. We find that the relations hardly change even with different D/T, even for the bulge-dominated regime D/T < 0.3.

THE SLOPE AND SCATTER IN THE SCALING RELATIONS
We have shown that the high-z relation is consistent with the locally measured ones (both in slope and normalization). However here we wish to investigate possible selection effects and/or physical parameters that may affect the slope and scatter in the scaling relations ( Figure 5 and Table 2).
In particular, we have found that there is a more significant scatter in the M • − σ relation (larger than in local measurements) than the M • − M ⋆ relation. The relatively large scatter in the M • − σ relation appears to be due to a significant amount of objects that lie below the main relation: galaxies with relative high σ compared to their relative low M • . Note that M ⋆ denotes M ⋆,total and σ denotes σ hm hereafter unless stated otherwise. effects at these high redshifts (which typically bias samples toward high SFR and M ⋆ objects) are likely to play an important role and lead to biased interpretations of evolutionary effects in these relations when compared to those seen at low-z. The scatter also tends to increase particularly when high SFR or high M ⋆ samples are selected (which populate the σ ≥ 200 km/s range).

5.2
The gas fraction: f = M gas /M ⋆ In Section 5.1, we have examined the dependence of the scaling relation on the various galaxy or AGN parameters; none of them helps to explain the relatively large intrinsic scatter in the M • − σ relation. Here, we test for the effects due to the large range and high value of gas fraction in the high-z galaxies and how that may affect the M • − σ relation. We use the gas-to-stellar ratio ( f ) to measure the gas fraction in our galaxies in BlueTides. Figure 8 shows the M • − σ relation at z = 8 color coded according to f . We find that the gas fraction in galaxies f has a significant impact on the scatter of the M • −σ relation. In particular, the scatter is significantly reduced for galaxies with a smaller value of f . We find that ǫ decreases significantly from 0.36 (all galaxies) to 0.28 (with f < 10) while the slope of the relation (β) decreases from 5.31 (all galaxies) to 4.64 (with f < 10). We further illustrate this trend of decreasing of ǫ and β with different limiting f , in the bottom panel in Figure 8. Lower gas fractions are indeed more representative of local galaxies, which have been used to measure the M • − σ.
The decrease of both ǫ and β with the lower limiting f implies that objects with higher σ but lower M • are those that have higher f . To look into the relations between f , M • , and σ, we show Figure 9, color coded according to σ and M • respectively. The top panel indicates a trend between a decreasing f at increasing M • . Gas fractions of galaxies decrease as BH masses are high/reach the relation, indicating that BH growth is quenched/self-regulated due to AGN feedback. The galaxies with higher σ but lower M • are indeed those that have larger f (the light blue area in the bottom panel), which results in that the M • − σ relation is more scattered if there is no limiting f applied. These are relatively sizeable galaxies where significant BH growth is still occurring (see also Figure 10 and objects are still moving toward the relation (feedback has not yet saturated the BH growth). Such an actively growing BH population is indeed rare among the sample of quiescent local BHs.

THE ASSEMBLY HISTORY
To further investigate how the M • − M ⋆ and M • − σ relations are established we trace the evolution of ∼ 200 black holes (from z ∼ 10 to z ∼ 8) and their hosts. In Figure 10, we show sample tracks on the M • − M ⋆ and M • − σ plane from z = 10 to z = 8. The solid lines shown in orange and blue show the average track in the evolution for higher mass (with M • 5 × 10 7 M ⊙ ) and lower mass (with M • 5 × 10 6 M ⊙ ) respectively. The tracks in the M • −σ plane suggest a slightly steeper growth in the low mass galaxies compared to the high mass galaxies. This likely indicates a fast growth of the black hole mass at approximately fixed values of σ up to the point where the galaxies reach the average relation.
To characterize the overall evolution in these planes, we parameterize the growth of M • , M ⋆ , and σ as where the exponents are γ • = −9.1 ± 1.0, γ ⋆ = −8.1 ± 2.2, and γ σ = −1.6 ± 0.8, (the error bars are standard deviation errors). Note that γ • /γ ⋆ ∼ 1.1 and γ • /γ σ ∼ 5.7, which is consistent with the slope of the scaling relation shown in Table 2 (if we use Eqs. (2) and eliminate (1 + z), then . This suggests that, on average, the redshift evolution of these black holes traces the overall scaling relation (at a given redshift), with commensurate growth in black hole mass and stellar mass. We now take a step further and look into two distinct mass regimes of the above sample of 200 BHs: 1) A low mass range of BHs with M • < 5 × 10 6 M ⊙ , and 2) A high mass range of BHs with M • > 1 × 10 7 M ⊙ . We note that, as discussed in the previous section, the lower mass subsample has more scatter in the M − σ relation (see Section 4). We find that γ • = −5.7±3.5, γ σ = −0.72±0.49 and γ • = −10±2.6, γ σ = −1.8 ± 0.5 for low mass and high mass subsamples respectively; therefore, the higher mass BHs have steeper increase in both M • and σ, as is clearly illustrated in Figure 10. Interestingly, we also find γ • /γ σ ∼ 7.9, which is higher than the slope of the overall fit of the M − σ relation in Section 4, suggesting that indeed these lower mass BHs tend to grow their black holes faster than their velocity dispersion and saturate their growth as they move closer to the mean relation. This behavior can then potentially lead to a decrease in the scatter in the relation at the low redshifts, where most black holes and galaxies have low gas fractions and have quenched their growth and star formation (particularly in the bulge dominated samples where the relations are typically measured). For the high mass sample, γ • /γ σ ∼ 5.55 which is consistent with the slope of the overall fit of the M − σ relation, so that high mass objects continue to move along the relation.

CONCLUSIONS
We investigate the global properties of supermassive black holes at high redshifts (z ∼ 8, 9, 10), which include scaling relations w.r.t properties of their host galaxies (σ, M ⋆ ) and their redshift evolution using the BlueTides simulation. The bolometric luminosity of BHs span more than two orders of magnitude around a mean of 0.3 of the Eddington luminosity. The BH mass functions in our simulation tend to have steeper slopes compared to the one inferred at z = 6 measured from optical quasars. While this may be due to obscuration, we find that it is consistent with the large range of luminosities spanned and the flux cuts implied by observations. We have also shown that the BH mass density and BH accretion rate are broadly consistent with current observational constraints at the highest redshifts (z ∼ 7).
The scaling relations, M • − M ⋆ , M • − σ predicted by BlueTides reveal that correlations between the growth of black holes and their host galaxies persist at high-z (z = 8 to z = 10), with the slopes and normalizations consistent with published relations at low-z. For the scatter, we find that the M • − σ relation has a significantly higher scatter compared to current measurements as well as the M • − M ⋆ relation. We further show that this large scatter can be primarily attributed to the gas-to-stellar ratio ( f = M gas /M ⋆ ), wherein we observe a significant decrease in the scatter (ǫ = 0.36 to ǫ = 0.28) upon exclusion of galaxies with f = M gas /M ⋆ > 10. Such high gas fraction systems have the largest star formation rates and black hole accretion rates indicating that these fast-growing systems have not yet converged to the relation.
We also find that the assembly history of the evolution of BHs on M • − M ⋆ and M • − σ planes is, on an average, consistent with the corresponding scaling relations; in other words, the average trajectory of the evolution of BHs traces the mean scaling relations well.