Anomalous Brownian motion of colloidal particle in a nematic environment: effect of the director fluctuations

As recently reported [Turiv T. et al., Science, 2013, Vol. 342, 1351], fluctuations in the orientation of the liquid crystal (LC) director can transfer momentum from the LC to a colloid, such that the diffusion of the colloid becomes anomalous on a short time scale. Using video microscopy and single particle tracking, we investigate random thermal motion of colloidal particles in a nematic liquid crystal for the time scales shorter than the expected time of director fluctuations. At long times, compared to the characteristic time of the nematic director relaxation we observe typical anisotropic Brownian motion with the mean square displacement (MSD) linear in time $\tau$ and inversly proportional to the effective viscosity of the nematic medium. At shorter times, however, the dynamics is markedly nonlinear with MSD growing more slowly (subdiffusion) or faster (superdiffusion) than $\tau$. These results are discussed in the context of coupling of colloidal particle's dynamics to the director fluctuation dynamics.


Introduction
Unlike the visible world around us, the atomic, molecular and small particle worlds are in a state of constant motion. This motion is widely occurring in nature and plays important role in physics, chemistry, biology and engineering [1]. The physical approach to this motion (Brownian motion [2]), developed early in the 20th century by Einstein [3,4], Smoluchowski [5], and Langevin [6], still forms the basis for our understanding of these stochastic dynamics, with the main result, derived by Einstein [3,4], being that the mean squared displacement (MSD) ∆r 2 of a particle undergoing Brownian motion in a Newtonian fluid increases linearly with time, ∆r 2 (τ) = 6Dτ, where D is the diffusion constant. For a spherical particle with hydrodynamic radius R in a fluid with viscosity η, the diffusion coefficient is given by the Stokes-Einstein relation D = k B T /ζ, where k B is the Boltzmann constant, T is the temperature and ζ is the viscous friction coefficient, which under no-slip conditions is given by ζ = 6πRη. The normal, linear diffusion regime for a particle of mass m is established on time scales long compared to the microscopic time m/ζ, which is typically in the nanosecond range. These results are valid for Brownian motion under the influence of two forces, the viscous frictional force linear in particle velocity, F = −ζv, and a random [8,9], in F-actin networks [10], in surfactant formed lyotropic liquid crystal [11] may exhibit a subdiffusive behavior with 0 < α < 1; superdiffusion, α > 1, was observed in concentrated suspensions of swimming bacteria [12,13] and in "living polymers" [14].
The behaviour of colloidal particles in a nematic liquid crystal (NLC) is in many respects more complicated as compared to the isotropic fluid host. First of all, the particle sets a certain director distortion around itself, due to the anisotropic nature of surface interactions (surface anchoring) [15]. These director distortions lead to long-range elastic interactions of the particle with the bounding walls. Second, the orientational nematic order leads to an anisotropy of the Stokes drag [16][17][18][19][20][21][22][23][24][25]. As a result, Brownian motion in a nematic host becomes anisotropic with two different diffusion constants D ∥ and D ⊥ , corresponding to the directions parallel and perpendicular to the nematic director n. Third, although the LC medium is homogeneous, the average axes of orientation fluctuate in time and space and thus influence the diffusive regimes [26,27]. In addition to normal diffusion, the particle experiences two anomalous regimes, with MSD growing slower (subdiffusion) and faster (superdiffusion) than τ. The anomalous diffusion occurs at time scales that correspond to the relaxation times of director fluctuations [26]. All three regimes of diffusion are anisotropic, with the MSD being larger for the motion along the director. Once the nematic is melted, the diffusion becomes normal and isotropic.
In this paper, we elaborate in detail the Brownian dynamics of colloidal particles in a nematic host on time scales shorter than the expected time of director fluctuations τ relax . Following the experiments described in [26], we first give a detailed description of the particle tracking experiments with an extensive analysis of the possible errors in the observed trajectories of the colloids. This is followed by an analysis of the diffusion of colloids, which turns out to be nonlinear with time for a certain time lag, as it is described in [26]. In the theoretical section, we present a model consideration of the diffusion of the colloidal particle in LC, where we explain the anomalous regime by coupling of the particle motion and director fluctuations.

