Geometric percolation of hard-sphere dispersions in shear flow

Abstract

We combine a heuristic theory of geometric percolation and the Smoluchowski theory of colloid dynamics to predict the impact of shear flow on the percolation threshold of hard spherical colloidal particles, and verify our findings by means of molecular dynamics simulations. It appears that the impact of shear flow is subtle and highly non-trivial, even in the absence of hydrodynamic interactions between the particles. The presence of shear flow can both increase and decrease the percolation threshold, depending on the criterion used for determining whether or not two particles are connected and on the Péclet number. Our approach opens up a route to quantitatively predict the percolation threshold in nanocomposite materials that, as a rule, are produced under non-equilibrium conditions, making comparison with equilibrium percolation theory tenuous. Our theory can be adapted straightforwardly for application in other types of flow field, and particles of different shape or interacting via other than hard-core potentials.

---

Introduction

The electrical conductivity of polymeric materials can be varied over ten orders of magnitude by the incorporation of a relatively small fraction of conductive nano-fillers, such as carbon black, graphene and metallic particles.^1–9 The intense interest in this topic, evidenced by a huge surge in the number of studies dealing with polymeric nano-composite materials in the last few decades, is perhaps not entirely surprising given their potential technological applications in, say, opto-electronics, photo-voltaics and electromagnetic interference shielding.^10–14 It turns out that the degree of homogeneity of the nano-particle dispersion in the host material is of crucial importance to the level of conduction of the composite achieved.^15–22 Hence, great care is taken in the manufacturing process to disperse the nano-fillers evenly when the host material is still in the fluid stages of the production. For this purpose, techniques are applied that include sonication, manual mixing and shear mixing, followed by casting and curing of the composite.^23–27 These methods contribute to the homogeneous dispersion of the particles and prevent their aggregation. Because the curing of the fluid is typically (but not always²⁸) done relatively quickly after the mixing so as to avoid re-aggregation, the out-of-equilibrium structure of the nano-fillers should be expected to be essentially frozen-in in the final, solid composite.

Consequently, if one attempts to predict the percolation threshold of composite materials using the standard tools of liquid state theory, as is usually done in connectedness percolation theory,²⁹ then it stands to reason that these predictions must be flawed. Indeed, connectedness percolation theory assumes the particle distribution to obey (equilibrium) Boltzmann statistics.^30,31 To remedy this for those conditions where the particle distribution does not obey equilibrium statistics, a quantitative, out-of-equilibrium continuum percolation theory is sorely needed. Unfortunately, no such theory is, as far as we are aware, currently available. One could envisage setting up a non-equilibrium version of connectedness percolation theory, based on the Smoluchowski equation for the steady-state pair correlation function under flow.³² However, in the dynamical theory there is no obvious way of separating the non-equilibrium equivalent of the so-called connectedness and blocking functions, as is possible in thermodynamic equilibrium.³⁰

An alternative that does not have this problem, is to take a more heuristic approach such as that we put forward in this work. Here, we use a simple geometric criterion for the percolation threshold that depends solely on the pair correlation function and a connectivity criterion. The method quantitatively describes results for different particle shapes from computer simulations under conditions of thermal equilibrium.³³ We combine this with the known steady-state solution of the Smoluchowski equation for the pair correlation function in shear flow, that we solve in the limit of low volume fractions but apply also to intermediate and high concentrations.³⁴ Together, these two ingredients provide us with a simple and tractable way to predict the percolation threshold in dispersions that are out-of-equilibrium, which we apply to the case of colloidal hard spheres in simple shear flow. We choose to study hard spheres over particles with more realistic colloidal interactions for reasons of simplicity. We note that the method is quite generic and allows for a straightforward extension to other types of flow field and other types of inter-particle interaction and particle shape. In future work, we intend to additionally study the effects of colloidal attractions. Simple shear flow, parameterised by the velocity field v ∝ (y, 0, 0), is among the most studied flow fields due to it being the simplest model for the more complex shear flows that are ubiquitous in industrial and experimental processes. Moreover, it is well known that the properties of colloidal materials change drastically when subjected to such flows.^35,36 Specifically, both shear-induced cluster formation and breaking-up has been reported in the literature,^37–40 making its influence on the percolation threshold highly unclear. With this work we aim to shed light on what impact a shear flow field may have on connected clusters, and in particular under what conditions these break up or grow to give rise to percolating particle networks.

As we shall see, while our theory loses its quantitative nature for high shear rates when compared with our Langevin dynamics simulations, it does overall describe the full phenomenology of how the percolation threshold depends on the connectivity range, that is, the geometric criterion defining whether two neighbouring particles are connected or not. This agreement is reached without including any free (fitting) parameters. In agreement with our simulations, we find that for sufficiently large connectivity range, fluid flow increases the percolation threshold and more so the larger the Péclet number, whilst for small connectivity range the opposite happens. This is caused by the delicate balance between the compression and extension that the flow field exerts on the pair structure. The impact of fluid flow on the percolation threshold remains modest, however, staying within 15% of the static case for Péclet numbers up to about ten. Our predictions are summarised in Fig. 1a and the results of our simulations in Fig. 1b.

Fig. 1 Theoretical (a) and simulation (b) results of the dimensionless percolation threshold of a dispersion of hard spherical particles subject to a simple shear flow as function of the ratio of the hard core diameter D and connectivity length λ. The strength of the shear flow is quantified with the Péclet number Pe. In (b), we indicate with the cross the literature value of ρπλ_c³/6 = 0.341889 valid in the case of non-interacting particles⁴¹ and we add lines between the obtained data points as guides to the eye. The insets show the relative change of the percolation threshold compared to the no-flow case. In all figures, the maximal value of the hard-core diameter is given by D/λ_c = 0.96.

