Modeling Early Clustering of Impact-induced Ejecta Particles Based on Laboratory and Numerical Experiments

The purpose of our impact experiment is to observe pattern formation in a granular ejecta curtain immediately after an impact. Imaging the early phase of crater formation using a video camera is generally hindered by self-emission from high-temperature vapor plumes. To avoid this issue, we use a monochromatic light source and a bandpass filter, which enable us to image an ejecta curtain as early as 10 μs after the hypervelocity impact of a millimeter-sized projectile. For comparison, Kadono et al. (2020) conducted similar experiments but were only able to image the ejecta curtains ∼10 ms after the impacts.

A schematic diagram of the experimental setup is shown in Figure 1. We used a vertical two-stage light gas gun at the Institute of Space and Astronautical Science of the Japan Aerospace Exploration Agency. A spherical polycarbonate projectile with a diameter of 4.8 mm was accelerated and impacted on the surface. We used zircon beads (FZB-400; the span of sizes is 30–63 μm) as granular targets. The median size of the beads was 50 μm. The impact velocity was set to 6.4 km s⁻¹. The experimental chamber was evacuated down to 5 Pa prior to the shot. To cut off the propellant gas mixing associated with projectile acceleration, we surrounded the targets with a gas shield with only a narrow opening of several centimeters for the projectile. A high-speed video camera (Shimadzu, HPV-X) and a 640 nm laser light source were used to observe the ejecta curtain. The effects of the self-luminous high-temperature plume on the curtain observation was minimized by mounting a bandpass filter corresponding to the laser wavelength in front of the camera lens. The spatial resolution was 110 μm pixel^–1. We obtained a total of 128 frames with a frame rate of 1 μs. The frame rate is close to the characteristic time for the projectile penetration. Thus, we could observe the curtain growth as early as ∼10 μs after the impact.

Zoom In Zoom Out Reset image size

Figure 2 shows the snapshots of the ejecta curtain within 100 μs after the impact. The vertical and radial ⁴ velocities of the ejecta in this early phase are estimated to be in the range of ∼100–1000 m s⁻¹. Because the ejecta immediately after the impact have a radial extent of ∼10 mm, they expand radially on a timescale of ∼10–100 μs. Below, we call the timescale defined in this way the expansion timescale.

Zoom In Zoom Out Reset image size

We find that the ejecta curtain already exhibits a mesh-like pattern 10 μs after the impact (see the upper right panel of Figure 2). Moreover, the snapshots at 76 and 91 μs after the impact show that the pattern simply expands with the curtain (see the red circles in the lower panels of Figure 2). Kadono et al. (2020) also observed the expansion of ejecta curtain patterns in their similar experiments, but at ∼10 ms after impact. Our experimental results presented here indicate that pattern formation occurs and is completed much earlier, on a timescale comparable to the expansion timescale of the ejecta curtain.

Our simulations use the open-source N-body code REBOUND (Rein & Liu 2012). We use IAS15 (Rein & Spiegel 2015) as an integrator. We approximate a portion of an ejecta curtain as a thin rectangular particle sheet and treat the motion of the ejecta particles in this portion in a two-dimensional Cartesian reference frame (x, y) comoving with the particles' center of mass, with the x- and y-directions corresponding to the directions parallel and perpendicular to the circumference of the cone-like curtain, respectively. We neglect the gravity in this comoving frame, assuming that the particles' center of mass follows a parabolic trajectory under the gravity. Because the particles in the laboratory frame move away from the impact point, the particle sheet expands in the y-direction. Including this expansion motion is essential here, as we aim to see how the expansion of the curtain terminates the growth of particle clusters at late times. In our simulations, we mimic the expanding motion by adding a velocity proportional to x to the x component of the initial velocity of each particle (see below). In reality, the particle sheet would also expand in the direction perpendicular to the curtain's circumference if particles ejected earlier have higher ejection speeds (see Section 5.2 for the potential implications for crater-ray formation). For those cases, the x-axis in our simulations should be regarded as the axis along which the particle system expands faster.

The particles are initially placed in a square area of side length L₀. Before we set the initial velocities of the particles, we allow the particles to repel elastically in a periodic box to remove initial particle overlap. Each particle (labeled with the subscript j) is given an initial velocity consisting of random and expansion components,

