A generalised drift-correcting time integration scheme for Brownian suspensions of rigid particles with arbitrary shape

The efficient computation of the overdamped, random motion of micron and nanometre scale particles in a viscous fluid requires novel methods to obtain the hydrodynamic interactions, random displacements and Brownian drift at minimal cost. Capturing Brownian drift is done most efficiently through a judiciously constructed time-integration scheme that automatically accounts for its contribution to particle motion. In this paper, we present a generalised drift-correcting (gDC) scheme that accounts for Brownian drift for suspensions of rigid particles with arbitrary shape. The scheme seamlessly integrates with fast methods for computing the hydrodynamic interactions and random increments and requires a single full mobility solve per time-step. As a result, the gDC provides increased computational efficiency when used in conjunction with grid-based methods that employ fluctuating hydrodynamics to obtain the random increments. Further, for these methods the additional computations that the scheme requires occur at the level of individual particles, and hence lend themselves naturally to parallel computation. We perform a series of simulations that demonstrate the gDC obtains similar levels of accuracy as compared with the existing state-of-the-art. In addition, these simulations illustrate the gDC's applicability to a wide array of relevant problems involving Brownian suspensions of non-spherical particles, such as the structure of liquid crystals and the rheology of complex fluids.


Introduction
Brownian motion, the random motion exhibited by colloidal particles due to thermal fluctuations in the surrounding viscous fluid, is a common and important feature of fluidic systems at the micron and nanometre scales [1,2,3].It plays a fundamental role in determining the distribution and configuration of polymers and particles that comprise complex fluids [4,5], which then in turn affect suspension rheology and bulk mechanical response to applied stresses [6,7,8].Further, Brownian motion and thermal fluctuations play a critical role in biological and cellular transport processes, often in conjunction with the cell's active mechanisms that utilise chemical energy [9].
To simulate these systems, one often considers the equations of motion in the Brownian dynamics [10], or overdamped [11], limit where momentum variables are taken to have reached thermal equilibrium and the velocity field in the surrounding fluid is governed by the steady Stokes equations.In the context of rigid particles, this means that the dynamics of the particle positions and orientations are provided by the overdamped Langevin equations, and accordingly, the dynamics of the corresponding probability distribution of the particle positions and orientations is described by Smoluchowski's equation.
The hydrodynamic interactions between the particles due to the surrounding fluid are embodied in the dense, position-dependent mobility matrix [12,13] that relates the generalised forces on the particles to the resulting particle velocities.The hydrodynamic interactions between the particles also impact the fluctuations that they experience.The fluctuation-dissipation theorem [14] stipulates that the covariance of the random particle velocities is proportional to the mobility matrix and hence so too is the particle diffusion matrix.This requires the random particle velocities to be proportional to the square root of the mobility matrix.Finally, due to hydrodynamic interactions, the overdamped limit produces a nontrivial thermal, or Brownian, drift term [10] in the equations of motion that is proportional to the divergence of the mobility matrix with respect to the particle positions and orientations.
Initially, the computations of these three quantities -the mobility matrix, random velocities, and Brownian drift -in approaches such as Stokesian Dynamics [15] utilised a variety of direct methods.The entries of the mobility matrix were determined pairwise from analytical expressions derived in relevant asymptotic limits and the random particle velocities were computed using a Cholesky decomposition of the resulting matrix.Finally, Brownian drift was incorporated directly by evaluating expressions for the derivatives of the mobility matrix entries.A reliance on these direct approaches limited system sizes, as well as particle shapes to spheres.
Since these first computations, there have been several key advances, some of which have been incorporated into Stokesian dynamics [16,17], including fast, matrix-free methods for determining the particle velocities given the forces applied to them.These methods provide the action of the mobility matrix on a vector in O(N log N ) operations, where N is the vector length, without ever computing the matrix itself.Such approaches are built around fast summation In this work, we generalise the DC scheme to simulate rigid bodies of arbitrary shape.Like the original DC, the generalised DC (gDC) scheme requires a single full mobility solve per time step.This significantly accelerates time integration for schemes built around grid-based solvers that can take advantage of the flows generated by a fluctuating stress.The gDC may also be used for matrix-based computations at a cost comparable to the existing state-of-the-art.The main idea behind the scheme is to advance to the mid-step using particle velocities obtained by orthogonally projecting the fluctuating velocities onto the space of rigid body motions.At the midstep, the full mobility is treated and important factors based on the divergence of the projected fluctuating velocity are incorporated into the update.The gDC has the nice property that many of the additional computations needed to capture Brownian drift occur at the level of individual particles, and hence naturally lend themselves to parallel computation and largerscale simulation.We show that the resulting scheme is weakly first-order accurate and, in practice, provides errors similar in magnitude to the T-S scheme.Finally, we demonstrate the applicability of the scheme for larger-scale simulation by considering confined suspensions of rod-like particles, as well as the rheology of Czech hedgehog particle suspensions.Consider a suspension of N rigid particles where the centre of mass position of particle p at time t is Y p (t), while the rotation relative to its initial orientation is given by the unit quaternion q p (t) (see fig. 1).We provide a detailed overview of using quaternions to represent rotations in Appendix A, but present some of the key facts here.Any vector b in the body frame of particle p is given by B(t) = R(q p (t))b in the lab frame at time t.The rotation matrix R is related to entries of the quaternion through

