Tensile Strength of Porous Dust Aggregates

We calculate interactions of each connection between two monomers using the theoretical model by Dominik & Tielens (1997) and Wada et al. (2007). Based on the JKR theory (Johnson et al. 1971) and the following studies by Dominik & Tielens (1995, 1996), Dominik & Tielens (1997) carried out two-dimensional simulations of monomer interactions. To expand into three-dimensional simulations, Wada et al. (2007) tested their recipe, and then Wada et al. (2008) conducted three-dimensional simulations of dust aggregate collisions. In the model, there are four kinds of interactions: normal (sticking and breaking), sliding, rolling, and twisting. The material parameters needed to describe the interactions are the monomer radius r₀, material density ρ₀, surface energy γ, Poisson's ratio ν, Young's modulus E, and the critical rolling displacement ξ_crit. These parameters of ice and silicate are listed in Table 1. To compare our results with those by Seizinger et al. (2013), we set the same values for silicate.

If a rolling displacement exceeds the critical one, ξ_crit, a monomer begins to roll inelastically. The critical rolling displacement has different values between the theoretical one (ξ_crit = 2 Å, Dominik & Tielens 1997) and the experimental one (ξ_crit = 32 Å, Heim et al. 1999). We adopt ξ_crit = 8 Å as a fiducial value and investigate the dependence of our results on ξ_crit in Section 3.3.

The rolling energy E_roll needed to rotate a monomer around its connection point by 90° is described as

$$\begin{eqnarray}\begin{array}{rcl}{E}_{\mathrm{roll}} & = & 12{\pi }^{2}\gamma R{\xi }_{\mathrm{crit}}\\ & = & 6{\pi }^{2}\gamma {r}_{0}{\xi }_{\mathrm{crit}}\\ & \sim & 4.7\times {10}^{-16}\left(\displaystyle \frac{\gamma }{100\,\mathrm{mJ}\,{{\rm{m}}}^{-2}}\right)\\ & & \times \,\left(\displaystyle \frac{{r}_{0}}{0.1\,\mu {\rm{m}}}\right)\left(\displaystyle \frac{{\xi }_{\mathrm{crit}}}{8\mathring{\rm A} }\right)\,{\rm{J}},\end{array}\end{eqnarray} \tag{ 1 }$$

where R is the reduced monomer radius (Wada et al. 2007). The reduced radius R of monomer radii r₁ and r₂ is defined as

$$\begin{eqnarray}&&\displaystyle \frac{1}{R}=\displaystyle \frac{1}{{r}_{1}}+\displaystyle \frac{1}{{r}_{2}}.\end{eqnarray} \tag{ 2 }$$

Here, the reduced monomer radius is R = r₀/2 because we assume no size distribution of monomers.

In our simulations, the maximum force needed to separate two sticking monomers (breaking) is

$$\begin{eqnarray}\begin{array}{rcl}{F}_{{\rm{c}}} & = & 3\pi \gamma R\\ & \sim & 4.7\times {10}^{-8}\left(\displaystyle \frac{\gamma }{100\,\mathrm{mJ}\,{{\rm{m}}}^{-2}}\right)\left(\displaystyle \frac{{r}_{0}}{0.1\,\mu {\rm{m}}}\right)\,{\rm{N}}.\end{array}\end{eqnarray} \tag{ 3 }$$

The force in the normal direction induces oscillation at each connection between two monomers. In reality, the oscillation would attenuate because of viscoelasticity or hysteresis of monomers (Greenwood & Johnson 2006; Tanaka et al. 2012). Therefore, we add an artificial damping force in the normal direction (Paszun & Dominik 2008; Suyama et al. 2008; Seizinger et al. 2012; Kataoka et al. 2013b). The dependence of our results on the damping force is investigated in Section 3.2.

We describe the damping force as follows. In the case in which two contacting monomers have their position vectors ${{\boldsymbol{x}}}_{1}$ and ${{\boldsymbol{x}}}_{2}$, and velocities ${{\boldsymbol{v}}}_{1}$ and ${{\boldsymbol{v}}}_{2}$, respectively, the contact unit vector ${{\boldsymbol{n}}}_{{\rm{c}}}$ is defined as

