Hydrodynamic interactions and the diffusivity of spheroidal particles

TL;DR

This paper investigates how hydrodynamic interactions affect the translational and rotational diffusivities of spheroidal particles with arbitrary aspect ratios in dilute suspensions, revealing small anisotropic effects relevant to various scientific fields.

Contribution

It provides the first calculation of hydrodynamic interaction effects on spheroidal particle diffusivity, extending prior spherical particle models to anisotropic shapes.

Findings

Hydrodynamic interactions cause a few percent anisotropy in diffusivity.

The stresslet field is the origin of anisotropic effects.

Results are applicable to colloidal suspensions and cytoplasm dynamics.

Abstract

It is intuitive that the diffusivity of an isolated particle differs from those in a monodisperse suspension, in which hydrodynamic interactions between the particles are operative. Batchelor (1976,1983) calculated how hydrodynamic interactions influenced the diffusivity of a dilute suspension of spherical particles and Russel et al.(1991), and Brady (1994) treated non-dilute (higher particle volume fraction) suspensions. Although most particles lack perfect sphericity, little is known about the effects of hydrodynamic interactions on the diffusivity of spheroidal particles, which are the simplest shapes that can be used to model anisotropic particles. Here, we calculate the effects of hydrodynamic interactions on the translational and rotational diffusivities of spheroidal particles of arbitrary aspect ratio, in dilute monodisperse suspensions. The origin of the hydrodynamic anisotropy is that found in the stresslet field for the induced-dipole induced-dipole interaction. However, in the dilute limit the anisotropy effects are at the level of a few percent. These effects have influence in a vast range of settings, from partially frozen colloidal suspensions to the dynamics of cytoplasm.

Equations

$${\mathbf{D}}_{\mathbf{tr}}=\frac{k_{b}T}{6\pi \mu L}\left[\frac{1}{X_A}{\mathbf{p}}_1{\mathbf{p}}_1+\frac{1}{Y_A}(\mathbf{I}-{\mathbf{p}}_1{\mathbf{p}}_1)\right]$$

$$D_r=\frac{k_{b}T}{8\pi \mu L^3}\frac{1}{Y_C}$$

$$M_{ij}{\dot{X}}_j-(I_s-I_t)X_i^I=-k_{b}T{D}_{ij}^{-1}{X}_j+{\alpha }_{ij}{g}_j,$$

$$\mathbf{X}={\left[\begin{array}{c}{v}_{x_1},v_{y_1},\dots ,v_{z_N},{\omega }_{x_1},{\omega }_{y_1},{\omega }_{z_1},\dots ,{\omega }_{y_N},{\omega }_{z_N}\end{array}\right]}^T,$$

$$\mathbf{g}={\left[\begin{array}{c}{F}_{x_1},F_{y_1},\dots ,F_{z_N},{\tau }_{x_1},{\tau }_{y_1},{\tau }_{z_1},\dots ,{\tau }_{y_N},{\tau }_{z_N}\end{array}\right]}^T.$$

$$⟨g_i(t)g_j(t^{\prime })⟩=2{\delta }_{ij}\delta (t-t^{\prime }).$$

$$\mathbf{X}=$$

$${\mathcal{F}}^h=$$

$$\mathbf{D}=$$

$$(k_{b}T{)}^2{D}_{ij}^{-1}={\alpha }_{ik}{\alpha }_{jk},$$

$$⟨\Delta r_p^1\Delta r_q^1⟩=2⟨D_{pq}^{tt}⟩t,$$

$$⟨\Delta {\theta }_p^1\Delta {\theta }_q^1⟩=2⟨D_{pq}^{rr}⟩t.$$

$$⟨D_{pq}^{tt}⟩=\frac{1}{4\pi V}{\int }_0^{\pi }{\int }_0^{2\pi }{\int }_0^{\pi }{\int }_0^{2\pi }{\int }_{R_{\text{restrict}}}^{\infty }{D}_{pq}^{tt}{R}_p^2\sin {\theta }_R$$

$$\sin {\theta }_{2}d{R}_{p}d{\varphi }_{2}d{\theta }_{R}d{\varphi }_{R}d{\theta }_2,$$

$$⟨D_{pq}^{rr}⟩=\frac{1}{4\pi V}{\int }_0^{\pi }{\int }_0^{2\pi }{\int }_0^{\pi }{\int }_0^{2\pi }{\int }_{R_{\text{restrict}}}^{\infty }{D}_{pq}^{rr}{R}_p^2\sin {\theta }_R$$

$$\sin {\theta }_{2}d{R}_{p}d{\varphi }_{2}d{\theta }_{R}d{\varphi }_{R}d{\theta }_2.$$

$${\mathbf{D}}_{\mathbf{tr}}^{\mathbf{f}}=$$

$$D_r^f=$$

$$D_a=\frac{D_{\parallel }^s+\varphi D_{\parallel }}{D_{\perp }^s+\varphi D_{\perp }},$$

$$\frac{d}{dt}⟨\mathbf{X}(t)\mathbf{X}(t)⟩=2⟨\mathbf{X}(t)\dot{\mathbf{X}}(t)⟩.$$

$$⟨{\dot{X}}_l(t)X_f(t)⟩-(I_s-I_t)M_{li}^{-1}⟨X_i^I(t)X_f(t)⟩=-k_{b}T{M}_{li}^{-1}{D}_{ij}^{-1}$$

$$⟨X_j(t)X_f(t)⟩+M_{li}^{-1}{\alpha }_{ij}⟨g_j(t)X_f(t)⟩.$$

$$k_{b}T{D}_{ij}^{-1}⟨X_j(t)X_f(t)⟩-{\alpha }_{ij}⟨g_j(t)X_f(t)⟩=0.$$

$$⟨g_j(t)F_f[\mathbf{g}]⟩={\int }_0^{t}d{t}^{\prime }⟨g_j(t)g_k(t^{\prime })⟩⟨\frac{\delta F_f[\mathbf{g}]}{\delta g_k(t^{\prime })}⟩.$$

$$⟨\frac{\delta F_f[\mathbf{g}]}{\delta g_k(t^{\prime })}⟩$$

$$(k_{b}T{)}^2{D}_{ij}^{-1}{M}_{jf}^{-1}=2{\alpha }_{ij}H(0){\alpha }_{pj}{M}_{pf}^{-1}.$$