The overdamped Langevin equation for rigid body motion
Suppose that the particles are subject to an external potential U (Y 1 , u 1 , . . ., Y N , u N ).The force on particle p is then −∂ Yp U , while the corresponding torque is −D − up ∂ up U .A derivation of this expression for the torque is presented in Appendix B. The matrix D −1 up is the 3 × 3 'dexpinv' matrix defined by where The 'dexpinv' matrix relates the angular velocity of particle p, Ω p , to the time derivative of u p through In the overdamped limit, the force and torque vectors are linearly related to their contributions to the particle velocity and angular velocity vectors via the 6N × 6N mobility matrix, N, whose entries depend on the particle configuration.The application of the mobility matrix is equivalent to solving the Stokes boundary value problem for a collection of rigid particles subject to applied forces and torques.Thus, in the absence of Brownian motion, the equations of motion are given by the system of differential equations where x = Y 1 , u 1 , . . ., Y N , u N describes the positions and orientations of all particles and N = ΦNΦ with When Brownian motion is present, the particles will also move randomly due to thermal fluctuations in the surrounding fluid.In the overdamped limit, the effects of Brownian motion are captured through the inclusion of random increments of the positions and Lie algebra elements, turning the equations of motion into a system of stochastic differential equations.Specifically, these increments are given by √ 2k B T N 1/2 dW , where k B T is the thermal energy, N 1/2 = ΦN 1/2 with N 1/2 being the matrix square root of the mobility matrix, and W is a 6N × 1 vector of independent Wiener processes.The dependence on the mobility matrix ensures that the fluctuation-dissipation theorem is satisfied, a necessary condition for the Boltzmann distribution to be recovered at equilibrium.
Finally, along with the random velocities, the overdamped equations of motion require the inclusion of the thermal drift term, k B T ∂ x • Ndt [10].This ensures that the stochastic differential equation with an Itô interpretation of the stochastic integral yields dynamics consistent with those of Smoluchowski's equation for the corresponding probability distribution.Putting these terms together, we arrive at the equations of motion for the particle positions and orientations,

Mobility matrix
The purpose of this paper is to develop time integration schemes for eq.( 8) that avoid the direct computation of the thermal drift term and naturally interface with matrix-free methods for computing the random increments and particle velocities arising from the applied forces.In doing so, we will assume that the particles are discretised into M total degrees of freedom, as depicted in fig. 2. This discretisation can involve the surface of the particles, as is done when considering the first-kind boundary integral representation of the Stokes equations [46,12], or can be related to a volume discretisation where elements of the particle volume are represented by regularised distributions of force, for example.Let r i denote the position of discrete degree of freedom i, v i ≡ dr i /dt be its velocity and λ i be the force it exerts on the fluid.The hydrodynamic interactions between the discrete degrees of freedom provide a linear relationship between v i and λ i such that where M ij is the matrix that relates the force on discrete degree of freedom j to the velocity of discrete degree of freedom i.Additionally, the requirement that all discrete degrees of freedom belonging to particle p move as a single rigid body (see fig. 2) stipulates that where Ẏp is the translational velocity of particle p, Ω p is its angular velocity and B p is the set of discrete degrees of freedom belonging to particle p. 9) and ( 10) for all particles p = 1, . . ., N can be equated to yield where is the 3M × 3M mobility matrix for the discrete degrees of freedom and K is a 3M × 6N block matrix that maps the velocity and angular velocity from the space of rigid body motions to the velocities of the discrete degrees of freedom.Specifically is the 3 × 6 block that relates the rigid body motion of particle p to the translational velocity of discrete degree of freedom i.Note that the specific form of M will depend on how the discretisation is performed, however, here we will assume that the discretisation results in a positive definite M, thereby preserving the structure of the underlying continuous equations.This is an important ingredient for the fluctuation-dissipation theorem to be satisfied exactly in the discrete setting.
In the absence of inertia, the conservation of linear and angular momentum reduce to the balance of force and torque, which can be written compactly for the N particles as where F is the total force and torque vector applied on the particles.The mobility problem given by eqs.( 11) and ( 14) forms a saddle point system After eliminating λ, one obtains and therefore the mobility matrix is given by As we describe later, our discretisation will follow from the rigid blob framework [47], but the gDC time integration scheme is applicable to any spatial discretisation provided that the resulting mobility matrix N can be related to the saddle point system in eq. ( 15).This particular decomposition has the advantage that the computation of the stochastic particle velocities can be readily incorporated into the system by including a random velocity for the discrete degrees of freedom [44].Specifically, we consider the system with the random velocity, where M 1/2 is the square root of M which satisfies M 1/2 (M 1/2 ) = M, and W (t) is a 3M × 1 vector of Wiener processes.After eliminating λ, one obtains where V = NK M −1 v.As a result the fluctuating body velocities are given by V where N 1/2 satisfies the fluctuation dissipation balance with N 1/2 (N 1/2 ) = N.
As shown by Balboa-Usabiaga et al. [48], the symmetric saddle point problem (18) can be efficiently solved with an iterative Krylov subspace method, such as GMRES, using a simple block diagonal preconditioner P .The matrix P is constructed by setting to zero the entries of M in (18) if they correspond to interactions between discrete degrees of freedom from different rigid bodies.As K is a block diagonal matrix, the cost of solving the linear system (18) will be related primarily to the costs of applying the matrix M and computing v.