$$\begin{eqnarray}&&{{\boldsymbol{n}}}_{{\rm{c}}}=\displaystyle \frac{{{\boldsymbol{x}}}_{1}-{{\boldsymbol{x}}}_{2}}{| {{\boldsymbol{x}}}_{1}-{{\boldsymbol{x}}}_{2}| }\end{eqnarray} \tag{ 4 }$$

(Wada et al. 2007). The damping force applied to each monomer is introduced as

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{\mathrm{damp}}=-{k}_{{\rm{n}}}\displaystyle \frac{{m}_{0}}{{t}_{{\rm{c}}}}({{\boldsymbol{n}}}_{{\rm{c}}}\cdot {{\boldsymbol{v}}}_{{\rm{r}}}){{\boldsymbol{n}}}_{{\rm{c}}},\end{eqnarray} \tag{ 5 }$$

where k_n is the damping coefficient, m₀ is the monomer mass, t_c is the characteristic time, and ${{\boldsymbol{v}}}_{{\rm{r}}}$ is the relative velocity (Kataoka et al. 2013b). When we calculate the damping force experienced by the monomer (${{\boldsymbol{x}}}_{1},{{\boldsymbol{v}}}_{1}$), ${{\boldsymbol{v}}}_{{\rm{r}}}={{\boldsymbol{v}}}_{2}-{{\boldsymbol{v}}}_{1}$ is the relative velocity of the other monomer (${{\boldsymbol{x}}}_{2},{{\boldsymbol{v}}}_{2}$). We adopt k_n = 0.01 as a fiducial value.

The characteristic time is given by Wada et al. (2007) as

$$\begin{eqnarray}&&{t}_{{\rm{c}}}=0.95\left(\displaystyle \frac{{r}_{0}^{7/6}{\rho }_{0}^{1/2}}{{\gamma }^{1/6}{{E}^{* }}^{1/3}}\right),\end{eqnarray} \tag{ 6 }$$

where E* is the reduced Young's modulus of monomers 1 and 2 defined as

$$\begin{eqnarray}&&\displaystyle \frac{1}{{E}^{* }}=\displaystyle \frac{1-{\nu }_{1}^{2}}{{E}_{1}}+\displaystyle \frac{1-{\nu }_{2}^{2}}{{E}_{2}}.\end{eqnarray} \tag{ 7 }$$

In our simulation, the Young's modulus E₁ = E₂ = E and the Poisson's ratio ν₁ = ν₂ = ν are uniform.

The initial dust aggregates are statically and isotropically compressed ballistic cluster–cluster aggregations (BCCAs) investigated by Kataoka et al. (2013b). We set these initial conditions to simulate the planetesimal formation mechanism. The calculation boundary is treated periodically, thus we do not have to consider the aggregate radius.

We set moving boundaries in the x-axis direction and fixed boundaries in the y- and z-axis directions to measure the tensile strength of dust aggregates. The initial calculation box is a cube whose length on each side is L₀. The length in the y- and z-axis directions does not change, while the length in the x-axis direction L_x increases. Therefore, the coordinates in the x-, y-, and z-axis directions are $-{L}_{x}/2\lt x\lt {L}_{x}/2$, $-{L}_{0}/2\lt y\lt {L}_{0}/2$, and $-{L}_{0}/2\lt z\lt {L}_{0}/2$, respectively.

The velocity at the boundary in the x-axis direction ${v}_{{\rm{b}}}\gt 0$ has to be constant and less than the effective sound speed of dust aggregates for static stretching. We investigate the dependence on the velocity in Section 3.2. The effective sound speed of dust aggregates ${c}_{{\rm{s}},\mathrm{eff}}$ is described as

$$\begin{eqnarray}&&{c}_{{\rm{s}},\mathrm{eff}}\sim \sqrt{\displaystyle \frac{P}{\rho }},\end{eqnarray} \tag{ 8 }$$

where P and ρ are the pressure and mean internal density of dust aggregates, respectively. Because the initial dust aggregates are statically and isotropically compressed, their pressure is given as

$$\begin{eqnarray}&&P=\displaystyle \frac{{E}_{\mathrm{roll}}}{{r}_{0}^{3}}{\left(\displaystyle \frac{\rho }{{\rho }_{0}}\right)}^{3}\end{eqnarray} \tag{ 9 }$$

(Kataoka et al. 2013b). From Equations (8) and (9), we can obtain the effective sound speed of dust aggregates as

