Circumstellar disks and outflows in turbulent molecular cloud cores: possible formation mechanism for misaligned systems

TL;DR

This study uses simulations to explore how magnetic fields and turbulence influence the formation, size, and alignment of circumstellar disks and outflows in turbulent molecular cloud cores, revealing mechanisms for misaligned systems.

Contribution

It introduces detailed numerical simulations showing the impact of magnetic fields and turbulence on disk formation, size, and misalignment with outflows in star-forming regions.

Findings

Magnetized disks are smaller than nonmagnetized ones but grow over time.

Stronger magnetic fields lead to smaller disks.

Disks tend to align with their angular momentum, but outflows can be misaligned due to magnetic and geometric effects.

Abstract

We investigate the formation of circumstellar disks and outflows subsequent to the collapse of molecular cloud cores with the magnetic field and turbulence. Numerical simulations are performed by using an adaptive mesh refinement to follow the evolution up to $\sim 1000$~yr after the formation of a protostar. In the simulations, circumstellar disks are formed around the protostars; those in magnetized models are considerably smaller than those in nonmagnetized models, but their size increases with time. The models with stronger magnetic field tends to produce smaller disks. During evolution in the magnetized models, the mass ratios of a disk to a protostar is approximately constant at $\sim 1-10$\%. The circumstellar disks are aligned according to their angular momentum, and the outflows accelerate along the magnetic field on the $10-100$~au scale; this produces a disk that is misaligned with the outflow. The outflows are classified into two types: a magneto-centrifugal wind and a spiral flow. In the latter, because of the geometry, the axis of rotation is misaligned with the magnetic field. The magnetic field has an internal structure in the cloud cores, which also causes misalignment between the outflows and the magnetic field on the scale of the cloud core. The distribution of the angular momentum vectors in a core also has a non-monotonic internal structure. This should create a time-dependent accretion of angular momenta onto the circumstellar disk. Therefore, the circumstellar disks are expected to change their orientation as well as their sizes in the long-term evolutions.

Tables

Table 1. Table 1 Model parameters

Model	$ℳ$	$\alpha$	${ℳ}_c$	$\mu$
M05B0	0.5	0.0	0.35	$\infty$
M1B0	1.0	0.0	0.71	$\infty$
M05B01	0.5	0.1	0.35	2.81
M1B01	1.0	0.1	0.71	2.81
M1B025	1.0	0.25	0.35	1.12
M05B025	0.5	0.25	0.71	1.12

Equations