Material
The refractive indices for ordinary (n o ) and extraordinary (n e ) modes of light propagation in a nematic are different from each other. The director n is also the local optic axis. Whenever the director fluctuates, birefringence ∆n = n e − n o and the associated "lens" effect of the distorted optic axis around the particle translate these fluctuations into phantom drifts of the image. To overcome the problem, we use the nematic IS-8200 synthesized at Merck (Germany) with ultra-low birefringence ∆n = 0.0015 (at light wavelength λ = 520 nm) [28]. IS-8200 is a thermotropic (solvent free, single-component) material with the nematic state in the temperature range 47-53°C, above which it melts into the isotropic fluid, figure 1. Small ∆n suppresses the image shifts caused by fluctuations.

Director alignment
Alignment of the director at both the particle's surface and the bounding plates influences the diffusion and thus needs to be controlled. We explored the alignment at the surface of spheres to be perpendicular to the surface. In his case, the particles were functionalized with dimethyl-octadecyl-[3(trimethoxysilyl)propyl] ammonium chloride (DMOAP) [29]. The overall uniform orientation of the nematic was set by two glass plates covered with rubbed polyimide PI-2555 (Nissan Chemicals) alignment layers that produce a uniaxial planar alignment n 0 = (1, 0, 0) = const in the cell. The locally distorted director around the spheres should smoothly transform into n 0 [22]. The resulting equilibrium director configuration is of a dipolar type (with a point defect -hyperbolic hedgehog accompanying the sphere, figure 2) for the normal anchoring. The director distortions around the particle [30] lead to a repulsion from the bounding substrates [31]. The particles levitate [32] in the bulk at some height determined by the balance of gravity and elastic forces.

Time scale
The director field n (r, t ) fluctuates in time and space [32]. The relaxation time of fluctuations is much slower than the relaxation time of vorticity ∼ ρl 2 /η ∼ 1 µs, which is the time needed by the perturbed fluid of density ρ and viscosity η ∼ 0.1 Pa · s to flow over the distance l . For a director perturbation with a wavevector q, this time is τ relax ∼ η eff /(q 2 K ), where η eff and K are the effective viscosity and elastic constant. A uniaxial nematic features five independent viscosities and at least four Frank elastic constants, thus η eff and K are complex combinations of these that depend on the director field [33]. In a flat cell of a finite thickness h, the fluctuation spectrum is restricted; in case of strong director anchoring at the boundaries, the minimum wavevector component in the z-direction perpendicular to the cell boundaries is q z min = π/h [34]. We are interested in wavevectors not much larger than q = π/d , where d is the sphere's diameter, since perturbations with wavelength much shorter than d will produce on average a negligible effect. The corresponding range of viscous relaxation times is thus τ d < τ relax < τ h , where τ h = η eff h 2 /π 2 K and τ d = η eff d 2 /π 2 K . To estimate these, we experimentally determined the quantity η eff /K in the splay Frederiks transition of IS-8200 and found η eff /K ≈ 10 11 s/m 2 , figure 3.

