Formation of Unipolar Outflow and Protostellar Rocket Effect in Magnetized Turbulent Molecular Cloud Cores

Observed protostellar outflows exhibit a variety of asymmetrical features, including remarkable unipolar outflows and bending outflows. Revealing the formation and early evolution of such asymmetrical protostellar outflows, especially the unipolar outflows, is essential for a better understanding of the star and planet formation because they can dramatically change the mass accretion and angular momentum transport to the protostars and protoplanetary disks. Here we perform three-dimensional nonideal magnetohydrodynamics simulations to investigate the formation and early evolution of the asymmetrical protostellar outflows in magnetized turbulent isolated molecular cloud cores. We find, for the first time to our knowledge, that the unipolar outflow forms even in the single low-mass protostellar system. The results show that the unipolar outflow is driven in the weakly magnetized cloud cores with the dimensionless mass-to-flux ratios of μ = 8 and 16. Furthermore, we find the protostellar rocket effect of the unipolar outflow, which is similar to the launch and propulsion of a rocket. The unipolar outflow ejects the protostellar system from the central dense region to the outer region of the parent cloud core, and the ram pressure caused by its ejection suppresses the driving of additional new outflows. In contrast, the bending bipolar outflow is driven in the moderately magnetized cloud core with μ = 4. The ratio of the magnetic to turbulent energies of a parent cloud core may play a key role in the formation of asymmetrical protostellar outflows.


INTRODUCTION
Protostellar outflows play essential roles in the formation process of protostars and planets.The protostellar outflows are one of the observable signs of the birth of protostars in the dense regions of molecular cloud cores.The mass accretion and angular momentum transport to protostars and protoplanetary disks are regulated by the driving and subsequent evolution of the protostellar outflows, contributing to the determination of the star formation efficiency in the isolated protostellar cores (Machida & Hosokawa 2013).Moreover, the grown dust with a size of about a centimeter in the inner region of the protoplanetary disk is entrained by the protostellar outflow and refluxed from the outflow onto the outer region of several tens of astronomical units of the disk (Tsukamoto et al. 2021b).The above process, which is referred to as the ashfall phenomenon, can circumvent the radial drift barrier that typically hinders planet formation.Therefore, the protostellar outflows contribute to crucial physical processes for the formation of protostars and planets.
One of the remarkable discoveries is that the molecular outflows driven by YSOs exhibit a variety of asymmetrical features.For instance, a statistical study of the properties of molecular outflows conducted by Wu et al. (2004) show that 50 sources (13 %) out of 391 targets present unipolar outflows, of which 28 are redshifted ones, indicating that the redshifted unipolar outflows are equally abundant compared to the blueshifted ones and the unipolar outflows are intrinsic.Recent highresolution observations of the protostellar environment in the Orion A molecular cloud using the Atacama Large Millimeter/submillimeter Array (ALMA) have more clearly detected the isolated protostars driving the unipolar outflow (Hsieh et al. 2023).The resolved unipolar cavities in the circumstellar envelopes are delineated by the near-infrared images of the scattered light around protostars, which typically trace the hot shocked gas regions in the outflows, obtained by the survey of protostellar outflow cavities in the Orion molecular clouds using the Hubble Space Telescope (Habel et al. 2021).Although the optical and near-infrared observations are inherently affected by the extinction effect, the result can provide complementary evidence for the existence of the unipolar outflows.Additionally, many detections of unipolar outflows in low-mass and high-mass YSOs are reported from recent ALMA observations (e.g., Aso et al. 2018;Louvet et al. 2018;Kong et al. 2019;Aso et al. 2019;de Valon et al. 2020;Garufi et al. 2020;Li et al. 2020;Dutta et al. 2020;Baug et al. 2021;Sato et al. 2023a;Dutta et al. 2023;Sato et al. 2023b).Moreover, observations have frequently detected bipolar molecular outflows with asymmetrical features, such as the bendings and different sizes of redshifted and blueshifted lobes (Arce et al. 2013;Yen et al. 2015Yen et al. , 2017a;;Aso et al. 2018Aso et al. , 2019;;Okoda et al. 2021;Hsieh et al. 2023;Kido et al. 2023).These observational results suggest that the asymmetrical features frequently emerge in the outflows.
Although the unipolar outflows have been observed, the formation and early evolution of the unipolar outflows driven by low-mass protostars remain unclear so far.Understanding the formation and early evolution of the outflows with the outstanding asymmetrical features, in particular the unipolar outflows, is crucial because they can strongly regulate the mass accretion and angular momentum transport to the protostars and protoplanetary disks.Observed molecular cloud cores are often threaded by the magnetic field (Crutcher 2012;Pattle et al. 2023), and in the context of low-mass star formation, the bipolar outflows magnetically driven by the first cores and the protoplanetary disks are successfully reproduced by many numerical simulations of the gravitational collapse of rotating magnetized cloud cores (e.g., Tomisaka 1998Tomisaka , 2002;;Matsumoto & Tomisaka 2004;Banerjee & Pudritz 2006;Machida et al. 2008;Tomida et al. 2010Tomida et al. , 2013;;Bate et al. 2014;Tomida et al. 2015;Tsukamoto et al. 2015aTsukamoto et al. , 2017;;Wurster et al. 2018a;Vaytet et al. 2018;Tsukamoto et al. 2018Tsukamoto et al. , 2020;;Hirano et al. 2020;Tsukamoto et al. 2021b;Wurster et al. 2021;Marchand et al. 2023).On the other hand, the isolated low-mass cloud cores are often associated with the turbulent velocity field as well (e.g., Goodman et al. 1993;Barranco & Goodman 1998;Burkert & Bodenheimer 2000;Misugi et al. 2019Misugi et al. , 2023)), which is typically subsonic or at most transonic (Ward-Thompson et al. 2007).The turbulence in the low-mass cloud cores can naturally produce the complex asymmetrical structures such as the warped protoplanetary disks and the filamentary infalling envelopes, and the rotation directions of the disks change dynamically as time proceeds due to the chaotic accretions (e.g., Tsukamoto & Machida 2013;Takaishi et al. 2020Takaishi et al. , 2021)).Thus, the turbulence is expected to generate the asymmetry of the magnetic field, leading to the unipolar outflows.However, few studies have performed that the asymmetrical bipolar outflows form in the realistic simulations of the collapse of magnetized turbulent low-mass cloud cores with/without coherent rotation (Matsumoto & Hanawa 2011;Joos et al. 2013;Matsumoto et al. 2017;Lewis & Bate 2018), and they have not yet identified the unipolar outflows driven by low-mass protostars.
In contrast to the above simulations, in a recent paper, Mignon-Risse et al. (2021a) report that the weak and transient unipolar outflow is driven by evolving binary massive protostars formed in the most turbulent case of their simulations, which calculate the gravitational collapse of magnetized turbulent massive cores of 100 M ⊙ by solving the magnetohydrodynamics (MHD) equations including the ambipolar diffusion and hybrid radiative transfer in the context of high-mass star formation (see also Mignon-Risse et al. (2021b) for details of their simulations).The results indicate that the initial turbulence of the parent cloud core, i.e. the environmental ram pressure, can strongly affect the outflow driving.Machida & Hosokawa (2020) also show that the ram pressure caused by the infalling envelope with the high accretion rate suppresses the outflow driving when the rigidly rotating cloud core initially has a weaker magnetic field.These results suggest that the ram pressure caused by the stronger turbulence and weaker magnetic field may be crucial for the formation of the unipolar outflows in the low-mass protostar formation as well as in the high-mass protostar formation.
This paper reports, for the first time to our knowledge, the formation of the unipolar outflows driven by the single low-mass protostellar systems.We perform the simulations of the gravitational collapse of magnetized turbulent low-mass cloud cores of 1 M ⊙ with strong, moderate and weak magnetic fields to investigate the formation and early evolution of the unipolar outflows, which have not yet been explored in the previous studies.In addition, this paper presents the subsequent evolution of the protostellar system driving the unipolar outflow, which is similar to the launch and propulsion of a rocket.
The rest of the paper consists of the following sections.Section 2 describes the numerical method and the initial conditions.Section 3 presents the numerical results.Finally, Section 4 summarizes and discusses the results and findings.

