Self-subdiffusion in solutions of star-shaped crowders: non-monotonic effects of inter-particle interactions

In figure 2, we present the behaviour of the translational and rotational MSDs for varying interparticle attraction strength ${\epsilon }_{A}$ at crowder packing fraction $\phi =0.12.$ The initial crowder diffusion is ballistic, stemming from the simulation of inertial particles, see also [18, 72, 79, 82, 86] for other simulations and experimental results. At intermediate time scales of ${\rm{\Delta }}\sim 0.1\ldots 10$ we observe a nonmonotonic behaviour of the time averaged MSD that we ascribe to the events of the first collision of a given crowder molecule with another crowder and the oscillations of the centre monomers connected by its springs. In the long lag time limit the translational and rotational MSDs grow linearly with Δ reflecting the Brownian behaviour of the crowder particles.

Zoom In Zoom Out Reset image size

The region of non-Brownian diffusion is centred at ${\rm{\Delta }}\sim 1$ or at about 1 ns in the physical time. Note that for larger tracers, as often used in single particle tracking experiments, the elementary time scale will increase correspondingly and the region of transient anomalous diffusion will shift to longer times. The time span of anomalous diffusion can be extended—as compared to the current self-diffusion in the sea of mobile crowders—if the obstacles/crowders are partly or fully immobilised, see the results of [61, 72].

In this limit the diffusion is ergodic, as we demonstrate in figure 3(A),B. This statement is not necessarily trivial: in many weakly non-ergodic systems the time averaged MSD turns out to be a linear function of the lag time Δ while the ensemble averaged MSD scales as a power-law or logarithmically in time t. This phenomenon was observed in a number of experiments [14, 54–56] and explained in terms of various stochastic processes [7, 29, 30, 111–119]. Figure 3(A) demonstrates both that to very good approximation ergodicity in the sense of the equality

$$\begin{eqnarray*}&&\langle {x}^{2}({\rm{\Delta }})\rangle =\bar{{\delta }^{2}({\rm{\Delta }})}\end{eqnarray*}$$

is fulfilled and that a very small amplitude scatter around the mean $\left\langle \bar{{\delta }^{2}({\rm{\Delta }})}\right\rangle$ exists and thus the individual time averages are reproducible quantities. We furthermore detail the dependence of the particle diffusivity in this Brownian limit versus the attraction strength and the filling fraction in figures 4 and 7, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Let us be more specific. Figure 3(A) illustrates the time averaged MSD for different lengths of the time series as well as the superimposed ensemble averaged MSD shown as the bold black line. As can be seen from the figure, the amplitude scatter of single traces $\bar{{\delta }^{2}}$ around their mean remains small along the entire trajectory, except when ${\rm{\Delta }}\sim T,$ as expected. This growing spread as ${\rm{\Delta }}\sim T$ is a standard feature of even canonical Brownian motion appearing due to progressively poorer statistics when taking the time average [7].

More importantly we observe that the amplitude spread of the time averaged MSD at a fixed lag time Δ decreases as the length T of the time traces increases. This property is ubiquitous for ergodic diffusion processes [7]. We note that the magnitude of the amplitude scatter that we observe for $\bar{{\delta }^{2}}$ for moderately adhering star-shaped crowders are similar to that of a tracer in a network of sticky spherical obstacles, compare figure 3(A) above and figure 7 in [72]. Computing the magnitude of the mean time averaged MSD for varying trace length T we observe that its magnitude stays nearly unchanged with T (figure 3(A)). In the short lag time regime the ballistic scaling is visible.

These observations as well as the explicit evaluation of the ergodicity breaking parameter EB in equation (10) indicate that even at moderate strengths of attraction the diffusion of our crowders stays ergodic. Namely, in figure 3(B) we demonstrate that—using the individual time averaged MSD traces from simulations—the ergodicity breaking parameter (9) follows at intermediate-to-long times the theoretical asymptote for the Brownian motion, equation (10). The EB grows linearly with the lag time Δ and decreases with the length of the trajectory as EB $(T)\sim 1/T,$ see figure 3(B). Some discrepancies from the theoretical prediction (10) we observe at short times stem likely from the ballistic and anomalous diffusion regimes of the time averaged MSD traces at small Δ values. Our process is thus fundamentally different from other anomalously diffusive systems such as those described by continuous time random walks [7] or heterogeneous diffusion processes [29]. For the latter a pronounced scatter of the time averaged MSD trajectories around their mean and a clear dependence of the amplitude of $\left\langle \bar{{\delta }^{2}({\rm{\Delta }})}\right\rangle$ on T at fixed Δ exist, that is, the system ages [7, 120, 121].

