# Collisional Cooling of a Charged Granular Medium

## I Introduction

The particles in most granular materials carry a net electrical charge. This charge emerges naturally due to contact electrification during transport or is artificially induced in industrial processes. It is well known , for instance, that particles always charge when transported through a pipe. In industry, contact electrification is used for dry separation of different plastic materials or salts , which tend to get oppositely charged and hence are deflected into opposite directions when falling through a condenser. Another application is powder varnishing, where uniformly charged pigment particles are blown towards the object to be painted, which is oppositely charged.

Whereas the dynamics of electrically neutral grains have been studied in great detail, little is known about what will change, if the grains are charged. In this paper we present the answer for collisional cooling, a basic phenomenon, which is responsible for many of the remarkable properties of dilute granular media. By collisional cooling one means that the relative motion of the grains, which lets them collide and can be compared to the thermal motion of molecules in a gas, becomes weaker with every collision, because energy is irreversibly transferred to the internal degrees of freedom of the grains.

In 1983 Haff  showed, that the rate, at which the kinetic energy of the relative motion of the grains is dissipated in a homogeneous granular gas, is proportional to $`T^{3/2}`$, where $`T`$ is the so called granular temperature. It is defined as the mean square fluctuation of the grain velocities divided by the space dimension:

$$T=\stackrel{}{v}^2\stackrel{}{v}^2/3.$$
(1)   
A consequence of this dissipation rate is that the granular temperature of a freely cooling granular gas decays with time as $`t^2`$. We shall discuss, how these laws change, if the particles are uniformly charged.

Due to the irreversible particle interactions large scale patterns form in granular media, such as planetary rings  or the cellular patterns in vertically vibrated granular layers . This happens even without external driving , where one can distinguish a kinetic, a shearing and a clustering regime. The regimes depend on the density, the system size and on the restitution coefficient $`e_\mathrm{n}=v_\mathrm{n}^{}/v_\mathrm{n}`$, which is the ratio of the normal components of the relative velocities before and after a collision between two spherical grains. The $`T^{3/2}`$ cooling law holds, provided the restitution coefficient may be regarded as independent of $`v_\mathrm{n}`$ , and the system remains approximately homogeneous . The latter condition defines the kinetic regime, which is observed for the highest values of the restitution coefficient, whereas the two other regimes are more complicated because of the inhomogeneities. Such inhomogeneities can only occur as transients, if all particles are equally charged, because the Coulomb repulsion will homogenise the system again.

In order to avoid additional dissipation mechanisms due to eddy currents within the grains we consider only insulating materials. Unfortunately, up to now, no consistent microscopic theory for contact electrification of insulators exists . In powder processing two types of charge distribution are observed : A bipolar charging, where the charges of the particles in the powder can have opposite sign and the whole powder is almost neutral. The other case is monopolar charging, for which the particles tend to carry charges of the same sign and the countercharge is transferred to the container walls. It depends largely on the type of processing, whether one observes bipolar or monopolar charging, which means, that the material of the container, the material of the powder and other more ambiguous things, like air humidity or room temperature are important .

The outline of this article is as follows: The next section specifies the model we are considering. A simple derivation of the dissipation rate in dilute charged granular media based on kinetic gas theory is given in section III. We find, that the dissipation rate is essentially the one known from uncharged granular media multiplied with a Boltzmann factor. Section IV compares the analytic results with computer simulations. We find that in non-dilute systems the Coulomb repulsion is effectively reduced. This reduction will be explained, and we determine its dependence on the solid fraction of the granular gas. In the appendix we discuss the new simulation method we developed for this investigation. It is a molecular dynamics method, that avoids the so called brake-failure .

## II The model

In this paper, we consider monopolar charging, which is the usual case if insulators are transported through a metal pipe . For simplicity we assume, that all particles have the same point charge $`q`$ centred in a sphere of diameter $`d`$ and mass $`m`$. No polarisation and no charge transfer during contact will be considered. The particle velocities are assumed to be much smaller than the velocity of light, so that relativistic effects (retardation and magnetic fields due to the particle motion) can be neglected. The electrodynamic interaction between the particles can then be approximated by the Coulomb potential:

$$\mathrm{\Phi }_{ij}=q^2/r_{ij},$$
(2)   
where $`r_{ij}`$ is the distance between the centers of particles $`i`$ and $`j`$.

We consider the collisions as being instantaneous, which is a good approximation for the dilute granular gas, where the time between collisions is much longer than the duration of the contact between two particles. As the incomplete restitution ($`e_\mathrm{n}<1`$) is the main dissipation mechanism in granular gases, Coulomb friction will be neglected in this paper. Also, the dependence of the restitution coefficient on the relative velocity  will be ignored, so that the constant $`e_\mathrm{n}`$ is the only material parameter in our model.

The particles are confined to a volume $`V=L^3`$ with periodic boundary conditions in all three directions. The periodic volume can be thought of as a sufficiently homogeneous subpart of a larger system, which is kept from expanding by reflecting walls. For vanishing particle diameter this model corresponds to the One Component Plasma (OCP) . In the OCP a classical plasma is modelled by positive point charges (the ions) acting via the Coulomb potential, whereas the electrons are considered to be homogeneously smeared out over the whole system. In the OCP the electron background cannot be polarised, i.e. Debye screening does not exist, as is the case in our model, too.

## III Analytical results for dilute systems

In this section we derive an approximate expression for the dissipation rate in a dilute system of equally charged granular spheres, neglecting particle correlations. Basically we apply the kinetic gas theory, but include inelastic collisions. Using the analytic form of the dissipation rate in the dilute limit derived here, we will discuss the dissipation in a non-dilute system, where correlations are important, in the next chapter.

We start with calculating the collision frequency of a fixed particle $`i`$ with any of the other particles $`j`$. If they were not charged, two particles would collide provided the relative velocity $`\stackrel{}{u}`$ points into the direction of the distance vector $`\stackrel{}{r}=\stackrel{}{r}_j\stackrel{}{r}_i`$ connecting the particle centers, $`\stackrel{}{u}\stackrel{}{r}>0`$, and the impact parameter $`b=|\stackrel{}{r}\times \stackrel{}{u}|/u`$ is smaller than the sum of the particle radii, $`bb_{\mathrm{max}}=d`$. If the particles carry the charge $`q`$, they repel each other and the maximum impact parameter $`b_{\mathrm{max}}`$ becomes smaller than $`d`$ (see Fig.1). By the conservation laws for angular momentum and for energy one gets:

$$b_{\mathrm{max}}^2=d^2\left(1\frac{2E_\mathrm{q}}{\mu u^2}\right)$$
(3)   
where $`\mu =m/2`$ is the reduced mass. $`E_\mathrm{q}=q^2/d`$ denotes the energy barrier which must be overcome to let two particles collide in the dilute limit. It is the difference of the potential energies at contact and when they are infinitely far apart. Eq. (3) is independent of the actual form of the potential, as long as it has radial symmetry. (Note that energy is conserved as long as the particles do not touch each other.)

Imagine a beam of particles, all having the same asymptotic velocity $`\stackrel{}{u}`$ far away from particle $`j`$. All particles within an asymptotic cylinder of radius $`b_{\mathrm{max}}`$ around the axis through the center of $`j`$ with the direction of $`\stackrel{}{u}`$ will collide with particle $`j`$. There will be $`\pi b_{\mathrm{max}}^2un`$ such collisions per unit time, where $`n=N/V`$ is the number density. Integrating over all relative velocities $`\stackrel{}{u}`$ gives the collision frequency of a single particle in the granular gas in mean field approximation:

$$f=\pi n\underset{uu_0}{}d^3uub_{\mathrm{max}}^2(u)p(u).$$
(4)   
$`u_0=\sqrt{2E_\mathrm{q}/\mu }`$ is the minimal relative velocity at infinity, for which a collision can occur overcoming the repulsive interaction. We assume that the particle velocity distribution is Gaussian with variance $`3T`$ (see (1)), so that the relative velocity will have a Gaussian distribution as well, with

