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

A projectile impact onto a granular target produces an ejecta curtain with the heterogeneous material distribution. Understanding how the heterogeneous pattern forms is potentially important for understanding how crater rays form. Previous studies predicted that the pattern formation is induced by inelastic collisions of ejecta particles in the early stages of crater formation and is terminated by the ejecta's expanding motion. In this study, we test this prediction based on a hyper-velocity impact experiment together with N-body simulations where the trajectories of inelastically colliding granular particles are calculated. Our laboratory experiment suggests that pattern formation is already completed on a timescale comparable to the geometrical expansion of the ejecta curtain, which is ~ 10 microseconds in our experiment. Our simulations confirm the previous prediction that the heterogeneous pattern grows through initial inelastic collisions of particle clusters and subsequent geometric expansion with no further cluster collisions. Furthermore, to better understand the two-stage evolution of the mesh pattern, we construct a simple analytical model that assumes perfect coalescence of particle clusters upon collision. The model shows that the pattern formation is completed on the timescale of the system's expansion independently of the initial conditions. The model also reproduces the final size of the clusters observed in our simulations as a function of the initial conditions. It is known that particles in the target are ejected at lower speeds with increased distance to the impact point. The difference in the ejection speed of the particles may result in the evolution of the mesh pattern into rays.


INTRODUCTION
Mutual collisions between planets, satellites, and minor bodies frequently occur in the solar system. A hypervelocity impact of a small body onto a larger body leads to the production of an excavation flow that eventually opens a crater. The materials in the excavation flow launch from the outer edge of a growing crater and form an inverted cone-like structure called an ejecta curtain.
Impact craters are often accompanied by radial streaks of ejecta called crater rays (e.g. Oberbeck 1971, Melosh 1989, Hawke et al. 2004. Crater rays cause degradation of preexisting craters, which may control crater equilibrium found for small lunar craters (Minton et al. 2019). Understanding how the morphology of crater rays is related to impact conditions may help us better understand how the size frequency distribution of small craters on planetary and lunar surfaces are determined.
While there are several studies characterizing crater rays on a variety of planetary bodies (e.g. Trask & Guest 1975, McEwen et al. 2005, only a few studies have addressed the physical mechanisms responsible for crater ray formation (e.g. Shuvalov 2012, Kadono et al. 2015, 2020, Sabuwala et al. 2018. Shuvalov (2012) propose that interaction of impact-induced shock waves with pre-existing craters can produce crater rays.
More recently, Kadono et al. (2015Kadono et al. ( , 2020 propose an alternative idea, which we focus on in this study, that crater rays are formed by an enhanced contrast in the particle number density in an ejecta curtain. They conducted impact experiments onto granular targets and found that the resulting ejecta curtains have mesh-like patterns. The observed mesh pattern simply expands with the ejecta curtain and suggested that the pattern formation had already been completed in an early stage of the ejecta curtain formation. They also performed numerical calculations and suggested that the mesh patterns likely result from inelastic collisions between gran-ular particles. However, the previous work did not address when exactly the pattern formation occurred in the experiments and, more fundamentally, how the final size of the particle clusters constituting the mesh pattern is related to the initial condition of the ejecta curtain formation. A better constraint on the timing of pattern formation is necessary for further investigating the processes of the pattern formation through future experiments. In this study, we focus on the very early phase of ejecta curtain formation and study in detail how particle clustering proceeds and how it is terminated in an expanding curtain. We present the results of a hyper-velocity impact experiment where we observed the development of a mesh-like pattern in the ejecta curtain on the timescale of the ejecta's expansion (Section 2). We also numerically simulate pattern formation of inelastically colliding granular particles with expanding systematic velocities (Section 3). To interpret the pattern formation observed in the simulations, we construct a simple physically-based model for clustering in an expanding particle system (Section 4). We discuss our model's validity and the evolution of the mesh pattern into rays in Section 5. A conclusion is presented in section 6.

Method
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 band pass filter, which enable us to image an ejecta curtain as early as 10 µs after the hyper-velocity 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 Institute of Space and Astronautical Science of 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 −1 . 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 selfluminous high-temperature plume on the curtain observation was minimized by mounting a band-pass filter corresponding to the laser wavelength in front of the camera lens. The spatial resolution was 110 µm per pixel. 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. Figure 2 shows the snapshots of the ejecta curtain within 100 µs after the impact. The vertical and radial 1 velocities of the ejecta in this early phase are estimated to be in the range of ∼ 100-1000 m s −1 . 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.

Results
We find that the ejecta curtain exhibits a mesh-like pattern already at 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 impacts. 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. In our laboratory experiment, we were unable to observe how the heterogeneous pattern developed during the very early phase of the ejecta curtain formation. To better understand this process, we carry out N-body simulations where we calculate the motion and mutual inelastic collisions of granular particles in an expanding curtain. We particularly focus on the evolution of the particles' velocity dispersion and density pattern. To quantify the properties of the particle pattern, we regard the pattern as a collection of particle clusters and employ a cluster analysis.

Method
Our simulations use the open-source N-body code RE-BOUND (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 twodimensional 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 its 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 0 . 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 by the subscript j) is given an initial velocity consisting of random and expansion components, 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 rand,max to v rand,max , where v 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 where v rel,n and v rel,t are the normal and tangential components of the relative velocity before the collision, respectively, v rel,n and v rel,t are those after the collision, and e is the coefficient of restitution. We ignore the rotation of the particles by assuming that the particles 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 with 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 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 −1 so that the expansion velocity of the system, ΩL 0 /2 = 200 m s −1 , 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 are 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) show 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.
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 rel r/v rand,max ≈ 1.3 µs, where v rel is the relative velocity of the particles. We therefore set ∆t = 0.15 µs.

Data Analysis
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 root-mean-squared (RMS) particle velocity as where the brackets represent an average over all particles (again, each particle is labeled by j). To reduce statistical errors, we also take the average of v RMS over 10 runs. The masses of 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) where a j and b j are the average distances between 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 is, 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 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 and therefore the cluster distribution within it is nonuniform.  Table 1 for the adopted parameter set)

Results
. 6 Nakazawa et al. visible already at t = 15 µs (upper 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. 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 to 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.
To illustrate how we determine the optimal number of clusters from the cluster analysis, we plot in Figure 5 the particleaveraged silhouette coefficient s j as a function of the number of clusters N at different times for a single run. Here, the values of s j are calculated from N = 1 to 30000 for every 100 increment in N. At all times, s j (N) has a single maximum, allowing us to uniquely define the optimal number of clusters as the number of clusters N( s j 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.   Table 2 lists the optimal cluster number N( s j max ) averaged over 10 simulation runs, the corresponding mean cluster size N = N tot /N( s j max ), and the standard deviation δN of the cluster size averaged over the 10 runs. Here, we search for N( s j max ) for individual runs in two steps, first obtaining a rough estimate of s j max by computing s j from N = 1 to 30000 for every 100 increment in N, and second obtaining a more precise value of s j max by computing s j for every 5 increment in N in the vicinity of N( s j max ) from the first step. Figure 7 plots N and δN versus t given in Table 2. One can see that 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. 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).

ANALYTIC MODEL
Our simulations presented in the previous section confirm the expectation by Kadono et al. (2015Kadono et al. ( , 2020) that the mesh pattern in ejecta evolves in two stages, initial growth through inelastic particle collisions and subsequent geometric expansion with no particle collision. However, it is yet to be addressed when the initial growth stage is terminated and what initial conditions set the final size of the clusters constituting the pattern. To address these questions, we here construct an analytic model for the growth of particle clusters in an expanding particle system.

Model Description
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 experience 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 the random velocity component dominates over the non-isotropic 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.

Mass
The assumptions made above gives a simple relationship between the cluster velocity dispersion and cluster mass. Let the velocities of two colliding clusters 1 and 2 be V 1 and V 2 and the velocity of the merged cluster be V . Using the conservation of momentum law and the assumption that the clusters are approximately equal in size, we have The assumption that the isotropic, random motion dominates the clusters' velocities gives 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, σ ≡ (V 2 ) 1/2 , can be written as Equation (8) shows that the velocity dispersion of the clusters decreases by a factor of 1/ √ 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 where N is the average number of particles per cluster (which is equal to the mean cluster mass in units of the monomer mass) and σ 0 and N 0 are the initial values of σ and N, respectively.

Time Evolution of the Cluster Mass
The time evolution of N through cluster mergers can be described by  Table 1. where τ is the mean collision time. For spherical clusters moving on a two-dimensional space, τ can be estimate as where R and n are the average radius and surface number density of the clusters (the factor 4R represents the collisional cross section for equal-sized disks on a plane). Assuming that each two-dimensional cluster is densely packed, R is related to N and the particle radius r as N = πR 2 /(πr 2 ), or Since the number of clusters in the system is ≈ N tot /N, the cluster surface number density n can be estimated as 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 constant expansion rate Ω, with L H 0 being the initial value of L H . Substituting Equations (9) and (11)-(14) into Equation (10), we obtain whereṄ 0 is a constant defined bẏ The solution to eq. (15) is where N 0 is the initial value of N. The set of Equations (9) and (15) predicts σ and N as a function of t under the assumption that the expansion motion of the system has negligible effect on the cluster collision velocity.

"Freezeout" of the Cluster Pattern
For late times where 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 then should "freeze out," i.e., should simply expand self-similarly with time as observed in our experiment and simulation. We expect the freezeout to occur when the collisional timescale τ becomes longer than the expansion timescale ∼ Ω −1 . We thus define the freezeout time t freeze by To account for the freezeout, we stop the evolution of σ and N stops at t = t freeze . The model confirms the expectation from our laboratory expriment that the cluster growth is completed on a timescale comparable to the expansion timescale 1/Ω. By definition, the freezeout time satisfies τ(t = t freeze ) = 1/Ω. Using Equations (12)-(14), we rewrite τ as For N 1, we have Substituting Equations (16) and (20) into Equation (19), we obtain τ ≈ (1 + Ωt) ln(1 + Ωt) Ω .

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 Substituting Equation (23) into Equation (9), we obtain the final velocity dispersion of the clusters as

Comparison with Simulations
Now we test the analytic model derived above against our numerical simulations. Figure 9 compare 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 freezeout of the clusters at t > t freeze . The parameters in the model are set to N 0 = 1, L V = L H 0 = 20 mm, and σ 0 = 32.6 m s −1 . The value of σ 0 is taken to be close to the initial RMS speed of the particles in the simulations. The assumed parameters givė N 0 = 0.49 µs −1 . The definition of the freezeout time (Equation (18)) gives t freeze = 35 µs for these parameter choices.
We find that the analytic model neglecting the freezeout 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 freezeout 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 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 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 freezeout 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, σ 0 = 32.6 m s −1 ) into Equations (23) and (24), we obtainN(t freeze ) = 14.9, σ(t freeze ) = 8.4 m s −1 . 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. 5. DISCUSSION

Validity and Limitations of the Model
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 ramains valid, we reran 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 freezeout. We find that the 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.  Figure 9).
(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 −1 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 equal-sized. Kadono et al. (2019) report 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.

Evolution of the Mesh Pattern to Crater Rays
In this section, we discuss a possible evolution path from the mesh-like pattern of an ejecta curtain to 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); (2) the sediments consisting of 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).
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 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.

SUMMARY AND CONCLUSIONS
We have conducted laboratory and numerical experiments to study pattern formation and evolution in an early phase of ejecta curtain formation. Our key findings are summarized as follows.
1. In our hyper-velocity impact experiments using a 10 mm-sized projectile, the ejecta curtain already exhibited mesh-like patterns at 10 µs after the impact ( Figure  2). This timescale is comparable to the timescale of the ejecta's expansion, suggesting that pattern formation is completed within the expansion timescale.
2. In the numerical simulations, we employed a cluster analysis to quantify cluster formation in an expanding particle system. At early times, we find that the particle random velocity decays and the mean cluster size increases with time, confirming the prediction by Kadono et al. (2015) that the pattern formation (particle clustering) is induced by inelastic particle collisions. At later times, the cluster pattern simply expands with time, consistent with the results of the previous laboratory experiments and ours.
3. We have constructed an analytic model for pattern formation that assumes perfect coalescence of clusters at early times and collisionless geometric expansion of the clusters at later times (Section 4.1). The model confirms the expectation from our laboratory experiment that the initial cluster growth stage is completed on a timescale comparable to the expansion timescale of the particle system (Section 4.1.3). Furthermore, the model reasonably reproduces the final mass of the clusters seen in our simulations as a function of initial conditions (Figure 9). Our model is valid as long as the coefficient of restitution of the particles is 0.8 or smaller ( Figure 10).
If the ejecta particles have a large positive vertical velocity gradient, the mesh pattern forming through inelastic particle collisions would evolve into a vertically elongated pattern. This may eventually form radially extended crater rays upon deposition (Section 5.2 and Figure 11). The analytic model presented in this study will allow us to better understand the pattern formation in this process. However, our current model is limited to targets made of monodisperse particles. Extension of the model to polydisperse granular targets should be done in future work.

Collisional Growth
Mesh-like pattern formation Pattern stretching Crater ray like structure Deposition Figure 11. Illustration of a possible evolution path from the mesh pattern to the observed ray. The mesh-like pattern formed by the inelastic collision is stretched in the radial direction by the velocity contrast according to the distance from the impact point, resulting in a crater ray-like pattern.