Formation of Intermediate-Mass Black Holes through Runaway Collisions in the First Star Clusters

TL;DR

This study investigates how runaway stellar collisions in early dense star clusters can lead to the formation of intermediate-mass black holes, linking early universe conditions to present-day black hole observations.

Contribution

It demonstrates that runaway stellar collisions in first star clusters can produce IMBHs with masses up to 1900 solar masses, influenced by host halo properties.

Findings

Runaway collisions occur within a few million years in most models.

Central massive stars reach 400-1900 solar masses.

IMBH mass to cluster mass ratios are comparable to observed black hole-bulge ratios.

Abstract

We study the formation of massive black holes in the first star clusters. We first locate star-forming gas clouds in proto-galactic haloes of $\gtrsim \!10^7\,{\rm M}_{\odot}$ in cosmological hydrodynamics simulations and use them to generate the initial conditions for star clusters with masses of $\sim \!10^5\,{\rm M}_{\odot}$. We then perform a series of direct-tree hybrid $N$-body simulations to follow runaway stellar collisions in the dense star clusters. In all the cluster models except one, runaway collisions occur within a few million years, and the mass of the central, most massive star reaches $\sim \!400-1900\,{\rm M}_{\odot}$. Such very massive stars collapse to leave intermediate-mass black holes (IMBHs). The diversity of the final masses may be attributed to the differences in a few basic properties of the host haloes such as mass, central gas velocity dispersion, and mean gas density of the central core. Finally, we derive the IMBH mass to cluster mass ratios, and compare them with the observed black hole to bulge mass ratios in the present-day Universe.

Tables

Table 1. Table 1: Properties of the host haloes and the main results of our star cluster simulations with a fiducial set of model parameters, α sfe = 6.32 × 10 − 4 subscript 𝛼 sfe 6.32 superscript 10 4 \alpha_{\rm sfe}=\!6.32\times 10^{-4} , m min = 3 ​ M ⊙ subscript 𝑚 min 3 subscript M direct-product m_{\rm min}=\!3\,{\rm M}_{\odot} , m max = 100 ​ M ⊙ subscript 𝑚 max 100 subscript M direct-product m_{\rm max}=\!100\,{\rm M}_{\odot} , β = 2.35 𝛽 2.35 \beta=\!2.35 , Q = 0.5 𝑄 0.5 Q=\!0.5 and m DM = 1.87 ​ M ⊙ subscript 𝑚 DM 1.87 subscript M direct-product m_{\rm DM}=\!1.87\,{\rm M}_{\odot} . Overlines indicate that the values are averaged over 3 realizations.

	z	$R_{\mathrm{vir}}$	$M_{\mathrm{vir}}$	${\overline{M}}_{\mathrm{cl}}$	$\overline{N}$	${\overline{r}}_c$	${\overline{\rho }}_c$	${\overline{t}}_{\mathrm{rh}}$	${\overline{t}}_{\mathrm{rc}}$	${\overline{ϵ}}_{\mathrm{sfe}}$	$M_{\mathrm{DM}}$	$N_{\mathrm{DM}}$	${\overline{m}}_{\max ,f}$	${\overline{N}}_{\mathrm{coll}}$
		($\mathrm{pc}$)	(${10}^7{M}_{\odot }$)	(${10}^4{M}_{\odot }$)	(${10}^3$)	($\mathrm{pc}$)	($M_{\odot }{\mathrm{pc}}^{-3}$)	($\mathrm{Myr}$)	($\mathrm{kyr}$)	(%)	(${10}^7{M}_{\odot }$)	(${10}^7$)	($M_{\odot }$)
A	19.7	281	4.03	16.4	19.9	0.401	$6.45\times {10}^5$	19.7	528	5.91	4.79	2.56	929	11.7
B	19.6	276	2.97	13.0	15.7	0.387	$5.82\times {10}^5$	12.6	783	6.12	3.78	2.02	409	4.67
C	19.7	208	2.03	12.1	14.7	0.380	$8.36\times {10}^5$	9.67	67.1	10.1	6.60	3.53	1330	18.3
D	14.9	321	2.60	11.7	14.1	0.357	$9.75\times {10}^5$	13.6	8.93	7.16	5.67	3.03	971	13.7
E	17.1	264	1.47	4.76	5.76	0.224	$1.16\times {10}^6$	4.42	2.98	8.15	3.25	1.74	773	9.67
F	16.5	312	2.01	9.00	10.8	0.662	$7.05\times {10}^5$	15.0	3.62	8.67	5.13	2.74	1100	14.0
G	16.9	242	1.99	12.5	15.0	0.353	$1.01\times {10}^6$	10.1	4.25	9.48	4.17	2.23	1660	25.0
H	11.7	541	4.22	7.70	9.32	0.276	$1.08\times {10}^6$	10.8	0.807	5.55	5.25	2.81	964	15.0

Table 2. Table 2: Parameters of the hybrid N-body simulations.

$\eta$	$\Delta t$	$L_{\mathrm{box}}$	$n_{\mathrm{crit}}$	$\theta$	${ϵ}_{\mathrm{cl}}$	${ϵ}_{\mathrm{DM}}$
	(yr)	(pc)			(pc)	(pc)
0.11	$1.16\times {10}^3$	1024/2048	512	0.5	0	0.0313

