Cooling of a granular gas mixture in microgravity

Granular gases are fascinating non-equilibrium systems with interesting features such as spontaneous clustering and non-Gaussian velocity distributions. Mixtures of different components represent a much more natural composition than monodisperse ensembles but attracted comparably little attention so far. We present the observation and characterization of a mixture of rod-like particles with different sizes and masses in a drop tower experiment. Kinetic energy decay rates during granular cooling and collision rates were determined and Haff’s law for homogeneous granular cooling was confirmed. Thereby, energy equipartition between the mixture components and between individual degrees of freedom is violated. Heavier particles keep a slightly higher average kinetic energy than lighter ones. Experimental results are supported by numerical simulations.


Introduction
Exploring the behavior of granular gases in microgravity environments holds an immense scientific and practical significance.This area of research advances our understanding of physics, engineering, and even space exploration.From the viewpoint of fundamental physics, ensembles of individual macroscopic particles colliding in a manner similar to gas molecules offer unique perspectives on the basic laws of multiparticle physics.Microgravity allows to observe and analyze pure granular interactions, eliminating the complex influence of gravity.This can lead to breakthroughs in understanding particle dynamics, energy dissipation, and entropy production.
In astronomy and cosmology, understanding behavior of granular systems sheds light on the formation and dynamics of celestial bodies, such as asteroids, comets, planetesimals and planetary rings.One can learn a lot about the way these objects evolve and interact.The study of granular gases can also provide valuable insights into energy dissipation as well as energy and heat transfer mechanisms, with manifold implications for the design of efficient applications on Earth and in space.It should be emphasized that microgravity experiments with granular gases can serve as captivating educational tools, inspiring students in physics, space science and computer vision.
Whereas the majority of previous experiments and theoretical investigations of granular gases was focused on monodisperse systems, more realistic studies have to take into account that such systems in general are polydisperse.Mixtures introduce an additional layer of complexity.One of the fundamental questions is the partition of kinetic energies among the constituents, and among their different degrees of freedom.Additionally, one may investigate the emergent behaviors and self-organization in long-term experiments with these mixtures.
This study is one of the first steps towards a generalization of granular gases to polydisperse systems.We report experiments and numerical simulations of bidisperse mixture and compare energy partition, dissipation (granular cooling) and collision statistics observed in microgravity experiments and numerical simulations.
Already in the 1980s, Peter Haff [30] proposed a scaling law for the kinetic energy loss of a homogeneously cooling dense granular gas of frictionless monodisperse spheres.He predicted a time dependence of the form which yields the scaling E kin ∝ t −2 for times t ≫ τ H . E 0 is the kinetic energy at time t = 0.The Haff time τ H (E 0 ) for a given initial state defines a time scale of energy loss.It depends on material properties and the shape of the grains, and on other system parameters.
Haff's law was confirmed for free cooling dilute ensembles of monodisperse rod-like particles [29], spheres [26] and oblate ellipsoids [27].Estimates for the mean kinetic energies and τ H were obtained.For mixtures, an important and unresolved question is whether and how the components differ in their kinetic and cooling parameters.Granular gas mixtures of spherical particles with different radii and masses were analysed in numerical simulations by Bodrova et al. [31,32].Heating and cooling were considered, yet particle rotations were disregarded for simplicity.Larger, heavier particles were found to have a higher granular temperature than lighter ones.The particular distribution of kinetic energies depends on the choice of the contact model in the simulation.Different granular temperatures for different components are quite common in granular matter, they were, e.g., also predicted for vibrated granular fluids [33].Experimental data for a confirmation of these theoretical predictions were not available so far.

Experimental setup and parameters
In the present paper, we analyze the dynamical properties of a two-component granular gas mixture.Our experimental setup follows the description given in Ref. [29] and is presented on Figure 1: Two cameras viewing along axes y and z were used for a stereoscopic observation of a granular gas consisting of 384 rods (192 of each of sort) in a container with dimensions L x × L y × L z = 11.2 × 8.0 × 8.0 cm 3 .The rods consist of insulated copper wire pieces of length ℓ = 10 mm and diameters d 1 = 0.75 mm (component 1) and d 2 = 1.35 mm (component 2), respectively.The mass of the thin rods is m 1 = 22 mg, their moment of inertia for rotations about the long axis is J ∥1 = 0.99 pN m s 2 and perpendicular to it J ⊥1 = 183 pN m s 2 .For the thicker particles, the mass is m 2 = 37.5 mg, J ∥2 = 4.6 pN m s 2 and J ⊥2 = 315 pN m s 2 .For easier detection and tracking, 48 rods of each type were colored, the rest serves as "thermal" background.The experiment was performed in the ZARM drop tower in Bremen, where microgravity (µg) is achieved for about 9.2 seconds.Initially, the system was excited by sinusoidal vibration of two side walls of the container (in x direction) at an amplitude of A = 0.24 cm and frequency f = 30 Hz, i. e. a maximum plate acceleration of ≈ 8 g.Then, vibration was stopped and the granular gas mixture was left without energy input in a granular cooling phase.
3D particle detection in the experiments, tracking and trajectory post-processing mainly follow the outlines of Refs.[8,29].The rod mixture is excited during the first two seconds of µg, then it undergoes granular cooling.

