Massive Outflows Driven by Magnetic Effects in Star Forming Clouds with High Mass Accretion Rates

TL;DR

This study uses 3D magnetohydrodynamics simulations to explore how magnetic effects influence massive outflows in star-forming clouds with high mass accretion rates, revealing the importance of magnetic strength.

Contribution

It demonstrates that strong magnetic fields are essential for launching massive outflows in high accretion rate environments, extending understanding of star formation mechanisms.

Findings

Massive outflows are driven by magnetic effects in high accretion rate clouds.

Strong magnetic energy leads to comparable mass ejection and accretion rates.

Simulation results align with observed properties of massive star outflows.

Abstract

The relation between the mass accretion rate onto the circumstellar disc and the rate of mass ejection by magnetically driven winds is investigated using three-dimensional magnetohydrodynamics simulations. Using a spherical cloud core with a varying ratio of thermal to gravitational energy, which determines the mass accretion rate onto the disc, to define the initial conditions, the outflow propagation for approximately 10^4 yr after protostar formation is then calculated for several cloud cores. The mass ejection rate and accretion rate are comparable only when the magnetic energy of the initial cloud core is comparable to the gravitational energy. Consequently, in strongly magnetised clouds a higher mass accretion rate naturally produces both massive protostars and massive outflows. The simulated outflow mass, momentum, kinetic energy and momentum flux agree well with observations, indicating that massive stars form through the same mechanism as low-mass stars but require a significantly strong magnetic field to launch massive outflows.

Tables

Table 1. Table 1: Model parameters. Column 1 lists model names. Columns 2–7 give the density enhancement factor f 𝑓 f , diffusion factor C η subscript 𝐶 𝜂 C_{\eta} , cloud mass M cl subscript 𝑀 cl M_{\rm cl} , cloud radius R cl subscript 𝑅 cl R_{\rm cl} , magnetic field strength B 0 subscript 𝐵 0 B_{0} and angular velocity Ω 0 subscript Ω 0 \Omega_{0} , respectively. Columns 8, 9 and 10 give the ratios of thermal α 0 subscript 𝛼 0 \alpha_{0} , rotational β 0 subscript 𝛽 0 \beta_{0} and magnetic energy γ 0 subscript 𝛾 0 \gamma_{0} , respectively, to the gravitational energy of the initial cloud core. Column 11 gives the initial mass-to-flux ratio normalised by the critical value μ 𝜇 \mu .

Model	$f$	$C_{\eta }$	$M_{\mathrm{cl}}$	$R_{\mathrm{cl}}$	$B_0$	${\Omega }_0$	${\alpha }_0$	${\beta }_0$	${\gamma }_0$	$\mu$
			[$M_{\odot }$]	[pc]	[$\mu$ G]	[${10}^{-14}$  s-1]
A	1.4	1	32	0.28	23	3.3	0.5	0.02	0.2	2
B	3.4		77	0.28	56	5.1	0.2
C	8.4		192		140	8.0	0.08
D	17		385		280	11	0.04
E	34		771		560	16	0.02
F	67		1542		1120	23	0.01
CE1	8.4	0.01	192	0.28	140	8.0	0.08	0.02	0.2	2
CE2		0.1
CE3		10
CE4		100
CW1	8.4	1	192	0.28	56	8.0	0.08	0.02	0.03	5
CW2	8.4		192		28	8.0	0.08		0.008	10
EW1	34		771		224	16	0.02		0.03	5
EW2	34		771		112	16	0.02		0.008	10

Equations

$$P=c_{s,0}^2\rho \left[1+{\left(\frac{\rho }{{\rho }_{\mathrm{cri}}}\right)}^{2/5}\right],$$

$$\eta$$

$$M_{\mathrm{out}}={\int }_{t^{\prime }}^{t^{\prime }+\Delta t}\frac{d{M}_{\mathrm{out}}}{dt}dt,$$

$$M_{\mathrm{in}}={\int }_{t^{\prime }}^{t^{\prime }+\Delta t}\frac{d{M}_{\mathrm{in}}}{dt}dt,$$

$$M_{\mathrm{out}}={\int }_{v_r>v_{\mathrm{cri}}}\rho dV,$$

$$P_{\rm out}=\int_{v_{r}>v_{\rm cri}}\,\rho\,|\mbox{\boldmath$v$}|\,dV,$$

$$E_{\rm out}=\int_{v_{r}>v_{\rm cri}}\,\rho\,\mbox{\boldmath$v$}\cdot\mbox{\boldmath$v$}\,dV,$$

$$F=P_{\mathrm{out}}/t_{\mathrm{out}},$$

$$\rho (r)=\frac{{\rho }_{\mathrm{c}}}{(1+(r/r_0{)}^2)},$$

Full text

1 Introduction

Stellar masses are distributed over four orders of magnitude from $\sim 0.01\thinspace M_{\odot}$ to \mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{>}}}100\thinspace M_{\odot}. Observations of the star forming regions in our galaxy reveal that the stellar initial mass function (IMF) peaks at $\sim 0.5\thinspace M_{\odot}$, indicating that the galaxy is mostly populated by low-mass stars (Kroupa, 2002; Chabrier, 2003). The low-mass star formation process has been thoroughly investigated and supported by numerous theoretical works (e.g. Larson, 1969; Masunaga & Inutsuka, 2000). Conversely, because massive stars (with masses \mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{>}}}8\thinspace M_{\odot}) are much rarer than low-mass stars (Kroupa, 2002) and form in regions far from the Sun, their formation process is less clearly understood. However, massive stars can greatly disturb the star forming environment through stellar radiative feedback, injecting kinetic and thermal energy as well as polluted gas into interstellar space through powerful stellar winds and supernova explosions (Zinnecker & Yorke, 2007; Beuther et al., 2007; Tan et al., 2014). Thus, despite being relatively low in number, massive stars significantly affect the evolution of the galaxy, and clarifying their formation process is an important aspect to astrophysical study.

Protostellar outflow offers useful clues for elucidating the star formation process. Such outflows are ubiquitously observed around young low-mass protostars (Bontemps et al., 1996; Cabrit & Bertout, 1992; Wu et al., 2004; Hatchell et al., 2007; Arce et al., 2010), and their driving mechanism, which has been investigated in many theoretical studies, can be briefly summarized as follows. Once a prestellar cloud core begins to collapse, gravitational energy is converted into magnetic and kinetic energy that contribute to the driving outflow either through magnetic pressure or magneto-centrifugal forces around the protostar or circumstellar disc (Blandford & Payne, 1982; Uchida & Shibata, 1985; Pudritz & Norman, 1986; Tomisaka, 2002; Lynden-Bell, 2003; Machida et al., 2008). Protostellar outflow significantly influences the low-mass star formation process through the transport of excess angular momentum from the protostar forming region or circumstellar disc into interstellar space, enabling the infalling gas to reach the protostar. As such, it promotes the mass accretion and protostellar growth during the early star formation phase (Pudritz & Norman, 1983; Konigl, 1989; Wardle & Koenigl, 1993; Tomisaka, 2002). Protostellar outflow also plays an important role in the low-mass star formation process by sweeping up much of the infalling gas in the collapsing cloud and expelling it into interstellar space. Thus, protostellar outflow driven by a low-mass protostar is thought to control star formation efficiency and determine the resultant stellar mass in the cloud (Nakano et al., 1995; Matzner & McKee, 2000; Machida & Hosokawa, 2013).

Recent developments in observational techniques and technologies have enabled the observation of massive star forming regions at high spatial resolutions and have led to the confirmation of many massive outflows believed to be driven by massive protostars (Beuther et al., 2002; Wu et al., 2004, 2005; Zhang et al., 2005; López-Sepulcre et al., 2009; de Villiers et al., 2014, 2015; Maud et al., 2015). Early observations hinted that outflows driven by massive protostars have wider opening angles than those driven by low-mass protostars; consequently, massive and low-mass outflows were originally regarded as different phenomena (Richer et al., 2000; Ridge & Moore, 2001). However, as the observations accumulated, it was confirmed that the structure and shape of the outflows driven by massive and low-mass protostars are very similar (Qiu et al., 2009; Beltrán et al., 2011; Beuther et al., 2010; Palau et al., 2013; de Villiers et al., 2014; Carrasco-González et al., 2015). Observations have also shown that the physical quantities of outflows driven by massive protostars are strongly correlated with protostellar luminosity or mass. In particular, the mass, momentum, momentum flux and kinematic energy of outflows increase with increasing protostellar mass. Consequently, a massive protostar has a massive protostellar outflow. The correlation between the physical quantities of the outflow and the bolometric luminosity or protostellar mass is seen over a very wide mass range from low- to very high-mass protostars (Beuther et al., 2002; Wu et al., 2005; Zhang et al., 2005; de Villiers et al., 2014; Maud et al., 2015). This implies that massive outflows are scaled-up versions of low-mass outflows. Although massive star forming regions are distant and it is not possible to resolve massive protostars themselves, the physical quantities of their massive outflows can be observationally determined because the outflows (\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{>}}}0.1-1 pc) have much larger dimensions than the protostars (\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}0.1 AU). Thus, these massive outflows can provide important clues for understanding the formation processes of massive stars.