Table 3. Table 3: Additional star cluster models with different parameters. The main difference from our fiducial model is described in the last column. The values of M DM subscript 𝑀 DM M_{\rm DM} and N DM subscript 𝑁 DM N_{\rm DM} , if DM is included, are the same as in the corresponding fiducial model and thus are not listed here. Note that the initial distribution of the stars for Model AnoDM and Arad are exactly the same as in the fiducial model A. Unless mentioned, the values are obtained by averaging over 3 realizations.

	${\overline{M}}_{\mathrm{cl},4}$	${\overline{N}}_3$	${\overline{r}}_c$	${\overline{\rho }}_{c,5}$	${\overline{t}}_{\mathrm{rh}}$	${\overline{t}}_{\mathrm{rc}}$	${\overline{ϵ}}_{\mathrm{sfe}}$	${\overline{m}}_{\max ,f}$	${\overline{N}}_{\mathrm{coll}}$	Notes
	(${10}^4{M}_{\odot }$)	(${10}^3$)	($\mathrm{pc}$)	($M_{\odot }{\mathrm{pc}}^{-3}$)	($\mathrm{Myr}$)	($\mathrm{kyr}$)	(%)	($M_{\odot }$)
A	16.4	19.9	0.401	$6.45\times {10}^5$	19.7	528	5.91	929	11.7	fiducial model (Table 1)
AnoDM	16.4	19.9	0.401	$6.45\times {10}^5$	19.7	528	5.91	915	10.7	no DM
Arad	16.4	20.0	0.400	$6.30\times {10}^5$	19.5	553	5.90	958	9.00	half radii for stars, 1 realization
Amax	16.5	18.6	0.394	$7.31\times {10}^5$	17.2	157	5.95	1510	11.0	$m_{\max }=200M_{\odot }$, 1 realization
Amin	16.6	53.3	0.401	$7.94\times {10}^5$	44.6	1040	5.99	980	12.0	$m_{\min }=1M_{\odot }$, no DM, 1 realization
Asfe1	28.5	34.5	0.509	$5.47\times {10}^5$	48.3	511	10.3	1060	14.5	${\alpha }_{\mathrm{sfe}}=1.26\times {10}^{-3}$, 2 realizations
Asfe2	8.56	10.4	0.359	$4.47\times {10}^5$	12.7	351	3.09	602	5.67	${\alpha }_{\mathrm{sfe}}=3.16\times {10}^{-4}$
Asfe3	1.75	2.10	0.361	$8.51\times {10}^4$	8.61	144	0.629	186	1.67	${\alpha }_{\mathrm{sfe}}=6.32\times {10}^{-5}$

Equations

$${ϵ}_{\mathrm{loc},i}=\max \left({\alpha }_{\mathrm{sfe}}\sqrt{\frac{n_{\mathrm{H},i}}{1\mathrm{c}{\mathrm{m}}^{-3}}}{\mathrm{e}}^{-r_i/R_{\mathrm{cl}}},1.0\right)\times \frac{m_{\mathrm{gas},i}}{{\overline{m}}_{\mathrm{s}}},$$

$${\mathbf{v}}_{\mathrm{star}}=\sqrt{\frac{Q}{T/|W|}}({\mathbf{v}}_{\mathrm{SPH}}-\overline{\mathbf{v}}),$$

$$t_{\mathrm{rh}}=\frac{0.651\mathrm{Gyr}}{\ln (\gamma N)}\frac{1{\mathrm{M}}_{\odot }}{{\overline{m}}_{\mathrm{s}}}{\left(\frac{M_{\mathrm{cl}}}{10^5{\mathrm{M}}_{\odot }}\right)}^{1/2}{\left(\frac{r_{\mathrm{h}}}{1\mathrm{pc}}\right)}^{3/2},$$

$$t_{\mathrm{rc}}=\frac{0.065{\sigma }_{\mathrm{c}}^3}{G^2{\overline{m}}_{\mathrm{s}}{\rho }_{\mathrm{c}}\ln (\gamma N)},$$

$$m=m_{\mathrm{seed}}+4\times 10^{-3}{M}_{\mathrm{cl}}{f}_{\mathrm{c}}\ln \Lambda \ln \left(\frac{t}{t_{\mathrm{cc}}}\right),$$

$$t_{\mathrm{df}}=\frac{0.186v^3}{G^{2}m\ln {\Lambda }^{\prime }}{\int }_0^{r_{\mathrm{i}}}\frac{\mathrm{d}r}{r\rho },$$

$$m_{\max ,\mathrm{f}}\propto M_{\mathrm{vir}}\ln (const.\times {\overline{\rho }}_{\mathrm{c},\mathrm{gas}}/{\sigma }_{\mathrm{c},\mathrm{gas}}^3),$$

$$m_{\max ,\mathrm{f}}=30+8\times 10^{-4}{M}_{\mathrm{cl}}\ln {\Lambda }_{\mathrm{cl}},$$

Full text

1 Introduction

The origin of supermassive black holes (SMBHs) at $z\gtrsim\!6$ remains largely unknown (e.g. Mortlock et al., 2011; Wu et al., 2015). Prompt formation of such massive objects appears implausible, and thus viable models posits formation of black hole ‘seeds’ and also their efficient growth. There are a few suggested formation mechanisms of the seed black holes. For example, very massive Population III stars, which are thought to be formed in primordial gas clouds (e.g. Susa, Hasegawa & Tominaga, 2014; Hirano et al., 2014), leave remnant BHs with masses greater than $\sim\!100\,{\rm M}_{\odot}$. While the Population III remnants appear as promising BH seeds, it is unlikely that the BH subsequently grow in mass efficiently to be as massive as the observed SMBHs by $z=\!6-7$ (e.g. Madau & Rees, 2001; Haiman & Loeb, 2001; Schneider et al., 2002; Volonteri, Haardt & Madau, 2003). Another popular physical model is the so-called direct-collapse of super-massive stars (e.g. Loeb & Rasio 1994; Oh & Haiman 2002; Bromm & Loeb 2003; Begelman, Volonteri & Rees 2006). The model is motivated to kick-start the formation of a SMBH from a BH with $\sim\!10^{5}\,{\rm M}_{\odot}$ at a sufficiently early epoch. The direct-collapse model, however, often resorts to a few peculiar conditions such as the existence of luminous galaxies nearby, and hence the overall formation rate of massive BHs remains rather uncertain.

