2.1. Hydrodynamics of an axisymmetric particle

We consider a general axisymmetric rigid particle suspended in a quiescent Newtonian fluid and undergoing Brownian motion, as illustrated in figure 1. Under Stokes flow conditions ( $Re \to 0$ ), hydrodynamic forces and torques depend linearly on the translational and angular velocities of a particle. An axisymmetric particle possesses a unique point where translational and rotational motions are decoupled (Kim & Karrila Reference Kim and Karrila1991): a force applied at this point (regardless of direction) does not produce rotational motion. This unique point, known as the hydrodynamic centre $\boldsymbol{R}_{{h}}$ , necessarily lies on the axis of symmetry. As long as the surface geometry of the particle is fixed, $\boldsymbol{R}_{{h}}$ does not vary along the axis and specifies the particle’s position. The frictional forces on particles arise from the viscosity of the surrounding fluid. Referring to $\boldsymbol{R}_{{h}}$ , the hydrodynamic force $\boldsymbol{F}_{{h}}$ and torque $\boldsymbol{T}_{{h}}$ acting on an axisymmetric particle can be expressed as

(2.1) \begin{equation} \begin{bmatrix} \boldsymbol{F}_{{h}} \\ \boldsymbol{T}_{{h}} \end{bmatrix} =- \begin{bmatrix} \boldsymbol{\zeta }_{{t}} & \boldsymbol{0} \\ \boldsymbol{0} & \boldsymbol{\zeta }_{{r}} \end{bmatrix} \boldsymbol{\cdot }\begin{bmatrix} \boldsymbol{u} \\ \boldsymbol{\omega } \end{bmatrix}, \end{equation}

where $\boldsymbol{u}$ and $\boldsymbol{\omega }$ represent the translational velocity of the hydrodynamic centre and the angular velocity of the particle, respectively. In terms of time derivatives, we have $\dot {\boldsymbol{R}}_{{h}}=\boldsymbol{u}$ and $\dot {\boldsymbol{n}}=\boldsymbol{\omega }\times \boldsymbol{n}$ , where $\boldsymbol{n}$ is the unit vector along the symmetry axis. Because hydrodynamic force depends on particle orientation, the translational ( $\boldsymbol{\zeta }_{{t}}$ ) and rotational ( $\boldsymbol{\zeta }_{{r}}$ ) resistance matrices are expressed through parallel and transverse components relative to the symmetry axis, respectively, i.e.

(2.2) \begin{equation} \boldsymbol{\zeta }_{{t}}=\zeta _{{t}}^\parallel \boldsymbol{n}\boldsymbol{n}+\zeta _{{t}}^\perp (\boldsymbol{\delta }-\boldsymbol{n}\boldsymbol{n}), \qquad \boldsymbol{\zeta }_{{r}}=\zeta _{{r}}^\parallel \boldsymbol{n}\boldsymbol{n}+\zeta _{{r}}^\perp (\boldsymbol{\delta }-\boldsymbol{n}\boldsymbol{n}). \end{equation}

The four scalar resistance coefficients, $\zeta _{{t}}^\parallel$ , $\zeta _{{t}}^\perp$ , $\zeta _{{r}}^\parallel$ and $\zeta _{{r}}^\perp$ , depend on the particle geometry (i.e. its size and shape) and on the fluid viscosity. Explicit expressions for the resistance coefficients for prolate and oblate spheroids are tabulated in Appendix D.

Figure 1.

Illustration of an axisymmetric rigid particle, with $\boldsymbol{R}_{{h}}$ and $\boldsymbol{R}_{{c}}$ denoting the hydrodynamic centre and force centre, respectively. Vector $\boldsymbol{n}$ is the unit vector along the axis of symmetry, and $l_{{c}}$ represents the centre offset. The force acting at the force centre is the sum of gravity ( $M\boldsymbol{g}$ ) and buoyancy ( $-M_{{b}}\boldsymbol{g}$ ).

Figure 1. Long description

A three-dimensional diagram of an axisymmetric rigid particle. The diagram includes labeled centers: the hydrodynamic center denoted as R_h and the force center denoted as R_c. The unit vector n is shown along the axis of symmetry, and l_c represents the center offset. The force acting at the force center is depicted as the sum of gravity (M - M_b)g and buoyancy. The diagram illustrates the relationship between these components and the orientation of the particle.

Before proceeding, let us consider the properties of gravity and buoyancy for an axisymmetric particle. Sedimentation arises from density difference between particle and fluid. Gravitational forces distributed over the particle are mechanically equivalent to a force acting at the centre of mass $\boldsymbol{R}_{{m}}$ , given by

(2.3) \begin{equation} \boldsymbol{R}_{{m}} = \frac {1}{M} \int _{V_{{p}}} {\rm d}\boldsymbol{r} \rho (\boldsymbol{r}) \boldsymbol{r}, \end{equation}

where $\rho (\boldsymbol{r})$ is the local density at position $\boldsymbol{r}$ within the particle, $M=\int _{V_{{p}}} {\rm d} \boldsymbol{r} \rho (\boldsymbol{r})$ is the total mass of the particle and $ V_{{p}}$ denotes the particle volume region. Similarly, the buoyancy forces are represented by a force acting at the centre of buoyancy $\boldsymbol{R}_{{b}}$ , given by

(2.4) \begin{equation} \boldsymbol{R}_{{b}} = \frac {1}{M_{{b}}} \int _{V_{{p}}} {\rm d}\boldsymbol{r} \rho _0 \boldsymbol{r}, \end{equation}

where $\rho _0$ is the density of the surrounding fluid and $M_{{b}}=\rho _0 \int _{V_{{p}}} {\rm d}\boldsymbol{r}$ is the total mass of the displaced fluid. Centres $\boldsymbol{R}_{{m}}$ and $\boldsymbol{R}_{{b}}$ must lie on the axis due to axisymmetry, but they are generally different because they depend on different mass distributions. Nevertheless, the two gravitational forces can be combined further. The total potential energy of the system is written as

(2.5) \begin{equation} U= -\left (M \boldsymbol{g}\boldsymbol{\cdot }\boldsymbol{R}_{{m}} - M_{{b}} \boldsymbol{g} \boldsymbol{\cdot }\boldsymbol{R}_{{b}} \right ) = -(M-M_{{b}}) \boldsymbol{g}\boldsymbol{\cdot }\boldsymbol{R}_{{c}}. \end{equation}

Substituting (2.3) and (2.4) into (2.5), we can combine the effects of gravity and buoyancy acting on the particle, as shown in the second equality, which is equivalently represented by the force $(M-M_{{b}}) \boldsymbol{g}$ exerted on a point referred to as the sedimentation force centre (or simply the force centre). Even though $\boldsymbol{R}_{{m}}$ and $\boldsymbol{R}_{{b}}$ are different, the position of the force centre is given by

(2.6) \begin{equation} \boldsymbol{R}_{{c}} = \frac {1}{M-M_{{b}}} \int _{V_{{p}}} {\rm d}\boldsymbol{r} [\rho (\boldsymbol{r})-\rho _0] \boldsymbol{r}. \end{equation}

The force centre is defined only for particles with $M \neq M_{{b}}$ . In the following, we only consider such particles, since when $M=M_{{b}}$ , the particle does not sediment and therefore shows no Taylor dispersion.

The explicit expressions for the dispersion coefficients have been obtained for the case where the sedimentation force centre coincides with the hydrodynamic centre, $\boldsymbol{R}_{{c}} = \boldsymbol{R}_{{h}}$ . As illustrated in figure 1, $\boldsymbol{R}_{{c}}$ generally differs from $\boldsymbol{R}_{{h}}$ , because $\boldsymbol{R}_{{h}}$ is exclusively determined by surface geometry, whereas $\boldsymbol{R}_{{c}}$ depends on the mass distribution within the particle. For an axisymmetric particle with $\boldsymbol{R}_{{c}} \neq \boldsymbol{R}_{{h}}$ , one has

(2.7) \begin{equation} \boldsymbol{R}_{{c}}=\boldsymbol{R}_{{h}}+l_{{c}}\boldsymbol{n}, \end{equation}

since the two centres are on the same axis, where $l_{{c}}$ is the centre offset. One can also conduct the analyses with the separate offsets of the centres of mass $l_{{m}}$ and buoyancy $l_{{b}}$ from $\boldsymbol{R}_{{h}}$ (defined by $\boldsymbol{R}_{{m}}=\boldsymbol{R}_{{h}}+l_{{m}}\boldsymbol{n}$ and $\boldsymbol{R}_{{b}}=\boldsymbol{R}_{{h}}+l_{{b}}\boldsymbol{n}$ ). The offsets $l_{{m}}$ and $l_{{b}}$ are all rolled up into $l_{{c}}$ , which is expressed as

(2.8) \begin{equation} l_{{c}} = \frac {l_{{m}} M - l_{{b}} M_{{b}}}{M - M_{{b}}}. \end{equation}

Note that $l_{{c}}$ is zero even if $l_{{m}}$ and $l_{{b}}$ are non-zero when $l_{{m}} M = l_{{b}} M_{{b}}$ . Adjustments to $l_{{c}}$ can be achieved by modifying either the particle’s asymmetric shape or its non-uniform mass distribution. While changes in shape also affect the friction coefficients, changes in mass distribution do not.

2.2. Smoluchowski equation under gravity