$${\dot{X}}_l+k_{b}T{M}_{li}^{-1}{D}_{ij}^{-1}{X}_j=(I_s-I_t)M_{li}^{-1}{X}_i^I+M_{li}^{-1}{\alpha }_{ij}{g}_j.$$

$$\frac{d}{dt}\left(e_{ml}^{Qt}{X}_l\right)=e_{ml}^{Qt}\left[(I_s-I_t)M_{li}^{-1}{X}_i^I+M_{li}^{-1}{\alpha }_{ij}{g}_j\right].$$

$$\Delta r_p^1=$$

$$⟨\Delta r_p^1\Delta r_q^1⟩=2⟨D_{pq}^{tt}⟩t,$$

$$⟨\Delta {\theta }_p^1\Delta {\theta }_q^1⟩=2⟨D_{pq}^{rr}⟩t.$$

$${\mathbf{U}}_1=$$

Full text

††preprint: AIP/123-QED

I Introduction

Confinement of single particles and particle-particle interactions in suspensions, share a common set of hydrodynamic phenomena through which the single particle Stokes-Einstein diffusivity is modified. The modification arises from the hydrodynamic interactions between suspended particles and/or particles in confined geometries. Central are the effects of the particle-generated disturbance velocity field on the motion of other particlesBatchelor (1976). For spherical particles, the influence of hydrodynamic interactions is known for a range of the particle volume fractions Russel, Saville, and Schowalter (1991). However, little is known about the diffusivity of anisotropic particles, in no small part due to the difficulty in modeling the inter-particle interactions.

This commonality of these hydrodynamic interactions naturally implies the existence of a swath of settings where they play out, ranging from confinement of single BrownianFaucheux and Libchaber (1994) and biopolymersBraun and Libchaber (2004), to suspensions in experimental cellsBrenner (1999) or proximity of particles due to phase changes of the solventWettlaufer (2019) and subsequent restricted mobilityPeppin, Spannuth, and Wettlaufer (2009). In this latter setting, freezing of colloidal suspensions is an important phenomenon that influences many natural Taber (1929); Hallet (1990) and technological Mazur (1970) processes. The phase-change boundary in a freezing colloidal suspension can become morphologically unstable due to the constitutional supercooling of the suspension close to the boundary Peppin, Elliott, and Worster (2006); Peppin, Worster, and Wettlaufer (2006); Peppin, Wettlaufer, and Worster (2008); You, Wang, and Worster (2018). The degree of constitutional supercooling is influenced by the diffusivity of particles in the suspension, which affects the transport of the particles away from the boundary. The diffusivity also plays a role in determining the thickness of ice lensesSchollick et al. (2016), and in the redistribution and aggregation of particles during the freezing of colloidal suspensions Spannuth et al. (2011).

I.1 Spherical suspensions

Prior to treating the diffusivity of spheroids, for context we briefly survey the effects of hydrodynamic interactions on the diffusivity of spherical particles. Dynamical light scattering experiments and the diffusivity of particles are intimately wed Rallison and Hinch (1986); Berne and Pecora (2000) and hence the key observations arise from this approach. The mean squared displacement of a single Brownian spherical particle of radius $a$ and mass $m^{s}$ suspended in a Newtonian fluid of viscosity $\mu$ grows diffusively for times, $t\gg t^{s}_{\text{m}}=\frac{m^{s}}{6\pi\mu a}$, where $t^{s}_{m}$ is the time scale for the relaxation of the momentum of the particle. The diffusivity of the particle is given by the Stokes-Einstein formula; $\mathbf{D}=\frac{k_{b}T}{6\pi\mu a}\mathbf{I}$, where $\mathbf{I}$ is the second-order isotropic identity tensor, $k_{b}$ is the Boltzmann constant, and $T$ the temperature of the fluid. BatchelorBatchelor (1976, 1983) calculated the effects of hydrodynamic interactions on the diffusivity of a dilute suspension of spherical particles, which act as a leading order correction in the particle volume fraction ($\phi_{p}$) to the Stokes-Einstein diffusivity in the Stokes limit–when the effect of inertial forces in the fluid are neglected. For any suspension, one can define two kinds of diffusivity; the gradient diffusivity, which quantifies the flux of particles down their concentration gradient, and the tracer diffusivity, which relates the mean squared displacement of a tracer particle in the suspension to the elapsed timeRallison and Hinch (1986). Batchelor showed that due to hydrodynamic interactions the gradient (tracer) diffusivity for the suspension increases (decreases) relative to the single particle diffusivity. Moreover, the tracer diffusivity takes on two values; one for short-times, $t^{s}_{m}\ll t\ll t^{s}_{c}=a^{2}\left(\frac{k_{b}T}{6\pi\mu a}\right)^{-1}$, and another for long-times, $t\gg t^{s}_{c}$, which is the time needed for a particle in the suspension to diffuse its own distance. Therefore, $t^{s}_{c}$ corresponds to the time scale at which there is a change in the configuration of the suspension. Batchelor’s results were subsequently verified experimentallyDe Kruif, Jansen, and Vrij (1987); Russel, Saville, and Schowalter (1991), using a Langevin approach Pusey and Tough (1983) and using Stokesian dynamics simulations Phillips, Brady, and Bossis (1988); Brady and Bossis (1988). In the non-dilute limit, that is with increasing particle volume fraction, the long-time tracer diffusivity decreases monotonically Brady (1994), whereas the gradient diffusivity does notRussel, Saville, and Schowalter (1991); Peppin, Worster, and Wettlaufer (2006).

I.2 Spheroidal suspensions

Here, we investigate the effects of hydrodynamic interactions on the diffusivity of anisotropic particles by modeling them as spheroids, which are axisymmetric ellipsoids, whose shapes are determined by their aspect ratio. The aspect ratio ($r$) is the ratio of a spheroid’s length along its axisymmetric axis to its maximum diameter perpendicular to the axis, and is thus a measure of its anisotropy. Spheroids with $r<1$ and $r>1$ are called oblate and prolate spheroids respectively; a sphere has $r=1$, a flat disk has $r\rightarrow 0$ and a slender fiber $r\rightarrow\infty$.