$$\begin{eqnarray}&&{({v}_{{jx}},{v}_{{jy}})}_{t=0}=({v}_{\mathrm{rand},{jx}}+{v}_{\exp ,j},{v}_{\mathrm{rand},{jy}}),\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{v}_{\exp ,j}={\rm{\Omega }}{x}_{j},\end{eqnarray} \tag{ 2 }$$

where t is time, v_rand,jx and v_rand,jy are random numbers, ${v}_{\exp ,j}$ is the initial expansion velocity, and Ω is the rate of expansion. With this initial condition, the particles expand in the x-direction and also mutually collide nearly isotropically. The random velocity components are sampled from a uniform distribution ranging from $-{v}_{\mathrm{rand},\max }$ to ${v}_{\mathrm{rand},\max }$, where ${v}_{\mathrm{rand},\max }$ is a constant. Whenever two particles collide, we decompose their velocities into the center-of-mass and relative velocity components and change the latter as

$$\begin{eqnarray}&&v{{\prime} }_{\mathrm{rel},{\rm{n}}}=-{{ev}}_{\mathrm{rel},{\rm{n}}},v{{\prime} }_{\mathrm{rel},{\rm{t}}}={v}_{\mathrm{rel},{\rm{t}}},\end{eqnarray} \tag{ 3 }$$

where v_rel,n and v_rel,t are the normal and tangential components of the relative velocity before the collision, respectively; $v{{\prime} }_{\mathrm{rel},{\rm{n}}}$ and $v{{\prime} }_{\mathrm{rel},{\rm{t}}}$ are those after the collision; and e is the coefficient of restitution. We ignore the rotation of the particles by assuming that they exert no torque on each other upon collision. In our simulations, we set e as a constant. We carry out 10 independent simulation runs with 10 different sets of initial particle random velocities but the same initial particle configuration.

The values of the parameters adopted in the simulations are listed in Table 1. All particles in the simulations are assumed to be spherical and equal-sized, with the radius r being equal to those of the zircon beads used in the experiment presented in Section 2. Because we exert no external force on individual particles, our simulation results are independent of the particle mass. For this reason, we adopt the code units in which the particle mass is taken to be unity. The expansion rate Ω is set to 0.02 μs⁻¹ so that the expansion velocity of the system, ΩL₀/2 = 200 m s⁻¹, is comparable to the horizontal velocity of the ejecta curtain pattern observed in our laboratory experiment. The maximum random velocity is taken to be 20% of the system expansion velocity. The coefficient of restitution of our zircon particles is unknown, so we arbitrarily take the default value of e to be 0.6 so that the pattern formed is visible enough to analyze. The validity of our choice of e is discussed in Section 5. Kadono et al. (2015) showed that particles of a lower restitution coefficient produce a clearer particle pattern through inelastic collisions. For completeness, we also present simulations for e = 0.1, 0.9, and 1.0.

Table 1. Parameters Adopted in the Simulations

Parameter	Value
Total number of particles Ntot	30,000
Particle radius r	50 μm
Expansion rate Ω	0.02 μs−1
Maximum random speed ${v}_{\mathrm{rand},\max }$	40 m s−1
Restitution coefficient e	0.6
Initial box size L0	20 mm

Download table as:  ASCIITypeset image

The time step Δt must be chosen to be smaller than the crossing times of the particles so that every particle collision can be detected. The particle crossing time is estimated as $2r/{v}_{\mathrm{rel}}\gtrsim r/{v}_{\mathrm{rand},\max }\approx 1.3\,\mu {\rm{s}}$, where v_rel is the relative velocity of the particles. We therefore set Δt = 0.15 μs.

Particles lose their random kinetic energies through inelastic collisions. Since the particles also form clusters through the collisions, the random energy of the system and the size of the clusters should be correlated with each other. Below, we describe how we quantify the two quantities.

As a diagnostic for the random motion of the particles, we define the rms particle velocity as

$$\begin{eqnarray}&&{v}_{{\rm{r}}{\rm{m}}{\rm{s}}}= \langle {({v}_{jx}-{v}_{\exp ,j})}^{2}+{v}_{jy}^{2}{ \rangle }^{1/2},\end{eqnarray} \tag{ 4 }$$

