<figure><img src="image_1.png" /></figure>

IMPACT OF DUST COOLING ON DIRECT-COLLAPSE BLACK HOLE FORMATION

M. A. LATIF¹,², K. OMUKAI³, M. HABOUZIT¹,², D. R. G. SCHLEICHER⁴, AND M. VOLONTERI¹,²

¹Sorbonne Universits, UPMC Univ Paris 06, UMR 7095, Institut d'Astrophysique de Paris, F-75014, Paris, France; latif@iap.fr

²CNRS, UMR 7095, Institut d'Astrophysique de Paris, F-75014, Paris, France

³Astronomical Institute, Tohoku University, Aoba, Sendai 980-8578, Japan

⁴Departamento de Astronoma, Facultad Ciencias Fsicas y Matemticas, Universidad de Concepcin, Av. Esteban Iturra s/n Barrio Universitario, Casilla 160-C, Chile
*Received 2015 November 13; accepted 2016 March 23; published 2016 May 20*

ABSTRACT

Observations of quasars at $z > 6$ suggest the presence of black holes with a few times $10^9 M_\odot$. Numerous models have been proposed to explain their existence, including a direct collapse, which provides massive seeds of $10^5 M_\odot$. The isothermal direct collapse requires a strong Lyman-Werner (LW) flux to quench $H_2$ formation in massive primordial halos. In this study, we explore the impact of trace amounts of metals and dust enrichment. We perform three-dimensional cosmological simulations for two halos of $>10^7 M_\odot$ with $Z/Z_\odot = 10^{-4}-10^{-6}$ illuminated by an intense LW flux of $J_{21} = 10^5$. Our results show that initially the collapse proceeds isothermally with $T \sim 8000$ K, but dust cooling becomes effective at densities of $10^8-10^{12} \text{ cm}^{-3}$ and brings the gas temperature down to a few 100-1000 K for $Z/Z_\odot \ge 10^{-6}$. No gravitationally bound clumps are found in the $Z/Z_\odot \le 10^{-5}$ cases by the end of our simulations, in contrast to the case with $Z/Z_\odot = 10^{-4}$. Large inflow rates of $\ge 0.1 M_\odot \text{ yr}^{-1}$ are observed for $Z/Z_\odot \le 10^{-5}$, similar to a zero-metallicity case, while for $Z/Z_\odot = 10^{-4}$ the inflow rate starts to decline earlier because of dust cooling and fragmentation. For given large inflow rates, a central star of $\sim 10^4 M_\odot$ may form for $Z/Z_\odot \le 10^{-5}$.

**Key words:** black hole physics – cosmology: theory – galaxies: formation – methods: numerical

1. INTRODUCTION

The discovery of high-redshift quasars at $z > 6$ reveals the presence of supermassive black holes (SMBHs) of about a few billion solar masses (Fan et al. 2006; Willott et al. 2010; Mortlock et al. 2011; Venemans et al. 2013, 2015; Wu et al. 2015). The formation of SMBHs a few hundred million years after the Big Bang presents a challenge for our understanding of early structure formation. The pathways for their assembly include the stellar-mass black holes forming from the collapse of Pop III stars (Abel et al. 2002; Yoshida et al. 2008; Latif et al. 2013c; Hirano et al. 2014), black hole seeds from the collapse of a dense nuclear cluster either because of a relativistic instability (Baumgarte & Shapiro 1999), via a stellar dynamical process (Omukai et al. 2008; Devecchi & Volonteri 2009; Katz et al. 2015; Yajima & Khochfar 2015), even the core collapse of a cluster consisting of stellar-mass black holes (Davies et al. 2011; Lupi et al. 2014), or a direct collapse of a protogalactic gas cloud (Loeb & Rasio 1994; Bromm & Loeb 2003; Begelman et al. 2006; Spaans & Silk 2006; Latif et al. 2013a, 2013d; Ferrara et al. 2014; Inayoshi et al. 2014; Choi et al. 2015; Shlosman et al. 2016). The details of these models are discussed in dedicated reviews (Volonteri 2010; Volonteri & Bellovary 2012; Haiman 2012).

The masses of seed black holes from the above-mentioned scenarios range from a few hundred to a thousand solar masses, except the direct-collapse model, which provides massive seeds of about $10^5-10^6 M_\odot$. The lighter seeds require continuous accretion close to the Eddington limit in order to reach a billion solar masses in a few hundred million years. Therefore, massive seeds forming via a direct collapse are favored for the assembly of high-redshift quasars. The formation of so-called direct-collapse black holes (DCBHs) requires large inflow rates of about $\ge 0.1 M_\odot \text{ yr}^{-1}$ (Begelman 2010; Hosokawa et al. 2013; Schleicher et al. 2013; Sakurai et al. 2015). These inflow rates can be obtained either via dynamical processes such as the "bars-

in-bars" instability (see Begelman et al. 2006; Begelman & Shlosman 2009) or thermodynamically because of a large speed of sound ($\dot{M} \sim c_s^3/G \sim 0.1 M_\odot \text{ yr}^{-1}(T/8000 \text{ K})^{3/2}$, where $c_s$ is the thermal sound speed). Large inflow rates can be more easily achieved in an isothermal direct collapse, which necessitates the presence of a strong Lyman-Werner (LW) flux to quench the formation of molecular hydrogen in metal-free halos (Omukai 2001; Schleicher et al. 2010; Shang et al. 2010; Latif et al. 2014). The critical value of the LW flux ($J_{crit}$) above which isothermal collapse occurs depends on the spectra of the stars (Omukai 2001; Shang et al. 2010; Latif et al. 2014, 2015; Regan et al. 2014b; Sugimura et al. 2014; Agarwal & Khochfar 2015). Latif et al. (2015) computed the strength of $J_{crit}$ for realistic Pop II spectra from three-dimensional (3D) cosmological simulations and found that it corresponds to a value of $J_{21}$ of a few times $10^4$ where $J_{21} = 1$ implies a flux of $10^{-21} \text{ erg cm}^{-2} \text{ s}^{-1} \text{ Hz}^{-1} \text{ sr}^{-1}$ below the Lyman limit. Another constraint for the isothermal direct collapse is that the halos should be metal-free because trace amounts of metals and dust can lead to fragmentation and star formation (Omukai et al. 2008).

Numerical simulations show that large inflow rates of about $0.1-1 M_\odot \text{ yr}^{-1}$ are available in the massive primordial halos of $10^7-10^8 M_\odot$ illuminated by a strong LW flux at $z = 15$, and massive objects of $10^5 M_\odot$ can be formed within about 100,000 years (Wise et al. 2008; Latif et al. 2013a, 2013d; Inayoshi et al. 2014; Regan et al. 2014a; Van Borm et al. 2014; Becerra et al. 2015). However, the value of $J_{crit}$ ($>10^4$ in terms of $J_{21}$) required for the formation of isothermal DCBHs is much higher than the background UV flux, and such a flux can only be achieved in the close vicinities of star-forming galaxies (Dijkstra et al. 2008, 2014; Agarwal et al. 2012; Visbal et al. 2014; Habouzit et al. 2016b; Latif et al. 2015). Supernova winds from nearby galaxies can also pollute halos with metals and make the sites for an isothermal
direct collapse even more rare (Johnson et al. 2013). In the context of second-generation star formation (Pop II), it has been found that above a critical value of metallicity ($Z_{crit} = 3 \times 10^{-4}Z_{⊙}$) fragmentation becomes inevitable and leads to the formation of Pop II stars (Omukai et al. 2005; Glover & Jappsen 2007; Wise et al. 2012; Bovino et al. 2014; Safranek-Shrader et al. 2014; Ritter et al. 2015). The critical value of the metallicity is further reduced by about two orders of magnitude in the presence of trace amounts of dust (Schneider et al. 2003; Omukai et al. 2005, 2008), and the $H_2$ formation on dust grains enhances the cooling (Omukai 2001; Cazaux & Spaans 2009; Latif et al. 2012). Numerical simulations show that low-mass star formation seems plausible because of the dust cooling occurring at high densities in minihalos with metallicities as low as a few times $10^{-5}$ solar (Tsuribe & Omukai 2006, 2008; Dopcke et al. 2013; Smith et al. 2015).