The configuration of the particle is entirely determined by independent variables $\boldsymbol{R}_{{h}}$ and $\boldsymbol{n}$ . Due to the Brownian motion, $\boldsymbol{R}_{{h}}$ and $\boldsymbol{n}$ evolve in a stochastic manner. Let $\psi ( \boldsymbol{R}_{{h}},\boldsymbol{n},t)$ denote the probability density function of finding the particle at position $\boldsymbol{R}_{{h}}$ and orientation $\boldsymbol{n}$ at time $t$ . The particle dynamics is governed by the hydrodynamic force (2.1) and gravitational force (gradient of (2.5)). Following Onsager’s variational principle (Appendix A), the probability density $\psi$ evolves according to the Smoluchowski equation:

(2.9) \begin{align} \frac {\partial \psi }{\partial t} &= D_{{r}} \frac {\partial }{\partial \boldsymbol{n}}\boldsymbol{\cdot }\left (\boldsymbol{\delta }-\boldsymbol{n}\boldsymbol{n}\right ) \boldsymbol{\cdot }\left (\frac {\partial \psi }{\partial \boldsymbol{n}} - \frac {(M-M_{{b}})l_{{c}}\boldsymbol{g}}{k_{{B}}\kern-1pt T}\psi \right ) \nonumber \\&\quad + D^\perp \frac {\partial }{\partial \boldsymbol{R}_{{h}}}\boldsymbol{\cdot }\left (\boldsymbol{\delta } + \frac {\zeta _{{t}}^\perp -\zeta _{{t}}^\parallel }{\zeta _{{t}}^\parallel } \boldsymbol{n}\boldsymbol{n}\right )\boldsymbol{\cdot }\left (\frac {\partial \psi }{\partial \boldsymbol{R}_{{h}}}-\frac {(M-M_{{b}})\boldsymbol{g}}{k_{{B}}\kern-1pt T}\psi \right )\!, \end{align}

where $D_{{r}} = k_{{B}}\kern-1pt T / \zeta _{{r}}^\perp$ is the rotational diffusion coefficient, $D^\perp = k_{{B}}\kern-1pt T / \zeta _{{t}}^\perp$ is the transverse diffusion coefficient, $k_{{B}}$ is the Boltzmann constant and $T$ is the temperature. We adopt the rotational relaxation time $\tau _{{r}}$ as the reference time scale:

(2.10) \begin{equation} \tau _{{r}} = \frac {1}{D_{{r}}} = \frac {\zeta _{{r}}^\perp }{ k_{{B}}\kern-1pt T}. \end{equation}

The rotational relaxation time arises from the rotational diffusion in the absence of gravitational torque. It is the time scale beyond which correlations in the rod orientation become negligible. In the presence of gravity, the magnitude of the gravitational force and the maximum magnitude of the gravitational torque exerted on particles are given by $F_{{g}}=(M - M_{{b}})g$ and $T_{{g}}=(M - M_{{b}})gl_{{c}}$ , respectively, where $g=|\boldsymbol{g}|$ denotes the magnitude of the gravitational acceleration. The corresponding characteristic magnitudes of the translational velocity and angular velocities are $u=F_{{g}}/\zeta _{{t}}^{\perp }$ and $\omega =T_{{g}}/\zeta _{{r}}^{\perp }$ , respectively. Therefore, the system has two other characteristic time scales induced by gravitational effects. One is the reorientation time defined by

(2.11) \begin{equation} \tau _{{o}} = \frac {1}{\omega } = \frac {\zeta _{{r}}^{\perp }}{ (M-M_{{b}})gl_{{c}} }, \end{equation}

which corresponds to the time scale associated with rotation caused by the gravitational torque in the absence of Brownian motion. The other is the sedimentation time, defined by

(2.12) \begin{equation} \tau _{{s}} = \frac {L}{u} = \frac {L\zeta _{{t}}^{\perp }}{ (M-M_{{b}}) g}, \end{equation}

where $L$ is the length scale of the particle, defined by its volume $V$ as $L=[3V/(4\pi )]^{1/3}$ . It corresponds to the time scale for translation due to the gravitational force in the absence of Brownian motion. The ratios of rotational relaxation time to gravity-induced time scales, $\tau _{{r}} / \tau _{{o}}$ and $\tau _{{r}} / \tau _{{s}}$ , are two key dimensionless parameters in this system. They are referred to as the reorientation Péclet number $\alpha$ (also referred to as the dimensionless Langevin parameter in Brenner (Reference Brenner1979)) and the sedimentation Péclet number $\beta$ , respectively:

(2.13) \begin{equation} \alpha = \frac {\tau _{{r}}}{\tau _{{o}}} =\frac { (M-M_{{b}})gl_{{c}}}{ k_{{B}}\kern-1pt T}, \qquad \beta = \frac {\tau _{{r}}}{\tau _{{s}}} = \frac {\zeta _{{r}}^{\perp }}{\zeta _{{t}}^{\perp }} \frac {(M-M_{{b}}) g}{L k_{{B}}\kern-1pt T}. \end{equation}

For a very large $\alpha$ , the gravitational torque overwhelms rotational Brownian fluctuations, and particle rotation is dominated by gravitational torque rather than rotational diffusion. Similarly, for a very large $\beta$ , particle translation is dominated by gravitational force.

We adopt $\tau _{{r}}$ and $L$ as the unit of time and length. Introducing dimensionless time $\tilde {t} = t / \tau _{{r}}$ and position $ \tilde {\boldsymbol{R}}_{{h}} = \boldsymbol{R}_{{h}} / L$ , the Smoluchowski equation (2.9) is non-dimensionalised to

(2.14a) \begin{equation} \frac {\partial \psi }{\partial \tilde {t}} = \tilde {\mathcal{L}}\psi , \end{equation}

with

(2.14b) \begin{equation} \tilde {\mathcal{L}} = \frac {\partial }{\partial \boldsymbol{n}}\boldsymbol{\cdot }\left (\boldsymbol{\delta }-\boldsymbol{n}\boldsymbol{n}\right )\boldsymbol{\cdot }\left (\frac {\partial }{\partial \boldsymbol{n}} - \alpha \hat {\boldsymbol{g}}\right ) + \frac {\partial }{\partial \tilde {\boldsymbol{R}}_{{h}}}\boldsymbol{\cdot }\left (\boldsymbol{\delta }+\chi \boldsymbol{n}\boldsymbol{n}\right )\boldsymbol{\cdot }\left (\tilde {D}^\perp \frac {\partial }{\partial \tilde {\boldsymbol{R}}_{{h}}}-\beta \hat {\boldsymbol{g}}\right )\!, \end{equation}

where $\hat {\boldsymbol{g}} = \boldsymbol{g} /g$ , and we have defined

(2.15) \begin{equation} \chi = \frac {\zeta _{{t}}^\perp -\zeta _{{t}}^\parallel }{ \zeta _{{t}}^\parallel }, \end{equation}

which characterises the hydrodynamic anisotropy of the particle. The dimensionless transverse diffusion coefficient $\tilde {D}^\perp$ is given by

(2.16) \begin{equation} \tilde {D}^\perp = \frac {D^\perp \tau _{{r}}}{L^2} = \frac {\zeta _{{r}}^\perp }{\zeta _{{t}}^\perp L^2} . \end{equation}

Using the results in Appendix D, the dimensionless friction-related quantities $\chi$ and $\tilde {D}^\perp$ depend solely on the aspect ratio $r$ .

Let $\tilde {\varOmega }=(\tilde {\boldsymbol{R}}_{{h}},\boldsymbol{n})$ denotes the particle configuration at time $\tilde {t}$ and $\tilde {\varOmega }^\prime=(\tilde {\boldsymbol{R}}_{{h}}^\prime,\boldsymbol{n}^\prime)$ denotes the configuration at $\tilde {t}^\prime$ . Given the initial condition $\left .\psi \right |_{\tilde {t}=\tilde {t}^\prime}=\delta (\tilde {\boldsymbol{R}}_{{h}}-\tilde {\boldsymbol{R}}_{{h}}^\prime)\delta (\boldsymbol{n}-\boldsymbol{n}^\prime)$ , the solution of the Smoluchowski equation (2.14) is the Green’s function, which is denoted by $\mathcal{G}(\tilde {\varOmega },\tilde {t};\tilde {\varOmega }^\prime,\tilde {t}^\prime)$ . The Green’s function represents the conditional probability density of finding the particle in $\tilde {\varOmega }$ at time $\tilde {t}$ , given that it was in $\tilde {\varOmega }^\prime$ at $\tilde {t}^\prime$ . Therefore, the ensemble average of a quantity $\boldsymbol{\mathcal{F}}(\tilde {\varOmega },\tilde {\varOmega }^\prime)$ can be generally obtained using the Green’s function and initial distribution $\psi _{\textit{in}}$ at $\tilde {t}^\prime$ , i.e.

(2.17) \begin{equation} \left\langle \boldsymbol{\mathcal{F}}\left(\tilde {t},\tilde {t}^\prime\right)\right\rangle = \int {\rm d}\tilde {\varOmega } \int {\rm d}\tilde {\varOmega }^\prime\mathcal{G}\left(\tilde {\varOmega },\tilde {t};\tilde {\varOmega }^\prime,\tilde {t}^\prime\right)\psi _{\textit{in}}\left(\tilde {\varOmega }^\prime,\tilde {t}^\prime\right)\boldsymbol{\mathcal{F}}\left(\tilde {\varOmega },\tilde {\varOmega }^\prime\right)\!, \end{equation}

where ${\rm d}\tilde {\varOmega }^\prime={\rm d}\tilde {\boldsymbol{R}}_{{h}}^\prime {\rm d}\boldsymbol{n}^\prime$ is the volume element in the dimensionless configuration space $\tilde {\varOmega }^\prime $ .