$$u^2=6T$$
(5)   
Hence, the total number of binary collisions per unit time and per unit volume is given by:

$$\dot{N}_g=1/2fn=2\sqrt{\pi }n^2d^2\sqrt{T}\mathrm{exp}\left(\frac{E_\mathrm{q}}{mT}\right)$$
(6)   
The factor $`1/2`$ avoids double counting of collisions. This corresponds to textbook physics for chemical reaction rates as can be found for example in Present.

Now we calculate the dissipation rate: The energy loss due to a single inelastic collision is:

$$\delta E(u,b)=\frac{\mu }{2}\left(1e_\mathrm{n}^2\right)u_{\mathrm{n}}^{}{}_{}{}^{2}$$
(7)   
where $`u_\mathrm{n}^{}`$ means the normal component of the relative velocity $`\stackrel{}{u}^{}`$ at the collision. It can be calculated easily from $`u_{\mathrm{n}}^{}{}_{}{}^{2}=u_{}^{}{}_{}{}^{2}u_{\mathrm{t}}^{}{}_{}{}^{2}`$: The tangential component is determined by angular momentum conservation,

$$\mu ub=\mu u_\mathrm{t}^{}d,$$
(8)   
and energy conservation gives

$$u_{}^{}{}_{}{}^{2}=u^2\left(1\frac{2E_\mathrm{q}}{\mu u^2}\right).$$
(9)   
This yields

$$u_{\mathrm{n}}^{}{}_{}{}^{2}=u^2\left(1\frac{b^2}{d^2}\frac{2E_\mathrm{q}}{\mu u^2}\right)$$
(10)   
The energy loss in one collision is therefore:

$$\delta E(u,b)=\frac{\mu }{2}\left(1e_\mathrm{n}^2\right)u^2\left(1\frac{b^2}{d^2}\frac{2E_\mathrm{q}}{\mu u^2}\right)$$
(11)   
Assuming a homogeneous distribution of particles, we eliminate the $`b`$-dependence by averaging over the area $`\pi b_{\mathrm{max}}^2`$ (see Fig.1):

$`\delta E(u)`$ $`=`$ $`{\displaystyle \frac{1}{\pi b_{\mathrm{max}}^2}}{\displaystyle \underset{0}{\overset{b_{\mathrm{max}}}{}}}𝑑b\mathrm{\hspace{0.17em}2}\pi b\delta E(u,b)`$ (12)   
$`=`$ $`{\displaystyle \frac{\mu }{4}}u^2\left(1e_\mathrm{n}^2\right)\left(1{\displaystyle \frac{2E_\mathrm{q}}{\mu u^2}}\right)`$ (13)   
The dissipated energy per unit time due to collisions with relative velocity $`u`$ is then the number of such collisions per unit volume, $`1/2n^2\pi b_{\mathrm{max}}^2u`$, times the energy loss $`\delta E`$, Eq. (13).

Finally we get the dissipation rate per unit volume in the dilute limit ($`\nu 0`$) by integration over the relative velocity distribution:

$`\gamma `$ $`=`$ $`{\displaystyle \frac{\pi }{2}}n^2{\displaystyle \underset{uu_0}{}}d^3ub_{\mathrm{max}}^2u\delta E(u)p(u)`$ (14)   
$`=`$ $`2\sqrt{\pi }n^2d^2m\left(1e_\mathrm{n}^2\right)T^{3/2}\mathrm{exp}\left({\displaystyle \frac{E_\mathrm{q}}{mT}}\right)`$ (15)   
The dissipation rate of an uncharged granular system in the dilute limit in the kinetic regime is given by :

$$\gamma _0=2\sqrt{\pi }n^2d^2m\left(1e_\mathrm{n}^2\right)T^{3/2}$$
(16)   
Thus the dissipation rate (15) in a monopolar charged dilute granular gas and the one for the uncharged case differ only by a Boltzmann factor, $`\gamma =\gamma _0\mathrm{exp}\left(E_\mathrm{q}/mT\right)`$. This is the main result of the analytic treatment in this section. It remains valid for any repulsive pair interaction between the grains that has rotational symmetry.