$$\begin{eqnarray}\begin{array}{rcl}{c}_{{\rm{s}},\mathrm{eff}} & \sim & \sqrt{\displaystyle \frac{{E}_{\mathrm{roll}}}{{\rho }_{0}{r}_{0}^{3}}}\displaystyle \frac{\rho }{{\rho }_{0}}\\ & \sim & 2.2\times {10}^{3}{\left(\displaystyle \frac{{r}_{0}}{0.1\mu {\rm{m}}}\right)}^{-1}{\left(\displaystyle \frac{{\rho }_{0}}{1.0{\rm{g}}{\mathrm{cm}}^{-3}}\right)}^{-1/2}\\ & & \times \,{\left(\displaystyle \frac{\gamma }{100\mathrm{mJ}{{\rm{m}}}^{-2}}\right)}^{1/2}{\left(\displaystyle \frac{{\xi }_{\mathrm{crit}}}{8\mathring{\rm A} }\right)}^{1/2}\phi \,\mathrm{cm}\,{{\rm{s}}}^{-1},\end{array}\end{eqnarray} \tag{ 10 }$$

where ϕ = ρ/ρ₀ is the volume filling factor of dust aggregates. Since v_b is independent of time t, the length in the x-axis direction L_x can be written as

$$\begin{eqnarray}&&{L}_{x}={L}_{0}+2{v}_{{\rm{b}}}t.\end{eqnarray} \tag{ 11 }$$

We treat the coordinates and velocity of a monomer across a periodic boundary as follows. When a monomer passes the moving periodic boundary at x = L_x/2, its position x and velocity v_x are converted as

$$\begin{eqnarray}&&x\to x-{L}_{x},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{v}_{x}\to {v}_{x}-2{v}_{{\rm{b}}}.\end{eqnarray} \tag{ 13 }$$

On the other hand, in the case of the moving periodic boundary at x = −L_x /2, its position and velocity are converted as

$$\begin{eqnarray}&&x\to x+{L}_{x},\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{v}_{x}\to {v}_{x}+2{v}_{{\rm{b}}}.\end{eqnarray} \tag{ 15 }$$

At y = ±L₀/2 and z = ±L₀/2, the coordinates of a monomer are converted similarly, but its velocity is not changed.

We calculate tensile stress in the same way as Kataoka et al. (2013b) because we have no walls. This is different from Seizinger et al. (2013), who measured the tensile stress considering the force exerted on walls.

The tensile stress is calculated only in the x-axis direction as follows. At first, we assume a virtual box, which is the same as the calculation box. The equation of motion of the monomer i in the x-axis direction is described as

$$\begin{eqnarray}&&{m}_{0}\displaystyle \frac{{d}^{2}{x}_{i}}{{{dt}}^{2}}={W}_{x,i}+{F}_{x,i},\end{eqnarray} \tag{ 16 }$$

where W_x,i is the force exerted from the walls of the virtual box on the monomer i and F_x,i is the total force from other monomers on the monomer i. We multiply Equation (16) by x_i and take a long-time average with the time interval τ. The left side of Equation (16) becomes

$$\begin{eqnarray}&&\displaystyle \frac{{m}_{0}}{\tau }{\int }_{0}^{\tau }{x}_{i}\displaystyle \frac{{d}^{2}{x}_{i}}{{{dt}}^{2}}{dt}=\displaystyle \frac{{m}_{0}}{\tau }{\left[{x}_{i}\displaystyle \frac{{{dx}}_{i}}{{dt}}\right]}_{0}^{\tau }-\displaystyle \frac{{m}_{0}}{\tau }{\int }_{0}^{\tau }\displaystyle \frac{{{dx}}_{i}}{{dt}}\displaystyle \frac{{{dx}}_{i}}{{dt}}{dt}.\end{eqnarray} \tag{ 17 }$$

The first term on the right side of Equation (17) becomes zero when $\tau \to \infty$. Here, we define the long-time average as $\langle {\rangle }_{t}$ and take a summation of all monomers of Equation (16). Then, Equation (16) can be written as