In the following, we first summarise the ingredients of our theory, consisting of the geometric percolation theory of Alon and collaborators,³³ and the analytical prediction of the correction of the pair structure by shear flow of Bławzdziewicz et al.³⁴ The latter we compare with our Langevin dynamics simulations. We subsequently integrate both ingredients and obtain a prediction of the percolation threshold, and end the paper with conclusions and an outlook. Details of our calculations and simulations are given at the end of this paper.

A heuristic approach to percolation

The theoretical framework of Alon et al. for predicting the percolation threshold is ideally suited for our problem of percolation far out of equilibrium, as it is geometric in nature and does not require thermodynamic equilibrium to hold.³³ We summarise it below for the case of spherical interacting particles but note that our treatment can be extended to non-spherical particles as well, see ref. 42 for an example in the case of non-interacting particles.

The theory presumes N homogeneously dispersed particles to be present in a volume V. We quantify the structure of this dispersion with the pair correlation function g(r): if a particle is placed at the origin, then the probability of finding another particle in volume element d³r at position r is equal to ρg(r)d³r, where ρ = N/V is the number density. We further assume that pairs of particles with a centre-to-centre distance |r| smaller than the connectivity length λ are directly connected to each other, implying that charge carriers can be transported efficiently between them. Physically, this connectivity range λ can e.g. be interpreted as a tunnelling length.^43,44

The question arises how to find the smallest connectivity length λ_c for which a connected cluster of particles exists that spans the entire material for a given density ρ. This is equivalent to asking what the critical density ρ_c is at which a percolating cluster appears for a given connectivity length λ.^45,46 The advantage of the former formulation is that finding the structure for a given density and finding the percolation threshold for a given structure are now completely decoupled problems, and no longer have to be solved self-consistently.

The argument of the theory is as follows. Clearly, a macroscopically connected cluster must comprise many so-called backbone particles with two or more connections, since particles with one or no connections at all cannot propagate connectivity and therefore do not contribute to percolation. Heuristically, Alon and coworkers argue that percolating networks only exist if the average distance L between such backbone particles is smaller than twice the average distance l between directly connected particles (for which |r| < λ). The percolation threshold may thus be found by requiring that L = 2l.³³ Estimates for both lengths l and L can be obtained relatively straightforwardly from the pair correlation function, g(r).

To start with the first and recalling that the probability of finding a particle in volume d³r at position r is ρg(r)d³r, we can use an appropriately normalised (statistical) moment of g(r), evaluated within the connectivity region, to estimate the mean distance l between connected particles: . We use the second moment rather than the first because it yields more accurate predictions for the percolation threshold.³³ In the expression for l², V_λ is the connectivity region, which in our case is the region within the sphere of radius λ.

As to the distance L between backbone particles with two or more direct connections, it is convenient to assume that the probability P_k that a particle has k direct connections is Poissonian, i.e., P_k = B^k exp(−B)/k!, in which is the average number of direct connections. This assumption is exact for non-interacting particles,^47,48 and from our molecular dynamics simulations we find that it remains a reasonable approximation for dispersions of hard spheres at low to intermediate densities (results shown in Fig. 3(g)). The mean number density of particles with at least two neighbours, ρ₂, can now straightforwardly be found from ρ₂ = ρ(1 − P₀ − P₁) = ρ(1 − (1 + B)exp(−B)). Assuming that the volume available to each particle is spherical, we have L = 2(4πρ₂/3)^−1/3.

In conclusion, the criterion for percolation, L = 2l, we find to be a function only of ρ, g(r), and λ, and for any given density ρ and corresponding pair correlation function g(r), we can use known numerical routines⁴⁹ to solve for the percolation threshold λ_c. Of course, we need to know g(r) too. It can be evaluated directly either from computer simulations, or, say, from the (numerical) integration of the Ornstein–Zernike equations with a suitable closure. For the case of hard particles, the Percus–Yevick closure is known to be highly accurate.⁵⁰

In Fig. 1a and b we also show our theoretical and simulation results for zero Péclet number, so in the absence of a shear field, the former within Percus–Yevick theory. For this, we numerically integrated the Ornstein–Zernike equation within the Percus–Yevick closure.⁵¹ The difference between theory and simulation is less than 6% for all values of the connectivity length λ, which incidentally outperforms approximations derived from rigorous liquid state theory significantly.⁵² Notice that the limit D/λ → 0 corresponds to the penetrable-sphere limit. For increasing values of the hard-core diameter D at fixed value of λ, D/λ increases and the percolation threshold initially decreases before it increases again. This is caused by the competition of an increasing contact value of the pair correlation function and a shrinking connectivity region. As we shall see, the presence of a flow field alters that competition. To study this, we first need to evaluate the impact of flow on the pair structure. A brief discussion of this we provide in the next section.

Pair correlations under shear flow

In order to find the percolation threshold in particle suspensions subject to a flow field, we must take into account the impact of the flow field on the pair correlation function, g(r). To do this, we approximate the pair correlation function using a steady-state two-particle Smoluchowski equation.

In the absence of hydrodynamic interactions and many-body correlations, this two-body Smoluchowski equation reads³²

∇·(Γrg − 2D₀(g∇V/k_BT + ∇g)) = 0, (1)

for an arbitrary velocity gradient tensor Γ and spherically symmetric interaction potential V/k_BT. Here, we introduced the thermal energy k_BT to non-dimensionalise the pair potential.

We restrict ourselves to the case of simple shear flow, with velocity field v = Γr = [small gamma, Greek, dot above]yx̂, where [small gamma, Greek, dot above] is the shear rate, x̂ the unit vector in the flow direction and y the coordinate along gradient direction of the flow field. To calculate the pair structure in such simply sheared liquids, we view it as a (not necessarily small) correction δg(r, Pe) = g(r, Pe) − g₀(r) of the equilibrium pair correlation function g₀(r). Here, Pe = [small gamma, Greek, dot above]D²/4D₀ is the Péclet number, with D the particle diameter and D₀ the self-diffusion constant.† For spherical particles, the equilibrium pair correlation function g₀(r) only depends on the radial distance r = |r| between the particles, for which reason it is usually referred to as the radial distribution function. Under shear, this is no longer true as the flow field breaks the spherical symmetry. We neglect the effect of hydrodynamic interactions to keep the theory concise. At the end of the next section, we shall discuss the validity of this approximation.

