Intermittent fragmentation and statistical variations during gas collapse in magnetised atomic cooling haloes

TL;DR

This study uses advanced cosmological magnetohydrodynamic simulations to investigate how magnetic fields and turbulence influence the formation of supermassive black hole seeds during gas collapse in atomic cooling haloes, revealing intermittent fragmentation and high accretion rates.

Contribution

First large eddy simulations incorporating unresolved turbulence models to explore magnetic field amplification and fragmentation during direct collapse in early universe haloes.

Findings

Magnetic fields are amplified to dynamically relevant levels.

Turbulence remains supersonic and super-Alfvénic throughout collapse.

Fragmentation is intermittent with high mass accretion rates.

Abstract

Observations reveal the presence of supermassive black holes (SMBH) as early as ~700 million years after the Big Bang. Their formation path is still subject to current debate. We explore the influence of magnetic fields, which are strongly amplified via the turbulent small-scale dynamo, on the formation of SMBH seeds within the direct collapse scenario. In this study, we perform for the first time cosmological magnetohydrodynamic large eddy simulations that employ a model for unresolved, compressible MHD turbulence. In total we perform 36 simulations for 9 haloes each with two different initial magnetic field strengths, and with and without employing the unresolved turbulence model. We make use of the adaptive mesh refinement approach to achieve an effective spatial resolution of less than one proper astronomical unit. We consider a regime where cooling is regulated by atomic hydrogen and the molecular hydrogen gets dissociated by a strong radiation field. Our main finding is that the majority of the gas properties in the haloes at the final output are predominantly determined by the run-away gravitational collapse. Turbulence is supersonic and super-Alfv\'enic in all cases, and magnetic fields are amplified to an approximately dynamically relevant regime. Finally, fragmentation during the collapse is intermittent and mass accretion rates range from 0.2-3 Msun/yr. This suggests that the presence of strongly amplified magnetic fields and turbulence provides additional pressure support on small scales and make the direct collapse a viable scenario for the formation of massive objects under the required ambient conditions.

Figures

Tables

Table 1. Table 1: Fragmentation results of all 36 simulations, i.e., 9 different haloes (with slightly different initial conditions) each with 4 different setups. The masses of the clumps are given in M ☉  and the distance between their centers of mass is proper au. The asterisk indicates whether the fragmentation was first identified at a peak density of 10 − 9 superscript 10 9 10^{-9}  g/cm 3 and not yet at a peak density 10 − 12 superscript 10 12 10^{-12}  g/cm 3 . The two snapshots are ≈ 150 absent 150 \approx 150 years apart (compared to a free-fall time of ≈ 50 , 000 absent 50 000 \approx 50,000 years at 1 pc). The fragments of halo 8 LES-B10 are visible in the 1 pc projections of Fig.  5 .

	ILES-B12	ILES-B10	LES-B12	LES-B10
halo 1	281 & 384 *	1335 & 1680 *	no	no
	4040 au	23126 au
halo 2	no	no	no	831 & 588
				17908 au
halo 3	no	no	no	no
halo 4	no	no	no	no
halo 5	no	no	no	356 & 156 *
				3227 au
halo 6	no	no	no	2 & 4 *
				64 au
halo 7	5589 & 5185	no	no	no
	56156 au
halo 8	no	28 & 42*	no	3833 & 1753
		672 au		76893 au
halo 9	no	no	128 & 121	no
			3352 au

Equations

$$\sqrt{⟨{\left(\frac{|{\rho }_1(\mathbf{x})-{\rho }_2(\mathbf{x})|}{0.5\left({\rho }_1(\mathbf{x})+{\rho }_2(\mathbf{x})\right)}\right)}^2⟩}=0.09.$$

$${\tau }_{ij}={\tau }_{ij}^{\mathrm{u}}-{\tau }_{ij}^{\mathrm{b}}+{\delta }_{ij}{E}_{\mathrm{sgs}}^{\mathrm{b}},$$

$${\tau }_{ij}^{\mathrm{u}}$$

$${\tau }_{ij}^{\mathrm{b}}$$

$$\mathcal{E}$$

$$M_{\mathrm{J}}(r)=\frac{4\pi c_s^3(r)}{3G^{3/2}{\rho }_0^{1/2}(r)}.$$

Full text

††pubyear: 2019††pagerange: Intermittent fragmentation and statistical variations during gas collapse in magnetised atomic cooling haloes–LABEL:lastpage

1 Introduction

