Ejection of Chondrules from Fluffy Matrices

The goal of this paper is to explore whether fluffy aggregates composed of chondrules and matrix grains can grow into chondrite parent bodies without chondrule losses. Therefore, it is necessary to calculate the structure of the dust aggregates, such as the radius and the bulk density, and the collision velocity. Here we assume that there are two types of aggregates, that is, chondrule-dust compound aggregates (hereinafter referred to as compound aggregates or CAs) and aggregates that are composed purely of matrix grains (hereinafter referred to as matrix aggregates or MAs). The formation processes of MAs and CAs are illustrated in Figure 1.

Zoom In Zoom Out Reset image size

In Stage I, the growth of MAs primarily occurs owing to Brownian motion, and they reduce their internal densities via fractal aggregation. We define Stage I as where the timescale of collisions between two MAs ${t}_{\mathrm{MA}\to \mathrm{MA}}$ is shorter than the timescale for MAs to accrete onto chondrules/CAs ${t}_{\mathrm{MA}\to \mathrm{CA}}$. Then, in Stage II, the masses of the CAs increase via the accretion of the grown MAs, and finally, in Stage III, the growth of CAs via the collision of two CAs with static compression takes place. We also define Stage III as where the growth timescale of CAs by the collision of two CAs ${t}_{\mathrm{grow}}^{\mathrm{CA}\leftarrow \mathrm{CA}}$ is shorter than the growth timescale of CAs by the accretion of the grown MAs ${t}_{\mathrm{grow}}^{\mathrm{CA}\leftarrow \mathrm{MA}}$. The growth timescale of CAs by the collision of two CAs ${t}_{\mathrm{grow}}^{\mathrm{CA}\leftarrow \mathrm{CA}}$ is the same as the timescale of collisions between two CAs ${t}_{\mathrm{CA}\to \mathrm{CA}}$. On the other hand, the growth timescale of CAs by the accretion of the grown MAs is ${t}_{\mathrm{grow}}^{\mathrm{CA}\leftarrow \mathrm{MA}}=({m}_{\mathrm{CA}}/{m}_{\mathrm{MA}}){t}_{\mathrm{CA}\to \mathrm{MA}}$, where ${m}_{\mathrm{MA}}$ and ${m}_{\mathrm{CA}}$ are the masses of MAs and CAs, respectively, and ${t}_{\mathrm{CA}\to \mathrm{MA}}$ is the timescale for CAs to collide with MAs. We can imagine that fine-grained matrices could help CAs grow in protoplanetary disks due to their large adhesive force. Forming highly porous matrices also encourages CAs to grow without significant radial drift because fluffy dust aggregates have large collisional cross-sections and can grow rapidly.

In this study, we focus on Stage III and only consider the existence and growth of CAs for simplicity. We calculate the density evolution of the CAs by considering the compression processes and examine whether the CAs can retain chondrules inside their fluffy matrices.

We assume that the disk structure is in steady state, i.e., the mass accretion rate $\dot{M}$ is constant with distance from the Sun R. Because a large amount of fine dust particles must have been present when planetesimals formed in the solar nebula, we also assume that the disk is optically thick and that the main heating source at the midplane is viscous dissipation. Then the local surface density Σ and the midplane temperature ${T}_{{\rm{m}}}$ are given by (e.g., Pringle 1981; Oka et al. 2011)

$$\begin{eqnarray}&&{\rm{\Sigma }}=\displaystyle \frac{\dot{M}}{3\pi {\nu }_{\mathrm{acc}}},\end{eqnarray} \tag{ 1 }$$

and

$$\begin{eqnarray}&&{T}_{{\rm{m}}}={\left(\displaystyle \frac{9\kappa {\rm{\Sigma }}}{32}\displaystyle \frac{{{GM}}_{\odot }\dot{M}}{4\pi {\sigma }_{\mathrm{SB}}{R}^{3}}\right)}^{1/4},\end{eqnarray} \tag{ 2 }$$

where ${\nu }_{\mathrm{acc}}$ is the effective viscosity of the accretion disk, κ is the Rosseland mean opacity of the disk medium, $G=6.67\,\times {10}^{-8}\ \mathrm{dyn}\ {\mathrm{cm}}^{2}\ {{\rm{g}}}^{-2}$ is the gravitational constant, and ${\sigma }_{\mathrm{SB}}\,=5.67\times {10}^{-5}\ \mathrm{erg}\ {\mathrm{cm}}^{-2}\ {{\rm{K}}}^{-4}\ {{\rm{s}}}^{-1}$ is the Stefan-Boltzmann constant. The solar mass M_⊙ is ${M}_{\odot }=1.99\times {10}^{33}\ {\rm{g}}$. The Kepler frequency ${{\rm{\Omega }}}_{{\rm{K}}}$ is given by ${{\rm{\Omega }}}_{{\rm{K}}}=\sqrt{{{GM}}_{\odot }/{R}^{3}}$ and the isothermal sound velocity ${c}_{{\rm{s}}}$ is given by ${c}_{{\rm{s}}}=\sqrt{{k}_{{\rm{B}}}{T}_{{\rm{m}}}/{m}_{{\rm{g}}}}$, where ${k}_{{\rm{B}}}=1.38\times {10}^{-16}\ \mathrm{erg}\ {{\rm{K}}}^{-1}$ is the Boltzmann constant and ${m}_{{\rm{g}}}$ is the mean molecular weight. We set the mean molecular weight to ${m}_{{\rm{g}}}=2.34{m}_{{\rm{H}}}$, where ${m}_{{\rm{H}}}=1.67\times {10}^{-24}\ {\rm{g}}$ is the mass of a hydrogen atom. Then, the effective viscosity ${\nu }_{\mathrm{acc}}$ is written as ${\nu }_{\mathrm{acc}}={\alpha }_{\mathrm{acc}}{{c}_{{\rm{s}}}}^{2}/{{\rm{\Omega }}}_{{\rm{K}}}$ (Shakura & Sunyaev 1973). Here, we assume ${\alpha }_{\mathrm{acc}}={10}^{-2}$ in this study (e.g., Hartmann et al. 1998).

