Mesoscale simulations of diffusion and sedimentation in shape-anisotropic nanoparticle suspensions

We determine the long-time self-diffusion coefficient and sedimentation coefficient for suspensions of nanoparticles with anisotropic shapes (octahedra, cubes, tetrahedra, and spherocylinders) as a function of nanoparticle concentration using mesoscale simulations. We use a discrete particle model for the nanoparticles, and we account for solvent-mediated hydrodynamic interactions between nanoparticles using the multiparticle collision dynamics method. Our simulations are compared to theoretical predictions and experimental data from existing literature, demonstrating good agreement in the majority of cases. Further, we find that the self-diffusion coefficient of the regular polyhedral shapes can be estimated from that of a sphere whose diameter is average of their inscribed and circumscribed sphere diameters.


I. INTRODUCTION
4][5] Many factors impact the motion of NPs, including their size, interactions with each other, and interactions with their surroundings. 6This work focuses specifically on the effect of shape, which has emerged as an important factor for modulating the properties and function of NPs in many practical applications. 7,8For example, shape-anisotropic iron-oxide-based magnetic NPs were shown to enhance contrast for magnetic resonance imaging compared to spherical NPs, 8 while quantum rods were shown to have enhanced diffusion compared to quantum dots in confined networks. 9Given that NPs with a variety of shapes can be readily synthesized [10][11][12] and that many naturally occurring NPs (e.g., the rod-like tobacco mosaic virus 13,14 and gibbsite platelets 15 ) also exhibit pronounced shape anisotropy, it is important to develop a fundamental understanding of the relationship between shape and transport properties, such as diffusion coefficients, in order to engineer NPs for practical applications.
Experimentally characterizing how NP dynamics depend on shape and concentration can be challenging.For example, dynamic light scattering is a common technique for measuring NP diffusion from fluctuations in scattered a) These authors contributed equally.b) Electronic mail: anikouba@ipfdd.dec) Electronic mail: mphoward@auburn.edulight intensity. 16However, knowledge about the distribution of NP sizes and/or shapes is needed to extract the diffusion coefficient from the raw measurement data, 17,18 and it is difficult to perform this analysis for non-dilute solutions. 19Camera-based tracking of tagged NPs is an alternative approach that allows for the direct measurement of the NP diffusion coefficient, 9,20 but this method has limited spatial and temporal resolution. 20NP properties may also be affected if labeling agents, such as fluorescent markers, are used. 21Further, it can be difficult to prepare NP suspensions with sufficiently low polydispersity and at high enough concentrations to accurately assess how the diffusion coefficient varies with both NP characteristics and concentration.
As a result, theory and simulations have proven to be useful approaches for studying the dynamics of NP suspensions.Early theories predominantly focused on spherical NPs, for which the single-particle translational and rotational diffusion coefficients can be calculated using the classical Stokes-Einstein and Stokes-Einstein-Debye relations, respectively.Theoretical predictions for the first-order concentration dependence of the long-time self-diffusion coefficient for suspensions of spherical NPs have also been derived. 22,235][26][27][28][29] At finite concentration, the diffusive motion of the rods becomes more complex but can be split qualitatively into three regimes: at dilute concentrations, rods have essentially unrestricted motion in all directions; at semi-dilute concentrations, their motion is slightly hindered perpendicular to the long axis of the rod; and at high concentrations, the perpendicular diffusive motion is entirely suppressed. 30However, predicting the dynamics of rodlike NPs with quantitative accuracy still remains challenging because their anisotropic shape can lead to com-arXiv:2403.03009v1[cond-mat.soft]5 Mar 2024 plex flow patterns around individual NPs and to nontrivial collective behavior such as nematic or smectic ordering.For more complicated NP shapes than rods, predicting even single-particle diffusion coefficients becomes challenging, and numerical approaches are often required. 31,32In general, fully analytic descriptions of NP dynamics in suspensions are challenging to construct due to the many-body hydrodynamic interactions (HIs) between NPs that are mediated by the solvent.
Computer simulations are highly useful tools for numerically investigating NP dynamics in suspensions.The main challenge is to construct models that capture the relevant physics while remaining computationally tractable.Explicitly resolving both the NPs and the solvent molecules they are suspended in using, e.g., classical molecular dynamics (MD) approaches, quickly becomes infeasible because NPs are typically much larger than solvent molecules.However, given the corresponding separation of time scales between the solvent dynamics and NP dynamics, it is often possible to overcome this difficulty using coarse-grained models having simplified or implicit treatments of the solvent. 33For example, Brownian dynamics (BD) is a well-known implicit-solvent technique that accounts for solvent-induced drag and fluctuating forces on the NPs, 34 but which neglects HIs between the NPs in its most basic form.HIs can be introduced to BD through appropriate mobility tensors, 35 such as the pairwise far-field Rotne-Prager-Yamakawa tensor for spherical particles. 36,37Stokesian dynamics, a gold-standard approach for simulating colloidal suspensions, additionally accounts for short-range lubrication forces between NPs within the BD framework. 38,39However, BD approaches that include HIs are often still computationally demanding to implement and require expressions for the mobility tensor, which may be difficult to obtain for complex NP shapes.
To circumvent issues determining inputs needed for a fully implicit treatment of the solvent, several mesoscale simulation methods, including multiparticle collision dynamics (MPCD), 33,40,41 dissipative particle dynamics, 42,43 and the lattice Boltzmann method, 44,45 use simplified particle-based solvent models that are less demanding to simulate than an atomistic model but still have properties resembling that of real solvents.In this work, we will use MPCD because we have recently shown that MPCD can reasonably reproduce expected results for the long-time self-diffusion coefficient and sedimentation coefficient for suspensions of spherical NPs over a range of NP concentrations, 46 and the same approach used to model the spherical NPs can be extended to NPs with other shapes.][48][49][50][51][52] The current state-of-the-art coupling scheme, first proposed by Poblete et al., uses a discrete particle model that represents an NP as a mesh of "vertex" particles interconnected via elastic springs. 47The solvent particles interact with the NPs only through stochastic collisions that are straightforward to compute.4][55][56][57][58] However, we are unaware of a systematic study using MPCD to characterize the long-time self-diffusion coefficients and sedimentation coefficients for suspensions of shape-anisotropic NPs at varying concentrations.
In this work, we use MPCD with a discrete particle model to study the long-time self-diffusion and sedimentation coefficients of octahedra, cubes, tetrahedra, and spherocylinders as a function of NP concentration.We investigate the effect of shape by comparing the results for the different NP shapes with each other and with spheres.We also assess the influence of solvent-mediated HIs by comparing the MPCD simulations with implicitsolvent Langevin dynamics simulations.

