Velocity Distribution of a Homogeneously Cooling Granular Gas.

In contrast to molecular gases, granular gases are characterized by inelastic collisions and require therefore permanent driving to maintain a constant kinetic energy. The kinetic theory of granular gases describes how the average velocity of the particles decreases after the driving is shut off. Moreover it predicts that the rescaled particle velocity distribution will approach a stationary state with overpopulated high-velocity tails as compared to the Maxwell-Boltzmann distribution. While this fundamental theoretical result was reproduced by numerical simulations, an experimental confirmation is still missing. Using a microgravity experiment which allows the spatially homogeneous excitation of spheres via magnetic fields, we confirm the theoretically predicted exponential decay of the tails of the velocity distribution.

In contrast to molecular gases, granular gases are characterized by inelastic collisions and require therefore permanent driving to maintain a constant kinetic energy. The kinetic theory of granular gases describes how the average velocity of the particles decreases after the driving is shut off. Moreover it predicts that the rescaled particle velocity distribution will approach a stationary state with overpopulated high-velocity tails as compared to the Maxwell-Boltzmann distribution. While this fundamental theoretical result was reproduced by numerical simulations, an experimental confirmation is still missing. Using a microgravity experiment which allows the spatially homogeneous excitation of spheres via magnetic fields, we confirm the theoretically predicted exponential decay of the tails of the velocity distribution.
Introduction.-Granular gases consist of macroscopic particles, i.e. their diameter d is micrometers or larger. From this property follow two important differences between granular and molecular gases. First, all collisions between granular particles are inelastic: a part of the mechanical energy of the relative motion of the particles is converted to heat [1]. Secondly, the energy scales involved in granular dynamics, such as the energy needed to lift a grain its own diameter, are 10 to 20 orders of magnitude larger than the thermal energy of the system. In consequence, maintaining a dynamic state of a granular gas does require the constant injection of energy.
The theoretical analysis of granular gases focuses on two stationary states: First the homogeneously driven gas where energy is injected in a spatially homogeneous way into the system in order to compensate the energy loss due to inelastic collisions. Kinetic theory predicts [35] that the distribution of individual particle velocities P (v) develops overpopulated tails when compared to the Maxwell-Boltzmann distribution: P (v) ∼ exp(−kv 3/2 ) for v > v / where is the coefficient of restitution and v is the average velocity of all particles. Second, the * peidong.yu@dlr.de † schroeter@science-consulting.info ‡ matthias.sperl@dlr.de homogeneous cooling state (HCS) where the system is not disturbed by external forces [35][36][37][38][39][40]. This results in a monotonous decrease of v with time t according to the so-called Haff's law [41] : τ is a material and density dependent characteristic time scale. If the coefficient of restitution of the particles can be considered to be velocity independent, γ equals one. If the velocity-dependence of for viscoelastic particles has to be taken into account, γ becomes 5/6. The reason that the HCS is qualified as a stationary state, even though v is a function of time, is that the distribution of the rescaled particle velocities c = v/v T , becomes stationary. v T is the velocity derived from granular temperature T of the particles: v T = 2T /m where m is the particle mass. As discussed in appendix B, this thermal velocity can be computed as v T = 2/ √ π · v . Kinetic theory [35,36] predicts that the tail of the velocity distribution P (c) becomes exponential, i.e. P (c) ∼ exp(−kc). The HCS ends with the onset of particle clustering which destroys its homogeneity [31].
Experimental confirmation of these kinetic theory results is difficult for two reasons. First, creating and maintaining a granular gas requires constant energy injection. Typically, this is done by vibrating the boundaries of the container, which however creates spatial and temporal inhomogeneities within the granular gas [25][26][27][28][43][44][45][46], thereby invalidating one of the preconditions of the kinetic theory approaches. A spatially homogeneous driving can be obtained by either restricting the experiment to a horizontally aligned two-dimensional system which is shaken vertically [47][48][49][50][51], or by using magnetic [30,[52][53][54][55] or electrostatic [56] forces to drive a three-dimensional system in the bulk.
Second, on earth the gravitational field will always collapse the granular gas into a granular solid where collisions are replaced by enduring contacts. This effect can again be counteracted by limiting the granular gas to a horizontally aligned two-dimensional system. This comes however at the price of a significantly larger number of particle-boundary than particle-particle collisions.
In consequence, experimental confirmations of the main results of the kinetic theory of granular gases have been scarce. A number of groups tried to measure the the velocity distribution of the homogeneously driven gas using boundary-driven, two-dimensional [47][48][49][64][65][66] and three-dimensional [67,68] systems. They reported high energy tails proportional to exp(−kv α ) with values of α in the range 0.8 to 2. Especially for twodimensional systems it has been shown that α is controlled by the friction coefficient of particle-sidewall collisions [69]. The only confirmations of α = 1.5 were obtained in an electro-statically driven three-dimensional system [56], and in a two-dimensional, horizontal system containing rotationally-driven disks [70].
For the homogeneous cooling state, the validity of Eq. (1) has been shown both in numerical simulations [31,[71][72][73] and experiments [30,60,62]. For the exponents α characterizing the tail of the velocity distribution values in the range 0.6 to 1.5 have been reported for boundarydriven, two-dimensional systems [49]. The theoretically predicted exponential tail, i.e. α = 1, has up to now only been seen in numerical simulations [74,75]; an experimental confirmation is missing.
In this letter, we describe a granular gas experiment which is both force-free, due to being performed in microgravity, and homogeneously heated via magnetic excitation. Our results confirm Haff's cooling law and show clear evidence of an exponential tail of the particle velocity distribution, thus directly confirming the validity of the kinetic theory of granular gases in the homogeneous cooling state.
Experimental Setup.-The experiment was performed in a sounding rocket of the MAPHEUS campaign of German Aerospace Center (DLR) which provided 375 s of microgravity with a remnant gravity on the order of 10 −5 times the earth acceleration [54]. The granular gas consists of 2796 spheres with a diameter d = 1.6 ± 0.02 mm. The particles are contained in a cubic container made of PMMA acrylic glass with an inner side length of 50 mm. Taking into account the slanted corners of the cell, the average volume fraction of the gas is φ = 0.05 and the number density n = 0.0234 mm -3 . The interior of the container is connected to the outside atmosphere, resulting in an air pressure < 0.01 Pa during the experiment. The resulting air drag reduces the particle velocity by 0.01% per second [54].
The spheres are made of ferro-magnetic MuMetall (Sekels GmbH). Energy is injected into all granular particles in the bulk by pairs of magnets which are aligned along the diagonals of the cube, cf. Fig. 1. The excitation protocol cycles through the four pairs, switching each one on for 20 ms and then pausing for 80 ms. The low coercivity of the particles ensures that the long-range interactions between the particles are negligible during cooling [54], while during heating, such long-range interactions actually contribute to a more isotropic excitation [55]. As shown in numerical simulations [55] and experiments [54], this type of excitation results in a homogeneously heated granular gas. It can also be seen in the movie of the experiment which is part of the supplementary material.
The granular dynamics inside the cell is observed with light field camera Raytrix R5 at 165 frames/s and a magnification of 70 µm/pixel. Due to its reduced precision, the depth information provided by the light field camera is not used in our analysis. Appendix A provides further information on the two-dimensional image analysis. . Each fit provides two parameters: v0, the average particle velocity at t = 0, and τ , the time by which the velocity has decreased to v0/2.
Freely cooling granular gas.-During the flight of the sounding rocket four individual experiments were performed. Each experiment consisted of a first phase of at least 40 s duration during which the granular gas was heated by switching the magnetic fields on and off. Then the excitation was stopped and a cooling period of at least 40 s duration followed. However, as shown in appendix C, after about 4 seconds the remnant gravitation starts to induce spatial inhomogeneities in the sample cell. We therefore limit the subsequent analysis to the first 3 seconds of the cooling phase; these were sufficient to reach a stationary state as shown in appendix F. Fig. 2 presents the evolution of the average particle velocity in the last second of the heating and the first three seconds of the cooling phase. Two features of Fig. 2 are important to point out. First, Haff's law (Eq. (1), displayed in red) provides a good fit to all four experiments. We have tested both exponents γ = 1 and 5/6; the differences are within our experimental uncertainty. This signifies that the velocity-dependence of plays only a minor role in the velocity range analyzed here. This is to be expected because for (g) = 1 − ζg 1/5 + h.o.t. and within the studied range of relative collision velocities g, the change of epsilon is small (ζ contains material parameters) [1]. We assume in our subsequent analysis a constant = 0.66, which is measured from lab calibrations (see appendix G for details).
Second, the reproducibility is high between the four experiments. The parameters v 0 and τ of the Haff fits agree within ±8%, and v 0 · τ , which is determined by system properties alone, agrees within ±3%. Additionally, the particle velocities averaged over the last two seconds of the heating phase differ only by ± 0.3 mm/s around the mean of the four experiments: 20.9 mm/s. While the exact value of the fit parameter v 0 depends on the stochastic fluctuations within the granular gas during the heating phase, the value of τ can be compared to its prediction from granular kinetic theory [1,35] Here χ(φ) is the contact value of the pair correlation function. We use the usual approximation for hard sphere systems [77]: . v T 0 is the thermal velocity of the particles at the beginning of the cooling: v T 0 = 2/ √ π · v 0 (see appendix B for details). Using the fitted value v 0 = 20.6 mm/s from our first experiment, Eq. (2) predicts τ k.t. = 1.34s, which is a factor of 2.6 times larger than the fitted value τ = 0.507s.
The fact that τ < τ k.t. matches with the limits of Eq. (2) which includes only the energy dissipation in the normal direction of the collisions, and not dissipation in tangential direction which can e.g. be caused by inter-particle friction. Furthermore, inclusion of tangential dissipation means the rotational degrees of freedom join the dissipation mechanism and couple with the translational motion [6]. We expect that theories for rough particles [7,71,78,79] provide predictions closer to the experimental results. However, such a comparison, would also require 3D particle tracking including the rotational motion of the particles. Therefore this question depicts the direction of future experimental improvements. Velocity Distribution. -Fig. 3 shows the distribution of the rescaled particle velocities P (c). The high reproducibility of the experiments allows us to improve the statistics by averaging the distributions of the four experiments; the four individual distributions can be found in appendix D. For values of c smaller than 1.5, P (c) can be approximated by a Maxwell-Boltzmann distribution. However, the main result of Fig. 3 is that the high energy tail of P (c) decays exponentially, as predicted by the kinetic theory of granular gases. A nonlinear fit of the experimental data using P (c) ∼ exp(−kc α ) results in α = 0.96 ± 0.19, which is again in good agreement with kinetic theory.
Due to the intermittent nature of the magnetic driving (i.e. 20 ms on, 80 ms off), the velocity distribution of the heated gas, which is shown in appendix E, is different from a Maxwell-Boltzmann distribution and it does also not display a high velocity tail proportional to exp(−kv 3/2 ) as predicted by kinetic theory for the homogeneously driven gas. A fit of the tail with P (c) ∼ exp(−kc α ) results in α = 0.72 ± 0.21.
The fact that the system still relaxes towards the distribution predicted by granular kinetic theory can be interpreted as the equivalent of Boltzmann's H-theorem for thermal systems [34]. It is essential for the comparison with theoretical predictions that the observed velocity histogram allows for the interpretation as a distribution function including the definition of a well-defined temperature.
Conclusions.-We have shown experimentally that the kinetic theory of granular gases predicts correctly the two most fundamental properties of a freely-cooling threedimensional gas: the evolution of the average particle velocity and the appearance of an exponential tail in the particle velocity distribution. Building on this foundation, we expect quantitative deviations between experiment and existing theory to stimulate further fruitful work.
We acknowledge partial funding from BMWi/DLR through project 50WM1651 and 50WM1945. We also thank an anonymous referee for providing valuable insights into the theoretical implications of our data.
Appendix A PARTICLE TRACKING Fig. 4 displays a typical raw images from which we track the motion of the visible subset of spheres. Usually, Particle Tracking Velocimetry (PTV) identifies in a first step in each image the positions of all visible particles. In a second step the lists of particle coordinates originating from subsequent images are used to construct self-consistent particle trajectories [80][81][82]. This strategy will not work for our raw data because the identification of the particles in an individual image is hampered by three experimental constraints: a) the intensity contrast between particles and background is small. b) the four illumination sources result in a spatially dependent number of reflections on each sphere which prohibits the use of template matching algorithms. And c) the relatively high particle density results in frequent partial occlusion of spheres.
In order to avoid these problems, we have decided to merge identification and tracking into a single step, i.e., to use multiple consecutive images to identify where our particles are. The basic idea can be seen in the left column of Fig. 5. The image in the top row is a cropped version of the image in Fig. 4. The image in the bottom row is a space-time plot which has been constructed as follows: take all 1000 images corresponding to a given experimental run (heating phase and then free cooling) and stack them on top of each other, effectively creating a 3D volume image. Then we cut a single vertical slice from this volume. Which shows how a specific row in one of our images evolves with time (i.e. going from top to bottom in that slice.) The bright streaks in this space-time image correspond to reflections on the sides of the spheres, the dark streaks are the dark areas at the centers of the spheres (because there is no light source at the camera side of the cell.) Streaks start and end because particles also have velocity components in the y-direction, i.e. perpendicular to the single line we are visualizing here. In some cases the ends of the streaks will also correspond to particles appearing or becoming hidden behind other particles. But the main point here is that the extended nature of this streaks allows us to identify moving entities together with their trajectory. While the streaks themselves will not necessarily correspond to the particle centers, the symmetric and diffuse illumination guarantees that we do not make a parallax error if we identify the streak velocity with the particle velocity. In consequence, this method is not suitable to measure pair correlation functions or mean square displacements, but it will produce reliable velocity distributions and mean velocities.
The actual image processing pipeline consist of two steps: preprocessing which creates pixel-thin representations of the streaks and track identification which detects the streaks and converts them into particle velocities. These steps are further described in the next two subsections.