To apply (1) to the case of a hard-sphere liquid, Bławzdziewicz and Szamel³⁴ have imposed no-flux boundary conditions at r = D and set V = 0 for r > D. This results in the boundary value problem for δg(r, Pe) given by

Since many-body correlations have been neglected here, the prediction for δg(r) does not depend on the particle density. This approximation is equivalent to setting the reference pair correlation function g₀(r) = exp(−V(r)/k_BT). We therefore expect the predicted flow-induced correction to be accurate only for sufficiently low volume fractions. In the next section, we combine this correction, strictly valid in the dilute limit, with a more realistic equilibrium structure at high densities.

We choose to use the theory of Bławzdziewicz and Szamel,³⁴ rather than more sophisticated approaches,^53–61 because even with this relatively simple model we obtain quite accurate results for the percolation threshold. An additional benefit is that the partial differential eqn (2), together with hard-sphere no-flux boundary condition (3), admit an analytical treatment that we summarise in Methods section.³⁴

We perform a direct test of the theory by comparing the prediction for δg(r, Pe) with results from our Langevin dynamics simulations in Fig. 2. The figure shows the shear-induced correction δg obtained from the theory and that from our simulations at two densities in the xy-plane, for two different Péclet numbers. We see that for low volume fractions φ = πρD³/6 ≪ 1, the simulation results quantitatively match the theoretical model, at least up to Pe = 10. (For Pe = 1 we take for our simulations φ = 0.1 rather than the φ = 0.01 that we used for Pe = 10 to improve the statistics.) For high volume fractions, however, we find that highly complex and long-ranged pair correlations are induced by the shear field, which differ not only quantitatively but in fact also qualitatively from those that our theoretical model predicts. This is not all that surprising, given the approximations of the model.

Fig. 2 Comparison between theory and simulations of the shear-induced correction of the pair correlation function δg(r, Pe) at low and high volume fractions of hard, spherical particles in the flow-gradient plane. The figures in the first column (a and d) correspond to the theoretical model, whereas the second (b and e) and third (c and f) column correspond to simulation results at low and high volume fractions φ. The first row (a–c) represents dispersions at Péclet number Pe = 1, and for the second we set Pe = 10. Each row shares a colour scheme that indicates the values of δg(r, Pe). For clarity, we added a thin black line at r = D indicating the theoretical excluded volume of hard, spherical particles.

At high densities, long-ranged correlations are in fact also present in the equilibrium pair correlation function, shown in Fig. 3(e), to which the shear flow couples. At low Péclet numbers, the peaks of the equilibrium pair correlation function are amplified in the compression quadrants and suppressed extensional quadrants. The opposite happens with the troughs of the equilibrium radial distribution function. In the compression quadrants, the flow increases the magnitude of structural correlations because the local density increases, whereas the opposite happens in the extensional quadrants. If the strength of the shear flow increases, this persists in the compression quadrant but leads to complex structural patterns in the extensional quadrant.

Fig. 3 Quantitative comparison of the predicted structure with molecular dynamics simulation data. (a–d) Shear-induced correction of the pair correlation function δg(r, Pe) of spherical particles along different curves through three dimensional space at volume fraction φ = 0.01 (blue circles) and φ = 0.4 (red diamonds) for Péclet number Pe = 10. Respectively, they correspond to the correction at constant radial distance r = 1.2D in the xy-plane, and along the curves ., , and r = rẑ. Indicated are also the theoretical predictions (drawn line). (e) Equilibrium radial distribution function at volume fractions φ = 0.01 (blue) and φ = 0.4 (red). (f) Flow-induced correction of the pair correlation function averaged over all solid angles for different Péclet numbers as predicted by the theory. Colour coding of the curves is given in the legend. (g) Probability P_k that a particle has k direct connections at the percolation threshold λ_c for Pe = 10. The simulation results are compared to Poisson distributions where the average number of neighbours B is determined from the simulation results, that is, B = 2.7, 1.8, 1.1 for ρD³ = 0.0005, 0.25, 0.88, respectively. (h) Percolation probability P_p(λ) in the absence of shear flow for D/λ = 0, determined for different particle number N indicated in the legend. The intersection of the curves is shown in the inset.

A more quantitative comparison we provide in Fig. 3(a–d), where we plot the same results for δg(r, Pe) presented in Fig. 2(d–f) for Pe = 10 along specific curves. In Fig. 3(a) we vary the azimuthal angle ϕ = arctan y/x at fixed distance in the xy-plane and in (b–d) the radial distance r for fixed angles θ and ϕ.

The figures confirm that the theoretical model reproduces very well the shear-induced correction of the pair correlation function δg at low volume fraction. At high volume fractions of particles, we again see that the theoretical model qualitatively fails to predict the features of δg(r, Pe). In fact, we notice the emergence of an additional peak in the correction of the pair correlation function at constant radius, see Fig. 3(a). In the constant-angle plots, (b–d), we obtain additional peaks and a large increase in the structural correlations at high particle density. In Fig. 3(d), we show the δg along the line perpendicular to the shear-plane. This correction is clearly nonzero, which is consistent with earlier findings from theory and simulations.^34,57,62,63 Notice also that the range of these correlations increases drastically at a volume fraction of φ = 0.4, which is close to the freezing transition in the absence of a flow field.