The surface densities of dust and gas, ${{\rm{\Sigma }}}_{{\rm{d}}}$ and ${{\rm{\Sigma }}}_{{\rm{g}}}$, respectively, are given by ${{\rm{\Sigma }}}_{{\rm{d}}}=Z{\rm{\Sigma }}$ and ${{\rm{\Sigma }}}_{{\rm{g}}}=(1-Z){\rm{\Sigma }}$, where Z is the vertically integrated dust-to-total mass ratio. The vertical structure of the gas is assumed to be in hydrostatic equilibrium, and the midplane gas density ${\rho }_{{\rm{g}}}$ is given by ${\rho }_{{\rm{g}}}={{\rm{\Sigma }}}_{{\rm{g}}}/(\sqrt{2\pi }{h}_{{\rm{g}}})$, where ${h}_{{\rm{g}}}={c}_{{\rm{s}}}/{{\rm{\Omega }}}_{{\rm{K}}}$ is the gas scale height. The mean-free path of gas molecules ${\lambda }_{{\rm{m}}}$ is given by ${\lambda }_{{\rm{m}}}\,={m}_{{\rm{g}}}/({\sigma }_{\mathrm{mol}}{\rho }_{{\rm{g}}})$, where ${\sigma }_{\mathrm{mol}}=2\times {10}^{-15}\ {\mathrm{cm}}^{2}$ is the collisional cross-section of the gas molecules (Okuzumi et al. 2012). The opacity of the disk κ is determined by the surface density of the fine dust, i.e., the matrix grains and small MAs. We assume that the mass fraction of matrix grains to all solids χ is equal in the disk and in the CAs. According to Scott (2007), the volume fraction of the matrices in ordinary chondrites is approximately 10%, and the mass fraction is also approximately 10%. Therefore, we set $\chi =0.1$ as the canonical value and the opacity of the disk is given by $\kappa =\chi Z{\kappa }_{{\rm{d}}}$, where ${\kappa }_{{\rm{d}}}\,={10}^{3}\ {\mathrm{cm}}^{2}\ {{\rm{g}}}^{-1}$ is the mass opacity of the fine silicate grains (e.g., Draine 1985). The mass opacity of highly porous aggregates composed of fine silicate grains is also on the order of ${10}^{3}\ {\mathrm{cm}}^{2}\ {{\rm{g}}}^{-1}$ (Kataoka et al. 2014).

We consider the Brownian motion, turbulence, and radial and azimuthal drift as sources of the relative velocity of the gas and dust particles. We write the relative velocity v as the root sum square of these contributions:

$$\begin{eqnarray}&&v=\sqrt{{v}_{{\rm{B}}}^{2}+{v}_{{\rm{t}}}^{2}+{v}_{{\rm{r}}}^{2}+{{v}_{\phi }}^{2}},\end{eqnarray} \tag{ 3 }$$

where ${v}_{{\rm{B}}},{v}_{{\rm{t}}},{v}_{{\rm{r}}}$, and ${v}_{\phi }$ are the relative velocities induced by Brownian motion, turbulence, radial drift, and azimuthal drift, respectively. In this study, we assume that the dust aggregates have no mass distribution and that aggregates that have the same mass have the same volume for simplicity. The collision velocity ${\rm{\Delta }}v$ is

$$\begin{eqnarray}&&{\rm{\Delta }}v=\sqrt{{({\rm{\Delta }}{v}_{{\rm{B}}})}^{2}+{({\rm{\Delta }}{v}_{{\rm{t}}})}^{2}},\end{eqnarray} \tag{ 4 }$$

where ${\rm{\Delta }}{v}_{{\rm{B}}}$ and ${\rm{\Delta }}{v}_{{\rm{t}}}$ are the collisional velocities induced by the Brownian motion and turbulence, respectively. The model of the dust dynamics is the same as that of Arakawa & Nakamoto (2016b).

2.3.1. Brownian Motions

The Brownian-motion-induced velocities ${v}_{{\rm{B}}}$ and ${\rm{\Delta }}{v}_{{\rm{B}}}$ are given by

$$\begin{eqnarray}&&{v}_{{\rm{B}}}=\sqrt{\displaystyle \frac{8{k}_{{\rm{B}}}T}{\pi m}},\end{eqnarray} \tag{ 5 }$$

and

$$\begin{eqnarray}&&{\rm{\Delta }}{v}_{{\rm{B}}}=\sqrt{\displaystyle \frac{16{k}_{{\rm{B}}}T}{\pi m}},\end{eqnarray} \tag{ 6 }$$

where m is the mass of an aggregate.

2.3.2. Turbulent Motions

The dynamics of dust particles in a disk is affected by disk gas turbulence. The turnover time and the root-mean-squared random velocity of the largest turbulent eddies, ${t}_{{\rm{L}}}$ and ${v}_{{\rm{L}}}$, respectively, are given by ${t}_{{\rm{L}}}={{{\rm{\Omega }}}_{{\rm{K}}}}^{-1}$ and ${v}_{{\rm{L}}}=\sqrt{{\alpha }_{{\rm{t}}}}{c}_{{\rm{s}}}$, where ${\alpha }_{{\rm{t}}}$ is the dimensionless turbulent parameter at the midplane (Cuzzi & Hogan 2003). The value of the dimensionless parameter ${\alpha }_{{\rm{t}}}$ for our solar nebula is unclear; however, measurements of CO emission lines reveal that some circumstellar disks, such as the disk around HD 163296, have a small turbulent parameter (${\alpha }_{{\rm{t}}}\lesssim {10}^{-3}$ Flaherty et al. 2015). The disk around HL Tau also has a small turbulent viscosity, which is equivalent to an ${\alpha }_{{\rm{t}}}$ at a few 10⁻⁴ when considering the dust settling suggested by observations of the gap structure (Pinte et al. 2016); meanwhile, ${\alpha }_{\mathrm{acc}}\sim {10}^{-2}$ is indicated by the disk mass distribution and the accretion rate ($8.7\,\times {10}^{-8}\ {M}_{\odot }\ {\mathrm{yr}}^{-1};$ Beck et al. 2010). If the magnetorotational instability is active, ${\alpha }_{{\rm{t}}}$ may be up to 10⁻³ (e.g., Balbus & Hawley 1991), whereas ${\alpha }_{{\rm{t}}}$ might be lower than 10⁻³ to 10⁻⁵ if we focus on the low-ionization-fraction regions (e.g., Gammie 1996; Mori & Okuzumi 2016). In this study, we assume that the turbulence is weak in the midplane of the inner disk and that ${\alpha }_{{\rm{t}}}={10}^{-4}$. The turbulent viscosity ${\nu }_{{\rm{t}}}$ is given by ${\nu }_{{\rm{t}}}={\alpha }_{{\rm{t}}}{c}_{{\rm{s}}}{h}_{{\rm{g}}}$. The extent of the gas turbulence is determined by the turbulent Reynolds number ${\mathrm{Re}}_{{\rm{t}}}$. The turbulent Reynolds number ${\mathrm{Re}}_{{\rm{t}}}$ is given by ${\mathrm{Re}}_{{\rm{t}}}\,=\,{\nu }_{{\rm{t}}}/{\nu }_{{\rm{m}}}$, where ${\nu }_{{\rm{m}}}={\lambda }_{{\rm{m}}}{v}_{\mathrm{th}}/2$ is the molecular viscosity, and ${v}_{\mathrm{th}}=\sqrt{8/\pi }{c}_{{\rm{s}}}$ is the thermal velocity of the gas molecules. The turnover time and the root-mean-squared random velocity of the smallest eddies, ${t}_{\eta }$ and ${v}_{\eta }$, respectively, are given by ${t}_{\eta }={\mathrm{Re}}_{{\rm{t}}}^{-1/2}{t}_{{\rm{L}}}$ and ${v}_{\eta }={{\mathrm{Re}}_{{\rm{t}}}}^{-1/4}{v}_{{\rm{L}}}$.

