We develop a simple model to probe the ion transport and mechanical properties of low volume fraction colloidal nanoparticle gels. Specifically, we study the influence of the morphology of gels on ion diffusion and the corresponding roles of affinity to and enhanced ion transport along nanoparticle surfaces. We employ kinetic Monte Carlo simulations to simulate ion transport in the colloidal gels, and we perform nonequilibrium molecular dynamics to study their viscoelastic behavior. Our results indicate that in the presence of enhanced diffusion pathways for ions along the particle surface, morphology has a significant influence on the diffusivity of ions. We demonstrate that some gel morphologies can exhibit simultaneously enhanced ion transport and mechanical properties, thus illustrating a strategy to decouple ion transport and mechanical strength in electrolytes.

Liquid electrolytes are attractive components to enable various energy storage technologies such as lithium-ion batteries, fuel cells, and supercapacitors.1–3 However, safety and stability issues still limit their application in commercial devices. Solid-state polymer electrolytes and mixed-phase gel polymer electrolytes, which often possess low volatility, higher thermal stability, and mechanical strength, have emerged as possible alternatives.4–7 Yet, the latter face the challenge of the trade-off between ionic conductivity and mechanical strength, in which increases in the latter are accompanied by concomitant decreases in the former.1,8 As a consequence, there is interest in strategies that enhance both conductivity and mechanical integrity.

In the above context, some experiments have demonstrated that the addition of ceramic nanoparticle fillers to polymer electrolytes can simultaneously improve the mechanical strength and ionic conductivity.8–19 In amorphous polymers, the resulting ionic mobilities were suggested to be correlated with the modification of polymer conformations and segmental dynamics introduced by the nanoparticle fillers.20,21 By contrast, for semicrystalline polymers, it has been suggested that adding nanoparticle fillers reduces the crystallinity of polymers leading to an increase in the ionic conductivity.6,17,22 Other studies have also suggested that polymer-nanocomposite electrolytes promote conduction pathways along the filler/polymer interface as a result of interaction of surface groups on ceramic particles with the polymer chains and the ionic species.6,11,17,23–23 For instance, Maranas and co-workers24 reported higher ionic conductivities by biasing the nanoparticle surface with more acidic surface sites, indicating a greater degree of decoupling of ion mobility from polymer dynamics due to the resulting enhanced surface conduction pathways. Similarly, Yang et al.25 highlighted preferential diffusion of lithium ions at the filler/polymer interface as well as importance of morphology of fillers in the matrix. An electron hopping mechanism along the surface of TiO2 nanoparticles dispersed in an agarose gel electrolyte, similar to Grotthuss transport, was proposed by Shih et al.26 to explain enhanced ionic conductivity due to fillers.

Despite the reports of enhanced conductivities in polymer-nanoparticle composites, the systems explored have been in the regime of low nanoparticle volume fraction where the enhancement has primarily been attributed to changes in the structure and dynamics of the polymer component. Colloidal gels with percolated nanoparticle networks at higher nanoparticle loadings offer opportunities to take further advantage of conductivity increases from enhanced surface diffusion via the ability to modify the connected pathways for conduction via the morphology of the colloidal particles.

Colloidal gel electrolytes with sufficient nanoparticle loadings can also display attractive mechanical characteristics. Some of these gel electrolytes probed for ionic conductivity include fumed silica and/or colloidal silica in sulphuric acid,27–32 nanoparticle salts in a solvent,33 a nonpolymeric colloidal suspension of quaternary ammonium salt and ethylene glycol with sol-gel properties,34 and a yttria-stabilized zirconia (YSZ) sol-gel composite.35 Since colloidal particles are usually nonconducting, the resulting gel electrolyte is a crowded medium in which the colloidal particles hinder ion transport. Hence, a design goal might be to assemble percolated, porous gel microstructures of particles that are at a low enough volume fraction to effectively utilize the advantage of the surface diffusion pathways while still preserving favorable mechanical properties.

Elongated, space-spanning assemblies have been achieved by balancing short-range attraction and long-ranged repulsive interaction between colloids. For instance, Campbell et al.36 reported a connected network of chainlike clusters on introducing a nonadsorbing polymer due to the interplay of short-ranged attraction and long-ranged electrostatic repulsion between colloids in which the nonadsorbing polymer was introduced to generate short-ranged depletion attractions. Sacanna and co-workers used a combination of short-range attractions from polymeric depletants and electrostatic repulsive interactions from charged colloidal particles to assemble chains of anisotropic shaped colloidal particles.37 Similar experimental38–40 and computational41–44 studies have demonstrated that colloids interacting by the SALR (short-range attraction and long-range repulsion) class of pair potentials can lead to elongated cluster morphologies even at a low volume fraction of particles. Varying the range of interaction of the repulsive potential in such systems has shown to influence the morphology of the percolated particle network.45 Recently, a wide range of stringlike as well as compact cluster particle morphologies were realized in simulation via manipulating the ranges of competing attractive and repulsive pair potential components.46 

Inspired by the above developments, in this work, we consider a simple model for ion transport in systems of nanoparticles interacting via an isotropic pair potential described in Sec. II A. This pair potential has been demonstrated to drive spontaneous assembly of open (i.e., porous) microstructures including percolated strings (PS), clusters, and Bernal spiral configurations (BSP) at low particle packing fractions.46 Considering the presence of interface-assisted pathways in the resulting structures, we employ on-lattice kinetic Monte Carlo (KMC) simulations to a coarse-grained representation of the morphologies to probe the influence of different system parameters including nanoparticle volume fraction, surface diffusivity of ions (relative to bulk), morphology of the particles, and the enthalpic affinity of the ions to the particle surface. We further study the rheological behavior of these model gel configurations in their linear regime to characterize their mechanical properties. Our results demonstrate that such gel configurations have the potential to lead to simultaneously enhanced ion conductivities and improved mechanical properties.

The rest of the paper is organized as follows. Section II outlines the model employed to assemble different morphologies along with a description of methodologies used to simulate ion transport and the rheological behavior in the assemblies. The KMC results for ion diffusion in the presence of enhanced surface transport and the results for rheological characterization are presented in Sec. III. A concluding summary and an outlook for future research directions are provided in Sec. IV.

We consider a system of N identical particles of diameter σ in a cubic box of size L with periodic boundary conditions. The various particle morphologies probed in this study are generated by molecular dynamics simulations using the following pair potential U(r) proposed by Banerjee et al.46 based on its similarity to those designed by inverse methods to form stringy colloidal structures:

(1)

where β=(kBT)1, kB is the Boltzmann constant, T is temperature, εwca = 1.5, α = 12, εR = 6.76, w1 = (1 + δA)σ, and w2 = (1 + δA + δR)σ. System sizes and simulation methods employed for the morphology generation are as outlined in the original study.46 Modulating the range of attraction (δA) and repulsion (δR) at fixed particle volume fraction ϕ = Nπσ3/6L3 generates a variety of different particle structures. To study ion transport through these different geometries considering surface diffusion pathways, we focus on the following three different morphologies:

  • Percolated strings (PS): A percolated, porous network of single-stranded chains.

  • Clusters (C): A porous, crystalline, and unpercolated structure comprising compact tetrameric clusters.

  • Bernal spiral configurations (BSP): A percolated structure of compact, parallel strands with each strand comprising three intertwined helical chains.

The parameter choices for the pair potential where these different structures form in the molecular dynamics simulations are shown in Fig. 1(a).

FIG. 1.