II. MODELS A. Multiparticle collision dynamics
In MPCD, the solvent consists of point particles that are propagated in alternating streaming and collision steps that occur at a regular time interval ∆t.During the streaming step, the solvent particles move according to Newton's equations of motion, where r i is the position, v i is the velocity, and m i is the mass of particle i.All solvent particles have the same mass m.Unlike standard MD particles, MPCD particles do not interact with each other by pairwise forces, but each particle may be acted on by a body force F i .For a constant F i , eq. ( 1) can be integrated analytically to give the standard equations of ballistic motion.
In the collision step, the solvent particles are sorted into cubic cells of edge length ℓ, then exchange momentum with particles in the same cell according to a collision scheme.Here, we use the stochastic rotation dynamics (SRD) scheme without angular momentum conservation. 40,41SRD updates the velocity of particle i in cell j according to: where u j is the mass-averaged velocity of the particles in cell j and Ω j is the rotation matrix for cell j.The matrix Ω j rotates about an axis randomly selected for cell j by a fixed angle α.At each collision step, the collision cells are shifted along each Cartesian direction by a random amount drawn uniformly from [−ℓ/2, +ℓ/2] to ensure Galilean invariance, 59,60 and a cell-level Maxwellian thermostat is used to maintain a constant temperature T . 61he natural units for MPCD simulations are the length ℓ of the collision cells, the mass m of the solvent particles, and the thermal energy k B T , where k B is the Boltzmann constant.The corresponding unit of time is τ = mℓ 2 β, where β = 1/(k B T ).We adopted the standard SRD parameters ∆t = 0.1 τ , α = 130 • , and average solvent number density 5 ℓ −3 , which give a liquid-like Newtonian fluid with dynamic viscosity η 0 = 3.95 k B T τ /ℓ 3 . 62

