Example of Numerical Simulation

Figure 4.15: Explanation of nested grid method.

We have described the cloud run-away collapse and succeeding accretion process In the latter phase, there is a possibility that the outflow is driven by the effect of magnetic Lorentz force. A restricted number of numerical simulation can simulate such evolution throughout. The size of the molecular cores $\sim 0.1{\rm pc}$ is more than $3\times 10^4$ times as large as the typical size of the first core $\sim 1{\rm AU}$. Therefore numerical scheme to calculate such process must have a large dynamic-range. In the finite-difference scheme, the spatial dynamic range is restricted by the number of the cells. Although at least $3\times 10^5 \times 3\times 10^5$ grid points (in 2D) seems necessary to resolve a factor of $3\times 10^4$, to solve the MHD equation with the Poisson equation of the self-gravity have not been done. This is done by the Eularean nested-grid simulation. This method uses a finite number of grid systems with different spatial resolutions are prepared. Coarse grid covers whole the cloud and is to see a global structure. Fine grid covers only the central part of the cloud and is to see the fine-scale structure appeared near the center. Figure 4.15 explains grid cells of the nested grid method. The size of the cells of the $n$-th level grid (L$n$) is taken equal to a half of those of the L$n-1$.

Figure 4.16: $\Omega_c={10^{-13} \rm rad\,s^{-1}}$, $B_c=10\mu{\rm G}$, $\rho_c = 10^4 {\rm H_2~cm^{-3}}$, and $c_s = 190{\rm m\,s^{-1}}(T/10{\rm K})$. Snapshots at the time of $0.6066\tau_{\rm ff}=1.061{\rm Myr}$ represented in different levels are shown: L1 (upper-left), L5 (upper-right), and L10 (bottom). The actual size of the frames of L5 (upper-right) and L10 (bottom) are, respectively, 1/16 and 1/512 smaller than that of L1 (upper-left). Magnetic field lines (dotted lines) and isodensity contours (solid lines) are presented as well as the velocity vectors by arrows.

A cylindrical cloud with coaxial magnetic field in hydrostatic balance is studied. Rotation vector $\mbox{\boldmath${\omega}$}$ and the magnetic field $\mbox{\boldmath${B}$}$ are parallel to the symmetric axis of the cylinder ($z$-axis). For the energy equation, a double-polytropic equation of state is adopted, since the radiation hydrodynamical calculation (for example Masunaga & Inutsuka 2000) predicts $\rho-T$ relations like shown in Figure 4.13, which consists of a number of power-laws for some density ranges.

$$p=\left\{\begin{array}{l} c_s^2\rho\ \ \ ({\rm for}\ \ \rho<\rho_... ...{\rm A}}\right)^{7/5}\ \ \ ({\rm for}\ \ \rho>\rho_{\rm A}), \end{array}\right.$$ (4.102)

where we take $\rho_A=10^{10}{\rm H_2\,cm^{-3}}$. Figure 4.16 is the model of $\Omega _c={10^{-13} \rm rad\,s^{-1}}$, $B_c=10\mu {\rm G}$, $\rho_c = 10^4 {\rm H_2~cm^{-3}}$, $c_s = 190{\rm m\,s^{-1}}(T/10{\rm K})$, and scale-length $H=0.341 {\rm pc} (c_s / 190{\rm m\,s^{-1}}) (\rho_c / 10^4 {\rm H_2~cm^{-3}})^{-1/2}$. This is a snapshot of $0.6066\tau _{\rm ff}=1.061{\rm Myr}$. The cylindrical cloud breaks into pieces with a length in $z$-direction equal to the most unstable wavelength of the gravitational instability (Matsumoto, Nakamura, & Hanawa 1994) as

$$\lambda_{\rm max} \simeq \frac{2\pi (1+\alpha/2+\beta)^{1/2}... ...ht)^{1/3}-0.6\right]^{1/2}}\frac{c_s}{(4 \pi G \rho_c)^{1/2}}.$$ (4.103)