Omukai et al. (2008) explored the thermodynamic properties of a gas irradiated by a strong LW flux above $J_{crit}$ and enriched by trace amounts of metals and dust. They found from a one-zone model that in such conditions dust cooling becomes effective for $Z/Z_{⊙} \ge 5 \times 10^{-6}$ and proposed that a dense stellar cluster can form under these conditions that may later collapse into a black hole of 100–1000 $M_{⊙}$.

In this study, we explore the impact of dust and metal line cooling in massive primordial halos illuminated by a strong background LW flux and polluted by trace amounts of metals and dust. To accomplish this goal, we perform high-resolution 3D cosmological simulations for two halos of a few times $10^7 M_{⊙}$ by turning on an LW flux of strength $10^5$ in units of $J_{21}$ at $z=30$ and presume that they are pre-enriched with $Z/Z_{⊙} = 10^{-6}-10^{-4}$. Our results suggest that dust cooling occurs even for $Z/Z_{⊙} = 10^{-6}$ at densities above $10^{10} \text{ cm}^{-3}$, but the conditions are still suitable for the formation of massive objects. This work has important implications for assessing the feasibility of black hole formation models.

This article is organized as follows. In Section 2, we describe numerical methods and the chemical model employed in this work. We present our results in Section 3 and discuss the implications for black hole formation in Section 4. A summary of our main results and our conclusions are discussed in Section 5.

## 2. NUMERICAL METHODS

We have performed simulations using the open-source code Enzo version 2.4 (Bryan et al. 2014). Enzo is an adaptive mesh refinement, parallel, grid-based cosmological simulation code that has been extensively used to perform high-resolution simulations. It uses the message-passing interface to achieve scalability and portability on various platforms. The piecewise parabolic method is employed to solve hydrodynamics, the particle-mesh technique is used to solve dark matter dynamics, and the multigrid Poisson solver is used for gravity.

### 2.1. Simulation Setup

Our simulations are started with cosmological initial conditions at $z=100$ generated with the “inits” package available with Enzo. The computational domain has a size of $1 \text{ Mpc h}^{-1}$ in comoving units, and periodic boundary conditions are adopted both for gravity and hydrodynamics. We use parameters from the WMAP 7-year data to generate the initial

conditions (Jarosik et al. 2011). Our computational box is centered on the most massive halo forming at $z=15$ selected from uniform-grid, dark matter (DM)-only simulations of resolution $128^3$ particles. We rerun the simulations with a top grid resolution of $128^3$ grid cells (same number of DM particles) and add two nested refinement levels, each with a resolution of $128^3$ grid cells and $128^3$ DM particles. In all, we employ 5,767,168 particles to solve the DM dynamics, which provides us with an effective DM resolution of about 600 $M_{⊙}$. We smooth the DM particles at about a 2 pc scale to avoid spurious numerical artifacts; see Latif & Volonteri (2015) for a detailed discussion. We further add 28 dynamical refinement levels during the course of the simulations in the central 62 kpc comoving region of the box. This approach enables us to resolve the gravitational collapse down to scales of about 0.1 au. We employ a fixed Jeans resolution of 32 cells throughout the simulations to resolve turbulent eddies and to fulfill the Truelove criterion, which requires a resolution of at least four cells per Jeans length (Truelove et al. 1997; Federrath et al. 2011; Latif et al. 2013b). Our dynamical refinement criterion further includes the baryonic overdensity and particle mass resolution. Grid cells are marked for refinement if the gas density exceeds four times the cosmic mean. Similarly, cells are flagged for refinement if the DM density is above 0.0625 times $\rho_{DM} r^{\ell\alpha}$, where $\rho_{DM}$ is the dark matter density, $r=2$ is the refinement factor, $\ell$ is the refinement level, and $\alpha = -0.3$ makes the refinement super-Lagrangian.

We select two halos of $5.6 \times 10^7 M_{⊙}$ and $3.25 \times 10^7 M_{⊙}$. The properties of these halos, such as their spins, collapse redshifts, and $J_{crit}$, are listed in Table 1 of Latif & Volonteri (2015) as halo A and halo C, respectively. We presume that the halos are enriched with trace amount of metals and dust by supernova winds from a nearby star-forming galaxy. The latter also provides a strong LW flux above the critical threshold required for the quenching of $H_2$ formation. The prime aim of this study is to explore the impact of dust cooling in the presence of an intense LW flux. We therefore select $J_{21} = 10^5$, for which a halo collapses isothermally in the absence of dust cooling. At densities above $10^{-16} \text{ g cm}^{-3}$, where dust cooling becomes important, the gas cloud is already self-shielded against external radiation, and therefore the value of $J_{crit}$ does not depend on the metal content of the halo. This scenario is particularly relevant for the halos that get polluted while being irradiated by a strong LW flux, and we specifically explore here the cases with $Z/Z_{⊙} = 10^{-6}, 10^{-5}$ and $10^{-4}$.

Because of the computational constraints, we stop our simulations when a peak density of $3 \times 10^{-11} \text{ g cm}^{-3}$ is reached. Omukai et al. (2008) have shown that dust cooling becomes optically thick at densities above $\sim 10^{-10} \text{ g cm}^{-3}$ for the metallicities explored here, and thereafter collapse is expected to proceed adiabatically. The transition to the adiabatic regime during the gravitational collapse of primordial gas has been explored in hydrodynamic simulations (Inayoshi et al. 2014; Van Borm et al. 2014; Becerra et al. 2015; Latif et al. 2016). These studies found that the adiabatic phase stabilizes the gravitational collapse on small scales and breaks its self-similarity, and a protostar begins to form at this stage. For dust cooling, Omukai et al. (2008) have shown that the transition to the adiabatic phase occurs at earlier stages of the collapse (at lower densities compared to the primordial gas). They also considered the sputtering of dust grains in their model, finding that the dust temperature remains below the
vaporization temperature in the optically thin regime, and once the adiabatic phase begins, it exceeds the evaporation temperature but the gas cloud is already optically thick. Therefore, the thermal evolution remains adiabatic, similar to the metal-free gas. Further details of our model are described in Omukai et al. (2008).

In this study, we explore three cases with metallicities of $Z/Z_{\odot} = 10^{-6}$, $10^{-5}$ and $10^{-4}$ for each halo and turn on a strong LW flux of strength $J_{21} = 10^5$, well above the $J_{crit}$ found in Latif et al. (2015). Further details of our chemical network, as well as heating and cooling processes, are described in the following section.

## 2.2. Chemical Model

We employ the KROME package (Grassi et al. 2014) to solve the chemical and thermal evolution of the gas in cosmological simulations. The rate equations of H, $H^+$, $H^-$, He, $He^+$, $He^{++}$, $H_2$, $H_2^+$, $e^-$ are solved to study their evolution during the gravitational collapse. We employ a uniform-background LW flux of strength $J_{21} = 10^5$ emitted from a Pop II galaxy, which can be mimicked by a radiation temperature of $2 \times 10^4$ K (see Sugimura et al. 2014; Latif et al. 2015). Our results do not depend on the choice of the radiation spectrum as long as the flux is well above the critical threshold to dissociate molecular hydrogen (derived by Omukai et al. 2008; Sugimura et al. 2014 for the one-zone model and by Latif et al. 2015 in 3D simulations). The list of chemical reactions is given in Table 1 in the appendix of Latif et al. (2015). Our model includes $H_2$ formation, $H_2$ and $H_2^+$ photodissociation, $H^-$ photodetachment, collision-induced emission, and $H_2$ collisional dissociation. In addition to this, we add a reaction for $H_2$ formation on dust grain surfaces taken from Omukai (2000). Various cooling and heating mechanisms, such as cooling that is due to the collisional excitation, collisional ionization, radiative recombination, Bremsstrahlung radiation, and chemical cooling/heating from three-body reactions, are included in the chemical model. We also employ the $H_2$ self-shielding fitting formula given in Wolcott-Green et al. (2011).