Additionally, one can multiply both sides of the Smoluchowski equation (2.14a ) by the function $\boldsymbol{\mathcal{F}}(\tilde {\boldsymbol{R}}_{{h}},\boldsymbol{n},\tilde {\boldsymbol{R}}_{{h}}^\prime,\boldsymbol{n}^\prime)$ and use the definition of the Green’s function to obtain the following evolution equation (Doi & Edwards Reference Doi and Edwards1986; Xiong, Seto & Doi Reference Xiong, Seto and Doi2024):

(2.18a) \begin{equation} \frac {\partial \left\langle \boldsymbol{\mathcal{F}}\left(\tilde {t},\tilde {t}^\prime\right)\right\rangle }{\partial \tilde {t}} = \left\langle \tilde {\mathcal{L}}^\dagger \boldsymbol{\mathcal{F}}\left(\tilde {\boldsymbol{R}}_{{h}},\boldsymbol{n},\tilde {\boldsymbol{R}}_{{h}}^\prime,\boldsymbol{n}^\prime\right)\right\rangle\! , \end{equation}

where $\tilde {\mathcal{L}}^\dagger$ is the conjugate operator of $\tilde {\mathcal{L}}$ defined by

(2.18b) \begin{equation} \tilde {\mathcal{L}}^\dagger = \left (\frac {\partial }{\partial \boldsymbol{n}} + \alpha \hat {\boldsymbol{g}}\right )\boldsymbol{\cdot }\left (\boldsymbol{\delta }-\boldsymbol{n}\boldsymbol{n}\right )\boldsymbol{\cdot }\frac {\partial }{\partial \boldsymbol{n}} + \left (\tilde {D}^\perp \frac {\partial }{\partial \tilde {\boldsymbol{R}}_{{h}}} + \beta \hat {\boldsymbol{g}}\right )\boldsymbol{\cdot }\left (\boldsymbol{\delta }+\chi \boldsymbol{n}\boldsymbol{n}\right ) \boldsymbol{\cdot }\frac {\partial }{\partial \tilde {\boldsymbol{R}}_{{h}}}. \end{equation}

The ensemble average could also be calculated via (2.18) without explicitly knowing the Green’s function.

3.1. Application of the Smoluchowski equation

We now examine the sedimentation behaviour of such a particle in a gravitational field. Here, we set the gravity direction as $\hat {\boldsymbol{g}}=-\boldsymbol{e}_z$ and consider the sedimentation process starting from the origin, $\tilde {\boldsymbol{R}}_{{h}}(0)=\boldsymbol{0}$ . The unit basis vectors along the $x$ , $y$ and $z$ axes are denoted by $\boldsymbol{e}_x$ , $\boldsymbol{e}_y$ and $\boldsymbol{e}_z$ , respectively. Since the anisotropy in this system arises solely from gravity, the sedimentation process with horizontal motion can be statistically characterised by the mean position coordinate along the gravity direction $\langle \tilde {z}_{{h}}(\tilde {t}) \rangle$ , and the MSD in the $x$ – $y$ plane $\langle \tilde {x}_{{h}}^2(\tilde {t}) + \tilde {y}_{{h}}^2(\tilde {t})\rangle$ . They are expressed as follows:

(3.1a) \begin{align}&\quad \tilde {z}_{{h}}\big(\tilde {t}\big) =\boldsymbol{e}_z\boldsymbol{\cdot }\tilde {\boldsymbol{R}}_{{h}}\big(\tilde {t}\big)= -\hat {\boldsymbol{g}}\boldsymbol{\cdot }\left[\tilde {\boldsymbol{R}}_{{h}}\big(\tilde {t}\big)-\tilde {\boldsymbol{R}}_{{h}}(0)\right]\!, \end{align}

(3.1b) \begin{align}& \tilde {x}_{{h}}^2\big(\tilde {t}\big)+\tilde {y}_{{h}}^2\big(\tilde {t}\big) = \left \{\big(\boldsymbol{\delta }-\hat {\boldsymbol{g}}\hat {\boldsymbol{g}}\big)\boldsymbol{\cdot }\left[\tilde {\boldsymbol{R}}_{{h}}\big(\tilde {t}\big)-\tilde {\boldsymbol{R}}_{{h}}(0)\right] \right \}^2. \end{align}

The mean position $\langle \tilde {z}_{{h}}(\tilde {t})\rangle$ (the drift displacement in the $z$ direction) during sedimentation is obtained by applying (2.18a ) to (3.1a ), yielding

(3.2) \begin{gather} \frac {\partial \langle \tilde {z}_{{h}}\rangle }{\partial \tilde {t}} = -\left \langle \tilde {\mathcal{L}}^\dagger \hat {\boldsymbol{g}}\boldsymbol{\cdot }\left[\tilde {\boldsymbol{R}}_{{h}}\bigl(\tilde {t}\bigr)-\tilde {\boldsymbol{R}}_{{h}}(0)\right] \right \rangle = -\beta \big ( 1+\chi \hat {\boldsymbol{g}}\boldsymbol{\cdot }\left \langle \boldsymbol{n}\boldsymbol{n} \right \rangle \boldsymbol{\cdot }\hat {\boldsymbol{g}} \big ). \end{gather}

Similarly, applying (2.18a ) to the square of (3.1a ) yields (see Appendix C for details)

(3.3) \begin{align} \frac {\partial \left\langle \left [\tilde {z}_{{h}}-\langle \tilde {z}_{{h}}\rangle \right ]^2\right\rangle }{\partial \tilde {t}} &= 2\tilde {D}^\perp \big ( 1+\chi \hat {\boldsymbol{g}}\boldsymbol{\cdot }\left \langle \boldsymbol{n}\boldsymbol{n} \right \rangle \boldsymbol{\cdot }\hat {\boldsymbol{g}} \big ) \notag \\ &\quad +2\left (\beta \chi \right )^2\int _0^{\tilde {t}}{\rm d}\tilde {t}^\prime\left [\left \langle \hat {\boldsymbol{g}}\boldsymbol{\cdot }\boldsymbol{n}\boldsymbol{n}\boldsymbol{\cdot }\hat {\boldsymbol{g}}\hat {\boldsymbol{g}}\boldsymbol{\cdot }\boldsymbol{n}^\prime\boldsymbol{n}^\prime \boldsymbol{\cdot }\hat {\boldsymbol{g}}\right \rangle -\hat {\boldsymbol{g}}\boldsymbol{\cdot }\left \langle \boldsymbol{n}\boldsymbol{n}\right \rangle \boldsymbol{\cdot }\hat {\boldsymbol{g}}\hat {\boldsymbol{g}}\boldsymbol{\cdot }\left \langle \boldsymbol{n}^\prime\boldsymbol{n}^\prime \right \rangle \boldsymbol{\cdot }\hat {\boldsymbol{g}}\right ]\!, \end{align}

where $ \langle [\tilde {z}_{{h}}-\langle \tilde {z}_{{h}}\rangle ]^2\rangle$ represents the $z$ -direction MSD with eliminating the drift displacement in the $z$ direction. Here $\boldsymbol{n}$ and $\boldsymbol{n}^\prime$ denote the orientation of the particle at time $\tilde {t}$ and $\tilde {t}^\prime$ , respectively. The MSD in the $x$ – $y$ plane is obtained by applying (2.18a ) to (3.1b ), yielding

(3.4) \begin{align} \frac {\partial \left \langle \tilde {x}_{{h}}^2+\tilde {y}_{{h}}^2 \right \rangle }{\partial \tilde {t}} &= 2\tilde {D}^\perp \left [2+\chi \big (1- \hat {\boldsymbol{g}}\boldsymbol{\cdot }\left \langle \boldsymbol{n}\boldsymbol{n} \right \rangle \boldsymbol{\cdot }\hat {\boldsymbol{g}} \big )\right ] \notag \\ &\quad +2\left (\beta \chi \right )^2\int _0^{\tilde {t}}{\rm d}\tilde {t}^\prime\left \langle \hat {\boldsymbol{g}}\boldsymbol{\cdot }\boldsymbol{n}\boldsymbol{n}\boldsymbol{\cdot }(\boldsymbol{\delta }-\hat {\boldsymbol{g}}\hat {\boldsymbol{g}})\boldsymbol{\cdot }\boldsymbol{n}^\prime\boldsymbol{n}^\prime\boldsymbol{\cdot }\hat {\boldsymbol{g}}\right \rangle . \end{align}

3.2. Orientational distribution under gravitational torque

Note that calculation of the sedimentation velocity (from (3.2)) and diffusion (from (3.3) and (3.4)) requires first evaluating the orientation dynamics, especially calculating quantities such as $\left \langle \boldsymbol{n}\boldsymbol{n} \right \rangle$ , $\left \langle \hat {\boldsymbol{g}}\boldsymbol{\cdot }\boldsymbol{n}\boldsymbol{n}\boldsymbol{\cdot }\hat {\boldsymbol{g}}\hat {\boldsymbol{g}}\boldsymbol{\cdot }\boldsymbol{n}^\prime\boldsymbol{n}^\prime \boldsymbol{\cdot }\hat {\boldsymbol{g}}\right \rangle$ and $\left \langle \hat {\boldsymbol{g}}\boldsymbol{\cdot }\boldsymbol{n}\boldsymbol{n}\boldsymbol{\cdot }(\boldsymbol{\delta }-\hat {\boldsymbol{g}}\hat {\boldsymbol{g}})\boldsymbol{\cdot }\boldsymbol{n}^\prime\boldsymbol{n}^\prime\boldsymbol{\cdot }\hat {\boldsymbol{g}}\right \rangle$ . Since these ensemble averages depend solely on orientation, integration over $\tilde {\boldsymbol{R}}_{{h}}$ can be performed directly. Defining the probability distribution of the orientation by

