Arrested coarsening and large density fluctuations in driven particle mixtures in two dimensions

In both simple fluids and in more complex, colloidal and polymer systems, excluded volume effects play an important role. A general model describing such systems is the nonadditive binary hard sphere mixture, in which spheres A and B with radii a_A and a_B , respectively, have a mutual collision diameter $(a_A+a_B) (1+\Delta)$, where $\Delta\gt0$ characterizes the nonadditivity of the mixture. This model describes, among many other examples, polymer-colloid interactions (the Asakura–Oosawa mixture [1–3]), phase separation in fluid mixtures at high pressures [4], and liquid–vapor transitions (the Widom–Rowlinson model [5]). Here, we consider a symmetric mixture of two particle types A and B in which $a_A = a_B \equiv a$, with $\Delta = 0.4$. Such a situation can be realized, for example, in mixtures of nanoparticles decorated by polymer brushes. These particles have applications to self-assembled materials, and generally have repulsive entropic interactions [6]. Various polymers may be grafted onto a nanoparticle surface, including immiscible ones such as polystyrene and poly methyl methacrylate (PMMA) [7]. Mixtures of this type are known to phase separate, with kinetics consistent with the Ising universality class [8]. Additionally, driven binary mixtures have also been shown to exhibit interesting phenomena including giant diffusion enhancement [9], segregation [10], enhanced density fluctuations [11], and lane formation [12], to name a few. Here we aim to establish some general features of the driven, nonadditive mixture immersed in a 'solvent' consisting of repulsive point particles.

Although the equilibrium phases of such mixtures are reasonably well-understood, not as much is known about mixtures when driven out of equilibrium. Even single-component driven systems, despite decades of study, continue to present surprises and challenges for both theory and experiment [13–15]. The Katz–Lebowitz–Spohn (KLS) model [16, 17], a lattice model of a single species gas with attractive interactions (exhibiting a liquid–gas transition), is perhaps the best-studied model with a driven, non-equilibrium steady state. In this case, an applied field orients dense liquid domains at low temperature along the field direction. More recent studies [18, 19] demonstrated that a lattice gas model of a binary mixture of particles with repulsive interactions exhibits striped order perpendicular to the drive direction. These models and associated classes of driven systems largely lack experimental realizations, although it is clear that driven colloidal systems are promising avenues for testing these theories [20].

In this work, we attempt to bridge the gap between the theoretical, lattice-model-based results on the driven, nonadditive binary mixtures and experimental observations by studying more realistic, molecular dynamics (MD) simulations which may, in principle, be practically realized. Specifically, we consider MD simulations of a (nonadditive) mixture of two colloids A and B immersed in a collection of solvent point particles S with repulsive interactions. All particles are confined to move in two dimensions, which could be realized experimentally by depositing the colloids on an aqueous interface. We apply a uniform force to the colloidal particles and achieve a steady velocity, mediated by an opposing friction force provided by the S particles. Such a drive may be created, for example, by an applied flow, an electric field for charged colloidal particles, or a magnetic field [21]. Selective snapshots of particle configurations are shown in figure 1 for low ($a,b$) and high ($c,d$) colloidal particle densities ρ (given as the area fraction of the simulation box occupied by the colloidal particles). At large particle densities, our mixture phase separates in the absence of a force. Similar phase separation has been observed experimentally for colloids of different sizes [22].

Zoom In Zoom Out Reset image size

As discussed in more detail below, we find that driving the mixture of particles has two major effects: first, the phase separation between colloidal species is arrested at a characteristic length scale that decreases with an increasing applied field magnitude. Second, the colloidal system develops large density fluctuations in the form of dense 'lanes' or anisotropic clusters separated by lower-density regions. This density distribution is reminiscent of the particle clustering observed in the driven single-component lattice gas (the KLS model). The paper is organized as follows: in the next section, we present the details of MD simulations. In section 3, we explore the structure of the domains of different colloidal species, focusing on what happens to the coarsening behavior in the presence of the applied force. In section 4, we consider the overall particle density distribution changes under an applied field. We discuss these results and conclude in section 5.

We consider two colloid particle types A and B. We apply a uniform force F to these particles and introduce a compensating friction force $- \gamma M \mathbf{v}$ to each particle (with v the particle velocity, M the particle mass, and γ the friction coefficient). This allows the system to establish a non-equilibrium steady state. A uniform force alone, however, would simply introduce a uniform velocity $\mathbf{v} = \mathbf{F}/\gamma M$ to all the particles. Such a uniform velocity would put the system in a moving reference frame, making it otherwise indistinguishable from the undriven case $(\mathbf{F} = 0)$. Thus, to break the Galilean invariance, we introduce a third particle type: a solvent point particle S with mass m which does not experience the applied force F. These solvent particles establish a preferred 'rest frame' and provide additional friction to A and B particles via the repulsive interactions. Therefore, our MD simulations contain three particle types: $\alpha = A, B, S$.

