The interactive dynamics of a binary asteroid and ejecta after medium kinetic impact

Kinetic impacts on a binary asteroid will produce ejecta debris as well as causing a disturbance to the spin-orbit motion of the binary components. Understanding the complex interactions between the ejecta and the binary system can be crucial to both theoretical study and mission design. In this paper, we develop a detailed model for the dynamical evolution of ejecta within a binary asteroid system, in which we combined various perturbations on the ejecta and considered the collisional process with the asteroid. Taking an example scenario like the DART impact, we calculated the fate of ejecta from size 0.01 mm to 10 m and analyzed the distribution of their final deposition locations. The results show that ejecta with diameter less than 0.1 mm escape quickly under the influence of solar radiation pressure, while ejecta of larger size survive longer, i.e., a few cycles in the binary system for more than 400 days. We analyzed the influence of the coefficient of restitution between ejecta and asteroid, and found more damped collisions would allow the ejecta to survive longer in orbit. We also studied the effect of the tumbling rotation of the secondary on the distribution of ejecta on the surface of the secondary, the results suggest that a tumbling rotational state causes obvious changes in Dimorphos’ surface slope, which, however, are not sufficient for triggering surface landslides. In addition, the influence of the internal structure of the primary on the evolution of the ejecta is studied, and we show it has little effect on the evolution of the ejecta but does affect the distribution of ejecta on the surface of the primary. This is due to the change of the distribution of relative acceleration on the surface of the primary, which changes the regions where the ejecta can remain.


Introduction
Binary asteroids account for approximately 16% of near-Earth asteroids (Margot et al. 2002) and are one of the popular targets for deep space exploration. Except for single asteroids that are confirmed to be active in the Solar system, such as (6478) Gault, (101955) Bennu, (136108) Haumea, etc, NASA's Double Asteroid Redirection Test (DART) mission, launched on November 24, 2021, plans to impact the secondary (Dimorphos) of binary asteroid system 65803 Didymos in late 26 September 2022 using a 563-kg impactor at 6.14 km/s, which is also expected to produce an artificial ejection event in this two-body system as well as to change their orbit-spin states Michel et al. 2016;Cheng et al. 2018;Rainey et al. 2019Rainey et al. , 2020Rivkin et al. 2021;Agrusa et al. 2022). Moreover, ESA's Hera mission will arrive in the Didymos system in 4 years later to observe the results of the impact (Agrusa et al. 2022). The evolution of the ejecta bears both theoretical and practical significances to the mission, which can help to understand the general laws of the geological evolution of binary asteroids and reveal the special mechanical behaviour of interplanetary materials in a complex environment. In addition to observation-based study, a theoretical analysis based on numerical models are also used to track and analyze the evo-  System configuration and 4 coordinate frames used for dynamics: H is Heliocentric ecliptic J2000; T is orbit translational frame; A is primary body-fixed frame; B is secondary body-fixed frame lution of ejecta near asteroids which helps to assess to understand the ejecta environment posterior to DART impact and provides safety assurance for the Hera mission.
There are still some intriguing questions worth discussing on the behaviours of ejecta. First, the ejecta that stay near orbit for a long time may pose a threat to future spacecraft observations. The size and velocity of the ejecta may significantly affect the evolution of ejecta. Studying whether there will be ejecta of certain sizes at certain velocities that can survive in orbit for a long time after impact is critical for predicting whether the spacecraft will be in danger. Second, Agrusa et al. (2021) found that Dimorphos is prone to unstable spins in its long-axis direction after DART's impact, which may change the dynamical environment of Dimorphos' surface. Third, the internal structure of the asteroid is not well understood, its possible heterogeneous internal structure may affect the distribution of the ejecta, which, in turn, allows for an indirect way to infer the internal mass distribution of the primary. In order to discuss these questions, we established a complete dynamic model, including considering the various forces on ejecta, using a gravitational field model that can describe asteroids with non-uniform internal mass distribution, establishing a fast collision detection method between ejecta and asteroid, and accounting for the uncertainty of rebound velocity of ejecta that re-impact the asteroid surface.
In this paper, we focus on model building to study the evolution of ejecta in the binary system under the complex dynamical perturbations of the Solar system and calculate its final distribution. The paper is organized as follows. In Sect. 2, we introduce the interactive dynamics modelling, which includes dynamic models of binary asteroids and ejecta, ejecting model and impact reflection model. In Sect. 3, we present the simulation schemes, results, and discussions, including the evolution of ejecta with different sizes and different coefficients of restitution of the surface of the asteroid, the influence of the tumbling rotational state of the secondary on its surface dynamic environment and the influence of the internal structure of the primary on the distribution of ejecta. Our conclusion is given in Sect. 4. These simulations are performed on a large-scale computing cluster using CUDA and OpenMP parallelism techniques.