The relation between the core mass function (CMF) and IMF may provide further insights into star formation processes. The similarity between CMF and IMF holds over a wide range of stellar masses (Motte et al., 1998; Enoch et al., 2006; Nutter & Ward-Thompson, 2007; Alves et al., 2007; Könyves et al., 2010; André et al., 2010) and indicates that only $\sim 30-50\%$ of the prestellar cloud mass can be converted into stellar mass (i.e. the star formation efficiency is $\varepsilon\sim 30-50\%$). In low-mass star formation, the remaining mass is considered to be ejected from the parent cloud into interstellar space by protostellar outflow. The similarity between CMF and IMF also holds in massive stars and clouds, indicating that massive stars also form from a fraction of the parent cloud mass (Nakano et al., 1995; Matzner & McKee, 2000). As low-mass and massive cloud cores have almost the same star formation efficiency, they likely form by the same fundamental mechanism, which can be surmise to function as follows: a massive cloud core collapses into a protostar, which then grows by mass accretion from the infalling envelope (the so-called core accretion scenario, e.g. Tan et al. 2014). In this case, massive star formation can be most simply characterised by two parameters: the initial cloud mass and the mass accretion rate. A massive cloud core provides the necessary mass reservoir for massive star formation and a large accretion rate enables the rapid formation of a massive star. The latter conditions are naturally satisfied in gravitationally highly unstable cloud cores (Bergin & Tafalla, 2007).

In low-mass star formation, the protostellar outflow emerges immediately after (or before) protostar formation and sustains itself during the mass accretion stage from the gravitational energy released by the infalling gas (Machida & Matsumoto, 2012). A similar evolutionary trend is expected for the massive star formation process, with the outflow from massive protostars also determining (at least partially) the resultant stellar mass. Note that radiative effects are also expected to significantly contribute to the determination of the final stellar mass (for details, see §4.4).

In the low-mass star formation process, outflow is driven by magnetic effects, which also collimate the outflow via a strong toroidal field and hoop stress (Pudritz et al., 2007). Prestellar clouds in low-mass star forming regions have strong magnetic fields, with magnetic energies comparable to their gravitational energies (Crutcher, 1999; Troland & Crutcher, 2008; Crutcher et al., 2010). Such fields are necessary for driving protostellar outflow (Hennebelle et al., 2011; Seifried et al., 2011; Machida & Hosokawa, 2013).

Recently, strong magnetic fields have also been observed in massive star forming regions (Falgarone et al., 2008). Such magnetic fields significantly affect the dynamical evolution of massive clumps and massive star formation (Girart et al., 2009; Tan et al., 2014; Li et al., 2015; Tan et al., 2016). Zhang et al. (2014) showed that, at cloud core scales of $\sim 0.01-0.1$ pc, the magnetic field (rather than turbulence) primarily controls massive star formation. This result is not surprising because the cloud collapses follows the magnetic field lines after turbulence dissipation, even in initially highly turbulent clouds. When the parent cloud of the massive star has a strong field, star formation can be controlled by the protostellar outflow, as is true in low-mass star formation (see §4.4 for radiation feedback effect). Note, however, that when the turbulence is super-Alfvénic, it may sustain for a long time (e.g. Burkhart et al., 2015), in which case it may greatly affect the massive star formation process. We need further observations of massive star forming regions to determine which effect primarily controls massive star formation.

Based on the core accretion scenario, many theoretical studies have investigated the outflow driven by magnetic effects in the low-mass star formation process (e.g. Tomisaka, 2002; Machida et al., 2008; Hennebelle & Fromang, 2008; Price et al., 2012; Tomida et al., 2013; Machida & Hosokawa, 2013; Machida, 2014; Bate et al., 2014; Tsukamoto et al., 2015; Tomida et al., 2015). By contrast, magnetically driven outflow in massive star formation has been investigated in only a few studies (Peters et al., 2011; Commerçon et al., 2011; Hennebelle et al., 2011; Seifried et al., 2012; Kuiper et al., 2015). Moreover, these studies have been limited to cloud cores with a specific cloud size, mass and accretion rate. Thus, the authors of these studies made no direct connections between low-mass (or intermediate-mass) and massive star formation processes.

In this study, we aim to derive the common mechanism of star formation processes using magnetohydrodynamics (MHD) simulations. We investigate the high-mass (or massive) star formation process and clarify the fundamental star formation process applicable to low/intermediate-mass and star formation. Under simple settings, we calculate the evolution of collapsing clouds at varying masses and stabilities (ratio of thermal to gravitational energy of the initial cloud), which control the mass accretion rate and characterise the resultant stellar mass (§3.1). As the outflow is closely related to the gravitational collapse and accretion process (the primary mechanisms in the core accretion scenario), we closely investigate the outflow properties and compare them with observations to validate the core accretion scenario. Note that, as described in §2, we set the prestellar cloud mass to $\geq 32\thinspace M_{\odot}$ and calculate the cloud evolution until the protostellar mass reaches \mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{>}}}2\thinspace M_{\odot}. Thus, this study focuses on intermediate or high-mass star formation. When discussing the star formation process from low- to high-mass stars, we refer to our previous study (Machida & Hosokawa, 2013) in addition to the results of this study.

The remainder of this paper is structured as follows. §2 describes the basic equations and numerical settings of out model and §3 presents the simulation results. In §4, we compare our results with observations and previous works. In §5, we summarise the relation between magnetically driven outflow and the star formation process.

2 Numerical Settings

We evolve gravitationally unstable cloud cores and outflow propagation using resistive MHD simulations that impose a sink around the centre of the computational domain following protostar formation (§2.3). In this way, it is possible to uniquely determine the protostellar mass and the mass accretion rate onto the protostar. However, because we ignore protostellar radiative feedback, which should affect the circumstellar environment, we evolve the system only up to a protostellar mass of $\sim 30\thinspace M_{\odot}$. During the corresponding period, the protostellar feedback is less significant than that at later stages (e.g. Krumholz et al., 2007, 2009). The effect of radiative feedback is discussed in §4.4.

2.1 Basic Equations

To evolve the star forming cloud, we use the resistive MHD nested grid method. The numerical methods and initial settings are nearly identical to those described in our previous studies (Machida & Matsumoto, 2011c, 2012; Machida & Hosokawa, 2013; Machida et al., 2014; Machida, 2014). The resistive MHD equations are given by Eqs. (1)–(4) in Machida & Matsumoto (2012). Instead of solving the energy equation or radiative transfer, we use a simple barotropic equation of state:

[TABLE]

where $c_{s,0}$ is the speed of sound in the initial cloud and $\rho_{\rm cri}$ is the critical density. We set $\rho_{\rm cri}=10^{4}\,\rho_{\rm c,0}$, where $\rho_{c,0}$ is the central density of the initial cloud. The gas behaves isothermally when $\rho<\rho_{\rm cri}$ and contracts adiabatically when $\rho>\rho_{\rm cri}$. The resistivity $\eta$ is that used in our previous works (Machida et al., 2007, 2011b; Machida & Hosokawa, 2013; Machida et al., 2014):

[TABLE]

where $\eta_{\rm fid}$ is determined by solving the chemical reactions of various neutral and charged particles (Nakano & Umebayashi, 1980, 1986; Nakano et al., 2002). To investigate the effect of the dependence of magnetic diffusivity on the calculation results, we change the diffusivity through the parameter $C_{\eta}$, where $C_{\eta}=1$ for fiducial models and $C_{\eta}=0.01-100$ in various other models (for details, see §4.1).

2.2 Initial Condition

A spherically symmetric cloud with a Bonnor–Ebert density profile is adopted as the initial state. We set the initial cloud radius as twice the critical Bonner–Ebert radius (Ebert, 1955; Bonnor, 1956). The Bonner-Ebert density profile is uniquely determined by the central density and isothermal temperature. To convert the non-dimensional parameters of the initial cloud into dimensional values, we assume a central density of $\rho_{\rm c,0}=3.8\times 10^{-19}$ g cm*-3* (or number density $n_{c,0}=10^{5}\,{\rm cm}^{-3}$) and an isothermal temperature of $T_{\rm iso}=40$ K. With these parameters, the initial cloud has a dimensional radius and mass of $R_{\rm cl}=0.28$ pc and $M_{\rm cl}=23\thinspace M_{\odot}$, respectively. After constructing the dimensional cloud with the Bonnor–Ebert density profile, we increase the cloud density by a factor of $f$. The parameter $f$ controls the cloud mass and initial cloud stability $\alpha_{0}$, where $\alpha_{0}$ is the ratio of thermal to gravitational energy. When the cloud density is multiplied by $f$, the mass changes as $M_{\rm cl}=23\thinspace M_{\odot}f$, and $\alpha_{0}$ also changes. Because the mass accretion rate is considered to be proportional to $\propto\alpha_{0}^{-3/2}$, the mass accretion rate is large in an initially thermally unstable cloud with larger $f$ or smaller $\alpha_{0}$. 111 The mass accretion rate can be described as $\dot{M}\propto Mt_{\rm ff}^{-1}$, where $M$ is the prestellar cloud mass and $t_{\rm ff}$ is the freefall timescale. Using the density enhancement factor $f$, the freefall timescale and cloud mass are defined as $t_{\rm ff}=(4\pi Gf\rho_{0})^{-1/2}$ and $M=fM_{\rm Jeans}$, respectively, where the Jeans mass $M_{\rm Jeans}$ is proportional to $c_{s,0}^{3}\rho_{0}^{-1/2}$. Thus, the mass accretion rate is proportional to $f^{3/2}c_{s,0}^{3}$. As the parameter $\alpha_{0}$ is proportional to $f^{-1}$ (or $M^{-1}$), the mass accretion rate is proportional to $\alpha_{0}^{-3/2}c_{\rm s,0}^{3}$.

