key: cord-0210107-v6b9hjnq
authors: Rosen, Anna L.; Krumholz, Mark R.
title: The Role of Outflows, Radiation Pressure, and Magnetic Fields in Massive Star Formation
date: 2020-06-08
journal: nan
DOI: nan
sha: 8ea08fdc3912c9d2e6cd0f8cce03386e69fb9858
doc_id: 210107
cord_uid: v6b9hjnq

Stellar feedback in the form of radiation pressure and magnetically-driven collimated outflows may limit the maximum mass that a star can achieve and affect the star-formation efficiency of massive pre-stellar cores. Here we present a series of 3D adaptive mesh refinement radiation-magnetohydrodynamic simulations of the collapse of initially turbulent, massive pre-stellar cores. Our simulations include radiative feedback from both the direct stellar and dust-reprocessed radiation fields, and collimated outflow feedback from the accreting stars. We find that protostellar outflows punches holes in the dusty circumstellar gas along the star's polar directions, thereby increasing the size of optically thin regions through which radiation can escape. Precession of the outflows as the star's spin axis changes due to the turbulent accretion flow further broadens the outflow, and causes more material to be entrained. Additionally, the presence of magnetic fields in the entrained material leads to broader entrained outflows that escape the core. We compare the injected and entrained outflow properties and find that the entrained outflow mass is a factor of $sim$3 larger than the injected mass and the momentum and energy contained in the entrained material are $sim$25% and $sim$5% of the injected momentum and energy, respectively. As a result, we find that, when one includes both outflows and radiation pressure, the former are a much more effective and important feedback mechanism, even for massive stars with significant radiative outputs.

Massive stars (M 8 M ) are rare, representing only ∼1% of the stellar population by number, yet they dominate energy injection into the interstellar medium (ISM) in star-forming galaxies because of their strong radiation fields, fast stellar winds, and supernova explosions. During their formation, massive protostars launch collimated bi-polar outflows that are qualitatively similar, but much more powerful, than those produced by low-mass protostars (Bontemps et al. 1996; Maud et al. 2015b; Bally 2016) . These similarities suggest that the driving mechanisms for such outflows originate from the same source, thereby providing evidence that low-and high-mass stars form in a similar fashion (Maud et al. 2015b) .

Massive stars form in dense (∼ 10 4 − 10 6 cm −3 ), cold (∼ 10 K), turbulent, and magnetized gas within giant molecular clouds and giant massive filaments (see detailed review by Rosen et al. 2020) . One of the key signatures of massive star formation, when the protostars are heavily embedded, are molecular outflows (e.g., Maud et al. 2015b; Pillai et al. 2019 ). These entrained outflows likely originate from collimated jets that are magnetically launched via the star-disk interaction (Shu et al. 1988; Pelletier & Pudritz 1992; Bontemps et al. 1996; Maud et al. 2015b; Kölligan & Kuiper 2018) . As these jets leave the star-disk system they encounter molecular material, which they sweep up, and this material can potentially be ejected from star-forming environments leading to low star formation efficiencies (SFEs, Cunningham et al. 2011; Maud et al. 2015a; Kuiper et al. 2016; Staff et al. 2019) . Additionally, the momentum injected by outflows may limit accretion onto the massive star, potentially setting an upper mass limit to the stellar initial mass function (IMF). Therefore, protostellar outflows are an important form of stellar feedback -the injection of momentum and energy into the surrounding ISM by young stars -during the star formation process.

Entrained molecular outflows from massive protostars are often observed via molecular line tracers (e.g., SiO, CS, or CO), methanol or water masers, and non-thermal synchrotron emission. They typically show complex orientations and morphologies with multiple shocks and broad density structures (e.g., de Villiers et al. 2014; Maud et al. 2015b; Brogan et al. 2018; Avison et al. 2019; Gieser et al. 2019; Sanna et al. 2019) . Observations of outflows from both low-and high-mass protostars indicate that the outflow and infall motions occur simultaneously and are closely linked. Typically, the outflow mass-loss rates are inferred to be 10-30% of the accretion rate onto the star. The outflow mass rate can also be highly variable, likely due to episodic accretion (Bachiller 1996; Matzner & McKee 2000; Burns et al. 2020; Nony et al. 2020 ).

Significant theoretical attention has focused primarily on the role that radiation pressure plays in massive star formation because massive stars have short Kelvin-Helmholtz time-scales (the time required for a star to radiate away its gravitational binding energy) and attain their main-sequence luminosities while they are still actively accreting (Palla & Stahler 1991 Behrend & Maeder 2001; Hosokawa & Omukai 2009; Krumholz et al. 2009; Kuiper et al. 2011; Rosen et al. 2016) . The radiation pressure associated with their high luminosities can oppose gravity and halt accretion. However, the rate of momentum deposition attributed to the direct radiation pressure from stars (i.e., first absorption of the stellar radiation by interstellar dust),ṗ rad = L c where L is the massive protostar's stellar luminosity, given bẏ

is typically lower than the rate of momentum deposition from protostellar outflows,ṗ OF =Ṁ OF v OF wherė M OF is the outflow mass-loss rate and v OF is the outflow velocity. This quantity is given bẏ 