A single spheroid executes diffusive motion for times, $t\gg t_{m}=\max\left[\frac{m}{6\pi\mu L\,\min(X_{A},Y_{A})},\frac{I_{t}}{8\pi\mu L^{3}Y_{c}}\right]$, where $L$ is the length of the semi-major axis of its generating ellipse, $m$ is its mass, $I_{t}$ is its moment of inertia transverse to the axis, $X_{A}$, $Y_{A}$ and $Y_{C}$ are functions of the eccentricityKim and Karrila (1991), $6\pi\mu LX_{A}$ ($6\pi\mu LY_{A}$) is the translational drag along (perpendicular) to the axis, and $8\pi\mu L^{3}Y_{c}$ is the rotational drag torque transverse to the axis. The eccentricity of an oblate (prolate) spheroid is defined as $e=\sqrt{1-r^{2}}$ ($e=\sqrt{1-1/r^{2}}$), ranging from [math] to $1$ as the shape ranges from a sphere to a slender fiber or a flat disk. The largest of the translational and the angular momentum relaxation time scales is given by $t_{m}$ above. The translational ($\mathbf{D_{tr}}$) and rotational (${D_{r}}$) diffusivities of a neutrally buoyant spheroid suspended in a Newtonian fluid of dynamic viscosity $\mu$ are given by

[TABLE]

and

[TABLE]

respectively, and are shown as a function of the spheroid eccentricity in figures 1 and 2, where $D_{\parallel}^{s}=X_{A}^{-1}$, $D_{\perp}^{s}=Y_{A}^{-1}$ and $D_{r}^{s}=Y_{C}^{-1}$. The orientation of the spheroid is given by $\mathbf{p}_{1}$, and $\mathbf{D_{tr}}$ is anisotropic for $X_{A}\neq Y_{A}$. The translational diffusivities parallel and perpendicular to $\mathbf{p}_{1}$, reduce to the Stokes-Einstein diffusivity, when $e=0$. In the limit that $e\rightarrow 1$, the diffusivities diverge as $\log{(1-e)}$ for prolate spheroids, and asymptote to finite values ($D_{\perp}^{s}=1.77$ and $D_{\parallel}^{s}=1.18$) for oblate spheroids. Furthermore, when $e\neq 0$ the diffusivity parallel to the axis is larger (smaller) than that perpendicular to the axis for prolate (oblate) spheroids. Figure 2 shows that ${D_{r}}(e=0)$ = $\frac{k_{b}T}{8\pi\mu L^{3}}$, the spherical value with dimension $s^{-1}$, and ${D_{r}}(e\rightarrow 1)$ diverges as $\log{(1-e)}$ for the prolate spheroid, and asymptotes to a finite value ($D_{r}^{s}=2.35$) for the oblate spheroid. On a time scale $t\gg D_{r}^{-1}$, the translational diffusivity of a spheroid given in Eq. (1) will become isotropic owing to this rotational diffusion, and thus unlike spherical suspensions, even for a single spheroid, there are two translational diffusivities, one at short-times ($t\ll D_{r}^{-1}$) and the other at long-times ($t\gg D_{r}^{-1}$).

Hydrodynamic interactions insure that both the gradient and tracer diffusivities of a spheroidal suspension differ from those given in Eqs. (1) and (2). Importantly, the tracer diffusivity should take different forms: (a) For short-times, $t_{m}\ll t\ll t_{c}=\min\left[D_{r}^{-1},L^{2}\left(\frac{k_{b}T}{6\pi\mu L\min({X_{A},Y_{A}})}\right)^{-1}\right]$, and (b) for long-times, $t\gg\max\left[D_{r}^{-1},L^{2}\left(\frac{k_{b}T}{6\pi\mu L\min({X_{A},Y_{A}})}\right)^{-1}\right]$. Here, $t_{c}$ is the time scale at which the configuration of the suspension changes. Before describing our calculation, we survey previous investigations of the diffusivity of spheroidal particles in monodisperse suspensions.

Treloar and Masters Treloar and Masters (1991) calculated the effects of hydrodynamic interactions on the short-time translational diffusivity of nearly spherical spheroids ($e\ll 1$), and hypothesized that the diffusivity for small $e$, might be a good approximation for all values of $e$. Claeys and Brady Claeys and Brady (1993a) used Stokesian dynamics to calculate the short-time translational and rotational diffusivities of a prolate spheroid with aspect ratio $6$ in unbounded monodisperse suspensions, and found that the diffusivity decreased monotonically with increasing particle volume fraction. Zheng and HanZheng and Han (2010) measured the translational diffusivity of a tracer spheroid in a monolayer of spheroids ($r=9.07$) near a flat wall (a two-dimensional suspension) as a function of spheroid area fraction. They found that as the particle area fraction increases, the anisotropy of the diffusivity decreases, reaching a minimum. Upon further increase in area fraction the anisotropy of the diffusivity increased. However, their system is two-dimensional and their results cannot be translated to unbounded three-dimensional suspensions.

In this paper, we calculate the short-time translational and rotational diffusivities of a tracer spheroid of an arbitrary aspect ratio in a dilute monodisperse suspension. The correction to the single particle diffusivities appears at O($\phi$), where $\phi$ is the hydrodynamic volume fraction defined as $nL^{3}$, with $n$ the number density of the spheroids in the suspension. We use a far-field approximation to model the hydrodynamic interactions, and a Langevin approach to calculate the diffusivities. In section II, we describe the calculation of the Langevin equation for $N$ spheroids to obtain the mean squared angular and translational displacements of one spheroid, and the diffusivities follow from averaging over the configuration space of the other spheroids. The results are discussed in section III before concluding in section IV. To facilitate continuity in the presentation four Appendixes contain a number of technical details.

II Analysis

Here we calculate the short-time diffusivities of a tracer spheroid in a monodisperse suspension of $N$ spheroids using a Langevin equation approach as follows. The diffusivities are related to the mean squared rotational and angular displacements of the tracer spheroid at time $t$, with $t_{m}\ll t\ll t_{c}$. The mean squared displacements are obtained by solving the Langevin equation for $N$ spheroids (see Ref. (Ermak and McCammon, 1978) for the case of $N$ spheres). For each spheroid in the suspension, we define a Cartesian coordinate system $x_{i}$, $y_{i}$, $z_{i}$, with $i=1,\dots N$} such that the orientation of the spheroid, ${\mathbf{p}_{i}}$, is aligned with the $z_{i}$ axis. The coordinate system does not rotate with the spheroid. Noting that the moments of inertia of any spheroid along the transverse axes ($x_{i}$ and $y_{i}$) are the same, one can write down the Langevin equation in index notation as

