Massive Star Formation via the Collapse of Subvirial and Virialized Turbulent Massive Cores
Anna L. Rosen Pak Shing Li Qizhou Zhang Blakesley Burkhart




Anna Rosen
anna.rosen@cfa.harvard.edu

0000-0003-4423-0660]Anna L. Rosen
NASA Einstein Fellow
ITC Fellow
Center for Astrophysics | Harvard & Smithsonian, 60 Garden St, Cambridge, MA 02138, USA


Astronomy Department, University of California, Berkeley, CA 94720, USA


Center for Astrophysics | Harvard & Smithsonian, 60 Garden St, Cambridge, MA 02138, USA


Center for Computational Astrophysics, Flatiron Institute,162 Fifth Avenue, New York, NY 10010, USA
Department of Physics and Astronomy, Rutgers, The State University of New Jersey,136 Frelinghuysen Rd, Piscataway, NJ 08854, USA




Similar to their low-mass counterparts, massive stars likely form via the collapse of pre-stellar molecular cores. Recent observations suggest that most massive cores are subvirial (i.e., not supported by turbulence) and therefore are likely unstable to gravitational collapse. Here we perform radiation hydrodynamic simulations to follow the collapse of turbulent massive pre-stellar cores with subvirial and virialized initial conditions to explore how their dynamic state affects the formation of massive stars and core fragmentation into companion stars. We find that subvirial cores undergo rapid monolithic collapse resulting in higher accretion rates at early times as compared to the collapse of virialized cores that have the same physical properties.

In contrast, we find that virialized cores undergo a slower, gradual collapse and significant turbulent fragmentation at early times resulting in numerous companion stars. In the absence of strong magnetic fields and protostellar outflows we find that the faster growth rate of massive stars that are born out of subvirial cores leads to an increase in the radiative heating of the core thereby further suppressing fragmentation at early times when turbulent fragmentation occurs for virialized cores. Regardless of initial condition, we find that the massive accretion disks that form around massive stars dominant the accretion flow onto the star at late times and eventually become gravitationally unstable and fragment to form companion stars at late times. 







§ INTRODUCTION
 

Massive stars play an essential role in the Universe. Their explosive deaths produce most of the heavy elements, enriching the interstellar medium (ISM) and future generations of stars. They are rare, yet the energy and momentum they inject into the ISM dwarfs the contribution of their more numerous low-mass cousins. This stellar feedback may set an upper limit on stellar masses, thereby affecting elemental abundances in the Universe. While the universal importance of massive stars is well understood the initial conditions, accretion history, and time span of their formation remains debated.

Massive stars form in dense (≳10^4-10^6 cm^-3), cold (∼10 K) turbulent gas within giant molecular clouds and giant massive filaments <cit.>. These condensations, commonly called clumps, have masses ranging from a ∼few tens to 10^5   M_⊙ and sizes of ≃ 0.5 - 2 pc in the Milky Way <cit.>. The most massive (i.e., those with M ≳ 10^3   M_⊙) and quiescent of these clumps have high surface densities with Σ≳ 0.05  g  cm^-2 and are likely the sites of massive star and protocluster formation <cit.>. 

These clumps host massive pre-stellar cores with typical sizes of 0.1 pc, which is a requirement for the Turbulent Core model for massive star formation <cit.>. In this model, massive stars form in a similar fashion to their low-mass counterparts via the monolithic collapse of massive pre-stellar cores that are supported by turbulence rather than thermal motions <cit.>. The stability of these cores can be described by their virial parameter given by α_ vir = 5 σ_ 1D^2 R_ c/G M_ c where σ_1D is the core's 1D velocity dispersion and M_ c and R_ c are the core mass and radius, respectively <cit.>. Neglecting external pressure and magnetic fields, cores with α_ vir≳ 1 are stable against gravitational collapse whereas those with α_ vir < 1 are subvirial and unstable to collapse. 

The Turbulent Core model requires that massive pre-stellar cores are in approximate virial equilibrium (i.e., α_ vir≃ 1). These cores then become marginally unstable to collapse to form a massive star or massive multiple system. The resulting formation time scale is several times the core free-fall time scale (t_ ff 10^5 yr) and the high degree of turbulence causes clumping, resulting in high accretion rates (Ṁ_ acc∼ 10^-4-10^-3 M_⊙  s^-1) that can overcome the radiation pressure associated with the star's large luminosity <cit.>. In agreement with this picture, observations show that massive cores live in highly pressurized environments and have non-thermal turbulent motions that dominate over thermal motions <cit.>. Additionally, due to the core's high angular momentum content, an optically thick accretion disk forms around the accreting massive star as the core collapses and delivers material at high rates via gravitational torques to the star <cit.>.

On smaller scales, massive clumps and cores fragment into pre-stellar cores via turbulent fragmentation that have masses larger then the mass and length scale dictated by thermal Jeans fragmentation <cit.>. Additionally, observations of massive protostellar cores show that they further undergo thermal Jeans fragmentation <cit.>. Such fragments may be the pre-cursers of low-mass pre-stellar cores that can accrete from the surrounding unbound gas to form massive stars as described by the Competitive Accretion model <cit.>. This model posits that a low-mass protostellar seed will accrete unbound gas within the clump as determined by its tidal limits and when it becomes massive enough it will then accrete at the Bondi-Hoyle accretion rate, Ṁ_⋆,   BH∝ v^-3 where v is the relative velocity of the gas. This model achieves high accretion rates onto the protostar under subvirial initial conditions, in contrast to the virialized conditions of the Turbulent Core model, since the gas velocity dispersion, and hence the level of turbulence, is low.

While both models can achieve the high accretion rates required for massive star formation, the specific initial conditions of the gas out of which they form will set the stage for the fragmentation properties of pre-stellar cores and the accretion history of massive stars. 

One way to determine the initial conditions of massive star formation is to study the demographics of massive starless clumps and cores. Most studies have found that massive cores are typically subvirial and should collapse within a gravitational free-fall time if the cores are not supported by magnetic fields <cit.>. In contrast, there have only been a few observational studies that have discovered virialized candidates <cit.>. If most massive cores are subvirial, then magnetic fields may provide additional support against collapse. Observations suggest that strong magnetic fields of the order of ∼1 mG are required for stabilizing massive pre-stellar cores. However, observations have demonstrated that magnetic pressure in dense, molecular gas is dynamically sub-dominant to gravity <cit.>. 






The dynamic state of massive cores will affect how mass is gathered to the star. Is most of the mass gathered via the direct global collapse of subvirial cores or is it slowly accumulated from the turbulent collapse of roughly virialized cores? The purpose of this work is to understand how the initial state of the pre-stellar core affects the growth of the resulting massive star and fragmentation of the core. However, massive stars are rare and form in highly embedded regions, therefore capturing the early moments of their formation and following how their mass is accreted onto the star is observationally challenging. They also have large luminosities throughout their formation since they contract to the main sequence while they are actively accreting and therefore radiation pressure is an important feedback mechanism during their formation <cit.>. Instead, we must turn to multi-dimensional radiation-hydrodynamic (RHD) numerical simulations to study the early formation of massive stars. 

Additionally, the accretion rate onto the massive star throughout the star formation process and the degree of core fragmentation into companions is likely dependent on the initial core properties such as the initial degree of turbulence and virial parameter.  We investigate this effect here by performing three-dimensional RHD numerical simulations of the collapse of subvirial and virialized turbulent massive pre-stellar cores. The simulations presented in this work are still highly idealized since we do not include magnetic fields or outflows, which will be addressed in future work. This paper is organized as follows: we describe our numerical methodology and simulation design in Section <ref>, we present and discuss our results in Sections <ref> and <ref>, respectively. Finally, we conclude in Section <ref>.









