Hydrodynamic simulations of sedimenting dilute particle suspensions under repulsive DLVO interactions

We present guidelines to estimate the effect of electrostatic repulsion in sedimenting dilute particle suspensions. Our results are based on combined Langevin dynamics and lattice Boltzmann simulations for a range of particle radii, Debye lengths and particle concentrations. They show a simple relationship between the slope $K$ of the concentration-dependent sedimentation velocity and the range $\chi$ of the electrostatic repulsion normalized by the average particle-particle distance. When $\chi \to 0$, the particles are too far away from each other to interact electrostatically and $K=6.55$ as predicted by the theory of Batchelor. As $\chi$ increases, $K$ likewise increases as if the particle radius increased in proportion to $\chi$ up to a maximum around $\chi=0.4$. Over the range $\chi=0.4-1$, $K$ relaxes exponentially to a concentration-dependent constant consistent with known results for ordered particle distributions. Meanwhile the radial distribution function transitions from a disordered gas-like to a liquid-like form. Power law fits to the concentration-dependent sedimentation velocity similarly yield a simple master curve for the exponent as a function of $\chi$, with a step-like transition from 1 to 1/3 centered around $\chi = 0.6$.


Introduction
The physics of sedimenting particles have proven to be surprisingly difficult to model despite many attempts over a large fraction of the 20th century. While a single particle slowly sedimenting in a sufficiently large container can be easily described by Stokes' law, the long-ranged nature of hydrodynamic interactions renders the dependence of the sedimentation speed on the particle concentration complicated to derive even in the dilute limit.
For the purpose of brevity we refer to the case of uncharged particles interacting only via hydrodynamic and hard sphere interactions as the case of non-interacting particles throughout this paper. The theory of non-interacting particles reached a major breakthrough when in 1972 Batchelor 1 derived the sedimentation velocity v at small particle volume fractions φ relative to the velocity v 0 at infinite dilution as with K = 6.55. Similarly, the sedimentation velocity is sometimes written as which is identical to Eq. (1) in the limit of small φ . The sedimena Helmholtz Institute Erlangen-Nürnberg for Renewable Energy, Forschungszentrum tation velocity remains positive for all φ following Eq. (2), unlike Eq. (1), though neither equation is accurate anywhere near the concentration φ ≈ 15% where Eq. (1) goes to zero. Beyond the dilute limit, the Rotne-Prager far-field approximation of hydrodynamic interactions was shown by Brady and Durlofsky in 1988 to be accurate for non-interacting spheres even up to a volume fraction of 50%. 2 Alternatively, via a Virial expansion, Cichoki et al. attempted to take into account three-particle contributions. 3 Experimentally, good agreement with Eq. (1) has been shown to require a Péclet number Pe < 1 due to its underlying assumption of a perfectly homogenous radial distribution function (RDF). 4 Especially for particles of nanometer scale neglecting any non-hydrodynamic interparticle interactions is a strict limitation though. Indeed, depending on the pH value, most types of colloidal particles tend to accumulate considerable surface charges when dissolved in water [5][6][7] . This leads to strong electrostatic interactions which typically decay over a Debye length of the order of 10 nm. The Debye length in water can in principle reach hundreds of nanometers, though this requires high degrees of purity that are in practice difficult to achieve. For this reason a majority of studies on the sedimentation of interacting particles focus on attractive potentials. 8,9 In organic solvents such as ethanol, however, Debye lengths of around 800 nm have been reached in experiments. 10 For particle suspensions with strong electrostatic interactions and weak screening (i.e. a Debye length λ D large enough to be comparable to the average particle-particle distance) a strongly nonlinear decrease of the sedimentation velocity with concentration has been both predicted theoretically and observed experimentally 10,11 even in the dilute limit where φ < 1%.
Early studies of electrostatic effects in particle sedimentation include the work of Booth 12 in 1954. They developed the dipole moment of sedimenting charged particles as a power series in terms of the particle charge or zeta potential and managed to calculate the first two coefficients in the series. The theory is thus appropriate for sufficiently low surface charges/zeta potentials, although this limitation was removed in a numerical extension of Booth's work by Stigter in 1980. 13 Both Booth's and Stigter's theories completely neglect hydrodynamic interactions between the particles and do not take changes in the RDF of the suspension into account.
A number of studies [14][15][16][17] of sedimentation under both electrostatic and hydrodynamic particle-particle interactions have been performed using methods based on geometric cells to obtain the hydrodynamic component, either with the free-surface boundary condition by Happel 18 or the zero vorticity condition by Kuwabara. 19 While experimental results confirm the cell models as adequate to calculate the sedimentation potential, 20 both the method by Happel and that of Kuwabara fail to correctly reproduce the sedimentation behavior of non-interacting particles in the dilute limit found by Batchelor about 14 years after the introduction of the method. 1,18,19 Furthermore, the methods based on geometric cells cannot take into account changes in the RDF of the sedimenting suspension induced by the electrostatic interactions and they assume an electrically neutral unit cell, which may be a too rough simplification if Debye layers overlap strongly. 21 Another promising approach in modeling charged particle sedimentation numerically was taken by Watzlawek and Nägele,22 though their approach is limited by the fact that it can only take into account pair-wise hydrodynamic interactions. Neglecting many-body hydrodynamic interactions was shown by Brady and Durlofsky 2 to lead to a significant error in the sedimentation rate at volume fractions as low as 5%, though the result could be improved considerably by additionally neglecting stresslet contributions as per the Rotne-Prager approximation. Approximate many-body hydrodynamic interactions can be taken into account using the Stokesian dynamics method 23 and advancements in recent years have improved its performance up to a linear scaling with the number of particles. 24 Nonetheless, the handling of hydrodynamic interactions remains fundamentally approximate in Stokesian dynamics due to a truncated expansion of the mobility matrix. Furthermore, the method is limited in terms of its extensibility to non-zero Reynolds numbers and polydisperse or non-spherical particles. Parallelized Stokesian dynamics implementations scale efficiently to up to a few hundred CPUs 25 and have been used to study the sedimentation of aggregates of thousands of polydisperse particles. 26 Banchio et al. and Gapinski et al., [27][28][29] have previously employed the Stokesian dynamics method to numerically study suspensions under repulsive interactions. They obtained the structure factor of the suspension and the so-called hydrodynamic function H(q) for selected values of salt and particle concentrations. Though their results are focused more on modelling diffusivity, the hydrodynamic function contains the relative sedimentation speed of the suspension under a spatially constant force for q = 0. Comparison of experiments with the hydrodynamic function for a given concentration and as a function of q requires measuring the static structure factor, e.g. via X-ray scattering, as well as the collective diffusion function, e.g. via dynamic light scattering. Our approach of quantifying the functional shape and the mean slope of the sedimentation velocity as a function of con-centration for a broad range of salt concentrations and different particle concentration ranges should lend itself to a more straightforward comparison to centrifugal sedimentation experiments. In fact we have recently applied an early version of our method described in this work in an experimental context. 30 As an alternative to Stokesian dynamics one can model the sedimentation of particles in a fluid by coupling the discrete element method for the dynamics of the suspended particles to a Stokes or Navier-Stokes level hydrodynamics solver. Many different methods have been used for the latter, such as directly solving the Navier-Stokes equation using the finite element method, 31 smoothed particle hydrodynamics, 32 or stochastic rotation dynamics. 33,34 In this work we employ the lattice Boltzmann method (LBM). It has been shown to be a viable tool to capture the full hydrodynamic interactions of large numbers of noninteracting sedimenting particles by Nguyen and Ladd in 2005, 35 though there is similar work by Ladd with smaller particle numbers dating back to 1994. 36 Later on the method has similarly been used to model particles with attractive interaction potentials. 37 Several different algorithms for coupling particles to the LBM fluid exist, the method is numerically efficient and is not limited to low Reynolds number flows. 38 For low Reynolds number flows the LBM has been found to give results consistent with the Stokesian dynamics method. 39,40 In this paper we numerically study the impact of electrostatic interactions modeled by DLVO theory on sedimenting suspensions under varied particle size, concentration and Debye length. By simulating the interactions of a large number of particles and the resulting changes in the RDF explicitely and by including full hydrodynamic interactions using the LBM we improve upon previous studies and contribute to a clearer picture of how electrostatic interactions influence particle sedimentation.