[TABLE]

where the terms on the left-hand side arise due to the inertia of the spheroids, and the first and the second terms on the right-hand side are the deterministic and the stochastic forces exerted by the fluid on the spheroids respectively. In Eq. (3), we have neglected all the direct forces between the spheroids.

In Eq. (3), $\mathbf{M}$ is a generalized mass/moment of inertia diagonal matrix of dimension $6N\times 6N$ , with $M_{ii}=m$ for $i=1$ to $3N$, $M_{ii}=I_{t}$ for $i=3N+j,3N+2j$ and $M_{ii}=I_{s}$ for $i=3N+3j$, with $j=1$ to $N$. The constants $I_{s}$ and $I_{t}$ are the components of the moment of inertia of any of the $N$ spheroids along and transverse to its orientation respectively, and $m$ is the spheroid mass. Additionally, $\mathbf{X}$ and $\mathbf{g}$ are $6N$-dimensional vectors. The first vector contains the components of the translational ($v$) and the angular ($\omega$) velocities of the $N$ spheroids along the $x_{i},y_{i}$ and $z_{i}$ axes, and the second vector contains the components of the translational ($F$) and the rotational ($\tau$) Brownian noises acting on the $N$ spheroids along these axes, and are given by

[TABLE]

and

[TABLE]

In Eq. (3), the term proportional to $\mathbf{X}^{I}$ is a $6N$-dimensional vector that describes the rate of change of the moment of inertia of the spheroid along the $x_{i},y_{i}$ and $z_{i}$ axes, on account of the rotation of the spheroidGoldstein, Poole, and Safko (2014). We have ${X}^{I}_{i}=0$, for $i=1$ to $3N$, ${X}^{I}_{i}=-{\omega}_{y_{j}}{\omega}_{z_{j}}$, for $i=3N+j$, ${X}^{I}_{i}={\omega}_{x_{j}}{\omega}_{z_{j}}$, for $i=3N+2j$ and ${X}^{I}_{i}=0$, for $i=3N+3j$, with $j=1$ to $N$. Note that the rate of change vanishes along the symmetry axis of the spheroid, and for a sphere, $\mathbf{X}^{I}$ is a null vector. We assume that the noise terms have zero mean, are Gaussian and $\delta$-autocorrelated in time as

[TABLE]

The angular brackets above denote an ensemble average over the rapidly fluctuating random forces and torques for a given configuration of $N$ spheroids. The quantity $(k_{b}T)^{-1}\mathbf{D}$ is the mobility matrix for $N$ spheroids, relating their translational and the angular velocities to the hydrodynamic torques ($\tau^{h}$) and forces ($F^{h}$) acting on them as

[TABLE]

where

[TABLE]

and

[TABLE]

The two superscripts of a matrix element indicate whether it multiplies a force (‘$t$’) or a torque (‘$r$’) component, to give a translational (‘$t$’) or an angular (‘$r$’) velocity component in Eq. (5). The mobility matrix is only a function of the configuration–positions and orientations–of the $N$ spheroids. In Eq. (3), $\bm{\alpha}$ is the strength of the stochastic noise, which is related to the mobility matrix through the fluctuation-dissipation theorem given by

[TABLE]

which is derived in Appendix A. Therefore, the strength of the fluctuations depends on the configuration of the spheroids through Eq. (8), so that the positions and the orientations of the spheroids are influenced by multiplicative noise.

In Appendix B, Eq. (3) is solved using Eqs. (4) and (8) to obtain the mean squared translational and rotational displacements of a tracer spheroid (say 1), which are given by

[TABLE]

and

[TABLE]

In Eqs. (9) and (10), $\Delta r_{p}^{1}$ and $\Delta\theta_{p}^{1}$ are the translational and the rotational displacements of the tracer spheroid ‘1’ from its initial position/orientation at $t=0$, with $p,q={1,2\mbox{ and }3}$. The angular brackets denote the ensemble average over all possible configurations of the other spheroids in the suspension, with the tracer spheroid fixed in space.

The mean squared displacements in Eqs. (9) and (10) are valid even in the non-dilute limit. However, in the dilute limit, they can be simplified further by neglecting the simultaneous interactions between three or more spheroids, since these interactions lead to corrections that are less than or equal to O($\phi^{2}$), as seen in the case for spheres Dhont (1996). Because the mobility matrix is pairwise additive Durlofsky, Brady, and Bossis (1987), the averaging procedure reduces to that over the configuration space of two spheroids with the tracer spheroid fixed in space. In Appendix C we use the method of reflectionsKim (1985a); Kim and Karrila (1991) to derive the mobility matrix for two prolate spheroids under the influence of forces and torques, and we truncate the mobility matrix at the level of force dipoles, or stresslets, and neglect higher multipole moments. mportantly, these higher multipole moments become as important as the lower moments when the interparticle distances is small compared to the particle length, in which case lubrication forces become operative and can be crucial in non-dilute suspensions.

The method of reflections must include at least two reflections to insure that the velocity (angular velocity) of the first spheroid has a non-trivial relation to the force (the torque) acting on it. The trivial relation is given by Eq. (29) [Eq. (30)] for the translational [angular] velocity, and the non-trivial relation given by Eq. (LABEL:eq:finmob3) [Eq. (LABEL:eq:finmob4)]. The averaging of the trivial relations over the configuration space recovers the single particle diffusivities given in Eqs. (1) and (2). For averaging, we use only the leading order term in the non-trivial relation, which decays as $1/R_{p}^{4}$ ($1/R_{p}^{6}$) for translational (angular) velocity, where $R_{p}$ is the interparticle distance (non-dimensionalized by ‘$L$’) between the centers of the spheroids.

The mobility matrix for two oblate spheroids is obtained from that of two prolate spheroids by using the transformation described in Appendix D. The mobility matrices are four-dimensional integrals, and the averaging in Eqs. (9) and (10) requires an additional integration over a five-dimensional space; two for the orientation of the second spheroid and three for the interparticle vector connecting the centers of the spheroids. The nine-dimensional integrals required to average the diffusivities are given by

[TABLE]

and

[TABLE]

