Formation of an embryonic supermassive star in the first galaxy
aa r X i v : . [ a s t r o - ph . GA ] S e p Mon. Not. R. Astron. Soc. , 000–000 (0000) Printed 17 October 2018 (MN L A TEX style file v2.2)
Formation of an embryonic supermassive star in the firstgalaxy
Kohei Inayoshi , ⋆ , Kazuyuki Omukai and Elizabeth Tasker , Department of Physics, Graduate School of Science, Kyoto University, Kyoto 606-8502, Japan Department of Astronomy, Columbia University, 550 West 120th Street, New York, NY 10027, USA Astronomical Institute, Tohoku University, Sendai 980-8578, Japan Department of Physics, Faculty of Science, Hokkaido University, Sapporo 060-0810, Japan
17 October 2018
ABSTRACT
We studied the gravitational collapse of a warm ( ∼ & M ⊙ ) using a three-dimensional hydrodynamical simulation including all the relevant cooling processes ofboth H and H, which can potentially induce cloud fragmentation. This is the firstsimulation of this kind to resolve protostar formation. We find that from a weaklyturbulent initial condition, the cloud undergoes runaway collapse without a majorepisode of fragmentation. Although the H fraction jumps by a large factor via thethree-body reaction at ∼ − g cm − , its cooling remains inefficient due to theoptical thickness, and the temperature remains & ∼ − g cm − , a hydrostaticprotostar with ≃ . ⊙ is formed. The protostar grows to the mass ≃ ⊙ and theradius ≃ ∼ ∼ ⊙ yr − , the protostar is expected to turn into an SMS withinits lifetime, eventually collapsing to a seed for the supermassive black hole observedin the early Universe at z ∼ Key words: cosmology: theory – dark ages, reionzation, first stars – stars: formation– quasars: supermassive black holes
Recent observations reveal the existence of supermassiveblack holes (SMBHs) with masses & M ⊙ as early as red-shift z & ∼
100 M ⊙ )exceeds the then Hubble time. As a solution to this co-nundrum, the formation from massive seed BHs originatingin the collapse of supermassive stars (SMSs; & M ⊙ )has been suggested (e.g., Begelman et al. 2006). Seed BHsformed in this way are expected to grow to & M ⊙ by z ∼ & K, provid-ing that H cooling is prohibited throughout the protostel-lar collapse. Without this latter constraint, the gas wouldrapidly cool via H and fragment into smaller pieces. H hasto therefore be dissociated e.g., via the photodissociation ⋆ E-mail: [email protected] by far-ultraviolet (FUV) radiation from nearby star-forminggalaxies (Omukai 2001; Bromm & Loeb 2003; Shang et al.2010; Inayoshi & Omukai 2011) or the collisional dissocia-tion in the shocked gas (Inayoshi & Omukai 2012). This inplace, the cloud can collapse almost isothermally at ∼ & . ⊙ yr − . Undersuch a high accretion rate, its growth is not hindered eitherby strong radiative feedback (Hosokawa et al. 2012, 2013) orby mass-loss due to stellar pulsations (Inayoshi et al. 2013).So far, however, most studies (e.g., Regan & Haehnelt2009; Latif et al. 2013b; Choi et al. 2013) have utilized sev-eral simplifying assumptions in studying the fragmentationprocess, e.g., turning off the H cooling, adopting optically-thin treatment of Ly α cooling, or using insufficient chemicalnetworks that neglect the H formation. Yet the efficiencyof fragmentation depends strongly on the thermal evolutiondetermined by the cooling processes and chemical reactions. c (cid:13) K. Inayoshi, K. Omukai and E. Tasker
As a result, it still remains unresolved whether the cloudfragments during the isothermal collapse at ∼ Letter , we use a three-dimensional hydrodynamical simulation to study the grav-itational collapse of a turbulent primordial-gas cloud withmass & M ⊙ . We follow the evolution until the formationof the protostar, including all relevant processes requiredfor examining possible fragmentation during this collapse.Among these, the H -line cooling is of primary importancedue to its ability to rapidly drop the temperature via thermalinstability if the gravitational collapse is delayed, a processpossible due to turbulence generated during the virializationof the halo. If the thermal instability occurs, the cloud canfragment into many smaller mass clumps instead of form-ing a single SMS. We therefore simulate the collapse to de-termine the likelihood of the outcome being a monolithiccollapse to a single star or fragmentation into a binary ormultiple member system. We performed a three-dimensional hydrodynamical simula-tion of the gravitational collapse of a primordial-gas cloudusing the adaptive mesh refinement code,
ENZO (Bryan et al.2014). Our main purpose is to investigate the gas dynam-ics over a wide range of the densities (10 − . ρ . − gcm − ). The cloud initially has a spherically symmetric den-sity profile enhanced by a factor f (= 1 .
6) above the criticalBonnor-Ebert (BE) distribution, an isothermal sphere em-bedded in a pressurized medium and supported in marginalhydrostatic equilibrium against gravitational collapse. Ac-cording to cosmological simulations (e.g., Wise et al. 2008),at the center of a first galaxy with virial temperature & K, forming in an environment where the H formation issuppressed, a warm ( T ∼ ∼ M ⊙ becomes gravitationally unstable at ρ ∼ − g cm − andcollapses. Based on this, we set the central density and tem-perature of the cloud to ρ c = 1 . × − g cm − and T = 8000 K, giving a mass and radius of 1 . × M ⊙ and 10 . is collisionally dissociated for ρ & − g cm − and T & andrefinement is controlled by insisting that one Jeans lengthis resolved by at least 64 grid cells (e.g., Turk et al. 2012).Under this condition, the simulation uses 23 out of the al-lowed 25 refinement levels, ensuring we are resolved by theabove criteria at all times and giving a limiting resolutionof . . ∼ − g cm − , the turbulence isstill subsonic in the cloud. To consider the density and veloc-ity perturbations due to the turbulence, we initially imposea subsonic velocity field (the root mean square of the veloc-ity is set to 0 . c s ) with power spectrum P ( k ) ∝ k − , whichcorresponds to the so-called Larson’s law for the contempo-rary star-forming regions (Larson 1981). To ensure that the -13-14-15-16-17-18 l og ρ ( g / c m ) -11-12-13-14-15-16-17 l og ρ ( g / c m ) -7-8-9-10-11-12-7-8-9-10-11 l og ρ ( g / c m ) l og ρ ( g / c m ) Figure 1.
Density distribution in the plane through the densitypeak for four spatial scales: from top-left, clockwise: the large-scale gas distribution ( ∼ − free-bound continuum cooling ( ∼ . ∼
100 AU), and the final protostar ( ∼
10 AU). turbulence is adequately resolved, we select the maximumk-mode value of 1/10 of the number of cells across the cloud.We consider the non-equilibrium primordial chemistryof 9 species (H, H , e − , H + , H +2 , H − , He, He + , and He ++ )and 13 hydrogen reactions selected to reproduce the cor-rect thermal/chemical evolution of the warm atomic-coolingcloud (reactions 3, 4, 7 −
10, 12, 15 −
18, 28, and 32 in table 2of Omukai 2001). We adopt the reaction rate coefficients up-dated by the following studies: 7 −
10 (Coppola et al. 2011),15 (Martin et al. 1996), 17 (Stibbe & Tennyson 1999), and28 (Ferland et al. 1992). The four helium reactions originallyincluded in
ENZO are also present, although they are notrelevant in our calculation. We initially assume a uniformdistribution of ionization degree with 10 − and H molecu-lar fraction with 10 − , respectively (e.g., Shang et al. 2010).At high density, the chemical reactions proceed faster thanthe cloud collapse and chemical equilibrium is achieved. Tosmoothly connect the non-equilibrium chemistry to that ofequilibrium, we solve the chemical network including boththe forward and reverse reactions for dominant processes. Tosolve the chemistry equations, we employ the piecewise ex-act solution method (Inoue & Inutsuka 2008) instead of theoriginal ENZO solver, which cannot follow the chemical evolu-tion with high enough density to reach the chemical equilib-rium. For the radiative cooling, we consider atomic cooling(H Ly α , two photon emission, and H − free-bound, free-freeemission) and H cooling (rovibrational line and collision-induced emission). We also include the suppression of thecooling rate in the optically thick case by using the opticaldepth estimated as ρκL c (e.g., Omukai 2001; Shang et al.2010), where κ includes the H -line opacity and the Rosse-land mean opacity considering the H Rayleigh scattering,the H collision-induced absorption, and the H − bound-freeand free-free absorption, and L c the size of the central core,which is approximately given by the Jeans length for thespherically symmetric cloud in the runaway collapse. Finally,note that we do not include the heating/cooling associatedwith the chemical reactions because their effect is negligibleduring the thermal evolution of the atomic-cooling clouds. c (cid:13) , 000–000 ormation of a supermassive star -6-8-10-12-14-16-18-20 l og ρ ( g / c m ) ∝ r (cid:16) log r (cm) T ( K )
11 12 13 14 15 16 17 18 19 20 (a)(b) -2-4-6-8-10 l og H fr ac ti on (c) Figure 2.
Mass-weighted radial profiles at different evolutionarystages of (a) mass density, (b) temperature, and (c) H fraction.The time sequences are indicated by numbers: (1) 8 . × yrafter the initial state of our simulation. (2) 9 . × yr after(1): the main coolants are the Ly α and two-photon emissions. (3)1 . × yr after (2): dominant cooling process shifts to the H − free-bound emission in the central core. (4) 1 . × yr after (3):H formation via the three-body reaction becomes active at thecentre. (5) 2 . × yr after (4): the cloud becomes optically thickto dominant H lines. (6) 1 . × yr after (5): the cloud becomesoptically thick to the continuum opacities and a hydrostatic coreis formed at the centre. (7) 1 . Fig. 1 shows the density distribution at the end of the sim-ulation, where the central density reaches ∼ − g cm − ,for four different spatial scales; from the top-left clockwise,large-scale gas distribution ( ∼ ∼ . ∼
100 AU region, and the proto-star formed at the center ( ∼
10 AU). The central portion ofthe cloud undergoes the runaway collapse. The turbulenceforms filamentary structures that channel material into thecentral region ( ρ ∼ − g cm − ), feeding the protostar. Theleft-bottom panel presents the density distribution aroundthe protostar. At the end of this simulation, the protostellarmass reaches ≃ ⊙ and its radius ≃ T e m p e r a t u r e ( K ) C e ll M a ss ( M s un ) (b) H fr ac ti on -2 -4 -6 -8 -10 -12 -20 -18 -16 -14 -12 -10 -8 -6 C e ll M a ss ( M s un ) Density (g/cm ) (a)
Figure 3.
Phase diagrams showing the distribution of (a)density-temperature and (b) density-H fraction of the collaps-ing cloud at the end of the simulation. The colors represent thetotal mass at the respective density and temperature. Fig. 2 shows the evolution of mass-weighted radial pro-files of (a) density, (b) temperature, and (c) H fraction.During collapse, the density profile obeys the self-similarsolution, which consists of the central core with flat densitydistribution and envelope with the ρ ∝ r − law (e.g., Larson1969). The central core collapses almost isothermally until ∼ − g cm − keeping the temperature at ∼ ρ < − g cm − , the cooling ismainly via the H Ly α and two-photon emission. At higherdensity, the dominant cooling process shifts to the H − free-bound emission (H + e − → H − + γ ). For & − g cm − ,photons from the H − free-bound emission are self-absorbed,as well as Rayleigh scattered by H. The gas cools further bythe H − free-free emission (H + e − → H + e − + γ ) until & − g cm − . Finally, the cloud becomes opaque to allthose continuum opacities at ∼ − g cm − and a hydro-static core, i.e., a protostar, with its mass ≃ . ⊙ andradius ≃ ∼ ∼ ⊙ ,the temperature inside the protostar adiabatically increasesto ∼ K.The mass-weighted H fraction initially approaches theequilibrium value ( ∼ − ), where the formation throughthe electron-catalyzed reaction (H + e − → H − + γ ; H − + H → H + H) and the collisional dissociation (H + H → ρ & − g cm − , the H fraction jumpsup to ∼ . → H + H) inthe inner region ( r . AU). However, neither the H linenor collision-induced-emission (CIE) cooling plays a signif-icant role in the thermal evolution: H lines are opticallythick for ρ & − g cm − and other continuum cooling c (cid:13) , 000–000 K. Inayoshi, K. Omukai and E. Tasker
30 km/s
200 AU
20 km/s T e m p e r a t u r e ( K ) Figure 4.
Temperature (color) and velocity field (arrows) in theplane through the density peak for two spatial scales. Cold re-gions are formed by the chemo-thermal instabilities due to theH formation by the electron-catalyzed reaction (left) and by thethree-body reaction (right). is more important than the H CIE cooling. After the pro-tostar formation, the H is dissociated inside owing to thehigh temperature.Fig. 3 presents the phase diagrams showing the distribu-tion of (a) temperature and (b) H fraction as a function ofthe density at the end of the simulation. The cloud consistsof two thermal phases of the gas, i.e., hot ( ∼ several 10 K)and cold ( ∼ K) components. As seen in Figs. 2(b) and3(a), most of the collapsing gas resides in the hot compo-nent, which ultimately forms a protostar at the center. Notethat since the density profile follows the self-similar formduring the runaway collapse and thus the radial positionhas one-to-one correspondence with the density (Fig. 2a),the density-temperature distribution of the hot componentis just a reflection of the temperature profile (Fig. 2b). TheH fraction in the hot component remains almost constantat ∼ − up to 10 − g cm − and then increases almostproportionally to the density for higher density by the three-body reaction until finally dissociated at & − g cm − asa result of the protostar formation.Meanwhile, the cold component exists over the widedensity range, 10 − . ρ . − g cm − . This gas is pro-duced by the thermal instability induced by the combina-tion of the adiabatic cooling due to the turbulent expansionand the subsequent H cooling. Once the temperature de-creases via adiabatic cooling, the H dissociation becomesinefficient, enhancing the H fraction and its cooling rate,causing the temperature to plummet. This process is knownas the chemo-thermal instability associated with the H for-mation/dissociation (Yoshii & Sabano 1979; Silk 1983).Fig. 4 presents the temperature distribution and veloc-ity fields for two different scales. In both panels, the coex-istence of the cold and hot components is clearly visible.Turbulence establishes a complex structure of interactingshocks and stagnation points. The cold components in thetwo scales are produced by the thermal instabilities due tothe H formation through the electron-catalyzed reaction( ∼ . ∼
200 AU), see alsoFig. 3(b). The cold components are not massive enough tobe gravitationally bound and have no influence on the evo-lution of the central collapsing region.Fig. 5 presents the profiles of the radial and tangentialvelocities. Also shown for comparison is half of the Kep-lerian velocity and the sound speed. Both the radial and V e l o c it y ( k m / s ) v rad v tan v kep /2 c s log r (cm)
11 12 13 14 15 16 17 18 19 20
Figure 5.
Profiles of the radial (solid) and tangential velocities(long-dashed) at the end of the simulation. For comparison, a halfof the Keplerian velocity (short-dashed) and sound speed (dotted)are also shown. All the quantities are spherically averaged. M ( M s un / y r) log M r (M sun ) -3 -2 -1 0 1 2 3 4 5 6-4 * ・ Figure 6.
Profile of the mass infall rate ( ˙ M ∗ = − πρr v rad )as a function of the enclosed mass at the end of the simulation.The horizontal line indicates the critical value 10 − M ⊙ yr − ,above which the protostar swells to supergiant and the radiativefeedback is strongly suppressed. (Hosokawa et al. 2012). tangential flows become supersonic with the Mach numberof 2 −
3. At the surface of the protostar, the radial flow isabruptly brought to a halt. In the accreting envelope, thetangential velocity is as large as half the Keplerian velocity,in accordance with previous studies of Pop III star forma-tion (e.g., Abel et al. 2002; Yoshida et al. 2008). It is knownthat the cloud can contract in the runaway fashion evenwith conserved angular momentum as long as the temper-ature does not increase with the density (e.g., Narita et al.1984; Saigo & Hanawa 1998). We note that the turbulentvelocity is accelerated to the well-known universal value of ∼ . v Kep soon after the gravitational collapse starts evenwith initially weak turbulence. Thus, the result seems un-likely to depend on the initial turbulent velocity. After proto-star formation, on the other hand, materials initially locatedin the outer radius and thus with higher specific angular mo-mentum begin to fall in, and the centrifugal radius increasesoutwards (Saigo et al. 2000). In our simulation, however, therotationally supported disk has not yet appeared becausethe centrifugal radius ( . . ∼ c (cid:13) , 000–000 ormation of a supermassive star the temporal evolution of the accretion rate after the pro-tostar formation. Note that the value inside the protostar( M r . ⊙ ) is not equivalent to the accretion rate. Theinfall rate becomes almost constant for the flat temperatureprofile because in the self-similar solution, the flat tempera-ture profile is proportional to c s /G , which depends only onthe temperature of the accreting envelope. The typical valueis as high as ∼ ⊙ yr − , which is consistent with the pre-vious simulations starting from the cosmological initial con-dition (e.g., Latif et al. 2013b). This infall rate is larger than20 c s /G for T = 8000 K and similar to the value found forthe runaway collapse starting from an initial condition notso far from the hydrostatic equilibrium. (Foster & Chevalier1993). The protostar is expected to grow via such rapidaccretion to an SMS within its lifetime ∼ yr. Whenthe stellar mass exceeds ∼ M ⊙ , the SMS is expectedto collapse to a BH by the general relativistic instability(Chandrasekhar 1964; Hosokawa et al. 2013). We simulated the collapse of a massive ( & M ⊙ ) andwarm ( ∼ ∼ . ⊙ is formed when the centralpart becomes optically thick to the continuum radiation at ρ & − g cm − and grows to the mass ≃ ⊙ and radius ≃ ∼ ⊙ yr − . Where the accretion rate is higherthan 10 − M ⊙ yr − (dashed line in Fig. 6), the protostar isknown to develop a giant-like structure with a bloated stel-lar envelope and contracting central core (Hosokawa et al.2012, 2013) and grows while avoiding the significant mass-loss due to the stellar pulsation (Inayoshi et al. 2013). Sincethe effective temperature of such a supergiant protostaris . K, the UV feedback is unlikely to prevent themass accretion on to the star. However, recent simulationby Regan et al. (2014) suggests the possibility of disk frag-mentation around the protostar. Further simulations of thedisk and the fragments with proper treatment of the H2-line cooling are needed to see whether such a high accretionrate continues to be maintained. After the protostar growsto an SMS ( & M ⊙ ) via rapid accretion, it eventuallycollapses through the general relativistic instability to turninto a seed of high- z SMBHs (e.g., Mortlock et al. 2011).In this
Letter , we started the calculation from the ini-tial condition of a critical BE sphere with turbulence, andfound a single protostar formed without a major episode offragmentation. Since at . . ρ ∝ r − and v tan ≃ . v Kep , independent of the ini-tial conditions (Figs. 2a and 5), we expect that our conclu-sions depend only weakly on the initial setup. To confirmthis speculation, we need to investigate the dependence onthe initial conditions. In particular, strong turbulence couldprompt the efficient fragmentation, instead of forming a sin- gle SMS (e.g., Clark et al. 2011). Similarly, SMS formationfrom the proper cosmological initial condition remains to beexplored for future studies.In this simulation, we have neglected the effect of mag-netic fields. Previous studies suggest that magnetic fieldstrength could rival that of the turbulent energy (e.g.,Federrath et al. 2011; Turk et al. 2012; Latif et al. 2013a)with the effect of either increasing accretion efficiency (viamagnetic breaking producing a more spherical flow) or de-creasing it via protostellar jets (Machida et al. 2006). Thisexploration will be left for future investigations.
ACKNOWLEDGEMENTS
We thank the Enzo and yt support teams, especially BrianO’Shea and Matthew Turk for their useful advice. Wealso thank Takashi Nakamura, Takashi Hosokawa, NaokiYoshida, Shu-ichiro Inutsuka, Tsuyoshi Inoue, and KeiTanaka for their fruitful discussions. The results are ana-lyzed using the visualization toolkit for astrophysical dataYT (Turk et al. 2011). Numerical computations were carriedout on Cray XC30 at the Center for Computational Astro-physics of the National Astronomical Observatory of Japan.This work is supported in part by the grants-in-aid by theMinistry of Education, Culture, and Science of Japan (KI23 · REFERENCES
Abel, T., Bryan, G. L., & Norman, M. L. 2002, Science, 295, 93Begelman, M. C., Volonteri, M., & Rees, M. J. 2006, MNRAS,370, 289Bromm, V., & Loeb, A. 2003, ApJ, 596, 34Bryan, G. L., Norman, M. L., O’Shea, B. W., et al. 2014, ApJS,211, 19Chandrasekhar, S. 1964, ApJ, 140, 417Choi, J.-H., Shlosman, I., & Begelman, M. C. 2013, ApJ, 774,149Clark, P. C., Glover, S. C. O., Klessen, R. S., & Bromm, V.2011, ApJ, 727, 110Coppola, C. M., Longo, S., Capitelli, M., Palla, F., & Galli, D.2011, ApJS, 193, 7Di Matteo, T., Khandai, N., DeGraf, C., et al. 2012, ApJL, 745,L29Fan, X. 2006, New Astron Rev., 50, 665Federrath, C., Sur, S., Schleicher, D. R. G., Banerjee, R., &Klessen, R. S. 2011, ApJ, 731, 62Ferland, G. J., Peterson, B. M., Horne, K., Welsh, W. F., &Nahar, S. N. 1992, ApJ, 387, 95Foster, P. N., & Chevalier, R. A. 1993, ApJ, 416, 303Greif, T. H., Johnson, J. L., Klessen, R. S., & Bromm, V. 2008,MNRAS, 387, 1021Hosokawa, T., Omukai, K., & Yorke, H. W. 2012, ApJ, 756, 93Hosokawa, T., Yorke, H. W., Inayoshi, K., Omukai, K., &Yoshida, N. 2013, ApJ, 778, 178Inayoshi, K., & Omukai, K. 2011, MNRAS, 416, 2748Inayoshi, K., & Omukai, K. 2012, MNRAS, 422, 2539Inayoshi, K., Hosokawa, T., & Omukai, K. 2013, MNRAS, 431,3036Inoue, T., & Inutsuka, S.-i. 2008, ApJ, 687, 303Larson, R. B. 1969, MNRAS, 145, 271Larson, R. B. 1981, MNRAS, 194, 809Latif, M. A., Schleicher, D. R. G., Schmidt, W., & Niemeyer, J.2013a, MNRAS, 432, 668c (cid:13) , 000–000
K. Inayoshi, K. Omukai and E. Tasker
Latif, M. A., Schleicher, D. R. G., Schmidt, W., & Niemeyer, J.2013b, MNRAS, 433, 1607Machida, M. N., Omukai, K., Matsumoto, T., & Inutsuka, S.-i.2006, ApJL, 647, L1McKee, C. F., & Tan, J. C. 2008, ApJ, 681, 771Martin, P. G., Schwarz, D. H., & Mandy, M. E. 1996, ApJ, 461,265Mortlock, D. J., et al. 2011, Nature, 474, 616Narita, S., Hayashi, C., & Miyama, S. M. 1984, Progress of The-oretical Physics, 72, 1118Omukai, K. 2001, ApJ, 546, 635Regan, J. A., & Haehnelt, M. G. 2009, MNRAS, 396, 343Regan, J. A., Johansson, P. H., & Haehnelt, M. G. 2014, MN-RAS, 439, 1160Saigo, K., & Hanawa, T. 1998, ApJ, 493, 342Saigo, K., Matsumoto, T., & Hanawa, T. 2000, ApJ, 531, 971Shang, C., Bryan, G. L., & Haiman, Z. 2010, MNRAS, 402, 1249Silk, J. 1983, MNRAS, 205, 705Stibbe, D. T., & Tennyson, J. 1999, ApJL, 513, L147Turk, M. J., Smith, B. D., Oishi, J. S., et al. 2011, ApJS, 192, 9Turk, M. J., Oishi, J. S., Abel, T., & Bryan, G. L. 2012, ApJ,745, 154Wise, J. H., & Abel, T. 2007, ApJ, 665, 899Wise, J. H., Turk, M. J., & Abel, T. 2008, ApJ, 682, 745Yoshida, N., Omukai, K., & Hernquist, L. 2008, Science, 321,669Yoshii, Y., & Sabano, Y. 1979, PASJ, 31, 505 c (cid:13)000