Particle-Based Simulations of Electrophoretic Deposition with Adaptive Physics Models

This work represents an extension of mesoscale particle-based modeling of electrophoretic deposition (EPD), which has relied exclusively on pairwise interparticle interactions described by Derjaguin-Landau-Verwey-Overbeek (DLVO) theory. With this standard treatment, particles continuously move and interact via excluded volume and electrostatic pair potentials under the influence of external fields throughout the EPD process. The physics imposed by DLVO theory may not be appropriate to describe all systems, considering the vast material, operational, and application space available to EPD. As such, we present three modifications to standard particle-based models, each rooted in the ability to dynamically change interparticle interactions as simulated deposition progresses. This approach allows simulations to capture charge transfer and/or irreversible adsorption based on tunable parameters. We evaluate and compare simulated deposits formed under new physical assumptions, demonstrating the range of systems that these adaptive physics models may capture.


Introduction
Electrophoretic deposition (EPD) is a well-studied process where colloidal particles, driven by an electric field, accumulate at an electrode to form a deposit. The accessible property space of these resulting deposits is enormous due to the large number of independently tunable parameters governing EPD, such as the electrode substrate, electric field strength and waveform, particle morphology and loading, supernatant characteristics, and so on. As such, EPD is found in a wide variety of applications, ranging from neural electrodes, 1,2 biomaterials processing, [3][4][5] high-contrast low-voltage 6,7 and quantum dot-based displays, 8 battery 9 and supercapacitor 10 manufacturing, fog harvesting, 11 superhydrophobic 12 and energetic 13 coatings, etc. Though efforts to optimize EPD deposits must be application-specific, the community typically relies on a general set of long-standing models for the stability of colloids and their deposition. For instance, particles are assumed to interact strictly according to Derjaguin-Landau-Verwey-Overbeek (DLVO) theory and adhere to a deposition rate that scales linearly with the strength of a uniformly applied field. 14 Fortunately, numerous models exist to accommodate the nuances of actual EPD systems.
Ferrarri, et. al.'s extensive review 15 lays out several analytical expressions that account for commonly encountered EPD scenarios, including flat 14 and cylindrical electrodes 16 with variable 17 and/or highly concentrated 18 particle loading. Empirically validated models exist that account for suspension resistivity variations alone 19 and in combination with solid loading. 20 Several models rely on finite element analysis (FEA) to simulate coating at the industrial-scale 21 and deposit thickness variations that result from electrode edge effects. 22 A recent multi-physics FEA model accounts for transient phenomena arising in a specialized EPD cell with cross flow across empirical time and length scales. It reproduces actual deposit morphologies over a wide range of inlet flow rates and applied voltages. 23 Finite element based approaches offer unique insights, specifically into how the deposition thickness evolves under (possibly) complicated electrode geometries and electric, flow, and particle concentration fields.
Particle-based models (PBMs) elucidate how deposit microstructure depends on the interplay between inter-particle and electric field forces during deposition. Specifically, PBMs compute the trajectory of all colloidal particles in a simulated control volume using Newton's equations of motion and assumed interaction potentials that model particle-particle interactions in a suspension bounded by a charged electrode. As such, PBMs reveal how particles accumulate at the electrode, though typically over shorter time scales than FEA-based simulations due to greater computational expense. PBM simulations accurately reproduce analytical models of bulk electrophoretic motion across an expansive experimentally relevant regime 24 and have elucidated how the degree of ordering within multi-layed deposits depends on both ionic concentration and field strength. 25 The same PBM applied to sub-monolayer deposits revealed correlations between electrochemical properties and composition of the supernatant with coating homogeneity. 2 PBM particle-spacing predictions were experimentally validated using a novel EPD cell that enables in situ characterization. Considering long-standing critiques of DLVO 26 applied to the wide range of available colloidal suspensions, it is fortunate that the underlying inter-particle potential computed by PBMs is fully customizable.
To date, PBMs exclusively assumed that interparticle interactions are of a fixed functional form that follows DLVO theory. Figure 1 shows representative simulation snapshots from a particle-based EPD simulation. Panel (a) shows the starting configuration, which consists of randomly distributed colloid particles. When the external electric field is applied (panel (b)), the particles begin to migrate in the −z direction, toward the electrode, assembling into a deposit. Panel (c) shows the representative simulation at equilibrium. We note that the position of each particle is stored and represented as a point, not a space-filling sphere, in these simulations and the data shown in this work. This representation allows us to clearly observe surface-induced ordering of the deposits onto the geometrically flat electrode in the snapshots by illustrating each particle as a gray point whose radius is smaller than the particle's physical radius. The apparent 'space' between the first deposition layer and the electrode (and the first layer and the bulk deposit) emerges from and is enhanced by this representation. Gray rectangles are included in this overview figure to indicate the position of the electrode. Subsequent simulation snapshots show only the particles. The density profile in panel (d) corresponds to the equilibrium configuration in (c) and quantifies the surface-induced ordering. We again note that this data captures the center points of the particles. Density profiles specifically report the density of the particles' centers, not a bulk density, which enhances our ability to observe ordering of particles in the deposit.
In this conventional PBM approach, interparticle interactions and colloidal motions described by an unchanging physics model persist at all times during deposition. Given the diversity of colloidal suspension and deposition conditions, we argue that many EPD systems of interest exhibit behavior not described by a static set of pairwise interactions defined a priori. Thus, we pose a PBM framework that allows for exploration of particles with physics models that may change on-the-fly according to a set of pre-defined heuristics. With these new models, we simulated parameter sweeps that vary the applied electric field and characteristic interparticle electrostatic repulsion. Specifically we simulate conditions that can occur in real systems in which particles can discharge and/or become immobile upon deposition e.g. when a metal nanoparticle deposits on a metal surface and forms a covalent or metallic bond at the interface. It is beyond the scope of this work to explore all parameter space made accessible by these new simulation protocols. Nevertheless, we simulate experimentally-accessible EPD conditions and compare particle packing density and degree of ordering resulting from the four model types. Further, the physical parameters used in these simulations may be adjusted to explore any EPD system of interest. Figure 1: Simulation snapshots show typical particle configurations for a conventional DLVO particle-based EPD simulation at (a) the start of the simulation, (b) as deposition begins, and (c) at equilibrium. (d) Equilibrium density profiles (with vertical axis aligned with the simulation z-coordinate) reveal a maximum density at the electrode followed by a few layers of decaying order and a region of deposited material whose density decreases linearly toward the supernatant region devoid of particles.