Optical particle tracking: experimental and data analysis procedures
Experimental setup consists of a CMOS high-speed video camera MotionBlitz EoSens mini1 (Microtron GmbH) mounted on an inverted microscope Nikon TE2000-U with a 100 × 1.3 N .A. immersion objective.
The camera is capable of a 5000 fps maximum frame rate (time resolution 0.2 ms); the images are grayscale with 256 gradations (bit depth 8). The pixel size is 14 × 14 microns, so that the magnification is 140 nm/pixel, with the Airy disk covering ∼ 100 pixels. For temperature control, we used a Linkam LTS120 heating stage (accuracy 0.3°C).
For a point light source, its image is diffraction-blurred into an Airy disk, whose intensity profile in the image plane has a root-mean-square (r.m.s.) size of s = 0.44λM /2 N .A., where λ is the wavelength, N .A. the numerical aperture, and M the transverse magnification of the optical system. N .A. should be maximized in order to decrease the diffraction blur. Meanwhile, since the image is pixelised as it is formed on a detector matrix, larger magnification M implies more pixels in the image, and thus more information and better ultimate position accuracy. In practice, N .A. and M are determined by the objective lens, typically with N .A. ∼ 1 and M ∼ 100, so that s ∼ 10 µm.
Particle trajectories are almost invariably analyzed in terms of mean square displacement [35], MSD, i.e., ∆x 2 (τ) and ∆y 2 (τ) , where τ is the time lag and angular brackets stand for ensemble average. In an isotropic medium ∆x 2 (τ) and ∆y 2 (τ) are equal, whereas in an anisotropic medium, such as NLC, they are different. Ensemble average cannot be obtained from a single trace. Instead, one can do a time average which is believed to be the same as the ensemble average in the limit of infinite averaging time. Specifically, one computes where angular brackets now stand for time average. Even though time averages equal the ensemble averages in an ergodic system, time averages over a trace of finite length and time step possess specific statistical errors [36].
The measurement errors (δx i , δy i ) (assumed to have zero mean) in particle coordinates (x i , y i ) yield a positive additive contribution to MSD computed through equation (1), as it is a quadratic form. Assuming where t = i ∆t , one gets for ∆x 2 (τ) (and similarly for ∆y 2 (τ) ) where the first term is the "true" MSD. Assuming that the errors are uncorrelated with the coordinates, 〈xδx〉 = 0, the second term in equation (2) averages to zero. Assuming further that the errors of different trace points are uncorrelated, 〈δx t +τ δx t 〉 = 0, and have the same variance δx 2 = ∆ 2 0 , the last term in equation (2)    so that measurement errors result in a constant additive background in MSD. The assumption of uncorrelated errors almost certainly holds true for static errors, but it may not be the case for dynamic errors. For instance, the errors due to birefringence fluctuations are likely correlated to some extent over the time interval corresponding to the time scales of the fluctuations. This will add a time lag dependence to the background so that, in general, it is a function of τ, i.e., ∆ 2 τ . The presence of the background should be taken into account in data analyses. We will experimentally estimate these errors and their time dependence by determining the apparent MSD of immobilized particles. Digital images of colloidal particles, captured at a maximum frame rate of 2400 fps (time resolution is 0.4 ms), were analyzed to find the coordinates (x, y) of the particle's center using the intensity-weighted algorithms. Specifically, we computed where x i and y i are coordinates of i -th pixel, I i is its intensity. On each frame with 8-bit gray scale (0-255 arbitrary units of intensity) we take into account only the pixels that have an intensity higher than certain threshold intensity, I i ,th . The different threshold intensity results in different apparent MSD for an immobilized particle. Images of the particle at different threshold levels are demonstrated in figure 4. Red area consists of pixels that are taken for the calculation of the particle's center coordinates. The optimal level of the threshold which gives minimum value of the apparent MSD of glued particle was extracted from the dependence of MSD vs intensity threshold, figure 4.

23001-5
To establish the limit of accuracy in the measurements of particle's positions that depends on birefringence, we used particles immobilized (by a Norland adhesive) at the bottom plate of the cell filled with three different fluids: two types of a nematic and water as isotropic fluid. The probing light beam traveled twice through the entire thickness of the cell. The apparent mean square displacement of the immobilized particles vs time lag is shown in figure 5. The apparent displacements represent a cumulative effect of errors in measuring the particle's position caused by the optical system of the microscope, vibrations and birefringence. It grows with birefringence of the material, being the largest for the nematic pentylcyanobiphenyl (5CB) with the highest birefringence (∼ 0.2). In all cases, the apparent MSD in the time range of interest was about 10 −16 m 2 or less; these values are about 100 times smaller than the MSD of free spheres experiencing on these timescales an anomalous diffusion described in the main text.