There is yet another mechanism of the formation of massive black holes in the early universe. It has been suggested that some star clusters in the present-day Universe harbor intermediate mass black holes (IMBHs) at their centres. Such IMBHs, if they exist, are likely remnants of very massive stars that could have been formed through runaway stellar collisions (Ebisuzaki et al. 2001; Vanbeveren et al. 2009; Devecchi et al. 2010; Katz, Sijacki & Haehnelt 2015). In the context of early structure formation, Omukai, Schneider & Haiman (2008) consider dense star clusters formed in the so-called atomic-cooling haloes with virial temperature of $\gtrsim\!10^{4}\,{\rm K}$. They argue that a metal-enriched gas cloud with metallicity greater than $Z_{\rm cr}\sim\!5\times 10^{-6}\,Z_{\odot}$ can undergo dust-induced cooling and fragmentation, to produce a star cluster. Devecchi et al. (2010) study formation of very massive stars and BH seeds based on a hierarchical galaxy formation model. Their semi-analytic model predicts a large BH mass density of $\rho_{\rm seed}\sim\!2000-4000\,{\rm M}_{\odot}\,{\rm Mpc}^{-3}$ at $z\gtrsim\!6$.

Star cluster formation at high redshift has been also studied by numerical simulations (Boley et al., 2009; Trenti, Padoan & Jimenez, 2015). Bromm & Clarke (2002) perform SPH simulations of early dwarf galaxy formation and show that globular clusters of $\sim\!10^{5}\,{\rm M}_{\odot}$ are formed in small-mass haloes. Kimm et al. (2016) perform cosmological radiation-hydrodynamics simulations of globular cluster formation at $z>10$. They find that gas fragmentation induced by metal cooling drives the formation of dense star clusters of $\sim\!6\times 10^{5}\,{\rm M}_{\odot}$ with half-light radius $\lesssim\!1\,{\rm pc}$ in about ten million years.

Overall, the formation of dense star clusters in the early universe appears plausible, but there are surprisingly few simulations to date that follow the evolution of a star cluster directly. Katz, Sijacki & Haehnelt (2015) use a cosmological simulation to study the dynamical evolution of a star cluster in a metal-enriched mini-halo. With the aid of direct $N$-body simulations, they show that runaway stellar collisions occur in the cluster, to yield a very massive star with mass of $\sim\!300-1000\,{\rm M}_{\odot}$. The result provide an interesting formation path of an IMBH. Unfortunately, the cluster model of Katz, Sijacki & Haehnelt (2015) is based on only one specific mini-halo, and thus it remains unclear how rare is such IMBH formation in the early universe.

In this study, we use direct-tree hybrid $N$-body simulations to follow the evolution of a number of dense star clusters. We generate realistic initial conditions by utilizing the output of cosmological hydrodynamics simulations. We show that all reasonable cluster models except one yield a very massive star of $\gtrsim\!1000\,{\rm M}_{\odot}$ as a product of runaway collisions in a few million years.

The rest of the present paper is organised as follows. In Section 2, we describe our cosmological simulations and direct-tree hybrid N-body simulations. There, we also explain generation of the initial conditions from the cosmological simulations. In Section 3, we show the results of the hybrid $N$-body simulations. Finally in Section 4, we give summary and concluding remarks.

2 Numerical Methods

Our simulations are performed in the following three-step manner. We first run cosmological hydrodynamics simulations of early galaxy formation. We use the output to locate a number of proto-galactic haloes that host a star-forming gas cloud. We then replace the gas cloud with a dense star cluster on the assumption that numerous stars with a wide range of masses are formed in the cloud. The positions and the velocities of the stars are configured by adopting a few simple models. The realizations generated in this manner serve as the initial conditions for star cluster evolution simulations. Finally, we run direct-tree hybrid $N$-body simulations and follow the stellar dynamics.

2.1 Cosmological simulations

We use the parallel $N$-body/Smoothed Particle Hydrodynamics (SPH) code Gadget-2 (Springel, 2005), suitably modified as in Hirano et al. (2014) so that the formation of primordial gas clouds can be followed (see also Yoshida et al. 2003; Yoshida, Omukai & Hernquist 2008). The simulations are initialized at $z_{\rm ini}=\!99$ with a box size $10\,h^{-1}\,{\rm Mpc}$ using the MUSIC software (Hahn & Abel, 2011). Cosmological parameters are adopted from the latest Planck data (Ade et al., 2016, last column of their table 4). The box size is chosen to be sufficiently large to locate about 10 haloes with virial mass $\sim\!10^{7}-10^{8}\,{\rm M}_{\odot}$ at redshift $z=\!10-20$ (Reed et al., 2007). We first run a dark matter only simulation with $N=\!512^{3}$, and run a friends-of-friends halo finder to identify dark haloes at $z=\!10$. Next, we perform zoom-in simulations for the selected target haloes with a high spatial resolution. With the multi-level zoom-in technique, we achieve a mass resolution of $m_{\rm DM}\sim\!1\,{\rm M}_{\odot}$. The resolution is designated by considering that, in our star cluster simulations (Section 2.3), the DM particle mass should be smaller than the minimum stellar mass that is $3\,{\rm M}_{\odot}$ in our fiducial models (Section 2.2). We perform the zoom-in simulations including SPH particles but with radiative cooling other than atomic hydrogen cooling switched off. Molecular hydrogen cooling is disabled in order to prevent gas cloud formation in early mini-haloes. The SPH simulations are stopped when the target haloes collapse gravitationally, and the central gas density reaches $n_{\rm H}\sim\!10^{7}\,{\rm cm^{-3}}$. We run simulations for a total of eight haloes. The basic halo properties are listed in Table 1.