In addition to the cooling and heating processes mentioned above, we include cooling that is due to the fine structure lines of C II and O I and cooling and heating by dust grains, which is relevant for $Z/Z_{\odot} = 10^{-6}$, $10^{-5}$, and $10^{-4}$. Our treatment for metal and dust cooling/heating processes comes from the model of Omukai et al. (2008); here we summarize its main features. We do not explicitly include the chemical reactions involving metals but instead assume that most of the oxygen and carbon are in O I and C II form because of the presence of a strong LW flux. Carbon has a lower ionization energy (11.26 eV) than hydrogen and therefore is in C II form, while O I remains in the atomic state because of the charge exchange with hydrogen atoms in a neutral medium. The dust grain composition and size distribution are assumed to be similar to the solar neighborhood, and the amount is scaled with the metallicity of the gas cloud. We assume a dust-to-gas mass ratio of 0.01 $Z/Z_{\odot}$. Dust grains in the supernova ejecta are more efficient in cooling and $H_2$ formation because of their smaller size and larger surface area (Schneider et al. 2003; Omukai et al. 2008). The dust is assumed to be in thermal equilibrium, and its temperature ($T_{dust}$) is evaluated separately from the gas temperature at each density and temperature. The following energy-balance equation is solved to compute the

$T_{dust}$ (Omukai et al. 2008):

$$4\pi \int \kappa_{a,\nu} B_{\nu}(T_{\text{dust}}) d\nu = \Lambda_{\text{gas}\to\text{dust}} + 4\pi \int \kappa_{a,\nu} B_{\nu}(T_{\text{CMB}}) d\nu. \quad (1)$$

Here $\Lambda_{gas\todust}$ is the net gas cooling rate per unit mass that is due to dust grain collisions, $\kappa_{a,\nu}$ is the dust absorption opacity, and $B_{\nu}$ is the blackbody spectrum. Also, $\Lambda_{gas\todust}$ is given as (Hollenbach & McKee 1979)

$$\Lambda_{\text{gas}\to\text{dust}} = 5.8 \times 10^{-8} n \left(\frac{T}{1000}\right)^{1/2} \times \left[1 - 0.8 \exp\left(\frac{-75K}{T}\right)\right] (T - T_{\text{dust}}) (Z/Z_{\odot}) \quad (2)$$

where $n$ is the gas number density and $T_{dust}$ is the dust temperature. We have set the temperature floor at $T_{CMB}$. Further details about the dust properties can be found in Section 2.2 of Omukai et al. (2008).

## 3. RESULTS

We here present our main results obtained both from a one-zone model and 3D cosmological simulations.

### 3.1. One-zone Model

To test the chemical model presented in the previous section, we performed a one-zone test with a LW flux of strength $J_{21} = 10^5$ (well above the $J_{crit}$ from the one-zone model with $T_{rad} = 2 \times 10^4$ K) and varied the metallicity in the range $Z/Z_{\odot} = 10^{-6}-10^{-4}$. We took an initial temperature of 160 K, a gas density of $\sim 10^{-23}$ g cm$^{-3}$, and species abundances of $2 \times 10^{-6}$ for $e^-$ and $2 \times 10^{-4}$ for $H_2$. In Figure 1, we show the temperature and the abundances of $H_2$, $H^-$, and $e^-$ as a function of density. In the presence of a strong LW flux, $H_2$ gets photodissociated, its abundance remains below $10^{-8}$, and it is unable to cool the gas. The temperature increases up to about $10^4$ K, where atomic line cooling becomes effective and cools the gas down to 8000 K up to densities of $10^{-17}$ g cm$^{-3}$. At higher densities for $Z/Z_{\odot} \ge 10^{-5}$, dust cooling comes into play and sharply cools the gas down to a few hundred kelvin. For the case of $Z/Z_{\odot} = 10^{-4}$, $H_2$ formation on dust grains becomes effective, and its abundance get boosted very quickly. The small increase in the temperature at a density of $10^{-13}$ g cm$^{-3}$ is due to the chemical heating produced by the formation of $H_2$ on grain surfaces. Dust cooling continues at high densities up to $10^{-10}$ g cm$^{-3}$, and after that the cloud becomes optically thick to dust self-absorption. Beyond this point, the thermal evolution proceeds adiabatically, and the inclusion of dust cooling breaks the self-similarity of an isothermal collapse. A protostar is expected to form at this stage with a central adiabatic core and is expected to grow rapidly by accretion from the protostellar envelope.

For $Z/Z_{\odot} = 10^{-5}$, the dust cooling becomes effective at densities of $10^{-14}$ g cm$^{-3}$, about two orders of magnitude higher than for the $Z/Z_{\odot} = 10^{-4}$ case. The temperature sharply declines down to a few hundred kelvin, and the $H_2$ fraction gets boosted. The cloud becomes optically thick to dust cooling, similar to the $Z/Z_{\odot} = 10^{-4}$ case at $10^{-10}$ g cm$^{-3}$, and the thermal evolution follows the adiabatic equation of state. At an even lower metallicity of $Z/Z_{\odot} = 10^{-6}$, the dust cooling becomes negligible, and the
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Temperature, H₂, H⁻, and e⁻ fractions plotted against the density from the one-zone test. The green, blue, and red lines represent Z/Z☉ = 10⁻⁶, 10⁻⁵, and 10⁻⁴, respectively.</figcaption></figure>

thermal evolution follows the atomic cooling track. The H₂ fraction is increased up to about a few times 10⁻⁴. The abundances of e⁻ and H⁻ remain very low and are similar at high densities for the Z/Z☉ = 10⁻⁵ and Z/Z☉ = 10⁻⁴ cases. The degree of ionization is higher for Z/Z☉ = 10⁻⁶. The impact of metal line cooling such as CII and OI is negligible for such metallicities. Overall, our results are in good agreement with Omukai et al. (2008).

## 3.2. 3D Simulations

In total, we have performed six cosmological simulations for two halos of 5.6 × 10⁷ M☉ (H1) and 3.25 × 10⁷ M☉ (H2) for metallicities of Z/Z☉ = 10⁻⁴, 10⁻⁵ and 10⁻⁶. The collapse redshifts of the halos are 10.4 and 10.8, respectively, and are independent of the dust content. We also performed one additional simulation with zero metallicity for the sake of comparison for H1. The gravitational collapse of the halos was followed from 1 Mpc down to sub-astronomical unit scales with the adaptive mesh refinement approach. We assumed that the halos are irradiated by a strong LW flux of strength 10⁵ in terms of J₂₁, well above the Jcrit found in Latif et al. (2015), and the LW flux as well as the fixed metallicities of the above-mentioned amount are turned on at z = 30.

### 3.2.1. Thermodynamical and Dynamical Properties

The density-temperature phase diagrams of both halos for Z/Z☉ = 10⁻⁴, 10⁻⁵, and 10⁻⁶ are shown in Figure 2. In the