$$\begin{eqnarray}&&{\left\langle \displaystyle \sum _{i=1}^{N}\displaystyle \frac{{m}_{0}}{2}{\left(\displaystyle \frac{{{dx}}_{i}}{{dt}}\right)}^{2}\right\rangle }_{t}=-\displaystyle \frac{1}{2}{\left\langle \displaystyle \sum _{i=1}^{N}{x}_{i}{W}_{x,i}\right\rangle }_{t}-\displaystyle \frac{1}{2}{\left\langle \displaystyle \sum _{i=1}^{N}{x}_{i}{F}_{x,i}\right\rangle }_{t}.\end{eqnarray} \tag{ 18 }$$

The left side of Equation (18) can be defined as

$$\begin{eqnarray}&&{K}_{x}={\left\langle \displaystyle \sum _{i=1}^{N}\displaystyle \frac{{m}_{0}}{2}{\left(\displaystyle \frac{{{dx}}_{i}}{{dt}}\right)}^{2}\right\rangle }_{t},\end{eqnarray} \tag{ 19 }$$

which is the time-averaged kinematic energy in the x-axis direction of all monomers. The first energy term on the right side of Equation (18) is related to the tensile stress in the x-axis direction P_x. Since the virtual box is the same as the calculation box, we obtain

$$\begin{eqnarray}&&{\left\langle \displaystyle \sum _{i=1}^{N}{x}_{i}{W}_{x,i}\right\rangle }_{t}={L}_{x}{P}_{x}{L}_{0}^{2}={P}_{x}V,\end{eqnarray} \tag{ 20 }$$

where $V={L}_{x}{L}_{0}^{2}$ is the volume of the calculation box. Therefore, Equation (18) gives an expression of the tensile stress P_x as

$$\begin{eqnarray}&&{P}_{x}=-\displaystyle \frac{2{K}_{x}}{V}-\displaystyle \frac{1}{V}{\left\langle \displaystyle \sum _{i=1}^{N}{x}_{i}{F}_{x,i}\right\rangle }_{t}.\end{eqnarray} \tag{ 21 }$$

The total force from other monomers on the monomer i can be described as

$$\begin{eqnarray}&&{F}_{x,i}=\displaystyle \sum _{j\ne i}{f}_{x,i,j},\end{eqnarray} \tag{ 22 }$$

where f_x,i,j is the force from the monomer j on the monomer i in the x-axis direction. Thus, Equation (21) can be written as

$$\begin{eqnarray}&&{P}_{x}=-\displaystyle \frac{2{K}_{x}}{V}-\displaystyle \frac{1}{V}{\left\langle \displaystyle \sum _{i\lt j}({x}_{i}-{x}_{j}){f}_{x,i,j}\right\rangle }_{t}\end{eqnarray} \tag{ 23 }$$

because of the relation that f_x,i,j = −f_x,j,i. Equation (23) is different from that of Kataoka et al. (2013b) because we consider the tensile stress only in the x-axis direction.

We take an average of the tensile stress P_x for every 10,000 time steps, at least. In some simulations, the tensile stress fluctuates, thus we take a longer time average to smooth it (see Section 3.1 for details). One time step in our simulation is $0.7{t}_{{\rm{c}}}=1.9\times {10}^{-11}$ s, and therefore 10,000 time steps corresponds to $1.9\times {10}^{-7}$ s.

The overview of our numerical simulations is as follows. First, we randomly create a BCCA to change the initial condition. Next, we compress it statically and isotropically (Kataoka et al. 2013b). The compression of a BCCA corresponds to the formation of a planetesimal. It is necessary for the BCCA to be attached to all boundaries so that we can stretch it. Then, we stop compression and define this volume filling factor as the initial one ϕ_init. Finally, we stretch it statically and one-dimensionally. Figure 1 shows the overview of our simulations.

Zoom In Zoom Out Reset image size

We measure the tensile stress of 10 runs for the fiducial parameter set. The fiducial values are N = 16384, ${v}_{{\rm{b}}}=10\,\mathrm{cm}\,{{\rm{s}}}^{-1}$, k_n = 0.01, ϕ_init = 0.1, r₀ = 0.1 μm, γ = 100 mJ m⁻², and ξ_crit = 8 Å. Figure 2 shows three snapshots of a fiducial run. Each particle represents a 0.1 μm radius ice monomer. The light gray monomers are in the calculation box with periodic boundaries, while the dark gray monomers are in the neighbor boxes. The box with white lines shows the final state of the calculation box. By stretching the dust aggregate, the chain-like structure appears.