Supermassive black holes (SMBH) with up to a billion solar masses are common in the centres of present-day galaxies (Kormendy & Richstone, 1995; Tremaine et al., 2002), and their presence has been confirmed as early as ${\sim}700$ million years after the Big Bang (Fan et al., 2003; Mortlock et al., 2011; Wu et al., 2015; Schleicher, 2018; Bañados et al., 2018). The formation of such massive black holes (BH) is still an enigma. It is important to obtain a physical understanding how they formed and how efficiently they grew in the first billion years. Different pathways to form supermassive black holes have been suggested (Rees, 1984; Volonteri, 2010; Latif & Ferrara, 2016). There are three main mechanisms that can potentially lead to the formation of SMBHs. They include the accretion and merging of Pop III remnants (Heger & Woosley, 2002; Whalen & Fryer, 2012; Madau & Rees, 2001), the relativistic instability in dense stellar cluster (Devecchi et al., 2010, 2012; Lupi et al., 2014; Katz et al., 2015; Boekholt et al., 2018; Reinoso et al., 2018; Schleicher et al., 2018; Sakurai et al., 2019) as well as the direct collapse of protogalactic gas clouds (Bromm & Loeb, 2003; Lodato & Natarajan, 2006; Wise et al., 2008; Latif et al., 2013c; Latif & Schleicher, 2015b; Becerra et al., 2015; Regan et al., 2014; Wise et al., 2019). With accretion at the Eddington limit, obtaining such masses during the available time is very difficult when starting from a stellar mass black hole (Shapiro, 2005). Given the difficulties with other scenarios like photo-evaporation of HII regions around the first stars, the direct collapse is the most plausible scenario where a metal-free protogalactic halo collapses directly to form a massive BH (Bromm & Loeb, 2003; Begelman et al., 2006; Volonteri & Rees, 2006; Latif et al., 2013c; Wise et al., 2019).

Molecular hydrogen is the only coolant in primordial haloes which can bring the temperature of the gas down to a few hundred Kelvin and may induce fragmentation by reducing the Jeans mass. The gas will show little fragmentation to form stars in zero metallicity haloes in the absence of H2 cooling (Li et al., 2003) which can be achieved via a strong photodissociating UV background in primordial haloes (Dijkstra et al., 2008; Shang et al., 2010; Schleicher et al., 2010). In the presence of only atomic line cooling, the Jeans mass becomes about two orders of larger (for $n=10^{4}\,\mathrm{cm}^{-3}$ and $T=8000$ K, $M_{J}\sim 10^{5}$ M*☉) and fragmentation to lower mass scales remains inhibited. Numerical simulations show that the collapse of metal-free haloes with $T_{vir}>10^{4}$ K irradiated by strong LW flux proceeds isothermally in the absence of H2* cooling and a massive object of $10^{5}-10^{6}$ M*☉* can be formed (Bromm & Loeb, 2003; Wise et al., 2008; Latif et al., 2014b).

Magnetic fields are expected to influence the formation of black holes by enhancing the Jeans mass due to the extra magnetic pressure and providing additional means for the transport of angular momentum by magnetic torques. The latter may become significant in the central accretion disc, implying the presence of strong rotation measures and enhanced accretion rates. In fact, the detection of strong rotation measures in quasars at z = 5.3 indicates the relevance of magnetic fields in the early universe (Hammond et al., 2012). The observations of nearby active galactic nuclei suggest that magnetic fields play a vital role in the transport of angular momentum (Beck et al., 1999, 2005).

Similarly, turbulence can affect the BH formation process by locally compressing the gas and regulating the transport of angular momentum. The hydrodynamical simulations employing a subgrid-scale (SGS) turbulence model to account for unresolved turbulence show that turbulence favours the formation of self-gravitating accretion discs, suppresses fragmentation (due to an additional unresolved turbulent pressure), and that a turbulent viscosity term contributes to the transport of angular momentum (Latif et al., 2013d, a). While fragmentation occasionally occurred, depending on the properties of the host halo, the formation of massive objects still proceeded.

Turbulence also amplifies weak magnetic fields by converting turbulent energy into magnetic energy a process known as the small scale dynamo. Strongly amplified magnetic fields may reach equipartition with the kinetic energy and become dynamically important. Latif et al. (2014a) have shown that strongly ($\mathcal{O}(10^{0}-10^{1})$ G) amplified magnetic fields provide a substantial amount of support on small scales, surpassing thermal pressure support by a factor of ${\sim}10$, and suppress fragmentation. However, these simulation lacked an MHD SGS model that treats the combined effect of unresolved turbulence and magnetic field dynamics.