B. Discrete particle model
A discrete particle model was used to represent the NPs and couple them to the solvent. 46,47The NP shapes we modeled were a sphere, an octahedron, a cube, a tetrahedron, and two spherocylinders with different aspect ratios (Fig. 1).Each NP consisted of N v vertex particles on the surface of the shape, and each vertex particle had mass 5 m.The vertex particles were bonded to their nearest neighbors with a harmonic potential, where r is the distance between two particles, r b is the distance required for the bond by the shape, and k b is the spring constant.To ensure that the NPs maintained their shapes, the vertex particles were also bonded to either an additional particle in the center of the NP (sphere, octahedron, cube, & tetrahedron) or their diametrically opposed vertex particle (spherocylinders).Excludedvolume interactions between NPs were modeled by applying the Weeks-Chandler-Andersen repulsive potential 63 between vertex particles All vertex particles (but not the central particle) were coupled to the MPCD solvent through the collision step eq.( 2). 47Between collision steps, the central and vertex particles moved according to eq. ( 1).Based on our prior work, 46 we used k b = 5000 ℓ −2 to make stiff bonds and σ = ℓ, and we integrated eq. ( 1) using the velocity Verlet algorithm with time step 0.005 τ .We visually confirmed that all NPs maintained a nearly rigid shape and that no NPs penetrated each other for the vertex-particle configurations chosen as described next.
Sphere.-We modeled a sphere having diameter d = 6 ℓ [Fig.1(a)] as a reference point.To create the vertex particles, we subdivided the triangular faces of a regular icosahedron twice and scaled the positions of all vertices to lie on the surface of the sphere.This process resulted in N v = 162 vertex particles with a typical nearest-neighbor distance between 0.83 ℓ and 0.97 ℓ.Note that this model differs from the one we used in ref. 46 in two ways: (1) the number of vertex particles is larger and (2) the excluded volume is handled through the vertex particles rather than through the central particle.These choices were made in this work so that the spheres would have a comparable surface density of vertex particles and the same style of excluded-volume interactions as the anisotropic NPs we studied.
Octahedron and tetrahedron.-Wemodeled a regular octahedron and a regular tetrahedron both having edge length a = 6 ℓ [Figs.1(b)  and 1(d)].Because the faces of these polyhedra are equilateral triangles, we first created a three-dimensional triangulated model of each shape using computer-aided design software, then subdivided the faces 3 times to create a triangular mesh of vertex particles.This process resulted in 9 vertex particles per edge and distance a/8 = 0.75 ℓ between all nearest-neighbor vertex particles for both shapes.The total number of vertex particles was N v = 258 for the octahedron and N v = 130 for the tetrahedron.
Cube.-We modeled a cube with edge length a = 6 ℓ [Fig.1(c)] using a square mesh with 8 vertex particles per edge.The total number of vertex particles was N v = 296, and the distance between nearest-neighbor vertex particles was a/7 ≈ 0.86 ℓ.This is the same vertex particle configuration as in ref. 46, and the description here corrects a typographical error for the number of vertex particles per edge.Unlike in ref. 46, though, a central particle was used to maintain rigidity to keep consistency with the sphere, octahedron, and tetrahedron.
Spherocylinder.-We modeled two types of spherocylinders: both had two hemispheres with diameter d = 6 ℓ, but one had a cylinder of length h = 6 ℓ while the other cylinder had the length h = 12 ℓ [Fig.1(e)].Thus, the spherocylinders had aspect ratios λ ≡ d/h = 1 and 2, respectively.This degree of anisotropy is much smaller than that of rods typically used in experiments, which usually are in the range of λ = 10 to 150. 65However, we found surprisingly little data on transport coefficients for rod-like NPs with these smaller aspect ratios and, as noted previously, theories are also hard to formulate in this regime. 66Hence, we chose to study spherocylinders with these smaller aspect ratios to begin to bridge the knowledge gap between spheres and long rod-like NPs.Discrete particle models for the spherocylinders were constructed through a multi-step process: First, a mesh of vertex particles for the hemispheres was created by slicing our discrete sphere model in half along a plane that exposed 20 evenly spaced vertex particles around its circumference and had 91 vertex particles in total.Then, vertex particles for the cylinder were generated from the ring of 20 exposed vertex particles by translating the ring by 0.75 ℓ and rotating it around the axis of the cylinder by 9 • to stagger the particles on consecutive rings.This process was repeated until the entire cylinder surface was covered with vertex particles.The total number of vertex particles per spherocylinder was N v = 322 for λ = 1 and N v = 482 for λ = 2, with the nearest-neighbor distance between vertex particles ranging from 0.83 ℓ to 0.97 ℓ.