What these results point at, is that in the xy-plane along the axis ϕ = π/4 (in the extensional quadrants) the contact value of g(r) goes down, whilst along the axis ϕ = 3π/4 (in the compressional quadrants) it increases very strongly. This, of course, has consequences for the structure, size and shape of (geometric) clusters of particles, as the contact value of g(r) informs on the likelihood of the presence of other particles near the surface of a test particle. In turn, this must have an impact on the percolation transition, which occurs when macroscopic clusters form. We investigate this next in more detail by making use of the prediction for δg(r, Pe) and the heuristic percolation criterion to investigate the percolation threshold for different Péclet numbers.

The percolation threshold

In the previous section we calculated the correction δg(r, Pe) = g(r, Pe) − g₀(r) of the equilibrium radial distribution function g₀(r), if we presume the latter to obey a simple Boltzmann weight. For hard particles this becomes a step function, and is accurate only for very low volume fractions φ ≪ 1. Since we are not necessarily only interested in results at low volume fractions, we use the g₀(r) obtained from the numerical integration of the Ornstein–Zernike equation together with the Percus–Yevick closure⁵¹ to obtain an ad hoc approximation for the full multi-body g(r, Pe). Clearly, it is only multi-body for the (reference) equilibrium structure, while the impact of flow is only described at the two-body level. Inserting this into the percolation criterion yields the percolation threshold, which we present in Fig. 1.

We show the results on the effect of shear flow on the percolation threshold according to our model in Fig. 1a. In Fig. 1b, we compare them to percolation thresholds obtained from the analysis of snapshots of our molecular dynamics simulations, on which more information can be found in the Methods section. We plot the critical volume fraction ρπλ_c/6 of connectivity regions as a function of the relative size of the hard-core particle diameters D/λ_c. If the hard-core diameter D goes to zero and particle interactions become negligible, we find that the effect of the shear flow on the percolation threshold does so too. In this case, our simulation results agree very accurately with the literature value of ρπλ_c/6 ≈ 0.341889.⁴¹ The fact that a flow field does not impact this number is intuitive, since the flow field cannot induce any structural changes if the particles do not interact, i.e., if D = 0. In the intermediate and high D/λ_c regime, both theory and simulations agree that an applied shear flow can both increase and decrease the percolation threshold of a fluid dispersion of spherical particles, depending on their hard-core diameter and strength of the flow field. As already alluded to, the effect of flow on the percolation threshold is modest with deviations of at most 15%. See also the insets of Fig. 1, that show relative changes compared to the no-flow case.

For all Péclet numbers studied, we find that there is a value of D/λ_c, below which the percolation threshold increases and above which it decreases if compared with the Pe = 0 case. Both the decrease and increase are explained by the theory. Since the orientationally averaged correction of the pair correlation function is positive for short distances, see Fig. 3(f), the flow field induces an increase in the average number of neighbours for small λ. This naturally translates to a lower percolation threshold for large D/λ. Conversely, for large separation distances r, the orientationally averaged correction becomes negative, meaning that the flow field decreases the number of neighbours for large λ, thereby increasing the percolation threshold for small D/λ.

We believe that the shear-induced decrease of the percolation threshold for large D/λ_c is closely related to the emergence of so-called shear-induced contact clusters,⁴⁰ which have also been found in earlier simulations.^39,64 Even at volume fractions where such clusters are finite, they might play a significant role in aiding long-range connectivity and therefore in decreasing the percolation threshold with respect to the equilibrium situation. In fact, from our simulations we find that the average cluster becomes progressively elongated as the Péclet number increases, in line with what has been found in simulations of such contact clusters.³⁹

The results of Fig. 1 confirm that there is good qualitative agreement between our theoretical model and simulation results of the percolation threshold. Our model correctly predicts the shear flow to induce an increase and subsequent decrease of the percolation threshold with decreasing connectivity range, and gives an accurate estimate of the location of the crossovers between these two regimes. It predicts the shift in the percolation threshold with almost quantitative accuracy as long as the material is sufficiently dilute, that is, as long as, say, D/λ_c < 0.8, which roughly corresponds to hard-core volume fractions φ < 0.2. However, as seen more clearly in the insets of Fig. 1, our theory does not quite capture the impact of shear flow on the percolation threshold for large values of D/λ_c. In that case, the hard-core volume fraction is high and the shear flow significantly affects many-body contributions to the pair structure, as we also show in Fig. 2. The reason why a flow field seems to only have a modest effect on the percolation threshold for all values of D/λ_c, even though it strongly impacts on the pair structure, is probably due to the circumstance that the orientational average of the many-body corrections remains small. This can in fact be deduced from Fig. 2(c and f), which show that corrections in the compression quadrants in part compensate for those those in the extensional quadrants, especially in the low Péclet number regime.

We expect that by taking many-body correlations in the distortion of the pair correlation function by the flow field explicitly into account, the accuracy of our predictions would improve slightly for values of D/λ_c approaching unity. Fortunately, for the largest part of our work, the hard-core volume fractions at the percolation threshold remain rather low, that is, lower than 0.2 as long as the hard-core diameter remains smaller than 80% of the connectivity range. At least for those cases, many-body correlations play only a subdominant role,⁶⁵ and should not be expected to influence the results much, as is evidenced by the agreement between the percolation thresholds obtained from our theory and simulations.

Even if many-body correlations had been included, they would not have remedied the slight disagreement between theory and simulation of the location of the minimum of the percolation threshold as a function of D/λ_c that the theory overestimates. As this is already the case for Pe = 0, this must in part be caused by the heuristic percolation theory we use. It would be interesting to see whether a novel continuum percolation theory that was recently proposed, nearest-neighbour connectedness theory,^66,67 improves the agreement. Although it is also based on geometric considerations, it is not obvious how to apply it to out-of-equilibrium percolation as it does not rely solely on the concept of a pair correlation function.