Small scale physical processes are often not captured by numerical simulations due to the limited spatial resolution. Large eddy simulations (LES) address this issue by incorporating unresolved small scale processes such as turbulent dissipation via an SGS model (Sagaut, 2006; Garnier et al., 2009; Schmidt, 2015). Simulations that treat turbulent dissipation implicitly by means of a shock capturing numerical method are usually referred to as implicit large eddy simulations (ILES, see, e.g., Grinstein et al., 2007). In addition to pure turbulent dissipation other small scale processes have been modelled by an SGS model in the context of astrophysics. These processes include turbulent deflagration in Type Ia supernovae (Ropke et al., 2007), thermal and turbulent feedback from supernovae in isolated disc galaxies (Braun et al., 2014), or turbulent pressure support in the intergalactic medium (Iapichino et al., 2011) and during black hole formation in the early universe(Latif et al., 2013c). All these examples focus on hydro and/or thermodynamic processes. Including magnetic fields and associated small scale processes such as magnetic reconnection, inverse (i.e., up-scale) cascades, and small scale dynamo action present challenges in SGS modelling (Miesch et al., 2015). The majority of MHD SGS models have been derived from their incompressible hydrodynamic counterparts, and only few have been tested in the compressible MHD regime (Miki & Menon, 2008; Chernyshov et al., 2014; Grete et al., 2015). Vlaykov et al. (2016) and Grete et al. (2016, 2017) have developed an SGS model that explicitly takes into account compressibility effects in MHD turbulence and demonstrated its applicability from the subsonic to the highly supersonic regime. Therefore, it is most suitable for the dynamics of a Direct Collapse scenario.

In this study, we perform for the first time cosmological simulations that employ a subgrid-scale model for unresolved MHD turbulence and investigate its impact on the formation of Direct Collapse BHs. We adopt the limit of a very strong radiation background, where molecular hydrogen is fully dissociated so that cooling is regulated only via atomic hydrogen lines. Such conditions tend to strongly suppress fragmentation (Latif et al., 2013c). We study the statistical properties of 9 haloes with different initial magnetic seed field strengths and the treatment of unresolved turbulence in the context of direct collapse scenario. We find that fragmentation is intermittent and the run-away collapse is predominantly determined by the microphysics of the atomic gas, while magnetic fields and turbulence play a secondary role.

The article is organised as follows. In Section 2 we present our simulation setup including details about employed turbulence and chemical models. In Section 3 we present our main results and we summarise our conclusions in Section 4.

2 Method

2.1 Simulation setup

We conduct cosmological zoom-in simulations with Enzo (Bryan et al., 2014) which is an open source multi-physics code111See http://enzo-project.org/. To perform cosmological zoom-in simulations, we first ran dark matter (DM) only simulations with $1024^{3}$ DM particles to identify haloes and reran them including baryons and additional physics. We employed the adaptive mesh refinement (AMR) approach to add dynamical refinement during the collapse.

All simulations start at a redshift of $z=100$ based on the $\Lambda$CDM-model by taking cosmological parameters from Planck 2015 data (Planck Collaboration et al., 2016a). The simulation domain is a cubic box with a side length of 1 Mpc/h and initial conditions are generated222Sample parameter files used for creating initial conditions and running both types of cosmological simulations are available as supplementary online material. with MUSIC (Hahn & Abel, 2011).

We use the Rockstar halo finder (Behroozi et al., 2013) to identify the most massive haloes at $z~{}=~{}12$ and traced back all particles that were within 4 virial radii of that halo to their original positions at $z~{}=~{}100$ to recreate nested grid initial conditions333See a script by John Wise get_halo_initial_extent.py at https://bitbucket.org/jwise77/enzo-mrp-music. The newly generated initial conditions have a root grid with $256^{3}$ cells (and an equal number of DM particles) with two nested grids of the same resolution that enclose all previously identified particles. In each nested grid the spatial resolution doubles so that the effective spatial resolution of the region of interest is initially identical to the previous DM only run. In total, we use 23,024,594 DM particles with an effective DM resolution of 99 M*☉*.

In order to follow the collapsing gas in the multi-physics simulations we use AMR that is triggered by one of the three following conditions: gas overdensity of a factor 4 (with a super-Lagrangian refinement exponent of $-0.3$), dark matter overdensity of a factor 4, and resolving the Jeans length by at least 64 grid cells throughout the simulations. The latter is required in order to resolve turbulence necessary for the small scale dynamo action (Latif et al., 2013b). DM particles are smoothed at level 12 to prevent artificial effects. We stop all simulations at a proper peak density of $10^{-9}$ g/cm3, corresponding to a number density of $\approx 4.9\times 10^{14}$ cm*-3*. This is reached at a redshift of $z\approx 11.8$. At this final time the central region of the halo is represented on a grid at the 29th level of refinement corresponding to an effective spatial resolution of 0.17 au (proper).

