The early growth of supermassive black holes in cosmological hydrodynamic simulations with constrained Gaussian realizations

TL;DR

This study uses constrained Gaussian realizations in cosmological hydrodynamic simulations to explore how initial conditions and seed masses influence the early growth of supermassive black holes, shedding light on their formation mechanisms.

Contribution

It introduces a novel constrained simulation approach to reproduce early SMBH growth and compares different seed mass scenarios to understand their impact on black hole evolution.

Findings

Low-tidal field and high-density peaks accelerate SMBH growth.

Small seed masses lead to more frequent early BH mergers.

All SMBHs converge to ~10^9 M_sun by z=6 except the smallest seed case.

Abstract

The paper examines the early growth of supermassive black holes (SMBHs) in cosmological hydrodynamic simulations with different BH seeding scenarios. Employing the constrained Gaussian realization, we reconstruct the initial conditions in the large-volume BlueTides simulation and run them to $z=6$ to cross-validate that the method reproduces the first quasars and their environments. Our constrained simulations in a volume of $(15\, h^{-1}{\rm Mpc})^3$ successfully recover the evolution of large-scale structure and the stellar and BH masses in the vicinity of a $\sim10^{12}\, M_{\odot}$ halo which we identified in BlueTides at $z\sim7$ hosting a $\sim10^9\, M_{\odot}$ SMBH. Among our constrained simulations, only the ones with a low-tidal field and high-density peak in the initial conditions induce the fastest BH growth required to explain the $z>6$ quasars. We run two sets of simulations with different BH seed masses of $5\times10^3$, $5\times10^4$, and $5\times10^5\, h^{-1}M_{\odot}$, (a) with the same ratio of halo to BH seed mass and (b) with the same halo threshold mass. At $z=6$, all the SMBHs converge in mass to $\sim10^9\, M_{\odot}$ except for the one with the smallest seed in (b) undergoing critical BH growth and reaching $10^8$ -- $10^9\, M_{\odot}$, albeit with most of the growth in (b) delayed compared to set (a). The finding of eight BH mergers in the small-seed scenario (four with masses $10^4$ -- $10^6\, M_{\odot}$ at $z>12$), six in the intermediate-seed scenario, and zero in the large-seed scenario suggests that the vast BHs in the small-seed scenario merge frequently during the early phases of the growth of SMBHs. The increased BH merger rate for the low-mass BH seed and halo threshold scenario provides an exciting prospect for discriminating BH formation mechanisms with the advent of multi-messenger astrophysics and next-generation gravitational wave facilities.

Tables

Table 1. Table 1 : Parameters adopted in our simulations.

$h$	0.697	${\Omega }_{\mathrm{matter}}$	0.2814	$M_{\mathrm{DM}}$	$1.7\times {10}^7{M}_{☉}$
${\sigma }_8$	0.820	${\Omega }_{\mathrm{baryon}}$	0.0464	$M_{\mathrm{gas}}$	$3.4\times {10}^6{M}_{☉}$
$n_s$	0.971	${\Omega }_{\Lambda }$	0.7186	$M_{\star }$	$8.4\times {10}^5{M}_{☉}$
$ϵ$	$1.5h^{-1}\mathrm{kpc}$	$N_{\mathrm{particle}}$	$2\times {264}^3$	$L_{\mathrm{box}}$	$15h^{-1}\mathrm{Mpc}$

Table 2. Table 2 : The sets, names, the BH seed mass M ∙ seed superscript subscript 𝑀 ∙ seed M_{\bullet}^{\rm seed} , and threshold halo mass M fof seed superscript subscript 𝑀 fof seed M_{\rm fof}^{\rm seed} in the simulations.

Set	Name	$M_{\bullet }^{\mathrm{seed}}$ [$h^{-1}{M}_{☉}$]	$M_{\mathrm{fof}}^{\mathrm{seed}}$ [$h^{-1}{M}_{☉}$]
A a,e	B3H8	$5\times {10}^3$	$5\times {10}^8$
A	B4H9	$5\times {10}^4$	$5\times {10}^9$
A, B	B5H10 c,d	$5\times {10}^5$	$5\times {10}^{10}$
B	B4H10	$5\times {10}^4$	$5\times {10}^{10}$
B b	B3H10	$5\times {10}^3$	$5\times {10}^{10}$

Table 3. Table 3 : The numbers of BHs in the simulations at different redshifts.

	B3H8	B4H9	B5H10	B4H10	B3H10
$z=10$	1678	36	1	1	1
$z=8$	4994	135	2	2	2
$z=6$	12210	488	11	11	11

Equations

$$\Gamma =\{C_i\equiv C_i[f;{\mathbf{x}}_i]=c_i;i=1,...,M\}$$

$$\bar{f}(\mathbf{x})=⟨f(\mathbf{x})|\Gamma ⟩={\xi }_i(\mathbf{x}){\xi }_{ij}^{-1}{c}_j$$

$$\mathcal{P}[F|{\Gamma }_p]=\mathcal{P}[F|{\Gamma }_q]$$

$$\begin{array}{rl}f(\mathbf{x})& =F(\mathbf{x})+\bar{f}(\mathbf{x})\\ & =\left(\tilde{f}(\mathbf{x})-{\xi }_i(\mathbf{x}){\xi }_{ij}^{-1}{\tilde{c}}_j\right)+{\xi }_i(\mathbf{x}){\xi }_{ij}^{-1}{c}_j\\ & =\tilde{f}(\mathbf{x})+{\xi }_i(\mathbf{x}){\xi }_{ij}^{-1}\left(c_j-{\tilde{c}}_j\right)\end{array}$$

Full text

††pubyear: 2019††pagerange: The early growth of supermassive black holes in cosmological hydrodynamic simulations with constrained Gaussian realizations–The early growth of supermassive black holes in cosmological hydrodynamic simulations with constrained Gaussian realizations

1 Introduction

The formation of the first supermassive black holes (SMBHs) remains challenging in our standard paradigm of structure formation. SMBHs, as massive as those in galaxies today, are known to exist in the early universe, even up to $z\sim 7.5$. Luminous, extremely rare, quasars at $z\sim 6$ were initially discovered in the Sloan Digital Sky Survey (Fan et al., 2006; Jiang et al., 2009) and, until recently, the highest redshift quasar known (Wu et al., 2015) at $z=7.09$ (Mortlock et al., 2011) has been surpassed by the discovery of a bright quasar at $z=7.54$ (Bañados et al., 2018), which is currently the record holder for known high redshift quasars. A further sample of two $z>7$ has also been recently discovered (Yang et al., 2018). The presence of luminous quasars observed within the first billion years of the Universe highlights that the BH seeds for the SMBH population must have assembled at the cosmic dawn, concurrently with the time of the formation of the first stars or galaxies.