All three kinds of particles interact via truncated Lennard-Jones potentials, approximating hard-core, repulsive interactions. The colloidal particles A and B are spheres, for which the associated potentials have been analytically determined previously [23, 24]. The solvent S is modeled as a collection of point particles with purely repulsive interactions. All particle interactions between particle types α and β (a center-to-center distance r apart) are captured by potentials $W_{\alpha \beta}(r)$ which, in units of $k_{\mathrm{B}}T$, read

$$\begin{align} W_{\alpha \beta }(r) &= - \frac{N_{\alpha \beta }}{6} \sum_{i=0}^1 \Bigg\{\frac{2 a_{\alpha} a_{\beta}}{r^{\,2}-(\Gamma_{\alpha \beta}^{i0})^{\,2}} +(-1)^i \ln [r^{\,2}-(\Gamma_{\alpha \beta}^{i0})^{\,2}]\nonumber \\ &\quad -\sum_{j=0}^1 \frac{(-1)^{i+j} \sigma^6}{6300\,r} \Bigg[ \frac{r^{\,2}+7r\Gamma_{\alpha \beta}^{ij}+6(a_{\alpha}^{\,2}+a_{\beta}^{\,2}+7(-1)^{i+j}a_{\alpha} a_{\beta})}{(r+\Gamma_{\alpha \beta}^{ij})^7} \Bigg]\Bigg\}, \end{align} \tag{ 1 }$$

where $\Gamma_{\alpha \beta}^{ij} \equiv (-1)^i a_{\alpha}+(-1)^j a_{\beta}$, and a_α is the particle radius, as illustrated in figure 2(a). Here, without loss of generality, the length σ is set to unity and all distances are determined in units of σ. The factors $N_{\alpha \beta}$ are appropriately chosen constants. For self-interactions between colloidal particles (W_AA and W_BB in equation (1)), we set $a_A = a_B = 5$ and $N_{AA} = N_{BB} = 39.478$. The repulsive interaction between colloids, W_AB , is modeled by setting $a_A = a_B = 7$ and $N_{AB} = 39.478$, so that the colloids appear as larger spheres to one another, introducing a 'nonadditivity' to the A and B particle mixture. The solvent-colloid interaction, $W_{AS} = W_{BS}$, has $a_{A,B} = 5$ and $a_S = 0$, with $N_{AS} = N_{BS} = 75.398$. Finally, the solvent self-interaction potential W_SS has $a_S = 0$ and $N_{SS} = 144$ in equation (1). Each interaction potential $W_{\alpha \beta}$ is truncated (held constant) for all $r \gt \sigma_{\alpha \beta}$, with $\sigma_{\alpha \beta}$ chosen to be the minimum of the interaction potential, so that all potentials $W_{\alpha \beta}(r)$ smoothly approach a constant value at $r = \sigma_{\alpha \beta}$. The truncation distances for the various interaction types read: $\sigma_{SS} = 2^{1/6}$, $\sigma_{AA} = \sigma_{BB} = 10.58113\ldots$, $\sigma_{AB} = 14.5771\ldots$, and $\sigma_{SA} = \sigma_{SB} = 5.85972\ldots$. Note that these truncation distances $\sigma_{\alpha \beta}$ are all slightly larger than $a_{\alpha}+a_{\beta}$, meaning that the repulsive interactions act over a narrow range of distances. The interaction potentials $W_{\alpha \beta}$ for the different choices of α and β are shown in figures 2(b) and (c). Note the sharp increase in the potentials $W_{\alpha \beta}$ at values of r near $a_{\alpha}+a_{\beta}$. This divergence scales as $W_{\alpha \beta} \sim (\Delta r)^{-12}$ (as in the repulsive portion of the Lennard-Jones potential), with $\Delta r = r-a_{\alpha}-a_{\beta}$ as shown in figure 2(c). These purely repulsive potentials, as derived from the Lennard-Jones potential, are also known as Weeks–Chandler–Anderson potentials [25].

Zoom In Zoom Out Reset image size

We implement the MD in LAMMPS, making use of parallelization on graphics processing units (GPUs) [26, 27]. In addition to the interaction potentials in equation (1), each particle experiences a random force and a friction due to a Langevin thermostat [28]. We also introduce solvent particles S which do not experience the applied force, but allow for the colloidal particles to dissipate the applied force via collisions. All particles move under constant volume and temperature, where the temperature is maintained by coupling the system to the thermostat, such that the motion of each particle $i$ can be described as

$$\begin{equation} m_i\frac{{\mathrm{d}}\mathbf{v}_{i}(t)}{{\mathrm{d}}t} = \mathbf{F}_i(t) - \gamma m_i \mathbf{v}_i(t) + \mathbf{F}_i^R(t), \end{equation} \tag{ 2 }$$

