# Hydrodynamic interactions in colloidal ferrofluids: A lattice Boltzmann study

## Abstract

We use lattice Boltzmann simulations, in conjunction with Ewald summation methods, to investigate the role of hydrodynamic interactions in colloidal suspensions of dipolar particles, such as ferrofluids. Our work addresses volume fractions of up to and dimensionless dipolar interaction parameters of up to . We compare quantitatively with Brownian dynamics simulations, in which many-body hydrodynamic interactions are absent. Monte Carlo data are also used to check the accuracy of static properties measured with the lattice Boltzmann technique. At equilibrium, hydrodynamic interactions slow down both the long-time and the short-time decays of the intermediate scattering function , for wavevectors close to the peak of the static structure factor , by a factor of roughly two. The long-time slowing is diminished at high interaction strengths whereas the short-time slowing (quantified via the hydrodynamic factor ) is less affected by the dipolar interactions, despite their strong effect on the pair distribution function arising from cluster formation. Cluster formation is also studied in transient data following a quench from ; hydrodynamic interactions slow the formation rate, again by a factor of roughly two.

## 1 Introduction

The tendency of ferromagnetic colloidal particles to form aggregated structures, stabilized by anisotropic dipole-dipole interactions, was insightfully discussed by Pierre-Gilles de Gennes and Philip Pincus in 1970 [1]. This tendency is central to the structure [2], and hence the phase equilibria [3, 4, 5] and dynamics [6, 7], of magnetic colloids in organic solvents (ferrofluids). The earliest direct observations of chain-like structures were obtained by Hess and Parker using electron microscopy [8] but remarkably the quantitative experimental study of strong pair correlations began only in 2003 with cryogenic transmission electron microscopy studies by Philipse and co-workers [9, 10, 11]. This delay in confirming the classic predictions of de Gennes and Pincus partly reflects the extreme opacity of most ferrofluids, which precludes both direct microscopy of bulk phases and light scattering as methods to elucidate structure. Alongside X-ray and neutron scattering, these two methods have been central to the widespread progress made in understanding other forms of colloidal aggregation since 1970 [12, 13, 14].

The same experimental difficulties have also made it hard to study structural relaxation, which in macroscopically isotropic fluids at equilibrium can be quantified by the time-dependent correlator (directly accessible in inelastic scattering experiments, where available) [15, 16]

(1) |

where . Difficulties with scattering methods mean that many relaxation studies on ferrofluids have been limited to strictly properties such as frequency dependent bulk magnetic susceptibilities [17, 18] or magnetoviscous and rheological properties [19].

It is increasingly possible for computer simulation methods to fill the gaps in our experimental knowledge of structural relaxation in complex fluids [20, 21, 22, 23, 24, 25]. However, for ferrofluids, such simulations have previously been oversimplified in their neglect of hydrodynamic interactions between particles. These interactions are mediated by the intervening solvent – an essentially incompressible, Newtonian fluid of viscosity and density . Previous simulations using molecular dynamics (MD) [26, 27, 28, 29], Monte Carlo (MC) [30, 31, 32], and Brownian dynamics (BD) [7, 33, 34] of dipolar fluids, while each capable of generating the correct equilibrium statistics governed by the Boltzmann distribution, are all compromised by either the neglect or incomplete treatment of thermal noise and many-body hydrodynamics. In particular BD, while capturing the overdamped, diffusive dynamics of an isolated particle (with diffusion constant , being the particle radius) fails to correlate the Brownian motion of two or more particles in the correct manner. As a result, will have the wrong time dependence.

In this work, we present simulation results for dipolar colloids generated using the lattice Boltzmann (LB) method. This approach treats full hydrodynamic interactions, at least for the colloids simulated here, which have a sufficiently repulsive soft-core potential to ensure that particles do not approach one another too closely. (The presence of such a potential avoids large hydrodynamic lubrication forces in thin fluid films between two solid particles in close contact, whose treatment within LB is possible, but costly [35].) Alongside LB, other methods that fully treat hydrodynamics include force methods [36], Stokesian dynamics (SD) and accelerated Stokesian dynamics (ASD) [37, 38], and stochastic rotation dynamics [39]. (Dissipative particle dynamics does so also, but without proper control of noise terms as required here [40].) Few of these methods have yet been used to treat systems with long-range interactions, although an SD method has been applied to ferrofluids in shear flow [41, 42], and ASD was recently used to address charge-stabilized colloids at low ionic strength [16]. We are aware of no systematic application of these methods to the equilibrium dynamics of dipolar fluids in three dimensions.

This paper is organized as follows. In section 2 we outline the numerical methodology, and in section 3 present results for equilibrium structure as calculated by BD, LB, and MC methods. (These should be, and nearly are, identical.) Then in section 4 we present results for and its orientational analogs, focusing on a like-for-like comparison between BD (no hydrodynamics) and LB (full hydrodynamics). In section 5 we present an analysis of transient behavior, addressing the time evolution of the static structure including an analysis of cluster statistics. Finally in section 6 we give our conclusions and discuss prospects for future work.