§ SIMULATION DETAILS


In this paper, we simulate the gravitational collapse of isolated turbulent massive pre-stellar cores using the  adaptive mesh refinement (AMR) code <cit.> to understand how their stability affects their fragmentation and the formation history of massive stars. We summarize the numerical methods, including the equations we solve in the  code, our refinement criteria, and how we treat the stellar radiation field in Section <ref>. We present the initial and boundary conditions for our simulations in Section <ref>. We refer the reader to  for more details on the numerics and overall algorithm (e.g., see Section 2.5 of ) since the physics included and the numerical methods used in this work are identical except where specified below. The initial conditions and numerical parameters are summarized in Table <ref>.



 §.§ Numerical Methods and Refinement Criteria


  uses the Chombo toolset to solve partial differential equations on block-structured AMR meshes <cit.> and solves the equations of gravito-radiation-hydrodynamics in the two-temperature, mixed-frame flux-limited diffusion (FLD) approximation on a Cartesian adaptive grid <cit.>. These equations are given by 

    ∂ρ/∂ t =  -∇·( ρ𝐯) - ∑_iṀ_i W(𝐱 - 𝐱_i)



    ∂(ρ𝐯)/∂ t   =    -∇·( ρ𝐯v) - ∇ P - ρ∇ϕ - λ∇ E_ R
          + ∑_i [𝐩̇_ rad,  i-𝐩̇_ iW(𝐱 - 𝐱_i)]



    ∂( ρ e)/∂ t   =    - ∇·[ (ρ e + P)v] - ρ𝐯·∇ϕ - κ_ 0Pρ(4π B - cE_ R ) 
          + λ( 2 κ_ 0P/κ_ 0R - 1 ) 𝐯·∇ E_ R - (ρ/m_ p)^2 Λ(T_ g) 
          + ∑_i [ε̇_ rad,  i-ε̇_i W(𝐱 - 𝐱_i) ]



    ∂ E_ R/∂ t   =   ∇·( c λ/κ_ 0Rρ∇ E_ R) + κ_ 0Pρ(4π B - c E_ R) 
          - λ( 2 κ_0  P/κ_0 R - 1) 𝐯·∇ E_ R 
    	- ∇·( 3 - R_2/2𝐯 E_ R) 
           +  (ρ/m_ p)^2 Λ(T_ g).



    dM_ i/dt = Ṁ_ i



    d x_i/dt = p_ i/M_ i



    dp_i/dt = -M_i ∇ϕ + ṗ_i



    ∇^2 ϕ = 4 π G [ ρ + ∑_i M_i δ(x - x_i)].

Here, ρ is the gas density,  ρv is the gas momentum, ρ e is the total internal plus kinetic gas energy density, and E_ R is the radiation energy density in the rest frame of the computational domain, κ_ 0P and κ_ 0R are the Planck and Rosseland mean opacities of the dust-plus-gas fluid, B is the blackbody function,  Λ is the rate of gas cooling that takes into account line and continuum processes for gas at temperatures above ≳ 10^3 K, λ is the flux-limiter, and R_2 is the Eddington factor. For more information on the flux-limiter, Eddington factor, hot gas cooling rate, and choice of dust opacities, we refer the reader to . 

In addition to the fluid,  contains Lagrangian sink particles used to represent radiating (proto)stars, indexed by subscript i, where each particle i is described with a mass m_ i, position x_i, and momenta p_i. Additionally, the particles accrete gas at a rate Ṁ_ i within four fine-level cells (i.e., the accretion around each star particle as weighted by the weighting kernel W(𝐱-𝐱_i) described in <cit.>).  We refer the reader to  for our description of our merging criterion for star particles when they pass within one accretion radius of each other. 

Equations <ref>-<ref> describe conservation of gas mass, gas momentum, gas total energy, and radiation total energy. Equations <ref>-<ref> describe how the star particles mass, velocity, and momenta are updated as they accrete mass and interact gravitationally with the surrounding fluid. Here ϕ is the gravitational potential that obeys the Poisson equation, given by Equation <ref>, that includes contributions from both the fluid and star particles. We also assume an ideal equation of state so that the gas pressure is 

    P=ρ k_ BT/μ m_ H = ( γ-1) ρ e_ T,

where T is the gas temperature, μ is the mean molecular weight,  γ is the ratio of specific heats, and e_ T is the thermal energy of the gas per unit mass. We take μ=2.33  and γ=5/3 that is appropriate for molecular gas of solar composition at temperatures too low to excite the rotational levels of H_2. The fluid is a mixture of gas and dust, and at the high densities that we are concerned with the dust is thermally coupled to the gas, allowing us to assume that the dust temperature is the same as the gas temperature.



Each star particle includes a model for protostellar evolution used to represent them as radiating protostars with frequency dependent stellar luminosities <cit.>. In order to properly model both the direct and indirect radiation pressure we use the  algorithm presented in <cit.>.  uses a frequency-dependent adaptive long-characteristics ray-tracing scheme that accurately solves the radiative transfer equation along rays to model the radial stellar radiation field from stars that is absorbed by dust and is coupled to a grey moment method, in our case FLD, which models the dust re-processed thermal radiation field intrinsic to the dusty fluid. The energy and momentum absorbed from the stellar radiation field, including contributions from the accretion luminosity,

    L_ acc = f_ rad G M_⋆Ṁ_⋆/R_⋆,

by the dust are added at a rate of ϵ̇_rad, i and ṗ_rad, i for particle i as shown in equations <ref>-<ref>. Here f_ rad is the fraction of the gravitational potential energy of the accretion flow that is converted to radiation and we take f_ rad=3/4 following <cit.> and M_⋆ and R_⋆ are the star's mass and radius, respectively. We use the frequency dependent dust opacities from <cit.> and divide the stellar spectrum and dust opacities into ten frequency bins (e.g., see Figure 1 in ). To reduce the computational cost, the rays adaptively split to conserve solid angle as they propagate radially away from their sources. We refer the reader to <cit.>, , and <cit.> for a complete description of our treatment of the direct and indirect radiation pressures modeled in this work.

For each simulation we begin with a base grid with volume (0.4 pc)^3 discretized by 128^3 cells and allow for five levels of refinement resulting in a maximum resolution of 20 AU.  As the simulation evolves, the AMR algorithm automatically adds and removes finer grids based on certain refinement criteria set by the user. Following , we refine cells if they meet at least one of the following criteria: (1) any cell on the base level (i.e., level 0) that has a density equal to or greater than the core's edge density so that the entire prestellar core is refined to level 1; (2) any cell where the density in the cell exceeds the Jeans density given by 

    ρ_ max,J  = π J^2_ max c_ s^2/G Δ x_l^2

where c_ s= √(k_B T/(μ m_ H)) is the isothermal sound speed, Δ x_l is the cell size on level l, and J_ max is the maximum allowed number of Jeans lengths per cell, which we set to 1/8  <cit.>; (3) any cell that is located within at least 8 cells of a sink particle; and (4) any cell within which the radiation energy density gradient exceeds |∇ E_ R| > 0.15 E_ R/Δ x_ l (i.e., where the radiation energy density changes by more than 15% over the size of a single cell). This procedure is repeated recursively on all levels after every two level updates.  Finally, a sink particle can only be created when the Jeans density is violated on the finest level in which we take J_ max=1/4 in Equation (<ref>) following the artificial fragmentation tests of <cit.> that found that values of J≲ 1/4 causes artificial fragmentation to be avoided in gravitational collapse simulations.

When a cell on the finest level exceeds the Jeans density we place a sink particle in that cell whose mass is the excess matter in that cell. The new sink particle will evolve according to equations <ref>-<ref>.



 §.§ Initial and Boundary Conditions


