Dynamics of Hard Colloidal Cuboids in Nematic Liquid Crystals

We perform Dynamic Monte Carlo simulations to investigate the equilibrium dynamics of hard colloidal cuboids in oblate and prolate nematic liquid crystals. In particular, we characterise the particles' diffusion along the nematic director and perpendicularly to it, and observe a structural relaxation decay that strongly depend on the particle anisotropy. To assess the Gaussianity of their dynamics and eventual occurrence of collective motion, we calculate two- and four-point correlation functions that incorporate the instantaneous values of the diffusion coefficients parallel and perpendicular to the nematic director. Our simulation results highlight the occurrence of Fickian and Gaussian dynamics at short and long times, locate the minimum diffusivity at the self-dual shape, the particle geometry that would preferentially stabilise biaxial nematics, and exclude the existence of dynamically correlated particles.


I. INTRODUCTION
Colloids are systems of macromolecules or nanoparticles dispersed through a continuous medium. If the dispersed phase consists of solid particles and the continuous phase is a liquid, like an ink or a paint, then they are referred to as colloidal suspensions. The particles' ability of remaining suspended in the liquid stems from a thermal energy that must overcome the gravitational potential energy, which, by contrast, promotes sedimentation. Colloidal suspensions of anisotropic particles are complex fluids that display a particularly rich phase behaviour, with liquid crystals (LCs) forming between the isotropic and the crystal phase. LCs are states of matter with intermediate features between those of a simple, disordered liquid and those of a crystalline solid. More specifically, they flow like liquids, but their building blocks are able to orient along one or more spatial directions, thus closely reproducing the internal ordering of molecular crystals and their optical properties.
A renewed interest has been recently devoted to colloidal suspensions of anisotropic particles with a biaxial geometry, such as board-like (cuboids) and bent-core particles, whose orientation in space is fully identified by two independent orthogonal axes. Biaxial particles are, in principle, excellent candidates to form biaxial nematic (N B ) LCs, whose existence, predicted by a generalized Maier-Saupe theory proposed by Freiser exactly half a century ago [1,2], is still the object of a fervent research interest among experimental, computational and theoretical groups [3][4][5][6][7][8][9][10][11][12][13][14][15][16]. This interest is rooted in the captivating idea of manufacturing biaxial LCDs, by many envisaged as the next-generation fast display technology, * acuemen@upo.es † alessandro.patti@manchester.ac.uk but still at an embryonic stage due to the difficulty of unambiguously identifying stable molecular N B phases. By contrast, stable colloidal N B phases have been observed for the first time approximately ten years ago in systems of polydisperse prolate goethite particles [3].
Despite such an intense research activity aiming at fully unveiling the phase and aggregation behaviour of biaxial particles [17], much less attention has been devoted to their dynamics. On one hand, this is not surprising as, apart from the especially successful case of the Gay-Berne fluid [18] and the purely orientational pair potential proposed by Straley [19], most of the model particles with a biaxial geometry (e.g. cuboids) lack a suitable interaction potential to compute time trajectories and are often modelled as arrays of spherical beads [20]. On the other hand, it is not an easy task to predict how and, especially, over what time scales the N B phase would respond to external stimuli and hence quantify the advantages of employing them in LCDs, if an insight into their dynamical behaviour is missing. For instance, it is relevant to know how the shape anisotropy of biaxial particles determines their ability of diffusing in a crowded and ordered mesophase, such as a nematic LC; it is also important to know whether these particles are dynamically correlated as this would influence the response time of the whole phase to applied forces. In this context, the Dynamic Monte Carlo (DMC) simulation method plays a very important role as it provides a theoretical framework to mimic the dynamics of colloidal particles of any possible shape without integrating stochastic or deterministic differential equations as done, respectively, in Brownian Dynamics (BD) or Molecular Dynamics (MD) simulations [21][22][23]. The DMC method, which shows excellent quantitative agreement with BD simulations and can be between 10 and 20 times faster than BD [23], is especially useful to mimic the dynamics of hard-core particles, which cannot be straightaway reproduced by MD or BD simulations. Consequently, the dynamics of biax-ial particles interacting via a mere steric repulsion can be analysed by DMC simulations. The main goal of the present work is therefore investigating the equilibrium dynamics of prolate and oblate hard board-like particles (HBPs) in nematic LCs and characterise their structural relaxation decay. This paper is organised as follows. In Section II, we introduce the model and simulation methodology, referring the reader to our previous contributions for specific details on the implementation of the DMC technique. The relevant correlation functions to assess the dynamics are also discussed in this section. In Section III, we discuss the main features of the equilibrium dynamics of cuboids and, finally, in Section IV, we draw our conclusions.