The particle diffusivities are defined in our simulations as

$$\begin{eqnarray}&&{D}_{x}=\displaystyle \frac{\left\langle \bar{{\delta }_{x}^{2}({\rm{\Delta }})}\right\rangle }{2{\rm{\Delta }}}\end{eqnarray} \tag{ 11 }$$

and

$$\begin{eqnarray}&&{D}_{r}=\displaystyle \frac{\left\langle \bar{{\delta }_{r}^{2}({\rm{\Delta }})}\right\rangle }{2{\rm{\Delta }}}\end{eqnarray} \tag{ 12 }$$

obtained from the linear behaviour of the mean time averaged MSD in the long time limit ${\rm{\Delta }}\gg 1.$ Figure 4 shows the values of D and D_r extracted from a linear fit of the translational and rotational time averaged MSDs in the range ${\rm{\Delta }}={10}^{3}\ldots {10}^{4}.$ We find that while the rotational diffusivity D_r decreases monotonically, the translational diffusivity exhibits a shallow yet significant maximum at ${\epsilon }_{A}^{*}\approx 1{k}_{B}T.$ This systematic trend persists for the variation of the crowder fraction in a quite broad range (figure 4). This implies that the self-diffusion of our star-like crowders can be facilitated by a weak interparticle attraction. This is one of the main conclusions of this study. One can rationalise this trend in the self-diffusion in terms of the concept of the effective crowder size that decreases for moderate attraction strengths ${\epsilon }_{A}\approx 1{k}_{B}T.$

Note that, as we work in the NVT-ensemble, the number of particles in the simulation box stays constant, whereas their effective size decreases with the strength of attraction ${\epsilon }_{A}$ between the monomers. Upon increase of ${\epsilon }_{A}$ the bond length of the stars decreases, namely about 0.8% for ${\epsilon }_{A}=1{k}_{B}T$ and 1.78% for ${\epsilon }_{A}=2{k}_{B}T.$ With this change of the bond length the effective crowding fraction gets smaller too for higher attraction strengths. If we account for this effect and keep the crowding fraction constant via adjusting the number of particles in the simulation box upon variations of ${\epsilon }_{A}$, the nonmonotonic behaviour of $D({\epsilon }_{A})$ indeed becomes slightly weaker, as demonstrated in figure 4(D). The optimal attraction strength, however, remains of the order of the thermal energy.

Figure 4 illustrates that for progressively stronger star–star attraction their mutual diffusivity decreases eventually to zero due to aggregate formation, see also figure 5 and its discussion below. Also note that the reduction of the rotational diffusion coefficient starts at smaller crowding fractions, see the blue squares in figure 4, as compared to the translational motion that might also feature a nonmonotonic behaviour. The relative reduction of the rotational diffusivity by crowding is also comparatively stronger, as physically expected, because of geometric frustration and overlapping with radially quite inflexible/non-responsive crowders.

We also detect a progressive aggregation of crowders at relatively large crowder–crowder attraction strengths ${\epsilon }_{A},$ as demonstrated in figure 1(C) and in the video files in the supplementary material. This is a well-known phenomenon, for instance, in the glass transitions of dense suspensions of sticky hard spheres [122]. Accordingly, the average diffusivity D of crowders as plotted in figure 4 decreases due to the averaging over an ensemble of particles that perform individual random motions. This average takes into account both particles forming transient aggregates as well as free particles. Roughly speaking the average diffusivity drops inversely proportionally to the number of particles in the cluster. The fraction of particles clustering in these aggregates increases with the mutual attraction strength. The average diffusivity therefore progressively decreases with ${\epsilon }_{A}$ due to a larger fraction of particles in transient aggregates.

At large attraction strength the majority of particles belong to big clusters (results not shown) which diffuse in a Brownian fashion as a whole, with the diffusivity correspondingly reduced by the cluster size. We observe that larger attraction strengths effectively reduce the temperature in the system of crowders, thus favouring aggregation. This temperature argument, however, gets reverted for the self-diffusion of weakly attractive stars, when at a fixed number of particles the crowders becomes slightly smaller and thus diffuse faster for stronger inter-particle attraction.

At a fixed cohesiveness ${\epsilon }_{A}$ of our stars, vicinal crowders create a rough energy landscape for the self-diffusion and the hopping of a given crowder particle. As the MMC fraction ϕ increases, the binding events give rise to a prolonged particle aggregation and reduced self-diffusivity. Above a critical MMC fraction, the barrier height exceeds the thermal energy, thus increasing the lifetime of crowder aggregates significantly. For stronger star–star attraction, the formation of essentially permanent aggregates sets in for less crowded systems, leading to an inhomogeneous, phase-separated spatial distribution, see figure 1(C).