(3.5) \begin{equation} \psi (\boldsymbol{n},\tilde {t})=\int {\rm d}\tilde {\boldsymbol{R}}_{{h}} \psi \left( \tilde {\boldsymbol{R}}_{{h}},\boldsymbol{n},\tilde {t}\right)\!, \end{equation}

with ${\rm d}\tilde {\boldsymbol{R}}_{{h}} = {\rm d}\tilde {x}_{{h}} {\rm d}\tilde {y}_{{h}} {\rm d}\tilde {z}_{{h}}$ , the integration of (2.17) gives

(3.6) \begin{equation} \left\langle \boldsymbol{\mathcal{F}}\big(\tilde {t}\big)\right\rangle = \int {\rm d}\tilde {\varOmega } \psi \left(\tilde {\boldsymbol{R}}_{{h}},\boldsymbol{n},\tilde {t}\right)\boldsymbol{\mathcal{F}}(\boldsymbol{n}) = \int {\rm d}\boldsymbol{n} \psi \left(\boldsymbol{n},\tilde {t}\right)\boldsymbol{\mathcal{F}}(\boldsymbol{n}). \end{equation}

The orientation vector $\boldsymbol{n}$ is expressed in spherical coordinates as $\boldsymbol{n}=\sin \theta \cos \phi \boldsymbol{e}_x +\sin \theta \sin \phi \boldsymbol{e}_y +\cos \theta \boldsymbol{e}_z$ with $\theta \in [0,\pi ]$ and $\phi \in [0,2\pi ]$ , illustrated in figure 1. Then the orientational distribution is expressed as $\psi (\boldsymbol{n},\tilde {t})=\psi (\theta ,\phi ,\tilde {t})$ , and ${\rm d}\boldsymbol{n}=\sin \theta {\rm d}\theta {\rm d}\phi$ . The equation governing $\psi (\theta ,\phi ,\tilde {t})$ is obtained by integrating both sides of (2.14a ) over $\tilde {\boldsymbol{R}}_{{h}}$ , which reads

(3.7) \begin{equation} \frac {\partial \psi }{\partial \tilde {t}}=\tilde {\mathcal{L}}_{{sp}}\psi = \frac {1}{\sin \theta }\frac {\partial }{\partial \theta }\left (\sin \theta \frac {\partial \psi }{\partial \theta }\right )+\frac {1}{\sin ^2\theta }\frac {\partial ^2\psi }{\partial \phi ^2}-\alpha \left (\sin \theta \frac {\partial \psi }{\partial \theta }+2\cos \theta \psi \right )\!, \end{equation}

with the periodic boundary condition $\left .\psi \right |_{\phi =0}=\left .\psi \right |_{\phi =2\pi }$ and boundedness conditions of $\psi$ at $\theta =0,\pi$ . The normalisation condition is $\int _0^{2\pi }{\rm d}\phi \int _0^\pi {\rm d}\theta \sin \theta \psi =1$ . The steady-state orientation distribution can be obtained analytically by solving $\tilde {\mathcal{L}}_{{sp}} \psi _{\textit{ss}} = 0$ , and the solution is

(3.8) \begin{equation} \psi _{\textit{ss}}(\theta ,\phi ) = \frac {\alpha }{4\pi \sinh \alpha }{\rm e}^{-\alpha \cos \theta } \end{equation}

(see also Brenner Reference Brenner1979; Dill & Brenner Reference Dill and Brenner1983). Figure 2 illustrates the distribution. The orientational distribution of a sedimenting particle in a steady-state is the same as the equilibrium distribution of a particle whose hydrodynamic centre is fixed in space. Due to the rotational symmetry about the $\hat {\boldsymbol{g}}$ direction, the steady-state orientation distribution is independent of $\phi$ . When the hydrodynamic centre coincides with the force centre ( $l_{{c}}=0$ ), $\alpha$ becomes equal to zero, and (3.8) reduces to the uniform distribution $\psi _{\textit{ss}}(\theta ,\phi ) = 1/(4\pi )$ (Frankel Reference Frankel1991), which represents a uniform distribution over a unit spherical surface.

Figure 2.

Steady-state orientation probability distribution $\psi _{\textit{ss}}(\theta , \phi )$ for an axisymmetric Brownian particle. Inset: the distribution weighted by $\sin \theta$ , accounting for the spherical area element associated with polar angle $\theta$ .

Figure 2. Long description

A line graph showing the steady-state orientation probability distribution for an axisymmetric Brownian particle. The main graph displays three lines representing different values of \alpha, with the x-axis labeled as theta ranging from 0 to pi and the y-axis labeled as psi subscript ss ranging from 0 to 0.3. The inset graph shows the distribution weighted by psi subscript ss sine theta, with the x-axis labeled as theta ranging from 0 to pi and the y-axis labeled as psi subscript ss sine theta ranging from 0 to 0.3. The lines in the main graph are color-coded: red for 1, blue for 2, and green for 10. The inset graph also features three lines corresponding to the same values. All values are approximated.

Using the steady-state distribution $\psi _{\textit{ss}}$ , the original (3.7) can be solved via the Green’s function $\mathcal{G}$ , determined through an eigenfunction expansion:

(3.9) \begin{gather} \mathcal{G}(\theta ,\phi ,\tilde {t};\theta ^\prime,\phi ^\prime) = \psi _{\textit{ss}}^{-1}(\theta ^\prime,\phi ^\prime)\sum _{p=0}^\infty \psi _p(\theta ,\phi )\psi _p(\theta ^\prime,\phi ^\prime){\rm e}^{-\lambda _p\tilde {t}}, \end{gather}

where $\lambda _p$ and $\psi _p$ ( $p=0,1,2, \dotsc$ ) are the eigenvalues and the eigenfunctions of the torque-dependent eigenvalue equation

(3.10) \begin{equation} \tilde {\mathcal{L}}_{{sp}}\psi _p=-\lambda _p\psi _p. \end{equation}

The eigenfunctions satisfy the orthonormality condition $\int _0^{2\pi }{\rm d}\phi \int _0^\pi {\rm d}\theta \sin \theta \psi _{\textit{ss}}^{-1}\psi _p\psi _q=\delta _{p,q}$ with weight $\psi _{\textit{ss}}^{-1}$ . It can generally be proved that the eigenvalues $\lambda _p$ are all non-negative and real (Doi & Edwards Reference Doi and Edwards1986). The sole zero eigenvalue $\lambda _0=0$ is included in (3.9), indicating that $\psi (\theta ,\phi ,\tilde {t})$ ultimately converges to the steady-state distribution $\psi _0=\psi _{\textit{ss}}(\theta ,\phi )$ over time. Other eigenfunctions $\psi _p$ ( $p=1,2,3,\dotsc$ ) are obtained numerically (see Appendix C), and they satisfy $\int _0^{2\pi }{\rm d}\phi \int _0^\pi {\rm d}\theta \sin \theta \psi _p=0$ for $p=1,2,3,\dotsc$ from the orthonormality. Therefore, using the Green’s function and the initial distribution of orientation $\psi _{\textit{in}}$ , the solution of (3.7) is

(3.11) \begin{align} \psi (\theta ,\phi ,\tilde {t}) &= \int _0^{2\pi }{\rm d}\phi ^\prime\int _0^\pi {\rm d}\theta ^\prime\sin \theta ^\prime\mathcal{G}(\theta ,\phi ,\tilde {t};\theta ^\prime,\phi ^\prime)\psi _{\textit{in}}(\theta ^\prime,\phi ^\prime) \notag \\ &= \psi _{\textit{ss}}(\theta ,\phi ) + \sum _{p=1}^\infty c_p \psi _p(\theta ,\phi ){\rm e}^{-\lambda _p\tilde {t}}, \end{align}

where $c_p = \int _0^{2\pi }{\rm d}\phi ^\prime\int _0^\pi {\rm d}\theta ^\prime\sin \theta ^\prime\psi _{\textit{ss}}^{-1}(\theta ^\prime,\phi ^\prime)\psi _{\textit{in}}(\theta ^\prime,\phi ^\prime)\psi _p(\theta ^\prime,\phi ^\prime)$ .

For example, ensemble averages $\langle n_z\rangle$ and $\langle n_z^2\rangle$ are calculated using the solution (3.11). The long-time behaviour can generally be obtained from the steady-state distribution function (3.8) as

(3.12a) \begin{align} \lim _{\tilde {t} \to \infty } \langle n_z\big(\tilde {t}\big)\rangle &= \langle n_z\rangle _{\textit{ss}} = \int _0^{2\pi }{\rm d}\phi \int _0^\pi {\rm d}\theta \sin \theta \psi _{\textit{ss}}(\theta ,\phi )\cos \theta =-\frac {\alpha \coth \alpha -1}{\alpha }, \end{align}

(3.12b) \begin{align} \lim _{\tilde {t} \to \infty } \langle n_z^2\big(\tilde {t}\big)\rangle &= \langle n_z^2\rangle _{\textit{ss}} = \int _0^{2\pi }{\rm d}\phi \int _0^\pi {\rm d}\theta \sin \theta \psi _{\textit{ss}}(\theta ,\phi )\cos ^2\theta =1-2\frac {\alpha \coth \alpha -1}{\alpha ^2}. \end{align}

It shows that the particle develops a preferred orientation along the direction of gravity when gravitational torque exists. Similarly, the steady-state ensemble average $\langle \boldsymbol{n}\boldsymbol{n}\rangle _{\textit{ss}}$ is obtained by

(3.13) \begin{equation} \langle \boldsymbol{n}\boldsymbol{n}\rangle _{\textit{ss}} = \frac {1}{3}\boldsymbol{\delta }+\left (1-3\frac {\alpha \coth \alpha -1}{\alpha ^2}\right )\left (\hat {\boldsymbol{g}}\hat {\boldsymbol{g}}-\frac {1}{3}\boldsymbol{\delta }\right )\!. \end{equation}