(a) Morphologies obtained by varying δA and δR at ϕ = 0.157 for (b) the model potential given by Eq. (1) for δR = 0.7, 1.0, and 2.0 at fixed δA = 0.2. Snapshots of particle configurations shown for (c) percolated strings, (d) clusters, and (e) Bernal spiral.

FIG. 1.

(a) Morphologies obtained by varying δA and δR at ϕ = 0.157 for (b) the model potential given by Eq. (1) for δR = 0.7, 1.0, and 2.0 at fixed δA = 0.2. Snapshots of particle configurations shown for (c) percolated strings, (d) clusters, and (e) Bernal spiral.

Close modal

Probing ion transport considering surface diffusion pathways and characterizing the viscoelastic properties for every state point corresponding to the above structures is a computationally arduous task. Therefore, to gain qualitative insights into the transport and mechanical properties of different structures, we restrict our study to three points in the morphological phase diagram, each corresponding to a distinct class of structures. The parameter combinations considered are (1) δA = 0.2, δR = 0.7, corresponding to percolated strings, (2) δA = 0.2, δR = 1, corresponding to clusters, and (3) δA = 0.2, δR = 2, corresponding to a Bernal spiral configuration. Figure 1(b) shows the interaction potentials for the above points, while the snapshots of the corresponding structures are depicted in Figs. 1(c)–1(e).

To study the influence of the surface diffusion pathways arising from nanoparticles on ionic conductivities, we adopt a coarse-grained (i.e., on-lattice) ion-diffusion model. In this representation, each colloid in the morphology occupies 179 sites on the lattice with its size and shape being as shown in Fig. 2 to roughly emulate a sphere on a cubic lattice (similar to the study of Vilaseca et al.47). Since the volume occupied by each particle in this cubic lattice representation is not exactly that of a sphere, the packing fraction of particles in the lattice representation scales by a factor of 0.643 relative to that of the off-lattice representation of the spherical particles. The volume fraction of particles reported in Sec. III A of this paper corresponds to the volume fraction in the lattice representation, while that in Sec. III B corresponds to that of the actual spherical particles, although they both represent the same number of particles in a fixed size box.

FIG. 2.

On-lattice mapping of colloidal particles arranged in a simple cubic lattice with tracer ions (green dots) placed randomly on the vacant sites. Yellow lines indicate the transport pathways along the surface of particles, blue points represent surface sites, and purple solid cubes represent blocked sites corresponding to a single particle as shown in the enlarged image.

FIG. 2.

On-lattice mapping of colloidal particles arranged in a simple cubic lattice with tracer ions (green dots) placed randomly on the vacant sites. Yellow lines indicate the transport pathways along the surface of particles, blue points represent surface sites, and purple solid cubes represent blocked sites corresponding to a single particle as shown in the enlarged image.

Close modal

This coarse-grained, on-lattice representation using periodic boundary conditions ensures that each of the N particles from the configurations generated in Sec. II C is mapped appropriately without overlap. Sites not occupied by particles (i.e., vacant sites) are considered to be part of an implicit solvent matrix. The vacant sites adjacent to outer sites of the particle as shown in Fig. 2 are referred to as surface sites in this study, while remaining vacant sites are referred to as bulk sites. The ions are modeled as tracers hopping on these vacant sites with a specified rate which depends on the spatial location relative to the particle. In this study, we assume a noninteracting infinitely dilute solution of ions and neglect ion-ion correlation or ion-pairing effects. We assumed a grid spacing (Δl) of length unity in the lattice. In each step, the tracer ions execute a unit hop to its neighboring site. We consider the case in which the particles are nonconducting, and hence, the tracer ion cannot enter the sites occupied by a particle. We do note that the relative size of the tracer and obstacle has been shown to influence the volume fraction dependence of tracer ion diffusivity.47,48 As expected, such effects become less important when the obstacles are much larger than the tracers. To minimize such effects and mimic the large ratio of particle to ion size, we consider point tracer ions. Furthermore, since the time scale for diffusion of colloidal particles is typically long compared to ions or polymers, we consider the colloidal particles to be immobile.

Kinetic Monte Carlo (KMC) simulations have been extensively used to model a variety of rate processes including mobility of defect vacancies,49–54 polymerization kinetics,55,56 mechanisms of catalytic activity,57–61 crystal growth,62 and traffic jams63 to mention a few. Furthermore, on-lattice KMC simulations have been utilized to simulate longer time scales than those accessible by more detailed dynamics simulations and hence capture behavior relevant to continuum models.51,64–67 Such simulations are based on a master equation68 in which the probability distribution P(σi, t) for being in state σi at time t evolves as follows:

(2)

where kij denotes the hopping rate from state σi to state σj. In the context of our system, a transition or a move corresponds to the hop of a tracer ion on the lattice sites. Thus, the set containing all possible states at any time t includes every surface and bulk site hosting an ion at that particular time. Since the possible states in our system to host an ion correspond to only surface sites and bulk sites, there exist four different processes with unique hopping rates:

  1. surface to surface hopping rate (kSS)

  2. bulk to bulk hopping rate (kBB)

  3. surface to bulk hopping rate (kSB)

  4. bulk to surface hopping rate (kBS).

We assume the colloidal particles to be nonconducting and that the ions cannot enter any particle site, thus making the hopping rate from the surface site to any site occupied by a particle to be zero. In some applications, Density Functional Theory (DFT) calculations are performed to probe the potential energy landscape and use them to deduce various microscopic rates kij, which are then used as input to KMC simulations.49,54,57–61,69–72 Here, we take a more general approach which aims to elucidate the qualitative influence of various rates of ion hopping in the various morphologies. Hence, we assume parametric values of kij to highlight the influence of different transport processes.

In this work, we implemented the direct KMC method,73–75 one of the widely used rejection-free algorithms to solve the master equation [Eq. (2)]. The crux of this method lies in accepting every selected process/move for a particular ion but advancing the time appropriately to reflect the different transport rates. The probability of selecting a process j at site i (pij) is proportional to the hopping rates of process j at site i and is

(3)

where ki,tot is the sum of all possible hopping rates (to a neighboring site in the ±x, ±y, and ±z direction) for an ion at site i and ktot is the total sum of hopping rates for the ions in the entire lattice. For every move, the system clock is then advanced by

(4)

where ζ is an uniformly distributed random number. On an average, the system clock is therefore advanced by ktot1 for every move.

The simulation box consists of ions placed randomly on vacant sites at t = 0. Every KMC simulation run lasted for 5 × 107 moves with the above steps being executed at each move. Using the resulting trajectories of the ions, we extract the mobility of the ions using the Einstein-Smoluchowski relation

(5)

where (ri(t)ri(0))2 is the mean squared displacement (MSD) of tracer ions at time t and D is the diffusion coefficient of the tracer ions. To emphasize the simultaneous effect of obstruction by colloidal particles in the presence of an enhanced transport rate along their interface, we define a dimensionless number, D* = D/Do, where Do is the diffusion coefficient of ions in the system without any colloidal particles. The values for D* reported in this study are averaged over 30 KMC simulation runs for each system with 1024 noninteracting tracer ions in the system for every run.

In this study, we probe the effect of the following parameters on tracer ion diffusion: (i) particle volume fraction (ϕ), (ii) the rate of surface to surface hopping relative to that in the bulk (kSS/kBB), (iii) morphology of colloidal particles, and (iv) relative affinity of a tracer ion to surface sites compared to bulk sites. The results for this section are organized as follows: In Sec. III A 1, we present the results for the obstruction effect of colloidal particles on the tracer ion diffusion coefficient. In Sec. III A 2, the effect of introducing a higher value of kSS relative to kBB is summarized. Results presented in Sec. III A 3 characterize the significance of spatial arrangement of colloidal particles. In Sec. III A 4, we present the influence of an enthalpic affinity between the tracer ions and the surface sites.