presence of a strong LW flux that dissociates H₂ molecules, the halos are unable to collapse via molecular hydrogen cooling and continue to grow via merging and accretion until they reach the atomic cooling limit. Both halos have recently gone through a major merger, and the merger history of these halos is described in Latif & Volonteri (2015). At densities of about ≥10⁻²⁴ g cm⁻³, atomic line cooling becomes effective and cools the gas to about 8000 K. The collapse proceeds isothermally mainly via atomic hydrogen cooling up to densities of 10⁻¹⁶ g cm⁻³, and the role of metal line cooling remains negligible. At densities >10⁻¹⁶ g cm⁻³, dust cooling comes into play and brings the gas temperature down to a few hundred kelvin depending on the amount of dust. For Z/Z☉ = 10⁻⁶, the collapse remains isothermal up to densities of about 10⁻¹² g cm⁻³. At higher densities, the dust cools the gas down to a temperature of ~1000 K for halo 2, while both warm and cold phases coexist for halo 1 at densities ≥10⁻¹² g cm⁻³, which is due to the local variation in the gas density and collapse velocity. In the intermediate case of Z/Z☉ = 10⁻⁵, the dust cooling becomes important earlier at densities of ~10⁻¹⁴ g cm⁻³ for both halos and decreases the gas temperature down to about 700 K.

For Z/Z☉ = 10⁻⁴, the dust grain cooling becomes strong at densities of ~10⁻¹⁶ g cm⁻³ and sharply brings the gas temperature down to a few hundred kelvin. The temperature stalls at about 1000 K at densities in the range ~10⁻¹⁴–10⁻¹² g cm⁻³ because of the chemical heating from
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 2.** Density–temperature phase diagram for $Z/Z_⊙ = 10^{-6}$, $10^{-5}$, and $10^{-4}$ from top to bottom, respectively. The left panel represents halo 1 and the right panel halo 2. The gas is initially heated up to a few times $10^4$ K and collapses isothermally up to densities of $10^{-16}$ g cm$^{-3}$ in the presence of a strong LW flux, and the dust cooling becomes effective at densities $\ge 10^{-16}$ g cm$^{-3}$, depending on the amount of metallicity. Color bars show the amount of gas in solar masses for a given density.</figcaption></figure>

the H₂ formation on grain surfaces, and again the dust cooling takes over at higher densities. These results are in agreement with the one-zone test presented in the previous section and also with Omukai et al. (2008). Small differences in the thermal

evolution are observed for both halos; the cooling at high densities is generally more effective for halo 2 in comparison with halo 1. We show the spherically averaged profiles of density, H₂ fraction, and temperature in Figure 3. The collapse
<figure><img src="image_9.png" /><figcaption>**Figure 3.** Radially averaged and spherically binned profiles of various quantities for halo 1 (H1) and halo 2 (H2) are shown here. The green, blue, and red lines represent $Z/Z_\odot = 10^{-4}$, $10^{-5}$, and $10^{-6}$, respectively. The solid lines represent halo 1, and the dashed lines show halo 2. The vertical black dashed line corresponds to the Jeans length. An isothermal case for halo 1 with zero metallicity is also plotted as a reference and is represented by the black solid line.</figcaption></figure>
remains isothermal down to the scales of 10,000 au irrespective of the metallicity explored in this work because of the strong LW flux. Below this scale, dust cooling becomes important and remains confined to the central 10 au for $Z/Z_⊙ = 10^{-6}$ and the inner ~100 au for $Z/Z_⊙ = 10^{-5}$, and it extends up to 3000 au for $Z/Z_⊙ = 10^{-4}$. The central temperature is about 300, 700, and 1000 K for $Z/Z_⊙ = 10^{-4} - 10^{-6}$, respectively. Similarly, the H₂ fraction remains low at scales above $10^5$ au and exceeds $10^{-3}$ (above which H₂ cooling becomes important) within the central part of the halo because of its formation on dust grains. For $Z/Z_⊙ = 10^{-4}$, the H₂ abundance gets boosted within the central few 1000 au and for $Z/Z_⊙ = 10^{-5}$ only in the central few hundred astronomical units, and for $Z/Z_⊙ = 10^{-6}$ it remains confined to the central 10 au.

The maximum density reached in our simulations is about $3 \times 10^{-11}$ g cm$^{-3}$, and the density increases with $R^{-1.8}$ above ~100 au, close to the expected profile from an isothermal collapse (i.e., $R^{-2.0}$). The density profile is shallower for $Z/Z_⊙ = 10^{-4}$ in the central 1000 au because of efficient dust cooling, while for the other cases it almost follows the isothermal profile. Small bumps in the profile indicate the formation of additional clumps. The turbulent velocity is about $20$ km s$^{-1}$ close to the viral radius and increases toward smaller scales up to about $30$ km s$^{-1}$. This increase in the turbulent velocity is due to the infall of gas toward the center and is higher for dust cooling cases compared to the metal-free case. To quantify whether a disk-like structure may form, we computed the ratio of the velocity dispersion (i.e., $\sqrt{c_s^2 + v_{\text{turb}}^2}$ where $c_s$ is the sound speed and $v_{\text{turb}}$ is the turbulent velocity) to the rotational velocity, corresponding to the vertical turbulent support within the disk. The latter provides an estimate for the ratio of vertical scale height to the disk radius (H/R). We found that this ratio remains about 1 for all cases, except for halo 2 with $Z/Z_⊙ = 10^{-5}$ where $V_{\text{rot}}/\sigma \sim 0.3$. This suggests that the rotational velocity is almost comparable to the velocity dispersion, and a very thick disk-like structure with an H/R of about 1 is formed in all cases. To further estimate the rotational support against gravity, we calculated the ratio of rotational velocity to the Keplerian velocity ($V_{\text{rot}}/V_{\text{Kep}}$). The ratio of $V_{\text{rot}}/V_{\text{Kep}}$ is shown in Figure 4 and is about 0.6 for all cases. It gets enhanced between 60 and 1000 au for dust cooling cases. A peak in $V_{\text{rot}}/V_{\text{Kep}}$ for the $Z/Z_⊙ = 10^{-5}$ case comes from the density structure inside the halo. Overall, the estimates of $V_{\text{rot}}/V_{\text{Kep}}$ suggest that the rotational support gets enhanced earlier for higher metallicity cases. However, enhanced rotational support can delay the runaway collapse but does not stop it. Moreover, the rotational support is only important in the envelope of the halo and decreases down within the central core.

### 3.2.2. Mass Inflow Rates

Our estimates for the mass inflow rates ($4\pi R^2 \rho v_{\text{rad}}$) are shown in Figure 3. Mass inflow rates of $0.1-1$ M$_{\odot}$ yr$^{-1}$ are observed for $Z/Z_⊙ \le 10^{-5}$ down to about 100 au, while for $Z/Z_⊙ = 10^{-4}$ the mass inflow rate starts to decline around $10^5$ au for halo 2. For comparison, we also show a zero-metallicity case where cooling is mainly due to atomic line cooling. We found that, apart from small differences, the inflow rates for $Z/Z_⊙ \le 10^{-5}$ almost follow the zero-metallicity case. Below 100 au the mass inflow rates decline in all cases. This comes from the fact that during the runaway collapse a core-

<figure><img src="image_10.png" /><figcaption>**Figure 4.** The ratio of rotational to Keplerian velocity for halo 1 (H1) and halo 2 (H2). The green, blue, and red lines represent $Z/Z_⊙ = 10^{-4}$, $10^{-5}$, and $10^{-6}$, respectively. The solid lines represent halo 1, and the dashed lines show halo 2. An isothermal case for halo 1 with zero metallicity is also plotted as a reference and is represented by the black solid line. The vertical black dashed line corresponds to the Jeans length.</figcaption></figure>