Model and simulation details Interaction potentials
Mesoscale particle dynamics simulations are based on the work of Giera and collaborators. 25 We refer the reader to this work for a full description of interaction potentials and model parameterization, but summarize the components most salient to the present work here.
The total energy acting on a colloid particle is a combination of the pairwise sum of solvent-colloid, colloid-colloid, and system-colloid interactions, where r is the distance between an i-j colloid pair and z is a colloid's distance from the electrode, here represented by a flat, uncharged plane in the x-y plane. The solvent is treated implicitly by adding pairwise dissipative lubrication and Brownian forces. [27][28][29][30][31] Pairwise colloidal interactions are treated with the DLVO potential, U DLVO (r) = U Yuk (r) + U steric (r) + U vdW-a (r) (2) which represents the sum of screened electrostatic repulsions by a variant of the Yukawa potential, 32 where is the relative permittivity of the solvent, 0 is the permittivity of free space, k B is Boltzmann's constant, and T is temperature. We note that this definition of λ D holds for colloidal suspensions in a symmetrical electrolyte with |q ± | = q, e is the elementary charge, and ρ ± is the bulk concentration of ions. As Eq. 3 indicates, the electrostatic repulsion between colloids is easy to tune empirically, especially in aqueous suspensions, e.g. λ D increases as ρ ± decreases via dilution, vice versa. For our particle-based simulations, the computational expense scales with the number of pairwise energies with U DLVO (r) ∼ k B T , which increases with λ D .
Interactions of the colloid with the electrode, U wall , and the electric field, U e-field are single body potentials, i.e., a function of the particles' z-position. The colloid-wall interaction prevents the particles from moving through the simulation cell wall in the −z-direction, resulting in simulated deposition. Electrophoretic motion is simulated by application of a uniform, external electric field where E field is the strength of the electric field and q col is the effective charge on a colloid particle in the direction perpendicular to the electrode, obtained from a force balance on a colloid particle experiencing steady-state hydrodynamics drag in an electric field. 36