A.1 Preprocessing
Before assembling all images of a given experiment into a 3D image stack, as shown in the bottom left panel of Fig. 5, we crop them to remove all pixels belonging to the cell boundaries.
The next step is to turn the spatio-temporal paths of the reflections and central shades (on the spheres) into separate objects. This is achieved with a double binarisation where we choose the upper and lower threshold such that the darkest 2.9 % and brightest 1.8 % of pixels are converted to white while the remaining background is set to black. The middle column of Fig. 5 provides an example. Changing these thresholds such that the number of pixels retained varies by a factor of 1.5 changes the τ values of the Haff's law fits by less than ±6%. Which is well within the variance between different experimental runs.
The final step is a reduction of the white objects to lines of thickness one pixel. This is achieved using a skeletonization algorithm with a 27 pixel neighborhood in 3D (Skeletonization works by making successive passes through the volume. On each pass, border pixels, i.e. white pixels with at least one black neighbor, are identified and removed on the condition that this does not break the connectivity of the skeleton). The right column of shows the raw images. The four corners of the image are cut out as they contain only the bevels in the corners of the experimental cell. The space-time view reveals the particle motion (in x direction) as inclined bright and dark streaks (originating from the reflections and shadows on the sphere surfaces). The smaller the angle between a streak and the vertical axis is, the slower the corresponding particle moves. The middle column shows the result of a double binarisation with an upper and a lower threshold: the streaks are now separated from the background. The right column corresponds to the middle column after the streaks have been reduced to lines of width one pixel using a 3D skeletonization. Images have been inverted for better visibility.