The transient behaviour before the system reaches the steady-state can also be calculated numerically via eigenfunction expansion of the Green’s function (3.9), or analytically via perturbation expansion for the moment $\langle (\hat {\boldsymbol{g}} \boldsymbol{\cdot }\boldsymbol{n})^2\rangle$ with respect to $\alpha$ (see (B1a )).

For small $\alpha$ , the latter method yields the following expression for $\langle n_z(\tilde {t})\rangle$ for a particle which is distributed isotropically at $\tilde t =0$ :

(3.14) \begin{align} \left\langle n_z\big(\tilde {t}\big)\right\rangle &= \boldsymbol{e}_z\boldsymbol{\cdot }\left\langle \boldsymbol{n}\big(\tilde {t}\big)\right\rangle = -\frac {\alpha }{3} \left (1-{\rm e}^{-2\tilde {t}} \right )+o\big(\alpha ^2\big). \end{align}

Equation (3.14) approaches the asymptotic value given by (3.12a ) as $\tilde {t} \to \infty$ .

For large $\alpha$ , the gravitational torque overwhelms rotational Brownian fluctuations, locking the particle into a nearly fixed orientation. If the particle is nearly in a gravity-aligned configuration, one has $n_x,n_y \ll 1$ , and $n_z=-[1-(n_x^2+n_y^2)]^{1/2}\approx -1+(n_x^2+n_y^2)/2$ . The equilibrium distribution (3.8) is reduced to $\psi _{\textit{ss}} \sim {\rm e}^{-\alpha n_z} \sim {\rm e}^{-\alpha (n_x^2+n_y^2)/2} \sim {\rm e}^{-U_\alpha /(k_{{B}}\kern-1pt T)}$ . This means that the particle experiences a restoring force for a small deviation from the fixed orientation. The restoring force can be evaluated from the gradient of a harmonic potential $U_\alpha =k_{{B}}\kern-1pt T\alpha (n_x^2+n_y^2)/2$ . Therefore, analytical calculation is possible since the dynamics of $\boldsymbol{n}(t)$ is the same as the Brownian motion of a particle trapped in a harmonic potential. The various time correlation functions for $\boldsymbol{n}(t)$ are obtained analytically as follows:

(3.15a) \begin{align} \left\langle n_x\big(\tilde {t}\big)n_z\big(\tilde {t}\big)n_x(0)n_z(0)\right\rangle &= \frac {1}{\alpha }{\rm e}^{-\alpha \tilde {t}}, \end{align}

(3.15b) \begin{align} \langle n_z^2\big(\tilde {t}\big)n_z^2(0)\rangle - \left\langle n_z^2\right\rangle ^2 &= \frac {4}{\alpha ^2}{\rm e}^{-2\alpha \tilde {t}}, \end{align}

where the initial equilibrium distribution has been used. These orientational correlation functions are useful for calculating the dispersion coefficients later.

3.3. Velocity of a sedimenting particle

We now turn to the motion along the $z$ direction during sedimentation. This is obtained by using (3.2). The sedimentation velocity $\tilde {u}_z$ is defined by

(3.16) \begin{equation} \tilde {u}_z = \lim _{\tilde {t} \to \infty }\frac {\left |\langle \tilde {z}_{{h}}\big(\tilde {t}\big)\rangle \right |}{\tilde {t}}. \end{equation}

The long-time behaviour is given by (3.2), which then yields

(3.17) \begin{align} \tilde {u}_z &= \beta \big ( 1+\chi \hat {\boldsymbol{g}}\boldsymbol{\cdot }\langle \boldsymbol{n}\boldsymbol{n}\rangle _{\textit{ss}}\boldsymbol{\cdot }\hat {\boldsymbol{g}} \big ) \notag \\ &=\beta \left [1+\chi \left (1-2\frac {\alpha \coth \alpha -1}{\alpha ^2}\right )\right ]\!, \end{align}

where (3.13) for $\langle \boldsymbol{n}\boldsymbol{n}\rangle _{\textit{ss}}$ has been used. The sedimentation velocity depends linearly on the sedimentation Péclet number $\beta$ . When $\alpha$ is large, the sedimentation velocity approaches $\lim _{\alpha \to \infty }\tilde {u}_z = \beta (1+\chi )$ , which corresponds to the sedimentation velocity of a rigid particle in a perfectly gravity-aligned configuration. The dimensional sedimentation velocity is given by $u_z = \tilde {u}_z L/\tau _{{r}}$ , which reads

(3.18a) \begin{align} u_z &= \frac {(M-M_{{b}}) g}{\zeta _{{t}}^{\perp }}\left [1+\chi \left (1-2\frac {\alpha \coth \alpha -1}{\alpha ^2}\right )\right ] \end{align}

(3.18b) \begin{align} &= \frac {(M-M_{{b}}) g}{\zeta _{{t}}^{\perp }} \left [1+\frac {\chi }{3}\left (1+\frac {2\alpha ^2}{15}\right )\right ]+o\big(\alpha ^2\big) . \end{align}

The second equality (3.18b ) holds when $\alpha$ is relatively small, indicating that the sedimentation velocity is increasing in $\alpha ^2$ for small gravitational torque.

The sedimentation velocity $u_z$ depends on various parameters: particle shape, mass and gravitational torque $\alpha$ . As an example, the sedimentation velocities of spheroids (defined by (D1)) of constant volume are shown in figure 3 for various aspect ratios $r$ and gravitational torques $\alpha$ . If the centrosymmetric spheroid density is uniform, the centre offset vanishes (thus gravitational torque $\alpha =0$ ). For such particles, it is known that spherical particles exhibit the maximum sedimentation velocity compared with prolate and oblate particles (see the black line in figure 3 a). This is because the average friction constant $[2(\zeta _{{t}}^\perp )^{-1}+(\zeta _{{t}}^\parallel )^{-1}]^{-1}$ of the spheroids is the smallest for a sphere. However, there is a very narrow band ( $1\lt r\lesssim 4$ ) where the sedimentation velocity of a prolate spheroid exceeds that of a sphere for non-vanishing $\alpha$ . The reason is that a non-vanishing $\alpha$ leads to an orientational preference along the gravity direction. Since the friction constant $\zeta _{{t}}^\parallel$ is smaller than $\zeta _{{t}}^\perp$ for a prolate spheroid, the sedimentation velocity could be larger than that of a sphere before the average friction constant effect becomes dominant. Figure 3(b) shows that with the increase of the gravitational torque $\alpha$ , the sedimentation velocity increases for prolate but decreases for oblate, which agrees with the early results of Brenner & Condiff (Reference Brenner and Condiff1972). This difference comes from the friction constant $\zeta _{{t}}^\parallel$ being smaller than $\zeta _{{t}}^\perp$ for prolate but being larger than that for oblate.

Figure 3.

Sedimentation velocity: (a) $u_z$ versus aspect ratio $r$ at various reorientation Péclet numbers $\alpha$ ; (b) $u_z$ versus reorientation Péclet number $\alpha$ at various aspect ratios $r$ .

Figure 3. Long description

The image contains two line graphs analyzing sedimentation velocity of non-spherical Brownian particles in viscous fluids. The first graph (a) plots sedimentation velocity against aspect ratio at various reorientation Peclet numbers, with curves for Peclet numbers 0, 5, 10, and 20. The second graph (b) plots sedimentation velocity against reorientation Peclet number at various aspect ratios, with curves for aspect ratios 0.01, 0.1, 1, 10, and 100. The graphs illustrate how sedimentation velocity changes with different aspect ratios and reorientation Peclet numbers, highlighting the complex behavior of non-spherical particles during sedimentation. The x-axes represent aspect ratio and reorientation Peclet number, while the y-axes represent normalized sedimentation velocity. The curves show how the velocity varies, indicating increased dispersion due to particle rotation and changes in velocity during sedimentation. All values are approximated.

To study the transient sedimentation behaviour, we compute the time evolution of $\langle \boldsymbol{n}\boldsymbol{n}\rangle$ . The ensemble average of $\boldsymbol{n}\boldsymbol{n}$ can be computed directly from definition (2.17) using solution (3.11) with $\boldsymbol{\mathcal{F}} = \boldsymbol{n}\boldsymbol{n}$ . Alternatively, using (2.18a ), the time evolution of the rotational correlation function is expressed as

(3.19) \begin{align} \frac {\partial \langle \boldsymbol{n}\boldsymbol{n}\rangle }{\partial \tilde {t}} = \left \langle \tilde {\mathcal{L}}^\dagger \boldsymbol{n}\boldsymbol{n} \right \rangle = -2\left (3\langle \boldsymbol{n}\boldsymbol{n}\rangle -\boldsymbol{\delta }\right ) +\alpha \left ( \hat {\boldsymbol{g}} \langle \boldsymbol{n}\rangle +\langle \boldsymbol{n}\rangle \hat {\boldsymbol{g}}-2\hat {\boldsymbol{g}} \boldsymbol{\cdot }\langle \boldsymbol{n}\boldsymbol{n}\boldsymbol{n} \rangle \right )\!. \end{align}

When the hydrodynamic centre and force centre coincide ( $\alpha =0$ ), the exact solution $\langle \boldsymbol{n}\boldsymbol{n}\rangle = \boldsymbol{\delta }/3$ holds for an initially isotropic state. For $\alpha \neq 0$ , we compute the time evolution of $\langle \boldsymbol{n}\boldsymbol{n}\rangle$ analytically by the perturbation method for small $\alpha$ (details appear in Appendix B). The perturbative expansion gives

