“Ash-fall” induced by molecular outflow in protostar evolution

We solve two-fluid magneto-hydrodynamics (MHD) equations for dust-gas mixture. The detail of the governing equations and numerical scheme are described in our previous paper (Tsukamoto et al., 2021). In Tsukamoto et al. (2021), we also showed that treating charged dust and neutral dust as a single fluid is justified for investigating the dynamics of the dust with size of $a_{d}\gtrsim 10{\rm\mu}{\rm m}$, which is the focus of this paper.

In addition to the governing equations presented in Tsukamoto et al. (2021), we further solve the equation for the dust-size evolution with single-size approximation (e.g., Sato et al., 2016; Tsukamoto et al., 2017) as follows. The equation for the time evolution of the dust-size $a_{d}$ is given by

where $t_{\rm growth}=1/(\pi a_{d}^{2}n_{d}\Delta v_{\rm dust})$, $n_{d}$ is the dust number density, $\Delta v_{\rm dust}$ is the collision velocity between the dust grains, $D_{d}/Dt=\partial/\partial t+v_{d}\cdot\nabla$, $v_{d}$ is the (mean) dust velocity, and

is a factor to model the collisional mass gain and loss (Okuzumi et al., 2016). The dust gains and loses a mass when $\Delta v_{\rm dust}<\Delta v_{\rm frag}$ and $\Delta v_{\rm dust}>\Delta v_{\rm frag}$, respectively, through dust collisions. To evaluate the dust relative velocity between the dust $\Delta v_{\rm dust}$, we consider the sub-grid scale turbulence and Brownian motion, and relative velocity is given as $\Delta v_{\rm dust}=\Delta v_{\rm turb}+\Delta v_{B}$ (for details, see Appendix A).

We model equation 2 so that it reproduces the results of the collision simulations of dust aggregates (Wada et al., 2013). In this study we set $\Delta v_{\rm frag}=30~{}{\rm m}~{}{\rm s}^{-1}$. We can rewrite $D_{d}/Dt$ as

where $\mathbf{v}$ is the barycentric velocity of the dust-gas mixture, $\rho_{g}$ is the gas density, $\mathbf{v}_{g}$ is the gas velocity, $\rho=\rho_{g}+\rho_{d}$ is the total density, $\rho_{d}$ is the dust density, $\epsilon=\rho_{d}/\rho$, and $\Delta\mathbf{v}=\mathbf{v}_{d}-\mathbf{v}_{g}$. Using this relation, we can rewrite the equation (1),

The SPH discretization form of this equation is given as

where the last term is an artificial viscosity and $v_{\rm sig,\Delta\mathbf{v}}=\frac{1}{2}(c_{s,i}+c_{s,j})+3\beta_{\rm visc}{\rm min}(0,(\Delta\mathbf{v}_{i}-\Delta\mathbf{v}_{j})\cdot\mathbf{e}_{ij})$ is the signal velocity where $c_{s,i}$ is the sound velocity of $i$ th particle and $\mathbf{e}_{ij}=(\mathbf{r}_{i}-\mathbf{r}_{j})/|\mathbf{r}_{i}-\mathbf{r}_{j}|$. We choose $\beta_{\rm visc}=2$.

The timestep for dust advection is chosen to be

where $h_{i}$ is the smoothing length of $i$ th particle.

During the simulations, $\epsilon_{i}$ of a few particles are found to become negative (in particular in the outflow regions) with $|\epsilon_{i}|\lesssim 10^{-7}$. When it happens, we correct $\epsilon_{i}$ to be $|\epsilon_{i}|$ and maintain it positive. We confirm that the error of the total dust mass introduced with this correction is negligible because the absolute value of the negative $\epsilon_{i}$, as well as the number of the particles with $\epsilon_{i}<0$, is very small.

Our numerical simulations consider non-ideal MHD effects, including the Ohmic and ambipolar diffusions, but ignore the Hall effect. For the resistivity model, we adopt the resistivity table with $a_{d}=0.035{\rm\mu}{\rm m}$ presented in Tsukamoto et al. (2020). Thus, the dust size for dust dynamics and that for the resistivity table are not consistent with each other in our present simulations.

A sink particle was dynamically introduced when the density exceeds $\rho_{\rm sink}=10^{-12}~{}{\rm g~{}cm}^{-3}$. The sink particle absorbs SPH particles with $\rho>\rho_{\rm sink}$ within $r_{\rm sink}<1$ AU.