Alternative deposition models
As described earlier, we introduce two new physical scenarios that may be employed individually or simultaneously to more accurately simulate a wide range of EPD applications.
In the first scenario, immobile, we consider systems where particle mobility decreases after deposition and model the extreme case of completely immobile deposited particles, e.g.
when practically irreversible chemical bonds between particle and wall or particle and particle form upon deposition. Our implementation consists of two parts. The first tracks each particle-electrode distance r i (equivalent to z i in these simulations). When a given particle is within a cutoff distance (r sticky ) of the electrode, the particle is flagged as deposited and prevented from moving by removing it from the molecular dynamics integrator. Though this renders the particle immobile, the particle continues to be accounted for in all interaction potential calculations; the particle is only rendered immobile. The second component of this implementation accounts for similar deposition beyond the first colloid layer, necessary since all deposited particles beyond the first layer would have z i > r sticky . To handle this deposition, we track the separation distances of all deposited-suspended i-j pairs during the simulation. If any mobile colloid particle is within r infect of a deposited particle, the mobile particle is flagged as deposited and similarly removed from the integrator. Another physically-relevant case not captured by existing PBMs is when a particle's charge changes upon deposition due to being in intimate contact with the electrode, e.g. Figure 2: Cartoon schematics of the three schemes used to simulate alternative deposition models. (a) In "sticky" a suspended particle is flagged as deposited when in close proximity to the electrode. (b) "Infect" permits the propagation of the 'deposited' state to an arbitrary distance away from the electrode. (c) "Cure" flags particles as resuspended when they separate from the deposit when a conductive contact between a metal particle and a metal electrode is established.
This case includes particles deposited directly onto the electrode and particles that come into contact with the accumulating deposit, which is in contact with the electrode. Our proofof-concept simulations capture this by making deposited particles assume the electrostatic interaction behavior of the electrode, in practice by zeroing the Yukawa potential between deposited and non-deposited particles. We refer to this scenario as decharge. This change to the newly-deposited particle's force field is made when the particle-electrode distance r i is less than a designated cutoff distance, similar in implementation to the particle immobilization detailed above. Also similar to the immobilization scheme, particles tagged as deposited can 'infect' other particles that move within a cutoff distance r infect of an altered particle. These other particles are then tagged as deposited and their force fields are modified accordingly. Source code for these schemes, named fix sticky, fix infect, and fix cure, is included in the Supporting Information.
In practice, combinations of the three implementations shown in Figure 2, sticky, infect, and cure, are applied to modify the physics models of colloid particles on-the-fly, enabling the immobile and/or decharge scenarios. Table 1 summarizes the changes to 'deposited' colloids when using any permutation of these alternative deposition models. For the remainder of this work, we use the nomenclature standard /decharge and mobile/immobile to indicate the off /on states of the two scenarios.

Simulation parameters
The proof-of-principle simulations presented in this work model a system of 10 nm diameter spherical platinum particles in an aqueous suspension, as in on our earlier studies of neural electrode coatings. 1,39 Table 2 summarizes the physical parameters used to specify these simulations. We note that two parameters of particular interest to EPD researchers, Debye

Results and discussion
Here we present an overview of the four model scenarios described in Table 1. Each scenario was run with three Debye lengths and under three electric fields for a total of 36 EPD simulations. This simulation set provides insights into various non-ideal deposition behaviour that can arise in real EPD systems, but are not accounted for in the standard DLVO particlebased model, termed here as standard/mobile. We present a quantitative overview of the four cases, discussing pertinent empirical systems that may exhibit such deposition phenomena.
We then walk through our analysis of density profiles and first-layer particle ordering and present quantitative comparisons between the EPD model types across our simulation set. additive molecular dynamics force field is too limited a framework to capture many EPD approaches. We present three additional EPD simulation scenarios with physics models that adapt to particles' local environments on-the-fly in an effort to explore new physical effects beyond conventional particle-based simulation methods.  charge density and zeta potential ζ. 42 In this example the decharge scenario would be applicable even though there was no externally applied electric field. Figure 3(d), characterized by the standard /immobile scenario, can occur when a gelation agent immobilizes particles at the electrode. Application drivers for this design include sensing, energy conversion, and catalysis. 43,44 Additionally, studying the differences between decharge/immobile and standard /immobile scenarios may provide insights into how the nature of particle-electrode bonding affects the emergent deposit.