C. Simulation details
We performed bulk simulations containing N NPs in a cubic simulation box with edge length L = 120 ℓ and periodic boundary conditions.We simulated a range of nominal NP volume fractions ϕ = N v/L 3 , where v is the nominal volume of each NP (Table I), by varying N .We created equilibrated configurations of NPs at the different volume fractions using Langevin dynamics (LD) simulations.LD simulations are faster to perform than MPCD simulations because they do not include HI, and we also chose the friction coefficient for the LD simulations to give faster NP dynamics than in the MPCD simulations in order to accelerate equilibration.Starting from these configurations, we measured the long-time self-diffusion coefficient as a function of ϕ using equilibrium simulations (Section III A) and the sedimentation coefficient as a function of ϕ using nonequilibrium simulations (Section III B).All simulations were conducted using HOOMDblue 67,68 (version 2.9.7) extended with azplugins 69 (version 0.12.0).
For the spheres and regular polyhedra, we performed one equilibrium simulation of length 2 × 10 5 τ and recorded the position of all central particles every 10 τ .We performed one nonequilibrium simulation consisting of a warmup period of 0.5 × 10 5 τ to reach steady state followed by a production period of length 1.5 × 10 5 τ in which we recorded the average velocity of the NPs every 0.105 τ and the average velocity of the solvent every 0.1 τ .To estimate error bars, we subdivided these trajectories into three blocks and computed the standard error between blocks.
For the spherocylinders, we performed eight equilibrium simulations of length 10 5 τ and recorded the position of enough vertex particles to reconstruct the center of mass of each NP every 2.5 τ .We performed three nonequilibrium simulations consisting of a 0.5 × 10 5 τ warmup period and 1 × 10 5 τ production period with the velocities sampled in the same way as for the other shapes.Error bars were estimated as the standard error of the independent simulations.
To improve statistics, we averaged ⟨∆r 2 ⟩ over NPs and time origins, and we extracted D from the time average of the long-time plateau of d⟨∆r 2 ⟩/ dt, which we fit in the time range 10 4 τ ≤ t < 2 × 10 4 τ for the spheres and regular polyhedra and in the range 3 × 10 4 τ ≤ t < 5 × 10 4 τ for the spherocylinders.Note that in defining D in this way, the long-time self-diffusion coefficient is a scalar quantity.For anisotropic NPs, the short-time motion is characterized by a diffusion tensor; this tensor is isotropic for the regular polyhedra we have studied, 70 but it is anisotropic for the spherocylinders. 71Hence, D reported in this work implies an orientational average at long times for the spherocylinders.3][74] For a cubic simulation box such as ours, the self-diffusion coefficient in an infinitely large box D ∞ is related to D in a finite box with edge length L by 73,74 where ξ ≈ 2.837297 and η is the suspension viscosity.Applying eq. ( 6) can be challenging in practice because it requires knowledge of η, which depends on the shape and volume fraction of the NPs.Analytic expressions for η exist for some NP shapes, 75,76 but they are typically only valid for small NP volume fractions; 77 hence, additional costly simulations are usually needed to accurately determine η.To avoid this step, we approximated η with a Stokes-Einstein-like proportionality,  46,74 where D 0 = k B T /γ 0 is the longtime self-diffusion coefficient at infinite dilution (i.e., the single-particle limit) and γ 0 is the corresponding hydrodynamic friction coefficient for the NP (again, orientationally averaged for the spherocylinders).Substituting for η in eq. ( 6) and solving for D ∞ yields We previously tested this approximate correction by computing D for spherical NPs in different box sizes L and confirming that D ∞ was independent of L within our measurement accuracy. 46o apply eq. ( 7), γ 0 must be determined for each NP shape.Experimental correlations 78 for γ 0 exist [e.g., eq. ( 9) below]; however, it is not guaranteed that the MPCD simulations are consistent with these.Instead, we noted that all diffusivities are corrected by the same factor in eq. ( 7) regardless of ϕ and that eq. ( 6) can be used directly when ϕ is sufficiently small that η ≈ η 0 .Accordingly, we linearly extrapolated our measured D to ϕ = 0, using the data from the smallest two values of ϕ that we simulated, to obtain a measured D 0 with finite-size effects.We then applied eq. ( 6) with η = η 0 to calculate D ∞ 0 from D 0 and used the ratio D ∞ 0 /D 0 as the finite-size correction factor for all D. In the rest of the paper, all diffusion coefficients have been corrected for finite-size effects in this way, but we will still refer to them as D and D 0 for notational simplicity.