## IV Dissipation rate for dense systems

In order to discuss the dissipation rate $`\gamma `$ in a non-dilute system of charged granular matter, let us recall the analytic form of $`\gamma `$ in an uncharged non-dilute system. The derivation is basically done by using the Enskog expansion of the velocity distribution function for dense gases . One gets the dissipation rate for a non dilute uncharged system:

$$\gamma =\gamma _0g_{\mathrm{hs}}(\nu )$$
(17)   
where $`\gamma _0`$ is given by Eq. (16) and $`g_{\mathrm{hs}}(\nu )`$ is the equilibrium pair distribution function of the non-dissipative hard-sphere fluid at contact. It only depends on the solid fraction $`\nu =\pi nd^3/6`$:

$$g_{\mathrm{hs}}(\nu )=\frac{2\nu }{2(1\nu )^3}$$
(18)   
(Carnahan and Starling , Jenkins and Richman ).

Our system consists of dissipative charged hard-spheres (CHS). The Boltzmann factor in Eq. (15) is just the equilibrium pair distribution function at contact in the dilute limit for a CHS-fluid, $`lim_{\nu 0}g_{\mathrm{chs}}(\nu ,q)=\mathrm{exp}\left(E_\mathrm{q}/mT\right)`$. So it is plausible, that the dissipation rate for a dense system of dissipative CHS is:

$$\gamma =\gamma _0g_{\mathrm{chs}}(\nu ,q)$$
(19)   
Unfortunately the literature is lacking a satisfying analytic expression for $`g_{\mathrm{chs}}`$. In 1972 Palmer and Weeks  did a mean spherical model for the CHS and derived an analytic expression for $`g_{\mathrm{chs}}`$, but this approximation is poor for low densities. Many methods  give $`g_{\mathrm{chs}}`$ as a result of integral equations, that can be solved numerically. We do not use those approximations, but make the following ansatz for $`g_{\mathrm{chs}}`$:

$$g_{\mathrm{chs}}(\nu ,q)g_{\mathrm{hs}}(\nu )\mathrm{exp}\left(\frac{E_{\mathrm{eff}}(\nu )}{mT}\right)$$
(20)   
As in the dilute case we assume that the long range Coulomb repulsion modifies the pair correlation function of the uncharged hard sphere gas by a Boltzmann factor. Note that the granular temperature enters the pair correlation function only through this Boltzmann factor. The hard core repulsion is not connected with any energy scale, so that the pair correlation function $`g_{\mathrm{hs}}`$ cannot depend on $`T`$. The effective energy barrier $`E_{\mathrm{eff}}`$ must approach $`E_\mathrm{q}`$ in the dilute limit. Hence the ansatz (20) contains both the uncharged and the dilute limit, (17) respectively (15).

In order to check the ansatz (20) we did computer simulations using the MD algorithm as described in the appendix. Test systems of varying solid fraction $`\nu `$ and particle number ranging from $`N=256`$ to $`N=1024`$ were prepared at a starting temperature $`T_0`$. As soon as the simulation starts, the granular temperature drops because of the inelastic collisions. We measured the dissipation rate $`\gamma `$ and the granular temperature during this evolution. According to Eq. (20) and Eq. (19) the dissipation rate is $`\gamma =\gamma _0g_{\mathrm{hs}}(\nu )\mathrm{exp}(E_{\mathrm{eff}}(\nu )/mT)`$. An Arrhenius plot ($`\mathrm{ln}(\gamma /\gamma _0g_{\mathrm{hs}})`$ versus $`E_\mathrm{q}/mT`$) should give a straight line whose negative slope is the effective energy barrier $`E_{\mathrm{eff}}`$. Fig. 2 shows two examples of these simulations. The Arrhenius plots are linear to a very good approximation. This confirms the ansatz (20). Systems with high densities show slight deviations from linearity.