A.2 Identifying tracks
In micro-gravity particles will move between collisions with constant speed. In our spatio-temporal image stack constant speed translates into straight lines; we therefore want to fit our streaks with straight line segments.
Identifying straight lines in an image is normally done with the so-called Hough transformation which introduces an additional accumulator space with axis which are parameters describing all possible lines. For 2D images these parameters are normally the pair of angle θ and distance r which constitute the Hesse Normal form of a straight line. Then for each white pixel the θ and r pairs of all lines possibly going through this pixel are computed and the corresponding bins in the accumulator space are incremented by one. Once all white pixels have been processed, the maxima in the accumulator space indicate the parameters of the lines on which many white pixels can be found.
Generalizing this idea to 3D requires a number of modifications and expansions. Luckily, most of this work has already been done by Dalitz and coworkers (2017). Their paper is also accompanied by a well documented, opensource implementation in C named houdh3dlines. Fig. 6 demonstrates how we can use this 3D Hough transformation to identify the line segments which form the backbone in the center of the streaks. The parameterization in 3D requires 4 parameters: two angles to describe the direction of the line and two components of a vector going from the origin of the coordinate system to the intersection point between the line and a plane which is perpendicular to that line (and also going through the origin). For a typical experimental run with an image stack of size 700 times 700 pixels times 1000 time steps this translates into a memory requirement of 10.1 GByte for the accumulator space.
After inserting the information contained in the white pixels of our 3D volume into this accumulator space, houdh3dlines returns the parameters of all maxima in that space in the sequence of decreasing bin height, i.e., decreasing line length. It also performs a conventional line fit to the points constituting each maximum in order to increase the accuracy of the four parameters describing the line.
We have slightly modified houdh3dlines for our purpose: we only retain line segments which are based on streaks with white pixels in at least six consecutive time steps. Moreover, we also record the earliest and latest frame contributing to that line segment. From the angle between the line direction vector and the time axis we can compute the x and y component of the velocity of the underlying particle. In the subsequent analysis this velocity is considered to be present between the start and end time of this line segment.

