{'fraction_non_alphanumeric': 0.038968776854769405, 'fraction_numerical': 0.020876539673446003, 'mean_word_length': 4.333740756584978, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Impact of Feedback During Massive Star Formation by Core Accretion

We study feedback during massive star formation using semi-analytic methods, considering the effects of disk winds, radiation pressure, photoevaporation and stellar winds, while following protostellar evolution in collapsing massive gas cores. We find that disk winds are the dominant feedback mechanism setting star formation efficiencies (SFEs) from initial cores of ~0.3-0.5. However, radiation pressure is also significant to widen the outflow cavity causing reductions of SFE compared to the disk-wind only case, especially for>100Msun star formation at clump mass surface densities Sigma<0.3g/cm2. Photoevaporation is of relatively minor importance due to dust attenuation of ionizing photons. Stellar winds have even smaller effects during the accretion stage. For core masses Mc~10-1000Msun and Sigma~0.1-3g/cm2, we find the overall SFE to be 0.31(Rc/0.1pc)^{-0.39}, potentially a useful sub-grid star-formation model in simulations that can resolve pre-stellar core radii, Rc=0.057(Mc/60Msun)^{1/2}(Sigma/g/cm2)^{-1/2}pc. The decline of SFE with Mc is gradual with no evidence for a maximum stellar-mass set by feedback processes up to stellar masses of ~300Msun. We thus conclude that the observed truncation of the high-mass end of the IMF is shaped mostly by the pre-stellar core mass function or internal stellar processes. To form massive stars with the observed maximum masses of ~150-300Msun, initial core masses need to be>500-1000Msun. We also apply our feedback model to zero-metallicity primordial star formation, showing that, in the absence of dust, photoevaporation staunches accretion at ~50Msun. Our model implies radiative feedback is most significant at metallicities ~10^{-2}Zsun, since both radiation pressure and photoevaporation are effective in this regime.


INTRODUCTION
Massive stars play important roles in a wide range of astrophysical settings. They are the sources of UV radiation, turbulent energy and heavy elements. Massive star close binaries are the likely progenitors of merging black hole systems that have been the first sources to be detected by their gravitational wave emission. However, massive star formation is relatively poorly understood compared to low-mass star formation (see Tan et al. 2014, for a recent review). One class of models of massive star formation is based on the Core Accretion scenario (e.g., the Turbulent Core Model of McKee & Tan 2003). These models are scaled-up versions of models of low-mass star formation from cores (e.g., Shu, Adams & Lizano 1987), invoking nonthermal forms of pressure support, i.e., turbulence and magnetic fields to help stabilize the initial massive pre-stellar core. However, there may also be significant differences compared to low-mass star formation due to the stronger feedback that is expected from massive protostars.
In low-mass star formation, the magnetohydrodynamic (MHD) outflow is thought to be the main feedback process, which may determine the star formation efficiency (SFE) from the pre-stellar core, i.e.,ε * f ≡ m * f /M c where m * f is the final mass that is achieved by the protostar at the end of its accretion and M c is the mass of the initial core. In relatively low-mass clusters that contain stars with masses up to ∼ 10M ⊙ , the core mass function (CMF) is reported to be similar in shape to the stellar initial mass function (IMF), but shifted to higher masses by a factor of a few (e.g., André et al. 2010;Könyves et al. 2010). One explanation for this is a nearly constant SFE as a function of core mass of aboutε * f ∼ 0.4. Matzner & McKee (2000) proposed that the accretion-powered, MHD-driven outflow sets the SFE from pre-stellar cores. They provided an analytic model and showed that the momentum injected by the disk wind sweeps up a certain fraction of material in the infalling envelope and sets a SFE of ∼ 0.3-0.5. The numerical simulation by Machida & Matsumoto (2012) confirmed this result obtaining a similar value of SFE. Therefore, in low-mass star formation, observations and theoretical models are in agreement that an individual star can be formed by collapse of a pre-stellar core with the MHD outflow setting a SFE of ∼ 0.4.
In massive star formation, additional feedback processes may become more significant than the MHD out-flow because of the high luminosities of massive stars. Especially, radiation pressure has been considered to be a potential barrier for massive star formation. In an idealized spherical geometry, radiation pressure acting on dust grains in an infalling envelope exceeds the gravitational force when the stellar mass reaches ∼ 10-20 M ⊙ preventing further mass accretion (Larson & Starrfield 1971;Wolfire & Cassinelli 1987). The fact that more massive stars exist tells us that the model of spherical infall is too simplified. Subsequent work on analytic and semi-analytic models (e.g., Nakano 1989;Jijina & Adams 1996;Tanaka & Nakamoto 2011) and numerical simulations (e.g., Yorke & Sonnhalter 2002;Krumholz et al. 2009;Kuiper et al. 2010;Rosen et al. 2016) of disk accretion found that mass infall and accretion can continue from behind the disk since this region is shielded from strong radiation pressure. The series of simulations by Kuiper and collaborators have shown that disk accretion continues while the direct stellar radiation sweeps up the material above the disk where the shielding effect is weak (Kuiper et al. 2010(Kuiper et al. , 2011(Kuiper et al. , 2012(Kuiper et al. , 2015. They found the SFE from 100M ⊙ -cores is about 0.5 in models without MHD disk wind feedback. The recent simulation with high resolution and moving sink particle method by Rosen et al. (2016) showed that the Rayleigh-Taylor (RT) instability strongly helps to bypass the radiation pressure barrier even above the disk. Also, if MHD outflow cavities exist before radiation pressure becomes significant, then radiation leaks away via these channels, i.e., enhancing the so-called "flashlight effect" (Yorke & Bodenheimer 1999;Yorke & Sonnhalter 2002;Krumholz et al. 2005;Kuiper et al. 2015Kuiper et al. , 2016. Thus, the radiation pressure barrier is not thought to be a catastrophic problem anymore for massive star formation. Rather the question now is what is its quantitative effect on the formation efficiency of massive stars from massive cores. Photoionization may also be a significant feedback process. When a massive protostar approaches the Zero-Age Main-Sequence (ZAMS), it contracts, increases its effective temperature and starts to emits significant fluxes of Lyman continnum photons with > 13.6eV that may ionize the infalling/accreting material. Such ionized gas has a high temperature of ∼ 10 4 K and its thermal pressure may drive mass-loss in a "photoevaporative" outflow. In the formation of primordial (Pop III) stars in the early universe, photoevaporation is thought to be significant, potentially stopping mass accretion at ∼ 50-100M ⊙ Hosokawa et al. 2011;Tanaka et al. 2013). Note that radiation pressure feedback is not very significant in primordial star formation since there are no dust grains. Coincidentally, the typical mass accretion rates in primordial star formation and in present-day massive star formation are expected to be similar, with values of ∼ 10 −3 M ⊙ yr −1 . Thus, one may speculate that photoevaporation also stops mass accretion in present-day massive star formation. However, the dependence of the photoevaporation rate on metallicity has not been studied very much and remains uncertain. The simulation of present-day massive star formation by Peters et al. (2010) suggested that photoionization feedback is not very significant, but a general theoretical framework of photoevaporation spanning the whole metallicity range from primordial to quasi-solar metallicities remains lacking.
Feedback by protostellar outflows, radiation pressure and photoevaporation act on the infalling/accreting material. Stellar winds launched from the protostellar surface could in principle also act against the accretion flow, but, as we will discuss below, they are expected to always be confined by the protostellar outflow and thus not have a direct impact on the accretion. However, massloss directly from these stellar winds could potentially become significant, especially for protostars at the highest masses and luminosities. The mass-loss by a stellar wind is certainly important during the later evolution of massive stars. For example, in the case of the R136a1 Wolf-Rayet (WR) star with current mass of 265M ⊙ , a stellar wind mass-loss rate of 5 × 10 −5 M ⊙ yr −1 is inferred, and its initial mass is evaluated to have been as high as 320M ⊙ (Crowther et al. 2010). The theoretical calculation by Vink et al. (2011) has found that the stellar wind mass-loss rate becomes extremely high if the Eddington factor to electron scattering is higher than 0.7. They interpreted this high mass-loss regime as leading to the observational appearance as WR stars and the lower mass-loss regime as O-type stars. However, the protostellar internal luminosity as a function of mass, and thus the Eddington factor, depends on the accretion history. Thus it is possible that in some circumstances the Eddington factor might potentially reach the extreme mass-loss regime even during the protostellar stage.
The feedback and mass-loss processes described above may impact the ability of very massive stars to form and thus reveal themselves in the observed distribution of the IMF, e.g., perhaps creating a break or turnover in the Salpeter (1955) power law that holds from lower masses ∼ 1 M ⊙ to at least ∼ 100 M ⊙ . In other words, feedback and mass-loss may imply there is a maximum stellar mass that can form. Observationally, Figer (2005) have reported the absence of stars with masses > 150M ⊙ in the Arches cluster near the Galactic center, whereas extension of the Salpeter mass function predicts there should be 18 of them. Thus Figer concluded there is an upper stellar mass limit of 150 M ⊙ . Later studies of the Tarantula nebula in the Large Magellanic Cloud (LMC), stars with initial masses of 200-300 M ⊙ were found (Crowther et al. 2010;Bestenlehner et al. 2011). Since LMC has lower metallicity (by about a factor of two) than the Galaxy, it can be speculated that the impact of feedback and/or mass-loss depends on the metallicity, which then affects the upper IMF in different environments.
A trend to a higher maximum stellar mass with decreasing metallicity is potentially supported by the theoretical studies of Hirano et al. (2014); Hosokawa et al. (2016), who found that Pop III stars may reach masses as high as 1000 M ⊙ . Unfortunately there are no direct observational constraints on the masses of Pop III stars. However, the chemical abundance patterns of Galactic metal-poor stars, which may be second generation stars polluted by Pop III supernova ejecta, have been interpreted as indicating that there were such very massive primordial stars (Keller et al. 2014;Aoki et al. 2014). Such conclusions, however, remain very tentative. Overall, a good theoretical understanding of how the stellar IMF depends on metallicity remains lacking.
Although there have been many studies concentrating on each radiative feedback mechanism in massive star formation, there has not yet been a study that has considered all the main processes together, including with the effects of an MHD-launched outflow. In this paper, we aim to carry out such a study with the goal of evaluating the SFE of pre-stellar cores of different masses and in different environments. A full numerical simulation with MHD and radiative effects that resolves the protostellar surface and the outer core scale and follows the full evolutionary growth of the protostar is computationally challenging and beyond current state-of-the-art capabilities. Here we present a semi-analytic model of this process that includes all the expected important physical processes and yet at the same time can be applied to large range of different conditions. This allows us to gain physical insight into the problem and can help guide future numerical simulation experiments. Our modeling builds upon our previous work that developed semi-analytic models for massive star formation (Zhang & Tan 2011;Zhang et al. 2013bZhang et al. , 2014Tanaka et al. 2016), but which did not yet include treatment of radiative feedback or stellar wind mass-loss. Additionally, we apply the same model to primordial star formation at zero metallicity, to gain insight into the metallicity dependence of massive star formation feedback. This paper is organized as follows. In §2 we review the basics of our model, and introduce the updates to include the effects of feedback processes. Next, in §3, we present our results: we show how the accretion rate and SFE are reduced by multiple feedback processes, and also reveal the differences caused by solar and zero metallicities. In §4, we discuss the relative importance of different feedback mechanisms, their impact on shaping the high-mass end of the IMF, and their dependence on metallicity. We conclude in §5.

METHODS
We calculate the accretion history of massive star formation including multiple feedback processes. The framework of our model has been constructed in a series of papers: Zhang & Tan (2011);Zhang et al. (2013bZhang et al. ( , 2014, and Tanaka et al. (2016). In these works, massive protostellar evolution with MHD disk wind feedback is calculated. We then estimate continuum emission from the protostar and disk, which is then followed in radiative transfer calculations, especially to predict infrared, sub-millimeter and cm-radio morphologies and spectral energy distributions. Now we extend this massive protostellar evolution model to include feedback by radiation pressure and photoevaporation, and stellar wind massloss. We note that this follows a similar methodology to that of , who considered the formation of primordial stars under the influence of multiple feedback processes.
We review the basics of our model that were developed in previous works (Zhang & Tan 2011;Zhang et al. 2013bZhang et al. , 2014 in §2.1, introduce the methods for each feedback process in §2.2, and explain how they are combined together in §2.3. Figure 1 shows the schematic view of our model. Although the main target of this study is present-day massive star formation, we also apply our model to primordial star formation for comparison of this Schematic view of massive star formation by core accretion including various feedback processes. The parameters for the initial conditions are core mass Mc, mass surface density of the ambient clump Σ cl , and the ratio of core initial rotational to gravitational energy βc. The model includes momentum feedback from a MHD disk wind and radiation pressure. It also follows massloss resulting from the MHD disk wind, photoevaporation and the stellar wind. different environment and for comparison with the previous results of . Thus in §2.4, we describe the modifications of methods that are used to apply the model to primordial star formation.

Evolution of infall rates, disks and protostars
Our model assumes a pre-stellar core collapses to form one massive star. This model should be a reasonable approximation even for multiple systems in which there is a single dominant protostar. The initial core is assumed to be spherical and close to virial equilibrium by the support of turbulence and/or magnetic fields (McKee & Tan 2003). The parameters to determine core properties are core mass M c , mass surface density of the ambient clump Σ cl , and the core's initial rotational to gravitational energy ratio β c . The core is assumed to be in pressure equilibrium with the ambient clump. If the clump is self-gravitating then this ambient pressure is related to its surface density Σ cl , which sets the pressure at the core surface, thus determining its size. The core radial density profile is assumed to be a power law, i.e., ρ ∝ r −kρ . Observations of dense cores in Infrared Dark Clouds find k ρ ≃ 1.3-1.6 (Butler & Tan 2012;Butler et al. 2014), and we adopt k ρ = 1.5 as a fiducial value, which is the same as the fiducial value used by McKee & Tan (2003) (also Zhang & Tan 2011;Zhang et al. 2013bZhang et al. , 2014. Then, the radius of a core is given as R c = 0.057(M c /60 M ⊙ ) 1/2 (Σ cl /g cm −2 ) −1/2 pc. The core radius is smaller for higher-Σ cl since the core is in pressure equilibrium with the ambient clump. The rotational parameter is fixed as β c = 0.02, i.e., similar to values derived from observations of lower-mass cores (Goodman et al. 1993;Li et al. 2012;Palau et al. 2013). In this study, we investigate the collapse of cores with M c = 10-3000 M ⊙ at Σ cl = 0.1-0.316 g cm −2 , 10-1000 M ⊙ at 1 g cm −2 , and 10-300 M ⊙ at 3.16 g cm −2 .
The inside-out collapse of a core that is a singular polytropic sphere is described by the self-similar solution (McLaughlin & Pudritz 1997;McKee & Tan 2003), which gives the infall rate onto the central protostar-disk system in the limit of no feedback: where M * d (t) = Ṁ * d dt is the collapsed mass, which indicates the mass of the protostar and disk if there was no feedback at all. A higher clump mass surface density leads to a more compact core and thus a shorter free-fall time and higher infall rate. Also, this formula indicates that the infall rate increases with time in the no-feedback case (set by the choice of k ρ = 1.5; a choice of k ρ = 2 would lead to a constant infall rate). To obtain the actual mass accretion rate, we need to calculate the effect of feedback processes, which will be described in §2.2. We note that, in all models in this study, the accretion rates are always smaller than the Eddington rate of ∼ 2 × 10 −2 (r * /10R ⊙ )M ⊙ yr −1 .
Since the initial core is rotating, a disk is assumed to form around the protostar. For simplicity, we only include the effect of rotation inside the sonic point where the infall becomes supersonic and assume that the ratio of the rotational to gravitational energy is constant at this location, β(< r) = β c . Based on the angular momentum conservation from the sonic point, the disk radius is given by (see §2.1 of Zhang et al. 2014). The protostellar disk is expected to be massive and self-gravitating due to high mass supply from the infalling envelope. The angular momentum is transported efficiently by torques in such a massive disk (e.g., Pérez et al. 2016), keeping the mass ratio of disk and protostar approximately constant at f d ≃ 1/3 (e.g, Kratter et al. 2008). We note that these density structures of rotating infall and protostellar disks were developed by Zhang & Tan (2011);Zhang et al. (2013bZhang et al. ( , 2014 and were then used in radiative transfer calculations for synthetic observations. However, in this study, we only focus on the accretion history of forming massive stars, and thus do not need the detailed structure of the envelopes and disks, except for the opening angle of the outflow cavity. Thus we do not expect the results to be very sensitive to the choice of β c as long as β c ≪ 1 so that the outer core structure is quasi spherical. The properties of the protstar, such as luminosity, radius, effective temperature, and their evolution are important to evaluate the strength of feedback. In our study, the protostellar evolution is calculated self-consistently, being adapted to the accretion rate using the model of Hosokawa & Omukai (2009a) and Hosokawa et al. (2010) (which is based on the method developed by Stahler et al. 1980;Palla & Stahler 1991). Since the typical mass-accretion rate in massive star formation is higher than that in low-mass star formation, the rate of entropy carried into the star is also high. This leads to a large protostellar radius of ∼ 100R ⊙ before Kelvin-Helmholtz (KH) contraction starts to be effective (Palla & Stahler 1991;Hosokawa & Omukai 2009a). This swelling causes lower effective temperatures and lower ionizing photon rates than those predicted by the ZAMS model at the same mass. The evolution also depends on the geometry of the accretion flow, i.e., spherical or disk accretion. The accretion geometry is quasispherical when the expected disk radius r d is smaller than the stellar radius r * . In this case, a shock front is produced when this flow hits the stellar surface and a fraction of the released gravitational energy is advected into the stellar interior, which is referred as the "hot" shock boundary. On the other hand, if r d > r * , the material accretes onto the stellar surface through a geometrically thin-disk. In disk accretion, much of the energy radiates away before the material settles onto the star. In the limiting case the entropy carried into star can be assumed to be the same as the gas in the stellar photosphere, which is referred as the "cold" photospheric boundary condition. In our model, the calculation starts from the hot shock boundary and switches to the cold photospheric boundary at r d = r * .
When the accreting material reaches the stellar surface, the accretion energy of L acc = Gm * ṁ * acc /(2r * ) (in the case of disk accretion) is released, whereṁ * acc † is the accretion rate onto the star. Following previous works, we treat this accretion luminosity and the intrinsic internal stellar luminosity as radiating isotropically with a single effective temperature: L * acc = L * + L acc = 4πr 2 * σT 4 * acc , where σ is the Stefan-Boltzmann constant. Following Tanaka et al. (2016), we adopt the stellar atmospheric model "Atlas" (Castelli & Kurucz 2004) to obtain the stellar spectrum L ν, * acc . Then the ionizing photon rate is evaluated as S * acc = ∞ νLy L ν, * acc (hν) −1 dν. Due to line absorption, the ionizing photon rate can be smaller by orders of magnitude than that simply evaluated by integrating over a blackbody spectrum, especially when T * acc 2 × 10 4 K.

Feedback processes
The accretion rate onto the star is smaller than the collapse rate given by equation (1) because of feedback. It is necessary to estimate the impact of feedback to obtain the final mass and the SFE. Here we explain how we evaluate the feedback processes and their effect on the accretion rate.
2.2.1. Outflow driven by momenta of MHD disk wind and radiation pressure † The accretion rate on to the star is described asṁ * in previous works (Zhang & Tan 2011;Zhang et al. 2013bZhang et al. , 2014Tanaka et al. 2016). However,ṁ * acc is adopted in this study since the actual mass growth rate of the star is smaller than this due to the mass loss by stellar wind, i.e.,ṁ * =ṁ * acc −ṁ * w.
The bipolar outflow sweeps up part of the core and thus helps to set the SFE. We calculate the opening angle of the outflow cavity θ esc considering momenta of the MHD disk wind and radiation pressure, i.e., p dw ‡ and p rad . Zhang et al. (2014) included MHD disk wind feedback using the model of Matzner & McKee (2000). In this model, if the outflow momentum is strong enough to accelerate the core material to its escape speed, the outflow extends in that direction. We simply extended this model including the additional term of the radiation pressure: the following equation is satisfied at the polar angle of θ = θ esc (t) where Ω is the solid angle, v esc = 2GM c /R c is the escape velocity from the core, and c g is a correction factor to account for the effects of gravity on the propagation of the shocked shell. Following Zhang et al. (2014), the angular distribution of the core mass is assumed to be isotropic: dM c /dΩ = M c /4π, even though in reality the core would be expected to flatten to some degree by rotation and/or large scale magnetic field support. Based on Appendix of Matzner & McKee (2000), we estimate c g = 2.63 for our core set up. The total MHD disk wind momentum p dw (t) is evaluated by integrating the momentum rate of the wind using a semi-analytic disk wind solution that is modified from the centrifugally driven MHD outflow model of Blandford & Payne (1982): where v K * = Gm * /r * is the Keplerian speed at the stellar radius, φ dw is the factor to measure the disk wind momentum in terms of m * acc v K * (Tan & McKee 2002), f dw is the mass loading rate of wind relative to the accretion rate onto the star (see Zhang et al. 2013bZhang et al. , 2014. We fix the mass loading rate as f dw = 0.1 as a typical value of disk winds (Königl & Pudritz 2000). According to results of our evolution calculation, we find the typical value of φ dw is 0.15-0.3. The angular distribution of the momentum of MHD disk wind is described as (Matzner & McKee 1999;Shu et al. 1995;Ostriker 1997) where θ 0 is a small angle which is estimated to be 0.01, and µ = cos θ (please note that 1 0 P dµ = 1). This angular distribution of P (µ) encapsulates the collimated nature of MHD disk winds. As a result of some trapping by the core, the actual disk wind mass-loss rate is smaller than f dwṁ * acc , which is the limiting value for a fully opened cavity. The fraction of the mass of the wind that can escape from the outflow cavity, f dw,esc , ‡ The subscript "w" was used to represent the MHD disk wind in previous works (Zhang & Tan 2011;Zhang et al. 2013bZhang et al. , 2014Tanaka et al. 2016). However, "dw" is adopted in this study to distinguish with the new component of stellar wind which is described by "*w".
is evaluated based on the fraction of the mass flow in the directions 0 ≤ θ ≤ θ esc . Zhang et al. (2014) where µ esc = cos θ esc . Then, we have the MHD disk wind rate asṁ Note that this value is the mass-loading rate from the disk, however, it is not the total mass-loss by the MHD disk wind from the core. The momentum by the MHD disk wind (and radiation pressure) sweeps up much large amount of gas from the envelope creating the outflow cavity. As we will see in §3, this outflow driven by the disk wind is the most significant feedback. In low-mass star formation, the MHD disk wind is the dominant source of momentum feedback. Radiation pressure becomes significant if the stellar mass reaches ∼ 20M ⊙ . The momentum from radiation pressure p rad is obtained by the integral of the radiation pressure momentum injection rate which is given bẏ where f trap is a trapping factor accounting for the increment of direct radiation pressure force by dust re-emission (Thompson et al. 2005;Murray et al. 2010Murray et al. , 2011. This radiation by dust re-emission should be reduced significantly by the pre-existing MHD outflow cavity (Krumholz et al. 2005;Kuiper et al. 2015Kuiper et al. , 2016 and/or the RT instability (Krumholz et al. 2009;Rosen et al. 2016). Therefore, in the implementation of our model in this paper the effect of dust re-emission is ignored and only direct stellar radiation is considered, i.e., f trap = 1. Therefore, since we are only considering direct stellar radiation, the angular distribution of the radiation pressure momentum is assumed to be isotropic: dp rad /dΩ = p rad /4π. Material in the envelope is swept-up by the momenta of MHD wind and radiation pressure as the opening-up of the outflow cavity. The mass-loss associated with this sweeping process can be evaluated aṡ where the negative sign is chosen to make the mass-loss rate positive sinceμ esc < 0 (see also §2.3). We note that it is not straightforward to clearly distinguish the separate mass-loss contributions here due to MHD disk wind and radiation pressure since the momentum from these two feedback mechanisms combine to open-up the outflow cavity. Below, we will compare the mass-loss from the system by this mechanism with that due to other feedback processes.
We also include effect of shielding by the inner disk. Since this inner disk shielding is efficient to overcome the direct stellar radiation pressure, infall can always continue from the disk shadow region (Tanaka & Nakamoto 2011;Kuiper et al. 2012). Therefore, we limit the maximum opening angle based on the aspect ratio of inner disk, i.e., θ esc,max = tan −1 (H/r). We calculate the inner disk structure with an α-disk model (Shakura & Sunyaev 1973) using the pseudo-viscosity α parameter, which depends on self-gravitational stability (Tanaka & Omukai 2014): α → 0.01 if the disk is stable with respect to self-gravity (the typical situation) implying turbulence is driven by the magneto-rotational instability; α → 1 if the disk is marginally gravitationally unstable implying angular momentum transport is governed by gravitational torques (but this case does not arise in our models for the inner disk region of interest). Note that angularmomentum transport by the disk wind is not explicitly considered here in this calculation (although it was accounted for in the larger scale disk structure calculations of Zhang et al. 2013bZhang et al. , 2014: one expects that its effects would be to change the effective value of α. However, the disk scale height is not very sensitive to α, i.e., H ∝ α 1/10 , and thus we consider that our estimate of the angular size of the shielded region is reasonably well estimated by this method. The aspect ratio is evaluated at the radius of r = 10r * following . Typically, the disk aspect ratio is about 0.1 and thus the maximum opening angle is about 84 • .

Dissipation of envelope and disk by photoevaporation
The ionizing photon rate S * acc increases dramatically after KH contraction starts, and the amount of ionizing photons creates an H II region even during the accretion phase. Ionized gas with high gas pressure can escape from the gravitational binding of the protostellar core, i.e., photoevaporation. We have derived a formula of the mass-loss rate by photoevaporation based on a ray-tracing radiative transfer calculation (Tanaka et al. 2013), which is the updated version of the classic analytic model by Hollenbach et al. (1994). However, Tanaka et al. (2013) did not consider the effect of the dust grains since it was applied to the case of primordial star formation. Here we extended the photoevaporation model including the effect of dust attenuation of ionizing photons.
The photoevaporation mass-loss rateṀ pe is evaluated contribution from both the upper and lower surfaces (Hollenbach et al. 1994;Tanaka et al. 2013), where r g is the gravitational radius inside which the ionized gas is gravitationally bound, r 0 (M * d ) is the collapse radius inside which the enclosing mass was originally equals to M * d (t), n 0 (r) is the base density at the ionization boundary, and c HII is the sound speed of the ionized gas. The gravitational radius in the dustfree case is determined as the escape velocity becomes comparable to c HII , r g,df = Gm * d (1 − Γ e ) /c 2 HII , where Γ e = 2.6 × 10 −5 (L * acc /L ⊙ )(m * /M ⊙ ) −1 is the Eddington factor for electron scattering (Hollenbach et al. 1994;). On the other hand, the gravitational radius in the dusty case can be evaluated as the dust-sublimation radius, r sub = κ * L * acc /4πσκ sub T 4 sub , where T sub is the dust sublimation temperature which we set as 1400 K, and κ * and κ sub are the dust opacity for the stellar radiation and at the dust sublimation temperature, respectively. This is because radiation pressure acting on dust grains assists the ionized gas to become unbounded from the stellar (and disk) gravity especially when photoevaporation occurs actively ( 20M ⊙ ). Therefore, we evaluate the gravitational radius as r g = min(r g,df , r sub ). The outer boundary of integration in equation (11) is chosen as the collapse radius, considering the evaporation not only from the protostellar disk but also from the infalling envelope.
The profile of the base density, n 0 (r), determines the total photoevaporation rate (eq. 11). In the dust-free case, the radiative transfer calculation by Tanaka et al. (2013) showed that the direct stellar radiation dominates at the ionization boundary, and derived an analytic formula of n 0 (r) in the dust-free case as, where α A is the recombination coefficient for all levels (so-called case A) and c pe ≃ 0.4 is the correction factor used to match numerical results. In the dusty region, we extend this formula including the absorption by dust grains as, where τ d is the optical depth caused by dust grains for ionizing photons evaluated from the dust sublimation radius, i.e., and σ a,d is the absorption cross sections of dust grains per H nucleon, which we fix at 10 −21 cm −2 from a typical value of the diffuse interstellar medium (Weingartner & Draine 2001) (however, note that the properties of dust in the upper layers of accretion disks around massive protostars are not well constrained). As we will see in §3, dust attenuation of ionizing photons is important for regulating the total photoevaporation rate. Using this base density profile (eqs. 12 and 13), we obtain the photoevaporation rateṀ pe integrating the equation (11). Please note that, we calculate the temperature of the ionized gas based on the protostellar spectrum and the gas density following Tanaka et al. (2016): it is typically close to 10, 000 K. We note that our model is not a fully self-consistent unification of MHD disk wind and photoevaporation feedback, since neither the magneto-centrifugal acceleration of the photoevaporation flow nor the photoionization mass-loading of the MHD disk wind are considered. The Alfvén speed decreases with distance as the Keplererian speed in the BP wind solution, while the ionized gas sound speed remains constant at ∼ 10 km s −1 . Therefore, the pure-MHD disk wind should dominate in the inner region of r ≪ r g . On the other hand, in the outer region where r ≫ r g , the pure-photoevaporative process is expected to be most important. Additionally, gas in the envelope rotates more slowly than Keplerian, so a magneto-centrifugal wind is not expected to be efficiently launched from this location. Thus, our model is expected to be appropriate at both of the extreme ends of inner and outer radii. Conventionally, those two flows are discussed separately. However, in reality, the mass-loss by thermo-and magneto-hydrodynamical processes occur together, and a unified model is necessary for a more accurate treatment, which we defer to a future paper. For more discussion about "magneto-photoevaporation," see Bai et al. (2016), who studied the MHD disk wind including far-UV/X-ray heating in protoplanetary disks.

Stellar wind mass-loss
The mass-loss via a stellar wind driven by radiative forces on spectral lines is also considered in our model. Vink et al. (2011) studied the stellar wind mass-loss rate up to m * = 300M ⊙ based on Monte Carlo radiative transfer models and dynamically consistent spherical structure. They found two regimes of stellar wind massloss: one is the normal O-type wind regime for Γ e < 0.7; the other is the extreme WR wind regime for Γ e > 0.7. The mass-loss rate dramatically increases with Γ e and they called this upturn at Γ e = 0.7 as the "kink." We adopt a stellar wind mass-loss rate as a function of stellar mass m * and luminosity L * acc based on the fiducial results of Vink et al. (2011): This mass-loss rate is evaluated based on a fixed effective temperature of 50, 000 K. Petrov et al. (2016) have suggested that the mass-loss rate would jump up about one order of magnitude if the effective temperature is lower than 25, 000K. However, our protostellar evolution calculation shows that the effective temperature is always higher than 35, 000K when the Eddington factor is higher than 0.4. Also the variation of mass-loss rate with effective temperature is less than a factor of a few in this high temperature range. Thus, the stellar wind mass-loss rate given by equation (15) is a reasonable approximation for our model, even ignoring the T * acc dependence. Indeed, we will show that the stellar wind mass-loss has only a minor effect compared to other feedback processes.
2.3. Net accretion rate onto stars with feedback We have introduced estimations of the impact of multiple feedback processes. We now evaluate the accretion rate of stars given the effects of these kinds of feedback. The total mass of the envelope at a certain moment is The first term on the right hand side is the sweeping rate by the opening-up of the outflow cavity created by the momenta of the MHD disk wind and radiation pressure (eq. 10). The second term represents the infall rate onto the star-disk system. From mass conservation in the infalling flow, we have µ escṀd * =ṁ * +ṁ * w +ṁ d +ṁ dw +Ṁ pe , (18) whereṁ d is the mass growth rate of the disk (see also Fig. 1). Note that, due to the stellar wind mass-loss, the net accretion rate, or the stellar-mass growth rate, is smaller than the accretion rate onto the star, Following our previous study, the mass ratio of disk and star is assumed to be constant at f d = m d /m * = 1/3 by the self-gravitational-torque regulation (e.g, Kratter et al. 2008), and thus the disk mass growth rate isṁ d = f dṁ * . Using also equation (8), the net mass growth rate of the star iṡ All quantities are time variable except f d and f dw . Eliminating the terms withṀ pe andṁ * w , this equation is identical to that in Zhang et al. (2014). Feedback by radiation pressure does not appear explicitly in equation (20), however, it increases the opening angle θ esc and escape fraction f dw,esc . We continue the protostellar evolution calculation as long asṁ * (t) > 0, i.e., the stellar mass increases, and determine the stellar mass at the moment ofṁ * = 0 as the final mass when it forms m * f . Note that, since the outflow opening angle has a limit set by the disk shielding effect, the outflow from the MHD disk wind and radiation pressure cannot stop mass accretion completely, i.e., µ esc > 0. Therefore, the accretion finishes when (1) mass-loss by photoevaporation and stellar wind is significant, or (2) the entire core collapses, i.e., M * d = M c . We define the instantaneous SFE as the ratio of net accretion rate to infall rate without feedback, i.e., ε * (t) ≡ṁ * (t)/Ṁ * d (t), and otherwise use "SFE" to refer to the ratio of the final stellar mass when the accretion stops to the initial core mass, i.e.,ε * f ≡ m * f /M c . The instantaneous SFE is important since it is in principle observable for individual protostellar systems. For example, Zhang et al. (2016) measured the detailed structure of the HH46/47 molecular outflow using Atacama Large Millimeter/sub-millimeter Array (ALMA), and reported the instantaneous SFE to be 1/4-1/3. However here, we focus mainly on the final SFE rather than the instantaneous SFE to discuss the relation between the CMF and IMF ( §4.2).
Note that the disk mass is so far ignored in the evaluation of the final stellar mass. However, some amount of disk accretion would still be able to continue even after the entire core collapses. The accretion rate is expected to decline as the disk mass to stellar mass ratio drops and self-gravitational torques become less effective. We expect that such a lower accretion rate disk would be more readily dissipated by photoevaporation and/or radiation pressure. However, the fraction of the disk mass that finally accretes onto the star is uncertain because the actual angular momentum transport processes are uncertain at this stage. Therefore, for simplicity, we ignore the disk mass in the SFE evaluation, and note that the actual SFE may be underestimated by up to a factor of 1 + f d → 4/3.

Primordial star formation
Although the main purpose of this paper is the study of feedback in massive star formation in present-day universe, we also apply the same feedback model to primordial star formation in the early universe for comparison and demonstration of our model. Here we describe modifications of the present-day massive star formation model for its application to primordial star formation. These modifications follow the methods of Tan & McKee (2004), Tan & Blackman (2004) and  for primordial star formation. Tan & McKee (2004) predicted the evolution of the mass infall rate, accretion disk structure, and protostellar evolution associated with primordial star formation. We use the results of Tan & McKee (2004) for the infall rate and disk evolution. The infall rate excluding effects of feedback is given bẏ Here K ′ is the entropy parameter of the polytropic equation of state of the cloud; larger values of K ′ correspond to denser gas cores. In this study, the entropy parameter is fixed at the fiducial value of K ′ = 1. The above infalling rate replaces that of equation (1) that is used to model present-day star formation. It is interesting that the typical infall rates in primordial star formation and present-day massive star formation are coincidentally of the same order. The high accretion rate in the primordial case is induced by the high gas temperature in the core due to inefficient cooling at zero metallicity, while for the present-day case the high turbulence and strong magnetic fields enhance the effective pressure of cores leading to their high accretion rates.
Based on the conservation of angular momentum from the sonic point to the outer radius of disk, the disk radius around Pop III protostars is evaluated as where f Kep is a angular momentum parameter of infalling gas, which is fixed at the fiducial value of 0.5 (Abel et al. 2002;Tan & McKee 2004).
As in the case of present-day massive star formation, protostellar evolution is calculated self-consistently adapted to the accretion rate using the code developed by Hosokawa & Omukai (2009a) and Hosokawa et al. (2010), except the opacity is modified for zero metallicity. As a result of similar accretion rates of 10 −3 M ⊙ yr −1 , the evolution of primordial protostars is expected to resemble that of present-day massive protostars (Omukai & Palla 2001. The main difference of protostellar evolution is that stellar radius in the main-sequence phase is smaller at zero metallicity than that at solar metallicity. This is because, due to the lack of C and N, the KH contraction continues until the temperature becomes high enough for small amounts of carbon to be produced by He burning, which then enables the operation of the CNO cycle. For more details of comparison of protostellar evolution in zero and solar-metallicities, see §3.4 of Hosokawa & Omukai (2009a).
For our modeling, we also update the stellar spectrum appropriate for the case of zero metallicity (Schaerer 2002) in order to evaluate the ionizing photon rate. This leads to higher ionizing photon luminosities for a given temperature due to a lack of metal line absorption.
The feedback model also needs some modifications to apply to primordial star formation. In the outflow feedback, we neglect the momentum by radiation pressure, i.e., p rad = 0 in equation (3), since there are no dust grains. Following Tan & Blackman (2004), we do include the MHD disk wind momentum, since the MHD disk wind could be driven by the disk-dynamo generated magnetic field. Due to the different density profile inside the core, the correction factor accounting for effects of gravity on shock propagation is c g ≃ 4.6. The escape velocity is evaluated by v esc,c = 3.22K ′5/7 (M * d /1000M ⊙ ) −1/7 km s −1 (Tan & Blackman 2004). In the photoevaporation feedback calculation, dust attenuation is set to zero, τ d = 0. This means photoevaporation is more efficient in primordial star formation. Finally, we neglect the the stellar wind mass-loss (i.e.,ṁ * w = 0), which is mainly driven by the metal lines.

RESULTS
In this section we first present the general evolution of massive formation by core accretion. Then, we examine details of individual feedback process. Next, we show the results of primordial star formation to demonstrate the effect of metallicity. Finally, we show the obtained SFE for various initial cores.
3.1. Accretion history and protostellar evolution Figure 2 shows results of our modeling of present-day protostars forming from cores with initial mass M c = 1000 M ⊙ . Three different clump mass surface densities are considered.
First, consider the case with Σ cl = 1 g cm −2 . As the infall rate increases (Fig. 2a), the net accretion rate also increases to 2 × 10 −3 M ⊙ yr −1 until the stellar mass reaches 170 M ⊙ (at t = 1.35 × 10 5 yr). Then the accretion rate drops, finally stopping at 285M ⊙ (t = 2.64×10 5 yr). The decline of accretion is mainly caused by the opening-up of the outflow cavity (Fig. 2b) given the increasing momentum of the disk wind and from radiation pressure (Fig. 2c), rather than by mass-loss by photoevaporation or via the stellar wind (Fig. 2d). It is clearly seen that the outflow sweeping rate is orders of magnitude higher than other mass-loss rates in Figure 2d, which indicates the MHD outflow, assisted by radiation pressure, is the most dominant feedback process. The MHD disk wind always dominates total momentum, however, the radiation pressure also can give significant assistance to open-up the outflow cavity (see §3.2.1 for more details). Mass-loss by photoevaporation quickly rises to ∼ 10 −5 M ⊙ yr −1 when the stellar mass is about 15M ⊙ , and then increases to ∼ 10 −4 M ⊙ gradually after that. However, the photoevaporation mass-loss rate reaches a maximum of about 10 −4 M ⊙ yr −1 , which is never enough to shut down accretion. Mass-loss by the stellar wind is even smaller than that from photoevaporation (thus, the accretion Figure 2. Evolution of protostellar and feedback properties as functions of the protostellar mass m * for initial core mass of Mc = 1000M ⊙ and Σ cl = 0.1, 0.316, and 1 g cm −2 (orange, green and blue lines, respectively). The panels show: (a) net accretion rate,ṁ * ; (b) outflow cavity opening angle, θesc; (c) momenta of MHD disk wind, p dw (solid), and radiation pressure, p rad (dashed); (d) outflow sweeping rate, Mswp (solid), photoevaporation mass-loss rate,Ṁpe (dashed), and stellar wind mass-loss rate,ṁsw (dotted); (e) stellar radius, r * ; (f) effective temperature, T * acc; (g) total luminosity, L * acc; (h) ionizing photon rate, S * acc. In the right-hand panels, the ZAMS properties are also plotted by black dash (Schaller et al. 1992) and dotted (Ekström et al. 2012) lines. rate onto the star is almost equal to the net accretion rate, i.e.,ṁ * acc =ṁ * ). Therefore, mass accretion only finishes when the entire initial core collapses. The SFE in this case isε * f = 285M ⊙ /1000M ⊙ = 0.285. Please remind that our evaluation of the final mass ignores the disk mass, and this SFE is the minimum estimation with the maximum error of ∆ε * f = 0.095 ( §2.3).
The evolution of protostellar properties are shown in the right panels of Figure 2. At around m * = 8M ⊙ , the stellar radius (Fig. 2e) suddenly increases by a factor of three which is due to the redistribution of entropy in the protostar (Hosokawa & Omukai 2009a;Hosokawa et al. 2010). The protostar reaches the local maximum radius of 30 R ⊙ at m * = 11 M ⊙ . Until this time the total luminosity ( Fig. 2f) is dominated by accretion luminosity, i.e., the total lumininosity is significantly larger than the ZAMS luminosity. The effective temperature of the protostar (Fig. 2g) is relatively low due to this large stel-lar radius. Therefore the ionizing photon rate (Fig. 2h), which is very sensitive to the effective temperature, is lower than that of ZAMS model by about five orders of magnitude.
At later times and greater masses the protostar undergoes KH contraction and approaches the main sequence structure. The effective temperature increases to 45, 000 K, and thus the ionizing photon rate also dramatically rises leading to the start of significant photoevaporation. The star evolves almost along the ZAMS line at the mass range of 30-100 M ⊙ . Then, the stellar radius again becomes slightly larger than the ZAMS model. This deviation is related to the metal opacity near the stellar surface and the high accretion rate of 10 −3 M ⊙ yr −1 (see Ishii et al. 1999;Gräfener et al. 2012). We note that this small deviation from the ZAMS has little effect on the significance of feedback and on the final stellar mass.
Next, we consider the cases of 1000 M ⊙ cores in lower mass surface density environments, i.e., protostars with lower accretion rates (see green and orange lines in Fig.  2). We find that the impact of radiation feedback becomes more significant. The momentum input due radiation pressure at a given stellar mass is higher in the lower Σ cl cases, while the MHD disk wind momentum is almost identical for all cases. Due to this higher radiation pressure momentum, the outflow opens up at lower masses in these lower Σ cl cases. Photoevaporation also has a larger impact since the accretion rate is lower. Especially in the case of Σ cl = 0.1 g cm −2 , photoevaporation shuts down mass accretion before the entire core collapses, i.e., M * d < M c . In this way, the relative importance of radiative feedback becomes higher and results in lower SFE, i.e.,ε * f = 0.29, 0.18 and 0.087 for Σ cl = 1, 0.316 and 0.1 g cm −2 , respectively.
The phases of protostellar evolution are shifted to lower stellar masses for the lower Σ cl cases because of their lower accretion rates, which thus mean it takes a longer time for the protostars to reach a given mass. However, after the KH contraction phase, the protostellar evolution following the ZAMS structure very closely in all cases, at least up to ∼ 100 M ⊙ .

Individual feedback processes
Here we describe results concerning each feedback process included in the modeling.

MHD disk wind and radiation pressure driven outflows
The outflow in our model is driven by the momenta of the MHD disk wind and by radiation pressure. As we have seen, the total momentum is dominated by the MHD disk wind. However, radiation pressure also plays a role in helping to open up the outflow cavities, since it acts more isotropically than the collimated disk wind. Figure 3a shows the evolution of net accretion rates as the protostars grow in mass for models with M c = 1000M ⊙ and Σ cl = 0.1, 0.316, and 1 g cm −2 . For comparison, the results with only MHD disk wind feedback and those with no feedback are also shown. For the case of Σ cl = 1 g cm −2 with no feedback, the SFE is unity and thus the stellar mass can reach the core mass of 1000 M ⊙ . Including MHD disk wind feedback, the accretion rate is lowered due to the outflow and accretion stops Figure 3. (a) Accretion rate history as a function of protostellar mass for full feedback models with Mc = 1000M ⊙ and Σ cl = 0.1, 0.316 and 1 g cm −2 (orange, green and blue solid lines, respectively). The results of these cases with no feedback (dotted) and those with only MHD disk wind feedback (dashed) are also shown. (b) Evolution of the mass-loss rate from the envelope due to outflow sweeping, i.e., opening-up of the outflow cavity, for the same full feedback and only MHD disk wind feedback models. (c) Evolution of the outflow opening angle θesc for the same full feedback and only MHD disk wind feedback models. (d) Evolution of the ratio of momenta from radiation pressure and from the MHD disk wind at θesc = 84 • , i.e., (dprp/ dΩ) / (dp dw / dΩ) | θesc=84 • , for the same full feedback models shown above. at m * = 470 M ⊙ . In the case of all feedback, accretion drops significantly after m * ≃ 200 M ⊙ and is finished by m * = 285M ⊙ . The plateau of the accretion rate around 250 M ⊙ is caused by disk shielding. The SFE including all feedback (ε * f = 0.285) is reduced compared to that resulting with only MHD disk wind feedback (ε * f = 0.470). As described above, the mass loss by photoevaporation and stellar wind is not very significant, and the decline of SFE is mainly due to the radiation pressure and its effect on opening up the outflow cavity. Figure 3b shows the evolution of the mass-loss rates from the envelope due to outflow sweeping, i.e., opening-up of the outflow cavity. When m * 100M ⊙ , the full feedback models are similar to the only MHD disk wind models. However, in the higher mass regime, the sweeping rate in the full feedback model becomes higher than that in the only MHD disk wind model and then quickly drops off. This phenomenon shows that the outflow cavity opening rate is enhanced by radiation pressure and more quickly reaches the limit set by disk shielding. This can be also seen in Figure 3c: the cavity opening in the full feedback models accelerates when θ esc reaches about 30 • .
The MHD disk wind dominates the total momentum, however, its angular distribution is highly collimated near the outflow axis (eq. 6). On the other hand, the radiation pressure is modeled as having an isotropic momentum distribution. Therefore, while the MHD disk wind initially creates the outflow cavity, radiation pressure has a significant impact in making it wider. Figure  3d shows the ratio of momenta due to radiation pressure and the MHD disk wind at an angle of θ = 84 • (close to the maximum angle allowed given disk shielding). In the case of Σ cl = 1 g cm −2 , the contribution of radiation pressure becomes similar to that of the disk wind at around 200 M ⊙ , which causes the accretion rate to start falling. This reduction of accretion rate also leads to the decline of the momentum rate by the disk wind since it is accretion powered (eq. 4), and the relative importance of radiation pressure increases even more. In the cases of lower Σ cl , the radiation pressure has a larger impact, starting to dominate at lower stellar masses. This is because the lower Σ cl leads to lower accretion rates and lower momentum injection rates from the disk wind, while the radiation pressure does not strongly depends on the accretion rate. In this way, radiation pressure has an important impact on the decline of SFE even though the MHD disk wind dominates the total outflow momentum.

Mass loss by photoevaporation
In the models shown in §3.1, photoevaporation is not a significant feedback in setting the SFE, even when m * > 100M ⊙ since the photoevaporation mass-loss rate is only ∼ 10 −4 M ⊙ yr −1 at its maximum, while the accretion rate is 10 −3 M ⊙ yr −1 . Here we show the importance of dust attenuation of ionizing photons in the reduction of photoevaporation feedback.
The calculation of the photoevaporation mass-loss rate in our model includes the effect of dust attenuation on the propagation of ionizing photons using the optical depth τ d (r) ( §2.2.2). To measure the effect of dust attenuation, we introduce a characteristic optical depth of the system asτ d = σ a,d n 0 (r sub )r sub . Note that, as we will see in §4.1, this characteristic optical depth is not exactly the same as the total optical depth of the flow τ d (∞), however it gives a good indication of the optically thin/thick boundary and well represents the effect of dust attenuation. In Figure 4 we show the characteristic optical depthτ d and the photoevaporation mass-loss rate as functions of ionizing photon production rate for models with M c = 1000 M ⊙ . It can be seen that the characteristic optical depth increases with S * acc andṀ pe , and reaches the optically thick regime wheṅ M pe ≃ 2 × 10 −6 M ⊙ yr −1 , which is much smaller than the typical infall rate. Even in the optically thick regime, the photoevaporation mass-loss rate still increases with S * acc , however it does not reach ∼ 10 −3 M ⊙ yr −1 , which is needed to be a significant feedback effect.
In the bottom panel of Figure 4, to illustrate the importance of dust attenuation on the photoevaporation massloss rate, we also plot hypothetical ratesṀ pe,τ d =0 , which are evaluated neglecting dust attenuation, i.e., assuming τ d (r) = 0. In the optically thin regime withτ d < 1, the rates with and without dust attenuation are similar. However, in the optically thick regime withτ d > 1, the actual photoevaporation mass-loss rate becomes much smaller thanṀ pe,τ d =0 . As we discuss in §4.1, we find that the reduction of photoevaporation mass-loss rate by dust attenuation can be approximately described aṡ M pe /Ṁ pe,τ d =0 ≃ 1/τ d in the case ofτ d ≫ 1. The reduction factor becomes more than one order of magnitude at high ionizing photon rates of 10 49 s −1 , when the hypothetical mass-loss rate without dust attenuation would exceed a few×10 −4 M ⊙ yr −1 . Thus, dust attenuation is very important in limiting the impact of photoevaporation feedback in present-day massive star formation.

Mass loss by stellar winds
Mass-loss via stellar winds is a minor effect compared with other processes. Figure 5 shows the stellar wind mass-loss rate and the Eddington factor Γ e as functions of protostellar mass from results of models with M c = 1000 M ⊙ . The stellar wind mass-loss rate is about 10 −5 M ⊙ yr −1 even at m * = 250 M ⊙ , which is much smaller than typical values of accretion rate and photoevaporation rate. The obtained Eddington factor with protostellar evolution calculation is slightly higher than that from the ZAMS model. However, it is not high enough to reach the critical "kink" value of Γ e = 0.7, above which the mass-loss rate dramatically increases . Therefore, we conclude that massloss by stellar winds is not a significant feedback effect for setting the SFE.
Note that the momentum input from the stellar wind has been ignored in our estimation of the outflow opening angle. The stellar wind momentum rate can be evaluated asṗ sw ≃ṁ sw v * esc , where v * esc = 2Gm * (1 − Γ e )/r * is the escape velocity from the stellar surface. We find that the stellar wind momentum is at most 10% of the radiation pressure component, and no more than about 1% of the MHD disk wind component. We thus expect that the stellar wind would anyway be confined and collimated by the MHD disk wind, so would not significantly impact the opening of the outflow cavity.

Primordial star formation
We apply our model also for the case of primordial star formation, which gives a limiting case for the effects of metallicity on massive star formation feedback ( §2.4). Figure 6 shows the comparison of primordial star formation (K ′ = 1) and present-day massive star formation (M c = 1000 M ⊙ and Σ cl = 1 g cm −2 ). Due to differences of core density structure, the evolution of accretion rate is different: the accretion rate decreases with time in primordial star formation, while it increases in present-day massive star formation (see Fig 6a and eqs. 1 and 21). Therefore, differences between the cases are not only due to metallicity. However, since the accretion rates are in fact quite similar at ∼ 10 −3 M ⊙ yr −1 when the protostar has ∼ 30 M ⊙ an approximate comparison to see the effects of metallicity is possible.
We find that accretion stops at 44.4 M ⊙ in the primordial case with K ′ = 1, which is a much lower mass than we found for the present-day case with M c = 1000 M ⊙ and Σ cl = 1 g cm −2 . The main reason for this is the high photoevaporation mass-loss rate from the primordial protostar (Fig. 6b). As described in §3.2.2, dust attenuation of ionizing photon strongly regulates the photoevaporation mass-loss rate to be 10 −4 M ⊙ yr −1 at solar metallicity. However, at zero metallicity the photoevaporation mass-loss rate can reach ∼ 10 −3 M ⊙ yr −1 .
One may suppose that this difference of photoevaporation mass-loss rate in the zero and solar metallicity cases is also related to differences in protostellar evolution and the stellar spectra. As described in §2.4, the primordial protostar contracts to a smaller ZAMS structure than the solar metallicity protostar because of its Figure 5. The stellar wind mass-loss rate (top) and the Eddington factor, Γe (bottom), as functions of protostellar mass, m * , for models with Mc = 1000M ⊙ and Σ cl = 0.1, 0.316 and 1g cm −2 (orange, green and blue lines, respectively). In the bottom panel, the Eddington factor evaluated for the ZAMS model (Ekström et al. 2012) is shown by the dotted line, and the critical value of 0.7 is indicated by the horizontal dashed line. initial lack of heavy elements. The evolution of stellar radii is shown in Figure 6c, showing that the primordial protostar contracts to a smaller size after m * ∼ 30 M ⊙ . This smaller radius causes higher effective temperatures and thus higher ionizing photon production rates (also aided by the lack of metal line absorption in the stellar atmosphere). However, these differences do not dramatically increase the ionizing photon rate at zero metallicity (Fig. 6d). When m * 10 M ⊙ , the ionizing photon rate is higher at zero metallicity than that at solar metallicity, because of higher accretion rates and luminosities in this earlier phase (Fig. 6a). At the higher mass range of m * 20 M ⊙ , the difference of ionizing photon rates by metallicity becomes modest (indeed, the solar metallicity case has even slightly higher ionizing photon production rates at 15-30 M ⊙ due to its smaller radius during this phase). The ionizing photon production rate difference is less than a factor of three at m * ∼ 40 M ⊙ , which is not enough to explain the one order of magnitude difference in photoevaporation mass-loss rate (see bottom panel of Fig. 4 and eq. 13). Therefore, we conclude that dust attenuation of ionizing photons is the most significant effect controlling the metallicity dependence of photoevaporation mass-loss rates.
Our model finds a smaller final stellar mass of 44 M ⊙ than the study of , who found 140M ⊙ for the K ′ = 1 case. The main difference is that we have included MHD disk wind feedback. As we have seen in the case of present-day massive star formation, Figure 6. Comparison of primordial star formation (K ′ = 1, orange lines) and present-day massive star formation (Mc = 1000 M ⊙ and Σ cl = 1 g cm −2 , blue lines): (a) net accretion rates,ṁ * ; (b) photoevaporation mass-loss rates,Ṁpe; (c) stellar radii, r * ; (d) ionizing photon production rates, S * acc. In panels c and d, the ZAMS models are also shown for reference: black dashed lines for zero metallicity (Schaerer 2002) and black dotted lines for solar metallicity (Ekström et al. 2012).
the MHD disk wind is the dominant feedback in the lowmass regime. This reduction of accretion rate at lower masses results in an earlier start of KH contraction and thus of effective photoevaporation feedback. Another difference also results from our updated protostellar evolution calculation and photoevaporation model compared with that of . For protostellar evolution, we use a detailed protostellar structure calculation code, which tends to predict a smaller protostellar size: e.g., at m * = 30 M ⊙ the model of  had r * ∼ 20 R ⊙ (see Fig.2 of  while we now estimate r * ≃ 10 R ⊙ . For the photoevaporation calculation, we adopt the model by Tanaka et al. (2013), who showed the importance of photoevaporation from the outer region based on an accurate radiative transfer calculation, while the analytic Hollenbach et al. (1994) model suggested the mass-loss rate is dominated by the region close to the inner gravitational radius r g . Including mass-loss from the outer region, including also the collapsing envelope that is exposed by the outflow cavity, the photoevaporation rate can be higher by a factor of (R c /r g ) 0.5 ∼ 10 than the Hollenbach et al. model (Tanaka et al. 2013). These differences lead to an enhancement of feedback compared to the study of  and result in a smaller final mass than found in their model.

Star formation efficiency
Now we return to the case of present-day massive star formation and explore how the SFE,ε * f , depends on core mass, M c , and ambient clump mass surface density, Σ cl . The left panel of Figure 7 shows SFEs for Σ cl = 0.1-3.16 g cm −2 as functions of final protostellar mass, m * f , i.e., when accretion stops. The SFE with only MHD disk wind feedback has a weak dependance on m * f , with values of ∼ 0.3-0.5, similar to the results of Matzner & McKee (2000) and Zhang et al. (2014). Note, that for the highest Σ case we have not run the MHD disk wind only cases since their very high accretion rates lead to protostellar structures that are difficult to model numerically with our adopted protostellar evolution code. On the other hand, the SFE in the models with radiation feedback decreases quite strongly with the final stellar mass for all Σ cl cases. The deviation from the MHD disk wind only case is small if the final stellar mass is less than 10M ⊙ . The SFE becomes much smaller as m * f increases, since radiative feedback grows strongly with stellar mass. The results of Σ cl = 0.1 g cm −2 shows the strongest impact of radiative feedback. In this case, the SFE is only 0.1 or less when forming 100 M ⊙ stars. On the other hand, for higher Σ cl , as we have seen in §3.1, the impact of radiative feedback is smaller due to the higher accretion rates. The dominant feedback mechanism for setting SFEs is the MHD disk wind for Σ cl 0.3 g cm −2 , even in the formation of very massive stars.
The right panel of Figure 7 shows the SFEs as functions of initial core radii, R c . We see that more compact cores result in higher SFE. Interestingly, all of our models with full feedback can be fitted by a single power law of within an error of 35%. This simple fitting formula is convenient analytic result that can be applied as a subgrid model to large scale simulations of star formation that resolve formation of massive pre-stellar cores (note, this result applies to cores from 10 M ⊙ to ∼ 10 3 M ⊙ ).

DISCUSSION
First, we summarize the relative importance of different feedback mechanisms in §4.1. Then we discuss the potential impact of radiation feedback on shaping the high-mass end of the IMF in §4.2. Next we consider the metallicity dependence of massive star formation in §4.3. Finally, we note the caveats and limitations of our modeling in §4.4.

Relative importance of feedback processes
We have studied multiple feedback mechanisms, i.e., MHD disk wind, radiation pressure, photoevaporation, and stellar wind, during star formation via core accretion. We find that for present-day massive star formation at solar metallicity the MHD disk wind plays a dominant role not only in low-mass star formation but also in massive star formation.
In simple spherical core collapse radiation pressure acting on dusty infall stops formation of massive star formation for m * 20M ⊙ . However, in non-spherical disk accretion, the optically thick inner region shields outer equatorial zone accretion. Additionally, the MHD disk wind outflow cavity effectively reduces the effects of radiation pressure by dust re-emission, i.e., f trap ≃ 1. Using equations (4) and (9), we can compare the momentum injection rates from the MHD disk wind and from radiation pressure: Here we use the luminosity and radius of the ZAMS model (Schaller et al. 1992) and adopt φ dw ≃ 0.2 from our results (see also Zhang et al. 2014). As seen from this evaluation, the MHD disk wind momentum injection rate is much higher than that from radiation pressure even for m * = 100 M ⊙ . However, as described in §3.2.1, the MHD disk wind is collimated while the stellar radiation acts isotropically. Considering the angular distribution of momenta (eq. 6), we obtain the following relation, It can be seen that the component of radiation pressure is not negligible at large angles θ away from the outflow axis. Also, these equations indicate that the contribution of radiation pressure is higher at lower accretion rates, i.e., lower Σ cl cases. The accretion rate also becomes smaller when the outflow cavity opens up, which enhances the importance of p rad . In this way, the MHD disk wind supplies a large measure of momentum to create the outflow, and the radiation pressure assists to open up the cavity and help set the SFE. As shown in §3.2.2, the photoevaporation mass-loss rate is regulated by dust attenuation of ionizing photons and is a relatively minor feedback process, unlike in the case of primordial star formation. Of course, dust attenuation only occurs in the region where dust survives, i.e., r > r sub . For this region, assuming a constant recombination rate α A , we can derive a simple differential equation from equations (13) and (14): where n sub = n 0 (r sub ) is the base density of the photoevaporation flow at the dust sublimation front, and x ≡ r/r d is a dimensionless radius. This equation has an analytic solution of The characteristic optical depth of the systemτ d , which also appears in §3.2.2, is evaluated aŝ from equation (12) and assuming an ionized gas temperature of 10 4 K. This solution is consistent with the dust-free case in the limit ofτ d → 0. Considering limits of a far distance of x ≫ 1, we see basic features of the effects of dust attenuation, i.e., τ d → 2 ln(1 +τ d ) and n 0 → n sub x −1.5 /(1 +τ d ). The photoevaporation flow reaches optically thick conditions, i.e., τ d = 1, when the characteristic optical depth reachesτ d ≃ 0.7. However, in the optically thick limit ofτ d ≫ 1, the total optical depth converges to 2 lnτ d , which is smaller thanτ d . Thus, the suppression of base density n 0 is not as strong as the exponent of e −τ d and involves only a factor ofτ d . Due to this suppression of photoionization, the total evaporation rate is also regulated tȯ M pe /Ṁ pe,τ d =0 ≃ 1/τ d since the evaporation rate is proportional to n 0 (eq. 11). In this manner, dust attenuation of ionizing photons regulates the mass-loss rate by photoevaporation whenτ d 1.
The stellar wind is found to be the weakest feedback process in our model. Please note that we have not explicitly considered the momentum injection from the wind, since it is always sub-dominant compared to the MHD disk wind and radiation pressure and would be confined along a narrow region of the outflow axis.
The Eddington factor Γ e evaluated by our protostellar calculation is typically higher than that of the ZAMS case, however it is still smaller than the critical value of 0.7 to initiate the extreme wind mass-loss regime (Fig.  5). Moreover, even assuming the maximum Eddington factor of Γ e = 1 in equation (15), the stellar wind massloss rate is lower than 10 −4 M ⊙ yr −1 at 100M ⊙ . Therefore, we conclude that, the stellar wind is not important during the protostellar accretion phase. Note, however, that during evolution after the mass accretion phase over timescales of ∼ Myr, the stellar wind has important effect leading to significant mass-loss.
To conclude, in massive star formation by core accretion, we find that the MHD disk wind is most important feedback mechanism, radiation pressure assists the opening-up of the outflow cavity to wide angles, photoevaporation is regulated by the dust attenuation and is thus of minor importance, and stellar wind mass-loss has a very minor effect during the accretion phase. In the sense that the bipolar MHD-driven outflow is the most significant feedback, massive star formation resembles low-mass star formation, however SFEs can be significantly reduced by the action of radiative feedback.

4.2.
Implications for the high-mass end of the IMF As we have seen, radiative feedback can significantly reduce SFE for the formation of very massive stars. Considering the stellar IMF to be the result of a multiplicative combination of the CMF and SFE, we can expect that the effects of radiative feedback may be seen in the high-mass end of the IMF. While the MHD disk wind only feedback sets a SFE, which depends only weakly on the initial core mass, the full model including radiative feedback yields smaller SFEs for higher core masses, M c (Fig. 7). Potentially this could induce a steepening of the IMF at the highest masses and if this is steep enough it may appear as an apparent truncation. Using the obtained SFE, we can relate the IMF to the CMF. We introduce the exponent of ε ′ ≡ d lnε * f /d ln M c . If the CMF is a power-law distribution of dN /d ln M c ∝ M −αc c and the exponent of ε ′ is constant, then the IMF would be (Nakano et al. 1995;Matzner & McKee 2000). The SFEs for massive cores with M c = 10-3000 M ⊙ and Σ cl = 0.1-3.16 g cm −2 obtained by our model are well within an error of 15%. Then, the power law exponent of the IMF at > 10M ⊙ is estimated as −1.13α c in clumps with Σ cl = 1 g cm −2 . Assuming an initial CMF slope of α c = 2.35, i.e., the same as the Salpeter IMF (Salpeter 1955), the total number of stars with 10-100 M ⊙ is smaller than than that simply expected from the CMF slope by about 11(Σ cl /g cm −2 ) −0.35 % due to MHD disk wind and radiative feedback. However, this reduction of massive stars is too small (the reduction factor is about 55% for the mass range 100-300 M ⊙ ) to explain the cut-off at 150 M ⊙ reported for the Arches cluster (Figer 2005). Thus we conclude that the highmass end of IMF, especially its potential truncation at masses ∼ 150-300 M ⊙ , is mainly determined by the prestellar core mass function rather than by feedback.
4.3. Radiation feedback in massive star formation at different metallicities In this paper, we mainly study the formation of massive stars at solar metallicity. However, we have also applied our model to the primordial star formation case ( §2.4 & 3.3) in order to compare to previous studies and to obtain a basic insight into the effect of metallicity. Our results show that the impact of radiative feedback depends strongly on such metallicity changes. Since massive stars are thought to have been important thoughout cosmic history as metallicities have evolved from primoridal, near-zero limits to ∼solar values and beyond, here we give some general discussion about the dependence of radiative feedback, especially radiation pressure and photoevaporation mass-loss (stellar wind feedback is weak at solar metallicity and would be even weaker at lower metallicities). However, we defer a detailed quantitative investigation of massive star formation at intermediate metallicities of 0 < Z < Z ⊙ to a future paper.
Radiation pressure is the strongest radiative feedback mechanism at solar metallicities. However, since it acts on dust grains in the infalling envelope, it must depend on metallicity. We have ignored dust re-emission since it is assumed to effectively escape from the outflow cavities. Then, the momentum injection by radiation pressure can be evaluated assuming f trap = 1 in equation (9), as long as the envelope is optically thick for direct stellar radiation. In other words, the effect of radiation pressure becomes weaker if the metallicity is low enough to keep the envelope transparent for direct stellar radiation. This transparency depends on the stellar spectrum and the grain components, but the typical opacity for direct stellar radiation is approximately evaluated as κ * ∼ 100(Z/Z ⊙ ) cm 2 g −1 for massive stars assuming the opacity is simply proportional to metallicity. Then, the trapping factor is approximately given by Assuming a typical massive core always forms in a clump with a mass surface density of 1 g cm −2 , the effect of radiation pressure would become weaker for metallicities of Z 10 −2 Z ⊙ . Photoevaporation is strongly suppressed at solar metallicities because of the dust attenuation of ionizing photons. As described §4.1, the photoevaporation rate with dust attenuation is about 1/τ d of that of the dustfree case ifτ d ≫ 1. Then, if the dust opacity for ionizing photons is simply proportional to the metallicity, we obtain the following relation oḟ Here we assumeτ d ∼ 100 at solar metallicity as a typical value for the high ionizing photon production rate of 10 50 s −1 at which photoevaporation mass-loss rate could be important (eq. 29). Therefore, the photoevaporation mass-loss rate could be as high as that at zero metallicity at metallicities of Z 10 −2 Z ⊙ . Considering both radiation pressure and photoevaporation, the critical metallicity for their transitions coincide at ∼ 10 −2 Z ⊙ . Dust absorption is efficient at higher metallicities than this critical value, which means that radiation pressure acts effectively. On the other hand, photoevaporation is suppressed at these higher metallicities. In the lower metallicity regime dust absorption is weak, which lessens the impact of radiation pressure, while photoevapration is more effective. These considerations suggest that the total effects of radiative feedback may be strongest at ∼ 10 −2 Z ⊙ .
However, note that we are not suggesting that massive star formation is necessarily rarer at metallicities of ∼ 10 −2 Z ⊙ . Only that SFE could be lower. The core mass function will also play an important role, along with the typical clump mass surface density. It is difficult to predict how these will vary with metallicity, especially since they may also be more strongly influenced by the degree of magnetization of the gas. Other processes, such as disk fragmentation (e.g. Tanaka & Omukai 2014), may also play a role.

Caveats
Even though our model predictions, including those from previous papers in this series, have some agreements with observations (Zhang et al. 2013a;Tanaka et al. 2016), this is a semi-analytic model that is still highly simplified and idealized. Ultimately, the predictions of the model need to be tested by full radiation-MHD numerical simulations, especially to study the interaction of the outflow with the infall envelope and establishment of the outflow cavity boundary. Below we discuss some additional caveats of our modeling.
We have considered only single star formation. The massive cores are expected to be supported mainly by non-thermal pressures, i.e., turbulence and magnetic fields, which keeps them from fragmenting to the thermal Jeans mass of ∼ 1M ⊙ (e.g., McKee & Tan 2003). Also the catastrophic fragmentation during collapse is expected to be suppressed by radiative heating by the high accretion luminosity and the efficient angular momentum transportation by magnetic breaking (Krumholz et al. 2007;Commerçon et al. 2011). However, a small amount of fragmentation may still occur, as seen simulations such as Krumholz et al. (2009). Indeed, it is observationally known that more than 70% of massive stars have close companions which eventually exchange their masses (Sana et al. 2012). We expect that our model is still quantitatively appropriate since the feedback is dominated by a single object as long as the total stellar mass is dominated by a primary star in the binary/multiple system.
On the other hand, our feedback model would need significantly modification if the system contains similar mass stars. Qualitatively, we expect that radiative feedback in the case of similar mass binaries would be weaker than that in the case of formation of a single massive star. This is because the stellar luminosity increases nonlinearly with the mass. If some amount of material is divided into two objects, the total luminosity is smaller than that of a single star with the same total mass. In contrast, the momentum rate from MHD disk winds is roughly proportional to the total accretion rate. Therefore, we expect that the conclusion that the MHD disk wind is the dominant feedback is correct also in the case of formation of a massive binary.
Next, we did not study the case with very high accretion rate ofṁ * 10 −2 M ⊙ yr −1 , which would arise in the collapse of very unstable cores of M c 1500(Σ cl /g cm −2 ) −1 M ⊙ (eq. 1). Even though the typical accretion rate of massive star formation is thought to be of the order of 10 −3 M ⊙ yr −1 , there may be cases with higher accretion rates that are especially important for formation of very massive stars. In the case of zero metallicity, Hosokawa et al. (2012) have found a new branch of protostellar evolution at 10 −2 M ⊙ yr −1 : the protostar balloons as r * ≃ 2.6 × 10 3 (m * /100 M ⊙ ) 1/2 R ⊙ without KH contraction. Such "supergiant" protostars have loweffective temperatures of about 5000 K and thus photoevaporation becomes ineffective, which is normally the most important feedback at zero metallicity. A similar phenomenon may appear also at solar metallicities, but it is non-trivial to calculate this evolution because of the presence of metals that alter the protostellar evolution (Hosokawa & Omukai 2009a,b). Due to this uncertainty of protostellar evolution and also the numerical difficulty of calculation of supergiant protostars, we have avoided models with parameter range of 10 −2 M ⊙ yr −1 in this paper. However, even without detailed calculations, we can expect that the MHD disk wind is still the most dominant feedback in the cases of such rapid accretion. This is because the momentum rate of MHD disk wind is proportional to the accretion rate, while the radiation pressure acts similarly as in lower-accretion rate case, and the photoevaporation becomes negligible in the supergiant protostar phase if it appears at solar metallicity. Future work on accurate protostellar evolution calculations with 10 −2 M ⊙ yr −1 at solar metallicities is needed to confirm this expectation.
We also did not consider short timescale variations of accretion rates, which may be induced by disk instabilties, e.g., due to self-gravity. Meyer et al. (2017) simulated the formation of a massive star and showed the accretion bursts occur repetitively. The accretion rate rapidly increases from 10 −4 -10 −3 M ⊙ yr −1 to 10 −1 M ⊙ yr −1 within a duration of 10 years and it recurs with the timescale of several kyr. The accretion burst has a significant impact on the observational aspects, since it results in luminosity outbursts similar to FU Orionis objects. Considering the evolution of the infalling envelope under the influence of feedback, however, we suspect that the accretion burst does not have too significant an impact. This is because the global evolution of the infalling envelope is affected by the accretion rate averaged over the accretion timescale of ∼ 10 4 (m * /10M ⊙ )(ṁ * /10 −3 M ⊙ yr −1 ) yr, which is longer than the expected durations and recurrence timescales of accretion bursts. Accretion bursts would also change the protostellar evolution, since, as described in the previous paragraph, such high-accretion rates can cause a supergiant phase. However, Sakurai et al. (2015) showed that, at least in zero-metallicity case, the supergiant phase cannot last as long as the recurrence timescale of ≥ 10 3 yr since the KH timescale is very short. To check these speculations, further study, deferred to a future paper, is needed to include accretion bursts self-consistently in our modeling.
In addition to further theoretical and numerical studies, better observational tests are needed to confirm the reliability of our theoretical model. We have applied the previous versions of our model to make predictions of observational features using radiative transfer calculations (Zhang & Tan 2011;Zhang et al. 2013bZhang et al. , 2014Tanaka et al. 2016). In a future paper we will model the radiative transfer predictions of the feedback models that we have presented here.

CONCLUSION
We have investigated the impact of multiple feedback mechanisms in massive star formation by core accretion, and calculated the star formation efficiency (SFE) from pre-stellar cores. Our model includes feedback by the outflows driven by momenta from MHD disk winds and radiation pressure, and the effects of mass-loss by photoevaporation and stellar winds. We found the MHD disk wind is the dominant feedback mechanism for all cases considered, while radiation pressure can cause a significant reduction in SFE at the highest masses and especially in lower mass surface density clumps. The obtained SFE can be fitted as 0.4(M c /100 M ⊙ ) −0.115 in the initial core mass range of M c = 10-1000 M ⊙ at the ambient clump mass surface density of 1 g cm −2 , which is a typical value for massive star formation. The gentle decline of M −0.115 c is caused by the radiative feedback which is stronger at higher masses. Therefore, we conclude that the shape of high-mass end of initial stellar mass function, especially potential truncation at m * ∼ 150-300 M ⊙ , is mainly determined by the pre-stellar core mass function and/or disk fragmentation rather than the effects of feedback.
The MHD disk wind provides the major portion ( 90%) of outflow momentum over all the considered mass range, and drives the outflow before the stellar mass reaches about 20M ⊙ , when radiation pressure acting on dust grains in a spherical envelope becomes stronger than the gravitational force. Such radiation pressure was once thought to be a potential barrier for massive star formation, but in more realistic disk accretion and outflow cavity geometries, the strong direct stellar radiation is shielded in the disk-shadowed region, and dust re-emission escapes along the cavities. Therefore, feedback by radiation pressure is not a catastrophic problem for massive star formation. Still, although the total momentum is dominated by the MHD disk wind, radiation pressure also assists to open the outflow cavity to wider angles, since it acts more isotropically than the collimated MHD outflow.
Mass-loss by photoevaporation is strongly suppressed by dust attenuation of ionizing photons. When the protostar starts the Kelvin-Helmholtz contraction at 10-20 M ⊙ , the ionizing photon rate increases with the effective temperature and the photoevaporation starts. However, as the mass-loss rate increases, the photoevaporation flow becomes opaque for the ionizing radiation due to dust opacity. Thus, the photoevaporation mass-loss rate is regulated to 10 −4 M ⊙ yr −1 , which is < 10% of the mass-loss rate without dust attenuation. Since the typical accretion rates of massive star formation are 10 −3 M ⊙ yr −1 , photoevaporation has only a minor impact for the SFE at solar metallicities.
The mass-loss by stellar wind is found to be 10 −5 M ⊙ yr −1 , and thus is not important to set the stellar mass. The stellar wind is, however, very important for the later evolution over timescales of several million years, i.e., the total mass-loss can be tens of M ⊙ .
We also applied our model to primordial, Pop III star formation. Due to the lack of dust grains at zero metallicity, the radiation pressure is negligible and also dust attenuation of ionizing photons does not occur. Therefore, photoevaporation is the major feedback effect in primordial star formation. In our fiducial model, the photoevaporation rate reaches ∼ 10 −3 M ⊙ yr −1 and stops the mass accretion at about 44M ⊙ . In this manner, radiation feedback depends on metallicity, mainly due to the dust absorption. We evaluated that the critical metallicity for two radiative feedback transitions is ∼ 10 −2 Z ⊙ . Dust absorption is effective at higher metallicities than this critical metallicity, which results in radiation pressure being strong while photoevaporation is suppressed. On the other hand, at lower metallicities dust absorption is weak and so radiation pressure eventually becomes negligible and the photoevaporation is more important. Since massive stars are thought to have been astrophysically important since the times of the first stars, more detailed studies are needed to investigate the quantitative effects of feedback as a function of metallicity.