Numerical simulations
In order to support the experimental analysis, we performed a numerical simulation.It provides the opportunity to extract particular properties that are not accessible in the experiment.At the same time, comparison of the simulation results with the experiment helps to choose a suitable collision model and realistic material parameters.We used a hybrid GPU-CPU implementation of discrete element modelling (DEM) [34][35][36][37][38], adapted to systems with moving walls.The collision detection and Hertz-Mindlin contact force model follows the previous simulations with spherocylinders [9,37,38].In addition, our program was modified to simulate the bidisperse mixture of particles with the same length and diameters as in the experiment.The energy loss was quantified using an effective restitution coefficient e n = 0.7 and friction coefficient µ = 0.4.These values provide reasonable agreement with the experiment, even though the computed Haff times are somewhat larger.The Young modulus was set to Y = 5 GPa, yet its particular value does not influence the results noticeably as long as it is large enough to restrict excessive particle overlap, and small enough so that the collisions are resolved in a sufficient number of simulation steps (δt = 5 × 10 −8 s) [39].Container and excitation parameters were the same as in the experiment.

Cooling rate and Haff time
Figure 2 shows the decay of the average kinetic energy per particle separately for the two components, thin (E 1 ) and thick (E 2 ) rods.Time is counted from the start of granular cooling after the excitation was switched off.Vertical lines represent the uncertainty which arises mainly from the limited ensemble size of evaluated rods [29].From the experimental data, two separate runs (drops) are presented.One observes satisfactory agreement between both experimental runs and simulation1 .In the subsequent figures, the data for both experimental runs are combined.The mean energies per particle at the end of the heating phase are E 1 ≈ 580 nJ and E 2 ≈ 660 nJ (see Fig. 2). Figure 3 shows that Haff's law, Eq. ( 1), well reproduces the loss of kinetic energy during granular cooling for both experimental and simulated data, except for the initial 0.2 s after heating was switched off.There, it fails because the heated state is too inhomogeneous [29], and this is the reason why the fit parameter E 0 does not represent the actual initial energy.Most importantly, fitted Haff times τ H differ only slightly for the two components, namely, τ H1 = 0.45 ± 0.02 s and τ H2 = 0.47 ± 0.02 s in the simulation, while τ H1 = 0.29 ± 0.03 s and τ H2 = 0.34 ± 0.03 s in the experiment.Within the statistical accuracy of our data, this is consistent with equal Haff times for both components.If a slight systematic deviation actually exists, it means that the heavier particles initially cool slower until an equilibrium is reached with equal Haff times and the energy distribution at later stages is slightly shifted further away from equipartition.The reported difference of Haff times between experiments and simulation is not problematic.Note that the Haff time has the property τ H (t 0 +t) = τ H (t 0 )+t, so that the difference of the initial Haff times in experiment and simulation appears as a simple time shift of 0.15 s between the experimental and simulated data (cf.Fig. 2).We first compare the shares of kinetic energies per individual degree of freedom (DOF) of each rod type.An excess of translational energy along the 'active' direction x is observed for both components during excitation (Fig. 4).The average energy for the directly excited DOF can reach more than half of the total kinetic energy.This is in accordance with previous findings [16,24,28,37,40,41] and arises mainly from 'hot' particles directly after collisions with vibrating walls.In addition, we find a granular temperature gradient towards the center of the container in the heated x-direction, similar to experiments with excited dense granular ensembles [42].The velocity distribution functions are non-Gaussian.The continuously heated phase was not the primary focus of this study.A thorough analysis of energy balance in continuously heated systems deserves a more dedicated experiment.