## 2 Simulation methods

We use the lattice Boltzmann method for a fluid incorporating spherical solid particles [35, 43]. In this method, the density, momentum, and stress in the fluid are associated with various moments of a kinetic distribution function , defined at each site on a 3D lattice and acting on a space of 19 discrete velocities that connect neighboring sites in one timestep (including the null velocity). Setting the lattice parameter, timestep, and mean fluid density all to unity defines a set of LB units that we use in results quoted below. (Of course, many physical quantities can be expressed in dimensionless forms, from which these units cancel out.)

Each of our spherical colloidal particles has a hard-core radius and resides off-lattice; its surface then cuts the lattice bonds at a set of links at which fluid and particle interact by momentum transfer. This transfer is achieved by a ‘bounce-back on links’ algorithm [35]. The force and torque on a colloidal particle are found by summing contributions across the boundary links; particle velocities and angular velocities are then updated via a standard molecular dynamics routine [35, 43]. The thermal noise, responsible for colloidal Brownian motion, is generated entirely within the fluid and is fully included in the description of the fluid momentum [44], using a method reported previously [45]. Momentum transport then causes the random forces (and torques) felt by different particles to become correlated, in accord with the fluctuation-dissipation theorem which relates random forces to the matrix of particle mobilities . This matrix obeys where is the velocity of particle in response to a force on particle , and are cartesian indices. In the presence of hydrodynamic interactions, is not diagonal in particle indices, but has long-range correlations. In contrast to SD-based methods, LB avoids explicit computation of .

Each colloidal particle has the same mass as the nominal volume of fluid that it displaces; is the volume of a colloidal particle. The fluid viscosity is set to and the temperature to ; these choices represent the best compromise we have found between numerical accuracy and efficiency for the systems under study here. (For a discussion of parameter optimization in LB see [22, 23].) Based on these parameters, the particle velocity relaxation time is , and the time scale for fluid momentum to equilibrate around a particle is . The single-particle diffusivity is , giving a diffusive timescale of . This offers a reasonably wide domain in which to study short time diffusion, even in the case of strongly interacting particles (as here) where the short-time diffusion time window closes off at with a typical surface-to-surface separation between adjacent particles. Our choices for other (reduced) simulation parameters – to be defined below – correspond to typical values for real ferrofluids. Time-dependent results are expressed in units of , which means that for , they will be directly relevant to almost any real ferrofluid system.

To avoid computationally expensive lubrication contacts between particles (as mentioned in section 1) we introduce a short-range, soft-core repulsion acting at separations beyond the hard-core radius . (The latter coincides with the hydrodynamic radius [35, 43].) We tested various options for and found that a satisfactory choice is

where , is the surface-to-surface separation for spheres whose centers are apart, and is a short-range cutoff. We set the cutoff separation , which in terms of the hard-core diameter is . This comprises a truncated and shifted inverse power law potential; note that the interaction force remains divergent at contact (). Although structural and magnetic properties of ferrofluids are known to be quite insensitive to the choice of short-ranged repulsive potential [29], our choice could qualitatively describe the effects of a screened electrostatic repulsion. The case of steric stabilisation is more complex since the polymer layer will have, in addition, a direct effect on the hydrodynamic forces.

The introduction of reduces discretization errors in the noise forces which become acute when fluid nodes are excluded from the space between particles. (Since absence of fluid entails absence of noise, particles that are too close to one another to have fluid nodes between them effectively feel a reduced temperature.) In fact this issue, rather than avoidance of lubrication contacts per se, prevents efficient use of a shorter-range soft-core potential than the one we have chosen; the effect of using such a potential is to exaggerate the peak in the pair distribution function for particles at close contact (as though particles become colder at close separations). Even with our choice of , this deviation remains visible in the data of section 3 for the largest dipole strength used, but this is considered acceptable. Indeed, in the current state of the art for LB, errors of several percent from this and other sources such as shape discretization remain unavoidable. To reduce these further is straightforward in principle: one simply increases . However, for a fixed number of particles at a given volume fraction , the system volume must then be increased as , with a further increase in the run time to reach . Thus a factor 2 increase in , as would be required to give a worthwhile reduction in discretization errors, entails a 32-fold increase in computational resource.

The total colloid-colloid pair potential in our simulations is the sum of short-range soft-core () and long-range dipolar () contributions:

(2) |

Here, is the center-center separation vector for particles and , denotes a unit vector pointing along the dipole of particle , and where . The dipole-dipole interaction potential is written

(3) |