envelope structure develops, where the core has flat density and the envelope has an ~$R^{-2}$ profile with a core-envelope boundary around 10–100 au. The radial velocity is high inside the envelope and decreases down in the core, and consequently a decline in the instantaneous inflow rate occurs toward the center. In the central core, the thermal pressure balances the gravity as the core is not yet gravitationally unstable, which results in a decline of the inflow rate. However, this may change at later stages when the central core further collapses. The fragmentation observed in the $Z/Z_⊙ = 10^{-4}$ cases may also partly contribute to the decline in the inflow rates. We also note that the mass inflow rate scales with the speed of sound cubed ($\dot{M} \propto c_s^3$) and therefore decreases with the gas temperature. The dust cooling reduces the gas temperature, resulting in lower inflow rates. The accreting matter inside the envelope will collapse into the protostar at later times, and high inflow rates are expected to be maintained. The increase in the mass inflow rates was observed in our previous simulations for a metal-free case (Latif et al. 2013a) as more mass collapsed into the core. We expect a similar trend for other cases until the feedback from the central star blows away the gas or rotational support inhibits the mass accretion.

A similar trend is observed in the radial infall velocity, which peaks around $15$ km s$^{-1}$ and brings large inflows into the center of the halo. The radial velocity temporarily becomes positive in a few places because of the presence of clumps, and the positive sign then describes the flow toward the other clump. This is particularly prominent in halo 2 for $Z/Z_⊙ = 10^{-4}$ and causes a decline in the mass inflow rate. Repeated peaks in the radial velocity indicate the onset of gravitational instabilities, which help in the transfer of angular momentum. The enclosed gas mass in the central 30 pc of the halo is a few times $10^6$ M$_{\odot}$. The mass profile increases with ~$R^2$ in the central 100 au and at larger scales linearly increases with radius. Some differences in the mass profile are observed between the metal-free and metal-poor cases. The enclosed mass is a factor of a few lower for $Z/Z_⊙ = 10^{-4}$ in comparison with the metal-free case and is probably a consequence of
<figure><img src="image_11.png" /><figcaption>**Figure 5.** Average gas density along the *x* axis for the central 1 pc of a halo. Each row represents a halo (halo 1 (H1) on top and halo 2 (H2) at the bottom), and each column represents a metallicity (increasing from left to right): H1Ze-6, H1Ze-5, and H1Ze-4 represent $Z/Z_\odot = 10^{-6}$, $10^{-5}$, and $10^{-4}$ for halo 1, respectively. Similarly, H2Ze-6, H2Ze-5, and H2Ze-4 represent $Z/Z_\odot = 10^{-6}$, $10^{-5}$, and $10^{-4}$ for halo 2.</figcaption></figure>

differences in the thermal structure as well as enhanced rotation for the $Z/Z_\odot = 10^{-4}$ cases.

Overall, the mass inflow rates at the end of our simulations seem to be sufficient for forming a supermassive star, at least for $Z/Z_\odot \le 10^{-5}$, while for higher metallicities enhanced rotation as well as fragmentation may limit accretion onto the central object. How long such accretion rates can be maintained depends on whether efficient fragmentation takes place or not. In the case of efficient fragmentation, radiation feedback from in situ star formation may shut down accretion onto the central star, and a supermassive star cannot form. Even in the absence of fragmentation, an enhanced rotation can limit mass accretion, as found for $Z/Z_\odot \le 10^{-4}$ in both halos. We discuss these possibilities in the following section.

### 3.2.3. Fragmentation

To quantify the fragmentation, we use the YT clump finder (Turk et al. 2011) to identify clumps. Our clump-finding algorithm locates topologically disconnected structures. The clumps are considered to be bound if the sum of kinetic and thermal energy is less than the potential energy; see Smith et al. (2009). The density structure within the central 1 pc of the halo is shown in Figure 5 for all cases. There is not much substructure found at this scale because the collapse is almost isothermal and dust cooling starts to become important around a few 1000 au. Zooming into the central 4000 au (see Figure 6), we observe significant differences in the morphology of the halos. Particularly for $Z/Z_\odot = 10^{-4}$, the structure is more elongated and collapses into a filament because of the more efficient cooling compared with other cases where central gas

clouds are more spherical. For $Z/Z_\odot = 10^{-4}$, two well-separated gravitationally bound clumps of about a solar mass are formed in halo 2, while in halo 1 the central clump is of similar mass and gravitationally bound, but the second clump is gravitationally unbound and has a subsolar mass. In contrast to this, the clouds are more spherical (although not completely) for $Z/Z_\odot = 10^{-6}$, very similar to a zero-metal case.

For $Z/Z_\odot = 10^{-5}$, there is more substructure in the central part of halo 1, with the possibility of multiple clump formation, while halo 2 looks similar to the metal-free case. Overall the clouds are more elongated compared to the $Z/Z_\odot = 10^{-6}$ case. We also ran a clump finder for these cases and found that all clumps in these cases are gravitationally unbound and mostly with subsolar masses, as shown in Figure 8. Figure 7 shows that only small pockets of gas are cooled by dust cooling. The $Z/Z_\odot = 10^{-6}$ case looks very similar to the metal-free case. The temperature of the dense clumps of gas is about a thousand kelvin for $Z/Z_\odot = 10^{-5}$, and it declines down to a few hundred kelvin for $Z/Z_\odot = 10^{-4}$. Similar differences in the density structure of a gas cloud are observed by Peters et al. (2012) employing an equation of state with polytropic index ($\gamma<1$).

We estimated the maximum scale of thermal instability (i.e., $l_{th} = c_s t_{cool}$) following the criterion given in Inoue & Omukai (2015) and found that it is comparable to the Jeans length. It shows that the resolution in our simulations is enough to resolve the thermal instability. Dust cooling triggers the Jeans instability, which results in the formation of multiple clumps. In a self-gravitating cloud, turbulence can also induce fragmentation by locally compressing the gas, as observed in an isothermal collapse; see Latif et al. (2013a). Overall, central
<figure><img src="image_12.png" /><figcaption>**Figure 6.** Average gas density along the *x* axis for the central 4000 au of a halo. Each row represents a halo (halo 1 on top and halo 2 at the bottom) and each column represents a metallicity (increasing from left to right).</figcaption></figure>

clumps forming for $Z/Z_\odot = 10^{-4}$ are gravitationally bound, and for low metallicities all of the clumps are gravitationally unbound. We caution the reader that gravitationally unbound clumps at the initial stages of the collapse do not reflect the true amount of fragmentation. To further estimate the fate of these clumps, we followed the collapse of halo 2 with $Z/Z_\odot = 10^{-5}$ to densities of $10^{-8} \text{ g cm}^{-3}$ by employing five additional refinement levels. In this case, clumps eventually collapse into the center of the halo. However, the possibility of fragmentation at later stages of the collapse cannot be ruled out, and it may vary from halo to halo. Evolving these simulations further in time becomes computationally extremely expensive because of the large dynamical range covered in these simulations, which results in very short time steps. In the following section, we discuss the possible implications of efficient fragmentation if it occurs at later times.

### 3.3. Fraction of Metal-enriched Halos

We estimate the fraction of metal-polluted halos with $Z/Z_\odot \le 10^{-4}$ from cosmological hydrodynamic simulations as in Habouzit et al. (2016a). The computational volume has a size of 10 Mpc, and the simulations include recipes for star formation and supernova feedback; a detailed discussion is provided in Habouzit et al. (2016a). We identified the halos with mass in the range $2 \times 10^7-10^8 M_\odot$ using the HaloMaker code (Tweed et al. 2009) at least resolved by 100 hundred DM particles, and we computed their mean metallicity. Our estimates for the fraction of halos with $Z/Z_\odot \le 10^{-4}$ are shown in Figure 9. The number of metal-polluted halos increases with decreasing redshift, as expected, and the fraction