In these integrals, all the configurations in which the second spheroid intersects the tracer spheroid are avoided, and $V$ corresponds to that volume in which one can place the second spheroid without intersecting the tracer spheroid. We non-dimensionalize all distances with ‘$L$’, and the inter-particle separation vector $\mathbf{R_{p}}$ is defined as ($R_{p}$,$\theta_{R}$,$\phi_{R}$) in a spherical coordinate system whose zenith direction is aligned with the orientation of the tracer spheroid. Because we avoid intersections described above, the inner limit of $R_{p}$, defined as $R_{\text{restrict}}$, will depend on $\theta_{R}$ and $\phi_{R}$. The orientation of the second spheroid is defined in the spherical coordinate system as ($1$,$\theta_{2}$,$\phi_{2}$). We make a further simplification by setting $R_{\text{restrict}}=2$, which insures that the spheroids don’t intersect. The nine-dimensional integrals are evaluated numerically using a Monte Carlo method with a sampling size of $10^{9}$ configurations. Additionally, we evaluated the integrals using the NIntegrate function in Mathematica as a cross check.

III Results

The short-time translational ($\mathbf{D^{f}_{tr}}$) and rotational (${D^{f}_{r}}$) diffusivities of the tracer spheroid that are obtained by numerically integrating the right-hand sides of Eqs. (11) and (12) are written as

[TABLE]

and

[TABLE]

in which $\mathbf{D_{tr}}$ and ${D_{r}}$ are the single particle diffusivities given in Eqs. (1) and (2). In Eq. (13), $D_{\parallel}$ and $D_{\perp}$ are the O($\phi$) correction to the translational diffusivities parallel and perpendicular to the spheroid respectively, and $D_{r}^{\text{hyd}}$ in Eq. (14) is the O($\phi$) correction to the rotational diffusivity.

The corrections $D_{\parallel}$ and $D_{\perp}$ are plotted as a function of the eccentricity of prolate and oblate spheroids in figure 3. Because these corrections are negative for all eccentricities, the hydrodynamic interactions reduce the translational diffusivity of the tracer spheroid. Physically, this reduction is due to the hydrodynamic hindrance created by the velocity fields of the other particles in the suspension, which on average suppress the displacement of the tracer spheroid. Spherical symmetry is exhibited in figures 3 (a) and (b), where for $e=0$, $D_{\parallel}$ and $D_{\perp}$ take the same value of $-1.87$ 111BatchelorBatchelor (1976) obtained a slightly different value of $-1.83$ for the O($\phi$) correction, by using a mobility matrix which was exact for any arbitrary separation between two spheres, whereas we use a mobility matrix which is approximate, since only two reflections are considered. However, even with the approximate mobility matrix, we obtain a good estimate for the correction at $e=0$.. Both the corrections decrease monotonically with increasing eccentricity because (a) the interactions weaken for a given $\phi$($nL^{3}=$ constant), and (b) in the limit of $e\rightarrow 1$, vanish as $\propto 1/\log(1-e)$ for prolate spheroids, and for oblate spheroids, $D_{\parallel}$ and $D_{\perp}$ asymptote to $-0.34$ and $-0.35$ respectively. Physically, in this limit the hydrodynamic interactions vanish for slender fibers whereas they survive for flat disks due to the finite area of the latter. For nearly spherical prolate spheroids ($e\ll 1$), $D_{\parallel}$ and $D_{\perp}$ asymptote to $-1.87+2.00\,e^{2}$ and $-1.87+2.06\,e^{2}$ respectively, and for nearly spherical oblate spheroids, $D_{\parallel}$ and $D_{\perp}$ asymptote to $-1.87+0.81e^{2}$ and $-1.87+0.75\,e^{2}$ respectively. Thus, we find that, as hypothesized by Treloar and Masters Treloar and Masters (1991), the nearly spherical asymptotes are not a uniformly valid approximation for all eccentricities.

In figures 4 (a) and 5 (a), the hydrodynamic corrections ($\phi D_{\parallel}$ and $\phi D_{\perp}$) for prolate and oblate spheroids are expressed as a percentage of their associated single particle diffusivities ($D_{\parallel}^{s}$ and $D_{\perp}^{s}$), and plotted as a function of their eccentricities for $\phi=0.1$. Clearly, the largest hydrodynamic reduction in diffusivity occurs for spheres, and for both prolate and oblate spheroids the reductions decrease monotonically as $e\rightarrow 1$. We see that for eccentricities greater than approximately $0.4$, there is a small difference between the reductions in the parallel and perpendicular diffusivities; for prolate (oblate) spheroids, the reduction is slightly larger for the perpendicular (parallel) diffusivity than for the parallel (perpendicular) diffusivity. The physical origin of the difference can be described using the case of the prolate spheroid. The hydrodynamic effect of a moving spheroid in the suspension is to induce stresslets (dipoles) on the other spheroids, which in turn induce the same back upon the spheroid. On average, because the gradient in the reflected velocity field is likely to be larger (smaller) parallel to (perpendicular to) the spheroidal axis, so too is the strength of the induced dipole222The net change in the reflected velocity field along a spheroidal axis is the integral over the spheroid of the gradient of the velocity field along that direction, and it is averaged over all possible configurations of the suspension with the test spheroid fixed. This underlies the strength of the induced dipole.. Now, the radial dipole flow field is towards (away from) the spheroid perpendicular to (parallel to) the axis Kim and Karrila (1991), thereby suppressing (enhancing) perpendicular (parallel) motion, as seen in figure 4 (a). This logic is immediately extended to the oblate spheroid, in which case the inward radial flow due to the induced dipole will be parallel to the spheroidal axis. Finally, as the eccentricity increases the anisotropy of the dipole strength increases, whereas the strength decreases, due to the reduction in the size of the spheroid, and the latter effect dominates the behavior as $e\rightarrow 1$.

To characterize the effects of the anistropy on $\mathbf{D^{f}_{tr}}$, we define a ratio of the diffusivities, $D_{a}$, as

[TABLE]

which is plotted for prolate (oblate) spheroids versus eccentricity for various volume fractions in figures 4 (b) and 5 (b). Clearly, there is a weak $\phi$-dependence, exhibited differently in figures 4 (c) and 5 (c) by scaling $D_{a}$ with its value at $\phi=0$. Here we see that the maximum change in the ratio for prolate (oblate) spheroids is about $0.8\%$ ($2.5\%$) for the the largest volume fraction considered. Therefore, at O($\phi$), the hydrodynamic interactions do not significantly affect the anisotropy of the diffusivity tensor.