In addition to the ion transport properties, we also compare viscoelastic properties of different gel structures in their linear regime. As a simplistic assumption, hydrodynamic effects are not accounted and only interparticle interactions are considered to explicitly identify the importance of particle networks and dynamics on viscoelastic properties. For this purpose, we use nonequilibrium molecular dynamics.

To carry out the above simulations, we start from a random initial configuration of particles and equilibrate the system for 3 × 108 steps using a box size of L = 30σ in the canonical ensemble at kBT/ε = 1 with time step dt=0.001ε/mσ2. Temperature is maintained using a Nosé-Hoover thermostat.76,77 After equilibration, the simulation box is deformed using the SLLOD equation of motion78 by applying Small Amplitude Oscillatory Shear (SAOS) strain γ with frequency ω, given by

(6)

The shear stress (σxy) response is of the same frequency but leads the shear strain function by angle δ, 0 ≤ δ ≤ 90°, given by

(7)

We calculate dynamic moduli, namely, (a) G′ (storage modulus) which determines the elastic characteristics and (b) G″ (loss modulus) which describes the viscous (liquidlike) behavior. We span a range of frequencies to calculate G′(ω) and G″(ω)

(8)
(9)

All the molecular dynamics simulations are executed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS)79 molecular dynamics package. Section III B reports the results for linear rheology of different gel structures considered in this study.

1. Particle obstruction effect

In this section, we present results for the long-time diffusivity of tracer ions as a function of particle volume fraction for the case where the surface hopping rates are identical to the bulk rates. Such a case isolates the obstruction effect on the diffusivity of ions arising from the presence of nonconducting particles. Previous computational studies have examined similar physics, and hence, we limit our discussion to only a few representative results to place them in the context of the results reported in Secs. III A 2–III A 4 of this paper. To probe such effects, we set the hopping rates to be unity, i.e., kBB = 1, kSS = 1, kBS = 1, and kSB = 1, to mimic random walk of tracers in obstructed media.

For the case of random arrangement of obstacles (i.e., equilibrium hard-sphere configurations), the results for dimensionless ion diffusivity D* are displayed in Fig. 3. We observe, as expected, that the effective ion diffusivity D* decreases with increasing ϕ. Since we are concerned with low volume fraction particle networks below their percolation threshold packing fraction, the long-time diffusion coefficient trend is found to be almost linear and matches well with the results reported in a similar study by Vilaseca et al.47 under similar conditions.

FIG. 3.

Circles denote dimensionless long-time diffusivity D* of tracer ions in configurations of hard particles (i.e., obstacles) at different packing fractions ϕ as computed from KMC simulations with equal surface and bulk hopping rates. Squares represent values reported in a similar study by Vilaseca et al.47 Dotted lines denote a linear fit for D*.

FIG. 3.

Circles denote dimensionless long-time diffusivity D* of tracer ions in configurations of hard particles (i.e., obstacles) at different packing fractions ϕ as computed from KMC simulations with equal surface and bulk hopping rates. Squares represent values reported in a similar study by Vilaseca et al.47 Dotted lines denote a linear fit for D*.

Close modal

2. Effect of enhanced surface hopping rate

Next, we present results for the long-time diffusion coefficient of tracer ions in a simple cubic lattice configuration of particles while accommodating enhanced transport along the nanoparticle-solvent interface. To probe the effect of enhanced surface transport, we vary the value of surface to surface hopping rate (kSS) while maintaining all other hopping rates to be unity. Since we limit our analysis to enhancement in surface transport (relative to the bulk), we report results only for the case of kSS ≥ 1.

From the results displayed in Fig. 4, for kSS > 1, we observe that the diffusion coefficients exhibit two classes of behaviors. In limit of high kSS, the ion diffusivities increase with increasing ϕ. By contrast, for low kSS, the ion diffusivities are seen to decrease with increasing ϕ. Such trends can be understood to arise from the competition between two factors which manifests with increasing particle volume fraction: (i) the greater obstruction effect from the impenetrability of ions into the nonconducting nanoparticles vs (ii) increased surface sites which present a faster hopping rate for ions. For relatively smaller values of kSS (kSS < 1.4), we see a decreasing trend for D* with respect to ϕ, reflecting the importance of the obstruction effect. By contrast, for higher values of kSS, a monotonically increasing trend is observed over the range of volume fraction considered due to the dominance of enhanced diffusion along the surface relative to the obstruction effect. Since the number of surface sites (relative to its volume) depends on the particle size, we also compare the dimensionless D* for two particle different sizes (Fig. 4). On increasing the particle size at a given volume fraction with exactly the same morphology and mesh size (i.e., effective ion size), the diffusion coefficient is seen to decrease. Such results can be understood as a consequence of the corresponding decrease in the number of ion-accessible surface sites per volume.

FIG. 4.

Dimensionless long-time diffusivity D* of tracer ions as a function of particle volume fraction ϕ with the enhanced hopping rate along particle surfaces. Solid symbols denote D* for the particle size considered in Fig. 2 (occupying 179 lattice sites), and open symbols represent D* for larger particles (occupying 493 sites with similar shape).

FIG. 4.

Dimensionless long-time diffusivity D* of tracer ions as a function of particle volume fraction ϕ with the enhanced hopping rate along particle surfaces. Solid symbols denote D* for the particle size considered in Fig. 2 (occupying 179 lattice sites), and open symbols represent D* for larger particles (occupying 493 sites with similar shape).

Close modal

3. Spatial configuration of colloidal particles

In this section, we present the results for the influence of morphology of colloidal particles on tracer ion diffusivity. The morphologies considered here include ordered structures, random arrangements of particles, and three different gel configurations obtained using the model discussed in Sec. II A. The gel configurations are generated by equilibrium molecular dynamics simulations as discussed in Sec. II C.

In Fig. 5(a), we present the results for ion diffusivities in different morphologies in the absence of enhanced hopping effects (all hopping rates equal to unity). We observe that for a specified volume fraction of particles, the ordered distribution of particles leads to a higher ion diffusivity compared to the random, disordered distribution. Such trends are in agreement with the results of Vilaseca et al.47 This trend can be explained by the fact that the ordered distribution of particles allows for a continuum path for tracers as opposed to random arrangement where the tracer might get trapped inside the pores spanned by the network of particles. Interestingly, we observe that the ion diffusivities for gel-like structures are also higher than those of random structures. Such results can be explained by the fact that gel configurations comprise an open, porous network of particles with more likelihood of connected pathways for ion transport.

FIG. 5.

(a) Dimensionless long-time diffusivities D* for tracer ions in different particle morphologies with kSS = 1, (b) D* for tracer ions in ordered and gel structures with kSS = 10 and kSS = 100, (c) D* for tracer ions in gel structures with kSS = 5, 50, and 500, and (d) MSDs for tracer ions in gel structures with restricted transport along the surface of particles at ϕ = 0.101 and ϕ = 0.1347.

FIG. 5.

(a) Dimensionless long-time diffusivities D* for tracer ions in different particle morphologies with kSS = 1, (b) D* for tracer ions in ordered and gel structures with kSS = 10 and kSS = 100, (c) D* for tracer ions in gel structures with kSS = 5, 50, and 500, and (d) MSDs for tracer ions in gel structures with restricted transport along the surface of particles at ϕ = 0.101 and ϕ = 0.1347.

Close modal

