"Ashfall" Induced by Molecular Outflow in Protostar Evolution

We solve two-fluid magnetohydrodynamics (MHD) equations for a dust–gas mixture. The details 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 and neutral dust as a single fluid is justified for investigating the dynamics of the dust with a size of a_d ≳ 10 μ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 , n_d is the dust number density, Δv_dust is the collision velocity between the dust grains, D_d /Dt = ∂/∂t + v_d · ∇, 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 mass when Δv_dust < Δv_frag and Δv_dust > Δv_frag, respectively, through dust collisions. To evaluate the dust relative velocity between the dust Δv_dust, we consider the subgrid scale turbulence and Brownian motion, and the relative velocity is given as Δv_dust = Δv_turb + Δ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 Δv_frag = 30 m s⁻¹. We can rewrite D_d /Dt as

where v is the barycentric velocity of the dust–gas mixture, ρ_g is the gas density, v _g is the gas velocity, ρ = ρ_g + ρ_d is the total density, ρ_d is the dust density, = ρ_d /ρ, and Δ v = v _d − v _g . Using this relation, we can rewrite Equation (1),

The smoothed particle hydrodynamics (SPH) discretization form of this equation is given as

where the last term is an artificial viscosity and is the signal velocity, where c_s,i is the sound velocity of the ith particle and e _ij = ( r _i − r _j )/∣ r _i − r _j ∣. We choose β_visc = 2.

The time step for dust advection is chosen to be

where h_i is the smoothing length of the ith particle.

During the simulations, the _i of a few particles are found to become negative (in particular in the outflow regions) with ∣_i ∣ ≲ 10⁻⁷. When this happens, we correct _i to be ∣_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 _i , as well as the number of the particles with _i < 0, is very small.

Our numerical simulations consider nonideal 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 μ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 is dynamically introduced when the density exceeds ρ_sink = 10⁻¹² g cm⁻³. The sink particle absorbs SPH particles with ρ > ρ_sink within r_sink < 1 au.

The initial condition of the cloud core is as follows. We adopt the density-enhanced Bonnor–Ebert sphere surrounded by a medium with a steep density profile of ρ ∝ r⁻⁴ as the initial condition, which is

with

where ξ_BE is a nondimensional 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. Here 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 ρ₀ = 7.3 × 10⁻¹⁸ g cm⁻³, ρ₀/ρ(R_c ) = 14, and f = 2.1. Then, the radius of the core is R_c = 4.8 × 10³ au, and the enclosed mass within R_c is M_c = 1 M_⊙. We adopt an angular velocity profile of with and Ω₀ = 2.3 × 10⁻¹³ s⁻¹. We adopt a constant magnetic field (B_x , B_y , B_z ) = (0, 0, 83 μG). The parameter α_therm ( ≡ E_therm/E_grav) is 0.4, where E_therm and E_grav are the thermal and gravitational energies of the central core (without the surrounding medium), respectively. The ratio of the rotational to gravitational energies β_rot within the core is β_rot ( ≡ E_rot/E_grav) = 0.03, where E_rot is the rotational energy of the core. The mass-to-flux ratio of the core is μ/μ_crit = 3. We resolve 1 M_⊙ with 3 × 10⁶ SPH particles. We adopt a dust density profile of , where f_dg = 10⁻² is the dust-to-gas mass ratio. The dust density profile has the same shape as the gas density profile in r ≲ 1.5R_c but is truncated at r ≥ 1.5R_c . The initial dust size is a_d = 0.1 μm.

We simulated the evolution of a protostar, a protoplanetary disk, and an outflow until 1.2 × 10⁴ yr after the epoch of protostar formation. Figure 1 shows the 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 the parent core, respectively.

Zoom In Zoom Out Reset image size

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