However, the precise SMBH seed formation mechanism remains unknown, nor is it clear if there is only one seed formation channel at play over the entire SMBH seed mass spectrum of models. The current scenarios suggest that seed BHs are (a) remnants of the first generation of stars (PopIII) (e.g. Madau & Rees, 2001; Abel et al., 2002; Johnson & Bromm, 2007) or (b) direct gas collapse within the first massive halos (e.g. Lodato & Natarajan, 2006; Begelman et al., 2006; Regan & Haehnelt, 2009; Ferrara et al., 2014; Latif et al., 2013) or (c) runaway collapse of dense nuclear star clusters (e.g. Begelman & Rees, 1978; Devecchi & Volonteri, 2009; Yajima & Khochfar, 2016; Katz et al., 2015). The seed BHs then range in mass from a few hundred for (a) to $10^{5}M_{\sun}$ for (b) and (c).

In large-volume cosmological simulations, a common and widely used sub-grid model for SMBHs and active galactic nuclei (AGN) feedback has been proposed in Di Matteo et al. (2005). Since the SMBH seed formation process is not resolved by cosmological simulations (see Regan & Haehnelt, 2009, for a review), it is assumed that every halo above a certain threshold mass hosts a central BH seed. Halos are selected for seeding by regularly running the ’Friends-of-Friends’ (FoF) halo finder on the dark matter distribution. The BH seed mass ($M_{\bullet}^{\rm seed}$) and the threshold halo mass ($M_{\rm fof}^{\rm seed}$) are the parameters in simulations. Although this is an ad-hoc seeding procedure, the initial seed BH mass subsequently grows in these simulations via mergers and accretion. Many simulations have adopted this or a similar scenario and gotten good agreements with observations such as the MassiveBlack simulation (Di Matteo et al., 2012), the Illustris simulation (Vogelsberger et al., 2014), the Evolution and Assembly of GaLaxies and their Environment (EAGLE) suite of SPH simulation (Schaye et al., 2015), the MassiveBlack II simulation (Khandai et al., 2015), and the BlueTides simulation (Feng et al., 2016a). Some recent studies have implemented different, physically motivated approaches where the BH seeding is based on gas properties such as Bellovary et al. (2011); Habouzit et al. (2017); Tremmel et al. (2017). However, it is worth noting that these models were adopted in much smaller volume simulations than, for example, MassiveBlack II or BlueTides. From those, it is not possible to validate the basic statistical properties of the BH population (e.g. luminosity functions or mass functions) against currently observed samples.

Taking BlueTides, a large-volume and high-resolution cosmological hydrodynamic simulation with $2\times 7040^{3}$ particles in a box of $400\,h^{-1}{\rm Mpc}$ on a side, as an example, SMBHs are modeled as follows. For each FoF halo with a mass above $M_{\rm fof}^{\rm seed}=5\times 10^{10}\,h^{-1}M_{\sun}$, a SMBH is seeded with an initial seed mass $M_{\bullet}^{\rm seed}=5\times 10^{5}\,h^{-1}M_{\sun}$ at the position of the local minimum potential if there is no SMBH in that halo. After being seeded, gas accretion proceeds according to Bondi (1952) while the BH accretion rate is limited to two times the Eddington rate. When SMBHs are accreting, we assume that some fraction of the radiated luminosity can couple thermally and isotropically to surrounding gas in the form of feedback energy (Springel, 2005; Di Matteo et al., 2005).

Adopting this SMBH model and appropriate sub-grid physics for the galaxy formation modeling, the BlueTides simulation has predicted various quantities in good agreements with current observational constraints in the high-$z$ universe such as UV luminosity functions (Feng et al., 2016a; Waters et al., 2016a, b; Wilkins et al., 2017), the first galaxies and the most massive quasars (Feng et al., 2015; Di Matteo et al., 2017; Tenneti et al., 2018), the Lyman continuum photon production efficiency (Wilkins et al., 2016, 2017), galaxy stellar mass functions (Wilkins et al., 2018), angular clustering amplitude (Bhowmick et al., 2018), BH-galaxy scaling relations (Huang et al., 2018), and gas outflows from the $z=7.54$ quasar (Ni et al., 2018). Important for our work here, BlueTides, with its large volume and appropriate resolution, is currently the only cosmological hydrodynamic simulation that makes direct contact with the rare, first quasar population at $z>7$.

However, an essential question for the SMBH sub-grid model is how different parameters (e.g. $M_{\bullet}^{\rm seed}$ or $M_{\rm fof}^{\rm seed}$) may affect the growth of SMBHs and the local environment in cosmological simulations. Changing the BH seed mass and re-running such a large-volume simulation multiple times is completely prohibitive because it is computationally expensive even on the largest current national facilities. To reduce the demand on computational resources, a common method is to run a "zoom-in re-simulation" with a higher resolution or different physical parameters from a certain region selected from a large-volume lower-resolution simulation. This allows people to focus on a specific environment numerically and has been applied to study SMBHs in simulations for various purposes (e.g. Li et al., 2007; Sijacki et al., 2009; Hopkins & Quataert, 2010; Bournaud et al., 2011; Romano-Diaz et al., 2011; Bellovary et al., 2013; Dubois et al., 2013; Anglés-Alcázar et al., 2014; Costa et al., 2014; Feng et al., 2014).

In this paper, we combine the technique of constrained Gaussian realization and cosmological hydrodynamic simulations to reduce the demand on computational resources. Hoffman & Ribak (1991) first introduced an optimal solution to the problem of the construction of constrained realizations of Gaussian fields by demonstrating how the algorithm generates constrained fields with a simple single-density peak. Later on, van de Weygaert & Bertschinger (1996) in addition developed an algorithm to set up initial Gaussian random density and velocity fields containing multiple constraints of arbitrary amplitudes and positions. Integrating the algorithm to cosmological hydrodynamic simulations has arisen in the past few years to explore dark matter halos and galaxy formation (Roth et al., 2016; Porciani, 2016; Pontzen et al., 2017).