(3.20) \begin{equation} \left\langle \boldsymbol{n}\big(\tilde {t}\big)\boldsymbol{n}\big(\tilde {t}\big)\right\rangle = \frac {1}{3}\boldsymbol{\delta }+\frac {\alpha ^2}{15}\left (1-\frac {3}{2}{\rm e}^{-2\tilde {t}}+\frac {1}{2}{\rm e}^{-6\tilde {t}}\right )\left (\hat {\boldsymbol{g}}\hat {\boldsymbol{g}}-\frac {1}{3}\boldsymbol{\delta }\right )+o\big(\alpha ^2\big), \end{equation}

for a system with initially isotropic distribution. This gives the following perturbation solution:

(3.21) \begin{align} \langle \tilde {z}_{{h}}\big(\tilde {t}\big)\rangle &= -\tilde {u}_z\tilde {t} + \frac {4\alpha ^2}{135} \beta \chi \left (1-\frac {9}{8}{\rm e}^{-2\tilde {t}} + \frac {1}{8}{\rm e}^{-6\tilde {t}}\right ) + o\big(\alpha ^2\big), \end{align}

where $\tilde {u}_z$ is now given by (3.18b ).

3.4. Dispersion of a sedimenting particle

We now calculate the diffusion in the $x$ – $y$ plane and $z$ direction. As the orientational distribution approaches the steady-state distribution, the ensemble averages in (3.4) and (3.3) converge to constant values. We therefore define the horizontal diffusion coefficient $\tilde {D}_{\text{$x$-$y$}}$ and the vertical diffusion coefficient $\tilde {D}_z$ by

(3.22a) \begin{align} \tilde {D}_{\text{$x$-$y$}} &= \lim _{\tilde {t} \to \infty } \frac {\left \langle \tilde {x}_{{h}}^2\big(\tilde {t}\big)+\tilde {y}_{{h}}^2\big(\tilde {t}\big)\right \rangle }{4 \tilde {t}}, \end{align}

(3.22b) \begin{align} \tilde {D}_z &= \lim _{\tilde {t} \to \infty } \frac {\left\langle \left [\tilde {z}_{{h}}\big(\tilde {t}\big)-\big\langle \tilde {z}_{{h}}\big(\tilde {t}\big)\big\rangle \right ]^2\right\rangle }{2 \tilde {t}}. \end{align}

Equations (3.4) and (3.3) then yield the following expressions for the horizontal and vertical diffusion coefficients (see Appendix C):

(3.23a) \begin{align} \tilde {D}_{\text{$x$-$y$}} &= \tilde {D}^\perp \left [1+\frac {\chi }{2}\big(1-\hat {\boldsymbol{g}}\boldsymbol{\cdot }\langle \boldsymbol{n}\boldsymbol{n}\rangle _{\textit{ss}}\boldsymbol{\cdot }\hat {\boldsymbol{g}}\big)\right ] +\left (\beta \chi \right )^2\lim _{\tilde {t}\to \infty }\varXi (\tilde {t};\alpha ) \notag \\ &= \tilde {D}^\perp \left (1+\chi \frac {\alpha \coth \alpha -1}{\alpha ^2}\right ) +\left (\beta \chi \right )^2\varXi _{\textit{ss}}(\alpha ) , \end{align}

(3.23b) \begin{align} \tilde {D}_z &= \tilde {D}^\perp \left (1+\chi \hat {\boldsymbol{g}}\boldsymbol{\cdot }\langle \boldsymbol{n}\boldsymbol{n}\rangle _{\textit{ss}}\boldsymbol{\cdot }\hat {\boldsymbol{g}}\right ) +\left (\beta \chi \right )^2\lim _{\tilde {t}\to \infty }\varTheta (\tilde {t};\alpha ) \notag \\ &= \tilde {D}^\perp \left (1+\chi -2\chi \frac {\alpha \coth \alpha -1}{\alpha ^2}\right ) +\left (\beta \chi \right )^2\varTheta _{\textit{ss}}(\alpha ) , \end{align}

where (3.13) for $\langle \boldsymbol{n}\boldsymbol{n}\rangle _{\textit{ss}}$ has been used. Here $\varXi _{\textit{ss}}(\alpha )$ and $\varTheta _{\textit{ss}}(\alpha )$ are defined by

(3.24a) \begin{align} \varXi _{\textit{ss}}(\alpha ) &= \int _0^{\infty }{\rm d}\tilde {t} \big\langle n_x\big(\tilde {t}\big)n_z\big(\tilde {t}\big)n_x(0)n_z(0) \big\rangle , \end{align}

(3.24b) \begin{align} \varTheta _{\textit{ss}}(\alpha ) &= \int _0^{\infty }{\rm d}\tilde {t} \left [\left\langle n_z^2\big(\tilde {t}\big)n_z^2(0) \right\rangle -\left\langle n_z^2 \right\rangle ^2\right ]\!, \end{align}

which can be calculated numerically by (C4a ) and (C4b ). These are the key quantities to characterise the orientational time correlation in steady-state for the dispersions. However, the analytical expressions of (3.24a ) and (3.24b ) are difficult to obtain in general and are only valid for small and large $\alpha$ . We can use the iterative method in Appendix C to obtain the expressions in series of $\alpha$ . For small $\alpha$ , the results up to order of $\alpha ^2$ are

(3.25) \begin{equation} \varXi _{\textit{ss}}(\alpha ) = \frac {1}{90}\left (1+\frac {59}{252}\alpha ^2\right )\!,\qquad \varTheta _{\textit{ss}}(\alpha ) = \frac {2}{135}\left (1+\frac {5}{14}\alpha ^2\right ) \qquad (\alpha \ll 1). \end{equation}

Notice that (3.25) reduces to the classical work of Brenner (Reference Brenner1979) in the limit of $\alpha \to 0$ .

For large $\alpha$ , we can use (3.15). The results are

(3.26) \begin{equation} \varXi _{\textit{ss}}(\alpha ) = \frac {1}{\alpha ^2}, \qquad \varTheta _{\textit{ss}}(\alpha ) = \frac {2}{\alpha ^3} \qquad (\alpha \gg 1). \end{equation}

Figure 4 illustrates the functions $\varXi _{\textit{ss}}(\alpha )$ and $\varTheta _{\textit{ss}}(\alpha )$ , where analytic solutions are shown in blue lines.

Figure 4.

Plots of (a) $\varXi _{\textit{ss}}(\alpha )$ defined in (3.24a ) and (b) $\varTheta _{\textit{ss}}(\alpha )$ defined in (3.24b ). The blue curves show the analytic solutions.

Figure 4. Long description

Two line graphs illustrate the analytic solutions for sedimentation of a non-spherical Brownian particle in viscous fluids. The left graph shows the relationship between a variable and alpha, with a peak around alpha equals 2 and a subsequent decline. The right graph depicts another variable against alpha, starting high and decreasing sharply before leveling off. Both graphs feature blue curves representing the analytic solutions, with specific mathematical expressions annotated at key points.

The apparent diffusivity $\tilde {D}_{\text{$x$-$y$}}$ (or $\tilde {D}_z$ ) originates from two distinct physical mechanisms. The first mechanism is normal diffusion due to Brownian motion, represented by the first term in (3.23a ) for $\tilde {D}_{\text{$x$-$y$}}$ or (3.23b ) for $\tilde {D}_z$ . The second mechanism is gravity-induced Taylor dispersion due to hydrodynamic anisotropy, corresponding to the second term in (3.23a ) for $\tilde {D}_{\text{$x$-$y$}}$ or (3.23b ) for $\tilde {D}_z$ . The torque acting on the particle affects both normal diffusion and Taylor dispersion. When $\alpha =0$ , the well-known Taylor dispersion of Brenner (Reference Brenner1979) is recovered: the diffusivities have terms proportional to $\beta ^2$ , which can be much larger than the terms of Brownian diffusion. For finite $\beta$ , the diffusivities approach limiting values $\lim _{\alpha \to \infty } \tilde {D}_{\text{$x$-$y$}}=\tilde {D}^\perp$ and $\lim _{\alpha \to \infty } \tilde {D}_z=\tilde {D}^\perp (1+\chi )$ . These limiting diffusivities correspond to the normal diffusion constants for a rigid particle in a perfectly gravity-aligned configuration.

Note that the original parameter $\beta$ in (2.13) depends on the aspect ratio $r$ , and here we define a new parameter $\beta _0=(M-M_{{b}}) gL/(k_{{B}}\kern-1pt T)$ to represent solely the magnitude of the external force. In addition, the original parameter $\alpha$ in (2.13) depends on $\beta _0$ as well, and here we define a new parameter $\epsilon =l_{{c}}/L$ to represent solely the dimensionless centre offset. Therefore, we have

(3.27) \begin{equation} \alpha =\beta _0\epsilon , \qquad \beta =\tilde {D}^\perp \beta _0. \end{equation}

For completeness, the diffusivities with dimension ( $D_{\text{$x$-$y$}} = \tilde {D}_{\text{$x$-$y$}} L^2/\tau _{{r}}$ , $D_z = \tilde {D}_z L^2/\tau _{{r}}$ ; see (2.16)) are given by

(3.28a) \begin{align} D_{\text{$x$-$y$}} &= D^\perp \left (1+\chi \frac {(\beta _0\epsilon )\coth (\beta _0\epsilon )-1}{(\beta _0\epsilon )^2} +\tilde {D}^\perp \chi ^2\beta _0^2\varXi _{\textit{ss}}(\beta _0\epsilon )\right )\!,\\[-6pt]\nonumber \end{align}

