Diffusion of a Brownian ellipsoid in a force field

# Diffusion of a Brownian ellipsoid in a force field

## Abstract

We calculate the effective long-term convective velocity and dispersive motion of an ellipsoidal Brownian particle in three dimensions when it is subjected to a constant external force. This long-term motion results as a “net” average behavior from the particle rotation and translation on short time scales. Accordingly, we apply a systematic multi-scale technique to derive the effective equations of motion valid on long times. We verify our theoretical results by comparing them to numerical simulations.

05.40.Jc
05.70.Ln
###### pacs:
05.40.-a

Brownian motion Nonequilibrium and irreversible thermodynamics Fluctuation phenomena, random processes, noise, and Brownian motion

## 1 Introduction

The Brownian motion [1, 2] of small particles suspended in an (aqueous) solvent, driven by the erratic impacts from the solvent molecules, is an ubiquitous phenomenon below micrometer length scales. Its theoretical foundations have been studied for over 100 years, with huge impact on the natural sciences in general and on physics in particular [3, 4]. Still, there are many puzzles of surprisingly fundamental nature which are not yet fully resolved. An important example is the effect of hydrodynamic coupling between the translation and rotation of particles with arbitrary, non-spherical shape, like colloidal particles, colloidal clusters, DNA, proteins, nanotubes etc., on their overall diffusive behavior. This problem has recently attracted considerable interest [5, 6, 7, 8, 9, 10, 11], presumably spurred by the developments in single particle tracking techniques which are able to record both, position and orientation, with high precision [5, 6, 12, 13, 14].

Brownian motion of non-spherical particles is characterized by a crossover from short-term anisotropic diffusion, dominated by the initial particle orientation, to effective “net” diffusion on very long times. For free Brownian motion without external forces, this behavior has been studied in some detail and is quite well understood [5, 7, 8, 9]. The situation is different, however, if the diffusive motion is driven by an external force. The only theoretical studies in that direction, that we are aware of, are a recent analysis of driven Brownian motion of an asymmetrical particle in the plane [10], and a work by Brenner published in 1981 and little noted in the physical literature [15]. Having in mind the settling of a dilute suspension of small ellipsoidal particles due to gravitation, Brenner studied this problem under the term “sedimentation dispersion”. However, for calculating this dispersion, i.e. the effective long-term diffusion coefficient, from the full description of particle translation and orientation, Brenner used what he himself calls an “ad-hoc approach”, and which in fact is not explained in any detail in [15].

In the present paper we analyze the effective long-term motion of a Brownian ellipsoidal particle in three dimensions, using a systematic multi-scale perturbation scheme. This scheme does not require an explicit parametrization of rotations in three dimensions, but can rather be performed on a relatively general and abstract level. We are thus able to fully recover the pioneering results by Brenner using an approach which over the last decades has become more established, and which opens up the prospect for further generalizations. We furthermore compare our results to numerical simulations of the coupled equations of motion for translation and rotation, as far as we know for the first time.

## 2 Model

We model the dynamics of the Brownian ellipsoid by the force and torque balance relations

 0 = −γ˙x+f+√2kBTγ1/2ξ(t), (1a) 0 = −ηω+√2kBTη1/2ζ(t). (1b)

The right-hand sides comprise all (non-inertial) forces and torques which add up to zero, because we neglect inertia effects [16, 17] (overdamped approximation). The term represents the viscous friction force acting on the particle center , given by the particles’ translational velocity multiplying the friction tensor . All the externally applied forces are collected in . Throughout this paper, we consider only constant external force fields , independent of particle position and orientation. The last term in (1a) models thermal fluctuation by unbiased Gaussian white noise sources with correlations . The strength of these fluctuations is characterized by the thermal energy ( is Boltzmann’s constant and the temperature) and the friction tensor. Since is symmetric and positive definite, its square root is uniquely defined, i.e. . The thermal environment is assumed to be homogeneous, so that and are constant in space. In (1b), all quantities represent the rotational counterparts to the ones from (1a), with the (positive definite and symmetric) rotational friction tensor (independent of position ), its square root , the angular velocity , and unbiased Gaussian noise sources which are independent of . We restrict ourselves to the case without externally applied torques so that only viscous friction and thermal fluctuations contribute to the torque balance.

