# Formation of massive seed black holes via collisions and accretion

Formation of massive seed black holes via collisions and accretion Abstract Models aiming to explain the formation of massive black hole seeds, and in particular the direct collapse scenario, face substantial difficulties. These are rooted in rather ad hoc and fine-tuned initial conditions, such as the simultaneous requirements of extremely low metallicities and strong radiation backgrounds. Here, we explore a modification of such scenarios where a massive primordial star cluster is initially produced. Subsequent stellar collisions give rise to the formation of massive (104−105 M⊙) objects. Our calculations demonstrate that the interplay among stellar dynamics, gas accretion, and protostellar evolution is particularly relevant. Gas accretion on to the protostars enhances their radii, resulting in an enhanced collisional cross-section. We show that the fraction of collisions can increase from 0.1 to 1 per cent of the initial population to about 10 per cent when compared to gas-free models or models of protostellar clusters in the local Universe. We conclude that very massive objects can form in spite of initial fragmentation, making the first massive protostellar clusters viable candidate birth places for observed supermassive black holes. methods: numerical, stars: black holes, stars: kinematics and dynamics, stars: Population III 1 INTRODUCTION More than 100 supermassive black holes have already been detected at z > 6 (Gallerani et al. 2017), with the number still continuously increasing, including 32 quasars recently discovered between z = 5.7 and 6.8 via Subaru (Matsuoka et al. 2017) and 8 z > 6 quasars via Spectral Energy Distribution (SED) model fitting of VISTA (Sutherland et al. (2015), WISE (Wright et al. 2010), and Dark Energy Survey (Flaugher et al. 2015) Year 1 observations (Reed et al. 2017). The currently most distant known quasars are at z = 7.085, i.e. 0.77 billion years after the big bang, with a mass of about 2 × 109 M⊙ (Mortlock et al. 2011), and at z = 7.54 with a mass of about 8 × 108 M⊙ (Bañados et al. 2018). Explaining the existence of these objects provides a significant challenge to our cosmological model (D’Amico et al. 2018), as the accretion at an Eddington rate requires initial seed masses of the order 104 M⊙, when realistic spin parameters and accretion disc models are taken into account (Shapiro 2005). The only solution is very massive seeds or extended periods of super-Eddington accretion or potentially also combinations of both during the formation and early growth of massive black holes. The possible scenarios leading to their formation were already outlined by Rees (1984), including the direct collapse of massive gas clouds either to a black hole or to a supermassive star, which later collapses to a black hole via general relativistic instabilities, or alternatively through the formation of a dense stellar cluster, which may either collapse into a black hole through relativistic instabilities or evolve due to stellar mergers. From these, the direct collapse model was often considered as the most promising scenario, as it can potentially produce the most massive seeds. Their formation has been proposed through low-angular momentum material in the first pre-galactic haloes (Koushiappas, Bullock & Dekel 2004), through efficient angular momentum transport by gravitational instabilities (Begelman, Volonteri & Rees 2006; Lodato & Natarajan 2006) or through bars-within-bars instabilities (Begelman & Shlosman 2009). Numerical simulations probing the formation of massive black holes via direct collapse were indicating that the latter are possible only if cooling is efficiently suppressed, for instance if the formation of molecular hydrogen is photodissociated via a strong radiation background (Bromm & Loeb 2003). Cosmological hydrodynamics simulations following such a collapse over many orders of magnitude indeed found signatures of self-gravitational instabilities in the regime of atomic hydrogen cooling (Wise, Turk & Abel 2008). An update of the pathways leading to massive black hole formation was given by Regan & Haehnelt (2009), considering the cosmological conditions that may help to keep the gas metal-free, thereby preventing strong fragmentation events. The first systematic fragmentation study in the atomic cooling regime has been pursued by Latif et al. (2013a), showing that fragmentation indeed can be strongly suppressed if cooling is only feasible via atomic hydrogen lines (see also Schleicher, Spaans & Glover 2010, for the role of the atomic line cooling transitions). The large accretion rates of 0.1 M⊙ yr−1 or more found in the simulations strongly favour protostellar models with cool atmospheres and highly extended envelopes (Hosokawa et al. 2013; Schleicher et al. 2013; Haemmerlé et al. 2018; Woods et al. 2017), implying that feedback is rather inefficient. Under such conditions, final black hole masses of 105 M⊙ can be reached based on the results of cosmological simulations (Latif et al. 2013b; Umeda et al. 2016). Such a scenario however represents the most optimistic case. In particular, one needs to investigate the strength of the radiation background to keep the gas atomic, which is typically described through the strength of the radiation background at 13.6 eV and parametrized via J21, with a value of 1 corresponding to 10−21 erg s−1 cm−3 Hz−1 sr−1. The first numerical investigations suggested a critical value of J21 ∼ 100 to prevent the formation of molecular hydrogen, while updated chemical networks and more realistic models for the radiation background (see e.g. Sugimura, Omukai & Inoue 2014; Agarwal & Khochfar 2015) have led to much larger critical values of the order 105 when applied in cosmological simulations (Latif et al. 2014, 2015, see also Agarwal et al. 2017 for a discussion of the impact of the spectral shape). Under the threshold, massive stars may potentially still form, even though the mass may be reduced by factors of 10–100 (Latif et al. 2014). In addition to molecular hydrogen line cooling, fragmentation can also be induced via metals or dust grains (Omukai et al. 2008). In the case of metal line cooling, a metallicity of 10−3 Z⊙ can already increase the cooling and trigger fragmentation within cosmological simulations (Bovino et al. 2014), while even lower metallicities of 10−5 Z⊙ are sufficient when dust cooling is considered (Schneider et al. 2006; Dopcke et al. 2013; Bovino et al. 2016; Latif et al. 2016). The need to have very strong radiation backgrounds, while keeping the gas metal-free, leads to a strong need of fine-tuning, which at best can be satisfied under very rare conditions (e.g. Agarwal et al. 2017), while in fact the need for large values of J21 provides a problem for the direct collapse scenario (Dijkstra, Ferrara & Mesinger 2014). Surprisingly, the alternative pathway of black hole formation through stellar clusters has been investigated to a lesser degree in the context of early Universe black hole formation. Analytical models by Devecchi & Volonteri (2009), Devecchi et al. (2012), and Lupi et al. (2014) predict black hole masses of 100–1000 M⊙ forming in the first stellar clusters. Katz, Sijacki & Haehnelt (2015) model the formation of a dense stellar cluster in a cosmological simulation, and show the subsequent formation of ∼1000 M⊙ black hole via N-body simulations. Similarly, Sakurai et al. (2017) showed the formation of black holes with 400–1900 M⊙ via a combination of cosmological and N-body simulations. So far, these simulations have not considered the enhanced protostellar radii in the presence of accretion (Hosokawa et al. 2013; Schleicher et al. 2013), and they also neglected the initial presence of gas in the cluster after the formation of the protostars. As we will show in this paper, their ongoing accretion events may however considerably favour collisions and mergers, and thus the formation of a central massive object. For present-day protostellar clusters, Baumgardt & Klessen (2011) have shown that 0.1–1 per cent of the protostars in the cluster may collide and help to form a particularly massive star within the cluster (see also Moeckel & Clarke 2011; Oh & Kroupa 2012; Fujii & Portegies Zwart 2013, for similar results). For star clusters at high redshift, we can however expect significantly higher collision rates. One reason is that these clusters are more dense, as the trace amount of metals will trigger cooling and fragmentation only when high densities of order 109 cm−3 are reached, thereby leading to the formation of very compact star clusters (Clark, Glover & Klessen 2008; Clark et al. 2011; Greif et al. 2011). Due to the low metallicity and the larger gas temperatures, the accretion rates will also be enhanced, thus favouring collisions through the accretion and mass growth itself, but also due to the protostellar radii that are enhanced as a result (Smith et al. 2012). We present here the first investigation which explores the formation of massive black hole seeds from a dense stellar cluster, where gas-phase effects like accretion as well as the resulting enhanced protostellar radii are taken into account. We describe our numerical methods in Section 2 and our experimental set-up in Section 3. The validation is described in Section 4. Our results are presented in Section 5, including both our reference model and an exploration of the parameter dependence. We summarize the main conclusions in Section 6. 2 NUMERICAL METHODS Modelling the early evolution of a Population III (Pop. III) protostar cluster is challenging. The main reason is the variety of physical processes that play a role, i.e. gravitational N-body dynamics, gravitational coupling between the stars and the gas, stellar growth in mass and size due to gas accretion, and stellar collisions. In this section, we describe each of these physical ingredients that go into our simulations, and how we couple them into one numerical model. We use the Astrophysical Multi-purpose Software Environment (AMUSE;1 Portegies Zwart, McMillan & Harfst 2009; Pelupessy et al. 2013; Portegies Zwart et al. 2013), which was designed for performing such multiphysics simulations as required for our study. Particularly, it has the flexibility to introduce new physical ingredients, such as a mass–radius parametrization for accreting Pop. III protostars and to couple it to existing N-body codes and background potentials (e.g. Martínez-Barbosa et al. 2016; Boekholt et al. 2017). 2.1 Initial conditions and dynamics Our astrophysical system under investigation consists of Pop. III protostars embedded in their natal gas cloud. We will assume a simplified initial condition: the protostars and the gas are distributed equally, and they both follow the commonly used Plummer distribution (Plummer 1911). In order to have a well-defined size of the cluster and to ensure that each protostar starts out within the gas cloud, we introduce a cut-off radius after which the density is set to zero. We set this radius to five times the Plummer radius so that the cluster remains stable. The parameters specifying the initial conditions are then the total gas mass, Mg; the cut-off radius, Rg; and the number of protostars, N. The initial mass of the protostars is set to m0 = 0.1 M⊙. In Appendix B, we confirm that the specific choice of initial mass does not change the results much, as long as we start out with a gas-dominated system. Complicating factors such as a flattened distribution, cluster rotation, and an initial binary fraction are not taken into account in the this study. Our main aim is to explore the complex interplay between dynamics and accretion that will lead to the collisional growth of a massive object. To model the star–star gravitational interactions, we use the N-body code ph4 (e.g. McMillan & Hut 1996, , section 3.2), which is a fourth-order Hermite algorithm in combination with the time-symmetric integration scheme of Hut, Makino & McMillan (1995). We simplify the gravitational dynamics of the gas cloud by an analytical background potential. This potential is coupled to the stars using the BRIDGE method (Fujii et al. 2007), so that the stars experience both the gravitational force from each other as well as from the gas. Especially at the start of the simulation, when the protostars are still low mass, their orbits will be completely determined by the dominant background potential. 2.2 Gas accretion models The initially low mass Pop. III protostars will gain mass by accreting from the gas reservoir. Since the accretion rate might vary with cluster environment and cluster evolution, we define six different accretion models based on: (in)finite gas reservoir, position-(in)dependent accretion rate, and time-(in)dependent accretion rate. An infinite gas reservoir resembles a system that is constantly being fed fresh gas. This prolongs the gas-dominated phase and the time-scale on which the protostars can accrete gas. This is contrary to the finite gas reservoir models, where the gas will eventually run out, and the protostars will stop accreting. For the position-dependent models, we set the accretion rate proportional to the local gas density. In this way, the protostars in the core accrete at a higher rate than protostars in the halo. Such a system naturally produces a range of stellar masses and radii. In order to compare the position-(in)dependent models, we make sure that the cumulative accretion rate is initially equal. A time dependence to the accretion rate is introduced if we set the accretion rate proportional to the gas density, which in turn decreases in time for the models with a finite gas reservoir. The main effect of a decreasing accretion rate is that the stars will migrate to lower mass–radius tracks, and thus have a decreasing collisional cross-section. The different combinations of the accretion properties described above define the six different accretion models presented in Table 1. An illustration of the time evolution of the accretion rate in the different models is presented in Appendix B. Table 1. Six different gas accretion models for Pop. III protostars embedded in their natal gas cloud. Model  Gas reservoir  Position-dependent  Time-dependent      accretion model  accretion model  1  Infinite  No  No  2  Infinite  Yes  No  3  Finite  No  No  4  Finite  Yes  No  5  Finite  No  Yes  6  Finite  Yes  Yes  Model  Gas reservoir  Position-dependent  Time-dependent      accretion model  accretion model  1  Infinite  No  No  2  Infinite  Yes  No  3  Finite  No  No  4  Finite  Yes  No  5  Finite  No  Yes  6  Finite  Yes  Yes  View Large 2.3 Mass–radius evolution The mass–radius evolution of accreting Pop. III protostars is very uncertain, and so we investigate two different sets of models to have a handle on the uncertainties introduced by different protostellar models. We base our mass–radius parametrizations on two different studies: Hosokawa et al. (2012, fig. 5) and Haemmerlé et al. (2018, fig. 2). Both of these studies performed detailed stellar evolution calculations using codes by Omukai & Palla (2003), Hosokawa & Omukai (2009), Hosokawa, Yorke & Omukai (2010), and GENEVA (Eggenberger et al. 2008), respectively. The radius of a protostar is completely determined by its mass, m, and accretion rate, $$\dot{m}$$. At every time-step in our simulation, we keep track of these two quantities and update the radius of the protostar. For values of the accretion rate in between the parametrized mass–radius tracks, we use interpolation between the two nearest tracks in log-space. In Fig. 1, we present our approximate parametrization of the mass–radius tracks of accreting Pop. III protostars, based on the detailed calculations of Hosokawa et al. (2012, fig. 5) and Haemmerlé et al. (2018, fig. 2). For the higher accretion rates, the models of Haemmerlé et al. (2018) produce somewhat smaller radii, but they show the same behaviour at large masses. We use both models to test the sensitivity of the formation of massive objects to the underlying mass–radius parametrization. The analytical form of the parametrization is given in Appendix A. Figure 1. View largeDownload slide Two parametrizations of the mass–radius evolution of accreting Pop. III protostars based on Hosokawa, Omukai & Yorke (2012, fig. 5) (left-hand panel) and Haemmerlé et al. (2018, fig. 2) (right-hand panel). Figure 1. View largeDownload slide Two parametrizations of the mass–radius evolution of accreting Pop. III protostars based on Hosokawa, Omukai & Yorke (2012, fig. 5) (left-hand panel) and Haemmerlé et al. (2018, fig. 2) (right-hand panel). 2.4 Stellar collisions We adopt the commonly used ‘sticky-sphere’ approximation to treat collisions between protostars. Whenever the distance between two protostars is less than the sum of their radii, we replace the two protostars by a single object at their centre of mass. We assume that during the collision the total mass is conserved, and the new radius is determined by the mass–radius parametrization. In case the accretion rate is very low, i.e. < 10 − 6 M⊙ yr − 1, then the new radius is determined by conserving the density of the primary star. 3 EXPERIMENTAL SET-UP The input parameters specifying a simulation and their range in values are as follows: Gas cloud mass Mg = 104, 3 × 104, 105, 3 × 105, 106 M⊙ Gas cloud radius Rg = 0.01, 0.03, 0.1, 0.3, 1.0 pc Number of protostars N = 64, 128, 256, 512, 1024 Average accretion rate $$\dot{m}$$ = 0.001, 0.003, 0.01, 0.03, 0.1, 0.3 $$\rm {M}_{\odot }\,\rm {yr}^{-1}$$. We also vary the accretion model (see Section 2.2) and the mass–radius parametrization (see Section 2.3). We note that the adopted protostellar accretion rates are high compared to present-day star formation. However, they are realistic for the primordial case. Similar to current star formation, the protostellar accretion process in Pop. III clusters is not regulated via Bondi–Hoyle–Littleton accretion, but rather it has been shown that gravitational instabilities in the gas are driving fragmentation as well as the accretion process on to the fragments. This leads to the typical range of accretion rates under Pop. III conditions that we adopted, see, e.g. Clark et al. (2011), Greif et al. (2011), Smith et al. (2012), Latif et al. (2013b), Latif & Schleicher (2015a), Latif & Schleicher (2015b), Hirano et al. (2015), Regan, Johansson & Wise (2016), Hosokawa et al. (2016), and Latif et al. (2016). We define a standard set of parameters as: Mg = 105 M⊙, Rg = 0.1 pc, N = 256, $$\dot{m} = 0.03\,\rm {M_{\odot }\,yr^{-1}}$$, and mass–radius parametrization based on Hosokawa et al. (2012). This choice of parameter values reflects that we are particularly interested in very massive Pop. III protostar clusters and the formation of very massive objects. The initial crossing time of our systems is given by   $$T_{\rm {cross,0}} = 853\,\left( \frac{10^5\,\rm {M}_{\odot }}{M_{\rm {g}}} \right)^{\frac{1}{2}} \left(\frac{R_{\rm {g}}}{0.1\,\rm { pc}}\right)^{\frac{3}{2}}\,\rm {yr},$$ (1)where we used Tcross,0 = 2Rv/σ, with Rv the virial radius given by Rv = 16Rg/15π (using the definition of our gas cloud truncation radius), σ the velocity dispersion estimated by $$\sigma = \sqrt{G M_{\rm {g}} / 2 R_{\rm {v}}}$$, and G the gravitational constant. We consider a simulation finished if most of the collisions have occurred. We determine this by keeping track of the average collision rate as   $$R_{\rm {av}}\left( t \right) = \frac{N_{\rm {col}}\left( t \right)}{t_{\rm {last\,collision}}},$$ (2)and an upper limit of the current collision rate as   $$R\left( t \right) = \frac{1}{t-t_{\rm {last\,collision}}}.$$ (3)If the ratio of R/Rav < 0.015, we stop the simulation. This criterion was chosen to make sure that the majority of collisions have occurred, while also limiting the duration of the simulation. We also stop the simulation if no collision has occurred in the last million years. For each simulation, we store regular snapshots of the full phase space information, allowing us to retrace the mass–radius evolution and calculate collision rates. In Tables 2 and 3, we provide an overview of our simulations and their input parameters, together with several statistics describing the outcome of the simulations, such as the maximum mass and maximum collision rate. An estimate of the measurement uncertainties can be found in Table 4. Table 2. Overview of the simulations. The input parameters are: gas cloud mass Mg, gas cloud radius Rg, number of protostars N, and the average accretion rate $$\dot{m}$$. Statistics describing the output are: final mass of the most massive object Mm and total number of collisions Nc. The number in the subscript denotes the accretion model (see Table 1). Mg  Rg  N  $$\dot{m}$$  Mm,1  Mm,2  Mm,3  Mm,4  Mm,5  Mm,6  Nc,1  Nc,2  Nc,3  Nc,4  Nc,5  Nc,6  M⊙  pc    M⊙ yr − 1  M⊙  M⊙  M⊙  M⊙  M⊙  M⊙              105  0.1  256  0.03  128 698  122 888  73 463  88 127  13 167  81 490  238  223  196  197  52  186  104  0.1  256  0.03  52 567  60 425  4069  7997  627  986  242  231  110  124  15  13  3 × 104  0.1  256  0.03  68 380  91 207  18 737  26 004  2240  14 670  239  234  162  182  19  81  3 × 105  0.1  256  0.03  165 820  133 309  138 241  179 454  111 502  131 137  231  196  216  198  166  205  106  0.1  256  0.03  179 746  172 681  178 408  172 964  158 720  142 053  230  191  219  191  209  187  105  0.01  256  0.03  10 707  9449  11 734  9343  12 557  12 406  239  203  239  205  239  210  105  0.03  256  0.03  48 319  43 226  32 766  33 458  36 281  31 088  239  211  226  204  209  191  105  0.3  256  0.03  253 063  362 374  53 793  79 250  5796  47 880  230  230  140  150  16  81  105  1.0  256  0.03  701 633  1 020 974  38 607  45 433  2345  4092  223  207  100  59  7  5  105  0.1  64  0.03  48 671  48 945  32 849  56 258  38 327  44 628  53  48  40  44  42  40  105  0.1  128  0.03  87 493  85 516  77 464  73 564  16 884  42 905  114  88  105  81  43  69  105  0.1  512  0.03  155 664  181 666  68 136  86 562  7507  70 803  500  439  371  362  52  319  105  0.1  1024  0.03  226 161  275 657  74 954  91 700  6101  57 747  1004  985  792  752  68  455  105  0.1  256  0.001  23 649  36 107  11 059  28 385  10 773  19 905  62  114  31  91  35  68  105  0.1  256  0.003  58 318  65 221  24 661  47 494  11 610  45 009  104  188  68  149  43  133  105  0.1  256  0.01  97 024  94 974  35 739  79 289  8204  54 539  174  210  103  191  26  171  105  0.1  256  0.1  195 809  289 809  78 644  88 052  20 661  79 274  245  245  206  192  77  170  105  0.1  256  0.3  244 229  309 201  53 819  89 040  13 292  73 790  245  250  149  174  51  136  Mg  Rg  N  $$\dot{m}$$  Mm,1  Mm,2  Mm,3  Mm,4  Mm,5  Mm,6  Nc,1  Nc,2  Nc,3  Nc,4  Nc,5  Nc,6  M⊙  pc    M⊙ yr − 1  M⊙  M⊙  M⊙  M⊙  M⊙  M⊙              105  0.1  256  0.03  128 698  122 888  73 463  88 127  13 167  81 490  238  223  196  197  52  186  104  0.1  256  0.03  52 567  60 425  4069  7997  627  986  242  231  110  124  15  13  3 × 104  0.1  256  0.03  68 380  91 207  18 737  26 004  2240  14 670  239  234  162  182  19  81  3 × 105  0.1  256  0.03  165 820  133 309  138 241  179 454  111 502  131 137  231  196  216  198  166  205  106  0.1  256  0.03  179 746  172 681  178 408  172 964  158 720  142 053  230  191  219  191  209  187  105  0.01  256  0.03  10 707  9449  11 734  9343  12 557  12 406  239  203  239  205  239  210  105  0.03  256  0.03  48 319  43 226  32 766  33 458  36 281  31 088  239  211  226  204  209  191  105  0.3  256  0.03  253 063  362 374  53 793  79 250  5796  47 880  230  230  140  150  16  81  105  1.0  256  0.03  701 633  1 020 974  38 607  45 433  2345  4092  223  207  100  59  7  5  105  0.1  64  0.03  48 671  48 945  32 849  56 258  38 327  44 628  53  48  40  44  42  40  105  0.1  128  0.03  87 493  85 516  77 464  73 564  16 884  42 905  114  88  105  81  43  69  105  0.1  512  0.03  155 664  181 666  68 136  86 562  7507  70 803  500  439  371  362  52  319  105  0.1  1024  0.03  226 161  275 657  74 954  91 700  6101  57 747  1004  985  792  752  68  455  105  0.1  256  0.001  23 649  36 107  11 059  28 385  10 773  19 905  62  114  31  91  35  68  105  0.1  256  0.003  58 318  65 221  24 661  47 494  11 610  45 009  104  188  68  149  43  133  105  0.1  256  0.01  97 024  94 974  35 739  79 289  8204  54 539  174  210  103  191  26  171  105  0.1  256  0.1  195 809  289 809  78 644  88 052  20 661  79 274  245  245  206  192  77  170  105  0.1  256  0.3  244 229  309 201  53 819  89 040  13 292  73 790  245  250  149  174  51  136  View Large Table 3. Overview of the simulations. Continuing Table 2 presenting the maximum collision rate Rm, and the average collision rate Rav. Mg  Rg  N  $$\dot{m}$$  Rm,1  Rm,2  Rm,3  Rm,4  Rm,5  Rm,6  Rav,1  Rav,2  Rav,3  Rav,4  Rav,5  Rav,6  M⊙  pc    M⊙ yr − 1  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  105  0.1  256  0.03  21  22  14  19  4  14  1.60  3.37  0.32  1.29  0.36  0.30  104  0.1  256  0.03  145  125  8  31  1  8  14.35  20.17  0.20  0.08  0.02  0.43  3 × 104  0.1  256  0.03  55  69  16  62  2  43  9.66  10.52  0.15  0.14  0.05  4.12  3 × 105  0.1  256  0.03  13  9  8  17  6  11  0.87  1.35  0.78  0.32  0.42  0.50  106  0.1  256  0.03  7  5  7  7  6  6  0.39  0.26  0.38  0.26  0.35  0.45  105  0.01  256  0.03  11  15  14  13  12  9  1.47  0.69  0.65  0.78  0.55  0.28  105  0.03  256  0.03  15  13  21  17  13  12  0.67  0.48  1.58  1.01  0.30  0.71  105  0.3  256  0.03  47  55  9  49  2  30  5.12  4.11  0.16  0.17  0.04  4.22  105  1.0  256  0.03  83  71  7  23  2  3  7.16  5.58  0.30  0.21  0.10  0.79  105  0.1  64  0.03  5  5  6  9  3  3  0.45  0.53  0.48  0.38  0.26  0.25  105  0.1  128  0.03  10  9  7  13  5  7  0.77  0.38  0.90  0.64  0.52  0.66  105  0.1  512  0.03  81  57  48  59  2  54  5.66  3.86  2.36  2.91  0.28  2.18  105  0.1  1024  0.03  190  185  52  232  3  123  14.71  7.53  2.06  3.10  0.42  4.86  105  0.1  256  0.001  2  4  2  3  1  4  0.07  0.13  0.04  0.11  0.01  0.11  105  0.1  256  0.003  3  7  4  5  3  6  0.13  0.33  0.01  0.43  0.07  0.14  105  0.1  256  0.01  7  12  7  12  3  8  0.57  0.68  0.08  0.57  0.03  0.59  105  0.1  256  0.1  42  62  28  44  12  22  6.12  1.93  0.58  1.78  0.66  1.09  105  0.1  256  0.3  80  114  23  72  10  47  14.11  17.55  1.93  0.92  0.77  7.44  Mg  Rg  N  $$\dot{m}$$  Rm,1  Rm,2  Rm,3  Rm,4  Rm,5  Rm,6  Rav,1  Rav,2  Rav,3  Rav,4  Rav,5  Rav,6  M⊙  pc    M⊙ yr − 1  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  105  0.1  256  0.03  21  22  14  19  4  14  1.60  3.37  0.32  1.29  0.36  0.30  104  0.1  256  0.03  145  125  8  31  1  8  14.35  20.17  0.20  0.08  0.02  0.43  3 × 104  0.1  256  0.03  55  69  16  62  2  43  9.66  10.52  0.15  0.14  0.05  4.12  3 × 105  0.1  256  0.03  13  9  8  17  6  11  0.87  1.35  0.78  0.32  0.42  0.50  106  0.1  256  0.03  7  5  7  7  6  6  0.39  0.26  0.38  0.26  0.35  0.45  105  0.01  256  0.03  11  15  14  13  12  9  1.47  0.69  0.65  0.78  0.55  0.28  105  0.03  256  0.03  15  13  21  17  13  12  0.67  0.48  1.58  1.01  0.30  0.71  105  0.3  256  0.03  47  55  9  49  2  30  5.12  4.11  0.16  0.17  0.04  4.22  105  1.0  256  0.03  83  71  7  23  2  3  7.16  5.58  0.30  0.21  0.10  0.79  105  0.1  64  0.03  5  5  6  9  3  3  0.45  0.53  0.48  0.38  0.26  0.25  105  0.1  128  0.03  10  9  7  13  5  7  0.77  0.38  0.90  0.64  0.52  0.66  105  0.1  512  0.03  81  57  48  59  2  54  5.66  3.86  2.36  2.91  0.28  2.18  105  0.1  1024  0.03  190  185  52  232  3  123  14.71  7.53  2.06  3.10  0.42  4.86  105  0.1  256  0.001  2  4  2  3  1  4  0.07  0.13  0.04  0.11  0.01  0.11  105  0.1  256  0.003  3  7  4  5  3  6  0.13  0.33  0.01  0.43  0.07  0.14  105  0.1  256  0.01  7  12  7  12  3  8  0.57  0.68  0.08  0.57  0.03  0.59  105  0.1  256  0.1  42  62  28  44  12  22  6.12  1.93  0.58  1.78  0.66  1.09  105  0.1  256  0.3  80  114  23  72  10  47  14.11  17.55  1.93  0.92  0.77  7.44  View Large Table 4. Comparing two different mass–radius parametrizations: Hosokawa et al. (2012) (subscript 1) and Haemmerlé et al. (2018) (subscript 2) (defined in Section 2.3). We compare their influence on the formation of a massive object: its final mass Mmax, number of collisions Ncol, and maximum collision rate Rmax. We compare runs with the exact same initial realization of the cluster, and calculate the average and standard deviation over 10 simulations. Models 2 and 4 give consistent results, which are the position-dependent accretion models where a dominant object is formed in the core. For the other models we observe a statistical difference due to the more conservative radii of Haemmerlé et al. (2018). This difference does not change the conclusion that a single massive object is formed due to a high collision rate. Model  Mmax,1  δMmax,1  Mmax,2  δMmax,2  Ncol,1  δNcol,1  Ncol,2  δNcol,2  Rmax,1  δRmax,1  Rmax,2  δRmax,2    M⊙  M⊙  M⊙  M⊙          $$\rm {t_{cross,0}^{-1}}$$  $$\rm {t_{cross,0}^{-1}}$$  $$\rm {t_{cross,0}^{-1}}$$  $$\rm {t_{cross,0}^{-1}}$$  1  120 529  13 127  173 130  13 856  235.3  6.67  182.5  6.17  23.3  5.48  14.2  3.01  2  139 439  16 503  160 317  16 480  228.6  5.30  225.4  3.60  28.9  4.12  25.7  6.60  3  74 613  3777  33 329  4637  201.5  6.00  87.4  10.60  16.8  3.68  5.4  1.43  4  84 285  4581  83 162  5962  189.1  13.74  181.1  15.88  21.7  3.47  19.6  4.58  5  13 381  1791  9238  1174  47.3  4.79  25.3  3.30  3.7  0.82  2.3  1.06  6  71 523  3729  56 040  2352  183.2  5.71  121.7  9.44  18.8  4.08  12.8  2.20  Model  Mmax,1  δMmax,1  Mmax,2  δMmax,2  Ncol,1  δNcol,1  Ncol,2  δNcol,2  Rmax,1  δRmax,1  Rmax,2  δRmax,2    M⊙  M⊙  M⊙  M⊙          $$\rm {t_{cross,0}^{-1}}$$  $$\rm {t_{cross,0}^{-1}}$$  $$\rm {t_{cross,0}^{-1}}$$  $$\rm {t_{cross,0}^{-1}}$$  1  120 529  13 127  173 130  13 856  235.3  6.67  182.5  6.17  23.3  5.48  14.2  3.01  2  139 439  16 503  160 317  16 480  228.6  5.30  225.4  3.60  28.9  4.12  25.7  6.60  3  74 613  3777  33 329  4637  201.5  6.00  87.4  10.60  16.8  3.68  5.4  1.43  4  84 285  4581  83 162  5962  189.1  13.74  181.1  15.88  21.7  3.47  19.6  4.58  5  13 381  1791  9238  1174  47.3  4.79  25.3  3.30  3.7  0.82  2.3  1.06  6  71 523  3729  56 040  2352  183.2  5.71  121.7  9.44  18.8  4.08  12.8  2.20  View Large In the next section, we provide several validation tests of our experimental set-up. We show that the initial condition is a stable Plummer sphere distribution (Fig. B1), and we validate both the correct implementation of the six different accretion models and the mass–radius parametrization based on Hosokawa et al. (2012). 4 VALIDATION OF THE NUMERICAL METHOD We describe the numerical methods used in this study in Section 2. Here, we present validation experiments of our simulation set-up and numerical implementation. We have performed a validation experiment of the initial conditions and dynamics by choosing the following parameters: Mg = 105 M⊙, Rg = 0.1 pc, and N = 256. We evolved the system for a time T = 105 yr and note that accretion and collisions are not included yet. In Fig. B1, we confirm that the 10, 50, and 90 per cent Lagrangian radii follow those of a Plummer sphere, with a slight discrepancy in the 90 per cent radius due to the truncation that we have introduced. We also note that all the protostars remain within the gas cloud, i.e. there are no escapers due to the initial condition, and the velocity dispersion of the protostars is close to the analytical value of ∼ 80 km s − 1. We have constructed a stable star cluster which will only alter its configuration due to gas accretion and stellar collisions. We have performed a numerical validation experiment by evolving the same initial condition as above, but this time only with gas accretion, i.e. dynamics is turned off. We set the initial accretion rate $$\dot{m} = 0.03\,\rm {M_{\odot }\,yr^{-1}}$$, i.e. here and in the rest of this study $$\dot{m}$$ refers to the average accretion rate per star over the whole cluster. In Fig. B2, we present the time evolution of the total star and gas mass (top row) and the time evolution of the average accretion rate (bottom row). We note that the position-dependent and -independent accretion models are consistent on average (their curves overlie). Together all these different models cover a variety of physical regimes. In the same validation experiment, we also kept track of the mass–radius evolution of the protostars, which we present in Fig. 2 (for the model of Hosokawa et al. 2012). There we clearly observe the difference between the position-independent (top row) and -dependent (bottom row) models. The latter produce a spectrum of masses and radii. We note that in the time-dependent models (5 and 6), the accretion rate decreases in time, and as a result the protostars will migrate to lower mass–radius tracks. After an initial phase of growth these protostars will eventually shrink, which is expected to decrease the collision rate. Figure 2. View largeDownload slide Validation of the mass–radius parametrization implementation. We consider a Plummer distribution of protostars, which are accreting gas, resulting in evolving masses and radii. The protostars do not move in this validation experiment, which explains the regularity in the distribution of mass–radius tracks. We confirm that for the position-independent accretion models (top row), the protostars follow the mass–radius parametrization based on Hosokawa et al. 2012. For the time-dependent accretion model (right-hand panel), the accretion rate decreases in time, so that the protostars migrate to mass–radius tracks of lower accretion rate, causing the protostars to eventually shrink again. For the position-dependent accretion models (bottom row), we confirm a correct interpolation of the mass–radius tracks in log-space, and the production of a range of masses and radii. Figure 2. View largeDownload slide Validation of the mass–radius parametrization implementation. We consider a Plummer distribution of protostars, which are accreting gas, resulting in evolving masses and radii. The protostars do not move in this validation experiment, which explains the regularity in the distribution of mass–radius tracks. We confirm that for the position-independent accretion models (top row), the protostars follow the mass–radius parametrization based on Hosokawa et al. 2012. For the time-dependent accretion model (right-hand panel), the accretion rate decreases in time, so that the protostars migrate to mass–radius tracks of lower accretion rate, causing the protostars to eventually shrink again. For the position-dependent accretion models (bottom row), we confirm a correct interpolation of the mass–radius tracks in log-space, and the production of a range of masses and radii. 5 RESULTS We start by analysing the simulations with the standard set of parameters defined in Section 3. We show a proof of concept that many accretion-induced collisions can lead to the formation of a single massive object. Next, we determine the sensitivity of this result to variations of the input parameters. 5.1 Simulations with the standard set of parameters Our simulations start with a cluster of single protostars, that is, there are no primordial binaries in the system. In time, collisions between protostars will occur and will result in collision products. We distinguish between the most massive collision product and other less massive ones. Finally, strong dynamical encounters among protostars can eject them from the cluster. We have thus defined four categories to which a protostar can belong: Single protostar, i.e. not a collision product Part of the most massive collision product Part of a less massive collision product Escaper. We note that binaries and higher order multiple systems form in our simulations, but we only categorize the individual stellar components. In Fig. 3 Figure 3. View largeDownload slide Time evolution of the fraction of stars belonging to four categories: (1) escapers (green); (2) stars that have collided with and are part of the most massive star in the system (blue, ‘Max. mass’); (3) stars that are part of other collision products (red, ‘Coll. prod.’); and (4) single stars (orange). We indicate the moment of transition to a stellar mass dominated system (see Section 2.2) with a vertical dashed line. Initially, all the stars are single stars in the system. As the stars accrete gas, they become larger and will dynamically interact with the other stars. We generally observe a decrease in the fraction of single stars, and an increase in collision products, and at late times a small fraction of escapers. Figure 3. View largeDownload slide Time evolution of the fraction of stars belonging to four categories: (1) escapers (green); (2) stars that have collided with and are part of the most massive star in the system (blue, ‘Max. mass’); (3) stars that are part of other collision products (red, ‘Coll. prod.’); and (4) single stars (orange). We indicate the moment of transition to a stellar mass dominated system (see Section 2.2) with a vertical dashed line. Initially, all the stars are single stars in the system. As the stars accrete gas, they become larger and will dynamically interact with the other stars. We generally observe a decrease in the fraction of single stars, and an increase in collision products, and at late times a small fraction of escapers. , we present the time evolution of the fraction of protostars belonging to each of these four categories. We observe that the majority of protostars end up in the final most massive object, which can be understood intuitively as follows. As a result of the gas accretion, the protostellar radii are substantially increasing for two main reasons: first, the overall increase of the protostellar mass, which naturally provides larger radii; and secondly, the rather high accretion rates in the primordial environment, which lead to rather extended envelopes that may exceed several 1000 solar radii. As a consequence, the collision probability in this environment is highly enhanced once accretion is switched on, and once the effect of accretion on to the protostellar radii is being considered. In particular for the infinite gas reservoir models (1 and 2), we find that 80–90 per cent of the protostars end up in a single object. For finite but large gas reservoirs (as expected in the large atomic cooling haloes), the fraction is somewhat reduced to 60–70 per cent in models 3, 4, and 6, but still quite significant. In model 5, the fraction reduces to about 10 per cent as due to the uniform accretion the central object in that cluster is less pronounced, and due to the time-dependent accretion, the protostellar radii shrink again with decreasing gas reservoir. The latter represents an enhancement by roughly a factor of 10 compared to the corresponding gas-free models or models that assumed smaller protostellar radii (Baumgardt & Klessen 2011; Reinoso et al. 2018). Comparing the radially independent accretion models (top row of Fig. 3) and the radially dependent models (bottom row), we observe a higher fraction of less massive collision products for the radially independent accretion models. Initially, this fraction is larger than the fraction of stars in the most massive collision product, implying that collisions tend to occur throughout the cluster between different pairs of stars. The radial-independent models also produce a larger fraction of escaping protostars, implying that these systems produce stronger dynamical encounters. For the radially dependent accretion models, most stars tend to immediately collide with the central protostar that is accreting fastest. Next, we aim to better understand the time evolution of the collision rate of protostars. In Figs 4 and 5, we present the collision rate (bottom panels) and correlate it with the total star and gas mass (top panels), the fraction of stars belonging to the four categories defined earlier (second panels) and the radius of the most massive protostar and the average radius of the remaining protostars (third panels). Figure 4. View largeDownload slide Correlating the time evolution of the collision and escape rate (bottom panel) with: total star and gas mass (top panel), fraction of stars belonging to the same four categories as in Fig. 3 (second panel), and maximum stellar radius and average stellar radius of the remaining stars (third panel). We show these results for the different accretion models (per column, see also next figure). Figure 4. View largeDownload slide Correlating the time evolution of the collision and escape rate (bottom panel) with: total star and gas mass (top panel), fraction of stars belonging to the same four categories as in Fig. 3 (second panel), and maximum stellar radius and average stellar radius of the remaining stars (third panel). We show these results for the different accretion models (per column, see also next figure). Figure 5. View largeDownload slide Same as the previous figure but for the position-dependent accretion models. Figure 5. View largeDownload slide Same as the previous figure but for the position-dependent accretion models. There is a small delay time for the collision rate to start to increase, which corresponds to the time it takes for the protostars to grow in mass and radius, and for the total stellar mass to become significant compared to the total gas mass. Then the collision rate increases rapidly and reaches a peak value. This rapid growth is fuelled by accretion, which results in the protostars reaching larger sizes and which makes the system stellar-mass dominated, resulting in dynamical encounters and collisions. Values for the maximum and average collision rates can be found in Table 3, which together specify a range of typical collision rates in our simulations. We find that the maximum collision rate can be very high, i.e. a collisional fraction per crossing time up to about 8 per cent (see left-hand panels in Figs 4 and 5), but that the average collision rate is more consistent with previous studies. For example, Baumgardt & Klessen (2011, fig. 4) measure a collisional fraction of 0.1–1 per cent, and this result is confirmed by Reinoso et al. (2018) for gas-free systems. Portegies Zwart et al. (2004) estimate the collision rate in the cluster MGG-11 to be 10–100 collisions within the first 3 Myr, or since the crossing time is about 105 yr (McCrady, Gilbert & Graham 2003), ∼0.3–3 collisions per crossing time, which is comparable to the average rates in Table 3. Our results show that accretion-induced collisions in massive Pop. III protostar clusters can increase the peak collision rate by an order of magnitude or even more. The saturation and eventual decrease of the collision rate have several causes. For the infinite gas reservoir models, the collision rate decreases again due to the sparsity of collision partners left in the system. For the finite gas reservoir models, the collision rate decreases due to a lack of gas, i.e. the stars do not expand anymore due to accretion. In this regime, we observe that the fraction of escaping stars starts to increase, which implies strong dynamical encounters. After the initial evolution with a high accretion-induced collision rate, the system experiences a transition into a regime with a dynamics-dominated collision rate, which is consequently much lower. For the time-dependent accretion models, the collision rate also decreases due to the shrinking of the protostars. We note however that even though the collision rate decreases, the collisions that do occur can be quite massive if both collision partners are collision products themselves. In Fig. 6, we plot the maximum mass in the system as a function of time. We observe that all models produce a very massive object with a mass of order 104 − 5 M⊙. The fastest growth in mass occurs between 5 and 20 crossing times, which is synchronous to the moment of highest collision rate (see Figs 4 and 5, bottom panels). For model 5, where accretion-induced collisions have the least effect, there is a longer delay time of about 30–50 crossing times, consistent with results from Reinoso et al. (2018, fig. 5). Figure 6. View largeDownload slide Time evolution of the maximum mass in the system. We show the results for the six different accretion models and standard set of parameters. Except for model 5, all models efficiently convert at least half of the initial gas mass into one single massive object. Figure 6. View largeDownload slide Time evolution of the maximum mass in the system. We show the results for the six different accretion models and standard set of parameters. Except for model 5, all models efficiently convert at least half of the initial gas mass into one single massive object. The results presented so far show that the formation of a central massive object is the natural outcome, in the regime of high-accretion rates and taking into account the large radii of accreting Pop. III protostars. It is favourable to have a position-dependent accretion rate, such that the formation of a massive object in the core is more likely, which fuels subsequent collisions with other stars. In model 5, such a formation channel is missing. This plus the fact that the stars will eventually shrink due to the time-dependent accretion rate causes the massive object to be an order of magnitude less massive . Model 5 comes closest to resembling the evolution of a gas-free, equal mass Plummer sphere in which collisions usually do not become important until core collapse. This implies that for systems in which accretion-induced collisions do play a significant role, the formation of a massive object is sped up. Next, we want to further investigate how the collisions occur in the cluster: do they preferably happen with a dominant central object, or are there collisions between different pairs of objects occurring throughout the cluster? In Fig. 7, we plot the number of collision products in the cluster as a function of the total number of collisions that have occurred. We observe that the position-independent accretion models (black curves) can produce factor 2–3 more collision products than the position-dependent models (grey curves). This is somewhat expected considering that the position-dependent accretion models rapidly produce a dominant mass in the core, which will have a much larger cross-section for collisions. We note that the mass ratio between the two most massive objects in the system ranges from an order of magnitude for model 5 up to two orders of magnitude for the other models. Figure 7. View largeDownload slide The number of collision products present in the cluster as a function of number of collisions that have occurred. This is a measure of whether collisions occur with predominantly one massive object or between many different pairs of stars throughout the cluster. We clearly observe a difference between the radially independent (thick curves) and radially dependent (thin curves) accretion models. As is somewhat expected, the radially dependent accretion models produce a massive object in the centre which is involved in most collisions. Figure 7. View largeDownload slide The number of collision products present in the cluster as a function of number of collisions that have occurred. This is a measure of whether collisions occur with predominantly one massive object or between many different pairs of stars throughout the cluster. We clearly observe a difference between the radially independent (thick curves) and radially dependent (thin curves) accretion models. As is somewhat expected, the radially dependent accretion models produce a massive object in the centre which is involved in most collisions. 5.2 Variation of input parameters We perform an ensemble of numerical experiments with varying input parameters to determine the robustness of the result in the previous section, i.e. that a single massive object of mass ∼ 104 − 5 M⊙ is formed due to a high collision rate. We vary the initial total gas mass, gas cloud radius, number of protostars, and average accretion rate per protostar. At the end of the simulation, we determine the maximum mass in the system and the total number of collisions. The results are presented in Fig. 8. Figure 8. View largeDownload slide Mass of the most massive object in the system at the end of the simulation, Mmax, and the total fraction of collisions, Nc/N0, as a function of the model input parameters: initial gas cloud mass (top left), gas cloud radius (top right), initial number of protostars (bottom left), and initial average accretion rate per protostar (bottom right). Figure 8. View largeDownload slide Mass of the most massive object in the system at the end of the simulation, Mmax, and the total fraction of collisions, Nc/N0, as a function of the model input parameters: initial gas cloud mass (top left), gas cloud radius (top right), initial number of protostars (bottom left), and initial average accretion rate per protostar (bottom right). Starting with the total gas mass, we observe that the maximum mass increases with higher gas mass. This makes sense intuitively as more gas is available to accrete from, and the resulting protostars will be more massive and larger in size. This increases the collisional cross-section, which is confirmed by the increasing number of collisions as a function of Mg. The trend seems to saturate at a gas mass of about 105.5 M⊙, and all the different models start to converge. In this regime, the protostars are able to reach such large radii that a runaway collision process is triggered, irrespective of the detailed accretion model. At the low mass end of the diagram, we observe that the time-dependent accretion models become inefficient, producing maximum masses ≤ 103 M⊙. With little mass to accrete from, the protostars in the time-dependent models will generally remain small, resulting in a low collision rate, which is consistent with gas-free star clusters in virial equilibrium. In the top right-hand panel of Fig. 8, we present the maximum mass and number of collisions with varying gas cloud radius. For the infinite gas reservoir models, we observe an increasing trend. For a larger gas cloud radius, the protostars need to accrete more gas before they reach the required radius for collisions to occur. As a result, the collision products will be more massive. For the finite gas reservoir models, there is a peak at Rg ∼ 10 − 1.5-10 − 1 pc after which the maximum mass decreases again. Intuitively, we can relate this decrease to the decrease in the number density, which is proportional to the collision rate. This is confirmed by the decrease of the number of collisions with Rg. For the densest configurations, i.e. the smallest gas cloud radii, we again observe a convergence of the different accretion models. With decreasing Rg, the ratio of the protostellar radii to Rg increases, resulting in an increased collision rate, which is confirmed by the data. The trend that Mmax increases with Rg, for the densest configurations, is explained by the fact that the protostars have to accrete for a longer time to reach a significant radius relative to Rg, thus resulting in more massive objects. The initial increase and the subsequent decrease of Mmax with increasing Rg are thus explained by the increase of the accretion time-scale and the decrease of the number density. In the bottom left-hand panel of Fig. 8, we vary the number of initial protostars in the cluster. Intuitively, we would expect that with a higher number of stars, the collision rate will increase. However, with fewer stars in the cluster, they will be able to accrete more and thus increase their collisional cross-section. For the infinite gas reservoir models, we again observe a steady increase of the maximum mass with N. The finite models produce massive objects between 104.5 − 5 M⊙ irrespective of the variation of N from 64 to 1024. The total collision fraction also remains roughly constant between 0.6 and 0.8. Therefore, the increasing number density and the decreasing cross-section with N tend to cancel out. Only for model 5, we observe a decrease of Mmax and Nc/N0 with N. The extra effect of the time-dependent accretion rate and subsequent shrinking of the protostellar radii causes the decreasing collisional cross-section to dominate. For the time-dependent accretion models, it is therefore favourable to have fewer stars in the system in order to produce a massive object. In the bottom right-hand panel of Fig. 8, we plot the maximum mass and total collision fraction as a function of average accretion rate per protostar. We observe an increasing collision fraction as a function of $$\dot{m}$$. This can be explained by the mass–radius tracks: a higher $$\dot{m}$$ produces larger protostellar radii. For the largest values of $$\dot{m}$$ however, the variation in the protostellar radii saturates (see Fig. 1), and we observe a flattening of the total collision fraction. The maximum mass depends weakly on $$\dot{m}$$, showing a relatively mild increase. For our standard set of parameter values defined in Section 3, we conclude that for $$\dot{m} \ge 10^{-3}\,\rm {M}_{\odot }\, \rm {yr^{-1}}$$ a massive object of M ≥ 104 M⊙ is formed. Finally, we measure the dependence of the formation of a massive object to a different underlying mass–radius parametrization. We compare our parametrization based on Hosokawa et al. (2012) to that of Haemmerlé et al. (2018). We use the standard set of parameters defined in Section 3 and create 10 initial realizations of the cluster. We integrate these initial conditions using both mass–radius parametrizations and calculate the average and standard deviation of the final maximum mass, number of collisions, and maximum collision rate (see Table 4). For models 2 and 4, we find the results to be mutually consistent. In these radial-dependent accretion models, a dominant central object is formed in the core. Its rapid growth causes many collisions to occur irrespective of the detailed mass–radius parametrization. For models 3, 5, and 6, we observe that the models based on Haemmerlé et al. (2018) produce somewhat fewer collisions and less massive objects. These models are sensitive to the more conservative radii of the Haemmerlé et al. (2018) parametrization because the accretion model is radially independent (so no dominant object in the core due to accretion) and/or because the accretion model is time dependent, in which protostars shrink to small sizes and differences in the radii will have a large effect on the dynamics-dominated collision rate. For model 1, we observe the counter-intuitive result that the parametrization by Haemmerlé et al. (2018) produces fewer collisions but a more massive object. This can be explained by the difference in time at which the stopping condition was fulfilled. The collisions rate is lower, and the collisions therefore more spread out in time. In this infinite gas reservoir model, this gives the protostars more time to accrete mass. 6 CONCLUSION We present the first numerical study of the formation of massive black hole seeds through the formation channel of accretion and collisions in a dense Pop. III protostellar cluster. We take into account accurate mass–radius parametrizations for accreting Pop. III protostars, which are based on detailed stellar evolution calculations performed by Hosokawa et al. (2012) and Haemmerlé et al. (2018). These studies have shown that for high accretion rates ($$\dot{m} \ge 10^{-3}\,\rm {M}_{\odot }\, \rm {yr^{-1}}$$), the protostellar radii can become of order 102 − 4 R⊙, thus significantly increasing the collisional cross-section. Using the amuse simulation framework, we perform a series of multiphysics simulations, including N-body dynamics, an analytical gas potential, six different accretion models, mass–radius parametrizations, and stellar collisions. In Tables 2 and 3, we present an overview of the simulation input parameters and outcome statistics. Our most conservative model (model 5), which most closely resembles the dynamics of ordinary Plummer spheres, confirms the fractional collision rate of 0.1–1 per cent per crossing time found in previous studies. Depending on the values of the input parameters, about 10 per cent of the initial protostellar population collides into a single massive object of order 104 M⊙ on a time-scale of 50 crossing times or longer. In our remaining accretion models, the effects of accretion-induced collisions are more prominent for two distinct but related reasons. First, in case of a more extensive gas reservoir the protostars will have more time to accrete gas and as a result be able to reach larger radii. Second, massive objects form in the centre of the cluster where gas volume densities, and thus the accretion rates and collision rates, are highest. Our results show that the maximum number of collisions per crossing time can be increased up to 1–10 per cent of the initial protostellar population. This increased collision rate leads to the formation of a single massive object with a mass of order M = 104 − 5 M⊙. The increase of the collision rate is fuelled by the high accretion rates of the Pop. III protostars and the resulting large radii. After the maximum collision rate is reached, the collision rate decreases again, either due to the sparsity of collision partners left in the system, or due to the lack of gas, which truncates the growth of the protostars and allows for dynamical ejections from the cluster. In Figs 4–6, we also demonstrate that the formation of a massive object can occur within 20 crossing times, which is about twice as fast compared to gas-free models (e.g. Reinoso et al. 2018). We have varied the initial gas cloud mass and gas cloud radius, and confirm the intuitive result that most collisions occur in the densest systems (see Fig. 8, top two panels). Also, the final mass of the seed black hole increases with higher gas cloud mass. Our calculations based on models with a limited gas reservoir show a characteristic gas cloud radius of Rg = 10 − 1.5-10 − 1 pc, at which the final mass of the most massive object is a maximum, i.e. M ∼ 104.5 − 5 M⊙. For smaller gas cloud radii, the system becomes more dense, which triggers a runaway collision between the protostars at earlier times. The less time there is for the protostars to accrete gas, the less massive will be the collision products. In the opposite regime of larger gas cloud radii, a runaway collision might not occur as the protostars cannot reach the required radii. These systems eventually become stellar-mass dominated, and since the collision rate is proportional to the number density, we would expect fewer collisions in larger systems. The final mass of the seed black hole depends mildly on the initial number of protostars (ranging from N = 64–1024) and average accretion rate (ranging from $$\dot{m} = 0.001 - 0.3\,\rm {M}_{\odot }\,\rm {yr^{-1}}$$), consistently producing massive objects of order M = 104 M⊙ or higher. The mass–radius evolution of accreting Pop. III protostars is very uncertain. In order to estimate the associated uncertainty and the robustness of the results mentioned above, we compare two mass–radius parametrizations: one based on Hosokawa et al. (2012) and another more conservative model based on Haemmerlé et al. (2018). Our calculations confirm that a parametrization with more conservative radii results in lower collision rates. However, the two different mass–radius parametrizations, each based on detailed stellar evolution calculations, produce no strong differences in the final mass of the seed black hole. We conclude that the mechanism of accretion-induced collisions in dense Pop. III protostellar systems is a viable mechanism for explaining the formation of the first massive black hole seeds. Therefore, this investigation warrants follow-up studies, which improve on the realism of the detailed implementation and coupling of the different physics ingredients. Acknowledgements TB acknowledges support from Fundação para a Ciência e a Tecnologia (grant SFRH/BPD/122325/2016), and support from Center for Research & Development in Mathematics and Applications (CIDMA) (strategic project UID/MAT/04106/2013), and from ENGAGE SKA, POCI-01-0145-FEDER-022217, funded by COMPETE 2020 and FCT, Portugal. DRGS and BR thank for funding through Fondecyt regular (project code 1161247). BR thanks Conicyt for financial support (CONICYT-PFCHA/MagísterNacional/2017-22171385). DRGS, MF, and AMS acknowledge funding through the ‘Concurso Proyectos Internacionales de Investigación, Convocatoria 2015’ (project code PII20150171) and the BASAL Centro de Astrofísica y Tecnologías Afines (CATA) PFB-06/2007. DRGS is further grateful for funding via ALMA-Conicyt (project code 31160001), Quimal (project code QUIMAL170001), and the Anillo program ‘Formation and Growth of Supermassive Black Holes’ (project code ACT172033). RSK acknowledges financial support from the Deutsche Forschungsgemeinschaft via SFB 881 ‘The Milky Way System’ (sub-projects B1, B2, and B8) and SPP 1573 ‘Physics of the Interstellar Medium’. LH and RSK were supported by the European Research Council under the European Community’s Seventh Framework Programme (FP7/2007–2013) via the ERC Advanced Grant ‘STARLIGHT: Formation of the First Stars’ (project number 339177). Footnotes 1 http://www.amusecode.org/ REFERENCES Agarwal B., Khochfar S., 2015, MNRAS , 446, 160 https://doi.org/10.1093/mnras/stu1973 CrossRef Search ADS   Agarwal B., Cullen F., Khochfar S., Klessen R. S., Glover S. C. O., Johnson J., 2017, MNRAS , 468, L82 https://doi.org/10.1093/mnrasl/slx028 CrossRef Search ADS   Agarwal B., Regan J., Klessen R. S., Downes T. P., Zackrisson E., 2017, MNRAS , 470, 4034 https://doi.org/10.1093/mnras/stx1528 CrossRef Search ADS   Bañados E. et al.  , 2018, Nature , 553, 473 CrossRef Search ADS PubMed  Baumgardt H., Klessen R. S., 2011, MNRAS , 413, 1810 https://doi.org/10.1111/j.1365-2966.2011.18258.x CrossRef Search ADS   Begelman M. C., Shlosman I., 2009, ApJL , 702, L5 https://doi.org/10.1088/0004-637X/702/1/L5 CrossRef Search ADS   Begelman M. C., Volonteri M., Rees M. J., 2006, MNRAS , 370, 289 https://doi.org/10.1111/j.1365-2966.2006.10467.x CrossRef Search ADS   Boekholt T. C. N., Stutz A. M., Fellhauer M., Schleicher D. R. G., Matus Carrillo D. R., 2017, MNRAS , 471, 3590 https://doi.org/10.1093/mnras/stx1821 CrossRef Search ADS   Bovino S., Grassi T., Schleicher D. R. G., Latif M. A., 2014, ApJL , 790, L35 https://doi.org/10.1088/2041-8205/790/2/L35 CrossRef Search ADS   Bovino S., Grassi T., Schleicher D. R. G., Banerjee R., 2016, ApJ , 832, 154 https://doi.org/10.3847/0004-637X/832/2/154 CrossRef Search ADS   Bromm V., Loeb A., 2003, ApJ , 596, 34 https://doi.org/10.1086/377529 CrossRef Search ADS   Clark P. C., Glover S. C. O., Klessen R. S., 2008, ApJ , 672, 757 https://doi.org/10.1086/524187 CrossRef Search ADS   Clark P. C., Glover S. C. O., Smith R. J., Greif T. H., Klessen R. S., Bromm V., 2011, Science , 331, 1040 https://doi.org/10.1126/science.1198027 CrossRef Search ADS PubMed  D’Amico G., Panci P., Lupi A., Bovino S., Silk J., 2018, MNRAS , 473, 328 https://doi.org/10.1093/mnras/stx2419 CrossRef Search ADS   Devecchi B., Volonteri M., 2009, ApJ , 694, 302 https://doi.org/10.1088/0004-637X/694/1/302 CrossRef Search ADS   Devecchi B., Volonteri M., Rossi E. M., Colpi M., Portegies Zwart S., 2012, MNRAS , 421, 1465 https://doi.org/10.1111/j.1365-2966.2012.20406.x CrossRef Search ADS   Dijkstra M., Ferrara A., Mesinger A., 2014, MNRAS , 442, 2036, https://doi.org/10.1093/mnras/stu1007 CrossRef Search ADS   Dopcke G., Glover S. C. O., Clark P. C., Klessen R. S., 2013, ApJ , 766, 103 https://doi.org/10.1088/0004-637X/766/2/103 CrossRef Search ADS   Eggenberger P., Meynet G., Maeder A., Hirschi R., Charbonnel C., Talon S., Ekström S., 2008, Astrophys. Space Sci. , 316, 43 https://doi.org/10.1007/s10509-007-9511-y CrossRef Search ADS   Flaugher B. et al.  , 2015, AJ , 150, 150 CrossRef Search ADS   Fujii M. S., Portegies Zwart S., 2013, MNRAS , 430, 1018 https://doi.org/10.1093/mnras/sts673 CrossRef Search ADS   Fujii M., Iwasawa M., Funato Y., Makino J., 2007, PASJ , 59, 1095 https://doi.org/10.1093/pasj/59.6.1095 CrossRef Search ADS   Gallerani S., Fan X., Maiolino R., Pacucci F., 2017, PASA , 34, e022 https://doi.org/10.1017/pasa.2017.14 CrossRef Search ADS   Greif T. H., Springel V., White S. D. M., Glover S. C. O., Clark P. C., Smith R. J., Klessen R. S., Bromm V., 2011, ApJ , 737, 75 https://doi.org/10.1088/0004-637X/737/2/75 CrossRef Search ADS   Haemmerlé L., Woods T. E., Klessen R. S., Heger A., Whalen D. J., 2018, MNRAS , 474, 2757 CrossRef Search ADS   Hirano S., Hosokawa T., Yoshida N., Omukai K., Yorke H. W., 2015, MNRAS , 448, 568 https://doi.org/10.1093/mnras/stv044 CrossRef Search ADS   Hosokawa T., Omukai K., 2009, ApJ , 691, 823 https://doi.org/10.1088/0004-637X/691/1/823 CrossRef Search ADS   Hosokawa T., Yorke H. W., Omukai K., 2010, ApJ , 721, 478 https://doi.org/10.1088/0004-637X/721/1/478 CrossRef Search ADS   Hosokawa T., Omukai K., Yorke H. W., 2012, ApJ , 756, 93 https://doi.org/10.1088/0004-637X/756/1/93 CrossRef Search ADS   Hosokawa T., Yorke H. W., Inayoshi K., Omukai K., Yoshida N., 2013, ApJ , 778, 178 https://doi.org/10.1088/0004-637X/778/2/178 CrossRef Search ADS   Hosokawa T., Hirano S., Kuiper R., Yorke H. W., Omukai K., Yoshida N., 2016, ApJ , 824, 119 https://doi.org/10.3847/0004-637X/824/2/119 CrossRef Search ADS   Hut P., Makino J., McMillan S., 1995, ApJL , 443, L93 https://doi.org/10.1086/187844 CrossRef Search ADS   Katz H., Sijacki D., Haehnelt M. G., 2015, MNRAS , 451, 2352 https://doi.org/10.1093/mnras/stv1048 CrossRef Search ADS   Koushiappas S. M., Bullock J. S., Dekel A., 2004, MNRAS , 354, 292 https://doi.org/10.1111/j.1365-2966.2004.08190.x CrossRef Search ADS   Latif M. A., Schleicher D. R. G., 2015a, MNRAS , 449, 77 https://doi.org/10.1093/mnras/stu2573 CrossRef Search ADS   Latif M. A., Schleicher D. R. G., 2015b, A&A , 578, A118 https://doi.org/10.1051/0004-6361/201525855 CrossRef Search ADS   Latif M. A., Schleicher D. R. G., Schmidt W., Niemeyer J., 2013a, MNRAS , 433, 1607 https://doi.org/10.1093/mnras/stt834 CrossRef Search ADS   Latif M. A., Schleicher D. R. G., Schmidt W., Niemeyer J. C., 2013b, MNRAS , 436, 2989 https://doi.org/10.1093/mnras/stt1786 CrossRef Search ADS   Latif M. A., Schleicher D. R. G., Bovino S., Grassi T., Spaans M., 2014, ApJ , 792, 78 https://doi.org/10.1088/0004-637X/792/1/78 CrossRef Search ADS   Latif M. A., Bovino S., Van Borm C., Grassi T., Schleicher D. R. G., Spaans M., 2014, MNRAS , 443, 1979, https://doi.org/10.1093/mnras/stu1230 CrossRef Search ADS   Latif M. A., Bovino S., Grassi T., Schleicher D. R. G., Spaans M., 2015, MNRAS , 446, 3163 https://doi.org/10.1093/mnras/stu2244 CrossRef Search ADS   Latif M. A., Omukai K., Habouzit M., Schleicher D. R. G., Volonteri M., 2016, ApJ , 823, 40 https://doi.org/10.3847/0004-637X/823/1/40 CrossRef Search ADS   Lodato G., Natarajan P., 2006, MNRAS , 371, 1813 https://doi.org/10.1111/j.1365-2966.2006.10801.x CrossRef Search ADS   Lupi A., Colpi M., Devecchi B., Galanti G., Volonteri M., 2014, MNRAS , 442, 3616 https://doi.org/10.1093/mnras/stu1120 CrossRef Search ADS   Martínez-Barbosa C. A., Brown A. G. A., Boekholt T., Portegies Zwart S., Antiche E., Antoja T., 2016, MNRAS , 457, 1062 https://doi.org/10.1093/mnras/stw006 CrossRef Search ADS   Matsuoka Y. et al.  , 2017, preprint (ArXiv e-prints:1704.05854) McCrady N., Gilbert A. M., Graham J. R., 2003, ApJ , 596, 240 https://doi.org/10.1086/377631 CrossRef Search ADS   McMillan S. L. W., Hut P., 1996, ApJ , 467, 348 https://doi.org/10.1086/177610 CrossRef Search ADS   Moeckel N., Clarke C. J., 2011, MNRAS , 410, 2799 https://doi.org/10.1111/j.1365-2966.2010.17659.x CrossRef Search ADS   Mortlock D. J. et al.  , 2011, Nature , 474, 616 https://doi.org/10.1038/nature10159 CrossRef Search ADS PubMed  Oh S., Kroupa P., 2012, MNRAS , 424, 65 https://doi.org/10.1111/j.1365-2966.2012.21152.x CrossRef Search ADS   Omukai K., Palla F., 2003, ApJ , 589, 677 https://doi.org/10.1086/374810 CrossRef Search ADS   Omukai K., Schneider R., Haiman Z., 2008, ApJ , 686, 801 https://doi.org/10.1086/591636 CrossRef Search ADS   Pelupessy F. I., van Elteren A., de Vries N., McMillan S. L. W., Drost N., Portegies Zwart S. F., 2013, A&A , 557, A84 https://doi.org/10.1051/0004-6361/201321252 CrossRef Search ADS   Plummer H. C., 1911, MNRAS , 71, 460 https://doi.org/10.1093/mnras/71.5.460 CrossRef Search ADS   Portegies Zwart S. F., Baumgardt H., Hut P., Makino J., McMillan S. L. W., 2004, Nature , 428, 724 https://doi.org/10.1038/nature02448 CrossRef Search ADS PubMed  Portegies Zwart S. F., McMillan S. L., van Elteren A., Pelupessy F. I., de Vries N., 2013, Comput. Phys. Commun. , 184, 456 https://doi.org/10.1016/j.cpc.2012.09.024 CrossRef Search ADS   Portegies Zwart S. et al.   2009, New Astron. , 14, 369 https://doi.org/10.1016/j.newast.2008.10.006 CrossRef Search ADS   Reed S. L. et al.  , 2017, MNRAS , 468, 4702 https://doi.org/10.1093/mnras/stx728 CrossRef Search ADS   Rees M. J., 1984, ARA&A , 22, 471 CrossRef Search ADS   Regan J. A., Haehnelt M. G., 2009, MNRAS , 396, 343 https://doi.org/10.1111/j.1365-2966.2009.14579.x CrossRef Search ADS   Regan J. A., Johansson P. H., Wise J. H., 2016, MNRAS , 459, 3377 https://doi.org/10.1093/mnras/stw899 CrossRef Search ADS   Reinoso B., Schleicher D., Fellhauer M., Klessen R., Boekholt T. C. N., 2018, preprint (arXiv:1801.05891) Sakurai Y., Yoshida N., Fujii M. S., Hirano S., 2017, MNRAS , 472, 1677 https://doi.org/10.1093/mnras/stx2044 CrossRef Search ADS   Schleicher D. R. G., Palla F., Ferrara A., Galli D., Latif M., 2013, A&A , 558, A59 CrossRef Search ADS   Schleicher D. R. G., Spaans M., Glover S. C. O., 2010, ApJL , 712, L69 https://doi.org/10.1088/2041-8205/712/1/L69 CrossRef Search ADS   Schneider R., Omukai K., Inoue A. K., Ferrara A., 2006, MNRAS , 369, 1437 https://doi.org/10.1111/j.1365-2966.2006.10391.x CrossRef Search ADS   Shapiro S. L., 2005, ApJ , 620, 59 https://doi.org/10.1086/427065 CrossRef Search ADS   Smith R. J., Hosokawa T., Omukai K., Glover S. C. O., Klessen R. S., 2012, MNRAS , 424, 457 https://doi.org/10.1111/j.1365-2966.2012.21211.x CrossRef Search ADS   Sugimura K., Omukai K., Inoue A. K., 2014, MNRAS , 445, 544 https://doi.org/10.1093/mnras/stu1778 CrossRef Search ADS   Sutherland W. et al.  , 2015, A&A , 575, A25 CrossRef Search ADS   Umeda H., Hosokawa T., Omukai K., Yoshida N., 2016, ApJL , 830, L34 https://doi.org/10.3847/2041-8205/830/2/L34 CrossRef Search ADS   Wise J. H., Turk M. J., Abel T., 2008, ApJ , 682, 745 https://doi.org/10.1086/588209 CrossRef Search ADS   Woods T. E., Heger A., Whalen D. J., Haemmerlé L., Klessen R. S., 2017, ApJL , 842, L6 https://doi.org/10.3847/2041-8213/aa7412 CrossRef Search ADS   Wright E. L. et al.  , 2010, AJ , 140, 1868 CrossRef Search ADS   APPENDIX A: ANALYTICAL FORM OF THE MASS–RADIUS RELATION The analytical power-law parametrization of the mass–radius relation from the models of Hosokawa et al. (2012) and Haemmerlé et al. (2018) is given in Tables A1 and A2, respectively. Table A1. Simplified parametrization of the mass–radius relation of accreting Pop. III protostars based on Hosokawa et al. (2012, fig. 5). $$\dot{\rm {m}}$$  m  R  $$\rm {M}_{\odot }\,\rm {yr}^{-1}$$  M⊙  R⊙  ≥ 1  < 150  $$\rm {250} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.4}}$$    < 450  $$\rm {1855} \left( {\rm {m} \over \rm {150}\,\rm {M}_{\odot }} \right)^{\rm {0.8}}$$    ≥ 450  $$\rm {4468} \left( {\rm {m} \over \rm {450}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\rm {0.3} \le \dot{m} < \rm {1}$$  < 80  $$\rm {190} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.4}}$$    < 140  $$\rm {1096} \left( {\rm {m} \over \rm {80}\,\rm {M}_{\odot }} \right)^{\rm {1.5}}$$    ≥ 140  $$\rm {2538} \left( {\rm {m} \over \rm {140}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\rm {0.1} \le \dot{m} < \rm {0.3}$$  < 40  $$\rm {140} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.4}}$$    < 90  $$\rm {612} \left( {\rm {m} \over \rm {40}\,\rm {M}_{\odot }} \right)^{\rm {1.5}}$$    ≥ 90  $$\rm {2066} \left( {\rm {m} \over \rm {90}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\rm {0.06} \le \dot{m} < \rm {0.1}$$  < 25  $$\rm {110} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.4}}$$    < 70  $$\rm {399} \left( {\rm {m} \over \rm {25}\,\rm {M}_{\odot }} \right)^{\rm {1.5}}$$    ≥ 70  $$\rm {1868} \left( {\rm {m} \over \rm {70}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\rm {0.03} \le \dot{m} < \rm {0.06}$$  < 20  $$\rm {90} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.4}}$$    < 70  $$\rm {298} \left( {\rm {m} \over \rm {20}\,\rm {M}_{\odot }} \right)^{\rm {1.5}}$$    ≥ 70  $$\rm {1953} \left( {\rm {m} \over \rm {70}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\rm {0.006} \le \dot{m} < \rm {0.03}$$  < 20  $$\rm {50} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.4}}$$    < 35  $$\rm {166} \left( {\rm {m} \over \rm {20}\,\rm {M}_{\odot }} \right)^{\rm {-1.5}}$$    ≥ 35  $$\rm {72} \left( {\rm {m} \over \rm {35}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\dot{m} < \rm {0.006}$$  < 9  $$\rm {25} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.4}}$$    < 50  $$\rm {60} \left( {\rm {m} \over \rm {9}\,\rm {M}_{\odot }} \right)^{\rm {-1.5}}$$    ≥ 50  $$\rm {5} \left( {\rm {m} \over \rm {50}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\dot{\rm {m}}$$  m  R  $$\rm {M}_{\odot }\,\rm {yr}^{-1}$$  M⊙  R⊙  ≥ 1  < 150  $$\rm {250} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.4}}$$    < 450  $$\rm {1855} \left( {\rm {m} \over \rm {150}\,\rm {M}_{\odot }} \right)^{\rm {0.8}}$$    ≥ 450  $$\rm {4468} \left( {\rm {m} \over \rm {450}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\rm {0.3} \le \dot{m} < \rm {1}$$  < 80  $$\rm {190} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.4}}$$    < 140  $$\rm {1096} \left( {\rm {m} \over \rm {80}\,\rm {M}_{\odot }} \right)^{\rm {1.5}}$$    ≥ 140  $$\rm {2538} \left( {\rm {m} \over \rm {140}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\rm {0.1} \le \dot{m} < \rm {0.3}$$  < 40  $$\rm {140} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.4}}$$    < 90  $$\rm {612} \left( {\rm {m} \over \rm {40}\,\rm {M}_{\odot }} \right)^{\rm {1.5}}$$    ≥ 90  $$\rm {2066} \left( {\rm {m} \over \rm {90}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\rm {0.06} \le \dot{m} < \rm {0.1}$$  < 25  $$\rm {110} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.4}}$$    < 70  $$\rm {399} \left( {\rm {m} \over \rm {25}\,\rm {M}_{\odot }} \right)^{\rm {1.5}}$$    ≥ 70  $$\rm {1868} \left( {\rm {m} \over \rm {70}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\rm {0.03} \le \dot{m} < \rm {0.06}$$  < 20  $$\rm {90} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.4}}$$    < 70  $$\rm {298} \left( {\rm {m} \over \rm {20}\,\rm {M}_{\odot }} \right)^{\rm {1.5}}$$    ≥ 70  $$\rm {1953} \left( {\rm {m} \over \rm {70}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\rm {0.006} \le \dot{m} < \rm {0.03}$$  < 20  $$\rm {50} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.4}}$$    < 35  $$\rm {166} \left( {\rm {m} \over \rm {20}\,\rm {M}_{\odot }} \right)^{\rm {-1.5}}$$    ≥ 35  $$\rm {72} \left( {\rm {m} \over \rm {35}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\dot{m} < \rm {0.006}$$  < 9  $$\rm {25} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.4}}$$    < 50  $$\rm {60} \left( {\rm {m} \over \rm {9}\,\rm {M}_{\odot }} \right)^{\rm {-1.5}}$$    ≥ 50  $$\rm {5} \left( {\rm {m} \over \rm {50}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  View Large Table A2. Simplified parametrization of the mass–radius relation of accreting Pop. III protostars based on Haemmerlé et al. (2018, fig. 2). $$\dot{\rm {m}}$$  m  R  $$\rm {M}_{\odot }\,\rm {yr}^{-1}$$  M⊙  R⊙  ≥ 10  < 10  $$\rm {200} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.0}}$$    < 50  $$\rm {200} \left( {\rm {m} \over \rm {10}\,\rm {M}_{\odot }} \right)^{\rm {-0.4}}$$    < 60  $$\rm {105} \left( {\rm {m} \over \rm {50}\,\rm {M}_{\odot }} \right)^{\rm {13}}$$    ≥ 60  $$\rm {1124} \left( {\rm {m} \over \rm {60}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\rm {1} \le \dot{m} < \rm {10}$$  < 10  $$\rm {200} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.0}}$$    < 35  $$\rm {200} \left( {\rm {m} \over \rm {10}\,\rm {M}_{\odot }} \right)^{\rm {-0.6}}$$    < 42  $$\rm {94} \left( {\rm {m} \over \rm {35}\,\rm {M}_{\odot }} \right)^{\rm {13}}$$    ≥ 42  $$\rm {1009} \left( {\rm {m} \over \rm {42}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\rm {0.1} \le \dot{m} < \rm {1}$$  < 10  $$\rm {200} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.0}}$$    < 25  $$\rm {200} \left( {\rm {m} \over \rm {10}\,\rm {M}_{\odot }} \right)^{\rm {-0.6}}$$    < 29.5  $$\rm {115} \left( {\rm {m} \over \rm {25}\,\rm {M}_{\odot }} \right)^{\rm {13}}$$    ≥ 29.5  $$\rm {993} \left( {\rm {m} \over \rm {29.5}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\rm {0.01} \le \dot{m} < \rm {0.1}$$  < 30  $$\rm {200} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.0}}$$    < 70  $$\rm {200} \left( {\rm {m} \over \rm {30}\,\rm {M}_{\odot }} \right)^{\rm {-2.5}}$$    ≥ 70  $$\rm {24} \left( {\rm {m} \over \rm {70}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\dot{m} < \rm {0.01}$$  < 10  $$\rm {200} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.0}}$$    < 50  $$\rm {200} \left( {\rm {m} \over \rm {10}\,\rm {M}_{\odot }} \right)^{\rm {-2.5}}$$    ≥ 50  $$\rm {3.6} \left( {\rm {m} \over \rm {50}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\dot{\rm {m}}$$  m  R  $$\rm {M}_{\odot }\,\rm {yr}^{-1}$$  M⊙  R⊙  ≥ 10  < 10  $$\rm {200} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.0}}$$    < 50  $$\rm {200} \left( {\rm {m} \over \rm {10}\,\rm {M}_{\odot }} \right)^{\rm {-0.4}}$$    < 60  $$\rm {105} \left( {\rm {m} \over \rm {50}\,\rm {M}_{\odot }} \right)^{\rm {13}}$$    ≥ 60  $$\rm {1124} \left( {\rm {m} \over \rm {60}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\rm {1} \le \dot{m} < \rm {10}$$  < 10  $$\rm {200} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.0}}$$    < 35  $$\rm {200} \left( {\rm {m} \over \rm {10}\,\rm {M}_{\odot }} \right)^{\rm {-0.6}}$$    < 42  $$\rm {94} \left( {\rm {m} \over \rm {35}\,\rm {M}_{\odot }} \right)^{\rm {13}}$$    ≥ 42  $$\rm {1009} \left( {\rm {m} \over \rm {42}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\rm {0.1} \le \dot{m} < \rm {1}$$  < 10  $$\rm {200} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.0}}$$    < 25  $$\rm {200} \left( {\rm {m} \over \rm {10}\,\rm {M}_{\odot }} \right)^{\rm {-0.6}}$$    < 29.5  $$\rm {115} \left( {\rm {m} \over \rm {25}\,\rm {M}_{\odot }} \right)^{\rm {13}}$$    ≥ 29.5  $$\rm {993} \left( {\rm {m} \over \rm {29.5}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\rm {0.01} \le \dot{m} < \rm {0.1}$$  < 30  $$\rm {200} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.0}}$$    < 70  $$\rm {200} \left( {\rm {m} \over \rm {30}\,\rm {M}_{\odot }} \right)^{\rm {-2.5}}$$    ≥ 70  $$\rm {24} \left( {\rm {m} \over \rm {70}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\dot{m} < \rm {0.01}$$  < 10  $$\rm {200} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.0}}$$    < 50  $$\rm {200} \left( {\rm {m} \over \rm {10}\,\rm {M}_{\odot }} \right)^{\rm {-2.5}}$$    ≥ 50  $$\rm {3.6} \left( {\rm {m} \over \rm {50}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  View Large APPENDIX B: VALIDATION OF THE NUMERICAL METHOD We present here additional plots for the validation experiments discussed in Section 4. In Fig. B1, we confirm that in the absence of accretion, the 10, 50, and 90 per cent Lagrangian radii follow those of a Plummer sphere, with a slight discrepancy in the 90 per cent radius due to the truncation that we have introduced. In Fig. B2, we present the time evolution of the total star and gas mass (top row) and the time evolution of the average accretion rate (bottom row) in a set-up where the N-body integrator has been switched off. The figure follows the expected behaviour of the accretion models. Finally, we show in Fig. B3 that the specific choice of the initial protostellar mass does not have a large influence on the time evolution of the most massive object in the system. The initial total stellar mass is given by Ms = Nm0, which becomes comparable to the initial gas cloud mass, Mg = 105 M⊙, and N = 256, if m0 = 390.625 M⊙. When m0 ≪ 390M⊙, the total mass (gas + stars) remains close to ∼ 105 M⊙ and, with a fixed collisional fraction of about 10 percent, produces a maximum mass of about 104M⊙. When m0 = 100 M⊙, the initial stellar mass Ms = 2.56 × 104M⊙, and so we have effectively increased the initial total mass by a quarter, which also allows for the formation of a more massive object. However, such a system would not start out as a gas-dominated system. Figure B1. View largeDownload slide Validation of the initial conditions and dynamics. We present the 10, 50, and 90 per cent Lagrangian radius of a Plummer sphere (see the text for the Plummer parameters), calculated analytically (dashed line, for an untruncated Plummer) and from our numerical model (solid line, truncated at 5 Plummer radii). We confirm that our cluster is initially stable and will only change its structure in the presence of accretion and collisions. Figure B1. View largeDownload slide Validation of the initial conditions and dynamics. We present the 10, 50, and 90 per cent Lagrangian radius of a Plummer sphere (see the text for the Plummer parameters), calculated analytically (dashed line, for an untruncated Plummer) and from our numerical model (solid line, truncated at 5 Plummer radii). We confirm that our cluster is initially stable and will only change its structure in the presence of accretion and collisions. Figure B2. View largeDownload slide Illustration and validation of the gas accretion models. In the top row, we show the total mass evolution of the stars and gas. We confirm that models 1 and 2 are consistent (top left-hand panel), and similar for models 3 and 4 (top middle panel) and models 5 and 6 (top right-hand panel). In the bottom row, we present the average accretion rates per star, i.e. the time derivative of the total stellar mass in the top row divided by the total number of stars. Figure B2. View largeDownload slide Illustration and validation of the gas accretion models. In the top row, we show the total mass evolution of the stars and gas. We confirm that models 1 and 2 are consistent (top left-hand panel), and similar for models 3 and 4 (top middle panel) and models 5 and 6 (top right-hand panel). In the bottom row, we present the average accretion rates per star, i.e. the time derivative of the total stellar mass in the top row divided by the total number of stars. Figure B3. View largeDownload slide Time evolution of the maximum mass in the system, for our most conservative accretion model (Model 5) with standard parameters (see Section 3). We vary the initial masses of the protostars and find that as long as the system starts out as a gas-dominated system, that is with m0 ≪ 390 M⊙, the final maximum mass does not change much. For m0 = 100 M⊙ the overall potential starts to become strongly influenced by the stars, and mass growth gets modified by Bondi–Hoyle–Littleton accretion. Figure B3. View largeDownload slide Time evolution of the maximum mass in the system, for our most conservative accretion model (Model 5) with standard parameters (see Section 3). We vary the initial masses of the protostars and find that as long as the system starts out as a gas-dominated system, that is with m0 ≪ 390 M⊙, the final maximum mass does not change much. For m0 = 100 M⊙ the overall potential starts to become strongly influenced by the stars, and mass growth gets modified by Bondi–Hoyle–Littleton accretion. © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Monthly Notices of the Royal Astronomical Society Oxford University Press

# Formation of massive seed black holes via collisions and accretion

, Volume 476 (1) – May 1, 2018
15 pages

/lp/ou_press/formation-of-massive-seed-black-holes-via-collisions-and-accretion-GOnsweZcCB
Publisher
The Royal Astronomical Society
ISSN
0035-8711
eISSN
1365-2966
D.O.I.
10.1093/mnras/sty208
Publisher site
See Article on Publisher Site

### Abstract

Abstract Models aiming to explain the formation of massive black hole seeds, and in particular the direct collapse scenario, face substantial difficulties. These are rooted in rather ad hoc and fine-tuned initial conditions, such as the simultaneous requirements of extremely low metallicities and strong radiation backgrounds. Here, we explore a modification of such scenarios where a massive primordial star cluster is initially produced. Subsequent stellar collisions give rise to the formation of massive (104−105 M⊙) objects. Our calculations demonstrate that the interplay among stellar dynamics, gas accretion, and protostellar evolution is particularly relevant. Gas accretion on to the protostars enhances their radii, resulting in an enhanced collisional cross-section. We show that the fraction of collisions can increase from 0.1 to 1 per cent of the initial population to about 10 per cent when compared to gas-free models or models of protostellar clusters in the local Universe. We conclude that very massive objects can form in spite of initial fragmentation, making the first massive protostellar clusters viable candidate birth places for observed supermassive black holes. methods: numerical, stars: black holes, stars: kinematics and dynamics, stars: Population III 1 INTRODUCTION More than 100 supermassive black holes have already been detected at z > 6 (Gallerani et al. 2017), with the number still continuously increasing, including 32 quasars recently discovered between z = 5.7 and 6.8 via Subaru (Matsuoka et al. 2017) and 8 z > 6 quasars via Spectral Energy Distribution (SED) model fitting of VISTA (Sutherland et al. (2015), WISE (Wright et al. 2010), and Dark Energy Survey (Flaugher et al. 2015) Year 1 observations (Reed et al. 2017). The currently most distant known quasars are at z = 7.085, i.e. 0.77 billion years after the big bang, with a mass of about 2 × 109 M⊙ (Mortlock et al. 2011), and at z = 7.54 with a mass of about 8 × 108 M⊙ (Bañados et al. 2018). Explaining the existence of these objects provides a significant challenge to our cosmological model (D’Amico et al. 2018), as the accretion at an Eddington rate requires initial seed masses of the order 104 M⊙, when realistic spin parameters and accretion disc models are taken into account (Shapiro 2005). The only solution is very massive seeds or extended periods of super-Eddington accretion or potentially also combinations of both during the formation and early growth of massive black holes. The possible scenarios leading to their formation were already outlined by Rees (1984), including the direct collapse of massive gas clouds either to a black hole or to a supermassive star, which later collapses to a black hole via general relativistic instabilities, or alternatively through the formation of a dense stellar cluster, which may either collapse into a black hole through relativistic instabilities or evolve due to stellar mergers. From these, the direct collapse model was often considered as the most promising scenario, as it can potentially produce the most massive seeds. Their formation has been proposed through low-angular momentum material in the first pre-galactic haloes (Koushiappas, Bullock & Dekel 2004), through efficient angular momentum transport by gravitational instabilities (Begelman, Volonteri & Rees 2006; Lodato & Natarajan 2006) or through bars-within-bars instabilities (Begelman & Shlosman 2009). Numerical simulations probing the formation of massive black holes via direct collapse were indicating that the latter are possible only if cooling is efficiently suppressed, for instance if the formation of molecular hydrogen is photodissociated via a strong radiation background (Bromm & Loeb 2003). Cosmological hydrodynamics simulations following such a collapse over many orders of magnitude indeed found signatures of self-gravitational instabilities in the regime of atomic hydrogen cooling (Wise, Turk & Abel 2008). An update of the pathways leading to massive black hole formation was given by Regan & Haehnelt (2009), considering the cosmological conditions that may help to keep the gas metal-free, thereby preventing strong fragmentation events. The first systematic fragmentation study in the atomic cooling regime has been pursued by Latif et al. (2013a), showing that fragmentation indeed can be strongly suppressed if cooling is only feasible via atomic hydrogen lines (see also Schleicher, Spaans & Glover 2010, for the role of the atomic line cooling transitions). The large accretion rates of 0.1 M⊙ yr−1 or more found in the simulations strongly favour protostellar models with cool atmospheres and highly extended envelopes (Hosokawa et al. 2013; Schleicher et al. 2013; Haemmerlé et al. 2018; Woods et al. 2017), implying that feedback is rather inefficient. Under such conditions, final black hole masses of 105 M⊙ can be reached based on the results of cosmological simulations (Latif et al. 2013b; Umeda et al. 2016). Such a scenario however represents the most optimistic case. In particular, one needs to investigate the strength of the radiation background to keep the gas atomic, which is typically described through the strength of the radiation background at 13.6 eV and parametrized via J21, with a value of 1 corresponding to 10−21 erg s−1 cm−3 Hz−1 sr−1. The first numerical investigations suggested a critical value of J21 ∼ 100 to prevent the formation of molecular hydrogen, while updated chemical networks and more realistic models for the radiation background (see e.g. Sugimura, Omukai & Inoue 2014; Agarwal & Khochfar 2015) have led to much larger critical values of the order 105 when applied in cosmological simulations (Latif et al. 2014, 2015, see also Agarwal et al. 2017 for a discussion of the impact of the spectral shape). Under the threshold, massive stars may potentially still form, even though the mass may be reduced by factors of 10–100 (Latif et al. 2014). In addition to molecular hydrogen line cooling, fragmentation can also be induced via metals or dust grains (Omukai et al. 2008). In the case of metal line cooling, a metallicity of 10−3 Z⊙ can already increase the cooling and trigger fragmentation within cosmological simulations (Bovino et al. 2014), while even lower metallicities of 10−5 Z⊙ are sufficient when dust cooling is considered (Schneider et al. 2006; Dopcke et al. 2013; Bovino et al. 2016; Latif et al. 2016). The need to have very strong radiation backgrounds, while keeping the gas metal-free, leads to a strong need of fine-tuning, which at best can be satisfied under very rare conditions (e.g. Agarwal et al. 2017), while in fact the need for large values of J21 provides a problem for the direct collapse scenario (Dijkstra, Ferrara & Mesinger 2014). Surprisingly, the alternative pathway of black hole formation through stellar clusters has been investigated to a lesser degree in the context of early Universe black hole formation. Analytical models by Devecchi & Volonteri (2009), Devecchi et al. (2012), and Lupi et al. (2014) predict black hole masses of 100–1000 M⊙ forming in the first stellar clusters. Katz, Sijacki & Haehnelt (2015) model the formation of a dense stellar cluster in a cosmological simulation, and show the subsequent formation of ∼1000 M⊙ black hole via N-body simulations. Similarly, Sakurai et al. (2017) showed the formation of black holes with 400–1900 M⊙ via a combination of cosmological and N-body simulations. So far, these simulations have not considered the enhanced protostellar radii in the presence of accretion (Hosokawa et al. 2013; Schleicher et al. 2013), and they also neglected the initial presence of gas in the cluster after the formation of the protostars. As we will show in this paper, their ongoing accretion events may however considerably favour collisions and mergers, and thus the formation of a central massive object. For present-day protostellar clusters, Baumgardt & Klessen (2011) have shown that 0.1–1 per cent of the protostars in the cluster may collide and help to form a particularly massive star within the cluster (see also Moeckel & Clarke 2011; Oh & Kroupa 2012; Fujii & Portegies Zwart 2013, for similar results). For star clusters at high redshift, we can however expect significantly higher collision rates. One reason is that these clusters are more dense, as the trace amount of metals will trigger cooling and fragmentation only when high densities of order 109 cm−3 are reached, thereby leading to the formation of very compact star clusters (Clark, Glover & Klessen 2008; Clark et al. 2011; Greif et al. 2011). Due to the low metallicity and the larger gas temperatures, the accretion rates will also be enhanced, thus favouring collisions through the accretion and mass growth itself, but also due to the protostellar radii that are enhanced as a result (Smith et al. 2012). We present here the first investigation which explores the formation of massive black hole seeds from a dense stellar cluster, where gas-phase effects like accretion as well as the resulting enhanced protostellar radii are taken into account. We describe our numerical methods in Section 2 and our experimental set-up in Section 3. The validation is described in Section 4. Our results are presented in Section 5, including both our reference model and an exploration of the parameter dependence. We summarize the main conclusions in Section 6. 2 NUMERICAL METHODS Modelling the early evolution of a Population III (Pop. III) protostar cluster is challenging. The main reason is the variety of physical processes that play a role, i.e. gravitational N-body dynamics, gravitational coupling between the stars and the gas, stellar growth in mass and size due to gas accretion, and stellar collisions. In this section, we describe each of these physical ingredients that go into our simulations, and how we couple them into one numerical model. We use the Astrophysical Multi-purpose Software Environment (AMUSE;1 Portegies Zwart, McMillan & Harfst 2009; Pelupessy et al. 2013; Portegies Zwart et al. 2013), which was designed for performing such multiphysics simulations as required for our study. Particularly, it has the flexibility to introduce new physical ingredients, such as a mass–radius parametrization for accreting Pop. III protostars and to couple it to existing N-body codes and background potentials (e.g. Martínez-Barbosa et al. 2016; Boekholt et al. 2017). 2.1 Initial conditions and dynamics Our astrophysical system under investigation consists of Pop. III protostars embedded in their natal gas cloud. We will assume a simplified initial condition: the protostars and the gas are distributed equally, and they both follow the commonly used Plummer distribution (Plummer 1911). In order to have a well-defined size of the cluster and to ensure that each protostar starts out within the gas cloud, we introduce a cut-off radius after which the density is set to zero. We set this radius to five times the Plummer radius so that the cluster remains stable. The parameters specifying the initial conditions are then the total gas mass, Mg; the cut-off radius, Rg; and the number of protostars, N. The initial mass of the protostars is set to m0 = 0.1 M⊙. In Appendix B, we confirm that the specific choice of initial mass does not change the results much, as long as we start out with a gas-dominated system. Complicating factors such as a flattened distribution, cluster rotation, and an initial binary fraction are not taken into account in the this study. Our main aim is to explore the complex interplay between dynamics and accretion that will lead to the collisional growth of a massive object. To model the star–star gravitational interactions, we use the N-body code ph4 (e.g. McMillan & Hut 1996, , section 3.2), which is a fourth-order Hermite algorithm in combination with the time-symmetric integration scheme of Hut, Makino & McMillan (1995). We simplify the gravitational dynamics of the gas cloud by an analytical background potential. This potential is coupled to the stars using the BRIDGE method (Fujii et al. 2007), so that the stars experience both the gravitational force from each other as well as from the gas. Especially at the start of the simulation, when the protostars are still low mass, their orbits will be completely determined by the dominant background potential. 2.2 Gas accretion models The initially low mass Pop. III protostars will gain mass by accreting from the gas reservoir. Since the accretion rate might vary with cluster environment and cluster evolution, we define six different accretion models based on: (in)finite gas reservoir, position-(in)dependent accretion rate, and time-(in)dependent accretion rate. An infinite gas reservoir resembles a system that is constantly being fed fresh gas. This prolongs the gas-dominated phase and the time-scale on which the protostars can accrete gas. This is contrary to the finite gas reservoir models, where the gas will eventually run out, and the protostars will stop accreting. For the position-dependent models, we set the accretion rate proportional to the local gas density. In this way, the protostars in the core accrete at a higher rate than protostars in the halo. Such a system naturally produces a range of stellar masses and radii. In order to compare the position-(in)dependent models, we make sure that the cumulative accretion rate is initially equal. A time dependence to the accretion rate is introduced if we set the accretion rate proportional to the gas density, which in turn decreases in time for the models with a finite gas reservoir. The main effect of a decreasing accretion rate is that the stars will migrate to lower mass–radius tracks, and thus have a decreasing collisional cross-section. The different combinations of the accretion properties described above define the six different accretion models presented in Table 1. An illustration of the time evolution of the accretion rate in the different models is presented in Appendix B. Table 1. Six different gas accretion models for Pop. III protostars embedded in their natal gas cloud. Model  Gas reservoir  Position-dependent  Time-dependent      accretion model  accretion model  1  Infinite  No  No  2  Infinite  Yes  No  3  Finite  No  No  4  Finite  Yes  No  5  Finite  No  Yes  6  Finite  Yes  Yes  Model  Gas reservoir  Position-dependent  Time-dependent      accretion model  accretion model  1  Infinite  No  No  2  Infinite  Yes  No  3  Finite  No  No  4  Finite  Yes  No  5  Finite  No  Yes  6  Finite  Yes  Yes  View Large 2.3 Mass–radius evolution The mass–radius evolution of accreting Pop. III protostars is very uncertain, and so we investigate two different sets of models to have a handle on the uncertainties introduced by different protostellar models. We base our mass–radius parametrizations on two different studies: Hosokawa et al. (2012, fig. 5) and Haemmerlé et al. (2018, fig. 2). Both of these studies performed detailed stellar evolution calculations using codes by Omukai & Palla (2003), Hosokawa & Omukai (2009), Hosokawa, Yorke & Omukai (2010), and GENEVA (Eggenberger et al. 2008), respectively. The radius of a protostar is completely determined by its mass, m, and accretion rate, $$\dot{m}$$. At every time-step in our simulation, we keep track of these two quantities and update the radius of the protostar. For values of the accretion rate in between the parametrized mass–radius tracks, we use interpolation between the two nearest tracks in log-space. In Fig. 1, we present our approximate parametrization of the mass–radius tracks of accreting Pop. III protostars, based on the detailed calculations of Hosokawa et al. (2012, fig. 5) and Haemmerlé et al. (2018, fig. 2). For the higher accretion rates, the models of Haemmerlé et al. (2018) produce somewhat smaller radii, but they show the same behaviour at large masses. We use both models to test the sensitivity of the formation of massive objects to the underlying mass–radius parametrization. The analytical form of the parametrization is given in Appendix A. Figure 1. View largeDownload slide Two parametrizations of the mass–radius evolution of accreting Pop. III protostars based on Hosokawa, Omukai & Yorke (2012, fig. 5) (left-hand panel) and Haemmerlé et al. (2018, fig. 2) (right-hand panel). Figure 1. View largeDownload slide Two parametrizations of the mass–radius evolution of accreting Pop. III protostars based on Hosokawa, Omukai & Yorke (2012, fig. 5) (left-hand panel) and Haemmerlé et al. (2018, fig. 2) (right-hand panel). 2.4 Stellar collisions We adopt the commonly used ‘sticky-sphere’ approximation to treat collisions between protostars. Whenever the distance between two protostars is less than the sum of their radii, we replace the two protostars by a single object at their centre of mass. We assume that during the collision the total mass is conserved, and the new radius is determined by the mass–radius parametrization. In case the accretion rate is very low, i.e. < 10 − 6 M⊙ yr − 1, then the new radius is determined by conserving the density of the primary star. 3 EXPERIMENTAL SET-UP The input parameters specifying a simulation and their range in values are as follows: Gas cloud mass Mg = 104, 3 × 104, 105, 3 × 105, 106 M⊙ Gas cloud radius Rg = 0.01, 0.03, 0.1, 0.3, 1.0 pc Number of protostars N = 64, 128, 256, 512, 1024 Average accretion rate $$\dot{m}$$ = 0.001, 0.003, 0.01, 0.03, 0.1, 0.3 $$\rm {M}_{\odot }\,\rm {yr}^{-1}$$. We also vary the accretion model (see Section 2.2) and the mass–radius parametrization (see Section 2.3). We note that the adopted protostellar accretion rates are high compared to present-day star formation. However, they are realistic for the primordial case. Similar to current star formation, the protostellar accretion process in Pop. III clusters is not regulated via Bondi–Hoyle–Littleton accretion, but rather it has been shown that gravitational instabilities in the gas are driving fragmentation as well as the accretion process on to the fragments. This leads to the typical range of accretion rates under Pop. III conditions that we adopted, see, e.g. Clark et al. (2011), Greif et al. (2011), Smith et al. (2012), Latif et al. (2013b), Latif & Schleicher (2015a), Latif & Schleicher (2015b), Hirano et al. (2015), Regan, Johansson & Wise (2016), Hosokawa et al. (2016), and Latif et al. (2016). We define a standard set of parameters as: Mg = 105 M⊙, Rg = 0.1 pc, N = 256, $$\dot{m} = 0.03\,\rm {M_{\odot }\,yr^{-1}}$$, and mass–radius parametrization based on Hosokawa et al. (2012). This choice of parameter values reflects that we are particularly interested in very massive Pop. III protostar clusters and the formation of very massive objects. The initial crossing time of our systems is given by   $$T_{\rm {cross,0}} = 853\,\left( \frac{10^5\,\rm {M}_{\odot }}{M_{\rm {g}}} \right)^{\frac{1}{2}} \left(\frac{R_{\rm {g}}}{0.1\,\rm { pc}}\right)^{\frac{3}{2}}\,\rm {yr},$$ (1)where we used Tcross,0 = 2Rv/σ, with Rv the virial radius given by Rv = 16Rg/15π (using the definition of our gas cloud truncation radius), σ the velocity dispersion estimated by $$\sigma = \sqrt{G M_{\rm {g}} / 2 R_{\rm {v}}}$$, and G the gravitational constant. We consider a simulation finished if most of the collisions have occurred. We determine this by keeping track of the average collision rate as   $$R_{\rm {av}}\left( t \right) = \frac{N_{\rm {col}}\left( t \right)}{t_{\rm {last\,collision}}},$$ (2)and an upper limit of the current collision rate as   $$R\left( t \right) = \frac{1}{t-t_{\rm {last\,collision}}}.$$ (3)If the ratio of R/Rav < 0.015, we stop the simulation. This criterion was chosen to make sure that the majority of collisions have occurred, while also limiting the duration of the simulation. We also stop the simulation if no collision has occurred in the last million years. For each simulation, we store regular snapshots of the full phase space information, allowing us to retrace the mass–radius evolution and calculate collision rates. In Tables 2 and 3, we provide an overview of our simulations and their input parameters, together with several statistics describing the outcome of the simulations, such as the maximum mass and maximum collision rate. An estimate of the measurement uncertainties can be found in Table 4. Table 2. Overview of the simulations. The input parameters are: gas cloud mass Mg, gas cloud radius Rg, number of protostars N, and the average accretion rate $$\dot{m}$$. Statistics describing the output are: final mass of the most massive object Mm and total number of collisions Nc. The number in the subscript denotes the accretion model (see Table 1). Mg  Rg  N  $$\dot{m}$$  Mm,1  Mm,2  Mm,3  Mm,4  Mm,5  Mm,6  Nc,1  Nc,2  Nc,3  Nc,4  Nc,5  Nc,6  M⊙  pc    M⊙ yr − 1  M⊙  M⊙  M⊙  M⊙  M⊙  M⊙              105  0.1  256  0.03  128 698  122 888  73 463  88 127  13 167  81 490  238  223  196  197  52  186  104  0.1  256  0.03  52 567  60 425  4069  7997  627  986  242  231  110  124  15  13  3 × 104  0.1  256  0.03  68 380  91 207  18 737  26 004  2240  14 670  239  234  162  182  19  81  3 × 105  0.1  256  0.03  165 820  133 309  138 241  179 454  111 502  131 137  231  196  216  198  166  205  106  0.1  256  0.03  179 746  172 681  178 408  172 964  158 720  142 053  230  191  219  191  209  187  105  0.01  256  0.03  10 707  9449  11 734  9343  12 557  12 406  239  203  239  205  239  210  105  0.03  256  0.03  48 319  43 226  32 766  33 458  36 281  31 088  239  211  226  204  209  191  105  0.3  256  0.03  253 063  362 374  53 793  79 250  5796  47 880  230  230  140  150  16  81  105  1.0  256  0.03  701 633  1 020 974  38 607  45 433  2345  4092  223  207  100  59  7  5  105  0.1  64  0.03  48 671  48 945  32 849  56 258  38 327  44 628  53  48  40  44  42  40  105  0.1  128  0.03  87 493  85 516  77 464  73 564  16 884  42 905  114  88  105  81  43  69  105  0.1  512  0.03  155 664  181 666  68 136  86 562  7507  70 803  500  439  371  362  52  319  105  0.1  1024  0.03  226 161  275 657  74 954  91 700  6101  57 747  1004  985  792  752  68  455  105  0.1  256  0.001  23 649  36 107  11 059  28 385  10 773  19 905  62  114  31  91  35  68  105  0.1  256  0.003  58 318  65 221  24 661  47 494  11 610  45 009  104  188  68  149  43  133  105  0.1  256  0.01  97 024  94 974  35 739  79 289  8204  54 539  174  210  103  191  26  171  105  0.1  256  0.1  195 809  289 809  78 644  88 052  20 661  79 274  245  245  206  192  77  170  105  0.1  256  0.3  244 229  309 201  53 819  89 040  13 292  73 790  245  250  149  174  51  136  Mg  Rg  N  $$\dot{m}$$  Mm,1  Mm,2  Mm,3  Mm,4  Mm,5  Mm,6  Nc,1  Nc,2  Nc,3  Nc,4  Nc,5  Nc,6  M⊙  pc    M⊙ yr − 1  M⊙  M⊙  M⊙  M⊙  M⊙  M⊙              105  0.1  256  0.03  128 698  122 888  73 463  88 127  13 167  81 490  238  223  196  197  52  186  104  0.1  256  0.03  52 567  60 425  4069  7997  627  986  242  231  110  124  15  13  3 × 104  0.1  256  0.03  68 380  91 207  18 737  26 004  2240  14 670  239  234  162  182  19  81  3 × 105  0.1  256  0.03  165 820  133 309  138 241  179 454  111 502  131 137  231  196  216  198  166  205  106  0.1  256  0.03  179 746  172 681  178 408  172 964  158 720  142 053  230  191  219  191  209  187  105  0.01  256  0.03  10 707  9449  11 734  9343  12 557  12 406  239  203  239  205  239  210  105  0.03  256  0.03  48 319  43 226  32 766  33 458  36 281  31 088  239  211  226  204  209  191  105  0.3  256  0.03  253 063  362 374  53 793  79 250  5796  47 880  230  230  140  150  16  81  105  1.0  256  0.03  701 633  1 020 974  38 607  45 433  2345  4092  223  207  100  59  7  5  105  0.1  64  0.03  48 671  48 945  32 849  56 258  38 327  44 628  53  48  40  44  42  40  105  0.1  128  0.03  87 493  85 516  77 464  73 564  16 884  42 905  114  88  105  81  43  69  105  0.1  512  0.03  155 664  181 666  68 136  86 562  7507  70 803  500  439  371  362  52  319  105  0.1  1024  0.03  226 161  275 657  74 954  91 700  6101  57 747  1004  985  792  752  68  455  105  0.1  256  0.001  23 649  36 107  11 059  28 385  10 773  19 905  62  114  31  91  35  68  105  0.1  256  0.003  58 318  65 221  24 661  47 494  11 610  45 009  104  188  68  149  43  133  105  0.1  256  0.01  97 024  94 974  35 739  79 289  8204  54 539  174  210  103  191  26  171  105  0.1  256  0.1  195 809  289 809  78 644  88 052  20 661  79 274  245  245  206  192  77  170  105  0.1  256  0.3  244 229  309 201  53 819  89 040  13 292  73 790  245  250  149  174  51  136  View Large Table 3. Overview of the simulations. Continuing Table 2 presenting the maximum collision rate Rm, and the average collision rate Rav. Mg  Rg  N  $$\dot{m}$$  Rm,1  Rm,2  Rm,3  Rm,4  Rm,5  Rm,6  Rav,1  Rav,2  Rav,3  Rav,4  Rav,5  Rav,6  M⊙  pc    M⊙ yr − 1  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  105  0.1  256  0.03  21  22  14  19  4  14  1.60  3.37  0.32  1.29  0.36  0.30  104  0.1  256  0.03  145  125  8  31  1  8  14.35  20.17  0.20  0.08  0.02  0.43  3 × 104  0.1  256  0.03  55  69  16  62  2  43  9.66  10.52  0.15  0.14  0.05  4.12  3 × 105  0.1  256  0.03  13  9  8  17  6  11  0.87  1.35  0.78  0.32  0.42  0.50  106  0.1  256  0.03  7  5  7  7  6  6  0.39  0.26  0.38  0.26  0.35  0.45  105  0.01  256  0.03  11  15  14  13  12  9  1.47  0.69  0.65  0.78  0.55  0.28  105  0.03  256  0.03  15  13  21  17  13  12  0.67  0.48  1.58  1.01  0.30  0.71  105  0.3  256  0.03  47  55  9  49  2  30  5.12  4.11  0.16  0.17  0.04  4.22  105  1.0  256  0.03  83  71  7  23  2  3  7.16  5.58  0.30  0.21  0.10  0.79  105  0.1  64  0.03  5  5  6  9  3  3  0.45  0.53  0.48  0.38  0.26  0.25  105  0.1  128  0.03  10  9  7  13  5  7  0.77  0.38  0.90  0.64  0.52  0.66  105  0.1  512  0.03  81  57  48  59  2  54  5.66  3.86  2.36  2.91  0.28  2.18  105  0.1  1024  0.03  190  185  52  232  3  123  14.71  7.53  2.06  3.10  0.42  4.86  105  0.1  256  0.001  2  4  2  3  1  4  0.07  0.13  0.04  0.11  0.01  0.11  105  0.1  256  0.003  3  7  4  5  3  6  0.13  0.33  0.01  0.43  0.07  0.14  105  0.1  256  0.01  7  12  7  12  3  8  0.57  0.68  0.08  0.57  0.03  0.59  105  0.1  256  0.1  42  62  28  44  12  22  6.12  1.93  0.58  1.78  0.66  1.09  105  0.1  256  0.3  80  114  23  72  10  47  14.11  17.55  1.93  0.92  0.77  7.44  Mg  Rg  N  $$\dot{m}$$  Rm,1  Rm,2  Rm,3  Rm,4  Rm,5  Rm,6  Rav,1  Rav,2  Rav,3  Rav,4  Rav,5  Rav,6  M⊙  pc    M⊙ yr − 1  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  $$\rm {T_{cr,0}^{-1}}$$  105  0.1  256  0.03  21  22  14  19  4  14  1.60  3.37  0.32  1.29  0.36  0.30  104  0.1  256  0.03  145  125  8  31  1  8  14.35  20.17  0.20  0.08  0.02  0.43  3 × 104  0.1  256  0.03  55  69  16  62  2  43  9.66  10.52  0.15  0.14  0.05  4.12  3 × 105  0.1  256  0.03  13  9  8  17  6  11  0.87  1.35  0.78  0.32  0.42  0.50  106  0.1  256  0.03  7  5  7  7  6  6  0.39  0.26  0.38  0.26  0.35  0.45  105  0.01  256  0.03  11  15  14  13  12  9  1.47  0.69  0.65  0.78  0.55  0.28  105  0.03  256  0.03  15  13  21  17  13  12  0.67  0.48  1.58  1.01  0.30  0.71  105  0.3  256  0.03  47  55  9  49  2  30  5.12  4.11  0.16  0.17  0.04  4.22  105  1.0  256  0.03  83  71  7  23  2  3  7.16  5.58  0.30  0.21  0.10  0.79  105  0.1  64  0.03  5  5  6  9  3  3  0.45  0.53  0.48  0.38  0.26  0.25  105  0.1  128  0.03  10  9  7  13  5  7  0.77  0.38  0.90  0.64  0.52  0.66  105  0.1  512  0.03  81  57  48  59  2  54  5.66  3.86  2.36  2.91  0.28  2.18  105  0.1  1024  0.03  190  185  52  232  3  123  14.71  7.53  2.06  3.10  0.42  4.86  105  0.1  256  0.001  2  4  2  3  1  4  0.07  0.13  0.04  0.11  0.01  0.11  105  0.1  256  0.003  3  7  4  5  3  6  0.13  0.33  0.01  0.43  0.07  0.14  105  0.1  256  0.01  7  12  7  12  3  8  0.57  0.68  0.08  0.57  0.03  0.59  105  0.1  256  0.1  42  62  28  44  12  22  6.12  1.93  0.58  1.78  0.66  1.09  105  0.1  256  0.3  80  114  23  72  10  47  14.11  17.55  1.93  0.92  0.77  7.44  View Large Table 4. Comparing two different mass–radius parametrizations: Hosokawa et al. (2012) (subscript 1) and Haemmerlé et al. (2018) (subscript 2) (defined in Section 2.3). We compare their influence on the formation of a massive object: its final mass Mmax, number of collisions Ncol, and maximum collision rate Rmax. We compare runs with the exact same initial realization of the cluster, and calculate the average and standard deviation over 10 simulations. Models 2 and 4 give consistent results, which are the position-dependent accretion models where a dominant object is formed in the core. For the other models we observe a statistical difference due to the more conservative radii of Haemmerlé et al. (2018). This difference does not change the conclusion that a single massive object is formed due to a high collision rate. Model  Mmax,1  δMmax,1  Mmax,2  δMmax,2  Ncol,1  δNcol,1  Ncol,2  δNcol,2  Rmax,1  δRmax,1  Rmax,2  δRmax,2    M⊙  M⊙  M⊙  M⊙          $$\rm {t_{cross,0}^{-1}}$$  $$\rm {t_{cross,0}^{-1}}$$  $$\rm {t_{cross,0}^{-1}}$$  $$\rm {t_{cross,0}^{-1}}$$  1  120 529  13 127  173 130  13 856  235.3  6.67  182.5  6.17  23.3  5.48  14.2  3.01  2  139 439  16 503  160 317  16 480  228.6  5.30  225.4  3.60  28.9  4.12  25.7  6.60  3  74 613  3777  33 329  4637  201.5  6.00  87.4  10.60  16.8  3.68  5.4  1.43  4  84 285  4581  83 162  5962  189.1  13.74  181.1  15.88  21.7  3.47  19.6  4.58  5  13 381  1791  9238  1174  47.3  4.79  25.3  3.30  3.7  0.82  2.3  1.06  6  71 523  3729  56 040  2352  183.2  5.71  121.7  9.44  18.8  4.08  12.8  2.20  Model  Mmax,1  δMmax,1  Mmax,2  δMmax,2  Ncol,1  δNcol,1  Ncol,2  δNcol,2  Rmax,1  δRmax,1  Rmax,2  δRmax,2    M⊙  M⊙  M⊙  M⊙          $$\rm {t_{cross,0}^{-1}}$$  $$\rm {t_{cross,0}^{-1}}$$  $$\rm {t_{cross,0}^{-1}}$$  $$\rm {t_{cross,0}^{-1}}$$  1  120 529  13 127  173 130  13 856  235.3  6.67  182.5  6.17  23.3  5.48  14.2  3.01  2  139 439  16 503  160 317  16 480  228.6  5.30  225.4  3.60  28.9  4.12  25.7  6.60  3  74 613  3777  33 329  4637  201.5  6.00  87.4  10.60  16.8  3.68  5.4  1.43  4  84 285  4581  83 162  5962  189.1  13.74  181.1  15.88  21.7  3.47  19.6  4.58  5  13 381  1791  9238  1174  47.3  4.79  25.3  3.30  3.7  0.82  2.3  1.06  6  71 523  3729  56 040  2352  183.2  5.71  121.7  9.44  18.8  4.08  12.8  2.20  View Large In the next section, we provide several validation tests of our experimental set-up. We show that the initial condition is a stable Plummer sphere distribution (Fig. B1), and we validate both the correct implementation of the six different accretion models and the mass–radius parametrization based on Hosokawa et al. (2012). 4 VALIDATION OF THE NUMERICAL METHOD We describe the numerical methods used in this study in Section 2. Here, we present validation experiments of our simulation set-up and numerical implementation. We have performed a validation experiment of the initial conditions and dynamics by choosing the following parameters: Mg = 105 M⊙, Rg = 0.1 pc, and N = 256. We evolved the system for a time T = 105 yr and note that accretion and collisions are not included yet. In Fig. B1, we confirm that the 10, 50, and 90 per cent Lagrangian radii follow those of a Plummer sphere, with a slight discrepancy in the 90 per cent radius due to the truncation that we have introduced. We also note that all the protostars remain within the gas cloud, i.e. there are no escapers due to the initial condition, and the velocity dispersion of the protostars is close to the analytical value of ∼ 80 km s − 1. We have constructed a stable star cluster which will only alter its configuration due to gas accretion and stellar collisions. We have performed a numerical validation experiment by evolving the same initial condition as above, but this time only with gas accretion, i.e. dynamics is turned off. We set the initial accretion rate $$\dot{m} = 0.03\,\rm {M_{\odot }\,yr^{-1}}$$, i.e. here and in the rest of this study $$\dot{m}$$ refers to the average accretion rate per star over the whole cluster. In Fig. B2, we present the time evolution of the total star and gas mass (top row) and the time evolution of the average accretion rate (bottom row). We note that the position-dependent and -independent accretion models are consistent on average (their curves overlie). Together all these different models cover a variety of physical regimes. In the same validation experiment, we also kept track of the mass–radius evolution of the protostars, which we present in Fig. 2 (for the model of Hosokawa et al. 2012). There we clearly observe the difference between the position-independent (top row) and -dependent (bottom row) models. The latter produce a spectrum of masses and radii. We note that in the time-dependent models (5 and 6), the accretion rate decreases in time, and as a result the protostars will migrate to lower mass–radius tracks. After an initial phase of growth these protostars will eventually shrink, which is expected to decrease the collision rate. Figure 2. View largeDownload slide Validation of the mass–radius parametrization implementation. We consider a Plummer distribution of protostars, which are accreting gas, resulting in evolving masses and radii. The protostars do not move in this validation experiment, which explains the regularity in the distribution of mass–radius tracks. We confirm that for the position-independent accretion models (top row), the protostars follow the mass–radius parametrization based on Hosokawa et al. 2012. For the time-dependent accretion model (right-hand panel), the accretion rate decreases in time, so that the protostars migrate to mass–radius tracks of lower accretion rate, causing the protostars to eventually shrink again. For the position-dependent accretion models (bottom row), we confirm a correct interpolation of the mass–radius tracks in log-space, and the production of a range of masses and radii. Figure 2. View largeDownload slide Validation of the mass–radius parametrization implementation. We consider a Plummer distribution of protostars, which are accreting gas, resulting in evolving masses and radii. The protostars do not move in this validation experiment, which explains the regularity in the distribution of mass–radius tracks. We confirm that for the position-independent accretion models (top row), the protostars follow the mass–radius parametrization based on Hosokawa et al. 2012. For the time-dependent accretion model (right-hand panel), the accretion rate decreases in time, so that the protostars migrate to mass–radius tracks of lower accretion rate, causing the protostars to eventually shrink again. For the position-dependent accretion models (bottom row), we confirm a correct interpolation of the mass–radius tracks in log-space, and the production of a range of masses and radii. 5 RESULTS We start by analysing the simulations with the standard set of parameters defined in Section 3. We show a proof of concept that many accretion-induced collisions can lead to the formation of a single massive object. Next, we determine the sensitivity of this result to variations of the input parameters. 5.1 Simulations with the standard set of parameters Our simulations start with a cluster of single protostars, that is, there are no primordial binaries in the system. In time, collisions between protostars will occur and will result in collision products. We distinguish between the most massive collision product and other less massive ones. Finally, strong dynamical encounters among protostars can eject them from the cluster. We have thus defined four categories to which a protostar can belong: Single protostar, i.e. not a collision product Part of the most massive collision product Part of a less massive collision product Escaper. We note that binaries and higher order multiple systems form in our simulations, but we only categorize the individual stellar components. In Fig. 3 Figure 3. View largeDownload slide Time evolution of the fraction of stars belonging to four categories: (1) escapers (green); (2) stars that have collided with and are part of the most massive star in the system (blue, ‘Max. mass’); (3) stars that are part of other collision products (red, ‘Coll. prod.’); and (4) single stars (orange). We indicate the moment of transition to a stellar mass dominated system (see Section 2.2) with a vertical dashed line. Initially, all the stars are single stars in the system. As the stars accrete gas, they become larger and will dynamically interact with the other stars. We generally observe a decrease in the fraction of single stars, and an increase in collision products, and at late times a small fraction of escapers. Figure 3. View largeDownload slide Time evolution of the fraction of stars belonging to four categories: (1) escapers (green); (2) stars that have collided with and are part of the most massive star in the system (blue, ‘Max. mass’); (3) stars that are part of other collision products (red, ‘Coll. prod.’); and (4) single stars (orange). We indicate the moment of transition to a stellar mass dominated system (see Section 2.2) with a vertical dashed line. Initially, all the stars are single stars in the system. As the stars accrete gas, they become larger and will dynamically interact with the other stars. We generally observe a decrease in the fraction of single stars, and an increase in collision products, and at late times a small fraction of escapers. , we present the time evolution of the fraction of protostars belonging to each of these four categories. We observe that the majority of protostars end up in the final most massive object, which can be understood intuitively as follows. As a result of the gas accretion, the protostellar radii are substantially increasing for two main reasons: first, the overall increase of the protostellar mass, which naturally provides larger radii; and secondly, the rather high accretion rates in the primordial environment, which lead to rather extended envelopes that may exceed several 1000 solar radii. As a consequence, the collision probability in this environment is highly enhanced once accretion is switched on, and once the effect of accretion on to the protostellar radii is being considered. In particular for the infinite gas reservoir models (1 and 2), we find that 80–90 per cent of the protostars end up in a single object. For finite but large gas reservoirs (as expected in the large atomic cooling haloes), the fraction is somewhat reduced to 60–70 per cent in models 3, 4, and 6, but still quite significant. In model 5, the fraction reduces to about 10 per cent as due to the uniform accretion the central object in that cluster is less pronounced, and due to the time-dependent accretion, the protostellar radii shrink again with decreasing gas reservoir. The latter represents an enhancement by roughly a factor of 10 compared to the corresponding gas-free models or models that assumed smaller protostellar radii (Baumgardt & Klessen 2011; Reinoso et al. 2018). Comparing the radially independent accretion models (top row of Fig. 3) and the radially dependent models (bottom row), we observe a higher fraction of less massive collision products for the radially independent accretion models. Initially, this fraction is larger than the fraction of stars in the most massive collision product, implying that collisions tend to occur throughout the cluster between different pairs of stars. The radial-independent models also produce a larger fraction of escaping protostars, implying that these systems produce stronger dynamical encounters. For the radially dependent accretion models, most stars tend to immediately collide with the central protostar that is accreting fastest. Next, we aim to better understand the time evolution of the collision rate of protostars. In Figs 4 and 5, we present the collision rate (bottom panels) and correlate it with the total star and gas mass (top panels), the fraction of stars belonging to the four categories defined earlier (second panels) and the radius of the most massive protostar and the average radius of the remaining protostars (third panels). Figure 4. View largeDownload slide Correlating the time evolution of the collision and escape rate (bottom panel) with: total star and gas mass (top panel), fraction of stars belonging to the same four categories as in Fig. 3 (second panel), and maximum stellar radius and average stellar radius of the remaining stars (third panel). We show these results for the different accretion models (per column, see also next figure). Figure 4. View largeDownload slide Correlating the time evolution of the collision and escape rate (bottom panel) with: total star and gas mass (top panel), fraction of stars belonging to the same four categories as in Fig. 3 (second panel), and maximum stellar radius and average stellar radius of the remaining stars (third panel). We show these results for the different accretion models (per column, see also next figure). Figure 5. View largeDownload slide Same as the previous figure but for the position-dependent accretion models. Figure 5. View largeDownload slide Same as the previous figure but for the position-dependent accretion models. There is a small delay time for the collision rate to start to increase, which corresponds to the time it takes for the protostars to grow in mass and radius, and for the total stellar mass to become significant compared to the total gas mass. Then the collision rate increases rapidly and reaches a peak value. This rapid growth is fuelled by accretion, which results in the protostars reaching larger sizes and which makes the system stellar-mass dominated, resulting in dynamical encounters and collisions. Values for the maximum and average collision rates can be found in Table 3, which together specify a range of typical collision rates in our simulations. We find that the maximum collision rate can be very high, i.e. a collisional fraction per crossing time up to about 8 per cent (see left-hand panels in Figs 4 and 5), but that the average collision rate is more consistent with previous studies. For example, Baumgardt & Klessen (2011, fig. 4) measure a collisional fraction of 0.1–1 per cent, and this result is confirmed by Reinoso et al. (2018) for gas-free systems. Portegies Zwart et al. (2004) estimate the collision rate in the cluster MGG-11 to be 10–100 collisions within the first 3 Myr, or since the crossing time is about 105 yr (McCrady, Gilbert & Graham 2003), ∼0.3–3 collisions per crossing time, which is comparable to the average rates in Table 3. Our results show that accretion-induced collisions in massive Pop. III protostar clusters can increase the peak collision rate by an order of magnitude or even more. The saturation and eventual decrease of the collision rate have several causes. For the infinite gas reservoir models, the collision rate decreases again due to the sparsity of collision partners left in the system. For the finite gas reservoir models, the collision rate decreases due to a lack of gas, i.e. the stars do not expand anymore due to accretion. In this regime, we observe that the fraction of escaping stars starts to increase, which implies strong dynamical encounters. After the initial evolution with a high accretion-induced collision rate, the system experiences a transition into a regime with a dynamics-dominated collision rate, which is consequently much lower. For the time-dependent accretion models, the collision rate also decreases due to the shrinking of the protostars. We note however that even though the collision rate decreases, the collisions that do occur can be quite massive if both collision partners are collision products themselves. In Fig. 6, we plot the maximum mass in the system as a function of time. We observe that all models produce a very massive object with a mass of order 104 − 5 M⊙. The fastest growth in mass occurs between 5 and 20 crossing times, which is synchronous to the moment of highest collision rate (see Figs 4 and 5, bottom panels). For model 5, where accretion-induced collisions have the least effect, there is a longer delay time of about 30–50 crossing times, consistent with results from Reinoso et al. (2018, fig. 5). Figure 6. View largeDownload slide Time evolution of the maximum mass in the system. We show the results for the six different accretion models and standard set of parameters. Except for model 5, all models efficiently convert at least half of the initial gas mass into one single massive object. Figure 6. View largeDownload slide Time evolution of the maximum mass in the system. We show the results for the six different accretion models and standard set of parameters. Except for model 5, all models efficiently convert at least half of the initial gas mass into one single massive object. The results presented so far show that the formation of a central massive object is the natural outcome, in the regime of high-accretion rates and taking into account the large radii of accreting Pop. III protostars. It is favourable to have a position-dependent accretion rate, such that the formation of a massive object in the core is more likely, which fuels subsequent collisions with other stars. In model 5, such a formation channel is missing. This plus the fact that the stars will eventually shrink due to the time-dependent accretion rate causes the massive object to be an order of magnitude less massive . Model 5 comes closest to resembling the evolution of a gas-free, equal mass Plummer sphere in which collisions usually do not become important until core collapse. This implies that for systems in which accretion-induced collisions do play a significant role, the formation of a massive object is sped up. Next, we want to further investigate how the collisions occur in the cluster: do they preferably happen with a dominant central object, or are there collisions between different pairs of objects occurring throughout the cluster? In Fig. 7, we plot the number of collision products in the cluster as a function of the total number of collisions that have occurred. We observe that the position-independent accretion models (black curves) can produce factor 2–3 more collision products than the position-dependent models (grey curves). This is somewhat expected considering that the position-dependent accretion models rapidly produce a dominant mass in the core, which will have a much larger cross-section for collisions. We note that the mass ratio between the two most massive objects in the system ranges from an order of magnitude for model 5 up to two orders of magnitude for the other models. Figure 7. View largeDownload slide The number of collision products present in the cluster as a function of number of collisions that have occurred. This is a measure of whether collisions occur with predominantly one massive object or between many different pairs of stars throughout the cluster. We clearly observe a difference between the radially independent (thick curves) and radially dependent (thin curves) accretion models. As is somewhat expected, the radially dependent accretion models produce a massive object in the centre which is involved in most collisions. Figure 7. View largeDownload slide The number of collision products present in the cluster as a function of number of collisions that have occurred. This is a measure of whether collisions occur with predominantly one massive object or between many different pairs of stars throughout the cluster. We clearly observe a difference between the radially independent (thick curves) and radially dependent (thin curves) accretion models. As is somewhat expected, the radially dependent accretion models produce a massive object in the centre which is involved in most collisions. 5.2 Variation of input parameters We perform an ensemble of numerical experiments with varying input parameters to determine the robustness of the result in the previous section, i.e. that a single massive object of mass ∼ 104 − 5 M⊙ is formed due to a high collision rate. We vary the initial total gas mass, gas cloud radius, number of protostars, and average accretion rate per protostar. At the end of the simulation, we determine the maximum mass in the system and the total number of collisions. The results are presented in Fig. 8. Figure 8. View largeDownload slide Mass of the most massive object in the system at the end of the simulation, Mmax, and the total fraction of collisions, Nc/N0, as a function of the model input parameters: initial gas cloud mass (top left), gas cloud radius (top right), initial number of protostars (bottom left), and initial average accretion rate per protostar (bottom right). Figure 8. View largeDownload slide Mass of the most massive object in the system at the end of the simulation, Mmax, and the total fraction of collisions, Nc/N0, as a function of the model input parameters: initial gas cloud mass (top left), gas cloud radius (top right), initial number of protostars (bottom left), and initial average accretion rate per protostar (bottom right). Starting with the total gas mass, we observe that the maximum mass increases with higher gas mass. This makes sense intuitively as more gas is available to accrete from, and the resulting protostars will be more massive and larger in size. This increases the collisional cross-section, which is confirmed by the increasing number of collisions as a function of Mg. The trend seems to saturate at a gas mass of about 105.5 M⊙, and all the different models start to converge. In this regime, the protostars are able to reach such large radii that a runaway collision process is triggered, irrespective of the detailed accretion model. At the low mass end of the diagram, we observe that the time-dependent accretion models become inefficient, producing maximum masses ≤ 103 M⊙. With little mass to accrete from, the protostars in the time-dependent models will generally remain small, resulting in a low collision rate, which is consistent with gas-free star clusters in virial equilibrium. In the top right-hand panel of Fig. 8, we present the maximum mass and number of collisions with varying gas cloud radius. For the infinite gas reservoir models, we observe an increasing trend. For a larger gas cloud radius, the protostars need to accrete more gas before they reach the required radius for collisions to occur. As a result, the collision products will be more massive. For the finite gas reservoir models, there is a peak at Rg ∼ 10 − 1.5-10 − 1 pc after which the maximum mass decreases again. Intuitively, we can relate this decrease to the decrease in the number density, which is proportional to the collision rate. This is confirmed by the decrease of the number of collisions with Rg. For the densest configurations, i.e. the smallest gas cloud radii, we again observe a convergence of the different accretion models. With decreasing Rg, the ratio of the protostellar radii to Rg increases, resulting in an increased collision rate, which is confirmed by the data. The trend that Mmax increases with Rg, for the densest configurations, is explained by the fact that the protostars have to accrete for a longer time to reach a significant radius relative to Rg, thus resulting in more massive objects. The initial increase and the subsequent decrease of Mmax with increasing Rg are thus explained by the increase of the accretion time-scale and the decrease of the number density. In the bottom left-hand panel of Fig. 8, we vary the number of initial protostars in the cluster. Intuitively, we would expect that with a higher number of stars, the collision rate will increase. However, with fewer stars in the cluster, they will be able to accrete more and thus increase their collisional cross-section. For the infinite gas reservoir models, we again observe a steady increase of the maximum mass with N. The finite models produce massive objects between 104.5 − 5 M⊙ irrespective of the variation of N from 64 to 1024. The total collision fraction also remains roughly constant between 0.6 and 0.8. Therefore, the increasing number density and the decreasing cross-section with N tend to cancel out. Only for model 5, we observe a decrease of Mmax and Nc/N0 with N. The extra effect of the time-dependent accretion rate and subsequent shrinking of the protostellar radii causes the decreasing collisional cross-section to dominate. For the time-dependent accretion models, it is therefore favourable to have fewer stars in the system in order to produce a massive object. In the bottom right-hand panel of Fig. 8, we plot the maximum mass and total collision fraction as a function of average accretion rate per protostar. We observe an increasing collision fraction as a function of $$\dot{m}$$. This can be explained by the mass–radius tracks: a higher $$\dot{m}$$ produces larger protostellar radii. For the largest values of $$\dot{m}$$ however, the variation in the protostellar radii saturates (see Fig. 1), and we observe a flattening of the total collision fraction. The maximum mass depends weakly on $$\dot{m}$$, showing a relatively mild increase. For our standard set of parameter values defined in Section 3, we conclude that for $$\dot{m} \ge 10^{-3}\,\rm {M}_{\odot }\, \rm {yr^{-1}}$$ a massive object of M ≥ 104 M⊙ is formed. Finally, we measure the dependence of the formation of a massive object to a different underlying mass–radius parametrization. We compare our parametrization based on Hosokawa et al. (2012) to that of Haemmerlé et al. (2018). We use the standard set of parameters defined in Section 3 and create 10 initial realizations of the cluster. We integrate these initial conditions using both mass–radius parametrizations and calculate the average and standard deviation of the final maximum mass, number of collisions, and maximum collision rate (see Table 4). For models 2 and 4, we find the results to be mutually consistent. In these radial-dependent accretion models, a dominant central object is formed in the core. Its rapid growth causes many collisions to occur irrespective of the detailed mass–radius parametrization. For models 3, 5, and 6, we observe that the models based on Haemmerlé et al. (2018) produce somewhat fewer collisions and less massive objects. These models are sensitive to the more conservative radii of the Haemmerlé et al. (2018) parametrization because the accretion model is radially independent (so no dominant object in the core due to accretion) and/or because the accretion model is time dependent, in which protostars shrink to small sizes and differences in the radii will have a large effect on the dynamics-dominated collision rate. For model 1, we observe the counter-intuitive result that the parametrization by Haemmerlé et al. (2018) produces fewer collisions but a more massive object. This can be explained by the difference in time at which the stopping condition was fulfilled. The collisions rate is lower, and the collisions therefore more spread out in time. In this infinite gas reservoir model, this gives the protostars more time to accrete mass. 6 CONCLUSION We present the first numerical study of the formation of massive black hole seeds through the formation channel of accretion and collisions in a dense Pop. III protostellar cluster. We take into account accurate mass–radius parametrizations for accreting Pop. III protostars, which are based on detailed stellar evolution calculations performed by Hosokawa et al. (2012) and Haemmerlé et al. (2018). These studies have shown that for high accretion rates ($$\dot{m} \ge 10^{-3}\,\rm {M}_{\odot }\, \rm {yr^{-1}}$$), the protostellar radii can become of order 102 − 4 R⊙, thus significantly increasing the collisional cross-section. Using the amuse simulation framework, we perform a series of multiphysics simulations, including N-body dynamics, an analytical gas potential, six different accretion models, mass–radius parametrizations, and stellar collisions. In Tables 2 and 3, we present an overview of the simulation input parameters and outcome statistics. Our most conservative model (model 5), which most closely resembles the dynamics of ordinary Plummer spheres, confirms the fractional collision rate of 0.1–1 per cent per crossing time found in previous studies. Depending on the values of the input parameters, about 10 per cent of the initial protostellar population collides into a single massive object of order 104 M⊙ on a time-scale of 50 crossing times or longer. In our remaining accretion models, the effects of accretion-induced collisions are more prominent for two distinct but related reasons. First, in case of a more extensive gas reservoir the protostars will have more time to accrete gas and as a result be able to reach larger radii. Second, massive objects form in the centre of the cluster where gas volume densities, and thus the accretion rates and collision rates, are highest. Our results show that the maximum number of collisions per crossing time can be increased up to 1–10 per cent of the initial protostellar population. This increased collision rate leads to the formation of a single massive object with a mass of order M = 104 − 5 M⊙. The increase of the collision rate is fuelled by the high accretion rates of the Pop. III protostars and the resulting large radii. After the maximum collision rate is reached, the collision rate decreases again, either due to the sparsity of collision partners left in the system, or due to the lack of gas, which truncates the growth of the protostars and allows for dynamical ejections from the cluster. In Figs 4–6, we also demonstrate that the formation of a massive object can occur within 20 crossing times, which is about twice as fast compared to gas-free models (e.g. Reinoso et al. 2018). We have varied the initial gas cloud mass and gas cloud radius, and confirm the intuitive result that most collisions occur in the densest systems (see Fig. 8, top two panels). Also, the final mass of the seed black hole increases with higher gas cloud mass. Our calculations based on models with a limited gas reservoir show a characteristic gas cloud radius of Rg = 10 − 1.5-10 − 1 pc, at which the final mass of the most massive object is a maximum, i.e. M ∼ 104.5 − 5 M⊙. For smaller gas cloud radii, the system becomes more dense, which triggers a runaway collision between the protostars at earlier times. The less time there is for the protostars to accrete gas, the less massive will be the collision products. In the opposite regime of larger gas cloud radii, a runaway collision might not occur as the protostars cannot reach the required radii. These systems eventually become stellar-mass dominated, and since the collision rate is proportional to the number density, we would expect fewer collisions in larger systems. The final mass of the seed black hole depends mildly on the initial number of protostars (ranging from N = 64–1024) and average accretion rate (ranging from $$\dot{m} = 0.001 - 0.3\,\rm {M}_{\odot }\,\rm {yr^{-1}}$$), consistently producing massive objects of order M = 104 M⊙ or higher. The mass–radius evolution of accreting Pop. III protostars is very uncertain. In order to estimate the associated uncertainty and the robustness of the results mentioned above, we compare two mass–radius parametrizations: one based on Hosokawa et al. (2012) and another more conservative model based on Haemmerlé et al. (2018). Our calculations confirm that a parametrization with more conservative radii results in lower collision rates. However, the two different mass–radius parametrizations, each based on detailed stellar evolution calculations, produce no strong differences in the final mass of the seed black hole. We conclude that the mechanism of accretion-induced collisions in dense Pop. III protostellar systems is a viable mechanism for explaining the formation of the first massive black hole seeds. Therefore, this investigation warrants follow-up studies, which improve on the realism of the detailed implementation and coupling of the different physics ingredients. Acknowledgements TB acknowledges support from Fundação para a Ciência e a Tecnologia (grant SFRH/BPD/122325/2016), and support from Center for Research & Development in Mathematics and Applications (CIDMA) (strategic project UID/MAT/04106/2013), and from ENGAGE SKA, POCI-01-0145-FEDER-022217, funded by COMPETE 2020 and FCT, Portugal. DRGS and BR thank for funding through Fondecyt regular (project code 1161247). BR thanks Conicyt for financial support (CONICYT-PFCHA/MagísterNacional/2017-22171385). DRGS, MF, and AMS acknowledge funding through the ‘Concurso Proyectos Internacionales de Investigación, Convocatoria 2015’ (project code PII20150171) and the BASAL Centro de Astrofísica y Tecnologías Afines (CATA) PFB-06/2007. DRGS is further grateful for funding via ALMA-Conicyt (project code 31160001), Quimal (project code QUIMAL170001), and the Anillo program ‘Formation and Growth of Supermassive Black Holes’ (project code ACT172033). RSK acknowledges financial support from the Deutsche Forschungsgemeinschaft via SFB 881 ‘The Milky Way System’ (sub-projects B1, B2, and B8) and SPP 1573 ‘Physics of the Interstellar Medium’. LH and RSK were supported by the European Research Council under the European Community’s Seventh Framework Programme (FP7/2007–2013) via the ERC Advanced Grant ‘STARLIGHT: Formation of the First Stars’ (project number 339177). Footnotes 1 http://www.amusecode.org/ REFERENCES Agarwal B., Khochfar S., 2015, MNRAS , 446, 160 https://doi.org/10.1093/mnras/stu1973 CrossRef Search ADS   Agarwal B., Cullen F., Khochfar S., Klessen R. S., Glover S. C. O., Johnson J., 2017, MNRAS , 468, L82 https://doi.org/10.1093/mnrasl/slx028 CrossRef Search ADS   Agarwal B., Regan J., Klessen R. S., Downes T. P., Zackrisson E., 2017, MNRAS , 470, 4034 https://doi.org/10.1093/mnras/stx1528 CrossRef Search ADS   Bañados E. et al.  , 2018, Nature , 553, 473 CrossRef Search ADS PubMed  Baumgardt H., Klessen R. S., 2011, MNRAS , 413, 1810 https://doi.org/10.1111/j.1365-2966.2011.18258.x CrossRef Search ADS   Begelman M. C., Shlosman I., 2009, ApJL , 702, L5 https://doi.org/10.1088/0004-637X/702/1/L5 CrossRef Search ADS   Begelman M. C., Volonteri M., Rees M. J., 2006, MNRAS , 370, 289 https://doi.org/10.1111/j.1365-2966.2006.10467.x CrossRef Search ADS   Boekholt T. C. N., Stutz A. M., Fellhauer M., Schleicher D. R. G., Matus Carrillo D. R., 2017, MNRAS , 471, 3590 https://doi.org/10.1093/mnras/stx1821 CrossRef Search ADS   Bovino S., Grassi T., Schleicher D. R. G., Latif M. A., 2014, ApJL , 790, L35 https://doi.org/10.1088/2041-8205/790/2/L35 CrossRef Search ADS   Bovino S., Grassi T., Schleicher D. R. G., Banerjee R., 2016, ApJ , 832, 154 https://doi.org/10.3847/0004-637X/832/2/154 CrossRef Search ADS   Bromm V., Loeb A., 2003, ApJ , 596, 34 https://doi.org/10.1086/377529 CrossRef Search ADS   Clark P. C., Glover S. C. O., Klessen R. S., 2008, ApJ , 672, 757 https://doi.org/10.1086/524187 CrossRef Search ADS   Clark P. C., Glover S. C. O., Smith R. J., Greif T. H., Klessen R. S., Bromm V., 2011, Science , 331, 1040 https://doi.org/10.1126/science.1198027 CrossRef Search ADS PubMed  D’Amico G., Panci P., Lupi A., Bovino S., Silk J., 2018, MNRAS , 473, 328 https://doi.org/10.1093/mnras/stx2419 CrossRef Search ADS   Devecchi B., Volonteri M., 2009, ApJ , 694, 302 https://doi.org/10.1088/0004-637X/694/1/302 CrossRef Search ADS   Devecchi B., Volonteri M., Rossi E. M., Colpi M., Portegies Zwart S., 2012, MNRAS , 421, 1465 https://doi.org/10.1111/j.1365-2966.2012.20406.x CrossRef Search ADS   Dijkstra M., Ferrara A., Mesinger A., 2014, MNRAS , 442, 2036, https://doi.org/10.1093/mnras/stu1007 CrossRef Search ADS   Dopcke G., Glover S. C. O., Clark P. C., Klessen R. S., 2013, ApJ , 766, 103 https://doi.org/10.1088/0004-637X/766/2/103 CrossRef Search ADS   Eggenberger P., Meynet G., Maeder A., Hirschi R., Charbonnel C., Talon S., Ekström S., 2008, Astrophys. Space Sci. , 316, 43 https://doi.org/10.1007/s10509-007-9511-y CrossRef Search ADS   Flaugher B. et al.  , 2015, AJ , 150, 150 CrossRef Search ADS   Fujii M. S., Portegies Zwart S., 2013, MNRAS , 430, 1018 https://doi.org/10.1093/mnras/sts673 CrossRef Search ADS   Fujii M., Iwasawa M., Funato Y., Makino J., 2007, PASJ , 59, 1095 https://doi.org/10.1093/pasj/59.6.1095 CrossRef Search ADS   Gallerani S., Fan X., Maiolino R., Pacucci F., 2017, PASA , 34, e022 https://doi.org/10.1017/pasa.2017.14 CrossRef Search ADS   Greif T. H., Springel V., White S. D. M., Glover S. C. O., Clark P. C., Smith R. J., Klessen R. S., Bromm V., 2011, ApJ , 737, 75 https://doi.org/10.1088/0004-637X/737/2/75 CrossRef Search ADS   Haemmerlé L., Woods T. E., Klessen R. S., Heger A., Whalen D. J., 2018, MNRAS , 474, 2757 CrossRef Search ADS   Hirano S., Hosokawa T., Yoshida N., Omukai K., Yorke H. W., 2015, MNRAS , 448, 568 https://doi.org/10.1093/mnras/stv044 CrossRef Search ADS   Hosokawa T., Omukai K., 2009, ApJ , 691, 823 https://doi.org/10.1088/0004-637X/691/1/823 CrossRef Search ADS   Hosokawa T., Yorke H. W., Omukai K., 2010, ApJ , 721, 478 https://doi.org/10.1088/0004-637X/721/1/478 CrossRef Search ADS   Hosokawa T., Omukai K., Yorke H. W., 2012, ApJ , 756, 93 https://doi.org/10.1088/0004-637X/756/1/93 CrossRef Search ADS   Hosokawa T., Yorke H. W., Inayoshi K., Omukai K., Yoshida N., 2013, ApJ , 778, 178 https://doi.org/10.1088/0004-637X/778/2/178 CrossRef Search ADS   Hosokawa T., Hirano S., Kuiper R., Yorke H. W., Omukai K., Yoshida N., 2016, ApJ , 824, 119 https://doi.org/10.3847/0004-637X/824/2/119 CrossRef Search ADS   Hut P., Makino J., McMillan S., 1995, ApJL , 443, L93 https://doi.org/10.1086/187844 CrossRef Search ADS   Katz H., Sijacki D., Haehnelt M. G., 2015, MNRAS , 451, 2352 https://doi.org/10.1093/mnras/stv1048 CrossRef Search ADS   Koushiappas S. M., Bullock J. S., Dekel A., 2004, MNRAS , 354, 292 https://doi.org/10.1111/j.1365-2966.2004.08190.x CrossRef Search ADS   Latif M. A., Schleicher D. R. G., 2015a, MNRAS , 449, 77 https://doi.org/10.1093/mnras/stu2573 CrossRef Search ADS   Latif M. A., Schleicher D. R. G., 2015b, A&A , 578, A118 https://doi.org/10.1051/0004-6361/201525855 CrossRef Search ADS   Latif M. A., Schleicher D. R. G., Schmidt W., Niemeyer J., 2013a, MNRAS , 433, 1607 https://doi.org/10.1093/mnras/stt834 CrossRef Search ADS   Latif M. A., Schleicher D. R. G., Schmidt W., Niemeyer J. C., 2013b, MNRAS , 436, 2989 https://doi.org/10.1093/mnras/stt1786 CrossRef Search ADS   Latif M. A., Schleicher D. R. G., Bovino S., Grassi T., Spaans M., 2014, ApJ , 792, 78 https://doi.org/10.1088/0004-637X/792/1/78 CrossRef Search ADS   Latif M. A., Bovino S., Van Borm C., Grassi T., Schleicher D. R. G., Spaans M., 2014, MNRAS , 443, 1979, https://doi.org/10.1093/mnras/stu1230 CrossRef Search ADS   Latif M. A., Bovino S., Grassi T., Schleicher D. R. G., Spaans M., 2015, MNRAS , 446, 3163 https://doi.org/10.1093/mnras/stu2244 CrossRef Search ADS   Latif M. A., Omukai K., Habouzit M., Schleicher D. R. G., Volonteri M., 2016, ApJ , 823, 40 https://doi.org/10.3847/0004-637X/823/1/40 CrossRef Search ADS   Lodato G., Natarajan P., 2006, MNRAS , 371, 1813 https://doi.org/10.1111/j.1365-2966.2006.10801.x CrossRef Search ADS   Lupi A., Colpi M., Devecchi B., Galanti G., Volonteri M., 2014, MNRAS , 442, 3616 https://doi.org/10.1093/mnras/stu1120 CrossRef Search ADS   Martínez-Barbosa C. A., Brown A. G. A., Boekholt T., Portegies Zwart S., Antiche E., Antoja T., 2016, MNRAS , 457, 1062 https://doi.org/10.1093/mnras/stw006 CrossRef Search ADS   Matsuoka Y. et al.  , 2017, preprint (ArXiv e-prints:1704.05854) McCrady N., Gilbert A. M., Graham J. R., 2003, ApJ , 596, 240 https://doi.org/10.1086/377631 CrossRef Search ADS   McMillan S. L. W., Hut P., 1996, ApJ , 467, 348 https://doi.org/10.1086/177610 CrossRef Search ADS   Moeckel N., Clarke C. J., 2011, MNRAS , 410, 2799 https://doi.org/10.1111/j.1365-2966.2010.17659.x CrossRef Search ADS   Mortlock D. J. et al.  , 2011, Nature , 474, 616 https://doi.org/10.1038/nature10159 CrossRef Search ADS PubMed  Oh S., Kroupa P., 2012, MNRAS , 424, 65 https://doi.org/10.1111/j.1365-2966.2012.21152.x CrossRef Search ADS   Omukai K., Palla F., 2003, ApJ , 589, 677 https://doi.org/10.1086/374810 CrossRef Search ADS   Omukai K., Schneider R., Haiman Z., 2008, ApJ , 686, 801 https://doi.org/10.1086/591636 CrossRef Search ADS   Pelupessy F. I., van Elteren A., de Vries N., McMillan S. L. W., Drost N., Portegies Zwart S. F., 2013, A&A , 557, A84 https://doi.org/10.1051/0004-6361/201321252 CrossRef Search ADS   Plummer H. C., 1911, MNRAS , 71, 460 https://doi.org/10.1093/mnras/71.5.460 CrossRef Search ADS   Portegies Zwart S. F., Baumgardt H., Hut P., Makino J., McMillan S. L. W., 2004, Nature , 428, 724 https://doi.org/10.1038/nature02448 CrossRef Search ADS PubMed  Portegies Zwart S. F., McMillan S. L., van Elteren A., Pelupessy F. I., de Vries N., 2013, Comput. Phys. Commun. , 184, 456 https://doi.org/10.1016/j.cpc.2012.09.024 CrossRef Search ADS   Portegies Zwart S. et al.   2009, New Astron. , 14, 369 https://doi.org/10.1016/j.newast.2008.10.006 CrossRef Search ADS   Reed S. L. et al.  , 2017, MNRAS , 468, 4702 https://doi.org/10.1093/mnras/stx728 CrossRef Search ADS   Rees M. J., 1984, ARA&A , 22, 471 CrossRef Search ADS   Regan J. A., Haehnelt M. G., 2009, MNRAS , 396, 343 https://doi.org/10.1111/j.1365-2966.2009.14579.x CrossRef Search ADS   Regan J. A., Johansson P. H., Wise J. H., 2016, MNRAS , 459, 3377 https://doi.org/10.1093/mnras/stw899 CrossRef Search ADS   Reinoso B., Schleicher D., Fellhauer M., Klessen R., Boekholt T. C. N., 2018, preprint (arXiv:1801.05891) Sakurai Y., Yoshida N., Fujii M. S., Hirano S., 2017, MNRAS , 472, 1677 https://doi.org/10.1093/mnras/stx2044 CrossRef Search ADS   Schleicher D. R. G., Palla F., Ferrara A., Galli D., Latif M., 2013, A&A , 558, A59 CrossRef Search ADS   Schleicher D. R. G., Spaans M., Glover S. C. O., 2010, ApJL , 712, L69 https://doi.org/10.1088/2041-8205/712/1/L69 CrossRef Search ADS   Schneider R., Omukai K., Inoue A. K., Ferrara A., 2006, MNRAS , 369, 1437 https://doi.org/10.1111/j.1365-2966.2006.10391.x CrossRef Search ADS   Shapiro S. L., 2005, ApJ , 620, 59 https://doi.org/10.1086/427065 CrossRef Search ADS   Smith R. J., Hosokawa T., Omukai K., Glover S. C. O., Klessen R. S., 2012, MNRAS , 424, 457 https://doi.org/10.1111/j.1365-2966.2012.21211.x CrossRef Search ADS   Sugimura K., Omukai K., Inoue A. K., 2014, MNRAS , 445, 544 https://doi.org/10.1093/mnras/stu1778 CrossRef Search ADS   Sutherland W. et al.  , 2015, A&A , 575, A25 CrossRef Search ADS   Umeda H., Hosokawa T., Omukai K., Yoshida N., 2016, ApJL , 830, L34 https://doi.org/10.3847/2041-8205/830/2/L34 CrossRef Search ADS   Wise J. H., Turk M. J., Abel T., 2008, ApJ , 682, 745 https://doi.org/10.1086/588209 CrossRef Search ADS   Woods T. E., Heger A., Whalen D. J., Haemmerlé L., Klessen R. S., 2017, ApJL , 842, L6 https://doi.org/10.3847/2041-8213/aa7412 CrossRef Search ADS   Wright E. L. et al.  , 2010, AJ , 140, 1868 CrossRef Search ADS   APPENDIX A: ANALYTICAL FORM OF THE MASS–RADIUS RELATION The analytical power-law parametrization of the mass–radius relation from the models of Hosokawa et al. (2012) and Haemmerlé et al. (2018) is given in Tables A1 and A2, respectively. Table A1. Simplified parametrization of the mass–radius relation of accreting Pop. III protostars based on Hosokawa et al. (2012, fig. 5). $$\dot{\rm {m}}$$  m  R  $$\rm {M}_{\odot }\,\rm {yr}^{-1}$$  M⊙  R⊙  ≥ 1  < 150  $$\rm {250} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.4}}$$    < 450  $$\rm {1855} \left( {\rm {m} \over \rm {150}\,\rm {M}_{\odot }} \right)^{\rm {0.8}}$$    ≥ 450  $$\rm {4468} \left( {\rm {m} \over \rm {450}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\rm {0.3} \le \dot{m} < \rm {1}$$  < 80  $$\rm {190} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.4}}$$    < 140  $$\rm {1096} \left( {\rm {m} \over \rm {80}\,\rm {M}_{\odot }} \right)^{\rm {1.5}}$$    ≥ 140  $$\rm {2538} \left( {\rm {m} \over \rm {140}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\rm {0.1} \le \dot{m} < \rm {0.3}$$  < 40  $$\rm {140} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.4}}$$    < 90  $$\rm {612} \left( {\rm {m} \over \rm {40}\,\rm {M}_{\odot }} \right)^{\rm {1.5}}$$    ≥ 90  $$\rm {2066} \left( {\rm {m} \over \rm {90}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\rm {0.06} \le \dot{m} < \rm {0.1}$$  < 25  $$\rm {110} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.4}}$$    < 70  $$\rm {399} \left( {\rm {m} \over \rm {25}\,\rm {M}_{\odot }} \right)^{\rm {1.5}}$$    ≥ 70  $$\rm {1868} \left( {\rm {m} \over \rm {70}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\rm {0.03} \le \dot{m} < \rm {0.06}$$  < 20  $$\rm {90} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.4}}$$    < 70  $$\rm {298} \left( {\rm {m} \over \rm {20}\,\rm {M}_{\odot }} \right)^{\rm {1.5}}$$    ≥ 70  $$\rm {1953} \left( {\rm {m} \over \rm {70}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\rm {0.006} \le \dot{m} < \rm {0.03}$$  < 20  $$\rm {50} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.4}}$$    < 35  $$\rm {166} \left( {\rm {m} \over \rm {20}\,\rm {M}_{\odot }} \right)^{\rm {-1.5}}$$    ≥ 35  $$\rm {72} \left( {\rm {m} \over \rm {35}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\dot{m} < \rm {0.006}$$  < 9  $$\rm {25} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.4}}$$    < 50  $$\rm {60} \left( {\rm {m} \over \rm {9}\,\rm {M}_{\odot }} \right)^{\rm {-1.5}}$$    ≥ 50  $$\rm {5} \left( {\rm {m} \over \rm {50}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\dot{\rm {m}}$$  m  R  $$\rm {M}_{\odot }\,\rm {yr}^{-1}$$  M⊙  R⊙  ≥ 1  < 150  $$\rm {250} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.4}}$$    < 450  $$\rm {1855} \left( {\rm {m} \over \rm {150}\,\rm {M}_{\odot }} \right)^{\rm {0.8}}$$    ≥ 450  $$\rm {4468} \left( {\rm {m} \over \rm {450}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\rm {0.3} \le \dot{m} < \rm {1}$$  < 80  $$\rm {190} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.4}}$$    < 140  $$\rm {1096} \left( {\rm {m} \over \rm {80}\,\rm {M}_{\odot }} \right)^{\rm {1.5}}$$    ≥ 140  $$\rm {2538} \left( {\rm {m} \over \rm {140}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\rm {0.1} \le \dot{m} < \rm {0.3}$$  < 40  $$\rm {140} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.4}}$$    < 90  $$\rm {612} \left( {\rm {m} \over \rm {40}\,\rm {M}_{\odot }} \right)^{\rm {1.5}}$$    ≥ 90  $$\rm {2066} \left( {\rm {m} \over \rm {90}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\rm {0.06} \le \dot{m} < \rm {0.1}$$  < 25  $$\rm {110} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.4}}$$    < 70  $$\rm {399} \left( {\rm {m} \over \rm {25}\,\rm {M}_{\odot }} \right)^{\rm {1.5}}$$    ≥ 70  $$\rm {1868} \left( {\rm {m} \over \rm {70}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\rm {0.03} \le \dot{m} < \rm {0.06}$$  < 20  $$\rm {90} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.4}}$$    < 70  $$\rm {298} \left( {\rm {m} \over \rm {20}\,\rm {M}_{\odot }} \right)^{\rm {1.5}}$$    ≥ 70  $$\rm {1953} \left( {\rm {m} \over \rm {70}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\rm {0.006} \le \dot{m} < \rm {0.03}$$  < 20  $$\rm {50} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.4}}$$    < 35  $$\rm {166} \left( {\rm {m} \over \rm {20}\,\rm {M}_{\odot }} \right)^{\rm {-1.5}}$$    ≥ 35  $$\rm {72} \left( {\rm {m} \over \rm {35}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\dot{m} < \rm {0.006}$$  < 9  $$\rm {25} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.4}}$$    < 50  $$\rm {60} \left( {\rm {m} \over \rm {9}\,\rm {M}_{\odot }} \right)^{\rm {-1.5}}$$    ≥ 50  $$\rm {5} \left( {\rm {m} \over \rm {50}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  View Large Table A2. Simplified parametrization of the mass–radius relation of accreting Pop. III protostars based on Haemmerlé et al. (2018, fig. 2). $$\dot{\rm {m}}$$  m  R  $$\rm {M}_{\odot }\,\rm {yr}^{-1}$$  M⊙  R⊙  ≥ 10  < 10  $$\rm {200} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.0}}$$    < 50  $$\rm {200} \left( {\rm {m} \over \rm {10}\,\rm {M}_{\odot }} \right)^{\rm {-0.4}}$$    < 60  $$\rm {105} \left( {\rm {m} \over \rm {50}\,\rm {M}_{\odot }} \right)^{\rm {13}}$$    ≥ 60  $$\rm {1124} \left( {\rm {m} \over \rm {60}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\rm {1} \le \dot{m} < \rm {10}$$  < 10  $$\rm {200} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.0}}$$    < 35  $$\rm {200} \left( {\rm {m} \over \rm {10}\,\rm {M}_{\odot }} \right)^{\rm {-0.6}}$$    < 42  $$\rm {94} \left( {\rm {m} \over \rm {35}\,\rm {M}_{\odot }} \right)^{\rm {13}}$$    ≥ 42  $$\rm {1009} \left( {\rm {m} \over \rm {42}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\rm {0.1} \le \dot{m} < \rm {1}$$  < 10  $$\rm {200} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.0}}$$    < 25  $$\rm {200} \left( {\rm {m} \over \rm {10}\,\rm {M}_{\odot }} \right)^{\rm {-0.6}}$$    < 29.5  $$\rm {115} \left( {\rm {m} \over \rm {25}\,\rm {M}_{\odot }} \right)^{\rm {13}}$$    ≥ 29.5  $$\rm {993} \left( {\rm {m} \over \rm {29.5}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\rm {0.01} \le \dot{m} < \rm {0.1}$$  < 30  $$\rm {200} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.0}}$$    < 70  $$\rm {200} \left( {\rm {m} \over \rm {30}\,\rm {M}_{\odot }} \right)^{\rm {-2.5}}$$    ≥ 70  $$\rm {24} \left( {\rm {m} \over \rm {70}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\dot{m} < \rm {0.01}$$  < 10  $$\rm {200} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.0}}$$    < 50  $$\rm {200} \left( {\rm {m} \over \rm {10}\,\rm {M}_{\odot }} \right)^{\rm {-2.5}}$$    ≥ 50  $$\rm {3.6} \left( {\rm {m} \over \rm {50}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\dot{\rm {m}}$$  m  R  $$\rm {M}_{\odot }\,\rm {yr}^{-1}$$  M⊙  R⊙  ≥ 10  < 10  $$\rm {200} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.0}}$$    < 50  $$\rm {200} \left( {\rm {m} \over \rm {10}\,\rm {M}_{\odot }} \right)^{\rm {-0.4}}$$    < 60  $$\rm {105} \left( {\rm {m} \over \rm {50}\,\rm {M}_{\odot }} \right)^{\rm {13}}$$    ≥ 60  $$\rm {1124} \left( {\rm {m} \over \rm {60}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\rm {1} \le \dot{m} < \rm {10}$$  < 10  $$\rm {200} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.0}}$$    < 35  $$\rm {200} \left( {\rm {m} \over \rm {10}\,\rm {M}_{\odot }} \right)^{\rm {-0.6}}$$    < 42  $$\rm {94} \left( {\rm {m} \over \rm {35}\,\rm {M}_{\odot }} \right)^{\rm {13}}$$    ≥ 42  $$\rm {1009} \left( {\rm {m} \over \rm {42}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\rm {0.1} \le \dot{m} < \rm {1}$$  < 10  $$\rm {200} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.0}}$$    < 25  $$\rm {200} \left( {\rm {m} \over \rm {10}\,\rm {M}_{\odot }} \right)^{\rm {-0.6}}$$    < 29.5  $$\rm {115} \left( {\rm {m} \over \rm {25}\,\rm {M}_{\odot }} \right)^{\rm {13}}$$    ≥ 29.5  $$\rm {993} \left( {\rm {m} \over \rm {29.5}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\rm {0.01} \le \dot{m} < \rm {0.1}$$  < 30  $$\rm {200} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.0}}$$    < 70  $$\rm {200} \left( {\rm {m} \over \rm {30}\,\rm {M}_{\odot }} \right)^{\rm {-2.5}}$$    ≥ 70  $$\rm {24} \left( {\rm {m} \over \rm {70}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  $$\dot{m} < \rm {0.01}$$  < 10  $$\rm {200} \left( {\rm {m} \over \rm {1}\,\rm {M}_{\odot }} \right)^{\rm {0.0}}$$    < 50  $$\rm {200} \left( {\rm {m} \over \rm {10}\,\rm {M}_{\odot }} \right)^{\rm {-2.5}}$$    ≥ 50  $$\rm {3.6} \left( {\rm {m} \over \rm {50}\,\rm {M}_{\odot }} \right)^{\rm {0.5}}$$  View Large APPENDIX B: VALIDATION OF THE NUMERICAL METHOD We present here additional plots for the validation experiments discussed in Section 4. In Fig. B1, we confirm that in the absence of accretion, the 10, 50, and 90 per cent Lagrangian radii follow those of a Plummer sphere, with a slight discrepancy in the 90 per cent radius due to the truncation that we have introduced. In Fig. B2, we present the time evolution of the total star and gas mass (top row) and the time evolution of the average accretion rate (bottom row) in a set-up where the N-body integrator has been switched off. The figure follows the expected behaviour of the accretion models. Finally, we show in Fig. B3 that the specific choice of the initial protostellar mass does not have a large influence on the time evolution of the most massive object in the system. The initial total stellar mass is given by Ms = Nm0, which becomes comparable to the initial gas cloud mass, Mg = 105 M⊙, and N = 256, if m0 = 390.625 M⊙. When m0 ≪ 390M⊙, the total mass (gas + stars) remains close to ∼ 105 M⊙ and, with a fixed collisional fraction of about 10 percent, produces a maximum mass of about 104M⊙. When m0 = 100 M⊙, the initial stellar mass Ms = 2.56 × 104M⊙, and so we have effectively increased the initial total mass by a quarter, which also allows for the formation of a more massive object. However, such a system would not start out as a gas-dominated system. Figure B1. View largeDownload slide Validation of the initial conditions and dynamics. We present the 10, 50, and 90 per cent Lagrangian radius of a Plummer sphere (see the text for the Plummer parameters), calculated analytically (dashed line, for an untruncated Plummer) and from our numerical model (solid line, truncated at 5 Plummer radii). We confirm that our cluster is initially stable and will only change its structure in the presence of accretion and collisions. Figure B1. View largeDownload slide Validation of the initial conditions and dynamics. We present the 10, 50, and 90 per cent Lagrangian radius of a Plummer sphere (see the text for the Plummer parameters), calculated analytically (dashed line, for an untruncated Plummer) and from our numerical model (solid line, truncated at 5 Plummer radii). We confirm that our cluster is initially stable and will only change its structure in the presence of accretion and collisions. Figure B2. View largeDownload slide Illustration and validation of the gas accretion models. In the top row, we show the total mass evolution of the stars and gas. We confirm that models 1 and 2 are consistent (top left-hand panel), and similar for models 3 and 4 (top middle panel) and models 5 and 6 (top right-hand panel). In the bottom row, we present the average accretion rates per star, i.e. the time derivative of the total stellar mass in the top row divided by the total number of stars. Figure B2. View largeDownload slide Illustration and validation of the gas accretion models. In the top row, we show the total mass evolution of the stars and gas. We confirm that models 1 and 2 are consistent (top left-hand panel), and similar for models 3 and 4 (top middle panel) and models 5 and 6 (top right-hand panel). In the bottom row, we present the average accretion rates per star, i.e. the time derivative of the total stellar mass in the top row divided by the total number of stars. Figure B3. View largeDownload slide Time evolution of the maximum mass in the system, for our most conservative accretion model (Model 5) with standard parameters (see Section 3). We vary the initial masses of the protostars and find that as long as the system starts out as a gas-dominated system, that is with m0 ≪ 390 M⊙, the final maximum mass does not change much. For m0 = 100 M⊙ the overall potential starts to become strongly influenced by the stars, and mass growth gets modified by Bondi–Hoyle–Littleton accretion. Figure B3. View largeDownload slide Time evolution of the maximum mass in the system, for our most conservative accretion model (Model 5) with standard parameters (see Section 3). We vary the initial masses of the protostars and find that as long as the system starts out as a gas-dominated system, that is with m0 ≪ 390 M⊙, the final maximum mass does not change much. For m0 = 100 M⊙ the overall potential starts to become strongly influenced by the stars, and mass growth gets modified by Bondi–Hoyle–Littleton accretion. © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society

### Journal

Monthly Notices of the Royal Astronomical SocietyOxford University Press

Published: May 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 12 million articles from more than
10,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Unlimited reading Read as many articles as you need. Full articles with original layout, charts and figures. Read online, from anywhere. ### Stay up to date Keep up with your field with Personalized Recommendations and Follow Journals to get automatic updates. ### Organize your research It’s easy to organize your research with our built-in tools. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. ### DeepDyve Freelancer ### DeepDyve Pro Price FREE$49/month

\$360/year
Save searches from Google Scholar, PubMed
Create lists to organize your research
Export lists, citations