We solve the cosmological magnetohydrodynamic equations using the MUSCL-Hancock framework with piecewise linear reconstruction and the HLL Riemann solver. For increased numerical stability, we set Theta_Limiter = 1.0, corresponding to a minmod flux limiter. The divergence constraint of the magnetic field is maintained by using hyperbolic divergence cleaning (Dedner et al., 2002). The chemical model employed in all simulation is described in more detail in subsection 2.3.

In total, we conducted and analysed 36 simulations that vary with respect to initial conditions, initial magnetic field strength, and usage of a subgrid-scale (SGS) model for unresolved compressible MHD turbulence, see the following subsection 2.2 for more details. The initial conditions were varied by using 9 different random seeds in MUSIC. Given that our random seeds varied by 1 the random number generator in MUSIC produced similar cosmological initial conditions444 See https://groups.google.com/forum/#!topic/cosmo_music/BGrUohOlbsE for a discussion.. Here, similar means that the large scales in the initial conditions are virtually identical

as illustrated in Fig. 1 for two initial conditions whose seed varies by 1. The root mean square of the relative density variation between those two initial conditions over the entire domain is

[TABLE]

Over the course of the evolution, the small variations in the initial conditions result in a sample of 9 similar haloes (e.g., with respect to virial radii and masses, see Sec. 3.1) that we use for statistical analysis. The similarity between the 9 haloes allows us to explore the effect of very small variations in the initial conditions on the final small scale halo properties. In other words, we are able to differentiate between robust statistics, which are not sensitive to small initial variations, and statistically varying results that are prone to small physical variations and/or numerical effects. All 9 haloes have been initialised with a proper uniform initial magnetic field strength of both $10^{-10}$ G and $10^{-12}$ G corresponding to comoving magnetic field strengths of $10^{-5}$ nG and $10^{-8}$ nG, respectively. This is well below current estimates of the upper limit of, e.g., $4.4$ nG (comoving at 1 Mpc scale) from the CMB power spectrum (Planck Collaboration et al., 2016b), see Subramanian (2016) for a recent review.In addition, all 9 haloes were simulated with an explicit nonlinear SGS model (the LES case) and without an explicit model (the ILES case). We will refer to the simulations with different physics throughout the paper as ILES-B10, ILES-B12, LES-B10, and LES-B12, respectively, see Table 1 for a list of all simulations.

2.2 Subgrid-scale turbulence model

To include the effect of unresolved, i.e., below the grid-scale, turbulence we employ the nonlinear SGS model for MHD turbulence presented by Vlaykov et al. (2016); Grete et al. (2016). The main features of this model are that it specifically takes compressibility effects into account and allows for energy transfer in both directions. In other words, energy can be transferred from large to small scales (e.g., the forward cascade and dissipative effects) and from small to large scales (e.g., the inverse cascade). The model has extensively been tested a priori (Grete et al., 2016) and a posteriori (Grete et al., 2017) in different turbulent regimes ranging from the subsonic to the highly supersonic regime. Thus, it suitable for application in simulations of dynamic direct collapse scenarios. The model is incorporated via two additional terms in the cosmological (comoving) MHD equations555See (Bryan et al., 2014) for the general set of equations solved..

First, the turbulent stress tensor,

[TABLE]

enters as a source term ($\ldots=-\frac{1}{a}\nabla\cdot\tau$) in the momentum equation with $a$ as the cosmological scale factor. It consists of the turbulent Reynolds stress,

[TABLE]

the turbulent Maxwell stress,

[TABLE]

and the turbulent magnetic pressure given via the turbulent magnetic energy $E^{\mathrm{b}}_{\mathrm{sgs}}$. In general, the turbulent kinetic and magnetic energies are given by the traces of the turbulent stresses via the identities $2E^{\mathrm{u}}_{\mathrm{sgs}}=\tau_{ii}^{\mathrm{u}}$ and $2E^{\mathrm{b}}_{\mathrm{sgs}}=\tau_{ii}^{\mathrm{b}}$. Einstein summation applies and $\Box_{i,j}$ designates the $j$th partial derivative of component $i$ of $\Box$. The filtered density $\overline{\rho}$, the filtered velocity $\widetilde{u}={\overline{\rho u}}/{\overline{\rho}}$, and the filtered magnetic field $\overline{B}$ are calculated via an explicit filter with a filter width of $\Delta=2.711\Delta_{x}$ grid cells (see Grete et al. (2017) for details). The explicit filtering ensures a scale decomposition between the LES scale (on which the SGS terms are calculated) and the grid scale, which is most affected by numerical dissipation.