A nonmonotonicity of the translational diffusivity D at similar strengths of the particle–crowder attraction was found in [64] for the tracer diffusion in dense suspensions of spherical Brownian particles. While we here detect that the attraction strength yielding the highest value of D is a function of the crowding fraction ϕ of the stars for the spherical particles, the stickiness facilitating the particle diffusivity was almost ϕ-independent in [64]. The nonmonotonic $D({\epsilon }_{A})$ dependence was interpreted in [64] in terms of the roughness of the free-energy landscape for the tracer diffusion using the concept of the chemical potential. Interestingly, the tracer diffusivity was also nonmonotonic in ϕ in a static regular array of sticky obstacles, as quantified in [72].

We checked the universality of the observed dependencies for $D({\epsilon }_{A})$ and $D(\phi )$ also for disk-like particles. Namely, we simulated just a single monomer of our star-shaped crowders with the given adhesive properties. The diffusivity was indeed found to reveal a maximum at ${\epsilon }_{A}^{\mathrm{opt}}\sim 0.5\ldots 1{k}_{B}T$ (not shown), indicating some universality of this a priori, counterintuitive, faster diffusion for a weak interparticle attraction [64]. Note also that for a polymer chain diffusing in an array of sticky obstacles, a weak chain-obstacle attraction can also substantially enhance the polymer diffusivity [64, 123].

To rationalise the observed behaviour of $D({\epsilon }_{A})$, we calculate in figure 5 the potential of the mean force between two crowders as

$$\begin{eqnarray}&&F(r)=-{k}_{B}T\mathrm{log}[\rho (r)].\end{eqnarray} \tag{ 13 }$$

In this reconstructed, approximate free energy (valid for dilute systems only) the quantity $\rho (r)$ is the average radial distribution function of the centre crowders monomer in the steady-state long time limit. As the mutual attraction strength ${\epsilon }_{A}$ increases, we observe that the potential well at the separation $r\approx 2\sigma$ becomes deeper, see the first well in figure 5. Concurrently, the distance at which F(r) sharply increases becomes shorter for larger ${\epsilon }_{A}.$ For a stronger star–star attraction the crowders feature a more organised appearance, resulting in measurable oscillations of $\rho (r)$ and $F(r),$ as evidenced in figure 5.

These trends indicate that the effective crowder radius gets smaller with increasing ${\epsilon }_{A},$ and, at an optimal value ${\epsilon }_{A}^{\mathrm{opt}}$, the crowders approach one another more closely, yet without sticking. This in turn might result in a faster average diffusivity D at ${\epsilon }_{A}\approx {\epsilon }_{A}^{\mathrm{opt}},$ as we indeed observe. An effective reduction of the crowder size at optimal attraction strength is one important cause—albeit possibly not the only one—for this facilitated diffusion. In the current system, the equilibrium distance of the outer monomers from the central monomer is reduced by about 2% for the inter-monomer attraction strength of $2{k}_{B}T.$ Even higher values of ${\epsilon }_{A}$ give rise to the formation of large clusters of crowders, see the supplementary material. As shown in figure 4, the corresponding diffusivity of an average particle drops dramatically with ${\epsilon }_{A}.$

In figure 6(A), we show the translational and rotational MSD for varying packing fraction of crowders ϕ. As expected—from a linear fit to the long time time averaged MSD—the diffusivity is a monotonically decreasing function of ϕ, as evidenced by figure 7. For more severely crowded systems, the tracer diffusion gets more obstructed and the magnitude of the corresponding mean time averaged MSD $\left\langle \bar{{\delta }^{2}}\right\rangle$ decreases. To elucidate these effects further, we evaluate from the time averaged MSD traces of figure 6(A) for the translational motion the local diffusion exponent [7, 124]

$$\begin{eqnarray}&&\beta ({\rm{\Delta }})=\displaystyle \frac{d\mathrm{log}\left(\left\langle \bar{{\delta }^{2}({\rm{\Delta }})}\right\rangle \right)}{d\mathrm{log}({\rm{\Delta }})}.\end{eqnarray} \tag{ 14 }$$

For the rotational motion the exponent ${\beta }_{r}(t)$ is defined analogously.