The generalised drifter-corrector (gDC)
The primary result of this paper is the following generalisation of the DC scheme to rigid bodies of arbitrary shape.The gDC updates the particle positions and Lie algebra elements while automatically accounting for the Brownian drift term that arises in the equations of motion.It has the advantage of requiring only a single full mobility solve per time step.For gridbased solvers, the gDC incurs nearly the same cost as applying the well-known Euler-Maruyama scheme, outlined in Appendix D, which does not capture the Brownian drift.For matrix-based approaches the cost will be similar to the existing state-of-the-art [44,45].
Below we present the gDC algorithm.The time-step is indicated by superscript, n, e.g.V n would correspond to V (x n ) = V (x(t n )) in the continuous setting.For the purposes of implementing the algorithm, it should be noted that the collection of Lie algebra elements contained in x n is set to zero at the start of each time step (see Appendix A).
2. Solve the following saddle point problem at the start of the time-step, yielding Note that this is also the solution of the least-squares problem min For a small parameter δ, the divergence can be calculated via1 (a) finite-differencing: where is the i-th component of the velocity of particle p, and e (p) i is a 6N × 1 vector whose nonzero entries, corresponding to the position and orientation of particle p, are given by the i-th Cartesian basis vector of R 6 ; or (b) random finite-differences: where W is a vector of independent N (0, 1) random variables which is also independent of W .By replacing the divergence in the definition of ν with the term inside the expectation in eq. ( 25), the appropriate value of ν will be achieved in expectation, which is ultimately all that is required to demonstrate weak-accuracy.
We find that option (a) is best partnered with grid-based mobility methods, whilst option (b) is best suited to matrix-based methods to reduce the computational cost.We discuss this in more detail in section 5.
4. Move the body positions and orientations to the mid-step according to 5. Solve the full mobility problem at the mid-step, for where F m contains the forces and torques acting on the bodies at the mid-step.Note that quantities are evaluated at the mid-step using the mid-step quaternions q m p = exp u m p • q n p .
6. Update the positions and orientations of the bodies according to with the new orientation quaternion of each particle p defined by q n+1 p = exp u n+1 p • q n p .
In words, the gDC first computes a set of particle velocities by solving a least squares problem to find the appropriate stochastic velocity, V, in the space of rigid body motions.This velocity is used to advance the particle positions to the mid-step where the full computation is performed.The positions are then updated using the velocity computed at the mid-step scaled by a factor based on the divergence of V.As shown in Appendix E, the scheme recovers the first and second moments of the increment ∆x = x n+1 − x n to first-order in expectation.

Applying M and computing v
As shown in section 2.2, hydrodynamic interactions between bodies are directly obtained from the matrix M. The action of M on a vector is required to compute the deterministic velocities V = NF , while the product M 1/2 W is necessary to obtain the Brownian velocities V .It is therefore essential to use and develop efficient methods to compute the action of M and its square root.
In this section, we outline two distinct strategies with similar resolution of hydrodynamic interactions between the discrete degrees of freedom.The first follows directly from an implementation of the rigid multiblob model [47] where the rigid particles are discretised into surface or volume elements that interact via the Rotne-Prager-Yamakawa (RPY) tensor.The other approach called the fluctuating Force Coupling Method (FCM) is matrix-free and simultaneously applies M and computes V by solving the fluctuating Stokes equations on a grid with a forcing term that accounts for the presence of the particles.

RPY tensor
The well-known RPY mobility matrix was originally developed to provide a symmetric positive definite, pairwise approximation of the mobility matrix for a collection of spherical particles of equal radii in an unbounded domain [49].Extensions of the RPY matrix for particles of different radii [50], in a background shear flow [49], and above a no-slip boundary [51] are available in the literature.Following the rigid multiblob model, the RPY tensor can be used to provide the hydrodynamic interactions between the discrete degrees of freedom making up the rigid particles.When discretising a particle surface, the rigid multiblob model provides a first-order accurate approximation to the particle mobility.
In general, the computational cost of these matrix-vector products scales quadratically with the number of discrete degrees of freedom.More sophisticated methods, such as the fast multipole methods (FMM) [52] and Ewald methods, achieve a linear scaling.For periodic domains, one can use the positively split Ewald method [22,23] to compute the action of M on a vector.In our computations below, we perform a direct pairwise evaluation of the wall-corrected RPY tensor [51].
The action of M 1/2 on the random vector W is obtained through the Lanczos algorithm [31].To achieve convergence, the Krylov subspace K is enriched iteratively with basis vectors that are linear combinations of the powers of the mobility matrix times the random vector: K = span W , MW , M 2 W , ... .The cost of the method therefore depends on the number of basis vectors, and thus mobility-vector products, required to reach a given tolerance .

Fluctuating Force Coupling Method
The other approach relies on a matrix-free method called fluctuating Force Coupling Method (FCM) [53], which combines FCM [54] with fluctuating hydrodynamics [32].With fluctuating FCM, the coupling between the discrete degrees of freedom and the fluid is achieved through a forcing term added to the fluctuating Stokes equations, where η is the dynamic fluid viscosity, u the fluid velocity and p is the pressure.The first term in the RHS of eq. ( 30) is the divergence of the fluctuating stress tensor, W . W is delta correlated in time and space with the following statistics, The second term on the RHS of eq. ( 30) is the forcing transferred to the fluid with a spreading operator S, where the finite size of the discrete degrees of freedom is accounted for with a Gaussian spreading envelope, The length scale σ functions in a similar way to the bead radius, a, in the RPY tensor.The diagonal entries of the FCM mobility matrices will match those of the RPY tensor if σ = a/ √ π [54].
The fluid velocity, obtained after solving eqs.( 30) and (31), is the sum of a deterministic part, due to the forcing λ, and a fluctuating term, stemming from the fluctuating stress.We may express the total fluid velocity as where S is the spreading operator in eq. ( 34), L −1 is the inverse Stokes operator (i.e. the fluid solver), and D the divergence operator applied to the fluctuating stress [36].
The velocities of the discrete degrees of freedom are then obtained from the fluid velocity using an averaging operator, J, such that where J = S is adjoint to the spreading operator.We see then that the FCM mobility matrix can be written as the composition of three linear operators: M F CM = JL −1 S. Additionally, as demonstrated in [39] the velocity v satisfies the fluctuationdissipation theorem with the covariance given by the FCM approximation of the mobility matrix, F CM = JL −1 D. Thus, we see that fluctuating FCM simultaneously computes the actions of M and M 1/2 on λ and W , respectively, by considering solutions to the forced fluctuating Stokes equations.
Following [39,40], we implement fluctuating FCM in periodic domains.This involves first evaluating the forcing, f (x), and the fluctuating stress, W , on a regular grid.Next, the Stokes equations are solved using a Fourier spectral method.Finally, the trapezoidal rule is used to numerically integrate eq. ( 38) to obtain the translational velocity for each of the discrete degrees of freedom making up the rigid bodies.