Results
The measured MSD vs. time lag τ dependencies for d = 5 µm silica spheres in IS-8200 are presented in figure 6 (a), for perpendicular anchoring. In the isotropic phase (elevated temperature T = 60°C) the diffusion is normal in the entire range of time lags, with the diffusion coefficient D = 9.2 · 10 −16 m 2 /s. In the nematic, at T = 50°C, the diffusion becomes anisotropic, with MSD being different when measured parallel and perpendicular to n 0 . At relatively long time scales, τ > (20 − 40) s, both MSD components grow linearly with τ, with the diffusion coefficients D ∥ = 1.9 · 10 −16 m 2 /s and D ⊥ = 1.4 · 10 −16 m 2 /s for the normally anchored spheres. At the times shorter than about (20-40) s, the MSD time dependence becomes markedly nonlinear, figure 6 (a). An apparent MSD of particles glued to the bottom of the cell (see figure 5) is much smaller than the MSD in the nonlinear range and is practically time-independent; it merely adds a constant background to the MSD of the moving particles. Thus, the non-linear behavior at short times is not a spurious effect due to the finite accuracy of measurements [37] or birefringence.
To obtain a better insight into the different diffusion regimes and the characteristic times limiting their borders, we calculated the velocity autocorrelation function C v∥ (τ) = 〈v x (τ)v x (0)〉 [38,39], where v x is the translational velocity of the particle along the x-axis, and a similar quantity C v⊥ (τ) = v y (τ)v y (0) for y-direction. In fact, one calculates the autocorrelation function of the mean velocity over finite time interval between the position of the particle (time lag) which is much shorter than the correlation time. Under this condition, the mean velocity autocorrelation serves as a good estimate of the velocity autocorrelation function [40].
For the diffusion in the isotropic fluid and for the diffusion at large time scales in the nematic, C v∥ (τ) and C v⊥ (τ) are close to zero, figure 6 (b), as it should be for the normal diffusion with MSD growing linearly with τ. However, both C v∥ (τ) and C v⊥ (τ) become negative in the nematic, when the time lag

Discussion
The anomalous Brownian motion is only observed in the nematic phase. As soon as the nematic is brought into isotropic state, the sub-and superdiffusive behavior changes to a normal linear MSD time dependence down to the shortest experimentally accessible times. Thus, anomalous dynamics must be related to the dynamics of additional degrees of freedom that exist in the nematic, but not in the isotropic phase, namely to the nematic director dynamics. The experiments demonstrate that the orientationally ordered environment influences the Brownian motion of a particle most profoundly, causing, in addition to anisotropy, anomalous super-and subdiffusion. The corresponding times, τ sub and τ sup , vary with the type of anchoring at the particle's surface, size, and displacement direction [26]. Above these time scales, the diffusion becomes normal (but still anisotropic). The anomalous character of the diffusion in the range τ < τ sub is evident not only in the nonlinear dependency of MSD on τ, figure 6 (a), but also in the behavior of velocity correlation functions, figure 6 (b).
The current models of Brownian motion in a nematic [16,20,24,25] consider the director field around the particle as being stationary. The predicted diffusion is always normal albeit anisotropic. Our results agree well with these models if τ is large, τ > τ sub . At τ < τ sub , however, the diffusion becomes anomalous; we attribute the effect to the director fluctuations.
The relaxation times of director fluctuations relevant to the Brownian motion are expected to be limited by the interval τ d < τ < τ h , where τ d = ηd 2 /π 2 K and τ h = ηh 2 /π 2 K are the characteristic times associated with the particle diameter and the cell thickness, respectively. At times τ > τ h , the fluctuations are suppressed by the surface anchoring at the cell plates. For τ τ d , the influence of perturbations with the wavelength much shorter than d averages to zero over the distance d . For d = 5 µm in an IS-8200 cell of thickness 50 µm, the limits are estimated as τ d ≈ 0.3 s and τ h ≈ 30 s. The time range τ d < τ < τ h thus embraces the experimentally determined τ sub = (20 − 42) s and τ sup = (4 − 10) s, see figure 6. The nematic is a viscoelastic medium, in which the director field n(r, t ) is coupled to the the velocity field v(r, t ). Both n(r, t ) and v(r, t ) are perturbed by the particle and by the director fluctuations. Translational motion is coupled to the orientational dynamics of n(r, t ). In its turn, director reorientations induce torques and forces that cause the material to flow (the so-called backflow effect [33,41,42]) and thus modify v(r, t ). The director fluctuations establish an intrinsic memory at the scales τ d < τ < τ h , which is typically much longer than the hydrodynamic memory time of the isotropic fluid [32]. The intrinsic memory is known to cause both superdiffusion and subdiffussion, sometimes just by varying the parameters of the very same system [43]. Below we present qualitative effects that help to understand the connection of director fluctuations to the intrinsic viscoelastic memory and anomalous diffusion in The equation of motion for the director fluctuations δn in a bulk nematic liquid crystal reduces to the torque equation for viscous and elastic torques [32], where γ 1 is the rotational viscosity and K is the effective elastic constant. Performing Fourier-Laplace transform of equation (5) yields the dispersion relation with purely imaginary frequency ω = −iK q 2 /γ 1 , from which it follows that the fluctuation modes are overdamped and thus purely relaxational with the relaxation time τ q = γ 1 /(K q 2 ). Their power spectrum is then with the corresponding correlation function Fluctuation amplitude is thus δn q If a nematic is confined in a flat cell, then, due to restrictions imposed by the boundary conditions and particle's size, only a discrete set of fluctuation wave vectors is allowed. Consequently, the relaxation times of different fluctuation modes depend on the size of the inclusion, cell thickness and anchoring conditions at the surfaces [36,44].
A coupling between the LC director and particle dynamics may result in a variety of scenarios of a particle movement. A dipolar inclusion behaves as an elastic dipole with the dipole moment P = aR 2 (a = 2.04, reference [30]) and, therefore, interacts with inhomogeneities of the director field that arise due to thermal fluctuations. The dipole energy is then U = −4πK P ∇n, so that the particle experiences a force F = −∇U = 4πK P ∇ (∇n). Neglecting the inertial effects, the particle then moves with a velocity v that is proportional to the force, Following the Fourier-transformation into reciprocal space, the director fluctuation component δn q with wave vector q results in a particle velocity component Particle velocity autocorrelation function is then For the fluctuation mode that involves splay deformation, δn is in the (n, q) plane and perpendicular to n [32], so that q·δn = qδn sin θ, where θ is the angle between q and n. Particle velocity autocorrelation