The model name, parameter $f$, cloud mass $M_{\rm cl}$, cloud radius $R_{\rm cl}$ and $\alpha_{0}$ for each model are listed in Table 1. The parameter $C_{\eta}$ in the third row of Table 1 controls magnetic diffusivity (for details, see §4.1). The parameters $f$ ($\alpha_{0}$) are varied as $f=1.4-67$ ($\alpha_{0}=0.01-0.5$). The mass accretion rate between models A ($\alpha_{0}=0.5$) and F ($\alpha_{0}=0.01$) should differ by approximately 350 ((0.5/0.01)3/2). The cloud masses in the models range from $32$ to $1542\thinspace M_{\odot}$.

A uniform magnetic field $B_{0}$ is imposed over the entire computational domain and the initial cloud is subjected to a rigid rotation $\Omega_{0}$. The initial magnetic field strength, $B_{0}$, and angular velocity, $\Omega_{0}$, are also listed in Table 1, along with the ratios of rotational ($\beta_{0}$) and magnetic ($\gamma_{0}$) energies to the gravitational energy. The last row of Table 1 gives the mass-to-flux ratio $\mu$ normalised by the critical value $(2\pi G^{1/2})^{-1}$. In models A–F (fiducial models) and CE1–CE4 (different diffusivity models, see §4.1), the magnetic field strength $B_{0}$ is adjusted to have the non-dimensional parameter $\mu=2$ ($\gamma_{0}=0.2$). Thus, models with different $f$ have different magnetic field strengths. As seen from Table 1, different angular velocities are used in models with different $f$ to fix the ratio of rotational to gravitational energy ($\beta_{0}=0.02$) among the models. To investigate the effects of magnetic field strength on outflow properties, we develop four weak magnetic field models, CW1, CW2, EW1 and EW2 (for details, see §4.2).

2.3 Numerical Method

To calculate considerably different spatial scales, we use a nested grid method in which each grid has the same number of cells (i, j, k) = (64, 64, 32) and impose mirror symmetry with respect to the $z=0$ plane (the nested grid method is detailed in Machida et al. 2005a, b). Six grid levels ($l=1-6$) are initially prepared and the initial cloud cores are embedded in the fourth grid level ($l=4$), with a box size of $L_{4}=0.56$ pc ($=2R_{\rm cl}$) and a cell width of $h_{4}=8.75\times 10^{-3}$ pc. To avoid artificial refraction of Alfén waves at the computational boundary, we set a large interstellar space outside the initial cloud core. In each model, the interstellar media outside the initial cloud core have a uniform density of $\rho_{\rm ISM}=7.6\times 10^{-21}$ g$\,{\rm cm}^{-3}$ ($n_{\rm ISM}=2\times 10^{3}\,{\rm cm}^{-3}$) and a uniform magnetic field $B_{0}$. Note that gravity (comprising gas self-gravity and protostellar gravity) is imposed only in the region $r<R_{\rm cl}$ (for details, see Machida et al., 2010a; Machida & Matsumoto, 2011c, 2012; Machida & Hosokawa, 2013). The outer computational boundary is imposed on the surface of the first grid level $l=1$. The box size of the first grid level ($L_{1}=4.4$ pc) is 16 times the initial cloud radius. The grid level is dynamically increased to maintain the Truelove condition (Truelove et al., 1997) in which the Jeans wavelength is resolved to at least 16 grid points. Both the box size and cell widths are halved with each increment of the grid level. The maximum grid level ($l=15$) has a box size of $L_{15}=56$ AU and a cell width of $h_{15}=0.8\,$AU.

To avoid a prohibitively short time step, we impose a sink at the centre of the computational domain. Protostar formation begins when the central density $n$ reaches the threshold density $n_{\rm thr}$ (set here to $n_{\rm thr}=10^{13}\,{\rm cm}^{-3}$). The threshold density satisfies the necessary condition for resolving the disc forming region (Machida et al., 2014). Once the protostar is formed, we remove the gas exceeding $n>n_{\rm thr}$ in the region $r<r_{\rm sink}$, where $r_{\rm sink}$ is the sink (accretion) radius set to $1$ AU in this study.

In models E and F, we evolve the cloud until the protostellar mass reaches $M_{\rm ps}\sim 30\thinspace M_{\odot}$. In other models, we stop the calculation when $M_{\rm ps}<30\thinspace M_{\odot}$, because the accretion rate in these models is too slow to reach this mass within a reasonable time frame (typically, the calculations require from two to four months of wall-clock time).

3 Results

3.1 Mass Accretion Rate

Figure 1 plots the mass accretion rate onto the protostar against the protostellar mass in models A–F. The figure indicates that the initially more unstable clouds have higher accretion rates. Within the range $M_{\rm ps}<0.1-1\thinspace M_{\odot}$, the mass accretion rate is very similar in all models because it depends on the properties of the first core (Larson, 1969; Masunaga & Inutsuka, 2000), which is supported by thermal pressure, at the very early stage. Just before protostar formation, the thermal pressure decreases inside the first core (Larson, 1969); just after protostar formation the central region of the first core rapidly falls onto the protostar with a higher accretion rate. Note that the outer region of the first core (or the first core remnant) is supported by centrifugal force and evolves into a rotationally supported disc (Bate, 1998; Walch et al., 2009; Machida et al., 2010a; Bate, 2011; Walch et al., 2012; Tsukamoto & Machida, 2013).

The mass accretion rate differs noticeably among the models once the protostellar mass exceeds $\sim 1\thinspace M_{\odot}$, with a mass accretion rate of two orders of magnitude higher in model F than in model A. After rotationally supported disc formation, the gas falls onto the disc and gradually moves inward towards the protostar. Eventually, it falls onto the protostar (or sink). Thus, the mass accretion rate onto the protostar is roughly determined by the mass accretion rate onto the disc, which is itself determined by the initial cloud stability (as described in §2.2). Figure 2 shows the averaged mass accretion rate for each model against the initial cloud stability $\alpha_{0}$, with the mass accretion rate averaged from just after protostar formation until the calculation ends (circle symbols) or until twice the freefall timescale elapses (cross symbols). The figure indicates that the averaged mass accretion rate in each model corresponds closely to the theoretical prediction $\sim\alpha_{0}^{-3/2}\,(c_{s}^{3}/G)$.

Figure 1 also shows that the mass accretion rate drops in the range of 0.4\thinspace M_{\odot}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}M_{\rm ps}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}1\thinspace M_{\odot}. Masunaga & Inutsuka (2000) showed that the mass accretion rate is high in the very early main accretion phase because the gas is accreted from the inner region of the first core (remnant), where the gas temperature is high (the mass accretion rate is proportional to $\propto c_{\rm s}^{3}$) and is primarily supported by the thermal pressure gradient force. Then, after the thermally supported region of the first core is accreted, the gas accretes from the infalling envelope through the rotationally supported disc and the mass accretion rate drops. The oscillation of the mass accretion rate seen in Figure 1 is caused by gravitational instability in the disc, as will be described later.

3.2 Evolution of Protostellar Outflow in Typical Models

Figures 3 and 4 show the time evolutions of the density and velocity distributions in models A, C and F. In model A (left column of both figures), the outflow already extends up to $\sim$2000 AU at the protostellar formation epoch. As the accretion rate in model A is relatively low (Fig. 1), the first core sustains over a long duration. Note that the mass accretion rate onto the first core determines its lifetime (Tomida et al., 2010). Initially, the wide-angle outflow is driven by the first core (Tomisaka, 2002). Once the protostar has formed, the outflow is driven by the rotationally supported disc that evolves from the first core or its remnant. As the transition from the first core to the rotationally supported disc is smooth (Bate, 2011; Machida & Matsumoto, 2011c), the outflow continues to be driven by the central region of the collapsing cloud after the first core formation, as seen in the left column of Figures 3 and 4.

At the protostar formation epoch, the outflow extends up to approximately 500 AU in model C (central column of Figs. 3 and 4) but is absent in model F (right column of Figs. 3 and 4). A higher accretion rate rapidly increases the mass (and hence the central density) of the first core, and the second collapse occurs shortly afterward. Therefore, a higher accretion rate shortens the lifetime of the first core (Saigo & Tomisaka, 2006; Saigo et al., 2008). In model F, the mass accretion rate is high and the lifetime of the first core is very short; consequently, the protostar forms before the magnetic field sufficiently twists and outflow driving occurs. Thus, in model F the outflow evolves only after protostar formation. We therefore attribute the differences in the outflow evolution or size at the protostar formation epoch (first row in Figs. 3 and 4) to differences in the first core lifetime, which is in turn related to the mass accretion rate and initial cloud stability, $\alpha_{0}$.