Energy partition
Figure 4 shows the energy partition between DOF associated with translational motion (E x , E y , E z ) as well as two rotational DOF about the short rod axes (E Rot⊥ ).After the onset of cooling, the share of energy associated with translation along x rapidly decreases.We obtain kinetic energies per DOF which are close to equipartition, except for slight residual dominance of translations along x.Energy partition for cooling monodisperse rod ensembles was studied in Refs.[28,29]: After a relatively short initial period, the particles were found to reach a steady partition of E  and E Rot⊥ .The average energy associated with the rotational DOF was found around 10-20% less than for translations.
The kinetic energy E Rot∥ contained in the 6 th degree of freedom, rotation around the long rod axis, is quite difficult to determine experimentally, and was not accessible in the present experiments.An approximate value was given in Ref. [28] where it was estimated that E Rot∥ is about one order of magnitude lower than the mean energies of the other DOF.In our simulations, all DOF are directly accessible.The energy partition extracted from the simulation (Fig. 4 c,d) yields energy levels for E Rot∥ contributing around 10% to the total kinetic energy at the begin of cooling, and even slightly growing afterwards.This is significantly more than the value obtained earlier experimentally [28].These rotations about the long axis are exclusively excited by frictional contacts.Trusting the experiment, we presume that the realization of these frictional contacts in the simulation is not yet satisfactory and needs refinement.
A crucial question for granular gas mixtures is the dependence of the mean kinetic energy on particle properties such as relative particle size and mass.In our experiments, the mass ratio is m 1 /m 2 = 0.59 (diameter ratio d 1 /d 2 = 0.56).An excess of average kinetic energy of the larger particles persists during the complete cooling phase, as seen from the ratio E 1 /E 2 ≈ 0.8 of the average total (observed) kinetic energies per particle for the mixture components in Fig. 5(a    Fig. 6 Cumulative statistics of particle-particle (blue markers) and particle-wall collisions (green markers).Solid lines are fits with logarithmic functions of the form given in Eq. ( 2).Double collisions of the rods with a wall [43] count as one.
The statistics of collisions as the elementary steps of the cooling process are compared in Figures 6 and 7. Figure 6 shows the cumulative number of collisions of each particle type with other particles and with the container walls, extracted directly from the simulations, as well as from an analysis of experimental particle trajectories.Following Haff's model, the average cumulative number of particle-particle collisions in the system can be approximated by [29,40,44] where time t starts with the beginning of cooling.Here, C p is a positive constant related to the mean energy loss in a single particle-particle collision [29,44].Since the number of collisions is the ratio of the mean speed of the particles and some timeindependent geometric parameter (i.e. the mean free path) both particle-particle and particle-wall collisions should obey the same logarithmic dependence as in Eq. ( 2), with different prefactors C p and C w .The four curves in Figure 6 were fitted by Eq. ( 2) with mean initial Haff times τ H = 0.31 s (experiment) and 0.46 s (simulation).While the simulation collision numbers are accurate and match the fits quite nicely, the experimental results differ significantly.The possible reason is that not all collisions could be correctly identified in the videos with our automated collision detection, both in the initial phase where particles move too fast and in the later stages of cooling where noise in the detected trajectories can be mistaken for collisions.A more accurate collision detection would be possible in experiments with somewhat lower packing fractions and an improved camera setup.Nevertheless, it is evident that particle-wall collisions are around 2.5 − 3 times less frequent than particle-particle ones for the thinner rods and around 3.5 − 4 times less frequent for the thicker rods.
For a comparison, we modified the formula for the collision cross-section of rods [29] to include cylinders of different diameters.A random particle orientation respective to the flight direction and a uniform, homogeneously mixed particle distribution are assumed for simplicity.The number of rods of each type be N 1 and N 2 , resp.Then, we obtain estimates of the mean free paths λ 1 and λ 2 as: where the scattering cross sections σ ij corresponding to collisions of particle types i and j can be approximated as: The estimated mean free paths are λ 1 = 1.92 cm and λ 2 = 1.66 cm, resp., in absence of wall collisions.A rough estimate of the mean free path of a sphere in a cubic container with sides ℓ yields λ w ≈ 0.583 ℓ, which for our experiment is roughly 2.7 λ 1 or 3.2 λ 2 .Thus, the predicted ratio of collisions with particles and with walls should be around 3, in fairly satisfactory agreement with the experiment.Another detail of the collision statistics is elucidated in Fig. 7.When the particleparticle collisions counted in the simulation are compared to the expected value from v/λ with the mean particle velocity v1,2 and the mean free paths λ 1,2 from Eqs. (3,4), a factor of ≈ 1.5 is evident.The mean free paths directly extracted from the simulations are also ≈ 1.5 times shorter than their theoretical estimates.A possible reason is that the simplified collision model considers only non-rotating rods.Fast rotations about the short axes obviously increase the scattering cross sections.

Summary
Summarizing, a first quantitative confirmation of the main features of the granular cooling of a bidisperse mixture was achieved here.Heavier particles have larger mean energies both during excitation and during cooling.The ratios E 1 /E 2 of average kinetic energies per particle were found to be around 0.8 both in experiment and simulation (with the mass ratio being 0.59 and the moment of inertia ratio 0.58 for the rotations about the short axes).The Haff time determined for heavier particles was found to be slightly larger than that of the lighter ones, which might indicate a slower cooling of the heavier particles until an equilibrium E 1 /E 2 is reached, that is smaller than the value reached within our observation time scale.In perspective, this study shall be extended to bidisperse mixtures of particles which are significantly different to each other [45], particularly with different shapes.The investigation of polydisperse mixtures is highly desired, but several problems including the automatic detection and distinction of components need to be solved first.

Methods
The ensemble was observed with two video cameras GoPro Hero 3 Ribcage.The image resolution is 1280 × 960 pixel 2 at a frame rate of 100 fps.Altogether, 8 different colors were available with 12 particles of each color.The remaining 144 non-tracked particles of each of the two types are black and serve as a 'thermal' background.The choice of 96 colored particles in each experiment is approximately at the limits of the detection and tracking feasibility with the current setup and packing fraction.
For particle detection and tracking, we followed the workflow which was described in [9].The program was significantly improved, transferred to Detectron 2 framework [46], and a custom GUI interface was added for preview and correction of 2D and 3D data.
A dataset which contains camera image frames together with the corresponding rod endpoints was assembled, first manually, later with an iterative procedure for corrected rods.The current dataset includes around 2500 images (around 300 000 object instances of colored rods) from both camera views from several independent runs of the experiments.Data processing scripts and pre-trained Detectron 2 network model files, as well as the GUI program for correction and annotation of data are being prepared for publication as an open-source software package.
In order to track the colored rods between the frames, rods endpoints are triangulated and matched by solving the 3D axial optimal assignment problem (also known as bipartite graph matching).We found that an optimization towards both reprojection error for rod endpoints and displacements of rod endpoints between frames is required for robust tracking of particles.
The trajectories of the rods were fitted taking into account constant translation velocities and rotation rates during the free flight phases of the particles between the collisions.For the rod centers of mass, the l 1 trend filtering optimization method [47,48] was used.It allows to fit the noisy particle center trajectories into a sequence of piecewise linear functions with kinks (bends of fitting function) in between.
We extracted angular velocities by differentiating the rod orientation quaternions [34,49,50].To smoothen the high-frequency noise due to the orientation measurement error, a moving average filter was applied.Then, the resulting angular velocity data were fitted with the similar trend filtering approach, only instead of the piecewise linear approximation we fit the angular velocities with piecewise constant (step) functions.
For translational as well as angular velocities, we employ the special case the l 1 trend filtering for vector time series as described in [47].The advantage to fitting the x, y, z center coordinates together is that the fitted coordinate components tends to show simultaneous trend changes at common kink points, which correspond to change of velocities due to collisions of the particle.One problem that arises when applying this procedure is that as particles slow down, the relative weight of the kinks in the optimization procedure decreases and the standard l 1 trend filtering procedure tends to overestimate the number of collisions.To improve the fitting, we have employed an iterative weighted heuristics from [47] which allows to optimize the fit towards the number of kinks instead of sum of their residual norm.
We assume that the optimization error which arises at the bending points of piecewise linear approximations of particle endpoint velocities, as well as jumps in angular velocities, signalize that the particle collided with another particle or the wall.Then, knowing the positions of the wall, collisions can be attributed to either particle-particle or particle-wall collision.This way, the collisions undergone by rods can be determined automatically, even though the absolute precision may not always be satisfactory.

Fig. 1 a
Fig. 1 a) Experimental setup: Inside the container (experiment box) is a mixture of elongated grains which can be excited by two vibrating walls.It is observed by two cameras though front and top transparent walls.b) Example of video frame which shows granular gas mixture during its cooling in microgravity.The image width is about 12 cm.