23001-8
Anomalous Brownian motion in a nematic medium function is then proportional to the director autocorrelation function of equation (7), where we substitutesd the director correlation function from equation (7) and absorbed all constant prefactors into A = (2P /3ηR) 2 k B T K . The obtained equation describes the contribution of thermal director fluctuations with wavevector q to the particle velocity autocorrelation function. To obtain a full velocity correlation function, the equation should be integrated over q. (Note that fluctuations with different q are uncorrelated.) Integration should only involve fluctuations occurring on length scales large compared to the particle size d = 2R, which corresponds to the wavenumbers smaller than q d = C /d , where C is a constant of the order one. (The high-q fluctuations with the correlation length smaller than the size of the particle will exert uncorrelated forces on different parts of the particle, averaging out to zero.) Thus, where dq is the volume element in q-space.
The MSD is expressed through the velocity autocorrelation function as follows [39,45]: Conversely, C v x (τ) = 1/2 d 2 dτ 2 ∆x 2 (τ) . The coupling mechanism discussed above evidently leads to C v (τ) > 0, and, therefore, to superdiffusion. Thus, direct dipole coupling to the director fluctuation dynamics may be responsible for the superdiffusion that we observe.
If the director field n(r) around the particle is bent by an external action or thermal agitation, it will exert a torque on the particle, proportional to the director rate of changeṅ(t ). In response, the particle will rotate with angular velocity ω ∝ṅ(t ). Neglecting the inertial effects, a particle rotating with an angular velocity ω will also move with a velocity v proportional to |ω|, and thus to |ṅ(t )|, see, e.g., references [25,46]. Thus, the particle velocity is coupled to the director fluctuations of the surrounding nematic, and thus the particle velocity autocorrelation function is proportional to the (q-dependent) director angular velocity correlation function, where c is a coupling constat. One needs to integrate over q, noting as before that only fluctuations with q < q d are relevant. The particle velocity correlation function then becomes Recalling that, for any mechanical property A that is a function on the phase space of a classical manyparticle system, there holds 〈Ȧ(0)Ȧ(t )〉 = − d 2 dt 2 〈A(0)A(t )〉 [45], the director angular velocity correlation function Cṅ ,q (τ) = ṅ −q (0)ṅ q (τ) = − d 2 dτ 2 C n,q (τ), where C n,q (τ) is given by equation (7), so that Clearly, the director angular velocity autocorrelation function is negative. This is easy to understand in view of the fact that director fluctuations are only small angular excursions from the mean, so that if the director rotates in certain direction at a given instant of time, at a later time it should be rotating back, which means a negative angular velocity autocorrelation. At short times, however, the relation equation (16) is not valid, as the initial value of an autocorrelation function must be positive. At short times, the