where m_i is the particle mass, $\mathbf{v}_i(t)$ is the particle velocity, and $\mathbf{F}_i(t)$ denotes the net deterministic force acting on the $i{\mathrm{th}}$ particle, as derived from the interaction potential in equation (1). The thermal fluctuations are implemented via a stochastic force $\mathbf{F}_i^R(t)$ with Gaussian, white noise components, zero average value $\langle \mathbf{F}_i^R(t)\rangle = 0$, and correlations $\langle \mathbf{F}_i^R(t) \cdot \mathbf{F}_i^R(t^{^{\prime}})\rangle = 4 k_{\mathrm{B}}T m_i\gamma \delta (t-t^{^{\prime}})$. The friction coefficient γ is set to γ = 1 τ⁻¹ where τ is the reduced time unit $\tau = \sigma(m/k_{\mathrm{B}}T)^{1/2}$. The mass of A and B particles is $M = 78.54m$, where $m$ is the solvent particle mass. An external force with a magnitude F ranging from 0 to 2.0 $k_{\mathrm{B}}T/\sigma$ is applied to A and B particles in the $\hat{\mathbf{x}}$ direction. Note that the thermostat is applied in both x and y directions for the solvent particles S but only in the $\hat{\mathbf{y}}$ direction for colloidal particles A and B, with the x component of the associated random force set to zero. A similar technique was implemented when performing MD under shear conditions [29, 30]. Removing the thermostat (both the friction and the random force) along the drive direction ensures that the applied force is uniform for all colloidal particles. Periodic boundary conditions in both directions are implemented.

The systems were initialized with a random distribution of A and B particles with locations initially on a square grid. The small solvent particles were then added in as a dilute suspension which was rapidly compressed to a fixed solvent number density. We used a total of 1600 colloidal particles (800 per species) and a variable number of solvents, depending on the overall simulation box size which was used to vary the colloidal density. For the low density systems in figures 1(a) and (b), we have 274 888 solvent particles and a box size L × L with L = 780, while the higher density systems in figures 1(c) and (d) have 82 418 solvent particles and L = 520. The velocity-Verlet algorithm with a time step of $\Delta t = 0.01~\tau$ was used for integrating the equations of motion in equation (2) for these initial conditions. The simulations proceeded up to times of at least $10^6\;\tau$ (10⁸ total time steps). We will now set $k_{\mathrm{B}}T = \tau = 1$ for convenience and report times in units of τ and energies (forces) in units of $k_{\mathrm{B}}T$ ($k_{\mathrm{B}}T/\sigma$).

We verified that the explicit solvent provides a friction force along the drive direction by calculating the colloidal particle velocities as a function of the applied force F. We found that the average velocity and pressure rapidly approach constant values, on a time scale shorter than 10⁴ (early in the simulation runs). We find a linear increase of the colloidal particle velocity with F, as shown in the inset of figure 3(a), so that $F = \gamma_{\mathrm{e}}M \bar{v}_x$, with $\gamma_{\mathrm{e}}$ an effective friction coefficient supplied by the small solvent particles S along the drive direction. We find that the friction generally increases with decreasing colloidal particle density. The relationship is shown in figure 3(a). This makes sense as, at lower densities, there are more solvent particles available to provide friction for each colloidal particle. At higher densities, the colloids crowd out the solvents. This is also evident in the velocity snapshots in figure 3(b), where we observe that large, dense clusters of colloidal particles move faster than the more dilute regions for both high and low density systems. We will explore these density variations in more detail section 4. For the lower density ρ ≈ 0.2065, we find a friction coefficient $\gamma_{\mathrm{e}} = 0.560 \pm 0.009$. For the higher density ρ ≈ 0.4647, we find $\gamma_{\mathrm{e}} = 0.40 \pm 0.01$. Note that these coefficients are about half of the friction coefficient γ = 1 supplied by the Langevin thermostat in the direction perpendicular to the drive. It would be possible to increase $\gamma_{\mathrm{e}}$ by supplying more solvent particles, which would cost significant computational resources.

Zoom In Zoom Out Reset image size

We have also explored the effects of removing the thermostat entirely from the colloidal particles A and B. We find no qualitative differences and have opted to keep the thermostat in the perpendicular direction throughout the rest of the paper. In addition, we used a momentum-conserving thermostat and found that, in this case, the applied force continuously supplies energy into the system so that both the colloids and solvent particles accelerate, never reaching a stable velocity. It would be interesting to explore the effects of modifying the thermostat in the future.

It is clear that our explicit solvent is compressible and has density variations, leading to the different effective frictions for the colloids as shown in figure 3(a). In an experiment, the colloidal mixture would likely have an incompressible solvent (e.g. water). In this case, hydrodynamic interactions will generate inter-particle interactions between the colloids. For example, moving colloids will generate hydrodynamic flows which can entrain adjacent colloidal particles [31], also effectively decreasing the friction of denser particle clusters. So, although our MD simulations are closer to experimental realizations than previous studies of this system using a lattice [18, 19], there may be important differences. It may be necessary, for example, to add small particles to the colloidal mixture in an experiment to generate the effects provided by our solvent particles S.

We begin by looking at colloidal mixtures with no applied force, which we expect to coarsen at sufficiently high density. It is easier to analyze such structures in Fourier space, so we calculate the static structure factors $S_{\pm}(\mathbf{q},t)$ associated with the 'charge' (−) and density (+) of the colloidal species. Specifically,