The magnetic fields prevent gas motions perpendicular to the field lines $\mbox{\boldmath${B}$}$. Further, the centrifugal force works in the perpendicular direction to $\mbox{\boldmath${\omega}$}$. A disk forms which running in the perpendicular direction to $\mbox{\boldmath${B}$}$ and $\mbox{\boldmath${\omega}$}$ and thus to the $z$-axis. Figure 4.16 illustrates a snapshot at the time of $0.6066\tau _{\rm ff}=1.061{\rm Myr}$. Level 1 grid (upper-left) captures the global structure of the contracting disk perpendicular to the $z$-axis. Close up view is shown in L5, which has 16 times finer spatial resolution (upper right). Since speed, $v_z$, of infalling gas is larger than the sound speed, a number of shock fronts are formed near $z\simeq 0.02H$ and $z \mbox{\raisebox{0.3ex}{$<$}\hspace{-1.1em} \raisebox{-0.7ex}{$\sim$}} 0.01H$. Further inner structure is represented in L10 grid (bottom) which has 32 times finer resolution than L5. Another discontinuity in density is now forming near $z\sim 0.0002H$. Since gas inside of this discontinuity has the density larger than the critical density $\rho > \rho_{\rm A}$, in this phase the first core (p.) begins to be formed. In the isothermal collapse phase, a characteristic power-law for the Larson-Penston self-similar solution is seen in the radial distributions of density as $\rho(r,z=0) \propto r^{-2}$, and $B_z(r,z=0)\propto r^{-1}$.

After $t > 0.6067\tau_{\rm ff}$ the first core is formed. Radial infall motion is accelerated toward the surface of the first core and the rotational motion $v_\phi$ also increases. This amplification promotes the toroidal magnetic field, $B_\phi$. Since the torque is proportional to $(\mbox{\boldmath${\nabla \times}$}B_\phi\mbox{\boldmath${e}$}_\phi) \mbox{\boldmath${\times B}$}_p$, the amplitude of the torque increases after the first core formation.

Structure just after $\sim$500 yr has passed from Figure 4.16 is shown in Figure 4.16 In $\tau \sim 500$ yr from the epoch of Figure 4.16, flow is completely changed and outflow is launched. Panel (b), which corresponds to a four-times close-up of panel (a), shows clearly that outflow is ejected from a region near $r \simeq 1.5\times 10^{-4}H$ and $z\simeq 0.5\times 10^{-4}H$. And the outflow vectors and the disk has an gngle of approximately 45 deg. The magneto-centrifugal wind model favors a small angle between the magnetic field and the disk, $\theta_{\rm mag}< 60$deg. This model predicts that the angle between the flow and the disk is also smaller than 60 deg. Numerical results is not inconsistent with this prediction of magneto-centrifugal wind model.

To see which force is driving such outflow, amplitudes of three forces are compared: the Lorentz force $\mbox{\boldmath${F}$}_m=(\mbox{\boldmath${\nabla\times B}$})\mbox{\boldmath${\times B}$}/4\pi$, the centrifugal force $\mbox{\boldmath${F}$}_c=(\rho v_\phi^2/r)\mbox{\boldmath${e}$}_r$, and the thermal pressure gradient $\mbox{\boldmath${F}$}_t = -\nabla p$. The parallel components to the poloidal magnetic field of the three forces are as