Figures 3 and 4 also show that the structures and shapes of the outflows are similar across the models. In each model, the opening angle increases with time over a small scale (or disc scale) around the disc (Fig. 3), while a well-collimated outflow structure appears at the cloud core scale of $\sim 5000$ AU (Fig. 4). The opening angle in Figures 3 and 4 strongly depends on both the evolutionary stage and the observed scale. Even when the outflow grows sufficiently (see bottom row of Fig. 4), it is recognizable as a very wide-opening angle flow over a very small scale around the protostar (see bottom row of Fig. 3). Note that, as the purpose of this study is to relate the accretion rate to the outflow rate (and not to investigate the evolutionary properties of outflows among models), we do not quantitatively compare the outflow shape and collimation degree among the models.

3.3 Evolution of Circumstellar Disc in Typical Models

Figure 5 shows the time sequence of the density and velocity distributions on the equatorial plane in models A, C and F, confirming that a higher accretion rate produces a larger disc. In model F, which has the largest accretion rate among the models (see Fig. 1), a large rotationally supported disc forms and then fragments.

To quantitatively compare the disc properties among the models, we plot the radii and masses of the rotationally supported disks in Figure 6. As a result of trial and error, we found that the disks contain all cells that fulfill the following criteria (Machida et al., 2010a; Joos et al., 2012):

the number density exceeds $n_{\rm disk}>10^{9}\,{\rm cm}^{-3}$; 2. 2.

the rotational velocity is five times larger than the infall (or radial) velocity ($|v_{\phi}|>5|v_{\rm r}|$) and larger than 90% of the Keplerian velocity ($v_{\phi}>0.9\,v_{\rm kep}$). The Keplerian velocity $v_{\rm kep}$ is defined as $v_{\rm kep}=(GM_{\rm tot}/r)^{1/2}$, where $M_{\rm tot}=M_{\rm gas}+M_{\rm ps}$ is the sum of the high density gas $M_{\rm gas}=\int^{n>10^{11}\,{\rm cm}^{-3}}\rho\,dV$ and the protostellar mass $M_{\rm ps}$; 3. 3.

the rotational energy ($\rho v_{\phi}^{2}/2$) is larger than the thermal pressure $P$ by a factor of two.

From Figure 6 it is seen that the models with larger accretion rates produce large massive disks in the early evolutionary stage. The oscillation of radius and mass seen in the figure is caused by the appearance and disappearance of spiral arms in the disc (see below). The rapid decrease of mass in model F is caused by the infalling of fragments to the protostar; as discussed below, fragmentation occurs only in model F. Later, we will show that the disc size also depends on the magnetic field strength and the magnetic dissipation process (§4.1).

When the magnetic field is absent or weak, a high accretion rate causes vigorous fragmentation (e.g. Hennebelle et al., 2011) because the disc surface density increases within a short time, while the Toomre’s $Q$ value ($\equiv c_{s}\kappa/(\pi G\Sigma)$,Toomre 1964) decreases. However, under a strong magnetic field, outflow and magnetic braking will effectively transport the angular momentum outward, promoting mass accretion from the disc to the protostar (Seifried et al., 2011; Joos et al., 2012, 2013; Machida et al., 2014). Therefore, the disc surface density decreases, the $Q$ value increases, and the disc consequently stabilises. In the massive star formation process, vigorous fragmentation within the disc (or over a small scale) should be suppressed because each fragment would be expected to evolve into a low- or intermediate-mass star without forming a (very) massive star (see also Hennebelle et al., 2011; Peters et al., 2011; Myers et al., 2013). Note that fragments tend to fall onto more massive protostars; therefore, a few massive protostars may continue to acquire mass and evolve into very massive stars (Machida et al., 2011a; Vorobyov & Basu, 2015; Machida & Nakamura, 2015; Sakurai et al., 2016). With an initially strong magnetic field (models A–F), fragmentation occurs when the mass accretion rate is as high as \dot{M}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{>}}}10^{-2}\thinspace M_{\odot} yr*-1* (model F), while fragmentation is suppressed with an accretion rate of \dot{M}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}10^{-2}\thinspace M_{\odot} yr*-1*(models A–E). However, there are few fragmentation events in model F, and therefore this system should evolve into a binary or a few-star multiple system rather than a small stellar cluster comprising several tens of stars. Note that, as the fragmentation condition and Toomre’s $Q$ parameter depend on the gas temperature, we should also carefully consider fragmentation with thermodynamics in addition to the mass accretion rate or disc surface density (see also §4.4).

Although no fragmentation occurs in models A–E, a higher accretion rate induces gravitational instability and establishes a spiral structure in the disc. Figure 7 shows the time sequence of density distributions on the equatorial plane in model D. Once the spiral structure develops (Fig. 7a), the gas in the disc rapidly falls onto the central protostar under the effective angular momentum transfer via the gravitational and pressure torques. The disc transiently stabilises against the self-gravity (Fig. 7b), until the high mass accretion rate onto the disc increases the disc surface density, and the spiral structure reforms (Fig. 7c). The gravitational instability of the disc explains the non-steady accretion seen in Figure 1. This phenomenon has been thoroughly investigated in low-mass star formation processes (e.g. Vorobyov & Basu, 2006; Tomida et al., 2017). Thus, both magnetic effects and gravitational instability contribute to the angular momentum transfer.

3.4 Mass Inflow and Outflow Rates

To investigate the efficiency of mass ejection, we plot the mass inflow and outflow rates in Figure 8. These rates are estimated on the cell surfaces of different grids (for details, see Tomisaka 2002). On the $l=9$ grid (Fig. 8a), the inflows and outflows cover 1/32 of the initial cloud (3500 AU), whereas on the level $l=12$ grid (Fig. 8d) they cover the entire disc region (450 AU). Clearly, the infall rate onto each grid depends on the cloud parameter $\alpha_{0}$ and strongly affects the mass outflow rate (the higher the accretion rate, the higher the outflow rate). Figure 8 also shows that the respective infall rate of the models do not significantly depend on the scale of the collapsing region (compare the infall rate of each model in panels (a)-(d)). Moreover, the outflow rate is time-variable in all models and is exaggerated on small scales (or disc scales). This indicates that the outflow is intermittently driven by the circumstellar disc. As shown in Figure 7, mass accretion on the protostar is promoted by a non-axisymmetric structure that develops in the disc. Thus, the disc time variability caused by gravitational instability also significantly affects outflow driving. On the disc scale (Fig. 8d), the outflow time variability is significant in models with higher accretion rates (models E and F). The high accretion rate gravitationally destabilises the disc over short time intervals, generating a highly time-variable outflow.

Figure 9 plots the ratios of outflow to inflow masses, averaged every 50 (left panels) and 1000 yr (right panels). To derive the ratio $M_{\rm out}/M_{\rm in}$, the inflow $M_{\rm in}$ and outflow $M_{\rm out}$ masses of each grid are, respectively, estimated as

[TABLE]

where $dM_{\rm out}/dt$ and $dM_{\rm in}/dt$ are the mass inflow and outflow rates of each grid (Fig. 8) and $\Delta t=$50 yr (left panels) or 1000 yr (right panels). At first, the protostar formation epoch $t_{\rm ps}=0$ is set to $t^{\prime}$ and the mass inflow and outflow rates are integrated over $\Delta t$. Then, $t^{\prime}=t^{\prime}_{\rm pre}+\Delta t$ is used as the next integration and the integration is repeated up to the end of the simulation data, where $t^{\prime}_{\rm pre}$ is the previous lower limit of integration. The outflow is driven near the disc. On the disc scale (Figs. 9c and d), the ratios range from approximately 0.2 to 0.4, indicating that about 20-40% of the infalling gas is ejected from the disc as outflow. The ejection rate estimated from our study corresponds closely to previously reported rates (Tomisaka, 1998, 2002; Hennebelle & Fromang, 2008; Seifried et al., 2012). In the large scale, the ratio of outflow to inflow mass increases with time because the outflow opening angle at small scales widens with time, sweeping a larger fraction of the infalling gas (see Figs. 3 and 4). The ratios reach or exceed 50% (Figs.  9a and b) on the large scales (3600 AU), implying that the protostellar outflow finally ejects over half of the infalling gas from the star-forming cloud into interstellar space.

3.5 Comparison with Observations

