Shibananda
Das†
a,
Jonas
Riest†
b,
Roland G.
Winkler
a,
Gerhard
Gompper
a,
Jan K. G.
Dhont
b and
Gerhard
Nägele
b
aTheoretical Soft Matter and Biophysics, Institute for Advanced Simulation and Institute of Complex Systems, Forschungszentrum Jülich, 52425 Jülich, Germany. E-mail: r.winkler@fz-juelich.de; g.gompper@fz-juelich.de
bSoft Matter, Institute of Complex Systems, Forschungszentrum Jülich, 52425 Jülich, Germany. E-mail: j.g.k.dhont@fz-juelich.de; g.nägele@fz-juelich.de
First published on 27th November 2017
Dispersions of particles with short-range attractive and long-range repulsive interactions exhibit rich equilibrium microstructures and a complex phase behavior. We present theoretical and simulation results for structural and, in particular, short-time diffusion properties of a colloidal model system with such interactions, both in the dispersed-fluid and equilibrium-cluster phase regions. The particle interactions are described by a generalized Lennard-Jones-Yukawa pair potential. For the theoretical-analytical description, we apply the hybrid Beenakker–Mazur pairwise additivity (BM-PA) scheme. The static structure factor input to this scheme is calculated self-consistently using the Zerah-Hansen integral equation theory approach. In the simulations, a hybrid simulation method is adopted, combing molecular dynamics simulations of colloids with the multiparticle collision dynamics approach for the fluid, which fully captures hydrodynamic interactions. The comparison of our theoretical and simulation results confirms the high accuracy of the BM-PA scheme for dispersed-fluid phase systems. For particle attraction strengths exceeding a critical value, our simulations yield an equilibrium cluster phase. Calculations of the mean lifetime of the appearing clusters and the comparison with the analytical prediction of the dissociation time of an isolated particle pair reveal quantitative differences pointing to the importance of many-particle hydrodynamic interactions for the cluster dynamics. The cluster lifetime in the equilibrium-cluster phase increases far stronger with increasing attraction strength than that in the dispersed-fluid phase. Moreover, significant changes in the cluster shapes are observed in the course of time. Hence, an equilibrium-cluster dispersion cannot be treated dynamically as a system of permanent rigid bodies.
In contrast to the multitude of studies on static properties of SALR systems,9,17 comparably little is known about their dynamical aspects. A major challenge for simulation and theoretical studies of the dynamics of these systems is to account for solvent-mediated hydrodynamic particle interactions (HIs) which decisively influence diffusion and rheology. HIs are long-ranged, and, in general, non-pairwise additive for larger particle concentrations. In earlier studies,2 Brownian dynamics simulations based on the generalized Lennard-Jones–Yukawa SALR model potential have been performed, where the salient HIs were completely disregarded. In these simulation studies, the dynamic structure factor, S(q,t), as a function of correlation time t and wavenumber q, and the wavenumber-dependent short-time diffusion function, D(q), have been determined.
The theoretical description of cluster-phase states is further hampered by the occurrence of additional time and length scales associated with the distributions of particle-cluster lifetimes, cluster sizes and shapes, and cluster electric charges. This renders a clear distinction between colloidal short-time and long-time regimes complicated, in contrast to homogeneous suspensions of individually diffusing monodisperse particles. An interesting experimental finding exemplifying the intriguing dynamics of SALR systems was obtained by Liu et al. in ref. 4. The short- and long-time self-diffusion coefficients for salt-free lysozyme protein solutions, deduced from neutron spin echo (NSE) measurements, were found to exhibit roughly the same concentration dependence. Lysozyme solutions, however, are special SALR systems due to the non-isotropic short-range patchy attractions between the lysozyme proteins.18,19 In this context, by using mesoscale hydrodynamics simulations, the strong influence of particle surface patchiness on short-time diffusion properties was shown recently.20 The patchy SA-interaction contribution gives rise to smaller values of the cage diffusion coefficient D(qm), where qm is the location of the particles next neighbor peak in S(q), than those observed for isotropic attraction of comparable strength and range.
A comprehensive theoretical study of short-time dynamic properties of SALR systems in the dispersed-fluid phase has been presented,18,21 where the analytic hybrid Beenakker–Mazur (plank) pairwise-additivity (BM-PA) scheme for calculating D(q) has been used. This scheme requires the static structure factor, or likewise the associated radial distribution function, as the only input. The BM-PA method gives overall good results for the short-time diffusion properties, both of neutral and charge-stabilized particle dispersions, which was shown in comparison with elaborate hydrodynamic force-multipole simulations and experimental data on various colloidal systems.22–24 Interestingly, a non-monotonic dependence of the sedimentation velocity of SALR dispersions on the interaction strength parameter is predicted by this method,21 revealing the intricate influence of short-range attractive and long-range repulsive pair potential contributions on the short-time dynamics.
In the present paper, we study the short-time dynamics of SALR systems both in the monomer-dominated dispersed-fluid and the equilibrium-cluster phases. We assess the particle dynamics by employing both a mesoscale hydrodynamic simulation method, the multiparticle collision dynamics (MPC) approach, and the aforementioned analytic BM-PA scheme. MPC is a particle-based simulation technique which incorporates thermal fluctuations, provides hydrodynamic correlations,25 and is easily coupled with other simulation techniques such as molecular dynamics simulations for embedded particles.26–28 It has been successfully applied to study equilibrium and non-equilibrium dynamical properties of colloids,29–39 polymers,40–47 blood cells,48,49 and biological microswimmers.50–52 Specifically, MPC simulations account for the many-particle HIs in the overdamped Brownian dynamics regime of friction-dominated particle motion addressed in the present work.
From the comparison with our MPC simulation results, we demonstrate the high accuracy of the hybrid BM-PA scheme for systems in the dispersed-fluid phase. In addition, we explore the inter- and intra-cluster dynamics in the equilibrium-cluster phase by simulations, and determine the characteristic cluster lifetime also in comparison with analytic first-passage time calculations. The intra-cluster dynamics in the equilibrium-cluster phase of SALR systems is of special interest, and here in particular the influence of the HIs. If the particles forming a cluster can move individually, i.e., if the cluster is non-rigid, there is no hydrodynamic screening inside the cluster. In fact, it is well known that HIs between an ensemble of interacting mobile particles is not screened owing to momentum conservation,53,54 and this is shown quantitatively in a recent Stokesian dynamics simulation study.55 The absence of HIs screening in turn implies that a polydispersity model describing the SALR dispersion as a mixture of differently sized rigid clusters is not appropriate, at least regarding the system dynamics. Our simulations clearly reveal internally very dynamic clusters, and a cluster lifetime which depends significantly on HIs. This emphasizes the importance of HIs in the dynamics of the particles also in the equilibrium-cluster phase.
The paper is structured as follows. Section 2 comprises the presentation of the modified-LJY model system of spherical Brownian particles under SALR conditions, a short description of the analytic BM-PA approach, and of the MPC simulation method. Section 3 presents results for structural properties and the phase behavior. Dynamical aspects are discussed in Section 4. Summary and conclusions are given in Section 5.
(1) |
The phase behavior of the LJY model system for fixed potential parameter values A = 2, ξ = 1.794, and varying values of ε has been studied by MD simulations and the thermodynamic Gibbs–Duhem integration method.6 Thereby, the selected reduced screening length ξ corresponds, for a typical particle diameter of σ ≈ 100 nm, to a 1:1 aqueous electrolyte solution at a concentration of about 3 μM.6 The shape of the LJY potential for these parameters, and various attraction strengths ε, is shown in Fig. 1. The minimum of the potential is located at xmin = 1.014, where βV(xmin) ≈ 2 − ε (cf.Fig. 1). The range of attraction strengths considered in this work is ε = 0–8. The second zero-crossing x0 > xmin of V(x) is x0 ≈ 1.032 for ε = 3 and x0 ≈ 1.06 for ε = 8. Here, x0 characterizes the radial position of the comparably small interval xmin < x < xmax of short-range attraction, where V′(x) = 0 at xmax in the interval xmin < x < ξ.
Fig. 1 The generalized Lennard-Jones–Yukawa (LJY) potential of eqn (1) for various values of the reduced attraction strength ε, as indicated, and the fixed values A = 2 and ξ = 1.794 for the strength and range of the long-range Yukawa-repulsion part, respectively. |
The short-time dynamics of Brownian particles can be probed experimentally by measuring the q-dependent dynamic structure factor, S(q,t), using dynamic light scattering (DLS) or neutron spin echo spectroscopy (NSE), depending on the particle size, concentration, and other system properties. For short times, S(q,t) decays single exponentially, and it can be expressed in terms of the short-time diffusion coefficient D(q) as
S(q,t ≪ τD) = S(q)exp(−q2D(q)t). | (2) |
(3) |
(4) |
We assess here, by the comparison with elaborate MPC simulations, the applicability of the BM-PA method to SALR systems for the first time quantitatively. Notice that it is not clear a priori whether the BM-PA method applies to SALR systems, where intermediate-range ordering is accompanied by a larger likelihood of near-contact configurations than for purely repulsive particle systems.
(5) |
The dynamics of a point particle is described by Newton's equation of motion, which is solved via the velocity Verlet algorithm with the forces derived from the above potentials.68
The simulations of ref. 37 demonstrate that this model closely resembles a hard-sphere colloidal particle with no-slip boundary conditions in a dilute dispersion. In particular, the flow field determined for an individually sedimenting MPC-colloidal sphere agrees very well with the theoretical expression for an isolated colloidal hard sphere. However, due to discretization, the hydrodynamic radius is somewhat larger than the radius a. The actual value depends on the radius of the colloid and the size of the simulation box. Explicit values are provided in Section 2.3.4.
ri(t + h) = ri(t) + hvi(t), | (6) |
vi(t + h) = vcm(t) + R(α)(vi(t) − vcm(t)). | (7) |
(8) |
(9) |
To maintain the temperature of the system at a desired value, we use a collision-cell level thermostat, referred to as Maxwell–Boltzmann scaling (MBS) method.70 This ensures that a canonical ensemble is simulated.70,71
In the following, we will compare simulation results with analytical calculations in the overdamped time regime. This requires that the characteristic time scale for Brownian motion τB = NsM/γs, momentum diffusion (hydrodynamics) τν = a2/ν, and diffusion of a dispersion particle across its own radius τD = γsa2/kBT obey the relation τB ≲ τν ≪ τD,31 where γs = 6πη0a is the Stokes friction coefficient of a sphere, and ν = η0ã3/m〈Nc〉 is the kinematic viscosity. With the above parameters, these relations are satisfied, specifically τν/τD < 5 × 10−3. In the subsequent chapters, we will consider time scales t/τν > 1 only.
Fig. 2 Static structure factor, S(q), of LJY systems at reduced particle concentration ρ* = ρσ3 = 0.1, for the indicated interaction strengths ε. Simulation results are displayed by symbols and theoretical ZH integral equation predictions by solid lines. (a) MPC simulation and theoretical results of the dispersed-fluid phase for ε = {2, 4, 6}. (b) MPC simulation results for dispersed-fluid phase systems with ε = {2, 4, 6} (open symbols), and equilibrium-cluster phase systems with ε = {7, 8} (filled symbols). The horizontal dashed line in (b) marks the critical value Scrit(qc) ≈ 2.7 for a first-order transition according to Godfrin et al. in ref. 5. |
Typical microstructural snapshots from simulations are displayed in Fig. 3. We identify here two phases, namely a dispersed-fluid phase at small attraction strengths and an equilibrium-cluster phase at larger ε. The emerging clusters in the latter case are clearly visible.
Insight into the phase behavior of the LJY system is gained from the cluster-size distribution function (CSD) defined as
(10) |
We focus here on the microstructure and dynamics of the dispersed-fluid and equilibrium-cluster phases. In the dispersed-fluid phase, N(s) is monotonically decaying with increasing s, as shown in Fig. 4, which indicates a monomer-dominated dispersion, where only transient small clusters are formed, corresponding to a mean cluster size of = 1–3. In contrast, the distribution function in the equilibrium-cluster phase exhibits a second maximum (peak) at a particular cluster size s* > 2, in addition to the monomer peak at s = 1. The second local peak reflects equilibrium clusters with a preferred size around s*, coexisting in thermodynamic equilibrium with the monomers.
Fig. 4 Distribution of cluster sizes, N(s), obtained from simulation configurations as function of the cluster size s, for A = 2, ξ = 1.794, and the indicated attraction strengths ε. The concentration is ρ* = 0.1 and in the inset ρ* = 0.2. For dispersed-fluid phase systems, open symbols connected by dashed lines are used, while equilibrium-cluster phase systems are marked by filled symbols connected by solid lines. Additionally, the inset shows a typical equilibrium cluster for ε = 7 and ρ* = 0.1. Two dispersion particles are part of a cluster for distances r ≤ rcluster = xmaxσ, where xmax > 1 is the location of the first maximum of V(x) (cf.Fig. 1). |
From the shape of the CSD functions at ρ* = 0.1, the systems with ε = 2 and 4 are identified as belonging to the dispersed-fluid phase, while the systems with ε = 7 and 8 are part of the equilibrium-cluster phase. This classification is in accord with the ε–ρ* LJY-phase diagram of ref. 6 obtained for the same potential parameters (A = 2 and ξ = 1.794) (cf.Fig. 5), as well as the phase transition criterion S(qc) ≳ 2.7 of ref. 5. Note that for ε = 6, a shallow peak at s* > 1 appears, but the monomer peak at s = 1 still dominates. A close inspection of the simulation configurations suggests that the system for ε = 6 is still in the dispersed-fluid phase. This conclusion is supported by the schematic phase diagram of ref. 6 (see also Fig. 5, where the system with ρ* = 0.1 and ε = 6 is inside the dispersed-fluid phase region.)
Fig. 5 Schematic ε–ρ* phase diagram deduced from our MPC-generated cluster size distribution functions for the eight considered LJY systems. Open circles and filled circles label systems in the dispersed-fluid and equilibrium-cluster phases, respectively, according to the CSD shapes. The two uniformly colored areas mark the dispersed-fluid (red) and equilibrium-cluster (blue) phase regions of the LJY system as predicted in our simulation studies (see also ref. 6). |
Pair correlation functions, g(r), obtained from ZH calculations and simulations are compared in Fig. 6 for ρ* = 0.1. For the dispersed-fluid phase systems, excellent agreement between analytical and simulation results is observed (cf.Fig. 6(a)). Note that the simulation results show two additional small peaks at x ≈ 1.67 and 1.75 for ε = 6. This again reflects that this system is close to the boundary separating the dispersed-fluid form the equilibrium-cluster phase. The g(r)'s of the systems in the equilibrium-cluster phase exhibit sharp peaks at various distinct positions indicating strong local ordering (Fig. 6(b)). The positions of the peaks can be attributed to preferred geometric particle configurations such as an in-line configuration of three particles (peak at x ≈ 2), an octahedral arrangement of four particles (peak at x ≈ 1.66), formation of two equilateral triangles with a common side (peak at x ≈ 1.76),13 and a cubic arrangement (peak at x ≈ 1.43).
Fig. 6 Radial distribution functions, g(r), for the dispersions of Fig. 2 with ρ* = 0.1 and indicated ε values. (a) MPC (symbols) and ZH (solid lines) results for three systems in the dispersed-fluid phase. (b) MPC results for two systems in the equilibrium-cluster phase. The depicted cluster configurations and the arrows indicate the cluster structures associated with the respective peaks in g(r). |
Similarly, good agreement between ZH and MPC results for g(r) and S(q) is observed for the dispersed-fluid phase systems at the larger concentration ρ* = 0.2 and for ε ≤ 5, as displayed in Fig. 7. The system with ε = 6 belongs to the equilibrium-cluster phase. This follows from the respective N(s) shown in the inset of Fig. 4, and is in accord with the phase diagram of Mani et al. in ref. 6. The ZH results for S(q) and g(r) at ε = 6 significantly differ from the simulation data, in particular around the IRO peak position at qc ≲ 2.5. Interestingly, g(r) for the equilibrium-cluster phase system with ε = 6 has a broad peak at larger x that starts roughly at xc = 2π/qc and becomes more pronounced with increasing S(qc) (see inset of Fig. 7(b)). However, as it was shown recently,10 such a broad peak in g(r), visible here in the range 3 < xc < 5, does not necessarily imply an IRO peak in the corresponding S(q).
A relation was recently established between the IRO wavenumber qc, or likewise the average center-to-center distance of clusters estimated by 2π/qc, the mean particle number of clusters, and the dispersion particles concentration ρ*, which reads17
(11) |
ρ* | ε | 2π/(qcσ) | (/ρ*)1/3 |
---|---|---|---|
0.1 | 6 | 3.57 | 3.34 |
7 | 4.17 | 4.11 | |
8 | 4.17 | 4.31 | |
0.2 | 5 | 2.78 | 3.35 |
6 | 3.33 | 4.17 |
Following ref. 72, we have also analyzed the thermal width, ξT, of the principal (cluster) peak of S(q), obtained from a quadratic Ornstein-Zernike fit of S(qc)/S(q) around qc. Consistent with the phase transition criterion of ref. 72, for ε > 6 our ξT is larger than the (reduced) long-range repulsion length ξ = 1.794 in eqn (1), while it is smaller for ε < 6. This reconfirms the classification of our two systems to be in the cluster-fluid and dispersed-fluid phase, respectively.
Fig. 8 MPC results for the inverse of the short-time diffusion function, D(q), expressed in units of the single-particle diffusion coefficient d0, for LJY systems at ρ* = 0.1 and the indicated ε values. Open (filled) symbols label systems in the dispersed-fluid (equilibrium-cluster) phase. The dashed lines are simulation results without HIs considered of d0/D(q)|no-HI = S(q), for ε ∈ {2, 6, 8}. Inset: Comparison of MPC (symbols) with BM-PA scheme results (solid lines), for systems in the dispersed-fluid phase. The BM-PA scheme uses the ZH S(q) and g(r) as input (cf.Fig. 2(a) and 6(a)). |
For the dispersed-fluid phase systems, the BM-PA results for D(q) are in excellent agreement with the simulation results. Note that the system with ρ* = 0.2 and ε = 6 is in the equilibrium-cluster phase (cf.Fig. 4 and 7(a)), hence no quantitative agreement can be expected between the simulation and BM-PA results for D(q). The quantitative agreement for ε < 6 with the simulation results establishes the hybrid BM-PA scheme as an easy-to-apply, precise, and fast tool for assessing the short-time diffusion in SALR systems in the dispersed-fluid phase state. This has been additionally verified recently by the comparison of BM-PA predictions for the two-Yukawa SALR model with NSE measurements of D(q) of lysozyme solutions at different temperatures and concentrations.18 The pronounced slowing down of collective diffusion in SALR systems due to HIs is illustrated by the comparison of the MPC results for d0/D(q) with d0/D(q)|no-HI = S(q), where HIs are disregarded. According to Fig. 9, neglecting HIs results in an overestimation of D(q) by a factor of about 2.
(12) |
For comparison, we estimate the cluster lifetime analytically in an approximate manner by considering only the dissociation of an isolated pair of particles. The associated dissociation or binary escape time, τdc, is the mean time required for a pair initially at the minimum distance xmin to unbind by increasing its mutual distance to a value xd > xmax larger than the potential barrier position, xmax. Theoretically, this characteristic time is obtained by applying the absorbing boundary condition f(x = xd,t|xmin) = 0 for the conditional pair-probability density function f(x,t|xmin) of finding the pair at a distance x = r/σ at time t, given its distance was xmin at time t = 0. Using the two-particle Smoluchowski equation, τdc is given by75
(13) |
G(r) = X11(r) − X12(r), | (14) |
Fig. 10 shows the cluster lifetimes obtained from the simulations and the two-particle Smoluchowski theory. The dissociation time τdc is calculated from eqn (13) using xmin(ε), xd = xmax(ε), and xl = 1 (cf.Fig. 1). We find τdc to be insensitive to the precise choice of xl ≤ 1, owing to the steep near-contact repulsion of the LJY potential in eqn (1). The upper boundary value xb = xmax is consistent with the cluster cutoff rcluster used in the simulation calculations of τc.
Fig. 10 Cluster lifetime as a function of attraction strength ε for dispersions with ρ* = 0.1. The squares indicate MPC simulation results for τc. The solid line is the pair dissociation time, τdc, according to eqn (13). The dashed line represents τdc/5, and the dashed-dotted blue line is the first passage time, τKc, according to Kramer's theory eqn (15). |
The cluster lifetime τc increases with increasing interaction strength ε as expected. Quite interestingly, the simulation results in the dispersed-fluid regime, where ε ≤ 5, exhibit qualitatively the same ε dependence as the binary dissociation time, quantitatively, however, τc ≈ τdc/5 (see dashed black line). The simulation values for τc increase more strongly than τdc in the equilibrium-cluster range of larger ε values. Here, for a large potential barrier of height Vmax − Vmin, the simulation values for τc are close to the binary escape time τKc derived from Kramer's theory for barrier crossing. The latter is given by
(15) |
The time evolution of a cluster of seven dispersion particles obtained in a MPC simulation is illustrated in Fig. 11. Evidently, the cluster undergoes substantial conformational changes during the considered time span of t ≈ 19τD. The particle labeled in green that finally detaches from the cluster is initially well bound to several of the cluster particles; as time progresses, it moves to a less stable place, where it is bound to only two particles before it dissociates from the cluster after the time t ≈ 19τD. The time sequence in Fig. 11 illustrates the pronounced relative mobilities of dispersion particles forming an equilibrium cluster. We like to emphasize that, due to rigid constraints, the hydrodynamic mobility tensors in a (fictitious) dispersion of rigid clusters with the same equilibrium cluster distribution as the considered LJY system are distinctly different from those of the LJY system of non-rigid, thermodynamically formed clusters. Transport properties such as D(q) and the dispersion viscosity are consequently expected to be different for the two systems. Hence, the dynamics of a SALR system cannot be adequately described by treating it as a polydisperse dispersion of rigid Brownian cluster objects. A quantitative assessment of the differences in the transport properties for rigid-cluster and non-rigid equilibrium-cluster dispersions can be the subject of a future simulation study. Furthermore, by the nonlinear hydrodynamic coupling of LJY cluster particles with sensitive dependence on the initial configuration and the stochastic Brownian forces, different MPC realizations of a tagged cluster sequence with identical initial configurations can develop quite differently in the course of time.
Fig. 11 MPC generated time sequence of a cluster consisting of 7 particles for ε = 7 and ρ* = 0.1, tracked until a particle (labeled in green) dissociates. |
We have studied the structural properties via g(r) and S(q), and the equilibrium phase behavior for varying values of the attraction strength (potential well depth) ε and two particle concentrations ρ*. Our phase-state mapping based on alternative criteria invoking the cluster-size distribution function and the IRO peak height and width of S(q), is in accord with a previous study by Mani et al. in ref. 6 on the same model system. A phase transition from the dispersed-fluid to the equilibrium-cluster phase is observed with increasing attraction strength. In the dispersed-fluid phase, the self-consistent ZH scheme gives results for S(q) and g(r) in excellent agreement with MPC simulation data. Likewise for the dispersed-fluid phase, we have shown the good accuracy of the analytic BM-PA results regarding the short-time diffusion function, D(q), by the comparison with corresponding MPC simulation predictions. Thus, the BM-PA scheme allows for the convenient assessment of short-time diffusion properties for dispersed-fluid phase SALR systems.
Using MPC simulations, we have studied the dynamics of LJY systems in the equilibrium-cluster phase, where for the first time the long-range HIs are fully accounted for. The simulations yield for this phase significantly smaller D(q) values as for dispersed-fluid systems. To explore the intra-cluster dynamics, we calculated the dissociation (first passage) time, τdc, for an isolated pair of LJY particles, and compared it with the MPC cluster lifetime, τc, for concentrated dispersions. This study is amended by the visualization of the shape changes of a cluster up to the event when a particle is leaving the cluster. The simulations reveal a highly mobil internal cluster dynamics, with pronounced conformational changes within a time span of a few τD values. Thus, both, the intra-cluster dynamics and the dynamics of the dispersion as a whole are significantly influenced by many-particle HIs, which are not screened inside flexible clusters. As a consequence, the dynamics of an equilibrium-cluster dispersion cannot be described simply as a polydisperse dispersion of rigid Brownian-cluster objects.
Our results contribute to an improved understanding of the dispersed-particles dynamics in SALR systems. Thereby, the results on the influence of HIs on cluster conformational changes and lifetimes are applicable, on a qualitative level, to general Brownian particle systems, where thermodynamic clusters are formed. As a future technological application, the obtained transport properties can be be used as salient input to a previously developed membrane ultrafiltration model,79 in its application to protein solutions under SALR conditions. Furthermore, using MPC simulations, we intend to study, in more detail, the intra-cluster structure and dynamics.
Footnote |
† These authors contributed equally. |
This journal is © The Royal Society of Chemistry 2018 |