Computational cost of the gDC
Due to the differences in their implementations, there is a distinct difference in cost between RPY and grid-based fluctuating FCM that manifests itself when computing the random velocity vector, v.For both implementations, the most expensive computation is applying the matrix M. With RPY, as we are directly handling M, any time the particle positions change, we must compute v = M 1/2 W using the Lanczos algorithm, for which each iteration requires a matrix-vector multiplication involving M.
For fluctuating FCM, however, we are instead working with a decomposition of M In this decomposition, the only matrix that depends on the particle configuration is J. Thus, provided that W remains the same, changing the particle positions simply requires reaveraging the fluid velocity ȗ at the new positions, and avoids having to recompute the velocity field itself.Cost of Step 2: Along with the computation of v, Step 2 of the gDC algorithm in section 5 involves solving a simplified mobility problem with M replaced by I.This is equivalent to solving the least squares problem min V KV − v , the solution of which is given by eq. ( 23) and involves the 6N × 6N symmetric, block-diagonal matrix K K.The 6 × 6 block associated with body p is given by where B p is the set of discrete degrees of freedom belonging to body p.Therefore, computing K K −1 just requires inverting N small 6 × 6 blocks using dense linear algebra, which is negligible in terms of computational cost and, additionally, naturally lends itself to parallel computation.
Thus, the primary cost associated with Step 2 involves computing v, which for fluctuating FCM incurs the cost of one Stokes solve and M averaging operations.The RPY implementation will require M Lanczos mobility-vector products with the Lanczos method.Cost of Step 3: The main cost of this step is computing the scalar term ν using finite differences (FD) or RFD to generate ∂ x • V. Finite differencing involves perturbing the positions and orientations of each rigid body and calculating the new fluctuating velocities for the discrete degrees of freedom at the disturbed states.With RFD, the new fluctuating velocities are only calculated once, after disturbing the particles' positions and orientations with a random vector W .
For fluctuating FCM, FD is straight forward and incurs minimal cost.Since the body velocities are uncorrelated in eq. ( 22), the finite-differencing step in eq. ( 24) requires 6M averaging operations.As discussed above, however, changing the particle position with the RPY-based implementation will require recomputing M 1/2 W with each displacement.Thus, we require 6N × M Lanczos mobility-vector products.This is a cost that could rise quickly as the particle number increases.On the other hand, using RFD in eq. ( 25) only requires M Lanczos mobilityvector products with matrix-based methods and M averaging operations with fluctuating FCM.Cost of Step 5: This step involves the costliest computation of solving the full mobility problem, eq. ( 18), along with the additional computation of v.
For the fluctuating FCM implementation, each GMRES iteration requires a full fluctuating FCM computation that involves spreading the force vector λ, solving the Stokes equations, and finally averaging the resulting flow field to obtain the velocities for the discrete degrees of freedom.A single additional averaging is needed to incorporate the fluctuating velocities, vm , on the RHS of eq. ( 27).
For matrix-based methods, each GMRES iteration will incur the cost of a matrix-vector product involving M.There will also be the cost of M Lanczos Lanczos interations needed to obtain vm .
For a tolerance of = 10 −3 , both the Lanczos and GMRES with block preconditioning require approximately M Lanczos = M GMRES = 5 interations on average.The tolerance of = 10 −3 has been shown [44] to be a suitable choice for simulations involving fluctuations.Adding these costs together assuming = 10 −3 , we see that the cost for one gDC-FD time iteration with fluctuating FCM is (1 + 5 =) 6 Stokes solves, 5M spreading operations and (1 + 6 + 1 + 5 =) 13M averaging operations.For the RPY based implementation, we have (5 + 6N × 5 + 5 + 5 =) 5 × (6N + 3) matrix-vector products involving M. With gDC-RFD the number of averaging operations reduces to 8M for fluctuating FCM, and the total number of matrix-vector products involving M with the RPY based implementation lowers to 20, becoming independent of the particle number N .
As shown in section 6.2, the gDC is more stable with FD than with RFD.For grid-based methods, the 5M additional averaging operations incurred by FD are negligible compared to the Stokes solve.This is not true for matrix based methods for which RFD represents a significant cost reduction.
Based on these costs, we see that the gDC with FD is well-suited for grid-based methods utilising fluctuating hydrodynamics such as fluctuating FCM [53,55], PSE [22,23] or the fluctuating and stochastic Immersed Boundaries methods [56,57,36,45], while gDC with RFD is the best compromise between numerical stability and computational cost for matrix-based methods.

Simulations
In this section, we demonstrate the performance of the gDC by performing simulations of particulate suspensions under both dynamic and equilibrium conditions.In Section 6.2, we investigate the accuracy of the gDC by examining the equilibrium distributions of the position and orientation of a boomerang-shaped particle in a gravitational field above a planar no-slip boundary.Owing to the small size of this system, we can integrate for long times and acquire many realisations to quantify the temporal accuracy of the gDC.We compare these results with those obtained using the state-of-the-art integrators developed in Sprinkle et al. [44].After this, we demonstrate the suitability of the gDC for large scale simulations of anisotropic Brownian particles.In Section 6.3, we study the equilibrium properties of confined suspensions of rodshaped particles, a model for a liquid crystal or stiff-polymer system.Finally, inspired by recent experimental work [58], in Section 6.4 we use the gDC to perform simulations to obtain the rheology of colloidal suspensions made from Czech hedgehog-shaped (CH) particles.

