Evolutionary strategy for inverse charge measurements of dielectric particles

Electrostatically charged granular particles are important in a wide variety of applications, ranging from particulate matter pollution to industrial handling of pharmaceutical products, food grains, and inks for printing and additive manufacturing, to name a few.^1,2 Granular dielectric particles often acquire charge through tribocharging or contact electrification;^3–10 the charges they carry can significantly affect their dynamics and their interactions with the surrounding environment. In order to better understand how such charged, polarizable particles interact, it is therefore of fundamental importance to make detailed measurements of their actual charge.¹¹ Recently developed experimental techniques have made attempts to determine the charges of individual particles in a vacuum environment using free-fall videography.^12–14 In those experiments, particles falling under the influence of gravity were filmed as they interacted in a vacuum tube. In one experiment, by accelerating charged particles in a horizontal electrical field and analyzing approximately ∼10⁴ trajectories,¹³ the average particle charges were estimated by relying on the relationship between acceleration, mass, and charge. In another experiment, by identifying the relative positions of two particles and fitting Kepler-like orbits¹⁵ to their motion, it was possible to determine their charges. The interactions that arise amongst polarizable particles are inherently many-body, and it is therefore essential that new approaches be developed that are capable of taking such effects into account. In this work, an approach is proposed that is capable of simultaneously measuring the charges of many individual particles from a single set of trajectories (i.e., the trajectories of the particles from a single experiment, as opposed to an ensemble of trajectories from many different experiments). The approach relies on two advances: (1) the availability of new numerical algorithms and new analytical expressions capable of describing polarizability effects on interacting particles^16–25 and (2) the availability of modern evolutionary computation strategies^26–29 that enable direct interpretation of experimental data from numerical computer experiments.

To computationally determine charges on granular particles from a given single set of target trajectories assembled on a time-sequence of N_f frames, we adopt an evolutionary optimization technique that seeks to minimize a fitness function. Here that function f is defined as the deviation between trial trajectories generated in each optimization step and the target set of experimental or computational trajectories,

where N_f excludes the initial configuration, N_p is the number of granular particles, and $ri,trial(k)$ and $ri,target(k)$ are the positions of the ith particle at kth frame in the trial and the target trajectories, respectively. If the masses of the particles are known, the trial trajectories can be obtained using simulations with a suitable force field, i.e., in this case the electrostatic interactions, which include pair-wise Coulombic forces and many-body dielectric polarization contributions. The bare Coulombic interaction can be attractive or repulsive, depending on the sign of the charges. Polarization effects, however, are purely attractive when the internal dielectric permittivity of materials is greater than that of the medium and are purely repulsive in the opposite case.²⁵ This polarization-induced attraction is summarized in Fig. 1 for two and three particles of equal-sign charge. The figure shows the conditions, ϵ_in/ϵ_out and Q₁/Q₂, under which two and three particles will form stable (attractive) aggregates. The boundary between the stable and unstable states is calculated using a recently developed analytical, perturbative theory²⁵ (see Sec. II B below). It is clear that dielectric polarization can strongly influence the nature of interactions between charged dielectric objects, particularly when sharp dielectric discontinuities are involved. In this work, we use the recently proposed analytical formalism (the image method) to calculate electrostatic interactions between polarizable granular dielectric particles, and we have implemented the resulting electrostatic force field into Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS, http://lammps.sandia.gov)³⁰ to simulate trajectories of the particles.