where is a unit vector pointing along the center-center separation vector. The interaction parameter is a dimensionless, dipolar coupling constant defined such that two colloids at hard core contact (), with dipoles mutually aligned and parallel to , have . This ‘nose-to-tail’ parallel conformation is of lowest energy; the other minimum-energy conformation is ‘side-by-side’ antiparallel, for which (with corresponding energy maxima on reversing the direction of one dipole). In experimental ferrofluids values up to are readily available; much larger ones can be achieved but are relatively unusual [9, 10, 11].

All of our simulations were performed in a cubic box of side with periodic boundary conditions for both fluid and particles. A standard Ewald summation technique was deployed to handle the long-range aspect of the dipolar interactions [46, 47]. The Ewald sum computes dipolar forces directly in real space for particle pairs with separation , a cutoff distance, and deals with the remainder in reciprocal space. To allow a parallel implementation using domain decomposition, which is important given the relatively high computational requirements associated with the LB fluid [43], we chose in all cases. This is small enough that some parallelism is possible, but not so small that the number of terms required in the reciprocal space sum for acceptable accuracy becomes unwieldy. A convergence factor [47] of was chosen, with wavevectors with and for and , respectively. These parameters were optimised using established methods [48]. Finally, following normal practice, the Ewald sum boundary condition at infinity was chosen to be “conducting”, i.e., the dipolar (magnetic) susceptibility of the surroundings is infinite. This removes a zero-wavevector, bulk-magnetization term from the Hamiltonian which, in the opposite extreme of vacuum boundary conditions, can lead to unphysically small polarised domains and hence slow down simulation convergence [49].

For comparison with our LB results, a BD algorithm was set up within the colloidal MD module of our LB code, deploying standard techniques to generate the required Langevin dynamics with independent noise acting directly on each particle. The inertia of the particles is retained but the many-body hydrodynamics are replaced by a Stokes drag that is independent of the location of other particles. By setting exactly the same value for , we thus create an algorithm which differs from LB solely in the omission of many-body hydrodynamics.

To further validate both codes, and to monitor the achievement of Boltzmann equilibrium, canonical Monte Carlo simulations were performed in a cubic simulation cell with periodic boundary conditions applied [47]. The long-range dipolar interactions were handled using the Ewald summation with conducting boundary conditions, a convergence factor , and wavevectors with . The maximum translational and orientational displacement parameters were adjusted independently to give acceptance rates of approximately and , respectively; it is efficient to employ low acceptance rates for translations of particles with hard cores, due to the possibility of rapidly identifying overlaps. For each state point considered, we performed equilibration runs of MC cycles, where one MC cycle consisted of, on average, one attempted translation or rotation per particle. Production runs consisted of MC cycles.

Our LB work was performed primarily on lattices of size or . The colloid volume fraction is given by , where is the number of colloids. For the larger system size at a volume fraction , a run of timesteps required 56 hours on 64 cores of a cluster of 3GHz Intel dual-core processors. For resource reasons most of our results on fully equilibrated samples (- at ) concern the smaller . BD and MC runtime requirements were modest in comparison.

## 3 Static structure

In this section we confirm that our LB algorithm generates, to acceptable accuracy, the Boltzmann distribution for thermal equilibrium properties, as does our BD code, and that both are in agreement with MC data. This is of course necessary if the dynamical data in subsequent sections are to be trusted. Such agreement is not automatic but is in fact a very demanding test of the ability of our LB algorithm to generate correlated noise forces satisfying the fluctuation dissipation condition. In particular, without adopting the methods of [45] (in which noise terms are applied to not only the hydrodynamic but also the local modes of the fluid degrees of freedom) we would not be confident of achieving such agreement. Even with these methods, LB parameter values must be carefully chosen to maintain acceptable performance. We have already mentioned that correct treatment of noise was the limiting factor in our choice of ; it also prevents us using a much larger value (which sets the intrinsic noise level in the LB fluid) and/or a much smaller viscosity . Either step could in principle drastically reduce the run time required to reach the basic timescale for colloidal diffusion. Thus the control of discretization errors within the noise sector currently remains the primary efficiency bottleneck in the application of large-scale LB simulations to colloidal diffusion.

### 3.1 Energy equilibration

Table 1 shows time-averaged energy data for various simulation runs. Those for have fully equilibrated (run times ) as judged by convergence of the energy parameters to their long-term averages. The results indicate excellent agreement between LB and the other methods at and , and adequate agreement (errors of less than 5%) at . Some runs with the larger system size, , are also reported in the Table. These have run times and while their energies appear to have saturated, those at are showing continued structural evolution by other measures (such as cluster statistics; see section 5). Accordingly, the reported energy discrepancies for these runs may include systematic errors arising from incomplete structural equilibration, and should not be taken as a guide to the relative accuracy of the LB algorithm.

Method | ||||||
---|---|---|---|---|---|---|

0 | 0.10 | 529 | LB | - | ||

529 | BD | - | ||||

529 | MC | - | ||||

4 | 0.10 | 529 | LB | |||

529 | BD | |||||

529 | MC | |||||

8 | 0.10 | 529 | LB | |||

529 | BD | |||||