Basic Equations and Numerical Method
The simulations solve the non-ideal MHD equations including self-gravity of the gas: where v is the gas velocity, ρ is the gas density, P is the gas pressure, B is the magnetic field, ϕ is the gravitational potential, and G is the gravitational constant.The unit vector along the magnetic field is denoted by B ≡ B/|B|.η O and η A are the resistivities for the Ohmic dissipation and ambipolar diffusion of the non-ideal MHD effects.We note that the Hall effect, which is another important factor of the non-ideal MHD effects, is currently ignored because the simulation incorporating the Hall effect requires highly small time steps in its calculation and thus it is computationally highly demanding for the long-term simulation up to ∼ 10 4 yr after the protostar formation.Instead of solving the radiation transfer, we adopt the barotropic equation of state that mimics the thermal evolution of the cloud core presented by the radiation hydrodynamics simulations (Masunaga & Inutsuka 2000): where c s,iso = 1.9 × 10 4 cm s −1 is the isothermal sound speed at the temperature of 10 K, and ρ crit = 4 × 10 −14 g cm −3 is the critical density at which the thermal evolution changes from the isothermal to adiabatic.The equations are solved with the smoothed particle hydrodynamics (SPH) method (Lucy 1977;Gingold & Monaghan 1977;Monaghan & Lattanzio 1985).The ideal MHD part of the equations is solved with the Godunov smoothed particle magnetohydrodynamics (GSPMHD) method (Iwasaki & Inutsuka 2011).In addition, we adopt the hyperbolic divergence cleaning method proposed by Iwasaki & Inutsuka (2013) so that the divergence-free condition of the magnetic field is satisfied.The Ohmic dissipation and ambipolar diffusion are calculated with the methods of Tsukamoto et al. (2013a) and Wurster et al. (2014).The processes of the Ohmic dissipation and ambipolar diffusion are accelerated by super-time-stepping (STS) method (Alexiades et al. 1996).The parameters of ν sts = 0.01 and N sts = 5 are adopted for STS (Tsukamoto et al. 2013a).The self-gravity of the gas is computed by the Barnes-Hut octree algorithm with an opening angle parameter of θ gravity = 0.5 (Barnes & Hut 1986).The spline interpolation for the gravitational softening is adopted with the technique of the adaptive softening length (Price & Monaghan 2007).The numerical code is parallelized with the Message Passing Interface (MPI).The numerical code has already been applied to a variety of problems (e.g., Tsukamoto & Machida 2011, 2013;Tsukamoto et al. 2013aTsukamoto et al. ,b, 2015cTsukamoto et al. ,b,a, 2017Tsukamoto et al. , 2018;;Takaishi et al. 2020;Tsukamoto et al. 2020Tsukamoto et al. , 2021a,b;,b;Takaishi et al. 2021;Tsukamoto et al. 2023a).