II. MODEL AND SIMULATIONS
The colloidal cuboids studied in this work are modelled as hard board-like particles (HBPs), whose main geometrical features are shown in Fig. 1. They are rectangular prisms of thickness T , length L = 12T and width T ≤ W ≤ L. The particle thickness, T is our unit length and is kept constant, while the particle width, W , is a simulation parameter that varies between T and L in order to reproduce prolate (rod-like) and oblate (disklike) particles, respectively. At W = √ LT , HBPs are not oblate nor prolate, but display an ambivalent nature that can produce nematic LCs with oblate, prolate or biaxial symmetry. This special geometry is referred to as self-dual shape. The width-to-thickness ratio, W * ≡ W/T , and the system packing fraction, η ≡ N p v 0 /V , completely determine the phase behaviour and dynamics of these systems. In the latter definition, V is the volume of the simulation box, N p = 2000 the number of HBPs, and v 0 = 12W T 2 the particle volume. The phase diagram of monodisperse HBPs, studied by Monte Carlo (MC) simulations, reveals the existence of oblate and prolate nematic and columnar liquid-crystalline phases as well as biaxial smectics [13]. More recently, simulations of especially large systems suggested that the prolate columnar phase might actu-ally be metastable with respect to either a smectic phase or a crystal phase, depending on the system density [24]. Over the last two years, the quest for the elusive biaxial nematic (N B ) phase has produced a number of theoretical [4], simulation [10,14,15] and experimental [8] studies that confirm the difficulty of obtaining a stable N B phase of cuboids. The focus of the present contribution is investigating and characterising the relaxation dynamics of prolate (N + ) and oblate (N − ) uniaxial nematic phases of HBPs. To this end, we first performed standard MC simulations in the canonical ensemble to equilibrate the systems at the packing fraction, η = 0.34, at which N + and N − phases are known to be stable over the full set of aspect ratios studied here, namely 1 ≤ W * ≤ √ L * ≈ 3.46 for N + and 3.46 ≤ W * ≤ 12 for N − [13]. To study the relaxation dynamics of these systems, we performed Dynamic Monte Carlo (DMC) simulations in the canonical ensemble. Details of the DMC method are available to the interested reader in our former works [21][22][23]. Here, we briefly summarise the main points that are essential to gain a general insight.
DMC simulations are performed in cubic boxes with periodic boundaries and, since our intention is to mimic the time evolution of physical systems, no unphysical moves, such as cluster moves or particle swaps, are allowed. One DMC cycle consists of N p attempts of simultaneously displacing and rotating randomly selected HBPs. These combined moves are accepted or rejected according to the standard Metropolis algorithm, that is with probability min[1, exp(−β∆E)], where ∆E is the energy difference resulting from the particle's move and β the inverse temperature. Because the interactions between HBPs are modelled by a hard-core potential, moves are always accepted unless an overlap is produced. Overlap tests are based on the separating axes algorithm proposed by Gottschalk et al. [25] and later adapted by John and Escobedo to investigate the phase behaviour of colloidal cuboids with square cross section [26]. The position of the center of mass of a given particle j is updated by decoupling the displacement into three contributions: δr j = X Tûj + X Wvj + X Lŵj , whereû j ,v j , andŵ j are unit vectors parallel to the particle thickness, width and length, respectively. The magnitude of displacements along the three particle axes is selected within uniform distributions that depend on the particle translational diffusivities at infinite dilution, D tra α,j , with α = T , W , or L. In particular: where δt MC is the timescale for a single DMC cycle. As far as the rotations are concerned, we attempt to reorient the particle's axesŵ j ,v j , andû j by three consecutive rigid rotations around L, W , and T , respectively. These rotations must satisfy the following condition: where D rot α,j are the particle rotational diffusivities at infinite dilution. The diffusion coefficients D tra α and D rot α have been estimated by employing HYDRO ++ , a program that calculates the solution properties of macromolecules and colloidal particles by approximating their shape and volume with an array of spherical beads of arbitrary size [27,28]. Numerical values of the translational and rotational diffusion coefficients used in this work are given in Table I in units of D 0 ≡ T 2 τ −1 , where τ is the time unit. In order to compare the dynamics of the systems listed in Table I, it is crucial to recover the actual timescale of the Brownian dynamics (BD), δt BD [21][22][23].
To this end, one needs to rescale the DMC time scale, which has been set to δt MC = 10 −2 τ for all the systems studied, via the acceptance rate A: This result has been applied to rescale the dynamical properties that characterise the structural relaxation of our systems, namely (i) the mean-square displacement; (ii) the self part of the van Hove function; (iii) the selfintermediate scattering function; and (iv ) the dynamic susceptibility.
The mean square displacement (MSD) is the ensemble average of the particles' displacement from their position over a given time window and is defined as where r j is the position vector of particle j and the angular brackets indicate ensemble average over all the HBPs and at least 100 independent trajectories. The self part of the van Hove function (s-VHF) provides the probability distribution of the HBPs' displacements at time t 0 +t, given their position at time t 0 and reads where the symbol δ is the Dirac delta. If the distribution of the displacements over time is Gaussian, then G s (r, t) can be approximated by a Gaussian function that reads where D t is the instantaneous diffusion coefficient and d the dimensionality of the system dynamics. In particular, G ′ s,1 and G ′ s,2 are the Gaussian approximations of the parallel and perpendicular s-VHFs. As far as G ′ s,3 is concerned, we have recently shown that it would not correctly determine the Gaussianity of the displacements in anisotropic systems, and thus proposed a Gaussian distribution that has an ellipsoidal, rather than spherical, symmetry [29]. More specifically, for prolate particles, the Gaussian approximation of the s-VHF function takes the following form: and for oblate particles where D t, and D t,⊥ are, respectively, the instantaneous values of the diffusion coefficients in the direction parallel and perpendicular ton, respectively, F (...) is the Dawson's integral, erf(...) the error function, The self-intermediate scattering function (s-ISF) is useful to gain an insight into the timescales associated to the system's structural relaxation occurring over a length 2π/|k|, where k is a wave vector that corresponds to the location of the relevant peaks in the static structure factor. It is defined as The fluctuations of the self-intermediate scattering function define the so-called dynamical susceptibility, χ 4 (k, t), which provides useful information on the extent of dynamic heterogeneities or, equivalently, on the number of dynamically correlated particles: where f s (k, t) = Np j=1 cos(k[(r j (t + t 0 ) − r j (t 0 ))]/N p is the real part of the instantaneous value of the s-ISF.

III. RESULTS
In this section, we analyse the main features of the equilibrium dynamics of HBPs in N − and N + phases by DMC simulations. All the results shown here have been rescaled with the acceptance rate to recover the time scale of the Brownian dynamics [21][22][23]. We stress that this is an essential step in order to correctly compare the dynamic properties of different systems. The typical behaviour detected in dense colloidal LCs generally consists of a short-time diffusive regime, where the MSD depends linearly on time, followed by an intermediate regime where the particles' motion is slowed down  [27,28].
By contrast, the perpendicular MSD, ∆r 2 ⊥ , reported in the right frame of Fig. 2, is not especially influenced by the particle anisotropy, although it can be observed that, at long time-scales, the MSD of oblate HBPs is slightly larger than that of prolate HBPs. To clarify the role of the shape anisotropy on the particles' ability of diffusing in nematic LCs, we have calculated the diffusion coefficients as D = lim t→∞ ∆r 2 /2dt, where d = 1, 2 or 3 for parallel, perpendicular or total MSD, respectively. The resulting diffusivities are reported in Fig.  3 as a function of W * . It is interesting to notice that at the prolate-oblate crossover, indicated by the vertical dashed line W * = √ L * ≈ 3.46, the parallel diffusivity, D , drops abruptly by a factor of 7, from approximately 6.5×10 −3 D 0 at W * = 3.46 to 10 −3 D 0 at W * = 4. An opposite tendency is detected in the perpendicular diffusivity, D ⊥ , which increases from 2.5 × 10 −3 D 0 at W * = 3.46 to 7×10 −3 D 0 at W * = 4. The combined effect of parallel and perpendicular diffusion results in a total diffusion coefficient, D T OT = (2D ⊥ + D )/3 in Fig. 3, that displays a minimum exactly at the self-dual shape. Despite being the geometry that would more easily stabilise the biaxial nematic phase, the self-dual shape delays diffusion, at least for the range of particle anisotropies explored in this work.
The effect of shape anisotropy is not limited to the single particle ability to diffuse, but extends to the whole system by contributing to determine the degree of Gaussianity of its dynamics as well as the decay of its structural relaxation. More specifically, eventual deviations from Gaussian dynamics are estimated by computing the self-van Hove correlation functions, while the long-time relaxation  sion; t/τ = 3.3 (middle frame) approximately indicates the cross-over from the short-time diffusion to the longtime diffusion, and matches the beginning of the caging effect for the parallel diffusion of oblate HBPs; finally, t/τ = 2400 (right frame) is a representative time within the fully-developed long-time diffusive regime.
The Gaussianity of the parallel and perpendicular s-VHFs has been estimated by comparing the simulation results with the Gaussian distribution of Eq. (6), with d = 1 and 2, respectively. By contrast, the Gaussianity of the total s-VHFs has been assessed with Eq. (7) (prolate HBPs) and (8) (oblate HBPs). The agreement between DMC simulations and theory is excellent, with almost insignificant deviations from a pure Gaussian distribution of displacements detected in systems of prolate HBPs. These tiny deviations, predominantly observed at long distances, suggest the existence of a number of fast particles that are able to displace longer distances than those normally expected. By contrast, at short distances, no significant deviations between simulation and theory are noticed, thus indicating an unambiguous Gaussian behaviour and the absence of particularly slow particles. The s-VHFs of HBPs with 1 < W * < 12 (not shown here) display a similar Gaussian behaviour, the main differences being the shift of the peak towards larger distances and the gradually lower probability of observing fast particles at increasing W * . It should also be noticed that the essentially Gaussian nature of the s-VHFs excludes, at both short and long times, where the MSD is a linear function of time, the occurrence of Fickian yet non-Gaussian (FNG) dynamics, as we had already observed in nematics of uniaxial particles [29]. This finding further confirms that not all soft-matter systems, now also including complex fluids of biaxial particles, necessarily exhibit FNG dynamics.
In the light of these considerations, we now discuss the decay of the density fluctuations in the N + and N − phases by examining the self-intermediate scattering functions, F s (k, t), shown in Fig. 6. The s-ISFs of oblate and prolate HBPs have been calculated by randomly selecting 5 × 10 4 wave vectors k = (2πi 1 /l, 2πi 2 /l, 2πi 3 /l), with i 1 , i 2 and i 3 integer numbers and l the side of the cubic simulation box, whose magnitude is |k| = 2π/T . All curves in Fig. 6 display a single-step decay, which can be closely reproduced by a stretched exponential function of the type exp [−(t/t r ) α )], with t r the relaxation time of the system's structural decay and α the stretching exponent, which ranges between 0.82 at W * = 1 and 0.88 at W * = 12. The former corresponds to the time the s-ISF takes to decays to e −1 and it is reported in the inset of Fig. 6 as a function of the particle reduced width. The latter provides an insight into the nature of  the long-time relaxation dynamics. In particular, α < 1 would indicate the emergence of separated domains of particles with different relaxation dynamics, usually referred to as dynamic heterogeneity. This behaviour, often detected in supercooled liquids, glasses and gels [40][41][42], it is generally attributed to particles that either (i) relax exponentially at different time rates or (ii) relax nonexponentially at the same time rates [43,44].
To assess the origin of dynamic heterogeneity in nematic LCs of HBPs, we calculated the dynamic susceptibility, χ 4 (k, t), which is proportional to the number of dynamically correlated particles over a time t and a distance ≈ 2π/k. This function hence ponders the presence of clusters that, if existed, would be formed by particles relaxing at the same time rate, thus corroborating scenario (ii). The dynamic susceptibility of N + and N − phases of HBPs, reported in Fig. 7, shows a peak at t/τ ≈ 3.3, at the crossover from short-time to long-time diffusive regime. The magnitude of this peak is not significant and, as expected by the relatively large value of the stretching exponent α, which is not much lower than 1, excludes the occurrence of HBPs that are dynamically correlated. In other words, in agreement with uniaxial disks and rods [39], there is no evidence of clusters of particles moving collectively and consequently the structural relaxation of both oblate and prolate nematics of HBPs is entirely determined by the single-particle dynamics.