The process of protostellar outflows in low-mass star formation has been well investigated in many observational and theoretical studies and is known to be driven by magnetic effects (e.g. Konigl & Pudritz, 2000; Pudritz et al., 2007; Arce et al., 2007; Frank et al., 2014). Although there have been only a limited number of studies focusing on massive outflow in the massive star formation process, these are also considered to be driven by magnetic effects (Hennebelle et al., 2011; Seifried et al., 2012). Recent observations also imply the importance of the magnetic field and protostellar outflows in massive star forming regions (Tan et al., 2013; Li et al., 2015; Tan et al., 2016). Therefore, to better clarify the massive star formation process, we need to further investigate the effects of the magnetic field on massive outflows. This study focuses on the relations between the outflow and accretion rates and the outflow properties, with results that can hopefully provide insights into massive star formation. For this reason, we compare our results with observed massive outflows in massive star forming regions.

Figure 10 shows the outflow mass $M_{\rm out}$, momentum $P_{\rm out}$ and kinetic energy $E_{\rm out}$ as functions of time (left panels) and protostellar mass (right panels). The results are presented as solid lines. In our simulations, these quantities are, respectively, estimated as

[TABLE]

[TABLE]

[TABLE]

where we integrate each quantity at each grid point (over the entire computational domain) with radial velocity $v_{r}>v_{\rm cri}$. Through trial-and-error, we determined $v_{\rm cri}$ to be $0.5$ km s*-1*. In this determination, we estimated the outflow physical quantities only without integrating other components such as the circumstellar disc. Note that, by setting $v_{\rm cri}=0.5$ km s*-1*, we somewhat underestimate the outflow physical quantities because we ignore the low-velocity outflow components ($0<v_{\rm out}<v_{\rm cri}$). As a reference, the left panels of Figure 10 plot the most recent observational data taken from de Villiers et al. (2014). These observations are plotted against the outflow dynamical timescale $t_{\rm dyn}$ (defined as the outflow length $L_{\rm out}$ divided by the typical outflow velocity $v_{\rm out,typ}$, namely $t_{\rm dyn}=L_{\rm out}/v_{\rm out,typ}$). We find that the simulated outflows and their maximum ($v_{\rm out,max}$) and typical ($v_{\rm typ}$) velocities strongly vary with time. Thus, we cannot easily understand the time evolution of the outflow physical properties in light of the dynamical timescale of the simulation because the dynamical timescale ($t_{\rm dyn}=L_{\rm out}/t_{\rm out,max}$ or $L_{\rm out}/t_{\rm out,typ}$) is not a simply increasing function of time. For this reason, the definition of time (the abscissa axis in the left panels of Fig. 10) differs between our calculations (elapsed time after the protostar formation) and observations (outflow dynamical timescale), and further effort may be therefore be needed to make a valid comparison.

As shown in Figure 10, the simulated quantities gradually increase as functions of both time (left panels) and protostellar mass (right panels). Models with higher accretion rates or smaller $\alpha_{0}$ increase more rapidly over time (Fig. 10, left panels). This result can be explained in the context of the ratio of the outflow to the inflow rate. As shown in Figure 9, this ratio is almost identical among the models and, therefore, increasing the accretion rate naturally generates a massive outflow with high momentum and kinetic energy. Because these simulated outflow quantities continue to increase, they will reach the observed values of de Villiers et al. (2014) at a further evolutionary stage.

On the other hand, each quantity seems to converge to a certain value in the right panels of Figure 10, in which the simulated quantities are plotted against the protostellar mass. The protostellar mass roughly corresponds to the mass accretion rate (or inflow rate) multiplied by the elapsed time after protostar formation ($M_{\rm ps}\sim\dot{M}\times t_{\rm ps}$). The convergence in these physical quantities means that the outflow properties are primarily determined by the mass falling onto the central region, which is an expected result as the outflow is powered by the release of gravitational energy, which is proportional to $\dot{M}\times t_{\rm ps}$.

3.5.1 Outflow Mass

We first compare the respective outflow masses produced by the simulations and observations. In a statistical study of outflow observations, Wu et al. (2005) showed that the outflow mass is widely distributed over the range 10^{-3}\thinspace M_{\odot}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}M_{\rm out}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}10^{3}\thinspace M_{\odot}, which includes both low- and high-mass stars. Beuther et al. (2002) observed 26 high-mass star forming regions and identified 21 outflows driven from protostars, with masses ranging from $0.04$ to $1423\thinspace M_{\odot}$. Nearly all of these outflows massed from $\sim 1$ to $100\thinspace M_{\odot}$ (average $101\thinspace M_{\odot}$), with the wider range attributable to a few extremely massive and extremely low-mass outflows. Zhang et al. (2005) identified 35 outflows among 69 luminous $IRAS$ sources. They showed that the masses of massive protostellar outflows were tens of $\thinspace M_{\odot}$, with an average of $M_{\rm out,ave}=20.6\thinspace M_{\odot}$. Recently, Maud et al. (2015) observed 99 massive young stellar objects and confirmed massive outflows with a mass range of $\sim 10-100\thinspace M_{\odot}$. de Villiers et al. (2015) listed 44 molecular outflows appearing in 13CO lines and 6.7-GHz methanol masers. They plotted the outflow mass distributions in their own and previous studies (Beuther et al., 2002; Wu et al., 2004; Zhang et al., 2005) and reported massive outflows with typical masses of $\sim 10-100\thinspace M_{\odot}$ and an overall mass range of 10\thinspace M_{\odot}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}M_{\rm out}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}10^{3}\thinspace M_{\odot} in massive protostars. The majority of their investigated outflows also massed $\sim 100\thinspace M_{\odot}$. Note that, as massive outflows are easily observed, they may have been preferentially found in these observations.

As already described, a massive protostar has a massive outflow. This tendency does not contradict the observed star formation efficiency of $\sim 30-50$% (e.g. André et al., 2010). Observations have also shown that the outflow mass increases with protostellar luminosity or protostellar mass (Wu et al., 2004; de Villiers et al., 2014; Maud et al., 2015). In the main accretion phase, protostellar luminosity is proportional to $\propto M_{\rm ps}\dot{M}$. Thus, the observations imply that massive outflows arise from massive protostars or high mass accretion rates. Although we terminate the calculation when the protostellar mass reaches $\sim 30\thinspace M_{\odot}$ in models with higher mass accretion rates, the simulated outflow masses ($\sim 20-30\,\thinspace M_{\odot}$) agree reasonably with the observations ($\sim 10-100\thinspace M_{\odot}$).

3.5.2 Outflow Momentum and Energy

To more precisely compare the simulated and observed results, we plot the momenta and energies of the outflows derived by simulations in Figures 10b, e, c and f. Figures 10b and c also show the observed outflow momenta and energies from de Villiers et al. (2014). In their statistical study of both low- and high-mass protostars, Wu et al. (2005) showed that the momenta of observed outflows range from $10^{-3}$ to $10^{4}\thinspace M_{\odot}\,{\rm km}\,{\rm s}^{-1}$. In our simulation, the outflow momentum depends on the time and on $\alpha_{0}$, but the massive protostars have outflow momenta of 3\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}P_{\rm out}/(\thinspace M_{\odot}\,{\rm km}\,{\rm s}^{-1})\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}1000 (Fig. 10b and e). Our derived outflow momenta and kinetic energies closely agree with other observations of massive outflows (Beuther et al., 2002; Zhang et al., 2005; Maud et al., 2015); however, they are somewhat smaller than the observed values. Nevertheless, the simulated outflows continue to be driven, indicating that the outflow momentum and kinetic energy would increase at a further evolutionary stage. Similarly, our derived massive outflow kinetic energies approach those in the observations (Fig. 10c) and would be expected to converge to the observed values under continuation of the simulations. In addition, the outflow kinetic energies derived from our simulations are comparable to observations by Beuther et al. (2002), Zhang et al. (2005) and Maud et al. (2015). Note that, to estimate the outflow kinetic energy in the late accretion stage, we need to consider the effect of radiative feedback. Note also that observational studies will preferentially detect more massive (or more evolved) outflows. Thus, for a proper comparison of simulated and observed massive outflows, we require longer-term simulations or more spatially resolved observations (which would identify younger outflows).

3.5.3 Outflow Momentum Flux

The outflow momentum flux $F$ can be used to understand the driving mechanism or engine of the outflow (Bally & Lada, 1983; Cabrit & Bertout, 1990, 1992; Bontemps et al., 1996). Observations have determined the outflow momentum flux as $F_{\rm out}=P/t_{\rm dyn}$ (note that the outflow momentum flux has the dimensions of force). The momentum flux in each of our simulated models, derived as

[TABLE]

is plotted in Figure 11. Here, $t_{\rm out}$ is the elapsed time after the outflow emerges. The left panel of Figure 11 also plots the momentum fluxes listed in Table 5 of de Villiers et al. (2014).

After the outflow emergence, the outflow momentum flux in each model remains nearly constant, with a slight increase when M_{\rm ps}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{>}}}1\thinspace M_{\odot} except in model A. Clearly, the momentum flux is larger in models with higher accretion rates (or smaller $\alpha_{0}$) than in models with lower accretion rates (or larger $\alpha_{0}$), which is to be expected because the gravitational energy is converted into magnetic and rotational (or kinetic) energies and the outflow is driven by the Lorentz and centrifugal forces (Blandford & Payne, 1982; Uchida & Shibata, 1985). Thus, as the accretion rate increases, a strong magnetic field and rapid rotation develop around the protostar or circumstellar disc, yielding a powerful outflow.