When the surface to surface hopping rate (kSS) is larger than the other hopping rates, we expect a stronger dependence of ion diffusivities on particle morphology. To probe such behavior, we fix the value of all hopping rates except kSS to be unity. As shown in Fig. 5(b), at relatively low values of kSS = 10, the ordered and gel structures exhibit very similar ion diffusion characteristics. By contrast, we observe clear differences in the diffusivity values between ordered and connected gel networks at the higher value of kSS = 100. Explicitly, all three gel configurations are seen to result in higher mobility of tracer ions compared to ordered configurations. Such trends are attributed to the fact that gel networks allow for surfaces of neighboring particles to be connected, allowing the tracer ions to directly hop from the surface of a particle to the surface of its neighboring particle and span larger distances. By contrast, the ordered structures lack interconnectivity of particles at a low volume fraction and cannot utilize the advantage of connected surface pathways.

Interestingly, among the different gel configurations, a difference in the ion diffusivities is seen with certain morphologies exhibiting higher diffusivities than others at the same volume fraction. As shown in Figs. 5(b) and 5(c), at higher values of the surface hopping rate corresponding to kSS = 50, 100, and 500, the order for D* in the volume fraction range, 0.091 < ϕ < 0.101 (where all the three gel configurations are observed), is Bernal spiral > percolated strings > clusters. By contrast, the trend for D* at ϕ = 0.1347 for kSS = 50, 100, and 500 is the following: percolated strings > Bernal spiral > clusters.

The above results can again be explained by invoking the extent of connectivity and percolation of the different morphologies. Explicitly, for kSS = 50, 100, and 500, we expect the transport of the tracer ion along the surface of the particle to be the dominant contributor to the overall diffusion coefficient, D*, in comparison to diffusion through the bulk. Within such a picture, we hypothesize that the different trends observed at different volume fractions can be rationalized based on the configurations which allow ions to span larger distances along the faster, connected surface pathways. To support such a hypothesis, we probed the MSDs of tracer ions for the case where they are restricted to hop just on the surface sites of particles and not allowed to jump back into the bulk. As shown in Fig. 5(d), the MSDs at ϕ = 0.101 for this case follow the order Bernal spiral > percolated strings > clusters, while at ϕ = 0.1347, the order is percolated strings > Bernal spiral > clusters. This trend for MSDs at these different volume fractions is exactly analogous to the trends observed for ion diffusivities (D*). Thus, we observe that the influence of the morphology on ion transport arises from the connectivity and percolation of surface pathways.

To provide a more physical explanation for the morphological effects observed above, we note that the coordination number for a perfect Bernal spiral configuration is six, which is maximum amongst the three gel networks considered in this study. Therefore, the Bernal spiral configuration of particles allows for tracer ions to easily hop from a surface site of a particle to the surface site of its neighbor and traverse larger distances per unit time due to the higher coordination number. Furthermore, each strand in the Bernal spiral configuration spans the entire length of the simulation box in one direction. By contrast, while the percolated strings similarly allow for a percolated 3D network, the average coordination number is only slightly larger than two, which is lower than that of the Bernal spiral configuration. The cluster configuration has higher coordination, but their lack of connectivity (i.e., percolated pathways) limits their conductivity.

As the volume fraction is increased from ϕ = 0.101 to ϕ = 0.1347, the particles are added to the pre-existing compact strands in the Bernal spiral configuration making them thicker and effectively exposing lesser surface sites as compared to the single-stranded percolated strings. As a consequence, we observe an inversion in the relative values of the diffusivity. This is reflected in the altered order between MSDs of percolated strings and Bernal spiral configurations at ϕ = 0.1347 as shown in Fig. 5(d) and consequently in overall ion diffusivities at higher values of kSS = 50, 100, and 500.

When the values of kSS > 1, but not too large (kSS = 5 and kSS = 10), the trends for D* are not in accordance with the surface-restricted MSDs shown in Fig. 5(d). For kSS = 5, the order for D* is percolated strings > clusters > Bernal spiral, while for kSS = 5, the order for D* is percolated strings > Bernal spiral > clusters. At such lower values of kSS, apart from the extent of connectivity and percolation, we expect significant contribution of bulk diffusion. Such characteristics are expected to depend on the distribution of pore sizes. In Fig. 6, we present the pore size distribution in different morphologies. It is seen that the average pore size for the Bernal spiral configuration is larger than that of percolated string and cluster morphologies. Thus, for an ion in the bulk, it is easier to find a surface site in its vicinity and utilize the faster diffusion pathways in percolated string and cluster morphologies than in the Bernal spiral configuration. Hence, contributions from both factors, namely, the pore size distribution and the enhanced transport along the particle surface, result in the above trend for D* for the different gel morphologies at relatively lower values of kSS = 5 and kSS = 10. Similarly, the slight differences in the D* values for different ordered structures at kSS = 10 and kSS = 100 as shown in Fig. 5(b) can be explained based on the effect of pore size distribution.

FIG. 6.

Pore-size distribution for different gel configurations at (a) ϕ = 0.101 and (b) ϕ = 0.135.

FIG. 6.

Pore-size distribution for different gel configurations at (a) ϕ = 0.101 and (b) ϕ = 0.135.

Close modal

In summary, when the transport rate along the particle surface is much higher in comparison to the transport rate in the bulk, the extent of connectivity and percolation of particles dictates the ion transport through different particle morphologies. However, we point out that that there exists no straightforward relation between the morphology and the ion diffusion properties. In this context, simulating surface-restricted diffusion of ions allows us to predict the performance of ion diffusion in different morphologies. By contrast, when the transport rate along the particle surface is not too large in comparison to the bulk transport rate, a significant effect of the average pore size in the network is observed along with the connectivity and percolation of particles.

4. Influence of affinity of ions to the particle surface

Since we are modeling ion transport in a solvent matrix, it is also of interest to examine the effect of affinity of ions (to the nanoparticle surface) on their transport properties. In realistic systems, such affinity differences arise from the distinct ion-polymer, ion-nanoparticle interactions. In this section, we probe the influence of preference to either particle surface or solvent matrix on the overall diffusivities in the presence of enhanced surface diffusion pathways.

To model the above physics, we assumed that the affinity to particle surface vs solvent matrix is quantified by a single parameter, ΔU, representing the difference in the affinity of the tracer ion to surface sites against that to the bulk sites. Based on the Boltzmann probability distribution of finding a tracer ion in these two states at thermodynamic equilibrium and the detailed balance criteria, we arrive at the following condition on the hopping rates (see Sec. S1 of the supplementary material):

(10)

We choose the value of kSB to be unity and fix the value of kBS based on the adopted bias, ΔU. For the results in this section, the volume fraction of particles in the lattice representation is fixed to be 0.101.

We probed two different values of the hopping rate along the surface, kSS = 10 and kSS = 100, while varying the affinity of tracer ions to the particle surface. The value of kBB is fixed to be unity. As shown in Fig. 7, when the tracer ions prefer bulk sites over the surface sites, i.e., ΔU > 1, the overall diffusion coefficients are similar for the different gel configurations. This result can be understood by noting that in such regimes, the ions exhibit an enthalpic preference to be in the bulk and hence do not exploit the faster transport pathways along the particle surface. However, upon biasing the preference of ions toward surface sites over the bulk sites, i.e., ΔU < 1, the trends for D* become dependent on the value of kSS. Explicitly, at comparatively low values of the surface hopping rate, kSS = 10, the value of D* is seen to decrease with increasing surface preference after a slight initial increase. For a larger value of kSS = 100, D* increases with an increase in the surface preference in all gel configurations except for clusters.