Regular polyhedra
We first investigated the shape-dependence of the longtime self-diffusion coefficient extrapolated to infinite dilution D 0 for the regular polyhedra we simulated (octahedron, cube, and tetrahedron).Pettyjohn and Christiansen experimentally measured the settling rates of particles with these shapes at low Reynolds number. 70,78hey found that the settling rate could be correlated with particle shape using the sphericity ψ, defined as the ratio of the surface area of a sphere having the same volume as the shape to the actual surface area A of the shape, The sphericities of our regular polyhedra are listed in Table I.Using the correlation for the settling velocity from ref. 78, a correlation for the hydrodynamic friction coefficient γ 0 is Note that the first term in parentheses is the diameter of an equivalent-volume sphere to the shape, so eq.( 9) gives γ 0 = 3πη 0 d for a sphere with diameter d as expected.
Based on eq. ( 9), a cube should diffuse more slowly than an octahedron, while an octahedron should diffuse more slowly than a tetrahedron when all have the same edge length a; a sphere with diameter d = a is predicted to have D 0 between that of the octahedron and the tetrahedron (Table II).Indeed, our simulation results for D 0 were qualitatively consistent with these predictions.Quantitatively, D 0 from the cube simulations was in excellent agreement with the value predicted using eq.( 9), but D 0 from the octahedron and tetrahedron simulations was 9% and 18% smaller, respectively.We calculated a similar deviation between the measured and predicted D 0 for tetrahedra in recent experiments by Hoffmann and coworkers, who fabricated tetrahedral clusters from four spherical polystyrene NPs with diameter 154 nm; they measured a self-diffusion coefficient of D 0 = 1.72 × 10 −12 m 2 /s in water, 79,80 which is 22% smaller than the predicted value of D 0 = 2.2×10 −12 m 2 /s when using an edge length of a = 308 nm in eq. ( 9).These clusters are, however, not true tetrahedra so it is unclear whether this deviation from eq. ( 9) should be expected in the MPCD simulations too.
We and others previously found that D 0 for a cube can also be reasonably well-approximated by D 0 for a sphere with diameter d = (d I + d C )/2, the arithmetic mean of the diameters d I and d C of the spheres that inscribe and circumscribe it, respectively. 32,46We carried out the same calculation for the octahedron and tetrahedron, and we again found good agreement with our TABLE II.Diffusion coefficient at infinite dilution D0 for the sphere and regular polyhedra calculated from our simulations, using eq.( 9), and using the Stokes-Einstein relationship for a sphere with mean diameter d given in Table I.All are in units of 10 −3 ℓ 2 /τ .simulation using eq.( 9  simulations (Table II).Thus, using d seems to provide a quick and reasonable estimate of D 0 for regular polyhedra as an alternative to eq. ( 9).We next investigated the volume-fraction dependence of D [Fig.2(a)].Given that the different NP shapes had different D 0 , we report D/D 0 to facilitate comparison between shapes [Fig.2(b)].The tetrahedra exhibited the strongest dependence on ϕ, the spheres exhibited the weakest dependence on ϕ, while both the cubes and octahedra exhibited a similar dependence on ϕ that was intermediate between the spheres and tetrahedra.In general, we expected D to decrease when ϕ increased because increased interactions between NPs usually slow their motion.At low NP volume fractions, long-ranged solventmediated HIs are important because short-ranged interactions are infrequent.Differences in the dependence of D/D 0 on ϕ seen in Fig. 2 when ϕ is small are then likely caused by differences in HIs between shapes.
At higher NP volume fractions, direct interactions between NPs become more frequent and significant, partic- for the transition to pentagonal dipyramids at 0.55aex and 0.75aex. 81,82arly those due to excluded-volume between NPs.Indeed, we expect that eventually D/D 0 → 0 when the NPs reach a freezing or jamming transition that essentially traps each NP in a local cage of surrounding NPs.2][83][84] However, we noted that the actual excluded volume v ex of the NPs (and hence excludedvolume fraction ϕ ex ) differs from the nominal volume v (and nominal volume fraction ϕ) because the vertex particles in our discrete model have a finite size.For example, the vertex particles on the surface of the cube [Fig.1(c)] protrude roughly σ/2, so the edge length of the volume excluded by the cube is roughly σ longer than the nominal edge length.In general, we can estimate the excluded-volume edge length a ex of our regular polyhedra as a ex = a(1+σ/d I ), such that the excluded-volume polyhedron contains all spheres with diameter σ on the surface of the nominal polyhedron.The ratio of the excluded volume to nominal volume is then v ex /v ≈ (1 + σ/d I ) 3 and ϕ ex is proportionally larger than ϕ by the same factor.This larger excluded size a ex is evident in the radial distribution function g(r) (Fig. 3) for all shapes.We attempted to assess the effect of this difference in nominal and excluded volume using the spherical NPs.We performed additional simulations where we implemented the excluded-volume interaction between spheres through a core-shifted Weeks-Chandler-Andersen potential between only their central particles, like in ref. 46.
As expected, we found that there was less structuring in the fluid, measured through g(r), at a given nominal volume fraction ϕ due to the smaller excluded volume of each sphere [Fig.4(a)].However, we found little difference in the diffusivity over the range of volume fractions investigated [Fig.4(b)].Moreover, the simulation data of D/D 0 agreed well with experimental data when plotted using ϕ.We observed similar agreement between MPCD simulations and experiments for cubes using the nominal volume fraction ϕ in our prior work 46 [see also Fig. 5(c)].Hence, at least over the range of volume fractions considered for the spheres, the nominal volume fraction ϕ seems to be a good description of the concentration.We note, though, that structural effects caused by differences in nominal and excluded volume may still become significant at sufficiently high excluded volume fractions, particularly if a phase transition is approached.The tetrahedron, which has the largest v ex /v ratio of our regular polyhedra, is an excellent example of this point.Previous simulations of hard tetrahedra 81,82 revealed a transition to a fluid consisting of pentagonal dipyramids when the volume fraction was 0.47.That study found that g(r) showed a distinct signature of this transition: at low volume fractions where dipyramids did not form, g(r) had its first peak at r = 0.75 a; whereas, at higher volume fractions where dipyramids formed, this original peak disappeared, and the first peak shifted to a much smaller distance r = 0.55 a.Our largest nominal volume fraction ϕ = 0.20 is well below the reported transition to dipyramids, but if we instead consider the excluded volume fraction (ϕ ex = 0.56), then the system should have surpassed this transition.When we computed g(r) for the tetrahedra [Fig.3(d)], we observed these signature peaks emerging at the reported r if a ex was used rather than a.Thus, the tetrahedra appear to undergo a transition to dipyramids that is not expected using only ϕ.The more dramatic slowing down of the tetrahedra dynamics with ϕ compared to the other shapes could be partially due to this transition.
Finally, we assessed the influence of HIs between the NPs on their long-time self-diffusion by performing complementary LD simulations that do not have these interactions (Fig. 5).Due to the lack of long-ranged solventmediated HIs, no finite-size corrections are needed.Qualitatively, D/D 0 had a similar dependence on ϕ both with and without HIs, with differences for the tetrahedra being most pronounced and differences for the cubes being least pronounced.However, there were clear quantitative differences between the MPCD simulations with HIs and the LD simulations without HIs.For all shapes, D/D 0 was smaller for a given ϕ (had a stronger ϕ dependence) in the LD simulations compared to the MPCD simulations.Taken together, these differences support the established picture that solvent-mediated HIs and excluded-volume interactions between NPs that determine their fluid structure both play a role in determining the NP dynamics.
As an aside, we remark that the agreement between our MPCD simulations and experiments 85,86 significantly improved for the spherical NPs compared to our previous study, 46 which is likely due to the higher surface density of vertex particles used in this work.The accuracy of discrete particle models typically improves with increasing surface density, 88 and the surface density of vertex particles on the sphere was roughly four times that of ref. 46.We note that Poblete et al. recommended an optimal surface density of 0.53/ℓ 2 for spheres in MPCD to balance discretization and inertia effects, 47 which lies between the value of 0.37/ℓ 2 used in ref. 46 and 1.43/ℓ 2 used here.The surface density of vertex particles used for the other regular polyhedra was comparable to that of the spheres.