Zoom In Zoom Out Reset image size

Figure 3 shows the time evolution of tensile stress of 10 fiducial runs averaged for every 10,000 time steps (left) and 200,000 time steps (right). The volume filling factor at each time step is calculated as

$$\begin{eqnarray}&&\phi =\displaystyle \frac{(4/3)\pi {r}_{0}^{3}N}{V}\end{eqnarray} \tag{ 24 }$$

and the tensile stress is calculated according to Equation (23). We choose the number of time steps when the rate of change of the volume filling factor does not exceed 10%. As tensile displacement increases, the volume filling factor ϕ decreases and the tensile stress P_x increases. The maximum value of tensile stress is called the tensile strength. To calculate the tensile strength for every parameter set, we find 10 maximum values of the obtained tensile stress and take an average of them.

Zoom In Zoom Out Reset image size

To investigate the dependence on the number of particles N, we plot tensile stress when N = 2¹⁰ = 1024, N = 2¹² = 4096, N = 2¹⁴ = 16384, and N = 2¹⁶ = 65536 in Figure 4(a). Changing N corresponds to changing the size of the calculation box. Obviously, the tensile strength has no dependence on N. Because of the smoothness of the tensile stress plot and calculation costs, we set N = 16384 as the fiducial value.

Zoom In Zoom Out Reset image size

Figure 4(b) shows tensile stress when boundary velocity v_b = 1, 10, and 100 cm s⁻¹. All boundary velocities are less than the effective sound speed of dust aggregates c_s,eff (Equation (10)). There is no difference among the three boundary velocities. Therefore, we can conclude that dust aggregates are stretched statically. We set v_b = 10 cm s⁻¹ as the fiducial value considering sampling rates of tensile stress and calculation costs.

Tensile stress with various damping coefficients is plotted in Figure 4(c). No damping force corresponds to k_n = 0. We change the strength of damping force from weak damping (k_n = 0.01) to strong damping (k_n = 1). Undoubtedly, there is no dependence on the damping force in this range. We use k_n = 0.01 for all the other simulations.

We measure the tensile strength of dust aggregates that have various initial volume filling factors as shown in Figure 5. The calculated tensile strength is proportional to ${\phi }_{\mathrm{init}}^{1.8}$ from the fitting. The tensile strength P_x,max can be described with the initial volume filling factor ϕ_init as

$$\begin{eqnarray}&&{P}_{x,\max }={P}_{0}{\phi }_{\mathrm{init}}^{1.8},\end{eqnarray} \tag{ 25 }$$

where ${P}_{0}\sim 6\times {10}^{5}\,\mathrm{Pa}$ in this case. The analytical interpretation of Equation (25) is discussed in Section 4.1.

Zoom In Zoom Out Reset image size

To investigate the dependence of tensile strength on the monomer radius, we perform simulations in the case of ice monomers whose radii are 0.3 and 0.9 μm. Figure 6 shows the summary of the monomer radius dependence. The plotted dashed lines are based on Equation (31), which is an analytical expression of tensile strength (see Section 4.1). It is confirmed that the tensile strength is in inverse proportion to the monomer radius.

Zoom In Zoom Out Reset image size

To clarify the surface energy dependence, we calculate the tensile strength when γ = 50 and 25 mJ m⁻² and plot it in Figure 7. The other parameters are the same as the fiducial values. The dashed lines represent Equation (31) (see Section 4.1). Obviously, the tensile strength is in proportion to the surface energy.

Zoom In Zoom Out Reset image size

Finally, we investigate the dependence of tensile stress on the critical rolling displacement ξ_crit in Figure 8. The critical rolling displacement is changed from ${\xi }_{\mathrm{crit}}=2\mathring{\rm A}$ (e.g., Dominik & Tielens 1997) to ${\xi }_{\mathrm{crit}}=32\mathring{\rm A}$ (e.g., Heim et al. 1999). Tensile stress has a marginal dependence on ξ_crit because the main mechanism of displacement is rolling (see Section 4.1). On the other hand, tensile strength, which is the maximum value of tensile stress, has no difference. We can conclude that the tensile strength is almost the same even if the critical rolling displacement has uncertainty. Therefore, we fix ${\xi }_{\mathrm{crit}}=8\mathring{\rm A}$ in our simulations.