IV. CONCLUSIONS
In summary, we have studied the equilibrium dynamics of prolate and oblate hard cuboids in nematic LCs. To this end, we have employed the DMC simulation technique, which can accurately reproduce the Brownian dynamics of colloids. In particular, we find Fickian and Gaussian dynamics at both short-time and long-time scales, which are separated by a very smooth crossover indicating an almost completely negligible caging effect. The cuboids' anisotropy (precisely, their width-tothickness ratio) plays a relevant role in the particle's ability to diffuse along the nematic director,n, and perpendicularly to it. Prolate and oblate cuboids diffuse faster in the direction parallel and perpendicular ton, respectively. Interestingly, despite being indicated as the most appropriate geometry for the formation of the elusive biaxial nematic phase, the self-dual shape reduces the particles' total diffusivity, which displays a minimum exactly at W * = √ L * . The analysis of the self-van Hove correlation functions reveals the full Gaussianity of the distribution of displacements over time. Very small deviations are actually observed in N + phases at large distances, suggesting the existence of few fast particles that are not found in N − phases. Finally, the structural relaxation decay, investigated with the self-intermediate scattering functions, exhibits a stretched exponential behaviour with a relaxation time that approximately increases by a factor 2 with increasing particle width-to-thickness ratio from W * = 1 to 12. The stretched exponential nature of the structural relaxation of dense systems, such as supercooled liquids, glasses and gels, usually suggests the occurrence of particles dynamically correlated and eventually displaying collective dynamics [40][41][42]. To assess the emergence of this phenomenon, we calculated the four-point dynamic susceptibility, which measures the magnitude of the dynamic correlations. However, as also noticed for uniaxial particles in nematic LCs [39], no sign of dynamic clusters and collective dynamics has been detected. Their occurrence in positionally ordered LC phases is currently under investigation.