For reasons of conciseness, neither our theory nor our simulations include hydrodynamic interactions between the colloids, e.g., in the form of short-ranged lubrication forces or long-ranged hydrodynamic many-body interactions. For low volume fractions and small Péclet numbers, we expect that changes in the perturbation of the structure due to hydrodynamic interactions are of quantitative nature only, and consequently this must also hold for their impact on the percolation threshold. In the limiting case of vanishing flow fields, indeed, the material is in equilibrium and hydrodynamic interactions cannot change the structure of the dispersion.⁶⁸ Consequently, they can neither have an effect on the percolation threshold. For strong shear cases, their exclusion might not be justified.

Conclusions and outlook

We have presented a theoretical framework that describes the geometric percolation of hard, spherical particles under shear flow within the free-draining approximation. The predictions of the theory compare favourably with results of our Langevin dynamics computer simulations. We find that the percolation threshold is determined, on the one hand, by the ratio of the connectivity range and the hard-core diameter of the particles and, on the other, by the Péclet number that measures how much diffusion is affected by the flow field. The Péclet number is proportional to the flow rate.

In the absence of a flow field, so at zero Péclet number, and at fixed connectivity range, the percolation threshold initially decreases with increasing hard-core diameter, to subsequently increase again when the hard-core diameter approaches the connectivity range. If we ramp up the strength of the flow field, the percolation threshold increases for small hard-core diameters smaller than some critical value, yet decreases if the diameter is larger than that. This critical value we find to depend on the Péclet number. Our calculations show that this is caused by the balance between compressional and extensional effects that the shear flow exerts on the local particle density field near a reference particle.

According to our simulations and theoretical predictions, the influence of simple shear flow on the percolation threshold of hard colloidal spheres remains modest for Péclet numbers below 10. This suggests that for many practical applications, where a high-precision prediction is not required, the subtle influence of shear flow on the percolation threshold may well be disregarded. In that case, calculations based on equilibrium connectedness percolation theory would suffice. Whether this conclusion extends to sticky colloids or to particles that are not isometric, remains to be seen.

The disagreement of our theoretical predictions with our computer simulations of the percolation threshold is at most 7% for Péclet numbers up to ten. This error depends only very weakly on the shear rate, and because of that must originate mainly from either the heuristic percolation criterion or the equilibrium pair correlation function. It becomes shear rate dependent only if the hard-core diameter approaches the connectivity range. In that case, there are obvious routes to improvement for the theory. A clear path forward involves the inclusion of higher-than-two-body correlations in the bare correlation that serves as input for the calculations. See, for example, ref. 50, 57 and 68–71 for more detailed discussions on how to deal with the inclusion of such higher order correlation functions.

Throughout this work, we have disregarded hydrodynamic interactions because they would render the molecular dynamics simulations very computationally expensive, and complicate the analytic treatment of the theory. As we argue in the main text, we expect that this does not introduce serious errors to the percolation threshold for low Péclet numbers. For stronger flow fields, however, hydrodynamic interactions should probably be taken into account to retain qualitative accuracy of the theory.^57,72 We note that this is highly non-trivial, and requires the inclusion of additional terms in eqn (1).^{32,36,62,69,72–75}

Because of its simplicity and generality, the framework outlined here may prove a useful tool for modelling percolation in more realistic and perhaps interesting systems or setups. Our theory can straightforwardly be extended to more complicated flow fields and interaction potentials by solving the two-body Smoluchowski eqn (1) either analytically or numerically, and inserting the result in the heuristic percolation criterion.

To apply the theory to the case of non-spherical particles, the pair correlation function depends also on the orientations of both particles and additional terms need to be added in eqn (1) to account for rotational diffusion.^68,76 Computer simulations and experiments indicate that such orientational degrees of freedom give rise to fundamentally different clustering behaviour than is the case for dispersions of spherical particles.^77–81 It appears this is due, in large part, to the aligning effect a flow field exerts on highly non-isometric particles.^82,83

To study the percolation of colloidal particles in polymeric hosts, which is the original motivation of this work, there are two obvious routes to improve the model. Firstly, the effects of colloidal attractions would have to be incorporated into the theory. We intend to pursue this in further research. Experiments suggest that these effects have a major influence on the percolation threshold of the nanocomposite.⁸⁴ Secondly, the non-Newtonian nature of the background fluid has to be incorporated too, especially in cases of strong flow fields. This can be achieved pragmatically by making the colloid self-diffusion coefficient shear-rate dependent in a way that faithfully models the shear thinning or thickening behaviour of the studied host material.

Methods

In this section we present technical details that are not essential the the main thesis of this work. We start by presenting a sketch of the derivation of the solution of the Smoluchowski equation, and present details on the Langevin dynamics simulations thereafter.

Analytical solution of the Smoluchowski Equation

We follow Bławzdziewicz and Szamel to solve the two-particle Smoluchowski boundary value problem given by eqn (2) and (3), using the method of induced multipole sources.^34,85 First, we note that the fundamental solutions of the time-dependent version of the Smoluchowski eqn (2)are given by⁸⁶

To be able to satisfy the boundary conditions (3), we place a linear combination of multipole sources at r = 0. The multipole source functions T_αβ(x, y, z, t) can be found by taking repeated spatial derivatives of the fundamental solutions

Since we are focused on finding the steady-state solution of eqn (4) rather than its time-dependent behaviour, we use the limit method⁸⁷ to find steady state multipole source functions T_αβ(x, y, z)

The full solution of the steady state boundary value problem is now

in which C^αβ are coefficients that are used to satisfy the boundary condition (3).

To fix the coefficients C^αβ, we insert (8) into (3) and expand both the left- and the right-hand side in terms of real-valued spherical harmonics Y_lm(x, y, z). Using the fact that , the boundary condition can now be formulated as

or, more compactly,