In figure 6, the orientational diffusivity correction $D_{r}^{\text{hyd}}$ of Eq. (14) is plotted versus the eccentricity of prolate and oblate spheroids. The correction is negative as in the translational case, and starts at the sphere value of $-0.31$ at $e=0$. The orientational diffusion does not affect the dynamics of a spherical particle due to its isotropy, and $D_{r}^{\text{hyd}}$ has a finite value just as in the single particle case (see figure 2). For prolate spheroids, in the limit of $e\rightarrow 1$, the correction decays as $1/\log(1-e)$ because of the weakened hydrodynamic interactions, and for oblate spheroids in the same limit the correction asymptotes to $-0.04$. In figure 7, the correction ($\phi D_{r}^{\text{hyd}}$) for a spheroid is expressed as a percentage of the associated single particle diffusivity ($D_{r}^{s}$), and plotted versus eccentricity. The percentage reduction in the orientational diffusivities of prolate and oblate spheroids are far smaller than that found for their translational counterparts in figures 4 (a) and 5 (a).

Claeys and BradyClaeys and Brady (1993a) used Stokesian dynamics to calculate the short-time diffusivities of a tracer spheroid of aspect ratio $6$ ($e=0.98$) in a monodisperse suspension. They modeled the hydrodynamic interactions by including lubrication forces, in addition to a far-field approximation of the mobility matrix. When averaging the diffusivity, they generated different configurations of $N$ spheroids (a maximum $N=64$) in a unit cell, with periodic boundary conditions, using a Monte Carlo procedure. They reported a linear combination of the translational diffusivities parallel and perpendicular to the symmetry axis of the tracer spheroid, which decreased with increasing volume fraction. In particular, for small volume fractions, their result is in excellent agreement with ours for the linear combination; $2.51-0.07\phi$ (when scaled with $\frac{k_{b}T}{6\pi\mu L}$) from figures 1 (a) and 3 (a).

IV Conclusions

We have calculated the hydrodynamic corrections to the short-time translational and rotational diffusivities of a tracer spheroid of arbitrary aspect ratio in a dilute monodisperse suspension to O($\phi$). The corrections reduce the diffusivities of the spheroid from that of an isolated spheroid of the same aspect ratio. The corrections are monotonically decreasing functions of the spheroid eccentricity, peaking for a sphere, vanishing for a slender fiber, and asymptoting to finite values for a flat disk. Moreover, the corrections do not significantly alter the anisotropy of the translational diffusion.

Further theoretical, numerical and experimental studies are needed to understand the behavior at higher particle volume fractions, where the lubrication effects neglected here become important. As $\phi$ increases, intuition suggests that the short-time diffusivities will be further reduced, since the reduction arises from the hydrodynamic hindrance created by the other particles in the suspension.

When one solidifies a colloidal suspension several processes relevant to the work described here are of interest. Firstly, as noted in §I, the solid/suspension interface can become morphologically unstable owing to constitutional supercooling, the degree of which is controlled by the diffusivity of particles in the suspension. Thus, our finding here of the suppression of the diffusivity implies that in natural settings, wherein the colloids are rarely spherical, morphological instability will be enhanced. Moreover, because one observes a van’t Hoff law in colloidal suspensionPeppin, Wettlaufer, and Worster (2008), our finding provides an underpinning for studies of its origin. Secondly, once particles are trapped in ice they are surrounded by a thin unfrozen (“premelted”) water film that allows for their subsequent motion by thermomolecular pressure gradientsWettlaufer (2019). The mobility of those particles can be treated within the framework of Onsager reciprocityPeppin, Spannuth, and Wettlaufer (2009), suggesting a range of interesting generalizations using spheroids, and in particular the question of thermodynamic buoyancyRempel, Wettlaufer, and Worster (2001). Moreover, because of the temperature dependence of the viscosity of supercooled water Dehaoui, Issenmann, and Caupin (2015), these dynamical processes will be influenced in a variety of ways, one of which concerns the temperature dependence of the ordering effect of water in confined supercooled filmsPramanik and Wettlaufer (2017).

Acknowledgements.

The authors acknowledge the support of the Swedish Research Council grant no. 638-2013-9243, conversation with their colleagues in Stockholm and Cristobal Arratia for access to his workstation on which the numerical integration code was executed.

Appendix A Derivation of the fluctuation dissipation theorem

In this section, we derive the fluctuation-dissipation theorem (FDT) for a suspension of $N$ spheroids. One starts with Eq. (3), albeit with an unknown relationship between the strength of the fluctuations and the mobility matrix, and noting that the time differentiation and ensemble average (over rapidly fluctuating random forces and torques for a given configuration of $N$ spheroids) commute, one can write,

[TABLE]

As opposed to Eq. (4), Eq. (16) does not hold when the noise is $\delta$-autocorrelated because $\langle\mathbf{X}(t)\mathbf{\dot{X}}(t)\rangle$ is ill-defined, although it does hold when the noise has a finite correlation time, a point to which we return in this derivation.

For $t>t_{m}$, the angular and translational velocities of the spheroids equilibrate, and hence $\frac{d}{dt}\langle\mathbf{X}(t)\mathbf{X}(t)\rangle=0$. Contracting Eq. (3) with ${M}^{-1}_{li}X_{f}$ and averaging gives,

[TABLE]

Using Eq. (16), the first term of Eq. (17) is zero, and the second term vanishes as it is an odd-moment in $\mathbf{X}$, and thus

[TABLE]

From the equipartition theorem Dhont (1996), the average in the first term above is equal to $k_{b}TM_{jf}^{-1}$. The second term is evaluated by noting that $\mathbf{X}$ is a functional of $\mathbf{g}$, namely $\mathbf{X}={{F}}[\mathbf{g}]$. Therefore, $\langle g_{j}(t)X_{f}(t)\rangle$=$\langle g_{j}(t)F_{f}[\mathbf{g}]\rangle$ and because $\mathbf{g}$ is Gaussian noise, one findsZinn-Justin (1996),

[TABLE]

The first average in the integrand above is given by Eq. (4), and the second average is evaluated by solving Eq. (3) to obtain the velocity vector (${\mathbf{X}}$), and then differentiating the vector with respect to $\mathbf{g}$. The average is given as

[TABLE]

where $H(t)$ is the Heaviside function. The $\delta$-autocorrelation of the noise in the integrand of Eq. (19) implies taking the limit of $t^{\prime}\rightarrow t$ in Eq. (20). In this limit, the contribution to Eq. (19) from the first term in Eq. (20) vanishes, and hence Eq. (18) becomes

