1. Introduction
Transport and mixing of solutes or particles in the presence of hydrodynamic flows are important for various biological and industrial processes. For micron-sized particles immersed in flows, the coupling between advection and diffusion can often lead to enhanced mass transport as compared with diffusion without flow. A classical example of such a coupling effect is the Taylor dispersion of Brownian solutes in pressure-driven channel flows (Taylor Reference Taylor1953, Reference Taylor1954a,Reference Taylorb; Aris Reference Aris1956). Brownian motion allows the solute particles to migrate across streamlines and then be advected downstream with different velocities. At long times, the coupling between transverse diffusion and longitudinal advection gives rise to diffusive transport of the solutes with an effective longitudinal dispersion coefficient that can be much larger than the bare diffusivity of the solute particle. Since the work of Taylor (Reference Taylor1953), a generalized Taylor dispersion (GTD) framework (Frankel & Brenner Reference Frankel and Brenner1989) has been developed to accommodate a wide range of transport problems including complex geometries, chemical reactions, spatial and/or time periodicity and active particles (Brenner Reference Brenner1980; Shapiro & Brenner Reference Shapiro and Brenner1986, Reference Shapiro and Brenner1987, Reference Shapiro and Brenner1990; Morris & Brady Reference Morris and Brady1996; Hill & Bees Reference Hill and Bees2002; Bearon Reference Bearon2003; Manela & Frankel Reference Manela and Frankel2003; Zia & Brady Reference Zia and Brady2010; Takatori & Brady Reference Takatori and Brady2014; Burkholder & Brady Reference Burkholder and Brady2017; Alonso-Matilla, Chakrabarti & Saintillan Reference Alonso-Matilla, Chakrabarti and Saintillan2019; Jiang & Chen Reference Jiang and Chen2019; Peng & Brady Reference Peng and Brady2020). More recently, longitudinal dispersion of elongated Brownian rods in a two-dimensional Poiseuille flow has been considered (Kumar et al. Reference Kumar, Thomson, Powers and Harris2021; Khair Reference Khair2022).
In contrast to the extensive study of the long-time effective transport of particles in position (linear) space, the dynamics and transport of particles in orientation space remains relatively less developed. For spherical or ‘point’ particles, the consideration of the orientational dynamics is often unnecessary. For anisotropic particles, their orientational dynamics plays a role in the overall dynamics and rheology of the suspension composed of the particles and the fluid (Leal & Hinch Reference Leal and Hinch1971; Hinch & Leal Reference Hinch and Leal1972; Khair Reference Khair2016). A typical example is the orientational dynamics of an isolated spheroid in simple shear. Under shear, the orientation of the spheroid undergoes complex dynamics governed by Jeffery's equation (Jeffery Reference Jeffery1922). As a result, a Brownian spheroid in simple shear experiences both rotational diffusion and angular advection that is non-uniform. An interesting question we wish to consider is: Does the coupling of advection and diffusion in orientation space lead to enhanced rotational transport?
Using experiments and particle-based simulations, Leahy et al. (Reference Leahy, Cheng, Ong, Liddell-Watson and Cohen2013) showed that advection–diffusion coupling indeed results in enhanced rotational diffusion at long times for an axisymmetric particle under shear. In a later paper (Leahy, Koch & Cohen Reference Leahy, Koch and Cohen2015), a continuum theory is developed to calculate the time-dependent orientation distribution for non-spherical axisymmetric particles confined to rotate in the velocity-gradient plane, in the limit of weak diffusion or large Péclet number (see § 3.1 for the definition). In this limit, a coordinate transformation is discovered and used to map the orientation dynamics to a diffusion equation, which ultimately allowed the calculation of the long-time rotational dispersion coefficient. Furthermore, a remarkably simple analytic expression is obtained for the dispersion coefficient in the large-Péclet-number limit. In comparison to the classical Taylor dispersion, Leahy et al. (Reference Leahy, Koch and Cohen2015) concluded that their theoretical consideration does not fit nicely under the canonical GTD framework.
In this work we show that the flow-enhanced rotational transport in the velocity-gradient plane can be treated as a GTD in orientation space. To setup the system for such a treatment, the key step is to consider the unbounded angular displacement 
$\varphi$ rather than the orientation angle 
$\phi$, which is bounded to an interval of length 
$2{\rm \pi}$. With this, one can then break down the unbounded displacement into an infinite sequence of cells, each of which has a length of 
$2{\rm \pi}$. In the language of GTD, one then identifies the cell index 
$j \in \mathbb {Z}$ as the global coordinate and 
$\phi$ as the local coordinate. The derived GTD theory works for generic linear flows and for arbitrary Péclet numbers. In the large-Péclet-number limit, we show that our asymptotic result agrees with that obtained by Leahy et al. (Reference Leahy, Koch and Cohen2015) for steady simple shear. Our results from the GTD theory is validated against Brownian dynamics (BD) simulations.
In § 2, starting from the Smoluchowski equation governing the orientational dynamics of a spheroidal particle, we develop the GTD formulation for generic linear flows. Following Leahy et al. (Reference Leahy, Koch and Cohen2015), the particle is constrained to rotate in the velocity-gradient plane. In § 3 we consider the long-time rotational transport in simple shear and extensional flows. The transport coefficients are calculated using perturbation expansions in both small- and large-Péclet-number limits. The results obtained in these asymptotic limits are compared with numerical solutions of the macrotransport equations and results from BD simulations. Lastly, we conclude the paper in § 4.
2. Problem formulation
2.1. The Smoluchowski equation
Consider a spheroidal particle immersed in a linear ambient flow field in an unbounded, incompressible Newtonian fluid. The particle is sufficiently small so that inertia effects are neglected and the fluid is in the Stokes regime. The particle is subject to rotational Brownian motion and no external torque is applied. In the absence of Brownian motion, the time evolution of the unit orientation vector 
$\boldsymbol {q}$ (
$|\boldsymbol {q}|=1$) of the particle is governed by Jeffery's equation (Jeffery Reference Jeffery1922):
\begin{equation} \frac{\text{d}\boldsymbol{q}}{\text{d} t} = \boldsymbol{\varOmega}\times \boldsymbol{q} \quad \text{and}\quad \boldsymbol{\varOmega} = \frac{1}{2}\boldsymbol{\omega} + B \boldsymbol{q} \times (\boldsymbol{E} \boldsymbol{\cdot} \boldsymbol{q}). \end{equation}
Here 
$\boldsymbol {\varOmega }$ is the instantaneous angular velocity, 
$\boldsymbol {\omega } = \boldsymbol {\nabla } \times \boldsymbol {u}$ is the vorticity vector, 
$\boldsymbol {E} = \frac {1}{2}( \boldsymbol {\nabla } \boldsymbol {u} + ( \boldsymbol {\nabla } \boldsymbol {u})^\intercal )$ is the rate-of-strain tensor, 
$\boldsymbol {u}$ is the ambient flow field and 
$B = (r^2-1)/(r^2+1) \in [0,1)$ is the Bretherton constant that characterizes the non-sphericity (Bretherton Reference Bretherton1962), with 
$r$ the aspect ratio of the spheroid. For a sphere, 
$r=1$ and 
$B=0$. For an infinitely thin rod, 
$r \to \infty$ and 
$B \to 1$.
With Brownian motion, a statistical mechanical description is required. To this end, we define the orientational probability density function 
$\varPsi (\boldsymbol {q}, t)$, which satisfies the total conservation condition 
$\int _{\mathbb {S}} \varPsi (\boldsymbol {q}, t) \text {d}\boldsymbol {q} = 1$ at (any) time 
$t$. Here, 
$\mathbb {S} = \{\boldsymbol {q}\, \lvert \, \boldsymbol {q}\boldsymbol {\cdot }\boldsymbol {q}=1 \}$ denotes the unit sphere of orientations. The orientational probability density function is governed by the Smoluchowski equation (Brenner & Condiff Reference Brenner and Condiff1974; Doi & Edwards Reference Doi and Edwards1988)
\begin{equation} \frac{\partial \varPsi }{\partial t} + \boldsymbol{\nabla}_R\boldsymbol{\cdot} \boldsymbol{j}_R =0, \end{equation}
where 
$\boldsymbol {j}_R = \boldsymbol {\varOmega } \varPsi - D_R \boldsymbol {\nabla }_R \varPsi$ is the rotational flux vector, 
$\boldsymbol {\nabla }_R = \boldsymbol {q} \times \partial /\partial \boldsymbol {q}$ is the rotational gradient operator and 
$D_R$ is the rotational diffusivity.
We note that (2.2) can be treated as a marginalization of the full probability density function 
$P(\boldsymbol {x}, \boldsymbol {q}, t)$ that describes the joint distribution of the particle in both position and orientation space, where 
$\boldsymbol {x}$ is the position vector of the particle centre. This full probability is governed by
\begin{equation} \frac{\partial P }{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{j}_T + \boldsymbol{\nabla}_R\boldsymbol{\cdot} \boldsymbol{j}^\prime_R =0, \end{equation}
where, for a Brownian particle, 
$\boldsymbol {j}_T =\boldsymbol {u} P - \boldsymbol {D}_T(\boldsymbol {q}) \boldsymbol {\cdot }\boldsymbol {\nabla } P$ and 
$\boldsymbol {j}^\prime _R = \boldsymbol {\varOmega } P - D_R \boldsymbol {\nabla } P$. Here, 
$\boldsymbol {D}_T$ is the translational diffusivity of the particle, which is a function of 
$\boldsymbol {q}$ for a spheroid. It is clear that 
$\varPsi (\boldsymbol {q},t) = \int P(\boldsymbol {x}, \boldsymbol {q}, t)\text {d}\kern0.06em \boldsymbol {x}$. For active (self-propelled) Brownian particles with a constant swim speed 
$U_s$, an additional term 
$U_s\boldsymbol {q} P$ would appear in the translational flux 
$\boldsymbol {j}_T$. However, this would not affect the resulting equation for 
$\varPsi$. In fact, any advective linear velocity is allowed provided that the translational flux vanishes at infinity. A difficulty in the marginalization would appear if the angular velocity 
$\boldsymbol {\varOmega }$ or the rotational diffusivity 
$D_R$ depends on 
$\boldsymbol {x}$. For linear flows as we consider here, the angular velocity is spatially uniform.
2.2. Rotational Taylor dispersion theory
It is cumbersome to work with the unit orientation vector 
$\boldsymbol {q}$ in the consideration of the long-time rotational dispersion because 
$\boldsymbol {q}$ is bounded to the unit sphere (Kämmerer, Kob & Schilling Reference Kämmerer, Kob and Schilling1997; De Michele & Leporini Reference De Michele and Leporini2001; Mazza et al. Reference Mazza, Giovambattista, Starr and Stanley2006, Reference Mazza, Giovambattista, Stanley and Starr2007; Hunter et al. Reference Hunter, Edmond, Elsesser and Weeks2011). From a micromechanical perspective in considering the stochastic trajectory of a particle, one needs to be able to track the unbounded or cumulative angular displacement. The particle orientation vector is constrained to rotate in the velocity-gradient plane (Leahy et al. Reference Leahy, Cheng, Ong, Liddell-Watson and Cohen2013, Reference Leahy, Koch and Cohen2015). In two dimensions, the orientation vector can be parameterized as 
$\boldsymbol {q} = \cos \phi \boldsymbol {e}_x + \sin \phi \boldsymbol {e}_y$, where 
$\boldsymbol {e}_x$ and 
$\boldsymbol {e}_y$ are the unit basis vectors of the Cartesian coordinate system 
$(x,y)$, and 
$\phi \in [0,2{\rm \pi} )$ is the orientation angle. The cumulative angular displacement 
$\varphi$ that is not bounded to the interval 
$[0, 2{\rm \pi} )$ can be defined via
\begin{equation} \varphi = 2 {\rm \pi}j + \phi, \end{equation}
where 
$j \in \mathbb {Z}$. Conversely, the bounded orientation angle 
$\phi$ is 
$\varphi$ modulo 
$2{\rm \pi}$. For a constant angular velocity 
$\boldsymbol {\varOmega } = \varOmega \boldsymbol {e}_z$ with 
$\boldsymbol {e}_z=\boldsymbol {e}_x\times \boldsymbol {e}_y$, we have 
$\varphi (t) - \varphi (0) = \int _0^t \varOmega \,\text {d} s = \varOmega t$, where 
$\varphi (0)$ is a reference value.
We remark that alternative methods exist to quantify the rotational dynamics. In particular, one may extract a long-time dispersion coefficient from an orientational correlation function as a function of time in a BD simulation of the orientational Langevin equation of motion (Dhont Reference Dhont1996; Zwanzig Reference Zwanzig2001; Leahy et al. Reference Leahy, Cheng, Ong, Liddell-Watson and Cohen2013). One can also directly keep track of the unbounded angular displacement 
$\varphi$ in a BD simulation and infer the long-time transport coefficients (Kämmerer et al. Reference Kämmerer, Kob and Schilling1997; De Michele & Leporini Reference De Michele and Leporini2001; Mazza et al. Reference Mazza, Giovambattista, Starr and Stanley2006, Reference Mazza, Giovambattista, Stanley and Starr2007; Hunter et al. Reference Hunter, Edmond, Elsesser and Weeks2011). Because our aim is to derive a GTD theory from a continuum (Smoluchowski) perspective, such methods are not pursued here. In Leahy et al. (Reference Leahy, Koch and Cohen2015) a coordinate transformation is discovered and used to map the orientational dynamics to a diffusion equation, which ultimately leads to a closed-form asymptotic solution to the long-time rotational diffusivity in the high shear rate limit.
In terms of the bounded orientation angle 
$\phi$, the Smoluchowski equation (2.2) is written as
\begin{equation} \frac{\partial \varPsi }{\partial t} + \frac{\partial }{\partial \phi} \left(\varOmega(\phi) \varPsi - D_R \frac{\partial \varPsi }{\partial \phi}\right)=0, \end{equation}
where the angular velocity 
$\varOmega (\phi )$ depends on the orientation angle. It is clear that (2.5) remains unchanged in terms of the unbounded coordinate 
$\varphi$. Noting that 
$\varPsi = \varPsi (\varphi, t) = \varPsi (\,j, \phi, t)$, one can rewrite 
$\varPsi$ in terms of the sequence 
$\{\varPsi _j(\phi, t),\ \forall j \in \mathbb {Z}\}$. In other words, to locate the particle in the unbounded orientation space, one can first identify the cell index, 
$j$, in which the particle resides, and then use the local angular position 
$\phi$ within this cell. In the language of GTD, 
$\phi$ is identified as the local coordinate whereas the cell index 
$j$ is the global coordinate (Frankel & Brenner Reference Frankel and Brenner1989; Brenner & Edwards Reference Brenner and Edwards1993). In other words, the long-time diffusive dynamics holds only when the particles have traversed many cells.
Because the 
$\varphi$ space is unbounded, it is more convenient to work in Fourier space. In the following, we make use of the flux-averaging approach of Brady and coworkers (Morris & Brady Reference Morris and Brady1996; Zia & Brady Reference Zia and Brady2010; Takatori & Brady Reference Takatori and Brady2014, Reference Takatori and Brady2017; Burkholder & Brady Reference Burkholder and Brady2017; Peng & Brady Reference Peng and Brady2020). We note that the original approach was developed for the transport of particles in unbounded domains where the global coordinate is continuous (e.g. the longitudinal coordinate along a flat channel). Recently, Barakat & Takatori (Reference Barakat and Takatori2023) extended this approach to accommodate the transport and dispersion in an oscillating array of harmonic traps, where the global coordinate is the discrete cell index of the infinite lattice. (An equivalent real-space approach may be taken where one makes use of the method of moments; see, e.g. Brenner Reference Brenner1980; Alonso-Matilla et al. Reference Alonso-Matilla, Chakrabarti and Saintillan2019.) In the current problem, we have a one-dimensional lattice of unit cells. Following Barakat & Takatori (Reference Barakat and Takatori2023), we introduce the semi-discrete Fourier transform (Trefethen Reference Trefethen2000)
\begin{equation} \hat{f}(k) = 2{\rm \pi} \sum_{j =-\infty }^\infty \exp({- {\rm i}k j 2{\rm \pi}})\, f_j, \end{equation}
where 
$k$ is the wavenumber and 
$i$ (
$i^2=-1$) is the imaginary unit. Note that the transform is from 
$j$ to 
$k$, and the local coordinate 
$\phi$ is unchanged. In Fourier space, the Smoluchowski equation becomes
\begin{equation} \frac{\partial \hat{\varPsi}}{\partial t} + \left({\rm i}k + \frac{\partial }{\partial \phi}\right) \left[ \varOmega \hat{\varPsi} - D_R \left({\rm i}k + \frac{\partial }{\partial \phi} \right) \hat{\varPsi}\right]=0, \end{equation}
where 
$\hat {\varPsi } = \hat {\varPsi }(k, \phi, t)$ is the Fourier transform of 
$\varPsi$. We note that in (2.7) the local gradient is written in real space whereas the global gradient appears in powers of 
$k$. The cell-averaged distribution,
\begin{equation} \langle \hat{\varPsi} \rangle(k,t) = \frac{1}{2{\rm \pi}}\int_0^{2{\rm \pi}} \hat{\varPsi}(k, \phi, t) \,\text{d}\phi, \end{equation}satisfies
\begin{equation} \frac{\partial \langle \hat{\varPsi} \rangle }{\partial t} + {\rm i}k \langle \varOmega \hat{\varPsi} \rangle + D_R k^2 \langle \hat{\varPsi} \rangle=0, \end{equation}
which is obtained by averaging (2.7). In writing (2.9) we have invoked the periodic boundary condition on the local coordinate 
$\phi$ (Barakat & Takatori Reference Barakat and Takatori2023). One can relate 
$\hat {\varPsi }$ to its average by defining the structure function 
$\hat {G}$ such that 
$\hat {\varPsi }(k,\phi, t) = \langle \hat {\varPsi } \rangle (k,t) \hat {G}(k,\phi,t)$. The structure function is normalized, i.e.
\begin{equation} \langle \hat{G} \rangle =1. \end{equation} To derive an effective advection-diffusion equation for 
$\langle \hat {\varPsi } \rangle$, we first take a small wavenumber expansion of 
$\hat {G}$ (Morris & Brady Reference Morris and Brady1996; Zia & Brady Reference Zia and Brady2010), giving
\begin{equation} \hat{G}(k,\phi,t) = g(\phi, t) + {\rm i}k b(\phi, t)+O(k^2), \end{equation}
where 
$g$ is the average (zero wavenumber) field and 
$b$ is the displacement field. Inserting the expansion 
$\langle \varOmega \hat {\varPsi } \rangle = \langle \hat {\varPsi } \rangle \langle \varOmega (g + {\rm i}k b) \rangle +O(k^2)$ into (2.9), we obtain
$$\begin{gather} \frac{\partial \langle \hat{\varPsi} \rangle }{\partial t} + {\rm i}k \varOmega^{eff} \langle \hat{\varPsi} \rangle + k^2 D^{eff} \langle \hat{\varPsi} \rangle=0, \end{gather}$$
$$\begin{gather}\varOmega^{eff} = \langle \varOmega g \rangle,\quad D^{eff} = D_R - \langle \varOmega b \rangle. \end{gather}$$Note that (2.12a) is an effective advection–diffusion equation in Fourier space, with the effective rotational drift and rotational dispersion coefficient given in (2.12b).
Subtracting (2.9) multiplied by 
$\hat {G}$ from (2.7), we obtain
\begin{equation} \frac{\partial \hat{G}}{\partial t}+ {\rm i}k \left[\left(\varOmega - \frac{\langle \varOmega \hat{\varPsi} \rangle}{\langle \hat{\varPsi} \rangle}\right)\hat{G} -D_R \frac{\partial \hat{G} }{\partial \phi}\right] + \frac{\partial}{\partial \phi}\left[ \varOmega \hat{G} - D_R \left({\rm i}k +\frac{\partial }{\partial \phi}\right)\hat{G}\right]=0. \end{equation}
Inserting the expansion (2.11) into (2.13), at 
$O(1)$ we obtain
\begin{equation} \frac{\partial g}{\partial t} + \frac{\partial}{\partial \phi}\left( \varOmega g - D_R \frac{\partial g}{\partial \phi}\right)=0. \end{equation}
At 
$O(k)$, we have
\begin{equation} \frac{\partial b}{\partial t} + \frac{\partial}{\partial \phi} \left(\varOmega b - D_R \frac{\partial b}{\partial \phi}\right)=2 D_R \frac{\partial g}{\partial \phi} + \left( \varOmega^{eff} - \varOmega\right)g. \end{equation}
The average displacement field vanishes: 
$\langle b \rangle =0$. Equations (2.12), (2.14) and (2.15) are the main results of this paper.
It follows that for a particle undergoing a steady rotation (
$\varOmega =const.$), 
$g=1$ and 
$b=0$, which implies that 
$\varOmega ^{eff}=\varOmega$ and 
$D^{eff} = D_R$. As a result, a non-uniform or 
$\boldsymbol {q}$-dependent angular velocity is required to achieve a long-time dispersivity that is potentially different from the bare diffusivity 
$D_R$. To calculate the average drift, one needs to solve (2.14) and then take the average of 
$\varOmega g$. With the solution of (2.14), one can solve (2.15) and then use (2.12b) to calculate the dispersion coefficient.
3. Results
3.1. Simple shear
Consider the simple shear flow given by 
$\boldsymbol {u} = \dot {\gamma } y \boldsymbol {e}_x$, where 
$\boldsymbol {e}_x$ is the unit basis vector in the 
$x$ direction of the Cartesian coordinate system 
$(x,y,z)$, 
$\dot {\gamma }$ is the shear rate. The problem is non-dimensionalized using the time scale 
$\tau _R = 1/D_R$. Two dimensionless groups dictate the behaviour of the problem. The first is a Péclet number, 
$Pe = \dot {\gamma }\tau _R$, which compares the time scale of the flow to that of rotational diffusion. The second parameter is the Bretherton constant that characterizes the aspect ratio of the spheroid (Bretherton Reference Bretherton1962). The non-dimensional (scaled by 
$\tau _R$) angular velocity is
\begin{equation} \varOmega \tau_R =-\tfrac{1}{2}Pe [1 - B \cos(2\phi)]. \end{equation} For spherical particles, 
$B =0$, and the angular velocity is a constant, 
$\varOmega \tau _R = -Pe/2$. In this case, 
$g = 1$ and 
$\varOmega ^{eff} = \varOmega$. From (2.15), we readily obtain 
$b =0$ and 
$D^{eff} = D_R$. Because the angular velocity is a constant, the average drift is simply the angular velocity of the flow and the flow does not affect the dispersion coefficient. Similar to the classical Taylor dispersion in which a spatially non-uniform advection is present, an orientation-dependent angular velocity is required to have potentially flow-enhanced dispersion (Leahy et al. Reference Leahy, Cheng, Ong, Liddell-Watson and Cohen2013, Reference Leahy, Koch and Cohen2015).
To probe the effect of non-uniform angular velocity on the long-time drift and dispersion, we seek a regular series solution by writing
\begin{equation} g= \sum_{n=0}^\infty B^n g_n \quad \mathrm{and}\quad b = \sum_{n=0}^\infty B^n b_n. \end{equation}The resulting drift and dispersion coefficient are written as, respectively,
\begin{equation} \varOmega^{eff} \tau_R = \sum_{n=0}^\infty B^n \varOmega^{eff}_n \quad \mathrm{and}\quad \frac{D^{eff}}{D_R} = \sum_{n=0}^\infty B^n D^{eff}_n. \end{equation}
We calculate the series solution up to 
$O(B^6)$ in Appendix A. The drift terms are given by
$$\begin{gather} \varOmega^{eff}_0 =-\frac{1}{2}Pe, \quad \varOmega^{eff}_2 = \frac{Pe^3}{4 Pe^2+64},\quad \varOmega^{eff}_4= \frac{ Pe^5 (Pe^2-80)}{16 (Pe^2+16)^2 (Pe^2+64)}, \end{gather}$$
$$\begin{gather}\varOmega^{eff}_6 = \frac{ (Pe-4) Pe^7 (Pe+4) (Pe^2-368)}{32 (Pe^2+16)^3 (Pe^2+64) (Pe^2+144)}, \end{gather}$$where the odd terms are zero. For the dispersion coefficient, we obtain
$$\begin{gather} D^{eff}_0 =1,\quad D^{eff}_2 = \frac{Pe^2 (3 Pe^2-16)}{2 (Pe^2+16)^2}, \end{gather}$$
$$\begin{gather}D^{eff}_4 = \frac{Pe^4 (3 Pe^6-124 Pe^4-12992 Pe^2+20480)}{2 (Pe^2+16)^3(Pe^2+64)^2}, \end{gather}$$
$$\begin{gather}D^{eff}_6 = \frac{3C Pe^6}{2 (Pe^2+16)^4 (Pe^4+208 Pe^2+9216)^2}, \end{gather}$$
where 
$C = -36175872 + 51011584 Pe^2 - 919040 Pe^4 - 47328 Pe^6 - 170 Pe^8 + Pe^{10}$ and odd terms vanish.
In figure 1 we plot the average drift scaled by the drift of a sphere and the dispersion coefficient as a function of 
$Pe$ for several values of 
$B$. The scaled drift is shown in figure 1(a) and the non-dimensional dispersion coefficient is plotted in figure 1(b). The truncated series solution is shown by dashed lines. The circles in figure 1 are results obtained by solving (2.14) and (2.15) at steady state using a Fourier collocation method. The squares are from BD simulations of the orientational Langevin equation. In two dimensions, the Langevin equation (dimensional) is written as 
$\text {d} \varphi /{\rm d} t = \varOmega + \sqrt {2D_R}\xi$, where 
$\xi$ is a white-noise process satisfying 
$\overline {\xi (t)}=0$ and 
$\overline {\xi (t)\xi (t^\prime )} = \delta (t-t^\prime )$. Here, the overhead bar denotes an ensemble average and 
$\delta$ is the delta function. We remark that in the Langevin equation, the unbounded angular coordinate 
$\varphi$ is used in order to calculate the mean and mean-squared angular displacements, from which the drift and dispersion coefficient can be obtained. The full numerical solutions (circles) of (2.14) and (2.15) agree with the results from BD (squares), which validates our theory.
Figure 1. (a) Plots of the average angular drift velocity scaled by the sphere result as a function of 
$Pe$ for different values of 
$B$. (b) Plots of the non-dimensional effective long-time dispersion coefficient as a function of 
$Pe$ for different values of 
$B$. For spheres (
$B=0$), the dispersion coefficient is not affected by the flow. For non-spherical particles, shear-enhanced dispersion is observed.
For spheres, 
$B=0$, and the drift is equal to the constant angular velocity. As 
$B$ increases, the drift decreases compared with that of the sphere because the alignment term due to the rate of strain becomes more important. In the limit 
$Pe \to 0$, the drift of non-spherical particles approaches that of spheres, 
$\varOmega ^{eff}/\varOmega ^{eff}_{B=0} \to 1$. The reduction of the scaled drift occurs at non-zero 
$Pe$ and is most prominent for large 
$Pe$ where the scaled drift asymptotes to a plateau as 
$Pe \to \infty$. For non-spherical particles, we observe a shear-enhanced angular dispersion as shown in figure 1(b). Similar to the classical Taylor dispersion in linear position space, the enhanced angular dispersion is a result of the coupling between non-uniform advection and diffusion. For spheres, the angular velocity is constant and no shear-enhanced dispersion is observed. When 
$B >0$, the dispersion coefficient increases monotonically as a function of 
$Pe$ until it asymptotes to a plateau at large 
$Pe$. In dimensional terms, we have 
$D^{eff}/D_R = O(Pe^0)$ as 
$Pe \to \infty$, which is different from the classical Taylor dispersion of Brownian solutes in Poiseuille flow where 
$D^{eff}/D_T \sim Pe^2$ as 
$Pe \to \infty$. For dispersion in position space, 
$D_T$ is the translational diffusivity and 
$Pe = U L/D_T$ with 
$U$ the characteristic fluid velocity and 
$L$ the characteristic width of the channel. We further note that Brownian particles in unbounded shear flows exhibit anomalous diffusion in position space (San Miguel & Sancho Reference San Miguel and Sancho1979; Foister & Van De Ven Reference Foister and Van De Ven1980; Krishnan & Leighton Reference Krishnan and Leighton1992; Katayama & Terauti Reference Katayama and Terauti1996; Orihara & Takikawa Reference Orihara and Takikawa2011; Takikawa & Orihara Reference Takikawa and Orihara2012).
To understand the behaviour of the system at large 
$Pe$, we consider a perturbation expansion in powers of 
$1/Pe$,
$$\begin{gather} g =g_0 + \frac{1}{Pe}g_1 + \frac{1}{Pe^2}g_2+\cdots, \end{gather}$$
$$\begin{gather}b =b_0 + \frac{1}{Pe}b_1 + \frac{1}{Pe^2}b_2 +\cdots. \end{gather}$$Inserting the expansion (3.6a) into (2.14), one can solve the resulting equations order by order. We derive after some algebra that
$$\begin{gather} g_0 = \frac{\sqrt{1 - B^2}}{1 - B \cos(2\phi)}, \end{gather}$$
$$\begin{gather}g_1 = \frac{4B \sqrt{1-B^2} \sin(2 \phi)}{[1-B\cos(2\phi)]^3}. \end{gather}$$
We note that in obtaining the preceding solutions, the conservation conditions 
$\langle g_0 \rangle =1$ and 
$\langle g_1 \rangle =0$ are enforced. From (3.7), we obtain
\begin{equation} \frac{\varOmega^{eff}}{\varOmega^{eff}_{B=0}} = \sqrt{1-B^2} +O(1/Pe^2)\quad\mathrm{as}\ Pe \to \infty. \end{equation}
As 
$Pe \to \infty$, the scaled drift approaches a finite value that depends on 
$B$. Due to symmetry, the function 
$g_1$ does not contribute to the drift. The correction is 
$O(1/Pe^2)$, which comes from 
$g_2$ in the expansion. In figure 2(a) the leading-order asymptotic solution 
$g_0$ is plotted as a function of 
$\phi$ for 
$B=0.6$. The full numerical solution of (2.14) at 
$Pe=1000$ and 
$B=0.6$ is also shown in figure 2(a). In figure 4(a) we plot the leading-order drift given in (3.8) as a function of 
$B$ and the numerical solution of the macrotransport equations for 
$Pe=1000$. Good agreement between the asymptotic and numerical solutions is observed.
Figure 2. Plots of numerical solutions (dashed) at 
$Pe=1000$ and the leading-order asymptotic solutions (dotted) at large 
$Pe$ for (a) the average field and (b) the displacement field. In the plots, 
$B=0.6$. The numerical and asymptotic solutions have excellent agreement and cannot be distinguished visually.
Substituting (3.6b) and (3.7) into (2.15), we obtain at 
$O(1)$,
\begin{equation} \frac{\partial}{\partial \phi}\left[-\frac{1}{2}(1-B\cos(2\phi))b_0\right] = \left[-\frac{1}{2} \sqrt{1-B^2} +\frac{1}{2}(1-B\cos(2\phi))\right]g_0, \end{equation}where we have used (3.8). Defining
\begin{equation} \tilde{b}_0(\phi) = \frac{\sqrt{1-B^2}}{1-B\cos(2\phi)} \left[\arctan\left(\frac{(1+B) \tan \phi }{\sqrt{1-B^2}}\right) - \phi\right],\quad 0 \leq \phi \leq {\rm \pi}/2, \end{equation}
one can write the solution at 
$O(1)$ as
\begin{equation} b_0(\phi) = \begin{cases} \tilde{b}_0(\phi), & 0 \leq \phi \leq {\rm \pi}/2, \\ - \tilde{b}_0({\rm \pi}-\phi), & {\rm \pi}/2 \leq \phi \leq {\rm \pi}, \\ \tilde{b}_0(\phi - {\rm \pi}), & {\rm \pi}\leq \phi \leq 3{\rm \pi}/2,\\ - \tilde{b}_0(2{\rm \pi}-\phi), & 3{\rm \pi}/2 \leq \phi \leq 2{\rm \pi}. \end{cases} \end{equation}
In figure 2(b) the leading-order asymptotic solution 
$b_0$ is plotted as a function of 
$\phi$ for 
$B=0.6$. The full numerical solution of (2.15) at 
$Pe=1000$ and 
$B=0.6$ is also shown in figure 2(b). Good agreement between the asymptotic and numerical solutions is observed. Because of the symmetry of 
$\varOmega$ and 
$b_0$, the average 
$\langle \varOmega b_0 \rangle$ vanishes.
To obtain the first non-zero term of shear-induced dispersion, we therefore need to consider the 
$O(1/Pe)$ solution to 
$b$. At 
$O(1/Pe)$, the displacement field equation is given by
\begin{equation} \frac{\partial }{\partial \phi}\left[ -\frac{1}{2}[1-B \cos(2\phi)] b_1 - \frac{\partial b_0}{\partial \phi} \right] = 2 \frac{\partial g_0}{\partial \phi}. \end{equation}Similar to (3.10), we define
\begin{equation} \tilde{b}_1(\phi) = \frac{4 B \sqrt{1-B^2} \sin (2 \phi )}{[1-B \cos (2 \phi)]^3} \left[\arctan\left(\frac{(1+B) \tan \phi }{\sqrt{1-B^2}}\right) - \phi\right] + \frac{3 B^2-3}{[1-B \cos (2 \phi )]^3}, \end{equation}
which is valid for 
$\phi \in [0, {\rm \pi}/2]$. Using (3.13), a particular solution to 
$b_1$ is constructed as
\begin{equation} b_1^p(\phi) = \begin{cases} \tilde{b}_1(\phi), & 0 \leq \phi \leq {\rm \pi}/2, \\ \tilde{b}_1({\rm \pi}-\phi), & {\rm \pi}/2 \leq \phi \leq {\rm \pi}, \\ \tilde{b}_1(\phi - {\rm \pi}), & {\rm \pi}\leq \phi \leq 3{\rm \pi}/2,\\ \tilde{b}_1(2{\rm \pi}-\phi), & 3{\rm \pi}/2 \leq \phi \leq 2{\rm \pi}. \end{cases} \end{equation}The full solution can be written as
\begin{equation} b_1(\phi) = \frac{B^2+\sqrt{1-B^2}+2}{1-B^2} \frac{1}{1-B\cos(2\phi)} + b_1^p(\phi), \end{equation}
where the first term on the right-hand side is the homogeneous solution. To compare the 
$O(1)$ asymptotic solutions with numerical solutions, we first construct a hybrid approximation to 
$g_1$, given by 
$g_1^{num} =Pe( g^{num} - g_0)$, where the superscript ‘
$num$’ denotes the numerical solution. In figure 3(a), 
$g_1^{num}$ is compared with the asymptotic solution 
$g_1$, and good agreement is observed. Similarly, we define 
$b_1^{num} =Pe(b^{num} - b_0)$. The results for 
$b_1^{num}$ and 
$b_1$ are plotted in figure 3(b).
Figure 3. Plots of numerical solutions (dashed) at 
$Pe=1000$ and the 
$O(1/Pe)$ asymptotic solutions (dotted) at large 
$Pe$ for (a) the average field and (b) the displacement field. In the plots, 
$B=0.6$. The numerical approximation to 
$g_1$ is obtained via 
$g_1^{num} =Pe(g^{num} - g_0)$, where the superscript ‘
$num$’ denotes a numerical solution and 
$g_0$ is given in (3.7a). Similarly, 
$b_1^{num} =Pe( b^{num} - b_0)$. The numerical and asymptotic solutions have excellent agreement and cannot be distinguished visually.
Using (3.15) and the expansion
\begin{equation} D^{eff} = D_0^{eff} + \frac{1}{Pe}D_1^{eff}+\cdots, \end{equation}we obtain
\begin{equation} \frac{D_0^{eff}}{D_R} = 1 + \left\langle \frac{1}{2}[1-B\cos(2\phi)] b_1 \right\rangle = \frac{2+B^2}{2(1-B^2)}. \end{equation}
In figure 4(b) we compare the asymptotic solution given in (3.17) with the full numerical solution of the macrotransport equations for 
$Pe=1000$. We note that (3.17) agrees with the result obtained by Leahy et al. (Reference Leahy, Koch and Cohen2015) where a coordinate transformation is employed to map the orientational dynamics to a diffusion equation in the large-
$Pe$ limit. In contrast to their approach, the current theory follows closely the GTD framework and allows us to calculate both the average drift and effective dispersion for arbitrary Péclet numbers in generic linear flows.
Figure 4. Plots of numerical solutions (circles) at 
$Pe=1000$ and the leading-order asymptotic solutions (dashed) at large 
$Pe$ for (a) the scaled drift and (b) the dispersion coefficient.
3.2. Extensional flow
As another case study, we now consider an extensional flow where the angular velocity is 
$\boldsymbol {\varOmega } = Pe B \cos (2\phi ) \boldsymbol {e}_z$. The extensional flow tends to align the particle orientation with the extensional axis (see figure 5). The particle has a non-zero angular velocity when the orientation deviates from the extensional axis. As shown in figure 5(b), the direction of rotation depends on the particle orientation relative to the extensional axis.
Figure 5. (a) Schematic of the extensional flow. (b) Schematic of the direction of rotation due to the extensional flow.
The average field in the extensional flow is governed by
\begin{equation} \frac{\partial}{\partial \phi}\left(B \,Pe \cos(2\phi) g - \frac{\partial g}{\partial \phi}\right)=0. \end{equation}
Because the particle can align with both directions of the extensional axis (
$\phi ={\rm \pi} /4$ and 
$5{\rm \pi} /4$ in figure 5b), the average field has a periodicity of 
${\rm \pi}$, 
$g(\phi +{\rm \pi} ) = g(\phi )$. Defining 
$\phi ^\prime = \phi -{\rm \pi} /4$, we consider 
$g$ in the interval 
$\phi ^\prime \in [-{\rm \pi} /2, {\rm \pi}/2]$. In terms of the shifted variable 
$\phi ^\prime$, we have
\begin{equation} \frac{\partial}{\partial \phi^\prime}\left(B \,Pe \sin(2\phi^\prime) g - \frac{\partial g}{\partial \phi^\prime}\right)=0,\quad \phi^\prime \in \left[-\frac{\rm \pi}{2}, \frac{\rm \pi}{2}\right]. \end{equation}
From (3.19), we see that 
$g$ is an even function of 
$\phi ^\prime$; this can also be understood by examining figure 5(b). The average angular drift 
$\langle \varOmega g \rangle$ vanishes because 
$\varOmega g$ is an odd function of 
$\phi ^\prime$ in the interval 
$\phi ^\prime \in [-{\rm \pi} /2, {\rm \pi}/2]$.
Integrating (3.18) once, we obtain
\begin{equation} B \,Pe \cos(2\phi) g - \frac{\partial g}{\partial \phi}= C_1. \end{equation}
Averaging the above equation over one period, we have 
$\langle \varOmega g \rangle = C_1$. Because 
$\langle \varOmega g \rangle =0$ as determined from symmetry, we must have 
$C_1=0$. With this, we obtain
\begin{equation} g(\phi) = \frac{1}{A_1}\exp\left(\frac{1}{2}B \,Pe \sin(2\phi)\right),\quad \phi \in [0, 2{\rm \pi}], \end{equation}
where 
$A_1$ is determined from 
$\langle g \rangle =1$.
The displacement field in the extensional flow satisfies
\begin{equation} \frac{\partial}{\partial \phi}\left(B \,Pe \cos(2\phi) b - \frac{\partial b}{\partial \phi}\right)=2 \frac{\partial g}{\partial \phi} - B \,Pe \cos(2\phi) g = \frac{\partial g}{\partial \phi}. \end{equation}
Because (3.22) does not admit a simple closed-form analytic solution, we instead seek pertubative solutions in the following two limits: (1) 
$Pe\,B\ll 1$ and (2) 
$Pe\, B \gg 1$.
In the small 
$Pe\, B$ limit, we write
$$\begin{gather} g(\phi) = 1 + Pe\,B g_1(\phi) + \cdots, \end{gather}$$
$$\begin{gather}b(\phi) = 0 + Pe\,B b_1(\phi) + \cdots, \end{gather}$$where
\begin{equation} g_1= \tfrac{1}{2}\sin(2\phi)\quad \mathrm{and}\quad b_1 =\tfrac{1}{4} \cos (2 \phi). \end{equation}From this, we determine the dispersion coefficient as
\begin{equation} \frac{D^{eff}}{D_R} = 1 - (Pe\,B)^2\langle \cos(2\phi)b_1 \rangle+O(Pe^4B^4) = 1 - \frac{1}{8}Pe^2B^2+O(Pe^4B^4). \end{equation} In figure 6 we plot the effective dispersion coefficient as a function of 
$Pe\,B$ in the extensional flow. The leading-order asymptotic solution for small 
$Pe\,B$ in (3.25) is denoted by the dashed line. In the small 
$Pe\,B$ regime, the asymptotic solution agrees well with both the full numerical solution (circles) of the macrotransport equations and results obtained from BD simulations (squares). In the extensional flow, rotational dispersion is hindered and the dispersion coefficient vanishes in the large 
$Pe\,B$ limit. Because the extensional flow acts to align the particle orientation with the extensional axis, in the strong flow limit random reorientations are suppressed.
Figure 6. Plots of the effective dispersion coefficient as a function of 
$Pe\,B$ in the extensional flow. Circles are results obtained from the numerical solution of the Taylor dispersion theory, squares are results calculated from BD simulations and the dashed line denotes the leading-order asymptotic solution for small 
$Pe\,B$ given in (3.25).
We now consider the probability distributions in the strong flow limit characterized by 
$\epsilon = 1/(Pe\,B) \ll 1$. In the limit 
$Pe\,B\to \infty$ or 
$\epsilon \to 0$, the orientational distribution becomes a delta function localized at 
$\phi _0 = {\rm \pi}/4$ due to strong alignment. For strong flow, the particle orientation is closely aligned with the extensional axis, which implies the existence of a boundary layer near 
$\phi _0 = {\rm \pi}/4$. A dominant balance reveals that the boundary layer thickness is 
$O(\sqrt {\epsilon })$. Introducing the stretched coordinate 
$\xi = (\phi - \phi _0)/\sqrt {\epsilon }$, we rewrite (3.18) in the boundary layer as
\begin{equation} \frac{\partial}{\partial \xi}\left(\cos (2 \phi_0 + 2\epsilon^{1/2} \xi) \tilde{g} - \sqrt{\epsilon} \frac{\partial \tilde{g}}{\partial \xi}\right)=0, \end{equation}
where 
$\tilde {g}(\xi ) = g(\phi )$. Noting that 
$\cos (2\phi _0)=0$ and 
$\cos (2 \phi _0 + 2\epsilon ^{1/2} \xi )=-2 \xi \sqrt {\epsilon }+\frac {4}{3} \xi ^3 \epsilon ^{3/2}-\frac {4}{15} \xi ^5 \epsilon ^{5/2} +O(\epsilon ^{7/2})$, we expect an expansion of the form
\begin{equation} \tilde{g} = \frac{1}{\sqrt{\epsilon}}\tilde{g}_0 + \cdots, \end{equation}
where the leading-order term is 
$O(1/\sqrt {\epsilon })$ due to the conservation of probability. In the bulk (outside the boundary layer), the orientational distribution is zero to algebraic orders of 
$\sqrt {\epsilon }$.
The leading-order orientational distribution is governed by
\begin{equation} \frac{\partial}{\partial \xi}\left(2 \xi \tilde{g}_0 + \frac{\partial \tilde{g}_0}{\partial \xi}\right)=0, \end{equation}
which admits a solution of the form 
$\tilde {g}_0 = A_0 e^{-\xi ^2}$, where 
$A_0$ remains to be determined. In writing the solution, we have made use of the matching condition 
$\tilde {g}\to 0$ as 
$\xi \to \infty$. From the conservation, 
$\langle g \rangle =1$, we obtain 
$A_0 = \sqrt {{\rm \pi} }$. Due to symmetry, the leading-order solution is valid as 
$\phi$ approaches 
$\phi _0$ from either side. We then construct a leading-order composite solution as
\begin{equation} g = \sqrt{\frac{\rm \pi}{\epsilon}} \exp\left(- \frac{(\phi-\phi_0)^2}{\epsilon}\right) +\cdots, \end{equation}
where the dots denote higher-order corrections. In figure 7(a) we compare the leading-order solution in (3.29) with the numerical solution of the full equation for 
$Pe\,B=20$. We note that indeed the leading-order solution in (3.29) approaches a delta function upon appropriate scaling in the limit 
$\epsilon \to 0$.
We now proceed to analyse the displacement field, which in the boundary layer is expanded as
\begin{equation} \tilde{b} = \tilde{b}_0 +\cdots, \end{equation}
where we remark that 
$\tilde {b}$ is 
$O(1)$ at leading order. At leading order, we obtain
\begin{equation} {-}\frac{\partial}{\partial \xi}\left(2 \xi \tilde{b}_0 + \frac{\partial \tilde{b}_0}{\partial \xi} \right)= 2\xi \tilde{g}_0+ 2 \frac{\partial \tilde{g}_0}{\partial \xi}. \end{equation}The equation is solved by
\begin{equation} \tilde{b}_0=(C_0 - \sqrt{\rm \pi} \xi)\ {\rm e}^{-\xi^2}. \end{equation}
Using the normalization condition 
$\langle b \rangle =0$, we must have 
$\int _{-\infty }^\infty \tilde {b}_0\,\text {d}\xi =0$, which gives 
$C_0=0$. In terms of 
$\phi$, we may write
\begin{equation} b_0(\phi) =-\sqrt{\rm \pi} \frac{\phi- \phi_0}{\sqrt{\epsilon}}\exp \left(-\frac{(\phi-\phi_0)^2}{\epsilon}\right). \end{equation}
In figure 7(b) we compare the leading-order solution of 
$b$ in (3.33) with the numerical solution of the full equation for 
$Pe\,B=20$.
Taking the limit 
$\epsilon \to 0$, we obtain
\begin{equation} \frac{D^{eff}}{D_R} =1 - \langle \varOmega b \rangle \to 1- \frac{1}{\rm \pi}\int_{-\infty}^\infty (-2) \xi \tilde{b}_0(\xi) \,\text{d}\xi =0. \end{equation}As a result, we have shown that the dispersion coefficient vanishes in the strong flow limit. This asymptotic result agrees with the full numerical solution and the BD simulation results (see figure 6).
3.3. General linear flows
We now consider a general two-dimensional linear flow of the form 
$\boldsymbol {u} = \dot {\gamma }y \boldsymbol {e}_x + \alpha \dot {\gamma }x \boldsymbol {e}_y$, where 
$\alpha$ is a dimensionless parameter that controls the characteristics of the external flow. As a function of 
$\alpha$, the flow varies from purely rotational (
$\alpha =-1$) to simple shear (
$\alpha =0$) and extensional flow (
$\alpha =1$). From (2.1a,b), we obtain the angular velocity of the particle as
\begin{equation} \varOmega \tau_R = \tfrac{1}{2}Pe [-1+\alpha +B(1+\alpha)\cos(2\phi)]. \end{equation} To probe the effect of 
$\alpha$ on the long-time rotational dispersion coefficient, we solve the macrotransport equations (2.14) and (2.15) numerically using a Fourier collocation method for 
$Pe=100$. In figure 8 the numerical results are plotted as circles, and the results obtained from BD simulations are marked by squares. The dispersion coefficient exhibits a non-monotonic dependence on the flow parameter 
$\alpha$. The maximum of the dispersion coefficient is obtained when 
$\alpha = (1-B)/(1+B)$. For a general linear flow, the angular velocity is a combination of a constant part (purely rotational) and a varying part (extensional flow). For 
$\alpha =-1$, the flow is purely rotational, and we have 
$D^{eff}=D_R$. As 
$\alpha$ increases, the varying part of the angular velocity gives enhanced dispersion. In the absence of the constant part (
$\alpha =1$), particles becomes strongly aligned. As a result, the dispersion is suppressed (see § 3.2). To achieve maximum enhancement in dispersion, the two terms need to balance, which occurs when 
$1-\alpha = B(1+\alpha )$ or 
$\alpha = (1-B)/(1+B)$. In figure 8 this value of 
$\alpha$ is marked by a vertical line for 
$B=0.6$. On the right side of the vertical line, the extensional flow becomes dominant, and the dispersion coefficient decays to zero as 
$\alpha$ increases and approaches that of the extensional flow (
$\alpha =1$).
Figure 8. Plots of the dispersion coefficient as a function of 
$\alpha$ for 
$Pe=100$ and 
$B=0.6$. Circles are results from numerical solutions of the macrotransport equations. Squares are results from BD simulations. The vertical line is given by the equation 
$\alpha =(1-B)/(1+B)$. The solid line is the asymptotic solution in (3.41).
In the large-
$Pe$ limit, we consider a perturbation expansion in powers of 
$1/Pe$,
$$\begin{gather} g = g_0 + \frac{1}{Pe}g_1+\cdots, \end{gather}$$
$$\begin{gather}b = b_0 + \frac{1}{Pe}b_1+\cdots. \end{gather}$$Following the notations and calculations in § 3.1, we obtain
$$\begin{gather} g_0= \frac{A_1}{1-\alpha -B(1+\alpha)\cos(2\phi)}, \end{gather}$$
$$\begin{gather}g_1 = \frac{4B(1+\alpha)A_1\sin(2\phi)}{[1-\alpha -B(1+\alpha)\cos(2\phi)]^3}, \end{gather}$$
where 
$A_1=\sqrt {(1-\alpha )^2-B^2(1+\alpha )^2}$. We note that (3.37a) and (3.37b) are only valid in the region 
$1-\alpha > B(1+\alpha )$, which is shown as the shaded region in figure 9. At the boundary, 
$1-\alpha = B(1+\alpha )$, the solutions of (3.37a) and (3.37b) become singular at 
$\phi _0$, where 
$\cos (2\phi _0)=1$ or 
$\phi _0\in \{0,{\rm \pi},2{\rm \pi} \}$. One may similarly obtain the displacement field as
$$\begin{gather} \tilde{b}_0=A_1 \frac{\arctan\left(\dfrac{1-\alpha +B(1+\alpha)}{A_1}\tan\phi\right) -\phi}{1- \alpha - B(1+\alpha)\cos(2\phi)}, \end{gather}$$
$$\begin{gather}\tilde{b}_1 = \frac{C_1\left[\arctan\left(\dfrac{(1-\alpha +B(1+\alpha))}{A_1}\tan\phi\right) -\phi\right]+C_2}{[1-\alpha - B(1+\alpha)\cos(2\phi)]^3}+ \frac{C_3}{1-\alpha - B(1+\alpha)\cos(2\phi)}, \end{gather}$$where
$$\begin{gather} C_1=4B(1+\alpha)A_1,\quad C_2=-3 A_1^2, \end{gather}$$
$$\begin{gather}C_3=\frac{B^2(1+\alpha)^2+(1-\alpha)(2(1-\alpha)+A_1)}{A_1^2}. \end{gather}$$
Figure 9. Plot of the parameter space 
$(\alpha, B)$. The solid line is given by 
$1-\alpha = B(1+\alpha )$. In the shaded region, we have 
$1-\alpha > B(1+\alpha )$.
The leading-order dispersion coefficient can then be shown to be
\begin{equation} \frac{D_0^{eff}}{D_R} = \frac{2(1-\alpha)^2+B^2(1+\alpha)^2}{2(1-\alpha)^2-2B^2(1+\alpha)^2}. \end{equation}
For simple shear, 
$\alpha =0$, we recover (3.17). In figure 8 the solution of (3.41) for 
$B=0.6$ is plotted as a solid line. The asymptotic solution agrees well with both the numerical and BD results. As expected, the asymptotic solution deviates from the numerical or BD results as 
$\alpha$ approaches the singular boundary where 
$1-\alpha = B(1+\alpha )$. Furthermore, the solution of (3.41) becomes singular when 
$1-\alpha = B(1+\alpha )$.
In Appendix B, using asymptotic theory, we analyse the macrotransport equations and characterize the long-time dispersion behaviour for the case of 
$1-\alpha = B(1+\alpha )$ in the strong flow limit. When 
$1-\alpha = B(1+\alpha )$, our analysis reveals that there exist boundary layers of thickness 
$O(\epsilon ^{1/3})$ at 
$\phi = n {\rm \pi}$, 
$n \in \mathbb {Z}$, where 
$\epsilon = 1/(Pe(1-\alpha ))$. We show that the dispersion coefficient is 
$O(\epsilon ^{-2/3})$ as 
$\epsilon \to 0$. For a fixed 
$\alpha$, this means that 
$D^{eff}/D_R = O(Pe^{2/3})$.
Our results show that depending on the parameters 
$\alpha$ and 
$B$, the rotational dispersion coefficient exhibits distinctly different scaling in the large-
$Pe$ limit. For 
$1-\alpha > B(1+\alpha )$, we have 
$D^{eff}/D_R = O(1)$ as 
$Pe \to \infty$. On the curve where 
$1-\alpha = B(1+\alpha )$, 
$\alpha \neq 1$, 
$D^{eff}/D_R =O(Pe^{2/3})$. Finally, for extensional flow, 
$D^{eff}/D_R \to 0$ as 
$Pe \to \infty$. We note that as a function of 
$\alpha$, the transition from the case of 
$1-\alpha =B(1+\alpha )$ to the vanishing dispersion at 
$\alpha =1$ occurs in a small region, which is shown by the rapid decay of 
$D^{eff}$ in figure 8. Within this small region, a perturbation analysis for a small deviation from the curve 
$1-\alpha =B(1+\alpha )$ may be useful. Such an analysis is left for future considerations.
4. Concluding remarks
In this paper we have developed a GTD theory that describes the long-time rotational dispersion of a spheroidal particle in linear flows that is constrained to rotate in the velocity-gradient plane. As is standard for Taylor dispersion, the average drift and the effective dispersion are treated in a unified framework by leveraging a flux-averaging method in Fourier space. The results obtained from the continuum theory are corroborated by BD simulations of the orientational equation of motion. Using asymptotic analysis in the strong flow limit, we have shown that a simple shear enhances while an extensional flow hinders the long-time rotational dispersion. More specifically, we showed that 
$D^{eff}/D_R=O(1)$ in simple shear and 
$D^{eff}/D_R \to 0$ in extensional flow as 
$Pe \to \infty$ for a given non-zero 
$B$. These results reveal that the long-time rotational dispersion depends qualitatively on the characteristics of the background flow.
While we focused on the long-time dispersion in steady flows, our GTD theory applies equally well to time-periodic flows such as oscillatory shear. In oscillatory shear one needs to first obtain the long-time (time-dependent) solutions to the average and displacement fields. In addition to the cell average employed for steady linear flows, a time average over the oscillation period is needed to obtain the long-time transport coefficients.
In linear flows as considered in this paper, the angular velocity of the particle due to the flow does not depend on the spatial position. This spatial uniformity allows us to consider the conservation equation in orientation space only. A difficulty would arise if one wishes to consider the rotational dispersion of a particle in flows in the presence of no-slip boundaries such as pressure-driven channel flows. In such a quadratic flow field, the angular velocity depends linearly on the transverse coordinate. Mathematically, the marginalization of 
$P(\boldsymbol {x}, \boldsymbol {q}, t)$ does not lead to a closed equation for 
$\varPsi (\boldsymbol {q}, t)$. In this case, one often needs to solve the full probability distribution 
$P(\boldsymbol {x}, \boldsymbol {q}, t)$ in order to calculate the long-time translational or rotational transport properties. As a concrete example, consider the rotational dynamics of a ‘point’ Brownian particle in a planar Poiseuille flow. The Smoluchowski equation can be written as
\begin{equation} \frac{\partial P}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{u} P - D_T \boldsymbol{\nabla} P) + \boldsymbol{\nabla}_R \boldsymbol{\cdot}(\boldsymbol{\varOmega} P - D_R \boldsymbol{\nabla} P)=0, \end{equation}
where 
$\boldsymbol {\varOmega }$ depends on 
$y$, which is the transverse coordinate. Following the consideration given in § 2, one can show that the average field is governed by
\begin{equation} \frac{\partial }{\partial y}\left(-D_T \frac{\partial g}{\partial y} \right) + \frac{\partial}{\partial \phi}\left( \varOmega(y) g -D_R \frac{\partial g}{\partial \phi}\right)=0. \end{equation}
At the channel walls (
$\kern0.09em y=\pm H$), the no-flux condition, 
$-\partial g/(\partial y)=0$, is imposed. The normalization involves both spatial and orientational integrals, which is given by 
$\langle \int _{-H}^H g(y, \phi ) \,\text {d} y \rangle =1$. Similarly, one may write down the governing equation for the displacement field.
The rotational Taylor dispersion theory can accommodate particles of different shapes such as chiral or ring-shaped particles (Singh, Koch & Stroock Reference Singh, Koch and Stroock2013). In the presence of hydrodynamic translation–rotation coupling, one again needs to consider the dynamics in both position and orientation space. In the case of general particle shapes, much like the case of a Brownian particle in Poiseuille flow, one performs the marginalization to obtain net orientational distributions after solving the average and displacement fields as functions of both 
$\boldsymbol {x}$ and 
$\boldsymbol {q}$.
The rotational Taylor dispersion theory is developed in two dimensions where the particle rotates about a single axis. To accommodate particles that rotate in three-dimensional space, one may consider the rotational dynamics about perpendicular axes. For a particular chosen axis of rotation, one can then utilize the current theory to characterize the long-time dynamics. We note that measuring rotational diffusion about perpendicular axes in three dimensions has been used in experiments (Rogers et al. Reference Rogers, Lisicki, Cichocki, Dhont and Lang2012; Carrasco-Fadanelli et al. Reference Carrasco-Fadanelli, Mao, Nakakomi, Xu, Yamamoto, Yanagishima and Buttinoni2024). To conclude, we also note that similar to classical Taylor dispersion, our theory is only valid at long times.
Funding
This work is supported by the Faculty of Engineering at the University of Alberta.
Declaration of interests
The author report no conflict of interest.
Appendix A. Asymptotic solution in simple shear
In this appendix we discuss the asymptotic solution for nearly spherical particles in simple shear in two dimensions. Inserting the expansion (3.2a,b) into (2.14), we obtain at steady state
\begin{equation} \frac{1}{2}Pe \frac{\partial g_n}{\partial \phi} + \frac{\partial^2 g_n}{\partial \phi^2} = \frac{1}{2}Pe \frac{\partial }{\partial \phi}[\cos(2\phi) g_{k-1}], \end{equation}
where 
$k=1,2,\ldots\,$. The normalization becomes 
$\langle g_0 \rangle =1$ and 
$\langle g_k \rangle =0$ for 
$k=1,2,\ldots\,$.
The first four terms in the expansion are written as
$$\begin{gather} g_0(\phi) =1,\quad g_1(\phi)=a_{11}\cos(2\phi) + a_{12}\sin(2\phi), \end{gather}$$
$$\begin{gather}g_2(\phi)=a_{21}\cos (4 \phi )+a_{22}\sin (4 \phi ), \end{gather}$$
$$\begin{gather}g_3(\phi) = a_{31} \cos(2\phi) + a_{32}\sin(2\phi) + a_{33}\cos(6\phi) + a_{34}\sin(6\phi), \end{gather}$$where
$$\begin{gather} a_{11}= \frac{Pe^2}{Pe^2+16}, \quad a_{12}= \frac{4 Pe }{Pe^2+16}, \end{gather}$$
$$\begin{gather}a_{21}=\frac{Pe^2 (Pe^2-32)}{2 (Pe^2+16) (Pe^2+64)},\quad a_{22}=\frac{6 Pe^3}{(Pe^2+16) (Pe^2+64)}, \end{gather}$$
$$\begin{gather}a_{31}= \frac{Pe^4(Pe^2-80)}{4 (Pe^2+16)^2 (Pe^2+64)}, \end{gather}$$
$$\begin{gather}a_{32}=-\frac{Pe^3 (-16 Pe^4-2176Pe^2+18432)}{4 (Pe^2+16)^2 (Pe^2+64) (Pe^2+144)}, \end{gather}$$
$$\begin{gather}a_{33}=-\frac{Pe^3 (2816 Pe-Pe^3 (Pe^2-160))}{4 (Pe^2+16)^2(Pe^2+64) (Pe^2+144)}, \end{gather}$$
$$\begin{gather}a_{34}=-\frac{Pe^3 (6144-24 Pe^4)}{4 (Pe^2+16)^2 (Pe^2+64) (Pe^2+144)}. \end{gather}$$ The solution at 
$O(B^4)$ is written as
\begin{equation} g_4(\phi) = a_{41} \cos(4\phi) + a_{42}\sin(4\phi) + a_{43}\cos(8\phi) + a_{44}\sin(8\phi), \end{equation}where
$$\begin{gather} a_{41} = \frac{Pe^2 (a_{31}+a_{33})-8Pe (a_{32}+a_{34})}{2 (Pe^2+64)}, \end{gather}$$
$$\begin{gather}a_{42} = \frac{8Pe (a_{31}+a_{33})+Pe^2 (a_{32}+a_{34})}{2 (Pe^2+64)}, \end{gather}$$
$$\begin{gather}a_{43} = \frac{Pe (Pe\,a_{33} -16 a_{34})}{2 (Pe^2+256)},\quad a_{44}= \frac{Pe (16 a_{33}+Pe\, a_{34})}{2 (Pe^2+256)}. \end{gather}$$
The solution at 
$O(B^5)$ is given by
\begin{align} g_5(\phi) &= a_{51} \cos (2 \phi )+a_{52} \sin (2 \phi )+a_{53} \cos (6 \phi )\nonumber\\ &\quad +a_{54} \sin (6 \phi )+a_{55} \cos (10 \phi )+a_{56} \sin (10 \phi ), \end{align}where
$$\begin{gather} a_{51} = \frac{Pe (a_{41} Pe-4 a_{42})}{2 (Pe^2+16)},\quad a_{52}=\frac{Pe (4 a_{41}+a_{42} Pe)}{2 (Pe^2+16)}, \end{gather}$$
$$\begin{gather}a_{53}= \frac{Pe^2 (a_{41}+a_{43})-12Pe (a_{42}+a_{44})}{2 (Pe^2+144)},\quad a_{54}=\frac{12Pe (a_{41}+a_{43})+Pe^2 (a_{42}+a_{44})}{2 (Pe^2+144)}, \end{gather}$$
$$\begin{gather}a_{55}=\frac{Pe (a_{43} Pe-20 a_{44})}{2 (Pe^2+400)},\quad a_{56}=\frac{Pe (20 a_{43}+a_{44} Pe)}{2 (Pe^2+400)}. \end{gather}$$
The solution at 
$O(B^6)$ is given by
\begin{align} g_6(\phi) &= a_{61} \cos (4 \phi )+a_{62} \sin (4 \phi )+a_{63} \cos (8 \phi )\nonumber\\ &\quad +a_{64} \sin (8 \phi )+a_{65} \cos (12 \phi )+a_{66} \sin (12 \phi ), \end{align}where
$$\begin{gather} a_{61}=\frac{Pe^2 (a_{51}+a_{53})-8Pe (a_{52}+a_{54})}{2 (Pe^2+64)},\quad a_{62}=\frac{8Pe (a_{51}+a_{53})+Pe^2 (a_{52}+a_{54})}{2 (Pe^2+64)}, \end{gather}$$
$$\begin{gather}a_{63}=\frac{Pe^2 (a_{53}+a_{55})-16Pe (a_{54}+a_{56})}{2 (Pe^2+256)},\quad a_{64}=\frac{16Pe (a_{53}+a_{55})+Pe^2 (a_{54}+a_{56})}{2 (Pe^2+256)}, \end{gather}$$
$$\begin{gather}a_{65}=\frac{Pe (a_{55} Pe-24 a_{56})}{2 (Pe^2+576)},\quad a_{66}=\frac{Pe (24 a_{55}+a_{56} Pe)}{2 (Pe^2+576)}. \end{gather}$$ Similarly, the displacement field can be solved order by order. At 
$O(B^n)$, (2.15) is given by
\begin{equation} \frac{\partial}{\partial \phi}\left(\varOmega_0 b_n + \varOmega_1 b_{n-1} - \frac{\partial b_n}{\partial \phi}\right) = 2 \frac{\partial g_n}{\partial \phi} + \sum_{i= n-j}\sum_{j} \varOmega_i^{eff} g_j - (\varOmega_0 g_n + \varOmega_1 g_{n-1}), \end{equation}
where 
$\varOmega _0 = -\frac {1}{2}Pe$, 
$\varOmega _1 = \frac {1}{2}Pe\cos (2\phi )$ and 
$n=0,1,\ldots\,$. In (A10), 
$b_k=0$ for 
$k<0$. At 
$O(1)$, we have
\begin{equation} \frac{\partial}{\partial \phi}\left(\varOmega_0 b_0 - \frac{\partial b_0}{\partial \phi}\right) = 2 \frac{\partial g_0}{\partial \phi} + (\varOmega_0^{eff}-\varOmega_0)g_0. \end{equation}
Since 
$g_0 =1$, we have 
$b_0=0$. At 
$O(B)$, we have
\begin{equation} \frac{\partial}{\partial \phi}\left(\varOmega_0 b_1 +\varOmega_1 b_0 - \frac{\partial b_1}{\partial \phi}\right) = 2 \frac{\partial g_1}{\partial \phi} + (\varOmega_0^{eff}-\varOmega_0)g_1+ (\varOmega_1^{eff}-\varOmega_1)g_0. \end{equation}The solution is given by
\begin{equation} b_1(\phi)=\frac{Pe}{2 (Pe^2+16)^2}\left[Pe (Pe^2-48) \sin (2 \phi )+(64-12 Pe^2) \cos (2 \phi)\right]. \end{equation}
At 
$O(B^2)$, we have
\begin{align} b_2(\phi) &= \frac{Pe^2}{8 (Pe^2+16)^2 (Pe^2+64)^2} [-4 Pe (17 Pe^4+208 Pe^2-19456) \cos (4 \phi ) \nonumber\\ &\quad + 3 (Pe^6-144 Pe^4-6656 Pe^2+32768) \sin (4 \phi)]. \end{align}Following this procedure, we have
$$\begin{gather} b_3(\phi) = c_{31} \cos(2\phi) + c_{32}\sin(2\phi) + c_{33}\cos(6\phi) + c_{34}\sin(6\phi), \end{gather}$$
$$\begin{gather}b_4(\phi) = c_{41} \cos(4\phi) + c_{42}\sin(4\phi) + c_{43}\cos(8\phi) + c_{44}\sin(8\phi), \end{gather}$$
\begin{gather} b_5(\phi) = c_{51} \cos (2 \phi )+c_{52} \sin (2 \phi )+c_{53} \cos (6 \phi) \nonumber\\ \,\qquad\qquad\qquad\qquad +\,c_{54} \sin (6 \phi )+c_{55} \cos (10 \phi )+c_{56} \sin (10 \phi), \end{gather}
\begin{gather} b_6(\phi) = c_{61} \cos (4 \phi )+c_{62} \sin (4 \phi )+c_{63} \cos (8 \phi) \nonumber\\ \,\qquad\qquad\qquad\qquad +\,c_{64} \sin (8 \phi )+c_{65} \cos (12 \phi )+c_{66} \sin (12 \phi), \end{gather}
where 
$c_{ij}$ can be readily obtained by inserting the expressions into (A10).
Appendix B. Asymptotic analysis for 
$1-\alpha = B(1+\alpha )$ in the strong flow limit
We develop an asymptotic solution in the large-
$Pe$ limit for a general linear flow under the constraint that 
$1-\alpha = B(1+\alpha )$. Under this condition, we have 
$\varOmega \tau _R = -Pe(1-\alpha ) \sin ^2\phi$. For convenience, we define 
$\epsilon = 1/(Pe(1-\alpha ))$ and consider the behaviour of the macrotransport equations when 
$\epsilon \ll 1$. We note that 
$\alpha \neq 1$ is assumed.
The dimensionless average field at steady state, 
$g=g(\phi )$, is governed by
\begin{equation} \frac{\text{d}}{\text{d}\phi}(\sin^2\phi g) +\epsilon \frac{\text{d}^2 g}{\text{d}\phi^2}=0. \end{equation}
Noting that 
$g(\phi +{\rm \pi} ) = g(\phi )$, we consider 
$\phi \in [-{\rm \pi} /2, {\rm \pi}/2]$. The normalization is given by
\begin{equation} \frac{1}{\rm \pi}\int_{-{\rm \pi}/2}^{{\rm \pi}/2} g(\phi)\,\text{d}\phi =1. \end{equation}
An asymptotic solution to (B1) subject to the normalization condition (B2) has been developed by Brenner (Reference Brenner1970). Since the solution to 
$g$ is required for the consideration of the displacement field 
$b$, we first repeat the calculations of Brenner (Reference Brenner1970) for 
$g$ and then proceed to consider 
$b$.
Integrating (B1) once, we obtain
\begin{equation} \sin^2\phi g +\epsilon \frac{\text{d} g}{\text{d}\phi}=C^\prime(\epsilon), \end{equation}
where we remark that 
$C^\prime$ only depends on 
$\epsilon$. In the bulk, we have 
$\sin ^2\phi g \sim C^\prime$ or 
$g\sim C^\prime /\sin ^2\phi$. This bulk result breaks down when 
$\sin ^2\phi g \sim \epsilon \,\text {d} g/(\text {d}\phi )$ or 
$\phi =O(\epsilon ^{1/3})$, which suggests the existence of a boundary layer with thickness 
$O(\epsilon ^{1/3})$ at 
$\phi =0$. In the boundary layer, from (B3), we obtain 
$g\sim C^\prime \epsilon ^{-2/3}$. Using the normalization condition (B2), one can show that 
$C^\prime = O(\epsilon ^{1/3})$. Following Brenner (Reference Brenner1970), we define 
$C^\prime = C \epsilon ^{1/3}$, where 
$C=O(1)$. In the bulk, we define 
$g(\phi ) = C\epsilon ^{1/3} H(\phi ; \epsilon )$. Introducing the stretched coordinate 
$s = \phi /\epsilon ^{1/3}$, we define 
$g(s) = C\epsilon ^{-1/3} Q(s;\epsilon )$ in the boundary layer. Expanding 
$H$ and 
$Q$, we write
\begin{equation} H = H_0 +\epsilon H_1+\cdots,\quad Q = Q_0 + \epsilon^{2/3}Q_1+\cdots, \end{equation}
where 
$H_0 = 1/\sin ^2\phi$ and 
$Q_0$ is governed by
\begin{equation} s^2 Q_0(s) + \frac{\text{d} Q_0}{\text{d} s} =1. \end{equation}
The solution to 
$Q_0$ is given by
\begin{equation} Q_0(s) = \exp({-s^3/3})\int_{-\infty}^s \exp({r^3/3})\,\text{d} r. \end{equation}To leading order, the normalization condition reduces to
\begin{equation} C \int_{-\infty}^{+\infty} Q_0(s) \,\text{d} s ={\rm \pi}, \end{equation}which allows us to obtain
\begin{equation} C = \frac{3^{5/6} {\rm \pi}}{6 \sqrt{3} (\varGamma (4/3))^2+2^{4/3} \sqrt{\rm \pi} \varGamma (7/6)} = \frac{3^{1/3}{\rm \pi}}{[\varGamma(1/3)]^2}, \end{equation}
where 
$\varGamma (z) = \int _0^\infty t^{z-1}\ {\rm e}^{-t}\,\text {d} t$ is the gamma function and 
$C\approx 0.6313$. We can construct a leading-order composite solution to 
$g$ as 
$C \epsilon ^{-1/3} Q_0(\phi /\epsilon ^{1/3})$. From this, we obtain the leading-order drift as 
$\varOmega ^{eff}\tau _R= -C \epsilon ^{-2/3} +o(\epsilon ^{-2/3})$. In figure 10(a) we compare the asymptotic result with the results from numerical solutions of the macrotransport equations for 
$B=1$ and 
$\alpha =0$ in the large-
$Pe$ regime.
Figure 10. (a) Plots of the numerical and asymptotic solutions to the average angular drift velocity as a function of 
$1/\epsilon$ for 
$B=1,\alpha =0$. (b) Plots of the numerical solutions to the non-dimensional effective long-time dispersion coefficient as a function of 
$1/\epsilon$ for 
$B=1,\alpha =0$. The solid line is given by 
$D^{eff}/D_R = \tilde {C} \epsilon ^{-2/3}$, where 
$\tilde {C}$ is obtained from numerical data.
The displacement field is governed by
\begin{equation} - \frac{\text{d}}{\text{d}\phi}(\sin^2\phi b) - \epsilon \frac{\text{d}^2 b}{\text{d}\phi^2} = 2 \epsilon \frac{\text{d} g }{\text{d}\phi} + \epsilon (\varOmega^{*{eff}} - \varOmega^*)g, \end{equation}
where 
$\varOmega ^{*{eff}} = \varOmega ^{eff}\tau _R$ and 
$\varOmega ^* = \varOmega \tau _R$. Noting that the dominant term on the right-hand side outside the boundary layer is 
$\epsilon \varOmega ^*g=O(\epsilon ^{1/3})$, we have 
$b = O(\epsilon ^{1/3})$ in the bulk. In the boundary layer, a dominant balance reveals that 
$b=O(\epsilon ^{-1/3})$. While the resulting equations are not analytically tractable, the above analysis allows us to show that 
$D^{eff}/D_R = O(\epsilon ^{-2/3})$ as 
$\epsilon \to 0$. In figure 10 we present the dispersion coefficient obtained from numerical solutions of the macrotransport equations for 
$B=1$ and 
$\alpha =0$ in the large-
$Pe$ regime. The solid line has the form 
$D^{eff}/D_R = \tilde {C} \epsilon ^{-2/3}$, where 
$\tilde {C}$ is a numerical constant that fits the data closely.