The second term of the model is the turbulent electromotive force,

[TABLE]

which enters the induction equation as a source term ($\ldots~{}=~{}\frac{1}{a}\nabla\times\bm{\mathcal{E}}$) with Levi-Civita symbol $\varepsilon_{ijk}$.

Finally, it should be noted that this SGS model falls in the category of instantaneous or zero-equation models. There is no additional dynamical equation for the turbulent energies as the energies are determined instantaneously via the turbulent stresses. This also means that there is no intermediate reservoir of energy, as, for example is used in Schmidt et al. (2014) in the context of cosmological hydrodynamics.

2.3 Chemical model

We use the chemical model described in Latif et al. (2016) and briefly summarise here its main features. We solve the chemical and thermal evolution of the gas using the KROME package (Grassi et al., 2014). In our simulations the rate equation of the following species, H, H-, H+, He, He+, He*++*, H2, H${}_{2}^{+}$, and e-, are solved for non-equilibrium and are coupled to MHD. The chemical network and reaction rates are listed in Appendix A of Latif & Schleicher (2015a). It includes all the relevant processes for the formation and dissociation of H2 as well as photo-detachment of H-. We ignore the species involving deuterium as they get dissociated for even weaker radiation fluxes. In this work we assume a background Lyman Werner (LW) flux of strength $10^{5}$ in units of $J_{21}=10^{-21}$ erg/s/cm2/Hz/sr with a fixed radiation temperature of $T_{\mathrm{rad}}=2\times 10^{4}$ K. Such a choice of a strong LW flux ensures that collapse proceeds isothermally as it is well above the critical value of LW flux estimated from 3D simulations.

Furthermore, our model includes cooling due to the collisional excitation, collisional ionization, radiative recombination and Bremsstrahlung radiation. Also the cooling due to the molecular hydrogen, collisionally induced emission and chemical heating and cooling due to the three-body rates is included in our model. At high densities, we also take into account $H^{-}$ cooling as well as employ realistic opacities for bound-free $H^{-}$ and Lyman alpha cooling. For further details about our chemical model, see Section 2.2 of Latif et al. (2016).

3 Results

3.1 Overview

Overall, we observe a self-similar gravitational collapse that results in virtually identical large-scale properties of the most massive halo in all simulations independent of initial magnetic field strength and the unresolved turbulence model. This is expected as the large-scale properties of the initial conditions are practically identical. In general, the most massive halo in every simulation, which we analyse in this paper, has a virial radii of $\approx 0.95$ kpc and a mass of $\approx 5.5\times 10^{7}$ M*☉* (dark matter) and $\approx 8.6\times 10^{6}$ M*☉* (gas). We stop all simulations at the same peak proper density of $10^{-9}$ g/cm3, corresponding to a number density of $\approx 4.9\times 10^{14}$ cm*-3*. The central region has an effective spatial resolution of 0.17 au. This allows us to differentiate between robust and statistically varying or intermittent properties on the small scales across haloes with different initial magnetic field strength and unresolved turbulence model.

Given the overall similarities, we discuss the gas dynamics for only one representative halo (i.e., the set of four simulations with different physics for halo 8 selected arbitrarily) and also discuss the significant differences found in the complete suite of 36 simulations.

3.2 Dynamical properties

The dynamical properties of simulated haloes are discussed in this section. We show the density-temperature phase diagram for a representative halo in Fig. 2. In the ubiquity of a strong LW radiation collapse proceeds isothermally with $T\approx~{}8000$ K. At densities above $10^{-17}$ g/cm3 H*-* cooling becomes efficient and lowers the gas temperature down to $\approx~{}5500$ K., see also the averaged temperature radial profile in Fig. 3(b). This behaviour is consistent with previous studies investigating the impact of H- cooling Latif et al. (2016). The median density is shown in Fig. 3(a) and follows roughly $r^{-2}$ power-law outside the central Jeans length of $\approx 10$ au and gets flattened inside the Jeans radius. The small bumps on the density profiles are due to intermittent fragmentation (see also the density projection in the left column of Fig. 5).

Within the central region, the isothermal collapse stops (the infall velocity is $\approx 0$ km/s, see Fig. 3(c), and the density profile is flattened) and further evolution will lead to the formation of an adiabatic core, see Latif et al. (2016). Outside of this region gas moves with a radial velocity of $\approx-10$ km/s indicating large inflows towards the halo centre. The tangential velocities are approximately constant throughout the inner parsec