In principle, the dynamics of the Brownian ellipsoid is not fully specified yet by eqs. (1). While the force balance (1a) indeed represents an equation of motion for the translational Brownian movement of the center of the ellipsoid (overdamped Langevin equation [2, 18]), the torque balance (1b) is a kinematic relation for the momentary angular velocity induced by the acting torques. For a full description, it has to be supplemented by a representation of the particle orientation and its equation of motion. As (1) is written in the laboratory frame, the friction tensors and then directly depend on the parameters specifying the particle orientation, i.e. the equations (1) are coupled. Common representations of rotations in three dimensions include Euler angles [19], quaternions [20] and unit vectors (“directors”) attached to the particle [15, 7, 21]. It turns out, surprisingly, that the analysis we are going to present in the following can be performed without choosing a specific representation of orientation, so that the description provided by eqs. (1) is sufficient for our purposes.

As already mentioned in the Introduction, the driven diffusive motion of the ellipsoid described by (1) shows a crossover from short-term anisotropic behavior, dictated by the particle’s initial orientation, to effective long-term convection and diffusion, after the initial orientation has been “forgotten”. This crossover is characterized by the time it takes the particle to diffusively perform a full rotation. We can estimate it as [11]

 τc=1/¯Q, (2)

with and the rotation diffusion tensor. Being the average of the eigenvalues of , the choice of as a “net” rotational diffusion is physically intuitive. The length scale associated with is given by the typical distance the particle covers by convection and diffusion during ,

 lc=max{τc|f|Tr(γ−1)/3,√τc¯D}, (3)

where we again used averaged friction and diffusion coefficients, in particular with the translation diffusion tensor .

On time and length scales and , much larger than and , the ellipsoid will perform an “effective” translational motion, incorporating its rotational diffusion in an averaged way. We can therefore define the small dimensionless parameter

 ε=τc/τL≪1 (4)

to quantify the separation of the small scales and from the large scales and . Note that the long scales do not directly appear in the model (1), but rather are introduced by the question we ask: What is the effective long-term dynamics of the ellipsoid after transients have died out? This is unlike other typical problems involving distinct scales, where these scales are inherent to the problem so that a small (or large) parameter explicitly appears already in the equations of motion [22, 23, 21]. Nevertheless, also in the present case the standard multi-scale or homogenization technique [24, 25, 26] can be applied as a systematic perturbation procedure to derive the sought effective equations.

## 3 Multi-scale analysis

The starting point of the multi-scale analysis is the (forward) Fokker-Planck equation for the probability density that the ellipsoid has a certain position and orientation at time [2, 27, 28]:

 ∂p∂t−(L†+M†)p=0, (5a) with L†=−∂∂xi[(γ−1f)i−∂∂xjDij]. (5b)

Here, we use index notation and the summation convention to sum over repeated indices; are the components of the diffusion tensor . The operator in (5a) represents the generator of rotary diffusion and thus involves the Laplace-operator in rotation space whose specific form depends on the choice of parametrization of orientation. It is clear that this Fokker-Planck equation resolves the probability density on all scales. The essential step to disentangle small and large scales is to explicitly introduce them as independent variables and to presume that is a function of all these variables. Using the symbol for collecting the three parameters representing orientation (without specifying them in any further detail), we define

 ~x=ε0x,X=ε1x,~α=ε0α, (6a) and θ=ε0t,ϑ=ε1t,τ=ε2t, (6b)

and require that

 p=p(θ,ϑ,τ,~x,X,~α). (7)

With these definitions, is of order one only for very large , and thus represents the large scale, which we are interested in. Likewise, the time variables and become of order one at large times , where we expect from their relative scaling with respect to that on scales long-term convective motion occurs and on scales long-term diffusive motion. The small scale variables and essentially correspond to the original variables and , but are restricted to small scales by imposing periodic boundary conditions for , eq. (7), in and . The spatial periodicity is assumed to be , whereas the small rotational scales are obviously periodic by definition. Relevant rotational motion occurs on small scales only. Accordingly, there is no large scale rotational variable defined in (6a), and the long-term effective equations of motion will not include rotational degrees of freedom explicitly.