529 | MC | |||||

4 | 0.20 | 8239 | LB | |||

8239 | BD | |||||

529 | MC | |||||

8 | 0.20 | 8239 | LB | |||

8239 | BD | |||||

529 | MC |

### 3.2 Radial distribution functions (RDFs)

In Figures 1 and 2 we plot the radial distribution function , and the projections of the ‘molecular’ pair distribution function onto rotational invariants [50, 51], as measured in various LB runs:

(4) | |||||

(5) | |||||

(6) | |||||

(7) |

Note that all the runs reported have . We do not report equilibrium structural data for here because, even though the energy can be equilibrated for , the modest number of particles combined with a very long autocorrelation time means that good statistics cannot be gathered for structural quantities even with millions of timesteps. For the statistics are better due to larger , but the maximum available runtime () is not long enough to guarantee equilibration as discussed above.

Each of the above RDFs could equally be plotted in Fourier space (with then transforming into the static structure factor ), but the real space versions offer the more sensitive tests of equilibration. This is because, for the reasons already discussed, any errors are likely to occur in a localized range of at close contact (). We find excellent agreement between LB, BD, and MC for at (data not shown). Figure 1 shows adequate agreement between all methods at although BD and LB both show slight discrepancies from the MC data (which should be the most accurate) in the neighborhood of the first peak for each correlator. (A slight discrepancy for LB is also detectable near the first minimum of .)

(a) | (b) |

(c) | (d) |

For – Figure 2 – there is a clear discrepancy between LB and the other two methods, overestimating by 10% or so the height of the first peak in all the RDFs. This error is consistent in sign and magnitude with the energy discrepancies in Table 1, and with the specific type of noise discretization error reported in section 2. Given the other sources of error in LB [22], we consider it acceptable.

(a) | (b) |

(c) | (d) |

Overall, the observed behavior of and is in good accord with that established in earlier studies [32, 51]. With , shows only short-range correlations (data not shown); with , the primary peak in becomes very strongly pronounced due to the high degree of particle association to form, at these low densities, chain-like aggregates. helps distinguish parallel from anti-parallel correlations between the dipole moments, contains (minus) the dipolar potential, and picks out ‘nematic’ orientational ordering. All of these functions show positive peaks at short range – at intervals close to – confirming the prevalence of the ‘nose-to-tail’ parallel conformations of nearby dipole moments within chains. As expected, the peaks become more pronounced with increasing . For a visual confirmation of the chaining, Figure 3 shows two snapshots from equilibrated LB simulations of -particle systems with and at . Each particle is color-coded to reflect the total number of particles in the cluster to which it belongs; monomers, clusters with - particles, and clusters with particles are given unique colors. (See section 5.2.) The chain-like structural motif is clearly visible in both systems. Figures 1 and 2 show that interparticle correlations are short-ranged (compared to the box dimensions), and hence the thermodynamic properties should not show any pronounced finite-size effects. Nonetheless, the cluster network in Figure 3 appears to span the simulation cell, and so we might anticipate some finite-size effects in the long-time, long-wavelength dynamics. We have simulated the largest possible system sizes throughout.

(a) | (b) |

## 4 Dynamic correlators in equilibrium

We now present results for the intermediate scattering function , and its orientational analogs. As for the static structure, we restrict attention to the case where we can combine complete equilibration with adequate statistical averaging. We examined relaxations in the Fourier components of the number density and magnetization density at wavevectors commensurate with the periodic boundary conditions; to improve the statistics, we averaged the appropriate correlators at wavevectors of the same magnitude .

### 4.1 Intermediate scattering function

Figure 4 shows (eq 1) as a function of for , , and at . In each case curves are plotted for three different wavevectors; one close to the peak () in the static structure factor , one larger, and one smaller. These linear-linear plots allow one to see clearly the long-time relaxation of the structure. (The short time dynamics is considered in section 4.3 below.)

(a) | (b) |

(c) | (d) |

In all cases, the effect of hydrodynamic interactions is to slow down the relaxation of ; this effect is most marked for wavenumbers well below (which in any case relax more slowly than those at the peak). The long-time relaxations are not far from exponential in all cases, and in particular show no sign of decomposition into separate and relaxation processes as expected in colloidal systems on approach to a glass transition [15]. Absence of the latter is confirmed by plotting the time on a logarithmic scale, as presented for in Figure 4(d). It is notable, however, that the slowing by hydrodynamic interactions of the long-time relaxation, at least for , is much diminished at strong dipolar interactions (). This might be taken as evidence that structural rearrangement in this case is controlled mainly by the energetics of aggregate rearrangement (breaking and reformation of dipolar contacts), and is no longer limited by the rate at which solvent can flow around the evolving structure – a state of affairs generally expected to hold for glassy colloids. Caution is needed before drawing such a conclusion, however; many authors would, on adopting this reasoning, expect BD and LB curves to superpose only after rescaling of time by the short-time diffusion constant at the peak, [52]. As discussed in section 4.3, hydrodynamic interactions continue to cause a factor of 2 change in this quantity even for .