The simulated momentum fluxes range from $10^{-4}$ to $0.3\thinspace M_{\odot}\,{\rm km}\,{\rm s}^{-1}\,{\rm yr}^{-1}$. The simulation results are consistent with the observations of de Villiers et al. (2014) (Fig. 11, left panel) and with those of Beuther et al. (2002), Zhang et al. (2005) and Maud et al. (2015). From these agreements, the driving mechanism of massive outflow may be explained as analogous to the magnetically driven outflow in the low-mass star formation process. Because the momentum flux does not (significantly) depend on the evolutionary stage, we may estimate the mass accretion rate from the observed momentum flux.

4 Discussion

4.1 Effects of Magnetic Dissipation on Outflow Driving

Whereas the protostellar outflow is considered to be driven by magnetic effects (magnetocentrifugal or magnetic pressure-driven winds), the magnetic field dissipates in the high-density gas region in which the ionisation degree is very low (e.g. Nakano & Umebayashi, 1980, 1986; Nakano et al., 1995, 2002; Tomida et al., 2013; Tsukamoto et al., 2015a, a; Tomida et al., 2015). Thus, to investigate the outflow driving process, we must carefully consider the magnetic dissipation. Here, we include ohmic dissipation and adopt the dissipation coefficient formulated in Machida et al. (2007) and Machida et al. (2008). The coefficient $\eta_{\rm fid}$ in equation (2) is essentially determined by the degree of gas ionisation, which significantly depends on the dust properties. For instance, the ionisation degree increases and the dissipation coefficient decreases with dust growth because the dust surface area decreases and charged particles do not efficiently stick to the dust grains. Thus, to correctly determine the magnetic dissipation coefficient, it is necessary to solve the dust growth and chemical reactions of charged particles in the MHD calculations. However, present computational resources are insufficient for this task.

To investigate the effect of magnetic dissipation on outflow driving, we vary magnetic diffusivity as $\eta=C_{\eta}\,\eta_{\rm fid}$, ranging parameter $C_{\eta}$ from 0.01 to 100 (models C, CE1, CE2, CE3, CE4) (see Table 1 and §2.1). Among the models, $C_{\rm\eta}$ varies by a factor of $10^{4}$. Because the magnetic Reynolds number $Re_{\rm m}$ is roughly proportional to $\rho^{-1}$ (see Fig. 1 and the corresponding description in Machida et al. 2007), the difference in the density of the magnetic dissipation region between models CE1 ($C_{\eta}=0.01$) and CE4 ($C_{\eta}=100$) is about $10^{4}$. Simply assuming constant temperature and the Jeans wavelength as the typical spatial scale, the spatial scales of the magnetic dissipation regions differ by approximately $10^{2}$. Note that in reality, the temperature varies and estimating the differences in the magnetic dissipation regions is non-trivial (Machida et al., 2007; Dapp & Basu, 2010).

We then calculate the cloud evolution by fixing the parameters to those of model C but varying the diffusivity $\eta$ (i.e. $C_{\eta}$). The outflow momentum for each model is plotted against the protostellar mass in the left panel of Figure 12. In this panel, the outflow momenta are larger in models CE3 and CE4 than in models C, CE1 and CE2 for M_{\rm ps}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}3\thinspace M_{\odot}. As models CE3 and CE4 have larger diffusivities, the magnetic field of the first core, which forms before protostar formation, is weak. Thus, in these models, angular momentum is not effectively transported by magnetic effects and the first core sustains over a long duration with high rotational energy (Saigo & Tomisaka, 2006; Saigo et al., 2008). The outflow, which appears just after the first core formation (Machida et al., 2004; Tomisaka, 2002; Tomida et al., 2013, 2015), travels a large distance during the long lifetime of the first core and possesses a large momentum prior to protostar formation. The right panel of Figure 12 plots the outflow momentum against the elapsed time following the beginning of the cloud collapse. As the outflow momenta are very similar among the models, we can attribute the differences in the left panel of Figure 12 to the first core lifetimes. In summary, the diffusivity affects the first core lifetime, but does not significantly influence the outflow momentum.

Figure 13 shows the density and velocity distributions when the protostellar mass reaches $M_{\rm ps}=3\thinspace M_{\odot}$ in models CE1 (left panels) and CE4 (right panels). It is seen that diffusivity affects the disc size: a larger diffusivity effectively dissipates the magnetic field in the disc, weakening the angular momentum transport by magnetic effects (e.g. Masson et al., 2016; Wurster et al., 2016; Tomida et al., 2015; Tsukamoto et al., 2015). Consequently, a large disc forms. As the outflow is driven by the rotationally supported disc, its foot is located near the disc surface; therefore, the disc size influences the outflow driving region (see lower panels of Fig. 13).

Although diffusivity affects magnetic field lines at the small scales around the disc region, if the initial magnetic fields are strong, diffusivity has little impact on the field lines at large scales around the cloud core and the outflow. Instead, the structure of the magnetic field lines at the cloud core scale is primarily determined by the collapse geometry. Note that, because the amplification of the magnetic field in a small scale determines the outflow properties, the outflow may not be determined by the global collapse geometry when the magnetic field of the prestellar cloud is weak (Tomisaka, 2002). In such cases, the degree of magnetic diffusivity may affect large-scale structure of the outflow (see §4.2 and A).

Whereas the outflow at large scales significantly contributes to the total outflow momentum, the outflow momentum at small scales has a negligible contribution. Therefore, because the outflow driving region is located at small scales near the disc region, the difference in diffusivities among the models has a negligible impact on the outflow momentum (Fig. 12) when the initial magnetic field is strong. In other words, the outflow momentum is determined by the opening angle or geometry of the large-scale magnetic field lines, which is independent of magnetic diffusivity.

4.2 Dependence of Magnetic Field Strength

The fiducial models A–F (Table 1) are simulated under a strong magnetic field with a mass-to-flux ratio of two. In observations of prestellar cores, the magnetic and gravitational energies are comparable and the most plausible mass-to-flux ratio is $\mu=2-3$ (Crutcher, 1999). Recent observations by Tan et al. (2013) imply that massive starless cores are primary supported by the strong magnetic fields that control the subsequent core accretion phase. Although clarifying the origin of the strong magnetic field in prestellar cores is beyond the scope of this study, relating the magnetic field strength to the outflow properties may be important. To this end, we prepared four models with weak magnetic fields, namely CW1 and EW1 (with $\mu=5$) and CW2 and EW2 (with $\mu=10$). CW and EW have the same cloud parameters as their C and E counterparts, differing only in their magnetic field strengths (Table 1). The cloud evolution of these models is conducted identically to that in the fiducial models.

Figure 14 shows the outflow momenta in models C (fiducial model), CW1, CW2, E (fiducial model), EW1 and EW2 versus protostellar mass. The outflow momenta in the weak (CW1, CW2, EW1 and EW2) and strong (C and E) magnetic field models are not significantly different up to $M_{\rm ps}=1-2\thinspace M_{\odot}$, but noticeable differences appear beyond this mass. In the weak magnetic field models, the outflow physical quantities (mass, momentum, kinetic energy and momentum flux) are smaller than those in the strong magnetic field models after protostar formation, especially in the range M_{\rm ps}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{>}}}1-2\thinspace M_{\odot}. When the protostellar mass exceeds $\sim 3\thinspace M_{\odot}$, the outflow momenta among the models differs by over one order of magnitude. In the weak magnetic field models, the outflow finally disappears. Further calculation in the weak magnetic field models leads to recurrent fragmentation.

Figure 15 shows the time evolution of the weak magnetic field model EW2 ($\mu=10$). The outflow appears after protostar formation (Fig. 15a) and reaches $\sim 300$ AU from the central protostar (Fig. 15b). Thereafter, the outflow weakens and shrinks (Fig. 15c). Eventually, the outflow component disappears and a torus or toroid-like structure, supported mainly by the magnetic field, appears. Conversely, in the strong magnetic field model E, the outflow ultimately reaches $\sim 5000$ AU. Thus, changing the initial magnetic field strength qualitatively changes the results.

Our results imply that a magnetic energy comparable to the gravitational energy may be necessary to reproduce the observed outflow properties around a massive protostar. However, it should be noted that the necessary magnetic field strength for outflow driving is still controversial. In the ideal MHD approximation, Hennebelle et al. (2011) showed that outflow appears even when the initial cloud has an extremely weak magnetic field of $\mu=120$. On the other hand, recent studies have shown that an outflow does not appear in clouds with $\mu=20$ under ideal MHD calculations (Lewis & Bate, 2017) and $\mu=10$ under non-ideal MHD calculations (Wurster et al., 2016). Extending our calculations using the parameter of mass-to-flux ratio may produce precisely the magnetic field strength required for the observed massive outflows, but doing so is time-prohibitive; currently, a single run requires two to four months of wall clock time. The dependence of magnetic field strength on outflow driving will be the focus of a future study. We discuss non-ideal MHD effects on outflow driving in Appendix §A.

4.3 Comparison with Previous Studies