ranging between 10-20 km/s as depicted in Fig. 3(d) with the exception of LES-B10. For comparison, the Keplerian velocities $v_{\mathrm{Kep}}~{}=~{}\sqrt{GM_{\mathrm{encl}}/r}$ are shown in Fig. 4. With 15-25 km/s between 10 au and 1 pc the Keplerian velocities are generally equal to or larger than the tangential velocities. In combination with infall velocities of $-10$ km/s this suggests that there is little support against collapse by coherent rotation in these haloes yet. Similarly, no morphological changes from (turbulent) spherical distributions towards disk like structures are observed, see projections at various length scales for one halo in Fig. 5. This equally applies to all other haloes in the sample. Virtually all motion in the halo is supersonic with sonic Mach numbers ranging between $\mathrm{M_{S}}=2-4$, see radial profiles in Fig. 3(f). Particularly the differences between ILES-B12 and LES-B12 in the inner region (e.g., $\mathrm{M_{S}}<2$ for ILES-B12 and $\mathrm{M_{S}}>4$ for LES-B12) may be attributed to the statistical variations, see Section 3.6.

3.3 Mass profiles and fragmentation

The enclosed baryon mass versus radius is illustrated in Fig. 6(a). At the end of the simulations the central core within $\approx 10$ au (proper) has a total mass a few M*☉. The entire unstable core with $M_{\mathrm{encl}}/M_{\mathrm{J}}>1$ extends to approximately 1 pc with $M_{\mathrm{J}}~{}\approx~{}10^{5}$M☉* as shown in Fig. 6(b) where the ratio of the enclosed gas mass $M_{\mathrm{encl}}$ to the Jeans mass $M_{\mathrm{J}}$ versus radius is illustrated. $M_{\mathrm{J}}$ is given by

[TABLE]

Here, $c_{s}(r)$ is the median sound speed at radius $r$ and $\rho_{0}(r)~{}=~{}3M_{\mathrm{encl}}(r)/4\pi r^{3}$ the mean density. For a gravitational collapse to occur this ratio $M_{\mathrm{encl}}/M_{\mathrm{J}}$ should be larger than one. Based on these estimates most of the gas within the central pc region is unstable to gravitational collapse. The typical mass accretion rates shown in Fig. 6(c) are between $0.2-3$ M*☉/yr. The mass accretion rate at the edge of the unstable core is about 0.5 M☉*/yr. If this mass accretion rate persists for 1 Myr, we expect the central core to reach $5\times 10^{5}$ solar masses within the 1 Myr. Thus, they are well above the minimum value required to form a direct collapse black hole (Schleicher et al., 2013; Hosokawa et al., 2012; Hosokawa et al., 2013; Haemmerlé et al., 2018). The spikes in the mass accretion are due to the substructure in the halo and no strong differences are observed between LES and ILES runs.

In order to quantify the fragmentation in our simulations we employ yt’s clump finder666 We consider density levels in factors of 2 and only track clumps that are gravitationally bound throughout the process. to find gravitationally bound clumps within the haloes at two different times: at a peak density of $10^{-12}$ g/cm3 and at a peak density of $10^{-9}$ g/cm3. The time difference between the two peak densities is $\approx 150$ yr (for a reference,the free-fall time at 1 pc is $\approx 50,000$ yr).

Table 1 lists the number of clumps and their masses in all simulations. Only in 9 out of 36 simulations more than one clump is formed and in those cases usually there are two clumps of similar masses.

However, these clumps are formed at small scales and are expected to merge within a dynamical time based on estimates of Latif & Schleicher (2015b). Overall, our results suggest that the process of fragmentation is stochastic and intermittent for our chosen setup and under the assumption of atomic hydrogen cooling. Given the large accretion rates and expected merging of clumps the formation of a central massive object is the most likely outcome.

3.4 Magnetic field dynamics

The magnetic field strength increases towards the centre of the halo, as shown in Fig. 7(a) for one reference halo. No significant differences in magnetic field amplification are observed between LES and ILES runs. Simulations with larger initial seed fields have larger magnetic field strength but the overall amplification is the same for both cases.