Zoom In Zoom Out Reset image size

The relationship between P_x,max and ϕ_init can be derived by considering the maximum force needed to separate two sticking monomers F_c and the radius of a dust aggregate r_agg. When the tensile stress has a maximum value, the force F_c is applied on a connection between two monomers of a dust aggregate. This means that

$$\begin{eqnarray}&&{P}_{x,\max }\propto \displaystyle \frac{{F}_{{\rm{c}}}}{{r}_{\mathrm{agg}}^{2}}.\end{eqnarray} \tag{ 26 }$$

The radius of a dust aggregate is given as

$$\begin{eqnarray}&&{r}_{\mathrm{agg}}\propto {N}_{\mathrm{agg}}^{1/D}{r}_{0},\end{eqnarray} \tag{ 27 }$$

where D and N_agg are the fractal dimension and the number of monomers of a dust aggregate, respectively. The initial volume filling factor is described as

$$\begin{eqnarray}&&{\phi }_{\mathrm{init}}={N}_{\mathrm{agg}}{\left(\displaystyle \frac{{r}_{0}}{{r}_{\mathrm{agg}}}\right)}^{3},\end{eqnarray} \tag{ 28 }$$

and then the radius of a dust aggregate is obtained as

$$\begin{eqnarray}&&{r}_{\mathrm{agg}}\propto {r}_{0}{\phi }_{\mathrm{init}}^{-1/(3-D)}.\end{eqnarray} \tag{ 29 }$$

From Equation (26), the tensile strength can be written as

$$\begin{eqnarray}&&{P}_{x,\max }\sim C\displaystyle \frac{{F}_{{\rm{c}}}}{{r}_{0}^{2}}{\phi }_{\mathrm{init}}^{2/(3-D)},\end{eqnarray} \tag{ 30 }$$

where C is a constant. The fractal dimension D of BCCAs is ∼1.9 (Mukai et al. 1992; Okuzumi et al. 2009).

Using the fitting result of Equations (3) and (25), we obtain

$$\begin{eqnarray}\begin{array}{rcl}{P}_{x,\max } & \sim & 6\times {10}^{5}\left(\displaystyle \frac{\gamma }{100\,\mathrm{mJ}\,{{\rm{m}}}^{-2}}\right)\\ & & \times \,{\left(\displaystyle \frac{{r}_{0}}{0.1\mu {\rm{m}}}\right)}^{-1}{\phi }_{\mathrm{init}}^{1.8}\,\mathrm{Pa}.\end{array}\end{eqnarray} \tag{ 31 }$$

In the case of an ice monomer with a radius of 0.1 μm, we find C = 0.12 ± 0.01.

We confirm Equation (31) from the perspective of energy dissipation. All energy dissipations, which are caused by the normal, sliding, rolling, twisting, and damping force in the normal direction, are plotted in Figure 9. The curves in Figure 9 run time-wise from right to left and arise during the stretching of a dust aggregate. The main energy dissipation mechanism is the rolling, which corresponds to the ξ_crit-dependence of tensile stress (Section 3.3). The energy dissipation by the normal arises when the tensile stress has a maximum value. This energy dissipation is caused by connection breaking between two contacting monomers. For this reason, tensile strength is determined by the connection breaking, i.e., F_c.

Zoom In Zoom Out Reset image size

To confirm that D ∼ 1.9 on a small scale of a dust aggregate in our simulations, we calculate the number of monomers inside the radius r_in for five snapshots of a fiducial run and plot it in Figure 10. We take the snapshots during the continuous strain of the dust aggregate. The parameters of the fiducial run are N = 16384, ${v}_{{\rm{b}}}=10\,\mathrm{cm}\,{{\rm{s}}}^{-1}$, k_n = 0.01, ϕ_init = 0.1, r₀ = 0.1 μm, γ = 100 mJ m⁻², and ξ_crit = 8 Å. The method to count the number of monomers N(r < r_in) is as follows. At first, we set a monomer in the calculation box as the center and count N(r < r_in) including monomers outside the periodic boundaries. Next, we take an average of N(r < r_in) for all monomers in the calculation box.

Zoom In Zoom Out Reset image size

Also, we plot N(r < r_in) as a function of ${r}_{\mathrm{in}}/{r}_{0}$ when D = 2 and D = 3 in Figure 10. This relationship is derived by Equation (27) as