of halos polluted with $Z/Z_\odot \le 10^{-5}$ is about a factor of 1.5 higher in comparison with metal-free halos.

The halos polluted by a trace amount of metals can be irradiated by a strong LW flux if they form in the close vicinity of a star-forming galaxy. Our simulations do not self-consistently take into account the local variations in the LW flux required for the direct-collapse sites, and the simulation volume is not large enough to capture rare sources that provide a strong LW flux of $10^5$ in terms of $J_{21}$. Assuming that the metallicity distribution is, however, uncorrelated with fluctuations in the LW background, if a direct collapse is feasible for metallicities up to $10^{-5}$, the resulting abundance of the halos is increased by a factor of 1.5 compared to the purely primordial case.

### 4. IMPLICATIONS FOR BLACK HOLE FORMATION

Stellar evolution calculations suggest that the main requirement for the formation of DCBHs are mass inflow rates of $\ge 0.1 M_\odot \text{ yr}^{-1}$. For such accretion rates, the stellar evolution differs from normal stars because the stellar radius increases with mass and no strong UV feedback is produced by the protostar. Under these conditions, a supermassive star of about $10^5 M_\odot$ can be formed that may later collapse into a black hole while retaining most of its mass (Hosokawa et al. 2012, 2013; Schleicher et al. 2013).

In this study, we explored the impact of dust cooling occurring at high densities in the presence of a strong LW flux for two halos of a few times $10^7 M_\odot$ and metallicities of $Z/Z_\odot = 10^{-4}-10^{-6}$. We resolved the collapse down to sub-astronomical unit scales and densities of about a few times $10^{-11} \text{ g cm}^{-3}$. The presence of large mass inflow rates of
<figure><img src="image_13.png" /><figcaption>**Figure 7.** Average temperature weighted by the gas density along the *x* axis for the central 4000 au of a halo. Each row represents a halo (halo 1 on top and halo 2 at the bottom), and each column represents a metallicity (increasing from left to right).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 8.** The ratio of clump mass to the Jeans mass is plotted against the clump masses for various metallicities. The left panel is for halo 1 and the right panel halo 2. The red symbols represent the clumps forming for $Z/Z_\odot = 10^{-6}$, the blue for $Z/Z_\odot = 10^{-5}$, and the green for $Z/Z_\odot = 10^{-4}$. The clumps for $Z/Z_\odot \le 10^{-5}$ are gravitationally unbound, while the central clumps for the $Z/Z_\odot = 10^{-4}$ cases are gravitationally bound and more massive. In general, the ratio of clump mass to the Jeans mass increases with clump masses.</figcaption></figure>

$\ge 0.1 \ M_\odot \ yr^{-1}$ suggests that trace amounts of metals and dust do not inhibit the formation of massive objects. At least for $Z/Z_\odot \le 10^{-5}$, the mass-accretion rates are comparable with an atomic cooling case, and no strong fragmentation is observed. It should be noted that the conditions found here are different from minihalos cooled by dust cooling in the absence of a strong LW background. Both the mass inflow rates and available gas mass are about two orders of magnitude higher compared to minihalos. Because of the extremely large spatial

range covered in our simulations, the time step becomes very short and does not allow us to evolve simulations long enough to assess fragmentation and inflow rates at later times. In the following, we consider two possible outcomes while assuming that efficient fragmentation occurs at later times.

In the first case, we consider that efficient fragmentation occurs at later times but clumps move inward because of the short migration timescale and merge with the central clump (Inayoshi & Haiman 2014; Latif & Schleicher 2015a, 2015b;
<figure><img src="image_15.png" /><figcaption>**Figure 9.** Fraction of halos with metallicities below a given value (see legend) as a function of redshift. In this figure, we only show the fraction of halos with masses in the range $2 \times 10^7-10^8 M_\odot$ in a computational box of size 10 Mpc for $Z/Z_\odot \le 10^{-4}$. The fraction of halos with $Z/Z_\odot \le 10^{-5}$ is about 1.5 times higher than for the metal-free halos.</figcaption></figure>

Schleicher et al. 2016). In the presence of the large inflow rates of $0.1 M_\odot \text{yr}^{-1}$ found here, viscous heating may become important in the interior of the disk: it will heat the gas up to a few thousand kelvin. Therefore, viscous heating may suppress fragmentation by increasing the thermal Jeans mass and may help in the formation of massive objects (Latif & Schleicher 2015b; Schleicher et al. 2016). Moreover, for such large inflow rates, the impact of UV feedback from stars is expected to be quite weak. Even if some of the clumps survive or get ejected via three-body processes and form low-mass stars, a massive central object is still expected to form because of the clumpy mode of accretion, and the stabilization of the disk is due to the viscous heating. We argue that, in such a scenario, massive seed black holes of $\sim 10^4 M_\odot$ can form for $Z/Z_\odot \le 10^{-5}$.

For the second case, we assume that efficient fragmentation takes place that is due to the dust cooling at high densities. Then the clumps are unable to migrate inward, and a stellar cluster forms. A dense stellar cluster is expected to form in such a scenario because dust cooling becomes effective only in the central few 100–1000 au surrounded by hot gas with a temperature around 8000 K. The minimum Jeans mass at high densities where dust cooling operates in our simulations is about $>0.1 M_\odot$, and the enclosed mass within a radius of a few hundred astronomical units is about 1000 $M_\odot$. As a consequence of gravitational collapse, a higher gas mass will collapse in the center over time and may reach up to $10^6 M_\odot$, which is about a few percent of the halo mass, as often observed in numerical simulations (Johnson et al. 2009; Latif et al. 2011). For a star-formation efficiency of about 10%, a stellar cluster of about $10^5 M_\odot$ may form. The expected radius of such a stellar cluster is about 1000 au. Various studies suggest that such dense clusters may collapse into a massive black hole either via relativistic instabilities or stellar dynamical processes (Baumgarte & Shapiro 1999; Portegies Zwart et al. 1999; Omukai et al. 2008; Devecchi & Volonteri 2009; Devecchi et al. 2012). We estimate the mass of a seed black

hole forming from the core collapse of a cluster in the following way (Portegies Zwart & McMillan 2002):

$$M_{\text{BH}} = m_{*} + 4 \times 10^{-3} f_{c} M_{c0} \gamma \ln \Lambda_{c}, \quad (3)$$

where $m_*$ is the mass of a massive star in the cluster that initiates runaway growth, $f_c$ is the effective fraction of binaries formed dynamically, $M_{c0}$ is the mass of the cluster at its birth, $\ln\Lambda_c$ is the Coulomb logarithm, and $\gamma \sim 1$ is the ratio of various timescales; see Portegies Zwart & McMillan (2002). For a star cluster of $M_{c0} = 10^4 M_\odot$, $m_* = 100 M_\odot$, $f_c = 0.2$, $\ln\Lambda_c = 10$, we get $M_{\text{BH}} \sim 180 M_\odot$. For $M_{c0} = 10^5 M_\odot$, the expected black hole mass is $\sim 900 M_\odot$. Therefore, depending on the initial mass of a cluster and the central star, massive black hole seeds of up to a thousand solar masses can be formed in metal-poor halos illuminated by the strong LW flux.

## 5. DISCUSSION AND CONCLUSIONS

We have performed high-resolution cosmological simulations to study the impact of trace amounts of metals and dust cooling in massive primordial halos irradiated by the strong LW flux. To accomplish this goal, we selected two massive halos of a few times $10^7 M_\odot$ at $z > 10$ with metallicities of $Z/Z_\odot = 10^{-4}-10^{-6}$ and turned on a strong LW radiation of strength $10^5$ in terms of $J_{21}$. Our simulations cover a large dynamical range by resolving the collapse starting from cosmological scales down to scales of sub–astronomical units. To take into account the effect of dust and metal line cooling, we extended our previous chemical model and included H$_2$ formation on dust grains, dust grain cooling and heating, and C II and O I metal line cooling following Omukai et al. (2008).