Resistivity model
The simulations use the tabulated resistivities for η O and η A that are presented as the single-sized dust model with the dust size of a d = 0.035 µm in Tsukamoto et al. (2020).The resistivities are generated by the chemical reaction network calculation using the method of Susa et al. (2015) The calculation also includes the neutral and singly charged dust grains, G 0 , G + , and G − .The calculation takes into account the cosmic-ray ionization, gas-phase and dust-surface recombination, and ion-neutral reactions.The indirect ionization by high-energy photons emitted by direct cosmicray ionization (described as CRPHOT in the UMIST database) is also considered in the calculation.The initial abundance and reaction rates are taken from the UMIST2012 database (McElroy et al. 2013).The grain-ion and grain-grain collision rates are calculated using the equations in Draine & Sutin (1987).The chemical reaction network calculation is conducted using the CVODE package (Hindmarsh et al. 2005) assuming the system is in the chemical equilibrium.The resistivities are calculated using the abundances of charged species in the equilibrium state.The momentum transfer rate between the charged and neutral species is calculated using the equations in Pinto & Galli (2008).The temperature for the chemical reaction network calculation is modeled as T chem = 10(1 + γ T (ρ/ρ crit ) (γ T −1) ) K, where γ T = 7/5.It is assumed that the dust has the internal density of ρ d = 2 g cm −3 and the size of a d = 0.035 µm.The dust-to-gas mass ratio is fixed to be 0.01 in the calculation.The cosmic ray ionization rate is assumed to be ξ CR = 10 −17 s −1 .

Initial conditions for the density profile
The simulations start from the collapse of isolated molecular cloud cores including the magnetic field and turbulence simultaneously.The initial cloud core follows the density profile presented by Tsukamoto et al. (2020): where f is the density enhancement factor, ρ c is the characteristic density, and R c = 6.45a is the radius of the initial cloud core.ϱ BE is the non-dimensional density profile of the Bonnor-Ebert sphere (Bonnor 1956;Ebert 1955), which is a pressure-confined and self-gravitating isothermal gas sphere in hydrostatic equilibrium state against the gravitational collapse and well fits the density distribution of observed isolated molecular cloud cores (e.g., Alves et al. 2001;Kandori et al. 2005).
The density enhancement factor f controls the strength of gravity.More specifically, f can be written as f = 0.84/α, where with E thm and E grav are the thermal and gravitational energies of the initial cloud core (without the surrounding medium of ρ(r) ∝ r −4 ); see Appendix A of Matsumoto & Hanawa (2011).The density profile of the initial cloud core with f = 1 in a region of r < R c corresponds to that of the critical Bonnor-Ebert sphere.The initial cloud core with f > 1 is gravitationally unstable.
The initial cloud core has the temperature of T c = 10 K, the radius of R c = 4.8 × 10 3 au, the mass of M c = 1 M ⊙ within r < R c (∼ 2.1 M ⊙ in the entire domain of r < 10R c ), and the ratio of α = 0.4, which are determined by specifying f = 2.1 and the central density of ρ 0 = f ρ c = 7.3 × 10 −18 g cm −3 .The free-fall time is calculated as t ff = (3π/(32Gρ 0 )) 1/2 = 2.5 × 10 4 yr.The simulations resolve M c = 1 M ⊙ with the number of SPH particles of N SPH = 10 6 , which corresponds to the mass resolution of M c /N SPH = 10 −6 M ⊙ .The number of all the SPH particles in the entire domain of

Turbulence and magnetic fields
The initial cloud core has the divergence-free turbulent velocity field with the velocity power spectrum of P v (k) ∝ k −4 (Burkert & Bodenheimer 2000), where k is the wavenumber.The amplitude of the turbulence is characterized by the mean sonic Mach number: or the ratio of the turbulent to gravitational energies: where σ v and E turb = 3σ 2 v M c /2 are the one-dimensional velocity dispersion and turbulent energy of the initial cloud core without the surrounding medium of ρ(r) ∝ r −4 .We adopt M s = 0.86, which corresponds to γ turb = 0.1 initially.Hence, the initial cloud core has a non-vanishing net angular momentum of |J c,net | = 4.4 × 10 53 g cm 2 s −1 due to the stochastic nature of the turbulent velocity field.The direction of J c,net is set as the z-axis of the simulations.
The turbulent velocity field is generated with the method of Tsukamoto & Machida (2013), which is also adopted in our previous studies (Takaishi et al. 2020(Takaishi et al. , 2021)).All the simulations use the same turbulent velocity field, which is assigned to the cloud core of r < R c alone, and vanishes for R c ≤ r < 10R c .We note that the initial velocity field consists only of the turbulent velocity field, and it is not the superposition of the turbulent and rotational velocity fields.
The initial cloud core has the axisymmetric magnetic field modeled by Tsukamoto et al. (2020).In cylindrical coordinates (R, φ, z), they are where B c denotes the strength of the central magnetic field.The magnetic field described by the equations has a constant uniform component for z-direction in the central region ( B R → 0 and B z → B c as r → 0) and becomes an hourglass-shaped structure with |B| ∝ r −2 in the outer region (except at the midplane) as r → ∞ (see Appendix A of Tsukamoto et al. ( 2020) for details).The simulations can avoid to emerge low-β plasma regions in the surrounding medium of ρ(r) ∝ r −4 , where β plasma is the plasma beta parameter.Such the hourglass-shaped structures of the magnetic field have been estimated by previous observations on protostellar cores (e.g., Girart et al. 2006;Kandori et al. 2020a) and starless cores (e.g., Kandori et al. 2017Kandori et al. , 2018Kandori et al. , 2020b)).
The initial cloud cores are parameterized with the strength of the central magnetic field B c , which can be expressed by using the dimensionless mass-to-flux ratio: where Φ mag (R) is the magnetic flux of the initial cloud core calculated as and (M c /Φ mag ) crit = (0.53/3π)(5/G) 1/2 is the critical mass-to-flux ratio on stability for uniform spheres (Mouschovias & Spitzer 1976).
The simulations are conducted with four different dimensionless mass-to-flux ratios, µ = 2, 4, 8, and 16, which correspond to B c = 252 µG, 126 µG, 63 µG, and 31 µG, respectively.The dimensionless mass-to-flux ratio with the constant value of the central magnetic field B c is defined as As noted in Tsukamoto et al. (2020), µ const would be a suitable indicator to compare the strength of the magnetic field of this study with those of previous studies because we focus on the time evolution of the central region of the cloud core.µ can be converted into the ratio of the magnetic to gravitational energies: where E mag is the magnetic energy of the initial cloud core without the surrounding medium of ρ(r) ∝ r −4 .The model names and corresponding parameters are summarized in Table 1.We note that the last column of Table 1 summarizes the results of driven outflow morphologies.3. RESULTS

