Self-charging of identical grains in the absence of an external field

We investigate the electrostatic charging of an agitated bed of identical grains using simulations, mathematical modeling, and experiments. We simulate charging with a discrete-element model including electrical multipoles and find that infinitesimally small initial charges can grow exponentially rapidly. We propose a mathematical Turing model that defines conditions for exponential charging to occur and provides insights into the mechanisms involved. Finally, we confirm the predicted exponential growth in experiments using vibrated grains under microgravity, and we describe novel predicted spatiotemporal states that merit further study.

In 1963, the volcanic island Surtsey, named after the legendary fire giant Surtr, rose out of the North Atlantic Ocean. True to its name, over the next year and a half, the island's volcanic debris cloud spat fire in the form of multimillion volt lightning displays 1 . Desert sandstorms similarly have long been known to generate lightning 2 .
How grains generate charge in volcanic plumes, sandstorms or industrial problems such as pharmaceutical mixing 3 or dust explosions 4 remains controversial 5 . Mechanisms have been proposed for charging granular materials that differ in shape, size, or chemical composition [6][7][8][9][10][11][12][13] . However, charging is also observed for materials that are absolutely identical, and experiments show that charging grows with repeated contacts 14,15 . Beyond the unexpected nature of these findings, the fundamentally surprising thing about charging of identical materials is that one appears to get something from nothing: charge that requires energy appears from materials that ought to discharge one another on contact.
Previous work has shown that one origin of this energy can be an external electric field that feeds electrification [16][17][18][19][20] . In the present work, we demonstrate that, remarkably, an external field is not needed: infinitesimal charges on grains themselves can induce charges on their neighbors, bootstrapping one another to grow exponentially rapidly in agitated beds. Unlike prior work 21 , we show that the energy for this charge growth can arise strictly conservatively, trading mechanical work for electrical energy. It is by no means self-evident whether or when this could cause a buildup of significant charges, but as we will show, this feedback mechanism can promote a rich set of macroscopic charge patterns within the bed and can generate very large electric fields.
The underlying mechanism that we explore here is that the electric field from charged particles induces a polarization on a neighbor given by, where χ e is the grain polarizability ranging from zero to unity, R d is an effective dipole radius, → E i is the electric field at the center of the i th grain due to surrounding permanent charges, and k e is Coulomb's constant. In a bed of colliding grains, we also assume that when the i th and the j th particles come into contact, charges q i,k and q j,l on the grain surfaces can partially neutralize 16,17 according to: Scientific RepoRts | 7:39996 | DOI: 10.1038/srep39996 where the prime denotes the constituent charges after collision, η is a neutralization efficiency, and k and l denote the two constituent charges of the i th and j th grains, respectively, that are in contact upon collision. We assume that particles initially have charges of the order of 10 5 elementary charges and we observe a growth by three orders of magnitude as we will show later. Details of the mechanism defined by Eqs (1) and (2) appear in Methods.