Our results show that, even in the presence of a trace amount of metals and dust, the collapse proceeds isothermally with temperatures around 8000 K up to densities of about $10^{-16} \text{ g cm}^{-3}$. Dust cooling becomes effective at densities in the range $10^{-16}-10^{-12} \text{ g cm}^{-3}$ and brings the gas temperature down to 100–1000 K for $Z/Z_\odot = 10^{-4}-10^{-6}$, respectively. As expected, dust cooling is more efficient and occurs at earlier stages of the collapse for higher metallicities in comparison with a lower metallicity case. In contrast to the isothermal case, H$_2$ formation on dust grains takes place in the core of the halo, and its fraction gets boosted. For the metal-poor cases studied here, metal line cooling is not important but could be effective at metallicities higher than explored here.

We found that large inflow rates of $0.1 M_\odot \text{yr}^{-1}$ are available for $Z/Z_\odot \le 10^{-5}$, while for the higher metallicity case the inflow rates start to decline earlier. The decline in accretion rate for higher metallicity cases is due to the lower gas temperature in the presence of dust cooling and due to the fragmentation. No gravitationally bound clumps are found for $Z/Z_\odot \le 10^{-5}$ by the end of our simulations, while for $Z/Z_\odot = 10^{-4}$ the clumps are gravitationally bound and more massive. The morphology of the central gas cloud is significantly different in $Z/Z_\odot = 10^{-4}$ because the cloud is more elongated compared to other cases, which is probably due to the more efficient dust cooling. Larson (1969) and Penston (1969) found a self-similar solution under isothermal conditions, and later Yahil (1983) showed that it changes with the equation of state by including a dependence of the collapse time on the equation-of-state parameter $\gamma = 1 + \frac{d \log \rho}{d \log T}$. The
latter suggests that the collapse dynamics will be changed by adding dust cooling. Overall, no strong fragmentation is observed for the $Z/Z_⊙ \le 10^{-5}$ cases, although the possibility of fragmentation at later stages of the collapse cannot be ruled out. Because of computational constraints, we were unable to evolve simulations for longer times and therefore cannot make strong statements about fragmentation occurring at later times.

The presence of large inflow rates suggests that a massive central star is expected to form for $Z/Z_⊙ \le 10^{-5}$ where dust cooling is confined to the central 10–100 au depending on the metallicity. Even if the cloud fragments at later stages of the collapse, the clumps are expected to migrate inward and merge in the center, as proposed by theoretical works (Inayoshi & Haiman 2014; Latif & Schleicher 2015b). Viscous heating, which becomes particularly important in the presence of large inflow rates and rapid rotation, may further help in suppressing the fragmentation by heating the gas and evaporating the dust grains (Latif & Schleicher 2015b; Schleicher et al. 2016). For $Z/Z_⊙ = 10^{-4}$, efficient fragmentation is expected to occur that may lead to the formation of a stellar cluster.

The strength of the LW flux required to keep the collapse completely isothermal in metal-free halos is about a few times $10^4$ (Latif et al. 2015), and in the present work we employed an even higher value of $10^5$ in units of $J_{21}$. The expected number density of DCBHs forming in primordial halos for $J_{crit} \ge 10^4$ is a few orders of magnitude lower than the observed abundance of quasars at $z=6$ (i.e., 1 per Gpc$^3$; Willott et al. 2010), extrapolating the results of Habouzit et al. (2016b). Even relaxing the constraint of primordial gas and allowing the DCBHs to form in halos with $Z/Z_⊙ \le 10^{-5}$ in the presence of a strong LW flux only increases their abundance by a factor of two. Even under these conditions, the formation sites for the DCBHs remain rare.

For metallicities higher than explored here, cooling due to the fine structure lines of C II and O I becomes important at low densities about $10^{-20} \text{ g cm}^{-3}$ even in the presence of a strong LW flux and cools the gas down to a few 10 K. In such cases, the thermal instability induces fragmentation and determines the mass spectrum of the resulting star cluster (Inoue & Omukai 2015).

This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 656428. The research leading to these results has also received funding from the European Research Council under the European Community's Seventh Framework Programme (FP7/2007-2013 Grant Agreement No. 614199, project "BLACK"). KO acknowledges the Grants-in-aid from the Ministry of Education, Culture, Sports, Science, and Technology (MEXT) of Japan (25287040). This work was granted access to the HPC resources of TGCC under the allocation x2015046955 made by GENCI. The simulation results are analyzed using the visualization toolkit for astrophysical data YT (Turk et al. 2011).

## REFERENCES

* Abel, T., Bryan, G. L., & Norman, M. L. 2002, *Sci*, 295, 93

* Agarwal, B., & Khochfar, S. 2015, *MNRAS*, 446, 160

* Agarwal, B., Khochfar, S., Johnson, J. L., et al. 2012, *MNRAS*, 425, 2854

* Baumgarte, T. W., & Shapiro, S. L. 1999, *ApJ*, 526, 941

* Becerra, F., Greif, T. H., Springel, V., & Hernquist, L. E. 2015, *MNRAS*, 446, 2380

* Begelman, M. C. 2010, *MNRAS*, 402, 673

* Begelman, M. C., & Shlosman, I. 2009, *ApJL*, 702, L5

* Begelman, M. C., Volonteri, M., & Rees, M. J. 2006, *MNRAS*, 370, 289

* Bovino, S., Grassi, T., Schleicher, D. R. G., & Latif, M. A. 2014, *ApJL*, 790, L35

* Bromm, V., & Loeb, A. 2003, *ApJ*, 596, 34

* Bryan, G. L., Norman, M. I., ÓShea, B. W., et al. 2014, *ApJS*, 211, 19

* Cazaux, S., & Spaans, M. 2009, *A&A*, 496, 365

* Choi, J.-H., Shlosman, I., & Begelman, M. C. 2015, *MNRAS*, 450, 4411

* Davies, M. B., Miller, M. C., & Bellovary, J. M. 2011, *ApJL*, 740, L42

* Devecchi, B., & Volonteri, M. 2009, *ApJ*, 694, 302

* Devecchi, B., Volonteri, M., Rossi, E. M., Colpi, M., & Portegies Zwart, S. 2012, *MNRAS*, 421, 1465

* Dijkstra, M., Ferrara, A., & Mesinger, A. 2014, *MNRAS*, 442, 2036

* Dijkstra, M., Haiman, Z., Mesinger, A., & Wyithe, J. S. B. 2008, *MNRAS*, 391, 1961

* Dopcke, G., Glover, S. C. O., Clark, P. C., & Klessen, R. S. 2013, *ApJ*, 766, 103

* Fan, X., Strauss, M. A., Richards, G. T., et al. 2006, *AJ*, 131, 1203

* Federrath, C., Sur, S., Schleicher, D. R. G., Banerjee, R., & Klessen, R. S. 2011, *ApJ*, 731, 62

* Ferrara, A., Salvadori, S., Yue, B., & Schleicher, D. 2014, *MNRAS*, 443, 2410

* Glover, S. C. O., & Jappsen, A.-K. 2007, *ApJ*, 666, 1

* Grassi, T., Bovino, S., Schleicher, D. R. G., et al. 2014, *MNRAS*, 439, 2386

* Habouzit, M., Volonteri, M., Latif, M., Dubois, Y., & Peirani, S. 2016a, *MNRAS*, submitted (arXiv:1601.00557)

* Habouzit, M., Volonteri, M., Latif, M., et al. 2016b, *MNRAS*, 456, 1901