2.2 Generation of star cluster/DM distributions

We generate the initial conditions for our stellar dynamics simulations directly from the outputs of cosmological zoom-in simulations. For each target halo, we dump a snapshot when the central gas density is $n_{\rm H}\sim\!10^{7}\,{\rm cm^{-3}}$. The density roughly corresponds to a critical density for fragmentation for a gas cloud with metallicity of $Z\sim\!10^{-4}\,Z_{\odot}$ (see fig. 5 of Omukai, Schneider & Haiman 2008). We expect the cloud is already gravitationally unstable, and likely yields multiple stars. However, our zoom-in simulations do not resolve the formation of individual stars, and thus we employ the following simplified model to place stars within the parent gas cloud.

We sample a fraction of the SPH particles as “stars”. A selected star particle is re-assigned its mass and velocities, while its position is kept the same as that of the original SPH particle. We use the following five physical parameters to determine the sample probability and to calculate the stellar mass and velocities: local star formation efficiency (SFE) $\alpha_{\rm sfe}$, the minimum stellar mass $m_{\rm min}$, the maximum stellar mass $m_{\rm max}$, index $\beta$ of a power-law initial mass function (IMF) ${\rm d}N/{\rm d}m\propto m^{-\beta}$, and virial ratio $Q$ which is the ratio of the total stellar kinetic energy to the total stellar potential energy. The probability of replacing a SPH particle $i$ with a star particle is calculated according to the local SFE (Fujii & Portegies Zwart, 2015) given by

[TABLE]

where $r_{i}$ is the distance of the SPH particle from the point of maximum density (cloud centre), $R_{\rm cl}$ is the radius where the enclosed gas mass is equal to that of DM, $m_{{\rm gas},i}$ is SPH particle mass and $\overline{m}_{\rm s}$ is the mean stellar mass for the specified IMF. The equation (1) is based on the star-formation law $\dot{\rho}_{\rm star}\propto t_{\rm ff}^{-1}$ (Schmidt, 1959; Kennicutt, 1998), where the star formation rate $\dot{\rho}_{\rm star}$ is inversely proportional to free-fall time $t_{\rm ff}\propto n_{\rm H}^{-1/2}$. The factor $m_{{\rm gas},i}/\overline{m}_{\rm s}$ guarantees mass conservation when the gas to star conversion is performed. The exponential cutoff is applied to set a finite cluster size, but the choice of the value of $R_{\rm cl}$ is unimportant. To assign the velocity to each star particle, we re-scale the velocity of the SPH particle as

[TABLE]

where $\bm{v}_{\rm SPH}$ is the SPH particle velocity after the cloud’s bulk velocity is subtracted, $\bm{\overline{v}}=\!\sum_{\rm replaced}m_{\rm star}\bm{v}_{\rm SPH}/\sum m_{\rm star}$, $m_{\rm star}$ is stellar mass, $T=\!\sum_{\rm replaced}m_{\rm star}\bm{v}_{\rm SPH}^{2}/2$, and $W$ is potential energy of the stars. With the fiducial virial ratio of $Q=\!0.5$, we can realize a marginally stable cluster. For each sample, we generate three realizations using different random number seeds to select star particles. We thus run a total of 24 simulations for our fiducial model. We further run additional simulations to examine the effect of model parameters, as will be discussed in Section 3.2.

The dark matter distribution is kept essentially the same as in the original cosmological simulation. In practice, we split DM particles so that all the DM particles have the same mass $m_{\rm DM}$ which is chosen as the minimum dark matter particle mass in the cosmological simulation. The splitting is necessary to avoid artificial mass segregation of DM particles in the hybrid $N$-body simulations (Section 2.3). When the splitting is done, the daughter particles are randomly distributed in a sphere whose radius is the mean separation of the DM particle, and retain the same velocity. The DM halo’s bulk velocity is also subtracted, and the particle velocities are re-scaled so that the DM virial ratio becomes 0.5, as in equation (2). By doing this, we prevent the outer part of the DM halo from evaporating during our star cluster simulations.

The fiducial parameters for the initial condition generation are set as $\alpha_{\rm sfe}=\!6.32\times 10^{-4}$, $m_{\rm min}=\!3\,{\rm M}_{\odot}$, $m_{\rm max}=\!100\,{\rm M}_{\odot}$, $\beta=\!2.35$, $Q=\!0.5$ and $m_{\rm DM}=\!1.87\,{\rm M}_{\odot}$. The resulting star cluster/DM global properties are shown in Table 1. The value of $\alpha_{\rm sfe}$ is chosen such that the particle number in Model A is about $\sim\!2\times 10^{4}$. With this choice, the global star formation efficiency (column 10 in Table 1) is $\epsilon_{\rm sfe}\sim\!0.06$-$0.1$, which is approximately consistent with the value of $\sim\!0.1$ found in the hydrodynamics simulations of Kimm et al. (2016). Note that $Q$ is set to 0.5, but it does not necessarily mean that the system is in virial equilibrium since not all the star particles are bound. The Salpeter mass function is assumed throughout the present paper. We note here that stars with lower masses of $m_{\rm min}<1\,{\rm M}_{\odot}$ may exist in real clusters. Setting a smaller $m_{\rm min}$ makes the number of star particles very large and our hybrid $N$-body simulations become virtually impractical. We study the effect of varying $m_{\rm min}$, as well as the other parameters, in Section 3.2.