Interactive dynamics modeling
In this paper, we take the example of the binary asteroid system Didymos to introduce the modelling of interaction dynamics. Its main purpose is the description and presentation of our model which includes five tasks, including (1) Build a dynamic model of the motion of a binary system, which is only affected by mutual gravity and the solar, since the ejecta has almost no effect on the motion of the binary system.
(2) Build a dynamic model of the motion of the ejecta, which is affected by multiple forces of the Solar System.
(3) Calculate the size and velocity of the ejecta formed after the impact event. (4) Build a secondary collision ejection model between ejecta and asteroid. (5) Set the stopping conditions of the simulation. The configuration of the system and coordinate is shown in Fig. 1. Four coordinate systems are adopted: Heliocentric ecliptic J2000 H; orbit translational frame T ; primary body-fixed frame A; secondary body-fixed frame B.

Binary dynamics modeling
The shape model of binary asteroid constructed by Benner et al. (2010) is adopted to describe primary (Didymos), and the shape model constructed by Michel et al. (2016) is adopted to describe secondary (Dimorphos). Table 1 gives some physical properties parameters of the Didymos system, and for more detailed data please refer to Naidu et al. (2020). We set the impact direction to be along the tangential direction of the secondary's orbit, with the impact location at the centre of the secondary. Here we assume a kinetic energy transfer efficiency of 1 (see Rivkin et al. 2021 for specific description), which corresponds to a velocity change of 4 mm/s to the asteroid. For the description of the gravitational field of small bodies, as we are studying the effects of the internal structure of asteroid, we need to choose models that are convenient to change the internal structure of asteroid in terms of gravitational field modelling. Therefore, we use the finite element method for calculating the full 2-body problem proposed by Yu et al. (2019).
The essence of this method is to divide two rigid bodies into multiple elements, derive the expression for the gravitational potential between every two elements, and then perform the summation. Assuming that both asteroids are rigid bodies and there are N 1 nodes in the primary and N 2 nodes in the secondary, for every two nodes α and β (α is one node of the primary, β is one node of the secondary), we need to calculate the gravitational attraction force f αβ on the nodes α, torque t α on the primary and torque t β on the secondary. Then the summation is performed and we get the gravitational attraction force F P S on the primary, torque T P on the primary and torque T S on the secondary. The specific formula can be found in Yu et al. (2019). The finite element method is used to generate two new rigid body models (by the shape models) constructed from tetrahedra, and we can reconstruct precisely the binary components' interior structures by defining the density of the element of model. The model has no restrictions on the shape and internal mass distribution of an asteroid and has acceptable accuracy in the mutual gravitational potential description between two gravitational bodies, therefore avoiding the problem of non-convergence caused by traditional series. Here f αβ , t α and t β are calculated N 1 × N 2 times respectively. We use GPU parallel technology to accelerate it and the relevant C++ code can be found in Gao et al. (2022). In addition, we added the solar tidal forces F ,P and F ,S acting on primary and secondary in the dynamic model, which are calculated by the following equations where R is the vector from the Sun to the binary mass centre, M is the mass of solar, M P and M S is the mass of the primary and the secondary, respectively. R P and R S indicate vectors from the system mass centre to the primary and secondary, respectively. Then we use the Newton-Euler equations to calculate the position and attitude of the two asteroids at each moment.