Appendix B RELATIONSHIP BETWEEN MEASURED AVERAGE VELOCITY v AND THE THERMAL VELOCITY vT
It is clear from the previous section that the particle velocity we measure is the magnitude of the velocity vector projected from a three-dimensional volume into a two-dimensional image plane. Therefore, the relevant Maxwell-Boltzmann (MB) distribution becomes the 2D version: . For the calculation of the Haff's time (Eq.2 of the paper) we need the thermal velocity v T 0 = v T | t=0 at the begin of the cooling. As described above, v T 0 can be computed from the mean velocity measured at the starting of the

Appendix C CLUSTERING AFTER 4 SECONDS
The spatial homogeneity of the particle distribution is a consequence of the magnetic driving. Once the gas enters the free cooling phase, the remnant gravitational acceleration will start to create inhomogeneities. Figures  7 and 8 demonstrate that denser cluster begin to form at the "top" and "right" edge of the sample cell after approximately 4 seconds.
Based on figures 7 and 8, we can try to estimate the remnant acceleration. If we assume that particles migrate about 1 cm in 8 s, the remnant acceleration would correspond to 3 × 10 −5 g, which is in agreement with the specification of the sounding rocket flight. Fig. 9 shows the rescaled velocity distributions of the four individual experiments. The average of the four distribution is shown in Fig. 3 of the paper.