Model and methodology
Each sedimentation simulation for a given set of concentration, particle size, and Debye length parameters consists of two major steps. First, a set of particle positions representative of an equilibrated bulk suspension of charged particles interacting via DLVO potentials is generated from a Langevin dynamics simulation. Second, the hydrodynamic interactions and the resulting sedimentation velocity under added constant acceleration (representing gravitational or centrifugal forces) are calculated for the particle positions obtained previously using the LBM. The final result is the particle velocity in the direction of the constant acceleration averaged over all particles.
While in the first step both particle positions and particle velocities evolve in time, only the velocities are updated in the last step while the positions remain fixed. In this way we neglect changes in the RDF induced by hydrodynamic interactions and greatly reduce the convergence time and numerical cost of our hydrodynamic simulations. We consider this simplification to be justified in the limit of small Péclet numbers, where particle advection plays a small role compared to diffusion and drift fluxes induced by strong DLVO interactions. As we keep the particle positions fixed, neither advection nor diffusion occur in our hy-drodynamic simulations so that the Péclet number is not obviously defined. However, by keeping the Reynolds number small (Re 5 × 10 −6 ) we can consider the fluid flow velocity around the particles to be arbitrarily small. It follows that the Péclet number Pe ∝ Re/D calculated using the diffusivity D of the particles in the preceding Langevin simulation is likewise vanishingly small.
Ongoing research into the possible causes of an observed slow decay of sedimentation velocity fluctuations has led to the widespread assumption that subtle changes in the RDF may be taking place in sedimenting suspensions over long time spans up to several hours, even in the limit of small Pe. 41,42 Reproducing this experimentally observed decay of velocity fluctuations accurately would require significantly longer simulation times 43 and the presence of confinement 44 with a geometry matching the experimental system. 45 To our knowledge, however, no corresponding long term evolution of the mean sedimentation velocity has been observed so far in monodisperse suspensions.