The key parameter that determines the motion of dust aggregates in the gas disk is the normalized stopping time, the Stokes number $\mathrm{St}$, given by

$$\begin{eqnarray}&&\mathrm{St}={{\rm{\Omega }}}_{{\rm{K}}}{t}_{{\rm{s}}},\end{eqnarray} \tag{ 7 }$$

where ${t}_{{\rm{s}}}$ is the stopping time of the dust aggregates. The stopping time depends on the radius of the aggregates a and the particle Reynolds number ${\mathrm{Re}}_{{\rm{p}}}\,=\,2{av}/{\nu }_{{\rm{m}}}$. The stopping time ${t}_{{\rm{s}}}$ is given by Probstein & Fassio (1970) and Weidenschilling (1977) as follows. If the radius of an aggregate a is smaller than $(9/4){\lambda }_{{\rm{m}}}$, then the stopping time of the dust aggregates is determined by Epstein's law, ${t}_{{\rm{s}}}^{\mathrm{Ep}}=3m/(4\pi {\rho }_{{\rm{g}}}{v}_{\mathrm{th}}{a}^{2})$. Conversely, when the radius a is larger than $(9/4){\lambda }_{{\rm{m}}}$, the stopping time depends on the particle Reynolds number. If the particle Reynolds number ${\mathrm{Re}}_{{\rm{p}}}$ is less than unity, the motion of the aggregate is determined by Stokes' law, ${t}_{{\rm{s}}}^{\mathrm{St}}=m/(6\pi {\rho }_{{\rm{g}}}{\nu }_{{\rm{m}}}a)$, while for the case of $1\,\leqslant {\mathrm{Re}}_{{\rm{p}}}\,\lt \,{54}^{5/3}\sim 800$, the stopping time is obtained from Allen's law, ${t}_{{\rm{s}}}^{\mathrm{Al}}=({2}^{3/5}m)/(12\pi {\rho }_{{\rm{g}}}{{\nu }_{{\rm{m}}}}^{3/5}{v}^{2/5}{a}^{7/5})$. Finally, for the case in which the aggregate moves fast and the particle Reynolds number is larger than ${54}^{5/3}$, we use Newton's law to obtain the stopping time, ${t}_{{\rm{s}}}^{\mathrm{Nw}}=(9m)/(2\pi {\rho }_{{\rm{g}}}{{va}}^{2})$.

We use the analytic formulae derived by Ormel & Cuzzi (2007) for the turbulence-driven velocity. If the stopping time ${t}_{{\rm{s}}}$ is shorter than the turnover time of the smallest eddies ${t}_{\eta }$, the relative velocity between the dust particle to disk gas ${v}_{{\rm{t}}}$ is given by ${v}_{{\rm{t}}}^{{\rm{S}}}$ such that

$$\begin{eqnarray}&&{v}_{{\rm{t}}}^{{\rm{S}}}={{\mathrm{Re}}_{{\rm{t}}}}^{1/4}\mathrm{St}{v}_{{\rm{L}}}.\end{eqnarray} \tag{ 8 }$$

For the case of ${t}_{\eta }\ll {t}_{{\rm{s}}}\ll {t}_{{\rm{L}}}$, the relative velocity ${v}_{{\rm{t}}}$ is given by ${v}_{{\rm{t}}}^{{\rm{M}}}$ such that

$$\begin{eqnarray}&&{v}_{{\rm{t}}}^{{\rm{M}}}=1.7{\mathrm{St}}^{1/2}{v}_{{\rm{L}}},\end{eqnarray} \tag{ 9 }$$

and for the case of ${t}_{{\rm{s}}}\gg {t}_{{\rm{L}}}$, the relative velocity ${v}_{{\rm{t}}}$ is given by ${v}_{{\rm{t}}}^{{\rm{L}}}$ such that

$$\begin{eqnarray}&&{v}_{{\rm{t}}}^{{\rm{L}}}={\left(1+\displaystyle \frac{1}{1+\mathrm{St}}\right)}^{1/2}{v}_{{\rm{L}}}.\end{eqnarray} \tag{ 10 }$$

Here, we assume that the turbulence-driven velocity ${v}_{{\rm{t}}}$ is given by the minimum of these three terms for continuity and better handling of ${v}_{{\rm{t}}}$:

$$\begin{eqnarray}&&{v}_{{\rm{t}}}=\min ({v}_{{\rm{t}}}^{{\rm{S}}},{v}_{{\rm{t}}}^{{\rm{M}}},{v}_{{\rm{t}}}^{{\rm{L}}}).\end{eqnarray} \tag{ 11 }$$

The turbulence-driven velocity between two dust particles was also obtained by Ormel & Cuzzi (2007). We assume the collision velocity of two particles with a stopping time $t\lt {t}_{\eta }$ to be

$$\begin{eqnarray}&&{\rm{\Delta }}{v}_{{\rm{t}}}^{{\rm{S}}}=\varepsilon {{\mathrm{Re}}_{{\rm{t}}}}^{1/4}\mathrm{St}{v}_{{\rm{L}}},\end{eqnarray} \tag{ 12 }$$

where $\varepsilon =0.1$ is a nondimensional value representing the variation in the density fluctuation between same-mass aggregates, which is neglected in this study (e.g., Okuzumi et al. 2011). Note that there is a large uncertainty in the estimation of the ε for CAs; the estimation by Okuzumi et al. (2011) is only valid for non-compressed fractal aggregates composed of 10–10⁶ monomers. For the intermediate region, we assume ${\rm{\Delta }}{v}_{{\rm{t}}}^{{\rm{M}}}$ such that

$$\begin{eqnarray}&&{\rm{\Delta }}{v}_{{\rm{t}}}^{{\rm{M}}}=1.4{\mathrm{St}}^{1/2}{v}_{{\rm{L}}},\end{eqnarray} \tag{ 13 }$$

and the relative velocity for high Stokes number aggregates ${\rm{\Delta }}{v}_{{\rm{t}}}^{{\rm{L}}}$ is

$$\begin{eqnarray}&&{\rm{\Delta }}{v}_{{\rm{t}}}^{{\rm{L}}}={\left(\displaystyle \frac{2}{1+\mathrm{St}}\right)}^{1/2}{v}_{{\rm{L}}}.\end{eqnarray} \tag{ 14 }$$

Here, we also assume that the turbulence-driven relative velocity ${\rm{\Delta }}{v}_{{\rm{t}}}$ is given by the minimum of these three terms:

$$\begin{eqnarray}&&{\rm{\Delta }}{v}_{{\rm{t}}}=\min ({\rm{\Delta }}{v}_{{\rm{t}}}^{{\rm{S}}},{\rm{\Delta }}{v}_{{\rm{t}}}^{{\rm{M}}},{\rm{\Delta }}{v}_{{\rm{t}}}^{{\rm{L}}}).\end{eqnarray} \tag{ 15 }$$

2.3.3. Systematic Motions

Dust aggregates in the gas disk have systematic velocities such as the radial drift velocity, the azimuthal velocity, and the settling velocity. In this study, we consider the radial drift velocity ${v}_{{\rm{r}}}$ and the azimuthal velocity ${v}_{\phi }$ given by (Weidenschilling 1977)

$$\begin{eqnarray}&&{v}_{{\rm{r}}}=-\displaystyle \frac{2\mathrm{St}}{1+{\mathrm{St}}^{2}}\eta {v}_{{\rm{K}}},\end{eqnarray} \tag{ 16 }$$

and

$$\begin{eqnarray}&&{v}_{\phi }=\left(1-\displaystyle \frac{1}{1+{\mathrm{St}}^{2}}\right)\eta {v}_{{\rm{K}}},\end{eqnarray} \tag{ 17 }$$

where $\eta {v}_{{\rm{K}}}$ is the difference between the Keplerian velocity and the gas rotational velocity. The Keplerian velocity ${v}_{{\rm{K}}}$ is given by ${v}_{{\rm{K}}}=R{{\rm{\Omega }}}_{{\rm{K}}}$, and the dimensionless coefficient η is given by (Nakagawa et al. 1986)

$$\begin{eqnarray}&&\eta =-\displaystyle \frac{1}{2}{\left(\displaystyle \frac{{c}_{{\rm{s}}}}{{v}_{{\rm{K}}}}\right)}^{2}\displaystyle \frac{\partial \mathrm{ln}({\rho }_{{\rm{g}}}{{c}_{{\rm{s}}}}^{2})}{\partial \mathrm{ln}R}.\end{eqnarray} \tag{ 18 }$$

The aerodynamic properties of the aggregates depend on their bulk densities. In this paper, we assume that all CAs have the same mass and density and that the matrix regions of the CAs have a uniform density. Then, the bulk filling factor of the CAs ${\phi }_{\mathrm{CA}}$ is given by

$$\begin{eqnarray}&&{\phi }_{\mathrm{CA}}={\left(\displaystyle \frac{\chi }{{\phi }_{\mathrm{mat}}}+\displaystyle \frac{1-\chi }{1}\right)}^{-1},\end{eqnarray} \tag{ 19 }$$

where ${\phi }_{\mathrm{mat}}$ is the filling factor of the matrix regions and χ is the mass fraction of the matrices in a CA. The bulk density of the CAs ${\rho }_{\mathrm{CA}}$ is

$$\begin{eqnarray}&&{\rho }_{\mathrm{CA}}={\rho }_{0}{\phi }_{\mathrm{CA}},\end{eqnarray} \tag{ 20 }$$

where ${\rho }_{0}=3\ {\rm{g}}\ {\mathrm{cm}}^{-3}$ is the material density of the chondrules and matrix grains.

There are several processes that affect the density of the dust aggregates: the hit-and-stick growth of two colliding aggregates without compaction (e.g., Meakin 1991; Okuzumi et al. 2009), the dynamical compression via high-speed collisions (e.g., Paszun & Dominik 2009; Suyama et al. 2012), and static compression via the ram pressure of the disk gas and/or via the self-gravity of large planetary bodies (e.g., Seizinger et al. 2012; Kataoka et al. 2013a).

Dynamical compression, however, does not work in our calculations, and the CAs are compressed by the static pressure. Note that Ormel et al. (2008) assumed that dust rims around chondrules are formed via direct accretion of monomer grains. Therefore, Ormel et al. (2008) assumed that the initial filling factor of the matrix regions of CAs is ${\phi }_{\mathrm{mat}}=0.15$, which is the filling factor obtained by the particle-cluster aggregation. By contrast, we assume chondrules obtain fluffy dust rims via the accretion of highly porous MAs, and ${\phi }_{\mathrm{mat}}$ is far smaller than 0.15 at the start of Stage III in our calculations. In addition, Ormel et al. (2008) employs the collisional compression model obtained by Ormel et al. (2007); however, Wada et al. (2008) revealed that Ormel et al. (2007) overestimates the density of fluffy aggregates at the maximum compression.

Kataoka et al. (2013a) revealed that the static compression of a highly porous aggregate depends on the pressure that works on the aggregate P and the rolling energy ${E}_{\mathrm{roll}}$. The rolling energy is given by (Dominik & Tielens 1995)

$$\begin{eqnarray}&&{E}_{\mathrm{roll}}=6{\pi }^{2}\gamma {a}_{0}\xi ,\end{eqnarray} \tag{ 21 }$$

where γ is the surface energy of the silicate, a₀ is the radius of the monomers (i.e., matrix grains), and ξ is the critical displacement for rolling. We assume that the radius of the matrix grains is ${a}_{0}=0.025\ \mu {\rm{m}}$ in this study, which is the typical radius of opaque grains in pristine chondrites such as Acfer 094 and QUE 99177 (Vaccaro et al. 2015). We employ the canonical value of the surface energy of the silicate, which is $\gamma =25\ \mathrm{erg}\ {\mathrm{cm}}^{-2}$ (Kendall et al. 1987). The lower limit of ξ is $0.2\ \mathrm{nm}$, which is equal to the interatomic distance (Dominik & Tielens 1997), and the upper limit is the radius of the contact surface area of two matrix grains ${a}_{\mathrm{cont}}$ (Heim et al. 1999). Krijt et al. (2014) showed that ξ is approximately 0.3–$1\ \mathrm{nm}$ for micron-sized amorphous silica spheres, and we assume that the critical displacement for rolling is $\xi =0.3\ \mathrm{nm}$ for submicron-sized monomers. Therefore, we can calculate the filling factor of a highly porous dust aggregate using the relation obtained by Kataoka et al. (2013a), and the filling factor of matrices in CAs is given by

$$\begin{eqnarray}&&{\phi }_{\mathrm{mat}}={\left(\displaystyle \frac{{{a}_{0}}^{3}P}{{E}_{\mathrm{roll}}}\right)}^{1/3}.\end{eqnarray} \tag{ 22 }$$

The sources of the pressure P are the ram pressure ${P}_{\mathrm{gas}}$ and the self-gravitational pressure ${P}_{\mathrm{grav}}$ (Kataoka et al. 2013b). For CAs whose radii are ${a}_{\mathrm{CA}}$ and masses are ${m}_{\mathrm{CA}}$, the ram pressure of the disk gas is estimated to be ${P}_{\mathrm{gas}}=({m}_{\mathrm{CA}}v)/(\pi {{a}_{\mathrm{CA}}}^{2}{t}_{{\rm{s}}})$, whereas the self-gravitational pressure is ${P}_{\mathrm{grav}}=(G{{m}_{\mathrm{CA}}}^{2})/(\pi {{a}_{\mathrm{CA}}}^{4})$. We calculate the filling factors and densities assuming that $P\,=\max ({P}_{\mathrm{gas}},{P}_{\mathrm{grav}})$.