Ejecta ballistic modeling
Ejecta from Dimorphos are influenced by multiple forces, including the gravitational attraction of the binary system, solar gravity, solar radiation pressure, Poynting-Robertson drag, solar wind drag, and gravitational perturbations from other planets. For the calculation of the gravitational force of a binary system on ejecta, we simplify the gravitational calculation method for a full two-body system as mentioned in Sect. 2.1 to calculate the force of a rigid body on a point mass, as shown in the following formula where F P E and F SE represent the gravitational forces of primary and secondary on ejecta, respectively, r P α is the position of node α in the primary in the world coordinate system T , m i and r i is the mass and position in the frame T of ejecta i, w P α is the nodal weights of node α in the primary, σ ρ P α is the density at node α in the primary, replacing the angle labels with S and β corresponds to the secondary.
Since the DART mission impact time has mainly been determined, we can easily calculate the tidal forces due to the Sun and other planets. We use the J2000 coordinate system to calculate the planetary tides of the major planets in the solar system on the Diymos system form Sep. 1, 2022 to Sep. 1, 2032. This is done as follows. First look up the orbital elements of the major planets in the solar system at the moment of the impact. Then update the positions of all the planets and use the gravity formula to calculate the tides of these planets to the Didymos system in the following decade. The orbital elements can be found in: https://ssd.jpl.nasa.gov/. Figure 2 shows the ratio of planetary tides (magnitude) to solar tidal forces in 10 years. The result shows planetary tidal perturbations are much smaller than solar tidal perturbations in 10 years by at least two orders of magnitude. The planetary tidal disturbance is greatest when the Didymos system is close to Earth (October 4th, 2022), and its ratio to the solar tide is 0.005746. Therefore, our calculations do not take into account the effects of planetary tidal perturbations on ejecta.
Regarding the solar radiation pressure, we introduce β to denote the ratio of the solar radiation pressure to the solar gravitational force (Burns et al. 1979 where is the reflection coefficient, we take 1.15 here. S = 1360 W/m 2 , d and ρ is diameter and density of ejecta. The effect of the solar radiation pressure is highly dependent on the area-to-mass ratio, and in our simulations we focus on the size interval of the ejecta (considered as the homogeneous spherical particle) from 10 −5 m to 10 m, so the value of this ratio is from 7 × 10 −5 to 70. For Poynting-Robertson drag and solar wind drag, the results in Yu et al. (2017) show they have little effect on ejecta, which are ignored in our works. Through the above, the dynamic equation of ejecta i in the coordinate system T is: where F ,T is the tidal force and F ,R is the radiation pressure of the solar on ejecta i in T .

Impact ejection modeling
Here an improved scheme proposed by Housen and Holsapple (2011), Holsapple and Housen (2012), Yu and Michel (2018) is applied to the discretization of the ejecta (for the specific formula, please refer to the appendix in Yu and Michel (2018)). Assuming a medium kinetic impact, for the DART mission, we got the results shown in Table 2. The radius of the impact crater is expected to be 11.412 m, the total mass of the ejecta is 9.567 × 10 5 kg, and the velocity distribution of ejecta is from 0 m/s to 247.387 m/s. Here we Mass of ejecta < 1.2 × v esp 1.398 × 10 5 kg Number of ejecta < 1.2 × v esp 9.875 × 10 8 assume ejecta sizes ranging between 1 mm -10 m and the total number of ejecta generated by the model is 6.759 × 10 9 . The escape velocity v esc at the point of impact is 0.245 m/s, which is lower than the ejecting velocity of the majority of ejecta, meaning that the majority of ejecta ejected will escape. For the simulation, we focus on ejecta with an initial ejection velocity of less than 1.2 × v esc , which corresponds to ejecta ejected between 10.827 m and 11.412 m from the centre of the crater. The total number of this part of the ejecta is 9.875 × 10 8 and the mass is 1.398 × 10 5 kg, representing 14.61% of the total mass of all the ejecta. Figure 3a shows the relationship between the ejecting velocity and the distance of the ejecting point from the centre of the crater. Red indicates that the ejecting velocity is less than 1.2 × v esc , which is the part that we will simulate. Green indicates the ejecting velocity greater than 1.2 × v esc , which will not be simulated because at this ejecting velocity, the ejecta would quickly escape from the binary system. Figure 3b gives a diagram of the impact ejecting, again in red for the simulated Red indicates that the ejecting velocity is less than 1.2 × v esc , which is the part that will be simulated. Green indicates the ejecting velocity greater than 1.2 × v esc , which will not be simulated.
Note that the minimum ejecting velocity is the orbital velocity of the secondary. (b) is the diagram of the impact ejecting, the red is the simulated part, and the green is the unsimulated part (we hide the green part where y < 0 for better expression) (d) is to detect whether the ejecta E is inside or outside the quadrangular pyramid OACD, where E 1 is inside the polyhedron, E 2 is outside the polyhedron part and in green for the non-simulated part. Only the ejecta at the edge of the crater has a small enough ejecting velocity that we will use as a sample for simulation. For impact events, the duration of ejecta escaping the crater is mostly a few hundred milliseconds, much smaller than our simulation time. Therefore, in our work, all the ejecta are considered to escape the crater at the same time.

Second collision modeling
Some of the ejecta will collide with binary asteroids during the motion, here we propose a quick collision detection model that can efficiently locate on which face of the polyhedron the ejecta are projected without using a traversal approach. A traversal approach is to detect whether the ejecta are in contact with each tetrahedron formed by each face of the polyhedral and the centre of mass of the polyhedral. This method requires a large-scale loop calculation, which is easy to implement but time consuming. So the studies for this usually reduces the size of the loop, such as the study of Jiménez and Segura (2008). Our collision detection model is suitable for that any ray from the centre of mass of the asteroid (called O here and after) has only one intersection with the asteroid model outline. Assuming that the line connecting the ejecta E and the centre of mass of the asteroid O is OE, and its intersection with the asteroid's outline is P , the principle of our collision detection model is to compare the lengths of OE and OP . First, the vertex coordinates of the polyhedral model are interpolated to generate a polyhedral model with uniform distribution of latitude and longitude, and the longitude lines and latitude lines are numbered ( Fig. 4a and Fig. 4b). Then convert the Cartesian coordinates of the position of ejecta to spherical coordinates (θ e , φ e , r e ), and use Eq. (5) to calculate which of the two lines of longitude and latitude the projection of the ejecta on the model surface lies between (Fig. 4c). Here we can determine that the ejecta are projected within ABCD, and the next task is to determine whether it is projected within ABD or ACD. Assume that the projection point is P . P falling within triangle ACD needs to satisfy: Assuming that the projection point is inside triangle ACD, we finally need to determine whether ejecta E is inside the quadrilateral cone OACD by where N is the normal vector of ACD. As shown in Fig. 4d, − − → AE 2 projects positive on N , outside the polyhedron; − − → AE 1 projects negative on N , inside the polyhedron.
We then processed the collision, and since the surface characteristics of the Didymos are not well known and the incident and reflected velocities of the ejecta are most likely not specular reflection relations, we considered the diffusive of the ejecta according to the Cercignani-Lampis-Lord (CLL) model (Cercignani and Lampis 1971). Assuming the incident velocity of the ejecta is (u i , v i , 0), where u i is the normal velocity with respect to the plane and v i is the corresponding tangential velocity. The reflected velocity consists of two components, specular velocity (u s , v s , 0) and diffusive velocity (u d , v d , w d ), the direction of w d is perpendicular to u d and v d . The specular velocity is calculated by where α n is the normal coefficient of restitution and α t is the tangential coefficient of restitution. The diffusive velocity is calculated by where r and θ is a random number between 0 and 1, and The final reflection velocity is the vector sum of specular velocity and diffusive velocity where μ ∈ [0, 1], when μ = 1, particles complete specular, when μ = 0, particles complete diffusive. Figure 5 shows the reflection of 500 particles impacting on a plane under the same incident conditions. The main reason for introducing the CLL collision reflection model is the lack of a complete understanding of the asteroid surface. If a specular reflection model is used, the surface of the asteroid is divided into many surfaces (the current model is of low resolution and has poor performance at characterizing the shape of the surface of the asteroid) and the reflection model provided by a particular surface is consistent, i.e., some ejecta impact on the same surface with the same output velocity magnitude and similar output velocity direction, and these ejecta reflections will all be in a very similar state. In practice this is rare over a large area, the CLL collision reflection model can provide a terraininduced uncertainty about the impact.

Simulation stop condition
The simulation stops when all ejecta escape the binary asteroids' gravitational influence sphere or stop on the surface of one of the binary asteroids. For the former condition, it is sufficient to check whether the distance of the ejecta relative to the centre of mass of the binary asteroid is greater than the radius of the gravitational influence sphere. Here we discuss in detail the conditions for the ejecta to stop on the surface of the asteroid.
Assuming that the ejecta collide with the asteroid, one of the conditions under which it can remain on the surface of the asteroid is that the reflection velocity in the normal direction is extremely small. Due to the low spin angular velocity of the secondary, the ejecta can stay at any position on the surface of the secondary and we can use only one condition to judge the situation from the surface of the secondary. However, the spin angular velocity of the primary is so fast that the ejecta cannot stay in some parts of surfaces, and we can determine whether the ejecta can stay locally by calculating the relative acceleration on the surface of the primary. The relative acceleration vector can be obtained by subtracting the implicated acceleration and Coriolis acceleration from the absolute acceleration. The following vectors are in the T coordinate system if not otherwise specified. Absolute acceleration can be calculated by where F ,E is the combined force of the solar gravitational force and solar radiation pressure on the ejecta. The implicated acceleration can be calculated in the following formula. The position of the ejecta can be expressed as where R P is the position of the primary, r e is the position of the ejecta relative to the primary. According to Eq. (13) where Q is the directional cosine matrix of the primary and ω is the angular velocity in the coordinate system A. Then a e =r =R P + Qω × (Qω × r e ) Fig. 6 The specific method of ejecta landing on the surface in our model wherer is the acceleration of the primary, α is the angular acceleration of the primary in coordinate system A, and ω is the antisymmetric matrix of ω (Q ωω = 0). Thus the relative acceleration where a C is Coriolis acceleration. Assume that the local normal vector is n (in the coordinate system T ), if the ejecta can stay locally, it is necessary to satisfy a n r = a r · n < 0 where a n r is the projection of the relative acceleration in the n direction.
The specific processing method is shown in Fig. 6, where the simulation step size is 1 s, the solid line outline indicates that the ejecta are outside the surface of the asteroid, and the dotted line outline indicates the inside of the asteroid surface.
Step 1 (Fig. 6a) shows the state of the ejecta flying in orbit before the collision; Step 2 ( Fig. 6b) shows that the ejecta collide with the surface of the asteroid, and here we assume that its reflection velocity is almost zero. At this time, we need to check whether the current surface can the ejecta stay (whether a n r < 0). If the result is yes, the ejecta are considered to have stopped on the asteroid, and the subsequent position relative to the asteroid surface remains unchanged (Fig. 6c); if the result is no, we artificially place the ejecta on the outside of the asteroid surface (due to the small step size, the artificially raised distance is small), and make its velocity to be 0 relative to the local surface, then it will continue to fly in the orbit (Fig. 6d).

Ballistic propagation of ejecta in terms of ejecta sizes
In this section, we perform seven simulations, including one general simulation (GS) and six special simulations (SS). GS1 follows the ejecta initialization model in Sect. 2.3, we set the simulation x (x is the distance from ejecting position to the centre of crater) to be 10.827 m to 11.412 m (in the case of x < 10.827, the initial velocities of the ejecta are greater than 1.2v esc , which will escape the binary system) and the ejecta size to be 1 mm to 10 m. We set both the normal coefficient of restitution and the tangential coefficient of restitution between ejecta and asteroid to 0.5. Due to computing resource constraints, here we use 1 × 10 5 particles.
In order to study the evolution of ejecta of different sizes (diameter), we performed six special simulations (noted as SS1-SS6). For each case, we use 10 5 particles with diameters of uniform distribution within 10 −5 − 10 −4 m, 10 −4 − 10 −3 m, 10 −3 − 10 −2 m, 10 −2 − 10 −1 m, 10 −1 − 10 0 m and 10 0 − 10 1 m respectively. The specific simulation scheme is shown in Table 3. Figure 7 shows the ejecta fate of the (GS1). At the initial moment, a part of the ejecta stays on the surface of Dimorphos because they are too close to the crater edge so that their initial velocities relative to the secondary are too small, while the other part enters the orbital space of the binary system. Here we define that ejecta are considered to have escaped if it moves out of the sphere of influence of the binary system. Part of the ejecta entering orbit will gradually escape from the binary system due to its excessive initial velocity or the influence of solar radiation pressure; the other part will experience some collisions with the asteroid and eventually stay on the two asteroids. In our general simulation, the ejecta all stopped moving within 13 days, and the final distribution is shown in Table 3. Figure 8 shows the fate of the ejecta of SS1-SS6, and a very clear feature is that as the ejector size becomes larger, Fig. 7 Ejecta fate of GS1, the lines of four colours represent the percentage of ejecta staying on primary, staying on secondary, flying in orbit and escaping out of the binary system the number of escapes decreases significantly and the survival time in orbit becomes longer, the final distribution is shown in Table 3. The size of the ejecta in SS1 is so small that most of them escape out of the binary system quickly under the influence of solar radiation pressure, except for those that stay on the surface of the secondary at the beginning, and no ejecta even land on the surface of the primary. For SS5 and SS6, some of the ejecta survive in orbit for a long time, even more than 400 days. In our simulation model, the ejecta size variation affects the magnitude of solar radiation pressure. Here we introduce an index d, which is defined as the critical ejecta diameter d for a certain position, at which the effect of the binary system gravitational force and the solar radiation pressure force is equal. If the ejecta size is smaller than d, the influence of solar radiation pressure dominates, otherwise, the gravitational force of the binary asteroids plays a major role. Figure 9a shows the distribution of d in the sphere of influence of the binary system (viewed from the normal direction of the orbit), and the farther away from the mass centre of the binary system, the stronger the effect of solar radiation pressure. Figure 9b shows the distribution of d near the centre of mass of the binary system. It can be seen that there is a larger d area near Dimorphos because the gravitational forces in this region from the two asteroids are in opposite directions. Therefore, the overall effect is weak, making the solar radiation pressure dominant. The index d also indicates the minimum size at which ejecta can survive in orbit for a long time, if the size of the ejecta are smaller than d, they will finally escape from the binary system under the influence of solar radiation pressure.
We randomly selected 5000 ejecta as a sample to observe the effect of initial ejecting velocity on the final distribution (Fig. 10). The green dots represent the ejecta escaping from the binary system, the red dots represent those landing on the surface of the primary, and the blue dots indicate those landing on the surface of the secondary. It can be clearly seen that ejecta with small enough initial ejecting velocity fall directly on the surface of the secondary, while ejecta with large enough initial ejecting velocity escape, and other ejecta will eventually land on the surface of the primary or secondary after orbital motions. The distribution of green dots can also prove that our initial ejecta velocity threshold of less than 1.2 × v esc is reasonable.

Coefficient of restitution influence
Uncertain mechanical properties of asteroid surface regolith would introduce different coefficients of restitution. With a  Fig. 11 Ejecta fate of GS2, the meanings of different colour lines are consistent with Fig. 7 small coefficient of restitution, ejecta lose more energy and are more likely to land on the asteroid's surface after collision compared to a larger coefficient of restitution. Under the simulation conditions of GS1 and SS1-SS6, we adjusted the tangential coefficient of restitution and the normal coefficient of restitution to be 0.05, and perform 7 simulations (denoted as GS2 and SS7-SS12). The final distribution of the ejecta for GS2 and SS7-SS12 is given in Table 4, and the ejecta fate of GS2 is shown in Fig. 11. The results show that when the coefficient of restitution becomes small, the proportion of ejecta landing on the asteroid is larger. In addition, a smaller coefficient of restitution leads to a lower reflection velocity of the ejecta, which further reduces the survival time of the ejecta in the orbit. Figure 12 shows the orbital survival times of ejecta of different sizes under these two collision coefficients of restitution. For ejecta of size 10 −5 − 10 −4 m, the coefficient of restitution has no effect on orbital survival time, because the ejecta are so small that they are quickly cleared by the solar radiation pressure and escape the binary system without colliding with the asteroid surface. For larger size ejecta, a larger coef-ficient of restitution increases the number of collisions and thus makes the orbital survival time longer. Therefore, the time of ejecta surviving in orbit is sensitive to the surface properties of the asteroid.

Dependence on the tumbling rotational state of Dimorphos
The present observations show that Dimorphos is in tidal self-locking with Didymos, i.e. the angular velocity of the rotation of Dimorphos is equal to the angular velocity of the revolution. The impact can destabilize this state, and it was found that Dimorphos is more prone to destabilization in a rotational state around its long axis (Agrusa et al. 2021). This tumbling rotational state is related to the shape of the secondary and the momentum transfer efficiency. In this section, we discuss the effect that the tumbling rotational state of the secondary would have on the ejecta on the surface. We set the momentum transfer efficiency to 3 here under the applicable model described above. We assume that the impact takes place at the extreme edge of the long axis of the secondary (which can cause the largest change in the angular velocity of the secondary) and that the change in angular velocity is where I is the impulse acting on the secondary, l is the arm of the force of the impact to the centre of mass of Dimorphos (here we take the longest semi-long axis of the secondary), J Z S is the rotational inertia around its axis of rotation, and U k is impact velocity of kinetic impactor and m k is the mass of the kinetic impactor. For the description of the attitude of the secondary, we use the 1-2-3 Euler angles, defined as φ, θ and ψ respectively, and the configuration of the system is shown in Fig. 13a. Using the parameters of Eq. (18), the impact may cause a 50% change in angular velocity, so we changed the ratio of the rotation angular velocity to the revolution angular velocity of the Dimorphos from 1:1 to 1.5:1. The variation of the Euler angle is shown in Fig. 13b, which shows that at about the 8th day the secondary librates more than 90°around its long axis (φ), this is a situation we call 'flip'. We use the simulation conditions of GS2 to calculate the evolution of the ejecta in this state, denoted as SS13. Table 5 gives the final distribution of the ejecta and their orbital survival time. The results show that the secondary tumbling rotational state does not significantly affect the final distribution of the ejecta, because the increase in the spin angular velocity caused by its tumbling rotation does not make its surface ejecta degree of shedding. However, this tumbling rotational state may trigger mass movement from on the surface of the secondary, which we can study by calculating the relative acceleration of the material on the surface of the secondary with respect to the local surface (please refer to Sect. 2.1 for the specific calculation formula).
The distribution of the relative acceleration from the surface of the secondary in the case of self-locking to the primary is given in Fig. 14a (note that here we have converted the relative acceleration a r to the coordinate system B), the arrows indicate the tendency of the matter to move and the colour represents how stable the ejecta are in the local area, with warmer colours indicating that the ejecta are more likely to stay in place; cooler colours indicating that the ejecta are more likely to slide. The specific concrete definitions of the colour values are as follows: where a t r is the projection of the relative acceleration in the local surface.
The results show that in the latitudinal direction there is a tendency for material to slip towards the equator; in the longitudinal direction there is a tendency for material to slip towards the two apexes of the long axes of secondary. The boundary between the eastward and westward slip of material occurs mainly at longitudes of ±90°, 0°and ±180°. The region with the greatest tendency to slip is distributed in the north direction and south direction at two apexes of the long axes. Figure 14b shows the relative acceleration from the surface of the secondary at a given moment in the 'flip' state. The results show a clear change in the distribution of λ, with the material on the surface of the secondary becoming more susceptible to sliding, and a change in the location of the most susceptible regions. There is also a change in the boundary between the eastward and westward sliding of the material, which would allow the surface shape to change over the long evolution of the secondary.

Dependence on the internal structure of primary
The internal structure of the primary of the Didymos system has so far been little studied, and its plausible effect may have on the deposition of the ejected material is unknown.
In this section, we focus on the influence of the internal mass distribution of the primary on the fate of the ejecta. We keep the mass constant and adjust the internal structure of the primary, using four models. (a) Internal density distribution is uniform. (b) The density decreases from the inside to the outside, with the centre density set to 3.2 g/cm 2 (equivalent to the density of granite) and the surface density set to 1.6 g/cm 2 (equivalent to the density of sand soil). We generate the model by setting the node density varying linearly with its distance from the centre of mass of the model, the closer to the centre of mass, the greater the node density.
(c) the internal structure is hollow. We generate the model by setting the node density to 0 for nodes that near the centre of mass of the model. (d) The internal structure contains some holes, we implement it by removing some elements and replacing them with voids; The four models correspond to Fig. 15a, b, c and d respectively.
We started SS14 -SS16 separately, using the same initial conditions for each simulation, to calculate the evolution of 100,000 ejecta with a size distribution ranging from 1 mm to 10 m. Table 6 gives the final distribution of the ejecta and the orbital survival time. The results show that the internal structure of the primary has no significant effect on the evolution of the ejecta. The number of ejecta landing on the surface of the primary is almost identical.
The primary of the Didymos system has a very fast spin rate, causing ejecta in some regions to be thrown off the surface and into orbit. The internal structure of the primary can have an effect on the distribution of these regions, resulting in a different distribution of ejecta on the surface of the primary. We analyzed the mechanical environment at the surface of the four primary models shown in Fig. 15 to study the regions where ejecta can rest on the surface of the primary. Using the method of Sect. 2.5, we calculated the local normal component a n r of the relative acceleration on the surface of the primary, which if positive indicates that the ejecta cannot stay here and will take off. The distribution of a n r for the four models of the primary described above is given in Fig. 15 The internal structure of the primary. (a) Uniform internal mass distribution; (b) Density decreases from inside to outside; (c) Hollow structure; (d) Internal randomly distributed holes Fig. 16 Distribution of a n r on the surface of the primary for the four internal structures corresponding to Fig. 15, where the regions surrounded by the black line indicate that the ejecta cannot stay local  Fig. 16, which are due to the gravitational field errors generated by the vertices of the polyhedron model. This will be corrected later when a more accurate model is obtained. The results show that ejecta close to the equator cannot stay on the surface. The distribution of regions differs under the four models, suggesting that the internal structure of the primary influences the deposition distribution of the ejecta on the surface of the primary.

Conclusions
We developed a detailed model for numeric simulation of the post-impact motion of a binary-ejecta system. For the interaction between ejecta and two asteroids, we propose a fast collision detection model, which can greatly improve the efficiency of collision detection. The CLL model is used for collision processing, taking into account the diffusiveness of ejecta. With large-scale simulations, we get the evolution of ejecta with diameters from 0.01 mm to 10 m. The ejecta with lower velocity will land directly on the surface of the secondary, while others will enter the orbit of the binary system and eventually escape from the binary system or land on the surface of the two asteroids. We found that the ejecta less than 0.1 mm in diameter are more susceptible to escape quickly under the influence of solar radiation pressure, while the larger size ejecta will survive in the orbit for a long time. We examined the survival time of the ejecta in orbit with two coefficients of restitution, and the results show that a large coefficient of restitution results in a longer survival time of the ejecta in orbit. We analyzed the effect of the tumbling rotational state of the secondary on the distribution of material on the surface of the secondary after the impact. This state results in a significant change in surface slope of the secondary but without causing the material to slip. Besides, we analyzed the effect of the internal structure of the primary on the distribution of the ejecta on its surface. The result shows the area in which the ejecta can remain varies with its internal structure. In the future, when finer asteroid surface models are obtained, we can add asteroid surface material migration models to our models to study the evolutionary behaviour of binary surface features after impact. We have uploaded the model to GitHub. 1