Generating particle positions
In the first step, we initialize about 10 000 spherical particles with random positions r i without overlap in a 3D rectangular system with periodic boundary conditions. The particle positions are evolved in time t in each spatial dimension ι according to the Langevin equation The particle mass m is set to reproduce a particle density of 1800 kg m −3 , which is a realistic value for e.g. SiO 2 nanoparticles. Pairwise particle interaction forces F i (r i , r j ) account for DLVO and hard sphere interactions with all surrounding particles up to a cutoff radius carefully selected depending on the range of the DLVO interactions. Stokes' law provides the translational friction coefficient γ = 6π µR based on the dynamic viscosity µ. The randomized force η represents thermal fluctuations and fulfills the fluctuation dissipation theorem in each spatial dimension, which is given by The Langevin equation is discretized in time and the particle positions are updated according to the leapfrog algorithm. Convergence of the Langevin dynamics simulations is determined based on the time evolution of the total DLVO interaction energy in the system. When the drift in the energy over the last 5000 time steps is smaller than the standard deviation of the energy due to thermal fluctuations, the simulation is stopped. The final particle positions are then transferred to a lattice Boltzmann (LB) simulation to determine the hydrodynamic interaction of the particles. For simulations with no DLVO interactions, the Langevin dynamics simulations are skipped and random particle positions are used in the LB simulation.