Chondrules contained in fluffy matrices move in the aggregates when fluffy aggregates collide at large velocities. We estimate that these CAs cannot retain the chondrules inside their matrices if the stopping length of the chondrules exceeds the radius of the CAs (Figure 2). Here, we assume that the stopping length is the penetration length for stopping via deceleration forces, ${L}_{\mathrm{pen}}$.

Zoom In Zoom Out Reset image size

The penetration of chondrules in fluffy matrices is the same process as that for the penetration of spherical projectiles in a silica aerogel (Niimi et al. 2011) and the intrusion of chondrule-analogs into matrix analogs (Machii et al. 2013). Therefore, the equation of motion is

$$\begin{eqnarray}&&{m}_{\mathrm{ch}}\displaystyle \frac{{{dv}}_{\mathrm{pen}}}{{dt}}={F}_{\mathrm{pen}},\end{eqnarray} \tag{ 23 }$$

where ${m}_{\mathrm{ch}}$ is the mass of the chondrules. In this study, we assume that the radius of the chondrules is ${a}_{\mathrm{ch}}=0.01\ \mathrm{cm}$, which is a typical size of chondrules in H ordinary chondrites and EH enstatite chondrites (Rubin 2000). The mass of the chondrules is ${m}_{\mathrm{ch}}=(4\pi /3){{a}_{\mathrm{ch}}}^{3}{\rho }_{0}$. Note that some other groups (e.g., L and LL ordinary chondrites) contain chondrules whose radii are a few times larger than that, which we assume here.

There are two mechanisms for decelerating chondrules in matrices: "hydrodynamic" drag for high-speed penetration and "crushing" drag for low-speed penetration (e.g., Niimi et al. 2011). In the case of high-speed penetration, the deceleration process is due to the "hydrodynamic" drag caused by the ram pressure, and the drag force ${F}_{\mathrm{pen}}$ is given by

$$\begin{eqnarray}&&{F}_{\mathrm{pen}}^{\mathrm{HD}}=-\displaystyle \frac{{C}_{{\rm{d}}}}{2}\pi {{a}_{\mathrm{ch}}}^{2}{\rho }_{\mathrm{mat}}{{v}_{\mathrm{pen}}}^{2},\end{eqnarray} \tag{ 24 }$$

where ${C}_{{\rm{d}}}=1$ is the drag coefficient obtained from impact experiments (e.g., Trucano & Grady 1995). Conversely, in the case of low-speed penetration, the deceleration process is due to "crushing" against the compressive strength of the matrices and the equation of motion is given by

$$\begin{eqnarray}&&{F}_{\mathrm{pen}}^{\mathrm{CR}}=-\pi {{a}_{\mathrm{ch}}}^{2}{P}_{{\rm{c}}},\end{eqnarray} \tag{ 25 }$$

where ${P}_{{\rm{c}}}$ is the compressive strength of the matrices (e.g., Domínguez et al. 2004). The compressive strength ${P}_{{\rm{c}}}$ has been studied in experiments (e.g., Blum & Schräpler 2004; Güttler et al. 2009) and in numerical simulations (e.g., Seizinger et al. 2012; Kataoka et al. 2013a). In this study, we calculate the ${P}_{{\rm{c}}}$ of highly porous matrices using the following equation (Kataoka et al. 2013a):

$$\begin{eqnarray}&&{P}_{{\rm{c}}}=\displaystyle \frac{{E}_{\mathrm{roll}}}{{{a}_{0}}^{3}}{{\phi }_{\mathrm{mat}}}^{3}.\end{eqnarray} \tag{ 26 }$$

Here, we can rewrite the equation of motion as

$$\begin{eqnarray}&&\displaystyle \frac{{{dv}}_{\mathrm{pen}}}{{dt}}=-\max (A{{v}_{\mathrm{pen}}}^{2},B),\end{eqnarray} \tag{ 27 }$$

where $A=3{C}_{{\rm{d}}}{\phi }_{\mathrm{mat}}/(8{a}_{\mathrm{ch}})$ and $B=3{P}_{{\rm{c}}}/(4{\rho }_{0}{a}_{\mathrm{ch}})$. Then we can obtain the penetration length ${L}_{\mathrm{pen}}$ as a function of the initial penetration velocity v₀. In the case of high-speed penetration (${v}_{0}\geqslant \sqrt{B/A}$), ${L}_{\mathrm{pen}}$ is given by

$$\begin{eqnarray}&&{L}_{\mathrm{pen}}=\displaystyle \frac{1}{A}\mathrm{ln}({v}_{0}\sqrt{A/B})+\displaystyle \frac{1}{2A},\end{eqnarray} \tag{ 28 }$$

and in the case of low-speed penetration (${v}_{0}\lt \sqrt{B/A}$), ${L}_{\mathrm{pen}}$ is given by

$$\begin{eqnarray}&&{L}_{\mathrm{pen}}=\displaystyle \frac{{{v}_{0}}^{2}}{2B}.\end{eqnarray} \tag{ 29 }$$

If two colliding CAs have the same mass and the matrices of the two colliding aggregates stop immediately after the collision, then the initial penetration velocity v₀ can be evaluated as ${\rm{\Delta }}v/2$. However, the motion of the matrices after the collision is not yet understood. We examine two cases in this study: ${v}_{0}={\rm{\Delta }}v/2$ and ${v}_{0}={\rm{\Delta }}v/8$. Note that, if the mass ratio of the projectile to the target is small, then v₀ of the target decreases while v₀ of the impactor increases when considering the motions relative to the center of mass.

We also estimate the critical velocity for the collisional growth of aggregates ${\rm{\Delta }}{v}_{\mathrm{crit}}$ from the scaling relation obtained by Wada et al. (2009). The impact energy of two colliding aggregates ${E}_{\mathrm{imp}}$ is ${E}_{\mathrm{imp}}=m{({\rm{\Delta }}v)}^{2}/4$, and the critical impact energy for growth ${E}_{\mathrm{imp},\mathrm{crit}}$ is ${E}_{\mathrm{imp},\mathrm{crit}}=30{N}_{\mathrm{total}}{E}_{\mathrm{break}}$, where ${N}_{\mathrm{total}}$ is the total number of monomers contained in the colliding aggregates, and ${E}_{\mathrm{break}}$ is the energy for breaking a single contact between two monomers (Wada et al. 2009). The energy for breaking a monomer–monomer contact ${E}_{\mathrm{break}}$ is given by Wada et al. (2007). Clearly, the total binding energy of a CA is dominated by the contacts between matrix grains. Therefore, the total number of monomers in the two colliding CAs is ${N}_{\mathrm{total}}=2\chi {m}_{\mathrm{CA}}/{m}_{0}$, where ${m}_{0}=(4\pi /3){a}_{0}^{3}{\rho }_{0}$ is the mass of monomers, and the critical velocity for the collisional growth of CAs ${\rm{\Delta }}{v}_{\mathrm{crit}}$ is given by