where the brackets represent an average over all particles (again, each particle is labeled with j). To reduce statistical errors, we also take the average of v_rms over 10 runs. The masses of the particle clusters are measured using a cluster analysis. In this analysis, all particles in a given snapshot of simulations are grouped into clusters. We adopt hierarchical clustering in which pairs of clusters (initially individual particles) are combined successively. Specifically, we use Ward's method (Joe & Ward 1963), a commonly used hierarchical clustering algorithm that combines clusters so that the variance of the particle positions within the clusters is minimized. However, this algorithm alone does not tell us the optimal (i.e., most likely) number of clusters. To determine it, we introduce a quantity called the silhouette coefficient. For each particle j, its silhouette coefficient s_j is defined as (Rousseeuw 1987)

$$\begin{eqnarray}&&{s}_{j}=\displaystyle \frac{{b}_{j}-{a}_{j}}{\max \{{a}_{j},{b}_{j}\}},\end{eqnarray} \tag{ 5 }$$

where a_j and b_j are the average distances between the particle j and the other particles inside and outside the the cluster, respectively, to which it belongs. The higher the average silhouette coefficient 〈s_j 〉, the better the clusters are separated. For each snapshot of particle distribution, we pick up the set of clusters that give the highest 〈s_j 〉 and calculate the mean and standard deviation of their masses. We also take the averages of the means and standard deviations over 10 simulation runs. The cluster analysis we employ is more naturally suited to our problem than Fourier analysis, as used by Kadono et al. (2015), because our particle system has an open, expanding boundary; therefore, the cluster distribution within it is nonuniform.

Figure 3 shows snapshots of the particle distribution obtained from a single simulation run. Particle clustering is already visible at t = 15 μs (top right panel). At t > 30 μs, the cluster pattern simply expands in the x-direction rather than growing further (middle right, bottom left, and bottom right panels), consistent with the result of our laboratory experiment. The pattern is mesh-like, demonstrating that inelastic collisions between particles at early times induce the formation of the mesh pattern in ejecta curtains observed in laboratory experiments, including ours (see Section 2). Because of the unidirectional expansion, the particle clusters constituting the mesh pattern are elongated along the x-axis. If the particles had higher expansion velocities in the y-direction, the resulting mesh pattern would be more elongated along the y-axis.

Zoom In Zoom Out Reset image size

Figure 4 shows the rms random velocity v_rms (Equation (4)) averaged over the 10 runs as a function of t. We find that v_rms decreases at t < 20 μs when the mesh pattern is growing and approaches a constant at later times when the pattern is expanding. This serves as another piece of evidence for particles' inelastic collisions being the cause of pattern formation.

Zoom In Zoom Out Reset image size

To illustrate how we determine the optimal number of clusters from the cluster analysis, we plot in Figure 5 the particle-averaged silhouette coefficient 〈s_j 〉 as a function of the number of clusters ${ \mathcal N }$ at different times for a single run. Here the values of 〈s_j 〉 are calculated from ${ \mathcal N }$ = 1 to 30,000 for every 100 increments in ${ \mathcal N }$. At all times, $\langle {s}_{j}\rangle ({ \mathcal N })$ has a single maximum, allowing us to uniquely define the optimal number of clusters as the number of clusters ${ \mathcal N }(\langle {s}_{j}{\rangle }_{\max })$ that gives the maximum silhouette coefficient. Figure 6 illustrates how our cluster analysis works for the particular example shown in Section 3 at t = 0 and 60 μ s.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2 lists the optimal cluster number ${ \mathcal N }(\langle {s}_{j}{\rangle }_{\max })$ averaged over 10 simulation runs, the corresponding mean cluster size $\overline{N}={N}_{\mathrm{tot}}/{ \mathcal N }(\langle {s}_{j}{\rangle }_{\max })$, and the standard deviation δ N of the cluster size averaged over the 10 runs. Here we search for ${ \mathcal N }(\langle {s}_{j}{\rangle }_{\max })$ for individual runs in two steps, first obtaining a rough estimate of $\langle {s}_{j}{\rangle }_{\max }$ by computing 〈s_j 〉 from ${ \mathcal N }$ = 1 to 30,000 for every 100 increment in ${ \mathcal N }$ and then obtaining a more precise value of $\langle {s}_{j}{\rangle }_{\max }$ by computing 〈s_j 〉 for every five increments in ${ \mathcal N }$ in the vicinity of ${ \mathcal N }(\langle {s}_{j}{\rangle }_{\max })$ from the first step. Figure 7 plots $\overline{N}$ and δ N versus t given in Table 2. One can see that $\overline{N}$ increases in the first 30 μ s, when the rms particle random velocity decreases (see Figure 4), indicating that our cluster analysis well captures pattern formation through inelastic collisions. At later times, the cluster mass plateaus out, reflecting the collisionless expansion of the system. Close inspection shows that the derived cluster mass gradually decreases after t = 60 μs. We interpret this as suggesting that our cluster analysis does not fully resolve the clusters at t ≲ 60 μ s.