The initial condition of the cloud core is as follows. We adopt the density-enhanced Bonnor-Ebert sphere surrounded by medium with a steep density profile of $\rho\propto r^{-4}$ as the initial condition, which is

with

where $\xi_{\rm BE}$ is a non-dimensional density profile of the critical Bonnor-Ebert sphere. The steep density profile in $r>R_{c}$ is introduced to reduce the particle number and unnecessary computational costs. $f$ is a numerical factor related to the strength of the gravity, and $R_{c}=6.45a$ is the radius of the cloud core. In this study, we adopt $\rho_{0}=7.3\times 10^{-18}~{}{\rm g~{}cm}^{-3}$, $\rho_{0}/\rho(R_{c})=14$, and $f=2.1$. Then, the radius of the core is $R_{c}=4.8\times 10^{3}$ AU and the enclosed mass within $R_{c}$ is $M_{c}=1\thinspace M_{\odot}$. We adopt an angular velocity profile of $\Omega(d)=\Omega_{0}/(\exp(10(d/(1.5R_{c}))-1)+1)$ with $d=\sqrt{x^{2}+y^{2}}$ and $\Omega_{0}=2.3\times 10^{-13}{\rm s^{-1}}$. We adopt a constant magnetic field $(B_{x},B_{y},B_{z})=(0,0,83\mu G)$. The parameter $\alpha_{\rm therm}$ ($\equiv E_{\rm therm}/E_{\rm grav}$) is $0.4$, where $E_{\rm therm}$ and $E_{\rm grav}$ are the thermal and gravitational energies of the central core (without surrounding medium), respectively. The ratio of the rotational to gravitational energies $\beta_{\rm rot}$ within the core is $\beta_{\rm rot}$ ($\equiv E_{\rm rot}/E_{\rm grav}$) $=0.03$, where $E_{\rm rot}$ is the rotational energy of the core. The mass-to-flux ratio of the core is $\mu/\mu_{\rm crit}=3$. We resolve 1 $\thinspace M_{\odot}$ with $3\times 10^{6}$ SPH particles. We adopt a dust density profile of $\rho_{d}(r)=f_{dg}\rho_{g}(r)/(\exp(10(r/(1.5R_{c}))-1)+1)$ where $f_{dg}=10^{-2}$ is the dust-to-gas mass ratio. The dust density profile has the same shape with the gas density profile in $r\lesssim 1.5R_{c}$ but is truncated at $r\geq 1.5R_{c}$. The initial dust size is $a_{d}=0.1{\rm\mu}{\rm m}$.

We simulated the evolution of a protostar, a protoplanetary disk, and a outflow until $1.2\times 10^{4}$ years after the epoch of protostar formation. Figure 1 shows time evolution at four epochs of the gas density, dust density, dust size, and dust-to-gas mass ratio on the $x$-$z$ plane where the $x$- and $z$-axes are perpendicular and parallel to the rotation axis of parent core, respectively.

At the epoch $t_{*}=3\times 10^{3}$ yr (panels 1-a to 1-d in Figure 1) where $t_{*}=0$ corresponds to the formation epoch of protostar, the bipolar outflow has been formed and gas and dust are blown out from the central region. At this stage, however, dust has not grown significantly with maximum size of $\lesssim 10{\rm\mu}{\rm m}$, which is found at the centre of the system. The grown dust is distributed only in vicinity of the disk (panel 1-c). Given the small maximum size of the dust, the dust and gas should be still fully coupled, which is supported by the facts that their density structures are identical (panels 1-a and 1-b) and that the dust-to-gas mass ratio remains at the initial value of $0.01$ (panel 1-d).

At $t_{*}=5\times 10^{3}$ yr (panels 2-a to 2-d), the maximum size of the dust reaches $\sim 100{\rm\mu}{\rm m}$ at the centre, and dust with size of $a_{d}\sim 10{\rm\mu}{\rm m}$ begins to be entrained by the outflow (panel 2-c). Panels 2-a, b and d indicate that the dust and gas begin to decouple in the outflow and dust spreads out of the outflow. As a result, the dust-to-gas mass ratio begins to decrease inside the outflow (panel 2-d).