$$\begin{eqnarray}&&{\rm{\Delta }}{v}_{\mathrm{crit}}=\sqrt{\displaystyle \frac{240\chi {E}_{\mathrm{break}}}{{m}_{0}}}=5.5\times {10}^{2}\ \mathrm{cm}\ {{\rm{s}}}^{-1}.\end{eqnarray} \tag{ 30 }$$

Note, however, that the critical collision velocity of dust aggregates is still open to argument (e.g., Tanaka et al. 2012; Krijt et al. 2013). Hence we do not consider fragmentation of CAs in this study.

First, we define the "Epstein-drag region" as the region in which the CA affects the Epstein drag and the chondrules are decelerated by the compressive stress of the matrices. In the upper panel in Figure 3, this region appears within the range between approximately 0.3 and $8\ \mathrm{cm}$ for the lenient case and between approximately 0.1 and $8\ \mathrm{cm}$ for the severe case. From our calculation, the collision velocity is given by Equation (12) and we can rewrite Equation (29) as

$$\begin{eqnarray}&&{L}_{\mathrm{pen}}\simeq \displaystyle \frac{{\pi }^{2}}{6}{\left(\displaystyle \frac{{v}_{0}}{{\rm{\Delta }}v}\right)}^{2}\displaystyle \frac{{\alpha }_{{\rm{t}}}{\varepsilon }^{2}{{\mathrm{Re}}_{{\rm{t}}}}^{1/2}{{c}_{{\rm{s}}}}^{2}}{{{{\rm{\Sigma }}}_{{\rm{g}}}}^{2}}\displaystyle \frac{{{\rho }_{0}}^{3}{{a}_{\mathrm{CA}}}^{2}}{{\chi }^{3}{\phi }_{\mathrm{CA}}}{\left(\displaystyle \frac{{E}_{\mathrm{roll}}}{{{a}_{0}}^{3}}\right)}^{-1}{a}_{\mathrm{ch}}.\end{eqnarray} \tag{ 32 }$$

Here, we use the fact that the volume of a fluffy CA is dominated by its matrices; therefore, the filling factor of the matrices ${\phi }_{\mathrm{mat}}$ is nearly equal to $\chi {\phi }_{\mathrm{CA}}$ according to Equation (19). Then, the critical bulk density for keeping chondrules in their matrices, obtained from the requirement ${L}_{\mathrm{pen}}\lesssim {a}_{\mathrm{CA}}$, is given by

$$\begin{eqnarray}&&{\rho }_{\mathrm{CA}}\gtrsim \displaystyle \frac{{\pi }^{2}}{6}{\left(\displaystyle \frac{{v}_{0}}{{\rm{\Delta }}v}\right)}^{2}\displaystyle \frac{{\alpha }_{{\rm{t}}}{\varepsilon }^{2}{{\mathrm{Re}}_{{\rm{t}}}}^{1/2}{{c}_{{\rm{s}}}}^{2}}{{{{\rm{\Sigma }}}_{{\rm{g}}}}^{2}}\displaystyle \frac{{{\rho }_{0}}^{4}{a}_{\mathrm{ch}}}{{\chi }^{3}}{\left(\displaystyle \frac{{E}_{\mathrm{roll}}}{{{a}_{0}}^{3}}\right)}^{-1}{a}_{\mathrm{CA}}.\end{eqnarray} \tag{ 33 }$$

As in Section 4.1, we introduce the "Stokes-drag region" as the region in which the CA affects the Stokes drag and the chondrules are decelerated by the compressive stress of the matrices. In the upper panel in Figure 3, this region appears within the range from approximately $8\ \mathrm{cm}$ to near the summit of the ejection barrier for both the lenient and severe cases. The collision velocity is given by Equation (12), the penetration length is

$$\begin{eqnarray}&&{L}_{\mathrm{pen}}\simeq \displaystyle \frac{16\pi }{243}{\left(\displaystyle \frac{{v}_{0}}{{\rm{\Delta }}v}\right)}^{2}\displaystyle \frac{{\varepsilon }^{2}{{\mathrm{Re}}_{{\rm{t}}}}^{5/2}{{{\rm{\Omega }}}_{{\rm{K}}}}^{2}}{{\alpha }_{{\rm{t}}}{{{\rm{\Sigma }}}_{{\rm{g}}}}^{2}}\displaystyle \frac{{{\rho }_{0}}^{3}{{a}_{\mathrm{CA}}}^{4}}{{\chi }^{3}{\phi }_{\mathrm{CA}}}{\left(\displaystyle \frac{{E}_{\mathrm{roll}}}{{{a}_{0}}^{3}}\right)}^{-1}{a}_{\mathrm{ch}},\end{eqnarray} \tag{ 34 }$$

and the critical bulk density for retaining the chondrules in their matrices is given by

$$\begin{eqnarray}&&{\rho }_{\mathrm{CA}}\gtrsim \displaystyle \frac{16\pi }{243}{\left(\displaystyle \frac{{v}_{0}}{{\rm{\Delta }}v}\right)}^{2}\displaystyle \frac{{\varepsilon }^{2}{{\mathrm{Re}}_{{\rm{t}}}}^{5/2}{{{\rm{\Omega }}}_{{\rm{K}}}}^{2}}{{\alpha }_{{\rm{t}}}{{{\rm{\Sigma }}}_{{\rm{g}}}}^{2}}\displaystyle \frac{{{\rho }_{0}}^{4}{a}_{\mathrm{ch}}}{{\chi }^{3}}{\left(\displaystyle \frac{{E}_{\mathrm{roll}}}{{{a}_{0}}^{3}}\right)}^{-1}{{a}_{\mathrm{CA}}}^{3}.\end{eqnarray} \tag{ 35 }$$

The Stokes-drag region does not appear in the severe cases in Figures 4(a)–(c). This is because the initial velocity of penetration v₀ already exceeds $\sqrt{B/A}$ before the CAs enter the Stokes-drag region when they grow along the critical ejection lines in radius–density space.