### 4.2 Orientational relaxation

Defining a wavevector-dependent dipole density we can construct orientational correlators from the longitudinal (L) and transverse (T) components and [53]:

(8) | |||||

(9) | |||||

(10) |

Data for and at are plotted in Figure 5 for , , and . For clarity, we omit since this is a simple average of the longitudinal and transverse parts. The results are broadly comparable to the relaxation of at similar wavevectors. We note that the longitudinal relaxations are slower than the transverse ones. This might be ascribed to the slow rotational diffusion of chain orientations with respect to the wavevector , as compared to faster librational motions of dipoles perpendicular to the local chain orientation. Figure 5(d) shows (for ) the dependence of . This is again comparable to that for the density relaxation. Note, though, that is not a conserved quantity and therefore, unlike the density, is not compelled to relax slowly for . The fact that it does so suggests that is enslaved to slow particle rearrangements, as would arise if the dipole moments inside a cluster were to adopt frozen orientations relative to the positions of the constituent particles over the cluster’s lifetime. With a dipolar bonding energy of for two linearly aligned dipoles, such behavior is quite plausible.

(a) | (b) |

(c) | (d) |

### 4.3 Short time diffusion

The shape of is partly characterized by a short-time collective diffusion constant

(11) |

where denotes a measurement taken at time scales long enough that a single particle indeed moves diffusively, but short enough that its average displacement remains small compared to (or, if smaller, the surface-to-surface separation from neighboring particles). For an isolated particle we therefore require where we recall that is the time scale for steady fluid motion to be established at the particle scale, is the velocity autocorrelation time of the particle (of mass ), and is the time for a particle to diffuse its own radius. For particles with caging or bonding at surface-to-surface separations , the requirement is replaced by . Defining

(12) |

one finds that in the absence of hydrodynamic interactions the ‘hydrodynamic factor’ is always unity for all , whereas experiments on, e.g., hard sphere colloids show values that are not only smaller but also -dependent [16, 54]. For example, in hard-sphere colloids, at and [16, 54], while at .

The numerical evaluation of , and hence of , carries significant difficulties associated with finite size corrections [36]. That is, the long range nature of the hydrodynamic interactions, in conjunction with periodic boundary conditions, makes very slow to converge with system size , or equivalently with particle number at fixed . For the case of hard spheres, at least, this can be brought under good control at a semi-empirical level by adopting the following correction [36, 54]

(13) |

where is the so-called high-frequency viscosity of the suspension and is the solvent viscosity.

To evaluate eq 13 we numerically measured the high frequency viscosities for fully equilibrated systems with the given repulsive short-range potential and dipolar long-range interaction, using the recipe by Ladd [44]. This amounts to calculating the integrated stress-stress correlation in a time window that is long enough to relax fluid degrees of freedom but too short for the colloids to move significantly. For , , and we found to be , , and , respectively. A similar procedure was used in [16] for the case of colloids with long range coulombic repulsions. However, since eq 13 was invented to account for the observed finite-size behavior of systems of hard spheres, its use in other systems remains empirically questionable. Below we therefore present data both for as actually measured and for as estimated via eq 13, but use the latter value to calculate .

The above caveat applies particularly when long-range (e.g., dipolar) interactions are present. Arguably such interactions should create their own finite size corrections, somewhat akin to those from hydrodynamics. In this case one might expect that, even with hydrodynamics switched off, the measured would show size-dependent deviations from unity. In the data reported below we indeed find values significantly less than unity for BD at large ; however, we know of no method to correct for this and make no attempt to do so.

Figure 6 shows representative () short time data for the three values of studied at . In accordance with expectation, the regime of short time diffusion is established beyond a few hundred timesteps, and for and there is thereafter a wide region of exponential decay within which can be measured easily. For this window is foreshortened – which is not surprising since the short time regime should end on the timescale of particle collisions. (For high interaction strengths, particles are bonded to neighbors with which they collide frequently.) Nonetheless a reasonable numerical estimate of the decay rate can be made. In practice this was done by first identifying by eye the time window for short-time diffusion and then fitting to the log-linear plot within this window at each . Finally the data for distinct values were binned (each bin containing roughly ten wavevectors) and the statistical error then estimated for the binned data.

(a) | (b) | (c) |

Figures 7, 8, and 9 show plots of , , and generated from our dynamic datasets for , , and , respectively. We also show, for comparison, the uncorrected curves; data generated from MC to check accuracy; and direct comparison with our BD results for and . For the BD data no finite-size correction was made; for we recover to simulation accuracy, with smaller values at larger presumably attributable to finite size effects in the thermodynamic sector, as discussed above. (It is possible that, were these to be corrected, the curves for LB could depend less strongly on than in the results shown here.)

(a) | (b) |