[TABLE]

The Heaviside function at the origin, $H(0)$, is however discontinuous (often defined as the integral of a Dirac $\delta$-function) and arises from the fact that $\langle\mathbf{X}(t)\mathbf{\dot{X}}(t)\rangle$ is ill-defined for $\delta$-autocorrelated noise. A solution to this problem is to instead use noise that has a finite autocorrelation time, viz., $\langle g_{i}(t)g_{j}(t^{\prime})\rangle=2\delta_{ij}\,\eta(t-t^{\prime})$, where $\eta$ is an even function about $t=0$ such that $\int_{-\infty}^{\infty}\eta(t)dt=1$, and repeat the steps in this section up to Eq. (21). The net effect of using this colored noise is equivalent to replacing $H(0)$ with $\frac{1}{2}$ in Eq. (21), which leads to the FDT given in Eq. (8). Note that, although we have derived the FDT for $N$ spheroids, this derivation can easily be extended to $N$ rigid bodies.

Appendix B Derivation of the mean squared displacements

In this section, starting from the Langevin equation (Eq. 3), we derive the translational and angular mean squared displacements of a tracer spheroid in a suspension of $N$ spheroids at short-times ($t_{m}\ll t\ll t_{c}$). The strength of the noise term in Eq. (3) depends on the configuration of the particle described by Eq. (8). Hence, the governing equation for the positions and the orientations of the spheroids will be driven by a multiplicative noise. Although, such a noise leads to the well-known Itô-Stratonovich dilemmaMoon and Wettlaufer (2014), the mean squared displacements at short-times ($\propto t$) required here, are unaffected by the calculus that is used. Eq. (3) can be rewritten as

[TABLE]

Defining $Q_{lj}=k_{b}T\,M^{-1}_{li}D_{ij}^{-1}$, and noting that for times such that $t_{m}\ll t\ll t_{c}$, $\mathbf{Q}$ is independent of time, Eq. (22) can be written as

[TABLE]

One can integrate Eq. (23) twice with respect to time to find the translational and angular displacements of any spheroid. Since finite angular displacements are not vectorsGoldstein, Poole, and Safko (2014), integrating the angular velocities with respect to time results in the correct angular displacements only when the displacements are infinitesimal, which is true here since $t\ll t_{c}$. The final translational displacement ($\Delta\mathbf{r}$) of a tracer spheroid (say 1) in the suspension is given as

[TABLE]

where the integral containing $\mathbf{g}$ leads to a mean squared displacement that is proportional to $t$. The other two integrals lead to mean squared displacements that are smaller for short times, and hence are not considered. Thus, after neglecting those two integrals, the final mean squared displacement of the tracer spheroid is obtained by averaging the outer product $\Delta r_{p}^{1}\Delta r_{q}^{1}$, first over an ensemble of realizations of random forces and torques, and then over an ensemble of various configurations of the spheroids in which the tracer particle is fixed in space. Using Eqs. (4) and (8), one finds

[TABLE]

for the translational mean squared displacement of the tracer spheroid. The angular mean squared displacement of the spheroid is obtained similarly and is given by

[TABLE]

Appendix C Derivation of the mobility matrix for two prolate spheroids

Using the method of reflectionsKim (1985a); Kim and Karrila (1991), we derive the mobility matrix for two prolate spheroids that are subjected to forces as well as torques. Consider two identical prolate spheroids that are subjected to torques and forces as shown in figure 8. The two spheroids will disturb the flow around them, and thereby interact with each other hydrodynamically. The aim is to calculate the relationship between the translational and the angular velocities of the spheroids and the forces and the torques. In the present calculation, the method of reflections has to be carried out up to the second reflection to obtain the leading order effects of hydrodynamic interactions. The translational ($\mathbf{U}_{1}$) and the angular velocities ($\bm{\omega}_{1}$) of the first spheroid (with orientation $\mathbf{p}_{1}$ in figure 8) are obtained by summing the contributions from different reflections that are given as

[TABLE]

and

[TABLE]

where the superscripts on the right-hand side indicate the number of reflections involved in obtaining that term.

At the zeroth reflection, the velocity field generated by a given spheroid in the figure is obtained using the method of singularities Chwang and Wu (1974, 1975). The translational and the angular velocities of the first spheroid as well as the disturbance velocity field ($\mathbf{v}_{1}^{\prime}$) at a point $\mathbf{x}$ generated by it Kim (1985a); Claeys and Brady (1993b) are given as

[TABLE]

[TABLE]

and

[TABLE]

Here, $c$ and $e$ are the semi interfocal distance and the eccentricity of the spheroid respectively, $\bm{\xi}_{1}={\xi}_{1}\mathbf{p}_{1}$ and $\mathbf{J}(\mathbf{x}-\bm{\xi}_{1})=\frac{I}{|\mathbf{x}-\mathbf{\xi}_{1}|}+\frac{(\mathbf{x}-\mathbf{\xi}_{1})(\mathbf{x}-\mathbf{\xi}_{1})}{|\mathbf{x}-\mathbf{\xi}_{1}|^{3}}$ is the Oseen tensorKim and Karrila (1991). The quantities $X_{A}$, $Y_{A}$, $X_{C}$ and $Y_{C}$ are functions of the eccentricity Kim and Karrila (1991) and are given at the end of this Appendix. The integrands in Eq. (LABEL:eq:mobility1) are evaluated at all points between the foci from $\bm{\xi}_{1}=-c\mathbf{p}_{1}$ to $c\mathbf{p}_{1}$. Eq. (LABEL:eq:mobility1) implies that the velocity field due to the spheroid acted upon by a torque and a force is equivalent to that generated by a line distribution of singularities (Stokeslets, dipoles and degenerate quadrapoles) along the aforementioned points.

Since the force and the torque on the second spheroid (with orientation $\mathbf{p}_{2}$ in figure 8) are specified, $\mathbf{v}_{1}^{\prime}$ cannot exert any additional force and torque upon it. Therefore, the spheroid will rotate and translate such that the additional force as well as the torque are zero. However, $\mathbf{v}_{1}^{\prime}$ will induce a stresslet, $\mathbf{S}_{2}^{(1)}$, on the spheroid that will in turn lead to a reflected disturbance velocity field, $\mathbf{v}_{12}^{\prime}$. The angular velocity, $\mathbf{\omega}_{2}^{(1)}$, and the stresslet are obtained using Faxen’s law Kim (1985b), and are given by