Several theoretical studies of massive star formation have simulated protostellar MHD outflow in the framework of the core accretion scenario (Hennebelle et al., 2011; Commerçon et al., 2011; Seifried et al., 2011, 2012; Peters et al., 2011, 2012; Myers et al., 2013). In this study, we evolved the cloud over a longer duration in resistive MHD simulations, ignoring the radiative effects. In this subsection, we compare our current work with the previous literature.

As the initial state, Hennebelle et al. (2011) assumed cores with 100 solar masses and $\alpha_{0}=0.12$. They ran ideal MHD equations, also without radiation effects, and computed the outflow evolution in collapsing clouds for $\sim 10^{4}-10^{5}$ yr at spatial resolutions of 2 and 8 AU. They achieved an outflow mass of $\sim\thinspace M_{\odot}$ at the end of their calculation. Their initial cloud conditions correspond to our models B ($\alpha_{0}=0.2$) or C ($\alpha_{0}=0.08$), which yielded an outflow mass of \mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}10\thinspace M_{\odot} (Figs. 10a and d). As Hennebelle et al. (2011) did not impose a sink, their protostellar mass was not uniquely determined. However, their results closely agree with our simulations at relatively small accretion rates.

Under similar conditions to Hennebelle et al. (2011), Commerçon et al. (2011) investigated the evolution of a massive star forming cloud with radiation effects and showed that fragmentation is suppressed in an initially strongly magnetised cloud, as was also observed in our study. In their calculation, the outflow reached $\sim 1000$ AU with a mass accretion rate of $\sim 10^{-4}M_{\odot}\,{\rm yr}^{-1}$. Quantitatively, this result is consistent with our calculations. However, their simulation results are difficult to compare against observations because they computed only the earliest phase of massive star formation.

Seifried et al. (2011) and Seifried et al. (2012) also investigated massive star formation under the core accretion scenario. Their sink radius ($r_{\rm sink}=12$ AU) and spatial resolution (4.7 AU) were quite large and neither study appears to have resolved the disc forming region (Machida et al., 2014). They included radiative effects but ignored the magnetic dissipation process. As the disc forming or magnetic dissipation region determines the outflow properties, it cannot be ascertained whether they correctly (or quantitatively) estimated these properties in their calculation. They evolved the cloud to a protostellar mass of $\sim 1-4\thinspace M_{\odot}$, where the mass and momentum of the outflow were 0.1\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}M_{\rm out}/(\thinspace M_{\odot})\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}1.1 and $0.2\leq P_{\rm out}/(\thinspace M_{\odot}\,{\rm km}\,{\rm s}^{-1})\leq 2$, respectively. Both quantities are substantially smaller than our results and observations, which may be attributed to their lower spatial resolution and larger sink radius.

Peters et al. (2011) investigated the collapse of massive clouds in a radiation-MHD simulation. Although they also reproduced the outflow in massive star forming regions, the scale of their region was very large, with a sink radius of 590 AU. As the outflow driving region was completely embedded in their sink, we do not quantitatively compare their results with ours. However, the suppression of fragmentation by the magnetic effects was confirmed in both studies.

Our study serves as an advancement on previous studies in two major aspect:

• •

We managed to compare the simulated outflow properties with observations using long-term time integration simulations that resolved the regions of magnetic dissipation and outflow driving.

• •

In addition, by varying the mass accretion rates or initial cloud stabilities, we could relate massive star formation dynamics to those of low (Machida & Hosokawa, 2013) and intermediate star formation.

As described above, this study focused not merely on massive star formation but also covered a wide mass range. Moreover, the outflow properties simulated at higher mass accretion rates agreed with those of previous studies. Therefore, the results of this study might provide important clues to the better understanding of the massive star formation process.

4.4 Effects of Radiative Feedback

As described in §2, this study ignored radiative feedback from the protostar, because our goal was simply to relate mass accretion (inflow) to mass outflow driven by magnetic effects. However, we compared our results with observations of massive outflows from massive protostars that are very luminous in the late accretion phase. Based on this comparison, in this subsection, we discuss the effects of radiative feedback from massive protostars.

In recent years, massive star formation has been studied using radiation-hydrodynamic simulations. Krumholz et al. (2007) investigated massive star formation through the use of radiation hydrodynamic simulation in which they primary focused on disc fragmentation. They showed that fragmentation is suppressed around massive protostars because such stars heat their surroundings (Toomre, 1964). Although the magnetic field can also suppress fragmentation, they ignored magnetic effects. Thus, Krumholz et al. (2007) produced results that are qualitatively close to those of this study, in which fragmentation is suppressed by either radiative or magnetic effects and only a one or a small number of protostars tend to appear around the centres of massive collapsing clouds. In our calculations, fragmentation occurs only in model F, which has the highest accretion rate and the most massive initial prestellar cloud core. When the radiative effects are included, the disc heating may further suppress fragmentation, even for model F.

In addition to fragmentation, the disc heating is also expected to influence the magnetic dissipation process. The ionisation rate increases with the disc temperature, causing magnetic dissipation to become ineffective (e.g. Machida et al., 2007). In such cases, a more powerful magnetic wind may be driven from the circumstellar disc (Tomisaka, 2002). On the other hand, a strong magnetic field may suppress disc growth through magnetic braking even when the mass accretion rate is significantly high. When the circumstellar disc does not grow sufficiently, a powerful magnetic outflow may not appear (Li et al., 2011). Therefore, radiative feedback (or disc heating) is expected to have conflicting effects on magnetically driven outflow. Thus, further investigation into the massive star formation process in the context of the magnetic field and its dissipation process and radiation feedback–especially after the protostar has grown sufficiently—will be necessary to develop a better quantitative comparison between simulation results and observations (see below).

We next discuss the effect of radiative feedback on mass accretion. Mass accretion onto the circumstellar disc is necessary in driving the magnetic outflow because the outflow is powered by the release of the gravitational energy of accreting matter. Radiative feedback from massive protostars is often discussed with regard to whether the mass accretion lasts until the formation of a very massive star (\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{>}}}50-100\thinspace M_{\odot}; Wolfire & Cassinelli, 1987). However, as very massive star formation is beyond the scope of this study, we discuss the protostellar mass at which mass accretion and circumstellar environment are affected by radiative feedback. Krumholz et al. (2007) showed that the radiation from the protostar does not significantly interrupt mass accretion until the protostellar mass exceeds $\sim 20\thinspace M_{\odot}$. Kuiper et al. (2010, 2015) also showed that radiative feedback begins to influence the mass accretion and the circumstellar region when the protostellar mass exceeds \mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{>}}}15-20\thinspace M_{\odot} (Kuiper et al., 2016). Although the final stellar mass at which the mass accretion is stopped by radiative feedback is still controversial (Kuiper et al., 2010), previous studies indicate that radiative feedback does not significantly affect mass accretion until the protostellar mass reaches $\sim 10-20\thinspace M_{\odot}$. Kuiper et al. (2016) also showed that, after the protostellar mass exceeds \mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{>}}}20-30\thinspace M_{\odot}, the stellar radiation force is comparable to gravity above and below the disc and strongly affects the circumstellar environment (see also Yorke & Sonnhalter, 2002).

The comparison of our simulation results with observations in §3.5 was limited to the range M_{\rm ps}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}20-30\thinspace M_{\odot}, in which we can conclude that radiative feedback does not significantly affect magnetically driven outflow. On the other hand, the effect of radiative feedback on protostellar outflow cannot be ignored for M_{\rm ps}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{>}}}20-30\thinspace M_{\odot}, because the radiative pressure also contributes to the outflow (radiation-driven outflow) from the central region (Kuiper et al., 2016).

5 Summary

In this study, we investigated the relation between mass accretion onto protostars (or circumstellar disks) and the magnetically driven outflow in gravitationally collapsing cloud cores. As an initial state, we assumed a gravitationally unstable prestellar cloud core with spherical symmetry. We found that the mass accretion rate onto the central region of a collapsing cloud core is governed by the ratio of the thermal to gravitational energy in the initial cloud, which can be parametrised as $\alpha_{0}$. To investigate the effects of the mass accretion rate on outflow driving, we varied $\alpha_{0}$ and calculated the cloud evolution and outflow propagation.

When the initial prestellar cloud has a strong magnetic field, the outflow is powerful at any mass accretion rate. The outflow is driven by the outer disc region, where the magnetic field and neutral gas are well coupled. In our simulations, the mass outflow rate was found to be proportional to the mass accretion rate. In addition, the ratio of the outflow to the accretion rate does not significantly depend on the mass inflow rate or the initial cloud stability.

Approximately 20–40% of the infalling gas is ejected by the protostellar outflow from the circumstellar region. Regardless of the accretion rate, the outflow has a very wide opening angle at the disc scale and a well-collimated structure at large scales. When the prestellar cloud has a strong magnetic field, the outflow structure is determined by the large-scale geometry of the magnetic field lines. The outflow travels into a dense infalling envelope and sweeps part of the infalling gas through the wide opening angle. On the cloud core scale, the mass ratio of the outflowing gas (including the gas swept up by the outflow) to the infalling gas exceeds 50%. Thus, over half of the infalling gas is finally ejected into interstellar space by the outflow.