(c) | (d) |

(a) | (b) |

(c) | (d) |

(a) | (b) |

(c) | (d) |

When hydrodynamics is switched on, we obtain values of for all three values that are comparable to previous data on hard sphere colloids at [16, 54]. However, for large , and both suggest a rising trend at small wavenumbers (). The rise in at low is consistent with the formation of large dipolar clusters (and specifically, chains [31]). Clustering should also reduce the hydrodynamic friction per particle – as is familiar from the fact that large clusters sediment more quickly under gravity. That is, the body force increases linearly with particle number whereas the viscous friction scales with hydrodynamic radius which is generally sublinear. This reduction is consistent with the observed rise in at low . Note however that this rise, although detectable beyond the scatter in the data itself, is sensitive to the treatment of finite-size corrections and is therefore provisional.

## 5 Transient dynamics and cluster formation

We now consider the evolution of the structure following an initial quench, at time zero, from an equilibrium state with to one with either or . We have studied such quenches at , and .

### 5.1 Transient dipolar energy

Figure 10 shows the relaxation of the dipolar energy for each volume fraction as a function of time following the quench; LB and BD data are directly compared. In all cases, the effect of hydrodynamic interactions is to slow the approach to equilibrium. However, the effect is quite modest, and comparable to that reported earlier for in equilibrium. The relaxation time is increased by no more than roughly a factor two, even for . Note that addition of hydrodynamics is by no means guaranteed to slow down, rather than speed up, the approach to equilibrium. A familiar counterexample is binary fluid phase separation, where fluid flow of the two species creates a less dissipative, and hence faster, phase-separation route than pure diffusion at intermediate and late times [20, 24, 25].

(a) | (b) |

### 5.2 Cluster statistics

We define two dipolar particles to be in a bonded configuration if their pair dipolar interaction energy from eq 3 obeys

(14) |

This definition is somewhat arbitrary: in principle once could choose either an energy-based or a geometric criterion. Our choice corresponds to an energy criterion set by an equipotential surface, in configuration space, of the dipolar part of the interaction. This is more suitable than a criterion based solely on : the latter would count as a bond any close encounter between dipoles even if their orientation was such as to create a strongly repulsive force. In addition, our choice is designed to capture end-to-end bonding but reject most encounters between antiparallel dipoles even when their orientation is such as to create a bond. (Because of the short-range repulsion, the energy minimum for such bonds lies above the threshold in eq 14.) The particular value of the energy threshold – which is intermediate between those used in earlier studies [55, 56] – gave cluster distributions in good accord with what was expected from visual inspection of simulation snapshots, and was sufficient for the current purpose of examining transient cluster formation.

Using eq 14 we partition each configuration of particles into a set of disjoint clusters, and monitor the fraction of particles that are assigned to clusters of size . The time evolution of gives information about the growth of clusters following the quench from at . Figures 11 and 12 show data for various at volume fractions and . Once again, BD data is included for comparison. (The data is binned timewise with a stride of 25 timesteps for and 50 timesteps thereafter. This choice offered the best compromise between smoothness and sensitivity. The actual numbers of clusters are low, and the relative fluctuations are high, so it is not easy to iron out the noise.) The transient dynamics is subject to a similar slowing by hydrodynamic interactions as was the energy transient. Other than this there are no obvious differences between the LB and BD data. For large , both show characteristically peaked plots for , , and as small clusters build up and are then subsumed into larger ones.

(a) | (b) |

(c) | (d) |

(a) | (b) |

(c) | (d) |

Finally, in Figure 13 we show the mean number of particles per cluster, , as a function of time for the same runs. The same hydrodynamic slowing is evident. Mean cluster sizes in the LB simulations are slightly higher than those in the BD simulations. This can be traced back to the small errors in dealing with particles close to contact, leading to more pronounced near-neighbour correlations, as discussed in section 3.

(a) | (b) |

## 6 Summary and Conclusions

In this paper we have presented results for the equilibrium and transient dynamics of dipolar colloids with many-body hydrodynamic interactions. These were gained by incorporating the Ewald summation for the long-range dipolar interactions into our existing lattice Boltzmann algorithms, which handle the hydrodynamic forces by explicit propagation of momentum across the fluid residing on a lattice. The colloidal particles themselves are off-lattice and undergo molecular dynamics; Brownian motion is caused by fluctuating momentum transfer from the surrounding solvent, which creates correlated noise of the kind demanded by the fluctuation-dissipation theorem for hydrodynamically interacting particles. The full, fluid-driven noise can easily be replaced by local noise with no correlation between particles, creating a BD code. The resulting comparison of the LB and BD results allows the effect of many-body hydrodynamics to be isolated.