As a consequence of (6) and (7), the time and spatial derivatives in (5) turn into

 ∂∂t=∂∂θ+ε∂∂ϑ+ε2∂∂τ,∂∂xi=∂∂~xi+ε∂∂Xi, (8)

while the generator of rotational diffusion remains unchanged, in particular it does not involve any terms containing . In view of the scale separation (4) we can treat as a small perturbative parameter and expand in powers of ,

 p=p(0)+εp(1)+ε2p(2)+…, (9)

where all a priori inherit the functional dependence (7) on the various variables. In these variables, is normalized to one, while all other with are normalized to zero. Plugging (8) and (9) into (5), and collecting terms of equal powers in in the resulting expression, we obtain a hierarchy of inhomogeneous Fokker-Planck like equations of which we list the first three (order , and ): {widetext}

 ∂p(0)∂θ−(~L†+~M†)p(0) = 0 (10a) ∂p(1)∂θ−(~L†+~M†)p(1) = −∂p(0)∂ϑ−∂∂Xivip(0)+2∂∂~xi∂∂XjDijp(0) (10b) ∂p(2)∂θ−(~L†+~M†)p(2) = −∂p(0)∂τ−∂p(1)∂ϑ−∂∂Xivip(1)+∂∂Xi∂∂XjDijp(0)+2∂∂~xi∂∂XjDijp(1) (10c)
{floatequation}

see eqs. (10). In (10), we introduce the velocity vector and the tilde over the operators and to indicate that they act on the small scale variables and , respectively. Note that (10a) is exactly the same equation as (5a), with the essential difference, however, that obeys periodic boundary conditions in the variables and .

For finding the solutions of the equation hierarchy (10) we largely follow the standard procedure detailed, e.g., in [26]. We are interested in solutions of (10) which are stationary on small scales after short-term transients have died out. Hence, the desired solutions do not depend on such that we can set for all . A further important observation is that and do not depend on the large scale variable . The solution to (10a) is therefore given by a product ansatz

 p(0)=w(~x,~α)ρ(0)(ϑ,τ,X), (11)

where has to be normalized over and , and over . Exploiting that in is constant in space and that dependencies on particle orientation enter only via (likewise in where the rotational diffusion coefficient depends on orientation via ), we find to be uniform for all and , with a constant value set by normalization: . This uniform distribution carries a probability current (see (5b)), corresponding to an averaged particle velocity

 V=∫d~xd~αvw=¯¯¯¯¯¯¯¯¯¯¯¯γ−1f=¯¯¯¯¯¯¯¯γ−1f, (12)

where the overbar denotes the average over the uniform orientational distribution.

Solving eqs. (10b) and (10c) is a little more involved because of the inhomogeneities on the right-hand sides. For a non-trivial solution to exist, they have to fulfill a so-called solvability condition [26], stating that these inhomogeneities have to be orthogonal to the null-space of the operator adjoint to . The nullspace of contains all constants (in and ), so that the solvability condition for (10b) reads . Inserting our above result (11) for , we find

 ∂ρ(0)∂ϑ+∂∂XiViρ(0)=0. (13)

With this equation (and using ) we can simplify (10b) to,

 (~M†+~L†)p(1)=w(vi−Vi)∂ρ(0)∂Xi. (14)

As is a function of the small scale variables while depends on large scales only, we can solve (14) again by a product ansatz. We thus set

 p(1)=wλij(~x,~α)fj∂ρ(0)∂Xi(ϑ,τ,X), (15)

where the factor has been introduced for later convenience. Likewise, also the specific expression for the unknown vector function of the small scales involving an auxiliary tensor proves to be very convenient.

The ansatz (15) solves (14) provided that the functions solve the auxiliary equation , where we used the definition and eq. (12) for rewriting the right-hand side inhomogeneity in the last step. Observing that this inhomogeneity, as well as all and appearing in , are independent of , we conclude that the do not depend on either, so that . As is supposed to be a solution for any (constant) , the above auxiliary equation thus simplifies to the tensor equation

 ~M†λij=(γ−1)ij−(¯¯¯¯¯¯¯¯γ−1)ij. (16)