$$\begin{equation} S_{\pm}(\mathbf{q},t)= \frac{1}{N} \sum_{i,j=1}^N \,w^{\pm}_i w^{\pm}_j {\mathrm{e}}^{i \mathbf{q} \cdot [\mathbf{r}_i(t)-\mathbf{r}_j(t)]}, \end{equation} \tag{ 3 }$$

where we sum over the N colloidal particles at locations r_i . Here, the weights for the charge structure factor are $w^{- }_i = \pm 1$ if colloidal particle i is of type A or B, respectively. Conversely, for the density structure factor, we have equal weights $w^{+}_i = 1$ for both types. It is possible to average over many times t to get better statistics for these structure factors, although one must be careful that the static structure factor does not evolve substantially over time. This works best for larger force values for which the coarsening is arrested, as described in the next section. The structure factors are calculated on a square grid in q-space, with spacing $\Delta q = 1/150$ (in units of $1/\sigma$).

Note that our definitions of $S_{\pm}$ in equation (3) have a constraint at $\mathbf{q} = 0$: $S_{\pm}(0,t) = N^{-1}\sum_{i,j}\,w^{\pm}_i w^{\pm}_j$. For equal numbers of A and B particles, this means $S_-(\mathbf{q},t) = 0$ and $S_+(\mathbf{q},t) = N$ exactly. This creates an issue as these constraints generate discontinuities in our structure factor measurements. We therefore replace the values of $S_{\pm}$ at the $\mathbf{q} = 0$ grid site with the value at a neighboring grid site. The structure factors near $\mathbf{q} = 0$, then, become smoothly-varying. Also, due to the computational difficulty in running these simulations, we were limited to using single simulation runs per condition (colloidal density, force, etc). We used system sizes of at least $500 \times 500$, in units of the Lennard-Jones distance σ. Structure factors calculated from single runs in this way are noisy in the q-space. Therefore, in addition to the time-averaging, we perform a smoothing of the structure factor results by convolving the structure factor with a Gaussian kernel with standard deviation $2\Delta q = 1/75$. The removal of the $\mathbf{q} = 0$ discontinuities in $S_{\pm}$ is also helpful in this regard. Finally, we checked that larger system sizes did not significantly modify our results by running some of the systems at twice the (linear) system size L. For example, we checked that the structure factor calculation yields compatible results for box sizes L × L with L = 520 and L = 1040 for the dense colloidal mixtures (ρ ≈ 0.4647).

First, we consider results for the 'charge' structure factor $S_-(\mathbf{q},t)$ at large particle densities and F = 0. Here, we would expect domains of $A$ and $B$ particles to coarsen, much like an Ising model with conserved dynamics at low temperature. Without an applied drive to specify a direction, we expect isotropic domains in our simulation and a circularly symmetric average structure factor $\langle S_-(\mathbf{q},t) \rangle$, which is what we find, apart from some fluctuations, in the color maps of figure 4(a). As shown in the bottom panel, the 'ring' of the structure factor at $|\mathbf{q}| = q^*$ gets smaller with increasing time and more difficult to resolve. This is expected for a coarsening system. We can track the position of the peak more systematically by performing a circular average $\overline{\langle S_-\rangle}(q,t)$, as shown in figure 4(b). We find that the shift in the peak is consistent with the Lifshitz–Slyozov–Wagner coarsening characteristic of conserved dynamics [32, 33], for which $q^* \sim t^{-1/3}$. We also show in the inset of figure 4(b) the characteristic length scale $\lambda^* \equiv 2\pi/q^*$ associated with the peak at $q = q^*$. The red dashed line indicates the $\lambda^* \sim t^{1/3}$ scaling. Note that at the larger times t, the scale $\lambda^*$ approaches values close to the system size, so we expect to see finite size effects at these larger times. At smaller densities, such as those shown in figure 1(a) where the colloidal species remain disordered, we find a similar behavior except that at long times t, the peak stops growing and shifting toward smaller q values. Instead, we approach a steady state with a fixed peak corresponding to the average domain size in the disordered phase. We may thus expect a phase transition in this system (from a phase with mixed A and B particles to a phase where the particles phase separate into distinct A and B domains) somewhere between ρ ≈ 0.2 and ρ ≈ 0.5.

Zoom In Zoom Out Reset image size

Adding an applied force to the system significantly changes the structure of the AB mixture: The applied force F drives an instability which creates large density fluctuations (see figure 3(b)). Moreover, we find that the phase separation between A and B species at high colloidal particle density ρ is arrested in the presence of the force F > 0. Instead, the system prefers to form anisotropic domains of A and B particles with a particular size, as can be seen by comparing figures 1(c) and (d). We can check this property in another way by starting with a fully phase separated initial condition. We find that under an applied force, the phase separated state breaks up into domains and eventually settles into a non-equilibrium steady state similar to the one for an initially well-mixed A and B particle configuration. The snapshots are shown in figure 5(a) where we consider two different fully phase separated initial states. Note that in both cases (top and bottom panels), the system settles into a state that looks similar to figure 1(d). Let us now quantify these properties by calculating the structure factor of the system under an applied drive for the mixed initial condition.