$$\begin{eqnarray}&&N(r\lt {r}_{\mathrm{in}})\propto {\left(\displaystyle \frac{{r}_{\mathrm{in}}}{{r}_{0}}\right)}^{D}.\end{eqnarray} \tag{ 32 }$$

On a small scale that $N(r\lt {r}_{\mathrm{in}})\lesssim 20$, all snapshot results mostly correspond to the relationship when D = 2. When the scale becomes large enough, dust aggregates have a fractal dimension of three.

To compare our results with previous experiments and numerical simulations, we perform simulations with the parameters of silicate listed in Table 1 and summarize the results in Figure 11. Both results by experiments (Blum & Schräpler 2004; Blum et al. 2006; Gundlach et al. 2018) and simulations (this work; Seizinger et al. 2013) correspond very well. This means that there is little influence from the artificial adhesion force introduced by Seizinger et al. (2013) and small aggregates whose sizes range from micrometers to millimeters. The right panel of Figure 11 is derived by Equation (31), i.e., ${P}_{x,\max }\propto {r}_{0}^{-1}$.

Zoom In Zoom Out Reset image size

Gundlach et al. (2018) measured the tensile strength of ice aggregates with the monomer radius of $2.38\pm 1.11\,\mu {\rm{m}}$, which is shown in Figure 12. We do not know which monomer radius determines the tensile strength of dust aggregates when the monomer radius has a size distribution. Therefore, we plot dashed lines derived with Equation (31) when r₀ = 1.27 μm, ${r}_{0}=2.38\,\mu {\rm{m}}$, and ${r}_{0}=3.49\,\mu {\rm{m}}$. The experimental value is lower than the theoretical lines. From their low value of the tensile strength, Gundlach et al. (2018) inferred that the specific surface energy of ice, ${\gamma }_{\mathrm{ice}}$, has a value of $=0.02\,{\rm{J}}\,{{\rm{m}}}^{-2}$ at low temperatures ($\lesssim 150$ K). However, Gundlach & Blum (2015) also estimated ${\gamma }_{\mathrm{ice}}=0.19\,{\rm{J}}\,{{\rm{m}}}^{-2}$ from the constant sticking threshold velocity of $\sim 10\,{\rm{m}}\,{{\rm{s}}}^{-1}$ for $T\lt 210$ K in their ice impact experiments. The reason for low tensile strength in Gundlach et al. (2018) is thus unknown. We also confirm that the experimental value is higher than the compressive strength theoretically expected by Kataoka et al. (2013b) (Equation (9)).

Zoom In Zoom Out Reset image size

We can estimate the monomer radius of comet 67P using Equation (31) as follows. From the exploration of 67P, it was found that the microporosity of 67P is 0.75–0.85 (Kofman et al. 2015), while the surface porosity is 0.87 (Fornasier et al. 2015). In other words, the volume filling factor of 67P is about 0.13–0.25. The averaged nucleus bulk density of 67P was estimated to be 0.533 ${\rm{g}}\,{\mathrm{cm}}^{-3}$ (Pätzold et al. 2016). In consideration of the high dust-to-water mass ratio of 67P (e.g., Fulle et al. 2016), its main material is not H₂O ice, but silicate. Assuming that the bulk density is 0.533 ${\rm{g}}\,{\mathrm{cm}}^{-3}$ and the material density is 2.65 ${\rm{g}}\,{\mathrm{cm}}^{-3}$ (Table 1), we can estimate that the volume filling factor of 67P is about 0.20, which is consistent with the values of 0.13–0.25. Substituting ${P}_{x,\max }=3-200$ Pa (see Section 1), $\gamma =20\,\mathrm{mJ}\,{{\rm{m}}}^{-2}$, and ${\phi }_{\mathrm{init}}=0.20$ in Equation (33), we find that the monomer radius of 67P has to be 3.3–220 μm. To explain the low tensile strength of 67P using only our model, we have to consider a larger radius than that of interstellar materials.

Another idea to decrease the tensile strength is assuming an aggregate of aggregates (e.g., Blum et al. 2014, 2017) instead of a simple aggregate of monomers as discussed in this paper. The aggregate-of-aggregate model may explain the low strength with small monomers, but this is beyond the scope of this paper.