In addition, all have to fulfill the “normalization”

 ∫d~αλij=¯¯¯¯¯¯λij=0 (17)

to guarantee overall normalization of in (7).

Analogously to the above procedure for (10b), we next analyse the solvability condition for the second order equation (10c). Plugging in the results we derived so far, namely (11), (13) and (15), we obtain in a straightforward way

 ∂ρ(0)∂τ−Deffij∂∂Xi∂∂Xjρ(0)=0, (18)

with the effective diffusion coefficient on large scales

 Deffij=kBT(¯¯¯¯¯¯¯¯γ−1)ij−¯(γ−1)ikλjlfkfl, (19)

where we used (17) to arrive at the given form of the second term.

As anticipated when introducing the scaling ansatz (6), the main results (13) and (18) of the multi-scale analysis consist in an equation, which describes convective motion with an effective velocity on time scales and large spatial scales (eq. (13)), and an equation, which governs diffusion with an effective diffusion coefficient on time scales and large spatial scales (eq. (18)). In order to combine them into one effective equation of motion, we switch back to the original variables and , and use that the marginal density in translation space only is obtained by integrating over orientational degrees of freedom, . From (7) we conclude that in lowest order . Using (8) (and ) we find in lowest order. Plugging in our results (13) and (18), and noticing that (see (6)) we finally arrive at

 ∂ρ∂t+∂∂xi(Vi−Deffij∂∂xj)ρ=0. (20)

Although it is given in the original variables and , we know from the way it has been derived that this effective forward Fokker-Planck equation is a valid description of the dynamics of our system (1) only in the long-term regime and . As usual in multi-scale schemes [26], the explicit expressions for the “effective coefficients” and are obtained from solving an auxiliary equation on the small scales, here eq. (16), and from averaging over these small scales (see eqs. (12) and (19)).

## 4 Calculation of the effective coefficients

The relevant averages over small scale rotational variables in (12), (16) and (19) are of the form and . Since these averages are performed over a uniform distribution of orientations, the resulting tensors should be invariant under rotation. We can thus use invariance theory to show that

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯(γ−1)ij = 13Tr(γ−1)δij, (21a) ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯(γ−1)ikλjl = [−115δikδjl+110(δijδkl+δilδjk)] (21b) ×Tr(γ−1λ).

In analogy to (21a) we also find , so that has to be traceless according to (17). The property has already been used to simplify (21b).

In order to find the solution to the auxiliary equation (16), we first rewrite its right-hand inhomogeneity using (21a),

 (^γ−1)ij=(γ−1)ij−13Tr(γ−1)δij, (22)

i.e.  denotes the traceless part of the inverse friction tensor. This tensor depends on the particle orientation via rotation matrices , i.e. , where the traceless inverse friction tensor represents an arbitrary reference configuration and thus is constant. Typically, it is chosen such that the principal axes of the ellipsoid are aligned with the coordinate axes of the laboratory frame, because then is diagonal. Assuming that the mobility tensor and the rotational diffusion tensor appearing in can be diagonalized in the same frame, this is also the configuration of the ellipsoid in which (with components ) is diagonal. We may therefore expect that even “inherits” that property such that it can be written as with being diagonal and constant.

To find the explicit form of the rotation matrices and, in particular, the rotary diffusion operator for a specific parametrization of orientation may be quite cumbersome. However, in [29] Rallison showed that we can write

 M†=ϵiklRkm∂∂RlmQijϵjpqRpn∂∂Rqn, (23)

independent of the specific representation of rotation, as only derivatives with respect to the -component of the rotation matrix appear ( is the Levi-Civita tensor). With this expression and the above ansatz for and , (16) reduces to an algebraic equation, which is easily inverted. After some lengthy but straightforward algebra, we find the solution

 λij=16Δ[Tr(Q^γ−1)δij−3(Q^γ−1)ij], (24a) with Δ=Q(1)Q(2)+Q(2)Q(3)+Q(3)Q(1), (24b)