Zoom In Zoom Out Reset image size

The circular symmetry of our structure factors $\langle S_-(\mathbf{q},t)\rangle$ is broken as the domains of A and B particles become anisotropic. In our simulations with F > 0, we average over the simulation trajectory from time $t = 9 \times 10^5$ to $t = 9.8 \times 10^5$, where we verify that the structure factor remains relatively time-independent. The averaged charge structure factors $\langle S_-(\mathbf{q},t)\rangle$ are shown in figures 5(b) and (c). We note two important differences compared to the F = 0 case: first, the peak at small $|\mathbf{q}|$ is no longer symmetric. Instead of a ring, we find two peaks arranged either approximately parallel (figure 5(b) at low densities) or perpendicular (figure 5(c) at high densities) to the drive direction $\hat{\mathbf{x}}$, depending on the force and particle density. Thus, we can no longer perform the circular average in q-space, as we did for F = 0. Second, the two peaks have a smaller amplitude for larger applied forces and do not grow with increasing time, indicating an arrested coarsening of the A and B colloidal domains. We will now study this more systematically.

We can check that the coarsening is indeed arrested by calculating $\langle S_-(\mathbf{q},t) \rangle$ at different times t. To do this, we can take a slice through the q-plane that passes through the origin and the two central peaks in $\langle S_-(\mathbf{q},t)\rangle$ (shown in the lighter colors in panels of figures 5(b) and (c)). The results are shown in figure 6 for the larger colloidal particle density and a relatively small applied force F = 0.5. We see in figure 6(a) that the peaks at small values of $|\delta q| \approx 0.1$ fluctuate but remain at roughly the same positions, as seen more clearly in the inset. The full q dependence near these central peaks is shown for the earliest and latest times in figure 6(b). We see that, overall, the $\langle S_-(\mathbf{q},t)\rangle$ peaks near $|\delta q| \approx 0.1$ remain at well-defined positions for different times t. One can compare this behavior to the large difference at the earliest and latest times for the F = 0 case in figure 4(b). There, the peak at small $|\mathbf{q}|$ continues to grow with time t and shifts rapidly toward smaller q (note the logarithmic scale). The peaks at $|\delta q| \approx 0.7$ in figure 6(a) do shift in time, moving toward larger $|\delta q|$ values with increasing time. This corresponds to a decrease in inter-particle spacing, consistent with the formation of dense particle clusters (see figure 1(d)). The shift in the position of this peak can be compared to the position of the peak at $|\mathbf{q}| = 0.6$ in figure 4(b), which remains stationary in time. In other words, along with arresting the coarsening of A and B domains, the applied force generates a change in local particle densities, separating the system into regions of high and low colloidal particle density. We will explore this in more detail in the next section.

Zoom In Zoom Out Reset image size

Having established that an applied force F arrests the coarsening of colloidal domains, we may look at the steady-state profiles of $S_-(\mathbf{q},t)$ by averaging over late simulation times t. The results for low and high particle densities are shown in figures 7(a) and (b), respectively, for varying applied force F, where we again evaluate the structure factors through a line that passes through the two central peaks (see right panels of figure 6). The application of a force at low densities decreases the characteristic domain size of the different colloidal species, as evidenced by the shifting of the central peaks toward larger $|\delta q|$ values in figure 7(a). At larger colloidal particle densities, there is a similar shift which saturates at sufficiently large forces $F \gtrsim 0.25$. These results are partially analogous to what is observed in a driven Widom–Rowlinson mixture (i.e. a binary mixture with short-range repulsive interactions) on a lattice, studied using Monte Carlo simulations in [18, 19]. However, on the lattice, the mixture (at sufficiently large particle densities) jams up into stripes of alternating $A$ and $B$ particle domains, oriented perpendicular to the drive direction. The stripe size decreases with increasing drive magnitude. The precise reason for this stripe formation remains poorly understood. In our more realistic MD simulations described here, we do not find the ordered stripes. Instead, the domains of A and B particles reach a certain characteristic size depending on the magnitude of the force F. This size decreases with increasing F, much like the stripe width in the lattice simulations. The MD results are more reminiscent of the behavior of the lattice model at lower densities [19], where disordered clusters (instead of stripes) form, with a characteristic width along the drive direction. This characteristic width yields two peaks in the (charge) static structure factor $\langle S_-(\mathbf{q},t)\rangle$. We also find these peaks in our MD simulations, although the peaks are not always aligned along the drive direction as they are in the lattice model case (see [19] for details). We shall now explore the peak position and orientation in more detail.

Zoom In Zoom Out Reset image size