Implementations
The simulations that follow rely on two implementations: one using direct evaluations of the wall-corrected RPY tensor and the other using fluctuating FCM.
• Python with RPY: for the single boomerang simulations we use the collaborative code, called "RigidMultiblobWall", that one of us co-developed [48] and used for large scale simulations of colloidal active particles [59,47].The code contains most of the time integration schemes recently developed for matrix-based approaches [60,47,44].It relies on the pairwise evaluation of the wall-corrected RPY tensor.The mobility computations can be accelerated in different ways: using C++ routines, Numba or GPUs via the PyCUDA interface.We implemented the gDC in the code to compare its performance with the most accurate scheme developed by Sprinkle et al.: the Trapezoidal Slip scheme (T-S).
• C++ with fluctuating FCM: for the confined liquid crystal and rheology simulations, we use an FFT-based fluid solver with fluctuating FCM in a periodic box.FFTs are parallelised with the scalable MPI library FFTW.MPI is also used to parallelise the spreading and averaging operations as well as the nearest neighbour search necessary to compute short-ranged, pairwise interactions between discrete degrees of freedom.At each grid point, the entries of the fluctuating stress are independent Gaussian random variables.Since the fluctuating stress W is symmetric, at each grid point we generate six random numbers based on a Gaussian distribution with zero mean and unit variance.We then multiply the off-diagonal entries by 2k B T η/(∆x 3 )∆t −1/2 and the diagonal ones by 2 k B T η/(∆x 3 )∆t −1/2 .We truncate the Gaussian envelopes by setting ∆ i (x) = 0 for x − r i > 3a, and the length scale of the envelope, σ, is related to the grid size, ∆x, through σ/∆x = 1.86.For more details, we refer the reader to our previous work [40].

Convergence study: boomerang above no-slip boundary
In this section we study the accuracy of the gDC scheme by examining the equilibrium distributions of a boomerang particle in a gravitational field above a no-slip boundary.Micron-size boomerangs have been used extensively to study the diffusion of colloidal anisotropic particles in experiments [61] and simulations [60,48].Following Delong et al. [60] and as shown in fig.3, we discretise the boomerang using 15 RPY-particles separated by distance a, the radius of the RPY-particle.The boomerang is subject to a gravitational force of magnitude mg = 0.18k B T /a and each RPY-particle interacts with the planar boundary through a short-ranged, repulsive potential, where h is the height of the RPY-particle's centre from the surface, V 0 = 23k B T sets the potential strength, and b = 0.5a is the range of the potential.Simulations are performed using both the T-S and gDC schemes for three different time-step sizes ∆t = 0.1, 0.2, 0.3τ D , where τ D = a 2 /D = 6πηa 3 /k B T is the typical diffusive timescale associated with an RPY-particle.The divergence term, ν, in Step 3 of the gDC algorithm (see section 3) is computed using both FD (24) and RFD (25).The solver tolerance for the GMRES and Lanczos algorithm is set to 10 −4 , and the finite difference parameter is δ = 10 −6 .
In order to obtain sufficient statistics, we run 36 different simulations for each value of ∆t.Each simulation is initialized using a Monte Carlo generated sample of the Gibbs-Boltzmann distribution and is run to a final time of 30000τ D , where the solution is recorded at every t = 0.3τ D .
Figure 4 compares the height distributions from the T-S and gDC simulations with the true distribution obtained using Markov Chain Monte Carlo (MCMC).We see that, as expected, both schemes converge to the MCMC distribution as ∆t decreases.The L 2 errors for the distributions shown in the right panel indicate that the gDC achieves a similar accuracy and convergence rate as the T-S.However, the RFD computation of the ν term in the gDC algorithm leads to larger errors than FD.This is particularly true for the largest time-step, ∆t = 0.3τ D , for which gDC-RFD exhibit numerical stability issues due to particle overlaps across the wall.Similar stability problems were also encountered with the T-S scheme, resulting in several simulations being discarded and restarted when ∆t = 0.3τ D .We suspect that the loss of stability with gDC-RFD is due to overestimations of (∂ x • V) n caused by the random sampling W , which leads to large values of the corrective term ν, therefore potentially moving the particle too far into the wall at the next time step.
The orientation distribution, shown in fig.5, displays a similar trend, though the gDC exhibits a smaller error which does not change considerably with the time-step size.
p y i p 1 j + A P n 8 w f j Z J E L < / l a t e x i t > p < l a t e x i t s h a 1 _ b a s e 6 4 = " 7 6 F 0 P e 8 K i S 0 n The angle θ is that between the vertical axis and the unit vector, p, that is orthogonal to both boomerang arms.L2-error between the simulated distributions and the true distribution generated using MCMC.

