Phase Transitions of Oppositely Charged Colloidal Particles Driven by Alternating Current Electric Field

We study systems containing oppositely charged colloidal particles under applied alternating current electric fields (AC fields) using overdamped Langevin dynamics simulations in three dimensions. We obtain jammed bands perpendicular to the field direction under intermediate frequencies and lanes parallel with the field under low frequencies. These structures also depend upon the particle charges. The pathway for generating jammed bands follows a stepwise mechanism, and intermediate bands are observed during lane formation in some systems. We investigate the component of the pressure tensors in the direction parallel to the field and observe that the jammed to lane transition occurs at a critical value for this pressure. We also find that the stable steady states appear to satisfy the principle of maximum entropy production. Our results may help to improve the understand of the underlying mechanisms for these types of dynamic phase transitions and the subsequent cooperative assemblies of colloidal particles under such non-equilibrium conditions.


Charged Colloidal Particle Model
We adopt a smooth repulsive Mie potential developed by Jover and co-workers, 1 to reflect the steric repulsion between colloidal particles, which aided our molecular dynamics (MD) simulations before. 2 The Mie potential is given by where λ r = 50 and λ a = 49. The potential is truncated and shifted at its minimum, which results in the purely repulsive pseudo hard sphere potential:  (2) in which ϵ is the interaction parameter, set at 1 k B T , and k B is Boltzmann constant. The quantity, σ, represents the effective particle diameter, and is set to 150 nm in real unit.
The electrostatic interaction between the colloidal particles is modeled using the screened Coulomb, or Yukawa potential, 3,4 where q i , q j are the charges of particles i and j, respectively.
Bjerrum length. This is set to 10 nm, which corresponds to the dielectric constant ϵ r ≈ 5.6 and mimics an organic solvent. We also have

Overdamped Langevin Dynamics
The overdamped Langevin dynamics method is implemented in LAMMPS, which neglects inertia of particles. The equation of motion for particle i is expressed as, where U represents the total particle interaction including Mie and Yukawa potentials, F ext = Eq = E 0 sin(2πωt) is the time-dependent external force by the electric field, γ i = 3πησ is the friction coefficient, according to Stokes' Law and η is the solvent viscosity, which we choose to be approximately 1.17 cP, ξ i is a Gaussian random noise with zero mean and unit variance. The Brownian time τ = σ 2 /D is used as time unit, where D = k B T /γ is the bare diffusion coefficient of a Brownian particle, the time-step for integrations ∆t = 3.73 × 10 −6 τ .

Heat Dissipation and Rate of Work
The averaged heat dissipation rate is described as J = ⟨ṙ i · (γṙ i − ξ i )⟩, it is also shown in the main article. We useṙ i in Supporting equation 4, to substituteṙ i in the equation, and neglect the change rate of averaged internal energy dU/dt = ⟨ṙ i · F c,i ⟩, as dU/dt is much smaller than the rate of work ⟨ẇ⟩ and dissipation J in our systems. Thus we obtain the form for J as 5 Theṙ i in Supporting equation 5 is replaced by its expression in Supporting equation 4 again, after that we calculate the heat dissipation over an oscillating period τ 0 , then we obtain the first term on the right-hand side corresponds to the contribution from external force, or the work of external force on ideal particle, the second term means the contribution from particle interactions. Then we substitute F ext,i with the expression F ext,i = F 0 sin(2πωt), and get a simplified expression for ⟨J ⟩ which is also shown as equation 4 in the main article. The expression is similar as the heat dissipation rate calculated by Tociu and co-workers, 5 but as we adopt the sine-wave external force versus time, the first term contributed by the external forces is slightly different.

System Size Effect
We check the influence of size effect on the phase transition pathways for the two systems shown in Figure 2, by elongating the simulation box along field direction for generating jammed band, and increasing the box width perpendicular to field for lane formation. Both of the two larger boxes are increased by a factor λ = 2 in z direction, and x, y directions, respectively. The similar pathways could also be observed, which are shown in Figure S4.

Density Profiles
We plot the number density profiles of positively charged particles along the electric field of the two systems in Figure 2, after the systems reach steady state (t ≥ 2 s). We choose the density profiles at 4 moments, when the field strength is E 0 , 0, −E 0 and 0 in different periods, corresponding the simulation times are (n + π/2)τ 0 , (n + π)τ 0 , (n + 3π/2)τ 0 and nτ 0 , respectively. The results are displayed in Figure S5. The local oscillations of jammed particles could also be observed via the different profiles of densities ( Figure S5a), and the lanes formed in ±20e charged particles under 10 Hz are not so perfect, according to the density fluctuations in Figure S5b. We also check the influence of system size effect on the density profiles, and find that the elongating simulation system under 100 Hz shows similar density profile ( Figure S5c) as Figure S5a, and the density profile seems to be flatter in the broadening system under 10 Hz ( Figure S5d). Then we calculate the densities at the end of each period in x-y dimension perpendicular to the field for the laned structure. The contour plot of the x-y densities is shown in Figure S5e, and the result for the broadening system under 10 Hz is shown in Figure S5f. The x-y densities show a flattened square profile in the smaller simulation system, and band profile in the broadening system. The system size effect has no influence on the jamming band formation, but influences the lane formation a little.  Figure S4a and b, respectively. (e) Density profile perpendicular to field (x-y) under 10 Hz. (f) Result from the broadening system of (e). Figure S6: The pressure profiles p(z) along AC field direction under ω = 10 Hz at different times, which correspond to the two minima and maxima of p z within an oscillating period.

Work and Heat Dissipation
The rates of work on the time steps dumped during simulations ⟨ẇ(t)⟩ = ⟨∇U i (t) · F ext,i (t)⟩ /γ i are shown in Figure S7, in order to provide more comprehensive view for the time evolutions of rates of work. Figure S7: The rate of work ⟨ẇ⟩ versus simulation times for ±20e charged particle systems.   In order to investigate the influence of collision strength between interfaces on the rate of work ⟨ẇ⟩, we simulate a system with extremely weak collision. A phase-separated configuration like a lane-structure is constructed initially (the left snapshot in Figure S9a). Then we apply a harmonic restraint on x direction, so the colloids will fluctuate in a very small region on x direction, and the initial laned configuration will not be changed too much (snapshot in the lower-right corner of Figure S9a). The ⟨ẇ⟩ values are much lower than the original system without constraint, since the collision between the two lanes is much weaker.
The systems with larger box width (x and y) by a factor λ = 2 are also investigated with the similar approach as above. And the results are shown in Figure S9b. Similarly, the system with constraints shows smaller ⟨ẇ⟩ values, because of very weak collisions. All the two ⟨ẇ⟩ values are smaller than those shown in Figure S9a. So the ⟨ẇ⟩ for lane formation, is smaller as the system size is larger, as ⟨ẇ⟩ is a per particle quantity.