We observe a ballistic regime with $\beta \approx 2$ in the particle diffusion at short times, figure 6(B). This ballistic regime is followed by a decrease and further increase of the scaling exponent at ${\rm{\Delta }}\sim 1.$ These nonmonotonic trends are also clearly visible from the behaviour of the time averaged MSD traces themselves as a function of the lag time Δ, see figure 6(A). We find that the variations of the scaling exponent for translational and rotational motions of the star-like crowders appear correlated, indicating a coupling of these diffusion modes [20]. In the plots for the scaling exponent $\beta (t)$ in figure 6(B) the significant spike-like signal at ${\rm{\Delta }}\sim 1$ is interpreted as an effect of the first collision of particles and the resulting onset of an effective confinement. We note that even in effective one-particle theories pronounced oscillations occur at the crossover point between the initial ballistic and the overdamped regime [125, 126].

With increasing crowder fraction ϕ we also observe a more pronounced range of anomalous diffusion for lag times of the order of ${\rm{\Delta }}\sim 1\ldots 100.$ This range appears strongly correlated for rotational and translational particle motion, as shown in figure 6(B). For rotational diffusion, the scaling exponent drops practically to zero for times longer than those of the initial ballistic growth, and the corresponding mean time averaged MSD trace $\left\langle \bar{{\delta }_{r}^{2}}\right\rangle$ exhibits a short plateau (figure 6(B)). Such a transient subdiffusion was observed for a number of systems [5–7], see also the Introduction. Especially in dense colloidal systems close to the glass transition at $\phi ={\phi }^{\star }$ this subdiffusion is accompanied by an exponential growth of the solution viscosity $\eta =\eta (\phi )$, which is divergent at $\phi \to {\phi }^{\star }$ [76]. The colloidal glasses also exhibit progressive particle localisation effects as discussed in [76, 77]. In the long lag time limit the exponent becomes Brownian $\beta \approx 1$ for the crowded systems far from the critical occupation ${\phi }^{\star }.$

Remarkably, the relative variation of the translational and rotational diffusivities with the crowding fraction of stars is quite similar. For comparison, we plot in figure 7 the theoretical prediction for dense suspensions of hard spheres [127, 128]

$$\begin{eqnarray}&&\displaystyle \frac{D(\phi )}{D(0)}={\left[1-\displaystyle \frac{\phi }{{\phi }^{\star }}\right]}^{2},\end{eqnarray} \tag{ 15 }$$

with the critical packing fraction for our system of ${\phi }^{\star }({\epsilon }_{A}=1{k}_{B}T)\sim 0.52$ that provides the best fit to the data. Above this value ${\phi }^{*}$ both translational and rotational diffusivities of the crowders essentially vanish. If the explicit dependence of the particle size on the attraction strength is taken into account, as in Panel D of figure 4, the critical fraction ${\phi }^{\star }$ also becomes a function of ${\epsilon }_{A}.$ At this critical crowding fraction the interparticle attraction becomes so strong that the self-diffusion is almost completely localised and the motion of particles corresponds more to a very restricted wiggling and jiggling. The translational MSD of individual crowders—after subtracting the diffusive motion of the entire cluster as a whole—saturates to a plateau with the scaling exponent $\beta \to 0$ (results not shown). As expected, when the star–star interactions become stronger, some aggregate formation sets in for less crowded systems, and thus the critical value ${\phi }^{\star }$ is diminished (not shown). Albeit the theory in [127, 128] is developed for three-dimensional suspensions in the presence of hydrodynamic interactions, it agrees remarkably well with our results, as shown in figure 7. The reader is also referred to [129] for experimental data of the crowding dependent diffusivity of colloidal particles and alternative theoretical predictions for the diffusivity $D(\phi ).$

We note that [130] suggest exponential rather than power-law forms for the particle diffusivity in crowded solutions. Also, different scalings with the crowding fraction were predicted for semi-dilute solutions of multi-arm polymers and hyper-stars [92], where $D(\phi )\sim {\phi }^{-1/2}.$ The polymeric nature of star arms plays, however, a dominant role in this scaling relation. Also, the decrease of the particle diffusivity in a system of polydisperse hard disks exhibiting a glass transition was shown to follow the relation [82, 85]

$$\begin{eqnarray}&&D(\phi )\sim {(\phi -{\phi }^{\star })}^{2.4},\end{eqnarray} \tag{ 16 }$$

in agreement with the results of the mode coupling theory, see, e.g., [83, 131]. In figure 7, we present both theoretical asymptotes with the scaling exponents of 2 and 2.4, equations (15) and (16). The latter indeed fits better the decrease of the translational diffusivity $D(\phi )$ of star-like crowders, but not the rotational diffusivity ${D}_{r}(\phi ).$ From the limited simulation data available we cannot determine the value of the scaling exponent more precisely.