(3.28b) \begin{align} D_z &= D^\perp \left (1+\chi -2\chi \frac {(\beta _0\epsilon )\coth (\beta _0\epsilon )-1}{(\beta _0\epsilon )^2} +\tilde {D}^\perp \chi ^2\beta _0^2\varTheta _{\textit{ss}}(\beta _0\epsilon )\right )\!. \end{align}

Figure 5.

(a) The horizontal diffusion coefficient $D_{\text{$x$-$y$}}$ and (b) the vertical diffusion coefficient $D_z$ as functions of the gravity $\beta _0=(M-M_{{b}}) gL/(k_{{B}}\kern-1pt T)$ with fixed dimensionless centre offset $\epsilon =l_{{c}}/L$ . The dashed lines are analytical solutions for large $\beta _0$ . Coefficients (c) $D_{\text{$x$-$y$}}$ and (d) $D_z$ as functions of $\epsilon$ with fixed $\beta _0$ .

Figure 5. Long description

The image contains four graphs. The first graph (a) shows the horizontal diffusion coefficient as a function of gravity with a fixed dimensionless center offset. The second graph (b) shows the vertical diffusion coefficient as a function of gravity with the same fixed dimensionless center offset. The dashed lines in both graphs represent analytical solutions for large values. The third graph (c) and the fourth graph (d) show the diffusion coefficients as functions of the dimensionless center offset with a fixed gravity value. Each graph uses a logarithmic scale for the x-axis and a linear scale for the y-axis. The colors of the lines in the graphs represent different values of the dimensionless center offset or gravity, as indicated by the labels on the lines.

For small torque $\beta _0\epsilon \ll 1$ , the first-order perturbations of diffusivities $D_{\text{$x$-$y$}}$ and $D_z$ with respect to $(\beta _0\epsilon )^2$ are given by

(3.29a) \begin{align} D_{\text{$x$-$y$}} &= D^\perp \left [1+\frac {\chi }{3}\left (1-\frac {(\beta _0\epsilon )^2}{15}\right ) +\frac {\tilde {D}^\perp \chi ^2\beta _0^2}{90}\left (1+\frac {59(\beta _0\epsilon )^2}{252}\right )\right ]\!,\\[-6pt]\nonumber \end{align}

(3.29b) \begin{align} D_z &= D^\perp \left [1+\frac {\chi }{3}\left (1+\frac {2(\beta _0\epsilon )^2}{15}\right ) +\frac {2\tilde {D}^\perp \chi ^2\beta _0^2}{135}\left (1+\frac {5(\beta _0\epsilon )^2}{14}\right )\right ]\!. \end{align}

They indicate that for small torque, the Taylor dispersion of both diffusivities increases with the dimensionless centre offset $\epsilon$ .

For large torque $\beta _0\epsilon \gg 1$ , the asymptotic behaviour of diffusivities $D_{\text{$x$-$y$}}$ and $D_z$ can be calculated using (3.26):

(3.30a) \begin{align} D_{\text{$x$-$y$}} &= D^\perp \left (1+\frac {\chi }{\beta _0\epsilon } +\frac {\tilde {D}^\perp \chi ^2}{ \epsilon ^2}\right )\!, \end{align}

(3.30b) \begin{align} D_z &= D^\perp \left (1+\chi -\frac {2\chi }{\beta _0\epsilon } +\frac {2\tilde {D}^\perp \chi ^2}{\beta _0\epsilon ^3}\right )\!. \end{align}

The limit behaviour of diffusivities as $\beta _0\to \infty$ with $\epsilon$ fixed can be understood as follows. The diffusivity $D_{\text{$x$-$y$}}$ is estimated by $\langle u_x^2\rangle \tau _{{o}}$ in the ballistic regime of MSD (see figure 7 below) since the Péclet number is considerably large, where $u_x$ is the velocity in the $x$ direction. One has $u_x=un_x\sim \beta _0n_x$ since the sedimentation velocity $u$ is proportional to the magnitude of gravity $\beta _0$ . Therefore, the diffusivity is given by $D_{\text{$x$-$y$}}\sim \beta _0^2\langle n_x^2\rangle \tau _{{r}}\alpha ^{-1} \sim \epsilon ^{-2}$ , where $\langle n_x^2\rangle \sim \alpha ^{-1}$ from the same procedure as in (3.12b ) and $\alpha =\beta _0\epsilon$ from (3.27) have been used. Similarly, $D_z$ is estimated by $\langle [u_z-\langle u_z\rangle ]^2\rangle \tau _{{o}}$ and one has $u_z=un_z\sim \beta _0n_z$ . Therefore, the diffusivity is given by $D_z\sim \beta _0^2\langle [n_z-\langle n_z\rangle ]^2\rangle \tau _{{r}}\alpha ^{-1} \sim \beta _0^{-1}\epsilon ^{-3}$ , where $\langle [n_z-\langle n_z\rangle ]^2\rangle \sim \alpha ^{-2}$ from the same procedure as in (3.12b ) has been used.

Figure 6.

(a) Horizontal ( $D_{x-y}$ ) and vertical ( $D_z$ ) diffusion coefficients versus aspect ratio $r$ at $\alpha =0$ (in the absence of gravitational torque). Here, $\beta _0=(M-M_{{b}}) gL/(k_{{B}}\kern-1pt T)$ . (b) Dependence of horizontal ( $D_{x-y}$ ) and vertical ( $D_z$ ) diffusion coefficients on reorientation Péclet number $\alpha$ , showing Brownian and gravity-induced contributions. The maximum diffusion coefficients $D_{{max}}$ are the peak values of the black lines in (b). Maximum of horizontal and vertical diffusion coefficients as functions of (c) aspect ratio $r$ and (d) sedimentation Péclet number $\beta _0$ .

Figure 6. Long description

The image contains four graphs analyzing diffusion coefficients. Graph (a) shows horizontal and vertical diffusion coefficients versus aspect ratio at zero gravitational torque, with solid and dashed lines representing different coefficients. Graph (b) illustrates the dependence of horizontal and vertical diffusion coefficients on the reorientation Peclet number, highlighting Brownian and gravity-induced contributions. Graph (c) presents the maximum diffusion coefficients as functions of aspect ratio, while graph (d) shows these coefficients as functions of the sedimentation Peclet number. Each graph uses solid and dashed lines to differentiate between horizontal and vertical diffusion coefficients.

Figures 5(a) and 5(b) illustrate, respectively, the horizontal diffusion coefficient $D_{\text{$x$-$y$}}$ in (3.28a ) and the vertical diffusion coefficient $D_z$ in (3.28b ) as functions of the sedimentation Péclet number $\beta _0$ with fixed dimensionless centre offset $\epsilon$ . The dashed lines are analytical solutions in (3.30a ) and (3.30b ) for large $\beta _0$ . The black curves show Brenner’s results for the torque-free case. If there is a centre offset with a large $\beta _0$ , a significant deviation from the $\beta _0^2$ scaling appears even for a small centre offset. For sufficiently large $\beta _0$ , the horizontal diffusion coefficients $D_{\text{$x$-$y$}}$ approach a constant value which depends on the offset, while the vertical diffusion coefficients $D_z$ all converge to the same value. This indicates that Taylor dispersion is suppressed when the gravitational torque dominates over rotational Brownian fluctuations. Figures 5(c) and 5(d) illustrate $D_{\text{$x$-$y$}}$ and $D_z$ as functions of $\epsilon$ with fixed $\beta _0$ , respectively. The results show that for small centre offset, the Taylor dispersion effect is enhanced compared with the torque-free case, owing to orientational preference induced by the torque. However, if the centre offset is too large, the Taylor dispersion is suppressed since the particle loses the orientation variability required to generate Taylor dispersion.

Figure 6(a) shows the results of Brenner (Reference Brenner1979) for the situation of $\alpha =0$ , i.e. the force centre is at the hydrodynamic centre. A similar graph was shown by Goren (Reference Goren1979). It is known that for particles of constant volume, the diffusivity $D_{\text{$x$-$y$}}$ (or $D_z$ ) of neutrally buoyant ( $\beta _0=0$ ) particles is maximised for spherical particles, since the average friction constant $[2(\zeta _{{t}}^\perp )^{-1}+(\zeta _{{t}}^\parallel )^{-1}]^{-1}$ is smallest for spherical particles. For large $\beta _0$ , gravity-induced Taylor dispersion dominates over ordinary diffusion, but no Taylor dispersion occurs for spheres. Since the coefficient $D^{\perp }/(D_{{r}}L^2)$ scales as $r^{4/3}$ for $r \gg 1$ and $r^{-2/3}$ for $r \ll 1$ , the diffusivity is minimised for spheres. Therefore, for sedimenting particles with intermediate $\beta _0$ , competition appears between the ordinary diffusion and the Taylor dispersion so that the minimum occurs slightly over into the prolate range. Both $D_{\text{$x$-$y$}}$ and $D_z$ exhibit similar behaviour.

In figure 6(b), the diffusivities of particles with gravitational torque ( $\alpha \neq 0$ ) normalised by that without torque are shown as a function of $\alpha$ , where aspect ratio $r=10$ and sedimentation Péclet number $\beta _0=10$ are chosen. Here contributions of the two mechanisms, the normal diffusion caused by thermal motion and the Taylor dispersion induced by gravity, are shown separately. It is seen that for small gravitational torque, the Taylor dispersion is enhanced, but as $\alpha$ increases, the Taylor dispersion starts to be suppressed. This occurs because the gravitational torque overwhelms rotational Brownian fluctuations, locking the particle into a nearly fixed orientation and thereby removing the orientational variability required to generate Taylor dispersion. As a consequence both $D_{\text{$x$-$y$}}$ and $D_z$ show maximum at a non-zero value of $\alpha$ . As $\alpha$ goes to infinity, the diffusivities approach their asymptotic values, i.e. $D_{\text{$x$-$y$}} \to D^\perp$ and $D_{\text{$z$}} \to D^\perp (1+\chi )$ .