With the constrained Gaussian realization, we can constrain the initial density field by adding a desirable height of a density peak when generating the initial condition such that a more massive halo can still form in a relatively small box compared to those large-volume ($\sim{\rm Gpc}$ per side) cosmological simulations (with uniform/unconstrained initial condition) . For instance, we can grow a halo with a mass $\sim 10^{12}\,M_{\sun}$ in a box of $15\,h^{-1}{\rm Mpc}$ on a side at $z=8$, whose mass is similar to the one hosting the most massive BH in BlueTides simulation ($400\,h^{-1}{\rm Mpc}$ on a side) under the same resolution. This reduces the computational demand by a factor of $(400/15)^{3}\sim 20000$. This approach is a more general way to study the growth of SMBH compared to the zoom-in re-simulation method because the goal of the latter is to exactly study a particular object/region (for example, a particular halo). However, our approach is aiming to study characteristic environments by creating one or more different realizations but with similar properties such as halo mass or tidal field to the object/region of our interest (which we extract, for comparison, in the uniform large volume simulations with the exact physics).

As we shall further demonstrate, besides the density constraint, we need another condition related to the ICs that induces the fastest growth for SMBHs in cosmological simulations. This is expected, as the observed population of quasar-like SMBHs at high redshifts is even rarer than the massive halos. For example, there is only one SMBH with a mass above $10^{8}\,M_{\sun}$ in a halo with a mass $\sim 10^{12}\,M_{\sun}$ while there are $>50$ halos more massive than that halo at $z=8$. An environmental property, also related to the ICs that has been found to be relevant to induce fast BH growth is the local tidal field strength (Di Matteo et al., 2017). In particular, using the large volume BlueTides simulations, Di Matteo et al. (2017) has shown that isolated regions of low tidal fields are key to the fast growth of the firs SMBHs. As a consequence, we also choose the realization with a lower tidal field around the local environment where the halo forms, which indeed helps more massive SMBHs grow in simulations (see Section 2.4).

After significantly decreasing the demand on computational resources with the constrained Gaussian realizations and a lower tidal field realization, we are finally able to examine how sensitive the SMBH growth is to the BH seed mass in the sub-grid model by running multiple cosmological simulations. According to the different hypotheses of the BH formation scenario, the BH seed mass has been suggested to range from $10^{2}$ to $10^{6}\,h^{-1}M_{\sun}$ (Haehnelt & Rees, 1993; Loeb & Rasio, 1994; Eisenstein & Loeb, 1995; Bromm & Loeb, 2003; Koushiappas et al., 2004; Begelman et al., 2006; Lodato & Natarajan, 2006; Zhang et al., 2008; Volonteri, 2010; Latif et al., 2013; Schleicher et al., 2013; Ferrara et al., 2014). Therefore we focus on the three different seed masses: $5\times 10^{3}$, $5\times 10^{4}$, and $5\times 10^{5}\,h^{-1}M_{\sun}$, with the same ratio of halo to BH seed mass and with the same halo threshold mass in this paper.

We organize the paper as the following. Section 2 describes the constrained Gaussian realizations, compares the constrained and unconstrained initial conditions and simulations, and discusses the effect of different tidal fields on the local environment of SMBHs. Section 3 demonstrates the results of the early growth of SMBHs and their hosts in the simulations with different BH seeding scenarios. Section 4 concludes the paper.

2 methodology

2.1 Constrained initial conditions

The first quasars are extremely rare and hence the massive halos hosting these first SMBHs also have to be commensurately rare. Traditionally one needs extremely large-volume simulations to simulate such objects. For example, there is only one SMBH with a mass above $10^{8}\,M_{\sun}$ in a halo with a mass of $\sim 10^{12}\,M_{\sun}$ at $z=8$ in the BlueTides simulation with $400\,h^{-1}{\rm Mpc}$ on a side of the cube. Using simulations of such scale to study the effect on BH growth due to subgrid prescriptions is prohibitively expensive. In this work, we work with much smaller simulation boxes ($15\,h^{-1}{\rm Mpc}$), but we add constraints to the initial condition (linear field) to ensure the existence of extreme density peaks. With high-density peaks in the linear field, we can guarantee the formation of massive halos at later times even in a smaller simulation box. These halos are then used to study the early growth of SMBHs in the massive halos in cosmological hydrodynamic simulations with different BH seeding parameters.

According to the linear growth theory, to generate a $10^{12}\,M_{\sun}$ halo at z=8, we need a $5\,\sigma$ peak in the underlying over-density field. To achieve that, we use FastPM, a particle mesh based quasi N-body solver (Feng et al., 2016b), to generate initial conditions with constrained Gaussian realizations for the simulations in the paper for the first time. Thanks to the contribution of Aslanyan et al. (2016), FastPM is capable of producing constrained Gaussian density fields to the initial condition, which we will describe the basic idea and the implementation in the following paragraphs.

First introduced by Hoffman & Ribak (1991), the constrained realization technique was then explained in a lot more detail in van de Weygaert & Bertschinger (1996) which we refer readers to. The goal is to construct a field $f(\mathbf{x})$ subject to a set of $M$ constraints:

[TABLE]

where the constraint $C_{i}$ can be viewed as a functional of $f(\mathbf{x})$ field (here in our specific case, the overdensity field) to have the specific value $c_{i}$ at the position $\mathbf{x}_{i}$.

To obtain a field $f(\mathbf{x})$ satisfying the constraint $\Gamma$, one can start with a random, unconstrained Gaussian realization $\tilde{f}(\mathbf{x})$ and impose on that an "ensemble mean field" $\bar{f}(\mathbf{x})$ corresponding to the desired constraint $\Gamma$. More specifically, the ensemble mean field $\bar{f}(\mathbf{x})$ can be written in the form of:

[TABLE]

where $\xi_{i}(\mathbf{x})=\langle f(\mathbf{x})\,C_{i}\rangle$ is the cross-correlation between the $f(\mathbf{x})$ field and the $i^{\rm th}$ constraint $C_{i}$, and $\xi^{-1}_{ij}$ is the $\left(i,j\right)$ element of the inverse of the constraint covariance matrix $\langle C_{i}\,C_{j}\rangle$. Note that the summation over repeated indices is used. The ensemble mean field $\bar{f}(\mathbf{x})$ can be interpreted as the "most likely" field subject to the set of constraints $\Gamma$.