In this work, we perform two simulations that are identical in every way except for their initial one-dimensional velocity dispersion (σ_1D) in order to determine how the dynamic state of the star-forming core (i.e., virialized versus subvirial) affects their fragmentation and the growth of massive stars. Following , we begin with an isolated pre-stellar core of molecular gas and dust,  where we assume a dust-to-gas ratio of 0.01, with mass M_ c=150   M_⊙, radius R_ c = 0.1 pc, and initial gas temperature of 20 K corresponding to a surface density of Σ = M_ c/π R^2_ c = 1  g   cm^-2 consistent with extreme massive star-forming environments such as W49, Sgr B2, W51, and Cygnus X1 and massive core densities <cit.>. The corresponding mean density of the core is ρ̅ = 2.4 × 10^-18 g   cm^-3 (1.2 × 10^6  H   nuclei   cm^-3) and its characteristic free-fall collapse time scale is t_ ff≈ 42.6  kyr.  The core follows a density profile ρ(r) ∝ r^- 3/2 in agreement with observations of massive cores at the ∼0.1 pc scale and clumps at the ∼1 pc scale that find values of κ_ρ = 1.5-2 <cit.>. Each core is placed in the center of a 0.4 pc box that is filled with hot, diffuse gas with density ρ_ amb = 0.01 ρ_ edge where ρ_ edge is the density at the core boundary and temperature T_ amb = 2000 K so that the core is isolated and is in thermal pressure balance with the ambient medium. 
We set the opacity of the ambient medium to zero (i.e., the ambient medium is dust-free). 


We include turbulence by seeding the initial gas velocities (v_x, v_y, and v_z) with a velocity power spectrum, P(k) ∝ k^-2, as is expected for supersonic turbulence <cit.>. We include modes between k_ min = 1 to k_ max = 256 and take the turbulence mixture of gas to be 1/3 compressive and 2/3 solenoidal, which is consistent for the natural mixture of a 3D fluid <cit.>. The onset of turbulence modifies the density distribution and we allow the turbulence to decay freely. For both simulations, we use the same velocity perturbation power spectrum at initialization but set the initial σ_ 1D amplitude to two different values: our virial run, , has σ_ 1D= 1.2 km/s corresponding to α_ vir = 1.1 and ℳ = ⟨σ_ 1D⟩_m/c_s = 4.4 and our subvirial run, , is initialized with σ_ 1D= 0.42 km/s yielding α_ vir = 0.14 and ℳ =1.6.[We note that run  is the same run as  from .] We note that allowing the turbulence to decay is somewhat unrealistic. However, this simplification should have little effect on our results since the decay timescale, ≈ D/σ_ 1D where D is the core diameter <cit.>, is ∼0.16 Myr and ∼0.46 Myr for runs  and , respectively; which are much longer than the runtime for both simulations presented in this work. Therefore, we expect that allowing the turbulence to decay should have a negligible effect on our results. 

Our boundary conditions for the radiation, gravity, and hydrodynamic solvers are as follows. For each radiation update, we impose Marshak boundary conditions that bathe the simulation volume with a blackbody radiation field equal to E_0 = 1.21 × 10^-9  erg   cm^-3 corresponding to a 20 K blackbody but we allow radiation generated within the simulation volume to escape freely <cit.>. We set the gravitational potential, ϕ, to zero at all boundaries when solving equation <ref> <cit.>. Since the core boundaries are far removed from the domain boundaries we do not expect this choice of boundary conditions to lead to any significant square artifacts near the domain boundaries. Finally, we impose outflow boundary conditions for the hydrodynamic update by setting the gradients of the hydrodynamic quantities (ρ,  ρv, ρ e ) to be zero at the domain when advancing equations <ref>-<ref> <cit.>.





§ RESULTS


In this section, we describe and compare our results for runs  and  described in Section <ref> and summarized in Table <ref>. Run  was run on the Hyades supercomputer located at UC Santa Cruz and run  was run on the NASA supercomputer Pleiades located at NASA Ames. We run each simulation to the point where the simulation takes too long to evolve because the majority of the bubble shells are refined to the finest level, severely increasing the computational cost of the simulation. We use the yt package <cit.> to produce all the figures and quantitative analysis shown below.



 §.§ Core Collapse Properties and  Influence of Radiative Feedback



We show a series of the density slices and the corresponding Eddington ratio slices, f_ Edd= F_ rad/F_ grav where F_ rad is the total radiation force from both the stellar and dust-reprocessed radiation fields and F_ grav is the gravitational force, for runs  and  in Figure <ref>. The two left hand columns of Figure <ref> show the density structure and collapse evolution along the yz-plane with velocity vectors over-plotted for runs  (left-most column) and  (center left column), respectively. These panels show that the pre-stellar core in run  undergoes a rapid, monolithic collapse whereas the core in run  undergoes a gradual collapse leading to a slower growth rate for the massive star. At the end of run , t=0.87   t_ ff, the primary star has a final mass of 61.66 M_⊙ whereas for run  the primary star is significantly lower in mass with a mass of 46.78 M_⊙ at this time. This result suggests that the collapse of subvirial cores, as compared to virialized cores, may lead to massive stars that have larger masses at birth.

The slower growth of the primary star in run  occurs because the core envelope has a larger velocity dispersion and angular momentum content resulting in greater support against direct gravitational collapse. This is seen in the velocity vectors over-plotted on the density slice plots that show a randomized orientation in the core envelope whereas in run  the velocity vectors point radially towards the massive star demonstrating that the core collapses monolithically. Instead, for run , the velocity vectors in the core envelope's low density regions are oriented towards the high density flattened filaments that form out of the core's turbulence at early times. This effect causes the envelope material to be accreted onto filaments rather than fall directly towards the primary star.



At high stellar masses the radiation pressure associated with the massive star's high luminosity launches radiation pressure dominated bubbles that expand away from the star, reducing the infall of material to the accreting primary star <cit.>. When this occurs an optically thick accretion disk has formed around the massive star due to conservation of angular momentum as the core collapses and the bubbles are launched above and below the accretion disk as material is funnelled to the massive star. These low-density expanding bubbles are seen in the density slice plots for both runs in Figure <ref> at late times. We find that they are not sustained until the primary stellar mass reaches ≳ 30 M_⊙ for both runs because the ram pressure from the collapsing core material quenches the bubbles at early times. However, as the stellar mass of the primary star, and therefore its luminosity, increases these bubbles expand and grow larger with time. These bubbles occur in regions where f_ Edd≫ 1  as shown in the two right-hand columns of Figure <ref>. These panels show that at early times (t=0.4   t_ ff), when the star is 19.59 M_⊙ (13.82 M_⊙) for run  () the inner regions of the core nearest the primary star are super-Eddington (i.e., f_ Edd>1) but that the radiation pressure dominated bubbles have not been launched yet for either run at this point. This is because the acceleration of the material is initially controlled by the gravitational collapse of the core material and is falling towards the star. Thus, it takes time for the radiation pressure to overturn the accretion flow onto the star even though the gas near the star is super-Eddington. 

As the primary star continues to grow in mass the region of the core that becomes super-Eddington increases in size. This is more apparent for run , which leads to the majority of the core to become super-Eddington at late times because the star is very massive and therefore extremely luminous with a zero-age main sequence luminosity of L_⋆ = 5.38 × 10^5  L_⊙ <cit.>. The f_ Edd panels shown in Figure <ref> also show that the dense filaments that are inherent to the turbulent structure of the core for both runs achieve f_ Edd≲ 1. Thus, these filaments are able to self-shield from the strong radiation pressure associated with the primary star and can be incorporated into the star-disk system as the core collapses. 