Spherocylinders
Having studied the self-diffusion of these regular polyhedra, we next investigated the long-time self-diffusion of spherocylinders.Bolhuis and Frenkel numerically studied the phase diagram of hard spherocylinders for a range of aspect ratios λ. 89 For λ ≤ 2, the spherocylinders exhibited only two phases-a low-density isotropic phase and a high-density crystal phase-with the transition between these occurring at volume fraction 0.58 and 0.53 for λ = 1 and 2, respectively.We therefore restricted our simulations to ϕ ≤ 0.30, which corresponds to ϕ ex < 0.44 for our spherocylinders (v ex /v = 1.45 and 1.41 for λ = 1 and 2), in order to focus our calculations on the isotropic phase.We confirmed this was the case by computing a global nematic order parameter, 90,91 finding it to be close to zero (0.02 and 0.03 for λ = 1 and 2 when ϕ = 0.30) as expected for an isotropic phase.
In the isotropic phase, the translational diffusion of rod-like objects is the orientational average of their parallel and normal components.The self-diffusion coefficient of rods in the infinite dilution limit can be estimated as 92 where the last three terms in the parenthesis correct for end effects. 28,29,93This equation gives D 0 = 2.94 × 10 −3 ℓ 2 /τ and 2.41 × 10 −3 ℓ 2 /τ for rods with λ = 1 and 2 in our MPCD solvent, respectively.Our simulated values D 0 = 3.09 × 10 −3 ℓ 2 /τ and 2.57 × 10 −3 ℓ 2 /τ were within 5% of eq. ( 10), showing the expected decrease of D 0 with λ.We also note that eq. ( 10) underpredicts the diffusivity of a sphere (λ = 0) by about 5% compared to the classical Stokes-Einstein relation.
The concentration dependence of D/D 0 with ϕ was similar for both spherocylinders (Fig. 6).Indeed, D/D 0 for the spherocylinders with λ = 1 was nearly indistinguishable from that for the spheres (λ = 0).The longer spherocylinders with λ = 2 showed some systematic differences, consistently having a slightly smaller value than for λ = 1 at a given ϕ.This result indicates that even a small amount of anisotropy may begin to have an effect on the diffusive dynamics, but the magnitude of this effect seems to be small.We also compared our simulation data to the parametric fit of ref. 94.Our simulations qualitatively agreed with the prediction that D/D 0 should be smaller for a larger λ at a given ϕ, but the simulations consistently had smaller values of D/D 0 than predicted.We note that ref. 94 used BD simulations that lacked HIs to develop this fit, so it is unclear to what extent we should expect agreement to simulations with HIs.