We show the positions $q_{\mathrm{max}}$ and orientations $\theta_{\mathrm{max}}$ of the peaks in the structure factor $\langle S_-(\mathbf{q},t) \rangle$ in figure 8, where panel (a) illustrates these two features. The $q_{\mathrm{max}}$ gives the distance between the origin in q space and the position of each peak (which are both equidistant from $|\mathbf{q}| = 0$). The angle $\theta_{\mathrm{max}}\in (-\pi/2,\pi/2]$ gives the orientation of the two peaks relative to the drive direction $\hat{\mathbf{x}}$. In panel (b) of figure 8 we see that the peak positions at $|\mathbf{q}| = q_{\mathrm{max}}$ increase roughly linearly with the applied force F for low densities, while saturating to an approximately constant value for the higher colloidal particle densities. The low density results are analogous to those found for the lattice gas, where a disordered 'microemulsion' phase develops, characterized by two peaks in the static structure factor at positions $q_x = \pm q_{\mathrm{max}} \propto F$ (along the drive direction with $\theta_{\mathrm{max}} = 0$) [19].

Zoom In Zoom Out Reset image size

The orientation of the peaks in $\langle S_-(\mathbf{q},t)\rangle$ at small $|\mathbf{q}|$ is subtle. We already observed in figure 5 that the orientation of the peaks can change depending on colloidal particle density. We also find that the magnitude of the force makes a difference. A larger force F > 1 rotates the peaks perpendicular to the drive direction ($\theta_{\mathrm{max}} \rightarrow \pi/2$). This is consistent with the formation of lanes of colloidal particles of opposite type forming parallel to the drive. Note that this behavior is in marked contrast to the lattice model results discussed in [19], where the peaks in $\langle S_-(\mathbf{q},t)\rangle$ always remain along the drive direction with $\theta_{\mathrm{max}} = 0$. A more detailed analysis of the peak orientations is shown in figure 8(c). At smaller forces F and low colloidal particle densities, the orientation of the peaks is consistent with domains lying perpendicular to the drive direction, $\theta_{\mathrm{max}} \approx 0$ (green points in figure 8(c)). Conversely, at higher colloidal densities and large forces, the peaks are oriented perpendicular to the drive direction ($|\theta_{\mathrm{max}}| \approx \pi/2$), corresponding to domains of A and B particles elongated along the drive direction (red points in figure 8(c)). Finally, note that although large particle densities ρ and large forces F both tend to orient the structure factor peaks perpendicular to the drive direction ($\theta_{\mathrm{max}} \approx \pi/2$), increasing ρ can either increase or decrease $q_{\mathrm{max}}$ (depending on the value of $F$), while increasing F tends to increase $q_{\mathrm{max}}$ (because the force tends to break up the coarsening domains). Thus, the parameters $\rho$ and $F$ play different roles.

It is interesting that our results here suggest that the orientation of the structure factor peaks depends on both the density ρ and the force F, by contrast to the lattice gas simulations in [18, 19] where the peaks are always displaced from the origin along the drive direction. We may speculate that the more realistic dynamics of our simulations here allow for the colloidal domains to rotate more freely, so that cluster formation is not locked into the drive direction as it is on the lattice. In addition, it is possible that the large particle density fluctuations play an important role, as the MD simulations exhibit very dense particle clusters, such as the locally crystalline patches shown in figure 1(d), which is not possible on the lattice. We also do not find ordered stripes at high densities, which robustly form in the lattice models [18, 19]. One possibility is that, at higher densities, the systems are kinetically arrested in our MD simulations, so that the ordered striped phases found in [18, 19] are restricted from developing.

We now consider the overall distribution of the colloidal particles. As observed in figures 1(b) and (d), the applied force generates large density fluctuations and local regions with liquid-like and gas-like regions. At higher densities, the liquid-like regions generally move faster than the gas-like regions, which are hampered by more solvent particles, as shown in figure 3(b). These fluctuations can be quantified by calculating the probability P(N) of observing N particles in some circular 'window' placed at random in the simulation box. To calculate P(N), we generated 50 random points and counted the number of particles within a certain distance (the window radius) from each point, taking into account the periodic boundary conditions. We calculated histograms from these 50 points chosen for 2000 different frames at the latest times in the simulation, over a total time interval of $\Delta t = 2 \times 10^5$. Examples of the resulting distributions P(N) for a window radius of 40 are shown in figure 9(a) for the larger colloidal particle density ρ ≈ 0.4647 (see figures 1(c) and (d)) for different force values F. Note that applying a force F > 0 significantly changes the distribution: it widens and becomes asymmetric as the force generates large density fluctuations in the system. Systems with smaller particle densities, such as those in figures 1(a) and (b), also exhibit this broadening of the P(N) distribution, but only for a sufficiently large force F.

Zoom In Zoom Out Reset image size

The density fluctuations can be further quantified by plotting the average $\langle N \rangle$ of each distribution P(N) versus the standard deviation σ_N for a variety of window sizes. For equilibrium systems with a single phase and average particle separation, we would expect that $\sigma_N \sim \sqrt{\langle N \rangle}$, in accordance with the central limit theorem. However, in a strongly driven system, we might expect long-range correlations to develop, generating giant density fluctuations which are characterized by $\sigma_N \sim \langle N \rangle^{\nu}$ with ν > 0.5. Such giant density fluctuations are a hallmark of active matter systems [34, 35]. We find an exponent ν ≈ 0.7 for our driven systems, as shown in figure 10(b). Note that a system undergoing a liquid-gas phase separation would also have particle density fluctuations consistent with ν = 1. Such a separation into dilute and dense regions occurs during mobility-induced phase separation (MIPS), the onset of which can often be interpreted as a liquid-gas transition [36]. A more thorough analysis would be required to see to what extent our system is similar to MIPS.