from which the expansion coefficients j^lm_αβ follow by invoking the orthogonality property

evaluated at r = D. Having found the expansion coefficients j^lm_αβ, we evaluate the coefficients C_αβ by inverting the linear system of equations given by (10). In agreement with Bławzdziewicz and Szamel, we find that the procedure converges if we take into account all coefficients such that α + β ≤ 10, at least for Pe < 10.

The spherical integrals we performed numerically using the Lebedev quadrature of order 47.⁸⁸ The t-integrals we performed using Simpson's rule on a logarithmic grid of 100 points spanning from t = 10⁻⁵ to t = 10⁵.⁸⁹ In order to calculate the radial integrals occurring in the expressions of the two length scales l and L, we use a trapezoidal rule, and subsequently solve the percolation criterion 2l = L using Broyden's rule,⁹⁰ in order to find the percolation threshold.

Molecular dynamics simulations

In order to verify our theory, we perform particle-resolved simulations using the LAMMPS software package,⁹¹ explicitly integrating the Langevin equationfor every particle.⁹² In the Langevin equation, we have introduced the mass m = 1, friction coefficient γ = 10, and the fluctuating force R(t) that has zero mean and a variance dictated by the fluctuation dissipation theoremin which R_i and R_j are components of the vector R. The potential energy is given by the sum over all pair potentials, for which we choose the Weeks–Chandler–Anderson form,⁹³if r_ij < 2^1/6σ and U_ij/k_BT = 0 otherwise. We choose to set σ = 2^−1/6D and ε = 100. This choice ensures that we accurately model hard-sphere behaviour.

We introduce shear flow in our simulations by continuously deforming our simulation box every time step such that the shear strain is equal to [small gamma, Greek, dot above] = 1 in dimensionless units. This requires us to remap the particle velocities when they cross the periodic boundaries in the y-direction. This method is equivalent to introducing Lees–Edwards boundary conditions.⁹⁴ In the friction term on the right hand side of eqn (12), we compute the friction force by multiplying the friction coefficient γ by the velocity of the particle relative to the background flow field. Together, these methods result in a simulation of simple shear flow with a linear average velocity profile. As a consistency check, we verified that the slope of this profile agrees with the one that we impose.

The Péclet number we now vary through the diffusion constant D₀ = k_BT/γ by adjusting the temperature. By changing the density, we effectively change D/λ_c, allowing us to probe the influence of hard-core interactions on the percolation threshold.

For densities ρD³ < 0.01 we set our time step equal to Δt = 10⁻³, which regime corresponds to the leftmost five data points in Fig. 1b. For all higher densities ρD³ > 0.01 we set our time step equal to Δt = 10⁻⁴ in order to ensure that we correctly integrate the equations of motion. This is especially necessary for the highest shear rate. We perform production runs of 10⁸ time steps, saving the particle positions every 10⁴ time steps. Before we initiate this production run, we equilibrate the system and ascertain that the pressure and total potential energy of the system have relaxed to a steady state.

By performing a full hierarchical clustering procedure for each saved simulation snapshot,^95,96 we determine the smallest connectivity length λ for which a cluster exists that connects to any image of itself through the periodic boundaries.⁹⁷ Using the set of all such connectivity lengths, we construct a percolation probability P(λ) for a given simulation. An example of such curves is given for four different system sizes in Fig. 3(h).

To obtain the percolation threshold λ_c, we locate of the intersection of the percolation probabilities of simulations of two different system sizes at the same particle density.⁹⁸ We find that finite-size effects are negligible if the sizes of the two systems are chosen such that N = 500 and N = 4000. In the absence of a shear flow [small gamma, Greek, dot above] = 0, our obtained percolation thresholds agree to within a few percent with those found in the literature.⁹⁹ In the presence of shear flow, the validity of this procedure depends on whether geometric percolation remains a critical phenomenon, which has been verified to be the case by Gallier and coworkers.³⁹ We probe the effect of the flow field on cluster shape by evaluating the eigenvalues of the gyration tensors of clusters of connected particles.¹⁰⁰

Author contributions

I. P., R. d. B., and P. v. d. S. designed the research and wrote the paper; I. P. performed the research.

Code availability

The code that is used to produce all data presented in this work is available in a permanent repository with DOI: https://doi.org/10.5281/zenodo.5786307.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

Ilian Pihlajamaa has been financially supported by the Dutch Research Council (NWO) through a Vidi grant, and René de Bruijn and Paul van der Schoot acknowledge funding by the Institute for Complex Molecular Systems at Eindhoven University of Technology.

Notes and references