B. Sedimentation
After investigating the long-time self-diffusion coefficients of our shape-anisotropic NPs, we characterized their sedimentation coefficients.This complementary dynamic property of a suspension is important for understanding, e.g., how NPs settle under gravity.We defined the sedimentation coefficient K from the linear propor- tionality between the average velocity u of an NP under a sufficiently small applied force F, To measure K in our simulations, we imposed a constant force F = f x x on all NPs, where x is the unit vector in the x direction, and measured their average velocity u x = u • x.The applied forces were f x = 0.5 k B T /ℓ and 1.0 k B T /ℓ per NP, which we distributed evenly among all the vertex and central particles in each NP.A balancing force was applied to the MPCD particles to ensure that the total force on the system was zero.We allowed the system to reach a steady state under the imposed forces, performed a production run where we measured u x , and extracted K from a linear regression of u x and f x .As for diffusion coefficients, the sedimentation coefficients from our MPCD simulations must be corrected for finite-size effects from periodic boundary conditions.The sedimentation coefficient measured in an infinitely large box K ∞ is related to the one measured in a finite box by 95 where S(0) is the static structure factor at zero wavevector.This structure factor is related to the isothermal compressibility and so can be computed from an equation of state.Here, we used the virial expansion of the pressure, which gives where Bn = B n /v n−1 and B n is the n-th virial coefficient.We used ϕ ex in eq. ( 13) because it should characterize the structure of the suspension better than ϕ (see discussion of Fig. 3).We used up to the 8th virial coefficient for the regular polyhedra 96 and up to the 5th virial coefficient for the spherocylinders. 97,98Like eq. ( 6), eq. ( 12) also includes the suspension viscosity η so we made the same Stokes-Einstein-like approximation to eliminate this dependency, We used the finite-size-corrected D/D 0 and computed γ 0 = k B T /D 0 from the finite-size-corrected D 0 .Note that eqs.( 12) and ( 14) fix an error in eqs.( 19) and ( 20) of ref. 46.All sedimentation coefficients are corrected in this way, but for notational simplicity, we will refer to them as K in the remaining discussion.MPCD conserves linear momentum, so the sedimentation coefficients calculated directly from the simulation are in a frame of reference where the mass-averaged velocity of the NPs and solvent is zero.However, it is a common practice to consider suspensions in the frame of reference where the volume-averaged velocity is zero, i.e., ⟨u⟩ = ϕu + (1 − ϕ)u 0 = 0 where u 0 is the solvent velocity.Shifting from the mass-averaged to volume-averaged frame of reference amounts to a rescaling of K, which we implemented as in our previous work. 46All values of K are presented in the frame of reference where the volume-averaged velocity is zero.
The sedimentation coefficients of the regular polyhedra [Fig.7(a)] exhibited a qualitatively similar dependence on shape and concentration as the self-diffusion coefficients did.We consistently found that the spheres had the largest K, the tetrahedra had the smallest K, while the octahedra and cubes had an intermediate K.Moreover, all sedimentation coefficients decreased with increasing concentration, as expected, with the tetrahedra having the strongest concentration dependence.The sedimentation coefficients of both spherocylinders were highly similar to each other and to that of the sphere [Fig.7(b)].These behaviors are qualitatively similar to the self-diffusion coefficients, so we will not repeat that discussion here for brevity.