23001-9
inertial effects cannot be neglected, and correct asymptotic behavior of C n,q (τ) is such that its initial second derivative is negative, so that Cṅ ,q (τ) is positive at short times, as it should be.
From equations (15) and (16) it follows that, if the particle motion couples to the director rotations, then the particle velocity autocorrelation function becomes where B is a constant. Thus, C v (τ) due to this coupling mechanism is negative, and corresponds to subdiffusion.
The coupling mechanism between the director bending modes and particle diffusion, discussed above, is indirect in the sense that director reorientations are directly coupled to particle rotations, and the latter may then couple to translations. Another mechanism that would couple the director bending modes with particle translation is through the backflow effects [47], whereby director fluctuations induce flows that affect the embedded particles. In the simplest realization in two dimensions, the force exerted on the fluid and that generates a backflow depends on the director angular velocityṅ, its gradient ∇ṅ, and the director field gradient [48]. This force generates a viscose flow with a velocity v proportional to the force, and thus depends onṅ and ∇ṅ. Assuming that an embedded particle just follows the flow, its velocity autocorrelation function will thus depend on the director angular velocity autocorrelation function, which leads to subdiffusion, as discussed in the preceding paragraph. This is most readily seen for one of the contributions to the backflow, which is proportional to the projection of the director angular velocity gradient onto the director n, i.e., n∇ṅ [48]. The corresponding contribution to the particle velocity is then v = kn∇ṅ, where k is a proportionality factor. Transforming to the reciprocal space, v q = k n · iq ṅ q = ikq cos θṅ q , (18) where q is the wavevector and θ is the angle between q and n. Velocity autocorrelation function is then Substituting ṅ −q (0)ṅ q (τ) from equation (16) and integrating over q to obtain C v (τ), one ends up with the same integral as in equation (17) (of course with a different pre-factor). Thus, the coupling to the director fluctuations through backflow effects leads to a negative particle velocity autocorrelation function of the form of equation (17), and thus to subdiffusion.
The described mechanisms and their interplay may lead to intricate scenarios of coupled particle dynamics, depending on the time scales and relative strength of the two coupling effects. The time scales may in fact be well separated. Indeed, there are two independent director fluctuation modes that involve different types of director distortions [32]. If the elastic constants pertaining to these deformations are much different, then these two modes will relax on different time scales. In particular, only one of the fluctuation modes involves a splay deformation, which is the only distortion that directly couples to elastic dipoles through the first mechanism discussed above.
The diffusion dynamics shown in figure 6 is a complex and multidimensional process that involves different factors and mechanisms. This complexity does not allow one to take into account all aspects of the diffusion process; however, our analysis gives a physical picture where the elastic interaction between the embedded particles and director fluctuations may describe the complex behavior of particle's dynamics in a NLC. Figure 7 shows the experimental data of IS-8200 for MSD parallel to the nematic director, together with theoretical curves. The data are the same as in figure 6 (a), with the background due to position determination errors (figure 5) subtracted. In order to generate the curves in the figure, we used the viscosity η = 2.5 Pa·s (manufacturer supplied data) and the elastic constants K 11 = K 33 = 30 pN and K 22 = 10 pN (our estimates), from which the average elastic constants are K 1 = 1 3 K 33 + 2 3 K 11 = 30 pN for splay (superdiffusion) and K 2 = 1 3 K 33 + 2 3 K 22 ≈ 17 pN for bend (subdiffusion) [32]. Contributions of the super-and subdiffusion to MSD were computed by integration, equation (13), of the velocity correlation functions of equation (12) and (17), respectively. We assumed that q d = C /d , where C is a constant of the order of one. To obtain the fit in the figure, we varied C and the magnitudes of the super-and subdiffusive contributions. Fit results are seen to agree well with the experimental data, both in the regions of super-(inset in figure 7) and subdiffusion ( figure 7).

23001-10
Anomalous Brownian motion in a nematic medium