2.3 Direct-tree hybrid N-body simulations

The stellar dynamics simulations are performed using the hybrid $N$-body code BRIDGE (Fujii et al., 2007). The code follow the orbits of the star particles by a direct method in a dynamically consistent manner with DM particles whose motions are calculated by a tree method. In the current version, the sixth-order Hermite integrator is applied for the direct integration (Nitadori & Makino, 2008). The NINJA scheme realizes better parallel efficiency (Nitadori, Makino & Abe, 2006). The phantom-GRAPE library is used (Tanikawa et al., 2013) to speed-up the gravitational force calculation.

A pair of stars can collide and merge when its separation $d\equiv|\bm{r}_{1}-\bm{r}_{2}|$ becomes less than the sum of the stellar radii $R_{*,1}+R_{*,2}$, i.e., $d<R_{*,1}+R_{*,2}$. This merger criterion is well tested by Gaburov, Lombardi & Portegies Zwart (2010), who show that the criterion gives sufficiently accurate results when compared to the results of direct hydrodynamics simulations of stellar three-body interactions.

The stellar radius in our simulations is given by the fitting formula of Tout et al. (1996) for the non-evolving zero-age main-sequence stars with $Z=\!0.02$. The fitting, which is valid for stellar mass smaller than $100\,{\rm M}_{\odot}$, is extrapolated to larger stellar mass and thus it might possibly cause underestimation of the radii, especially for very massive stars. For example, the radii of stars with mass $100,200,500,1000\,{\rm M}_{\odot}$ from the formula of Tout et al. (1996) are $17,28,54,87\,R_{\odot}$ respectively, while those from the interior structure calculations of Ishii, Ueno & Kato (1999) for massive stars are $18,40,160,3000\,R_{\odot}$ respectively and those of Yungelson et al. (2008) are $14,27,66,129\,R_{\odot}$ respectively. Note that the stellar radius is generally smaller for stars with lower metallicity (e.g. Baraffe & El Eid, 1991; Baraffe, Heger & Woosley, 2001). Although we could technically construct a model of stellar radii for $Z=\!10^{-4}\,Z_{\odot}$ as in Katz, Sijacki & Haehnelt (2015), we have decided to use Tout’s fitting formula of the solar metallicity for simplicity. We study the effect of adopting a different model of stellar radius in (Section 3.2).

Parameters of the $N$-body simulations are presented in Table 2, where $\eta$ is an integration accuracy parameter (equation 16 in Nitadori & Makino 2008), $\Delta t$ is a tree timestep, $L_{\rm box}$ is a box size of the simulations, $n_{\rm crit}$ is the maximum group size for GRAPE calculation (Makino, 1991), $\theta$ is the tree opening angle, and $\epsilon_{\rm cl}$ and $\epsilon_{\rm DM}$ are softening lengths for star clusters and DM, respectively. The accuracy parameter $\eta$ is chosen so that the total energy errors during the simulations do not exceed $0.04\%$. The simulation boxsize $L_{\rm box}$ is set to be 1024 pc or 2048 pc, for all the particles to be encompassed during the simulations. Convergence of the results has been checked for $\Delta t$, $L_{\rm box}$ and $\epsilon_{\rm DM}$. The simulations are stopped at $t=\!3\,{\rm Myr}$, when the central massive star is supposed to end its life.

3 Results

The properties of the eight haloes are summarized in Table 1. The virial masses are $\sim\!(1.5$ - $4)\times 10^{7}\,{\rm M}_{\odot}$, and the formation epochs are $z=\!11-20$. Fig. 1 shows the gas distribution for our Halo B (top panels) and G (bottom panels), respectively. Halo B shows the features of turbulent motions in the gas distribution, whereas Halo G appears dynamically relaxed, and is approximately spherical. The generated stellar distributions are also shown in the right panels of Fig. 1. Most of the stars are located within the central few parsec region but a small number of stars are also found in the outer region.

3.1 Fiducial models

With our fiducial model parameters, the initial stellar mass in the cluster ranges from $\sim\!5\times 10^{4}\,{\rm M}_{\odot}$ to $\sim\!1.6\times 10^{5}\,{\rm M}_{\odot}$ (see Table 1).

Fig. 2 shows the time evolution of the cluster core radii. We calculate the core radii following the procedure of Casertano & Hut (1985). The core shrinks quickly within about a million years because of the short crossing time $t_{\rm cr}\equiv r_{\rm h}/\sigma\lesssim\!0.1\,{\rm Myr}$ of the cluster systems, where the half-mass radius $r_{\rm h}$ is typically $\sim\!1\,{\rm pc}$ and three-dimensional velocity dispersion $\sigma$ is $\sim\!10\,{\rm km\,s^{-1}}$.

Fig. 2 also shows the evolution of $r_{\rm h}$. In contrast to that of the core radii, $r_{\rm h}$ remains roughly constant with time. This suggests that the runaway collision is initiated by core collapse and not driven by cold collapse of the whole cluster due to the initially non-equilibrium configuration.

Core collapse time $t_{\rm cc}$, as defined by the time when the core radius has its minimum, ranges from $\lesssim\!0.1\,{\rm Myr}$ to $\sim\!2.7\,{\rm Myr}$ for our 24 simulations. Rewriting in terms of the half-mass relaxation time $t_{\rm rh}$ (Table 1) given by

[TABLE]

where $M_{\rm cl}$ is cluster mass and $\gamma$ is $\sim\!0.015$ in a system with a wide mass spectrum (Giersz & Heggie, 1996; Gürkan, Freitag & Rasio, 2004), we find that $t_{\rm cc}$ ranges from $\lesssim\!0.01t_{\rm rh}$ to $\sim\!0.2t_{\rm rh}$.