We further introduce the "residual field" $F(\mathbf{x})\equiv f(\mathbf{x})-\bar{f}(\mathbf{x})$ as the difference between an arbitrary Gaussian realization $f(\mathbf{x})$ satisfying the constraint set $\Gamma$ and the ensemble mean field $\bar{f}(\mathbf{x})$. The crucial idea of the constrained realization construction method is based on the fact that, the complete probability distribution $\mathcal{P}[F\,|\,\Gamma]$ of the residual field $F(\mathbf{x})$ is independent of numerical values $c_{i}$ of the constraints in $\Gamma$ (c.f. Hoffman & Ribak, 1991; van de Weygaert & Bertschinger, 1996, for the detailed derivations). That is, for any $\Gamma_{p}$ and $\Gamma_{q}$ where $p\neq q$,

[TABLE]

Therefore, we can construct the desired realization under a constraint set $\Gamma$ by properly sampling the residual field $F(\mathbf{x})$ from a random, unconstrained realization $\tilde{f}(\mathbf{x})$ and then add that $F(\mathbf{x})$ to the ensemble field $\bar{f}(\mathbf{x})$ corresponding to $\Gamma$. The formalism can be written as

[TABLE]

In other words, we treat the original $\tilde{f}(\mathbf{x})$ as a field subject to a constraint set $\tilde{\Gamma}$ with $\tilde{c}_{j}=C_{j}[\tilde{f};\mathbf{x}_{j}]$ where $\tilde{c}_{j}$ is the original value of the unconstrained field. We then have the ensemble mean field corresponding to $\tilde{\Gamma}$ as $\tilde{\bar{f}}(\mathbf{x})=\xi_{i}(\mathbf{x})\,\xi^{-1}_{ij}\,\tilde{c}_{j}$. Getting the residual field $F(\mathbf{x})$ from a random unconstrained realization by $\tilde{f}(\mathbf{x})-\tilde{\bar{f}}(\mathbf{x})$, we then add $F(\mathbf{x})$ to $\bar{f}(\mathbf{x})$ to obtain the field $f(\mathbf{x})$ satisfying the constraint $\Gamma$. It is well established in van de Weygaert & Bertschinger (1996) that the $f(\mathbf{x})$ field constructed in this way is a properly sampled realization subject to the desired constraint $\Gamma$. This is what we implemented in our code.

We note that, however, one limitation in our implementation is that the super-sampling variance (DC mode) is missing. Super-sampling variance is the effect of the coupling to modes at scales larger than the box size (Li et al., 2014). In our simulation, we assume that the overdensity of the whole simulation box to be zero, i.e., the DC mode is zero. The DC mode can be incorporated by the so-called separate universe technique that absorbs the overdensity of the simulation volume into a modified cosmology. (see, e.g., Sirko, 2005; Gnedin et al., 2011; Wagner et al., 2015; Li et al., 2014, 2018, for more details). However, we estimate that with our simulation box size ($15\,h^{-1}{\rm Mpc}$ per side) at these high redshifts ($z>6$), this effect accounts only about 10 percent.

Here we demonstrate the constrained realization generated via FastPM with a single $5\,\sigma$ density peak. Figure 1 shows examples of the density field with and without a constrained density peak and its associated residual map in a dense region at the center in the domain with a box size of $15\,h^{-1}{\rm Mpc}$. As expected, we find that the density increases in the region where we put the constraint without changing the overall pattern of the density field. According to the residual map, 0.04 percent of the pixels have greater than 10 percent residuals; 2 percent of the pixels exceed 5 percent residuals; none of the pixels have residuals less than $-5$ percent.

2.2 Simulation setup

We use the massively parallel cosmological smoothed particle hydrodynamic (SPH) simulation software, MP-Gadget (Feng et al., 2016a), to run all the simulations in this paper. Its hydrodynamics solver adopts the new pressure-entropy formulation of SPH (Hopkins, 2013). The main sub-grid models in MP-Gadget are

• •

star formation based on a multiphase star formation model (Springel & Hernquist, 2003) with modifications following Vogelsberger et al. (2013),

• •

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 (Nelson et al., 2015),

• •

SMBH growth and AGN feedback (Di Matteo et al., 2005).

All the new constrained simulations in the paper are run with periodic boundary conditions from $z=99$ to $z=6$ (as we are mostly interested in the seed mass and the early growth of SMBHs). Each simulation contains $2\times 264^{3}$ particles in a cube with the box size $L_{\rm box}=15\,h^{-1}{\rm Mpc}$ (the choice of $L_{\rm box}$ size will be further discussed later in detail). We adopt the cosmological parameters based on the Nine-Year Wilkinson Microwave Anisotropy Probe Observations (Hinshaw et al., 2013). All the simulations have the same cosmology and resolution as in BlueTides, so that we can always use the direct large volume simulations to assess the validity of the new simulations. Table 1 summarizes all the basic parameters of our new runs. Note that a star particle has a mass of $M_{\star}=\frac{1}{4}\,M_{\rm gas}$ and that the gravitational smoothing length $\epsilon$ is the same for all kinds of particles.

Constraining a high-density peak in the initial density field to get a massive halo allows us to study rare objects in a small simulation box rather than in large-volume cosmological hydrodynamic simulations with the box size of a few hundred Mpc on the side. In order to set a box size for our constrained runs, we need to make sure that the growth history of the halo needs to be well converged. In particular, we would like to make sure that all the particles that make into the halo and all the way into the central BH are captured. In particular, here, we want to track the growth of SMBHs at the center and their host galaxies, so we look into how far the particles in the halo hosting the most massive BH in BlueTides have traveled from $z=99$ to $z=8$ in Figure 2. The mean of the travel distance is $2\,h^{-1}{\rm Mpc}$ with a standard deviation of $0.5\,h^{-1}{\rm Mpc}$, indicating that at least a box size of $L_{\rm box}\geq 2\,h^{-1}{\rm Mpc}$ is necessary to contain the halo up to $z=8$. As we will run our new simulations beyond $z=8$ here we also make use of a dark matter only realization of BlueTides, the BTMassTracer (Tenneti et al., 2018). We track all the dark matter particles in the halo down to $z=5$. We show that the particles typically travel a mean of $4.6\,h^{-1}{\rm Mpc}$ and a standard deviation of $0.2\,h^{-1}{\rm Mpc}$. This implies an absolute minimum box size of $L_{\rm box}\sim 5\,h^{-1}{\rm Mpc}$. To be rather conservative and make sure we have the appropriate growth history of the halo and its black holes, we choose a size of $15\,h^{-1}{\rm Mpc}$ and stick to this for all the simulations in the paper.