At the volume fractions () and interaction strengths (, ) studied here, these effects are easily measurable but remain relatively modest. Quantitative shifts in both short-time and long-time diffusion were observed for wavevectors near and below the peak in the static structure factor. Likewise, we found shifts in transient relaxation rates for the cluster size distribution on approach to steady state following a quench from . In all cases, the system with hydrodynamic interactions relaxes more slowly than the equivalent BD system. However, for the range of volume fraction and interaction strengths studied, this slowing down was rarely by more than a factor of two.

Although it is possible that stronger hydrodynamic effects would be observed in the dynamics of a quiescent system at larger and , their effects on long-time relaxation appear already to be decreasing for the largest values studied here. This could be a precursor to entering a glassy regime in which the crossing of local energy barriers limits relaxation rates; within this regime, conventional wisdom holds that hydrodynamics affects relaxational dynamics only through a scale factor [52]. However, despite the observation of slow transients and the difficulty of attaining full equilibration, even at and , we find no direct evidence for a glassy regime; specifically we see no separate and relaxation processes. This is not surprising as our simulations run for at most , and a truly glassy system would certainly not approach equilibrium, as ours do, on this time scale.

Within our LB framework, there are serious obstacles to achieving much longer physical timescales using reasonable computational resources. One bottleneck remains the accurate treatment of noise; currently this requires a very large separation (order in our runs) between the simulation timestep and . Future algorithmic work will, we hope, partially address this issue.

Our simulations were engineered to avoid the very large hydrodynamic forces that arise when hard colloidal particles come into lubrication contact. This was done by including a soft-core repulsion to maintain adequate separation between particle surfaces even when strongly bonded by the dipolar interactions. It is possible that these lubrication effects could enhance the relative role of hydrodynamic interactions, by further slowing the timescale for bond breakage and re-formation. To address this effect specifically (in a bulk periodic system) an algorithm such as ASD might be more suitable than LB. Note that within LB one can include a routine to address lubrication forces via an SD-like algorithm, but the computational scaling becomes bad when there are large clusters of particles in mutual lubrication contact. We do not know how well ASD would perform under such conditions, as compared to the purely repulsive interactions studied in [16].

Even without lubrication, the effects of many-body hydrodynamics could, of course, also become much more pronounced in various nonequilibrium situations. These include the rheological response to steady and/or time-dependent shearing, and perhaps the nonlinear response to large orienting fields. We hope to address one or more of these topics in future work.

### Acknowledgments

This work was funded in part under EPSRC Grants GR/S10377/01 and EP/C536452/1 (RealityGrid). We thank ECDF (Edinburgh Compute and Data Facility) for computational resources. EK thanks SUPA and ORS for a studentship. MEC holds a Royal Society Research Professorship.

### References