Zoom In Zoom Out Reset image size

Table 2. Properties of the Optimal Clusters at Different Times t (Averaged over 10 Runs)

t (μ s)	${ \mathcal N }(\langle {s}_{j}{\rangle }_{\max })$	$\overline{N}$	δ N
0	10,285	2.9	1.0
30	1265	23.7	8.7
60	1195	25.1	10.8
90	1575	19.0	8.9
120	1800	16.7	8.3

Download table as:  ASCIITypeset image

Figure 8 compares the particle distribution at t = 60 μ s from simulations of different values of e. We see that particle clusters become more appreciable as e decreases, again supporting the idea that the pattern formation is driven by inelastic collisions. The observed trend is also qualitatively consistent with the simulation results by Kadono et al. (2015, their Figure 6).

Zoom In Zoom Out Reset image size

As discussed in the previous section, it is clear that inelastic particle collisions are the cause of the cluster formation seen in our simulations. Because particles constituting clusters inelastically collide with each other frequently, we can expect that colliding clusters lose their initial collision energy quickly. Based on this consideration, our analytic model assumes that clusters coalesce perfectly whenever they collide. As discussed in Section 5, this assumption is valid if the coefficient of restitution of the particles is sufficiently low.

For the moment, we also assume that the random velocity component dominates over the nonisotropic velocity component arising from the system's expansion motion. This assumption breaks down at late times where the expansion motion terminates the collisional evolution of the system. We account for this effect at the end of this subsection. To make our model as simple as possible, we approximate all clusters as equal-sized spheres and only consider their one-dimensional motion.

4.1.1. Relationship between Cluster Velocity Dispersion and Mass

The assumptions made above give a simple relationship between the cluster velocity dispersion and mass. Let the velocities of two colliding clusters 1 and 2 be V₁ and V₂ and the velocity of the merged cluster be $V^{\prime}$. Using the conservation of momentum law and the assumption that the clusters are approximately equal in size, we have

$$\begin{eqnarray}&&V^{\prime} =\displaystyle \frac{1}{2}({V}_{1}+{V}_{2}).\end{eqnarray} \tag{ 6 }$$

The assumption that the isotropic random motion dominates the clusters' velocities gives