The negative slopes $`E_{\mathrm{eff}}/E_\mathrm{q}`$ in Fig. 2 are smaller than $`1`$, which means, that the effective energy barrier is smaller than in the dilute system. The explanation is that two particles which are about to collide not only repel each other but are also pushed together by being repelled from all the other charged particles in the system.

For dimensional reasons the effective energy barrier to be overcome, when two particles collide, must be of the form

$$E_{\mathrm{eff}}=\frac{q^2}{d}\frac{q^2}{\mathrm{}}f(d/\mathrm{}),$$
(21)   
where $`\mathrm{}>d`$ is the typical distance between the charged particles and $`f`$ is a dimensionless function. The first term is the Coulomb interaction $`E_\mathrm{q}`$ of the collision partners at contact. The second term takes the interaction with all other particles in the system into account. It is negative, because the energy barrier for the collision is reduced in dense systems.

Obviously, for a dense packing, $`\mathrm{}d`$, the energy barrier for a collision must vanish, i.e. $`E_{\mathrm{eff}}|_{d=\mathrm{}}=0`$. Moreover, if one takes a dense packing and reduces the radii of all particles infinitesimally, keeping their centers in place, all particles should be force free for symmetry reasons. Therefore, the energy barrier must vanish at least quadratically in $`(\mathrm{}d)`$, i.e. $`E_{\mathrm{eff}}/d|_{d=\mathrm{}}=0`$. For the function $`f`$ this implies

$$f(1)=1\text{and}\frac{\mathrm{d}f(x)}{\mathrm{d}x}|_{x=1}=1.$$
(22)   
If the particle diameter $`d`$ is much smaller than the typical distance $`\mathrm{}`$ between the particles, the function $`f(d/\mathrm{})`$ may be expanded to linear order,

$$f(x)=c_0+c_1x+\mathrm{}$$
(23)   
In linear approximation the coefficients are determined by (22): $`c_0=2`$ and $`c_1=1`$. This determines the energy barrier (21).

In 1969 Salpeter and Van Horn  pointed out, that inside a strongly coupled OCP a short-range body centered cubic (BCC) ordering will emerge. In the BCC lattice the nearest neighbour distance $`\mathrm{}`$ is related to the volume fraction $`\nu `$ by

$$\frac{d}{\mathrm{}}=\frac{2}{\sqrt{3}}\left(\frac{3}{\pi }\nu \right)^{1/3}1.14\nu ^{1/3}.$$
(24)   
Assuming a BCC structure and using the linear approximation for $`f(x)`$ in (21), the effective energy barrier is therefore given by

$$E_{\mathrm{eff}}=E_\mathrm{q}\left(12.27\nu ^{1/3}+1.29\nu ^{2/3}\right)$$
(25)   
To test Eq. (25) we simulated systems with densities ranging from $`\nu =0.001`$ to $`\nu =0.216`$ and determined the ratio $`E_{\mathrm{eff}}(\nu )/E_\mathrm{q}`$ as in Fig. 2. The results are plotted in Fig. 3. The agreement of the theoretical formula (25) with the simulations is excellent. One can see, that in the dilute limit the effective energy barrier extrapolates to $`E_\mathrm{q}`$. We cannot simulate low density systems, because collisions are too unlikely.

For the highest densities one cannot expect that the linear approximation (23) remains valid. Also, the dense packing of spheres is achieved with an FCC (face centered cubic) rather than a BCC ordering. This may be responsible for the systematic slight deviation from the theoretical curve in Fig. 3.

The reduction of the Coulomb repulsion was also found in the OCP, when it was applied to dense stars . There the analogue of the second term in (21) is called the “screening potential” (somewhat misleadingly, as there is no polarizable counter charge and hence no screening). Monte Carlo simulations  of the OCP were interpreted in terms of a linear “screening potential” , which corresponds to (23), and the analogue of the conditions (22) also occurs in the plasma context , although based on a different physical reasoning. Corrections to the linear approximation are the subject of current research . However, applying these more sophisticated forms of the “screening potential” of the OCP model to dense charged granular gases seems arguable as for higher densities the influence of the hard spheres become more and more important and so the analogy to the OCP model, which uses point charges, does no longer hold.