2.3 Constrained versus unconstrained simulations

To illustrate the basic features of the constrained simulations, we first run the constrained and unconstrained initial conditions in Figure 1 down to $z=6$ while keeping all the other simulation parameters. Figure 3 shows the density fields of the constrained and unconstrained simulations at $z=6$ color-coded by temperature as well. As expected, the density around where we put the constrained peak in the constrained simulation is higher than the unconstrained one while the overall structure maintains. So is the temperature. Figure 4 shows the halo and stellar mass functions ($\Phi_{\rm halo}$ and $\Phi_{\rm\star}$) at $z=6$, 8, and 10, compared with BlueTides. In particular, there is one halo in the very massive end of these functions in the constrained simulation due to the constrained high-density peak. Aside from the massive objects, the consistency of both mass functions with each other and with the ones in BlueTides indicates that the constrained simulation appropriately captures the growth of halo and stellar mass function statistically.

We then investigate the growth history of the most massive BHs ($M_{\bullet}$) and their hosts (halo mass $M_{\rm halo}$ and stellar mass $M_{\star}$) in the two simulations compared with that of the BlueTides simulation in Figure 5. With a proper density peak, the growth history of the three masses in the constrained simulation converges to the ones in BlueTides (note that a total convergence is not expected as this a new constrained simulation but not a zoom-in simulation). On the other hand, the halo mass of the unconstrained simulation at $z=6$ is an order of magnitude less massive than the one in the constrained simulation; the stellar mass is around two orders of magnitude less massive; the BH mass is three orders of magnitude less massive.

2.4 Tidal fields of the SMBHs

While a highly biased region (as in our constrained simulations) is a necessary condition for growing a massive BH, it is not sufficient. For example in the BlueTides simulation, only one of the 50 most massive halos of mass similar or greater than the one hosting the most massive BH has a BH more massive than $10^{8}\,M_{\sun}$. As we shall further show, not all of the constrained Gaussian realizations that can grow massive halos guarantee to grow SMBHs in them as well. Directly related to the density field in the initial conditions, the local tidal field has been identified as the environmental property that is the most strongly correlated to the growth of the first quasars in Di Matteo et al. (2017). In these findings, the extreme early growth depends on the early interplay of high gas densities and the tidal field that shapes the mode of accretion in those halos.

The tidal field is characterized by the three eigenvalues ($t_{1}$, $t_{2}$, $t_{3}$) of the local tidal tensor $T_{ij}\equiv S_{ij}-\frac{1}{3}\sum_{i}S_{ii}$, where the strain tensor is the second derivative of the potential, $S_{ij}\equiv\nabla_{i}\nabla_{j}\phi$. According to Dalal et al. (2008), $S_{ij}$ is calculated in Fourier space as $\hat{S}_{ij}=\frac{k^{2}}{k_{i}k_{j}}\hat{\delta}$. The three eigenvalues are by definition $t_{1}>t_{2}>t_{3}$ and satisfy $t_{1}+t_{2}+t_{3}=0$ so that $t_{1}$ is always positive and $t_{3}$ is negative. Thus, the tidal field stretches material along with ${\bf t_{1}}$ and compresses material along with ${\bf t_{3}}$, where (${\bf t_{1}}$, ${\bf t_{2}}$, ${\bf t_{3}}$) are the corresponding eigenvectors. To use $t_{1}$ as the indicator of the local tidal field strength following the standard usage, we calculate $t_{1}$ numerically using nbodykit (Hand & Feng, 2015). We read all the particle data from a snapshot into a mesh object weighted by the particle mass to get $\hat{\delta}$; transform them to Fourier space; apply a kernel of $\frac{k^{2}}{k_{i}k_{j}}$ to get $\hat{S}_{ij}$; transform them back to the real space and evaluate $T_{ij}$ at the position of the SMBH.

To further evaluate the role of the tidal field in the growth of the first massive SMBHs, we generate a number of constrained realizations; select the ones with the minimum, intermediate, and maximum $t_{1}$ around the density peak as the initial conditions; run them from $z=99$ to $z=6$. The left panel of Figure 6 shows $t_{1}$ at the position of the most massive BHs in the three simulations. $t_{1}$ measured at $z=6$ on different scales of 1 – $5\,h^{-1}{\rm Mpc}$ in the inset suggests that the simulation with a lower or higher $t_{1}$ is always lower or higher across the scales. The evolution of $t_{1}$ measured at the scale of $1\,h^{-1}{\rm Mpc}$ in the main panel shows that a lower or higher $t_{1}$ environment tends to maintain a lower or higher $t_{1}$ as time goes.

The right panel of Figure 6 shows the growth history of the most massive BHs and their hosts in the three simulations. Several interesting results we find include that, the masses of the halos $M_{\rm halo}$ always differ by a factor less than 10 among the three simulations; the stellar mass $M_{\star}$ of the low $t_{1}$ simulation is around an order of magnitude higher than the others at an earlier stage; the BH masses $M_{\bullet}$ differ by a factor of 10 – 100. The three simulations suggest that the tidal field has a larger impact on the growth of SMBHs: a SMBH can grow more or less massive when it is in a lower or higher $t_{1}$ surrounding environment. Besides, the growth history of the most massive BH and its host galaxy and halo in the low $t_{1}$ simulation also converges better to the BlueTides simulation.

The fact that a lower tidal field environment helps a more massive growth of the SMBHs in our simulations strengthens the findings in Di Matteo et al. (2017) that the local tidal field is strongly correlated to the growth of the first quasars. Moreover, we utilize the constrained realization that provides the lowest $t_{1}$ environment as the initial condition for the study of the early growth of SMBHs with different BH seeding parameters in the following sections.

3 Results: different SMBH seeding scenarios

The main objective of this work is to study the effect of different BH seeding parameters on the early growth of SMBHs using a set of constrained cosmological simulations that statistically reproduce the environments for early growth in the large-volume BlueTides simulation. Here, we conduct two sets of new constrained simulations and investigate the growth history of SMBHs starting from different BH seed masses. In particular, we use three different BH seed masses $M_{\bullet}^{\rm seed}=5\times 10^{3}$, $5\times 10^{4}$, and $5\times 10^{5}\,h^{-1}M_{\sun}$ and perform two sets of simulations.