The core collapse time $t_{\rm cc}$ is expected to be proportional to the central relaxation time-scale (Gürkan, Freitag & Rasio, 2004; Fujii & Portegies Zwart, 2014, see Table 1), which is

[TABLE]

where $\sigma_{\rm c}$ is the three-dimensional velocity dispersion at a cluster centre, which we calculate using central 20 stars. The time $t_{\rm cc}$ ranges from $\sim\!2t_{\rm rc}$ to $200t_{\rm rc}$. Note that, in contrast to the result of Fujii & Portegies Zwart (2014), we find significant scatter in the ratio of $t_{\rm cc}/t_{\rm rc}$ even though $m_{\rm max}/\overline{m}_{\rm s}$ is fixed. This is likely due to the variation in the compactness of the initial stellar distribution.

The final stellar mass $m_{\rm max,f}$ and the number of collisions $N_{\rm coll}$ by $t=\!3\,{\rm Myr}$ are listed in Table 1. In all our runs, runaway collision occurs with $N_{\rm coll}\sim\!5$ - $25$ and the final stellar mass exceeds the threshold mass $\sim\!300\,{\rm M}_{\odot}$ for gravitational collapse (Heger et al. 2003; Yoon, Dierks & Langer 2012; Spera & Mapelli 2017); IMBHs are likely to be left in the clusters. Fig. 3 shows the mass of the central star that undergoes runaway collisions. We also show the mass evolution calculated based on the analytic model of (Portegies Zwart & McMillan, 2002):

[TABLE]

where $m_{\rm seed}$ is the seed mass of a star which commences runaway process, $f_{\rm c}$ is the fraction of binaries that contribute to collisions, $\ln\Lambda$ is the Coulomb logarithm, and $t_{\rm cc}$ is the core collapse time. The equation is derived by integrating the mass growth rate which is essentially the average mass increase per collision times the mean collision rate of binaries. The rate of mass accumulation through collisions in our simulations is well described by this model.

3.2 Dependence on model parameters

It is important to examine how the initial configuration and the values of our model parameters (see Section 2.2) affect the main result. To this end, we consider models with different parameters or with different set-up: models without no DM particles (AnoDM), adopting smaller stellar radii (Arad), with different mass limits for the same IMF (Amax and Amin), and with different values of the SFE parameter $\alpha_{\rm sfe}$ (Asfe1, Asfe2, and Asfe3) The assigned parameters and the resulting initial properties of the star clusters are summarized in Table 3.

Let us first compare Model A and AnoDM. Clearly, the existence of a dark matter halo has little impact on the runaway growth of the stars. This is expected because the characteristic time of dynamical friction between stars and DM is very long. We consider frictional force on a star with mass $m$ and velocity $v$ in a matter field with a density profile $\rho(r)$ (Binney & Tremaine, 2008), and use the equation of angular momentum change (Fujii & Portegies Zwart, 2014), to derive the dynamical friction time-scale $t_{\rm df}$

[TABLE]

where $\Lambda^{\prime}\simeq\!0.1N$ (Giersz & Heggie, 1994) and $r_{\rm i}$ is an initial radial position of the star. Using the inner part of the Navarro-Frenk-White dark halo density profile $\rho\propto r^{-1}$ with the normalization $\rho_{0}\sim\!7\times 10^{-20}\,{\rm g\,cm^{-3}}$ at $1\,{\rm pc}$ which is derived directly from our simulation, and using $v=\!10\,{\rm km\,s^{-1}}$, $r_{\rm i}=\!1\,{\rm pc}$, $m=\!100\,{\rm M}_{\odot}$ and $N=\!3\times 10^{7}$, we obtain $t_{\rm df}\sim\!7\,{\rm Myr}$ which is several dozen to hundreds times larger than $t_{\rm cc}$. Hence the dynamical friction by dark matter is unimportant.

Next, we compare Model A and Arad. The difference is in the model of stellar radius calculation. Model Arad assumes a smaller stellar radius for a given stellar mass. The difference causes only modest effect on the results. Specifically, the number of collisions in Model Arad is slightly smaller, but the final stellar mass $m_{\rm max,f}$ is almost the same.

A more direct effect is seen in Model Amax, where we set a larger maximum mass limit for the power-law IMF $m_{\rm max}$. The runaway growth of the central massive star is accelerated in this case. Contrastingly, setting a smaller $m_{\rm min}$ (Model Amin) does not significantly affect the runaway process. Altogether, these results are reasonable because mainly massive stars contribute to the runaway collisions through mass segregation / concentration toward the cluster centre. Note, however, that $m_{\rm min}$ is critically important in determining the total stellar mass formed within a gas cloud, hence the relative abundance of massive stars, especially when one adopts a steep or bottom-heavy IMF.

Varying the SFE parameter $\alpha_{\rm sfe}$ notably changes the results. With a higher SFE (Model Asfe1), runaway collision is accelerated because an initially denser cluster is made, while setting a smaller $\alpha_{\rm sfe}$ results in smaller $m_{\rm max,f}$ and $N_{\rm coll}$ (Asfe2 and Asfe3). In fact, in Model Asfe3, there are few collisions, and the runaway collision did not occur because of the small cluster mass (note the small value of the second term in the equation 5).

4 Summary and discussions