Particle-fluid coupling
LB simulations are performed using our in-house code LB3D. [46][47][48] In our LB implementation, fluid properties are calculated on a regular cubic lattice in three dimensions. On each lattice site 19 scalar populations f i are defined and the lattice constant is referred to as ∆x. Each population is proportional to the fraction of fluid flowing in the velocity direction c i toward a neighboring lattice site or remaining at rest (c 19 = 0) during that time step. In each time step ∆t they are shifted to a neighboring lattice site according to and then relaxed toward an equilibrium distribution f eq i in the collision step This approach, using a single (dimensionless) relaxation time τ r , is known as the BGK scheme, after Bhatnagar, Gross and Krook. 49 . f eq i is a truncated Maxwell-Boltzmann distribution for the discretized set of possible velocities along the velocity directions c i . 50,51 The source term S i stems from the action of body forces. A number of different schemes to calculate S i have been shown to produce physically accurate results. For this work we choose a scheme by Kupershtokh. 52,53 Letting ∆m be the unit mass, the fluid mass density ρ f and velocity v f at a given lattice site are calculated as and the dynamic viscosity is given as We performed all our simulations with τ r = 1 for reasons of numerical simplicity. Changing the relaxation time and thus the viscosity only changes the overall time scale of the system. Particles are coupled to the interpolated fluid velocity via a linear friction force in an approach based on work by Ahlrichs and Dünweg. 54 According to Stokes' law, the friction force experienced by a single particle of velocity v p inserted into a fluid flowing with velocity v f is The same friction force, with the opposite sign, also acts on the fluid following Newton's third law. The effect of each particle on the fluid flow is limited to the effect of this point-like friction force, yielding the Stokeslet approximation of its full flow field. As a result, particle rotation is not included in the model. Applying Eq. (10) in the numerical model using γ = 6π µR and setting v f equal to the fluid velocity from the LBM interpolated to the particle position results in steady-state velocities v p = |v p | that are higher than the expected result from Stokes' theory. This is because v f in Eq. (10) represents the fluid velocity without the Stokeslet contribution from the considered particle according to Stokes' law. We chose a particle radius equal to the LB grid spacing in order to have a relatively large radius while ensuring that the particle geometry remains comfortably within the extent of its stencil surrounding it. Furthermore we verified that the friction force densities remain smaller by more than an order of magnitude at all times compared to values deemed problematic in the lattice Boltzmann method. Fortunately, the contribution of the particle to its surrounding flow field can be easily subtracted by rescaling the friction coefficient, as shown by Ollila et al. 55 The correction factor γ s can be determined analytically in principle, 56 but it depends non-trivially on the stencil used to interpolate the fluid velocity as well as other details of the numerical implementation. Instead we choose the simpler approach of deriving γ s from fits to a series of numerical measurements of the steady-state single particle velocity as a function of the input friction coefficient. 55 Using a cubic stencil with a side length of four lattice discretization lengths ∆x and a weighting function derived by Peskin, 57 we obtain γ s ≈ 4.91 in simulation units. The corresponding fit is shown in Fig. ?? in the ESI † . From here on, γ always refers to the corrected friction coefficient according to the substitution in Eq. (11).
The sedimentation of particles with mass m in our LB simulations is triggered by a constant force F g representing gravitational or centrifugal acceleration as well as the counteracting buoyancy. The same force F g with opposite sign is distributed homogenously among all fluid sites in the system. This ensures global momentum conservation and mimics the backflow of displaced fluid occuring during sedimentation in a closed cell.
Assuming a constant v f , the particle velocity update by one time step due to the friction force alone can be written as If ∆t γ m < 1, v p approaches v f via an exponential decay. If 1 < ∆t γ m < 2, v p oscillates around v f due to discretization errors, but |v ∆ | still decays to zero in time. If, however, ∆t γ m > 2, then |v ∆ | diverges to infinity in an oscillating manner. The easiest way to avoid these discretization effects would be to choose the time step such that ∆t < m/γ, or, at least, ∆t < 2m/γ. However, large values of both γ and ∆t are desirable when simulating a suspension at low Reynolds number. In order to avoid this issue, we analytically integrate the friction force F s under the assumption of a constant v f but continuously varying v p and F s over one time step, add the constant F g , and calculate the average total force F T ∆t as Because we use the leapfrog algorithm to generate particle tra- ) shifted by half a time step with respect to positions and forces is readily available. The fluid velocity in Eq. (13) is ) because we require the fluid velocity averaged over the time step from t i− 1 2 to t i+ 1 2 . In the overdamped limit, when m/γ ∆t, Eq. (13) gives the same acceleration from F g as predicted by Brownian dynamics, plus advection by v f .
The averaged friction force acting on the fluid can be identified as −( F T ∆t − F g ) and it is distributed to the fluid sites surrounding the particle on the same stencil on which the interpolation of v f takes place.
The LB simulations are considered converged when the slope of the sedimentation velocity over time relative to the velocity at infinite dilution and averaged over all particles and the last 1000 time steps falls below a threshold value of 5 × 10 −8 . This procedure usually requires between 5000 and 20 000 LB time steps. We find that letting some simulations run up to about thirty times longer changes the final sedimentation velocity by less than 0.01%.