overview on 3D structure of the driven outflows and magnetic field
First, we show an overview of the simulations, which clearly show the very different morphologies of the driven outflows.Figure 1 shows the three-dimensional views of the density distribution and magnetic field structure for the three driven outflows: the no outflow (model MF2, top green box), bipolar outflow (model MF4, middle blue box), and unipolar outflow (model MF8, bottom magenta box).The simulations are conducted from the protostellar collapse to t p ∼ 10 4 yr, where t p denotes the elapsed time after the protostar formation epoch defined at the time when the central density becomes higher than 1.0 × 10 −11 g cm −3 .
The top green box of Figure 1 shows that no outflow is driven in the model MF2 in which the initial cloud core has a strong magnetic field of µ = 2.The hourglass-shaped structure of the magnetic field is formed and kept until the end of the simulation.The protoplanetary disk at the central region does not form in a timescale of t p ∼ 1.5 × 10 4 yr whereas the disk-shaped flattened infalling envelope, which is so-called the pseudo-disk structure (Galli & Shu 1993a,b), forms in this model.This result indicates that the relatively strong magnetic field rapidly extracts the angular momentum caused by the turbulent accretion from the central region via magnetic braking and suppresses the formation of a Keplerian rotating disk at the central region.
The middle blue box of Figure 1 shows that the bipolar outflow is driven in the model MF4 in which the initial cloud core has the magnetic field of µ = 4.The bipolar outflow is mainly driven by the magneto-centrifugal wind model (Blandford & Payne 1982) whereas the spiral-flow model (Matsumoto & Hanawa 2011;Matsumoto et al. 2017) also contributes to driving it.The magnetic field gradually becomes twisted as time proceeds, and the driven bipolar outflow evolves to be a large size of ∼ 5 × 10 2 au in the lower region and ∼ 10 3 au in the upper region, indicating that the outflow size of the upper region is greater than that of the lower region.Furthermore, the results show that the driven bipolar outflow is slightly bending and its driving directions in the upper and lower regions are not perfectly antiparallel with each other.The driving directions of the bipolar outflow are also misaligned with the direction of the large-scale global magnetic field (roughly the z-axis  Yellow isosurfaces show the density of 3.2×10 −18 g cm −3 , 10 −17 g cm −3 , 3.2 × 10 −17 g cm −3 , and 10 −16 g cm −3 with the radial velocity of v r > 0, representing the outflows.White lines show the magnetic field lines.Blue and green isosurfaces show the density of 3.2 × 10 −17 g cm −3 and 10 −16 g cm −3 with v r < 0, representing the infalling envelopes.Cut-plane densities on x − y (z = 0), x − z (y = 0), and y − z (x = 0) planes are projected for each panel.t p notes the elapsed time after the protostar formation epoch defined at the time when the central density becomes higher than 1.0 × 10 −11 g cm −3 .The scale of the box is ∼ 2, 000 au.The origin of a coordinate system is shifted to the center of mass of the system.direction).The warped structures of the infalling envelopes are emerged by the perturbation of the turbulent accretion.Our results suggest that the bending and misaligned structures of the bipolar outflow are created by the turbulent accretion via the surrounded infalling envelopes.Therefore, the turbulence of molecular cloud cores can naturally explain the observed asymmetrical features of bipolar molecular outflows such as the bendings and different sizes of redshifted and blueshifted lobes (Arce et al. 2013;Yen et al. 2015Yen et al. , 2017a;;Aso et al. 2018Aso et al. , 2019;;Okoda et al. 2021;Hsieh et al. 2023;Kido et al. 2023).The results perform that the bipolar outflows can be formed in not only the rotating cloud cores but also the turbulent cloud cores, which is consistent with previous theoretical studies (Matsumoto et al. 2017).
The bottom magenta box of Figure 1 shows that the unipolar outflow is driven in the model MF8 in which the initial cloud core has the weak magnetic field of µ = 8.The unipolar outflow evolves to be a large size of ∼ 10 3 au in the lower region and has the large opening angle in the bottom right panel of Figure 1.The magnetic field lines gradually becomes twisted only in the lower region as time proceeds.In contrast, the magnetic field lines in the upper region have extended and spread out structures without twisting.The results indicate that the unipolar outflow is mainly driven by getting accelerated with the magnetic pressure gradient force although the magneto-centrifugal force (Blandford & Payne 1982) slightly contributes to driving and accelerating it.Tomisaka (2002) indeed reports that the outflow driven by the magnetic pressure gradient force appears in the case of the weak magnetic field while the outflow driven by the magneto-centrifugal force appears in the case of the strong magnetic field because the toroidal components of the magnetic field can be amplified more easily by the disk rotation comparing to the case of the strong magnetic field.The results indicate that the strong turbulent accretions can locally amplify the toroidal components of the magnetic field relative to the poloidal ones, in contrasts to the bipolar outflow formed in model MF4.The results also show that the infalling envelopes have warped and elongated filamentary structures by the turbulent accretion.As shown in Figure 2 of the next subsection 3.2, these structures generated by the turbulent accretion can naturally explain the observed arc-like structures on the scale of ∼ 10 3 au (e.g., Tokuda et al. 2014).
The unipolar outflow is also driven in the model MF16 in which the initial cloud core has the very weak magnetic field of µ = 16.The outflows and initial parameters of the cloud cores are summarized in Table 1.It can be seen from Table 1 that the unipolar outflows are driven with E mag /E turb < 1, whereas the bipolar outflow is driven with E mag /E turb ∼ 1.The results of Figure 1 and Table 1 suggest that the ratio of the magnetic and turbulent energies of the parent cloud core, as represented by E mag /E turb , may play a key role in the driven outflow morphologies.