Results
We describe a particle-based simulation of charging of agitated grains subject to Eqs (1) and (2) followed by a continuum mathematical model, with the goal of establishing whether induced polarization combined with contact neutralization can predictably amplify small initial charges. We conclude with an experiment confirming that this amplification occurs as predicted.
Discrete element simulation. We begin by performing discrete-element method (DEM) 22 simulations of mechanical and electrical interactions between insulating particles. Typical qualitative charges on grains are shown in Fig.1, and details are presented Methods. In but in summary, we use standard procedures to evaluate forces and torques on particles including both mechanical and electrostatic interactions for particles subject to Eqs (1) and (2). We calculate polarizations and mechanical interactions by embedding within each particle three pairs of orthogonally placed constituent charges (see Methods, Fig. 8) at fixed separations 2R d , where R d is 2/3 of the mean grain radius. These constituent charges take part in neutralization events defined by Eq. (2). Our simulations are 3D and use horizontally periodic boundary conditions and a bottom surface that injects energy by vertically kicking impinging particles. For every grain, we set the initial charge of each of its domains randomly chosen from a uniform distribution between [− 10 −2 , 10 −2 ] pC. We enforce energy conservation to compensate for both the work needed to separate induced charges within each particle and the energy associated with a dipole inducing secondary dipole moments on its neighbors. For the first energy compensation term, we apply a force on particles such that the work done by the force compensates exactly for the energy necessary to induce the corresponding dipoles, and for the second, we reduce dipole moments to account for the energy associated with surrounding field effects. Both terms are detailed in Methods.
To quantify charge growth, we evaluate the evolution of the absolute value of all charges averaged over all grains: as a function of model parameters χ e and η, as shown in Fig. 2 using N = 10 3 particles (see Methods for initialization procedure). Figure 2(a) shows that for large constant polarizability, χ e , q typically does grow roughly exponentially following an initial transient and continuing up to an asymptote that we discuss shortly. We plot q here, but remark that polarizations, and charges of both signs, also grow with the same exponential rate. For fixed neutralization, η, q also exhibits an exponential growth period, however for small χ e , q decreases in time, as shown in Fig. 2(b). Figure 2(c) collects growth rates obtained from the slopes of least-squares fits in linear regions of q log ( ) 10 vs. time plots for χ e and η ranging from zero to one. Evidently, the growth rate increases with both χ e and η. Figure 2 contains several features that we discuss next. First, exponential growth is only seen for sufficiently large χ e and η: Fig. 2(c) suggests that there is an onset criterion for exponential growth. Below this onset, q shrinks monotonically. Second, for most parameter values, exponential growth is preceded by a transient during which q drops. Third, q reaches an asymptotic value for long times. And fourth, Fig. 2(c) indicates an oscillatory "breathing" state that we will describe.
To understand the first two of these observations, we will define a mathematical model that captures the problem's essential dynamics. This will involve some analysis, so we first discuss the simpler asymptotic and breathing states.
We begin with the asymptotic behavior, which provides insight into the bed's charging dynamics. The origin of this behavior can be established by comparing the magnitudes of typical Coulomb and gravitational forces. In Fig. 3, we plot two representative cases, one with moderate charging, (χ e = 0.6, η = 1) and one with rapid charging (χ e = η = 1). Figure 3(a-c) show color-coded charge densities, and Fig. 3(b-d) show corresponding bed charges alongside the ratio R cg between characteristic Coulomb and gravitational forces (defined in the figure caption). From Fig. 3(b), we see that for moderate charging, q reaches a noisy asymptote after about 5 seconds, at which point R cg averaged over the entire bed approaches 20%. Figure 3(a), by comparison, shows that the bed charge is dominated by grains in the middle of the bed: we measure that 70% of the charge is contained in the 14 th through 16 th layers. If we evaluate R cg in these central layers that dominate bed charging, we find that the charge saturates when R cg reaches one, as shown in Fig. 3(b). We conclude that the charging asymptote coincides with R cg approaching one in the fastest charging region. At this point, grains levitate or stick together -either of which will prevent the collisional charging mechanism that we have described.
As for the breathing state identified in Fig. 2(c) (see also Supplemental video 1), the same behavior occurs, but throughout the entire bed. Figure 3(c) shows that for χ e = η = 1, where breathing occurs, the highly charged region extends over most of the bed. In this case, Fig. 3(d) shows that R cg averaged over the entire bed reaches one so the whole bed must levitate or stick together. Indeed, Fig. 3(d) shows that when R cg exceeds one, charging stops, bed charge reduces, and simultaneously the bed contracts. This of course causes densities and collision rates to increase, which in turn must increase charging rates.
We confirm the link between charge oscillations and mechanical breathing in Fig. 4, where we compare power spectra of charge and bed expansion in breathing (χ e = η = 1) and non-breathing (χ e = 0.6, η = 1) states. To evaluate bed oscillations, we average displacements of grain heights, z i (t) from the center of mass height, z c (t): As shown in the main plots of Fig. 4(a,b), for the breathing state, both charges and average displacements of grains oscillate at the same frequency, while as shown in the insets, the non-breathing state exhibits broad spectrum noise, with no dominant frequency and much smaller peaks.
Therefore we propose that the cause of both asymptotic charge and breathing oscillations is that the region of the bed that dominates charging reaches R cg = 1, at which point particles cannot collide and so cannot charge. For moderate charging, this occurs over a limited bed height that we presume cannot lift overlying particles; for strong charging, this occurs over the entire bed, which appears to cause global oscillations.
The breathing state indicates that there can be substantial temporal and spatial variations in charging: moreover our DEM simulations also suggest spatial variations in charge density, as shown in Fig. 5(b-c) (left panels) for two values of χ e and η. In these plots, to identify spatial patterns we enlarge our simulations to 4 × 10 3 grains on a horizontal domain of 20 × 20 mean grain diameters squared.
Mathematical model. To better understand these spatial variations as well as the transient charge growth, we provide a mathematical model for bed charging and discharging. This model, in the long tradition of dynamical systems analysis, is simplified but captures the essential physics that we have described so far. That physics consists of an iterative growth in polarization, defined by Eq. (1), combined with a neutralization in charge, Eq. (2). Accordingly, let us define functions → P x t ( , ) and → x t C( , ) to represent the polarization in units of polarization per length and charge as functions of position, → x , and time, t: , where r min is the mean of the distance to each grain's nearest neighbor, m is the mean grain mass, and g is gravity. (a) Spatiotemporal evolution of bed charges using χ e = 0.6 and η = 1. Inset shows the same plot over longer time, (to 14 sec.). (b) Time evolution of q and R cg (blue: averaged over the entire bed, green: averaged over the layers 14-16), using χ e = 0.6 and η = 1. Note that although R cg is only 20% when averaged over the entire bed, in the fastest charging region, around height = 15, R cg = 1. (c) Spatiotemporal evolution of q using χ e = η = 1: again inset shows same plot over longer time, (to 14 sec.), highlighting periodic oscillations in charge. Notice that for χ e = η = 1, the entire bed becomes highly charged. (d) Here, R cg averaged over the entire bed exceeds one, at which point the charge starts to plummet. Data in (a) and (c) obtained by dividing the bed into one-mean-grain-diameter slices and calculating the sum of the absolute values of charges in each slice.
Scientific RepoRts | 7:39996 | DOI: 10.1038/srep39996 Here, the rate of polarization grows iteratively with electric field in proportion to a single parameter, A, and the charge is reduced in proportion to another parameter, B. We have included diffusivity terms D p and D c to account for migration of the agitated grains, and we anticipate that since polarized grains tend to align and attract one another while charged grains tend to repel, D p should be smaller than D c . We acknowledge that this model is simplified in several respects: it represents polarization as a scalar and neglects nonlinear interactions, particle motion, and Coulomb forces. Nevertheless, as we will show, Eq. (4) reproduces much of the charging dynamics in the bed shown in Fig. 2, it allows us to cleanly explain the charging transient that we have described, it predicts conditions for the onset of exponential charging, and it permits us to analyze spatial and temporal variations in a transparent and comprehensible way. Moreover, Eq. (4) is nothing more than the linear reaction-diffusion equation -arguably the simplest of Turing models, whose features have been fully analyzed elsewhere 23 .
We begin our discussion of this mathematical model by reiterating that the parameter A codifies the exponential rate of growth of polarization, as can be seen in the first line of Eq. (4). Additionally, considering the equations together, we note that -B causes a reduction in charge, so difference A-B determines a net exponential growth. Moreover, neutralization is a physical mechanism that produces charge transfer, which allows particles to gain or lose a net charge. The physical model, Fig. 8, relies on neutralization to permit charge transfer between particles, and the equivalent term in the Turing model, Eq. (4), that allows charge growth is A − B. Accordingly, we will associate A with polarization growth -analogous to χ e in our DEM simulations -and A-B with net growth -analogous to η.
This correspondence is not exact, but as shown in Fig. 6(c), plotting the exponential growth rates of the charge amplitude averaged over all grains, C, as a function of A-B and A produces substantially similar behavior to that seen in Fig. 2(c). Moreover in detail, if we plot the growth in C as a function of time, we obtain panels Fig. 6(a,b) for the values indicated. Again the behavior is similar to the growth seen in DEM simulations, Fig. 2(a,b). Here we We can go further with Eq. (4) and evaluate when exponential growth should be seen and determine the cause of the initial transient observed both in DEM and the Turing simulations.
It is well known that the Turing instability (at which point patterned growth emerges) occurs for Eq. (4) when A · D c − B · D p > 0 23 . In Fig. 6, we set D p = 1 · 10 −5 and D c = 4 · 10 −5 , as we mentioned to mimic the fact that charged grains are expected to diffuse faster than polarized ones. So we expect the onset of growing patterns to occur when A > B/4. We indicate parameters for which this inequality fails -and so growth is not expected -as a violet shaded region in Fig. 6(c).
So from stability analysis of this model, we predict that the onset of charging in the DEM simulations is associated with a simple Turing mechanism -i.e. that polarization builds faster than diffusion or neutralization can destroy it. On the other hand, parameters that do not generate charging are associated with charge transfer (mediated by neutralization) that is too slow to act before grains diffuse away.
The transient in charging is also amenable to analysis using Eq. (4). Our finite difference solutions in Fig. 6 are initialized with small random charges (uniformly distributed on [− 0.025, 0.025]), and since diffusion constants are O(10 −5 ), gradients that would trigger diffusion are negligible. If we remove the diffusive terms from Eq. (4), we obtain a pair of coupled ordinary differential equations (ODEs) with off-diagonal terms (i.e. A · C in the  P equation, and − B · P in the  C equation) that are of the same magnitude as the diagonal terms. This is a recipe for non-normal growth 24 in which transient contraction can occur in a system whose eigenvalues indicate growth, or vice versa.
A simple geometrical interpretation of non-normal growth can be obtained by considering a block of clay that is uniformly expanded in x and y (corresponding to positive diagonal terms in a set of ODEs), but is also sheared (corresponding to nonzero off-diagonal terms). In the long term, the expansion will cause points to move away from one another, while in the short term the shearing can temporarily bring points together.
In the case of Eq. (4), non-normal growth occurs when random orientations of vectors (P, C) re-orient along the expanding eigendirection, approaching smaller values as they do so. This can be confirmed by taking a cluster Apparently the transients shown in Fig. 6(a,b), and presumably as well in Fig. 2(a,b), are associated with a simple mathematical behavior that occurs while small random vectors align with their ultimate eigendirections. Once alignment occurs, exponential growth can proceed, causing strong local heterogeneities, and at this point, the approximation of neglecting diffusion can no longer be made.
To study these charge heterogeneities, we consider patterns produced both by DEM simulations and by the reaction diffusion model Eq. (4), summarized in Fig. 5.
In Fig. 5(a), we show a comparison between transient growth for the DEM (from Fig. 2(a)) and non-normal growth as described previously. In Fig. 5(b), we show a top view of color-coded charges from a DEM simulation that resemble horizontal stripes, along with a comparable Turing pattern from Eq. (4). In all of these plots, red indicates positive charge and blue indicates negative. We emphasize, however, that the dominant patterns seen in solutions of the Turing model consist of irregular spots, as shown in Fig. 5(c). Alongside each of these simulations, we plot color-coded charges from the Turing simulations in 3D.
These plots reveal that for aspect ratio 1:2 (i.e. width = 2 · height), the charge patterns extend throughout the thickness of the bed.
This occurs for all parameter values shown in Fig. 6, although we note that this is a property that results from the D c > D p assumption, which inhibits spatial growth of patterns. As with other reaction-diffusion models, much richer behaviors are also possible: for example relaxing the D c > D p constraint leads to temporally oscillatory charges, strong vertical gradients in charge, and other states such as tubes, labyrinths etc ref. 25. Similarly, changing boundary conditions has a strong and well documented 26 effect on the patterns expressed -for example modeling a tall thin column of grains rather than a broader bed produces horizontal charge striations 21 .
Our primary goal for presenting a simplified Turing model has been to show that Eq. (4) can reproduce charging seen in detailed DEM simulations, as evidenced by comparisons shown in Figs 2, 5 and 6, and we defer comprehensive analysis of this model for future studies. The essential thing that the Turing model appears to reveal is that granular charging follows a stereotypical evolution due to straightforward mathematical causes, including a transient decrease in charge followed by exponential growth, a well-defined onset criterion for growth of charging, and emergence of patterned states (chiefly vertical columns of like charge and polarization). The robustness of the mathematics underlying these behaviors suggests that these findings should be common in practical problems, and calls for experimental verification in laboratory and field experiments.