Note that the relative velocity ${\rm{\Delta }}v$ changes ${\rm{\Delta }}{v}_{{\rm{t}}}^{{\rm{S}}}=\varepsilon {{\mathrm{Re}}_{{\rm{t}}}}^{1/4}\mathrm{St}{v}_{{\rm{L}}}$ (Equation (12)) into ${\rm{\Delta }}{v}_{{\rm{t}}}^{{\rm{M}}}=1.4{\mathrm{St}}^{1/2}{v}_{{\rm{L}}}$ (Equation (13)) when the Stokes number reaches $\mathrm{St}={(\varepsilon /1.4)}^{-2}{{\mathrm{Re}}_{{\rm{t}}}}^{-1/2}$, and the critical Stokes number for $\varepsilon =0.1$ is on the order of ${10}^{-2};$ however, the critical Stokes number strongly depends on the uncertain parameter ε. If ε is close to unity, the relative velocity of the CAs is given by Equation (13), and the penetration length is given by

$$\begin{eqnarray}&&{L}_{\mathrm{pen}}\simeq 0.73{\left(\displaystyle \frac{{v}_{0}}{{\rm{\Delta }}v}\right)}^{2}\displaystyle \frac{{\mathrm{Re}}_{{\rm{t}}}{c}_{{\rm{s}}}{{\rm{\Omega }}}_{{\rm{K}}}}{{{\rm{\Sigma }}}_{{\rm{g}}}}\displaystyle \frac{{{\rho }_{0}}^{2}{{a}_{\mathrm{CA}}}^{2}}{{\chi }^{3}{{\phi }_{\mathrm{CA}}}^{2}}{\left(\displaystyle \frac{{E}_{\mathrm{roll}}}{{{a}_{0}}^{3}}\right)}^{-1}{a}_{\mathrm{ch}}.\end{eqnarray} \tag{ 36 }$$

In this case, the critical bulk density is

$$\begin{eqnarray}&&{\rho }_{\mathrm{CA}}\gtrsim 0.85\displaystyle \frac{{v}_{0}}{{\rm{\Delta }}v}\sqrt{\displaystyle \frac{{\mathrm{Re}}_{{\rm{t}}}{c}_{{\rm{s}}}{{\rm{\Omega }}}_{{\rm{K}}}}{{{\rm{\Sigma }}}_{{\rm{g}}}}}\displaystyle \frac{{{\rho }_{0}}^{2}{{a}_{\mathrm{ch}}}^{1/2}}{{\chi }^{3/2}}{\left(\displaystyle \frac{{E}_{\mathrm{roll}}}{{{a}_{0}}^{3}}\right)}^{-1/2}{{a}_{\mathrm{CA}}}^{1/2}.\end{eqnarray} \tag{ 37 }$$

Finally, we define the "hydrodynamic region" as the region in which chondrules in a CA are decelerated by the ram pressure caused by high-speed penetration. The penetration length is given by Equation (28), and we assume that the term $\mathrm{ln}({v}_{0}\sqrt{A/B})$ is on the order of unity; then, we can rewrite Equation (28) as

$$\begin{eqnarray}&&{L}_{\mathrm{pen}}\sim \displaystyle \frac{8\chi }{3{\phi }_{\mathrm{CA}}}{a}_{\mathrm{ch}}.\end{eqnarray} \tag{ 38 }$$

The critical bulk density for retaining chondrules is

$$\begin{eqnarray}&&{\rho }_{\mathrm{CA}}\gtrsim \displaystyle \frac{8\chi {\rho }_{0}{a}_{\mathrm{ch}}}{3}{{a}_{\mathrm{CA}}}^{-1}.\end{eqnarray} \tag{ 39 }$$

The critical bulk density depends little on the properties of matrix grains in the hydrodynamic region, while the properties of matrix grains (i.e., ${E}_{\mathrm{roll}}/{{a}_{0}}^{3}$) strongly affect the critical bulk density in the Epstein-drag and Stokes-drag regions.

In all regions, the penetration length is proportional to the radius of the chondrules. We assume ${a}_{\mathrm{ch}}=0.01\ \mathrm{cm}$ in this study, however, the sizes of chondrules vary between different types of chondrites. The mean radius of chondrules in H ordinary chondrites and EH enstatite chondrites is approximately $0.01\ \mathrm{cm}$, whereas the chondrules in L and LL ordinary chondrites and EL ordinary chondrites are approximately $0.03\ \mathrm{cm}$ in radii (Rubin 2000, and references therein). By contrast, the spread of sizes in individual groups is relatively small. Approximately 70% of chondrules in Qingzhen, Kota-Kota, and Allan Hills A77156 EH3 chondrites have similar radii within a factor of two (Rubin & Grossman 1987). It might be possible to derive some constraints on the accretion environment of CAs from the differences in chondrule size in the different groups of chondrites. We will discuss this topic in the future.

Recent high-precision measurements of chondrules, matrices, and bulk chondrites reveal that chondrites have chemical (e.g., Palme et al. 2015; Ebel et al. 2016) and isotopic (e.g., Budde et al. 2016a, 2016b) complementarities. These complementarities are therefore a strong constraint, indicating that the chondrules and the matrices must have formed from a single reservoir and that after their formation, neither the chondrules nor the matrix grains were lost. The simplest way to explain the chondrule-matrix complementarity is that the chondrules and the matrix grains were formed in the same heating events and that some parts of the matrix grains are condensates of evaporated dust. Miura et al. (2010) showed that a shock wave heating model for chondrule formation can predict the evaporation and formation of fine dust grains via chondrule formation; in addition, Miura & Nakamoto (2005) revealed that the minimum size of chondrules that avoids the complete evaporation is consistent with the observations (e.g., Eisenhour 1996; Nelson & Rubin 2002).

There are many models for chondrule formation, and these models are classified into two groups, that is, impact origin models (e.g., Sanders & Scott 2012; Dullemond et al. 2014; Wakita et al. 2017) and nebular origin models (e.g., Muranushi 2010; Boley et al. 2013; McNally et al. 2013). To determine the formation process(es) of chondrules, the key constraints are the heating and cooling timescales. According to Tachibana & Huss (2005), the chondrule-forming events must be flash heating events, i.e., a heating rate of 10⁴–${10}^{6}\ {\rm{K}}\ {{\rm{h}}}^{-1}$ would be required to avoid the isotopic fractionation of sulfur in chondrule precursors. In addition, although there are huge uncertainties, several studies suggest that the cooling rate of chondrules is as large as 10⁵–${10}^{6}\ {\rm{K}}\ {{\rm{h}}}^{-1}$ from observations of iron–magnesium and oxygen isotopic exchange (e.g., Yurimoto & Wasson 2002) and the crystal growth of olivine phenocrysts (e.g., Wasson & Rubin 2003; Miura & Yamamoto 2014). These estimates are larger by several orders of magnitude than the cooling rates of 10–${10}^{3}\ {\rm{K}}\ {{\rm{h}}}^{-1}$ obtained from furnace-faced experiments (Desch et al. 2012, and references therein).