$$\rho(r)=\left\{\begin{array}[]{ll}\rho_{0}\varrho_{\rm BE}(r/a)&{\rm for}\;r<R_{c}\\ \rho_{0}\varrho_{\rm BE}(R_{c}/a)&{\rm for}\;r\geq R_{c}\end{array}\right.\;,$$

$$a=c_s{\left(\frac{f}{4\pi G{\rho }_0}\right)}^{1/2},$$

$${\cal M}=\left(\frac{1}{V_{\mathrm{cd}}}\int_{V_{\mathrm{cd}}}|\mbox{\boldmath$v$}^{2}|\,dV\right)^{1/2},$$

$$I=\left(\begin{array}[]{ccc}I_{xx}&I_{xy}&I_{xz}\\ I_{yx}&I_{yy}&I_{yz}\\ I_{zx}&I_{zy}&I_{zz}\\ \end{array}\right)\;,$$

$$I_{xy}=\frac{{\int }_{V_d}\left(x-x_p\right)\left(y-y_p\right)dV}{{\int }_{V_d}dV},$$

$$R_d=2{\lambda }_2^{1/2},$$

$$M_d={\int }_{V_d}\rho dV,$$

$$\displaystyle\overline{\mbox{\boldmath$B$}}(r)$$

$$\displaystyle\overline{\mbox{\boldmath$j$}}(r)$$

$$M(r)={\int }_{V_s(r)}\rho dV,$$

$$V_{s}(r)=\{\mbox{\boldmath$r$}\in\mathbb{R}^{3}\mid|\mbox{\boldmath$r$}-\mbox{\boldmath$r$}_{p}|\leq r\}.$$

$$I_{xy}=\frac{{\int }_{\rho \ge {\rho }_d}\rho \left(x-x_p\right)\left(y-y_p\right)dV}{{\int }_{\rho \ge {\rho }_d}\rho dV},$$

$$\overline{\mbox{\boldmath$v$}}_{\mathrm{of}}(r)=\frac{1}{V_{\mathrm{of}}(r)}\int_{V_{\mathrm{of}}(r)}\mbox{\boldmath$v$}\,\mathrm{sign}(\mbox{\boldmath$v$}\cdot\overline{\mbox{\boldmath$B$}})\,dV,$$

$$V_{\mathrm{of}}(r)=\{\mbox{\boldmath$r$}\in V_{s}(r)\mid v_{r}(\mbox{\boldmath$r$})\geq 2c_{s}\}.$$

Full text

††slugcomment: Accepted by ApJ: March 27, 2017

1. Introduction

The magnetic field and turbulence play important roles in the early phase of star formation. Observations of the magnetic field have indicated that molecular clouds and molecular cloud cores have a large amount of magnetic energy, which is approximately equal to the kinetic energy (Crutcher, 1999). The magnetic field therefore has the potential to control the gravitational collapse of cloud cores. Molecular clouds exhibit broad molecular lines, which are interpreted as supersonic turbulence (Zuckerman & Evans, 1974). The turbulence seems to have a scaling relation such that a smaller scale has a smaller velocity dispersion (Larson, 1981). For high-density cloud cores, weak turbulence is suggested by the narrow molecular lines (e.g., Onishi et al., 1998). As shown by Burkert & Bodenheimer (2000), such turbulence reproduces the observed rotational properties of cloud cores. This indicates that the turbulence contributes angular momentum to the cloud cores, and it is the origin of the rotation of circumstellar disks and protostars.

The magnetic field extracts angular momentum from the circumstellar disks and the infalling envelopes, due to magnetic braking and outflows. As a consequence, the protostars accrete gas from the disks. In the past decade, the existence of the so-called “magnetic braking catastrophe” has been debated (e.g., Mellon & Li, 2008); in this phenomenon, the magnetic field prevents the formation of the circumstellar disks around the protostars. Axisymmetric models have been investigated intensively, and they show that circumstellar disks are formed around the protostar. The size of these disks increases with time, but they are considerably smaller than nonmagnetized disks (e.g., Machida et al., 2011). Models in which the magnetic field and rotation axes are misaligned have also been investigated (Matsumoto & Tomisaka, 2004; Li et al., 2013). Few studies have investigated the role of turbulence in disk formation (e.g., Seifried et al., 2012, 2013). Matsumoto & Hanawa (2011) performed numerical simulations of the collapse of magnetized turbulent cloud cores, but they only followed the evolution up to the formation of the first core. Their simulations reproduced the formation of the outflow but not that of disks, because the period simulated was too short.

Recent observations have revealed misaligned young stellar objects. Hull et al. (2013) reported that for 16 Class 0 and Class I sources, the magnetic field in protostellar cores of $\sim 1000$ au scale is not tightly aligned with outflows. It has been suggested that in the Class I source L1489 IRS, the central Keplerian disk is inclined with respect to a flattened infalling envelope (Brinch et al., 2007a, b). Moreover, it has been suggested that the class I binary source IRS 43 exhibits a misalignment between the orbit of the binary and the circumbinary disk (Brinch et al., 2016). The formation mechanism for such misaligned systems is not yet known.

In this paper, we use high-resolution numerical simulations to investigate the collapse of cloud cores to form protostars, and we include the effects of the magnetic field and turbulence of the cloud core. The formation of circumstellar disks and outflows is also investigated. In the simulations, it is expected that misaligned protostars will be formed, because a priori rotation axes are not assumed in the initial conditions. We focus on the early phase of protostar formation, because recent observations have provided high-resolution images of very young protostars and circumstellar structures (e.g., Tokuda et al., 2014).

This paper is organized as follows. Section 2 presents the model. Section 3 discusses the simulation methods, the results are presented in Section 4, and they are discussed in Section 5. Our conclusions are presented in Section 6.

2. Models

As the initial state of a molecular cloud core, we consider a turbulent, spherical cloud threaded by a uniform magnetic field. The cloud is confined by a uniform ambient gas. This is similar to the initial state considered by Matsumoto & Hanawa (2011). It is further specified by the initial strength of the turbulence and the magnetic field strength, as summarized in Table 1.

As a template for a molecular cloud core, we consider a cloud for which the density profile is that of the critical Bonnor-Ebert sphere (Bonnor, 1956; Ebert, 1955). We let $\varrho_{\rm BE}(\xi)$ denote the nondimensional density profile of the critical Bonnor-Ebert sphere (see Chandrasekhar, 1939), and then the initial density distribution is given by

[TABLE]

and

[TABLE]

where $r$, $G$, $c_{s}$, and $\rho_{0}$ denote the radius, gravitational constant, isothermal sound speed, and initial central density, respectively. The gas temperature is assumed to be 10 K ($c_{s}=0.19~{}\mathrm{km}~{}\mathrm{s}^{-1}$) . The initial central density is set at $\rho_{0}=10^{-18}\,{\rm g}\,{\rm cm}^{-3}$, which corresponds to a number density of $n_{0}=2.61\times 10^{5}\,{\rm cm}^{-3}$ for an assumed mean molecular weight of 2.3. The parameter $f$ denotes the nondimensional density enhancement factor, and the critical Bonnor-Ebert sphere is obtained when $f=1$. For a given central density, an increase in density by a factor of $f$ is equivalent to an enlargement of the spatial scale by a factor of $f^{1/2}$. We adopt $f=2$ in this paper. The radius of the cloud is set to be $R_{c}=6.45a=0.0434f^{1/2}\,\mathrm{pc}=0.0614\,\mathrm{pc}$, where the factor 6.45 comes from the nondimensional radius of the critical Bonnor-Ebert sphere. The density contrast of the initial cloud is $\rho(0)/\rho(R_{c})=14.0$. The initial freefall timescale at the center of the cloud is thus $t_{\rm ff}\equiv(3\pi/32G\rho_{0})^{1/2}=6.66\times 10^{4}\,{\rm yr}$. The mass of the cloud core ($r\leq R_{c}$) is $M_{c}=0.89f^{3/2}M_{\sun}=2.51\,M_{\sun}$. The ratio of the thermal energy to the gravitational energy is estimated as $E_{\mathrm{th}}/|E_{\mathrm{grav}}|=0.84f^{-1}=0.42$. The spherical cloud described above is located at the center of the computational domain of $x,y,z\in[-2R_{c},2R_{c}]^{3}$.

The turbulence is determined by the initial velocity field, and it is not driven during the course of the simulations; that is, we assume free decay of the turbulence. The initial velocity field is incompressible, with a power spectrum of $P(k)\propto k^{-4}$, and it is generated in accordance with that in Dubinski et al. (1995), where $k$ is the wavenumber. This power spectrum results in a velocity dispersion of $\sigma(\lambda)\propto\lambda^{1/2}$, which is in agreement with the Larson scaling relations (Larson, 1981), where $\lambda$ denotes the length scale. The root-mean-square (rms) Mach number in the computational domain is specified by a model parameter, ${\cal M}$;

[TABLE]

where $V_{\mathrm{cd}}$ denotes the volume of the computational domain. We utilized a common velocity field for generating the turbulence in all the models. We changed the amplitude of the turbulence by changing the Mach number ${\cal M}$. Note that even if we assume ${\cal M}=1$, the rms Mach number on a 0.06 pc scale (the cloud core scale) is estimated to be ${\cal M}_{c}={\cal M}(R_{c}/2R_{c})^{1/2}=0.71$, according to the scaling relations, and therefore the turbulence is subsonic on the cloud core scale. When ${\cal M}=0.5$, the rms Mach number is estimated to be ${\cal M}_{c}=0.35$. Note that weak turbulence was suggested by the narrow molecular lines in the dense cores in Taurus (Onishi et al., 1998). The rms Mach number on the cloud core scale ${\cal M}_{c}$ for each model is listed in Table 1.

The initial magnetic field is uniform in the $z$-direction. The field strength is given by $B_{z}=\alpha B_{\mathrm{cr}}$, where $\alpha$ denotes the nondimensional flux-to-mass ratio (see Table 1), and $B_{\mathrm{cr}}$ denotes the critical field strength for the center of the cloud core, given by $B_{\mathrm{cr}}=2\pi G^{1/2}\Sigma_{0}$ (Nakano & Nakamura, 1978; Tomisaka et al., 1988). The central column density $\Sigma_{0}$ is calculated by $\Sigma_{0}=\int_{-R_{c}}^{R_{c}}\rho dz=5.38\rho_{0}a$, where the integral is performed along a line passing through the center of the cloud core. In this paper, we consider a magnetically supercritical core ($\alpha<1$). The initial field strength is estimated to be $B_{z}=181\alpha f^{1/2}\,\mu\mathrm{G}=256\alpha\,\mu\mathrm{G}$; note this uses the model parameters $\alpha$ and $f$. Note that the model parameter $\alpha$ is inversely proportional to the dimensionless mass-to-flux ratio $\mu$, which is defined as $\mu=(M_{c}/\Phi)/(M_{c}/\Phi)_{\mathrm{cr}}$; its value for each model is listed in Table 1. The magnetic flux is defined by $\Phi=\pi R_{c}^{2}B_{z}$, and the mass of the cloud core is defined by $M_{c}=4\pi\int_{0}^{R_{c}}\rho r^{2}dr$. The critical value is $(M_{c}/\Phi)_{\mathrm{cr}}=(2\pi G^{1/2})^{-1}$. The parameter $\mu$ is the mass-to-flux ratio in the entire cloud core, while the parameter $\alpha$ reflects the mass-to-flux ratio only at the central axis.

The barotropic equation of state is assumed as $P(\rho)=c_{s}^{2}\rho+\kappa\rho^{7/5}$, where $\kappa=c_{s}^{2}\rho_{\mathrm{cr}}^{-2/5}$; the critical density is set at $\rho_{\mathrm{cr}}=10^{-13}\,\mathrm{g}\,\mathrm{cm}^{-3}$ (the corresponding number density is $n_{\mathrm{cr}}=2.62\times 10^{10}\,\mathrm{cm}^{-3}$), which is taken from the numerical results of Masunaga, Miyama, & Inutsuka (1998).

The ohmic dissipation is also considered. The resistivity $\eta$ is quantitatively estimated according to equations (9) and (10) of Machida et al. (2007). The magnetic Reynolds number is estimated as $Re_{m}=v_{f}\lambda_{J}\eta^{-1}$, where $v_{f}=[(4/3)\pi G\lambda_{J}^{2}\rho]^{1/2}$ is the free-fall velocity, $\lambda_{J}=(\pi c_{p}^{2}/G\rho)^{1/2}$ is the Jeans length, and $c_{p}=(dP/d\rho)^{1/2}$ is the sound speed for the barotropic equation of state. The magnetic Reynolds number is less than unity ($Re_{m}<1$) in the region where $n\gtrsim 2\times 10^{12}\,\mathrm{cm}^{-3}$, and thus the magnetic field is dissipative in that region. This indicates that the inner portion of the first core is magnetically dissipative.

We allowed the models listed in Table 1 to develop until $t_{p}\simeq 10^{3}$ yr for the magnetized models ($\alpha\neq 0$) and $t_{p}\simeq 10^{4}$ yr for the nonmagnetized models ($\alpha=0$), where $t_{p}$ denotes the elapsed time following formation of a sink particle (a model of a protostar).

The short elapsed times in the magnetized models are due to short time steps of the simulations, which are caused by an extremely fast Alfvén speed around the sink particles. In this study, we therefore focus on the early evaluational stages of low-mass star formation. The recent observational studies have investigated young stars in such early evolutionary stages, e.g., the first hydrostatic core candidates (e.g., Belloche et al., 2006; Pineda et al., 2011).

3. Methods

Gravitational collapse of the cloud cores was calculated using the three-dimensional adaptive mesh refinement (AMR) code SFUMATO (Matsumoto, 2007). The magnetohydrodynamics (MHD) scheme has second-order accuracy in space and time. The computational domain was resolved on a base grid of $l=0$ with $256^{3}$ cells. The maximum grid level was set at $l_{\mathrm{max}}=9$. The cell width was $\Delta x_{\mathrm{min}}=0.386$ au on the finest grid, $l=9$, and $\Delta x_{\mathrm{max}}=198$ au on the base grid of $l=0$. The Jeans condition was employed as a refinement criterion: blocks were refined when the Jeans length was shorter than eight times the cell width, i.e., $\lambda_{J}<8\Delta x$, where $\lambda_{J}$ is the Jeans length.

This condition is twice as strict as that originally proposed by Truelove et al. (1997), and tested in Matsumoto (2007). An even more strict refinement criterion was proposed by Federrath et al. (2011) to capture dynamo amplification of the very weak magnetic field in the gravitational collapse of primordial clouds. Although it should be better to use such a criterion to describe the small scale dynamo action in very high beta plasma, a very weak turbulent magnetic field does not appear anywhere in the present models.

For the MHD scheme, we adopted the HLLD Riemann solver (Miyoshi & Kusano, 2005), though the Roe-type MHD Riemann solver (Fukuda & Hanawa, 1999) was implemented on the original SFUMATO (Matsumoto, 2007). The HLLD Riemann solver is known to be more robust than the Roe-type scheme, because it preserves positivity. We did not apply any density floors, even when the plasma beta became low in the low-density regions.

For the sub-grid model of a protostar, we used sink particles. The details of the implementation of the sink particles is shown in Matsumoto et al. (2015a). The critical density for sink particle formation is set at $\rho_{\mathrm{sink}}=1\times 10^{-10}\,\mathrm{g}\,\mathrm{cm}^{-3}$ ($n_{\mathrm{sink}}=2.62\times 10^{13}\,\mathrm{cm}^{-3}$), and the sink radius is set at $r_{\mathrm{sink}}=4\Delta x_{\mathrm{min}}=1.55~{}\mathrm{au}$.

The minimum cell width $\Delta x_{\mathrm{min}}$ is determined so that the Jeans length of $\rho_{\mathrm{sink}}$ is resolved by the cells width $\Delta x_{\mathrm{min}}$ according to the Jeans condition.

The ohmic dissipation was calculated in accordance with Matsumoto (2011). The induction equation for the magnetic field was split into a hyperbolic term and a parabolic term. The former corresponds to the ideal MHD and was solved explicitly; the latter corresponds to the ohmic dissipation and was solved implicitly with the multigrid AMR. This provides second-order spatial accuracy and first-order temporal accuracy.

4. Results

4.1. Overall structures

Figures 1–3 show snapshots of the models at $t_{p}\simeq 10^{3}$ yr on the cloud core scale, the infalling envelope scale, and the circumstellar disk scale, respectively. Hereafter, we will refer to a circumstellar disk simply as a disk.

The cloud cores collapse on the timescale of $\sim 10^{5}$ yr. In general, a stronger initial turbulence and/or a stronger initial magnetic field retards the formation epoch of the sink particle. The sink particle forms at $t=1.2\times 10^{5}$ yr for model M05B0 (the earliest formation) and $t=1.9\times 10^{5}$ yr for model M1B025 (the latest formation).

On the cloud core scale (Figure 1), the cloud cores are deformed from a spherical shape. The nonmagnetized models and the weak magnetic field models (the left and middle panels of Figure 1) produce complex shapes for the cloud cores because of the turbulence, while the strong magnetic field models (the right panels of Figure 1) produce flattened cloud cores, which are oriented perpendicular to the mean magnetic field (the $z$-direction). The upper and lower panels in Figure 1 show the models with weak turbulence (${\cal M}=0.5$) and moderate turbulence (${\cal M}=1$), respectively. Models with moderate turbulence result in shapes that are more disturbed than those resulting from models with weak turbulence; this is consistent with what was reported by Matsumoto & Hanawa (2011).

Figure 2 shows the column density and velocity distributions on the envelope scale of $400~{}\mathrm{au}$. The nonmagnetized models and weak magnetic field models (left and middle panels of Figure 2) result in rotating infalling envelopes. The rotation velocities are comparable to the infall velocities. The angular momentum of the infalling envelopes comes from the turbulent velocities imposed on the initial conditions. The density distributions in the envelopes show spiral structures, and these are caused by the inhomogeneity of the rotation velocities. On the other hand, the strong magnetic field models with $\alpha=0.25$ (the right panels of Figure 2) produce flattened infalling envelopes that are perpendicular to the direction of the mean magnetic field. The arrows in these models reflect the silhouettes of the bipolar outflows along the $z$-axis. Outflow formation is discussed in Section 4.3.

Figure 3 shows the column density and velocity distribution for the central $(50~{}\mathrm{au})^{3}$ cubes. The disks show almost face-on views for all the models because a common seed was adopted for the initial turbulence, which provides angular momentum to the cloud cores. In each panel, the measured size of the disk is indicated by the green circle; the following two criteria were used to define a disk: (1) the density is higher than $\rho_{\mathrm{cr}}$, and (2) the rotation speed is considerably faster than the infall velocity. The measurement of the disk radii is shown in Appendix A.

The size of each disk is correlated with the strength of the magnetic field; in models with a stronger magnetic field, the disks are smaller, as shown in Figure 3. The strong magnetic field model M05B025 has a small disk with a radius of $R_{d}=3.36$ au at $t_{p}=10^{3}$ yr, though this radius is larger than the sink radius, $r_{\mathrm{sink}}=1.55~{}\mathrm{au}$. The disk is oriented nearly perpendicular to the $y$-axis (face-on view in Figure 3), and it is also perpendicular to the flattened infalling envelopes. In the disk, the rotation velocity is considerably greater than the infall velocity.

Model M1B025 has also a small disk, with a radius of $3.08$ au at $t_{p}=10^{3}$ yr; note that it is larger than the sink radius. The disk is elongated toward the upper left in Figure 3. The high-density portion of the disk, where $\rho\geq 10^{6}\rho_{0}=10^{-12}\,{\rm g}\,{\rm cm}^{-3}$, is approximately axisymmetric, while the edge of the disk with $\rho\sim 10^{-13}\,{\rm g}\,{\rm cm}^{-3}$ has a cometary shape.

For the weak field models with $\alpha=0.1$ (models M05B01 and M1B01), the disks have the face-on views shown in Figure 3. Each disk is divided into inner and outer parts. The inner disk has an approximately axisymmetric shape with a high density, and its radius is estimated to be $R_{d}\sim 5~{}\mathrm{au}$ at $t_{p}\sim 10^{3}$ yr, which is considerably larger than the sink radius. The outer disk has spiral arms that wind around the disk; this occurs because the rotation velocity is greater than the infall velocity.

We found that the spiral feature is also caused by the inhomogeneity of the ratio between the thermal pressure and magnetic pressure, and the inhomogeneity of the rotation velocity along the azimuthal direction. The outer disk exhibits the $Q$-value (Toomre, 1964) larger than $\sim 2$, indicating that it is gravitationally stable.

The Atacama Large Millimeter/submillimeter Array (ALMA) recently observed a circumstellar disk with a similar dual structure in ELias 2-27, a Class II star (Pérez et al., 2016).

The spiral arms of this object is probably caused by the gravitational instability (see also Tomida et al., 2017), in contrast to the models here.

The nonmagnetized models, M05B0 and M1B0, have disks with a large radius of $R_{d}\sim(20-30)$ au at $t_{p}\sim 10^{3}$ yr. For model M1B0, the disk fragments three times: at $t_{p}=7.1\times 10^{3}~{}\mathrm{yr}$, $7.8\times 10^{3}~{}\mathrm{yr}$, and $8.6\times 10^{3}~{}\mathrm{yr}$. Two of these fragments merge at $t_{p}=8.9\times 10^{3}~{}\mathrm{yr}$, and at the end of the calculation period, $t_{p}\simeq 10^{4}~{}\mathrm{yr}$, there are still three sink particles. Similar fragmentation is also seen when a filamentary cloud is used as the initial condition (Matsumoto et al., 2015b). For model M05B0, the disk has not fragmented by the end of the calculation period, $t_{p}\simeq 10^{4}~{}\mathrm{yr}$.

For all the models, the thicknesses of the disks are resolved by more than four cells. This is consistent with the fact that a self-gravitational disk has the scale hight $H=c_{p}/(2\pi\gamma G\rho)^{1/2}=\lambda_{J}/(\sqrt{2\gamma}\pi)=0.19\lambda_{J}$ (e.g., Larson, 1985), and the Jeans length adopted here requires $2H\gtrsim 3\Delta x_{\mathrm{min}}$.

4.2. Disk formation

Figure 4 shows the rotation velocity profile as a function of the distance form the sink particle for each of the models. The rotation velocity profile is obtained as follows. The orientation of the disk axis is determined according to the total angular momentum for a volume of $\rho\geq\rho_{\mathrm{cr}}$. According to the disk orientation, the azimuth velocity $v_{\varphi}$ is calculated with respect to the sink particle, and it is averaged with a density weight along the vertical direction of the disk for a volume of $\rho\geq\rho_{\mathrm{cr}}$. Finally, the density-weighted $v_{\varphi}$ is azimuthally averaged to obtain the rotation velocity profile shown in Figure 4. The nonmagnetized models (M05B0 and M1B0) exhibit rotation velocity profiles faster than those of the Keplerian rotation because of the massive disks. The self-gravity of the disk increases the rotation velocity. For the weak field models (M05B01 and M1B01), the inner parts of the disks ($r\lesssim(5-7)$ au) exhibit the Keplerian rotation. These regions correspond to those showing an approximately axisymmetric shape in Figure 3. In the outer parts of the disks, rotation velocity is slower than that of the Keplerian rotation, and the infall motion is observed there as shown in Figure 3. For the strong field models (M05B025 and M1B025), the disk sizes are as small as $\sim 3$ au, and the rotation velocities are close to those of the Keplerian rotation near the sink radius. We confirmed that the centrifugal force is a dominant force against the gravity in the disks, and the rotation velocity is considerably faster than the infall velocity. Note that, for all the models, the regions within the sink radius exhibit slower rotation than the Keplerian rotation because of the softening of the gravity therein.

Figure 5 shows the increase in the radius of the disks as a function of time following the formation of sink particles for each of the models. The strong magnetized models ($\alpha=0.25$) have disks that are smaller than those for the nonmagnetized models ($\alpha=0$). The disk radius for the strong field models remain at $R_{d}\sim 3$ au during the simulation period of $6\times 10^{2}\,\mathrm{yr}\lesssim t_{p}\lesssim 10^{3}\,\mathrm{yr}$. In the weak magnetic field models ($\alpha=0.1$), the disk radius increases with considerable undulations, which are caused by the dynamical changes in the spiral arms. The radius of the disks is sensitive to changes in the spiral arms associated with the outer parts of the disks.

Figure 6 shows the mass of the disks as a function of time following the formation of sink particles. The stronger magnetic field models show slower growth of the disk mass, and this is roughly independent of the Mach number ${\cal M}$. The models with $\alpha=0.25$ have a disk mass of $M_{d}\sim 10^{-3}\,M_{\sun}$ at $t_{p}=10^{3}$ yr, while models with $\alpha=0.1$ have a disk mass of $M_{d}\sim 8\times 10^{-3}\,M_{\sun}$.

Figure 7 shows the mass of the sink particles as a function of time following the formation of sink particles. The growth of the sink particles exhibits a clear tendency: magnetic field models have higher accretion rates than do the nonmagnetized models; all the magnetized models have accretion rates of approximately $(4-5)\times 10^{-5}\,M_{\sun}\,\mathrm{yr}^{-1}$, while the nonmagnetized models have rates of $\sim 1\times 10^{-5}\,M_{\sun}\,\mathrm{yr}^{-1}$.

Note that we follow only the first $\sim 10^{3}$ yr of evolution of the sink particles, and thus we cannot determine the final stellar masses in our models.

Figure 8 shows the ratio of the disk mass to the sink particle mass for each model. The strong magnetic field models ($\alpha=0.25$) exhibit a low ratio of $\sim 0.02-0.03$, indicating that a strong magnetic braking efficiently extracts angular momentum from the disks. The weak magnetic field models ($\alpha=0.1$) have a mass ratio of $\sim 0.1$, which is an order of magnitude larger than that of the strong field models. The nonmagnetized models ($\alpha=0$) have ratios that are on the order of unity, and such large ratios induce fragmentation of the disk (model M05B0).

4.3. Outflow formation

Figure 9 shows outflows at $t_{p}=700$ yr for the magnetized models. The strong magnetic field models (M05B025 and M1B025) have outflows that extend further than do those for the weak magnetic field models (M05B01 and M1B01). The envelopes for the weak magnetic field models are disturbed by the turbulence, as shown on the 800 au scale, while the strong magnetic field models have disk-shaped infalling envelopes. For all the models, the outflows extend in approximately the mean direction of the magnetic field on that scale. The direction of the magnetic field depends on its initial strength, and when the field is strong, the initial direction is approximately maintained ($z$-direction). In the weak magnetic field models, the magnetic field lines become steeply inclined with respect to the initial direction, and this has been reported by Matsumoto & Hanawa (2011).

The outflows are not always aligned with the disks, for which the orientations are shown in the lower panels of Figure 9. For models M05B025 and M1B025, the lower panels of Figure 9 are almost edge-on views of the disks, and the outflows extend nearly vertically in the upper and middle panels of Figure 9. This indicates that the outflows are roughly perpendicular to the disk axes for models M05B025 and M1B025. In Figure 9, the face-on disk of Model M05B01 is shown in the lower panel, and the outflow extends horizontally, as shown in the middle panel. This indicates that there is a misalignment between the disk and the outflow. On the other hand, model M1B01 produces a disk that is roughly aligned with the outflow. Note that, for each of the models, the infalling envelope on the $\gtrsim 100$ au scale is aligned with the outflow, as shown in Figure 9. The envelope is oriented perpendicular to the magnetic field on that scale, and the outflow extends along the magnetic field, as reported by Matsumoto & Tomisaka (2004).

Figure 10 shows the length of the outflow as a function of time following the formation of sink particles. The length is measured from the sink particle to the maximum distance of the outflow region, where the outflow region is defined by a volume of $v_{r}\geq 2c_{s}$. This plot indicates that the strong magnetic field models result in more rapid growth than do those with a weak magnetic field. The growth rates are $\sim 2\,\mathrm{km\,s}^{-1}$ and $\sim 1\,\mathrm{km\,s}^{-1}$ for the strong and weak magnetic field models, respectively. These rates are consistent with typical gas velocities in the outflow (Figure 11).

4.4. Cavity formation

One of the prominent features in the strong field models is a cavity structure in the envelope. Figure 11 shows the outflow and envelope structures for model M1B025 on the scales of 400 au and 100 au. On the 400 au scale, the flattened envelope is perpendicular to the outflow, which is associated with the helical magnetic field lines. On the 100 au scale, the flattened envelope has a cavity in which the magnetic field lines are straight. The bipolar outflow is not associated with the straight magnetic field lines of the cavity; instead, it is associated with the helical magnetic field lines, which thread the disk around the sink particles. The cavity is created beside the disk, as shown in Figure 12. The radius of the cavity increases with time, and it increases up to $\sim 50$ au by $t_{p}\simeq 10^{3}$ yr.

The cavity is caused by the magnetic pressure in the following way. The magnetic pressure is higher inside the cavity than it is in the other regions of the envelope, and thus it pushes the gas away from the flattened envelope. The gas accumulates on the rim of the cavity, which then has a higher density than the other regions, as shown in Figure 12. Inside the cavity, the gas moves outward at a velocity of $\sim 0.2-0.5\,\mathrm{km\,s}^{-1}$. The velocity at which the rim extends is typically $0.2\,\mathrm{km\,s}^{-1}\,(=50~{}\mathrm{au}/10^{3}~{}\mathrm{yr})$, which is roughly equal to the speed of sound.

Model M05B025 also results in cavities in the envelope, and the size of the cavities ($\sim 10$ au) is smaller than that for model M1B025 when compared at the same time ($t_{p}=10^{3}$ yr). In the upper right panel of Figure 3, cavities can be seen at both the upper and lower sides of the disk, but their density contrast is less than that for the cavities seen in model M1B025. The weak magnetic field models M05B01 and M1B01 do not produce cavities in the envelopes.

Similar cavity formation has been reported in recent MHD simulations (e.g., Zhao et al., 2011; Joos et al., 2012; Krasnopolsky et al., 2012; Machida et al., 2014). In these simulations, the cavities are formed whether or not sink particles are implemented, and whether or not magnetic diffusion is considered. Our simulations suggest that the formation of a cavity depends on the magnetic field strength. The cavities reproduced in these simulations likely correspond to the magnetic wall that has been predicted by theoretical studies (Li & McKee, 1996; Tassis & Mouschovias, 2005).

5. Discussion

5.1. Disk growth

Our simulations indicate that a cloud core with a stronger magnetic field produces a disk with a smaller radius. On the timescale of $t_{p}\sim 10^{3}$ yr, this is roughly consistent with long-term simulations of magnetized rotating cloud cores performed by Machida et al. (2011) (see also Machida & Hosokawa, 2013). However, their simulations tend to result in larger disk radii than those produced by the models considered here. In their simulations, the disk radii exceed 10 au by $t_{p}=10^{3}$ when $\mu=1$, and reach 20 au when $\mu=3$ (see Machida et al., 2011, Figure 9a) Meanwhile, the models here had a disk radius of $3-4$ au when $\mu=1.12$ and $10-20$ au when $\mu=2.81$.

The difference in the disk radii is due to the difference in the initial distribution of the angular momentum. Machida et al. (2011) assumed a uniform rotation, which provides a specific angular momentum distributed as $j\propto r^{2}$. On the other hand, the turbulence produces a velocity distribution of $\Delta v\propto r^{1/2}$, due to the Larson scaling relations assumed here (Larson, 1981), and the distribution of the specific angular momentum is expected to be $j\propto r^{3/2}$. The differences seen in the initial angular momentum distributions are therefore larger for larger radii. This has a greater effect on the disk radius in the later stages of the accretion phase, because the angular momentum of the infalling gas has a strong impact on disk evolution (e.g., Vorobyov et al., 2015). This suggests that, for the models here, the disk radius is expected to increase with time, e.g., $t_{p}=10^{5}$ yr, but it will still be smaller than that produced by a model with uniform rotation ($\sim 100$ au). Small disks with a radius of less than 100 au are expected to be produced by the magnetized turbulent model.

Recent high-resolution observations have revealed Keplerian disks around Class 0 and Class I protostars (e.g., Jørgensen et al., 2009; Tobin et al., 2012). Ohashi et al. (2014) recently suggested that the Class 0 protostar L1527 IRS has a Keplerian disk with a small radius of 54 au. A more extended disk has been suggested in the Class 1 source TMC 1A, which has an estimated radius of 100 au (Aso et al., 2015). Disks that have a radius that exceeds 100 au have also been suggested in the Class 0 source VLA1623A (Murillo et al., 2013) and the Class I source L1489 IRS (Hogerheijde, 2001; Yen et al., 2014). Such a variety of Keplerian disks may be responsible for the variety of magnetic fields and turbulence in natal cloud cores.

5.2. Alignment between an outflow, a disk, and an envelope

The various models examined here have shown that some disks are misaligned with the outflow, as described in Section 4.3. Matsumoto & Tomisaka (2004) investigated the direction of outflows when the initial magnetic field is misaligned with the rotation axis on the cloud core scale (see also Hennebelle & Ciardi, 2009; Ciardi & Hennebelle, 2010). Their simulations indicate that the outflows are extended in the direction parallel to the local magnetic field, even when this direction is not aligned with the magnetic field on the cloud core scale. This causes the outflow on the $\sim 100$ au scale to be misaligned with the magnetic field on the cloud core scale.

Figures 14 and 14 show the directions of the magnetic field, the rotation axes, the disk-like structures, and outflows for the representative models M05B01 and M1B025. Measurement of these directions is described in Appendix B.

For model M05B01, the mean magnetic field \overline{\mbox{\boldmathB}} on the scale of $10^{4}$ au is antiparallel to that on the 1 au scale, indicating that the magnetic field rotates up to $\sim 180\arcdeg$ around the vector \overline{\mbox{\boldmathj}}. The disk normal vector \overline{\mbox{\boldmathn}} is associated with \overline{\mbox{\boldmathB}} on scales larger than 10 au, indicating that the flattened envelope is perpendicular to the local magnetic field. On scales smaller than 10 au, the disk normal vector \overline{\mbox{\boldmathn}} is associated with the mean angular momentum \overline{\mbox{\boldmathj}}, indicating that the disk is aligned with the rotation axis. The outflow is accelerated on the 10 au scale, and it extends along the magnetic field, up to the 100 au scale. The outflow is therefore misaligned with the rotating disk, but it is aligned with the flattened envelope.

The ejection mechanism for this outflow is different from the ordinal magnetocentrifugal wind (Blandford & Payne, 1982; Pudritz & Norman, 1986); this outflow mechanism is a spiral flow (see Figure 15b). The model here demonstrates that the spiral flow reproduces a bipolar outflow. Similar spiral flows have been observed (see Figure 21 of Matsumoto & Hanawa, 2011). To confirm that the spiral flow mechanism continues to drive the outflow on the timescale of $10^{4}$ yr, further long-term simulations are necessary (e.g., Seifried et al., 2012, 2013).

For model M1B025, the outflow is driven by the magnetocentrifugal wind on the $\sim 100$ au scale (see Figure 15a), and the flow direction is aligned with the local magnetic field \overline{\mbox{\boldmathB}}, as shown in Figures 14 and 14. On this scale, the flattened envelope is also aligned with the magnetic field. On the scale of $\lesssim 10$ au, the disk normal vector \overline{\mbox{\boldmathn}} is aligned with the mean angular momentum \overline{\mbox{\boldmathj}}, indicating that the disk is aligned with the rotation axis. The flow direction on this scale is contaminated by the velocity associated with the formation of a cavity (see the right bottom panel of Figure 9).

The direction of the angular momentum vector depends on the radius of the cloud core, as shown in Figure 14. This indicates that the disk accretes lumps of gas for which the angular momentum vectors have highly nonuniform directions. The angular momentum of the infalling gas greatly affects the evolution of the disk, as shown in Vorobyov et al. (2015). As the disk further evolves, its orientation and size are both expected to change.

Observations of misalignment between outflows, magnetic field, circumstellar disks, and flattened envelopes have been reported. Hull et al. (2013) showed that outflows are misaligned or randomly aligned relative to the magnetic field on the $\sim 1000$ au scale for Class 0 and Class I objects, and this is consistent with the results of the models considered here, which predicted that the direction of the magnetic field depends on the scale length (see, e.g., Figure 14), though the outflow accelerates along the magnetic field on the $10-100$ au scale. If the outflow is aligned with the acceleration, it is expected to be misaligned with the magnetic field on the $\sim 1000$ au scale.

Observations of the Class I source L1489 IRS have suggested that the central Keplerian disk is inclined with respect to the flattened infalling envelope (Brinch et al., 2007a, b). This misalignment between the central disk and the flattened envelope was reproduced in all the models we considered, as shown in Figure 14. In each model, the disk normal vector \overline{\mbox{\boldmathn}} (the green line) drifts on the $\theta_{x}-\theta_{y}$ plane, indicating that with the flattened density structure, the inclination depends on the radius. Moreover, \overline{\mbox{\boldmathn}} is aligned with \overline{\mbox{\boldmathB}} when $r\gtrsim 100$ au, suggesting that the inclination of the flattened envelope is guided by the magnetic field. In other words, the flattened density structure of the envelope is caused by the magnetic field.

The Class I binary source in the Ophiuchus star-forming region, IRS 43, is a complex misaligned system. The most curious misalignment in this object is that of the circumbinary disk and the orbit of the binary. According to Brinch et al. (2016), IRS 43 has a circumbinary disk of which inclination is nearly edge-on, while the orbit of the binary is close to being in the plane of the sky. Such misalignment is possibly produced by a non-monotonic distribution of the angular momentum on the scale of the cloud core as shown in Figure 14. On the cloud scale, the infalling gas has misaligned angular momentum, leading the time-dependent accretion of angular momenta onto the circumbinary disk and circumstellar disks.

Such a misaligned system cannot be reproduced by axisymmetric models (e.g., Machida et al., 2008). Even if weak turbulence is assumed, misaligned systems are reproduced.

5.3. Cavity and arc-like structure

The rim of the cavity has a higher column density than does the envelope, as shown in Figure 12, and it may be observed as an arc-like structure. Tokuda et al. (2014) reported that the ALMA Cycle 0 observations reveal an arc-like structure at the center of the high-density molecular cloud core MC27 or L1521F. The arc-like structure was extended to a length of $\sim 1000$ au, and they proposed that it was caused by a dynamical interaction between the dense gas condensation and the envelope. Matsumoto et al. (2015b) performed hydrodynamical simulations to determine the origin of the arc-like structure, and they demonstrated that gravitational torque due to the orbiting protostars produces arc-like structures extending up to 1000 au.

The typical length of the cavity is consistent with the observations of MC27/L1521F. The rim of the cavity obtained with the model M1B025 was extended to $\sim 50$ au at $t_{p}=1000$ yr. If we assume that the cavity continues to extend at a constant velocity of $0.2\,\mathrm{km\,s}^{-1}$ (see Section 4.4), it takes $2\times 10^{4}$ yr for the rim to extended to 1000 au. This timescale agrees with that of the arc structures reproduced by Matsumoto et al. (2015b), and it is also consistent with the timescale of the protostar (Spizer source) in MC27/L1521F.

As shown in the literature (e.g., Zhao et al., 2011; Joos et al., 2012; Krasnopolsky et al., 2012; Machida et al., 2014), the detailed structure of the cavity and the rim seems to be sensitive to the simulation settings. Comparison between the models and the observations should be performed in terms of typical values, e.g., a typical length, as shown here.

The rim of the cavity can account for the dense gas condensations in MC27/L1521F. Our simulation results indicate that the number density of the rim is $\sim 10^{10}\,\mathrm{cm}^{-3}$ on the 50 au scale. When the cavity is expanded up to 1000 au, the number density of the rim is expected to be $2.5\times 10^{7}\,\mathrm{cm}^{-3}$, assuming that the density is distributed as $\rho\propto r^{-2}$. The number densities of the dense gas condensations, MMS-2 and MMS-3, are estimated to be $10^{7}\,\mathrm{cm}^{-3}$ and $10^{6-7}\,\mathrm{cm}^{-3}$, respectively (Tokuda et al., 2014). Therefore, MMS-2 and MMS-3 can be explained by the dense portions of the rim rather than by the fragments.

On the other hand, between the rim and the cavity, the column density differs by a factor of $\sim 30$; the column densities are $\sim 5\times 10^{24}\,\mathrm{cm}^{-3}$ on the rim and $\sim 2\times 10^{23}\,\mathrm{cm}^{-3}$ in the cavity (Figure 12). Such a high contrast in the column density has not been seen in observations of the ALMA Cycle 1 (see Figure 3b of Tokuda et al., 2016). Observations of the magnetic field of the cloud core will be of key importance in determining which model best accounts for the origin of the arc-like structure.

6. Summary

Gravitational collapse of molecular cloud cores and the formation of circumstellar disks and outflows around protostars were investigated by performing AMR simulations; the effects of both turbulence and the magnetic field were considered. Ohmic dissipation was considered in the MHD simulations. We allowed the system to evolve for $\sim 1000$ yr following the formation of a protostar. The main outcomes are summarized as follows.

In each of the magnetized models, the cloud core collapses to form a protostar surrounded by a circumstellar disk. The star-disk system is surrounded by an infalling envelope, and bipolar outflows are ejected. The nonmagnetized models produce massive circumstellar disks, one of which undergoes fragmentation at $\sim 10^{4}$ yr following the formation of a protostar. 2. 2.

The radius of the circumstellar disk depends on the initial strength of the magnetic field. Models with a stronger magnetic field produce a circumstellar disk with a smaller radius. The mass of the disk shows a similar dependence on the magnetic field, where a stronger field produces a less massive disk. The ratio of disk mass to stellar mass remains roughly constant at about $\sim 1-10$%, depending on the strength of the magnetic field. 3. 3.

The magnetized models reproduce the outflow, which can be classified into two types: a magnetocentrifugal wind and a spiral flow. In the latter, the outflow is not aligned with the rotational axis of the disk. In both the cases, the outflow and the flattened envelope are aligned with the magnetic field on that scale. In some models, the outflow is misaligned with the circumstellar disk. Similarly, the flattened envelope may be misaligned with the circumstellar disk. 4. 4.

The internal distribution of angular momentum in the cloud cores is nonuniform. After long-term evolution, the disk accretes lumps of gas in which the direction of the angular momentum vectors is highly nonuniform; hence, the disk is expected to change its orientation and size. This means that a planet formed during a later phase may have an orbital angular momentum that is highly misaligned with the angular momentum of the central star. 5. 5.

A strong magnetic field tends to produce a cavity in the infalling envelope; this is due to the strong magnetic pressure, and the gas accumulates on the rim. Thus, the rim can account for the arc-like structure and dense gas condensation observed in the high-density molecular cloud core MC27/L1521F, though it has not been verified by observation that there is a high contrast between the column density of the cavity and that of the rim.

We would like to thank K. Tokuda, T. Onishi, K. Tomida, R. Kawabe, and N. Ohashi for fruitful discussions. Numerical computations were carried out on the Cray XC30 and XT4 at the Center for Computational Astrophysics, National Astronomical Observatory of Japan, and the HITACHI HA8000 Clustre System (T2K-Todai) in the Information Technology Center, The University of Tokyo. This research was supported by JSPS KAKENHI Grant Numbers 16H02160, 15K05032, 26400233, 26287030, 25400232, 24244017, 23540270, and 23244027.

Appendix A Measurement of the Disk Radius

The disk radius for each model was estimated using the density and velocity distributions, as follows. The volume of the disk $V_{d}$ was defined by two criteria: $\rho\geq\rho_{\mathrm{cr}}$ and $(v_{\varphi}^{2}+v_{\theta}^{2})^{1/2}/(v_{r}^{2}+c_{s}^{2})^{1/2}\geq 3$. The former indicates that the disk has a density higher than the critical density of the equation of state. The latter indicates that the velocity of rotation is greater than the radial velocity. In spherical coordinates, the velocity is $(v_{r},v_{\theta},v_{\varphi})$, which was calculated by a transformation from Cartesian coordinates with the origin set at the location of the sink particle. Because the disk orientation is not aligned to any of the coordinate axes, the tangential velocity of $(v_{\varphi}^{2}+v_{\theta}^{2})^{1/2}$ was adopted as the rotational velocity. The value of three on the right-hand side of the second criterion was determined empirically.

The disk radius was obtained from the inertia tensor of the volume, $V_{d}$. The inertia tensor is calculated as

[TABLE]

and each element of the matrix is calculated by a moment of the coordinates, e.g.,

[TABLE]

where $(x,y,z)$ are the coordinates of a cell, and $(x_{p},y_{p},z_{p})$ is the position vector for a sink particle. The volume integrals in equation (A2) were performed by summing over the cells within the volume, $V_{d}$. The matrix $I$ yields three eigenvalues, $\lambda_{1}>\lambda_{2}>\lambda_{3}$, and the square root of each of the eigenvalues corresponds to the length of a principal axis; $\lambda_{1}^{1/2}$, $\lambda_{2}^{1/2}$, and $\lambda_{3}^{1/2}$ correspond to the semi-major axis, the semi-minor axis, and the thickness of the disk, respectively. The semi-minor axis $\lambda_{2}^{1/2}$ was adopted as the disk radius so that this method could be applied to a highly elongated disk, such as the disk of model M1B025 (the lower right panel of Figure 3). The disk radius is defined as

[TABLE]

where the factor of two comes from the inertial tensor for a uniform thin disk with a radius of $R_{d}$, i.e., $I_{xx}+I_{yy}=(1/2)R_{d}^{2}$.

We also calculated the mass of the disks using

[TABLE]

where $V_{d}$ is the disk volume.

Appendix B Measurement of the direction of the axes

In Figures 14 and 14, the direction of the magnetic field, angular momentum, disk, and outflow are each shown as a function of the radius. These axes were measured as follows.

In order to define the direction of the magnetic field and angular momentum, we measured the mean magnetic field and the angular momentum:

[TABLE]

where

[TABLE]

The volume $V_{s}(r)$ is that of a sphere with radius $r$, and the center coincides with the position of the sink particle:

[TABLE]

The vectors \mbox{\boldmathr}_{p} and \mbox{\boldmathv}_{p}, respectively, denote the position and velocity of the sink particle.

The orientation of the disk was calculated using the eigenvector of the inertia tensor, which is similar to $I$ in equation (A1). The elements of the inertia tensor are obtained as follows.

[TABLE]

where the integration is performed inside the region in which $\rho\geq\rho_{d}$, for a given threshold $\rho_{d}$. The eigenvector associated with the smallest eigenvalue represents a normal vector for the flattened disk-like structure. The radius of the disk-like structure is defined as the maximum extent of the region in which $\rho\geq\rho_{d}$, as measured from the position of the sink particle. Thus, we obtain the direction of the disk normal vector \overline{\mbox{\boldmathn}} as a function of the radius $r$. Note that the disk measured here corresponds to the circumstellar disk for a small radius, e.g., $r\lesssim 10$ au, and to the flattened infalling envelope for a large radius.

The direction of the outflow was obtained from the flow of the gas inside the outflow region. Because the outflow is bipolar and the gas flows roughly parallel to the magnetic field, the direction of the outflow is calculated as follows.

[TABLE]

where $V_{\mathrm{of}}(r)$ denotes the region of the outflow, defined as

[TABLE]