At t_* = 5 × 10³ yr (panels (2-a)–(2-d)), the maximum size of the dust reaches ∼100 μm at the center, and dust with a size of a_d ∼ 10 μm begins to be entrained by the outflow (panel (2-c)). Panels (2-a), (2-b), and (2-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 × 10³ yr (panels (3-a)–(3-d)), the dust with a size of a_d ∼ 100 μm is distributed over a large area of a few hundred au across. The streamlines 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 streamlines do not diverge but take a closed path, meaning that some of the large dust is refluxed from the outflow to the envelope. The ejection further decreases the amount of 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 × 10⁴ yr (panels (4-a)–(4-d)), the dust with a size of a_d ≳ 1 mm is distributed over a large area of several hundred au across (panel (4-c)). The streamlines and velocity field of the dust clearly indicate that the dust moves from the outflow back to the envelope (panels (4-a), (4-b), and (4-d)), while the streamlines and velocity field of the gas show that the gas is released to the interstellar medium (ISM). At this epoch, the gas and dust density distributions 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 has been refluxed to the envelope.

Figure 2 shows three-dimensional streamlines of the gas and dust at t_* = 1.1 × 10⁴ yr. The gas streamlines show that the gas is rotating and outflowing from the surface of the disk and blown away. By contrast, the dust streamlines 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 streamline of the gas, the centrifugal force drives the dust toward the radial direction. Gravity then pulls it back toward the disk, and the grown dust accretes to the edge of the disk.

Zoom In Zoom Out Reset image size

The overall process is schematically shown in Figure 3 and is very similar to ashfall from volcanic eruptions, where the 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 the "ashfall" phenomenon in protostar evolution.

Zoom In Zoom Out Reset image size

The reflux of the large dust due to the ashfall phenomenon increases the amount of 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 ∫ρ a_d dV/∫ρ dV in the disk, where the integration is performed within the disk. Here we define the disk as a region with the gas density ρ_g > 10⁻¹³ g cm⁻³. We find that, in an early phase of t_* < 5 × 10³ yr, where the mean dust size is a_d ≲ 10² μm, the dust-to-gas mass ratio stays constant at M_dust/M_gas = 10⁻², maintaining the initial value. This fact implies that the amount of refluxed dust is negligibly small, presumably because the outflowing dust is well coupled with gas and is not ejected from the outflow.

Zoom In Zoom Out Reset image size

Once the dust has grown to a size of a_d ≳ 10² μm at t_* ∼ 4 × 10³ yr, the dust-to-gas mass ratio begins to increase as a result of the increase of the amount of dust in the accretion flow from the envelope due to the dust reflux. The dust-to-gas mass ratio keeps increasing and becomes ∼ 0.01 at the end epoch of the simulation (t_* ∼ 1.2 × 10⁴ yr), implying that ∼10% of the dust in the disk is the refluxed large dust.

Note that the increase of the 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), it is not a consequence of dust drift in the prestellar collapse phase. We start from a dust size of 0.1 μ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⁻² during the prestellar collapse phase (see Figures (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 to the dust density may be minor (at least with an MRN-like size distribution). Therefore, we expect that an ∼10% increase in the 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 of 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 × 10⁴ 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 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 a roughly 50% increase from that of the typical ISM.

Zoom In Zoom Out Reset image size

A similar U-shaped dust-enhanced region has been observed in some class 0/I young stellar objects (YSOs), where the dust thermal emission was not detected, or was much weaker than the surrounding regions, inside the outflow (such as L1527, HOPS 136, B335, and IRAS 15398–3559; Yen et al. 2010; Fischer et al. 2014; Yen et al. 2017; Harris et al. 2018; Maury et al. 2018). The reason for the 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, this implies that the U-shaped dust enhancement is an indicator of the growth of dust in the disk and outflow.

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

• (i)  
  dust growth in the disk where the density is sufficiently high and the dust growth timescale t_growth is as short as t_growth ∼ 10³ yr (see Appendix B and Figure 4);
• (ii)  
  large dust with a_d > 1 mm being entrained by the powerful magnetized outflow launched from the inner region of the disk;
• (iii)  
  decoupling of the gas and dust in the outflow and dust ejection from the outflow due to centrifugal force, leading to a U-shaped dust-enhanced region around the outflow; and
• (iv)  
  mixing of the grown dust into the envelope accretion flow and the reflux of the large dust to the disk outer edge.

The ashfall supplies a nonnegligible amount of large dust to the outer edge of the disk even shortly after protostar formation (t ∼ 10⁴ yr) and keeps increasing (∼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 St is inversely proportional to the gas surface density Σ_g , given by St = π ρ_mat a_d /(2Σ_g ), where ρ_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 St < 1 at the inner region of the disk can have 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) connected to the envelope with an accretion shock (e.g., Miura et al. 2017). The envelope surface density at the connection with the edge of the disk (upstream of the accretion shock) is estimated as , where . Then, the disk surface density at the edge of the disk (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 is the mass accretion rate. Note that St can be even larger by a factor of with our parameters) when the accretion shock is adiabatic and the density increase from 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 ≳ 1 cm can already be St ≳ 1 when it is refluxed to the disk edge.

As such, the ashfall phenomenon can provide a pathway to circumvent the radial drift barrier and may allow dust to grow into planetesimals in the outer region of the 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 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 ∼ 10⁴ yr after protostar formation (i.e., investigating the evolution of class 0/I YSOs), we can expect that the ashfall mechanism also plays a role even in class II YSOs (or possibly more importantly because the Stokes number of the refluxed dust at the disk edge increases as envelope accretion diminishes; see Equation (10)).

In addition, a significant amount of large dust whose Stokes number is close to unity refluxed to the disk provides an ideal condition for the growth of the secular gravitational instability (SGI; Takahashi & Inutsuka 2014, 2016; Tominaga et al. 2018, 2020), which is expected to be a powerful mechanism for creating rings 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 ashfall and subsequent SGI can be an explanation for the multiple ringlike structures observed in protoplanetary disks (ALMA Partnership et al. 2015; Fedele et al. 2017, 2018; Andrews 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 jetlike structures. The coexistence of outflows/jets and very ordered axisymmetric multiple ringlike structures in the disk has been puzzling us since its discovery in 2014. Now we are proposing that these two processes might actually be cooperating in such an object.

Another important implication of ashfall is that the reflux of the large dust into the envelope of ∼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 β 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 a nonideal 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 a three-dimensional nonideal MHD simulation with dust growth, that the chain of steps (i)–(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 step (i), discussed the importance of dust growth in the discussion section (Section 7.5) but did not include dust growth in their simulations. We improve this point by actually solving the dust growth in the simulation. In the context of step (ii), they showed that the dust size with ∼100 μm can be entrained by outflow. On the other hand, we show that much larger dust of a_d > 1 mm (which has a more than 10 times longer stopping time) can be entrained. In the context of step (iii), from their results, whether the dust decoupling and ejection from the outflow occurs is not clear. They showed that the dust is enriched (rather than depleted) in the outflow, suggesting that the dust and gas are relatively well coupled. This may be reasonable because the dust stopping time in the outflow in their simulations is more than 10 times smaller than in ours. Note, however, that they set a numerical upper limit of 1 km s⁻¹ on the dust–gas relative velocity, which is smaller than the relative velocity realized in our simulation (see arrows of Figures (3-a) and (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 improve this point by explicitly solving the time dependence of the relative velocity without a velocity limit. In the context of step (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 ashfall phenomenon to the long-standing problem of chondrule and calcium-aluminum-rich inclusion (CAI) formation. By showing that dust reflux can happen, our simulation supports the formation model of CAIs and chondrules 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 chondrules will be gained with future studies of two-fluid radiation MHD simulations of the dusty gas with the resolution to the vicinity of the central protostar.