Experiment
The underlying hypothesis of both DEM and Turing models is that an iterative process leads to exponentially rapid growth of initially infinitesimally small charges, so we close by testing this hypothesis in experiments. In these experiments, shown in Fig. 7, hollow glass spheres are vibrated on a grounded metal plate at 2 kHz by a piezoelectric transducer. The thickness of the particle bed is close to that used in our DEM simulations (under 1 mm, or about 9 particle diameters). These experiments are performed under microgravity (see figure caption), yet as shown in Fig. 7(a), particles return to the plate along curved trajectories. Crucially, since gravity is essentially absent, the only known force that can act at a distance in this way is electrostatic. Moreover, the heights of particle flights diminish with time as can be seen in Fig. 7(a) and in Supplemental video 2, which we can use to obtain a quantitative evaluation of our DEM and Turing predictions, as follows. The maximum height, h, of a particle ballistically ejected from the bed is simply its kinetic energy, KE, divided by the force, F, attracting the particle to the bed: h ~ KE/F. We have seen that our model predicts exponential growth in charges of both signs, so the force, F, associated with these charges must also grow exponentially in time. Since h ~ 1/F, we predict that h will decrease exponentially in time: h ~ e −a·time where a defines the charging rate shown in Figs 2 and 6.
We assess this prediction by evaluating heights reached by particles near the top of the bed. As shown in Fig. 7(b), we horizontally sum the grayscales of pixels from successive 200 ms superpositions (as in Fig. 7(a)). The transducer could drift away from the agitated particle bed and change the particle number density that would affect the brightness of the images. However, the timescales for particle losses qualitatively are on the order of 10 s while the timescales for charging are a fraction of a second. Keeping these two different timescales in mind, we consider the brightness as a suitable measure of the charging.
Our procedure is as follows. High flying outlier particles produce noisy variations in grayscale, so we exclude the noisy region identified in Fig. 7(b), and select a moderate grayscale that shows little noise but provides the largest available height discrimination between superpositions. This grayscale is boxed in the main plot and enlarged in the inset to this figure. We evaluate the grayscale at the center of this region (broken line in the inset),  33 . Gravity is about 10 −6 g and pressure = 1 mbar. Spheres have density 0.14 g/cm 3  which we plot in Fig. 7(c), along with a least-squares fit to the predicted exponential, h = h 0 ± h 1 e −a·time . We find that a fit can be made using a = 1.31 ± 0.03 sec −1 with correlation coefficient, r 2 = 0.997. We repeat the experiment with the metal plate covered with insulating tape, and obtain the height vs. time plot shown in the inset to Fig. 7(c): here we obtain a = 1.40 ± 0.03 sec −1 and r 2 = 0.998. Both of these fits have growth rates, a, in the range expected from Fig. 2(c).
We are aware that other dependencies can also be fitted. However, this would require an alternative competitive idea to be tested, which does not currently exist. Consequently we find that the experimental results are consistent with the theoretical prediction and avoid drawing conclusions about the functional dependence just from the analysis of the experimental data.
These results seem to confirm our predictions of exponential charging, however other possibilities deserve mention. First, most bouncing particles return to the bed, yet some particles near the edges escape (see Supplemental video 2). It might be argued that loss of particles could account for the decrease in bed height, however we note that fewer bed particles would cause particles above the bed both to be less strongly attracted to the bed and to rebound more elastically, both of which would increase, rather than decrease, the measured heights shown in Fig. 7.
Second, it is possible that particles have been tribocharged by the vibrating plate. Although tribocharging doubtless occurs, we remark that (1) particles ejected from the bed are attracted back to the bed, so particles cannot simply be tribocharged with the same sign, which would cause repulsion; (2) spheres landing on grounded metal and on insulating tape produce nearly indistinguishable results; and (3) charging appears to occur exponentially in time. None of these results are consistent with tribocharging as it is traditionally understood 27 . It remains conceivable that particles near and far from the vibrating plate could acquire opposite charges 12 , however this would not explain the apparent exponential charging.
Third, several groups have described charging models for particles differing in size [8][9][10] , and indeed our hollow spheres range from 125 μm to 150 μm. Again, size-dependent charging doubtless does occur, however such mechanisms invariably produce monotonically decreasing charging rates, rather than the exponential growth that we observe.
Based on these considerations, we conclude that our agitated granular bed appears to produce exponential growth in charges of both signs, which to our knowledge our model is unique in predicting.