To quantify the amplification of magnetic fields by the small scale dynamo, we have divided the magnetic field strength by $\rho^{2/3}$, the maximum possible contribution from flux freezing under spherical symmetry, and the result is shown in Fig. 7(b). At the 0.1-1 pc scale the magnetic field amplification by the small scale dynamo action is a factor of a few and increases up to a factor of 10 in the centre. These results are in agreement with previous studies showing a critical resolution of at least 32-64 cells per Jeans length for dynamo action to occur (Sur et al., 2010; Federrath et al., 2011; Latif et al., 2013b). In principle, we note that a larger amplification would be expected at even higher resolutions, as the Reynolds number of the flow is effectively limited by the numerical resolution. Finally, for this particular halo the differences in the inner region (e.g., amplification of ${\lesssim}10$ for ILES-B12 and ${\sim}20$ for LES-B12) may indicate differences between the runs (e.g., between with and without SGS model), but our sample of different haloes allows us to attribute this to statistical variations as detailed in Section 3.6.

The Alfvénic Mach number $\mathrm{M_{A}}$, i.e., the ratio of gas velocity to local Alfvén velocity $v_{\mathrm{A}}=B/\sqrt{\rho}$, is shown in Fig. 7(c) and is a proxy for the ratio of kinetic to magnetic energy density. Kinetic motion clearly dominates the dynamics outside of the central region with $\mathrm{M_{A}}\approx 10^{3}$ (for ILES-B10 and LES-B10) and $\mathrm{M_{A}}\approx 10^{4}$ (for ILES-B12 and LES-B12) at $1$ pc. $\mathrm{M_{A}}$ decreases towards the centre in agreement with an active small-scale dynamo that converts kinetic energy to magnetic energy. For weaker initial seed fields the Alfvénic Mach number is 100 in the central region while for stronger initial fields the Alfvénic Mach number drops below 10 in the core but remains super-Alfvénic throughout. This suggests that even in the strong field case motions remain primarily driven by kinetic dynamics rather than magnetic field dynamics.

3.5 Pressure support

Apart from thermal pressure, the presence of magnetic fields introduces an additional pressure, the so-called magnetic pressure, that may support the gas against fragmentation. Furthermore, in our LES two additional pressure terms, the turbulent kinetic and turbulent magnetic pressure, are given by $E^{\mathrm{u}}_{\mathrm{sgs}}$ and $E^{\mathrm{b}}_{\mathrm{sgs}}$, respectively. These terms take into account the cumulative effect of numerically unresolved velocity and magnetic field fluctuations. In the following, turbulent kinetic and magnetic pressure always refer to subgrid scales as opposed to the resolved kinetic and magnetic pressure, which encompasses all (turbulent and non-turbulent) contributions from scales larger than the grid resolution scale. The ratio of the three non-ideal hydrodynamic pressures to the thermal pressure for all LES is illustrated in Fig. 8.

In the case of weak initial seed fields (LES-B12, left panel), both resolved and turbulent magnetic pressures are weak (<$10^{-4}$ relative to the thermal pressure). For strong initial seed field cases the resolved magnetic pressure reaches close to equipartition with the thermal pressure within the central core. This suggests that apart from some variations from halo to halo the presence of magnetic fields in the central region becomes significant. In all cases the turbulent magnetic pressure is approximately 10-20% of the resolved magnetic pressure indicating that the overall magnetic field dynamics is just barely resolved. The ratio of turbulent kinetic pressure to thermal pressure remains constant ($\approx 10^{-2}$) for both cases.

3.6 Statistical variations

The majority of quantities discussed in previous sections are in agreement between the 9 different haloes and across the 4 configurations. A few exceptions, e.g., the tangential velocity $v_{\mathrm{tang}}$ in Fig. 3 (d), the sonic Mach number $\mathrm{M_{S}}$ in Fig. 3 (f), or the magnetic field amplification $B/\rho^{2/3}$ in Fig. 7 (b), exhibit variations between the configurations for the specific halo discussed before. Fig. 9 illustrates the variations in the radial profiles of these quantities for all haloes and configurations in our sample. Overall, the variations in different haloes for each configuration individually is significant ($1\lesssim\mathrm{M_{S}}\lesssim 5$, and between one and two orders of magnitude for $v_{\mathrm{tang}}$ and $B/\rho^{2/3}$, respectively) as shown in the left panels. The differences observed between configurations for an individual halo are negligible compared to the statistical variation of the entire sample of haloes. This is illustrated in the right panels of Fig. 9 where the mean and standard deviations of the median profiles (of all haloes) are drawn on top of each other. All profiles (of different configurations) lie within the variations of each other. Thus, there is no significant statistical difference between the configurations.

In general, all profiles shown including the ones discussed in previous sections are well converged from 10 pc outwards to the virial radius (0.95 kpc). This highlights the similarities between the large-scale properties of haloes resulting from the similar initial conditions used in combination with the self-similar behavior of isothermal collapse.

4 Conclusions and discussion