DLVO interactions
The total force acting on a particle in the Langevin simulations is calculated with the same averaging of the friction force introduced in Eq. (13) but setting the fluid velocity v f to zero and exchanging F g with a sum of DLVO and hard sphere pair poten- . Each particle interacts only with particles within a numerical cutoff distance chosen according to the DLVO parameters of the simulation. The DLVO interactions consist of an attractive contribution stemming from van der Waals interactions and a repulsive contribution stemming from Coulomb repulsion screened by counterions: The van der Waals force of two spheres of equal radius R at a surface to surface distanceŝ = s/R in multiples of R is the derivative of the potential 58 We model van der Waals forces using an effective Hamaker constant of A H = 2 × 10 −20 J. This value is similar to that measured by Fielden et al. 59 for a silica particle interacting with a partially oxidized silicon wafer. Valmacco et al. 60 measured substantially lower values for pairs of silica particles in water, probably due to a high surface roughness. A Hamaker constant of the order of 1 × 10 −20 J is to be expected for interactions between polystyrene particles in water. 61 As shown by example in Fig. ?? of the ESI † , the strength of the repulsive component of the DLVO interaction in the parameter space of large Debye lengths studied by us renders the van der Waals interaction almost irrelevant for most of our simulations. We include van der Waals forces anyway for the sake of completeness. The repulsive component consists of a Coulomb interaction between like-charged spheres with an electrostatic potential ζ at the hydrodynamic slipping plane, which is exponentially screened over a decay length λ D by the presence of dissolved ions in a solvent of dielectric permittivity ε 58 A comparison of the resulting total DLVO potential E DLVO = E vdw + E coul with E coul alone for R = 300 nm, ζ = 50 mV and dif- ferent values of λ D is shown in Fig. ?? in the ESI † .
The simplified pair-wise interactions of DLVO theory are computationally efficient and allow us to reach large particle numbers with acceptable computational effort. However, this approach neglects the deformation of the Debye layer in the presence of a hydrodynamic flow. While taking this deformation into account could be achieved by coupling the solver for the fluid and particle dynamics to a solver for the Nernst-Planck equation, 62,63 the influence of such ion advection effects becomes negligible when the ions' Péclet number λ D v/D i is small. 64 As established in section 2, we are concerned in this work with systems of small particle Péclet number and Debye lengths comparable in size to the particle radius. The ions' Péclet number can be considered to be smaller still, owing to the smaller size and therewith larger diffusivity D i of the ions as compared to the particles. A fully resolved double layer would furthermore yield a reduction of the sedimentation velocity due to the restoring dipole force acting on the particle when it is accelerated by F g out of the center of its ionic atmosphere. 13,65,66 Because this so-called primary charge effect is also present in the sedimentation of a single particle, we assume its effect on the relative sedimentation speed v/v 0 to be negligible.
In order to avoid strongly overlapping particles due to the divergence of E vdw at contact when λ D is small and thermal fluctuations allow particles to cross the potential barrier posed by E coul , a hard sphere repulsion term of the form based on Hertzian contact theory 67 is applied to particles at center-to-center distances c < 2R. The stiffness k is chosen empirically based on the conditions that it needs to be sufficiently large to avoid significant particle overlap but small enough to not lead to excessive particle acceleration due to time discretization.

Validation
In order to check the accuracy of the particle-fluid coupling, we compare our simulations with known results for the sedimentation behavior of non-interacting particles. First we compute the velocity of a pair of neighboring particles under constant acceleration in Stokes flow. Two particles with a radius equal to the length of discretization of the LB solver are initialized in a fully periodic system. As described in section 2.2, F g is applied to each particle in the same direction and −2 F g is spread homogenously over all fluid lattice sites. The component of the final sedimentation velocity in direction of F g and relative to the velocity of a single particle can be written as where λ 1 and λ 2 as such are given in tabulated form as a function of the interparticle distance by Batchelor, 1 albeit the original computations were performed by Stimson and Jeffery 68 for λ 1 , and Goldman et al. 69 for λ 2 . Here, Θ is the angle between the connecting line of the particle centers and the direction of F g . As shown in Fig. 1(a), very good agreement with Eq. (17) is obtained even when the interparticle distance from center to center is less than 3 discretization lengths. This is remarkable, as the limitation of the fluid-particle coupling to the Stokeslet level means that hydrodynamic interactions are strictly accurate only in the far-field. Next, we benchmark Eq. (1) for the sedimentation velocity of non-interacting particles in bulk by simulating about 10 000 sedimenting particles in the same way as in the previous test. The corresponding results in Fig. 1(b) also show good agreement with Eq. (1), with a measured K = 6.10 ± 1.24. Fig. 1(c) shows results from an identical set of sedimentation simulations with radius R = 100 nm and a Debye length λ D = 950 nm. The long-ranged repulsive interactions in these simulations lead to a functional form of v(φ ) ∝ 3 √ φ , in agreement with theoretical expectations. 10