* Haiman, Z. 2012, arXiv:1203.6075

* Hirano, S., Hosokawa, T., Yoshida, N., et al. 2014, *ApJ*, 781, 60

* Hollenbach, D., & McKee, C. F. 1979, *ApJS*, 41, 555

* Hosokawa, T., Omukai, K., & Yorke, H. W. 2012, *ApJ*, 756, 93

* Hosokawa, T., Yorke, H. W., Inayoshi, K., Omukai, K., & Yoshida, N. 2013, *ApJ*, 778, 178

* Inayoshi, K., & Haiman, Z. 2014, *MNRAS*, 445, 1549

* Inayoshi, K., Omukai, K., & Tasker, E. 2014, *MNRAS*, 445, L109

* Inoue, T., & Omukai, K. 2015, *ApJ*, 805, 73

* Jarosik, N., Bennett, C. L., Dunkley, J., et al. 2011, *ApJS*, 192, 14

* Johnson, J. L., Dalla Vecchia, C., & Khochfar, S. 2013, *MNRAS*, 428, 1857

* Johnson, J. L., Greif, T. H., Bromm, V., Klessen, R. S., & Ippolito, J. 2009, *MNRAS*, 399, 37

* Katz, H., Sijacki, D., & Haehnelt, M. G. 2015, *MNRAS*, 451, 2352

* Larson, R. B. 1969, *MNRAS*, 145, 271

* Latif, M. A., Bovino, S., Grassi, T., Schleicher, D. R. G., & Spaans, M. 2015, *MNRAS*, 446, 3163

* Latif, M. A., Bovino, S., Van Borm, C., et al. 2014, *MNRAS*, 443, 1979

* Latif, M. A., & Schleicher, D. R. G. 2015a, *MNRAS*, 449, 77

* Latif, M. A., & Schleicher, D. R. G. 2015b, *A&A*, 578, A118

* Latif, M. A., Schleicher, D. R. G., & Hartwig, T. 2016, *MNRAS*, 458, 233

* Latif, M. A., Schleicher, D. R. G., Schmidt, W., & Niemeyer, J. 2013a, *MNRAS*, 433, 1607

* Latif, M. A., Schleicher, D. R. G., Schmidt, W., & Niemeyer, J. 2013b, *MNRAS*, 430, 588

* Latif, M. A., Schleicher, D. R. G., Schmidt, W., & Niemeyer, J. 2013c, *ApJL*, 772, L3

* Latif, M. A., Schleicher, D. R. G., Schmidt, W., & Niemeyer, J. 2013d, *MNRAS*, 436, 2989

* Latif, M. A., Schleicher, D. R. G., & Spaans, M. 2012, *A&A*, 540, A101

* Latif, M. A., & Volonteri, M. 2015, *MNRAS*, 452, 1026

* Latif, M. A., Zaroubi, S., & Spaans, M. 2011, *MNRAS*, 411, 1659

* Loeb, A., & Rasio, F. A. 1994, *ApJ*, 432, 52

* Lupi, A., Colpi, M., Devecchi, B., Galanti, G., & Volonteri, M. 2014, *MNRAS*, 442, 3616

* Mortlock, D. J., Warren, S. J., B. P., Venemans, et al. 2011, *Natur*, 474, 616

* Omukai, K. 2000, *ApJ*, 534, 809

* Omukai, K. 2001, *ApJ*, 546, 635

* Omukai, K., Schneider, R., & Haiman, Z. 2008, *ApJ*, 686, 801

* Omukai, K., Tsuribe, T., Schneider, R., & Ferrara, A. 2005, *ApJ*, 626, 627

* Penston, M. V. 1969, *MNRAS*, 144, 425

* Peters, T., Schleicher, D. R. G., Klessen, R. S., et al. 2012, *ApJL*, 760, L28

* Portegies Zwart, S. F., Makino, J., McMillan, S. L. W., & Hut, P. 1999, *A&A*, 348, 117

* Portegies Zwart, S. F., & McMillan, S. L. W. 2002, *ApJ*, 576, 899

* Regan, J. A., Johansson, P. H., & Haehnelt, M. G. 2014a, *MNRAS*, 439, 1160

* Regan, J. A., Johansson, P. H., & Wise, J. H. 2014b, *ApJ*, 795, 137

* Ritter, J. S., Sluder, A., Safranek-Shrader, C., Milosavljević, M., & Bromm, V. 2015, *MNRAS*, 451, 1190
Safranek-Shrader, C., Milosavljević, M., & Bromm, V. 2014, MNRAS, 438, 1669

Sakurai, Y., Hosokawa, T., Yoshida, N., & Yorke, H. W. 2015, MNRAS, 452, 755

Schleicher, D. R. G., Bovino, S., Latif, M. A., Ferrara, A., & Grassi, T. 2016, A&A, 585, A11

Schleicher, D. R. G., Palla, F., Ferrara, A., Galli, D., & Latif, M. 2013, A&A, 558, A59

Schleicher, D. R. G., Spaans, M., & Glover, S. C. O. 2010, ApJL, 712, L69

Schneider, R., Ferrara, A., Salvaterra, R., Omukai, K., & Bromm, V. 2003, Natur, 422, 869

Shang, C., Bryan, G. L., & Haiman, Z. 2010, MNRAS, 402, 1249

Shlosman, I., Choi, J.-H., Begelman, M. C., & Nagamine, K. 2016, MNRAS, 456, 500

Smith, B. D., Turk, M. J., Sigurdsson, S., O'Shea, B. W., & Norman, M. L. 2009, ApJ, 691, 441

Smith, B. D., Wise, J. H., O'Shea, B. W., Norman, M. L., & Khochfar, S. 2015, MNRAS, 452, 2822

Spaans, M., & Silk, J. 2006, ApJ, 652, 902

Sugimura, K., Omukai, K., & Inoue, A. K. 2014, MNRAS, 445, 544

Truelove, J. K., Klein, R. I., McKee, C. F., et al. 1997, ApJL, 489, L179

Tsuribe, T., & Omukai, K. 2006, ApJL, 642, L61

Tsuribe, T., & Omukai, K. 2008, ApJL, 676, L45

Turk, M. J., Smith, B. D., Oishi, J. S., et al. 2011, ApJS, 192, 9

Tweed, D., Devriendt, J., Blaizot, J., Colombi, S., & Slyz, A. 2009, A&A, 506, 647

Van Borm, C., Bovino, S., Latif, M. A., et al. 2014, A&A, 572, A22

Venemans, B. P., Findlay, J. R., Sutherland, W. J., et al. 2013, ApJ, 779, 24

Venemans, B. P., Verdoes Kleijn, G. A., Mwebaze, J., et al. 2015, MNRAS, 453, 2259

Visbal, E., Haiman, Z., & Bryan, G. L. 2014, MNRAS, 445, 1056

Volonteri, M. 2010, A&AR, 18, 279

Volonteri, M., & Bellovary, J. 2012, RPPh, 75, 124901

Willott, C. J., Delorme, P., Reylé, C., et al. 2010, AJ, 139, 906

Wise, J. H., Turk, M. J., & Abel, T. 2008, ApJ, 682, 745

Wise, J. H., Turk, M. J., Norman, M. L., & Abel, T. 2012, ApJ, 745, 50

Wolcott-Green, J., Haiman, Z., & Bryan, G. L. 2011, MNRAS, 418, 838

Wu, X.-B., Wang, F., Fan, X., et al. 2015, Natur, 518, 512

Yahil, A. 1983, ApJ, 265, 1047

Yajima, H., & Khochfar, S. 2015, arXiv:1507.06701

Yoshida, N., Omukai, K., & Hernquist, L. 2008, Sci, 321, 669