We have performed a suite of simulations of star cluster evolution and explored the formation of IMBHs through runaway collisions. The star cluster formation sites are located in realistic cosmological simulations. Runaway stellar collisions occur quickly within $3\,{\rm Myr}$, and very massive stars with mass of $m_{\rm max,f}\sim\!400-1900\,{\rm M}_{\odot}$ (Fig. 3) are formed, which is consistent with the results in Katz, Sijacki & Haehnelt (2015). In all our fiducial models, $m_{\rm max,f}$ exceeds $\sim\!300\,{\rm M}_{\odot}$, and thus IMBHs are expected to be left in the clusters (Heger et al., 2003). An interesting implication is that the IMBHs can seed the formation of the SMBHs observed at $z\gtrsim\!6$. The SMBHs are thought to have grown via gas accretion or merger. If near-Eddington accretion is sustained with accretion efficiency of $10\%$, a $1000\,{\rm M}_{\odot}$ seed BH grows to become as massive as $\sim\!10^{9}\,{\rm M}_{\odot}$ in $\sim\!0.6\,{\rm Gyr}$ (e.g. Sakurai, Inayoshi & Haiman, 2016). We conclude that IMBH formation in dense star clusters offers a viable mechanism for seeding SMBHs.

4.1 Correlation between the final mass and halo properties

It is intriguing to study the origin of the diversity of $m_{\rm max,f}$ (see Table 1) found in our simulations. To this end, we examine correlations between $m_{\rm max,f}$ and several halo properties. Possibly relevant physical quantities are, the halo virial mass $M_{\rm vir}$, central velocity dispersion of gas $\sigma_{\rm c,gas}$, and mean gas density of central core region $\overline{\rho}_{\rm c,gas}$. Let us naively expect the cluster mass scales with the halo mass, $M_{\rm cl}\propto M_{\rm vir}$, and the stellar density scales with the central gas density, $\rho_{\rm c}\propto\overline{\rho}_{\rm c,gas}$. From dynamical consideration, we also expect $\sigma_{\rm c}\propto\sigma_{\rm c,gas}$; we assume essentially that the bulk properties of the star cluster are determined from those of the parent gas cloud. Then, because the core collapse time $t_{\rm cc}$ scales with the relaxation time $t_{\rm rc}$, equation (5) can be rewritten in terms of the halo properties as

[TABLE]

with $t=\!3\,{\rm Myr}$. The log dependence of $\sigma_{\rm c}$ and $\rho_{\rm c}$ comes from integration of the rate of mass growth by collision (scales as $\propto t^{-1}$) over a time interval of $t_{\rm cc}\propto t_{\rm rc}\propto\sigma_{\rm c}^{3}/\rho_{\rm c}$ and $t=\!3\,{\rm Myr}$ (Portegies Zwart & McMillan, 2002). Fig. 4 shows the correlation between $m_{\rm max,f}$ and $M_{\rm cl}\ln(\rho_{\rm c}/\sigma_{\rm c}^{3})$ for our fiducial models. Although with substantial scatter, we find $m_{\rm max}$ is well correlated with the halo properties as in equation (7).

We can further derive a simple redshift dependence as follows. We use the cosmological scaling of halo properties and redshift as in, for example, the equations (18)-(20) of Ahn & Shapiro (2007); $M_{\rm vir}\propto(1+z)^{-3/2}$ with constant $T_{\rm vir}\sim\!8000\,{\rm K}$, $\overline{\rho}_{\rm c,gas}\propto(1+z)^{3}$, and $\sigma_{\rm c,gas}\propto\sigma_{\rm gas}\sim\!(GM_{\rm vir}/r_{\rm t})^{1/2}\propto M_{\rm vir}^{1/3}(1+z)^{1/2}$. Then we obtain $m_{\rm max,f}\propto(1+z)^{-3/2}\ln(1+z)$. Fig. 5 shows $m_{\rm max,f}$ against the cluster formation redshift $z$. Again, we find that the actual dependence is consistent with the above scaling (dashed line) but with substantial scatter. Despite of these interesting correlations, we argue that simulations with a significantly larger number of samples are needed to determine more accurate correlations if there is any.

4.2 Cluster mass-IMBH mass relation

There is a well-known SMBH mass-bulge mass relation (Magorrian et al., 1998; Merritt & Ferrarese, 2001) for galaxies and their central black holes. Although the objects and mass scales are very different, it is intriguing to compare the IMBH mass-cluster mass relation in our simulations. In Fig. 6, we plot the final mass $\overline{m}_{\rm max,f}$ of the runaway collision stars which are expected to collapse and leave IMBHs without significant mass loss against cluster mass $\overline{M}_{\rm cl}$. Note that our samples are located in the left-bottom portion of the figure. The lower line indicates the well-known Magorrian relation for SMBHs (equation 10 of Kormendy & Ho 2013), whereas the upper line is the black hole mass-cluster mass relation (see equation 16 and figure 3 of Portegies Zwart & McMillan, 2002). More specifically, it is given by

[TABLE]

where $\Lambda_{\rm cl}=\!\min(M_{\rm cl}/{\rm M}_{\odot},10^{6})$. Our results are indeed consistent with this. In other words, about one percent of the cluster mass contributes to the mass of the central very massive star. Overall, the first star clusters formed in early atomic-cooling haloes may be quite similar to present-day star clusters, but with smaller masses. The central black hole mass to galaxy (cluster) mass ratio is slightly higher than that for local SMBHs, perhaps reflecting the different formation mechanism.

4.3 Model uncertainties