where we parameterize the outflow mass flux and velocity asṀ OF = f wṀacc and v OF = f k v kep , whereṀ acc is the accretion rate, f w is the fraction of the accreted mass that is lost to outflows, v kep = GM /R is the Keplerian velocity of a star with mass M and radius R , and f k is the fraction of the Keplerian velocity with which the outflow is launched (Matzner & McKee 2000) . The values to which we have scaled in Equation 2 are typical for massive protostars (M , R , andṀ ) and magnetocentrifugal wind models (f w and f k ). Therefore, momentum feedback from protostellar outflows is a non-negligible feedback mechanism in massive star formation that must be included when studying how feedback limits accretion onto massive stars. Previous numerical studies have studied the effect outflows play in massive star formation in addition to radiation pressure. Cunningham et al. (2011) performed a series of 3D adaptive mesh refinement (AMR) radiationhydrodynamic (RHD) simulations with radiative and outflow feedback, where they modeled the dust reprocessed radiation field and used an ad-hoc prescription to treat the stellar radiation field from stars. In agreement with Krumholz et al. (2005) , they found that outflows evacuate polar cavities of reduced optical depth through the turbulent, ambient core. This effect enhances the radiative flux in the poleward direction so that radiative heating and the outward radiative force is diminished. Likewise, Kuiper et al. (2015 Kuiper et al. ( , 2016 performed cylindrically symmetric RHD simulations in which the star is held fixed (i.e., outflow launching location and outflow direction remained constant) but they include a hybrid treatment of radiation pressure that properly accounts for both the stellar and dust reprocessed radiation fields inherent to massive star formation. They found that outflows open a bipolar cavity extending to the outer edge of the protostellar core from which the simulation begins. The opening angles of the outflows and the amount of mass entrained and ejected from the core both increase with time. Additionally, they find that the importance of feedback from outflows depends on the amount of mass injected into the outflows at the point where they are launched, suggesting that outflows with a larger mass flux yields a star formation efficiency (SFE) from 50% in the case of very weak outflows to as low as 20% for very strong outflows.

These studies concluded that protostellar outflows makes radiation pressure less significant in massive star formation and causes molecular material to be ejected from the core leading to low SFEs. However, the simulations by Kuiper et al. (2015 Kuiper et al. ( , 2016 model the collapse of a laminar core with the star held fixed, so the outflow launching direction remains fixed. In reality, the accre-tion flow onto a massive protostar is chaotic because the accreted core material is turbulent, which in turn will cause the star's spin axis to precess. This effect should make the outflow launching direction highly variable, possibly leading to multiple outflows as are commonly observed in massive star forming regions (e.g., Gieser et al. 2019; Avison et al. 2019) . Additionally, magnetic fields are also not included in the Kuiper et al. 2015 , 2016 and Cunningham et al. 2011 simulations, and they will likely affect the outflow structure and ejection of material.

In this paper, we investigate these effects by performing 3D radiation-magnetohydrodynamics (RMHD) numerical simulations of the collapse of magnetized and unmagnetized turbulent massive prestellar cores into massive stellar systems, including both radiative and outflow feedback. This paper is organized as follows: we describe our numerical methodology and simulation design in Section 2; we present and discuss our results in Sections 3 and 4, respectively; finally, we conclude in Section 5.

In this paper, we simulate the gravitational collapse of isolated magnetized and non-magnetized turbulent massive pre-stellar cores with the Orion2 adaptive mesh refinement (AMR) code. Orion2 includes MHD (Li et al. 2012) , radiative transfer (Krumholz et al. 2007a; Shestakov & Offner 2008; Rosen et al. 2017) , self-gravity (Truelove et al. 1998) , and Lagrangian accreting sink particles (Krumholz et al. 2004 ) that include a protostellar evolution model used to represent them as radiating protostars ) coupled to a sub-grid prescription to model stellar feedback from protostellar outflows (Cunningham et al. 2011) . We describe the equations solved by our code in Section 2.1, our stellar radiation and outflow feedback prescriptions in Section 2.2, and the initial and boundary conditions, including our refinement and sink creation requirements, for our simulations in Section 2.3.

The full gravito-RMHD equations solved by Orion2 for the simulations that describe the dynamics of the fluid-sink (star) particle system presented in this work are

Here ρ is the density, ρv is the momentum density, ρe is the total internal plus kinetic gas energy density, E R is the radiation energy density in the rest frame of the computational domain, B is the magnetic field, and φ is the gravitational potential. Equations 3-6 describe conservation of gas mass, gas momentum, gas total energy, and radiation total energy, respectively. They include terms describing the exchange of these quantities with the star particles, and exchange of energy and momenta between radiation, magnetic fields, and gas. Equation 7 is the induction equation that describes the time evolution of the magnetic field in the ideal limit that assumes the magnetic field and fluid are well-coupled. Orion2 includes MHD using a constrained transport scheme (Li et al. 2012 ) that maintains ∇ · B = 0 to machine accuracy.

We assume an ideal equation of state so that the gas pressure is

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 . We assume the fluid is a mixture of gas and dust with a dust-to-gas ratio of 0.01. At the high densities that we are concerned with the dust will be thermally coupled to the gas, allowing us to assume that the dust and gas temperatures are the same. The radiation-specific quantities in Equations 5-6 are the Planck-and Rosseland-mean opacities κ 0P and κ 0R computed in the frame co-moving with the gas, the blackbody function B P = ca R T 4 /(4π), a dimensionless number λ called the flux limiter, and the Eddington factor R 2 . The last two quantities originate from the (gray) flux limited diffusion (FLD) approximation (see Krumholz et al. (2007a) and Rosen et al. (2016) for more detail). Lastly, Λ(T ) describes the cooling rate by atomic lines and the continuum, which only becomes significant when T 10 3 K when dust begins to sublime (Cunningham et al. 2011) .

We evolve the radiating (proto)star particles via Equations 8-10, indexed by the subscript i in the above equations. These particles accrete nearby gas and interact with the fluid via gravity, radiation, and protostellar outflows. We describe the modeling of their feedback (i.e., the momentum and energy injected into the fluid) associated with their radiation fields and outflows in Section 2.2, but note that the radiation and outflow specific terms in Equations 3-5 affiliated with star particles are denoted with the rad and o subscripts, respectively. The star particles are characterized by their position x i , momentum p i , mass M i , angular momentum that describes the particle's spin axis J i , and luminosityε rad,i (e.g., integrated stellar spectrum from Lejeune et al. 1997) , as determined by the protostellar evolution model described in Offner et al. (2009) . They accrete mass, momentum, and energy from the computational grid at ratesṀ a,i ,ṗ a,i , andε a,i . The distribution of these quantities over cells in the computational grid is described by a weighting kernel W a (x−x i ), which is nonzero only within 4 computational zones of each particle, following the algorithm of Krumholz et al. (2004) . We update the star particles' angular momentum and spin axis directions via the prescription described in Fielding et al. (2015) . The gravitational potential of the gas is advanced by solving Poisson's equation given by Equation 11, which includes contributions from both the fluid and star particles.

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. 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

where c s = kT /µm p is the isothermal sound speed, ∆x l is the cell size on level l, β = 8πρc 2 s /B 2 , and J max is the maximum allowed number of Jeans lengths per cell, which we set to 1/8, following the MHD Truelove Criterion (Myers et al. 2013);  (3) any cell that is located within at least eight cells of a sink particle; and (4) any cell within which the radiation energy density gradient exceeds ∇E R > 0.15E R /∆x l .

Star particles form when the Jeans condition for a Jeans number of N J = 0.25 is exceeded on the maximum AMR level following the resolution tests of Truelove et al. (1997) . They are allowed to merge when two star particles pass within one accretion radius of each other if the smaller particle has a mass less than 0.04 M , the threshold that corresponds to the largest plausible mass at which second collapse occurs for the protostar (Masunaga et al. 1998; Masunaga & Inutsuka 2000) . At masses lower than this value the protostar represents a hydrostatic core that is several au in size and will likely be accreted by the more massive star whereas larger mass protostars will have collapsed down to sizes of roughly several R and will unlikely merge with the nearby protostar.

Each star produces a (direct) stellar radiation field and collimated protostellar outflows that injects energy (ε) and momentum (p) into the fluid at a rate per unit volumeε rad,i ,ε o,i ,ṗ rad,i , andṗ o,i , where quantities subscripted by rad and o denote feedback from radiation and outflows, respectively. We use the multi-frequency Hybrid Adaptive Ray-Moment Method (HARM 2 ) described in Rosen et al. (2016) and Rosen et al. (2017) to treat both the direct (stellar) and indirect (dustreprocessed) radiation fields. This method combines direct solution of the frequency-dependent radiative transfer equation along long characteristics launched from stars to treat the direct stellar radiation field, including contributions from the accretion luminosity,

with a gray FLD method to treat the (indirect) radiation field produced by thermal emission from dust (Krumholz et al. 2007a; Rosen et al. 2016 Rosen et al. , 2017 . 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 ), and M and R are the star's mass and radius, respectively. We use the frequency-dependent stellar spectra and dust opacities from Lejeune et al. (1997) and Weingartner & Draine (2001) , respectively, and divide the stellar spectrum and dust opacities into 10 frequency bins (e.g., see Figure 1 in Rosen et al. (2016) ). We refer the reader to Krumholz et al. (2007a) and Rosen et al. (2016 Rosen et al. ( , 2017 for a detailed description of our treatment of the direct and indirect radiation pressures modeled in this work. Proper modeling of the magnetic launching of collimated protostellar outflows requires sufficiently highresolution (e.g., sub-au; Kölligan & Kuiper 2018) , which is prohibitively computationally expensive for the sim-ulations presented in this work. Instead, we adopt a sub-grid prescription for launching outflows from stars based on the protostellar outflow model of Matzner & McKee (2000) , first implemented by Cunningham et al. (2011) . In this model, the outflows are described by a collimation angle, θ c , and launching fraction, f w . This algorithm has the advantage that it can represent either a X-wind (Shu et al. 1988) or disk wind (Pelletier & Pudritz 1992) model. For the simulations presented in this work we adopt θ c = 0.01 and f w = 0.21, that assumes 21% of the accreted material is ejected in the outflows, and inject the outflows in the 8 nearest zones to the star with the weighting kernel W o,i (x − x i ) described in Cunningham et al. (2011) . The outflows are launched along the star's angular momentum (spin) axis at a fraction f k = 0.3 of the Keplerian velocity, such that their velocity is v o = f k GM /R . These parameter values are chosen to match observations of outflow momentum observed in low-and high-mass star formation (Cunningham et al. 2011) .

The outflows inject mass,

into the surrounding gas. To trace the outflow material we add a passively advected scalar to represent the outflow gas that is injected and we set T o equal to the protostar's surface temperature and µ w = 1.27, which is the mean molecular weight for a neutral gas of solar composition, since observations have shown that outflows from intermediate and massive protostars are predominately neutral (Reiter et al. 2016; Cesaroni et al. 2018; Fedriani et al. 2019) . When the massive (proto)star reaches a surface temperature 10 4 K we then assume the outflow is ionized and set T w = 10 4 K.

In this work, we perform three simulations of the collapse of turbulent, massive prestellar cores with feedback from stellar radiation and collimated outflows to determine how these feedback mechanisms effect the formation of massive stars: runs TurbRad, TurbRad+OF, and TurbRad+OFB. Run TurbRad 1 only includes radiative feedback, TurbRad+OF includes radiative feedback and collimated protostellar outflows, and TurbRad+OFB is identical to TurbRad+OF except that we also include magnetic fields. We use these simulations to compare how magnetic fields and feedback from collimated out-flows affect massive star formation and the growth rate of massive stars.

Following Rosen et al. (2016 Rosen et al. ( , 2019 , 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 massive prestellar core densities and radii in extreme massive star forming environments (e.g., Galván-Madrid et al. 2013; Battersby et al. 2014; Ginsburg et al. 2015 Ginsburg et al. , 2018 Contreras et al. 2018; Cao et al. 2019) . 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 (e.g., Caselli & Myers 1995; Beuther et al. 2002; Mueller et al. 2002; Beuther et al. 2007; Zhang et al. 2009; Longmore et al. 2011; Butler & Tan 2012; Battersby et al. 2014; Stutz & Gould 2016) . 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 in thermal pressure balance with the ambient dust-free medium and we set the opacity of the ambient medium to zero.

We explore the influence of magnetic fields on the collapse and outflow properties in run TurbRad+OFB. In this run, we set the initial magnetic field to be uniform in the z direction: B = B 0ẑ . We choose B 0 = 0.81 mG by selecting a mass-to-flux ratio

c B 0 is the magnetic flux through the core, consistent with observed values of Φ 2-3 (Crutcher 2012).

Observations of massive prestellar cores find that they contain supersonic turbulence and therefore we include supersonic turbulence by seeding the initial gas velocities (v x , v y , and v z ) with a velocity power spectrum, P (k) ∝ k −2 (Padoan & Nordlund 1999; Boldyrev 2002; Cho & Lazarian 2003; Kowal et al. 2007) . 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, consistent with the natural mixture of a 3D fluid (Kowal et al. 2007; Kowal & Lazarian 2010) . The onset of turbulence modifies the density and magnetic field distribution and we allow the turbulence to decay freely. For all runs, we use the same velocity perturbation power spectrum at initialization and a velocity dispersion of σ 1D = 1.2 km/s corresponding to a α vir = 5σ 2 1D R c /GM c = 1.1 so that the core is roughly virialized. We allow the turbulence to decay, which 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 (Goldreich & Sridhar 1995), is ∼0.16 Myr, which is much longer than the runtime for the simulations presented in this work. Our boundary conditions for the hydrodynamic, gravity, and radiation solvers are as follows. 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 3-5 (Cunningham et al. 2011; Myers et al. 2013; Rosen et al. 2016 Rosen et al. , 2019 and set the gravitational potential, φ, to zero at all boundaries when solving Equation 11 (Myers et al. 2013 ). We do not expect this choice of boundary conditions for the gravitational potential to lead to any significant square artifacts near the domain boundaries since the core boundaries are far removed from the domain boundaries. Finally, 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 Cunningham et al. 2011; Myers et al. 2013; Rosen et al. 2016 Rosen et al. , 2019 .

Here, we summarize the main results of our calculations. These simulations were run on the NASA supercomputer Pleiades located at NASA Ames. We use the yt package (Turk et al. 2011) to produce all the figures and quantitative analysis shown below.

3.1. Density Structure Figure 1 shows a series of density slices for runs TurbRad (top row), TurbRad+OF (middle row), and TurbRad+OFB (bottom row) at the same primary (most massive) stellar mass. Each panel is oriented so that the primary star's angular momentum (spin) axis is pointing up and the center of each panel corresponds to the location of the primary star. Each slice covers an area of (0.4 pc) 2 with the primary star at the center. We choose this orientation for the density slices because the collimated outflows are injected along the direction of the stellar spin axis, and thus the overall outflow structure is predominantly along this axis. Comparison of the overall core density structure for runs TurbRad+OFB and TurbRad+OF show that the entrained outflows break out of the core when the star reaches a TurbRad TurbRad+OF TurbRad+OFB Figure 2 . Zoom-in density slices for runs TurbRad (top row), TurbRad+OF (middle row), and TurbRad+OFB (bottom row). The most massive star is located at the center of each panel, as marked by the gray circle, and the slice is oriented such that the mass-weighted angular momentum axis of the gas within a radius of 250 au from the primary star points up in order to highlight the radiation pressure dominated bubbles that are perpendicular to the circumstellar gas due to the "flashlight effect" as described in the text. The primary stellar mass and the simulation time, in units of t ff , are shown in the bottom and top left corners of each panel, respectively. Each panel is (0.1 pc) 2 . mass of ∼ 30 M and that the outflows at breakout for run TurbRad+OF are more collimated than those in run TurbRad+OFB. The outflows in run TurbRad+OFB are surrounded by a lower density envelope of material. We discuss the outflow structure and energetics in more detail in Section 3.5. Comparison of these simulations with run TurbRad shows that the escaping outflows are a product of feedback from jets alone, because we do not see any outflow or bubble breakout from the core for run TurbRad.

Once the primary star in our simulations exceeds ≈ 30 M , feedback from radiation pressure drives lowdensity radiation pressure dominated bubbles that expand away from the massive star. We show this in Figure 2 , which shows zoom-in density slice plots for all three runs that cover an area of (0.1 pc) 2 . This effect is commonly referred to as the "flashlight effect", in which optically thick circumstellar material pinches the radiation field and beams it into the polar directions, driving low-density cavities that expand away from the star (e.g., Yorke & Sonnhalter 2002) . At equal primary stellar mass, the radiation pressure-dominated bubbles are more prominent in run TurbRad+OF (middle row) as compared to run TurbRad because outflows carve out low-density regions allowing the direct radiation pressure to be more effective at launching material at greater distances from the star. We find that these radiation pressure dominated low-density cavities are the least prominent in run TurbRad+OFB. However, at equal times, the radiation pressure dominated bubbles are most prominent in run TurbRad due to the faster mass growth and resulting larger luminosity of the primary star. Additionally, we find that the low-density cavities that develop in run TurbRad+OF are more collimated compared to the cavities in run TurbRad at equal primary mass.

The expanding radiation pressure-dominated regions develop because material near the star becomes super-Eddington, as shown in Figure 3 . As the central star gains mass, its luminosity increases, and the size of the super-Eddington regions near the star typically increase as well. We note that the super-Eddington region in run TurbRad+OF eventually decreases at late times (e.g., see last panel of the middle row in Figure 3 ). This result for run TurbRad+OF is likely due to the influence of outflows, which drive low-density channels near the star through which radiation can vent, an effect also seen by TurbRad TurbRad+OF TurbRad+OFB Figure 3 .

Slices of the Eddington ratio (f Edd = f rad /fgrav) for runs TurbRad (top row), TurbRad+OF (middle row), and TurbRad+OFB (bottom row). The most massive star is located at the center of each panel and the slice is oriented such that the mass-weighted angular momentum axis of the gas within a radius of 250 au from the primary star points up. The primary stellar mass and the simulation time, in units of t ff , are shown in the bottom and top left corners of each panel, respectively. Each panel is (0.25 pc) 2 . Cunningham et al. (2011) , and predicted theoretically by Krumholz et al. (2005) .

We also find that radiative feedback is less important when magnetic fields are present because magnetic confinement inhibits the expansion of the radiation pressure dominated bubbles: as the cavities expand their dense shells sweep up magnetic flux, amplifying the field strength in the shells (Krumholz et al. 2007b) . The increase in magnetic tension at the shells suppresses expansion, as shown in Figure 4 , which zooms in on the second to last panel of run TurbRad+OFB in Figure 2 and has the magnetic field vectors over-plotted to highlight the increased field strength in the bubble shells. The magnetic field structure also shows that the field lines are predominantly parallel to the shells, thereby opposing the radiation-pressure dominated bubble expansion.

In addition to radiation pressure, runs TurbRad+OF and TurbRad+OFB show that feedback associated with collimated protostellar outflows, which are present throughout the star formation process (i.e., from lowto high-masses), drive high-velocity entrained outflows whose opening angle and extent increase with time. The outflows eventually break out of the core and eject material when the primary star reaches ∼ 30 M . This behavior is more apparent in Figure 5 , which shows thin Density slice for the second to last snapshot in run TurbRad+OFB shown in Figure 2 with magnetic field vectors over-plotted. The area shown is (5,000 au) 2 . The massive star is located at the center of the panel, marked by the gray circle, and the slice is oriented so that the angular momentum axis of the gas within a radius of 250 au from the primary star points up. density projections of the radial momentum, p r = ρv r , with respect to the primary star. Negative values of p r denote material that is falling towards the star whereas positive values denote material that is moving away from the star. Comparison of Figures 2 and 5 show that the entrained outflowing material due to feedback from protostellar outflows is not necessarily aligned with the low-density cavities associated with the expanding radiation pressure dominated bubbles.

We show the properties of the primary protostar as a function of simulation time for runs TurbRad (teal solid lines), TurbRad+OF (pink dashed lines), and TurbRad+OFB (purple dot-dashed lines) in Figure 6 . The far left column shows the accretion rate (top panel) and primary stellar mass (bottom panel), demonstrating that the mass growth rate is fastest for run TurbRad since accreted material is not lost to outflows, as happens in runs TurbRad+OF and TurbRad+OFB. The growth rate is slowest for run TurbRad+OFB because, in addition mass loss by outflows, magnetic pressure slows down the gravitational collapse of the core. However, the overall accretion history for all three runs are similar in shape. The one difference we see between the runs is that the accretion rate decreases at late times for run TurbRad+OFB, which we run for the longest time, because feedback (aided by magnetic pressure) becomes sufficient to begin dispersing the core, thereby reducing the infall of material that can be accreted by the primary star. The accretion rate may also decrease because magnetic braking inhibits the formation of an optically thick accretion disk in run TurbRad+OFB, a possibility that we discuss in more detail in Section 3.4.

The accretion luminosity, stellar luminosity, and ratio of these two quantities are shown in the top middle left panel, bottom middle left panel, and bottom left panel, respectively. Early in all runs, when the star is less than ∼ several M , the accretion luminosity is larger than the stellar luminosity due to the high accretion rates inherent to massive star formation. As the primary star increases in mass and contracts to the zero age main sequence (ZAMS), as shown in the radial evolution plot in the top middle right panel, the stellar luminosity becomes larger than the accretion luminosity. Once this transition occurs, the accretion luminosity varies between ≈10%-100% of the stellar luminosity, suggesting that the accretion luminosity is non-negligible in massive star formation until late times when the accretion rate tapers off.

We show the outflow velocity and ratio of the rate of momentum deposition from outflows,ṗ OF =Ṁ OF v OF , to the rate of momentum deposition from direct radiation,ṗ rad = (L + L acc )/c, in the top far right and bottom far right panels, respectively. The outflow velocity weakly increases until t ∼ 0.5t ff because the primary protostar's radius and surface escape speed gradually increase. However, once the star begins to contract to the ZAMS, the outflow velocity increases rapidly. Throughout the majority of the accretion history the rate of momentum deposition by outflows is much larger than the rate of momentum deposition by radiation. At early times this ratio is of the order of ∼ few ×100, and it eventually decreases to a factor of a few by the end of run TurbRad+OF and a factor of ∼ 1 by the end of run TurbRad+OFB. Therefore, we find that the momentum input by outflows dominates over the momentum input by radiation throughout the majority of the star formation process.

We show the evolution of the primary star's position and spin axis in Figure 7 as a function of primary stellar mass, respectively. The top row of this figure shows the evolution of the primary star's spin axis in spherical coordinates: θ j = arctan (ĵ y /ĵ x ) (left panel) and φ j = arccos (ĵ z ) (right panel) whereĵ = (ĵ x ,ĵ y ,ĵ z ) is the unit vector describing the direction of the primary star's spin axis in Cartesian coordinates. This figure demonstrates that the momentum and angular momentum accreted by the star cause the star to move, and lead the primary star's spin axis to precess, especially at early times when the star is low in mass (i.e., M 10M ). We describe the impact this has on the outflow structure in more detail in Section 3.5.

The bottom left panel of Figure 7 show the rate of change of the angle traced by the primary star's spin axis defined as

where we have divided this quantity by 2π to convert from radians to revolutions where one revolution corresponds to the spin axis precessing by 360 • . This panel in Figure 7 shows that the precession decreases and dψ j /dt flattens once the star has a mass 10 M . Thereafter the rate of spin precession stays relatively constant and low until accretion begins to taper off, at which point precession stops as well. Figure 8 shows a series of density slices of the accretion disk that forms around the massive star with velocity streamlines over-plotted. Here, we compare each TurbRad TurbRad+OF TurbRad+OFB Figure 5 . Slices of the gas radial momentum with respect to the primary star with a diverging color scale chosen to highlight material that is moving toward (negative values of ρvr) and away from (positive values of ρvr) the star for runs TurbRad (top panels) and TurbRad+OF (bottom panels) at different times. Each panel is (0.4 pc) 2 in area and the gray star at the center of each panel denotes the location of the primary massive star. The slices are oriented so that the angular momentum axis of the protostar points up to highlight the radial momentum from the outflows. simulation at equal times since the accretion disk structure depends on the angular momentum content of the collapsing core. This figure shows that a noticeable highdensity accretion disk (i.e., a resolved accretion disk with a radius larger than the 80 au accretion zone radius of the sink particle) forms around the primary stars in runs TurbRad and TurbRad+OF. However, we do not see a resolved accretion disk in run TurbRad+OFB. This difference is almost certainly a result of magnetic braking. The accretion disk grows in size as runs TurbRad and TurbRad+OF progress owing to conservation of angular momentum: as time advances, the material reaching the vicinity of the primary star began its infall from greater distances, and thus has larger net angular momentum and therefore will be circularized at a distance farther from the star. Magnetic fields suppress this effect by carrying angular momentum away from the infalling material and farther out into the surrounding core suppressing disk growth (Seifried et al. 2011; Commerçon et al. 2011; Myers et al. 2013) . Nonideal effects, which remains an active field of study in protostellar disk formation (e.g., see detailed reviews by Wurster & Li 2018; Zhao et al. 2020) , may reduce the effect of magnetic braking leading to small accretion disks or toroids or mitigate the magnetic braking catastrophe entirely (Kölligan & Kuiper 2018; Wurster et al. 2019 ). However, we do not include these effects, and, even if they were included, it is still possible that any resulting disk would be unresolved in our simulation.

We do note that at early times for run TurbRad+OFB in Figure 8 a small disk-like structure forms around the massive star, but this disappears later in the simulation. Instead, radiation pressure near the star yields a magnetically confined low-density bubble surrounding the star that causes the accretion rate to drop at late times. We show this structure in Figure 9 , which shows the density structure near the primary star near the end of run TurbRad+OFB with magnetic field vectors over-plotted.

Thus far we have shown that outflows have a decisive effect: they deliver far more momentum to a protostellar core than radiation pressure, they entrain a significant amount of mass, and they reduce the influence of radiation pressure by providing channels through which radiation can escape. In this section we explore the structure and properties of the outflows in more detail.

We show projections of the material entrained by outflows for run TurbRad+OF (top row) and run TurbRad+OFB (bottom row) in Figure 10 . We define entrained material as consisting of all cells whose mass contains at least 5% of the launched material (i.e., cells where f t = ρ OF /ρ ≥ 0.05); recall that we add a passively advected scalar to the outflow material we inject, which allows us to measure ρ OF precisely for each cell. We only include gas that has a positive radial velocity, v r > 0, with respect to the primary star and subtract the launched material from the total density so that we only include contributions of entrained material. We note that our definition does include contributions from the low-mass companion stars as well as the primary star since we are not able to trace the ejected quantities for individual stars. However, this should only be a minor effect since most of the injected outflow mass is from the primary star. Figure 10 shows that the outflow structure is not steady. At early times there appear to be multiple outflows present. In reality the primary star dominates outflow production at all times, but, as discussed in Section 3.3, the star's angular momentum axis changes rapidly when its mass is low, causing outflows to be launched in multiple directions. The precession eventually decreases as the star increases in mass leading to a steadier outflow. Eventually the opening angles of the TurbRad TurbRad+OF TurbRad+OFB Figure 8 . Density slices in runs TurbRad (top row), TurbRad+OF (middle row), TurbRad+OFB (bottom row) shown at equal times; the slices have been oriented so that the angular momentum axis of the material within 250 au of the primary star points out of the page, in order to highlight the accretion disk. In each panel, the most massive star is at the center, and the region shown around it is (2000 au) 2 in size. Velocity streamlines and companion stars with masses greater than 0.04 M are over-plotted on all panels. The color of the star indicates its mass, as shown in the colorbar. Density g cm 3 Figure 9 . Similar as the bottom right panel of Figure 8 but for the snapshot shown in Figure 4 , and with magnetic field vectors rather than velocity streamlines over-plotted.

entrained outflows broaden such that the multiple outflows from the massive star merge. This effect leads to wider outflows as shown in the middle column of Figure 10 ; the outflow broadening is much larger in run TurbRad+OFB as compared to run TurbRad+OF. Hence, magnetic fields produce lower density and broader, lesscollimated outflows.

We quantify the outflow broadening in Figure 11 , which shows the volume filling fraction of the outflows and entrained gas for runs TurbRad+OF (pink dashed lines) and TurbRad+OFB (purple dot-dashed lines) as a function of time (top panel) and primary stellar mass (bottom panel). We calculate this value by summing over the volume of all cells whose mass contain launched material normalized to the initial core volume: f OF, V = i dV OF,i /V core, init where V core, init = 4 3 πR 3 c . We find that, as a function of simulation time, the volume filling fractions are roughly the same for both simulations and they slowly increase up to t ≈ 0.8 t ff . After this time, the entrained volume fraction increases much more rapidly for both runs. The rapid increase corresponds to the point where outflows begin to break out of the the initial core and expand freely into the low-density medium outside it. However, it is interesting to notice that, despite their similarity in time, there is a significant difference between the runs when we study the behavior in terms of primary star mass: the upturn in the volume filling fraction occurs when the star reaches ∼ 20M in run TurbRad+OFB as compared to ∼ 30M in run TurbRad+OF. This suggests that less momentum is required for outflow breakout in run TurbRad+OFB, likely because magnetic levitation reduces gravitational confinement of the core (Shu et al. 2004 ). We show the launched outflow mass versus the entrained outflow mass in Figure 12 for runs TurbRad (pink dashed line) and TurbRad+OFB (purple dashed line), in which we consider that material is entrained if the cell contains at least 5% of the launched outflow material. Entrained outflow versus injected mass for runs TurbRad+OF (pink dashed lines) and TurbRad+OFB (purple dot-dashed lines). The shaded pink and purple regions denote the entrained outflow mass where we consider cells whose density contain between 1-10% of outflow material for runs TurbRad+OF and run TurbRad+OFB, respectively. The gray line has a slope of three (i.e., MOF, ent = 3MOF, inj) .

The shaded regions denote the spread in the total entrained outflow mass considering cells that contain 1-10% of the launched outflow material. The over-plotted gray solid line has a slope of 3 showing that the entrained outflow mass is a factor of ∼ 2−3 larger than the launched outflows mass. Hence, we find that feedback from collimated outflows and radiation pressure have a mass-loading factor of ∼ 2 − 3 in massive star formation, so the total amount of mass directly removed is ∼ 50% of the core (since the outflow at launch contains ∼ 20% of the accreted mass). Additionally, Figure 12 shows that the entrained mass is quantitatively similar regardless if the core is magnetized or not. The shaded regions in this Figure in total entrained mass occur near the end of the simulation when we consider a low (e.g., 1%) fraction of outflow to entrained material, likely due to advection of the launched material. However, a larger fraction (e.g., 10%) yields quantitatively similar results of the total mass entrained as compared to when we consider cells that contain ≥ 5% of launched material.

We also compare the entrained outflow momentum and kinetic energy to the injected outflow momentum and kinetic energy in Figure 13 . We define the entrained outflow momentum and kinetic energy as

and

respectively, where ρ ent, i is the entrained outflow mass density, dV i is the cell volume, and v i is the gas velocity magnitude of cell i. Again, we remind the reader that we count a cell as entrained if it is at least 5% outflow ma-terial by mass, and if the radial velocity, with respect to the primary star, is moving away from the primary star. Also note that these outflow quantities include contributions from both the primary and companion stars. However, this should make little difference since the injected and entrained outflow momentum and energy is dominated by the primary star regardless. For comparison, we compute the time-integrated injected outflow momentum and kinetic energy from the stellar properties (see Figure 6 ) up to time t as p OF, inj (t) = t 0Ṁ OF v OF dt and E OF, inj (t) = t 0 1 2Ṁ OF v 2 OF dt, respectively. We find that the entrained material contains ∼ 25% of the injected momentum and ∼ 5% of the injected kinetic energy. At late times, we also find that the momentum and kinetic energy contained within the entrained outflows are quantitatively independent of the core's magnetic field strength. The reduction in the entrained outflow momentum compared to that which was originally injected is likely due to mixing between the outflow material and ambient gas that has a negative radial momentum as a result of gravitational infall. Similarly, the reduction in energy is due to inelastic entrainment of material coupled with fast radiative cooling, since the cooling time is always short compared to any mechanical timescale. The fact that the reduction in energy is larger than the reduction in momentum suggests that outflows should be considered a momentum-driven rather than an energy-driven feedback mechanism.

Companion stars form via turbulent fragmentation, in which over-densities in the surrounding core envelope collapse to form stars; and via disk fragmentation, in which the accretion disk becomes gravitationally unstable and fragments (Kratter & Matzner 2006; Rosen et al. 2019) . We show the total companion stellar mass (top row) and number of companion stars (bottom row) as a function of simulation time (left column) and stellar mass (right column) for runs TurbRad (solid teal lines), TurbRad+OF (dashed pink lines), and TurbRad+OFB (dot dashed purple lines) in Figure 14 . This Figure shows that run TurbRad+OFB does not form any companion stars because magnetic pressure suppresses turbulent fragmentation and additionally a resolved gravitationally unstable accretion disk does not form around the primary star (e.g., Commerçon et al. 2011) .

We find that from t ≈ 0.35 − 0.6 t ff companion stars form via turbulent fragmentation and at late times they form via disk fragmentation when the star is very massive for runs TurbRad and TurbRad+OF (Rosen et al. 2019) . For run TurbRad+OF more companion stars form via turbulent fragmentation as compared to run TurbRad because the primary star is less luminous due to its slower growth, thereby reducing radiative heating of the core material. In addition, outflows allow venting of the radiation field making radiative heating less effective. At late times we see that run TurbRad forms more companion stars via disk fragmentation because the star is much more massive than the primary star in run TurbRad+OF and becomes gravitationally unstable (e.g., see Figure 10 from Rosen et al. (2019) ). Given that we do not form companion stars via turbulent fragmentation in run TurbRad+OFB, but we form many in run TurbRad+OF we expect that weaker magnetic fields (i.e., cores with µ Φ > 2) should lead to a weaker degree of turbulent fragmentation (e.g., Palau et al. 2013; Myers et al. 2013; Fontani et al. 2018 ).

The purpose of this work is to understand how the interplay between magnetic fields and feedback from radiation pressure and collimated outflows effect the formation of massive stellar systems and affects the resulting entrained outflow structure from massive stars. Most notably, we find that momentum feedback from collimated outflows dominates over radiation pressure in massive star formation and ejects a significant fraction of molecular material from the core leading to low star formation efficiencies (SFEs; = M , tot /M c where M , tot is the total stellar mass and M c is the initial core mass).

Additionally, the presence of magnetic fields slows the growth rate of massive stars, reduces core fragmentation, and leads to broader and lower density outflows as compared to the outflows that emanate from unmagnetized cores.

In what follows, we discuss how magnetic fields and feedback from outflows and radiation pressure affects the growth rate of massive stars and the overall SFEs of massive cores in Section 4.1. Next, we discuss how feedback from outflows and radiation pressure affect the energetics of massive star forming cores in Section 4.2.

Star formation is an inefficient process, likely as a result of stellar feedback, turbulence, and magnetic fields (e.g., Matzner & McKee 2000; Offner et al. 2014; Louvet et al. 2014) . Observations of the stellar initial mass function (IMF) and prestellar core mass function (CMF) in low-mass star formation studies find a direct mapping between these two quantities in which the CMF is similar in shape to the IMF but offset by a factor ∼ 3; yielding a core SFE of ∼ 33%. This finding suggests that magnetic fields and feedback from outflows may be responsible for the mass offset since radiation pressure is unimportant in low-mass star formation (Alves et al. 2007; Offner & Chaban 2017) .

Given the rarity of and large distances to massive star forming regions, and the fact that low-mass star forma- The gray dotted line in both panels denotes where = 0.33 and = 0.4 as determined by the core mass function for low-mass star formation by Alves et al. (2007) and the SFE estimate for high-mass protocluster formation from Maud et al. (2015a), respectively. tion likely occurs coevally with high-mass star formation in massive cores, it remains uncertain what SFEs are expected for massive star formation and which physical properties set the SFE of massive star forming cores (e.g., Motte et al. 2018; Pillai et al. 2019) . Maud et al. (2015a) studied the core properties and SFEs of a sample of 89 distance limited (d 6 kpc) cores that hosted massive young stellar objects (MYSOs) and compact H ii regions. They found that the mass-luminosity plane of these sources is consistent with the luminosity expected from a proto-cluster that hosts at least one high-mass source and forms with a ∼40% SFE, slightly larger than the value inferred from low-mass studies.

Here, we use our simulations to determine the role that radiative and outflow feedback play in setting the SFE of magnetized and unmagnetized massive star forming cores and compare them to the observed values de-scribed above. We show the total SFE for runs TurbRad (teal solid line), TurbRad+OF (pink dashed line), and TurbRad+OFB (purple dot dashed line) in the top panel of Figure 15 and the SFE considering only the primary star (i.e., = M , p /M c ) in the bottom panel. Comparison of these two SFEs (total versus primary) allow us to determine how magnetic fields and feedback from radiation pressure and/or outflows affects the mass growth of the primary star and subsequent growth of the companion stars. We note that the SFE plotted for run TurbRad+OFB in Figure 15 is identical in both panels since this run does not form companion stars throughout the simulation run time.

The top panel in Figure 15 show that the SFE is largest for run TurbRad, reaching ∼ 0.45 at t = 0.95 t ff ; this is larger than the ≈ 0.33 and ≈ 0.4 (marked by the dashed gray lines) expected for low mass star formation (Alves et al. 2007 ) and high-mass protocluster formation (Maud et al. 2015a) , respectively. This high SFE occurs because accreted mass is not lost to outflows, while radiation pressure alone is not sufficient at ejecting material from the core, at least for the simulation run time considered here. Additionally, the SFE for run TurbRad continues to increase rapidly at late times suggesting it should further increase. These findings suggest that radiation pressure alone is not responsible for the SFEs observed in massive star forming environments.

Next, we find that the total SFE in the runs containing outflows are significantly lower than the SFE in run TurbRad but that run TurbRad+OF is larger than run TurbRad+OFB, reaching a value of ∼ 0.3 as compared to ∼ 0.2 for run TurbRad+OFB at t = 0.95 t ff . The higher SFE for the non-magnetic core is attributed to the increase of stellar mass for the companion stars that are still accreting from core material at the end of run TurbRad+OF and the slightly higher growth rate of the primary star. Hence, we find that outflows are required to set the SFEs expected for massive star formation and that the presence of strong magnetic fields leads to even lower SFEs. However, weaker magnetic fields and/or non-ideal MHD effects such as Ohmic resistivity, ambipolar diffusion, and the Hall effect, which we do not include in this work could lead to a higher degree of fragmentation and therefore higher SFEs in magnetized massive cores as seen in the numerical study of Fontani et al. (2018) and the observational study of Palau et al. (2013) . The importance of these non-ideal effects remains substantially uncertain. Radiative effects should increase the ionisation fraction near massive protostars, and thus reduce the importance of non-ideal effects. If non-ideal effects do allow greater fragmentation than we find in run TurbRad+OFB, this would likely lead to increase SFEs, perhaps close to those found by Maud et al. (2015a) .

When we only take into account the primary mass for calculating the SFE, as shown in the bottom panel, we find the SFE continues to increase at an accelerating rate by the end of run TurbRad but tapers off at the end of run TurbRad+OF and TurbRad+OFB, suggesting that outflows, rather than radiation pressure, are required to shut off accretion onto the massive star at late times. Comparison with the top panel, which shows that the SFE for the total system in runs TurbRad and TurbRad+OF continues to increase when these simulations end, demonstrates that accretion onto the low mass companion stars leads to increased SFEs even when feedback from outflows cuts off the accretion flow onto the primary star. Our result agrees with the observational study presented in Pillai et al. (2019) , which studied the protostellar content via outflow signatures in two infrared dark clouds, that found that low-mass protostars likely form co-evally at the earliest phase of high-mass star formation (i.e., t 50, 000 yrs).

As we demonstrated throughout this paper, collimated outflows that originate from accreting protostars dominate the momentum budget for stellar feedback in massive star formation, while radiation pressure is secondary. The momentum from these outflows entrains nearby molecular core material that may eventually be ejected from the core, leading to low SFEs. Studies have observed entrained outflows emanating from MYSOs via methanol masers and molecular line tracers in order to to determine the outflow energetics (e.g., the entrained outflow mass, momenta, energy, and force) (e.g., de Villiers et al. 2014; Maud et al. 2015b ). These studies have concluded that these quantities scale with the luminosity of the central driving source, indicating that the observed outflows are powered by a common, scalable driving mechanism similar to outflows observed in low-mass star formation (Bontemps et al. 1996; Maud et al. 2015b; Bally 2016) .

The simulations presented in this work allow us to follow the evolution of the the outflow energetics and determine how they depend on source luminosity, and how the relationship between these two quantities in our simulations compares to the observed one. In light of this goal, we show the entrained outflow mass, momenta, energy, and force (i.e., momentum flux) as a function of the massive protostar's luminosity taken directly from runs TurbRad+OF (pink dashed lines) and TurbRad+OFB (purple dot dashed lines), respectively, in Figure 16 . This figure is constructed to be similar to that presented in Maud et al. (2015b) . We find that the outflow mass, momentum, energy, and force are roughly similar as a function of source luminosity regardless of whether the core is magnetized or not, and that the relationship between outflow force and source luminosity in our simulations agrees well with that determined by Maud et al. (2015b) .

In the bottom right panel, where we show the outflow force as a function of source luminosity, we overplot the power law fits derived from observations of outflows from low-mass protostars, log 10 F = −5.6 + 0.9 × log 10 L(L ) (Bontemps et al. 1996 , gray dotted line), and massive protostars, log 10 F = −4.8 + 0.61 × log 10 L(L ) (Maud et al. 2015b , gray solid line). We find that the outflow force inferred from our simulations agree well with the fit derived by Maud et al. (2015b) . In their work, Maud et al. (2015b) conclude that the most massive protostars in the clusters are responsible for the energetics of the observed outflows, thereby leading to the shallower slope in their derived fit for the outflow force as compared to the steeper fit determined by Bontemps et al. (1996) for outflows from low-mass protostars. Our results are consistent with this conclusion.

In this work, we performed a series of 3D RHD and RMHD simulations of the gravitational collapse of isolated dense massive prestellar cores to determine how magnetic fields, turbulence, and radiative and collimated outflow feedback affects the formation of massive stellar systems. This is one of the first studies of massive star formation to include all of these effects, in the context of a realistic, turbulent medium without artificially-imposed geometries or symmetries (e.g., stars that are fixed at the origin of a symmetric cloud). By following the material launched by the outflows we have investigated outflow entrainment and the momentum and energy injection by outflows in massive star formation with and without magnetic fields. We reach the following conclusions:

1. Feedback from outflows dominates over momentum injection by radiation pressure in massive star formation.

2. The accretion and stellar luminosities are comparable throughout the star formation process, therefore accretion luminosity should not be neglected in massive star formation studies.

3. Strong magnetic fields suppress fragmentation and preferentially lead to monolithic collapse to single stars. Weaker magnetic fields and non-ideal effects Maud et al. (2015b) , and for low-mass protostars by Bontemps et al. (1996) , respectively.

might weaken this effect, and lead to the formation of a small number of companion stars.

4. The presence of magnetic fields leads to broader and lower density entrained outflows as compared to outflows in non-magnetic cores.

5. Magnetic fields and outflows are required to produce the low SFEs commonly observed in star formation. We find that radiative feedback alone does not lead to SFEs as low as those commonly observed in massive star formation.

6. We find that accretion continues onto the low-mass companion stars even after outflow feedback has significantly reduced accretion onto the massive primary star.

7. By tracking the launched and entrained outflows we calculated the mass-loading factor and momentum-and energy-scaling factors associated with feedback from protostellar outflows. We find that the mass-loading factor roughly falls between ∼ 2 − 3, but that the amounts of momentum and energy carried by the outflow are smaller than those injected by the stars, by factors of ∼ 0.25 and ∼ 0.05, respectively. Our simulations reproduce the relationship between source luminosity outflow force observed for massive protostars by Maud et al. (2015b) .

8. Our results for the momentum and energy contained in the entrained outflows suggests that outflows should be considered a momentum-driven rather than an energy-driven feedback mechanism since the reduction in energy is larger than the reduction in momentum.

Council through its Future Fellowship and Discovery Projects funding schemes (awards FT180100375 and DP190101258), and from the Australian National Computational Infrastructure (NCI, project jh2). A.L.R. would like to thank Stella Offner, Andrew Cunningham, Qizhou Zhang, and Alyssa Goodman for insightful conversations regarding this work. A.L.R. and M.R.K. would also like to thank their "coworkers," Nova Rosen and Jeff Anderson-Krumholz, for "insightful conversations" and unwavering support at home while this paper was being written during the Coronavirus pandemic of 2020. Their contributions were not sufficient to warrant co-authorship due to excessive napping. 2

2 Because they are cats

Protostars and Planets VI

The authors thank the anonymous referee for their advice and suggestions which 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. M.R.K. acknowledges support from the Australian Research