The $\alpha$ dependence of the diffusivities changes when the other parameters change. The maximum diffusion coefficients $D_{{max}}$ are defined as the peak values of the black lines in figure 6(b). The maximum diffusivity takes place around the same value of $\alpha \sim 2$ for both prolate and oblate spheroids, as well as for both horizontal and vertical diffusivities. Note that the maximum diffusivity at non-vanishing $\alpha$ exceeds the classical Taylor dispersion (at $\alpha =0$ ) due to the orientational preference of the particles. The maximum diffusivity itself varies with aspect ratio $r$ and sedimentation Péclet number $\beta _0$ . Figures 6(c) and 6(d) show the maximum diffusivity as a function of $r$ and of $\beta _0$ , respectively. The enhanced Taylor dispersion is larger than the Taylor dispersion without a centre offset and saturates at approximately $150\,\%$ of the dispersion without centre offset. This implies that one can further enhance particle dispersion by tuning $\alpha$ to $\sim 2$ , for example, by adjusting the mass distribution of the particles.

3.5. Transient behaviour of mean square displacement

We now compute the transient MSD to characterise dynamic evolution. First we consider the case of $\alpha =0$ . Since the orientational distribution is isotropic when $\alpha =0$ , the right-hand sides of (3.4) and (3.3) can be computed exactly for a system with initially isotropic distribution. The results of the horizontal and vertical MSD are

(3.31a) \begin{align} \left\langle \tilde {x}_{{h}}^2\big(\tilde {t}\big)+\tilde {y}_{{h}}^2\big(\tilde {t}\big) \right\rangle &= 4\left [\tilde {D}^\perp \left (1+\frac {\chi }{3}\right ) +\frac {\left (\beta \chi \right )^2}{90}\right ] \tilde {t}- \frac {1}{15} \left (\frac {\beta \chi }{3}\right )^2\left (1-{\rm e}^{-6\tilde {t}}\right )\!,\\[-12pt]\nonumber \end{align}

(3.31b) \begin{align} \left\langle \left [\tilde {z}_{{h}}\big(\tilde {t}\big)-\big\langle \tilde {z}_{{h}}\big(\tilde {t}\big)\big\rangle \right ]^2\right\rangle &= 2\left [\tilde {D}^\perp \left (1+\frac {\chi }{3}\right ) +\frac {2\left (\beta \chi \right )^2}{135} \right ] \tilde {t}- \frac {2}{45} \left (\frac {\beta \chi }{3}\right )^2\left (1-{\rm e}^{-6\tilde {t}}\right )\!. \end{align}

The transient MSD for $\alpha \neq 0$ is obtained from (C1a ) (or (C1b )) by perturbation up to $\alpha ^2$ using Laplace transform methods (C17).

Figure 7.

(a) The time evolution of the horizontal (solid line) and vertical (dashed line) MSD for varying sedimentation Péclet number $\beta$ ( $\alpha =0$ , $r=10$ ). The crossover time $\tau _{{\textit{cross}}}$ is defined by inverse of the smallest non-zero eigenvalue of (C6). (b) The scaled crossover time (by $\tau _{{r}}$ ) depends on $\alpha$ for both the horizontal (solid line) and vertical (dashed line) MSD. The analytical results for large $\alpha$ (from (3.15)) are represented by blue lines. Note that $1/\alpha =\tau _{{o}}/\tau _{{r}}$ (see (2.13)).

Figure 7. Long description

The image contains two graphs. The first graph on the left shows the time evolution of the horizontal and vertical mean squared displacement (MSD) for varying sedimentation Peclet numbers. The horizontal MSD is represented by solid lines, while the vertical MSD is represented by dashed lines. The graph includes different colors to denote varying Peclet numbers, with labels such as 0, 100, 300, and 1000. The x-axis is labeled with a tilde t, and the y-axis is labeled with MSD. The crossover time is defined by the inverse of the smallest non-zero eigenvalue of the eigenvalue. The second graph on the right shows the scaled crossover time divided by tau r as a function of alpha for both the horizontal and vertical MSD. The horizontal MSD is represented by a solid line, and the vertical MSD is represented by a dashed line. The analytical results for large alpha are represented by blue lines. The x-axis is labeled with alpha, and the y-axis is labeled with tau cross over tau r. The graph includes annotations indicating the relationships 1 over alpha and 1 over 2 alpha.

Figure 7(a) shows the MSD of a centrosymmetric ( $\alpha =0$ ) spheroid of $r=10$ , calculated by (3.31) and the formula for the friction constants in Appendix D. For $\beta =0$ , the MSD increases linearly with $t$ , showing the usual diffusion behaviour. As $\beta$ increases (e.g. $100$ ), the MSD starts to show distinct transient behaviour: the MSD first increases linearly with $\tilde {t}$ , then increases in proportion to $\tilde {t}^2$ , and finally shows the diffusion behaviour. Such behaviour can be understood from (3.31a ). For small $\tilde {t}$ , the right-hand side of (3.31a ) can be expanded with respect to $\tilde {t}$ :

(3.32) \begin{equation} \left\langle \tilde {x}_{{h}}^2\big(\tilde {t}\big)+\tilde {y}_{{h}}^2\big(\tilde {t}\big) \right\rangle = 4\tilde {D}^\perp \left (1+\frac {\chi }{3}\right )\tilde {t} + \frac {2(\beta \chi )^2}{15}\tilde {t}^2 + o\big(\tilde {t}^2\big). \end{equation}

This gives the linear dependence for $\tilde {t} \lt \beta ^{-2}$ , and the quadratic dependence for $\tilde {t} \gt \beta ^{-2}$ . Therefore, the MSD crosses over from the usual diffusion scaling ( $\sim \tilde {t}$ ) to the ballistic scaling ( $\sim \tilde {t}^2$ ) around $\tilde {t}= \beta ^{-2}$ . For large time $\tilde {t}$ , the right-hand side of (3.31a ) exhibits diffusive behaviour driven by both thermal diffusion and Taylor dispersion, i.e. $4 [\tilde {D}^\perp (1+\chi /3 ) + (\beta \chi )^2/90 ] \tilde {t}$ . There is a crossover time $\tau _{{\textit{cross}}}$ from the ballistic behaviour to the final diffusion behaviour. The dimensionless crossover time $\tau _{{\textit{cross}}}$ (scaled by $\tau _{{r}}$ ) is defined by the inverse of the smallest non-zero eigenvalue of (C6) with $m=1$ for horizontal dispersion and $m=0$ for vertical dispersion.

Figure 7(b) shows the dimensionless crossover time as a function of $\alpha$ . The crossover time indicates the dynamic transition controlled by the orientation of particles, while the orientation dynamics is affected by both rotational diffusion and gravitational torque. For large $\alpha$ , the rotation is governed by gravitational torque, and the crossover time $\tau _{{\textit{cross}}}$ approaches the reorientation time $\tau _{{o}}$ . As can be observed in figure 7(b), $\tau _{{\textit{cross}}}/\tau _{{r}}$ approaches $1/\alpha$ and $1/(2\alpha )$ for the horizontal and vertical MSD, respectively. Equation (3.15) clearly indicates the relevant time scales ( $1/\alpha$ and $1/(2\alpha )$ ) for large $\alpha$ (note that $1/\alpha =\tau _{{o}}/\tau _{{r}}$ ; see (2.13)). In the absence of Brownian motion, $\tau _{{o}}$ diverges as $\alpha \to 0$ , which explains why the settling behaviour of non-Brownian particles is sensitive to centre offsets. If Brownian motion is present, however, for small $\alpha$ the rotation is governed by rotational diffusion and the crossover time $\tau _{{\textit{cross}}}$ tends to $\tau _{{r}}/2$ . The divergence is smoothed out by rotational diffusion, which means that the settling behaviour of Brownian particles is no longer sensitive to centre offsets.

As $\beta$ increases to infinity, we can recover the non-Brownian limit. To analyse the situation, we define $\tilde {t}_{{s}} = t / \tau _{{s}}$ (i.e. $\tilde {t}_{{s}} = \mathrm{\beta } \, \tilde {t}$ ). Equations (C1a ) and (C8) in the limit $\beta \to \infty$ are

(3.33) \begin{equation} \frac {\partial \left \langle \tilde {x}_{{h}}^2+\tilde {y}_{{h}}^2 \right \rangle }{\partial \tilde {t}_{{s}}} = 2\chi \mathcal{H}^1, \qquad \frac {\partial \mathcal{H}^1}{\partial \tilde {t}_{{s}}} = \frac {2\chi }{15}, \end{equation}

where the initial condition with the isotropic orientation distribution has been used. Solving (3.33) yields

(3.34) \begin{equation} \lim _{\beta \to \infty } \left \langle \tilde {x}_{{h}}^2 + \tilde {y}_{{h}}^2 \right \rangle = \frac {2\chi ^2}{15} \tilde {t}_{{s}}^2. \end{equation}

Both short-time and long-time diffusion behaviours in figure 7(a) vanish. The result demonstrates $\left \langle \tilde {x}_{{h}}^2+\tilde {y}_{{h}}^2 \right \rangle \sim \left \langle \tilde {z}_{{h}}^2 \right \rangle \sim \tilde {t}^2$ , and is valid even for a non-zero value of $\alpha$ . This occurs because rotational Brownian motion is too weak to disrupt the gliding motion of the particle, yielding ballistic $\tilde {t}^2$ scaling.