Luis Ruiz
Pestana
a,
Natalie
Minnetian
d,
Laura Nielsen
Lammers
bc and
Teresa
Head-Gordon
*adef
aChemical Sciences Division, Lawrence Berkeley National Laboratory, University of California, Berkeley, USA. E-mail: thg@berkeley.edu
bEarth and Environmental Science Area, Lawrence Berkeley National Laboratory, University of California, Berkeley, USA
cDepartment of Environmental Science, Policy, and Management, University of California, Berkeley, USA
dDepartment of Chemistry, University of California, Berkeley, USA
eDepartment of Bioengineering, University of California, Berkeley, USA
fDepartment of Chemical and Biomolecular Engineering, University of California, Berkeley, USA
First published on 2nd January 2018
When driven out of equilibrium, many diverse systems can form complex spatial and dynamical patterns, even in the absence of attractive interactions. Using kinetic Monte Carlo simulations, we investigate the phase behavior of a binary system of particles of dissimilar size confined between semiflexible planar surfaces, in which the nanoconfinement introduces a non-local coupling between particles, which we model as an activation energy barrier to diffusion that decreases with the local fraction of the larger particle. The system autonomously reaches a cyclical non-equilibrium state characterized by the formation and dissolution of metastable micelle-like clusters with the small particles in the core and the large ones in the surrounding corona. The power spectrum of the fluctuations in the aggregation number exhibits 1/f noise reminiscent of self-organized critical systems. We suggest that the dynamical metastability of the micellar structures arises from an inversion of the energy landscape, in which the relaxation dynamics of one of the species induces a metastable phase for the other species.
In this work, we investigate the phase behavior of a model inspired by a binary mixture of ionic particles (K and Cs) that have different size, rK < rCs, and are under semiflexible planar confinement (Fig. 1a). The critical feature of the system is that the activation energy barrier to diffusion of the larger (Cs) particles is higher than that of the smaller (K) particles, given the same local environment, and the activation energy barrier of all the particles decreases with the local concentration of the larger particles because they mechanically expand the semiflexible slit locally (Fig. 1b). Non-equilibrium effects will arise when the particles are confined between two semiflexible interfaces, which introduces a non-local coupling between the particles and their diffusion energy barriers. These abstracted features of the model have been derived from a recent study of the kinetics of ion migration in anhydrous clay interlayers,12 which we investigate here in the limit of infinite screening (i.e. the electrostatic interactions between the ions are negligible).
We find that for a broad range of conditions, quantified below, the model reaches a state characterized by the cyclical formation and dissolution of metastable micelle-like clusters, with smaller K particles in the core and larger Cs particles in the corona. We call this a cyclical self-organized (CSO) state. We show that the cyclical nature of the dynamics arises from the inversion of the energy landscape13,14 of one particle type induced by the aggregation of particles of the other type, but by a mechanism in which the relaxation dynamics of the K particles is arrested by the Cs particles. We denote this novel mechanism dynamical inversion of the energy landscape (DIEL). Furthermore, the power spectrum of the fluctuations of the aggregation number in the CSO state exhibits a power law tail, reminiscent of self-organized critical systems.15
In what follows, we discuss first the phenomenology of the model, and then quantitatively analyze the phase behavior and dynamics of the system and its dependence on the relevant model parameters. We also derive a simple analytical form to predict the transition of the system to the CSO state and summarize our findings in a phase diagram that could guide experimental efforts towards reproducing the anomalous phase behavior observed here.
At each step of the KMC simulation we assign to each particle a characteristic waiting time determined using by eqn (1) and (2) in the main text. For each particle, a waiting time is randomly generated from an exponential distribution, effectively assuming a Poisson process. Then, we attempt to move a particle with the shortest waiting time to a neighboring lattice site selected at random. If the selected site is occupied by another particle, new waiting times are generated for all the particles, and a new move corresponding to the particle with new shortest waiting time is attempted. For each step we iterate this procedure until a move is successfully realized. The time step is calculated as the sum of the waiting times of all the attempted moves. We output the configuration of the lattice every 103 steps and run each simulation for a total of 3.5 × 108 steps.
(1) |
Ea = Ai − BflocCs | (2) |
The parameter Ai is the activation energy barrier of particle type i in a local environment where flocCs = 0, whereas B determines the strength of the dependence on the fraction of Cs within a local neighborhood of radius rloc. In what follows, we show that the transition to the non-equilibrium cyclical state is governed by the dimensionless parameter ΔA/B, where ΔA = ACs − AK. Further details of the KMC model and governing equations are described in detail in the Computational methods section.
In the trivial case where B = 0 and ΔA < kBT (i.e. the activation barriers of Cs and K are approximately equivalent, and there is no dependence of the activation energy barriers on the local environment) the entropy of mixing is maximized and the particles exhibit diffusive dynamics, as expected from a system without pairwise interactions in thermodynamic equilibrium (Fig. S1–S3†). Increasing ΔA alone has the effect of separating the timescales at which Cs and K diffuse, but the dynamics of both remain diffusive. Because increasing ΔA makes the Cs particles to be exponentially less mobile than K, at high occupancy densities and ΔA ≫ kBT, Cs particles act as fixed obstacles to the much faster K particles, which leads to sub-diffusive behavior of K. This phenomenon is equivalent to crowding-induced sub-diffusion observed in a variety of systems.16,17 Because in our system the particles aggregate into clusters with high local concentrations, sub-diffusive dynamics of individual particles is manifested even at moderately low occupancy densities.
When the dimensionless parameter ΔA/B takes finite values, we observe a phase transition from the well-mixed diffusive state (ΔA/B → ∞) to a state where K particles aggregate into clusters. In an initial well-mixed state, the local environment of every particle is the same on average, (flocCs)K = (flocCs)Cs = 0.5, which makes the activation energy barriers of K particles, EKa = AK − B(flocCs)K, lower than those of Cs, ECsa = ACs − B(flocCs)Cs, by ΔA. Therefore, the well-mixed configuration represents an unstable state for the K particles such that they drift towards regions with high concentrations of K (i.e. clusters of K) where they diffuse much more slowly. If we examine the scaling of the fluctuations in the number of particles, ΔN, averaged over time in a subsystem as a function of subsystem size, N0, i.e. ΔN ∼ N0α, we see that K exhibits out-of-equilibrium fluctuations (i.e. α > 0.5) for all finite values of ΔA/B (Fig. 2). We show in the ESI† that purely diffusive systems exhibit equilibrium fluctuations (Fig. S2†). The mechanism of aggregation of K particles, where the particles accumulate in low-diffusivity regions and diffuse slowly in the regions where they accumulate, is analogous to the positive feedback responsible for clustering in models of active matter where the mobility of the particles decreases with their concentration.9,18,19
For Cs particles to start diffusing once K has aggregated into clusters, the activation energy barrier of Cs dispersed in the lattice needs to be lowered such that ECsa = EKa. This condition can be written as:
(ΔA/B)c = ΔflocCs | (3) |
The micelle-like clusters are not stable, however, but kinetically arrested in a metastable state. The aggregation of Cs has the effect of increasing the local concentration of Cs for K particles at the boundaries of each cluster, thus enhancing the mobility of the K particles, which become unstable with respect to their aggregated current state. However, their relaxation path is hindered by the Cs particles in the corona. Eventually, a finite-size fluctuation occurs (i.e. the Cs corona becomes significantly disrupted) and the metastable micelle-like cluster disperses in the lattice, which brings the system back to a state where both species are mixed in the lattice, before the self-organization resumes the next cycle. The destabilizing effect that the dynamical relaxation of one particle type has on the other can be rationalized considering a recently developed concept: the inversion of the energy landscape (IEL). Originally, IEL was invoked to explain the anomalous formation of metastable states through spinodal decomposition, triggered by the variation of a thermodynamic external parameter.13,14 In our case, the inversion of the energy landscape of each particle type is dynamically driven by the aggregation of particles of the other type. This dynamical inversion of the landscape (DIEL) leads to the cyclical behavior illustrated in Fig. 3. (Movies of simulations showing this phenomenon are provided in the ESI†). We note that the phenomena we observe for finite ΔA/B is insensitive to initial conditions by showing the same features are exhibited when the binary system is initially phase-separated (Fig. S4†). Because the barrier to escape the transient state can be quite large (particularly for high values of B, high occupancy density of the lattice, and/or initial phase separated conditions), the system can take a long simulation time to manifest the cyclic behavior. Because our goal is to investigate the cyclical non-equilibrium behavior rather than the initial transient, we decided to use well-mixed starting configurations to more straightforwardly quantify the result.
To show that the DIEL phenomena is robust over a range of system parameters, and in order to gain quantitative insight into the order parameters that drive the phase behavior and dynamics of the system to make predictive outcomes, we analyze the effect of varying ΔA/B as well as its dependence on occupancy density, ρ, which is the fraction of lattice sites occupied by particles. Fig. 4a and c show the probability distribution of the aggregation fraction, fiagg ∈ [0,1], defined as the number of particles of type i in an aggregated state divided by the total number of particles of that type in the system, for different values of ΔA/B. For both types of particles the probability distributions display Gaussian profiles with well-defined averages, indicating that the average fraction, fiagg, could be a reasonable order parameter for the system. In Fig. 4b and d, fiagg is plotted as a function of ΔA/B. We find that K aggregates into clusters across all finite ΔA/B values (Fig. 4b). The nature of the K clustering is classifiable as non-equilibrium aggregation (Fig. S5a and b† shows that the aggregation fraction remains approximately constant through the simulation and that the average cluster size is inversely proportional to the instantaneous number of clusters). The aggregation of Cs, on the other hand, exhibits a relatively sharp transition at ΔA/B ≈ 0.7 (Fig. 4d). For ΔA/B > 0.7, Cs exhibits stochastic aggregation due to equilibrium fluctuations, where the average cluster size is small and independent on the instantaneous number of clusters (Fig. S5c and d†). For ΔA/B < 0.7 we observe the non-equilibrium aggregation behavior of Cs. This analysis provides quantitative evidence of the non-equilibrium and metastable phenomenology described earlier, namely, the aggregation first of K into clusters, followed by aggregation of Cs for values of ΔA/B ≤ (ΔA/B)c.
Fig. 4 Aggregation behavior as a function of ΔA/B for ρ = 0.3. Panels (a) and (c) show the probability distributions of the aggregation fraction, fiagg, for K and Cs, respectively. Panels (b) and (d) show the average aggregation fraction, 〈fiagg〉, as a function of ΔA/B for K and Cs, respectively. The results shown in the figure are for ρ = 0.3. It is interesting that the average stochastic aggregation fraction for K decreases slightly, but consistently, for the lower values of ΔA/B, which can be explained by the fact that, the lower ΔA/B is, the earlier in the clustering process Cs would arrest the aggregation of K by surrounding the existing clusters. Because the values of fiagg have been corrected for background values of aggregation under diffusive conditions, in the diffusive case, when ΔA/B → ∞, fiagg = 0. It is worth noting that the values of fiagg have been corrected according to , where fiagg,eq is the background aggregation (Fig. S3b†) and is the uncorrected value. |
Because it is the aggregation of Cs that triggers the metastability in the micelle-like clusters, we can use fCsagg as an order parameter to identify additional features of the transition to the CSO state (Fig. 4d). In the CSO state, the system adopts a dynamical behavior whereby micelles form and dissolve as a function of time, which we hypothesize, originates from the dynamical inversion of the energy landscape, and where many micelles at different stages of their lifetimes coexist. A particularly important question is whether the dynamics in the CSO state exhibit any particular time-scale or periodicity. To answer that question, we calculate in Fig. 5 the power spectrum, S(ω), of the time signal of the order parameter, fCsagg. For all cases, we observe that the power spectrum exhibits a power law tail, which indicates the absence of a characteristic frequency governing the dynamical behavior. We note that the inverse power law behavior that arises from the dynamical metastability of the micelles due to DIEL, is similar to that observed in self-organized critical systems.21
We also investigate the behavior of the system in the CSO state for values of the occupancy density, ρ, ranging between 0.2 and 0.7 (Fig. 6a). We observe that the K clusters become more compact as the occupancy density increases. What is conserved across systems with different occupancy density is the thickness of the compact K phase, RK ≈ rloc (Fig. S6a†). For example, for lower ρ values the large clusters that form display hollow cores that preserve RK ≈ rloc. Using simple geometry and assuming that the micellar clusters are circular with a K core of radius RK and a corona of Cs of thickness DCs, the condition that each micelle has the same amount of Cs than K leads to which agrees reasonably well with the numerical observations. Given this relation between RK and DCs, the question remains as to why the system tends to form micelles of size RK ≈ rloc? If the system evolves such that the activation energy barriers of Cs particles are maximized, or equivalently such that the Cs diffusivity is minimized, the clusters should evolve to a state where the total overlap area between the local neighborhood of Cs in the corona and the K particles in the core, , is maximized. In fact, this analysis, shown if Fig. S6b,† reveals that corresponds to RK/rloc ≈ 1, which agrees with the simulation data histogrammed in Fig. S6a.†
Finally, in order to provide a simple analytical model for the threshold of the adimensional parameter (ΔA/B)c under which transitions to the CSO state are observed as a function of density, we return to the condition of metastability defined in eqn (3) by analyzing ΔflocCs in the situation when K is aggregated into clusters of radius RK = rloc and Cs is still dispersed in the lattice. For mathematical simplicity, we again assume the clusters are circular, and that the unit area is a lattice site. In this configuration (flocCs)Cs ≈ 1, and the K particles with the lowest energy barriers in the system are at the boundary of the cluster. The local fraction of Cs, (flocCs)K, is by definition:
(4) |
NlocK = RK2(θRK − sinθRKcosθRK) + rloc2(θrloc − sinθrloccosθrloc) | (5) |
NlocCs = AlocCsρlocCs | (6) |
(7) |
Using eqn (4)–(7), we can calculate (ΔA/B)c as a function of fundamental parameters of the model ΔA, B, f0Cs, f0K, ρ, and rloc. The dashed black line in Fig. 6b shows the relationship (ΔA/B)cvs. ρ, for the simulation conditions rloc = 4 and f0Cs = f0K = 0.5, as well as the results from the simulations for different densities. This result suggests that the transition to the CSO state is facilitated at lower occupancy densities.
Footnote |
† Electronic supplementary information (ESI) available. See DOI: 10.1039/c7sc03524a |
This journal is © The Royal Society of Chemistry 2018 |