$$\mbox{\boldmath${F}$}_m\cdot \frac{\mbox{\boldmath${B}$}_p}{\vert\mbox{\boldmath${B}$}_p\vert} = \frac{\mbox{\boldmath${\nabla\times B}$})\mbox{\boldmath${\times ... ...}{4\pi}\cdot \frac{\mbox{\boldmath${B}$}_p}{\vert\mbox{\boldmath${B}$}_p\vert},$$

$$= -\frac{1}{8\pi r^2}\frac{\mbox{\boldmath${B}$}_p}{\mbox{\boldmath${B}$}_p\vert}\cdot \nabla(rB_\phi)^2$$ (4.104)

(Ustyugova et al. 1999)

$$\mbox{\boldmath${F}$}_c\cdot \frac{\mbox{\boldmath${B}$}_p}{... ...\phi^2}{r}\cdot \frac{B_r}{\vert\mbox{\boldmath${B}$}_p\vert},$$ (4.105)

and

$$\mbox{\boldmath${F}$}_p\cdot \frac{\mbox{\boldmath${B}$}_p}{... ...c{\mbox{\boldmath${B}$}_p}{\vert\mbox{\boldmath${B}$}_p\vert}.$$ (4.106)

Figure 4.17 illustrates the largest force at each grid point. Centrifugal force dominated region (region C), which is shown by asterisk (*), overlaps the outflow region around point P1 ($z=10^{-4}H$, $2.5\times 10^{-4}H$). Just radially exterior to this region C, there is another region near point P2 ( $z=2\times 10^{-4}H$, $3.5\times 10^{-4}H$) filled with plus signs where the magnetic force is dominated (region M). Comparing with the flow vectors in the left, the strong outflow coincides with this region M and the above region C. Therefore, the outflow is driven by the centrifugal force and the Lorentz force (pressure gradient by toroidal magnetic fields).

Difference between models with different magnetic field strength but the same rotation speed is illustrated in Figure 4.18. Magnetic to thermal pressure ratio $\alpha$ is taken 1 [panels (a), and (d)], 0.1 [panels (b), and (e)], and 0.01 [panels (c), and (f)]. Panels (a), (b), (c) in the left column are snapshots when the central density reaches $\rho_{\rm A}$ while panels (d), (e), (f) in the right column show the structure after $\tau=4.5\times 10^{-3}\tau_{\rm ff}=8000{\rm yr}$ passed. In the model with $\alpha = 1$, the outflow gas flows through a region whose shape resembles a capital letter U. This is similar to that observed in model ($\alpha = 1$ and $\Omega=5$) shown in Figure 4.17. Decreasing the magnetic field strength $\alpha=1\rightarrow 0.1$ and $\rightarrow 0.1$, the shape of disk changes from a flat disk to a round ellipsoid. The model with neither magnetic field nor rotation shows a spherical contracting core. Outflow is also affected by decrease in magnetic field strength. Globally folded magnetic field lines are seen in the outflow region in panel (e). Smaller-scale structure folding the magnetic field lines develops in the outflow region in panel (f). The strucure seen in the polidal magnetic field is well correlated to the distribution of toroidal magnetic field strength. That is, a packet of strong $B_\phi$ pinch the gas and poloidal magnetic field lines radially inwardly. This means that the toroidal magnetic field is amplified by the rotation motion and becomes predominant over the poloidal one in the case with meak initial magnetic field. Thus, the whoop stress, which is expressed by the first term of the Lorentz force of equation(B.13) and comes from the tension of toroidal magnetic field, pinches the plasma contained inside of the loop of $B_\phi$. For this case, the magnetic force dominant region is widely spread compared with the centrifugal force dominant region, which is far from the model with strong magnetic field. In this case, the magnetic pressure gradient plays an important role. Gas of outflow is ejected perpendicularly from the disk and forms flow similar to the capital letter I. In conclusion, U-type and I-type outflows are generated for strong magnetic model and weak magnetic model, respectively.

Figure 4.17: Continuation of the evolution shown in Fig.4.16. Snapshots at the time of $0.6069\tau_{\rm ff}=1.066{\rm Myr}$ just $\sim 500$ yr after from Fig.4.16. Panel (a) plots the structure of L10 which is similar to Fig.4.16 (c). Spherical region inside $\sim 1.2\times 10^{-3} H$ has been swept by the outflow. Panel (b) is a snapshot of L12 grid which has 4-times finer spatial resolution than L10 (a). In panel (c), we plotted magnetic force dominated region with `+', centrifugal force dominated with `*', and thermal pressure force dominated with blank ' '.

Figure 4.18: Difference in models with the same rotation speed $\Omega=2\times 10^{-14} \rm rad\,s^{-1}$ but different magnetic field strength ((a) and (d): $B_c=10\mu{\rm G}$, (b) and (e): $3\mu{\rm G}$, and (c) and (f): $1\mu{\rm G}$).