FIG. 7.

Long-time diffusivity D* of tracer ions as a function of relative stability of the ion on surface compared to bulk for (a) kSS = 10 and (b) kSS = 100. The results pertain to structures with particle volume fraction ϕ = 0.101.

FIG. 7.

Long-time diffusivity D* of tracer ions as a function of relative stability of the ion on surface compared to bulk for (a) kSS = 10 and (b) kSS = 100. The results pertain to structures with particle volume fraction ϕ = 0.101.

Close modal

The above trends can be explained by the interplay of two competing factors in the presence of bias toward (faster conducting) surface sites: (a) restricted transport on a 2D surface formed by connected surface sites which reduces the MSDs in comparison to that in a 3D space and (b) faster hops on the connected transport pathways along the particle surface in the gel networks. At comparatively small values of kSS, the restriction on moving over a 2D plane dominates the advantage of faster surface conduction pathways, resulting in lowered D* as observed in Fig. 7. However, at larger values of kSS, the relatively fast conduction along the surface becomes an important factor leading to an increase in diffusivities on increasing the preference toward surface sites.

On increasing the preference for surface over bulk, for all values of kSS, the order for D* is Bernal spiral > percolated strings > clusters at ϕ = 0.101. This is in agreement with the trends shown in Fig. 5(d). Such results highlight the fact that upon biasing the preference of tracer ions to be on the surface, the transport along the surface dominates the overall diffusion coefficient and effectively renders the contribution from the transport in bulk to become less important.

Together, the results discussed in this section demonstrate that on biasing the preference of ion toward surface sites, certain morphologies exhibit better diffusivities than others at a given volume fraction irrespective of the enhanced transport rate along the particle surface. The extent of connectivity and percolation of particles governs the ion diffusivity trend between different morphologies as discussed in Sec. III A 3. Furthermore, we observe that biasing the preference of the ion toward surface sites without a sufficiently higher transport rate along the particle surface (relative to bulk transport) is not advantageous.

Battery electrolytes are required to possess both high ionic conductivity and sufficient mechanical strength. While Sec. III A focused on ion transport, in this section, we characterized the rheological properties of the gel morphologies to probe whether the mechanical characteristics are decoupled from the ion transport properties in the presence of enhanced diffusion pathways along the particle surface. Since our results indicated that certain gels outperform others in the presence of enhanced surface transport, it is of interest to probe the rheological response to identify the importance of percolation and connectivity of particles on the mechanical strength. To quantify such features, we perform nonequilibrium molecular dynamics simulations (as discussed in Sec. II C) to probe the viscoelastic behavior of these gel networks.

In Fig. 8, we present the frequency sweep results for G′ and G at ϕ = 0.157 (corresponding to ϕ = 0.101 in lattice representations) for different gel configurations considered in this study. For Bernal spiral and cluster configurations, G′ is seen to exceed G″ over all the probed frequencies, suggesting a solidlike behavior. Percolated strings in contrast exhibit a viscoelastic response with a crossover frequency characterizing the transition from solidlike to liquidlike behavior.

FIG. 8.

SAOS frequency sweeps for different gel configurations. Solid symbols represent G′, and open symbols represent G″. The log-log slope (β) for G′ and G″ as a function of frequency(ω) is shown along the black dashed lines for percolated strings.

FIG. 8.

SAOS frequency sweeps for different gel configurations. Solid symbols represent G′, and open symbols represent G″. The log-log slope (β) for G′ and G″ as a function of frequency(ω) is shown along the black dashed lines for percolated strings.

Close modal

Recently, Sciortino and co-workers discussed the occurrence of arrested, ordered, and solidlike configurations of clusters at a low volume fraction with a specific SALR type potential.43 Such occurrences were attributed to the effective cluster-cluster repulsive interactions. For the SALR potentials considered in our study, we similarly hypothesize that for large enough repulsive ranges, i.e., of the order of particle diameter, the solidlike behavior observed in Bernal spiral and cluster configurations is a result of effective cluster-cluster interactions leading to arrested dynamics (see Sec. S2 of the supplementary material). By contrast, when the range of repulsive interactions is relatively small, as in the case of percolated strings, the cluster-cluster interactions are insignificant. Hence, percolated strings exhibit a transient, diffusive behavior (Fig. S1) with continuous reorganization of strings which rationalizes the viscoelastic response. Similarly, for the morphologies considered in this study, we observe that viscoelastic behavior is not related to percolation of particles in the network. The cluster configuration exhibits better mechanical properties than percolated strings despite the former structure lacking a percolating network. Such trends are in agreement with the observations reported in the literature wherein no direct correlation was shown to exist between percolation and dynamic arrest for colloidal gels interacting with SALR type potentials.41,80

For a fixed ω, the trend for elastic modulus (G′) is Bernal spiral > clusters > percolated strings as shown in Fig. 8. Thus, at ϕ = 0.157, the Bernal spiral configuration exhibits better mechanical strength along with optimal transport properties in the presence of surface diffusion pathways in comparison to other gel networks, illustrating a strategy to decouple ion transport and mechanical strength in colloidal gels.

To summarize, the cluster-cluster repulsive interactions seem to dictate the linear regime viscoelastic behavior of morphologies interacting with the SALR potential considered in this study. We observe no direct correlation between percolation of particles in the network and its viscoelastic behavior. From the morphologies considered in this study, we demonstrate that the mechanical properties of different structures are decoupled from their ion transport characteristics in the presence of enhanced diffusion pathways along the particle surface.

In summary, we studied the ion transport and mechanical properties of low volume fraction nanoparticle gel systems considering enhanced ion diffusion pathways along the particle surface. We considered a model isotropic potential to generate different gel microstructures at a low volume fraction of colloidal particles. The effects of parameters such as particle volume fraction, enhanced transport rate along the particle surface, morphology of particles, and enthalpic affinity to the particle surface on diffusivity of ions were probed. We probed the rheological characteristics of the gel networks in their linear viscoelastic regime.

In the limit of low volume fraction of particles (below the percolation threshold packing fraction), an increasing trend for diffusivity of ions with particle volume fraction was observed for higher transport rates along the particle surface. Thus, better diffusivity of ions was observed in comparison with that without the presence of nanoparticles. Furthermore, the diffusivity of ions increased with decreasing particle size. In the presence of enhanced transport along the particle surface, a stronger dependence of ion transport on the morphology of particles was observed. We observed that the extent of connectivity between the particles and network percolation significantly influences the diffusivity of ions at higher transport rates along the particle surface. The affinity of the tracer ion to the particle surface correlated with the diffusivity of ions at higher transport rates along the particle surfaces. Moreover, we observed that at ϕ = 0.101, the Bernal spiral configuration exhibits better mechanical properties along with optimal transport properties in the presence of enhanced transport along the particle surface in comparison with other gel networks.

The results presented in this work demonstrate the possibility of ion transport enhancement in colloidal gel electrolytes by realizing enhanced transport for ions along the colloidal surface. Furthermore, tuning the structure of particles might present a potential opportunity to enhance ion transport in colloidal gel electrolytes by exploiting the extent of percolation and connectivity between particles. Colloidal gel electrolytes provide an opportunity to decouple ion transport from mechanical strength in the presence of enhanced conduction pathways along the surface of colloidal particles by realizing a well connected network spanning large lengths.

See supplementary material for conditions on hopping rates to introduce affinity for surface/bulk sites and mean-squared displacements of particles for different morphologies during equilibration.