Zoom In Zoom Out Reset image size

One way our system is clearly different from active matter and MIPS is that our particles are not self-propelled and we have strong anisotropies induced by the applied drive: The particles cluster more perpendicular to the drive direction (forming dense 'lanes'). This can be evaluated by looking at the standard deviations $\sigma_N^{(x,y)}$ along the two directions (the drive direction $\hat{\mathbf{x}}$ and the perpendicular direction $\hat{\mathbf{y}}$) separately. The deviations $\sigma_N^{(x,y)}$ are measured by using long rectangular windows (spanning the system in the y and x direction, respectively) of different sizes. The results are shown in figure 9(c). Note that the blue curves, corresponding to particle density fluctuations $\sigma^{(x)}_N$ along the drive direction, are smaller than the perpendicular fluctuations $\sigma_N^{(y)}$ given by the black curves. At large force, however, we see large density fluctuations in both directions (solid lines in figure 9(c)).

The large density fluctuations are accompanied by changes in the density structure factor $S_+(\mathbf{q},t)$. The results are shown in figure 10 for the large particle density case. The key feature corresponding to large-scale density fluctuations is the appearance of a sharp peak at $|\mathbf{q}|=0$ in the presence of an applied force $F>0$. Also, we find that at larger $|\mathbf{q}|$, the 'ring' corresponding to local particle arrangements (light green circle for the F = 0 panel in figure 10(a)) becomes deformed at larger forces F, with higher intensity along the drive direction. This can be seen by comparing the intensity maps around the light green circles in the panels of figure 10(a). Such an asymmetry indicates that the driven systems tend to have local particle arrangements with orientations sensitive to the drive direction. This can be seen in figure 1(d), where the small crystalline domains of different colloidal species are largely oriented such that rows of particles align along the drive. A similar result was found for the smaller particle densities at a sufficiently large force (see figure 1(b)).

For smaller particle densities such as those in figures 1(a) and (b), we also find a peak at $|\mathbf{q}| = 0$, but only for sufficiently large forces $(F\gt1)$ at which we observe the large density fluctuations. It appears that these fluctuations are accompanied by a redistribution of the solvent particles away from dense colloidal particle clusters, but a detailed analysis of this possibility would require more simulations at different solvent concentrations. The appearance of the $|\mathbf{q}| = 0$ peak in $\langle S_+(\mathbf{q},t)\rangle$ for applied drives is also found in the Monte Carlo simulations of the lattice gas model [19]. There, the peak appears even at small forces and particle densities. As mentioned previously, it is possible that our MD simulations have kinetic barriers which prevent the system from settling into certain non-equilibrium steady states. We thus might be prevented from seeing the development of the small $|\mathbf{q}|$ peak at small forces and particle densities. We would expect the same to be true of experimental realizations of such driven systems.

In summary, we presented MD simulation results on a colloidal mixture of two particle types A and B with repulsive interactions, driven by a uniform applied field and immersed in a collection of solvent particles S. We found that the applied field generates two major effects: first, the coarsening of A and B particle domains (exhibited at zero force and sufficiently large colloidal particle density ρ) is arrested when the domains reach a certain characteristic size. Initially phase separated systems also break up into domains of this characteristic size. Second, the system exhibits large density fluctuations and separates into low density (gas-like) and high density (liquid-like or locally-crystalline) regions. These behaviors should be observable in an experimentally-realized colloidal system.

In addition, we verified that the coarsening at zero applied force (at sufficiently large ρ) is consistent with a conserved dynamics, where the characteristic domain size grows as $\lambda^* \sim t^{1/3}$. We also found that the characteristic domain size at non-zero applied forces gets smaller with increasing applied force (especially at low particle densities), much like in a lattice model of the driven Widom–Rowlinson mixture [18, 19]. We also found that at high colloidal particle densities, the system forms lanes of dense, quickly moving particles parallel to the driven direction. Conversely, the driven Widom–Rowlinson lattice gas studied in [18, 19] jams up into ordered stripes perpendicular to the drive direction. The effects of the applied field on the particle distributions in our simulations were quantified by calculating the 'density' and 'charge' structure factors $\langle S_{\pm}(\mathbf{q},t)\rangle$ that take into account the overall particle positions and the relative distribution of the two colloidal species (A and B), respectively. It would be interesting to understand precisely why our simulations do not produce the ordered stripes (perpendicular to the drive direction) found in the lattice gas simulations in [18]. It is likely that the MD simulations presented here are prone to kinetic arrest (perhaps due to the formation of rigid, locally crystalline clusters), so that this kind of long-range ordering is prevented. We would expect to see similar behavior in experiments.