At $t_{*}=7\times 10^{3}$ yr (panels 3-a to 3-d), the dust with a size of $a_{d}\sim 100{\rm\mu}{\rm m}$ is distributed over a large area of a few 100 AU across. The stream-lines and arrows indicate that the dust moves parallel to the $x$-axis inside the outflow while the gas moves parallel to the $z$-axis (panels 3-a, 3-b and 3-d). This indicates that the dust decoupled from the gas is ejected from the rotating outflow due to centrifugal force. Furthermore, some dust stream-lines do not diverge but take a closed path, meaning that some of large dust is refluxed from the outflow to the envelope. The ejection further decreases the amount of the dust in the outflow. In addition, the ejected dust gathers around the outflow, forming a U-shaped dust enhanced region, which is clearly visible in the dust-gas mass ratio map (panel 3-d).

At $t_{*}=1.1\times 10^{4}$ yr (panels 4-a to 4-d), the dust with a size of $a_{d}\gtrsim 1{\rm mm}$ is distributed over a large area of several 100 AU across (panel 4-c). The stream-lines and velocity field of the dust clearly indicate that the dust moves from the outflow back to the envelope (panels 4-a, b, d), while the stream lines and velocity field of gas shows that the gas is released to the interstellar medium. At this epoch, the gas and dust density distribution in the outflow are markedly different. This is caused by the dust-gas decoupling in the outflow. The majority of the large dust in the outflow have been refluxed to the envelope.

Figure 2 shows three-dimensional stream-lines of the gas and dust at $t_{*}=1.1\times 10^{4}$ yr. The gas stream-lines show that the gas is rotating and outflowing from the surface of the disk and is blown away. By contrast, the dust stream-lines indicate different dusty dynamics from that of the gas. The dust also outflows from the inner part of the disk. Then it immediately decouples from the gas in the outflow. Since the outflow rotates as shown in the stream line of gas, the centrifugal force drives the dust towards the radial direction. Gravity then pulls it back towards the disk and the grown dust accretes to the edge of the disk.

The overall process is schematically shown in the figure 3, and is very similar to ash-fall from volcanic eruptions where dust-gas mixture is ejected from the crater and then the gas and dust decouple in the atmosphere, causing the dust (or ash) to fall selectively. Hence we name this phenomenon as “ash-fall” phenomenon in protostar evolution.

The reflux of the large dust due to the ash-fall phenomenon increases the amount of the dust in the accretion flow from the envelope and thus the dust-to-gas mass ratio in the disk. Figure 4 shows the time evolution of the mean dust-to-gas mass ratio of the disk $M_{d,disk}/M_{g,disk}$ and the mean dust size $\int\rho a_{d}dV/\int\rho dV$ in the disk, where the integration is performed within the disk. Here, we define the disk as a region with the gas density $\rho_{g}>10^{-13}~{}{\rm g~{}cm}^{-3}$. We found that, in an early phase of $t_{*}<5\times 10^{3}$ yr, where the mean dust size is $a_{d}\lesssim 10^{2}{\rm\mu}{\rm m}$, the dust-to-gas mass ratio stays constant at $M_{\rm dust}/M_{\rm gas}=10^{-2}$, maintaining the initial value. This fact implies that the amount of the refluxed dust is negligibly small, presumably because the outflowing dust is well coupled with gas and is not ejected from the outflow.

Once dust has grown to a size of $a_{d}\gtrsim 10^{2}{\rm\mu}{\rm m}$ at $t_{*}\sim 4\times 10^{3}$ yr, the dust-to-gas mass ratio begins to increase as a result of the increase of the amount of the dust in the accretion flow from the envelope due to the dust reflux. The dust-to-gas mass ratio keep increasing and becomes $\epsilon\sim 0.01$ at the end epoch of the simulation ($t_{*}\sim 1.2\times 10^{4}$ yr), implying that $\sim 10$ % of the dust in the disk is the refluxed large dust.

Note that the increase of dust-to-gas mass ratio in the disk is a consequence of the dust growth in the disk and the dust reflux, and hence, in contrast to Lebreuilly et al. (2020), is not a consequence of the dust drift in the prestellar collapse phase. We start from the dust size with $0.1{\rm\mu}{\rm m}$ and, with this dust size, the dust and gas are essentially completely coupled and maintain the initial dust-to-gas mass ratio of $10^{-2}$ during prestellar collapse phase (see figure 1-d and 4 at $t_{*}=0$).

Note also that $a_{d}$ obtained in the single-size approximation represents the peak-mass dust size (Sato et al., 2016) and the contribution of smaller dust on the dust density may be minor (at least with MRN like size distribution). Therefore, we expect $\sim 10$% increase in dust-to-gas mass ratio in the disk shown in figure 4 may be realized even by considering a more realistic dust size distribution. Future studies with dust size distribution are imperative to confirm this prediction.