IV. CONCLUSIONS
We investigated the long-time self-diffusion and sedimentation of NPs with anisotropic shapes.
The anisotropic shapes we studied were an octahedron, a cube, a tetrahedron, and a spherocylinder.The NPs were represented with a discrete particle model and were hydrodynamically coupled to each other using the multiparticle collision dynamics method.Simulations were conducted across a range of volume fractions for each shape where the NPs remained in a fluid/isotropic phase.Our modeling approach can be easily extended to explore the dynamics of other NP shapes, e.g., irregular polyhedra and non-convex shapes. 99or regular polyhedra equal edge lengths, shape had a clear influence on their transport properties.Octahedra and cubes were slower diffusing than spheres with diameter equal to their edge length for all investigated volume fractions [Fig.2(a)].Tetrahedra diffused the fastest at small volume fractions but the slowest at larger volume fractions, which we partially attributed to the formation of pentagonal dipyramids.The simulated self-diffusion coefficients of all investigated NP shapes at infinite dilution were in good agreement with a correlation based on sphericity and also with an approximation using the mean diameter of the spheres that inscribed and circumscribed the shapes.After accounting for differences due to shape at infinite dilution [Fig.2(b)], the self-diffusion coefficient of the spheres showed the weakest volume-fraction dependence, that of the tetrahedra showed the strongest volume-fraction dependence, while the octahedra and cubes showed intermediate behavior.Similar trends were found for the dependence of the sedimentation coefficients on volume fraction [Fig.7(a)].
For small-aspect-ratio spherocylinders (λ = 1 and 2), the diffusion coefficients at infinite dilution showed a dependence on aspect ratio that was consistent with theoretical expectation, meaning that the spherocylinders diffused more slowly as aspect ratio increased [Fig.6(a)].However, after accounting for shape effects at infinite dilution, the self-diffusion coefficient [Fig.6(b)] had a volume-fraction dependence that closely followed that of spheres having diameter equal to the spherocylinders, with only minor differences for the spherocylinder with λ = 2.The sedimentation coefficient [Fig.7(b)] had essentially the same volume-fraction dependence for the spheres and both spherocylinders.We expect that the dynamics of spherocylinders should deviate more significantly from spheres as λ increases, and in principle, we can expand our spherocylinder model to study this regime.However, doing so incurs higher computational cost due to a substantial increase in the number of vertex particles per spherocylinder.Further, we would need larger simulation boxes to accommodate these spherocylinders and gather good statistics, thereby also increasing the number of solvent particles required.To mitigate these computational challenges, an alternative approach is to represent the spherocylinders as linear rods comprised of partially overlapping particles. 100However, establishing a connection between this model, our spherocylinder model, and experiments is still an open question, which we plan to explore.470113688, and 509039598.This work was completed with resources provided by the Auburn University Easley Cluster, the supercomputer Mogon at Johannes Gutenberg University Mainz (www.hpc.uni-mainz.de), and the Texas Advanced Computing Center (TACC) at The University of Texas at Austin.

FIG. 1 .
FIG. 1. Discrete particle model for (a) sphere, (b) octahedron, (c) cube, (d) tetrahedron and (e) two spherocylinders (aspect ratios λ = 1 and 2).The Nv vertex particles are shown in orange, and the bonds to their nearest neighbors are shown in blue.To improve the clarity of these renderings, the size of the vertex particles has been decreased, and central particles and additional bonds used to maintain the shape have been omitted.These snapshots were rendered using Visual Molecular Dynamics (version 1.9.3)64 III. RESULTS AND DISCUSSION A. Long-time self-diffusion coefficient We computed the long-time self-diffusion coefficient D of the NPs from the time derivative of the average mean squared displacement ⟨∆r 2 ⟩ of each NP, 34 D = lim t→∞ 1 6 d⟨∆r 2 ⟩ dt .

FIG. 2 .
FIG. 2. (a) Long-time self-diffusion coefficient D of spheres, octahedra, cubes, and tetrahedra as a function of nominal volume fraction ϕ.(b) D normalized by its value linearly extrapolated to infinite dilution D0.

86 FIG. 4 .
FIG. 4. (a)Radial distribution functions for spheres with excluded volume handled through either vertex particles (solid lines) or a central particle (dashed lines).(b) Long-time selfdiffusion coefficient D (normalized by D0) for the same systems compared with experimental data.85,86

FIG. 5 .
FIG. 5. Comparison of long-time self-diffusion coefficient D (normalized by D0) of (a) spheres, (b) octahedra, (c) cubes, and (d) tetrahedra as functions of volume fraction ϕ from MPCD and LD simulations.Experimental data is included in (a) and (c) from multiple sources.The experimental values of D for the spheres85,86 were scaled by the Stokes-Einstein prediction for D0, while the experimental values of D for the cubes87 were scaled such that D/D0 = 1 for the lowest-concentration point in that data set (ϕ ≈ 0).

FIG. 6 .
FIG. 6.(a) Long-time self-diffusion coefficient D for spherocylinders with aspect ratios λ = 1 and 2 as functions of volume fraction ϕ.The sphere data from Fig. 2 is included as a reference point with λ = 0. (b) D normalized by its value extrapolated to infinite dilution D0.The dashed curves in (b) are D/D0 predicted from the fit of ref. 94.

FIG. 7 .
FIG.7.Sedimentation coefficient K of (a) spheres, octahedra, cubes, and tetrahedra and (b) spherocylinders as a function of volume fraction ϕ.The frame of reference used to define K is the one where the volume-averaged velocity is zero.

TABLE I .
Geometric properties of the regular polyhedra investigated.General formulae are given in terms of the edge length a, with the specific value for a = 6 ℓ (the edge length for all our polyhedral NPs) quoted in parentheses.The properties are the volume v, surface area A, sphericity ψ, inscribed-sphere diameter dI, circumscribed-sphere diameter dC, and mean of inscribed-sphere and circumscribed-sphere diameters d.