$$\begin{eqnarray}&&\left\{\begin{array}{l}\overline{{V}_{1}}=\overline{{V}_{2}}=0,\\ \overline{{V}_{1}^{2}}=\overline{{V}_{2}^{2}}={\sigma }^{2},\\ \overline{{V}_{1}{V}_{2}}=0,\end{array}\right.\end{eqnarray} \tag{ 7 }$$

where σ is the velocity dispersion of the clusters, and the overlines denote ensemble averages over all clusters. Using Equations (6) and (7), the velocity dispersion of the clusters after a single collision, $\sigma ^{\prime} \equiv {(\overline{V{{\prime} }^{2}})}^{1/2}$, can be written as

$$\begin{eqnarray}&&\sigma ^{\prime} =\displaystyle \frac{1}{\sqrt{2}}\sigma \,.\end{eqnarray} \tag{ 8 }$$

Equation (8) shows that the velocity dispersion of the clusters decreases by a factor of $1/\sqrt{2}$ after every merging collision. Because the cluster mass increases by a factor of 2 after every collision, one can relate σ to the cluster mass as

$$\begin{eqnarray}&&\sigma (\overline{N})={\sigma }_{0}{\left(\displaystyle \frac{{N}_{0}}{\overline{N}}\right)}^{1/2},\end{eqnarray} \tag{ 9 }$$

where $\overline{N}$ is the average number of particles per cluster (which is equal to the mean cluster mass in units of the monomer mass), and σ₀ and N₀ are the initial values of σ and $\overline{N}$, respectively.

4.1.2. Time Evolution of the Cluster Mass

The time evolution of $\overline{N}$ through cluster mergers can be described by

$$\begin{eqnarray}&&\displaystyle \frac{d\overline{N}}{{dt}}=\displaystyle \frac{\overline{N}}{\tau },\end{eqnarray} \tag{ 10 }$$

where τ is the mean collision time. For spherical clusters moving in a two-dimensional space, τ can be estimated as

$$\begin{eqnarray}&&\tau =\displaystyle \frac{1}{4\overline{R}n\sigma },\end{eqnarray} \tag{ 11 }$$

where $\overline{R}$ and n are the average radius and surface number density of the clusters (the factor $4\overline{R}$ represents the collisional cross section for equal-sized disks on a plane). Assuming that each two-dimensional cluster is densely packed, $\overline{R}$ is related to $\overline{N}$ and the particle radius r as $\overline{N}=\pi {\overline{R}}^{2}/(\pi {r}^{2})$, or

$$\begin{eqnarray}&&\overline{R}=r{\overline{N}}^{1/2}.\end{eqnarray} \tag{ 12 }$$

Since the number of clusters in the system is $\approx {N}_{\mathrm{tot}}/\overline{N}$, the cluster surface number density n can be estimated as

$$\begin{eqnarray}&&n=\displaystyle \frac{{N}_{\mathrm{tot}}}{{L}_{V}{L}_{H}\overline{N}},\end{eqnarray} \tag{ 13 }$$

where L_H and L_V are the horizontal and vertical length of the system. We assume that the vertical system size is approximately constant, whereas the horizontal system size increases at a constant expansion rate Ω,

$$\begin{eqnarray}&&{L}_{H}={L}_{{H}_{0}}(1+{\rm{\Omega }}t),\end{eqnarray} \tag{ 14 }$$

with ${L}_{{H}_{0}}$ being the initial value of L_H . Substituting Equations (9) and (11)–(14) into Equation (10), we obtain

$$\begin{eqnarray}&&\displaystyle \frac{d\overline{N}}{{dt}}=\displaystyle \frac{{\dot{N}}_{0}}{1+{\rm{\Omega }}t},\end{eqnarray} \tag{ 15 }$$

where ${\dot{N}}_{0}$ is a constant defined by

$$\begin{eqnarray}&&{\dot{N}}_{0}=\displaystyle \frac{4{{rN}}_{0}^{1/2}{N}_{\mathrm{tot}}{\sigma }_{0}}{{L}_{V}{L}_{{H}_{0}}}\,.\end{eqnarray} \tag{ 16 }$$

The solution to Equation (15) is

$$\begin{eqnarray}&&\overline{N}(t)={N}_{0}+\displaystyle \frac{{\dot{N}}_{0}}{{\rm{\Omega }}}\mathrm{ln}(1+{\rm{\Omega }}t),\end{eqnarray} \tag{ 17 }$$

where N₀ is the initial value of $\overline{N}$. The set of Equations (9) and (15) predicts σ and $\overline{N}$ as a function of t under the assumption that the expansion motion of the system has a negligible effect on the cluster collision velocity.

4.1.3. "Freeze-out" of the Cluster Pattern

For late times when the expansion motion of the system dominates over the random motion of individual clusters, we expect that the clusters can no longer grow collisionally. The cluster pattern should then "freeze out," i.e., simply expand self-similarly with time as observed in our experiment and simulation. We expect the freeze-out to occur when the collisional timescale τ becomes longer than the expansion timescale ∼Ω⁻¹. We thus define the freeze-out time t_freeze by

$$\begin{eqnarray}&&\tau (t\gt {t}_{\mathrm{freeze}})\gt \displaystyle \frac{1}{{\rm{\Omega }}}\,.\end{eqnarray} \tag{ 18 }$$

To account for the freeze-out, we stop the evolution of σ, and $\overline{N}$ stops at t = t_freeze.

The model confirms the expectation from our laboratory experiment that the cluster growth is completed on a timescale comparable to the expansion timescale 1/Ω. By definition, the freeze-out time satisfies τ(t = t_freeze) = 1/Ω. Using Equations (12)–(14), we rewrite τ as

$$\begin{eqnarray}&&\tau =\displaystyle \frac{{L}_{V}{L}_{{H}_{0}}(1+{\rm{\Omega }}t)}{4{{rN}}_{\mathrm{tot}}{\sigma }_{0}{N}_{0}^{1/2}}\overline{N}.\end{eqnarray} \tag{ 19 }$$

For $\overline{N}\gg 1$, we have

$$\begin{eqnarray}&&\overline{N}\approx \displaystyle \frac{{\dot{N}}_{0}}{{\rm{\Omega }}}\mathrm{ln}(1+{\rm{\Omega }}t).\end{eqnarray} \tag{ 20 }$$

Substituting Equations (16) and (20) into Equation (19), we obtain

$$\begin{eqnarray}&&\tau \approx \displaystyle \frac{(1+{\rm{\Omega }}t)\mathrm{ln}(1+{\rm{\Omega }}t)}{{\rm{\Omega }}}.\end{eqnarray} \tag{ 21 }$$

Therefore, the condition τ(t = t_freeze) = 1/Ω implies

$$\begin{eqnarray}&&(1+{\rm{\Omega }}{t}_{\mathrm{freeze}})\mathrm{ln}(1+{\rm{\Omega }}{t}_{\mathrm{freeze}})\approx 1,\end{eqnarray} \tag{ 22 }$$

from which we find t_freeze ≈ 0.76/Ω. We note that the freeze-out time depends only on Ω.

4.1.4. Final Cluster Mass and Velocity Dispersion as a Function of Initial Conditions

From Equation (17) and t_freeze ≈ 0.76/Ω, we can write the final cluster size as a function of the initial cluster size and random velocity dispersion as

$$\begin{eqnarray}&&\overline{N}({t}_{\mathrm{freeze}})\approx {N}_{0}+0.57\displaystyle \frac{{\dot{N}}_{0}}{{\rm{\Omega }}}.\end{eqnarray} \tag{ 23 }$$

Substituting Equation (23) into Equation (9), we obtain the final velocity dispersion of the clusters as

$$\begin{eqnarray}&&\sigma ({t}_{\mathrm{freeze}})\approx {\sigma }_{0}{\left(\displaystyle \frac{{N}_{0}}{{N}_{0}+0.57{\dot{N}}_{0}/{\rm{\Omega }}}\right)}^{1/2}.\end{eqnarray} \tag{ 24 }$$

Now we test the analytic model derived above against our numerical simulations. Figure 9 compares the cluster random velocity and mean cluster size measured in the simulations with those predicted from the analytic model (Equations (9) and (15), respectively) with and without the freeze-out of the clusters at t > t_freeze. The parameters in the model are set to N₀ = 1, ${L}_{V}={L}_{{H}_{0}}=20\,\mathrm{mm}$, and σ₀ = 32.6 m s⁻¹. The value of σ₀ is taken to be close to the initial rms speed of the particles in the simulations. The assumed parameters give ${\dot{N}}_{0}=0.49\,\mu {{\rm{s}}}^{-1}$. The definition of the freeze-out time (Equation (18)) gives t_freeze = 35 μs for these parameter choices.

Zoom In Zoom Out Reset image size

We find that the analytic model neglecting the freeze-out well reproduces the evolution of v_rms, at least at early times t ≲ 30 μ s. This indicates that energy dissipation of particle clusters through perfect coalescence explains the decrease of the random velocity. At t ≳ 30 μ s, the random velocity from the simulation decreases more slowly than predicted from the model, suggesting that this model overestimates the frequency of cluster collisions at late times.

The analytic model neglecting the freeze-out reproduces the mean cluster size in the simulations to within a factor of 2. The agreement is less good than for the random velocity, in particular at t ≤ 60 μs, where $\overline{N}$ from the simulation is considerably higher than predicted by the model. This is likely because the clusters are so close that our cluster analysis using the silhouette coefficient s_j (Equation (5)) does not fully resolve them. Density-based cluster analysis algorithms could mitigate this issue but would introduce additional free parameters (for example, the DBSCAN algorithm (Ester et al. 1996) requires the detection radius and the minimum number of cluster members). We defer further investigation of this issue to future work. At t ≳ 80 μ s, the model predicts higher $\overline{N}$ than observed in the simulation. This is another indication that cluster collisions at late times are less frequent than expected by the model.

Accounting for the freeze-out of the clusters at t > t_freeze significantly improves the model predictions at late times (see the dotted lines in Figure 9). Substituting the initial conditions of the simulation (Table 1, σ₀ = 32.6 m s⁻¹) into Equations (23) and (24), we obtain $\bar{N}({t}_{\mathrm{freeze}})=14.9$, σ(t_freeze) = 8.4 m s⁻¹. These model predictions are consistent with the final values from the simulation.

To summarize this section, we have shown that the initial growth of particle clusters in an expanding particle system terminates when the collisional timescale exceeds the expansion timescale 1/Ω. This occurs at t_freeze ∼ 1/Ω after the clusters start growing. Our analytic model quantitatively reproduces the evolution of the cluster mass and velocity dispersion as a function of the initial conditions.

In this section, we discuss the validity and limitations of our analytic model presented in Section 4.

The analytic model assumes perfect coalescence of the clusters. It is natural to expect that the model would not be applicable to granular materials with a high coefficient of restitution e. To quantify the range of e where the model remains valid, we reran the simulations presented in Section 3 but adopting different values of e. Figure 10 shows the time evolution of the random velocity obtained from these simulations compared with the prediction from the model including the cluster freeze-out. We find that the simulation results for e < 0.8 deviate from the model prediction, suggesting that the assumption of perfect coalescence is only valid for e ≲ 0.8.

Zoom In Zoom Out Reset image size

The next question is whether the values of e of real granular materials fall into the range where our model remains valid. Measurements show that e = 0.97 for 3 mm glass beads and e = 0.87 for 6 mm cellulose (Foerster et al. 1994), indicating that perfect cluster coalescence may not be a good approximation for these millimeter-sized particles. However, for smaller particles, cohesive forces can lead to a smaller coefficient of restitution (Royer et al. 2009). There is also experimental evidence that collision velocities exceeding 10 m s⁻¹ cause a reduction of e to below 0.9 (Labous et al. 1997). Therefore, we expect that the assumption of perfect cluster coalescence is valid for a high-velocity impact onto a target made of submillimeter-sized granules.

Another strong assumption in our model is that all ejecta particles are of equal size. Kadono et al. (2019) reported that a target made of different-sized granules produces ejecta curtains with filament patterns. Extending our analytic model to polydisperse granular targets will be the subject of future work.

In this section, we discuss a possible evolution path from the mesh-like pattern of an ejecta curtain to the crater rays frequently observed around flesh craters.

Crater rays on planetary bodies have two important features: (1) the rays' azimuthal pattern has a characteristic angle scale (Kadono et al. 2015), and (2) the sediments that construct the rays are stretched in the radial direction from the host craters' center (e.g., McEwen et al. 2005; Kadono et al. 2015). As already discussed by Kadono et al. (2015) and also demonstrated in our experiment and simulations, mesh pattern formation through inelastic particle collisions followed by the pattern's geometric expansion naturally explains the first feature of the crater rays.

The second feature may also be explained if ejecta particles launched earlier have higher launch speeds. As we mentioned in Section 3.3, if the expansion associated with this velocity difference dominates over the expansion of the curtain's circumference, the particle mesh pattern should be stretched vertically. Indeed, it is known that the launch speed of the ejecta particles decreases with increasing distance from the impact point (e.g., Housen et al. 1983). Therefore, we can envision the scenario in which the mesh pattern forming in an ejecta curtain through initial particle clustering evolves into a vertically elongated, crater ray–like pattern by the vertical expanding motion as schematically illustrated in Figure 11. In our experiment, the pattern evolution envisioned above was not observed due to the limited time window of the high-speed imaging. On the other hand, such pattern evolution was observed in the previous impact experiments conducted by Kadono et al. (2015).

Zoom In Zoom Out Reset image size

In the case of natural impacts occurring under low-gravity environments, such as the lunar surface, the mesh pattern is more likely to evolve into rays than in laboratory environments under the Earth's gravity due to longer flight times. Nevertheless, not all observed craters have rays around them. We expect that this is because of degradation. On the Moon, the lifetime of the rays is estimated to be 1 Gyr, and most of the formation ages of lunar craters are older than 1 Ga (McEwen et al. 1997). Therefore, rays are often not observed in actual craters.