Appendix E PROPERTIES OF THE HEATED GRANULAR GAS
It has been pointed out, based on numerical simulations, that the evolution of the velocity distribution from a normal Maxwell-Boltzmann distribution towards the solution predicted by kinetic theory can take up to the order of ten times τ [40].
We find in our experiment that the fully developed velocity distribution can already be found after three times τ . This is a beneficial side effect of the 1:4 duty cycle of our driving protocol; during the 80 ms where the magnets are switched off, the system can already start to relax towards the HCS. This becomes evident from Fig. 10 where the velocity distribution in the heated gas is already approaching the HCS when compared to the Maxwell-Boltzmann distribution.
However, there are still significant differences between the heated and cooling velocity distribution; especially in the shape of the high velocity tail: Figure 11 shows that the tail of the heated distribution can not be described by an exponential decay.

Appendix F THE COOLING GAS REACHES A STATIONARY STATE
Given the limited time frame of cooling before clustering sets in, the question could be raised if we measure indeed the asymptotic distribution. This question is answered in figure 12 where we have split the time range shown in figure 3 in the paper (1.8 to 3 s) into two halves. Figure 12 then demonstrates that the velocity distribution does not change between these two intervals which corroborates the existence of a steady state. The heated case has a higher probability of particles with velocities in the range 1 < c < 1.5 and also a different exponent of the tail, as shown in figure 11. However, the heated velocity distribution is clearly more similar to the cooling distribution than a Maxwell-Boltzmann distribution.

Appendix G THE COOLING TIME SCALE FOR VISCOELASTIC PARTICLES
As mentioned in the manuscript, in addition to the Haff's law for constant , that for = (g), i.e., v(t) = v 0 /(1 + t/τ ) 5/6 for viscoelastic particles, has also been used to fit the experimental cooling curves (see Fig. 13). The resulting time scale τ can then be compared with theories.
Eq. 3.22, 9.14 and 9.29 of [1] establishes how τ depends on the viscoelastic collision property = (g). We rewrite them using the conventions in the main manuscript.
As illustrated in Fig. 14, used in the main manuscript is retrieved from a simple particle dropping experiment. Two MuMetall particles from the same batch used in the 0g experiment are prepared, with one completely fixed by glue to a solid surface below, and another temporarily attached to a nozzle above (a). The nozzle's position is finely adjusted, so that when the particle above drops (b), it collides with the one below almost  vertically (c). This close-to-vertical collision minimizes the rotation after the particle above is bounced back (d).
The translational velocities before and after the collision are calculated from image analysis and lead to the result of = 0.66 for a collision velocity g = 368mm/s. For viscoelastic particles, (g = 368mm/s) = 0.66 fully determines Aκ 2/5 in the first line of Eq. 3. The second and the third lines then predict, for our low-gravity cooling, a time scale of τ v.s. = 0.982s.
Comparing with the viscoelastic fit result shown in Fig. 13, we find that τ v.s. is again more than 2.5 times larger than the experimental value (e.g., τ = 0.344s from the 1st cooling). Therefore, viscoelastic collision theory does not interpret our experiments quantitatively better compared to the simpler kinetic theory assuming to be constant.