and the eigenvalues of the rotary diffusion tensor .

Plugging the expressions (21) and (24a) into (12) and (19), we obtain the final results for the effective coefficients,

 Vi = 13Tr(γ−1)fi, (25a) Deffij = ¯Dδij+κij, (25b) where we re-used the abbreviation ¯D=Tr(D)/3=kBTTr(γ−1)/3 in (25b) and defined κij=12ΔTr(^γ−1Q^γ−1)(130fifj+110f2δij). (25c)

For the specific choice , the effective diffusion tensor becomes diagonal with

 κ11=115ΔTr(^γ−1Q^γ−1)f2,κ22=κ33=34κ11. (26)

The -component quantifies effective diffusion parallel, and the - and -component perpendicular to the direction of the external force .

## 5 Numerical simulations

We compare our main results (25) to numerical simulations of the original equations of motion (1) for a representative ellipsoid with ratio between its three semi-axes. For the simulations, we choose quaternions as a concrete parametrization of rotation, i.e. the model (1) is supplemented by the equation of motion for the quaternion ,

 ˙q=12ω∘q, (27)

where the symbol denotes a quaternion product [20] evaluated in the Stratonovich sense [2, 27, 28], and where is the angular velocity in the laboratory frame from (1b), represented as a pure quaternion [20]. We solve (1) and (27) using the Euler algorithm with time-step 0.1 . The explicit values for the (eigenvalues of the) friction tensors and of the ellipsoid are calculated from the exact expressions given in [30]. For the simulations we choose an ellipsoid with semi-axis lengths 0.3 μ, 0.6 μ and 0.9 μ. The results of the simulations averaged over 10000 realizations of the noise sources (with identical initial conditions) are shown in Fig. 1. The crossover from short-term to long-term diffusion is clearly visible to occur at a time around the order of a second, well comparable to the estimation from (2). The long-term coefficients obtained from these simulations are in perfect agreement with the theoretical predictions (25).

## 6 Discussion and Conclusions

The effective diffusion coefficient (25b) is composed of two contributions. The first one, , is the purely thermal diffusion of the center of the ellipsoid, averaged over all particle orientations. The second one, , is a dispersion effect stemming from the variations of the translational velocity with diffusive changes of the particle orientation relative to the direction of the external force . This contribution is anisotropic, the dispersion effect is stronger in the direction of the force than perpendicular to it. Remarkably, the ratio between the parallel and perpendicular component (see (26)) is completely independent of the specific shape of the ellipsoid. Furthermore, this anisotropy is only present in three dimensions. For two dimension, the effective diffusion tensor analogous to (25c) turnes out to be isotropic; it has been obtained in [10] by formally solving the Langevin equations of motion for translation and rotation in the plane.

It is straightforward to verify that (25b) and (26) agree with the results obtained by Brenner in [15]. However, unlike [15] we have here shown that (20) and (25) follow from a rigorous perturbative multi-scale technique which has been widely used in physics and applied mathematics over the last 10-20 years. Mean convective velocity and effective diffusion tensor are both computed by solving an auxiliary equation (see eq. (16)), and by performing averages over the stationary rotational distribution of the original model (see eqs. (21)). In the case studied here, where there is no torque, this procedure is pretty straightforward: the stationary rotational distribution is uniform, and the solution of the auxiliary equation can be obtained by a simple ansatz. In more general settings, the stationary distribution as well as the solution of the auxiliary equation would have to be computed numerically, in analogy to the description of scalar transport in compressible flows presented in [24]; we will leave this for future work. Specifically, such more general settings may include external torques [31, 32], or other types of external “forces”, like hydrodynamic flows or phoretic mechanisms which drive particle motion. We finally remark that experiments for measuring the long-term convective and diffusive motion of ellipsoidal particles in three dimensions, and for verifying our theoretical predictions in (25), could be performed along the lines of [12].

###### Acknowledgements.
Financial support by the Swedish Science Council (Vetenskapsrådet) under the grants 621-2012-2982 and 621-2013-3956 is acknowledged.

### References