Quantifying simulated EPD
Our quantitative analysis of the 36 EPD simulation set focuses on two main features: the density profile and the spatial arrangement of the first layer of colloid particles at the planar electrode. Figure 4 shows representative data for the case of standard /mobile, λ D =5 nm, and E z =250 kV/m.   N 1 , ..., N p , each particle N i has a corresponding Voronoi cell, V i , which contains every point in the plane whose euclidean distance to N i is less than the distance to any other particle at the electrode. Each Voronoi cell represents one particle.
In this work, Voronoi diagrams representing the first deposition layer are calculated and visualized using an in-house customization of the freud Python library, where the color of each Voronoi cell corresponds to its area. 45,46 The inset Figure 4(c) is a histogram of the Voronoi cell sizes, which represents the electrode surface area occupied by the deposited colloid particles. Voronoi analysis has previously been applied to analyze scanning electron microscopy images of EPD monolayer films where the particles were characterized by the number of sides on its Voronoi cell, corresponding to the number of nearest neighbor particles. 47 Voronoi analysis presents a straightforward way survey surface coverage by EPD and to directly and quantitatively compare particle-based simulations and experiment.

Behavior of deposition models
The density profiles in Figure 5 provide direct, quantitative insight into the simulated deposits. In (a), the standard /mobile simulation, three layers of decreasing ordering are visible before the particle density decays to a smooth, monotonically decreasing profile, ending at ∼   presented in earlier work. 25 Implementation of the decharge model, Figure 6(b), results in a more densely packed layer of particles at the electrode, agreeing with data from the density profiles in Figure 5. Histograms of these Voronoi cell areas are overlaid and shown in Figure   7, where the standard /mobile distribution is centered at 150 nm 2 . The decharge/mobile particles are more closely packed, with a considerably smaller variance and an average area of 108 nm 2 . When projecting the first deposition layer onto a plane, this corresponds to 73% of the particles' area covering the electrode surface versus 51% in the standard /mobile case and 90.7% for the case of a hexagonal lattice: the ideal, highest-density configuration for packing circles in a plane. 48 The immobile cases, shown in Figure 7(c) and (d), show significantly less dense packing at the electrode. Although the electrostatics of the deposits formed in the decharge/immobile (c) and standard /immobile (d) cases differ substantially, the average Voronoi cell areas are within ∼4%. This indicates that the separation and ordering of the particles at the electrode in these simulations is a strong function of particle deposition parameters, r sticky and r infect , and not dominated by electrostatic repulsion as in the mobile cases. This similar behavior at the electrode is due to the suspended particles descending with identical physics until they interact with the electrode, then becoming 'deposited' and made immobile. In the first layer changes in the systems' electrostatics (i.e. U Yuk ) are relatively small. As the deposited layers grow, differences in the deposits' electrostatics increase which results in the differences observed in the immobile density profiles in Figure 5(c) and (d). The standard simulation, which retains DLVO behavior throughout the deposit, results in additional suspension-deposit particle-particle repulsion and a slightly more well-ordered and   Figure 6. Histograms are normalized to unit total area. The inset image magnifies the standard/immobile and decharge/immobile histograms and shows the bars as semitransparent to reveal differences between the two results.
Varying electric field for fixed Debye length of 5 nm particles adjacent to the particles tagged as 'deposited,' as seen in Figure 3(b). This layer of particles, electrostatically perturbed by contact with the electrode but suspended in solution, appears to act as a passivation layer, particularly at lower fields, preventing suspended particles in the bulk from reaching the electrode (or deposited particles in intimate contact with the electrode) and losing their electrostatic character. This effect may be made more evident by removing the 'resuspension' scheme (see Figure 2(c)). With fix cure disabled, deposited particles that diffuse away from the electrode move through the bulk suspended colloid, removing the electrostatic repulsion from all suspended particles they encounter and all decharge/mobile simulations result in a dense-packed deposit. Using the terminology of our LAMMPS implementation, fix cure is a first-order approximation to prevent this nonphysical runaway electrostatic deactivation by fix infect.
The immobile cases in panels (c) and (d) appear similar to each other when compared to the mobile cases. Upon close inspection decharge/immobile and standard /immobile density profiles show differences that result from the deposits' dissimilar electrostatic behavior.
Since particles are made immobile upon deposition these cases isolate the effect of depositsuspension repulsion. This repulsion affects the ordering of the deposits in the standard case, resulting in slightly higher density packing near the electrode. In the decharge case the lack of deposit-suspension interaction permits a stochastic deposition; this lack of ordering results in the observed lower densities. During immobile in silico deposition, the fraction of particles tagged as 'deposited,' varies more dramatically in simulations where depositsuspension electrostatic repulsion becomes stronger than the applied field, E z < 250k V/m.
The density profiles in (c) and (d) corresponding to E z = 50k V/m (orange curves) show more interesting differences. In (c), the decharge case, the rapid density decrease between 400 and 500 nm indicates the upper boundary of the deposit. A similar decrease is not observed in (d), although a slight decrease at around 250-300 nm occurs before the density plateaus. The inset images in panels (c) and (d), simulation snapshots corresponding to E z = 50k V/m, reveal the origin of these density features: deposit-suspension repulsion in (d) limits the extent of deposition, even in our immobilize protocol, which implements the most liberal case of immobile particle deposition.   a high-density first layer is present, followed by region of near-zero density that extends to ∼ 180 nm from the electrode surface. A small density peak is present after this excluded region, followed by a low-density colloid suspension.
All three density profiles in panel (b) look identical because the decharge/mobile case removes interparticle electrostatics, thus removing the λ D differences introduced in this parameter study. Deposited particles in this case effectively behave as λ D = 0. This explains the similarity between the λ D = 0.5 nm curve in panel (a) and the curves in (b) since the range of electrostatic repulsion in the λ D simulations is significantly smaller than the particles' radii.
The shorter λ D immobile cases in (c) and (d) look similar aside from the slightly more dense deposit in the standard /immobile case due to additional ordering of the descending particles by particle-particle repulsion. Conversely, the deposits resulting from the λ D = 50 nm cases look very different. In the decharge case, the extended electrostatics in the suspension results in more ordering during the deposition and a more dense deposit than the λ D = 0.5 and 5 nm cases, a trend that suggests a design strategy for higher-density EPD deposits within this physical regime. In panel (d), the standard /immobile case, the λ D = 50 nm deposit resembles the λ D = 50 nm density profile in (a), where the first layer of deposit repels the suspended bulk colloid. In the immobile case, this results in a suspended colloid regime separated from the electrode at that same distance as in the analogous mobile simulation. Figure 10: Deposit density profiles for the four simulation scenarios. Each panel shows the resulting deposits for the given scenario and each colored curve corresponds to a different value of λ D : 0.5 (blue), 5.0 (orange), and 50 (green) nm. E z = 250 kV/m for all simulations.