The low colloidal particle density results suggest that the non-equilibrium 'microemulsion' phase [19], characterized by two peaks in the 'charge' structure factor (indicating anisotropic domains with a characteristic 'width') and a single peak at $|\mathbf{q}| = 0$ in the density structure factor, can be realized in realistic systems. A preliminary theory for such a non-equilibrium microemulsion phase, based on a lattice model, is presented in [19]. The basic idea is that charge and density fluctuations propagate at different velocities in the system due to the repulsive interactions between colloidal species. This velocity difference, when combined with thermal fluctuations, yields a peak in the (averaged) structure factor $\langle S_-(\mathbf{q},t)\rangle$ at a characteristic $\mathbf{q} = (q^*,0)$ (along the drive direction $\hat{\mathbf{x}}$), with $q^*$ proportional to the drive strength. In our more realistic simulations, we find that the peaks in $\langle S_-(\mathbf{q},t)\rangle$ (and the corresponding anisotropic disordered domain geometries) may not always have a distinct relationship to the drive direction, with different orientations depending on the particle density and drive strength. Generally, we find that the domains are more elongated along the drive direction for increasing forces and densities. It would be interesting in future work to evaluate the dynamic charge/density structure factors for this model in order to see if there is a difference in the velocities and diffusivities of the charge and density fluctuations. This would require more robust statistics for our simulation runs, which may be computationally prohibitive. We expect that our work inspires such an endeavor in a real colloidal suspension.

Our results also illustrate the importance of the colloid-solvent interactions in these strongly driven systems. As a MD simulation integrates the equations of motion for each particle simultaneously, any uniform applied field with an implicit solvent (modeled by a Langevin thermostat, say) will simply uniformly accelerate all of the particles, with no effect on their relative dynamics or structure. Even introducing a large opposing friction force to mitigate the acceleration does not work because each particle will move with the same fixed velocity as set by the friction coefficient and applied force magnitude. Instead, here we have opted for a more computationally intensive approach and used an explicit solvent which does not experience the applied force directly. Collisions between colloid and solvent particles provide a friction for the colloids along the drive direction, without introducing an artificial, uniform velocity or acceleration. In addition, the choice of thermostat proved to be important: Our use of a Langevin thermostat for the solvent particles allowed for dissipation of the energy injected via the applied drive. Removing the thermostat or using momentum-conserving thermostats leads to acceleration of the entire system of colloids and solvents, preventing the establishment of a non-equilibrium steady state with a uniform average velocity.

There is clearly much more to explore in the future. We only considered a few simulation runs in detail, starting from either well-mixed or completely phase separated initial conditions. It would be important to check our results against larger system sizes and multiple runs. This would allow us to calculate more subtle statistical quantities, such as the dynamic structure factors. Given the similarities and differences between our results and the lattice simulations of a related system, it would be interesting to build an off-lattice Monte Carlo simulation to see if the ordered striped phases are possible. We did not observe them here, possibly due to kinetic arrest of the system into a more disordered phase. One potentially interesting avenue for future research would be to modify the initial conditions to see if striped phases could be stabilized under the applied drive.

MOL thanks R K P Zia, R Dickman, and I Schmidt for helpful discussions. MOL and DB are supported through the National Science Foundation, Division of Molecular and Cellular Biosciences (MCB), under Contract No. 2219289. Computational support was provided by the University of Tennessee. MD simulations research was supported by the Center for Nanophase Materials Sciences (CNMS), which is a US Department of Energy, Office of Science User Facility at Oak Ridge National Laboratory. This research used resources of the Oak Ridge Leadership Computing Facility (OLCF) at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725.

The data that support the findings of this study are available upon reasonable request from the authors.

As a control, it is useful to see what happens when we consider a single colloidal species (particle A, say) in the presence of the solvent particles S. This is easily accomplished by removing the larger repulsion between A and B particles so that all the repulsive interactions become additive. In this case, as there is no difference between A and B particles, the only relevant structure factor is for the overall density, $\langle S_+(\mathbf{q},t)\rangle$. We may also calculate the particle clustering properties. The results are shown in figure 11.

Zoom In Zoom Out Reset image size

For the system with a single colloidal particle type, we find a similar increase in the particle distribution fluctuations σ_N at non-zero forces F > 0 (figure 11(a)), and a corresponding development of a peak in $\langle S_+(\mathbf{q},t)\rangle$ at $\mathbf{q}=0$ (figure 11(b)). Interestingly, the large $|\mathbf{q}|$ behavior remains rotationally symmetric when a force is applied, with the outer 'ring' only inflating slightly compared to the F = 0 case, as can be seen by comparing the ring to the black dashed line in the inset of figure 11(b). The inflating ring suggests that there are more dense particle arrangements when the force is applied. However, unlike the case of the binary mixture of colloids, this system does not exhibit the localized crystallization and corresponding anisotropies in $\langle S_+(\mathbf{q},t)\rangle$. Also, the standard deviations $\sigma_N$ increase more graduatlly with the applied force $F$ in this single colloidal species case.