To further illustrate the gravitational collapse of the core and the launching of the low-density bubbles driven by radiative feedback from the primary star we show a series of thin projection plots of the radial momentum, ρ v_ r, of gas with respect to the primary star for runs  (top row) and  (bottom row) in the yz-plane in Figure <ref>. Each projection is taken over a depth of 1000 AU with the primary star at the center of each panel. Values of positive ρ v_ r denote gas that is moving away from the primary star whereas negative values denote gas that is falling towards it. These snapshots show the growth and expansion of the radiation pressure dominated bubbles is perpendicular to the accretion disk. The accretion disk is larger in size in run  since the pre-stellar core has a larger angular momentum content due to its larger initial velocity dispersion (e.g., see Table <ref>). This `flashlight' effect allows material to be funnelled to the star by the accretion disk, as demonstrated by the negative values of ρ v_ r along the equatorial plane of the star, while the radiative flux escapes along directions above and below the disk <cit.>. We describe the accretion disk structure and evolution in more detail in Section <ref>.

For run , Figure <ref> shows that the core undergoes a more significant collapse than run  in agreement with the density slice plots in Figure <ref> and that the full core, at early times, is collapsing overall.  This leads to an overall higher density for the core in run  than the core in run  as the simulations progress. In contrast, the radial momentum plots for run  show that while the majority of the core, especially the inner parts, are collapsing towards the star a significant portion of the outer envelope is expanding. This is due to the higher degree of turbulence providing support against collapse.  As the envelope expands, its velocity and density decreases with time. Eventually the majority of the envelope begins to collapse towards the primary star at late times. This Figure also shows that the turbulence produces dense, optically thick filaments as the core collapses. Since the core in run  has a low velocity dispersion and collapses immediately the filaments are typically wider than those that are produced in run . 

The evolution of the core expansion is shown in Figure <ref>, which shows the spherically averaged density profiles of the core as a function of radius  for run  (top left panel) and run  (top right panel). The gray dashed line in each panel denotes the initial density distribution of the core before the turbulent velocity structure is added to the gas velocities (i.e., at t=0), described in Section <ref>. These panels present the density in log scale and the radius in linear scale to show the overall density evolution of the cores and how turbulence can modify the density distribution. As these panels demonstrate the onset of turbulence causes the cores to expand and the core gas at the core boundary mixes with the ambient medium. While this effect is a limitation due to treating our core as an isolated object in pressure balance with a low-density ambient medium it should have little effect on our results since the cores undergo inside-out collapse <cit.>. For run  the expansion due to the onset of turbulence is small, only increasing by ∼3 kAU (15% increase) in size. However, the onset of turbulence for run  causes the overall core to expand to ∼30 kAU (50% increase) at early times and continues to expand as the simulation progresses. This is because the initial global velocity dispersion of the core in run  is a factor of ∼3 larger than run . As these panels demonstrate, the core in run  undergoes a less rapid inside-out collapse than run  because of the greater support from turbulence. This causes the inner core collapse in run  to be more gradual leading to a slower growth rate for the primary star. 

The bottom panels in Figure <ref> show a zoom-in of the inner density distribution of the cores in run  (bottom left panel) and run  (bottom right panel) but now both panels are in log-log scale to highlight the density distribution near the primary star.  The inner regions of the core collapse for both runs and this gas is incorporated into the massive star or is blown away due to the strong radiation pressure from the massive star. We find that the inner density profile of the core for run  is significantly larger in density at t=0.60   t_ ff as compared to run  for distances ≲ 300 AU from the massive star. This is because the radiation pressure dominated bubbles have not been launched yet at this time for run  but they have been launched in run  due to the larger stellar mass attained in a shorter amount of time (e.g., M_⋆=36.6   M_⊙ for run  as compared to M_⋆=26.3   M_⊙ for run ). 

At late times (e.g., t=0.80   t_ ff and t=0.87   t_ ff) the average density of the radiation pressure dominated bubbles is significantly larger in run   than the bubble density of run  in agreement with the density slice plots in Figure <ref>. This results in a larger overall radial momentum magnitude in the bubbles at late times for run  as compared to run  as shown in Figure <ref>. To further show this we plot the same radial momentum thin projection plots (top row) and the mass-weighted density thin projection plots (bottom row) along the yz-plane for both runs but at the same stellar mass at the end of run , which ran to t=0.95   t_ ff at which point the primary star has a mass of ∼52   M_⊙ in Figure <ref>. Comparing these snapshots for both runs shows that the radiation pressure dominated bubbles that are above the star are approximately the same size and that the bubbles below the primary star differ in size, with the bottom bubble in run  being smaller than that in run . This difference in bubble sizes is due to the turbulent structure and larger ram pressure of the infalling material that surrounds the bubble in run . Furthermore, the radial momentum in the radiation pressure bubbles in run  are larger in magnitude than the bubbles in run  because the bubbles in run  have a larger overall density as shown in the bottom row of Figure <ref>. Since the overall core density in run  is smaller in magnitude than run  as the cores evolve, radiation pressure is more efficient at launching the gas leading to lower densities in the bubbles.












We show thin projections of the mass-weighted density and temperature for runs  (density: left hand column, temperature: middle right column) and  (density: center left column, temperature: right hand column) in Figure <ref> that show the growth of the radiation pressure dominated bubbles and shell structure as each simulation progresses.  Each panel is (20,000 AU)^2 in size and the depth of the projection is 1000 AU with the massive star at the center of each panel. The low density radiation pressure dominated bubbles eventually become transparent to the stellar radiation field since the bubbles exceed the dust sublimation temperature of 1500 K due to radiative heating from the massive star <cit.>. When this occurs the majority of the stellar radiation is absorbed in the dense bubble shells and the momentum and energythereby accelerating the shells' expansion.

The second (third) panel for run  () show that the shells of the low-density radiation pressure dominated bubbles develop small scale turbulence as they expand. This turbulence may be due to the growth of radiative Rayleigh Taylor (RT) instabilities that develop in the dense shells. 
These instabilities occur at the interface between two fluids of different densities in which the fluid at the bottom of the interface is a lower-density radiation pressure dominated fluid (e.g., the rarefied radiation pressure dominated bubbles) that pushes on a higher density, less radiatively dominated medium <cit.>.

In classical RT instability theory <cit.>, the amplitude η of linear perturbations that develop due to corrugations or asymmetries that form at the interface will grow with time as η(t) ∝expω t where

    ω = √(g k ρ_2 - ρ_1/ρ_2+ρ_1)


and ρ_1 and ρ_2 are the densities at the interface (i.e., the shell and radiation dominated bubble, respectively), g is the magnitude of the gravitational acceleration, k=2π/λ is the wave number, and λ is the wavelength of the perturbation <cit.>. In the limit where ρ_2 ≫ρ_1, such as at the interfaces between the dense shells and low-density radiatively driven bubbles presented in this work, Equation <ref> reduces to ω≈√(-g k). Hence, the interface becomes unstable when ρ_2 < ρ_1.

As discussed in , the initial turbulence and anisotropic absorption of the stellar radiation field seeds the development of these perturbations that can grow exponentially with time as the rarefied bubbles expand. These instabilities may grow into RT fingers, as shown in the right top panel Figure <ref> that shows the radial momentum with respect to the star and shows two fingers penetrating the top bubble. Comparison with the last bottom panel in Figure <ref> show that these fingers likely originate from the two corrugations present in the top bubble at this time. Eventually, these fingers may be able to fall to the star-disk system and deliver material to the star potentially aiding accretion unless radiation pressure from the star overturns their infall. However, given the computational expense we are unable to follow the further evolution of the growth of these RT fingers to determine if they will be incorporated into the star-disk system.





 






























 §.§ Accretion History and Accretion Disk Evolution


Figure <ref> shows the stellar mass (top panel) and accretion rate (bottom panel) as a function of time for runs  and . The accretion rate for the massive star in run  is larger by a factor of ∼2 as compared to the accretion rate onto the massive star in run  up to t ≈ 0.55  t_ ff, resulting in a faster stellar growth. The higher accretion rate achieved in run  is a result of the core's global collapse as compared to run , which undergoes a slower collapse due to the additional turbulent support. At t ≈ 0.55   t_ ff the accretion rate onto the massive star in run  increases and becomes comparable to the accretion rate attained in run  and they remain comparable until t ≈ 0.75   t_ ff.

Collapsing optically thick filaments inherent to the core's turbulent structure can lead to variable accretion onto the star and also overcome the radiation pressure barrier, delivering material to the star that would otherwise be blown away. In this scenario, the radiative flux escapes through low density channels while the denser filaments can accrete onto the star <cit.>. As discussed in the previous Section, the filaments in run  are wider than those that form in run  and begin to collapse instantly. This results in an accretion rate onto the primary star that grows in time and is relatively smooth up to t ≈ 0.55  t_ ff. For run , the turbulence has more time to develop thin, dense filaments before they collapse onto the star as compared to run  and the accretion of these filaments causes the variability of the accretion rate onto the star at early times (i.e., between t=0.3 - 0.55   t_ ff). This variability can increase the accretion rate by up to a factor ∼2 for t ≲ 0.55   t_ ff. 

At later times the accretion rate becomes much more highly variable and can quickly increase by up to an order of magnitude before dropping to lower values for both runs. This increase in accretion variability is due to accretion of infalling filaments and/or the inflow of mass from gravitational torques induced in the accretion disk. In order to determine which process dominates the accretion variability we examine the evolution and structure of the accretion disk for runs  and . Figure <ref> shows a series of density slices of the accretion disk that forms around the massive star for run  (top row) and run  (bottom row) with the velocity streamlines over-plotted showing that this over-dense structure eventually undergoes circular motion around the central star. The massive star is at the center of each panel, denoted by the yellow circle and companion stars that are color-coded by mass are over-plotted. These stars are formed via turbulent fragmentation or via disk fragmentation, which we will describe in more detail in Section <ref>. As these panels show a noticeable high density accretion disk (i.e., an accretion disk with a radius larger than the 80 AU accretion zone radius of the sink particle that has a circular structure and undergoes circular motion around the massive star) begins to form around the massive star in run   at t ≈ 0.4   t_ ff whereas in run  a noticeable accretion disk does not begin to form until t ≈ 0.65   t_ ff when the star has a mass of 41.1   M_⊙. The accretion disk grows in size as both simulations progress due to conservation of angular momentum: as the core collapses, material that is farther away from the massive star has a larger net angular momentum and therefore will be circularized at a distance farther from the star. The disk can be disrupted due to the passage of companion stars as seen in the third snapshot of the accretion disk projection for run  in Figure <ref>. However, after this occurs and the core continues to collapse, infalling material is circularized rebuilding the disk.

As the accretion disk grows in size and mass, it becomes gravitationally unstable and develops spiral structure. The gravitational torques induced delivers material to the star <cit.>. This, in turn, leads to a highly variable accretion rate as shown in the bottom panel of Figure <ref>. Indeed, we see that the accretion rate, while relatively smooth throughout most of the simulation becomes highly variable for run  for t ≳ 0.7 t_ ff and the magnitude of the variability increases with time as the both the accretion disk and massive star increase in mass. This variability in run  occurs because the star and disk are orbiting around the star-disk system's center of mass and this further drives gravitational torques in the disk that eventually lead to disk fragmentation into companion stars. Likewise, we see the same behavior for run  but this behavior begins at earlier times, t ∼ 0.5 t_ ff, since the accretion disk forms earlier than the disk in run . Hence, we find that at late times for both runs most of the accretion variability is due to how mass is driven from the disk to the star via gravitational torques. 











While most of the accretion variability is due to gravitational torques in the accretion disk driving material to the central massive star at late times, accretion variability is also a result of infalling filaments and RT instabilities that accrete onto the star-disk system. To determine this and differentiate it from accretion variability attributed to the accretion disk we plot the ram pressure of the inflowing core material in Figure <ref> for material that is 1000 AU away from the star as a function of primary stellar mass since the accretion disk radius is < 1000 AU throughout both simulations. To compute the inflow ram pressure we take a sphere 1000 AU in radius centered on the most massive star for each run. We compute the area-weighted mean ram pressure (solid lines) and the mass-flux-weighted mean ram pressure (dashed lines) for inflowing material for runs  (teal lines) and  (pink lines). These quantities are defined by 

    ⟨ P_ ram⟩_w = ∫ρ v_ r^2 wdA/∫ w dA

where v_ r is the radial velocity, and the weighting function w is unity for the area-weighted average(⟨ P_ ram⟩), and we take w=ρ v_ r for the mass-flux-weighted average (⟨ P_ ram⟩_Ṁ) where we only include contributions from gas that is inflowing to the star (i.e., v_ r>0). We include the mass-flux-weighted mean ram pressure because it is a better representation of the ram pressure of the material that can be accreted onto the star. 

Throughout the simulation, the ram pressure of the inflowing material is larger for run  except at points when the accretion rate onto the star jumps significantly for run . These jumps in ram pressure are much more pronounced in the mass-flux-weighted mean ram pressure and can cause it to jump by up to ∼ 3 orders of magnitude. By comparing these ram pressure spikes in run  with the accretion rate shown in the bottom panel of Figure <ref>, we find that these spikes correspond to large jumps in the accretion rate. Hence, we find that the turbulent structure of the core can lead to highly variable accretion onto the massive star for run  but is negligible for run . This is likely because the core in run  develops a more filamentary structure than the core in run  because the core does not undergo instantaneous collapse and therefore has more time for supersonic shocks to collide to form more thin, dense filaments.



 §.§ Core and Disk Fragmentation into Companion Stars

 
Companion stars can form via turbulent fragmentation or disk fragmentation. Turbulent fragmentation occurs in molecular clouds when supersonic shocks collide and form high density thin sheets and elongated structures surrounded by low-density voids (i.e., filaments). Dense parts of these structures fragment into dense clumps and cores when they become Jeans' unstable <cit.>. The hierarchical fragmentation continues until stars form inside the dense cores. Since the filaments are due to the collision of supersonic shocks in the ISM, the degree of fragmentation depends on the velocity dispersion and therefore the virial state of the core. We find that the core in run  undergoes significant turbulent fragmentation but the core in run  does not undergo any turbulent fragmentation as shown in Figure <ref> that shows the total mass and number of companion stars as a function of time for runs  (teal lines) and  (pink lines). The pre-stellar core in run  has a larger initial velocity dispersion than the core in run  and therefore provides support against collapse allowing time for the filaments to form and then undergo longitudinal collapse (i.e., collapse modes that are parallel to the symmetry axis of the filament) due to their self-gravity <cit.>.  In contrast, the core in run  has less supersonic shocks due to its weaker turbulence and lower ℳ.  Because of these effects and aided by the fact that the core in run  instantly collapses,  high-density filaments that form do not have enough time to develop thin, elongated structures that can then collapse into stars since their motions are dominated by infall towards the primary star. Therefore, they do not become Jeans unstable and fragment into low-mass companion stars. 



As Figure <ref> shows, the core in run  starts to fragment into companion stars via turbulent fragmentation at t ≈ 0.35   t_ ff, forming eight companion stars by t ≈ 0.60   t_ ff. After that, turbulent fragmentation into companion stars for run  is suppressed. This is likely due to the enhanced radiative heating of the core material by the massive star when the star reaches ∼ 25 M_⊙, which we describe in more detail in the next subsection. At t ≈ 0.87   t_ ff the formation rate of companion stars is greatly enhanced due to disk fragmentation, nearly doubling the number of companion stars by a factor of two by t ≈ 0.9   t_ ff. This occurs because the disk becomes gravitationally unstable because the massive star and disk are orbiting their common center of mass. The gravitational torques induced in the disk drive high density spiral waves that can become Jeans' unstable and collapse to form disk borne stars <cit.>. Referring back to Figure <ref> and Figure <ref>, we see that the accretion disk in run  begins to fragment at the end of the simulation, forming two companion stars with M_⋆ > 0.04   M_⊙ (i.e., the mass that we initialize our protostellar model). 

To demonstrate how the accretion disk becomes unstable and fragments to form close-in companion stars we show the disk surface density (left panels) and disk Toomre Q parameter (right panels) for the final snapshots of runs  (top row) and  (bottom row) in Figure <ref>. The Toomre Q parameter is given by

    Q = c_s κ/π G Σ

where c_ s is the sound speed, κ is the epicyclic frequency which is equal to the angular velocity Ω for a Keplerian disk, and Σ is the disk surface density <cit.>. Values of Q<1 denote locations where the disk is unstable to gravitational collapse whereas regions with Q>1 are stable and not prone to collapse. For run , the majority of the disk has Q<1 and the disk has fragmented into several low-mass companion stars by this time. We also note that the final disks in both simulations are relatively massive with masses of ∼ 10.5   M_⊙ for run   and ∼9.4   M_⊙ for run  and they are still being fed by the collapsing core at this time.

Most of the disk-borne stars in run  are very low-mass, as shown in the histogram in Figure <ref> that shows the distribution of stellar masses for the companion stars before disk fragmentation (top panel) compared to the final distribution of companion star masses at the end of the simulation (bottom panel). The bottom panel shows that the newly formed disk borne stars have masses ≲ 0.2   M_⊙ whereas several of the stars that formed via turbulent fragmentation have attained masses ≳ 1   M_⊙ by the end of the simulation. At the end of both simulations we find that they consist of a hierarchical system that contains a massive primary and a series of low-mass companions instead of a massive binary system. This may also be attributed to our treatment of merging sink particles because we do not merge sink particles once their mass exceeds 0.04   M_⊙ and therefore our strict merging criteria may lead to an excess of low-mass companions, which we discuss further in Section <ref>.







 §.§ Evolution of the Core Properties


We show the average mass-weighted core velocity dispersion (top row) and temperature (bottom row) as a function of primary stellar mass (left column) and time (right column) for runs  (teal solid lines) and  (pink dashed lines) in Figure <ref>, respectively. We plot these values as a function of both simulation time and primary stellar mass because time is a good proxy to follow these properties as the core collapses whereas primary stellar mass traces the impact of stellar feedback. As <cit.> noted infall motions leads to an increase in the velocity dispersion in collapsing clouds. In agreement, we see that the rapid collapse of the core in run  causes the core's mass-weighted velocity dispersion to become larger than the mass-weighted velocity dispersion in run  beginning at t ≈ 0.2   t_ ff when the primary star is ∼few M_⊙ such that radiative feedback is unimportant at this time. The velocity dispersion for run  continues to increase and remains larger than the velocity dispersion in run  until the star has a reached a mass of ∼25   M_⊙. At this stellar mass for both runs radiative feedback from the massive star becomes super-Eddington and launches radiation pressure dominated bubbles along the polar directions of the star. As these bubbles expand they become RT unstable and drive turbulence in the shells as shown in the thin mass-weighted density projections for runs  (left most column) and  (center left column) in Figure <ref> <cit.>, thereby increasing the overall core's velocity dispersion in both runs. Hence, we find that both gravitational infall and radiative feedback can lead to an increase in the overall velocity dispersion of the cores.




Similarly, gravitational collapse and radiative feedback also increases the global core temperature as shown in the bottom panels of Figure <ref>. As the cores in runs  and  collapse and the primary stars grow in mass their overall mass-weighted temperature increases with time, however we see that this effect is much more pronounced for run  up to t ≈ 0.6   t_ ff (corresponding to M_⋆∼ 25   M_⊙ for the primary star in run ) in agreement with the behavior seen in the core's velocity dispersion. As described in Section <ref>, turbulent fragmentation for run  occurs during t ≈ 0.35-0.6   t_ ff but no turbulent fragmentation occurs for run  during this time. Radiative heating reduces fragmentation and therefore we find that both radiative heating from the more massive star in run  and enhanced heating due to collapse also suppresses fragmentation of the core in run  as compared to run  <cit.>. As the primary star becomes more massive, making radiative heating more important, the core temperature continues to increase for both runs. However, radiative feedback is more efficient at heating the core for run  when the star has a mass of ∼25 to ∼37 M_⊙. At higher stellar masses the temperature is comparable for both simulations. Hence, we find that at early times gravitational collapse can lead to higher core temperatures whereas at late times and high stellar masses the overall increase in the core temperature is due to stellar feedback.







§ DISCUSSION


The purpose of this work is to understand how the initial dynamic state of the core affects the accretion history of massive stars, core fragmentation, and formation of companion stars. 

Most notably we find that highly subvirial pre-stellar cores undergo rapid, monolithic collapse and no turbulent fragmentation whereas virialized cores undergo a slower, more gradual collapse and significant fragmentation thus forming a small cluster of stars. In what follows, we discuss how the pre-stellar core's virial state affects the growth rate of massive stars and core fragmentation in Section <ref> and how disk fragmentation can produce close-in companion stars in Section <ref>. The simulations presented here do not include magnetic fields which are ubiquitous  in massive star forming regions and are therefore an important ingredient in the star formation process. In light of this, we discuss our motivation for neglecting them, what effect we expect magnetic fields to play in the collapse of massive cores into massive stellar systems, and our plans to include them and other forms of stellar feedback such as collimated outflows in future work in Section <ref>.



 §.§ Fragmentation and Multiplicity



  §.§.§ Formation of Isolated Massive Stars versus Protoclusters via Turbulent Fragmentation


In this work we showed that the dynamic state of massive pre-stellar cores, as described by their virial parameter that compares the core's kinetic energy to its gravitational potential energy, affects the accretion history of massive stars and core fragmentation into companion stars. Specifically, we find that virialized cores for the initial conditions chosen in this work undergo hierarchical turbulent fragmentation up until t∼ 0.6   t_ ff, when radiative heating from the massive star becomes significant. This results in a small cluster that contains one massive star with several low-mass companions. In contrast, we find that highly subvirial cores undergo rapid, monolithic collapse that yields higher accretion rates onto massive stars and no turbulent core fragmentation into companion stars. When the star is sufficiently massive the enhanced radiative heating, in comparison to run , heats up the core material and therefore suppresses fragmentation. Therefore, we conclude that subvirial cores are more likely to form massive stars with higher accretion rates at the onset of core collapse than virialized cores that have similar physical properties (e.g., surface density, mass, and density profile) and that the lower degree of turbulence in subvirial cores likely result in less turbulent fragmentation for subvirial cores. Our results suggest that star clusters that host at least one massive star may have originated from cores that have a greater degree of turbulence at the onset of collapse and that wide companion stars form either coevally or at a later stage after the birth of massive protostars in agreement with observations <cit.>.

Similarly, <cit.> studied the early collapse of massive supercritical magnetized cores with different initial virial parameters and magnetic field strengths by performing a suite of radiation-magnetohydrodynamic simulations and found that fragmentation is inhibited for clumps with low turbulence (ℳ≲ 3) regardless of the magnetic field strength in agreement with run  presented in this work (ℳ = 1.7). They conclude that magnetized cores with a high initial turbulence (ℳ∼ 6) undergo fragmentation but that the number of fragments that form is larger for weakly magnetized cores. They conclude that the turbulent state of the core is more important than magnetic fields in determining whether the core will fragment into a small cluster or not but that the degree of fragmentation and distribution of fragments does depend on the magnetic field strength. Similarly, in the absence of magnetic fields we find that significant fragmentation also occurs in massive cores with weaker initial turbulence (ℳ∼ 4.5, run  presented in this work) leading to the formation of a small cluster in which ∼ 9 low-mass companion stars are formed via turbulent fragmentation. 

The simulations presented here and in <cit.> only explore a small parameter space of massive pre-stellar cores, however observed massive dense cores span a wide range of physical properties. Based on our results that the initial turbulence of these cores affect fragmentation we expect that weakly subvirial cores with α_ > 0.14, the initial value used for run , should undergo fragmentation but to a lesser degree than virialized cores. In agreement, observations of massive star-forming cores show a diverse population of fragments. <cit.> studied 20 evolved high-mass star-forming regions that host at least one young massive stellar object (MYSO; i.e., M_⋆≳ 8   M_⊙) and found that these clumps have diverse fragmentation morphologies ranging from regions that are dominated by single evolved high-mass cores to those that fragment up to 20 cores. They conclude that the diversity in fragmentation of evolved massive-star forming clumps depends on the clump's dynamic state, which depends on the interplay of self-gravity, turbulence, magnetic fields, and heating from stellar feedback that can suppress fragmentation. Hence, one way to possibly identify massive stars that form from highly subvirial cores would be to find young, embedded massive stars that have no long distance companions but close-in companions due to disk fragmentation (<cit.>, see Section <ref>). 
 


  §.§.§ Formation of Close-In Companion Stars via Disk Fragmentation



Turbulent fragmentation is likely responsible for wide binaries but ∼80 % of O-stars (M_⋆≳ 16   M_⊙) are found in close-in multiple binaries with separations ≲ 700 AU  and a flat mass-ratio, q=M_ 2/M_ 1, distribution where M_ 1 (M_ 2) is the primary (secondary) stellar mass <cit.>. This large binary fraction may originate from the formation process via disk fragmentation rather than by direct capture <cit.>. In agreement with this scenario, we find close-in companion stars form via disk fragmentation at late times regardless of the initial virial state of the core.  We find that the companion stars formed via disk fragmentation at the end of the simulations presented in this work have M_ 2≪ M_ 1 resulting in extreme mass ratios. 

In agreement with our results, recent observational studies have now found asymmetric, fragmented Keplerian accretion disks around massive proto-O stars. <cit.> found a binary system of compact dusty objects that are separated by ≈300 AU in projection at the center of an asymmetric accretion disk of size ≈ 1200 AU. They conclude that the primary massive proto-O star has a stellar mass of ∼ 20   M_⊙ and its low-mass companion likely has a mass within 0.2-2   M_⊙ yielding a mass ratio of q ≈ 0.01-0.1. Similarly, <cit.> observed the G11.92-0.61 system with ALMA  and found a fragmented disk around a ∼34 ± 5   M_⊙ proto-O star with a low-mass companion star with mass ≲ 0.6   M_⊙ yielding a mass ratio of q ∼ 0.015. These studies are the first to observe the formation of a binary star via disk fragmentation around a young proto O-star. 

The extreme mass ratio observed in these objects suggest that disk fragmentation into low-mass companions can occur when the primary star is substantially massive in agreement with the simulations presented in this work and numerical work by <cit.>. As these systems evolve, the secondary can become more massive by accreting disk material or merge with other fragments so that it becomes massive enough to resist being dragged inward by the primary star. Further accretion onto the secondary, if it begins to grow faster than the primary, can then push the mass-ratio from an extremely low-value at the onset of disk fragmentation to q ≈ 1 at late times after the core has been exhausted <cit.> consistent with observations <cit.>. For example, <cit.> find that binary systems with initial extreme mass ratios, q ≪ 1 can be pushed to higher values of q as the binary system accretes material from the disk because the accretion onto the secondary can be as much as a factor of 10 larger than the primary star's accretion rate pushing the system to a higher mass ratio. However, we are unable to follow the long-term growth of the low-mass companion stars or further disk fragmentation into additional companion stars due to the high computational cost[The simulations presented in this work each took ≈ few ×10^5 CPU hours to run, with a significant fraction of the time spent on the final few kyr. This is because at late times the intense luminosity from the primary star yields such strong radiative accelerations such that the timestep, set by the Courant condition, decreases drastically thereby making advancing the simulation prohibitively expensive.] of the simulations presented in this work. 

Additionally, the growth of disk-borne stars is also sensitive to our merging criteria for star particles. In this work we only allow two sink particles to merge when the lower mass particle has a mass less than 0.04 M_⊙. Hence, this strict criterion for particle merging may lead to an over abundance of very low-mass companion stars as discussed in . This assumption in merging sink particles may lead us to predict a higher multiplicity for the massive system formed in run  as compared to  because we see that in Figure <ref>, many of the low-mass companion stars that are formed in the disk at late times would merge if our merging criteria had no mass dependence. A more lenient prescription of our merging criteria such as allowing sinks to merge if they pass within one accretion radius of each other (80 AU), as was done in <cit.>, would lead to a lower number of higher mass companion stars that may be more consistent with observations of massive multiple stellar systems. In light of these limitations, we are unable to conclude if such massive multiple stellar systems are formed quickly via disk fragmentation and accretion or via capture due to dynamical interactions in clustered environments. 

In the simulations presented here we treat the stars as accreting sink particles with an 80 AU accretion zone and therefore we are not able to follow if such objects are disrupted via shear motions on smaller scales in the accretion disk. For comparison, <cit.> performed similar simulations of disk fragmentation around massive protostars to study disk fragmentation with a much higher spatial resolution. By using a spherically symmetric domain in which the resolution increases logarithmically from the origin they attain sub-au resolution closest to the massive star and < 10 AU  resolution in the ∼ 500 AU region of their domain. They find that fragments that form via gravitational instability are not disrupted by shear motions in the accretion disk at these size scales and continue to grow via accretion like the companion stars formed in this work. Hence, we expect that the companion stars that form via disk fragmentation will likely survive the shear forces in the accretion disk unless the fragments, before sink creation, is significantly close to the star where the shear is high. 





 §.§ Caveats and Future Work


In this work, we have omitted two potentially important effects. The first is that massive cores are magnetized and therefore magnetic pressure may be dynamically important during core collapse. The second is that we neglect other feedback mechanisms such as magnetically launched collimated outflows and photoionization, which will also affect the growth rate of massive stars and core fragmentation. We also note that simulations presented in this work focus on the collapse of very dense massive star-forming cores (Σ = 1  g   cm^-3) and are therefore more representative of cores that sample the high density end of massive star formation since massive stars form in cores and clumps with surface densities that range from Σ∼ 0.1-1  g   cm^-3 <cit.>. Hence, our choice to neglect magnetic fields and chosen initial conditions are biased towards faster collapse time scales and higher accretion rates than those typically observed in most galactic massive star forming regions that have accretions rates of ∼ 10^-5 - 10^-4  M_⊙ yr^-1  and collapse time scales of ∼ 10^5 yr thereby yielding higher star formation efficiencies (SFEs) per free-fall time than those typically observed <cit.>. 

Magnetic support in cores can slow down their collapse, reduce fragmentation, and result in lower accretion rates onto stars thereby reducing the SFE of massive pre-stellar cores <cit.>. Observations of molecular clouds find that they are super-critical (μ_ϕ≳ 1 where μ_ϕ is the mass-to-flux ratio) concluding that gravity dominates over magnetic pressure in molecular clouds <cit.>. In addition, most observed massive dense cores are highly super-critical (μ_ϕ≫ 1). For example, <cit.> studied the evolved massive dense core DR 21(OH) and found μ_ϕ≈ 6 whereas its parent filament has a lower value of μ_ϕ≈ 3.4. Additionally, <cit.> studied several more cores in this filament and found that all cores are roughly or highly super-critical with μ_ϕ≈ 1-4.3). Hence, their results suggest that as filaments and clumps collapse to form dense cores gravity can become even more important making magnetic regulated collapse subdominant. 

In addition to gravity, turbulence can also dominate over magnetic pressure in massive cores and clumps. Both numerical simulations and observations have shown that on large scales (i.e., l ≳ 0.1-100 pc) the ISM is largely sub-Alfvénic whereas on smaller size scales, such as the size scales of dense cores (e.g., 0.01 ≲ l  ≲ 0.1 pc), molecular gas becomes super-Alfvénic with typical values of ℳ_ A≈ 2-3 suggesting that turbulence may dominate over magnetic fields in most massive star forming cores <cit.>. However, there are a few cases in which the massive core is super-critical but sub-Alfvénic. For example, <cit.> studied the massive c8 core (∼200   M_⊙) in G035.39-00.33 and found that the magnetic field structure of c8 has an hour glass morphology that arises from the field lines being dragged in as the core forms due to filament collapse. This effect increases the magnetic pressure in the inner core causing it to have an infall velocity smaller than the core envelope. Likewise, <cit.>  found that the magnetic field orientation in massive star forming regions ranges from ordered hour-glass configurations likely owing to strong magnetic fields to more chaotic distributions that are a result of weak magnetic fields. They also found that the magnetic field throughout the clump shows preferential alignment and aligns from the clump to the core scale. Their results suggest that at least qualitatively the field may be dynamically important. However, a more complete statistical sample is required to determine if cores are preferentially sub- or super-Alfvénic <cit.>. Strong magnetic fields will also affect clump fragmentation and protocluster formation. <cit.> concluded that cores that show fragments distributed in a filamentary-like structure are likely characterized by a strong magnetic field. Hence, future work requires including magnetic fields and simulating the collapse of sub- and super-Alfvénic cores to determine how the magnetic field and its relative importance to gravity and turbulence affects the growth rate of massive stars and the fragmentation properties of collapsing massive dense cores.



The simulations presented in this work also neglected how feedback from collimated outflows and photoionization affects the accretion history of massive stars and fragmentation of massive cores. Outflows are magnetically launched and present during the accretion phase for both low-mass and high-mass star formation with typical mass-loss rates of ∼10-30% of the accretion rate <cit.>. The presence of outflows will reduce the growth history of accreting massive stars yielding lower SFEs  for massive pre-stellar cores and therefore reduce the evolutionary effect of radiative heating in the collapsing core since the presence of outflows provides an avenue for the the radiative flux to escape <cit.>. Such an effect will lead to a higher degree of thermal fragmentation, potentially resulting in more companion stars to the massive star regardless of the initial virial state of the pre-stellar core.

Photoionization, on the other hand, will only become important when the star is sufficiently massive and compact enough to generate a large amount of photoionizing radiation to produce an ultra compact HII region and this occurs at M_⋆∼25-35 M_⊙ <cit.>. Recent RHD simulations of massive star formation by <cit.> that includes radiation pressure, outflows, and photoionization found that the larger thermal pressure from the resulting ionized gas (T = 10^4 K) can actually enhance accretion onto the star when photoionization first becomes important. They also show that photoionization feedback broadens the outflow cavities but does not limit disk accretion onto the star. Instead, they conclude that radiation feedback on dust at late times is responsible for reducing the gravitational infall of core material onto the circumstellar disk and eventually cuts off accretion. Hence, outflow and photoionization feedback likely does not shut off mass accretion onto massive stars but these feedback components will aid in reducing the accretion rate onto massive stars.


Finally, we note that modeling the collapsing cores as isolated objects is an oversimplification since observed massive cores and clumps are typically found in embedded, dynamical environments that may feed the cores' material as they collapse <cit.>. This type of environment may affect the fragmentation and collapse of such cores. However, the simulations presented in this work provide a crucial first step in understanding how the dynamic state of massive pre-stellar cores and stellar feedback affects the growth rate of massive stars and core fragmentation. In future work we will also address the concerns presented above by including magnetic fields and feedback from protostellar outflows in addition to radiation pressure to further understand how the dynamic state of the pre-stellar core affects the formation history and multiplicity of massive stars. 



§ CONCLUSIONS


In this work we performed 3D RHD simulations of the collapse of turbulent massive pre-stellar cores that are identical in every way except for their initial virial state to determine how the dynamic state of massive pre-stellar cores affects the formation history of massive stars, core fragmentation, and the multiplicity of massive stellar systems. To address these questions we modeled the collapse of a highly subvirial core with α_ vir = 0.14 and a roughly virialized core with α_ vir = 1.1. Our simulations included radiation feedback from both the direct stellar and dust-reprocessed radiation fields inherent to massive star formation owing to the massive star's large luminosity as it is actively accreting. We note that the simulations presented in this work neglect magnetic fields and outflow feedback that are ubiquitous in massive star forming regions. However, we will address these physical processes in future work to determine how these processes further affect the growth of massive stars and fragmentation of massive cores.

Our main results are summarized as follows. We find that subvirial cores undergo a fast global collapse leading to higher accretion rates, by up to a factor of ∼2, onto the massive star as compared to the accretion rates obtained from the gradual inside-out collapse of virialized cores.

Additionally, we find that virialized cores undergo significant turbulent fragmentation into companion stars at early times owing to their greater degree of turbulence (i.e., higher ℳ). The faster growth rate of massive stars that originate from less turbulent, subvirial cores leads to higher stellar luminosities at earlier times and therefore enhances the radiative heating from both the direct stellar and indirect dust-reprocessed radiation fields, thereby suppressing turbulent fragmentation. We also find that the global gravitational collapse of subvirial cores and gradual inside out collapse of virialized cores increases the velocity dispersion of massive cores in agreement with observations <cit.>. Furthermore, radiative feedback from massive stars also increases the averaged mass-weighted velocity dispersion and temperature of massive cores.  

Regardless of the core's initial virial state, we find that an optically thick accretion disk forms around the massive star and supplies material to the star, especially at late times when radiation pressure has driven radiation pressure dominated bubbles above and below the accretion disk. We also note that the accretion disk forms around the massive star earlier for virialized initial conditions because of  the core's larger angular momentum content and the disk is much larger in size at late times as compared to the accretion disk that forms around the massive star in our subvirial run. At late times, the accretion disk undergoes disk fragmentation forming close-in low mass companions regardless of the core's initial virial state yielding multiple systems. 

yt, <cit.>,  <cit.>,  <cit.>



 §.§ Acknowledgements

The authors thank the anonymous referee for their advice and suggestions which greatly improved  the manuscript. A.L.R. acknowledges support from NASA through Einstein Postdoctoral Fellowship grant number PF7-180166 awarded by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-03060. PSL acknowledges support by NASA through a NASA ATP grant NNX17AK39G. B.B. acknowledges support from the Simons Foundation. A.L.R would like to thank Alyssa Goodman, Phil Myers, and Alessio Traficante for insightful conversations regarding this work.
apj