1. \NameDuplantier B. in Einstein, 1905-2005: Poincaré Seminar 2005, Progress in Mathematical Physics Vol. 47 (Birkhäuser Verlag, Basel, 2005), pp. 201-293.
2. \NameMazo R. M. \BookBrownian Motion: Fluctuations, Dynamics and Applications \PublOxford University Press, Oxford \Year2002.
3. \NameFrey E. Kroy K. \REVIEWAnn. Phys. (Leipzig)14200520.
4. \NameHänggi P. Marchesoni F. \REVIEWChaos152005026101.
5. \NameHan Y. Alsayed A. M. Nobili M. Zhang J. Lubensky T. C. Yodh A. G. \REVIEWScience3142006626.
6. \NameHan Y. Alsayed A. Nobili M. Yodh A. G. \REVIEWPhys. Rev. E802009011403.
7. \NameGonzalez O. Li J. \REVIEWSIAM J. Appl. Math.7020102627.
8. \NameCichocki B. Ekiel-Jezewska M. L. Wajnryb E. \REVIEWJ. Chem. Phys.1362012071102.
9. \NameCichocki B. Ekiel-Jezewska M. L. Wajnryb E. \REVIEWJ. Chem. Phys.1422015214902.
10. \NameGrima R. Yaliraki S. N. \REVIEWJ. Chem. Phys.1272007084511.
11. \NameRibrault C. Triller A. Sekimoto K. \REVIEWPhys. Rev. E752007021112.
12. \NameFung J. Manoharan V. N. \REVIEWPhys. Rev. E882013020302(R).
13. \NameChakrabarty A. Konya A. Wang F. Jonathan V. Selinger J. V. Sun K. Wei Q.-H. \REVIEWPhys. Rev. Lett.1112013160603.
14. \NameChakrabarty A. Konya A. Wang F. Jonathan V. Selinger J. V. Sun K. Wei Q.-H. \REVIEWLangmuir30201413844.
15. \NameBrenner H. \REVIEWJ. Colloid. Interface Sci.801981548.
16. \NamePurcell E. M. \REVIEWAm. J. Phys.4519773.
17. \NameDusenbery D. B. \BookLiving at Micro Scale \PublHarvard University Press, Cambridge, MA \Year2011.
18. \NameSnook I. \BookThe Langevin and Generalised Langevin Approach to the Dynamics of Atomic, Polymeric and Colloidal Systems \PublElsevier, Amsterdam \Year2007.
19. \NameGoldstein H. \BookClassical Mechanics \PublAddison-Wesley, New York \Year1980.
20. \NameCoutsias E. A. Romero L. The Quaternions with Applications to Rigid Body Dynamics; Sandia National Laboratories, Technical Report SAND2004-0153 (2004).
21. \NameMarino R. Eichhorn R. Aurell E. \REVIEWPhys. Rev. E932016012132.
22. \NamePagitsas M. Nadim A. Brenner H. \REVIEWPhysica135A1986533.
23. \NameMazzino A. Musacchio S. Vulpiani A. \REVIEWPhys. Rev. E712005011113.
24. \NameVergassola M. Avellaneda M. \REVIEWPhysica D1061997148.
25. \NameBender C. Orszag S. A. \BookAdvanced Mathematical Methods For Scientists and Engineers: Asymptotic Methods and Perturbation Theory \PublSpringer, New York \Year1999.
26. \NamePavliotis G. A. Stuart A. M. \BookMultiscale Methods: Averaging and Homogenization \PublSpringer, New York \Year2008.
27. \NameGardiner C. W. \BookHandbook of Stochastic Methods \PublSpringer, Berlin \Year1983.
28. \Namevan Kampen N. G. \BookStochastic Processes in Physics and Chemistry \PublNorth Holland, Amsterdam \Year1987.
29. \NameRallison J. M. \REVIEWJ. Fluid Mech.841978237.
30. \NameBrenner H. \REVIEWJ. Colloid Interface Sci.231967407.
31. \NameGüell O. Tierno P. Sagués F. \REVIEWEur. Phys. J. Special Topics187201015.
32. \NameSandoval M. \REVIEWPhys. Rev. E872013032708.
104905