The most significant uncertainty in our study lies in the process of generating the initial conditions for star clusters (Section 2.2). First of all, it is not trivial when a whole gas cloud should be replaced by a star cluster. We assume that a cluster is formed when the gas cloud density reaches $n_{\rm H}\sim\!10^{7}\,{\rm cm^{-3}}$. We adopt the value by noting that efficient gas cooling and fragmentation by OH cooling occur in a gas cloud with metallicity $\gtrsim\!10^{-4}\,Z_{\odot}$ (Chiaki, Yoshida & Hirano, 2016). Although protostars are formed at much higher density of $\gtrsim\!10^{20}\,{\rm cm^{-3}}$, we set a lower density threshold for “star cluster” formation to trace the global structure of the parent cloud. Related to this assumption, we also note that starburst is not strictly instantaneous but can last for over a million years (Kimm et al., 2016). Continuous star formation may promote runaway collisions by supplying newly formed stars to the cluster, or slower star formation process can actually hamper the growth of the central star. Also, a newly born star cluster does not necessarily have the virial ratio $Q=\!0.5$, and could have significantly anisotropic velocity structure. Finally, the star formation efficiency can be regulated by stellar evolution itself through various feedback effects, especially when a different IMF is assumed. Self-consistent treatment of the formation of star clusters and their evolution is beyond the scope of the present paper, but direct cosmological simulations that couple star formation and stellar dynamics as well as assembly of dark matter haloes will clarify many of the above issues.

There are also several physical processes that are not included in our hybrid $N$-body simulations. First, stellar evolution likely affects the collision rate through the evolution of the stellar radius and mass loss. As already examined in Section 3.2, the overall uncertainty in the stellar radii of main-sequence stars gives modest impact to our results. If post-main sequence evolution is considered, the collision rate can be increased during the giant phases of the stars. Stellar wind and mass loss are not expected to be significant over a short time of $\sim\!3\,{\rm Myr}$ if the stellar metallicity is as low as $\lesssim\!10^{-4}\,Z_{\odot}$ (Baraffe, Heger & Woosley 2001; Muijres et al. 2012; Katz, Sijacki & Haehnelt 2015. Second, binary evolution enhances mergers via tidal interaction (Hurley, Tout & Pols, 2002), and thus may promote runaway collisions. If this is the case, the formation of super-massive stars can occur even faster than in our simulations. Third, primordial binaries can delay core collapse of clusters by binary heating (Rasio, Fregeau & Joshi, 2001). They can also accelerate the core collapse by causing mass segregation due to the effectively increased mass (Heggie, Trenti & Hut, 2006). Finally, the actual collision conditions and outcome may not be as simple as in our model (Section 2.3). Tidal effects may enhance the rate of close encounters (Fregeau et al. 2004, and references therein), while mass loss during binary collision likely reduces the mass of the collision product (e.g. Glebbeek et al., 2009). Stellar rejuvenation by collisional mixing makes the lifetime of the collision products longer. Encountering two stars with high velocities that satisfy our merger condition could in reality just pass through each other depending on the impact parameter and the thickness of the stellar envelopes. Most of such stars in our simulations are in tight binary orbits, and thus they may eventually merge at later time. Interaction of stars and gas may enhance core collapse if dynamical friction or mass augmentation by accretion is effective. A dynamical friction time-scale of equation (6) for gas is, assuming an isothermal sphere profile $\rho=\!\rho_{0}(r/r_{0})^{-2}$ with $\rho_{0}\sim\!2\times 10^{-18}\,{\rm g\,cm^{-3}}$ at $r_{0}=\!1\,{\rm pc}$ derived from the halo data, and using $m=\!100\,{\rm M}_{\odot}$, $v=\!10\,{\rm km\,s^{-1}}$, $r_{\rm i}=\!1\,{\rm pc}$ and $\ln\Lambda^{\prime}\sim\!10$, $t_{\rm df}\sim\!0.2\,{\rm Myr}$. An accretion time-scale is defined by $t_{\rm acc}=\!m/\dot{m}_{\rm B}$, where $m$ is a typical stellar mass and $\dot{m}_{\rm B}$ is the Bondi accretion rate (e.g. Ryu et al., 2016). For $m=\!10\,{\rm M}_{\odot}$, using temperature $\sim\!10^{4}\,{\rm K}$ and density $n_{\rm H}\sim\!10^{7}\,{\rm cm^{-3}}$ as indicated in Kimm et al. (2016) for calculation of $\dot{m}_{\rm B}$, $t_{\rm acc}\sim\!2\,{\rm Myr}$. Both time-scales are within $3\,{\rm Myr}$ and the presence of gas can accelerate the core collapse of the clusters.

4.4 Fate of IMBHs in star clusters

We have successfully shown the formation of super-massive stars in dense star clusters, which undergo direct gravitational collapse to leave IMBHs. We thus suggest a promising initial process of seeding the formation of SMBHs. An important question still remains, however. It is unclear whether an IMBH in a star cluster can grow efficiently to be a SMBH in about a billion years. Either lack of gas supply or the BH’s radiation feedback can hamper its growth (e.g. Milosavljević, Couch & Bromm, 2009; Park & Ricotti, 2012). Some or even many IMBH might be left in star clusters or within galaxies in the present-day Universe (e.g. Maccarone et al., 2007; Pasham, Strohmayer & Mushotzky, 2014; Kızıltan, Baumgardt & Loeb, 2017). Further studies are warranted to examine the fate of the early IMBHs and their observational signatures.

Acknowledgements

YS is grateful to Junichiro Makino, Masaki Iwasawa, Kazumi Kashiyama, Tilman Hartwig, Kohei Inayoshi and Sunmyon Chon for fruitful discussions. The study is partly supported by Advanced Leading Graduate Course for Photon Science (YS), by Grant-in-Aid for JSPS Research Fellows (15J08816: YS), by JSPS KAKENHI Grant Number 26800108 (MSF), and by Grant-in Aid for JSPS Postdoctoral Fellowships for Research Abroad (SH). The cosmological simulations and hybrid N-body simulations are performed on Cray XC30 at Center for Computational Astrophysics, National Astronomy Observatory of Japan and XC40 at YITP in Kyoto University.