Results
We perform simulations for a range of Debye lengths from 5 nm to 950 nm with particle radii R set to 100, 200, 300, 450, and 600 nm. A Debye length of λ D = 5 nm corresponds to a monovalent salt concentration of about 3.7 mM in water at room temperature, whereas 950 nm is close to the Debye length obtained in perfectly pure water solely by self-dissociation at a pH value of 7. We keep the zeta potential fixed at ζ = 50 mV regardless of the particle size.
For the smallest Debye length of 5 nm combined with the largest particles of R = 600 nm, strong aggregation occurs, leading to neg-  ative values of K, as shown in Fig. 2. For smaller particles at the same Debye length, we observe K ≈ 6.55 and almost no aggregation. To understand this, first note that the van der Waals potential in Eq. (14) does not depend on R for a givenŝ. The repulsive potential in Eq. (15) on the other hand can be shown in a simple mathematical exercise to always decrease when R is increased as long as s > λ D . The proof can be found in section ?? of the ESI † . Thus, the attractive potential at distances beyond one Debye length is relatively stronger than the repulsion for larger particles. We exclude simulations showing extensive aggregation from further analysis.
For small Debye lengths around 10 nm, the sedimentation velocity is predicted well by Eq. (1) with K ≈ 6.55. As λ D increases, the slope increases rapidly, meaning that mutual hindrance is increased. While particles close to each other sediment faster than a single particle, as shown in Fig. 1(a), at larger interparticle distances the effect of fluid backflow dominates and particles mainly slow down each other. An increase in λ D leads directly to an increase in the mean distance between next neighbors due to a longer range of the repulsive potential.
When λ D is sufficiently large so that the particles cannot fully escape the repulsive potential of their neighbors, v(φ ) becomes distinctly nonlinear. The nonlinear regime begins to show in Fig. 2 for λ D = 950 nm. According to calculations by Thies et al., 10 the sedimentation behavior in the limit of large λ D should follow with ς ≈ 1.8 and ω ≈ 3. This functional form is quite general for ordered particle arrays. 10,22,70 The slope obtained from a linear fit over the concentration range [φ 1 , φ 2 ] in a system described by Eq. (18) can be predicted as While the more general Eq. (18) can give a better fit in the limit of large λ D , a simpler linear fit may often be preferable, particularly when working with data exhibiting significant statistical errors and when varying φ across a relatively narrow range.
When λ D is small enough for v(φ ) to remain in the linear regime, K can be approximated for general interaction potentials Φ following Batchelor and Wen 71 as K = 6.55 − 0.44α with If the potential Φ in Eq. (20) is set equal to E coul + E vdw + E hs , the resulting K depends on R, ζ , λ D , and k in a non-trivial way.
One might be tempted to identify ξ = λ D R . However, as both the radius R and the zeta potential ζ influence the strength of the repulsive potential at a given distance s/R significantly (see Fig. ?? in the ESI † ), we instead define ξ as the surface-to-surface distance in multiples of the radius at which E DLVO first exceeds k B T coming from infinity, thus constituting a significant potential barrier as compared to thermal energy. The exact choice of the threshold value makes relatively little difference in the resulting value of ξ due to the fast exponential decay of the repulsive potential. Because this measure depends on the charge state of the particle as well as its size, it encodes more information than the Debye length alone. In the dilute regime that we are focused on here, colloids are typically far apart. Hence, we can disregard the van der Waals contribution, Eq. (14), and we can approximate the DLVO potential by the solely electrostatic interaction, Eq. (15). Accordingly, the surface-surface dimensionless distance, ξ = ξ 0 , at which the colloid-colloid DLVO interaction is comparable to the thermal energy is obtained by numerically solving Fig.3a shows the dependence of ξ 0 on λ D for diverse particle radii. In particular, Fig.3a shows that ξ 0 can be comparable or even larger than the particle size.
Eq. (21) models the effect of the repulsive potential as an excluded volume around otherwise non-interacting and thus randomly distributed particles. A similar approach of modelling short-ranged DLVO interactions as an excluded volume, or alternatively an effective particle concentration, has been used previously for example by Gilleland et al. 11 or Antonopoulou et al. 72 The impact of the repulsive potential barrier atŝ = ξ 0 on the final particle distribution of course depends on the average particle-particle spacing, which in turn depends on the particle concentration φ . In order to account for this we furthermore in- troduce the naively calculated average interparticle spacinĝ using the particle volume V p . It is formulated in multiples of the radius and measured from surface to surface, just like ξ . Normalizing ξ 0 as we obtain a useful dimensionless measure for the range of the repulsive DLVO force relative to the average interparticle distance.
Interestingly, Fig.3b shows that, for the values of the parameters under scrutiny, ξ 0 attains quite large values as compared to both λ D /R. This implies that the relevant distance at which the colloids experience the mutual DLVO interaction can be quire larger than the Debye length hence supporting our definition of χ 0 in Eq. 24.
We have used the surface-surface distance ξ 0 (see Eq. (22)), and the associated value of χ 0 (see Eq. (24)), as the effective particle size in the hard-sphere model, Eq. (21). However, the agreement is qualitatively good yet we admit some quantitative discrepancies. To address the role of the softness of the DLVO potential at distance ξ 0 , as compared to the hard-sphere interaction, we define ξ (and the associated χ) as the distance at which the DLVO potential 10k B T by numerically solving Indeed, at such distance the DLVO potential is stiff and therefore it may resemble the hard-sphere interaction. Interestingly, the slope K as a function of χ in Fig. 4(a) neatly collapses onto a single curve when calculated over a fixed range of concentrations hence showing that χ is the dimensionless variable that captures the "dynamics".
The error bars in Fig. 4(a) account for variations due to the randomness involved in initial particle placement and the subsequent equilibration of particle distributions under thermal fluctuations. To estimate the error bars we repeat simulations up to 6 times at selected parameter combinations spanning the whole range of χ with different random number seeds and calculate the standard deviation of the resulting velocities as described in section ?? of the ESI † . The error bars strongly depend on χ and are largest for non-interacting particles. Knowing χ and the concentration at which K is measured, we can predict the value of K to a decent accuracy both for small (via Eq. (21)) and large χ (via Eq. (19)). At intermediate χ an interpolating fit drawn in dashed lines in Fig. 4(a) matches the observed trend well. This interpolating fit captures the transition from K Φ (ξ ) (Eq. (21)) to the constant value K ω given by Eq. (19) via a sigmoid function We remark that for χ → χ m Eq. (27) gives K K Φ (ξ ) whereas for χ → ∞ we get K K ω . We determine the two fit parameters χ m ≈ 0.38 -roughly corresponding to the position of the maximum, and δ K ≈ 0.096 -giving the scale in χ over which the transition to a locally ordered suspension occurs, from a combined fit to the complete data set for both φ = 0.2 − 0.8% and φ = 1 − 1.4%. The dashed lines in Fig. 4(a) both follow Eq. (27) using the same values of χ m and δ K . The two lines differ only due to the different values of φ1 / 2 used in calculating K ω and the different values of ξ = χŝ φ inserted in K Φ for a given χ, again due to the different volume fractions φ . In calculating K ω we set ω = 3 and ς ≈ 1.71 regardless of φ . We remark that the value ς ≈ 1.71 has been obtained by averaging over all data points from nonlinear fits to Eq. (18) at χ > 1.2, as indicated in Fig. 4(b).
One can reformulate the fitted K(χ) from Eq. (27) as a function of φ for fixed ξ and perform numerical integration to reconstruct the hindrance function v(φ )/v 0 as shown in section ?? in the ESI † . As shown in Fig. 4(a) we recover the case of non-interacting par- ticles for χ → 0 and K → 6.55 as in Eq. (1). Up to χ ≈ 0.3, K is well-approximated by Eq. (21), which is shown as full lines in Fig. 4(a). Eq. (21) fails as a valid approximation when the particle distribution cannot be approximated as homogenous in space, i.e. when the RDF deviates from the step function expected for dilute hard spheres with an effective radius increased by Rξ /2. Fig. 4(b) shows the obtained parameters ς and ω from nonlinear fits to Eq. (18). While the nonlinear fit works well in the locally ordered regime at χ 0.4 and ω ≈ 1 is correctly reproduced even for χ → 0, there is a large uncertainty in ς at χ 0.4. As shown in Fig. 1(b) and (c), the uncertainty in the velocity in the disordered regime is much larger than in the locally ordered regime. According to Eq. (18), v/v 0 depends much more sensitively on the exponent ω than on the prefactor ς , in particular when φ is small. This can be seen from the ratio of the derivatives (∂ v/∂ ς )/(∂ v/∂ ω) = −ω 2 /(ς ln φ ), which goes to zero for small φ . Accordingly, uncertainties in the velocity translate into much larger uncertainties in the values of ς than of ω, leading to the large spread in the obtained values of ς at low χ in Fig. 4(b). The exponent ω can be predicted well from χ via a fitted sigmoid function with the fit parameters χ 0 ≈ 0.63 and δ ω ≈ 0.13 again obtained via a single fit to the combined data for all φ as in Eq. (27). Fig. 5(a) shows the changes in the RDF leading in turn to the changes in K and ω. For small χ the RDF is a simple step function as expected for randomly distributed hard spheres. The first deviation from this idealized form is visible for χ = 0.4 in the form of a pronounced primary maximum next to the exclusion zone. The transition from a disordered to a locally ordered particle distribution is accompanied by an oscillatory component in the RDF, which becomes visible at χ ≈ 0.5. The length scale over which the oscillations decay can be interpreted as the length scale over which particle positions are correlated. Predictably, this length scale increases markedly as χ increases, with the RDF for χ = 1.7 showing visible correlation at distances well beyond 25 particle radii. Our results here qualitatively agree well with the RDF of charged sphere suspensions obtained in other works. 10,11,28,29 Going to higher values of φ or χ than those we studied should lead to crystalline bcc or fcc particle distributions. 73,74 Simulating crystalline systems would require great care though, because their very long-ranged order may be strongly affected by finite system sizes and take a long time to equilibrate. 75 Furthermore, DLVO force models may be ill-suited for such systems, as they fail to properly model the experimentally observed coexistence of colloidal crystals with disordered phases in dilute suspensions at small salt concentration. 76 In Fig. 5(b) we can see how changes in χ affect the average surface-to-surface interparticle distance s nn as well as its standard deviation ∆s nn within the particle configuration used in a given hydrodynamic simulation. Due to normalization by the average interparticle spacing s φ the results for different φ again collapse rather well on a single curve. Unlike in Figs. 4(a) and (b), where each data point represents a group of simulations at different φ and otherwise identical parameters (with χ being averaged over φ ), each data point here corresponds to a single simulation.
The full and dashed lines in Fig. 5(b) compare the simulation results with a homogenous suspension of particles interacting only via the step potential Φ = E 0 Θ(ξ −ŝ) with E 0 → ∞ at s = ξ -like hard spheres with a radius enlarged in proportion to χ. Because the particle distribution (derived in section ?? of the ESI † ) neglects particle-particle correlations beyond the range of the step potential, it is neccessarily inaccurate when either φ or χ are large.
A substantial discrepancy between the simulation results for s nn and the hard sphere distribution develops starting around χ = 0.3. This is consistent with the observation that the enlarged hard sphere model from Eq. (21) predicts K(χ) well only up to Distance to next-neighbor-particle s nn and its standard deviation ∆s nn = Var(s nn ) normalized by s φ . Black lines give theoretical results for homogenous particle distributions with an infinite step potential effectively extending the hard sphere radius in proportion to χ by ξ R/2. χ ≈ 0.3, as shown in Fig. 4(a). Interestingly, ∆s nn seems to diverge from the hard sphere distribution much faster though, showing that the microstructure of the DLVO suspension does differ noticeably from the hard sphere suspension for χ < 0.3, despite affecting sedimentation in much the same way as an increased hard sphere radius.
For χ 1, s nn approaches the maximal average interparticle distance s φ and ∆s nn indicates a narrow distribution of nextneighbor distances as expected for a locally ordered particle distribution.