1. C. Y. Li, L. Li, W. Cai, S. L. Kodjie and K. K. Tenneti, Adv. Mater., 2005, 17, 1198–1202.
2. M. Moniruzzaman and K. I. Winey, Macromolecules, 2006, 39, 5194–5205.
3. T. Wu, E. Chen, Y. Lin, M. Chiang and G. Chang, Polym. Eng. Sci., 2008, 48, 1369–1375.
4. K. H. Kim and W. H. Jo, Carbon, 2009, 47, 1126–1134.
5. J. A. King, M. D. Via, J. A. Caspary, M. M. Jubinski, I. Miskioglu, O. P. Mills and G. R. Bogucki, J. Appl. Polym. Sci., 2010, 118, 2512–2520.
6. B. Hornbostel, P. Pötschke, J. Kotz and S. Roth, Phys. Status Solidi B, 2006, 243, 3445–3451.
7. B. K. Satapathy, R. Weidisch, P. Pötschke and A. Janke, Compos. Sci. Technol., 2007, 67, 867–879.
8. S. Abbasi, P. J. Carreau, A. Derdouri and M. Moan, Rheol. Acta, 2009, 48, 943.
9. S. Stankovich, D. A. Dikin, G. H.-B. Dommett, K. M. Kohlhaas, E. J. Zimney, E. A. Stach, R. D. Piner, S. T. Nguyen and R. S. Ruoff, Nature, 2006, 442, 282–286.
10. M. C. Hermant, Manipulating the Percolation Threshold of Carbon Nanotubes in Polymeric Composites, PhD thesis, Eindhoven University of Technology, Eindhoven, 2009.
11. Z. Zeng, M. Chen, H. Jin, W. Li, X. Xue, L. Zhou, Y. Pei, H. Zhang and Z. Zhang, Carbon, 2016, 96, 768–777.
12. P. Avouris, M. Freitag and V. Perebeinos, Nat. Photonics, 2008, 2, 341–350.
13. E. Kymakis and G. A.-J. Amaratunga, Appl. Phys. Lett., 2002, 80, 112–114.
14. E. Kymakis, I. Alexandrou and G. A.-J. Amaratunga, J. Appl. Phys., 2003, 93, 1764–1768.
15. M. Sumita, K. Sakata, S. Asai, K. Miyasaka and H. Nakagawa, Polym. Bull., 1991, 25, 265–271.
16. V. Choudhary and A. Gupta, Carbon Nanotubes: Polym. Nanocompos., 2011, 65–90.
17. Y. Li and H. Shimizu, Polymer, 2007, 48, 2203–2207.
18. R. Ramasubramaniam, J. Chen and H. Liu, Appl. Phys. Lett., 2003, 83, 2928–2930.
19. L. Mathieu, P.-E. Bourban and J.-A. E. Månson, Compos. Sci. Technol., 2006, 66, 1606–1614.
20. X. Huang, Y. Li, F. Liu, P. Jiang, T. Iizuka, K. Tatsumi and T. Tanaka, IEEE Trans. Dielectr. Electr. Insul., 2014, 21, 1516–1528.
21. R. Ou, R. A. Gerhardt, C. Marrett, A. Moulart and J. S. Colton, Composites, Part B, 2003, 34, 607–614.
22. M. Majidian, C. Grimaldi, L. Forró and A. Magrez, Sci. Rep., 2017, 7, 1–9.
23. M. Kang, C.-H. Lim and J.-H. Han, Toxicol. Res., 2013, 29, 121–127.
24. Z. Spitalsky, D. Tasis, K. Papagelis and C. Galiotis, Prog. Polym. Sci., 2010, 35, 357–401.
25. F. Tiarks, K. Landfester and M. Antonietti, Macromol. Chem. Phys., 2001, 202, 51–60.
26. G. R. Cotten, Rubber Chem. Technol., 1984, 57, 118–133.
27. S.-P. Rwei, I. Manas-Zloczower and D. L. Feke, Polym. Eng. Sci., 1992, 32, 130–135.
28. N. Grossiord, P. J.-J. Kivit, J. Loos, J. Meuldijk, A. V. Kyrylyuk, P. van der Schoot and C. E. Koning, Polymer, 2008, 49, 2866–2872.
29. S. Torquato and H. W. Haslach Jr, Appl. Mech. Rev., 2002, 55, B62–B63.
30. A. Coniglio, U. De Angelis and A. Forlani, J. Phys. A: Math. Gen., 1977, 10, 1123.
31. J. Xu and G. Stell, J. Chem. Phys., 1988, 89, 1101–1111.
32. J. K.-G. Dhont, J. Fluid Mech., 1989, 204, 421–431.
33. U. Alon, I. Balberg and A. Drory, Phys. Rev. Lett., 1991, 66, 2879.
34. J. Bławzdziewicz and G. Szamel, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 1993, 48, 4632.
35. J. Bergenholtz, Curr. Opin. Colloid Interface Sci., 2001, 6, 484–488.
36. J. Vermant and M. J. Solomon, J. Phys.: Condens. Matter, 2005, 17, 187.
37. P. Butler, Curr. Opin. Colloid Interface Sci., 1999, 4, 214–221.
38. N. Koumakis, E. Moghimi, R. Besseling, W. C.-K. Poon, J. F. Brady and G. Petekidis, Soft Matter, 2015, 11, 4640–4648.
39. S. Gallier, E. Lemaire, F. Peters and L. Lobry, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 92, 020301.
40. P. G. de Gennes, J. Phys., 1979, 40, 783–787.
41. C. D. Lorenz and R. M. Ziff, J. Chem. Phys., 2001, 114, 3659–3661.
42. J. F. Thovert, V. V. Mourzenko and P. M. Adler, Phys. Rev. E, 2017, 95, 042112.
43. T. Hu and B. I. Shklovskii, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 74, 174201.
44. A. V. Kyrylyuk and P. van der Schoot, Proc. Natl. Acad. Sci. U. S. A., 2008, 105, 8221–8226.
45. G. Ambrosetti, C. Grimaldi, I. Balberg, T. Maeder, A. Danani and P. Ryser, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 155434.
46. B. Nigro, D. Ambrosetti, C. Grimaldi, T. Maeder and P. Ryser, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 83, 064203.
47. N. Wax, Selected papers on noise and stochastic processes, Courier Dover Publications, Mineola, 1954.
48. L. E. Reichl, A modern course in statistical physics, American Association of Physics Teachers, Maryland, 1999.
49. R. P. Brent, Algorithms for minimization without derivatives, Courier Corporation, Mineola, 2013.
50. J. P. Hansen and I. R. McDonald, Theory of simple liquids, Elsevier, Amsterdam, 1990.
51. H. H.-H. Homeier, S. Rast and H. Krienke, Comput. Phys. Commun., 1995, 92, 188–202.
52. T. DeSimone, S. Demoulini and R. M. Stratt, J. Chem. Phys., 1986, 85, 391–400.
53. R. Abe, Prog. Theor. Phys., 1959, 21, 421–430.
54. J. K.-G. Dhont and H. Verduin, J. Chem. Phys., 1994, 101, 6193–6205.
55. T. Ohtsuki, Phys. A, 1981, 108, 441–458.
56. S. A. Rice and J. Lekner, J. Chem. Phys., 1965, 42, 3559–3565.
57. E. Nazockdast and J. F. Morris, J. Fluid Mech., 2012, 713, 420.
58. J. F. Schwarzl and S. Hess, Phys. Rev. A: At., Mol., Opt. Phys., 1986, 33, 4277.
59. G. Szamel, J. Chem. Phys., 2001, 114, 8708–8717.
60. L. Banetta and A. Zaccone, Phys. Rev. E, 2019, 99, 052606.
61. N. J. Wagner and W. B. Russel, Phys. A, 1989, 155, 475–518.
62. J. F. Morris and B. Katyal, Phys. Fluids, 2002, 14, 1920–1937.
63. E. Nazockdast and J. F. Morris, Phys. Fluids, 2013, 25, 420–452.
64. K. Thøgersen, M. Dabrowski and A. Malthe-Sørenssen, Phys. Rev. E, 2016, 93, 022611.
65. Y. Yurkovetsky and J. F. Morris, J. Chem. Phys., 2006, 124, 204908.
66. F. Coupette, A. Härtel and T. Schilling, Phys. Rev. E, 2020, 101, 062126.
67. F. Coupette, R. de Bruijn, P. Bult, S. P. Finner, M. A. Miller, P. van der Schoot and T. Schilling, Phys. Rev. E, 2021, 103, 042115.
68. J. K.-G. Dhont, An introduction to dynamics of colloids, Elsevier, Amsterdam, 1996.
69. R. A. Lionberger and W. B. Russel, J. Rheol., 1997, 41, 399–425.
70. G. Szamel and J. A. Leegwater, Phys. Rev. A: At., Mol., Opt. Phys., 1992, 46, 5012.
71. N. J. Wagner, J. Non-Newtonian Fluid Mech., 2001, 96, 177–201.
72. J. F. Morris, Rheol. Acta, 2009, 48, 909–923.
73. J. F. Brady, Chem. Eng. Sci., 2001, 56, 2921–2926.
74. A. Tomilov, A. Videcoq, M. Cerbelaud, M. A. Piechowiak, T. Chartier, T. Ala-Nissila, D. Bochicchio and R. Ferrando, J. Phys. Chem. B, 2013, 117, 14509–14517.
75. R. C. Ball and J. R. Melrose, Adv. Colloid Interface Sci., 1995, 59, 19–30.
76. Y.-G. Tao, W. K. den Otter, J. K.-G. Dhont and W. J. Briels, J. Chem. Phys., 2006, 124, 134906.
77. G. Kwon, Y. Heo, K. Shin and B. J. Sung, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 011143.
78. A. E. Eken, E. J. Tozzi, D. J. Klingenberg and W. Bauhofer, J. Appl. Phys., 2011, 109, 084342.
79. J. Xu, W. Florkowski, R. Gerhardt, K.-S. Moon and C.-P. Wong, J. Phys. Chem. B, 2006, 110, 12289–12292.
80. I. Alig, T. Skipa, M. Engel, D. Lellinger, S. Pegel and P. Pötschke, Phys. Status Solidi B, 2007, 244, 4223–4226.
81. W. Bauhofer and J. Z. Kovacs, Compos. Sci. Technol., 2009, 69, 1486–1498.
82. J. K.-G. Dhont and W. J. Briels, Colloids Surf., A, 2003, 213, 131–156.
83. M. Ripoll, P. Holmqvist, R. G. Winkler, G. Gompper, J. K.-G. Dhont and M. P. Lettinga, Phys. Rev. Lett., 2008, 101, 168302.
84. R. Schueler, J. Petermann, K. Schulte and H. P. Wentzel, J. Appl. Polym. Sci., 1997, 63, 1741–1746.
85. D. G. Duffy, Green's functions with applications, CRC Press, Boca Raton, 2015.
86. D. E. Elrick, Aust. J. Phys., 1962, 15, 283–288.
87. K. Cole, J. Beck, A. Haji-Sheikh and B. Litkouhi, Heat conduction using Greens functions, Taylor & Francis, New York, 2010.
88. V. I. Lebedev, USSR Comput. Math. Math. Phys., 1976, 16, 10–24.
89. T. Pang, An introduction to computational physics, Cambridge University Press, Cambridge, 2012.
90. C. G. Broyden, Math. Comput., 1965, 19, 577–593.
91. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in't Veld, A. Kohlmeyer, S. G. Moore and T. D. Nguyen, et al. , Comput. Phys. Commun., 2021, 271, 108171.
92. M. Doi, Soft matter physics, Oxford University Press, Oxford, 2013.
93. J. D. Weeks, D. Chandler and H. C. Andersen, J. Chem. Phys., 1971, 54, 5237–5247.
94. A. W. Lees and S. F. Edwards, J. Phys.: Condens. Matter, 1972, 5, 1921.
95. S. C. Johnson, Psychometrika, 1967, 32, 241–254.
96. D. Müllner, et al. , J. Stat. Softw., 2013, 53, 1–18.
97. D. C. Rapaport, The art of molecular dynamics simulation, Cambridge University Press, Cambridge, 2004.
98. J. Škvor, I. Nezbeda, I. Brovchenko and A. Oleinikova, Phys. Rev. Lett., 2007, 99, 127801.
99. M. A. Miller, J. Chem. Phys., 2009, 131, 066101.
100. R. Sanchez and P. Bartlett, J. Phys.: Condens. Matter, 2005, 17, S3551.

---

Footnote

† By self-diffusion constant, we mean the diffusion constant that a particle would have if it is isolated from all other solute particles in the absence of a flow field. It can be estimated by the Stokes-Einstein relation , in which η₀is the viscosity of the medium.

---