formation and evolution of unipolar outflow
Next, we focus on the formation and subsequent evolution of the unipolar outflow.Figure 2 shows the evolution of the surface density distributions along the y-direction for model MF8.We note that the origin of a coordinate system is shifted to the center of mass of the system.
Panel (b) of Figure 2 shows that the unipolar outflow is driven around the protostar at t p ∼ 8×10 3 yr after the protostar formation.Subsequent evolution from panels (b) to (c), the unipolar outflow grows up with the outflow speed of 1 − 5 km s −1 and expands to the scale of ∼ 10 3 au in a timescale of ∼ 5 × 10 3 yr.Panels (c) to (f) show that the unipolar outflow has a highly wide opening angle.The velocity of the unipolar outflow gradually becomes large during the evolution from panels (b) to (f).As time proceeds from panels (b) to (f), the velocity of the central protostar v p evolves to be from subsonic (v p ∼ 0.5c s,iso ) to supersonic (v p ∼ 1.6c s,iso ).
All the panel of Figure 2 shows that arc-like structures appear from ∼ 10 2 au to ∼ 10 3 au scales in the x−z plane.The arc-like structures in the panels of Figure 2 are infalling by the turbulent accretion and not outflowing.Tokuda et al. (2014) have reported an arc-like structure on ∼ 10 3 au scale around a very low-luminosity protostar in the dense cloud core MC27/L1521F by ALMA observations.Many Figure 2. Evolution of the surface density distributions along the y-direction for model MF8 in which the unipolar outflow is driven around the protostar.Panels are labeled by t p that notes the elapsed time after the protostar formation epoch defined at the time when the central density becomes higher than 1.0 × 10 −11 g cm −3 .v p is the velocity of the protostar.White arrows of the panels show the cut-plane velocity at y = 0.The reference arrow plotted on the top right corresponds to 1 km s −1 .Sky-blue solid lines show the velocity contours at v r = 0, tracing the front lines of the outflowing gas, where v r is the radial velocity of the gas.The origin of a coordinate system is shifted to the center of mass of the system.interferometric observations have also detected similar arc-like structures, which have been recently called the accretion streamers, in a wide range from the cloud core scale of ∼ 10 4 au (e.g., Pineda et al. 2020) to the disk and envelope scales of 10 2 −10 3 au (e.g., Yen et al. 2014Yen et al. , 2017b;;Akiyama et al. 2019;Yen et al. 2019;Thieme et al. 2022;Garufi et al. 2022;Kido et al. 2023;Aso et al. 2023).Our results suggest that the accretion streamers on the scales of 10 2 − 10 3 au can be naturally explained by the filamentary envelope accretions caused by the turbulence of the parent cloud cores, as emerged in many previous numerical simulations of the collapse of self-gravitating low-mass cloud cores with the turbulence (e.g., Matsumoto & Hanawa 2011;Tsukamoto & Machida 2013;Matsumoto et al. 2017;Takaishi et al. 2020).
Figure 3 shows the evolution of the ratio of the ram pressure P ram and magnetic pressure P mag at the cut-plane of y = 0 for model MF8, following the analysis of Machida & Hosokawa (2020).Panels (a) to (f) of Figure 3, the ram pressure P ram is lower than the magnetic pressure P mag inside the lower region of the driven unipolar outflow (inside the sky-blue line).The result indicates that the 1.0 0.5 0.0 0.5 1.0 log (P ram /P mag ) (= log( plasma / ram )) magnetic pressure gradient force caused by the twisted magnetic field as shown in panels of model MF8 in Figure 1 gradually becomes large and continuously enhances the unipolar outflow getting accelerated until the end of the simulation of t p ∼ 2.3 × 10 4 yr.In the upper region, the ram pressure P ram , however, is always higher than the magnetic pressure P mag from panels (a) to (f) of Figure 3, suggesting that the outflow driving is suppressed and failed by the ram pressure due to the infalling envelopes.
In order to see the detailed evolution of the ram and magnetic pressures separately, we introduce the ram beta parameter as β ram = P thm /P ram , where P thm is the thermal pressure of the gas.Figure 4 shows the evolution of the ram beta parameter β ram outside the outflow (v r < 0, outside the sky-blue line) and the plasma beta parameter β plasma = P thm /P mag inside the outflow (v r > 0, inside the sky-blue line) at the cut-plane of y = 0 for model MF8.
Panel (a) of Figure 4 shows that β ram (color map outside the outflow) is asymmetrically distributed and its value in the upper region is smaller than that in the lower region before the unipolar outflow is 2.0 1.5 1.0 log ram (v r < 0), log plasma (v r > 0) driven.The results indicate that the asymmetric accretion due to the turbulence of the initial cloud core generates such different ram pressure distributions at the initial accretion stage if the turbulence dominates.Therefore, driving outflow is initially delayed and suppressed in the region with large values of the ram pressure, and the outflow is driven earlier in the region with a smaller value of the ram pressure.
Figure 4 shows that the low-β ram region gradually expands to ∼ 10 3 au in the upper region from panels (b) to (f), indicating that the ram pressure by the infalling envelopes increases and suppresses the outflow driving as time proceeds.Figure 4 also shows that low-β plasma region expands inside the region of the driven unipolar outflow.Therefore, it is expected that the ram pressure P ram keeps overcoming the magnetic pressure P mag in the further subsequent evolution and no outflow driving is sustained in the upper low-β ram region.