Conclusions
By simulating the hydrodynamic and DLVO interactions of large ensembles of particles we found seemingly universal trends in the sedimentation behavior for a wide range of Debye lengths and particle sizes. We quantified the effect of particle interactions depending on the range χ of electrostatic repulsion in our results either via the slope K or the exponent 1/ω of fits to the sedimentation velocity v(φ )/v 0 across different ranges of φ .
K is the slope extracted from a linear fit to v(φ )/v 0 and appears to be described well by our fit to Eq. (27) for any φ in the dilute limit. Eq. (27) predicts K(χ) assuming that the electrostatic repulsion at χ 0.3 acts merely like an increase of the effective hard sphere radius, whereas at χ 1 sedimentation follows the known solution Eq. (18) with ς ≈ 1.71 and ω = 3 for ordered particle arrays. The transition from one solution to another is approximated in Eq. (27) using a simple sigmoid function.
Applying non-linear fits following Eq. (18) to v(φ )/v 0 instead we find clear nonlinearity (ω > 1) commencing around χ = 0.4, where K(χ) reaches its maximum. Near this transition point from linear to nonlinear the RDF shows a transition from a disordered gas-like state to a liquid-like state. ω(χ) is likewise describable by a sigmoid function, with a smoothened step-like transition from ω = 1 at χ → 0 to ω = 3 at χ 1. This coincides with the point where the average next-neighbor distance reaches its maximum possible value.
Both K(χ) following Eq. (27) and ω(χ) following Eq. (28) offer themselves as a potentially useful gauge to estimate the extent of electrostatic interactions (encoded by χ) in a suspension directly from experimental measurements of the sedimentation velocity under varied particle concentration. The approach of estimating χ via ω has the advantage that ω increases monotonously with χ and hence can in principle be inverted to obtain a mapping χ(ω). The downside of this approach is that the nonlinear fits tend to be more sensitive to noise in v(φ ) than linear fits.
We note in conclusion that our results, while obtained under the assumption of a strongly repulsive DLVO potential at ζ = 50 mV, are in fact generally valid for any repulsive potential with a steep potential barrier at distance ξ = χŝ φ . The van der Waals interactions are strongly subdued in most of our parameter regime and the models we used to predict K at both small and large χ are not specific to details of the DLVO interaction.
In future work we aim to reproduce long-ranged electrostatic interactions in sedimentation velocity experiments for a broad parameter range and compare the experimental data directly to our simulations. In the experimental setup we wish to study model nanoparticle systems including a controlled degree of polydispersity. Other possible avenues of future research might include nonspherical, in particular rod-like charged particles, where orientation and rotation become important in addition to translational ordering.

Conflicts of interest
There are no conflicts to declare.