Even though it is true that furnace-faced experiments have succeeded in reproducing some textural features of chondrules (e.g., Hewins & Fox 2004; Tsuchiyama et al. 2004), the solidification processes might be affected by the contact with a furnace wall (Nagashima et al. 2006). Furthermore, rapid cooling of levitating supercooled precursors also succeeds in reproducing some textures of chondrules (e.g., Srivastava et al. 2010). The supercooling of chondrule melts is preferable to explain the textural features of compound chondrules, which might form via the collision of chondrule precursors (Arakawa & Nakamoto 2016a). The rapid cooling of the chondrules indicates that the evaporated dust would also be cooled rapidly, and the cooling rate of the dust vapor affects the size and the morphology of the condensates (e.g., Miura et al. 2010).

Understanding of the physical properties of the matrix grains is important when considering the growth process of CAs. As shown in Section 4, the penetration length ${L}_{\mathrm{pen}}$ is inversely proportional to ${E}_{\mathrm{roll}}/{{a}_{0}}^{3}$, and the rolling energy of two sticking matrix grains ${E}_{\mathrm{roll}}$ is ${E}_{\mathrm{roll}}=6{\pi }^{2}\gamma {a}_{0}\xi$. Therefore, we can obtain the relationship between the surface energy, the radius, and the penetration length, ${L}_{\mathrm{pen}}\propto {\gamma }^{-1}{{a}_{0}}^{2}$, and large surface energies and small grain sizes are preferable when the bulk density is fixed. We append that if the bulk density of the CAs is determined by static compression, ${\rho }_{\mathrm{CA}}$ also depends on ${E}_{\mathrm{roll}}/{a}_{0}^{3};$ however, the dependence is weak (e.g., Kataoka et al. 2013b; Arakawa & Nakamoto 2016b).

Here we mention the fact that the exact value of ${E}_{\mathrm{roll}}$ might differ from $6{\pi }^{2}\gamma {a}_{0}\xi$. This is because the equation ${E}_{\mathrm{roll}}=6{\pi }^{2}\gamma {a}_{0}\xi$ obtained by Dominik & Tielens (1995) is based on the so-called JKR contact model (Johnson et al. 1971), and the JKR model is valid only for large values of the Tabor parameter μ (Tabor 1977):

$$\begin{eqnarray}&&\mu ={\left(\displaystyle \frac{2{a}_{0}{\gamma }^{2}{(1-\nu )}^{2}}{{Y}^{2}{{z}_{0}}^{3}}\right)}^{1/3}=0.11\ {\left(\displaystyle \frac{{a}_{0}}{0.025\mu {\rm{m}}}\right)}^{1/3},\end{eqnarray} \tag{ 40 }$$

where $\nu =0.17$ is the Poisson's ratio, $Y=5.4\times {10}^{11}\ \mathrm{dyn}\ {\mathrm{cm}}^{-2}$ is the Young's modulus, and ${z}_{0}=0.2\ \mathrm{nm}$ is the interatomic distance. The JKR model can be used only for $\mu \gtrsim 5$ (Johnson & Greenwood 1997), and thus other contact models should be used for submicron-sized silicate spheres. In addition, the dependence on the sphere radius is weak, thus we should not adopt the JKR theory not only for submicron-sized spheres but also micron- or millimeter-sized silicate grains. Even if we consider chondrule-sized grains, the Tabor parameter is still small; the Tabor parameter is $\mu =1.7$ when the sphere radius is ${a}_{0}=0.01\ \mathrm{cm}$. The Maugis-Dugdale model (Maugis 1992) might be better than the JKR model where the Tabor parameter is $0.1\lesssim \mu \lesssim 5$, and several contact models such as the so-called DMT model (Derjaguin et al. 1975) and the semi-rigid sphere model (Greenwood 2007) are suggested for $\mu \lesssim 0.1$. However, the maximum force needed to separate the two contact spheres ${F}_{{\rm{c}}}$ is hardly dependent on the selection of contact models: ${F}_{{\rm{c}}}=(3/2)\pi \gamma {a}_{0}$ for the JKR model and ${F}_{{\rm{c}}}=2\pi \gamma {a}_{0}$ for the DMT model. Additionally, the rolling energy ${E}_{\mathrm{roll}}$ is proportional to ${F}_{{\rm{c}}}\xi$ (e.g., ${E}_{\mathrm{roll}}=4\pi {F}_{{\rm{c}}}\xi$ for the JKR model), therefore the rolling energy might be also subequal among contact models. Either way, future studies of the rolling resistance between two stiff spheres are necessary.

The morphology of the matrix grains also influences both the penetration process and the density evolution because the irregular shape of the monomers affects the rolling energy. In addition, the irregular shape affects not only the density and the penetration process but also the critical velocity for collisional growth (Poppe et al. 2000). It has been suggested that at least some parts of the matrix grains are condensates of evaporated dust and that the size of the condensate ranges from nanometer to micrometer (e.g., Yamamoto & Hasegawa 1977; Miura et al. 2010) depending on the cooling rate and the mass density of the dust vapor. Ishizuka et al. (2015) experimentally revealed that there are two types of nanometer-sized silicate condensates obtained from the re-condensation of magnesium and silicon monoxide vapor, that is, spherical forsterite particles and cubic MgO particles. However, there is no reliable contact theory that can be applied to non-spherical monomers. Therefore, improvements on the contact mechanical theory for non-spherical monomers are also necessary.

In this study, we assume the uniform density of the matrix regions throughout the CA for simplicity. However, the gradient and fluctuation of the density of the matrix regions of a CA might exist to some extent. Considering the accretion process in Stage II, the initial growth of chondrules/CAs is driven by the accretion of small and dense MAs. If the dense dust rims around chondrules can be maintained when CAs collide at a large velocity, dust-rimmed chondrules in the fluffy CA are subjected to a large drag force ${F}_{\mathrm{pen}}$ and the penetration length ${L}_{\mathrm{pen}}$ of dust-rimmed chondrules might be smaller than the penetration length of naked chondrules. Although it is still unclear whether chondrules can maintain the dense dust rims or not, we should consider the formation and influence of the dust rims around chondrules.

In this paper, we only consider the drag force for penetration ${F}_{\mathrm{pen}}$ as a decelerating mechanism. Even though the narrow size distribution of chondrules (e.g., Nelson & Rubin 2002) might prevent chondrules from colliding unless there is a large variation in the initial velocity of penetration, chondrules in a large CA are likely affected by chondrule collisions when they move in the CA. In addition, the recapture of chondrules by large CAs is likely affected by this process because there is a large velocity difference between intruding chondrules and chondrules contained in CAs. We will examine this effect in the future.

As shown in Section 3, CAs would lose their chondrules if the bulk density of these aggregates were determined by static compression processes. Therefore, we need to consider alternative scenarios for chondritic planetesimal formation. Here, we propose two scenarios: changing the density evolution pathway of CAs to retain chondrules in CAs and escaping from the ejection barrier of planetesimal formation via the streaming instability before CAs lose their chondrules.

5.5.1. Modification of the Density Evolution Pathway

5.5.2. Jumping over the Ejection Barrier