Direct Collapse to Supermassive Black Hole Seeds: Comparing the AMR and SPH Approaches
MMon. Not. R. Astron. Soc. , 000–000 (0000) Printed 26 September 2018 (MN L A TEX style file v2.2)
Direct Collapse to Supermassive Black Hole Seeds:Comparing the AMR and SPH Approaches
Yang Luo (cid:63) , Kentaro Nagamine , † , Isaac Shlosman , ‡ Theoretical Astrophysics, Department of Earth & Space Science, Osaka University, 1-1 Machikaneyama, Toyonaka, Osaka 560-0043, Japan Department of Physics & Astronomy, University of Nevada Las Vegas, 4505 S. Maryland Pkwy, Las Vegas, NV 89154-4002, USA Department of Physics & Astronomy, University of Kentucky, Lexington, KY 40506-0055, USA
Accepted ?; Received ??; in original form ???
ABSTRACT
We provide detailed comparison between the AMR code Enzo-2.4 and the SPH/ N -body code GADGET-3 in the context of isolated or cosmological direct baryonic col-lapse within dark matter (DM) halos to form supermassive black holes. Gas flow isexamined by following evolution of basic parameters of accretion flows. Both codesshow an overall agreement in the general features of the collapse, however, many sub-tle differences exist. For isolated models, the codes increase their spatial and massresolutions at different pace, which leads to substantially earlier collapse in SPH thanin AMR cases due to higher gravitational resolution in GADGET-3. In cosmologicalruns, the AMR develops a slightly higher baryonic resolution than SPH during halogrowth via cold accretion permeated by mergers. Still, both codes agree in the buildupof DM and baryonic structures. However, with the onset of collapse, this differencein mass and spatial resolution is amplified, so evolution of SPH models begins to lagbehind. Such a delay can have effect on formation/destruction rate of H due to UVbackground, and on basic properties of host halos. Finally, isolated non-cosmologicalmodels in spinning halos, with spin parameter λ ∼ . − .
07, show delayed col-lapse for greater λ , but pace of this increase is faster for AMR. Within our simulationsetup, GADGET-3 requires significantly larger computational resources than Enzo-2.4 during collapse, and needs similar resources, during the pre-collapse, cosmologicalstructure formation phase. Yet it benefits from substantially higher gravitational forceand hydrodynamic resolutions, except at the end of collapse. Key words: methods: numerical — galaxies: formation — galaxies: high-redshift —cosmology: theory — cosmology: dark ages, reionization, first stars
The origin of supermassive black holes (SMBHs) in thegalactic centers remains an unresolved problem in astro-physics. Unless the SMBHs are primordial, a number of al-ternative scenarios exists. First, they could grow from relicsof the Population III stars (e.g., Madau & Rees 2001; Abelet al. 2002; Bromm & Larson 2004; Volonteri & Rees 2005)— the largely pre-galactic objects whose masses appear tobe comparable or slightly higher than normal OB stars (e.g.,Turk et al. 2009; Hosokawa et al. 2011; Wise et al. 2012; Hi-rano et al. 2014; Fraser et al. 2015, see also Bromm 2013 forreview). Second, they could form from the runaway collapse (cid:63)
E-mail:[email protected] † E-mail:[email protected] ‡ E-mail: [email protected] of compact stellar clusters, subject to general relativistic ef-fects (e.g., Zel’dovich & Podurets 1965; Shapiro & Teukolsky1985), or stellar/gasdynamics evolution of stellar clusters(e.g., Begelman & Rees 1978). Most appealing, however, isthe third alternative which involves a direct baryonic col-lapse to form massive black hole seeds of ∼ − M (cid:12) ,(e.g., Haehnelt & Rees 1993; Loeb & Rasio 1994; Bromm &Loeb 2003; Koushiappas et al. 2004; Begelman et al. 2006;Volonteri & Rees 2006; Begelman et al. 2008; Begelman &Shlosman 2009; Milosavljevi´c et al. 2009; Mayer et al. 2010;Schleicher et al. 2010; Johnson et al. 2011; Choi et al. 2013,2015; Latif et al. 2013; Prieto et al. 2013; Shlosman et al.2016) which experience a substantial growth subsequently.Depending on the chemical composition of the collapsinggas, this can take place when suitable dark matter (DM)halos appear in the universe. For a primordial composition,hydrogen can be either atomic, with a cooling floor of just c (cid:13) a r X i v : . [ a s t r o - ph . GA ] M a r Yang Luo, Kentaro Nagamine, and Isaac Shlosman below 10 K, or molecular, which is able to cool down byadditional ∼ ∼ M (cid:12) . In the latter case, the critical halomasses will be smaller.Direct collapse scenario involves gravitational collapsefrom typical spatial scales of ∼ ∼ − N -body code (Springel 2005).Previous comparison between the grid and Lagrangiancodes have both achieved results within ∼ SPH particles,compared to 512 cells required for the same purpose in AMR. Unfortunately, this comparison differs fundamentallyfrom the one performed in the present work: it is entirelynon-gravitational. As a result, the density stratification wasmuch weaker than one observes during gravitational col-lapse, and hence the evolution in their simulation did notgo as far as in our current work.It remains unclear, whether SPH or AMR codes aremore accurate in general simulation regimes. Both meth-ods have their shortcomings. For example, in AMR, if thefluid moves rapidly across the mesh, the truncation errorslead to significant departure from the same simulation ifthere is no bulk velocity (Tasker et al. 2008). Hopkins (2013)have studied a number of SPH modifications, in particularalternative SPH equations of motion that guarantee con-servation and improved treatment of fluid contact disconti-nuities (e.g., Kelvin-Helmholz instability). Good agreementhas been reported for supersonic flows and associated strongshocks, but the SPH suppresses mixing in subsonic, thermalpressure-dominated regimes. Saitoh & Makino (2013) havealso provided a new density-independent formulation of SPHwhich removes the artificial surface tension. Their methodwas more formally treated in the Lagrangian formalism byHopkins (2013).Another important issue for classical AMR codes is theangular momentum conservation. Unless we can anticipatethe geometry of the flow and adjust accordingly the mesh sothat the fluid velocity is perpendicular to the boundary, theaccuracy is seriously compromised and subject to artificiallyenhanced diffusion. Finally, the refinement criteria is some-what arbitrary and can produce various artefacts. At therefinement boundaries, there is a significant loss of accuracyas the refinement is necessarily discontinuous.On the other hand, the SPH in 3-D suffers from a poorshock resolution, noise on the scale of the smoothing ker-nel, and low-order accuracy for the treatment of contactdiscontinuities. Hydrodynamic instabilities like the Kelvin-Helmholtz can be suppressed (Agertz et al. 2007), althoughrecent modifications put forth by Hopkins (2013) has beenclaimed to correct these problems. In real fluids, the entropyis raised in shocks because particle collisions randomize theirvelocities which increase the heat and the entropy. The SPHdoes not capture shocks properly because entropy must begenerated locally to dissipate the small scale velocities. Inorder to mimic this process, and to prevent interpenetra-tion of the SPH particles at the shock locations, artificialviscosity is required. The introduction of the latter usuallyresults in an overly viscous fluid away from shocked regions.To overcome this, the Balsara (1995) switch or Cullen &Dehnen (2010) formulation have been developed.The aim of this paper is to compare the ability of AMRand SPH codes in following the direct collapse of gas with aprimordial composition, i.e., involving identical atomic hy-drogen cooling. We compute and analyze models of gas col-lapse inside DM halos in isolated and cosmological settings.Our models span more than 7 orders of magnitude in radius,from 1 kpc down to 10 − pc. Of course, the main question iswhether both codes can describe the evolution that leads tothe same end product. However, specific details are impor-tant as well, namely, are various physical quantities (e.g.,gas density, temperature, tangential and radial velocity, andangular momentum distributions) similar at various timesof the collapse? Does it take the same time to collapse, and c (cid:13) , 000–000 omparing AMR and SPH Approaches in Direct Collapse Modeling how do accretion rates and their temporal and radial distri-butions compare?This paper is structured as follows. In Section 2, webriefly describe the AMR and SPH codes used. Section 3provides the results of our comparison runs, and in the lastsection, we discuss our results and summarize them. The AMR code Enzo-2.4 has been tested extensively andis publicly available. It uses a multi-grid particle-mesh N -body method to calculate the gravitational dynamics, in-cluding collisionless DM particles, and a second-order piece-wise parabolic method (PPM, Colella & Woodward 1984;Bryan et al. 1995) to solve hydrodynamics. ZEUS hydro isalso available, but we use PPM here.The structured AMR used in Enzo places no fundamen-tal restrictions on the number of rectangular grids used tocover some region of space at a given level of refinement, oron a number of refinement levels (Berger & Colella 1989).A grid is refined by a factor of 2 in lengthscale, if either thegas or DM density become greater than ρ N (cid:96) , where ρ isthe cosmic mean density for the gas or DM, respectively.The refinement factor is N = 2, and (cid:96) refers to the AMRrefinement level.The Jeans length has been resolved by at least 4 cellsin these simulations to satisfy the Truelove et al. (1997)requirement for resolution. For the SPH code, an equivalentcriterion is to resolve the local Jeans mass by having atleast twice the number of particles in an SPH kernel (Bate& Burkert 1997).In GADGET-3, the gas and the DM are represented byparticles. The N -body solver is essentially the same as inGADGET-2 (Springel 2005), with some improvements foroptimization purposes, such as the domain decomposition.To achieve the spatial adaptivity at a moderate computa-tional cost, it uses a hierarchical multipole expansion, i.e.,the tree algorithm (Barnes & Hut 1986). The gravitationalpotential is softened below a spatial scale specified by thegravitational softening length. In principle, this kernel (andthe associated softening length, (cid:15) , which represents the forceresolution) can differ from the smoothing length used in thehydrodynamical method, as we describe below. But mostlythe SPH uses the same cubic-spline kernel, with a fixed (cid:15) evaluated as the mean initial interparticle distance, dividedby a fudge factor of ∼ − α = 0 .
2, compared to usualconstant artificial viscosity. With this method, the gas be-comes virtually inviscid away from the shocks, while main-taining particle order. The viscosity decay length is taken as3.73, and the source scaling equal to 2.77.For the cooling package, we use Grackle version 2.0 forGADGET-3, and Grackle version 1.1 for Enzo. These ver- sions for the SPH and AMR codes have identical chem-istry and cooling (Bryan et al. 2014b; Kim et al. 2014,https://grackle.readthedocs.org/). This library was bornout of the chemistry and cooling routines of the Enzo simula-tion code. The Grackle solves for radiative cooling and inter-nal energy, calculates cooling time, temperature, pressure,and ratio of specific heats. It uses a non-equilibrium primor-dial chemistry network for atomic H and He (Abel et al.1997; Anninos et al. 1997), H and HD, Compton coolingoff the cosmic microwave background, tabulated metal cool-ing and photo-heating rates calculated with Cloudy (Fer-land et al. 2013). In addition, Grackle provides a look-uptable for equilibrium cooling. The gas is assumed to be dust-free and optically-thin and the metals are assumed to be inionization equilibrium. The cooling rate for a parcel of gaswith a given density, temperature, and metallicity, that isphotoionized by incident radiation of known spectral shapeand intensity can be pre-calculated. As we focus on gravita-tional collapse models at z >
10, and limit the runs to theoptically-thin regime, the UV background is not included,and we neglect the radiative transfer of ionizing photons inthe present work.
For isolated models, we adopt the WMAP5 cosmological pa-rameters (Komatsu et al. 2009), namely, Ω m = 0 . , Ω b =0 . , h = 0 . h is the Hubble constant in unitsof 100 km s − Mpc − . We set up the details of an isolatedDM halo that is consistent with the cosmological contextthat we work with. Therefore, some of the halo parametersare specified with units that include the Hubble parameter,although we use physical quantities (not comoving) in allisolated case calculations. In both codes, we define a DMhalo as having the density equal to the critical density ∆ c times the mean density of the universe, ρ b , which dependson the redshift z and the cosmological model. The top-hatmodel is used to calculate ∆ c ( z ), and the density is calcu-lated within a virial radius, R vir . The halo virial mass is M vir ( z ) = (4 π/ c ( z ) ρ b R vir .For isolated models we work in physical coordinates.We assume that the gas fraction in the model is equal to theuniversal ratio, and simulate the gas evolution within DMhalos of a virial mass of M vir = 2 × h − M (cid:12) and a virialradius R vir = 945 h − pc. The initial temperature of the gasis taken to be T = 3 . × K. The simulation domain isa box with a size L box = 6 kpc centered on the halo. Inthe following, we abbreviate the spherical radii with R andcylindrical ones with r .The initial DM and gas density profiles ρ DM ( R ) and ρ g ( R ) are those of an nonsingular isothermal spheres witha flat density core of R c , DM = 0 .
34 pc for DM and R c , g =142 .
65 pc for the gas, respectively (Fig. 1): ρ g ( R ) = ρ g , RR c , g ) ( R (cid:54) R vir )1(1 + R vir R c , g ) R R ( R > R vir ) , (1) ρ DM ( R ) = ρ DM , R , DM R ( R (cid:54) R vir ) , (2)where ρ DM , = 7 . × − g cm − , and ρ g , = 1 . × c (cid:13)000
65 pc for the gas, respectively (Fig. 1): ρ g ( R ) = ρ g , RR c , g ) ( R (cid:54) R vir )1(1 + R vir R c , g ) R R ( R > R vir ) , (1) ρ DM ( R ) = ρ DM , R , DM R ( R (cid:54) R vir ) , (2)where ρ DM , = 7 . × − g cm − , and ρ g , = 1 . × c (cid:13)000 , 000–000 Yang Luo, Kentaro Nagamine, and Isaac Shlosman
Table 1.
Simulation parameters for the isolated halo models withEnzo AMR and GADGET-3 SPH. Name of the run and the val-ues of spin parameter λ is listed. “I” in the model abbreviationstands for isolated runs, and the number λ ×
100 for DM halos.Cosmological models (not shown in the Table), are abbreviatedby “C” and by λ × λ AMR-I1 0.01AMR-I3 0.03AMR-I5 0.05AMR-I7 0.07SPH-I1 0.01SPH-I3 0.03SPH-I5 0.05SPH-I7 0.07
Figure 1.
Initial density profiles in the isolated halo model forthe gas (solid line) and DM (dashed line) as a function of radius. − g cm − . Such DM halos are similar to the NFW halos(Navarro et al. 1997), and are simple to construct. The DMhalo core size is actually set by the gravitational softeninglength for GADGET-3 SPH.The DM halo rotation is defined in terms of the cosmo-logical spin parameter λ , λ = J √ M vir R vir v c (3)where J is the total angular momentum of the DM halo, M vir the halo virial mass, and v c the circular velocity at R vir ,which is constant with radius for an isothermal sphere massdistribution. The mean DM spin is λ ∼ . ± .
005 (e.g.,Bullock et al. 2001). We explore the range of λ (cid:39) − . λ , we sample the sameinitial density distribution using four different random num-ber sequences. The properties of the collapsing baryons havebeen averaged over these sequences. Models having identi-cal mass distribution and differing only in the initial randomsequence, are abbreviated similarly.To produce the DM halos with a pre-specified λ for iso- lated halo models, we follow the the prescription of Longet al. (2014). In short, we set the isotropic distribution ofvelocity dispersion in DM. Then reverse tangential veloci-ties for a fraction of DM particles to reproduce the specificangular momentum described above, with λ equal to the re-quired value. For both SPH particles and the gas in AMRgrid cells, we calculate the average tangential velocities ofthe background DM in cylindrical shells, accounting for thedependence along the (rotation) z -axis. The rotational ve-locities for the gas in the equatorial plane of the DM haloare given by v t ( r ) = v × (cid:26) r/R c , g if r (cid:54) R c , g , R c , g (cid:54) r (cid:54) R vir , (4)where v is tangential velocity defined by λ . The radial veloc-ity dispersion in the center is given by σ = (cid:112) GM vir / R vir =21 . − , where G is the gravitational constant. The iso-lated halo models have been listed in Table 1 and are identi-fied with “I” and the value of the spin parameter multipliedby 100. In cosmological models, we use the same notation,but replace “I” with “C”.In the isolated models of GADGET-3, the DM res-olution is set by the fixed gravitational softening length, (cid:15) DM = 0 .
37 pc. This is done in order to match the conditionin Enzo runs. For Enzo, it corresponds to an initial root gridof 64 in a 6 kpc region with a maximal refinement level of8 allowed for gravity, (cid:15) DM , min = 6000 / / = 0 .
37 pc.For the gas, the gravitational softening is (cid:15) g = 10 − pc,which is a fixed number for GADGET-3 and serves as a min-imal value for Enzo. In other words, we use the initial res-olution of 512 SPH particles or cells for baryons, and 100 particles or particles-in-mesh for the DM. All other param-eters in Enzo and GADGET-3 runs are similar. The forceresolution in adaptive PM codes is twice the minimal cellsize (e.g., Kravtsov et al. 1997). Isolated and cosmologicalruns and their parameters are summarized in Tables 1 and2. We use the models AMR-I5 and SPH-I5 as representativeisolated models. For further comparison, we use zoom-in cosmological sim-ulations with the same composition as in isolated modelsto follow up the gas evolution within DM halos. The ini-tial conditions (ICs) are generated using WMAP5 cosmol-ogy: Ω Λ = 0 . m = 0 . b = 0 . h = 0 . σ = 0 . n = 0 . z = 200.The same cooling packages have been used for cosmologicalmodels as for isolated ones (Section 2.1).For the initial setup, we use the MUSIC algorithm(Hahn & Abel 2011) to generate cosmological zoom-in ICs.MUSIC uses a real-space convolution approach in conjunc-tion with an adaptive multi-grid Poisson solver to generatehighly accurate nested density, particle displacement, andvelocity fields suitable for multi-scale zoom-in simulationsof structure formation in the universe.Generating a set of zoom-in ICs is a two-step process.First, we generate 3 h − Mpc comoving 256 DM-only ICsfor the pathfinder simulation and run it without AMR until z = 10. Using the HOP group finder (Eisenstein & Hut1998), we select an appropriate DM halo, whose mass is > ∼ h − M (cid:12) at z = 10. Second, we generate 0 . h − Mpc c (cid:13) , 000–000 omparing AMR and SPH Approaches in Direct Collapse Modeling Table 2.
Summary of simulation setup for Enzo-2.4 (AMR) and GADGET-3 (SPH) simulations.Parameters Isolated model Cosmological modelSimulation Volume 6 kpc box 3 h − Mpc with a zoom-in region of 0 . h − MpcInitial baryonic resolution 64 root-grid with 3 levels of refinement; 256 root-grid and SPH particles;Effective 512 zoom grid and SPH particles Effective 2048 zoom grid and SPH particlesNumber of DM particles 100 for outer region;Effective 2048 for zoom region † Mass resolution 532 . M (cid:12) for DM particles 245 . h − M (cid:12) for DM particles0 . M (cid:12) for gas particles 44 . h − M (cid:12) for gas particlesGravitational softening (cid:15) DM , min = 0 .
37 pc for DM in Enzo (cid:15) DM , min = 0 . h − pc for DM in Enzo (cid:15) DM = 0 .
37 pc for DM in GADGET (cid:15) DM = 0 . h − pc for DM in GADGET (cid:15) g , min = smallest cell size for gas in Enzo (cid:15) g , min = smallest cell size for gas in Enzo (cid:15) g = 10 − pc for gas in GADGET (cid:15) g = 10 − h − pc for gas in GADGETMinimum gas smoothing length η min − pc 10 − h − pc † The size of the zoom-in region is about one-tenth of the outer region, therefore the actual grid/particle number in the zoom region isinitially ∼ for both baryons and DM particles. Maximum refinement level for DM gravity is set to 8 levels for Enzo in isolated models, i.e., (cid:15) DM , min = 6000 / / = 0 .
37 pc. Theisolated models do not have Hubble expansion in space, so there is no factor of h − and all numbers are physical. Maximum refinement level for DM gravity is set to 15 levels for Enzo in cosmological models, i.e., (cid:15) DM , min = 3 × / / = 0 . h − pc. Maximum refinement levels for the gas in Enzo are set to 17 levels in isolated models and 20 levels in cosmological models.
ICs with 2048 effective resolution in DM and gas, embeddedin the lower resolution outer region. Since we use the samerandom seeds for these ICs as the first step, the phases ofboth ICs are identical. The zoom-in region is centered on theselected halo position and is set to be large enough to coverthe initial positions of all selected halo particles. We performthe zoom-in procedure for each of the targeted halos.In the cosmological runs, the gravitational softening is10 − pc in the comoving coordinates, in the SPH run. Wehave measured the spin parameter of many halos in therange of λ ∼ . − .
07 in our zoom region at z ∼ λ ∼ .
04 as a representative one for both AMR andSPH runs.
We aim at comparing results for direct baryonic collapsewithin DM halos by Enzo-2.4 AMR and GADGET-3 SPHcodes. As a first step, we compare the isolated models. Thecosmological halos have a range of spin parameter, as dis-cussed in § λ (Table 1 and § § λ on the dynamics of the collapse inisolated and cosmological models ( §§ Direct collapse of isolated models with atomic gas has beenperformed and analysed by Choi et al. (2013) for λ = 0 . t = 0 at the start of the run, and stop therun when the central collapse reaches the final resolution of R fin ∼ − pc, in both Enzo and GADGET-3, in order toallow a meaningful comparison, because the SPH run would have difficulty reaching a higher resolution without invokingcomplicated procedures, such as SPH particle splitting.The DM in isolated models has isothermal profiles byconstruction, i.e., ρ DM ∝ R − . They remain nearly identi-cal to ICs. As the collapse proceeds and reaches the centralregion, the DM is dragged inward by the gas, and is com-pressed adiabatically, forming new cusps with slopes ∼ − ∼ − ρ ∝ R − (Larson 1969; Penston 1969). But significant differencesare obvious, as angular momentum is present and modifiesthe solution, and the DM dominates the gravitational po-tential everywhere during the initial stage.When collapsing gas reaches the centrifugal barrier, ashock forms, where the density jumps and the radial velocitydrops abruptly (top panels of Figure 2). The gas accumulatesbehind the shock and forms a disk which grows in mass. Asthe gas density in the inner disk surpasses that of the DM,the 2nd stage of the collapse ensues, and proceeds with an in-creasingly short timescale. Basically, most of the gas behindthe shock collapses and reaches the prescribed resolution of ∼ − pc (e.g., Choi et al. 2013), although here we com-pare the final distributions when R fin ∼ − pc has beenreached.The current AMR simulations with various λ follow ba-sically the same outline. The initial positions of shocks inall runs are located at ∼ λ , as the centrifugal barrier is reachedat progressively larger radii. Figure 3 displays the shock lo-cations as function of time, from their formation time tothe end of the simulations. Note, that for a specific λ , bothAMR and SPH shocks strictly follow each other. Moreover, c (cid:13)000
We aim at comparing results for direct baryonic collapsewithin DM halos by Enzo-2.4 AMR and GADGET-3 SPHcodes. As a first step, we compare the isolated models. Thecosmological halos have a range of spin parameter, as dis-cussed in § λ (Table 1 and § § λ on the dynamics of the collapse inisolated and cosmological models ( §§ Direct collapse of isolated models with atomic gas has beenperformed and analysed by Choi et al. (2013) for λ = 0 . t = 0 at the start of the run, and stop therun when the central collapse reaches the final resolution of R fin ∼ − pc, in both Enzo and GADGET-3, in order toallow a meaningful comparison, because the SPH run would have difficulty reaching a higher resolution without invokingcomplicated procedures, such as SPH particle splitting.The DM in isolated models has isothermal profiles byconstruction, i.e., ρ DM ∝ R − . They remain nearly identi-cal to ICs. As the collapse proceeds and reaches the centralregion, the DM is dragged inward by the gas, and is com-pressed adiabatically, forming new cusps with slopes ∼ − ∼ − ρ ∝ R − (Larson 1969; Penston 1969). But significant differencesare obvious, as angular momentum is present and modifiesthe solution, and the DM dominates the gravitational po-tential everywhere during the initial stage.When collapsing gas reaches the centrifugal barrier, ashock forms, where the density jumps and the radial velocitydrops abruptly (top panels of Figure 2). The gas accumulatesbehind the shock and forms a disk which grows in mass. Asthe gas density in the inner disk surpasses that of the DM,the 2nd stage of the collapse ensues, and proceeds with an in-creasingly short timescale. Basically, most of the gas behindthe shock collapses and reaches the prescribed resolution of ∼ − pc (e.g., Choi et al. 2013), although here we com-pare the final distributions when R fin ∼ − pc has beenreached.The current AMR simulations with various λ follow ba-sically the same outline. The initial positions of shocks inall runs are located at ∼ λ , as the centrifugal barrier is reachedat progressively larger radii. Figure 3 displays the shock lo-cations as function of time, from their formation time tothe end of the simulations. Note, that for a specific λ , bothAMR and SPH shocks strictly follow each other. Moreover, c (cid:13)000 , 000–000 Yang Luo, Kentaro Nagamine, and Isaac Shlosman
Figure 2.
From top to bottom, radial profiles of gas density,radial velocity, tangential velocity, and mass accretion rate, atthe end of the isolated runs for models AMR-I5 (solid lines) andSPH-I5 (dashed lines). All axes are logarithmic. The density andaccretion rates have been averaged over spherical shells, and thevelocities averaged over cylindrical shells. consistently, the AMR shocks reach larger radii, as shown inFigure 3 for all λ . The explanation for this difference comesfrom the fact that the isolated AMR models collapse laterthan the isolated SPH models, as we discuss below, andshocks have more time to advance.As a first step, we compare the radial profiles of a gasdensity, radial and tangential velocities, and accretion ratesat the end of the runs, i.e., when the collapse has reached R fin (Fig. 2). All values are given after averaging in cylindri-cal or spherical shells and at the run end. The gas densityand accretion rates have been averaged over spherical shells,and the velocities averaged over cylindrical shells. Shownare the representative AMR-I5 and SPH-I5 runs. In AMR-I5 run, the gas goes through a strong accretion shock andforms a disk, ending the first stage of the collapse. By theend of the runs, the SPH-I5 shock is positioned at a radiusabout 3 times smaller than in the AMR-I5, which is con-firmed by Figure 3. It appears also somewhat weaker thanthe AMR shock. In all other respects, the gas density profilesare very similar, albeit it takes another decade in radius forthe density in SPH-I5 to catch up with the AMR-I5 density Figure 3.
Evolution of the radial position of shocks, r s , in iso-lated models in Enzo (solid lines) and GADGET-3 (dashed lines)runs for various λ . Note that both AMR and SPH runs with thesame λ exhibit the same r s ( t ). The maximal r s correspond to theend of the simulation. just inside the shock. The final density is ρ ∼ − g cm − at R fin = 10 − pc.The radial velocity profiles, v r , are very similar in theouter few hundred pc. In the preshock region, both v r in-crease toward the center, then fall down substantially toa deep minimum in the post-shock region, at the respec-tive positions of the shocks in the AMR and SPH runs. Atsmaller radii of r < − pc, the curves are very similar andincrease by an order of magnitude toward R fin . The tangen-tial velocities, v t , reach maximal values in the post-shockregion of both runs, as expected, then decrease toward thecenter. The behavior of v r and v t with r is associated witha variable degree of rotational support with radius, as wediscuss below.The mass accretion rate profiles show a similar be-havior, but ˙ M differs by a factor of 3 in the preshockregion, where ˙ M ∼ . M (cid:12) yr − for the AMR run and ∼ . M (cid:12) yr − for the SPH run. This rate drops to below ∼ . M (cid:12) yr − at R <
10 pc for both runs. The 2nd stageof the collapse inside ∼ . r , reaching ˙ M ∼ . M (cid:12) yr − at R fin . The small spike inthe SPH run at r ∼ − pc is caused by insufficient massresolution there.The projected gas volume density shown in Fig. 4 em-phasizes both similarities and differences, in comparisonwith the 1-D radial density profiles. The most significantdifference is that of the central disk size in the top frames— disk radius of ∼
20 pc in the AMR-I5, and only ∼ m = 2 (gaseous bar) is seen within the central 0.3 pc,and m = 3 mode is driving a triple spiral on scales > ∼ . c (cid:13) , 000–000 omparing AMR and SPH Approaches in Direct Collapse Modeling Figure 4.
Mass-weighted gas volume density in the x - y equatorial plane (top and middle frames) and the x - z plane (bottom) of isolatedmodels AMR-I5 (left frames) and SPH-I5 (right frames), at the end of the runs. Shown are the snapshots of decreasing spatial scales(physical): 50 pc (top) and 5 pc (middle and bottom) on each side, respectively. The color palette reflects the range from minimum tomaximum densities in each panel. The face-on and edge-on disks can be seen clearly in both runs.c (cid:13)000
Mass-weighted gas volume density in the x - y equatorial plane (top and middle frames) and the x - z plane (bottom) of isolatedmodels AMR-I5 (left frames) and SPH-I5 (right frames), at the end of the runs. Shown are the snapshots of decreasing spatial scales(physical): 50 pc (top) and 5 pc (middle and bottom) on each side, respectively. The color palette reflects the range from minimum tomaximum densities in each panel. The face-on and edge-on disks can be seen clearly in both runs.c (cid:13)000 , 000–000 Yang Luo, Kentaro Nagamine, and Isaac Shlosman
Figure 5.
Specific angular momentum profile in the gas as a function of cylindrical radius for AMR-I5 (left panel) and SPH-I5 (rightpanel). Circular angular momentum profile is shown for comparison (dashed line).
The SPH-I5 run displays comparable m = 2 and 3 modeson the same scales.Redistribution of the angular momentum is an impor-tant ingredient of gravitational collapse. We compare spe-cific angular momenta, j z , in the gas at the end of the simu-lations, in AMR-I5 and SPH-I5 (Fig. 5). Two important dif-ferences can be observed. Positions of radial shocks (Fig. 2)corresponds to the radii where the specific angular momentacurves come close to the circular angular momenta, j z ∼ j c .At smaller radii, j z moves away from j c gradually. Still, thereis a non-negligible rotational support at R fin , in both runs,albeit smaller for the SPH-I5 run. The bottom frames of thisFigure provide the quantitative measure of rotational sup-port in the form of j z /j c radial profiles. The AMR exhibits adisk at larger radii than the SPH run, because it takes moretime for the AMR model to collapse, and, therefore, the ac-cretion shock has more time to propagate outwards. In the2nd stage of the collapse, i.e., inside ∼ j z /j c ratiodeclines inwards, and showing a plateau for < ∼ − pc, atabout 0.5. They demonstrate that, in the isolated models,rotational support for the collapsing gas never falls below ∼
50% in the 2nd stage, in both runs.The T − ρ g relation is shown in Figure 6 at the endof the representative AMR-I5 and SPH-I5 runs. While theupper envelope of T is similar in both models, there ap-pears to be much more gas at lower T at low densities, ρ < − g cm − , in the AMR-I5 run. Given our refinementcriteria, this is an anticipated property of mesh codes whichare able to follow the low-density gas, while SPH codes tend to resolve better the higher density regions. On the otherhand, a larger amount of low T gas appears to be presentat high densities, ∼ − − − g cm − , in the SPH-I5run. We note that AMR run has larger mass accumulated(at ρ ∼ − g cm − ) around the shock, due to a largershock radius, as evident from Figure 3. Since we use identicalcooling packages for all runs, this difference can arise fromthe ability of numerical schemes to resolve the shocks andthe cooling of a post-shock gas. In the AMR-I5, the radialshock is strong and advances to larger radii correspondingto lower gas densities of ∼ − g cm − . The low T gas isformed from the post-shock gas which cools down more effi-ciently. On the other hand, the SPH-I5 shock is weaker andhappens at ∼ − g cm − , where particles at lower T al-ready are much more common. Hence the shocked particles’contribution is hardly visible on the phase diagram.Figure 7 compares various parameters which affect theresolution in AMR and SPH runs. The left panel showsevolution of maximal density in the gas, ρ max , with time.Clearly, the collapse in the isolated model SPH-I5 proceedsmuch faster than in AMR-I5 already during the 1st stage.At the onset of 2nd stage, the density increases slowly, ex-hibiting kind of a plateau, then accelerates, and at the endof this stage ρ max increases rapidly. Both runs display anidentical behavior at the end.The evolution of ρ max allows us to quantify the separa-tion of the two stages of the collapse (see also Choi et al.2013). The plateau between the stages provides an indepen-dent quantitative support to the existence of these stages. c (cid:13) , 000–000 omparing AMR and SPH Approaches in Direct Collapse Modeling Figure 6.
Phase diagram of temperature vs. density of gas in isolated runs AMR-I5 (left panel) and SPH-I5 (right panel) at the end ofeach run. The color represents the amount of mass in each pixel.
Figure 7.
Comparing various resolution parameters for isolated models AMR-I5 (solid lines) and SPH-I5 (dashed lines).
Left panel:
Evolution of maximum gas density ρ max as a function of time. Middle panel:
Evolution of the minimal gas smoothing length for theSPH and minimum cell size for the AMR run, η min , as functions of time. Right panel:
Radial profile of resolution given by the smoothinglength for SPH and the minimum cell size at each radius for AMR, at the end of the runs. The 1st stage of the collapse which startsat t = 0 is identical for both models. The 2nd stage of the collapse is triggered when the gas-to-DM density ratio in the inner few pcexceeds unity, i.e., ∼ . ∼ . η min (middle panel). They are also observed in the evolution of the minimumgas smoothing length for SPH and minimum cell size forAMR, η min , shown in the middle frame of Figure 7 — thisframe describes the same process of ρ max but in terms ofthe gas minimal smoothing length η min , which is the mini-mal interparticle distance Following Choi et al. (2013, 2015) and Shlosman et al. (2016), we choose to define the be-ginning of the 2nd stage when the gas-to-DM density ra-tio becomes larger than unity within the central few pc.This happens at t ∼ . ∼ . c (cid:13)000
Radial profile of resolution given by the smoothinglength for SPH and the minimum cell size at each radius for AMR, at the end of the runs. The 1st stage of the collapse which startsat t = 0 is identical for both models. The 2nd stage of the collapse is triggered when the gas-to-DM density ratio in the inner few pcexceeds unity, i.e., ∼ . ∼ . η min (middle panel). They are also observed in the evolution of the minimumgas smoothing length for SPH and minimum cell size forAMR, η min , shown in the middle frame of Figure 7 — thisframe describes the same process of ρ max but in terms ofthe gas minimal smoothing length η min , which is the mini-mal interparticle distance Following Choi et al. (2013, 2015) and Shlosman et al. (2016), we choose to define the be-ginning of the 2nd stage when the gas-to-DM density ra-tio becomes larger than unity within the central few pc.This happens at t ∼ . ∼ . c (cid:13)000 , 000–000 Yang Luo, Kentaro Nagamine, and Isaac Shlosman
Table 3.
DM halo virial parameters in AMR and SPH cosmo-logical models at the time of direct collapse, i.e., at z coll .Models z coll M vir [ h − M (cid:12) ] R vir [ h − pc]AMR-C1.5 17.2 2 . × . × . × . × . × . × η min ∼ few × − pc at ∼ . ∼ . r fin .The final radial profile of η min is shown in the rightpanel of Figure 7. Both AMR and SPH curves for η min aresimilar, but still the AMR one is systematically lower by afactor of 2 − For cosmological models, we stop the runs when the AMRand SPH resolutions, measured by η min , reach R fin ∼ − pc. The reason for this is that high resolution of iso-lated models is more difficult to reach in cosmological runs.Therefore, we reduce the resolution demands in both AMRand SPH cosmological models by a factor of 10 compared tothe isolated ones. The time is measured from the Big Bang.As the DM halos form and grow with time in the cos-mological models, the DM density profiles tend to the NFWshape (Navarro et al. 1997). When the collapse reaches thecentral regions, the baryon drag in the DM and form a cusp,similar to that in isolated models. Note, that our isolatedmodels follow the isothermal sphere density profiles and notthe NFW ones. We shall return to this issue in § λ ∼ .
042 —very similar to the representative isolated cases (Fig. 8). Wehave also measured the DM halo virial masses and virialradii for Enzo and GADGET cosmological models at thetime of direct collapse, z coll (Table 3). Overall, the AMRmodels collapse slightly earlier than the SPH ones. Conse-quently, the DM halos appear slightly less massive and theirvirial radii slightly smaller. So a pronounced trend existsbetween the AMR and SPH cosmological models.The top panel displays the final density profiles whichare remarkably similar for the two runs. While density at R fin ∼ − pc is nearly identical to that in the correspond-ing isolated runs, the overall shape of ρ ( R ) somewhat dif-fers. In the isolated runs, inside the radial shock, the slopeof ρ ( R ) was nearly constant, while in the cosmological runsit exhibits the trend of becoming shallower with R . Thischange happens at R ∼ . − r , but end up similar withina factor of 2 at the center. The tangential velocity is consis-tently higher in AMR at r (cid:38) . r (cid:46) . M . The curves differ by a factor ∼ − < ∼ r ∼ −
100 pc and stay about constant, ˙ M ∼ − (cid:12) yr − inward. The reason for this behavior is increased rotationalsupport they experience in this region, i.e., v r and v t havelocal minima and maxima there. On the other hand, theisolated models display growing ˙ M toward r fin and reach amaximum there which is greater than the cosmological caseby a factor of 2. The cosmological ˙ M is higher by a factorof a few everywhere except at the very center. Lastly, at thecenter, for R < ∼ . λ = 0 . R fin ∼ .
01 pc. The top panels show the large scale distribu-tion of the DM on scales of 100 h − kpc (comoving). Somedifferences are noticeable in these frames, and typically canbe attributed to individual DM halos. Indeed, we focus onthe representative DM halos in the second row within (ap-proximately) their virial radii, with the box size of 5 h − kpc(comoving). Again, while overall a nice degree of similarityexists between the two runs, differences are clear as well,especially if one focuses on the substructure on this scale.The bottom panels show the corresponding gas distributionwithin these halos, on scales of 200 h − pc (comoving). Whilethe AMR panel displays the characteristic filamentary struc-ture, the SPH panel shows a rather centrally-concentrated,elongated gas distribution, although the central density isthe same in both runs. The SPH density distribution dis-plays less small-scale structure in general. The AGORA col-laboration (Kim et al. 2014) has also reported difference insubstructure within the similar large-scale structure.Next, we show the radial distribution of the specific an-gular momentum profile at the end of the simulations, aswell as j z /j c ratio (Fig. 10). Clearly, both runs display muchless angular momentum support at all radii compared to theisolated models. On top of this, the SPH run shows even lesssupport than the AMR one, and proceeds basically in thefree-fall fashion until it has reached r ∼ . T − ρ mapsfor the above models in Figure 11. The overall feature of c (cid:13) , 000–000 omparing AMR and SPH Approaches in Direct Collapse Modeling Figure 9.
Projections of final DM density in x − y plane in cosmological models AMR-C4.2 (left column) and SPH-C4.2 (right column).The spatial scale of each panel decreases from top to bottom: DM in 100 kpc (comoving), DM in 5 kpc (comoving), and mass-weighted gasvolume density in 200 pc (comoving), respectively. The depth of each frame is 10 kpc (top), and 0.5 kpc (middle) in comoving coordinates.The color palette reflects from minimum to maximum projection densities in each panel. the phase diagrams are very similar between the two runs,with the exception of the mass involved at particular T and ρ values. At high densities, the AMR run shows a some-what broader distribution in T -range with very small masscontained in each grid-cell, which are completely absent inthe SPH run. This is understandable, as the two codes havea different resolution at these densities. Although we start from identical initial conditions in Enzo and GADGET-3,the refinement history differs between the codes in detail.We have tested how sensitive is the evolution with respectto the refinement condition in Enzo and discuss it below.To specify the characteristic times of gravitational col-lapse in the cosmological models, we follow the definitionfrom isolated models. The onset of the 1st stage of the col- c (cid:13)000
Projections of final DM density in x − y plane in cosmological models AMR-C4.2 (left column) and SPH-C4.2 (right column).The spatial scale of each panel decreases from top to bottom: DM in 100 kpc (comoving), DM in 5 kpc (comoving), and mass-weighted gasvolume density in 200 pc (comoving), respectively. The depth of each frame is 10 kpc (top), and 0.5 kpc (middle) in comoving coordinates.The color palette reflects from minimum to maximum projection densities in each panel. the phase diagrams are very similar between the two runs,with the exception of the mass involved at particular T and ρ values. At high densities, the AMR run shows a some-what broader distribution in T -range with very small masscontained in each grid-cell, which are completely absent inthe SPH run. This is understandable, as the two codes havea different resolution at these densities. Although we start from identical initial conditions in Enzo and GADGET-3,the refinement history differs between the codes in detail.We have tested how sensitive is the evolution with respectto the refinement condition in Enzo and discuss it below.To specify the characteristic times of gravitational col-lapse in the cosmological models, we follow the definitionfrom isolated models. The onset of the 1st stage of the col- c (cid:13)000 , 000–000 Yang Luo, Kentaro Nagamine, and Isaac Shlosman
Figure 10.
Specific angular momentum profile of gas as function of cylindrical radius for Enzo (left panel) and GADGET3 (right panel)in the cosmological runs AMR-C4.2 and SPH-C4.2. Circular angular momentum profile is shown for comparison (dashed line). lapse is inferred from the sudden rise in refinement level inEnzo and from η min ( t ) in GADGET. The 2nd stage is trig-gered when the gas-to-DM density ratio exceeds unity inboth Enzo and GADGET.The left panel of Figure 12 shows the time evolutionof ρ max in cosmological simulations. The time prior to t ∼
178 Myr corresponds to DM structure formation. It involvesthe initial expansion and collapse of DM shells, forming asmall DM halo, which we target in the zoom-in simulations.This halo grows via accretion of cold gas and DM — a pro-cess permeated by merger events. As the halo approachesthe critical mass of ∼ M (cid:12) , direct baryonic collapse fol-lows. The gas collapse starts around 178 Myr and proceedsrapidly (but slower than in the isolated models), being trig-gered nearly simultaneously in both runs.The middle frame of Figure 12 displays the evolution ofthe minimal smoothing length, η min , in the gas. Again, un-til the onset of the gravitational collapse, the AMR curvedisplays a better resolution compared to the SPH one, bya factor of 2 −
3, but both codes are resolving basically thesame densities at this time (see § t ∼
178 Myr, this difference in η min gets amplified andthe collapse proceeds increasingly faster in the AMR run.This is a reverse situation with respect to isolated models,and shows up explicitly in the radial profiles of η min . Weexpect that the run with a smaller η min will collapse first.We note that the duration of the baryonic collapse here, ∼ −
20 Myr, is substantially longer than in respective iso-lated models. The reason for this difference comes from thecollapse inside NFW halos which have a substantially lowercentral DM densities than in similar isothermal spheres usedby us for isolated models.The right frame of Figure 12 displays the radial profileof η min ( R ) at the end of cosmological simulations AMR-C4.2and SPH-C4.2. Note, the bifurcation in η min at around fewpc from the center of the DM halo. Whereas the outer haloregion is equally resolved with both runs, the inner one isbetter resolved with the AMR. This explains why the 2ndstage of the collapse proceeds faster with Enzo. As additional and a broad test for model comparison, wemeasure the collapse time in each model and contrast theEnzo-2.4 and GADGET-3 runs in isolated models. Effect ofgas rotation on gravitational collapse has been also modeledby Jappsen et al. (2009) in the context of Pop III star for-mation. However, they assumed a rigid spherical DM halowhich, therefore, cannot absorb angular momentum and/orproduce gravitational torques on the gas.To analyze the collapse time for isolated models, wemeasure the time from the start of the run, using modelswith a range in the spin parameter λ , for AMR and SPHruns. In other words, we define the collapse time, ∆ t c , as c (cid:13) , 000–000 omparing AMR and SPH Approaches in Direct Collapse Modeling Figure 11.
Phase diagram of temperature vs. density in AMR-C4.2 (left panel) and SPH-C4.2 (right panel) at the end of the cosmologicalmodel simulations. The color represents the mass in each pixels.
Figure 12.
Same as Figure 7, but for cosmological models AMR-C4.2 (solid lines) and SPH-C4.2 (dashed lines). Note, that the collapsestarts around t ∼
178 Myr for AMR and SPH runs, when ρ max and η min are slightly in favor of the AMR, compared to being identicalin the isolated models. η min ( R ) shows that the innermost resolution is higher in Enzo inside few central pc, while it is identical in theouter DM halo. the time from t = 0 till the collapse reaches R fin = 10 − pc( § t c , is given in Figure 13.For isolated models, two trends can be observed. First,both AMR and SPH runs exhibit a monotonic increase of∆ t c with λ , up to the measured spin of λ ∼ .
07 (Fig. 13).Second, the AMR runs display a longer collapse time thanthe SPH runs. As a corollary, the difference between the collapse time of AMR and SPH runs is increasing monoton-ically with λ .The first trend can be easily explained. With increasingDM halo spin parameter, the gas must overcome a centrifu-gal barrier positioned at progressively larger radii, becausethe initial conditions for the isolated runs require the gas tohave the same specific angular momentum j ( r ) distribution c (cid:13)000
07 (Fig. 13).Second, the AMR runs display a longer collapse time thanthe SPH runs. As a corollary, the difference between the collapse time of AMR and SPH runs is increasing monoton-ically with λ .The first trend can be easily explained. With increasingDM halo spin parameter, the gas must overcome a centrifu-gal barrier positioned at progressively larger radii, becausethe initial conditions for the isolated runs require the gas tohave the same specific angular momentum j ( r ) distribution c (cid:13)000 , 000–000 Yang Luo, Kentaro Nagamine, and Isaac Shlosman
Figure 8.
From top to bottom, radial profiles of gas density, ra-dial velocity, tangential velocity, and mass accretion rate, at theend of the cosmological runs for models AMR-C4.2 (solid lines)and SPH-C4.2 (dashed lines). All axes are logarithmic. The den-sity and accretion rate have been averaged over spherical shells,and velocities averaged over cylindrical shells. as the DM. The gas evolution is complicated by the intri-cacies of gravitational collapse: gas at small radii collapsesfirst and it has lower j than the gas at larger radii. So theinitial circularization determines the position of the accre-tion shock which forms at ∼ j circularizes at larger r progressively later, the shock moves out, and does it fasterfor larger λ . Hence, the size of the forming disk behind theaccretion shock is increasing with λ and with time (e.g.,Fig. 3).For the second stage of the collapse to ensue, the gasmust decouple from the background DM potential, i.e., itsdensity must increase above the DM density at small radii(e.g., § Figure 13.
Duration of gravitational collapse time for isolatedmodels, ∆ t c , in Enzo-2.4 and GADGET-3 runs as function ofthe halo spin parameter λ . The grey region represents ± σ errorsusing the Poisson statistics. lation of baryons is required when decoupling happens atlarger radii, explaining the increased ∆ t c for larger λ . In-deeed, the disk sizes and, therefore, masses increase with λ . Hence, this simple physics explains the relatively mod-est increase in ∆ t c with λ . Remarkably, both Enzo andGADGET-3 runs agree in their ∆ t c for λ = 0 .
01. But whatabout the increasing delay in the collapse time of the AMRruns with λ ? To understand this, one should look carefullyinto different resolution pattern of both codes. While ini-tial conditions are identical in this respect, the SPH grav-itational softening is constant with time, while that of theAMR code are refined as the collapse proceeds. This by itselfdoes not have a direct effect on the dynamics, but there arecaveats. Probably the most important issue is the possibleeffect that refinement has on the development on degree ofturbulence in the collapsing gas.In Figure 14a, we show evolution of η min for three repre-sentative λ for Enzo and GADGET-3. For λ = 0 . η min ( t )displays little difference between Enzo-2.4 and GADGET-3.For higher spin, we observe an increasing delay in the 2ndstage of the collapse, especially for Enzo runs. Checking theassociated increase in the refinement levels, the latter ap-pears to stagnate before the onset of the 2nd stage of thecollapse. In other words, η min , which reflects the best spatialresolution at each time, as mentioned above, “hesitates” todecrease for a longer time in Enzo compared to GADGET-3.Therefore, we observe a longer plateau for Enzo in Figure 14. In the previous section, we have measured the collapsetimescales for isolated models and contrasted them betweenEnzo-2.4 and GADGET-3 runs. Similar exercise in cosmo-logical models is more complicated, because it involves cos-mological evolution of DM halos, which depends on struc- c (cid:13) , 000–000 omparing AMR and SPH Approaches in Direct Collapse Modeling ture formation in the universe and results in halos having dif-ferent shapes and other properties. Hence, while a detailedcomparison of the collapse time for cosmological models isoutside the scope of this work, we do compare some numer-ical aspects of their evolution with various λ .We define the cosmological collapse timescale, ∆ t c byfollowing the evolution of the baryonic refinement level, themaximal hydrodynamic resolution, η min , and the maximaldensity, ρ max — all defined as in isolated models. The be-ginning of the direct collapse in the cosmological DM halos, t , is determined when these parameters exhibit a sharp in-crease or decrease, correspondingly. For example, η min de-creases sharply at t ∼
178 Myr, which is considered as theonset of the collapse. The end point of the collapse is taken t fin , when it reaches R fin , and ∆ t c = t fin − t . To detecta sharp increase/decrease in ρ max and η min , we follow theirtime derivatives after ∼
140 Myr and test all the occurrencesfor the increase/decrease by a factor of 2 or more.We repeat the λ -test discussed in the previous section,for cosmological models. Naively, one expects to observe thesame correlation between ∆ t c and λ shown in Table 4 andin Figure 13. However, as exhibited by Figure 14b, this isnot the case. For three representative cosmological modelsin AMR and SPH, the one-to-one correlation between ∆ t c and λ does not exist here. Why does this happen?A number of physical reasons in cosmological modelswould ‘mess up’ this correlation. First, the host DM halosin isolated models are axisymmetric. This means that col-lapsing gas has difficulty to lose angular momentum to DM,and loses it mostly to gas. On the other hand, the DM haloshapes in cosmological models differ profoundly from ax-isymmetric shapes. Therefore, the gas will lose its angularmomentum continuously and efficiently to the DM, whichwill move the gas away from the centrifugal barrier at eachradius, as in fact is shown in Figure 10.Secondly, because the direct collapse takes about 10 −
20 Myrs in cosmological models, the basic properties of DMhalos can change during this time. They can undergo minor,intermediate and major mergers. In other words, the hosthalos will have enough time to interact with the substructure— which does not happen in isolated models.These and additional factors are expected to affect thedirect collapse onset and duration in cosmological simula-tions in a kind of unpredicatble way, which we indeed ob-serve in our simulations.For isolated and cosmological SPH and AMR models,the free-fall time for the DM halos at the onset of direct col-lapse is t ff ∼
35 Myr. In comparison with the actual collapsetimescales given in Table 4, the collapse times are shorterthan the free-fall times. This underlines an important pointthat direct collapse involves only the inner part of the gasinitially residing in DM halos.
We have performed comparison simulations of gas collapsein DM halos in the context of direct collapse scenario toform SMBH seeds using AMR and SPH codes. This prob-lem involves large dynamic range, from kpc down to AUscales — such a comparison has not been attempted before.Using Enzo-2.4 AMR and GADGET-3 SPH codes, we ex-
Table 4.
Values of ∆ t c for isolated and cosmological modelsshown in Figure 14. The definitions of ∆ t c are given in §§ t c (Myr)AMR-I1 1.7AMR-I5 6.3AMR-I7 8.3SPH-I1 1.7SPH-I5 2.3SPH-I7 3.4AMR-C1.5 14.6AMR-C2.8 17.5AMR-C4.2 14.3SPH-C1.5 20.6SPH-C2.8 19.8SPH-C4.2 20.6 amined the systematic differences which buildup when thesenumerical schemes are invoked. We made significant effortsto match the initial setup of the two simulations. We ad-dressed the associated problems using isolated and cosmo-logical models of gravitational collapse in DM halos witha range of cosmological spin parameter λ . To simplify thistask, we ignored the effect of molecular hydrogen and UVbackground radiation, and assumed an optically-thin atomiccooling, as described in §§ • The isolated models generally follow similar evolution-ary path and agree with Choi et al. (2013), but exhibitmany subtle differences among themselves. Specifically, thecollapse proceeds in two stages, first reaching the angularmomentum barrier, going through an accretion shock, andforming a disk behind it. In the second stage, which beginsafter the gas density in the disk region surpasses that of theDM density, the gas decouples from the background grav-itational potential of the DM and experiences a runawaycollapse, dragging some of the DM inward in an adiabaticcompression. We find that, although the AMR and SPHmodels start with identical initial conditions, the pace of in-crease in spatial resolution differs in both codes. This leadsto a substantially earlier collapse in the SPH models for theisolated models. • Cosmological models do not exhibit ‘standing’ accre-tion shocks, unlike their isolated model counterparts, and,consequently, do not form rotationally-supported disks. Thereason for this is the continuous loss of angular momentumby the collapsing gas to the background triaxial potential ofthe host DM halos, contrary to isolated models which staylargely axisymmetric. Nevertheless, one can distinguish be-tween first and second stages of the collapse, when the gas-to-DM density ratio crosses unity, within the central few pc.The cosmological evolution of DM halos and baryons ap-pears to be similar among the AMR and SPH runs, beingonly slightly in favor of the AMR models in baryonic spa-tial and mass resolutions. But a pronounced trend exists be-tween the AMR and SPH cosmological models, with the for-mer exhibiting direct collapse at higher redshifts, and, con-sequently, the collapse happens in less massive and smallerDM halos. With the onset of baryonic gravitational collapsewithin DM halos, the pace of the resolution change differs c (cid:13)000
Values of ∆ t c for isolated and cosmological modelsshown in Figure 14. The definitions of ∆ t c are given in §§ t c (Myr)AMR-I1 1.7AMR-I5 6.3AMR-I7 8.3SPH-I1 1.7SPH-I5 2.3SPH-I7 3.4AMR-C1.5 14.6AMR-C2.8 17.5AMR-C4.2 14.3SPH-C1.5 20.6SPH-C2.8 19.8SPH-C4.2 20.6 amined the systematic differences which buildup when thesenumerical schemes are invoked. We made significant effortsto match the initial setup of the two simulations. We ad-dressed the associated problems using isolated and cosmo-logical models of gravitational collapse in DM halos witha range of cosmological spin parameter λ . To simplify thistask, we ignored the effect of molecular hydrogen and UVbackground radiation, and assumed an optically-thin atomiccooling, as described in §§ • The isolated models generally follow similar evolution-ary path and agree with Choi et al. (2013), but exhibitmany subtle differences among themselves. Specifically, thecollapse proceeds in two stages, first reaching the angularmomentum barrier, going through an accretion shock, andforming a disk behind it. In the second stage, which beginsafter the gas density in the disk region surpasses that of theDM density, the gas decouples from the background grav-itational potential of the DM and experiences a runawaycollapse, dragging some of the DM inward in an adiabaticcompression. We find that, although the AMR and SPHmodels start with identical initial conditions, the pace of in-crease in spatial resolution differs in both codes. This leadsto a substantially earlier collapse in the SPH models for theisolated models. • Cosmological models do not exhibit ‘standing’ accre-tion shocks, unlike their isolated model counterparts, and,consequently, do not form rotationally-supported disks. Thereason for this is the continuous loss of angular momentumby the collapsing gas to the background triaxial potential ofthe host DM halos, contrary to isolated models which staylargely axisymmetric. Nevertheless, one can distinguish be-tween first and second stages of the collapse, when the gas-to-DM density ratio crosses unity, within the central few pc.The cosmological evolution of DM halos and baryons ap-pears to be similar among the AMR and SPH runs, beingonly slightly in favor of the AMR models in baryonic spa-tial and mass resolutions. But a pronounced trend exists be-tween the AMR and SPH cosmological models, with the for-mer exhibiting direct collapse at higher redshifts, and, con-sequently, the collapse happens in less massive and smallerDM halos. With the onset of baryonic gravitational collapsewithin DM halos, the pace of the resolution change differs c (cid:13)000 , 000–000 Yang Luo, Kentaro Nagamine, and Isaac Shlosman
Figure 14.
Comparison of the evolution of ρ max ( t ) for (a) three isolated models with λ = 0 . , .
05 and 0.07 for AMR-I1, -I5 and -I7(solid lines) and SPH-I1, -I5 and -I7 (dashed lines) runs. (b)
Same for three cosmological models, AMR-C1.5, - C2.7 and -C4.2 (solidlines) and SPH-C1.5, -C2.7, and -C4.2 (dashed lines). The collapse timescales, ∆ t c , for these models are listed in Table 4. among the codes and the resolution increases faster in theAMR runs. Consequently, the trend in the time-lag is re-versed here with respect to the isolated models, and theAMR models collapse before the SPH models. Thus, a smalldifference in the pre-collapse resolution is amplified substan-tially during the collapse. • The collapse time ∆ t c increases with increasing λ in iso-lated models. Greater centrifugal forces divert the infallinggas to larger radii, i.e., circularization happens further awayfrom the rotation axis. This is translated to larger expansionvelocities of the accretion shocks and larger accretion disksizes. Shock in the AMR runs propagate further out than inthe SPH runs, because the final collapse time is later for theAMR isolated models. • For cosmological models, the baryonic collapse time ismuch longer than in isolated models, because the averagedensity in the central regions of the NFW halos is lower thanthat of the isothermal spheres by a factor of ∼ ∼
10. Additionalfactors discussed above affect ∆ t c , unlike in isolated models.In this work we have focused on numerical aspects ofgravitational collapse to the SMBH seeds, and so avoid dis-cussing the associated physics, except when unavoidable.The main reasons for dissimilar performance of the AMRand SPH codes lies in the mismatch between the pace ofincrease in spatial and mass resolutions when comparingperformance of the AMR and SPH codes. These differencesappear to be profound and not easily reconcilable.Both numerical codes under comparison here, AMREnzo-2.4 and SPH GADGET-3, follow collisionless cold DMand dissipative fluid component. Each of these componentsis important in order to understand evolution of gravita-tional collapse. The DM plays a crucial role in diluting thegravitational interactions within the gas, and, therefore, in-creasing its Jeans mass. Together with development of virial supersonic turbulence in the gas, it reduces the fragmenta-tion in the collapsing flow, and allows the collapse to proceedvia large dynamic range.Most important is the resolution of the gravitationalforce. Enzo uses a multi-grid PM method, and therefore, itsresolution is governed by twice the minimum cell size ( § ∼
10 pc with three initial nestedgrid levels. At the same time, GADGET has had a constantsoftening length of (cid:15) g = 10 − pc for the SPH particles and (cid:15) DM = 0 .
37 pc for the DM, from the beginning (see Table 2).Of course the particles are not as close as these distancesfrom the start of the run, but the forces are computed ac-curately with the tree method down to (cid:15) g and (cid:15) DM . On theother hand, Enzo gravity was allowed to refine, but it did notrefine fast enough to catch up with GADGET in the isolatedmodels, judging from the evolution of the hydrodynamicalresolution parameter, η min . For periods of time, Enzo refine-ment level has been constant, then changed abruptly. Thus,Enzo resolution always lagged behind GADGET in the iso-lated models. Due to these differences, Enzo was unable tofollow the steep density cusp of the DM isothermal spherein the center.The above situation is similar to the one found in thecontext of cosmological simulation by O’Shea et al. (2005), c (cid:13) , 000–000 omparing AMR and SPH Approaches in Direct Collapse Modeling where Enzo and GADGET codes have been compared forcosmological structure formation. In this case, Enzo had toinvest in ∼ root-grid in order to reproduce the sameDM halo mass function as the GADGET run with 128 par-ticles. When simply using the root-grid of 128 , Enzo wasunable to follow the growth of early density perturbationsdue to poorer force resolution, and it resulted in an under-estimate of halo substructures and low-mass halos at earlytimes. This means that Enzo AMR code requires signifi-cantly larger computing resources than GADGET SPH codein order to obtain similar structures on small scales, partic-ularly in the context of cosmological structure formation.In the present work, the context and the spatial scaleof the problem is very different from that of O’Shea et al.(2005). Enzo’s main advantage is the refinement methodwhich allows it to obtain a superior hydrodynamic or grav-itational force resolution, albeit toward the end of the cos-mological simulation and in a very small volume. Therefore,during the pre-collapse phase, when high resolution is notrequired, the difference in computational resources used byboth codes are comparable, within a factor of two, in theCPU time. On the other hand, during the direct collapsein isolated and cosmological models, GADGET requires 5– 10 times more CPU time than Enzo. Compared to Enzo,GADGET over-resolves in the initial stage of the collapse,and under-resolves in the final stages, unless SPH particlesplitting method is used or the SPH has a dynamic grav-itational softening, as in GIZMO (e.g., Heller & Shlosman1994; Hopkins 2013).The bottom line is that in our simulations of directbaryonic collapse in isolated and cosmological framework,GADGET requires greater investment in computational re-sources, but benefits from better hydrodynamic and gravita-tional force resolution in the computational box, while Enzoachieves a superior resolution in a very small volume.Note that our cosmological models of direct collapseevolve similarly in the AMR and SPH cases. The left andmiddle panels of Figure 12 show that direct collapse happensnearly simultaneously in AMR and SPH runs, as exhibitedby the behavior of ρ max and η min . But these two parametersdiffer by a factor of a few until the onset of the direct collapseat t ∼
178 Myr. As the collapse develops, we observe thatSPH lags behind AMR, unlike in isolated models, where theSPH models collapse first. The reason for this behavior canbe observed in the right panel of Figure 12. Here the AMRhas higher mass and spatial resolution than SPH in the innerfew pc at the final phase of the collapse.The delay in the collapse time can have a significant im-pact on the cosmological evolution of direct collapse, in thepresence of an external UV background. The latter can dis-sociate H molecules, which is considered to be a critical fac-tor for the success of this scenario (e.g., Omukai 2001; Latifet al. 2011; Inayoshi et al. 2014; Sugimura et al. 2014). If thecollapse is delayed significantly, it might appear favorably forthe direct collapse scenario because of the higher UVB inten-sity from nearby star formation, which should be increasingwith time at redshifts z ∼
20 to 10. In our current simula-tions, Enzo-2.4 and GADGET-3 give mixed results on thecosmological collapse timescales, as we have discussed above.This happens because the cosmological collapse depends onrealistic properties of DM halos, such as their shapes. In amore controlled environment of isolated models, GADGET- 3 exhibits earlier direct collapse than Enzo-2.4, due to thedifference in the early force, mass and hydrodynamic reso-lution between the two codes.In summary, we have performed a comparison betweenthe AMR code Enzo-2.4 and the SPH code GADGET-3 inthe framework of direct collapse to SMBH seeds, using bothisolated and fully cosmological models of the collapse. Wehave followed the model evolution into a strongly nonlinearregime, when small initial differences have been substan-tially amplified. Although the overall model evolution hasbeen found similar at large, substantial differences also ex-ist at the end of the runs. We find that the main cause forthese differences lies in different evolutionary pace of massand spatial resolution with time and space. This results indifferent abilities of the numerical schemes to resolve hydro-dynamics and gravitational interactions, such as shocks andangular momentum transfer by gravitational torques.
ACKNOWLEDGMENTS
We thank the Enzo and YT support team for help. Allanalysis has been conducted using YT (Turk et al. 2011,http://yt-project.org/). We have used the Grackle chem-istry and cooling library (Bryan et al. 2014b; Kim et al.2014, https://grackle.readthedocs.org/). We are grateful toVolker Springel for providing us with the original versionof GADGET-3, and to Jun-Hwan Choi and Long Do Caofor their help in the early phase of this project. K.N. ac-knowledges the partial support by JSPS KAKENHI GrantNumber 26247022. I.S. acknowledges partial support fromSTScI grant AR-12639.01-A. I.S. and K.N. are grateful tothe support from the International Joint Research Promo-tion Program at Osaka University. Support for HST/STScIAR-12639.01-A was provided by NASA through a grantfrom the STScI, which is operated by the AURA, Inc., underNASA contract NAS5-26555. Numerical simulations were inpart carried out on XC30 at the Center for ComputationalAstrophysics, National Astronomical Observatory of Japan,as well as the VCC at the Cybermedia Center at Osaka Uni-versity.
REFERENCES
Abel T., Anninos P., Zhang Y., Norman M. L., 1997, NewAstronomy, 2, 181Abel T., Bryan G. L., Norman M. L., 2002, Science, 295,93Agertz O., Moore B., Stadel J., Potter D., Miniati F., ReadJ., Mayer L., Gawryszczak A., Kravtsov A., Nordlund˚A., Pearce F., Quilis V., Rudd D., Springel V., Stone J.,Tasker E., Teyssier R., Wadsley J., Walder R., 2007, MN-RAS, 380, 963Anninos P., Zhang Y., Abel T., Norman M. L., 1997, NewAstronomy, 2, 209Ascasibar Y., Yepes G., M¨uller V., Gottl¨ober S., 2003, MN-RAS, 346, 731Balsara D. S., 1995, Journal of Computational Physics, 121,357Barnes J., Hut P., 1986, Nature, 324, 446Bate M. R., Burkert A., 1997, MNRAS, 288, 1060 c (cid:13)000
Abel T., Anninos P., Zhang Y., Norman M. L., 1997, NewAstronomy, 2, 181Abel T., Bryan G. L., Norman M. L., 2002, Science, 295,93Agertz O., Moore B., Stadel J., Potter D., Miniati F., ReadJ., Mayer L., Gawryszczak A., Kravtsov A., Nordlund˚A., Pearce F., Quilis V., Rudd D., Springel V., Stone J.,Tasker E., Teyssier R., Wadsley J., Walder R., 2007, MN-RAS, 380, 963Anninos P., Zhang Y., Abel T., Norman M. L., 1997, NewAstronomy, 2, 209Ascasibar Y., Yepes G., M¨uller V., Gottl¨ober S., 2003, MN-RAS, 346, 731Balsara D. S., 1995, Journal of Computational Physics, 121,357Barnes J., Hut P., 1986, Nature, 324, 446Bate M. R., Burkert A., 1997, MNRAS, 288, 1060 c (cid:13)000 , 000–000 Yang Luo, Kentaro Nagamine, and Isaac Shlosman
Begelman M. C., Rees M. J., 1978, MNRAS, 185, 847Begelman M. C., Rossi E. M., Armitage P. J., 2008, MN-RAS, 387, 1649Begelman M. C., Shlosman I., 2009, ApJL, 702, L5Begelman M. C., Volonteri M., Rees M. J., 2006, MNRAS,370, 289Berger M. J., Colella P., 1989, Journal of ComputationalPhysics, 82, 64Bromm V., 2013, Reports on Progress in Physics, 76,112901Bromm V., Larson R. B., 2004, ARA&A, 42, 79Bromm V., Loeb A., 2003, ApJ, 596, 34Bryan G. L., Norman M. L., 1997, in Astronomical Societyof the Pacific Conference Series, Vol. 123, ComputationalAstrophysics; 12th Kingston Meeting on Theoretical As-trophysics, Clarke D. A., West M. J., eds., p. 363Bryan G. L., Norman M. L., O’Shea B. W., Abel T., WiseJ. H., Turk M. J., Reynolds D. R., Collins D. C., WangP., Skillman S. W., Smith B., Harkness R. P., Bordner J.,Kim J.-h., Kuhlen M., Xu H., Goldbaum N., Hummels C.,Kritsuk A. G., Tasker E., Skory S., Simpson C. M., HahnO., Oishi J. S., So G. C., Zhao F., Cen R., Li Y., EnzoCollaboration, 2014a, ApJS, 211, 19—, 2014b, ApJS, 211, 19Bryan G. L., Norman M. L., Stone J. M., Cen R., OstrikerJ. P., 1995, Computer Physics Communications, 89, 149Bullock J. S., Dekel A., Kolatt T. S., Kravtsov A. V.,Klypin A. A., Porciani C., Primack J. R., 2001, ApJ, 555,240Choi J.-H., Shlosman I., Begelman M. C., 2013, ApJ, 774,149—, 2015, MNRAS, 450, 4411Colella P., Woodward P. R., 1984, Journal of Computa-tional Physics, 54, 174Cullen L., Dehnen W., 2010, MNRAS, 408, 669Eisenstein D. J., Hut P., 1998, ApJ, 498, 137Ferland G. J., Porter R. L., van Hoof P. A. M., WilliamsR. J. R., Abel N. P., Lykins M. L., Shaw G., Henney W. J.,Stancil P. C., 2013, Rev. Mex. de Astronom. y Astrofs.,49, 137Fraser M., Casey A. R., Gilmore G., Heger A., Chan C.,2015, ArXiv:1511.03428Frenk C. S., White S. D. M., Bode P., Bond J. R., BryanG. L., Cen R., Couchman H. M. P., Evrard A. E., GnedinN., Jenkins A., Khokhlov A. M., Klypin A., Navarro J. F.,Norman M. L., Ostriker J. P., Owen J. M., Pearce F. R.,Pen U.-L., Steinmetz M., Thomas P. A., Villumsen J. V.,Wadsley J. W., Warren M. S., Xu G., Yepes G., 1999,ApJ, 525, 554Haehnelt M. G., Rees M. J., 1993, MNRAS, 263, 168Hahn O., Abel T., 2011, MNRAS, 415, 2101Heller C. H., Shlosman I., 1994, ApJ, 424, 84Hirano S., Hosokawa T., Yoshida N., Umeda H., OmukaiK., Chiaki G., Yorke H. W., 2014, ApJ, 781, 60Hopkins P. F., 2013, MNRAS, 428, 2840Hosokawa T., Omukai K., Yoshida N., Yorke H. W., 2011,Science, 334, 1250Inayoshi K., Omukai K., Tasker E., 2014, MNRAS, 445,L109Jappsen A.-K., Mac Low M.-M., Glover S. C. O., KlessenR. S., Kitsionas S., 2009, ApJ, 694, 1161Johnson J. L., Khochfar S., Greif T. H., Durier F., 2011, MNRAS, 410, 919Kim J.-H., Abel T., Agertz O., Bryan G. L., CeverinoD., Christensen C., Conroy C., Dekel A., Gnedin N. Y.,Goldbaum N. J., Guedes J., Hahn O., Hobbs A., Hop-kins P. F., Hummels C. B., Iannuzzi F., Keres D., KlypinA., Kravtsov A. V., Krumholz M. R., Kuhlen M., Leit-ner S. N., Madau P., Mayer L., Moody C. E., NagamineK., Norman M. L., Onorbe J., O’Shea B. W., PillepichA., Primack J. R., Quinn T., Read J. I., Robertson B. E.,Rocha M., Rudd D. H., Shen S., Smith B. D., Szalay A. S.,Teyssier R., Thompson R., Todoroki K., Turk M. J., Wad-sley J. W., Wise J. H., Zolotov A., AGORA Collaboration,2014, ApJS, 210, 14Komatsu E., Dunkley J., Nolta M. R., Bennett C. L., GoldB., Hinshaw G., Jarosik N., Larson D., Limon M., Page L.,Spergel D. N., Halpern M., Hill R. S., Kogut A., MeyerS. S., Tucker G. S., Weiland J. L., Wollack E., WrightE. L., 2009, ApJS, 180, 330Koushiappas S. M., Bullock J. S., Dekel A., 2004, MNRAS,354, 292Kravtsov A. V., Klypin A. A., Khokhlov A. M., 1997, ApJS,111, 73Larson R. B., 1969, MNRAS, 145, 405Latif M. A., Schleicher D. R. G., Schmidt W., Niemeyer J.,2013, MNRAS, 433, 1607Latif M. A., Zaroubi S., Spaans M., 2011, MNRAS, 411,1659Loeb A., Rasio F. A., 1994, ApJ, 432, 52Long S., Shlosman I., Heller C., 2014, ApJL, 783, L18Madau P., Rees M. J., 2001, ApJL, 551, L27Mayer L., Kazantzidis S., Escala A., Callegari S., 2010,Nature, 466, 1082Milosavljevi´c M., Bromm V., Couch S. M., Oh S. P., 2009,ApJ, 698, 766Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490,493Norman M. L., Bryan G. L., 1999, in Astrophysics andSpace Science Library, Vol. 240, Numerical Astrophysics,Miyama S. M., Tomisaka K., Hanawa T., eds., p. 19Omukai K., 2001, ApJ, 546, 635O’Shea B. W., Nagamine K., Springel V., Hernquist L.,Norman M. L., 2005, ApJS, 160, 1Penston M. V., 1969, MNRAS, 144, 425Price D. J., Federrath C., 2010, MNRAS, 406, 1659Prieto J., Jimenez R., Haiman Z., 2013, MNRAS, 436, 2301Saitoh T. R., Makino J., 2013, ApJ, 768, 44Schleicher D. R. G., Spaans M., Glover S. C. O., 2010,ApJL, 712, L69Shapiro S. L., Teukolsky S. A., 1985, ApJ, 298, 58Shlosman I., Choi J.-H., Begelman M. C., Nagamine K.,2016, MNRAS, 456, 500Springel V., 2005, MNRAS, 364, 1105Springel V., Hernquist L., 2003, MNRAS, 339, 289Sugimura K., Omukai K., Inoue A. K., 2014, MNRAS, 445,544Tasker E. J., Brunino R., Mitchell N. L., Michielsen D.,Hopton S., Pearce F. R., Bryan G. L., Theuns T., 2008,MNRAS, 390, 1267Truelove J. K., Klein R. I., McKee C. F., Holliman II J. H.,Howell L. H., Greenough J. A., 1997, ApJL, 489, L179Turk M. J., Abel T., O’Shea B., 2009, Science, 325, 601Turk M. J., Smith B. D., Oishi J. S., Skory S., Skillman c (cid:13) , 000–000 omparing AMR and SPH Approaches in Direct Collapse Modeling S. W., Abel T., Norman M. L., 2011, ApJS, 192, 9Vogelsberger M., Sijacki D., Kereˇs D., Springel V., Hern-quist L., 2012, MNRAS, 425, 3024Volonteri M., Rees M. J., 2005, ApJ, 633, 624—, 2006, ApJ, 650, 669Wise J. H., Abel T., Turk M. J., Norman M. L., SmithB. D., 2012, MNRAS, 427, 311Zel’dovich Y. B., Podurets M. A., 1965, Astron.Zh., 42, 963 c (cid:13)000