Fig. 2
Fig.2Total kinetic energy for the two mixture components from experiment and simulation in comparison.Time refers to the stop of excitation (vertical dashed line), 2 s after entry into the microgravity phase.

Fig. 3
Fig.3Total kinetic energy for the two mixture components in double logarithmic scale.Both are fitted well for t > 0.2 s with Eq. (1) and the parameters given in the graphs.

Fig. 4
Fig. 4 Partition of the mean particle energies per DOF, (a,b) experiment, (c,d) simulation.Here and in following figures, vertical dashed lines mark the start of cooling.

Fig. 5 a
Fig.5a) Ratio of the total kinetic energies E 1 , E 2 per particle for the mixture components 1 and 2. The red dashed line marks the start of cooling.b) Ratios of translational (E T 1 , E T 2 ) and rotational (E R1 , E R2 ) kinetic energies per particle for the mixture components.

Figure 5 (
Figure 5(b)  shows separately the ratios of translational (E T 1 /E T 2 ) and rotational (E R1 /E R2 ) energies.The data are more noisy than for the total energy due to the permanent energy exchange between translational and rotational DOF.Nevertheless, we observe satisfactory statistical agreement between experiment and simulation.

Fig. 7
Fig.7Average particle-particle collision frequencies ν p1,2 .Orange line shows the ratio of collision frequency observed in simulation to one obtained from mean velocity and theoretical mean free path estimate.