Figure 5 shows the gas density, dust density, and dust-gas ratio for a 1500 AU scale at $t_{*}=1.1\times 10^{4}$ yr. The gas density map shows that the gas is relatively dense inside the outflow while the dust density is significantly depleted inside the outflow and is enriched around the outflow, forming a U-shaped dust enhancement. In the dust-enhanced region, the dust-to-gas mass ratio is as high as 0.015, which is roughly a 50% increase from that of the typical interstellar medium.

A similar U-shaped dust-enhanced region has been observed in some Class 0/I YSOs, where the dust thermal emission was not detected, or much weaker than the surrounding regions, inside the outflow (such as L1527, HOPS136, B335, IRAS 15398-3559 Yen et al., 2010; Fischer et al., 2014; Yen et al., 2017; Maury et al., 2018; Harris et al., 2018). The reason of lack of detection of the dust thermal emission inside the outflow while it was detected in the surrounding regions, has been thought to be the low gas density in the outflow. However, our simulation suggests that dust cavities are formed due to the ejection of large dust from the outflow by centrifugal force. Then, it implies that the U-shaped dust enhancement is an indicator of growth of dust in the disk and outflow.

In this paper, we propose a new type of dust dynamics i.e., the “ash-fall” phenomenon in protostar evolution, in which the grown dust in the disk is blown away by outflows and is then refluxed on to the outer edge of the disk. The process consists of the following four steps (see also the schematic drawing in figure 3):

1. 1.

   dust growth in the disk where the density is sufficiently high and the dust growth timescale $t_{\rm growth}$ is as short as $t_{\rm growth}\sim 10^{3}$ yr (see Appendix B and figure 4),

2. 2.

   large dust with $a_{d}>1{\rm mm}$ can be entrained by the powerful magnetized outflow launched from the inner region of the disk,

3. 3.

   decoupling of the gas and dust in the outflow and dust ejection from the outflow due to centrifugal force, leading U-shaped dust-enhanced region around the outflow,

4. 4.

   mixing of the grown dust into the envelope accretion flow and the reflux of the large dust to the disk outer edge.

The ash-fall supplies a non-negligible amount of large dust to the outer-edge of the disk even at shortly after protostar formation ($t\sim 10^{4}$ yr) and keep increasing ($\sim 10$ % of the total amount of the accreting dust from the envelope; see figure 4).

The outer edge of a disk has a lower surface density than its inner region and the Stokes number $~{}{\rm St}$ is inversely proportional to the gas surface density $\Sigma_{g}$, given by $~{}{\rm St}=\pi\rho_{\rm mat}a_{d}/(2\Sigma_{g})$, where $\rho_{\rm mat}$ is the dust internal density. Therefore, the Stokes number of dust increases when the dust is refluxed to the disk outer edge. As a result, the dust of $~{}{\rm St}<1$ at the inner region of the disk can have $~{}{\rm St}>1$ when it is refluxed to the disk edge.

More quantitatively, we make the following estimates. We assume the disk has an exponential tail (as suggested from observations Williams & Cieza, 2011) and connected to envelope with accretion shock (e.g., Miura et al., 2017). The envelope surface density at the connection with the edge of the disk (at the upstream of accretion shock) is estimated as $\Sigma_{g}=\dot{M}_{\rm gas}/(2\pi rv_{r})$ where $v_{r}=\sqrt{2GM_{\rm star}/r}$. Then the disk surface density at the edge of the disk (at the downstream of accretion shock) is roughly given as

by assuming the shock is isothermal. Hence the Stokes number at the disk outer edge is

where $\dot{M}_{gas}$ is the mass accretion rate. Note that $~{}{\rm St}$ can be even larger by a factor of $(v_{r}/c_{s})^{2}(\sim 40$ with our parameters) when the accretion shock is adiabatic and density increase by the shock compression is small. Note also that the downstream density also decreases when the shock is oblique. This estimate suggests that the dust of $a_{d}\gtrsim 1{\rm cm}$ can be already $~{}{\rm St}\gtrsim 1$ when it is refluxed to disk edge.

As such, the ash-fall phenomenon can provide a pathway to circumvent the radial drift barrier and may allow dust to grow to planetesimals in the outer region of disk. Thus, the large dust could be the embryo of the planetesimal at the outer region of the disk.

Note that the radial drift barrier is conventionally regarded as a problem in Class II/III YSOs. Actually, the radial drift universally occurs even in Class 0/I YSOs, once the rotationally supported protostellar disk is created and dust growth starts. For example, Tsukamoto et al. (2017) showed that the dust growth and radial drift also happens in a disk of Class 0/I YSOs. In this sense, the occurrence of the radial drift barrier is not limited to Class II/III YSOs.