This work was primarily supported by the NSF through the Center for Dynamics and Control of Materials: an NSF Materials Research Science and Engineering Center (MRSEC) under Cooperative Agreement No. DMR-1720595. The authors acknowledge Texas Advanced Computing Center (TACC) for providing computing resources that have contributed to the research results reported within this paper. We also acknowledge the Welch Foundation (Grant Nos. F-1599 and F-1696) for support.

1.
K.
Xu
, “
Electrolytes and interphases in Li-ion batteries and beyond
,”
Chem. Rev.
114
,
11503
11618
(
2014
).
2.
M.
Armand
,
F.
Endres
,
D. R.
MacFarlane
,
H.
Ohno
, and
B.
Scrosati
, “
Ionic-liquid materials for the electrochemical challenges of the future
,”
Nat. Mater.
8
,
621
629
(
2009
).
3.
A.
Lewandowski
,
A.
Olejniczak
,
M.
Galinski
, and
I.
Stepniak
, “
Performance of carbon–carbon supercapacitors based on organic, aqueous and ionic liquid electrolytes
,”
J. Power Sources
195
,
5814
5819
(
2010
).
4.
M.
Armand
, “
Issues and challenges facing rechargeable lithium batteries
,”
Nature
414
,
359
367
(
2001
).
5.
Q.
Li
,
J.
Chen
,
L.
Fan
,
X.
Kong
, and
Y.
Lu
, “
Progress in electrolytes for rechargeable Li-based batteries and beyond
,”
Green Energy Environ.
1
,
18
42
(
2016
).
6.
J. L.
Schaefer
,
Y.
Lu
,
S. S.
Moganty
,
P.
Agarwal
,
N.
Jayaprakash
, and
L. A.
Archer
, “
Electrolytes for high-energy lithium batteries
,”
Appl. Nanosci.
2
,
91
109
(
2012
).
7.
J. Y.
Song
,
Y. Y.
Wang
, and
C. C.
Wan
, “
Review of gel-type polymer electrolytes for lithium-ion batteries
,”
J. Power Sources
77
,
183
197
(
1999
).
8.
E.
Quartarone
and
P.
Mustarelli
, “
Electrolytes for solid-state lithium rechargeable batteries: Recent advances and perspectives
,”
Chem. Soc. Rev.
40
,
2525
2540
(
2011
).
9.
W.
Krawiec
,
L.
Scanlon
,
J.
Fellner
,
R.
Vaia
,
S.
Vasudevan
, and
E.
Giannelis
, “
Polymer nanocomposites: A new strategy for synthesizing solid electrolytes for rechargeable lithium batteries
,”
J. Power Sources
54
,
310
315
(
1995
).
10.
X.
Huang
,
X.
Ma
,
R.
Wang
,
L.
Zhang
, and
Z.
Deng
, “
Combined effect of surface-charged latex nanoparticle AHPS and Al2O3 nano-fillers on electrochemical performance of the anionic gel polymer electrolytes PVA/P (MA-co-AHPS)
,”
Solid State Ionics
267
,
54
60
(
2014
).
11.
C.
Capiglia
,
P.
Mustarelli
,
E.
Quartarone
,
C.
Tomasi
, and
A.
Magistris
, “
Effects of nanoscale SiO on the thermal and transport properties 2 of solvent-free, poly(ethylene oxide) (PEO)-based polymer electrolytes
,”
Solid State Ionics
118
,
73
79
(
1999
).
12.
F.
Croce
,
S.
Sacchetti
, and
B.
Scrosati
, “
Advanced, high-performance composite polymer electrolytes for lithium batteries
,”
J. Power Sources
161
,
560
564
(
2006
).
13.
J. D.
Jeon
,
M. J.
Kim
, and
S. Y.
Kwak
, “
Effects of addition of TiO2 nanoparticles on mechanical properties and ionic conductivity of solvent-free polymer electrolytes based on porous P(VdF-HFP)/P(EO-EC) membranes
,”
J. Power Sources
162
,
1304
1311
(
2006
).
14.
F.
Croce
,
G. B.
Appetecchi
, and
B.
Scrosati
, “
Nanocomposite polymer electrolytes for lithium batteries
,”
Nature
394
,
456
458
(
1998
).
15.
M. M.
Jacob
,
E.
Hackett
, and
E. P.
Giannelis
, “
From nanocomposite to nanogel polymer electrolytes
,”
J. Mater. Chem.
13
,
1
5
(
2003
).
16.
P. G.
Bruce
,
B.
Scrosati
, and
J. M.
Tarascon
, “
Nanomaterials for rechargeable lithium batteries
,”
Angew. Chem., Int. Ed.
47
,
2930
2946
(
2008
).
17.
H. M.
Xiong
,
X.
Zhao
, and
J. S.
Chen
, “
New polymer-inorganic nanocomposites: PEO-ZnO and PEO-ZnO-LiClO4 films
,”
J. Phys. Chem. B
105
,
10169
10174
(
2001
).
18.
R.
Kumar
and
S. A.
Suthanthiraraj
, “
Segmental mobility and relaxation processes of Fe2O3 nanoparticle-loaded fast ionic transport nanocomposite gel polymer electrolyte
,”
J. Solid State Electrochem.
18
,
1647
1656
(
2014
).
19.
P. K.
Singh
,
B.
Bhattacharya
, and
R. K.
Nagarale
, “
Effect of nano-TiO2 dispersion on PEO polymer electrolyte property
,”
J. Appl. Polym. Sci.
118
,
2976
2980
(
2010
).
20.
B.
Hanson
,
V.
Pryamitsyn
, and
V.
Ganesan
, “
Mechanisms underlying ionic mobilities in nanocomposite polymer electrolytes
,”
ACS Macro Lett.
2
,
1001
1005
(
2013
).
21.
S.
Mogurampelly
,
V.
Sethuraman
,
V.
Pryamitsyn
, and
V.
Ganesan
, “
Influence of nanoparticle-ion and nanoparticle-polymer interactions on ion transport and viscoelastic properties of polymer electrolytes
,”
J. Chem. Phys.
144
,
154905
(
2016
).
22.
S. H.
Chung
,
Y.
Wang
,
L.
Persi
,
F.
Croce
,
S. G.
Greenbaum
,
B.
Scrosati
, and
E.
Plichta
, “
Enhancement of ion transport in polymer electrolytes by addition of nanoscale inorganic oxides
,”
J. Power Sources
97-98
,
644
648
(
2001
).
23.
O.
Borodin
,
G. D.
Smith
,
R.
Bandyopadhyaya
,
P.
Redfern
, and
L. A.
Curtiss
, “
Molecular dynamics study of nanocomposite polymer electrolyte based on poly(ethylene oxide)/LiBF4
,”
Modell. Simul. Mater. Sci. Eng.
12
,
S73
S89
(
2004
).
24.
L. V. N. R.
Ganapatibhotla
and
J. K.
Maranas
, “
Interplay of surface chemistry and ion content in nanoparticle-filled solid polymer electrolytes
,”
Macromolecules
47
,
3625
3634
(
2014
).
25.
T.
Yang
,
J.
Zheng
,
Q.
Cheng
,
Y. Y.
Hu
, and
C. K.
Chan
, “
Composite polymer electrolytes with Li7La3Zr2O12 garnet-type nanowires as ceramic fillers: Mechanism of conductivity enhancement and role of doping and morphology
,”
ACS Appl. Mater. Interfaces
9
,
21773
21780
(
2017
).
26.
C. M.
Shih
,
Y. L.
Wu
,
Y. C.
Wang
,
S. R.
Kumar
,
Y. L.
Tung
,
C. C.
Yang
, and
S. J.
Lue
, “
Ionic transport and interfacial interaction of iodide/iodine redox mechanism in agarose electrolyte containing colloidal titanium dioxide nanoparticles
,”
J. Photochem. Photobiol., A
356
,
565
572
(
2018
).
27.
M. P.
Vinod
,
K.
Vijayamohanan
, and
S. N.
Joshi
, “
Effect of silicate and phosphate additives on the kinetics of the oxygen evolution reaction in valve-regulated lead/acid batteries
,”
J. Power Sources
70
,
103
105
(
1998
).
28.
D. W. H.
Lambert
,
P. H. J.
Greenwood
, and
M. C.
Reed
, “
Advances in gelled-electrolyte technology for valve-regulated lead-acid batteries
,”
J. Power Sources
107
,
173
179
(
2002
).
29.
X.
Sun
,
J.
Zhao
,
T.
Chen
, and
X.
Liu
, “
Colloidal particle size of fumed silica dispersed in solution and the particle size effect on silica gelation and some electrochemical behaviour in gelled electrolyte
,”
J. Solid State Electrochem.
20
,
657
664
(
2016
).
30.
M. Q.
Chen
,
H. Y.
Chen
,
D.
Shu
,
A. J.
Li
, and
D. E.
Finlow
, “
Effects of preparation condition and particle size distribution on fumed silica gel valve-regulated lead-acid batteries performance
,”
J. Power Sources
181
,
161
171
(
2008
).
31.
Y. A.
Sang
,
D. J.
Euh
,
S. W.
Mi
, and
B. S.
Yoon
, “
The optimization of gel electrolytes on performance of valve regulated lead acid batteries
,”
Bull. Korean Chem. Soc.
29
,
998
1002
(
2008
).
32.
K.
Pan
,
G.
Shi
,
A.
Li
,
H.
Li
,
R.
Zhao
,
F.
Wang
,
W.
Zhang
,
Q.
Chen
,
H.
Chen
,
Z.
Xiong
, and
D.
Finlow
, “
The performance of a silica-based mixed gel electrolyte in lead acid batteries
,”
J. Power Sources
209
,
262
268
(
2012
).
33.
J. L.
Schaefer
,
D. A.
Yanga
, and
L. A.
Archer
, “
High lithium transference number electrolytes via creation of 3-dimensional, charged, nanoporous networks from dense functionalized nanoparticle composites
,”
Chem. Mater.
25
,
834
839
(
2013
).
34.
A.
Renjith
and
V.
Lakshminarayanan
, “
A novel colloidal suspension of TBA+BF4–EG and its applications as a soft solid electrolyte
,”
RSC Adv.
5
,
87956
87962
(
2015
).
35.
E.
Courtin
,
P.
Boy
,
T.
Piquero
,
J.
Vulliet
,
N.
Poirot
, and
C.
Laberty-Robert
, “
A composite sol-gel process to prepare a YSZ electrolyte for solid oxide fuel cells
,”
J. Power Sources
206
,
77
83
(
2012
).
36.
A. I.
Campbell
,
V. J.
Anderson
,
J. S.
Van Duijneveldt
, and
P.
Bartlett
, “
Dynamical arrest in attractive colloids: The effect of long-range repulsion
,”
Phys. Rev. Lett.
94
,
208301
(
2005
); e-print arXiv:0412108 [cond-mat].
37.
S.
Sacanna
,
W. T. M.
Irvine
,
P. M.
Chaikin
, and
D. J.
Pine
, “
Lock and key colloids
,”
Nature
464
,
575
(
2010
).
38.
C. A.
Saez Cabezas
,
G. K.
Ong
,
R. B.
Jadrich
,
B. A.
Lindquist
,
A.
Agrawal
,
T. M.
Truskett
, and
D. J.
Milliron
, “
Gelation of plasmonic metal oxide nanocrystals by polymer-induced depletion attractions
,”
Proc. Natl. Acad. Sci. U. S. A.
115
,
8925
8930
(
2018
).
39.
B.
Ruzicka
,
L.
Zulian
, and
G.
Ruocco
, “
Ergodic to non-ergodic transition in low concentration, Laponite
,”
J. Phys.: Condens. Matter
16
,
S4993
(
2004
).
40.
R. P.
Sear
,
S. W.
Chung
,
G.
Markovich
,
W. M.
Gelbart
, and
J. R.
Heath
, “
Spontaneous patterning of quantum dots at the air-water interface
,”
Phys. Rev. E
59
,
R6255
R6258
(
1999
).
41.
D.
Richard
,
C. P.
Royall
, and
T.
Speck
, “
Communication: Is directed percolation in colloid-polymer mixtures linked to dynamic arrest?
,”
J. Chem. Phys.
148
,
241101
(
2018
).
42.
N. E.
Valadez-Pérez
,
R.
Castañeda-Priego
, and
Y.
Liu
, “
Percolation in colloidal systems with competing interactions: The role of long-range repulsion
,”
RSC Adv.
3
,
25110
25119
(
2013
).
43.
F.
Sciortino
,
S.
Mossa
,
E.
Zaccarelli
, and
P.
Tartaglia
, “
Equilibrium cluster phases and low-density arrested disordered states: The role of short-range attraction and long-range repulsion
,”
Phys. Rev. Lett.
93
,
055701
(
2004
).
44.
F.
Sciortino
,
P.
Tartaglia
, and
E.
Zaccarelli
, “
One-dimensional cluster growth and branching gels in colloidal systems with short-range depletion attraction and screened electrostatic repulsion
,”
J. Phys. Chem. B
109
,
21942
21953
(
2005
).
45.
S.
Mossa
,
F.
Sciortino
,
P.
Tartaglia
, and
E.
Zaccarelli
, “
Ground-state clusters for short-range attractive and long-range repulsive potentials
,”
Langmuir
20
,
10756
10763
(
2004
).
46.
D.
Banerjee
,
B. A.
Lindquist
,
R. B.
Jadrich
, and
T. M.
Truskett
, “
Assembly of particle strings via isotropic potentials
,”
J. Chem. Phys.
150
,
124903
(
2019
).
47.
E.
Vilaseca
,
A.
Isvoran
,
S.
Madurga
,
I.
Pastor
,
J. L.
Garcés
, and
F.
Mas
, “
New insights into diffusion in 3D crowded media by Monte Carlo simulations: Effect of size, mobility and spatial distribution of obstacles
,”
Phys. Chem. Chem. Phys.
13
,
7396
7407
(
2011
).
48.
P. A.
Netz
and
T.
Dorfmüller
, “
Computer simulation studies of diffusion in gels: Model structures
,”
J. Chem. Phys.
107
,
9221
9233
(
1997
).
49.
B.
Sadigh
,
T.
Lenosky
,
S.
Theiss
,
M.-J.
Caturla
,
T.
Diaz de la Rubia
, and
M.
Foad
, “
Mechanism of boron diffusion in silicon: An ab initio and kinetic Monte Carlo study
,”
Phys. Rev. Lett.
83
,
4341
4344
(
1999
).
50.
S. M.
Gordon
,
S. D.
Kenny
, and
R.
Smith
, “
Diffusion dynamics of defects in Fe and Fe-P systems
,”
Phys. Rev. B
72
,
214104
(
2005
).
51.
J.
Dai
,
J. M.
Kanter
,
S. S.
Kapur
,
W. D.
Seider
, and
T.
Sinno
, “
On-lattice kinetic Monte Carlo simulations of point defect aggregation in entropically influenced crystalline systems
,”
Phys. Rev. B
72
,
134102
(
2005
).
52.
T.
Noda
, “
Modeling of indium diffusion and end-of-range defects in silicon using a kinetic Monte Carlo simulation
,”
J. Appl. Phys.
94
,
6396
6400
(
2003
).
53.
N.
Soneda
and
T. D.
de la Rubia
, “
Defect production, annealing kinetics and damage evolution in α-Fe: An atomic-scale computer simulation
,”
Philos. Mag. A
78
,
995
1019
(
1998
).
54.
M.
Sun
,
Z.
Li
,
G.-Z.
Zhu
,
W.-Q.
Liu
,
S.-H.
Liu
, and
C.-Y.
Wang
, “
Diffusion in Ni-based single crystal superalloys with density functional theory and kinetic Monte Carlo method
,”
Commun. Comput. Phys.
20
,
603
618
(
2016
).
55.
T.
Tongtummachat
,
S.
Anantawaraskul
, and
J. B. P.
Soares
, “
Dynamic Monte Carlo simulation of olefin block copolymers (OBCs) produced via chain-shuttling polymerization: Effect of kinetic rate constants on chain microstructure
,”
Macromol. React. Eng.
12
,
1800021
(
2018
).
56.
H.
Gao
,
A.
Waechter
,
I. A.
Konstantinov
,
S. G.
Arturo
, and
L. J.
Broadbelt
, “
Application and comparison of derivative-free optimization algorithms to control and optimize free radical polymerization simulated using the kinetic Monte Carlo method
,”
Comput. Chem. Eng.
108
,
268
275
(
2018
).
57.
Z.
Lian
,
S.
Ali
,
T.
Liu
,
C.
Si
,
B.
Li
, and
D. S.
Su
, “
Revealing the Janus character of the coke precursor in the propane direct dehydrogenation on Pt catalysts from a kMC simulation
,”
ACS Catal.
8
,
4694
4704
(
2018
).
58.
T. B.
Rawal
,
S. R.
Acharya
,
S.
Hong
,
D.
Le
,
Y.
Tang
,
F. F.
Tao
, and
T. S.
Rahman
, “
High catalytic activity of Pd1/ZnO(1010) toward methanol partial oxidation: A DFT+KMC study
,”
ACS Catal.
8
,
5553
5569
(
2018
).
59.
Y.
Cai
,
S.
Zhang
,
Y.
Xue
,
J.
Jiang
, and
Z. X.
Chen
, “
Density functional and kinetic Monte Carlo study of Cu-catalyzed cross-dehydrogenative coupling reaction of thiazoles with THF
,”
J. Org. Chem.
81
,
1806
1812
(
2016
).
60.
Z.-J.
Zuo
,
X.-Y.
Gao
,
P.-D.
Han
,
S.-Z.
Liu
, and
W.
Huang
, “
Density functional theory (DFT) and kinetic Monte Carlo (KMC) study of the reaction mechanism of hydrogen production from methanol on ZnCu(111)
,”
J. Phys. Chem. C
120
,
27500
27508
(
2016
).
61.
H. J.
Chun
,
V.
Apaja
,
A.
Clayborne
,
K.
Honkala
, and
J.
Greeley
, “
Atomistic insights into nitrogen-cycle electrochemistry: A combined DFT and kinetic Monte Carlo analysis of NO electrochemical reduction on Pt(100)
,”
ACS Catal.
7
,
3869
3882
(
2017
).
62.
G. H.
Gil
, “
Computer models of crystal growth
,”
Science
208
,
355
(
1980
).
63.
Y.
Sun
, “
Kinetic Monte Carlo simulations of two-dimensional pedestrian flow models
,”
Physica A
505
,
836
847
(
2018
).
64.
C. C.
Battaile
,
D. J.
Srolovitz
, and
J. E.
Butler
, “
A kinetic Monte Carlo method for the atomic-scale simulation of chemical vapor deposition: Application to diamond
,”
J. Appl. Phys.
82
,
6293
6300
(
1997
).
65.
A.
Magna
,
S.
Coffa
, and
L.
Colombo
, “
A lattice kinetic Monte Carlo code for the description of vacancy diffusion and self-organization in Si
,”
Nucl. Instrum. Methods Phys. Res., Sect. B
148
,
262
267
(
1999
).
66.
U.
Von Toussaint
,
T.
Schwarz-Selinger
, and
K.
Schmid
, “
First-passage kinetic Monte Carlo on lattices: Hydrogen transport in lattices with traps
,”
J. Nucl. Mater.
463
,
1075
1079
(
2015
).
67.
Y.
Zhang
,
C.
Jiang
, and
X.
Bai
, “
Anisotropic hydrogen diffusion in α-Zr and Zircaloy predicted by accelerated kinetic Monte Carlo simulations
,”
Sci. Rep.
7
,
41033
(
2017
).
68.
N. G.
Van Kampen
,
Stochastic Processes in Physics and Chemistry
, North-Holland Personal Library (
Elsevier Science
,
2011
).
69.
S.
Hong
,
A.
Karim
,
T. S.
Rahman
,
K.
Jacobi
, and
G.
Ertl
, “
Selective oxidation of ammonia on RuO2(110): A combined DFT and KMC study
,”
J. Catal.
276
,
371
381
(
2010
).
70.
M. A.
Albao
and
A. A. B.
Padama
, “
CO adsorption on W(100) during temperature-programmed desorption: A combined density functional theory and kinetic Monte Carlo study
,”
Appl. Surf. Sci.
396
,
1282
1288
(
2017
).
71.
S. I.
Shah
,
S.
Hong
, and
T. S.
Rahman
, “
Combined density functional theory and kinetic Monte Carlo study of selective oxidation of NH3 on rutile RuO2 (110) at ambient pressures
,”
J. Phys. Chem. C
118
,
5226
5238
(
2014
).
72.
P.
Xiao
and
G.
Henkelman
, “
Kinetic Monte Carlo study of Li intercalation in LiFePO4
,”
ACS Nano
12
,
844
851
(
2018
).
73.
D. T.
Gillespie
, “
A general method for numerically simulating the stochastic time evolution of coupled chemical reactions
,”
J. Comput. Phys.
22
,
403
434
(
1976
).
74.
D. T.
Gillespie
, “
Exact stochastic simulation of coupled chemical reactions
,”
J. Phys. Chem.
81
,
2340
2361
(
1977
).
75.
A.
Chatterjee
and
D. G.
Vlachos
, “
An overview of spatial microscopic and accelerated kinetic Monte Carlo methods
,”
J. Comput.-Aided Mater. Des.
14
,
253
308
(
2007
).
76.
S.
Nosé
, “
A unified formulation of the constant temperature molecular dynamics methods
,”
J. Chem. Phys.
81
,
511
519
(
1984
).
77.
W. G.
Hoover
, “
Canonical dynamics: Equilibrium phase-space distributions
,”
Phys. Rev. A
31
,
1695
1697
(
1985
).
78.
D. J.
Evans
,
G. P.
Morriss
,
D. P.
Craig
, and
R.
McWeeny
,
Statistical Mechanics of Nonequilibrium Liquids, Theoretical Chemistry
(
Elsevier Science
,
2013
).
79.
S.
Plimpton
, “
Fast parallel algorithms for short–range molecular dynamics
,”surfaces. Moreover, we observed that
J. Comput. Phys.
117
,
1
19
(
1995
).
80.
E.
Zaccarelli
, “
Colloidal gels: Equilibrium and non-equilibrium routes
,”
J. Phys.: Condens. Matter
19
,
323101
(
2007
).

Supplementary Material