The Impact of Feedback in Massive Star Formation. II. Lower Star Formation Efficiency at Lower Metallicity

We conduct a theoretical study of the formation of massive stars over a wide range of metallicities from 1e-5 to 1Zsun and evaluate the star formation efficiencies (SFEs) from prestellar cloud cores taking into account multiple feedback processes. Unlike for simple spherical accretion, in the case of disk accretion feedback processes do not set upper limits on stellar masses. At solar metallicity, launching of magneto-centrifugally-driven outflows is the dominant feedback process to set SFEs, while radiation pressure, which has been regarded to be pivotal, has only minor contribution even in the formation of over-100Msun stars. Photoevaporation becomes significant in over-20Msun star formation at low metallicities of<1e-2Zsun, where dust absorption of ionizing photons is inefficient. We conclude that if initial prestellar core properties are similar, then massive stars are rarer in extremely metal-poor environments of 1e-5 - 1e-3Zsun. Our results give new insight into the high-mass end of the initial mass function and its potential variation with galactic and cosmological environments.


INTRODUCTION
Massive stars are the main sources of UV radiation, turbulent energy, and heavy elements. Massive closebinaries are the progenitors of merging black holes which have been detected by their gravitational wave emission. Even though they play crucial roles across a wide range of astrophysics, the formation of massive stars is still not fully comprehended. Especially, it is important to understand how the massive star formation process depends on galactic environmental conditions, since this shapes the high-mass end of the initial mass function (IMF) and affects how the IMF may vary through cosmic history.
To address this topic, we have developed a model of feedback during massive star formation relevant for solar metallicity conditions, especially star formation efficiencies from prestellar gas cores (Tanaka et al. 2017b, hereafter Paper I). Here we apply this model to a wide range of metallicities that are expected to be relevant to galactic environments across most of the evolution of the universe.
Radiative feedback has been considered to be critical in setting the mass of massive stars at their birth. Especially, radiation pressure acting on dust grains has been modelled to be a potential barrier to the formation of present-day massive stars. For the idealized case of spherical accretion, radiation pressure acting on the dusty envelope exceeds gravitational attraction when the stellar mass reaches about 20M ⊙ preventing further mass accretion (Larson & Starrfield 1971;Wolfire & Cassinelli 1987). The existence of more massive stars indicates that a non-spherical accretion geometry, i.e., involving accretion disks, is important. Subsequent studies via (semi-)analytic models (Nakano 1989;Jijina & Adams 1996;Tanaka & Nakamoto 2011) and numerical simulations (Yorke & Sonnhalter 2002;Krumholz et al. 2009;Kuiper et al. 2010;Rosen et al. 2016) showed that the gas behind the disk, which is expected to be optically thick, is shielded from radiation pressure and thus accretion can continue to high masses, potentially 100 M ⊙ , depending on the initial condition of the core. Radiation predominantly escapes via lower density cavities above and below the disk. Such low-density cavities may be opened by the stellar radiation itself (Kuiper et al. 2010(Kuiper et al. , 2011, including via radiative Rayleigh-Taylor instabilities (Krumholz et al. 2009;Rosen et al. 2016) or, more likely, by magneto-centifugally-driven outflows (Krumholz et al. 2005;Kuiper et al. 2015Kuiper et al. , 2016Matsushita et al. 2017).
As a result of the above studies, radiation pressure is no longer regarded as a feedback mechanism that is catastrophic for massive star formation. Still, there is a Figure 1. Overview of the evolutionary stages of massive star formation in our model, which is based on the Core Accretion paradigm. (a): The initial prestellar cloud core is spherical and close to virial equilibrium. The structure is characterized by three main parameters: core mass, Mc; mass surface density of ambient clump, Σ cl ; and the ratio of core's initial rotational to gravitational energy βc (McKee & Tan 2003). Here we assume metallicity, Z, may alter feedback effects, but not core structure and accretion properties. (b): In the main accretion phase, the infalling envelope accretes onto the central protostar through the disk. The outflow cavity is opened up by the momentum of an MHD disk wind, with later contributions from radiation pressure, leading to reduction of the solid angle of the region that is able to infall. Additionally, mass loss by the MHD disk wind and photoevaporation reduces the accretion rate onto the star. (c): When infall from the envelope is finished, the disk starts to dissipate by mass accretion onto the star and mass loss caused by the MHD disk wind and photoevaporation. The stellar birth mass, in the approximate limit of formation of a single dominant star, is set when the remnant disk has finally dissipated.
more general remaining question about the quantitative effects of feedback mechanisms in setting star formation efficiencies from gas cores, potentially shaping the stellar IMF and its variation with metallicity.
In the formation of primordial (Pop III) stars, i.e., the limit of zero metallicity, radiation pressure is not expected to be significant because there are no dust grains. Instead of radiation pressure, photoevaporation is thought to be a critical feedback process for setting the mass of Pop III stars. As a massive primordial protostar approaches the Zero-Age Main-Sequence (ZAMS), it starts to emit vast amounts of Lyman continuum photons, i.e., with > 13.6 eV that would ionize infalling and accreting material. The thermal pressure of such ionized gas with 10 4 K drives a photoevaporative flow (Hollenbach et al. 1994), which staunches mass accretion at stellar masses of ∼ 50-100 M ⊙ (McKee & Tan 2008;Hosokawa et al. 2011;Tanaka et al. 2013Tanaka et al. , 2017b. Photoevaporation may also be important in the formation of massive stars in non-zero metallicity environments. Recently, Nakatani et al. (2018) performed radiative hydrodynamical simulations showing the metallicity dependence of the photoevaporation rate. However, their focus is on the dissipation of protoplanetary disks around lowmass protostars with 0.5 M ⊙ . Although there are some similarities, their model is not applicable to our study because the luminosity and the spectrum are quite different between low-and high-mass stars. It is still uncertain how photoevaporation feedback during massive star formation depends on metallicity.
Non-radiative feedback, namely magneto-centifugallydriven outflows, may also be important. In the mass range lower than 10 M ⊙ and in local Milky Way environments, the observed core mass function (CMF) is reported to be similar in shape to the stellar IMF, but with a shift to higher masses by a factor of a few (e.g., André et al. 2010;Könyves et al. 2010). One promising explanation for this is that SFEs from prestellar cores may be ∼ 0.4 for both low-and intermediate-mass star for-mation. Theoretical and numerical studies of low-mass star formation proposed this SFE value is set by outflow feedback that is driven by the momentum of a magnetohydrodynamic (MHD) disk wind (Matzner & McKee 2000;Machida & Matsumoto 2012;Zhang & Tan 2015;Offner & Chaban 2017). In the formation of massive stars, on the other hand, theoretical studies have paid most attention to radiative feedback because of their enormous luminosities. However, observations suggest that the structures of the outflows around low-and highmass protostars are similar (e.g., Qiu et al. 2009;Zhang et al. 2013a;De Buizer et al. 2017;Hirota et al. 2017;McLeod et al. 2018). The models of Zhang et al. (2013bZhang et al. ( , 2014; Zhang & Tan (2018) have considered the formation of massive stars from cores with the only feedback effect included being that due to MHD outflows. Scaling up the assumptions of the model of Matzner & Mc-Kee (2000), they find similar SFEs from massive cores of ∼ 0.5. Matsushita et al. (2017) recently performed MHD simulations of the collapse of massive magnetized cloud cores, ignoring radiative feedback. They showed that an MHD outflow is launched in a similar way to the case of low-mass star formation, but is more powerful due to the higher accretion rate and deeper gravitational potential. Hence MHD outflow feedback is expected to also play an important role in massive star formation.
In reality, massive stars are formed under the influence of all of these feedback processes. Paper I studied the impact of multiple feedback processes in massive star formation using semi-analytic methods and found that MHD disk wind feedback is more important compared to radiative feedback. In this sense and under the assumptions of the modeling via Core Accretion (McKee & Tan 2003) the formation of massive stars is similar to those of low-mass stars. Recently  also studied the combination of multiple feedback processes (disk winds, radiation pressure and photoionization) by radiative-hydrodynamical simulations, which well agreed with our semi-analytic work of Paper I. How-ever, in this paper we focused mainly on present-day massive star formation assuming solar metallicity. In this paper, to investigate how the formation processes of massive stars change with galactic environment and over cosmic history, we extend our model to lower metallicities and evaluate the impact of feedback and SFEs from given prestellar cloud cores. We note that, with a similar conceptual framework, Hosokawa & Omukai (2009b) and Fukushima et al. (2018) have studied radiative feedback and estimated the maximum stellar mass as a function of metallicity. However, they assumed spherical accretion geometry and ignored the MHD disk wind. As described above, disk accretion is a key factor to circumvent the radiation pressure barrier, and MHD disk winds likely have considerable contributions in massive star formation. We adopt an axisymmetric model allowing treatment of these processes, which we will see leads to completely different outcomes from the idealized spherical calculations.
This paper is organized as follows. In §2 we review the basics of our model and introduce the updates needed to treat the effects of feedback processes at a range of metallicities. Then, in §3 we present our results: we show the metallicity dependencies of the feedback processes and the SFE. In §4 we discuss the potential implications for IMF variation based on the obtained SFE model. We conclude in §5.

METHODS
We calculate the accretion history onto massive protostars including effects of several feedback processes. Figure 1 presents a schematic overview of our model. In Paper I, we focused on massive star formation at solar metallicity. Our model is developed under the paradigm of the Turbulent Core Model (McKee & Tan 2003), which is a Core Accretion model scaled-up from those developed for low-mass star formation. We now extend this model in the metallicity range from 10 −5 Z ⊙ to 1 Z ⊙ . The main update from Paper I is taking into account how metallicity influences protostellar evolution and radiative feedback. Here we review the basics components of the model and these updates.

Prestellar Cloud Cores
Our model assumes single star formation from a collapsing prestellar cloud core ( Figure 1a). The initial core is spherical and close to virial equilibrium, being supported by non-thermal pressure components, i.e., turbulence and magnetic fields. The core properties in this model are characterized by three fundamental parameters: core mass, M c ; mass surface density of the ambient clump, Σ cl ; and ratio of the core's initial rotational to gravitational energies, β c . The size of the core is set as R c = 0.057(M c /60 M ⊙ ) 1/2 (Σ cl /1 g cm −2 ) −1/2 pc under the assumption of pressure equilibrium of the core surface with the ambient clump medium, and with the normalization factor for a power-law internal core density profile ρ ∝ r −1.5 (observations of dense cores in Infrared Dark Clouds indicate a density power law index of ≃ 1.3-1.6, Butler & Tan 2012;Butler et al. 2014). In this study, we investigate core masses in the range of M c = 10-1000M ⊙ , fixing the clump mass surface density and the rotational parameter at fiducial values, i.e., Σ cl = 1 g cm −2 (Plume et al. 1997;McKee & Tan 2003;Tan et al. 2014) and β c = 0.02 (Goodman et al. 1993;Li et al. 2012;Palau et al. 2013). We note that we assume that metallicity does not alter the properties of prestellar cores in order to clarify the influence of metallicity on the feedback processes. We will discuss how core properties may vary at various metallicities in §4.

The Main Accretion Phase
The prestellar core undergoes gravitational collapse forming a protostar at its center. Material accretes onto the star through the disk under the influence of multiple feedback process (Figure 1b). In this work, we consider feedback due to MHD disk winds, radiation pressure and photoevaporation, as explained later in this section. In Paper I, we also considered mass loss by stellar winds using a wind model for hot stars with 30, 000-50, 000 K (Vink et al. 2011) and found that stellar wind rates are several orders of magnitude smaller than those by other processes at the solar metallicity. The metal-line-driven stellar winds are expected to be even weaker at lower metallicities. Recently, Vink (2018) showed that, in the case of inflated very massive stars with a relatively low effective temperatures of ∼ 15, 000 K, stellar wind mass loss rates can reach as high as 10 −3 M ⊙ yr −1 at stellar masses > 800(Z/Z ⊙ ) −0.35 M ⊙ . However, in our model setup, stars above 100 M ⊙ always have a high effective temperature of 30, 000 K. Therefore, in this study, the mass loss from stellar winds is not considered.

Infall, Disks and Protostars
The infall of the core is described by the self-similar solution (McLaughlin & Pudritz 1997;McKee & Tan 2003), which gives the infall rate onto the protostar-disk system in the limit of no feedback: where M * d = Ṁ * d dt is the collapsed mass, which is the total mass of the protostar and disk in the limit of no feedback. The obtained infall rate is orders of magnitude higher compared to the typical accretion rates of ∼ 10 −6 M ⊙ yr −1 that are considered to be characteristic of low-mass star formation. However, the actual accretion rate onto the protostar is reduced from the value in Equation (1) due to feedback processes ( §2.2.4).
A disk is formed around the protostar because the initial core is rotating. Assuming angular momentum conservation of infalling gas from the sonic point, where the infall velocity reaches the sound velocity, the disk radius is given by (see §2.1 of Zhang et al. 2014). The disk is massive and self-gravitating due to high supply rate from the envelope. The angular momentum in the disk is efficiently transported by torques in such a massive disk, keeping the mass ratio of the disk to protostar approximately constant at ∼ 1/3 (e.g., Kratter et al. 2008). Therefore, this ratio is fixed at f d = 1/3 during the main accretion phase following the assumption adopted in our previous series of papers (Zhang & Tan 2011;Zhang et al. 2013bZhang et al. , 2014Zhang & Tan 2018, Paper I).
To evaluate the strength of feedback, the properties of the protostar, such as luminosity, radius and effective temperature, along with their evolution, are important. Therefore, we calculate the protostellar evolution selfconsistently given the accretion rate using the model of Hosokawa & Omukai (2009a) and Hosokawa et al. (2010), which solves the basic stellar structure equations, i.e., continuity, hydrostatic equilibrium, energy conservation, and transport (Stahler et al. 1980;Palla & Stahler 1991). This model has also performed successfully at all metallicities in the range of Z = 0-1Z ⊙ (Hosokawa & Omukai 2009b;Fukushima et al. 2018). The stellar boundary condition is adopted from two types depending on the accretion geometry, i.e., spherical or disk accretion. In the earliest stages, the expected disk radius, r d , is smaller than the stellar radius r * and the accretion flow is quasispherical. In this case, a shock front is produced at the stellar surface and a fraction of released flow energy is advected into the interior, which is referred to as the hot shock boundary. On the other hand, if r d > r * , gas accretes onto the stellar surface through a geometrically thin disk. Then much of the energy is radiated away before the material settles onto the star. Thus, in this case, the cold photospheric boundary condition is adopted, in which the specific entropy carried into the star is assumed to be the same as the gas at the stellar photosphere.

MHD Disk Wind and Radiation Pressure
A bipolar outflow sweeps up part of the core reducing the amount of gas that can accrete onto the star. We evaluate the opening angle of the outflow cavity θ esc considering the total momenta of the MHD disk wind and radiation pressure. Matzner & McKee (2000) developed a basic model of outflows driven by the momentum of an MHD disk wind, applied in the context of lowmass star formation. Zhang et al. (2013bZhang et al. ( , 2014 applied this model to the case of high-mass star formation, finding MHD outflow feedback creates outflow cavities that open-up during the course of star formation and set formation efficiencies from the core of ∼ 50%. In Paper I, we introduced the contribution of radiation pressure to this momentum-driven outflow model. The outflow cavity extends to a certain angle if the outflow momentum is strong enough to accelerate the core material to its escape speed in that direction: the following equation is satisfied at the polar angle of θ = θ esc (t) where p dw and p rad are the momenta of the MHD disk wind and radiation pressure, respectively, Ω is the solid angle, v esc = 2GM c /R c is the escape speed from the core. The correction factor c g was introduced by Matzner & McKee (2000) to account for the effect of gravity and the propagation of the shocked shell, which is evaluated as 2.63 for our core setup. As the momenta from the MHD disk wind and the radiation pressure keep accu-mulating, the opening angle of the outflow increases with time until it reaches the maximum angle that is limited by the disk aspect ratio, i.e., θ esc,max = tan −1 (H/r), where H is the disk scale height. Infall can always continue from the equatorial region in the disk shadow because shielding by the inner disk region is efficient at overcoming radiation pressure (Tanaka & Nakamoto 2011;Kuiper et al. 2012). The inner disk aspect ratio is evaluated at the radius of r = 10r * following McKee & Tan (2008). The typical value of the aspect ratio is ∼ 0.1, corresponding to a maximum opening angle of ∼ 84 • . Disk winds driven magneto-centrifugally (Blandford & Payne 1982) are adopted in our study. The mass loading fraction of the wind relative to the accretion rate onto the star is assumed to be f dw = 0.1 as a typical value of disk winds (Königl & Pudritz 2000). Note that, due to the trapping by the core, the actual mass-loss rate in the disk wind is smaller than f dwṁ * , which is the value in the limit of a fully opened outflow cavity. Taking this into account, we describe the disk wind mass-loss rate asṁ where f dw,esc is the fraction of the mass of the wind that can escape from the outflow cavity, which is evaluated based on the the mass flow in the direction θ ≤ θ esc (Zhang et al. 2014). Then, the disk wind momentum p dw is evaluated by integrating the momentum rate of the disk wind,ṗ whereṁ * is the accretion rate onto the star and v K * is the Keplerian speed at the stellar radius. The factor of φ dw is introduced to measure the disk wind momentum toṁ * v K * (Tan & McKee 2002). For our disk wind model, the value of φ dw is about 0.15-0.3. We note that the momentum of the MHD disk wind obtained from our analytic model agrees well with the recent MHD simulation of massive star formation by Matsushita et al. (2017). The angular distribution of the momentum of MHD disk wind is described as (Shu et al. 1995;Ostriker 1997;Matzner & McKee 1999) where θ 0 is a small angle that is estimated to be 0.01 and µ = cos θ (note that 1 0 P dµ = 1). This angular distribution of P (µ) encapsulates the collimated nature of MHD disk winds. Higuchi et al. (2018) followed evolution of magnetic fields in collapsing star-forming clouds using non-ideal simulations before the main accretion phase starts. However, metallicity dependences of MHD disk winds during the accretion phase are still uncertain. Therefore, for simplicity, the MHD disk wind is assumed to be independent of the metallicity in this study. For more details about the disk wind momentum, see also §2.3 of Zhang et al. (2014) and §2.2.1 of Paper I. The momentum from radiation pressure, p rad , is obtained by the integral of the radiation pressure momentum injection rate,ṗ  where f trap is a trapping factor accounting for the optical depth to the stellar radiation and the dust re-emission. In the spherical limit at solar metallicity, the infalling envelope is optically thick not only for the direct stellar radiation but also to the infrared radiation re-emitted from dust grains. Then, the trapping factor is larger than unity, i.e., f trap ≃ τ IR ≫ 1 ( Thompson et al. 2005), boosting the contribution of radiation pressure feedback. However, in non-spherical accretion, this radiation pressure by dust re-emission is reduced significantly by the pre-existing MHD outflow cavity (Krumholz et al. 2005;Kuiper et al. 2015Kuiper et al. , 2016. Therefore, in our models, the effect of dust re-emission is ignored and only direct stellar radiation absorbed by dust grains is considered. In Paper I, we assumed the envelope is optically thick to the stellar radiation and f trap = 1 because we were interested in the solar metallicity case. To treat the low metallicity case properly, we modify f trap as where τ env is the optical depth of the infalling envelope to the direct stellar radiation, κ * acc is the Planck mean opacity at the stellar effective temperature of T * acc † , ρ env is the envelope density, and Σ env is the mass surface density of the envelope from the dust sublimation front, r sub , to the core radius, R c . The envelope density is evaluated from the self-similar solution by McLaughlin † We use the subscript of " * acc" to indicate that we take into account both contributions from the accretion powered radiation and the intrinsic radiation, summed together as the effective total stellar radiation. & Pudritz (1997). The dust-sublimation front can be evaluated as r sub = κ * acc L * acc /[4πσκ sub T 4 sub ], where κ sub is the Planck mean opacity at the sublimation temperature of T sub = 1400 K. The Planck mean opacity at solar metallicity is evaluated using the opacity table by Semenov et al. (2003), and it is simply assumed to be proportional to the metallicity at lower metallicities. At some low level of metallicity, the envelope eventually becomes optically thin reducing the impact of radiation pressure, i.e., f trap ≃ τ env < 1. The angular distribution of the radiation pressure momentum is treated as spherical because only direct stellar radiation is considered, i.e., dp rad /dΩ = p rad /[4π].
The momenta of the MHD disk wind and radiation pressure sweep up the envelope. The sweeping rateṀ swp , i.e., the rate for the envelope material swept up into outflow, is obtained as, where µesc = cos θ esc . Note that the opening angle of the outflow cavity, θ esc , monotonically increases with time, and thus the factor ofμ esc is negative. The mass loss by outflow sweeping is the dominant feedback mechanism in the solar metallicity case (Paper I).

Photoevaporation
Ionizing photons from the central massive protostar irradiate the surface of the disk and the outflow-cavityexposed regions of the infall envelope creating ionized regions. The ionized gas evaporates from the surface if its sound speed is strong enough to escape from the gravitational potential of the protostar. The ionized gas is gravitationally bounded in the inner region, while it evaporates away from the system in the outer region. The critical radius of this transition is named as the gravitational radius, which is usually evaluated as r g = Gm * d /c 2 HII where c HII is the sound speed of the ionized gas (Hollenbach et al. 1994;Tanaka et al. 2013). To take into account the repulsion by radiation pressure, we update this formula (following McKee & Tan 2008), where r g,e+d is the tentative gravitational radius, Γ e+d = (κ T + κ * acc )L * acc /4πcGm * d is the Eddington factor considering both electrons and dust, and κ T is the opacity due to Thomson scattering. As a result of the dust opacity, the tentative gravitational radius can be negative especially in higher metallicity cases. However, r g,e+d should not be smaller than r sub because the dust opacity is included in its estimation. Therefore, the sublimation radius is set as the minimum value of the gravitational radius (Eq. 12). Considering the evaporation speed is the sound speed of the ionized gas, the total photoevaporation rate can be described as, where r 0 (M * d ) is the collapse radius inside which the enclosing mass was originally equal to M * d , and n 0 (r) is the base density at the ionization boundary (Hollenbach et al. 1994). Based on an accurate radiative transfer calculation, Tanaka et al. (2013) provided a base density model in the dust-free case, where S * acc is the ionizing photon rate from the central star, α A is the recombination coefficient for all levels (socalled case A), and c pe ≃ 0.4 is the the correction factor used to match numerical results. In Paper I, we have extended this formula including the absorption by dust grains as where τ pe is the optical depth of the photoevaporation flow due to dust grains, and σ a,d is the absorption cross section of dust grains per H nucleon. For this cross section, we adopt the typical value of σ a,d = 10 −21 cm 2 from the diffuse interstellar medium for the solar metallicity case (Weingartner & Draine 2001), and reduce it by a factor of Z/Z ⊙ for the lower metallicity cases. In the evaluations of the ionizing photon rates, S * acc , and the sound speed of ionizing gas, we take into account of metallicity dependence by using the stellar atmosphere model ATLAS (Kurucz 1991;Castelli & Kurucz 2004) and the spectral synthesis code CLOUDY (Ferland et al. 2013). Following Paper I, we introduce the characteristic optical depth of the photoevaporation flow, which is evaluated from the physical values at the dust sublimation front:τ d = σ a,d r sub n 0 (r sub ).
As will be seen in later sections, this characteristic optical depth is a good measure of the strength of the photoevaporation.

Net Accretion Rates
As described above, several feedback processes combine to reduce the accretion rate below the value it would take in the limit of no feedback (Eq. 1). This reduced accretion rate is expressed aṡ whereṁ d is the mass growth rate of the disk. Since the mass ratio of the disk to the protostar is fixed as f d = 1/3 during the main accretion phase, the mass growth rate of the disk isṁ d =ṁ * /3. The first term of the right hand side represents the infall rate, which is reduced by a factor of µ esc from its no-feedback limit (Eq. 1). Diversion of mass into the outflow by direction injection in the disk wind and by photoevaporation is accounted for by the final two terms. The evolution of the protostar and the accretion structure is solved until accretion is cut off by the feedback processes, i.e.,ṁ * = 0 or the entire natal core collapses, i.e., M * d = M c . In the former case, the stellar mass at its birth m * f is set at that moment. In the latter case, while the main accretion phase has now finished, accretion still continues from the remnant disk until it dissipates.

The Disk Dissipation Phase
The final stage is the disk dissipation phase in which the remnant disk accretes onto the star and/or is blown away by feedback (Figure 1c). In Paper I this phase was ignored, and thus the calculated final stellar masses and the SFEs were minimum values, with maximum fractional errors of 1/3. Now in this paper, we extend the evolutionary calculation until the remnant disk dissipates.
In this phase, the disk mass decreases monotonically because supply from the core infall envelope has ended. The rate of change of the disk mass can be written as, We use Equation (13) with the maximum of the integral range is the disk radius r d to evaluate the photoevaporation from the diskṀ pe . We ignore radiation pressure in this phase since the self-shielding of the disk is expected to be efficient (Tanaka & Nakamoto 2011;Kuiper et al. 2012). To evaluate the accretion rate onto the star, the α-disk model of Shakura & Sunyaev (1973) is introduced, which describes the viscosity as ν vis = αc s H, where c s is the sound speed. Following Kuiper et al. (2010), the viscous parameter and the aspect ratio are fixed in the estimation of the viscosity as α = 0.3 and H/r = 0.1, which is equivalent to the so-called β-viscosity model for self-gravitating disks with β ≃ 0.003 (Duschl et al. 2000). This assumption is reasonable while the remnant disk has non-negligible mass compared to the stellar mass, i.e., m d /m * 0.1. Then, the viscous accretion rate is evaluated asṁ vis = m d r d /ν vis . Viscous accretion powers the MHD disk wind, with its mass flux still assumed to be a fraction f dw of the accretion onto the star, i.e., m vis =ṁ * +ṁ dw = (1 + f dw )ṁ * . Therefore, the mass accretion rate onto the star iṡ Note that the total mass loss rate from the star-disk system is the sum of the MHD disk wind and the photoevaporation, i.e.,ṁ * d = −ṁ dw −Ṁ pe . We solve the evolutionary sequence of the star and the disk until the end of the accretion from the remnant disk, i.e., m d = 0, and finally obtain the stellar mass at its birth m * f . We evaluate the SFEs from prestellar cloud cores with M c = 10-1000M ⊙ , Σ cl = 1 g cm −2 and Z = 10 −5 -1Z ⊙ . Following Paper I, we define the "instantaneous SFE" as the ratio of the accretion rate to the infall rate without feedback, i.e., ε(t) ≡ṁ * /Ṁ * d . The instantaneous SFE is important especially as it can be observable in individual protostellar systems (e.g., Zhang et al. 2016). However, in this series of papers, we focus mainly on the final SFE, rather than the instantaneous SFE, to investigate the relation between CMF to IMF. Therefore, we use "SFE" to refer to the ratio of the final stellar mass to the initial core mass, i.e.,ε * f ≡ m * f /M c .  Figure 2 shows the accretion histories as functions of protostellar mass and time at various metallicities for the initial core mass of M c = 1000 M ⊙ embedded in a Σ cl = 1 g cm −2 clump environment. The accretion rate in the no feedback limit is also shown for reference.

Accretion and Mass-Loss Histories at Various Metallicities
At solar metallicity, the accretion rate increases as the infall rate increases (Eq. 1). When the stellar mass reaches around 30 M ⊙ at a time of 5 × 10 4 yr, the accretion rate starts to deviate significantly below that of the no-feedback limit. However, the accretion rate still continues to rise until m * ≃ 175 M ⊙ (t ≃ 1.4 × 10 5 yr), where the peak accretion rate is about 2 × 10 −3 M ⊙ yr −1 . It then decreases as feedback becomes ever stronger at higher protostellar masses. The decline of accretion is mainly caused by the opening-up of the outflow cavity as was seen in Paper I. The plateau starting at m * ∼ 250M ⊙ (1.8 × 10 5 yr) appears because the opening angle reaches its maximum, limited by the disk thickness. Infall from The evolution of the trapping factor ftrap (top) and the momentum from the MHD disk wind p dw and from radiation pressure p rad (bottom) during the main accretion phase at various metallicities as indicated (for the fiducial cores with Mc = 1000 M ⊙ ). At metallicities ≤ 10 −4 Z ⊙ , the trapping factor becomes lower than unity and thus the radiation pressure momentum is weaker than in higher metallicity cases. Note, the MHD disk wind is the main source of momentum at all metallicities.
the envelope finishes at 290 M ⊙ , which is the end of the main accretion phase (t ≃ 2.6 × 10 5 yr, as indicated by the vertical line in the right panel). In the subsequent disk dissipation phase, the accretion rate decreases as the remnant disk dissipates and a final stellar mass of m * f ≃ 359 M ⊙ is achieved by 6.6 × 10 5 yr. The SFE from the core is thenε * f = 359 M ⊙ /1000 M ⊙ ≃ 0.36. The final stellar mass is higher than that obtained in Paper I (290 M ⊙ ), because accretion during the disk dissipation phase is newly included in this paper.
At lower metallicities, the accretion rate and the final stellar mass become smaller, although the other initial conditions are the same. This is because the impact of the total feedback becomes higher at lower metallicities, which is a trend that is opposite from that in the classic view with idealized spherical accretion. At metallicities lower than 10 −3 Z ⊙ , the accretion history is almost identical. In those cases, the accretion rate starts to drop at m * ≃ 15 M ⊙ (4 × 10 4 yr) and the main accretion phase finishes at m * ≃ 120M ⊙ (2.6 × 10 5 yr). The mass accreted in the disk dissipation phase is negligible unlike in the solar metallicity case. The SFE at 10 −5 Z ⊙ is ε * f = 120 M ⊙ /1000 M ⊙ ≃ 0.12 which is lower than that at Z ⊙ by a factor of three. Figure 3 shows the mass-loss histories at Z ⊙ and 10 −5 Z ⊙ with same initial core mass of M c = 1000 M ⊙ . As presented in Paper I, the dominant mass-loss mechanism at Z ⊙ is outflow sweeping driven by the momentum of the MHD disk wind and radiation pressure. The The evolution of the characteristic optical depth of the photoevaporation flowτ d (Eq. 17) (top) and the photoevaporation rateṀpe (bottom) at during the main accretion phase various metallicities as indicated (for the fiducial cores with Mc = 1000 M ⊙ ). At higher metallicities of 10 −1 Z ⊙ , the increasing photoevaporation rate slows down as the flow becomes optically thickτ d > 1 (shaded region in the top panel). At lower metallicities of 10 −3 Z ⊙ , the photoevaporation evolution converges to a low metallicity limit.
sweeping rate increases and becomes higher than the accretion rate at 30 M ⊙ . At m * ∼ 250 M ⊙ , the sweeping rate drops off as the opening angle reaches close to its maximum limit set by the disk thickness. The photoevaporation rate quickly increases from m * ∼ 10 M ⊙ , however, its rate of increase reduces above 20 M ⊙ . This regulation of the photoevaporation rate is mainly caused by dust absorption of the ionizing photons (Paper I).
On the other hand, in the case of the low metallicity of 10 −5 Z ⊙ , the photoevaporation rate is not limited by dust to be under 10 −4 M ⊙ yr −1 . Instead it overtakes the sweeping rate when the stellar mass reaches 20 M ⊙ . Although the rate of increase becomes smaller at this point, the photoevaporation mass loss rate does still keep rising. The sweeping rate, on the contrary, starts to decreases earlier than is seen in the solar metallicity case. This is because the momentum of the disk wind is powered by mass accretion (Eq. 4), which declines due to the efficient photoevaporation. In this manner, the metallicity changes which is the dominant feedback mechanism: i.e., MHD outflow sweeping at Z ⊙ and photoevaporation at 10 −5 Z ⊙ . The total impact of feedback is thus also altered.
To reveal the causes of the metallicity dependence of the feedback processes, we now discuss how the basic properties of protostars and the flows driven by feedback depend on Z, using the results from initial conditions of M c = 1000 M ⊙ and covering the range Z = 10 −5 -1 Z ⊙ (the corresponding accretion histories were shown in Figure 2). Figure 4 presents the evolution of stellar radii, r * , luminosities, L * acc , and ionizing photon rate S * acc at various metallicities. In the top panel, the basic evolution of the stellar radius is same for all metallicity cases: the protostar swells from 5-8 M ⊙ , then returns via Kelvin-Helmholtz (KH) contraction, and reaches the ZAMS phase at ∼ 30 M ⊙ . Additionally, there are some apparent metallicity dependences. The swelling phase starts earlier at lower metallicity. This is because the swelling occurs when the opacity becomes low enough to redistribute the interior entropy and the interior opacity is lower at lower Z (Hosokawa & Omukai 2009a). Another difference is the radius in the main-sequence phase is smaller at lower metallicity. This is due to the lower abundances of C, N and O atoms. The KH contraction continues until sufficient energy is produced by the CNO cycle, which requires higher temperatures for lower CNO abundances.
Differences that are seen in the radius evolution also influence radiation properties of L * acc and S * acc . In general a smaller radius will lead to a greater accretion luminosity and a hotter photosphere and so more H-ionizing photons per patch of the stellar surface. Still, as presented in Figure 4 middle and lower panels, the relative differences for different metallicities become quite small at higher protostellar masses (although note the larger dynamic range of these panels). Especially, the deviation is smaller when the stellar mass is higher than 20 M ⊙ when the radiative feedback would be significant. In other words, the metallicity dependence of the protostellar evolution does not significantly affect the strength of the radiative feedback in our model calculations.
The top panel of Figure 5 shows the trapping factor, f trap , (Eq. 8). As expected, the trapping is less efficient at lower metallicity: the trapping factor is always f trap = 1 for Z > 10 −3 Z ⊙ , while it becomes less than unity for Z ≤ 10 −4 Z ⊙ . As a result, the momentum feedback due to radiation pressure is less in the lowest metallicity cases than that in higher metallicity cases (bottom of Figure 5). This trend of the trapping factor has significant importance in the classical spherical models. However, in our axisymmetric model, the MHD disk wind is always the dominant source of the momentum that drives the opening of the outflow cavity at all metallicities (bottom of Figure 5). Therefore, the metallicity dependence of the radiation pressure momentum has a minor impact on the total feedback strength and the SFEs.
The top panel of Figure 6 shows the characteristic optical depth of the photoevaporation flowτ d (Eq. 17) and the photoevaporation rateṀ pe as a function of the ionizing photon production rate. In the solar metallicity case, the characteristic optical depth reaches the optically thick regime when S * acc 10 46 s −1 , and thus the increasing rate ofṀ pe is suppressed. On the other hand, in the cases with Z ≤ 10 −3 Z ⊙ , the characteristic optical depth is alwaysτ d ≪ 1. As a result, the photoevaporation rate increases smoothly and reaches higher than 10 −3 M ⊙ yr −1 , which is the typical value of the accretion rate (Figure 2). In this way, photoevaporation has a more and more significant impact as metallicity is lowered, as seen in Figure 3 (see also the analytic argument in §4.3 of Paper I). Figure 7 shows the mass fractions of the final stellar mass, m * f , the outflow mass, M out , and the photoevaporated mass, M pe , compared to the initial core mass at the end of each model calculation. Results are shown for initial conditions with M c = 10, 100, 1000 M ⊙ (all for Σ cl = 1 g cm −2 ) and log Z/Z ⊙ = −5, −4, −3, −2, −1, 0. The mass fraction of m * f (pink bars) is equivalent to the SFE,ε * f . The outflow mass M out is the sum of the timeintegral of the sweeping rate,Ṁ swp , and the mass-loss rate by the MHD wind,ṁ dw , while the photoevaporated mass, M pe , is the time-integral of the photoevaporation rate,Ṁ pe . As seen above, in the case with M c = 1000M ⊙ (right panel), the outflow is the dominant feedback effect at solar metallicity, while photoevaporation becomes significant at Z 10 −2 Z ⊙ reducing the SFE. Similar to Figure 2, it can be seen that all mass fractions are similar in the low metallicity regime, Z 10 −3 Z ⊙ . This is because photoevaporation becomes optically thin to dust opacity at these metallicities ( Figure 6). As the initial core mass decreases, however, the above trend becomes weaker. In particular, the results are almost identical in the case of M c = 10 M ⊙ : the outflow is the dominant feedback mechanism and photoevaporation is negligible. In this lower mass case, the stellar mass is too low to have significant radiative feedback, and thus the effective feedback is only from the MHD disk wind, which is assumed to be independent of the metallicity. Note that the SFE at lower masses are higher than that in Paper I, because the mass accreted during the disk dissipation phase is now included.

Star Formation Efficiencies at Various
Metallicities We have seen that the impact of feedback depends on metallicity and also on initial core mass. To show these trends more clearly, the SFEs at various metallicities are plotted as a function of the final achieved stellar mass m * f in Figure 8.
First, the SFE decreases with the final stellar mass at all metallicities, because radiative feedback becomes stronger in higher-mass star formation. This trend is true even in the solar metallicity case in which the MHD disk wind is the dominant feedback rather than radiative feedback (Paper I). Second, the SFE for m * f 10 M ⊙ is nearly independent of metallicity, since only MHD disk winds are effective feedback in this low-mass regime (left panel of Figure 7). Finally and most importantly in this paper, the SFE for m * f 20 M ⊙ is lower for the lower metallicity cases, and approaches converged results that are independent of metallicty once Z 10 −3 Z ⊙ . This is mainly caused by the metallicity dependence of the photoevaporation rate at higher metallicities, where it becomes suppressed by the presence of dust ( §3.2). This trend of lower SFE at lower metallicity may have potential important implications for systematic variation of the high-mass end of the IMF with galactic environment.
The obtained SFE can be fitted bȳ which agrees with our numerical results within a maximum error of (∆ε * f ) max = 0.03 over the wide ranges of parameters of M c = 10-1000 M ⊙ (all for Σ cl = 1 g cm −2 ) and Z = 10 −5 -1Z ⊙ . This simple fitting formula (and potential generalizations for different Σ cl ) can be applied as a sub-grid model to large-scale simulations of star formation that resolve formation of massive prestellar cores. We note one more interesting finding from our model calculations. Although the SFE decreases with stellar mass at all metallicities, its decline does not show any abrupt cut-off up to about 300 M ⊙ . In other words, for our adopted initial conditions, there is no evidence for an upper limit to the birth mass of stars being caused by feedback, in contrast to spherical models (Hosokawa & Omukai 2009b;Fukushima et al. 2018).

Implications for IMF Variations
Massive stars are short-lived and thus constitute the main source of heavy elements injected into the interstellar and intergalactic media, especially in the early universe. Additionally, their feedback by strong radiation and supernovae affects the dynamical and chemical evolution of galaxies. Therefore, the stellar IMF must be known to predict element production, the impact of feedback, and the formation rate of black holes. However, the universality of the IMF is still under investigation, and it is uncertain if it depends on environment and metallicity (e.g., Bastian et al. 2010).
We thus discuss the importance of feedback processes and their metallicity dependence to IMF variation based on our model calculations. Our model shows that SFE from a core is lower when forming stars of higher mass and under lower metallicity conditions ( §3.2). Considering the IMF to be the multiplicative product of the combination of CMF and SFE, the shape of the highmass end of the IMF is then expected to deviate from the CMF shape in contrast to present-day low-mass star formation Könyves et al. 2010;Cheng et al. 2018).
We can quantitatively link the IMF and CMF based on the obtained SFEs (Nakano et al. 1995;Matzner & McKee 2000). Assuming the CMF and SFE are power- where variables with a subscript of "0" indicate the normalized values, the IMF can be written as, . (23) As the exponent of ε ′ is a negative value of ∼ −0.11-−0.36 (Equation 22) the upper-end IMF slope is steeper than the CMF slope by a factor of ∼ 1.1-1.6 depending on the metallicity. Thus, the number of massive stars is lower than the simple estimation with a constant SFE.
As an example, assuming a CMF slope of α c = 1.35, similar in shape to the Salpeter (1955) IMF from ∼ 1 to ∼ 10 M ⊙ , we evaluate the upper-end IMF at various metallicities based on the fitting of our SFEs ( Figure  9). Here the initial CMF is normalized at 10 M ⊙ , i.e., dN /d ln M c = (M c /10 M ⊙ ) −1.35 . In the solar metallicity case, the upper-end IMF slope is then 1.53, which is a little steeper than the assumed CMF slope. The IMF slope becomes steeper as metallicity decreases, and it converges to α c /(1 − 0.36) ∼ 2.1 at around 10 −3 Z ⊙ as a result of the metallicity dependence of feedback processes ( §3.2). Due to this difference of the IMF slope, the number of stars with 30-100 M ⊙ at 10 −5 Z ⊙ is 2.2 times smaller than the number at Z ⊙ , and that factor is 4.6 for the mass range of 100-300 M ⊙ (assuming the same initial CMF is applied). In this manner, our model predicts that massive stars are relatively harder to form at lower metallicity, especially 10 −3 Z ⊙ .
An interesting observed feature of the high-mass end of the IMF is the maximum stellar mass in young massive clusters. Figer (2005) suggested the upper-mass limit of 150 M ⊙ based on the Arches cluster near the Galactic Center. However, recent observations of the 30 Doradus star-forming region in the Large Magellanic Cloud (LMC) reported very massive stars whose initial masses are estimated to be as high as 200-300 M ⊙ (Crowther et al. 2010(Crowther et al. , 2016. Schneider et al. (2018) reported that the IMF in 30 Doradus has an excess of massive stars  with the slope of 0.90 +0.37 −0.26 , which is shallower than the Salpeter value of 1.35. Since the metallicity in the LMC is about 0.4 Z ⊙ , one might suppose this difference of the maximum stellar mass may come from the effect of the metallicity dependence of the feedback. However, our model showed that the impact of total feedback is stronger at lower metallicity leading to the opposite trend. To reconcile with model results, we therefore speculate that the initial CMFs in these regions may have been different, i.e., there may have been a CMF cut-off in the Arches, and the CMF slope in 30 Doradus may have been shallower than the observed IMF. Our model predicts the CMF slope of the 30 Doradus is as shallow as ∼ 0.77 +0.32 −0.22 assuming a metallicity of 0.4 Z ⊙ (Equations 22 and 23). This speculation suggests that the CMF depends on environmental properties (e.g., Cheng et al. 2018;Motte et al. 2018;Liu et al. 2018), or that other mechanisms, like a small number of protostellar mergers, might influence the formation of the most massive stars.
While the CMF is not the main question in this study, the metallicity dependence of the CMF is crucial to fully understand potential IMF variation in different galactic environments and thus over cosmic history. In primordial star formation, the typical core mass is as high as 1000 M ⊙ or more due to the lack of an efficient coolant (e.g., Bromm & Larson 2004). The analytic model calculation by Omukai et al. (2005) showed that, at the metallicity of 10 −5 Z ⊙ , cloud fragmentation induced by dust cooling mainly creates low-mass fragments with 1 M ⊙ , rather than massive ones (see also simulations by Dopcke et al. 2013). Still, a full understanding of the CMF likely requires accounting for nonthermal processes, such as turbulence and/or magnetic fields (e.g, Padoan et al. 2017;Li et al. 2018).
Another important process to determine the stellar birth mass is disk fragmentation. Although massive cores have many thermal Jeans masses at solar metallicity, catastrophic fragmentation is suppressed by radiative heating by high accretion luminosity and efficient angular-momentum transport by magnetic breaking (e.g., Krumholz et al. 2007;Commerçon et al. 2011). However, small amounts of fragmentation may still occur forming binaries/multiple systems. Indeed, Sana et al. (2012) showed that the more than 70% of observed massive stars have close companions that eventually exchange mass. Tanaka & Omukai (2014) analytically studied the metallicity dependence of the self-gravitational instability of protostellar disks. They found that the protostellar disk is strongly unstable due to efficient dust cooling at 10 −5 -10 −3 Z ⊙ with typical accretion rates of massive star formation, i.e., 10 −4 -10 −3 M ⊙ yr −1 . However, this analysis did not allow for the effects of magnetic fields on disk fragmentation.

Caveats
We have adopted a semi-analytic model that is still highly simplified and idealized, even though it already has some agreements with observations (Zhang et al. 2013a;Tanaka et al. 2016;De Buizer et al. 2017). Below, we discuss some caveats of our modeling.
As described above, we considered only single star formation not allowing for fragmentation. In the formation of present-day massive stars (Krumholz et al. 2009;Rosen et al. 2016) and primordial stars (Stacy et al. 2010;Susa et al. 2014), three-dimensional simulations suggest that fragmentation of protostellar disks leads to the formation of multiple systems. However, magnetic fields are expected to suppress fragmentation (Machida et al. 2008;Commerçon et al. 2011), and so our model may apply in this limit. We expect that our model is still quantitatively appropriate as long as the total stellar mass is dominated by that of the most massive star. On the other hand, if the system contains two or more similar mass stars, our model would need modifications. The momentum rate from MHD disk winds is roughly proportional to the total accretion rate (Equation 5), and thus the number of stars would not significantly alter the MHD wind feedback. In contrast, the total radiative feedback would become weaker in multiple systems, because the luminosity and the ionizing photon rate increase nonlinearly with the stellar mass at least for 100 M ⊙ . Therefore, the total feedback could be somewhat weaker. We adopted the same dust model as at solar neighborhood even for lower metallicity cases, while the dust properties are thought to be different in the early universe. For example, dust grains in the early universe are considered to be produced in supernovae and affected by reverse shocks. These are thought to reduce the fraction of metals in the dust phase destroying especially smallersize grains (Nozawa et al. 2007). The metal fraction in dust would increase to the solar-neighborhood value during the prestellar collapse phase (Chiaki et al. 2013), which may justify our assumption in this work. However, if the dust distribution tends to be biased to large sizes, then the opacity for the ionizing photon would be smaller than our assumed value. Thus, the metallicity at which the photoevaporation becomes significant could be somewhat higher than 10 −2 Z ⊙ compared to our model result. However, the feedback impact at ≤ 10 −3 Z ⊙ of our model would not enhanced by this fact, because it already reaches the saturation level of the low metallicity limit.
Finally, we note that observational tests are needed to confirm the reliability of complex theoretical models. We have applied the previous versions of our models to make predictions on observations of massive protostars at infrared (Zhang & Tan 2011;Zhang et al. 2013bZhang et al. , 2014Zhang & Tan 2018) and at radio (Tanaka et al. 2016(Tanaka et al. , 2017a, and also compared them to observations (Zhang et al. 2013a;De Buizer et al. 2017). We will perform the radiative transfer predictions of the feedback models that we have presented here, and test (at least near solar metallicity cases) with current and future observations, including with Stratospheric Observatory For Infrared Astronomy (SOFIA), Very Large Array (VLA) and Atacama Large Millimeter/ submillimeter Array (ALMA).

CONCLUSIONS
Massive stars are thought to have been astrophysically important since the times of the first stars. Thus we have investigated the impact of several feedback mechanisms in massive star formation and evaluated, by semianalytic methods, the star-formation-efficiencies (SFEs) from prestellar cloud cores. Previously we focused on the formation of present-day massive stars at solar metallicity in Tanaka et al. (2017b) (Paper I). Here we have extended the model to cases with various metallicities of Z = 10 −5 -1 Z ⊙ , as one measure of the effects of galactic environment and cosmic evolution.
We found that the total impact of feedback and which process dominates depends on metallicity. Radiation pressure, which has been regarded as the crucial barrier for present-day massive star formation, has a relatively minor impact over all the metallicity range. At solar metallicities, the MHD disk wind is the dominant mechanism providing a major portion ( 90%) of the outflow momentum. As the metallicity decreases, photoevaporation becomes stronger and reduces the SFE, because dust attenuation of ionizing photons is inefficient. This metallicity dependence saturates at around 10 −3 Z ⊙ .
The obtained SFE from a given core decreases in the formation of higher-mass stars at all metallicities because their feedback is stronger. Moreover, this SFE decline is steeper at lower metallicities due to more efficient photoevaporation (Figure 8). If the initial CMF is described with the Salpeter index of 1.35, our model predicts that the number fraction of stars with 30-100 M ⊙ (100-300M ⊙ ) at 10 −5 Z ⊙ would be 2.2 (4.6) times smaller than that at Z ⊙ . We note that our modeling does not show any clear truncation of SFE at the highest masses. This means that the upper mass limit of stars (if it exists) is not determined by feedback processes and that this applies for all the metallicities we have explored.