On the other hand, some molecular outflows seem to last longer and remain even in Class II YSOs (e.g., Hartigan et al., 1995). Thus, although our simulation stops at $t\sim 10^{4}$ yr after protostar formation (i.e., investigating the evolution of Class 0/I YSOs), we can expect ash-fall mechanism also plays a role even in Class II YSOs (or possibly more important because Stokes number of refluxed dust at the disk edge increases as envelope accretion diminishes; see, equation (4.2)).

In addition, a significant amount of large dust whose Stokes number is close to unity refluxed to disk provides an ideal condition for the growth of the secular gravitational instability (SGI) (Takahashi & Inutsuka, 2014, 2016; Tominaga et al., 2018, 2020) that are expected to be a powerful mechanism for creating ring and planetesimals in the outer region of the disk. Since SGI creates many axisymmetric dusty ring structures in the disk and may subsequently create planetesimals, a combination of the ash-fall and subsequent SGI can be an explanation for multiple ring-like structure observed in protoplanetary disks (ALMA Partnership et al., 2015; Fedele et al., 2017; Andrews et al., 2018; Fedele et al., 2018; Long et al., 2018). It is interesting to note that the very first example of those disks, HL-Tau, is known to possess bipolar molecular outflow and collimated jet-like structures. The co-existence of outflows/jets and very ordered axisymmetric multiple ring-like structure in the disk has been puzzling us since its discovery in 2014. Now we are proposing that actually these two processes might be cooperating in such an object.

Another important implication of ash-fall is that the reflux of the large dust into the envelope of $\sim 1000$ AU may explain the presence of grown dust in the envelope suggested by the observations. As shown in, e.g., Kwon et al. (2009) and Galametz et al. (2019), the spectral index of dust opacity $\beta$ in the envelope of some Class 0 YSOs is lower than that of the ISM. We expect that this is explained by large dust supplied from the disk to the envelope by the outflow.

Recently, Lebreuilly et al. (2020) conducted the non-ideal MHD simulation with dust fluids including dust-size distribution for the first time. Here we summarize the novel points of our work in comparison with Lebreuilly et al. (2020).

The essential points of our study are that we show, for the first time, with 3D non-ideal MHD simulation with dust growth, that the chain of (i) to (iv) (see above) self-consistently occurs. We also suggest that this mechanism can be a theoretical breakthrough to overcome the radial drift barrier.

On the other hand, Lebreuilly et al. (2020), in the context of (i), discuss the importance of dust growth in the discussion section (section 7.5, ”Caveat: coagulation and fragmentation during the collapse”) but do not include dust growth in their simulations. We improve this point by actually solving the dust growth in the simulation. In the context of (ii), they show that dust size with $\sim 100{\rm\mu}{\rm m}$ can be entrained by outflow. On the other hand, we show that much larger dust of $a_{d}>1{\rm mm}$ (which has more than ten times longer stopping time) can be entrained. In the context of (iii), from their results, whether the dust decoupling and dust ejection from outflow occur is not clear. They show that the dust is enriched (rather depleted) in outflow, suggesting that the dust and gas are relatively well coupled. This may be reasonable because the dust stopping time in outflow in their simulations is more than ten times smaller than in ours. Note, however, that they set a numerical upper limit of $1~{}{\rm km}~{}{\rm s}^{-1}$ on the dust-gas relative velocity, which is smaller than the relative velocity realized in our simulaiton (see arrows of figure 4-a or 4-b). We think this treatment may restrict the range of relative motion between gas and dust and affect the dynamics of dust grains such as ejection. We improved this point by explicitly solving the time dependence of the relative velocity without velocity limit. In the context of (iv), the reflux to the disk outer edge, seems to be missing in their results or presentation.

Finally, we briefly comment on the relation of our ash-fall phenomenon to the long-standing problem of chondrule and calcium-aluminum-rich inclusions (CAIs) formation. By showing that dust reflux can happen, our simulation supports the formation model of CAIs and chondrule with a bipolar outflow (e.g., Shu et al., 2001; Haugbølle et al., 2019), where the reflux is analytically discussed to occur from the inner to outer parts of a disk. New insights into the formation process of CAIs and chondrule will be gained with future studies of two-fluid radiation magnetohydrodynamics simulation of the dusty gas with the resolution to the vicinity of the central protostar.