• •

In Set A: we lower the halo mass threshold to $M_{\rm fof}^{\rm seed}=5\times 10^{8}$, $5\times 10^{9}$, and $5\times 10^{10}\,h^{-1}M_{\sun}$ commensurate with keeping the ratio of $M_{\bullet}^{\rm seed}/M_{\rm fof}^{\rm seed}$ constant. This is motivated by the physical models for BH seed formation implying that smaller BHs may form earlier in the first, molecular cooling halos which have smaller mass (see, e.g., Johnson & Bromm, 2007). Hence in this set, BH seeds with smaller masses (than the canonical $5\times 10^{5}\,h^{-1}M_{\sun}$) are seeded at earlier times.

• •

In Set B: we fix the threshold halo mass at $M_{\rm fof}^{\rm seed}=5\times 10^{10}\,h^{-1}M_{\sun}$ and study the effect of changing the BH seed mass in a given fixed halo mass. All the BH seeds in this set are seeded at the same time and in the same halos but simply with different BH seed masses.

Table 2 summarizes the sets of simulations, the adopted naming and their respective BH and halo seeding parameters. We emphasize that all of the simulations have the same constrained initial condition. Also, in particular, B5H10 has the same BH seeding parameters as that of BlueTides with the canonical/reference choice of BH seed mass and halo mass.

To illustrate the results of our simulations, we start by showing the environments of the BHs at $z=6$, 8, and 10 in Figure 7. In particular, we show the projected gas density field color-coded by the gas temperature in each of the simulations. Each panel is $6\,h^{-1}{\rm Mpc}$ per side with the most massive BH residing at the center. The green circles mark out all the BHs, where the size of the circles scales with the BH masses. The relatively hot region of gas around the BH results mostly from the effects of AGN feedback. Here we specifically show the results of the simulations in Set A (B3H8, B4H9, and B5H10). We note that the density field/environments of B3H10 and B4H10 in Set B will be similar to those shown for B5H10 except for the gas temperatures in the region around the central BH (which would typically be less affected by AGN feedback and associated heating, as we will discuss later).

Figure 7 highlights two major points. First, the effects of BH feedback, represented by the heated gas phase (reddish colors) are more prominent as the seed BH grows and as in the case for a single, larger BH seed. Second, the BH populations in the simulations B3H8, B4H9, and B5H10 are different. This is due to the adopted values for the threshold halo mass in each of the respective simulations. In particular, the lower the threshold halo mass is, the more BH seeds are placed in a simulation, resulting in more BHs in B3H8 than in B5H10. Table 3 summarizes the numbers of BHs seeded and growing in each of the simulations.