- de Gennes, P. G.; Pincus, P. A. Phys. Kondens. Materie 1970, 11, 189-198.
- Holm, C.; Weis, J.-J. Curr. Opin. Colloid Interface Sci. 2005, 10, 133-140.
- Camp, P. J.; Shelley, J. C.; Patey, G. N. Phys. Rev. Lett. 2000, 84, 115-118.
- Tlusty, T.; Safran, S. A. Science 2000, 290, 1328-1331.
- Ganzenmüller, G.; Camp, P. J. J. Chem. Phys. 2007, 126, 191104.
- Murashov, V. V.; Camp, P. J.; Patey, G. N. J. Chem. Phys. 2002, 116, 6731-6737.
- Duncan, P. D.; Camp, P. J. Phys. Rev. Lett. 2006, 97, 107202.
- Hess, P. H.; Parker Jr., P. H. J. Appl. Polymer Sci. 1966, 10, 1915-1927.
- Butter, K.; Bomans, P. H. H.; Frederik, P. M.; Vroege, G. J.; Philipse, A. P. Nature Materials 2003, 2, 88-91.
- Butter, K.; Bomans, P. H.; Frederik, P. M.; Vroege, G. J.; Philipse, A. P. J. Phys.: Condens. Matter 2003, 15, S1451-S1470.
- Klokkenburg, M.; Dullens, R. P. A.; Kegel, W. K.; Erné, B. H.; Philipse, A. P. Phys. Rev. Lett. 2006, 96, 037203.
- Poon, W. C. K. J. Phys.: Condens. Matter 2002, 14, R859-R880.
- Scheffold, F.; Schurtenberger, P. Soft Materials 2003, 1, 139-165.
- Jenkins, M. C.; Egelhaaf, S. U. Adv. Colloid Interfac. 2008, 136, 65-92.
- Pusey, P. N. Colloidal Suspensions. In Liquids, Freezing and Glass Transition, Les Houches Session LI, Volume II; Hansen, J.-P.; Levesque, D.; Zinn-Justin, J., Eds.; North Holland: Amsterdam, 1991.
- Banchio, A. J.; Nägele, G. J. Chem. Phys. 2008, 128, 104903.
- Erné, B. H.; Butter, K.; Kuipers, B. W. M.; Vroege, G. J. Langmuir 2003, 19, 8218-8225.
- Huke, B.; Lücke, M. Rep. Prog. Phys. 2004, 67, 1731-1768.
- Odenbach, S. J. Phys.: Condens. Matter 2004, 16, R1135-R1150.
- Kendon, V. M.; Cates, M. E.; Pagonabarraga, I.; Desplat, J.-C.; Bladon, P. J. Fluid. Mech. 2001, 440, 147-203.
- Kremer, K. Macromol. Chem. Phys. 2003, 204, 257-264.
- Cates, M. E.; Stratford, K.; Adhikari, R.; Stansell, P.; Desplat, J.-C.; Pagonabarraga, I.; Wagner, A. J. J. Phys.: Condens. Matter 2004, 16, S3903-S3915.
- Stratford, K.; Adhikari, R.; Pagonabarraga, I.; Desplat, J.-C.; Cates, M. E. Science 2005, 309, 2198-2201.
- Stansell, P.; Stratford, K.; Desplat, J.-C.; Adhikari, R.; Cates, M. E. Phys. Rev. Lett. 2006, 96, 085701.
- Stratford, K.; Desplat, J.-C.; Stansell, P.; Cates, M. E. Phys. Rev. E 2007, 76, 030501(R).
- Wang, Z.; Holm, C.; Müller, H. W. Phys. Rev. E 2002, 66, 021405.
- Wang, Z.; Holm, C. Phys. Rev. E 2003, 68, 041401.
- Huang, J. P.; Wang, Z. W.; Holm, C. Phys. Rev. E 2005, 71, 061203.
- Ivanov, A. O.; Kantorovich, S. S.; Reznikov, E. N.; Holm, C.; Pshenichnikov, A. F.; Lebedev, A. V.; Chremos, A.; Camp, P. J. Phys. Rev. E 2007, 75, 061405.
- Weis, J. J.; Levesque, D. Phys. Rev. Lett. 1993, 71, 2729-2732.
- Camp, P. J.; Patey, G. N. Phys. Rev. E 2000, 62, 5403-5408.
- Levesque, D.; Weis, J. J. Phys. Rev. E 1994, 49, 5131-5140.
- Mériguet, G.; Jardat, M.; Turq, P. J. Chem. Phys. 2004, 121, 6078-6085.
- Mériguet, G.; Jardat, M.; Turq, P. J. Chem. Phys. 2005, 123, 144915.
- Nguyen, N.-Q.; Ladd, A. J. C. Phys. Rev. E 2002, 66, 046708.
- Ladd, A. J. C. J. Chem. Phys. 1990, 93, 3484-3494.
- Phung, T. N.; Brady, J. F.; Bossis, G. J. Fluid Mech. 1996, 313, 181-207.
- Banchio, A. J.; Brady, J. F. J. Chem. Phys. 2003, 118, 10323-10332.
- Malevanets, A.; Kapral, R. J. Chem. Phys. 1999, 110, 8605-8613.
- Groot, R. D.; Warren, P. B. J. Chem. Phys. 1997, 107, 4423-4435.
- Satoh, A.; Chantrell, R. W.; Coverdale, G. N.; Kamiyama, S. J. Colloid Int. Sci. 1998, 203, 233-248.
- Satoh, A.; Chantrell, R. W.; Coverdale, G. N. J. Colloid Int. Sci. 1999, 209, 44-59.
- Stratford, K.; Pagonabarraga, I. Comput. Math. Appl. 2008, 55, 1585-1593.
- Ladd, A. J. C. J. Fluid Mech. 1994, 271, 285-309.
- Adhikari, R.; Stratford, K.; Cates, M. E.; Wagner, A. J. Europhys. Lett. 2005, 71, 473-479.
- de Leeuw, S. W.; Perram, J. W.; Smith, E. R. Proc. R. Soc. Lond. A 1980, 373, 27-56.
- Allen, M. P.; Tildesley, D. J. Computer simulation of liquids; Clarendon Press: Oxford, 1987.
- Fincham, D. Information Newsletter for Computer Simulation of Condensed Phases 1993, 38, 17-24.
- Wei, D.; Patey, G. N. Phys. Rev. Lett. 1992, 68, 2043-2045.
- Hansen, J.-P.; McDonald, I. R. Theory of simple liquids; Academic Press: London, 1986.
- Weis, J. J.; Levesque, D. Phys. Rev. E 1993, 48, 3728-3740.
- Fuchs, M.; Ballauff, M. J. Chem. Phys. 2005, 122, 094707.
- Felderhof, B. U.; Jones, R. B. Phys. Rev. E 1993, 48, 1084-1090.
- Segrè, P. N.; Behrend, O. P.; Pusey, P. N. Phys. Rev. E 1995, 52, 5070-5083.
- Stevens, M. J.; Grest, G. S. Phys. Rev. Lett. 1994, 72, 3686-3689.
- Tavares, J. M.; Weis, J. J.; Telo da Gama, M. M. Phys. Rev. E 1999, 59, 4388-4395.