Protostellar rocket effect
Figure 5 shows the projected protostellar trajectories on the x − z plane relative to the center of the initial cloud core for all models.As shown in Figure 5, the protostellar systems driving the unipolar outflow in models MF8 and MF16 move from the inner to the outer regions of their parent cloud cores at a distance of approximately 5 × 10 2 au in the timescale of t p ∼ 2 × 10 4 yr because the unipolar outflow ejects the protostellar system due to the linear momentum transport from the unipolar outflow to the protostellar system.Figure 5 also shows that the projected protostellar velocities with the unipolar outflow highly increase as time proceeds, suggesting that the protostellar systems are accelerated by driving the unipolar outflow.This phenomenon is similar to the launch and propulsion of a rocket.In the following, we refer to the protostellar systems with the acceleration by the linear momentum transport of the driven unipolar outflow as the protostellar rocket.
The protostellar rocket amplifies the relative velocity of the infalling envelopes towards the central protostar and protoplanetary disk.Therefore, the results suggest that the increase of the ram pressure by the infalling envelopes, in the upper region as shown in panels (b) to (f) in Figures 3 and 4, is caused by the combination with such the protostellar rocket and the infalling envelopes.
One interesting finding is that the driving additional new outflows to the different directions against the unipolar outflow are prevented once the protostellar rocket forms, suggesting that the unipolar outflow is sustained via the above phenomenon.Furthermore, the protostellar rocket feedbacks itself to evolve and enhance the unipolar outflow driving further as shown in Figure 2, resulting that the feedback system performs to be the instability state.Hereafter, we refer to the instability state as the protostellar rocket effect.
We note that the protostellar systems driving no outflow and the bipolar outflow also move shorter distances of ∼ 10 2 au comparing those driving the unipolar outflows due to the linear momentum transport from initial chaotic accretions by the turbulence.The velocities of the protostellar systems driving the bipolar outflow and no outflow remain subsonic throughout the simulations.

observational signature of unipolar outflows
Although the unipolar outflows have been detected in many observations, it is possible that the bipolar outflow is observed as the unipolar outflow due to extinction effects by the geometry and/or the surrounding gas.We propose the expected linear polarization maps of the thermal emission from dust grains aligned with the magnetic field as the observational signature to identify whether the unipolar outflow is real or not.
Figure 6 shows the expected polarization maps of the thermal emission from dust grains aligned with the magnetic field for model MF8 in which the unipolar outflow is driven.The polarization maps are calculated from the relative Stokes parameters q and u by using the method described in Tomisaka (2011) (see also Appendix A for details).Panel (a) of Figure 6 shows that the expected polarization degree in both upper and lower regions has a large value of ∼ 10 % at the scale of ∼ 10 3 au.The result suggests that the expected polarization in the upper region is roughly the same as that in the lower region before driving the unipolar outflow.However, panels (b) to (c) of Figure 6 show that the expected polarization degree inside the region of the driven unipolar outflow (inside the sky-blue line) considerably decreases to ∼ 0−5 %.As shown and pointed out in Tomisaka (2011), the protostellar outflow has a low polarization degree because the toroidal components cancel out  : 1 c s, iso : 1 c s, iso : 1 c s, iso : 1 c s, iso : 1 c s, iso : 1 c s, iso : 1 c s, iso : 1 c s, iso : 1 c s, iso : 1 c s, iso : 1 c s, iso : 1 c s, iso : 1 c s, iso : 1 c s, iso : 1 c s, iso : 1 c s, iso : 1 c s, iso = 2 (no) = 4 (bipolar) = 8 (unipolar) = 16 (unipolar) Figure 5. Projected protostellar trajectories on the x − z plane relative to the center of the initial cloud core (gray point) for all models.t p is the elapsed time after the protostar formation.v p shows the velocity of the protostar.Gray arrows indicate the projected protostellar velocity on the x − z plane.The reference arrow plotted on the top right corresponds to 1 c s,iso = 1.9 × 10 4 cm s −1 .The symbols "S" and "E" denote the positions of the protostar at its formation (t p = 0) and at the end of the simulation.The numbers beside the symbols "S" and "E" mean the values of the dimensionless mass-to-flux ratio µ.A schematic illustration on the protostellar rocket induced by the unipolar outflow is plotted near the positions of "E8" and "E16".Gray circles indicate the same distances from the center of the initial cloud core.each other inside the outflow region.The depolarization along the line of sight therefore emerges by the configuration of the twisted magnetic field around the unipolar outflow.
During the subsequent evolution from panels (d) to (f) of Figure 6, the expected polarization in the upper region gradually increases to ∼ 14 % and has a large value relative to that in the lower region whereas the polarization values becomes slightly large inside the unipolar outflow.The results indicate that the expected polarization is quite different between the lower and upper regions once after driving the unipolar outflow.Thus, the unipolar outflow can be identified by the observation of the thermal dust polarization, although the detectable sensitivity should be investigated in our subsequent work.It should be pointed out in Figure 6 that the accretion streamers caused by the turbulent accretions of the infalling envelope have a low expected polarization degree.