[TABLE]

and

[TABLE]

where $\bm{\xi}_{2}={\xi}_{2}\mathbf{p}_{2}$ changes from $-c\mathbf{p}_{2}$ to $c\mathbf{p}_{2}$ along the symmetry axis of the second spheroid. The tensor $\mathbf{e}_{1}^{\prime}(\bm{\xi}_{2})$ is the strain rate tensor due to $\mathbf{v}_{1}^{\prime}(\bm{\xi}_{2})$. The fourth order tensor $\mathbf{W}^{(2)}_{4}$ and the third order tensor $\mathbf{Y}^{(2)}_{3}$ are given byKim and Karrila (1991)

[TABLE]

and

[TABLE]

where

[TABLE]

and

[TABLE]

The quantities $X_{M}$, $Y_{M}$, $Z_{M}$ and $Y_{H}$ are functions of the eccentricity and are given at the end of the section. The disturbance velocity field reflected by the second spheroid is given by

[TABLE]

where $\bm{\xi}_{2}^{\prime}={\xi}_{2}^{\prime}\mathbf{p}_{2}$ changes from $-c\mathbf{p}_{2}$ to $c\mathbf{p}_{2}$ along the symmetry axis of the second spheroid and the velocity field is proportional to $\mathbf{F}_{1}$ and $\mathbf{T}_{1}$, through $\mathbf{S}_{2}^{(1)}$. The velocity field cannot cause any additional torque and force on the first spheroid, and therefore the spheroid has to move with the translational and the angular velocities given by

[TABLE]

and

[TABLE]

where $\bm{\xi}_{1}^{\prime}={\xi}_{1}^{\prime}\mathbf{p}_{1}$ changes from $-c\mathbf{p}_{1}$ to $c\mathbf{p}_{1}$ along the symmetry axis of the first spheroid, and $\mathbf{e}_{12}^{\prime}$ is the strain rate tensor due to $\mathbf{v}_{12}^{\prime}$. The Eqs. (29), (30), (LABEL:eq:finmob3) and (LABEL:eq:finmob4) are substituted into Eqs.  (27)and (28) to obtain the elements of the mobility matrix that multiplies $\mathbf{F}_{1}$ and $\mathbf{T}_{1}$ to give $\mathbf{U}_{1}$ and $\bm{\omega}_{1}$ respectively. The velocities from the first reflection, $\mathbf{U}_{1}^{(1)}$ and $\bm{\omega}_{1}^{(1)}$ are not needed in the present calculation, as they are proportional to $\mathbf{F}_{2}$ and $\mathbf{T}_{2}$. Therefore, for a given configuration of the two spheroids in figure 8, the elements of the mobility matrix are 4-dimensional integrals over $\xi_{1}$, $\xi_{2}$, $\xi_{1}^{\prime}$ and $\xi_{2}^{\prime}$.

The quantities needed to evaluate the mobility matrix are given by Kim and Karrila (1991)

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

Appendix D Mobility matrix for two oblate spheroids

Here we derive the mobility matrix for two oblate spheroids. The velocity field due to an oblate spheroid acted upon by a force and a torque is equivalent to that generated by a surface distribution of singularities (Stokeslets, dipoles and degenerate quadrapoles), the surface being a circle whose diameter is equal to the inter-focal length of the spheroid, and is on the equatorial plane of the spheroid with its centre being the spheroid’s centre of massKim and Karrila (1991). While one could derive the mobility matrix using the aforementioned surface distribution, in order to simplify the diffusivity calculation, the matrix is obtained by transforming that of the prolate spheroid derived in appendix C. The transformation involves replacing the semi-inter-focal distance $c$ in the final expression for the velocities given in Eqs. (27) and (28) (after substituting for the velocities from Eqs. (29), (30), (LABEL:eq:finmob3) and (LABEL:eq:finmob4)) with $-i\,c$, and replacing the eccentricity $e$ with $-ie/\sqrt{1-e^{2}}$, where the latter substitution means that the quantities $X_{A},Y_{A},Y_{c},X_{m},Y_{m},Y_{c},Y_{H}$ will be replaced by their oblate counterparts given at the end of this Appendix. The arguments leading to the transformation is given below.

The mobility matrix for the prolate spheroid was derived in Appendix C using the singularity representation and Faxen’s law for prolate spheroids. It is already known that the transformation, when applied to the singularity representation of the velocity field generated by a prolate spheroid of eccentricity $e$ in a particular flow (rotation, translation, straining, etc.,), gives the velocity field generated by an oblate spheroid of the same eccentricity in the same flowShatz (2004); Dabade, Marath, and Subramanian (2015, 2016). Now, the singularity representation can be interpreted as follows; the velocity field generated by a spheroid is obtained by contracting a translational velocity (or angular velocity or strain rate) with a tensor obtained by operating a functional on $\frac{\mathbf{J}}{8\pi\mu}$, where $\mathbf{J}$ is the Oseen tensor. For example, from Eq. (C) the functional for translation would be $[6\pi\mu X_{A}\mathbf{p}_{1}\mathbf{p}_{1}+6\pi\mu Y_{A}(\mathbf{I}-\mathbf{p}_{1}\mathbf{p}_{1})]({2c})^{-1}\displaystyle\int_{-c}^{c}d\xi\left[1+(c^{2}-\xi^{2})(1-e^{2})/4e^{2}\nabla^{2}\right]$ , where the term in the first square bracket is the resistance tensor for a translating prolate spheroid. Applying the transformation on the singularity representation of the prolate spheroid is equivalent to applying it on this functional. The Lorentz reciprocal theorem can be used to show that the same functional, when used to operate on an ambient flow around a prolate spheroid, gives the driving force (or torque or stresslet)Kim and Karrila (1991) due to the flow on the spheroid, leading to the Faxen law for prolate spheroidsKim (1985b). Therefore, Faxen’s law for oblate spheroids can be obtained from that of prolate spheroids by applying the transformation. Hence, by applying this transformation to the mobility matrix of a prolate spheroid of a given eccentricity, which was derived based on Faxen’s law and the singularity representation, one can obtain that of an oblate spheroid of the same eccentricity.

The quantities needed to evaluate the mobility matrix for oblate spheroids are given by Kim and Karrila (1991).

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]