Discussion
We have performed simulations, modeling, and experiments of charged grains in an agitated bed. The simulations show that grains can charge exponentially rapidly by feeding back their electric fields through their neighbors. The Turing model provides a simple framework to understand the exponential growth in polarization and charge as well as more detailed predictions such as an onset criterion and non-normal charging transients. The microgravity experiments confirm that charging of agitated beds of insulating grains does appear to grow exponentially. Finally, the simulations and the Turing model predict a previously unreported oscillatory state and complex spatiotemporal charging dynamics, both of which merit further study.
We propose that our findings of exponential growth of charging may account for the generation of multi-million volt potentials observed in nature, and may contribute to improved understanding of electrical charging in mining 4 , and industrial powder handling 3,5,12 .
Beyond these findings, we remark that a Turing model contains a natural mechanism for producing "something from nothing," which we mentioned at the start characterizes the emergence of strong electrical effects by doing no more than agitating initially nearly neutral grains. Indeed, historically the appearance of complex dynamics from seemingly benign reactions was rejected as being "impossible" 28 for this reason. In our problem, the "reaction" comes in the form of known electrical effects, polarization and neutralization, but at its core these are not mathematically different from reactions between enzymes or autocatalysts -which also appear to produce something from nothing.

Methods
To account for a heterogeneous permanent charge distribution on the surface of insulating grains, we considered a multipole expansion and truncated after the second-order term. We simulate this by embedding in each grain six independent, orthogonally placed charge domains as shown in Fig. 8. As in previous work 17 , we track translational and rotational motions of each grain by evaluating its dynamics due to mechanical and electrical forces and torques acting on it. We then solve the equations of motion for each grain by means of the DEM on a domain that is periodic along the horizontal directions. We achieve periodicity by surrounding the computational domain by 8 copies in the horizontal plane, which we use for field and energy calculations.
The DEM algorithm itself uses spherical grains with restitution coefficient 0.935 and kinetic friction coefficient of 0.4. We use the model of Walton and Braun 29 for the elastic force and fix the two elastic coefficients, k l = 0.07 and k u = 0.08, to achieve the restitution coefficient, = .
. ≈ . k k / 007/0 08 0 935 l u . We use a time step of 50 msec., which produces over 10 2 steps per collision for the fastest moving grains. We use polydisperse grain sizes to prevent crystallization: the radius of each grain, R i , is Gaussian distributed with standard deviation 10% of the mean radius, = .
R 0 75 mm, and each grain has the density of glass, ρ g = 2.4 g/cm 3 . The top of the computational domain is free, and the bottom is fixed. Any grain that hits the bottom acquires additional upward energy defined by a kick velocity, → = .V g R z 2 7 2 , which maintains the granular bed in a collisional state.
Particles polarize according to Eq. (1), where we emphasize that E&vec; i is the electric field at the center of the i th grain due to all pre-existing permanent charges in the system. The distinction between permanent and induced charges is significant because induced charges are slaved to the external field, and cannot themselves do work. For example, induced charge dipoles always point in the direction of an external field and so cannot exert torque on a grain. Permanent charges, on the contrary, are fixed on a grain and exert forces on other charges 30 . An induced polarization, given as a superposition of multiple orders of electrostatic interaction, could point at any arbitrary direction. The effects of higher order polarizations, even for particles in close proximity, are appropriately captured in this scheme. It is only when particles collide that charges resulting from induced dipoles are mapped onto the domains to become permanent (described next).
At the instant when two grains collide, neutralization is imposed between contacting charge domains of colliding grains i and j according to Eq. (2). This permits charge transfer between grains, so for η = 0 all charges remain unchanged, and grains increasingly transfer charges as η grows. Equation (2) is applied during binary collisions, and whenever a grain contacts multiple neighbors during a single time step, we perform this operation for all pairs in random order.
Both permanent and induced charges take part in neutralization events, and to keep accounts straight we add exactly the fraction of induced charges needed to conserve charge to the permanent charges. That is, if an induced charge Δ q is added to one domain of a grain due to Eq. (2), then − Δ q will be made permanent on the opposing domain of the same grain to guarantee overall charge neutrality. Finally, we prevent spurious repetition of charging by only applying neutralization and induction operations at the moment when two grains first touch one another.
We enforce energy conservation in two ways. First, we compensate for the energy associated with assembling the induced dipole moment prescribed by Eq. (1) by integrating the work needed to bring the induced charges to their positions from infinity 31 . We then evaluate the gradient of this energy, which gives us a mechanical force that we apply to each particle. This is the force that must be exerted to polarize the particle, and the spatial integral of this gradient is mechanical work that exactly equals the required electrical energy.
Second, we note that an induced dipole moment changes the electric field of neighboring grains, and this change in turn induces secondary dipole moments according to Eq. (1). To conserve energy, we account for secondary dipoles by reducing each primary moment by exactly the energy associated with every secondary moment. This process feeds back iteratively, so that every secondary moment in turn induces another moment on the originating dipole. We have numerically confirmed that this feedback converges rapidly, and after two iterations the error in neglecting higher order terms is less than 0.8%. Consequently, in our simulations we perform two iterations of inducing new additive dipole moments based on this feedback process.
Neutralization events also involve energy considerations: when a particle of charge q i fully neutralizes (η = 1) during collision with a neighbor, each particle will leave the collision with a charge of up to q i /2. This produces repulsion between the particles that was not present prior to the collision. This repulsion is a real physical effect that is seen in experiments 32 , and so we include the repulsion in our simulations.
We initialize each simulation by dropping 10 3 grains onto the fixed bottom, of area 10R × 10R. No energy is injected (through kicks by the bottom surface) while particles settle, and grains are all initially neutral. We wait five seconds until grain velocities become negligibly small (kick velocity/10 3 ). We then add charges uniformly distributed on [− 10 −2 , 10 −2 ] pC to all six domains of all grains and thereafter kick particles at the bottom surface.