We performed 36 cosmological magnetohydrodynamic large eddy simulations (implicit and explicit) with adaptive mesh refinement, and followed the collapse of 9 similar metal-free massive haloes ($\approx 5.5\times 10^{7}$ M*☉* at $z\approx 11.8$) down to sub-au scales.All simulations employed a realistic chemical model for 9 species and assume a super-critical LW background flux from a nearby star-forming galaxy. In this study, we particularly explored the influence of magnetic fields and turbulence on the formation of supermassive black hole seeds within the direct collapse scenario. Thus, each of the 9 haloes was simulated using 4 different setups: a weak and a strong initial magnetic field field, and with or without a subgrid-scale model for unresolved MHD turbulence.

We find that, by and large, the self-similar run-away gravitational collapse is the main driver in determining the majority of the halo properties. The presence of (explicitly modelled) turbulence or strongly amplified magnetic fields is secondary. For example, all simulations agree with an isothermal collapse at $\approx 8000$ K exhibiting an $r^{-2}$ density profile. Within the central parsec, all velocities are supersonic ($\mathrm{M_{S}}\approx 2-4$) with infall velocities of $\approx 10$ km/s and tangential velocities of $\approx 15$ km/s. No significant coherent rotational support against the collapse is observed in any simulation yet. The resulting mass accretion rates vary between $0.2-3$ M☉/yr, and, thus, are large enough to support the direct collapse scenario.

Overall, fragmentation is quite intermittent in our simulations and appears as a stochastic process as it occurred mildly in 9 out of 36 simulations. Even clumps forming in these cases are expected to merge within a dynamical time scale.

In all cases, we observe magnetic field amplification beyond simple compression of frozen-in magnetic fields, indicating small-scale dynamo action. In the strong magnetic field case the final magnetic fluxes reach dynamically relevant strengths. Given that the seed field in our strong field case is still $>10^{5}$ times weaker than the upper limit from CMB power spectra (Planck Collaboration et al., 2016b) a further increase of the initial field strength is still realistic. These stronger seed fields may be amplified to dynamically relevant strengths at earlier time during the collapse and may potentially affect gas dynamics. We also note that our results on the magnetic field amplification correspond only to lower limits, as the Reynolds number in the simulations is limited by the numerical resolution, and also the analysis through the MHD SGS model shows that the MHD turbulence is currently under-resolved.

Systematic differences between ILES and LES are marginal in averaged radial profiles, although the turbulent magnetic pressure on the grid scale contributes about 20% to the total magnetic pressure – a significant fraction in the context of the fundamental assumption of explicit LES. While it has been found that the Jeans length needs to be resolved by at least $32-64$ cells to observe the onset of small-scale dynamo action (Sur et al., 2010; Federrath et al., 2011; Latif et al., 2013b), our results suggest that in some cases even a higher resolution is required to fully resolve magnetic fields. Strongly amplified fields are expected to stabilise collapse and help in transferring angular momentum making the Direct Collapse scenario scenario more feasible.

Finally, our simulated sample of haloes is limited. All 9 different initial conditions are virtually identical on the large scales and, thus, result in a very similar large scale evolution of the simulations, e.g., with respect to the virial properties of the halo during collapse. This allowed us to study the intrinsic variations of a run-away collapse, cf., the butterfly effect in cosmological simulations (Genel et al., 2019), in one particular environment. We demonstrated that some quantities, such as the magnetic field amplification factor within the core of the halo, vary by almost two orders of magnitude for a fixed configuration. Moreover, the large intrinsic variation even for haloes with similar properties outshines statistical changes introduced by the LES approach. Eventually, larger samples of similar haloes in combination with larger samples of different haloes (or environments) are necessary, requiring a carefully chosen ensemble of cosmological simulations in order to make precise quantitative statements.

Acknowledgements

The authors thank Brian W. O’Shea for useful discussions. PG acknowledges funding by NASA Astrophysics Theory Program grant #NNX15AP39G. ML thanks the funding from United Arab Emirates University via startup grant No 31S372. DRGS thanks for funding through Fondecyt regular (project code 1161247), Conicyt PIA ACT172033, Conicyt Programa de Astronomia Fondo Quimal QUIMAL170001, Conicyt Programa de Astronomia Fondo ALMA 31160001 and the “Concurso Proyectos Internacionales de Investigación, Convocatoria 2015” (project code PII20150171). The simulations were performed and analyzed with the HLRN-III facilities of the North-German Supercomputing Alliance under Grant No. hhp00039. Enzo and yt are developed by a large number of independent researchers from numerous institutions around the world. Their commitment to open science has helped make this work possible.