Figure 8 shows the halo mass functions $\Phi_{\rm halo}$ and the stellar mass functions $\Phi_{\star}$ in all our simulations (both Set A and Set B). The consistency of the mass functions among the simulations suggests that the choice of BH seeding parameters does not affect the global halo and galaxy population in significant ways. In particular, the lower-mass end is virtually unaffected in both halo and stellar mass functions. The high-mass end of the stellar mass function shows some differences. As we shall show later, this is a reflection of the different SFR histories (for $z<10$ in the different seed models which are modulated by different amounts of AGN feedback in different BH seed models.

In the following sections, we aim to explore in more detail the growth histories of most massive BHs and their hosts in the simulations including their masses and mass assembly rates. In particular, in Section 3.1 we will show results for the simulations with different BH seed masses and different halo thresholds, B3H8, B4H9, and B5H10 while in Section 3.2 the simulations with different BH seed masses at fixed halo mass, B3H10, B4H10, and B5H10.

3.1 Set A: different BH seed masses and halo mass thresholds

Here we describe the results from the simulations B3H8, B4H9, and B5H10. Those are the ones in which halo thresholds for BH seeding are adjusted such that the ratio $M_{\bullet}^{\rm seed}/M_{\rm fof}^{\rm seed}$ is fixed. The left-top panel of Figure 9 shows the growth history of host halo and galaxy ($M_{\rm halo}$ and $M_{\star}$) for the most massive BHs in the simulations. The halos show the same growth history among the simulations and their masses at $z=6$ are $3\times 10^{12}\,h^{-1}M_{\sun}$. The galaxy stellar mass also have a similar growth history except for B4H9 in which the host stellar mass ends up being larger than the others after $z=8$ but less than a factor of two different at $z=6$ according to the stellar mass ratio in the left-bottom panel. The stellar mass ratio $f_{M_{\star}}=\frac{M_{\star}}{M_{\star}^{\rm B5H10}}$ compares the galaxy mass in each of the simulation to $M_{\star}$ in B5H10. At $z=6$, $M_{\star}=10^{11}\,h^{-1}M_{\sun}$ for B3H8 and B5H10 and $M_{\star}=6\times 10^{10}\,h^{-1}M_{\sun}$ for B4H9. Our simulations suggest that the choice of BH seeding parameters does not affect the growth of the hosts for more than a factor of two.

Of our interest, the solid curves in the right panel of Figure 9 show the growth history of the most massive BHs in the simulations B3H8, B4H9, and B5H10. As apparent in the figure, the BHs are seeded with different masses according to $M_{\bullet}^{\rm seed}$ and at different times according to $M_{\rm fof}^{\rm seed}$; that is, a smaller BH seed emerges in a lower-mass halo at an earlier time than a corresponding higher-mass one. For example, when a BH seed of $5\times 10^{3}\,h^{-1}M_{\sun}$ is placed in halos of mass $5\times 10^{8}\,h^{-1}M_{\sun}$ it can be seeded at $z>20$. As halos of this mass are not so rare at high redshifts, a lot of BH seeds emerge rather than just one seed does as in the case of $M_{\bullet}^{\rm seed}=5\times 10^{5}\,h^{-1}M_{\sun}$. The most massive BHs in the simulations start to converge in mass at $z\sim 8$ and reach a mass of $2\times 10^{9}\,h^{-1}M_{\sun}$ by $z=6$ even though the seed population and the total number of BHs are different. This suggests that the choice of different pairs of BH seed mass and threshold halo mass does not significantly affect the growth of the most massive BHs, at least in the expected range of $M_{\bullet}^{\rm seed}=10^{3}$ – $10^{6}\,M_{\sun}$. However, the early growth of SMBHs in the three simulations can be faster or slower at $z>10$; that is, the small and large BH seeds start more massively while the intermediate one remains at its seed mass the longest but catch up drastically once it starts growing. Moreover, the discrete jumps in the growth history of B3H8 at early times indicate a lot of mergers occur even at such high redshifts compared to the others.

Figure 10 shows the evolution of BH accretion rates $\dot{M}_{\bullet}$ of the most massive BHs in the simulations. The solid curves show the $\dot{M}_{\bullet}$ and the shaded bands indicate the regime between the Eddington rate and two times the Eddington rate as the upper limit for $\dot{M}_{\bullet}$ in the simulations. We find that $\dot{M}_{\bullet}$ has the same overall evolution: starting with an initial low accretion phase, followed by a close to exponential Eddington growth, and ending with a final quenched feedback-regulated phase. It is noticeable that $\dot{M}_{\bullet}$ gradually converges at $z>10$, the time when the three SMBHs enter the feedback-regulated phase where their growth saturates and $M_{\bullet}$ starts converging. However, in the early phases during $10<z<13$, the BH accretion rates have different trajectories as the BHs experiencing exponential growth but constrained by the upper limit in the simulations.

The right panel of Figure 10 shows the star formation history of the host galaxies for the simulations. The evolution of the star formation rates (SFRs) appears similar at $z>10$ but diverges in the later phase because the AGN feedback starts to regulate the star formation rate by coupling significant energy to the star forming gas. Therefore, $\dot{M}_{\bullet}$ and SFR start to couple with each other after a significant BH growth phase ($z<10$) though the physical scales of the two quantities are quite different (SFR is in the galactic scale of tens of $h^{-1}{\rm kpc}$, whereas $\dot{M}_{\bullet}$ is determined by local gas density at the scale of $h^{-1}{\rm kpc}$).

As a further comparison, Feng et al. (2014) has examined the exact pairs of $M_{\bullet}^{\rm seed}$ and $M_{\rm fof}^{\rm seed}$ in Set A (see Table 2) using zoom-in simulations. With the zoom-in technique, they re-simulated a high-redshift ($z>5.5$) halo hosting a $10^{9}\,M_{\sun}$ BH from the $\rm\sim Gpc$ volume, MassiveBlack cosmological hydrodynamic simulation. They reported that regardless of the BH seed mass, the BH masses converged to $M_{\bullet}=10^{9}\,M_{\sun}$ at $z=6$, the BHs underwent a similar history of BH accretion, and the evolution of SFRs was the same at the earlier times before AGN feedback starts to regulate the SFR. It is interesting that both their findings and ours are fully consistent with completely different methods, the zoomed-in simulation and the constrained simulation. This strengthens the conclusion that the choice of BH seed mass - when keeping the ratio $M_{\bullet}^{\rm seed}/M_{\rm fof}^{\rm seed}$ fixed - does not affect the growth of SMBHs significantly.

3.2 Set B: different BH seed masses at fixed host halo mass

We then move on to investigate the growth history of host halo and galaxy ($M_{\rm halo}$ and $M_{\star}$) of the most massive BHs in the simulations B3H10, B4H10, and B5H10 in Figure 9. The halos show the same growth history among the simulations and their masses at $z=6$ are $3\times 10^{12}\,h^{-1}M_{\sun}$. On the other hand, the galaxies seem to have a very similar growth history at $z>10$ but then their masses start to differ. According to the stellar mass ratio $f_{M_{\star}}$, the galaxy in B3H10 is three times more massive than the one in B5H10 at $z=6$. This can be inferred through the evolution of SFR in Figure 10; that is, the galaxy in B3H10 undergoes a rather bursty star formation history at $z<10$ with variations in SFR up to an order of magnitude compared to the others. The reason is that the BH mass in B3H10 is smaller than the others and therefore the corresponding AGN feedback is not strong enough to bring sufficient suppression on the local star formation at $6<z<10$. Our simulations suggest that the choice of BH seed mass does not affect the growth of halo by $z=6$; does not affect the growth of galaxy by $z=10$; do affect the growth of galaxy after $z=10$ but less than a factor of three by $z=6$. There is another noticeable trend that a galaxy grows more when the BH seed mass is smaller.

The growth of the most massive BHs in Figure 9 starts with different masses and then converges to $\sim 10^{9}\,M_{\sun}$ by $z=6$ except for the one in B3H10, which is less than an order of magnitude difference. The evolution of BH accretion rates of the BHs in Figure 10 shows that the BH in B3H10 is still experiencing Eddington exponential growth at $z<8$. This indicating that the BH in B3H10 is still catching up in mass and will probably converge to the others at later times. Our simulation results hint that the growth of the most massive BHs will converge at the later times regardless of the choice of the SMBH seeding parameters in cosmological simulations.

3.3 BH-BH mergers

Figure 7 and Table 3 indicate that there is a major difference in the BH populations of the simulations B3H8, B4H9 and B5H10. Particularly in B3H8, there is a vast BH population in the environment of the most massive BH since we also adjust the minimum halo, implying that BH mergers are more likely to happen in the early times. The step-like feature in the BH mass assembly history in Figure 9 then infers that mergers occur. These pieces of evidence motivate the following investigation of BH mergers in the simulations. In MP-Gadget, a BH-BH merger occurs when the distance between two BHs is smaller than their SPH smoothing kernel and the relative velocity of the two BHs is smaller than $\frac{1}{2}\,c_{\rm s}$, where $c_{\rm s}$ is the local sound speed of the gas. Convection is also applied to the BH dynamics by repositioning the BH particle to the local minimum potential at every time step.

To investigate the BH merger history of the simulations B3H8, B4H9, and B5H10 (as well as B3H10 and B4H10), we plot the merger tree of BH projected positions in Figure 11. The left and right panels show the position of the most massive BHs in B3H8 and B4H9 respectively compared with B5H10. The size of the data points scales with the BH masses; the data points are color-coded by the ID of BH particles; the black diamonds mark where and when the BH mergers occur. By $z=6$, there are eight and six BH mergers that happen in B3H8 and B4H9 respectively whereas there is no merger in B5H10 despite mergers likely to occur below this redshift as the closest BH below 100 kpc distance from the most massive one. Besides, B3H10 and B4H10 have no merger as well as B5H10 since they contain the same number and position of BHs but with smaller BH masses.

At $z>12$, interestingly, four mergers happen in B3H8 whereas none in the others. This explains why the SMBH in B3H8 grows faster than the others during the earlier phase and further implies that mergers dominate the early growth of SMBHs in small BH seeding scenarios. Our simulations suggest that if less-massive BH seeds are more common, SMBHs can still grow via mergers at the early times even though they may be expected to grow slower. In other words, a different BH merger history results in a different growth of the SMBHs particularly at the early times. Despite the fact that the high halo occupation fraction of SMBHs in cosmological simulations will increase the number of BH-BH mergers since these simulations are implemented with simple merger models without considering the BH dynamics that could make BH-BH mergers more difficult between low-mass BHs, it is still interesting that different seeding scenarios are expected to produce different BH populations and associated merger rates that can discriminate the different scenarios at early times while the final BH mass converges to a similar value.

In contrast, at $6<z<10$, SMBHs seeded with a mass of $10^{3}$ – $10^{4}\,M_{\sun}$ undergo a few BH mergers whereas the one in the largest seed models ($\sim 10^{5}\,M_{\sun}$) does not experience any merger until $z<6$. These different predictions of the merger history for the first massive BHs constitute an interesting prospect for constraining BH seed masses or models for the first quasars that will become within reach with the planned LISA mission (Amaro-Seoane et al., 2017).

4 Conclusion

In this paper, we have investigated new constrained cosmological simulations designed to reproduce the environments and large-scale structures relevant for the growth of the first quasars at $z\geq 6$. In particular, we have focused on the effects of different choices of BH seeding scenarios (different parameters in the SMBH sub-grid model) on the growth of SMBHs at the early times. Employing the technique of constrained Gaussian realizations (Hoffman & Ribak, 1991; van de Weygaert & Bertschinger, 1996), we have reconstructed the initial conditions to reproduce the large-scale structure and the local environment of the most massive BH in the BlueTides simulation. BlueTides has been the only cosmological hydrodynamic simulation that directly predicted the rare-observed first quasars (Di Matteo et al., 2017; Ni et al., 2018; Tenneti et al., 2018) thanks to its sufficiently high resolution and large volume. The first quasars are extremely rare such that there were only four SMBHs with mass $\sim 10^{9}\,M_{\sun}$ by $z=7$ in BlueTides with $L_{\rm box}=400\,h^{-1}{\rm Mpc}$.

We have compared the new constrained simulations with the BlueTides simulation to validate this method by running the constrained initial conditions forward in time until $z=6$. Our new simulations in boxes of $15\,h^{-1}{\rm Mpc}$ on a side have successfully recovered the evolution of the large-scale structure, the mass functions, as well as the growth history of the most massive BHs and their hosts at the high redshifts of interests. At $z=8$, the most massive BH and its hosts had a halo mass of $\sim 10^{12}\,M_{\sun}$; a stellar mass of $\sim 4\times 10^{10}\,M_{\sun}$; a BH mass of $\sim 4\times 10^{8}\,M_{\sun}$. This is consistent with BlueTides within a factor of 1.5 in mass while keeping the resolution. More importantly, the demand on computational resources has decreased significantly by a factor of $(400/15)^{3}\sim 20000$.

By running a set of different realizations such that each of them has a different local tidal field, we have further shown that a low-tidal field environment is crucial for the growth of the earliest and most massive SMBHs. This is consistent with the finding from BlueTides in Di Matteo et al. (2017). For our highest tidal field realization, the mass of the most massive BH was only $\sim 2\times 10^{7}\,M_{\sun}$ at $z=7$ which was two orders of magnitude lower than the SMBH in the lowest tidal field realization. Among the simulations, the SMBH in the lowest tidal field environment had a mass an order of magnitude more than the one in the highest tidal field environment at $z=6$.

After selecting the initial conditions that best recovered the original quasar environment in BlueTides, we have run other simulations to investigate the effects of the choice of BH seeding parameters on the growth of these first massive objects. In BlueTides simulation, the BH seed mass has been chosen to be $5\times 10^{5}\,h^{-1}M_{\sun}$, which is at the high end of the predicted mass for SMBH seeds in theories. With the same constrained initial conditions, we have conducted two sets of simulations with different SMBH seed masses $M_{\bullet}^{\rm seed}=5\times 10^{3}$, $5\times 10^{4}$, and $5\times 10^{5}\,h^{-1}M_{\sun}$. Set A with different threshold halo masses such that the ratio $M_{\bullet}^{\rm seed}/M_{\rm fof}^{\rm seed}$ is fixed, while set B has a fixed halo threshold. Our simulations have suggested that the final mass of the SMBH is insensitive to the initial seed mass regardless of the choice of BH seeding parameters; the mass of SMBH in our constrained simulations has converged to $\sim 10^{9}\,M_{\sun}$ at $z=6$. In the early times at $z>10$, the growth of SMBHs varies among the simulations with different seeding scenarios; less massive seed models tend to grow slower initially unless they are seeded in more common but less massive halos so that they can merge frequently. A significant fraction of the early growth occurs in this mode in a low mass seed scenario in set A, effectively allowing the SMBH growth to catch up with that of a more massive seed. There were four SMBH mergers at $z\gtrsim 12$ for the most massive SMBH with the lowest seed mass while no mergers happened for the other two runs, suggesting that the smallest seed grows faster at earlier times when seeded in less massive halos.

The significant differences in the early merger rates provide an interesting discriminating feature for small versus large BH seed models at the early time. The space-based gravitational wave telescope LISA will open up new investigations into the dynamical processes involving SMBHs and new exciting prospects for tracing the origin, merger history of SMBHs across cosmic ages.

Acknowledgements

We acknowledge the funding from NSF ACI-1614853, NSF AST-1517593, NSF AST-1616168, NASA ATP 80NSSC18K1015 and NASA ATP 17-0123. The BlueTides simulation is run on the BlueWaters facility at the National Center for Supercomputing Applications. Some of the simulations in this work are carried out on the Stampede2 supercomputing cluster. The authors acknowledge the Texas Advanced Computing Center (TACC) at the University of Texas at Austin for providing HPC resources that have contributed to the research results reported within this paper. URL: http://www.tacc.utexas.edu.

This work has made use of the following software as well: numpy (van der Walt et al., 2011), scipy (Jones et al., 2001), matplotlib (Hunter, 2007), seaborn (Waskom et al., 2016), gaepsi2 (Feng, 2018), nbodykit (Hand & Feng, 2015).