## V Discussion

We derived the dissipation rate of a charged granular gas in the dilute limit. Compared to the uncharged case it is exponentially suppressed by a Boltzmann factor depending on the ratio between the Coulomb barrier and the granular temperature. This result was obtained assuming a Gaussian velocity distribution, although it is known that in the uncharged case deviations from a Gaussian behaviour emerge due to the inelastic collisions . These deviations, however, were shown to have little effect on the dissipation rate . As the system becomes less dissipative in our case, it is reasonable to expect that the effect of deviations from a Gaussian velocity distribution will be even weaker. One may say that a dilute granular gas with monopolar charging is more similar to a hard sphere gas in thermal equilibrium than a neutral one.

In a dense system particle correlations enter the collision statistics and hence the dissipation rate in two ways: First there is the well known Enskog correction as in the uncharged case. It describes that the excluded volume of the other particles enhances the probability that two particles are in contact. Second the Coulomb barrier which colliding particles must overcome is reduced and will vanish in the limit of a dense packing.

In our simulation we did not observe the shearing or the clustering instability, probably because our systems were rather small. It is reasonable to expect, however, that shearing or clustering instabilities may at most exist as transients in the presence of monopolar charging, because the Coulomb repulsion will homogenise the system in the long run.

## Acknowledgements

We thank Lothar Brendel, Alexander C. Schindler and Hendrik Meyer for useful comments. We gratefully acknowledge support by the Deutsche Forschungsgemeinschaft through grant No. Wo 577/1-2.

## A Computer Simulation Method

Distinct element (or molecular dynamics (MD)) simulations  are usually done with time step driven or event driven algorithms . None of them is well suited for investigating a charged granular gas. Therefore we developed a new simulation scheme, which combines the virtues of both and will be described in this section.

We use a “brute force” MD algorithm, which is simple and sufficient for our problem. More sophisticated ways of dealing with the long range interactions, such as the multipolar expansion , the particle-particle-particle-mesh  or the hypersystolic algorithms  should be used, if larger systems need to be studied.

The event driven method for simulating the motion of all particles in the granular gas can be applied, whenever the particle trajectories between collisions can be calculated analytically, so that the time interval between one collision event and the next can be skipped in the simulation. Obviously this is impossible in a system with long range Coulomb interactions. However, the idea to avoid the detailed resolution of a collision event in time is still applicable. So the velocities of the collision partners are simply changed instantaneously to the new values predicted by momentum and angular momentum conservation and an energy loss determined by the restitution coefficient. We shall keep this feature of event driven simulations.

In the time step driven simulation method the equations of motion of all particles in the granular gas are discretized using a fixed time step, which is small compared to the duration of a collision. Hence each collision, which is modelled as an overlap between particles, is temporally resolved in detail. This has the advantage, that the formation of long lasting contacts between particles can in principle be simulated realistically. If the particles carry equal charges, however, this will not happen, so that the collisions may be approximated as being instantaneous like in event driven simulations. Apart from being more efficient, this automatically avoids the so called brake-failure artifact , which hampers time-step driven molecular dynamics simulations with rapid relative motion. On the other hand, we need a time discretization of the particle trajectories between collisions, in order to take the changing electrostatic interactions properly into account.

Because of the long-range nature of the Coulomb potential, we have to include the interactions with the periodic images of the particles in the basic cell. One way to do this is by Ewald summation. Details of this method can be found in . Another method is the minimum image convention: Only the nearest periodic image is taken into account for the calculation of the interactions. The minimum image method has the advantage, that it is much faster than the Ewald summation. We checked the validity of the minimum image method compared to the Ewald summation and found, that as long as $`E_\mathrm{q}/mT<10`$ both methods yield indistinguishable results. This upper limit for the coupling has been found before in Monte-Carlo simulations of the OCP . As our systems all satisfy this condition, we used the minimum image convention.