SUMMARY AND DISCUSSION
Recent observations on star-forming regions have revealed that the protostellar outflows driven by the YSOs exhibit a variety of asymmetrical features, in particular the unipolar outflows.Although the observations have reported the unipolar outflows, the formation and early evolution of the unipolar outflows driven by low-mass protostars are yet unclear.This study investigates the formation and early evolution of the protostellar outflows with the asymmetrical features by using the threedimensional non-ideal MHD simulations of the gravitational collapse of magnetized turbulent isolated low-mass cloud cores.This paper presents, for the first time to our knowledge, the formation of the unipolar outflow in the early evolution of the low-mass protostellar system.Our results and findings are summarized as follows.
1.The unipolar outflows are driven by the protostellar systems formed in the weakly magnetized cloud cores with the dimensionless mass-to-flux ratios of µ = 8 and 16.In contrast, the bending bipolar outflow is driven by the protostellar system formed in the moderately magnetized cloud core with µ = 4. Furthermore, no outflow is driven by the protostellar system formed in the strongly magnetized cloud core with µ = 2.The results explain the observed asymmetrical features of the protostellar outflows.
2. The protostellar system is ejected by the unipolar outflow from the central dense region to the outer region of the parent cloud core.As a result, the protostellar system driving the unipolar outflow is gradually accelerated as time passes.The protostellar velocity evolves to be from subsonic to supersonic by the acceleration.This is very similar to the launch and propulsion of a rocket, and so we call the protostellar system ejected and accelerated by the unipolar outflow the protostellar rocket.
3. We find that the subsequent additional new outflows cannot be driven by the protostellar rocket until the end of the simulation.This is because they are suppressed by the ram pressure of the infalling envelopes, which is enhanced by the acceleration of the protostellar rocket itself.In the context of low-mass star formation, the unipolar outflows can set the protostellar rocket to be the instability state, and we call this phenomenon the protostellar rocket effect.
The remaining question is how the outflow morphologies change with the different turbulent energies E turb from the one adopted in this study.Our results show that the unipolar outflows are driven with E mag /E turb < 1, while the bipolar outflow is driven with E mag /E turb ∼ 1.This implies that the outflow morphologies may depend on the ratio of the magnetic and turbulent energies of the parent cloud core, E mag /E turb .
We suggest that the outflow morphologies can be used as a new tracer to indirectly estimate the magnetic field strengths of the parent molecular cloud cores.The measurements of the magnetic field strengths using the Zeeman effect show that the observed typical cloud cores are slightly supercritical with the dimensionless mass-to-flux ratios of µ obs ∼ 2 (e.g., Troland & Crutcher 2008;Falgarone et al. 2008;Crutcher 2012), which means that the magnetic pressure is not sufficient to prevent the gravitational collapse of the cloud cores.The measurements using the Davis-Chandrasekhar-Fermi (DCF) method also suggest that the cloud cores are magnetized with µ obs ∼ 2 − 3 (e.g., Kirk et al. 2006;Karoly et al. 2023).The cloud cores may actually form from somewhat subcritical initial conditions (e.g., Pattle et al. 2017;Karoly et al. 2020;Yin et al. 2021;Priestley et al. 2022;Ching et al. 2022;Karoly et al. 2023).However, the measurements of µ obs suffer fundamentally from their statistical and systematic uncertainties such as the selection bias towards sources with strong magnetic fields for the Zeeman effect observations and the overestimates of the magnetic field strengths in the DCF method (Liu et al. 2021).
In our simulations, a bipolar outflow forms in a cloud core with µ = 4 (µ const = 2) and M s = 0.86.Typical cloud cores have the turbulence of M s ≲ 1 (e.g., Ward-Thompson et al. 2007).Therefore, our results suggest that the formation of bipolar outflows in turbulent cloud cores requires the relatively strong magnetic field of the parent cloud cores.In contrast, the unipolar morphology may form with the relatively weak magnetic field of the parent cloud cores.Wu et al. (2004) show that bipolar outflows are observed more frequently than unipolar outflows.This suggests that the typical molecular cloud cores tend to have relatively strong magnetic fields.Note, however, that the turbulence strength is fixed in our current study and the impact of different turbulence strengths remains unclear.We will investigate the relation between the outflow morphologies, magnetic field strengths (E mag ), and turbulence strengths (E turb ) in our future work.
The protostellar rocket effect would drive shock waves into the envelope around the protostar, which may possibly be detected by some chemical shock tracers in the molecular line emissions.However, the detectability of the shock waves depends on the chemical species of the shock tracers and their lifetimes in the post-shock waves.In addition, even just accreting gas, without the protostellar rocket effect, is capable of driving the accretion shocks.Therefore, the detectability of the shock waves in the ambient cloud material needs to be investigated in detail to identify the probable chemical shock tracers with their characteristics.
This study ignores the Hall effect due to the computational cost to include it.Here, we briefly discuss possible impact of the Hall effect on our conclusion.The magnetic braking is strengthened or weakened by the Hall effect when the rotation and magnetic field vectors are aligned or anti-aligned, and the Hall effect introduces some interesting phenomena in the formation and early evolution of the protostars and protoplanetary disks in collapsing cloud cores (e.g., Wardle & Ng 1999;Krasnopolsky et al. 2011;Li et al. 2011;Braiding & Wardle 2012;Tsukamoto et al. 2015b;Wurster et al. 2016;Tsukamoto et al. 2017;Wurster et al. 2018aWurster et al. ,b, 2021)).If the Hall effect is strong enough even in our simulation environments, we speculate that when the rotation and magnetic field vectors are aligned, the magnetic field may be twisted more strongly and the outflow is stronger, making it harder for the unipolar outflow to form.On the other hand, when the rotation and magnetic field vectors are anti-aligned, the magnetic field may be twisted more weakly and the outflow is less likely to overcome the ram pressure, which may lead to the formation of the unipolar outflow.However, it should be noted that the magnetic resistivity depends on the dust model, cosmic ray ionization rate, and also on the strength of the magnetic field (e.g., Wardle & Ng 1999;Nakano et al. 2002;Padovani et al. 2014;Marchand et al. 2016;Wurster et al. 2018c;Koga et al. 2019;Tsukamoto & Okuzumi 2022;Kawasaki et al. 2022;Kobayashi et al. 2023;Tsukamoto et al. 2023b).These predictions will be verified by future simulations including the Hall effect.
How long the protostellar rocket effect continues in subsequent evolution is still an open question.More evolved protostellar rockets are expected to be distributed outside the region of the initial cloud core as long as the unipolar outflow is maintained.Thus, the unipolar outflow and the protostellar rocket may be observationally identified by combining the proper motions, the launch direction of the unipolar outflow on the plane of the sky, and the morphologies of the infalling envelope cavity using, for example, the Global Astrometric Interferometer for Astrophysics (GAIA), the James Webb Space Telescope (JWST), ALMA, and/or other methods.However, even if peculiar proper motions are detected, it may be difficult to distinguish between protostellar rockets being its origin and, for example, stellar encounters being its origin.The long-term evolution and observability of protostellar rockets will be considered in our future studies.