Confined liquid crystal
Inspired by the alignment of slender molecules in liquid crystals, this section examines the equilibrium distributions of confined suspensions of rod-shaped particles.In this set of simulations, each rod is constructed from 22 particles of radius a as illustrated in fig.6.The simulations are performed using fluctuating FCM in periodic computational domains with dimensions [0, Confinement is introduced by first including slip boundaries at y = 0 and y = L y through boundary conditions u • e y = 0 and I − e y e y ∇u = 0 on the flow field.These conditions are enforced using an image system [40] for both the forces on the bodies, as well as the fluctuating stress, in the y > L y half of the computational domain.Second, the rods are kept away from the channel walls by including a repulsive harmonic potential on each FCM-particle making up the rod.Here, r is the distance between an FCM-particle centre and the wall, b = 2.2a is the cut-off distance, and K = 40k B T ab .Along with this potential, overlap between FCM-particles is also discouraged using the soft potential with r now denoting the centre-to-centre distance between 2 FCM-particles, U 0 = 20k B T and b = 0.5a.In the following simulations, the domain length is set to L x = L z = L = 77.65a.We consider three different channel widths, namely L y = L/8, L y = L/4 and L y = L/2, to explore how confinement impacts the rod distributions.The number of rods in these simulations is taken to be 64, 128 and 256 respectively to maintain a constant volume fraction of 10%.A representative simulation snapshot shown in fig. 7.Each simulation was run until a final time of 50τ D /3 and the time step size was set to ∆t = τ D /600.Finally, the initial conditions for the gDC simulations were generated using samples from independent MCMC runs.
To measure the distribution of the rods within the channel, we compute the distribution of the y-coordinates of geometric centre of the rods.Additionally, we compute the distribution of the angle θ between the rod-axis and the xz-plane, such that θ = 0 corresponds to the rod being parallel to the channel boundaries.Figure 8 compares the y-coordinate distributions from the gDC simulations with the Gibbs-Boltzmann distributions calculated using MCMC.For all channel widths that were considered, good agreement is observed between these distributions.Similarly close agreement is found for the θ-distributions in fig.9; the shape of the distribution is well-preserved even when viewed on a logarithmic scale.All of the y-coordinate distributions exhibit the expected symmetry about the midpoint of the channel at y = L y /2, and have peaks in probability density near the slip boundaries.The position of the peaks depends on the width of the slip channel because of the finite thickness of the rods.It is also clear that the narrower the channel, the higher the probability of the y-positions associated with these peaks.The θdistributions all decrease from a maximum value at θ ≈ 0, although for the narrowest channel, the distribution exhibits a second local maximum just below θ = π/8.The wider the channel, the more slowly the probability density decays, and hence the larger the range of likely orientations.
For the narrowest channel, the channel width, L y = L/8 ≈ 9.7a, is slightly larger than twice the diameter of the rod cross-section (2 × 4a = 8a), but less than the rod length of 16.14a.As a result, the harmonic potential which keeps rods away from the slip surfaces will also induce a torque on the rods even for small deviations from θ = 0. Thus, the rods are most likely to be found with an orientation θ ≈ 0 in one of the two layers displaced from the channel centre.The distance between these layers balances the repulsive potentials between different rods and between particles and the walls.The smaller peak at the channel centre represents the small number of rods that get 'jammed' between the two layers.Additionally, as there is insufficient space for all rods to have θ ≈ 0, they also tend to span the channel diagonally with their ends trapped between the channel walls, leading to the secondary peak in the θ-distribution below θ = π/8.
As the channel width increases when we have L y = L/4 ≈ 19.4a, the central peak in the y-coordinate distribution becomes more pronounced, as shown in fig.8. Additionally, we observe a wider range of possible θ-value for rods near the channel centre, although we note that a small, secondary peak in the θ distribution remains as seen in the log-scale plot in fig.9. Note that this peak occurs at a larger θ value than for the narrowest channel due to the increased channel width.Rods in the outer layers between this central layer and the slip boundaries are still most likely to have θ ≈ 0. As a consequence, the rods in the middle layer are also likely to align with θ = 0 with the two outer rod layers acting like the channel walls.Therefore, we see that the angle distribution retains its maximum at low θ.
For the largest width, the θ distribution still exhibits its peak value at θ ≈ 0 as the rods nearest the slip boundaries continue to align parallel to the boundaries.We again observe peaks in the y-coordinate distribution close to the channel boundaries, though their magnitude is now reduced as shown in fig.8.The central peak observed in the two previous cases, however, has now spread out to form a series of peaks of decreasing magnitude in the interior of the channel.Note that this persistent pattern of density peaks which are largest near the walls and which decrease towards the channel centre is consistent with existing observations on the layering of colloidal particles near repulsive boundaries [62,63,64,65,66].Rods appear to spend comparable amounts of time between and in these internal layers, seemingly making layers away from the slip surfaces more transient than those near the boundaries where rods tend to spend significant amounts of time.Finally, rods in the channel interior achieve a broader distribution of θ values with no discernible secondary peak in the θ distribution.

Rheology of suspensions of Czech hedgehog colloids
In this section, we use fluctuating FCM with the gDC to perform nonequilibrium simulations that examine the rheology of Brownian Czech hedgehog-shaped (CH) particles, depicted in fig.10.This set of simulations is inspired by the recent experimental work of Bourrianne et al. [67], where suspensions of dendritic, silica particles, either hydrophobic or hydrophilic, were found to exhibit interesting rheological behaviour, including discontinuous shear thickening, depending on the particle interactions, as well as the relative strengths of shear and particle diffusion (Peclet number) and volume fraction.
In our simulations, the CH particles are formed of 19 FCM-particles of radius a, as shown in fig.10.This shape was chosen to reproduce, at some level, the complex structure and high specific surface area of the silica particles in the experiments.We consider two types of interactions between the CH particles.For the first type, the FCM-particles comprising the CH particles repel at short-range through a soft potential, while attracting at longer range.Specifically, defining for A = 5k B T and λ = 0.5a, the 'repel-attract' potential is given by where U (r) is the soft-sphere potential defined in eq. ( 43).For the second type of interaction, we have the 'repel-attract-repel' potential    where a barrier is introduced in the potential just beyond the well.The resulting shapes of these potentials are shown in fig.10.Note that strictly these definitions are only r ≥ 2a, with the potentials continuously extended by linear functions for r < 2a to ensure that the force is constant for blobs which overlap, just as for the original soft-sphere potential U (r).
The simulations are performed in a triply-periodic computational domain with dimensions L x = L y = L z = 77.65a.In order ensure that a strong rheological response is observed, a high volume fraction of 20% was used for all simulations.This corresponds to 1177 CH particles, yielding a 1177×19 = 22363 FCM-particles in each simulation.Simulations were run with a time step length of ∆t = τ D /600.The initial conditions for the full simulations were constructed by running simulations that ignore thermal fluctuations and hydrodynamic interactions, i.e. using a diagonal mobility matrix, but retain interactions due to the potentials.In practice, we find that by allowing the CH particles to settle in this way before the shear is applied, the 'long-time' velocity profile is realised more quickly.The simulations were run to final times between 35τ D /3 and 50τ D and in all cases, a regular oscillatory velocity profile emerges long before the end of the simulation.
To measure the suspension viscosity from this emergent velocity field, we adopt the approach described by Vázquez et al. [68].Here, the background fluid is forced using the periodic force density f (x, y, z) = f 0 sin (ky) e x (47) where f 0 is the force magnitude and k = 2π/L y is the wavenumber.The x-component of the resulting velocity field of the suspension is averaged over time and in the x-and z-directions, resulting in a sinusoidal velocity profile in y with magnitude v x .The effective viscosity of the suspension can then be defined as In fig.11, the ratio of the effective viscosity η ef f to the underlying fluid viscosity η is compared to the shear rate γ for suspensions with both types of particle interactions.In this work, γ is taken as the slope of the sinusoidal velocity profile.We find that both the 'repel-attract' and 'repel-attract-repel' exhibit shear at lower shear rates, in the original experiments of Bourrianne et al.We do not observe, shear thickening for the 'repel-attract-repel' CH particle suspensions whose interactions are meant to be similar to those of the hydrophilic particles in the experiments.Bourrianne et al. suggest that solid friction is required for the onset of shear thickening and that hydrogen bonding is additionally required for DST to be observed.Both of these interactions are absent from our simulations.The bottom right plot is displaying the same data as the bottom left plot but over a smaller r range.The colour of the curves fades from black to white as f0 decreases.
Given that this is a non-equilibrium problem, it is of particular interest to examine distributions of the particle states since they cannot be sampled using equilibrium procedures such as MCMC. Figure 11 shows the radial distribution function g(r) where r is the centre-to-centre distance of r for the CH particles with 'repel-attract' interactions; the distribution for 'repelattract-repel' interactions is similar.Simulations with higher values of f 0 correspond to darker curves in the figure.At the smallest f 0 , where the particle motion is dominated by the Brownian motion, g(r) increases from zero to a peak value at separations between r/a ≈ 8 and r/a ≈ 10.At larger separations, g(r) gradually decreases to approach values for a uniform distribution.As f 0 , and hence the shear rate, is increased, this peak diminishes and the distribution approaches a monotonically increasing function.

Summary and conclusions
In this paper, we developed the gDC scheme for the efficient time integration of hydrodynamically interacting rigid Brownian particles.The gDC automatically accounts for Brownian drift in advancing particle positions and orientations while retaining the need to perform only a single, full mobility problem at each time step, while achieving accuracy comparable to the current state-of-the-art.The gDC has been designed to be used with fast methods for applying the mobility matrix and generating random increments with the correct covariance and is ideally suited for grid-based computations that take advantage of fluctuating hydrodynamics.Additionally, many of the additional computations that the gDC requires are local to each particle, enabling parallelisation and large-scale computation.Our simulations demonstrate each of these features, as well as the applicability of gDC to enable simulations that complement experiments in modern suspension mechanics and colloidal science.eq.(C.1) and hence that the rotational dynamics in our approach are equivalent to those in the work of Delong et al. [60].Given our result in Appendix B, namely that the torque on the body can be expressed using the Lie algebra element as T = −D − u ∂ u U , we propose that the appropriate Itô SDE is the direct analogue of eq.(C.1).In the remainder of this section, we demonstrate that this does indeed result in the correct Itô SDE for the unit quaternion.Loosely speaking, the approach here is to convert the Itô SDE for u to the Stratonovich interpretation, for which many standard calculus results are recovered, then to mirror the deterministic result "du/dt = D −1 u Ω =⇒ dq/dt = Ψ q Ω" and finally to recover the Itô interpretation.To that end, recall that if y satisfies the Itô SDE bk /∂u j .For u = 0, we can simply differentiate eq. ( 4) and observe that ∂D −1 jb /∂u j ∝ [u×] 2 u b = 0.The case u = 0 can be shown easily using the limit definition of the partial derivative.Thus for the first moment in expectation.
To simplify notation, for the remainder of this section all quantities are evaluated at time n unless explicitly indicated otherwise.We start by expanding V m about time t n to obtain Recalling that ν = 1 + ∆t 2 ∂ j V j we obtain (E.7) We recognize the second and fourth terms on the right hand side as a product rule expansion and re-write the above as (E.8) All but the second term vanish in expectation, leaving us with (E.9) Note that if we are calculating the divergence using RFD instead, eqs.(E.7) and (E.8)only hold in expectation and are obtained by exploiting the independence of W and W to separate terms.Nevertheless, we obtain eq.(E.9).Now, recalling the definitions of V and V we have Using W n W s = δ ns it is easy to verify that Substituting this into eq.(E.9) and multiplying through by ∆t reveals that eq.(E.5) is satisfied to first-order.To show the same for eq.(E.3), recall eq.(E.8) and observe that (E.12) The second term vanishes in expectation and we have ∆x i ∆x a = ∆t 2 V i V a + O(∆t 2 ).(E.13) If we are using RFD to calculate the divergence, then the expansion in eq.(E.12)only holds in expectation but we obtain eq.(E.13) nonetheless.Again recalling the definition of V from section 3 we find that In a similar way to before, substituting W n W e = δ ne and summing over repeated indices yields giving us the desired result.

Figure 1 :
Figure 1: Illustration of a rigid body with position Yp and orientation qp.

Figure 2 :
Figure 2: Discretisation of the rigid body surface with discrete degrees of freedom.

✓
< l a t e x i t s h a 1 _ b a s e 6 4 = "J q E n Y v V 6 P t s K B J Y m B V w E p j I M A N w = " > A A A B 7 X i c b V D L S g N B E J y N r x h f U Y9 e B o P g K e y K o M e g F 4 8 R z A O S J c x O O s m Y 2 Z l l p l c I S / 7 B i w d F v P o / 3 v w b J 8 k e N L G g o a j q p r s r S q S w 6 9 5 e V 2 k 0 e R 5 G c k F N y T g J y R W r k j t R J g 3 D y S J 7 J K 3 n z t P f i v X s f i 9 a C l 8 8 c k z / w P n 8 A o / e P K A = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " J q E n Y v V 6 P t s K B J Y m B V w E p j I M A N w = " > A A A B 7 X i c b V D L S g N B E J y N r x h f U Y 9 e B o P g K e y K o M e g F 4 8 R z A O S J c x O O s m Y 2 Z l l p l c I S / 7 B i w d F v P o / 3 v w b J 8 k e N L G g o a j q p r s r S q S w 6 9 5 e V 2 k 0 e R 5 G c k F N y T g J y R W r k j t R J g 3 D y S J 7 J K 3 n z t P f i v X s f i 9 a C l 8 8 c k z / w P n 8 A o / e P K A = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " J q E n Y v V 6 P t s K B J Y m B V w E p j I M A N w = " > A A A B 7 X i c b V D L S g N B E J y N r x h f U Y 9 e B o P g K e y K o M e g F 4 8 R z A O S J c x O O s m Y 2 Z l l p l c I S / 7 B i w d F v P o / 3 v w b J 8 k e N L G g o a j q p r s r S q S w 6 9 5 e V 2 k 0 e R 5 G c k F N y T g J y R W r k j t R J g 3 D y S J 7 J K 3 n z t P f i v X s f i 9 a C l 8 8 c k z / w P n 8 A o / e P K A = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " J q E n Y v V 6 P t s K B J Y m B V w E p j I M A N w = " > A A A B 7 X i c b V D L S g N B E J y N r x h f U Y 9 e B o P g K e y K o M e g F 4 8 R z A O S J c x O O s m Y 2 Z l l p l c I S / 7 B i w d F v P o / 3 v w b J 8 k e N L G g o a j q p r s r S q S w 6 9 5 e V 2 k 0 e R 5 G c k F N y T g J y R W r k j t R J g 3 D y S J 7 J K 3 n z t P f i v X s f i 9 a C l 8 8 c k z / w P n 8 A o / e P K A = = < / l a t e x i t > h < l a t e x i t s h a 1 _ b a s e 6 4 = " S s YK O O t 1 j T f i N k y H K 7 x 6 G F F E b x M = " > A A A B 6 H i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l E q M e i F4 8 t 2 A 9 o Q 9 l s J + 3 a z S b s b o Q S + g u 8 e F D E q z / J m / / G b Z u D t j 4 Y e L w 3 w 8 y 8 I B F c G 9 f 9 d g o b m 1 v b O 8 X d 0 t 7 + w e F R + f i k r e N U M W y 7 w 5 j 8 6 L 8 + 5 8 L F s L T j 5 z C n / g f P 4 A z f m M 7 A = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " S s YK O O t 1 j T f i N k y H K 7 x 6 G F F E b x M = " > A A A B 6 H i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l E q M e i F4 8 t 2 A 9 o Q 9 l s J + 3 a z S b s b o Q S + g u 8 e F D E q z / J m / / G b Z u D t j 4 Y e L w 3 w 8 y 8 I B F c G 9 f 9 d g o b m 1 v b O 8 X d 0 t 7 + w e F R + f i k r e N U M W y 7 w 5 j 8 6 L 8 + 5 8 L F s L T j 5 z C n / g f P 4 A z f m M 7 A = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " S s YK O O t 1 j T f i N k y H K 7 x 6 G F F E b x M = " > A A A B 6 H i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l E q M e i F4 8 t 2 A 9 o Q 9 l s J + 3 a z S b s b o Q S + g u 8 e F D E q z / J m / / G b Z u D t j 4 Y e L w 3 w 8 y 8 I B F c G 9 f 9 d g o b m 1 v b O 8 X d 0 t 7 + w e F R + f i k r e N U M W y 7 w 5 j 8 6 L 8 + 5 8 L F s L T j 5 z C n / g f P 4 A z f m M 7 A = = < / l a t e x i t > g < l a t e x i t s h a 1 _ b a s e 6 4 = " v m P s / L 4 D T

Figure 3 :
Figure 3: Colloidal boomerang made of 15 RPY-particles above a no-slip boundary.The height, h, is the distance between the centre of the boomerang-point RPY-particle and the no-slip surface.The angle θ is that between the vertical axis and the unit vector, p, that is orthogonal to both boomerang arms.

Figure 4 :
Figure 4: (left panel) Height distribution of the boomerang above a no-slip surface (right panel) L2-error between the simulated distributions and the true distribution generated using MCMC.

Figure 5 :
Figure 5: (left panel) Orientational distribution of the boomerang above a no-slip surface (right panel)L2-error between the simulated distributions and the true distribution generated using MCMC.

Figure 7 :
Figure 7: A snapshot from a rod simulation with channel width Ly = L/2.Each rod has a randomlyselected greyscale colouring and the semi-transparent grey boundaries depict the locations of the slip boundaries.

Figure 9 :
Figure 9: The equilibrium distributions of rod orientations in slip channels of various widths, shown in both (left) linear and (right) logarithmic scales.

Figure 10 :
Figure 10: An illustration of a CH particle formed of 19 FCM-particles of radius a (left), and the different interaction potentials for suspended CH particles (right).

Figure 11 :
Figure 11: The effective viscosity of both types of suspension for different dimensionless shear rates (top), and the radial distribution functions for the different 'repel-attract' simulations (bottom).The bottom right plot is displaying the same data as the bottom left plot but over a smaller r range.The colour of the curves fades from black to white as f0 decreases.

2 ak D − 1 jb∂N 1 / 2
to cast eq.(C.2) in Stratonovich form asdu i = k B T D −1ia N Theorem 5.1 of Malham and Wiese[71] (recognizing the ξ as the columns of N 1/2 ) we find that the corresponding orientation quaternion satisfies the Stratonovich equationdq i = Ψ ia k B T N 1/bk ∂u j − N ab D −1 jb ∂U ∂u j dt + 2k B T N 1/2 aj • dW j , (C.8)since Φ n kl = δ kl .In Appendix C we show that ∂ j D −1 u jb = 0 and hence the first term vanishes.Furthermore, applying the definition of partial differentiation yields ∂ j D −1 u ia n = − 1 2 ija , and hence the components of the third term which don't vanish in the differentiation vanish as the inner product of a symmetric and skew-symmetric matrix.Thus we obtain the simpler expression ∆x n = k B T ∆t (∂ x • N) n (E.5)