Conclusions
The particle-based modeling approaches presented here represent extensions beyond standard DLVO-based implementations for EPD simulation. As such, they pave the way for exploring the interplay between electrode and colloidal suspension properties and EPD operating conditions that may result in non-DLVO phenomena. Specifically, we model four types of deposition in which particles (1) keep or (2) lose their charge upon deposition and/or (3) remain mobile or (4) lose mobility upon deposition. We simulate parameter sweeps of these scenarios and compare resulting particle density profiles throughout the deposit and ordering at the electrode. Compared to the standard/mobile case, decharge/mobile deposits pack more densely in the absence of electrostatic repulsion between deposited and (re)suspended particles. This reveals the impact of Debye length on the structure and density of EPD deposits: If λ D is too large, multilayer deposits may not form due to repulsion from deposited particles. Deposit density increases with electric field but is, in all cases, less dense than deposits formed in 'mobile' particle scenarios.
This extensible and versatile EPD modeling framework is an initial step toward expanding particle-based EPD simulation to capture many more physical systems. Modifications to any of the four EPD model types shown here are straightforward to implement. For instance, as with previous particle-based models, the functional form of the interaction potentials among particles and electric field are tunable and, with this work, dynamic. Furthermore, the cutoff distances used to specify whether particles are deposited or suspended can be modified accordingly to match or approximate the behaviour of any physical system of interest. Simulations can be configured to account for multiple deposition scenarios, e.g.
DLVO particles in which a fraction irreversibly adsorb to the surface. Since this framework allows for conditional changes to any property, it is possible to simulate colloids that partially decharge or exhibit reduced mobility near the electrode or deposit.
Though our proof-of-concept simulation set focuses on the extremes of EPD physical phenomena, it is possible to interpolate insights pertaining to more complicated, non-ideal systems commonly encountered in the EPD research community. In cases where multiple deposition effects may be at play and/or for EPD conditions not represented here, the model readily can be tailored to simulate a wide set of experimental systems. Thus, it offers new ways of guiding the experimental design of particle deposition beyond what was possible with previous modeling capabilities. Overall, the model provides an unprecedented quantitative glimpse into the EPD process for both measurable and (presently) empirically inaccessible EPD phenomena, e.g. particle migration within deposits. Considering the various falsifiable aspects of the model, e.g. bulk particle trajectories, surface characterization, etc., model validation and complementary model improvement strategies can be achieved through a tight coupling between experiment and theory-based research.