Figure 1 .
Figure1.Three-dimensional views of the density distribution and magnetic field structure for models of the no outflow (model MF2, top green box), bipolar outflow (model MF4, middle blue box), and unipolar outflow (model MF8, bottom magenta box).Yellow isosurfaces show the density of 3.2×10 −18 g cm −3 , 10 −17 g cm −3 , 3.2 × 10 −17 g cm −3 , and 10 −16 g cm −3 with the radial velocity of v r > 0, representing the outflows.White lines show the magnetic field lines.Blue and green isosurfaces show the density of 3.2 × 10 −17 g cm −3 and 10 −16 g cm −3 with v r < 0, representing the infalling envelopes.Cut-plane densities on x − y (z = 0), x − z (y = 0), and y − z (x = 0) planes are projected for each panel.t p notes the elapsed time after the protostar formation epoch defined at the time when the central density becomes higher than 1.0 × 10 −11 g cm −3 .The scale of the box is ∼ 2, 000 au.The origin of a coordinate system is shifted to the center of mass of the system.

Figure 3 .
Figure 3. Evolution of the ratio of the ram pressure P ram and magnetic pressure P mag at the cut-plane of y = 0 for model MF8.Black lines show the contours of the ratio.t p of each panel corresponds to that in Figure 2. v p is the velocity of the protostar.White lines show the contours of the surface density in Figure 2 ranging from −0.75 to 1.25 in 0.25 steps on a logarithmic scale.Black arrows of the panels show the cut-plane velocity at y = 0.The reference arrow plotted on the top right corresponds to 1 km s −1 .Sky-blue solid lines show the velocity contours at v r = 0, tracing the front lines of the outflowing gas.The origin of a coordinate system is shifted to the center of mass of the system.

Figure 4 .
Figure 4. Evolution of the ram beta parameter β ram = P thm /P ram outside the outflow (v r < 0) and the plasma beta parameter β plasma = P thm /P mag inside the outflow (v r > 0) at the cut-plane of y = 0 for model MF8, where P thm is the thermal pressure of the gas.Black lines show the contours of them.Sky-blue solid lines show the velocity contours at v r = 0, indicating the boundary between outside and inside of the outflow.t p of each panel corresponds to that in Figure 2. v p is the velocity of the protostar.White lines show the contours of the surface density in Figure 2 ranging from −0.75 to 1.25 in 0.25 steps on a logarithmic scale.White arrows of the panels show the cut-plane velocity at y = 0.The reference arrow plotted on the top right corresponds to 1 km s −1 .The origin of a coordinate system is shifted to the center of mass of the system.

Figure 6 .
Figure 6.Expected linear polarization along the y-direction for the model of the unipolar outflow (model MF8).The color maps in each panel show the polarization degree.The polarization degree vectors are plotted by red-brown bars in each panel.The elapsed times t p of each panel are the same as in Figure 2. White lines show the contours of the surface density in Figure ranging from −0.75 to 1.25 in 0.25 steps on a logarithmic scale.White of the panels show the cut-plane velocity at y = 0.The references of the polarization degree vector and cut-plane velocity are plotted on the top right.Sky-blue solid lines show the velocity contours at v r = 0, tracing the front lines of the outflowing gas.The origin of a coordinate system is shifted to the center of mass of the system.

Table 1 .
The model names and parameters Model µ B c (µG) µ const γ mag E mag /E turb Outflow Note-µ is the dimensionless mass-to-flux ratio.B c is the strength of the central magnetic field.µ const is the dimensionless mass-to-flux ratio with the constant value of the central magnetic field.γ mag = E mag /|E grav | is the ratio of the magnetic to gravitational energies.E mag /E turb = γ mag /γ turb is the ratio of the magnetic to turbulent energies.The last column indicates the morphology of the driven outflow.