To validate the importance of magnetically driven outflow, we quantitatively compared our simulated outflow mass, momentum, kinetic energy and momentum flux with their observed counterparts. The physical quantities derived from our simulations closely agree with those observed in massive outflows (Fig. 10). The simulated momentum fluxes are roughly distributed between the maximum and minimum observed momentum fluxes (Fig. 11), suggesting that the massive outflows from massive protostars are magnetically driven, as true in low-mass star formation. On the other hand, the outflow is subdued or absent in cloud cores with an initially weak magnetic field. In addition, weakly magnetised clouds fragment multiple times, preventing the formation of a single massive star. Because a strong magnetic field sufficiently explains the massive outflow and suppressed vigorous fragmentation, it may be a requirement for massive star formation.

By combining the results of this study with those of Machida & Hosokawa (2013), we can state that both low- and high-mass stars form through a common fundamental mechanism (i.e. the core accretion scenario). Once a protostar has formed in a gravitationally collapsing cloud, a circumstellar disc grows and establishes a surrounding outflow powered by gas accretion or the gravitational energy released by infalling gas; consequently, the outflow momentum and energy increase with the accretion rate. In this manner, the accretion rate determines both the physical quantities of the outflow and the resultant stellar mass. Thus, the massive outflows observed around massive protostars are evidence of these star’s evolution according to the core accretion scenario. This study only investigated the relation between accretion and outflow rate, and the comparisons between our simulated and previously observed massive outflows are rather tentative. In particular, we ignored turbulence, radiation effects and the surrounding environment outside the star forming cloud, all of which will affect massive star formation. These issues were not discussed here in detail because doing so would involve broadening our investigation beyond massive star formation. However, to more fully understand the formation processes of low- to high-mass stars, it will ultimately be necessary to consider these effects.

Acknowledgements

This study has benefited greatly from discussions with  T. Hirota,  K. Tomida,  K. Motogi,  M. Sekiya, and  T. Tsuribe. We are very grateful to an anonymous reviewer for a number of useful suggestions and comments. The authors would like to thank Enago (www.enago.jp) for the English language review. This research used computational resources of the HPCI system provided by (Cyber Sciencecenter, Tohoku University; Cybermedia Center, Osaka University, Earth Simulator, JAMSTEC) through the HPCI System Research Project (Project ID:hp150092, hp160079). This work was supported by Grants-in-Aid from MEXT (25400232, 26103707).

Appendix A Difference in Outflow Driving between Ideal and non-ideal MHD models

In §4.2, we showed that the outflow driven by the circumstellar disc weakens and disappears quickly in some weakly magnetised models; by contrast, it seems to persist for a long time in previous numerical studies on massive star formation (Hennebelle et al., 2011; Seifried et al., 2012). There are many physical and numerical differences between this and previous studies, including the initial condition, the sink treatment and the spatial resolution. We contend that the most important difference is whether magnetic dissipation is taken into account. Whereas Hennebelle et al. (2011) and Seifried et al. (2012) calculated cloud evolution using ideal MHD equations, we calculated it using non-ideal MHD equations.

To investigate the effect of the inclusion of magnetic dissipation on outflow driving, we re-calculated the cloud evolution of model EW2 using ideal MHD equations (hereafter we will call the amended model EW2I). The numerical settings were identical to those of the models listed in Table 1 except for the inclusion (in EW2, the non-ideal MHD model) or exclusion (in EW2I, the ideal MHD model) of the magnetic field dissipation term. As seen from Figure 15, in the non-ideal MHD model EW2, after the outflow reaches $\sim 400$ AU, it gradually weakens and finally disappears within $\sim 8000$ yr after protostar formation. On the other hand, in the ideal MHD model EW2I (Fig. 16), outflow reaches $\sim 5000$ AU and does not weaken within $\sim 8000$ yr after protostar formation. Thus, magnetic dissipation causes a qualitative difference in the outflow driving process.

The reason for the weakening and disappearance of outflows in the magnetic dissipation models is simple: because the neutral gas is not well coupled to the magnetic field in a magnetically inactive region in which the magnetic field dissipates, the outflow, which is mainly driven from the circumstellar disc by the Lorentz force, does not appear in such regions (Machida et al., 2007). Ohmic resistivity depends on the density (and gas temperature) and a high-density region in the range n\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}10^{15}\,{\rm cm}^{-3} has a lower ionisation rate or a larger resistivity (e.g. Nakano et al., 2002). In the massive star formation process, the disc surface density becomes high when the initial magnetic field strength is weak because magnetic braking is not effective even in the magnetically active region where the neutral gas is coupled to the magnetic field. Thus, the magnetically inactive region is expected to expand with time and the entire region of the disc might become magnetically inactive. In such cases, the outflow does not appear because the magnetic field lines spin free and the magnetic field, which is not coupled to the neutral gas, cannot drive the flow.

To compare the size of the disc and the magnetically inactive region, Figure 17 plots the density and velocity distributions at density contours of $n=10^{12}$ and $10^{13}\,{\rm cm}^{-3}$ on the $y=0$ plane for models E (left), EW2 (middle) and EW2I (right) at a protostellar mass of $\sim 4\thinspace M_{\odot}$. In the high density region $n>10^{12}\,{\rm cm}^{-3}$, which corresponds to the region inside the dotted black contour in Figure 17, the magnetic field dissipates through ohmic dissipation. For model E (Fig. 17 left), the outflow is driven by the outer edge of the rotationally supported disc, where the magnetic field is coupled to the neutral gas. On the other hand, in model EW2 (Fig. 17 middle) the disc and its surrounding envelope have a high density of \mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{>}}}10^{12}\,{\rm cm}^{-3}, and the overall disc region becomes magnetically inactive. Although the outflow is just barely driven at this point, it gradually weakens and finally disappears as shown in Figure 15 because the magnetic field dissipates and gradually weakens in the overall disc region. Model EW2I (and EW) has an initially weak magnetic field and a high mass accretion rate. The former weakens the magnetic braking while the latter promotes the formation and growth of a massive disc having a high surface density. Correspondingly, model EW2 produces a large magnetically inactive region in which outflow is not efficiently driven. Even when the magnetic field is weak, the outflow continues to be driven in the ideal MHD model, as shown in Figures 16 and 17 (right panel) because the magnetic field continues to be amplified by the disc rotation in the ideal MHD approximation.

Finally, to confirm outflow driving in an initially extremely weak magnetised cloud, we prepared the same initial condition and adopted the same numerical settings as in Hennebelle et al. (2011) and calculated the cloud evolution using ideal MHD equations and without a sink (hereafter, we call this model HI). The initial cloud has a density distribution of

[TABLE]

where $\rho_{c}=1.4\times 10^{-20}\,{\rm g}\,{\rm cm}^{-3}$ and $r_{0}=0.22$ pc. The initial cloud has a mass of $100\thinspace M_{\odot}$ and a radius of $1.35$ pc. A uniform magnetic field ($B_{0}=3\times 10^{-7}$ G) is adopted and adjusted to produce a mass-to-flux ratio $\mu=120$ in the overall cloud. Whereas an internal turbulent velocity dispersion is given as the source of cloud rotation in Hennebelle et al. (2011), a rigid rotation of $\Omega_{0}=10^{-14}$ $\rm{s}^{-1}$ is adopted for model HI.

Figure 18 shows the density and velocity distributions on the $y=0$ plane at $t=5.8\times 10^{5}$ yr after the cloud begins to collapse. Although the initial cloud has an extremely weak magnetic field, the outflow continues and does not weaken by the end of the calculation (see also §4.3). This indicates that, in the ideal MHD approximation, the outflow can be driven by the circumstellar disc even when the initial magnetic field is extremely weak. Using an ideal MHD simulation, Tomisaka (2002) also showed that magnetic field lines are highly twisted by the rotational motion of the central region and that the amplified magnetic field drives outflow. As we reproduced a long-lived outflow in a cloud with an extremely weak magnetic field using ideal MHD calculations, our result does not contradict previous studies using ideal MHD calculations (Tomisaka, 2002; Hennebelle et al., 2011; Seifried et al., 2012). On the other hand, in the non-ideal MHD model, since the magnetic field dissipates when the density of the disc and the surrounding medium is high, the magnetic field lines spin free and the magnetic field is not amplified. Therefore, even when outflow appears just after protostar or disc formation, it can occasionally weaken or disappear for short periods of time, as shown in Figure 15.

Observations of star forming clouds indicate that their magnetic fields are strong with magnetic energies comparable to their gravitational energies (Crutcher et al., 2010). Thus, the assumption of a star forming cloud with an extremely weak magnetic field is not very realistic. This study focused only on the relation between the mass accretion and ejection rates, and clarifying the outflow driving condition in high-mass star formation is beyond its scope. In addition, much more computationally expensive calculation would be required to accurately determine the outflow driving conditions at different magnetic field strengths. However, as it is important to comprehensively understand the outflow driving condition in the high-mass star formation process, we will focus on this issue in future work.