A Covariance Matrix Adaption Evolution Strategy (CMA-ES) is adopted, and we rely on the open source library libcmaes (https://github.com/beniz/libcmaes)³¹ to extract the charges through an iterative optimization process. We address the inverse problem under two scenarios. In the first, the initial velocities of the particles are known, but their charges are not. The search variables are therefore the N_p charges. In the second scenario, both the initial velocities and charges of the particles are unknown, so there are in total 4N_p search variables [N_p charges and 3N_p velocities in 3-dimensional (3D) space]. As our results demonstrate, the proposed strategy is able to determine the charges of the particles under both scenarios. It is difficult to know how many local optima may be encountered in the fitness landscape. In our simulations, we find that the optimizer can be trapped in one of several optima during a series of consecutive optimizations for the test cases in Secs. III B and III C, which suggests that multiple optima are generally accessible. To help the optimizer escape a local optimum, and move toward the global optimum, we restart the optimization and rescale all search variables as shown in Secs. III B and III C. Mathematically, the global optimum is found when the fitness function decays to zero; numerically, the global optimum is found when the fitness function is smaller than a tolerance. The value of the tolerance depends on the underlying errors in the experimental trajectories as well as on the numerical approximations involved in the particle simulations.

The image method is an analytical method capable of describing polarizability effects on interacting particles. We consider N spherical particles, with a radius a and dielectric permittivity ϵ_in, embedded in a continuum with dielectric permittivity ϵ_out. In principle, our approach may be used for polydispersed particles with different ϵ_in and a, but we will limit the discussion in this work to monodispersed systems of equally sized spheres. The ith particle carries a point charge Q_i(x_i) = z_ie at the center, where z_i is the valence and e is the elementary charge, which implies a homogeneous free surface charge density on the particle. Note, however, that in future work, it should also be possible to pursue the proposed inverse calculations using numerical methods—such as those proposed in our recent work—that take inhomogeneous free charge distributions into account.²¹ In this work, we use the numerical method to calculate the induced surface charges on particles and found that the induced surface charges are inhomogeneous and their distribution changes as particles move along the trajectory as shown in Fig. 5.

The main interactions between the particles are the usual pairwise Coulombic interactions. However, when the particles are in close proximity, they induce surface charges that give rise to additional interactions. We recently developed a systematic multiple-scattering formalism^24,25 to describe this polarization interaction. In this formalism, the polarization energy is grouped in terms according to the number of interacting particles. The lowest-order of the polarization energy, i.e., the three-body terms E₃, is contributed to by 3 particles. The higher-order terms E₄ and E₅ involve four-body and five-body interactions; the two-body terms are reserved for the normal pairwise Coulombic interaction. Symbolically, the total electrostatic energy, E_E, for an ensemble of dielectric spheres may then be written as E_E = E₂ + E₃ + E₄ + ⋯, where each term in such an expansion involves a summation over all possible two-body, three-body, four-body, and so on permutations. The key point to note about this multi-body expansion for E_E is that all interaction terms only depend on the particle positions. The references to surface charges are avoided by replacing the induced charges by the gradient of the electrostatic potential, therefore the degrees of freedom are greatly reduced. Furthermore, the forces on the particles can be computed via differentiation with respect to the particle positions, which enables N-body particle simulations. Three terms are preserved for the electrostatic multibody potential in the particle simulations, which have been shown to be essential in describing particle interactions in the presence of polarization.^21,24,25

CMA-ES is one type of evolutionary optimization algorithm³² that does not require derivative information of the fitness function, so it enables minimizing a broad range of fitness functions that have no analytical forms. In general, the idea of evolutionary optimization is to first generate a sample of random search variables every generation following a Gaussian distribution, then select the best search variables that produce the most optimized value of the fitness function. This process is then iterated until the fitness function is within a target convergence criterion. In CMA-ES, the mean and covariance matrix of the search variables as well as the step size are updated every generation to achieve fast and successful optimization. CMA-ES has found many applications in various materials design problems.^33–38 

We consider spherical granular particles in a 3D vacuum environment, where only electrostatic interactions are the driving force for their motions. Recent experiments¹⁵ in a similar setting have observed striking phenomena of aggregation and motion of charged granular particles, from which we adopt the particles’ parameters for our simulations. Specifically, the particles are monodisperse and have relative dielectric constants of 15, diameters of 260 μm, and mass densities of 3800 kg/m³. No thermal fluctuation or Brownian motion of the particles is included; the particles are in a vacuum and have diameters of hundreds of microns. The system may become chaotic when it is evolved over time scales longer than those used in this work. However, a short-time trajectory without chaotic behavior (maximum of 25 ms in this work) is sufficient to successfully extract the charges of particles using our proposed strategy. We also neglect the particles’ rotational motion and only account for their translational motion because particles studied in this work have spherical shapes. The charges on every particle are assumed to be uniformly distributed and the electrostatic interactions that include both the Coulombic interaction and the polarization effect are calculated using the aforementioned image charge method.^23–25 To simulate the trajectories of granular particles, Newton’s equation of motion is integrated by the velocity-Verlet algorithm in LAMMPS with a time step of 1 μs. For the electrostatic interaction, the boundary condition is such that the electrical potential decays to zero at infinity. The particles are simulated in the NVE ensemble without periodic boundary conditions. In test problems, the initial positions of the particles are randomly generated while ensuring that there are no overlaps between any two particles. For the test problems with ten and 30 particles in this work, the target trajectories are generated by simulations with initial velocities all set to zero. For the problem with ten particles, the particles’ charges are ±1, ±2, ±3, ±4, and ± 5 pC, respectively. After t = 0, particles start to move under the influence of the electrostatic forces. For the problems with ten and 30 particles, the dynamic simulations are run for 2000 steps (2 ms), and we sample and store the trajectories with a frequency of 1 frame per 10 steps (10 μs). There are no collisions between any two particles in the trajectories, so tribocharging phenomena are avoided. The trajectories generated by these simulation are then imported as the target trajectories to the optimization program for the inverse calculation of charges on the ten and 30 particles, respectively.

As alluded to earlier, in this work we assume that it is sufficient to assume a uniform charge distribution on the particles’ surface, which is equivalent to placing a point charge in the center of the spherical particle according to Gauss’ law.^24,25 In reality, the free charges on particles’ surface may be distributed inhomogeneously.³⁹ According to Ref. 39, immediately after the contact, the surface charge distributions are inhomogeneous in the contact areas; after 2.2 h, they become uniform. In our experiments, the particles are placed and stabilized in the chamber for more than 3 h before their free fall begins, such that the surface charge distributions are uniform. During free fall, there are collisions and contacts between particles and charges are expected to be inhomogeneous in contact areas. Because the contact area is much smaller than the total surface area of each particle, however, the charges can be reasonably assumed to be mostly uniform on the particles’ surfaces. As a result, the charge non-uniformity has a minor effect on the particles’ trajectories. Ideally, we would include charge non-uniformity in our simulations and estimate charges on individual patches on all particles; this could be done by relying on numerical methods such as those introduced in our earlier work²¹ but at greater computational expense. If there are many patches of charges, we conjecture that the number of solutions for the charges that can match experimental trajectories will still be one (with nonzero external electric field), or two (without external electric field, where the signs of the charges will be opposite in the two solutions). As the number of patches increases, the fitness function landscape becomes rougher, and the number of iterations to converge and the associated computational cost increase. Thus, using an efficient electrostatic solver²¹ would help reduce the computational cost for this problem. However, in view of the lack of systematic experimental data on how the charges are really distributed on the particles during their free fall, we find it difficult to assume a certain pattern of non-uniform surface charge distributions on the particles. By instead assuming a uniform distribution, we are still able to produce simulated trajectories that agree very well with experimental trajectories, serving to validate our assumptions. Note that while the free surface charges are assumed to be homogeneous, the induced surface charges on particles are inhomogeneous. These are calculated using the numerical method of Ref. 21; the distribution of the induced surface charges evolves as particles move as shown in Fig. 5.

We first study the inverse problem when the charges of the particles are the only unknowns; the known particles’ initial positions and velocities are used to start the simulation. In our test case, we first generate trajectories of ten charges with assigned charges and we use a single final frame at t = 2 ms with the final positions of the particles, so N_f in Eq. (1) is 1. The true (assigned) charges of the particles are ±1, ±2, ±3, ±4, and ±5 pC, respectively, and the aim of this first example is to demonstrate that the evolutionary optimization process can correctly recover those charges from the knowledge of the particles’ masses, initial positions, and velocities. We use the following three parameters in the CMA-ES evolutionary optimization. The initial values for charges are set to zero; the initial search step is 2, and the number of offsprings is 10. Figure 2(a) shows the velocities of three representative particles as a function of time throughout the simulated trajectory. One can see that the trajectory generated by the simulation using inversely calculated charges (the dotted line) agrees well with the target trajectory (the solid line). Figures 2(b) and 2(c) show the evolution of the fitness function and the estimated charges of ten particles as a function of the number of fitness function evaluations. The fitness function decreases as the number of evaluations increases, and it converges to about 2 × 10⁻² after 110 generations (∼1100 fitness function evaluations). The deviation between trial and target particle trajectories becomes smaller as the optimization proceeds, and the trial charges on the particles gradually evolve to their target values. When the fitness function reaches a plateau, the particles’ charges stabilize at the correct values of ±1, ±2, ±3, ±4, and ±5 pC, respectively.

To examine how the initial guess for the particle charges affects convergence, we chose four different starting values, i.e., 0 pC, 1 pC, 10 pC, and 100 pC, and examined the evolution of our optimization while keeping all other CMA-ES parameters constant. For an initial guess of 0 pC, 1 pC, and 10 pC, all three optimizations yielded correct estimates of the charges, but the number of generations to reach convergence increased as values of particle charges in the initial guess increased. When the initial guess is 100 pC, the optimizer cannot reach convergence in a single optimization, i.e., the estimated charges at the end of the first optimization are of a different order of magnitude. These findings serve to illustrate that it is critical to choose values for the initial particle’s charges that are of the same order of magnitude as the target value, otherwise, rescaling the search variables multiple times is necessary to achieve convergence, as shown in Sec. III B.

In this test case, we also find that using only the Coulombic interactions (i.e., neglecting polarization interactions) lead to 10% error in the inversely calculated charges. We note that the importance of polarization depends on (i) the ratio between the dielectric constant of the particles and that of their surrounding environment and (ii) the charge ratio between interacting particles. The polarization effect becomes more important when the dielectric ratio or charge ratio increases. In the first problem considered here, the dielectric ratio is 15 and the maximum charge ratio is 5. According to Fig. 1, with this combination of parameters, the polarization effect is not too strong, which is consistent with the 10% error that is observed when polarization is neglected. However, for particles with dielectric ratios that are larger than 15 and charge ratios that are above 5, neglecting polarization effects leads to errors that are much larger than 10%; in those cases, polarization effects must be taken into account. Moreover, attractions and adhesions between like-charged particles were observed in Ref. 15, where polarization was shown to play a central role.

Accurately recording both positions and velocities of granular particles using videography is challenging. In most experiments, only the particles’ positions are recorded. It is of course possible to approximate velocities at every frame using a finite difference approximation, but that may lead to a loss of accuracy. It is therefore of interest to explore the use of only information about the positions of the particles for inverse determination of their charge. We find that, when the particles’ initial velocities are unknown, using multiple frames from the particles’ trajectories to evaluate the fitness function can enable such inverse charge calculation. Specifically, 20 consecutive frames are selected from the initial stage of the simulated trajectories for ten particles (having the same charges as above) at an interval of 10 μs; N_f in Eq. (1) is set to 20. The parameters for performing the CMA-ES evolutionary optimization are as follows: the charges and initial velocities are set to zero for all particles, the initial search step is 2, and the number of offsprings is 10. Figure 3 shows the evolution of the charges and initial velocities as a function of the number of fitness function evaluations. As the optimization progresses, the initial velocities are estimated to be close to their true values, as shown in Fig. 3(b). Unfortunately, however, there are significant deviations between the estimated charges and their true values, even when the charges stabilize at the end of the optimization process (about 200 generations), see Fig. 3(a). The converged value of the fitness function for Figs. 3(a) and 3(b) is 1.28 μm.

The estimated initial velocities at the end of the optimization process in Fig. 3(b) are on the order of 10⁻³–10⁻⁴ m/s, while the charges are on the order of 10⁰ pC. Thus there is a difference of about 3 to 4 orders of magnitude between the numerical values of the estimated charges and the initial velocities. This points to a numerical artifact of the optimization process as originally implemented; for successful optimization using CMA-ES, it is essential that variables be rescaled in order to ensure that all numerical values of the search parameters are of the same order of magnitude.⁴⁰ Then a second optimization simulation is started by setting initial values for the search parameters to those corresponding to the last step of the previous simulation; in the subsequent simulation, we re-scale the search variables for initial velocities by 10⁻⁴, i.e., if one of the initial velocities in the optimizer is 1, then its value fed to the dynamical simulation is 10⁻⁴. The results of this second optimization process are shown in Fig. 3(c). One can appreciate that the charges evolve rapidly toward their true values after about 30 generations (300 fitness function evaluations). At the end of the optimization process, the charges agree very well with their true values and the converged value of the fitness function for Fig. 3(c) is 3.69 × 10⁻⁴ μm, serving to demonstrate that the inverse calculation process can accurately estimate charges from known trajectories, even if the particles’ initial velocities are unknown.

To further test the robustness and applicability of our proposed strategy, a third test was performed on a system of 30 particles with randomly assigned charges, drawn from a Gaussian distribution with zero mean and a standard deviation of 5 pC. Their values are represented by red dots in Fig. 4. True trajectories were then generated from simulations using these random charges and zero initial velocities as inputs. We applied the evolutionary optimization strategy to this problem assuming unknown charges and unknown initial velocities. Multiple optimization simulations were performed consecutively to rescale the search variables properly. As can be seen in Fig. 4, upon convergence, most of the estimated charges (blue dots) do not agree particularly well with their true values (red dots). In fact, for the 13-th to 30-th particles, the signs of the estimated charges are just the opposite of their true values. On the other hand, the absolute values of the estimated charges agree well with the true values. This happens because of the symmetry of the system, i.e., a trajectory remains the same when the particles’ charges all have reverse signs (the mass and shape for all particles are the same).

To resolve the sign problem, an external electric field was applied to break the charge symmetry during the generation of model trajectories. The evolutionary optimization strategy was then used on these new trajectories. Figure 4 shows that after optimization with the applied field, the charges obtained through the optimization process (black dots) agree very well with the target values. This result shows that by breaking the charge symmetry by applying an electric field, the evolutionary optimization strategy can correctly recover the charges of the particles. We have also varied the magnitude of the applied field from 0.1 to 100 V/μm, and found that all values lead to the correct sign of the charges. The difference between various magnitudes of the applied field is that larger electric fields can generate a trajectory that is different from that without a field within a shorter amount of time. Since the only good of applying a field is to determine the sign of a charge (and not its magnitude), such differences have no influence on our results.

Finally, we apply the evolutionary optimization strategy to experimental data to (i) calculate charges on granular particles in experiments and (ii) reproduce the experimental trajectories using simulations. A set of trajectories for three granular particles and another set of trajectories for four granular particles are chosen from the experimental data. The data set, which consists of 25 frames, covers a span of 25 ms. Data were captured from videography; a particle tracking technique was applied to extract the coordinates of all particles. The time interval between consecutive frames is 1 ms. In the trajectory with three particles, two particles are always in contact with each other, while the third particle moves freely around the other two; in the trajectories with four particles, particles are always in contact with their neighbors. A bond is formed between the two sticking particles by short-range cohesive forces, including van der Waals forces or capillary forces due to absorbed molecular layers.^12,41,42 These short range interactions are strong enough to hold the two particles together without relative translational and rotational momentum between each other. To reproduce this behavior, a rigid bond is implemented in the simulations between the sticking particles. The length of the rigid bond is set to the diameter of one granular particle and is maintained in every simulation step using the SHAKE algorithm.⁴³ The charge on each particle is kept constant in particle simulations, because charge transfer between particles is negligible during the short trajectory time of 25 ms.¹⁵ The coordinates of three and four particles extracted from experiments are then fed as the target trajectory into the optimizer, and the optimization strategy is applied. Note that in this case the signs of the particles’ velocities can be inferred from experimental data, and we constrain the sign of the search variables for initial velocities in the optimization program to increase efficiency.

Figures 5(a)–5(d) show the evolution of charges and initial velocities as the optimization proceeds, for three and four particles, respectively. The results presented here are from the last optimization simulation in a set of consecutive calculations. In the first few simulations, the fitness function decreases to a plateau at the end of each simulation. However, the fitness function is still large and on the order of 10². Moreover, there is a difference of about three orders of magnitude between some of the search variables. These consecutive simulations are merely used to properly rescale all search variables; we then feed the final charges and initial velocities to the particle simulation to generate a simulated set of trajectories of three and four particles. Figures 5(e) and 5(f) show six snapshots of three particles from experiments and simulation, respectively. The induced surface charge density (σ_pol) on every particle is calculated using a parallel O(N) numerical solver for electrostatic polarization available in the COPSS suite of codes (miccomcodes.org).²¹ We obtain excellent agreement between the simulated and experimental trajectories of the particles. The charges on the three particles are obtained as 1.785, −1.338, and 2.0 pC, respectively [from left to right in the first snapshot in Fig. 5(f)]. They could also be −1.785, 1.338, and −2.0 pC because of the symmetry of the system. The range of calculated charges is consistent with that inferred in previous experimental results;¹⁵ the experimental charge distributions P(q) have tails up to several million electron charges (10⁶ e ≈ 0.16 pC). It is also found that the two bound particles carry opposite charge, and the Coulombic attraction force helps bind them together. For the set of trajectories comprising four dielectric granular particles, Figs. 5(g) and 5(h) show six snapshots from experiments and simulations, respectively. Excellent agreement is again found between both, serving to demonstrate the applicability of our proposed optimization strategy for inverse charge calculations. The charges on the four particles are obtained as −9.97 × 10⁻³, −9.37 × 10⁻³, −5.12 × 10⁻³, and −1.37 × 10⁻² pC, respectively [from right to left in the first snapshot in Fig. 5(h)]. They could also be of positive sign because of the symmetry of the system. The range of calculated charges is consistent with that inferred from previous experimental measurements.¹⁵ Note that the magnitudes of the charges for four particles are much smaller than those for three particles, and the signs of the four particles’ charges are the same, indicating that polarizability is essential for describing the physics of the particles considered here. The converged values of the fitness function for Figs. 5(a) and 5(c) are 13.66 μm and 48.5 μm, respectively, which are larger than those for Fig. 3. The converged values in Fig. 3 are relatively small because (i) the target trajectory is generated by simulation and the trajectory data are an exact representation of the simulated particles’ motion; (ii) there is no error from the numerical simulation results since the models used in generating the trajectory and in optimizations are identical. The relatively large converged values in Fig. 5 when fitting the experimental data are likely due to errors from experimental measurements and approximations used in the numerical model. The experimental errors are from vibrations of the camera and the particles’ position tracking process; note that the resulting errors in the particles’ positions are on the order of 10 μm. The approximations in the numerical model include using rigid bonds to connect sticking particles and assuming a homogeneous free surface charge density on each particle. Movies for the above sets of trajectories from the experiment and simulation can be downloaded from the supplementary material.

In summary, we have combined the evolutionary optimization strategy CMA-ES with a particle dynamics simulator to obtain the charges on granular polarizable particles based on a given set of experimental trajectories. The availability of a polarizable force field for electrostatically interacting charged granular particles is central to the particle dynamics simulator; electrostatic polarization and Coulombic interactions can in some cases have opposite signs and lead to trajectories that are very different from those observed in the absence of polarizability effects. The proposed strategy was demonstrated in the context of several problems. In the first problem, the initial position and velocities of all particles were given and the algorithms were used to infer the particles’ charges. In the second and third problems, both the particles’ charges and initial velocities were unknown, and it was shown that the evolutionary optimization can be used to successfully determine the particles’ charges and their initial velocities. In the fourth problem, the evolutionary optimization strategy was applied to extract the charges from experimentally observed trajectories and the charges were found to be within the ranges reported in previous experimental measurements from the literature.

The proposed strategy could be extended to more complex systems containing electrostatically charged granular particles. For example, using a recently developed fast parallel O(N) numerical solver for electrostatic polarization interactions among arbitrary-shaped particles,²¹ which is made available under the COPSS suite of codes (miccomcodes.org), the evolutionary optimization strategy could be applied to determine charges not only of spherical particles but also of arbitrarily shaped particles with uniform or nonuniform surface charge distributions, including rotational motion. The proposed strategy could also be used to determine the charges of particles in micro- or nano-fluid environments by coupling the strategy with a recently developed parallel O(N) Stokes’ solver for hydrodynamically interacting objects in general geometries which is also made available under the COPSS code.⁴⁴ We envision that our proposed strategy could find applications in material property measurements and material designs.

See supplementary material for the experimental and simulated trajectories of dielectric granular particles.

The work was supported by MICCoM, as part of the Computational Materials Sciences Program funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division. The generation and analysis of experimental trajectories of interacting colloidal particles was supported by the NSF MRSEC program under No. DMR-1420709. V.L. acknowledges support from NSF No. DMR-1309611. We gratefully acknowledge the computing resources provided on Blues and Bebop, a high-performance computing cluster operated by the Laboratory Computing Resource Center at Argonne National Laboratory and the University of Chicago Research Computing Center for support of this work.