Active microrheology of dilute colloidal dispersions at large Péclet numbers

Abstract

We revisit the model problem of Squires & Brady (Phys. Fluids, vol. 17, 2005, 073101), where a Brownian probe is dragged through a dilute dispersion of Brownian bath particles. In this problem, the microrheology due to excluded-volume interactions is represented by an effective viscosity, with the nonlinearity in the driving force entering via the dependence of the viscosity increment (relative to the viscosity of a pure solvent) upon the deformation of the dispersion microstructure. Our interest is in the limit of large Péclet numbers, $ P{\kern-1pt}e\gg 1$, where the microstructural deformation adopts the form of a boundary layer about the upstream hemisphere of the probe. We show that the boundary-layer solution breaks down at the equator of the probe and identify a transition region about the equator, connecting the layer to a downstream wake. The microstructural deformation in this region is governed by a universal boundary-value problem in a semi-bounded two-dimensional domain. The equatorial region continues downstream as a transition layer, which separates the wake of the probe from the undisturbed ambient; in that layer, the microstructure is governed by a one-dimensional heat-like equation. Accounting for the combined contributions from the respective asymptotic provinces we find the approximation $ ({1}/{2})[1+ (\ln P{\kern-1pt}e + 1.046)/ P{\kern-1pt}e]$ for the ratio of the large-$ P{\kern-1pt}e$ viscosity increment to the corresponding linear-response increment. Our asymptotic approximation is in excellent agreement with the increment predicted by a finite-difference numerical calculation of the microstructure deformation, tailored to the large-$ P{\kern-1pt}e$ topology.

1. Introduction

Active microrheology, describing the non-equilibrium mechanics of a probe motion in a colloidal dispersion, is a vibrant research field (Squires & Mason 2010; Puertas & Voigtmann 2014; Zia 2017). The pioneering theoretical paradigm for arbitrary magnitude of the external force, motivated by early experiments (Furst 2003; Habdas et al. 2004), was introduced by Squires & Brady (2005) (denoted hereafter by SB05). The model problem considered by SB05 can be traced back to theoretical investigations of non-Brownian dispersions (Davis & Hill 1992; Almog & Brenner 1997). It consists of Brownian bath particles dispersed in a viscous liquid (the ‘solvent’), with a Brownian probe being dragged through the dispersion. The probe experiences forces due to both the liquid and the bath particles. Assuming hydrodynamic radii small compared with those representing hard-sphere collisions, SB05 neglected hydrodynamic interactions. The particles therefore interact only through excluded volumes.

Assuming a statistically homogeneous dispersion, SB05 employed an averaging technique, considering both the case where the probe is dragged by a constant force and that where it is dragged at a constant velocity. In the fixed-force problem, they obtained the mean velocity of the probe, while in the fixed-velocity problem they calculated the mean force on the probe. Both averages are functionals of the pair-distribution function $g$ – the conditional probability for finding a bath particle in the vicinity of the probe. By interpreting these averaged relations in the dilute limit via a Stokes law template, SB05 defined an effective viscosity. Its dependence upon the deformed bath-particle microstructure, as expressed via $g$ , introduces an implicit nonlinear dependence upon the external forcing. In the dilute limit, the microstructure is governed by an advection–diffusion problem, where the forcing is manifested in a uniform ‘probability advection’, and the excluded-volume interactions are represented by a no-flux condition on an ‘exclusion sphere’, whose radius is the sum of the probe radius and the radius of a bath particle.

In a dimensionless formulation, the microstructure problem is governed by a single parameter, the Péclet number $ P{\kern-1pt}e$ , expressing the magnitude of the external forcing. SB05 solved the leading-order problem for both small and large $ P{\kern-1pt}e$ . For small $ P{\kern-1pt}e$ , they perturbed about the uniform equilibrium distribution. For large $ P{\kern-1pt}e$ , they identified a boundary layer about the upstream portion of the exclusion sphere, where $g=O( P{\kern-1pt}e)$ ; outside that layer $g\approx 1$ , except in a downstream wake where $g\approx 0$ . SB05 found that the viscosity increment, relative to that of the solvent, approaches finite values in both limits, with the large- $ P{\kern-1pt}e$ limit smaller by a factor of two. SB05 verified their ‘force-thinning’ prediction using a semi-numerical solution of the advection–diffusion problem, constructed using eigenfunction expansions.

The force thinning predicted by SB05 is in agreement with the Brownian simulations of Carpen & Brady (2005) and is qualitatively supported by experiments (Meyer et al. 2006; Wilson et al. 2009; Sriram, Meyer & Furst 2010) carried out under conditions where the hard-sphere model is applicable. The paradigm of SB05 was later generalised to incorporate hydrodynamic interactions for both the fixed-force (Khair & Brady 2006) and fixed-velocity (Swan & Zia 2013) models. The related calculation of probe diffusion under an external force was carried out by Zia & Brady (2010).

Our interest is in illuminating the approach to the Newtonian plateau at large $ P{\kern-1pt}e$ . To that end, we revisit the problem of SB05, focusing upon the asymptotic limit $ P{\kern-1pt}e\to \infty$ . Going beyond the leading-order calculation of SB05, we calculate the $\textrm {ord}(1)$ correction to the $\textrm {ord}( P{\kern-1pt}e)$ approximation for $g$ in the boundary layer, observing a divergence at the equator. This breakdown indicates the formation of an equatorial transition region. Continuing downstream, this region connects to a transition layer surrounding the wake. Following asymptotic analysis in the separate asymptotic regions, we obtain a refined large- $ P{\kern-1pt}e$ approximation to the effective viscosity increment. It is corroborated by a finite-difference numerical solution of the microstructure problem.

The paper is arranged as follows. In the next section we review the physical problem addressed by SB05, their key approximations and the averaged expressions they derived for the effective viscosity increment. The advective–diffusive problem governing the deformed microstructure is discussed in § 3. In these two sections we generally follow the notation of SB05, although certain deviations were found beneficial. Section 4 describes the force thinning discovered by SB05 and outlines the objective of the present contribution in that context. The boundary-layer problem is analysed in § 5. The equatorial-region scalings are identified in § 6, where an appropriate universal problem is formulated and solved. The effective viscosity increment is calculated in § 7. The transition layer surrounding the wake is discussed in § 8. In § 9 we describe a finite-difference scheme for the solution of the microstructure problem, tailored to the boundary-layer structure at large $ P{\kern-1pt}e$ . We conclude in § 10, proposing possible future extensions.

2. Background: effective viscosity

The system considered by SB05 consists of a colloidal ‘probe’ particle and colloidal ‘bath’ particles (average number density $n$ ), dispersed in a liquid bath of viscosity $\eta$ . The self-diffusivities of the probe and bath particles are denoted by $D_a$ and $D_b$ , respectively. The interaction (e.g. electrostatic or steric) between the particles is modelled via excluded volumes, with the hard-sphere radii of the probe and bath particles denoted by $a$ and $b$ , respectively. These are assumed sufficiently larger than the hydrodynamic radii, so the hydrodynamic interactions can be neglected. The hydrodynamic radii, which determine the hydrodynamic mobilities, do not play an explicit role in the following description.

The probe is dragged by an external force $F$ in a fixed direction in space, specified by the unit vector $\hat {\boldsymbol {{k}}}$ . By averaging over all possible configurations of the bath particles, it was shown by SB05 that the average probe velocity $\langle \boldsymbol {U} \rangle$ is given by

(2.1) \begin{align} \langle \boldsymbol {U} \rangle = \frac {D_a}{k_BT}\left [\hat {\boldsymbol {{k}}} F - nk_BT\oint _{|\boldsymbol {x}|=a+b} \hat {\boldsymbol {n}} g(\boldsymbol {x}/\ell )\, \mathrm {d} S\right ]. \end{align}

Here, $k_BT$ is the thermal energy and $g(\boldsymbol {x}/\ell )$ is the pair-distribution function – the conditional probability density for finding a bath particle at position $\boldsymbol {x}$ relative to the probe particle. The integration is carried out over an ‘exclusion sphere’, of radius $a+b$ , centred about the probe. The normalising length scale $\ell$ is to be specified later. Note that the pre-factor $D_a/k_BT$ is the Stokes–Einstein mobility in free solution.

By symmetry $\langle \boldsymbol {U} \rangle$ is in the direction of $\hat {\boldsymbol {{k}}}$ . SB05 introduced a Stokes drag interpretation for the probe motion,

(2.2) \begin{align} \hat {\boldsymbol {{k}}} F = \frac {\eta _{\textit{eff}}}{\eta } \frac {k_BT}{D_a}\langle \boldsymbol {U} \rangle , \end{align}

which serves to define an effective viscosity $\eta _{\textit{eff}}$ of the bath. Substitution into (2.1) yields, upon forming the dot product with $\hat {\boldsymbol {{k}}}$ ,

(2.3) \begin{align} \frac {\eta }{\eta _{\textit{eff}}} = 1-\frac {nk_BT}{F}\hat {\boldsymbol {{k}}}\boldsymbol {\cdot }\oint _{|\boldsymbol {x}|=a+b} \hat {\boldsymbol {n}} g(\boldsymbol {x}/\ell )\, \mathrm {d} S. \end{align}

SB05 also considered the companion problem where the probe is dragged at a fixed speed, say $U$ , in a fixed direction specified by the unit vector $\hat {\boldsymbol {{k}}}$ . Here they found that (2.1) is replaced by

(2.4) \begin{align} \hat {\boldsymbol {{k}}} U = \frac {D_a}{k_BT}\left [\langle \boldsymbol {F} \rangle - nk_BT\oint _{|\boldsymbol {x}|=a+b} \hat {\boldsymbol {n}} g(\boldsymbol {x}/\ell )\, \mathrm {d} S\right ], \end{align}

wherein $\langle \boldsymbol {F} \rangle$ is the average external force required to drag the probe, which by symmetry is in the direction of $\hat {\boldsymbol {{k}}}$ . The Stokes law interpretation used in this fixed-velocity case is

(2.5) \begin{align} \langle \boldsymbol {F} \rangle = \frac {\eta _{\textit{eff}}}{\eta }\frac {k_BT}{D_a} U \hat {\boldsymbol {{k}}}. \end{align}

Substitution into (2.4) yields, upon forming the dot product with $\hat {\boldsymbol {{k}}}$ ,

(2.6) \begin{align} \frac {\eta _{\textit{eff}}} {\eta }= 1+\frac {nD_a}{U}\hat {\boldsymbol {{k}}}\boldsymbol {\cdot }\oint _{|\boldsymbol {x}|=a+b} \hat {\boldsymbol {n}} g(\boldsymbol {x}/\ell )\, \mathrm {d} S. \end{align}

The rheological quantity of interest is the viscosity increment $\eta _{\textit{eff}}-\eta$ . For a dilute dispersion, the fractional increment for a fixed force is found from (2.3) as

(2.7) \begin{align} \frac {\eta _{\textit{eff}}-\eta } {\eta }=\frac {nk_BT}{F}\hat {\boldsymbol {{k}}}\boldsymbol {\cdot }\oint _{|\boldsymbol {x}|=a+b} \hat {\boldsymbol {n}} g(\boldsymbol {x}/\ell )\, \mathrm {d} S. \end{align}

Similarly, the fractional increment for a fixed velocity is found from (2.6) as

(2.8) \begin{align} \frac {\eta _{\textit{eff}}-\eta } {\eta }=\frac {nD_a}{U}\hat {\boldsymbol {{k}}}\boldsymbol {\cdot }\oint _{|\boldsymbol {x}|=a+b} \hat {\boldsymbol {n}} g(\boldsymbol {x}/\ell )\, \mathrm {d} S \end{align}

(which is not restricted to dilute systems). By defining the velocity scale

(2.9) \begin{align} U \stackrel {\text {def}}{=} \frac {D_a}{k_BT}F \end{align}

for the fixed-force problem, the fixed-force expression (2.7) coincides with the fixed-velocity expression (2.8). We shall use the latter.

At this stage a dimensionless notation is beneficial, using the position vector $\boldsymbol {r}=\boldsymbol {x}/\ell$ and the differential areal element $\mathrm {d}{A}=\mathrm {d} S / \ell ^2$ . The Péclet number is defined as

(2.10) \begin{align} P{\kern-1pt}e = \frac {U\ell }{D}, \end{align}

where $D$ is the diffusivity of a bath particle relative to the probe. Upon choosing $\ell$ as dimensional radius of the exclusion sphere,

(2.11) \begin{align} \ell =a+b, \end{align}

approximation (2.8) reads

(2.12) \begin{align} \frac {\eta _{\textit{eff}}-\eta } {\eta }=\frac {3\phi }{4\pi P{\kern-1pt}e}\frac {D_a}{D} \left(\frac {a+b}{b}\right)^3\hat {\boldsymbol {{k}}}\boldsymbol {\cdot }\oint _{|\boldsymbol {r}|=1} \hat {\boldsymbol {n}} g(\boldsymbol {r})\, \mathrm {d}{A}. \end{align}

Here $\phi =(4\pi b^3/3)n$ is the average volume fraction of the bath particles. With definition (2.9), (2.12) applies for both the fixed-force and fixed-velocity scenarios. As explained by SB05, $D=D_a+D_b$ for the fixed-force problem and $D=D_b$ for the fixed-velocity problem.

Defining

(2.13) \begin{align} {V} = \frac {3}{2\pi P{\kern-1pt}e}\hat {\boldsymbol {{k}}}\boldsymbol {\cdot }\oint _{|\boldsymbol {r}|=1} \hat {\boldsymbol {n}} g(\boldsymbol {r})\, \mathrm {d}{A}, \end{align}

the fractional increment (2.12) is factored as

(2.14) \begin{align} \frac {\eta _{\textit{eff}}-\eta } {\eta }=\frac {\phi }{2}\frac {D_a}{D}\left (\frac {a+b}{b}\right )^3 {V}. \end{align}

The quantity $V$ isolates the effect of external forcing upon the fractional viscosity increment. Its pre-factor in (2.14) is independent of that forcing.

3. Microstructure

The evaluation of the viscosity increment using (2.13)–(2.14) amounts to the calculation of the pair-distribution function $g$ . Using a reference frame moving with the probe, SB05 derived the problem governing $g$ in the dilute limit. It consists of the advection–diffusion equation

(3.1) \begin{align} \boldsymbol {\nabla } \boldsymbol {\cdot } (\boldsymbol {\nabla } g+ P{\kern-1pt}e \,\hat {\boldsymbol {{k}}} \,g) = 0, \end{align}

the no-flux condition

(3.2) \begin{align} \hat {\boldsymbol {n}} \boldsymbol {\cdot } ( \boldsymbol {\nabla } g+ P{\kern-1pt}e \,\hat {\boldsymbol {{k}}} \,g) = 0 \quad \text {at}\; |\boldsymbol {r}|=1 \end{align}

and the far-field condition

(3.3) \begin{align} \lim _{|\boldsymbol {r}|\to \infty }g=1. \end{align}

The above linear problem defines $g(\boldsymbol {r})$ as a (nonlinear) function of the single parameter $ P{\kern-1pt}e$ .

We employ here cylindrical coordinates $(\varpi ,\varphi ,z)$ , with the $z$ axis pointing in the direction $\hat {\boldsymbol {{k}}}$ and the origin at $\boldsymbol {r}=\boldsymbol {0}$ , as well as spherical polar coordinates $(r,\theta ,\varphi )$ with the polar angle $\theta$ measured from the positive $z$ axis. The exclusion sphere now coincides with $r=1$ . By axial symmetry, $g$ is independent of the azimuthal angle $\varphi$ . The advection–diffusion equation (3.1) is alternatively written

(3.4) \begin{align} \nabla ^2 g + P{\kern-1pt}e \frac {\partial g}{\partial z} = 0. \end{align}

Since the boundary condition (3.2) applies at $r=1$ , where $\hat {\boldsymbol {n}}=\hat {\boldsymbol {e}}_r$ , it simplifies to

(3.5) \begin{align} \frac {\partial g}{\partial r} + P{\kern-1pt}e \,g \cos \theta = 0\quad \text {at} \; r=1. \end{align}

Also, condition (3.3) becomes

(3.6) \begin{align} \lim _{r\to \infty }g=1. \end{align}

When using spherical coordinates, we write (with a slight abuse of notation) $g(r,\theta ; P{\kern-1pt}e)$ instead of $g(\boldsymbol {r})$ . Equation (3.4) then reads

(3.7) \begin{align} \frac {\partial ^2g}{\partial r^2} + \frac {2}{r}\frac {\partial g}{\partial r} + \frac {1}{r^2\sin \theta } \frac {\partial }{\partial \theta } \left (\sin \theta \, \frac {\partial g}{\partial \theta }\right )+ P{\kern-1pt}e\left (\cos \theta \frac {\partial g}{\partial r} - \frac {\sin \theta }{r}\frac {\partial g}{\partial \theta }\right ) =0. \end{align}

Equations (3.5)–(3.7) define a linear problem governing $g$ . With $\mathrm {d}{A}=2\pi \sin \theta \,\mathrm {d}\theta$ at $r=1$ , expression (2.13) simplifies to

(3.8) \begin{align} {V}( P{\kern-1pt}e) = \frac {3}{ P{\kern-1pt}e}\int _0^\pi g(1,\theta ; P{\kern-1pt}e)\sin \theta \cos \theta \, \mathrm {d} \theta . \end{align}

For future reference, we also write the governing equations in the cylindrical coordinates, where (with another abuse of notation) the pair-distribution function is referred to as $g(\varpi ,z; P{\kern-1pt}e)$ . Equation (3.4) then reads

(3.9) \begin{align} \frac {\partial ^2g}{\partial \varpi ^2} + \frac {1}{\varpi }\frac {\partial g}{\partial \varpi } + \frac {\partial ^2g}{\partial z^2} + P{\kern-1pt}e\, \frac {\partial g}{\partial z} =0, \end{align}

while condition (3.5) becomes

(3.10) \begin{align} \varpi \frac {\partial g}{\partial \varpi } + z\frac {\partial g}{\partial z} + P{\kern-1pt}e \,g z = 0\quad \text {at} \ \varpi ^2+z^2=1. \end{align}

4. From small to large $ P{\kern-1pt}e$

SB05 considered both small and large $ P{\kern-1pt}e$ . For small $ P{\kern-1pt}e$ , perturbing upon a uniform distribution of $g$ , they obtained

(4.1) \begin{align} g-1 \approx \frac { P{\kern-1pt}e}{2}\cos \theta , \end{align}

from which they found that

(4.2) \begin{align} \lim _{ P{\kern-1pt}e\to 0}{V}( P{\kern-1pt}e) =1. \end{align}

At the limit of large $ P{\kern-1pt}e$ , SB05 identified three asymptotic regions: a boundary layer of $O( P{\kern-1pt}e^{-1})$ width, where $g=O( P{\kern-1pt}e)$ , about the upstream portion of the exclusion sphere; an advection-dominated region outside it, where $g\approx 1$ ; and a wake downstream of the sphere, where $g\approx 0$ , bounded by the cylindrical surface $\varpi =1$ . Using the leading-order boundary-layer solution, SB05 found that

(4.3) \begin{align} \lim _{ P{\kern-1pt}e\to \infty }{V}( P{\kern-1pt}e) =\frac {1}{2}. \end{align}

To illustrate the force thinning from (4.2) to (4.3), SB05 also addressed the problem of arbitrary $ P{\kern-1pt}e$ . Thus, by transforming (3.4) to a modified Helmholtz equation, SB05 obtained a semi-numerical solution of $g$ using an eigenfunction expansion.

Our interest is in the approach to the large- $ P{\kern-1pt}e$ Newtonian plateau (4.3). As a prelude to our analysis, we illustrate in figure 1 the numerical solution of the problem (3.4)–(3.6) for $ P{\kern-1pt}e=100$ , using a finite-difference scheme. The scalings indicated in the figure hint towards the asymptotic investigation in the next sections. The numerical scheme itself is discussed in § 9.

Figure 1. Top: diagram of the microstructure problem (3.4)–(3.6) for the pair-distribution function $g$ outside the exclusion sphere $r=1$ , with arrows indicating the streaming flow in the probe-fixed reference frame. The $g$ contours (calculated numerically using the finite-difference method described in § 9 with $ P{\kern-1pt}e = 100$ ) illustrate the typical structure of the solution for large $ P{\kern-1pt}e$ , with the asymptotic estimates of $g$ derived in §§ 5, 6 and 8. Bottom: schematic illustrating the two-part non-uniform numerical grid. The dashed segment is the equator $\theta =\pi /2$ .

5. The limit $ P{\kern-1pt}e\gg 1$ : boundary-layer analysis

Henceforth, we consider the limit

(5.1) \begin{align} P{\kern-1pt}e\gg 1. \end{align}

At leading order, it follows from the governing equation (3.4) that $g$ is uniform on ‘streamlines’ parallel to $\hat {\boldsymbol {{k}}}$ , i.e.

(5.2) \begin{align} g \text { is a function of $\varpi $ alone}. \end{align}

It then follows from the far-field condition (3.6) that

(5.3) \begin{align} g\approx 1. \end{align}

This solution, however, cannot satisfy the no-flux condition (3.5) at leading order. A boundary layer is therefore formed about the exclusion sphere $r=1$ . In the terminology of matched asymptotic expansions (Hinch 1991) the solution (5.3) is regarded as an ‘outer’ one. We note that it is an exact solution of (3.4) and (3.6) at all algebraic orders. We therefore conclude that

(5.4) \begin{align} g=1+\exp , \end{align}

wherein ‘exp’ stands for exponentially small terms – that is, terms that are asymptotically smaller than any negative powers of $ P{\kern-1pt}e$ – which are triggered by asymptotic matching.

Following SB05, we examine a boundary layer of $O( P{\kern-1pt}e^{-1})$ width where $g$ is $O( P{\kern-1pt}e)$ . We employ the rescaled radial coordinate

(5.5) \begin{align} r = 1 + P{\kern-1pt}e^{-1}\zeta , \end{align}

and define

(5.6) \begin{align} g(r,\theta ; P{\kern-1pt}e) = {p}(\zeta ,\theta ; P{\kern-1pt}e). \end{align}

Upon using (5.5)–(5.6), the advection–diffusion equation (3.7) becomes

(5.7) \begin{align} \frac {\partial ^2{p}}{\partial \zeta ^2} & + \frac {\partial {p}}{\partial \zeta } \cos \theta + P{\kern-1pt}e^{-1}\left (2\frac {\partial {p}}{\partial \zeta } - \frac {\partial {p}}{\partial \theta } \sin \theta \right ) \nonumber \\ &+ P{\kern-1pt}e^{-2}\left [-2\zeta \frac {\partial {p}}{\partial \zeta } + \frac {1}{\sin \theta }\frac {\partial }{\partial \theta }\left (\sin \theta \,\frac {\partial {p}}{\partial \theta }\right ) + \zeta \sin \theta \,\frac {\partial {p}}{\partial \theta }\right ] + O( P{\kern-1pt}e^{-3}{p}) = 0, \end{align}

while the no-flux condition (3.5) reads

(5.8) \begin{align} \frac {\partial {p}}{\partial \zeta } + {p} \cos \theta = 0 \quad \text {at} \ \zeta = 0. \end{align}

In addition, $p$ must also satisfy

(5.9) \begin{align} \lim _{\zeta \to \infty } {p} = 1, \end{align}

representing asymptotic matching with (5.4).

The governing equation (5.7) is integrated once using condition (5.8) to yield

(5.10) \begin{align} \frac {\partial {p}}{\partial \zeta } + {p} \cos \theta = & P{\kern-1pt}e^{-1} \int _0^\zeta \left (-2\frac {\partial {p}}{\partial \zeta } + \frac {\partial {p}}{\partial \theta } \sin \theta \right ) \, \mathrm {d}\zeta \nonumber \\ + & P{\kern-1pt}e^{-2}\int _0^\zeta \left [2\zeta \frac {\partial {p}}{\partial \zeta } - \frac {1}{\sin \theta }\frac {\partial }{\partial \theta }\left (\sin \theta \,\frac {\partial {p}}{\partial \theta }\right ) - \zeta \sin \theta \,\frac {\partial {p}}{\partial \theta } \right ] \mathrm {d}{\zeta } + O( P{\kern-1pt}e^{-3}{p}) . \end{align}

In what follows, we solve (5.10) subject to (5.9). Since $p=O( P{\kern-1pt}e)$ , we expand it as

(5.11) \begin{align} {p}(\zeta ,\theta ; P{\kern-1pt}e)= P{\kern-1pt}e\, {p}_{-1}(\zeta ,\theta ) + {p}_0(\zeta ,\theta ) + P{\kern-1pt}e^{-1} {p}_1(\zeta ,\theta ) + O( P{\kern-1pt}e^{-2}), \end{align}

and solve order by order.

5.1. Evaluation of ${p}_{-1}$

At $\textrm {ord}( P{\kern-1pt}e)$ the problem posed by (5.10) and (5.9) reads

(5.12a,b) \begin{align} \frac {\partial {p}_{-1}}{\partial \zeta } + {p}_{-1} \cos \theta = 0, \quad \lim _{\zeta \to \infty } {p}_{-1} = 0. \end{align}

The solution of (5.12 a) is

(5.13) \begin{align} {p}_{-1}(\zeta ,\theta ) = f_{-1}(\theta )\mathrm {e}^{-\zeta \cos \theta }. \end{align}

To satisfy condition (5.12 b), the boundary layer must be restricted to the upstream hemisphere:

(5.14) \begin{align} 0 \lt \theta \lt \pi /2, \end{align}

whereby (5.12 b) is trivially satisfied. We note that

(5.15) \begin{align} \int _0^\infty {p}_{-1}(\zeta ,\theta ) \, \mathrm {d}{\zeta } = \frac {f_{-1}(\theta )}{\cos \theta }. \end{align}

5.2. Wake

In the absence of a boundary layer at the ‘back’ hemisphere, the outer solution downstream must satisfy the no-flux condition (3.10). The only possible solution of (5.2) which is compatible with that condition at leading order is the trivial one. It therefore follows that $g$ must vanish in a wake downstream of the exclusion sphere. Since the trivial solution satisfies (3.9) and (3.10) exactly, we conclude that (cf. (5.4))

(5.16) \begin{align} g = \exp \quad \text {for}\ \varpi \lt 1 \ `\text {behind' the exclusion sphere}. \end{align}

5.3. Solvability

We determine $f_{-1}(\theta )$ by a solvability condition which we derive at the next asymptotic order. At $\textrm {ord}(1)$ , equation (5.10) and condition (5.9) give, respectively,

(5.17) \begin{align} \frac {\partial {p}_0}{\partial \zeta } + {p}_0 \cos \theta = -2 \left .{p}_{-1}\right |_0^\zeta + \sin \theta \frac {\partial }{\partial \theta } \int _0^\zeta {p}_{-1}\,\mathrm {d}{\zeta }, \end{align}

(5.18) \begin{align} \lim _{\zeta \to \infty } {p}_0 = 1. \end{align}

We can derive the requisite solvability condition by evaluating (5.17) at $\zeta =\infty$ using (5.13)–(5.15) and (5.18) to obtain the first-order ordinary differential equation

(5.19) \begin{align} 1 = 2\frac {f_{-1}(\theta )}{\cos \theta } + \tan \theta \,\frac {\mathrm {d} }{\mathrm {d} \theta }\left [\frac {f_{-1}(\theta )}{\cos \theta }\right ], \end{align}

whose general solution is

(5.20) \begin{align} \frac {f_{-1}(\theta )}{\cos \theta } = \frac {1}{2} + \frac {C_{-1}}{\sin ^2\theta }. \end{align}

Requiring regularity at $\theta =0$ necessitates $C_{-1}=0$ . Thus, $f_{-1}(\theta )= (1/2)\cos \theta $ . Substitution into (5.13) gives

(5.21) \begin{align} {p}_{-1}(\zeta ,\theta ) = \frac {\cos \theta }{2}\,\mathrm {e}^{-\zeta \cos \theta } . \end{align}

We note that

(5.22) \begin{align} \int _0^\zeta {p}_{-1}(\zeta ,\theta ) \,\mathrm {d}{\zeta } = \frac {1-\mathrm {e}^{-\zeta \cos \theta }}{2} \end{align}

and

(5.23) \begin{align} \int _0^\infty \zeta {p}_{-1}(\zeta ,\theta ) \,\mathrm {d}{\zeta } = \frac {1}{2\cos \theta }. \end{align}

Substitution into (3.8) of (5.11) and (5.21) reproduces (4.3), already obtained by SB05. (SB05 did not employ a solvability approach; rather, they obtained the counterpart of the present $f_{-1}$ using an integral flux argument. Note that (4.3) was derived in Appendix D of Khair & Brady (2006) using a solvability approach.) To go beyond this leading-order result we need the $\textrm {ord}(1)$ term ${p}_0$ .

5.4. Evaluation of ${p}_0$

We can now actually solve for $ {p}_0$ . Substituting (5.21)–(5.22) into (5.17) yields

(5.24) \begin{gather} \frac {\partial {p}_0}{\partial \zeta } + {p}_0 \cos \theta = (1-\mathrm {e}^{-\zeta \cos \theta })\cos \theta - \frac {\zeta }{2}\mathrm {e}^{-\zeta \cos \theta }\sin ^2\theta . \end{gather}

Using the integrating factor $\mathrm {e}^{\zeta \cos \theta }$ we find that the general solution of this first-order equation is

(5.25) \begin{gather} {p}_0(\zeta ,\theta ) = f_0(\theta )\mathrm {e}^{-\zeta \cos \theta } + 1 - \zeta \mathrm {e}^{-\zeta \cos \theta }\cos \theta - \frac {\zeta ^2}{4}\mathrm {e}^{-\zeta \cos \theta }\sin ^2\theta . \end{gather}

Given (5.14), condition (5.18) is trivially satisfied.

To derive a solvability condition on $f_0$ , we continue to the next asymptotic order. At $\textrm {ord}( P{\kern-1pt}e^{-1})$ , we find from (5.10) that

(5.26) \begin{align} \frac {\partial {p}_1}{\partial \zeta } + {p}_1 \cos \theta & = -2\left .{p}_0\right |_0^\zeta + \sin \theta \frac {\partial }{\partial \theta } \int _0^\zeta {p}_0\,\mathrm {d}\zeta \nonumber \\ & \quad + 2 \int _0^\zeta \zeta \frac {\partial {p}_{-1}}{\partial \zeta } \,\mathrm {d}\zeta - \frac {1}{\sin \theta }\frac {\partial }{\partial \theta }\left (\sin \theta \,\frac {\partial }{\partial \theta } \int _0^\zeta {p}_{-1} \mathrm {d}\zeta \right ) \nonumber \\ &\quad - \sin \theta \,\frac {\partial }{\partial \theta } \int _0^\zeta \zeta {p}_{-1}\, \mathrm {d}{\zeta }. \end{align}

This equation is subject to

(5.27) \begin{align} \lim _{\zeta \to \infty } {p}_1 = 0, \end{align}

obtained from (5.9).

We now apply the limit $\zeta \to \infty$ to (5.26), where its left-hand side vanishes by virtue of (5.27). The various terms on the right-hand side of (5.26) are evaluated in that limit using (5.21)–(5.22) and (5.25). The first term is simply $2f_0(\theta )$ . In evaluating the second term, interchanging the order of integration and differentiation requires some care. Making use of

(5.28) \begin{align} \int _0^\infty [{p}_0(\zeta ,\theta ) - 1]\,\mathrm {d}\zeta = \frac {f_0(\theta )}{\cos \theta } - \frac {1}{\cos \theta } - \frac {\sin ^2\theta }{2\cos ^3\theta }, \end{align}

this term is obtained as

(5.29) \begin{align} \sin \theta \,\frac {\mathrm {d} }{\mathrm {d} \theta } \left [\frac {f_0(\theta )}{\cos \theta } - \frac {1}{\cos \theta } - \frac {\sin ^2\theta }{2\cos ^3\theta }\right ]. \end{align}

Integration by parts of the third term gives $-1$ . The fourth term vanishes, while the fifth, using (5.23), is $ ({1}/{2})\tan ^2\theta$ . Adding these up we obtain the first-order ordinary differential equation

(5.30) \begin{align} \sin \theta \,\frac {\mathrm {d} }{\mathrm {d} \theta } \left [\frac {f_0(\theta )}{\cos \theta }\right ] + 2f_0(\theta ) = \frac {1}{\cos ^2\theta } + \frac {3\sin ^2\theta }{2\cos ^4 \theta }, \end{align}

whose general solution is

(5.31) \begin{align} \frac {f_0(\theta )}{\cos \theta } = \frac {C_0}{\sin ^2\theta } + \frac {1}{2\cos ^3\theta }. \end{align}

Regularity at $\theta =0$ yields $C_0 = 0$ . Substitution back into (5.25) then furnishes the result

(5.32) \begin{align} {p}_0(\zeta ,\theta )=1 + \mathrm {e}^{-\zeta \cos \theta }\left (\frac {1}{2\cos ^2\theta } - \zeta \cos \theta - \frac {\zeta ^2}{4}\sin ^2\theta \right ). \end{align}

Thus, the boundary-layer solution is

(5.33) \begin{align} {p}(\zeta ,\theta ; P{\kern-1pt}e) & = P{\kern-1pt}e\,\frac {\cos \theta }{2}\mathrm {e}^{-\zeta \cos \theta } \nonumber \\ & \qquad + \left \{1 + \mathrm {e}^{-\zeta \cos \theta }\left (\frac {1}{2\cos ^2\theta } - \zeta \cos \theta - \frac {\zeta ^2}{4}\sin ^2\theta \right )\right \} + O( P{\kern-1pt}e^{-1}). \end{align}

In particular,

(5.34) \begin{align} {p}(0,\theta ; P{\kern-1pt}e) = P{\kern-1pt}e \frac {\cos \theta }{2} + \left (1 + \frac {1}{2\cos ^2\theta }\right ) + O( P{\kern-1pt}e^{-1}). \end{align}

5.5. Boundary-layer breakdown

The solution (30) proposed by SB05 for the boundary-layer problem, while not originally posed as an asymptotic expansion, is eventually recast in (35) as the sum of an $\textrm {ord}( P{\kern-1pt}e)$ term and an $\textrm {ord}(1)$ term. Ignoring slight typos, the first term agrees with the present (5.21). (The exponent appearing in equation (30) of SB05 is dimensional; apparently, the factor $U/D$ there needs to be removed. The exponent appearing in equation (35) of SB05, $-Uz\cos\theta/D$ , while dimensionless, does not agree with the correct one $-\zeta\cos\theta$ . For it to agree, $\cos\theta$ must be removed and $z$ needs to be replaced with the projection $z-(a+b)\cos\theta$ . In any event, these typos do not affect equation (36) of SB05).

The uniform $\textrm {ord}(1)$ term in equation (35) of SB05 is tantamount, in the present notation, to ${p}_0=1$ . This term does not satisfy condition (5.8) at $\textrm {ord}(1)$ , and is therefore erroneous. The error is carried to equation (36) of SB05 – the counterpart of the present (5.34). The oversight of SB05 (propagated into the microrheology literature – see (26) of Puertas & Voigtmann (2014)) may be traced to their original (30), which seems like an attempt to superpose upon the boundary-layer solution a constant term that satisfies the approach to the uniform outer solution.

The oversight did not affect the remaining analysis of SB05, who were only interested in the leading-order viscosity increment. In going beyond that leading-order calculation, however, the difference between equation (36) of SB05 and the present (5.34) is striking. Indeed, a naive attempt to substitute (5.34) into (3.8) so as to obtain the asymptotic correction to (4.3) would fail: due to the $1/\cos ^2\theta$ term in (5.34), the integral appearing in (3.8) diverges logarithmically at $\textrm {ord}(1)$ . This failure indicates a breakdown of expansion (5.11), and with it the very boundary-layer structure, as $\theta \nearrow \pi /2$ . It suggests the asymptotic formation of a transition region about the equator $\theta = \pi /2$ , which we address next.

6. Equatorial region

With the first term in the parentheses of (5.32) diverging as $\theta \nearrow \pi /2$ , the boundary-layer solution (5.11) breaks down near the equator. Introducing for convenience the latitude angle

(6.1) \begin{gather} \vartheta = \frac {\pi }{2} - \theta \ll 1, \end{gather}

we observe that the two asymptotic terms in (5.34) behave there as $ ({1}/{2}) P{\kern-1pt}e\,\vartheta$ and $ ({1}/{2})\vartheta ^{-2}$ , respectively. Hence, these terms become comparable at

(6.2) \begin{align} \vartheta = O( P{\kern-1pt}e^{-1/3}), \end{align}

for which

(6.3) \begin{align} g = O( P{\kern-1pt}e^{2/3}). \end{align}

Considering now (5.21) and noting that $\zeta \cos \theta \sim \zeta \vartheta$ near the equator, we see that the boundary layer ‘expands’ in the radial direction to $\zeta =O( P{\kern-1pt}e^{1/3})$ , or, equivalently (recall (5.5)), to

(6.4) \begin{align} r-1 = O( P{\kern-1pt}e^{-2/3}). \end{align}

We therefore investigate an equatorial region, of $O( P{\kern-1pt}e^{-1/3})$ angular extent about the equator and $O( P{\kern-1pt}e^{-2/3})$ radial extent about $r=1$ , where the pair-distribution function is $O( P{\kern-1pt}e^{2/3})$ (see figure 1).

6.1. Stretched cylindrical coordinates

In order to gain physical intuition for the transition region, we consider first stretched cylindrical coordinates on the scales (6.2) and (6.4):

(6.5) \begin{align} t = - P{\kern-1pt}e^{1/3} z, \quad y = P{\kern-1pt}e^{2/3}(\varpi -1). \end{align}

Defining

(6.6) \begin{align} g(\varpi ,z; P{\kern-1pt}e) = {G}(y,t; P{\kern-1pt}e), \end{align}

equation (3.9) becomes

(6.7) \begin{align} \frac {\partial {G}}{\partial t}=\frac {\partial ^2{G}}{\partial y^2 } + P{\kern-1pt}e^{-2/3} \left (\frac {\partial {G}}{\partial y} + \frac {\partial ^2{G}}{\partial t^2 } \right )+ O \big( P{\kern-1pt}e^{-4/3} \big). \end{align}

We interpret this as a perturbed one-dimensional heat equation, with the downstream coordinate being a time-like variable. With the surface $\varpi ^2+z^2=1$ becoming

(6.8) \begin{align} y = -\frac {t^2}{2} - P{\kern-1pt}e^{-2/3} \frac {t^4}{8} + O \big( P{\kern-1pt}e^{-4/3} \big), \end{align}

condition (3.10) reads, on that surface,

(6.9) \begin{align} \frac {\partial {G}}{\partial y} - t {G} + P{\kern-1pt}e^{-2/3} \left (y \frac {\partial {G}}{\partial y} + t \frac {\partial {G}}{\partial t}\right ) = 0. \end{align}

We interpret this as a ‘parabolically’ moving no-flux boundary, with perturbations. Equations (6.7)–(6.9) are subject to the condition

(6.10) \begin{align} \lim _{y \to \infty } {G} = 1, \end{align}

representing asymptotic matching with (5.4).

Following (6.3), we pose the asymptotic expansion

(6.11) \begin{align} {G}(y,t; P{\kern-1pt}e) = P{\kern-1pt}e^{2/3} {G}_{-2}(y,t) + O(1). \end{align}

At $O( P{\kern-1pt}e^{2/3})$ , we then find from (6.7)–(6.10) that

(6.12) \begin{align} \frac {\partial {G}_{-2}}{\partial t} = \frac {\partial ^2{G}_{-2}}{\partial y^2} \quad \text {for}\ y \gt -\frac {t^2}{2}, \\[-24pt] \nonumber \end{align}

(6.13) \begin{align} \frac {\partial {G}_{-2}}{\partial y} = t {G}_{-2} \quad \text {at}\ y = -\frac {t^2}{2}, \\[-24pt] \nonumber \end{align}

(6.14) \begin{align} \lim _{t \to \infty } {G}_{-2} = 0. \end{align}

The partial differential equation (6.12) is parabolic with a time-like coordinate $t$ . Hence, we require an ‘initial condition’, which is obtained from matching with the boundary-layer solution (5.11) in the limit $t\to -\infty$ . Using the leading-order boundary-layer solution (5.21) in conjunction with (6.26), this requirement yields the limiting behaviour

(6.15) \begin{align} {G}_{-2} \sim -\frac {t}{2} \mathrm {e}^{ty+t^3/2} \quad \text {as} \ t\to -\infty \ \text {with} \ \text {$ty+\frac {t^3}{2}$ fixed}. \end{align}

Further, we can derive a conservation law for ${G}_{-2}$ as follows. Integration of (6.12) from $y=-t^2/2$ to $\infty$ yields, upon making use of (6.13)–(6.14),

(6.16) \begin{align} \int _{-t^2/2}^\infty \frac {\partial {G}_{-2}}{\partial t} \, \mathrm {d} y + t {G}_{-2}(-t^2/2,t) = 0, \end{align}

or, making use of Leibniz’s rule,

(6.17) \begin{align} \frac {\mathrm {d} }{\mathrm {d} t} \int _{-t^2/2}^\infty {G}_{-2} (y,t) \, \mathrm {d} y = 0, \end{align}

so that $\int _{-t^2/2}^\infty {G}_{-2}(y,t) \, \mathrm {d} y$ is independent of $t$ ; substitution of (6.15) then yields

(6.18) \begin{align} \int _{-t^2/2}^\infty {G}_{-2}(y,t)\, \mathrm {d} y = \frac {1}{2} \quad \text {for all }t. \end{align}

The leading-order equatorial transport (6.12)–(6.13) can be interpreted as a one-dimensional heat-like equation with a no-flux boundary that first ( $t\lt 0$ ) moves in and accumulates the pair-distribution probability in its way, and then ( $t\gt 0$ ) moves out again, leaving behind a diffusing peak. Thus, the ‘local’ behaviour of the transition-region solution far downstream, unaffected by the boundary condition that has moved far away, is that of the unbounded heat equation:

(6.19a,b) \begin{align} {G}_{-2}(y,t) \sim \frac {Q}{(4\pi t)^{1/2}} \mathrm {e}^{-{y^2}/{4t}} \quad \text {as} \ t\to \infty \ \text {with} \ y^2/t \ \text{fixed}. \end{align}

From (6.18) we then find

(6.20) \begin{align} Q=\frac {1}{2}. \end{align}

6.2. Stretched spherical coordinates

With the differential equation (6.12) adopting a simple form, the formulation in terms of the $(y,t)$ coordinates is physically illuminating. For numerical and asymptotic calculations, however, it is preferable for the exclusion-sphere boundary to coincide with a coordinate surface. We therefore employ another local formulation in term of stretched spherical coordinates $(X,Y)$ , defined by

(6.21) \begin{align} r = 1 + P{\kern-1pt}e^{-2/3} Y, \quad \theta = \frac {\pi }{2} - P{\kern-1pt}e^{-1/3} X, \end{align}

and replace (6.6) with

(6.22) \begin{align} g(r,\theta ; P{\kern-1pt}e) = \tilde G(X,Y; P{\kern-1pt}e). \end{align}

Introducing (6.21)–(6.22) into (3.7) and (3.5) gives

(6.23) \begin{align} \frac {\partial ^2\tilde G}{\partial Y^2} + X \frac {\partial \tilde G}{\partial Y} + \frac {\partial \tilde G}{\partial X} +O( P{\kern-1pt}e^{-2/3}\tilde G) = 0 \quad \text {for} \ Y \gt 0 , \\[-24pt] \nonumber \end{align}

(6.24) \begin{align} \frac {\partial \tilde G}{\partial Y} + X\tilde G + O( P{\kern-1pt}e^{-2/3}\tilde G) = 0 \quad \text {at} \ Y = 0. \end{align}

These equations are subject to the condition

(6.25) \begin{align} \lim _{Y \to \infty } \tilde G = 1, \end{align}

representing asymptotic matching with (5.4).

Similarly to (6.11), we pose the asymptotic expansion

(6.26) \begin{align} \tilde G(X,Y; P{\kern-1pt}e) = P{\kern-1pt}e^{2/3} \tilde G_{-2}(X,Y) + O(1). \end{align}

At $O( P{\kern-1pt}e^{2/3})$ , we then find from (6.23)–(6.25) that

(6.27) \begin{align} \frac {\partial ^2\tilde G_{-2}}{\partial Y^2 } + X \frac {\partial \tilde G_{-2}}{\partial Y} + \frac {\partial \tilde G_{-2}}{\partial X} = 0 \quad \text {for}\ Y \gt 0, \end{align}

(6.28) \begin{align} \frac {\partial \tilde G_{-2}}{\partial Y} + X\tilde G_{-2} = 0 \quad \text {at}\ Y = 0, \\[-24pt] \nonumber \end{align}

(6.29) \begin{align} \lim _{Y \to \infty } \tilde G_{-2} = 0. \end{align}

As before, the ‘initial condition’ is obtained from matching with the boundary-layer solution (5.11) in the limit $X \to \infty$ . This gives

(6.30) \begin{align} \tilde G_{-2} \sim \frac {X}{2}\mathrm {e}^{-XY} \quad \text {as} \ X \to \infty \ \text {with} \ \text {$XY$ fixed}. \end{align}

We again derive a conservation law. Integration of (6.27) from $Y=0$ to $\infty$ yields, upon making use of (6.28)–(6.29),

(6.31) \begin{align} \quad \int _0^\infty \frac {\partial \tilde G_{-2}}{\partial X}\,\mathrm {d} Y=0 \quad \text {for all }X. \end{align}

Interchanging the order of integration and imposing the initial condition (6.30) reveals that

(6.32) \begin{align} \int _0^\infty \tilde G_{-2}(X,Y)\, \mathrm {d} Y = \frac {1}{2} \quad \text {for all }X. \end{align}

We note that the problem (6.27)–(6.29) is recovered from (6.12)–(6.14) by the transformation

(6.33) \begin{align} t \mapsto -X, \quad y \mapsto Y - \frac {X^2}{2}, \quad G_{-2} \mapsto \tilde {G}_{-2}. \end{align}

Thus, (6.30) and (6.32) directly follow from (6.15) and (6.18), respectively.

Performing a local analysis in the limit $X \to \infty$ with $XY$ fixed, we find using (6.27)–(6.29) that (6.30) is refined to

(6.34) \begin{align} \tilde G_{-2} = \left [\frac {X}{2} + \left (\frac {1}{2X^2} - \frac {Y^2}{4}\right ) + o(X^{-2})\right ]\mathrm {e}^{-XY} \quad \text {as} \ X \to \infty \ \text {with} \ \text {$XY$ fixed}. \end{align}

Note that the second term in the brackets matches the boundary-layer correction (5.32).

We solve the universal problem (6.27)–(6.30) numerically using a finite-difference method. The $Y$ direction is discretised using an exponentially stretched grid with relative resolution 0.1 %, and the solution is evolved in the time-like coordinate X from the initial condition (6.34) at $X = 10$ in steps of size $0.001$ to $X = -10$ . The resulting variation of $\tilde G_{-2}$ with $Y$ at certain values of $X$ in that range is displayed in figure 2.

Figure 2. Variation of the leading-order equatorial deformation $\tilde G_{-2}$ with $Y$ for the indicated positive values of $X$ , as obtained from the numerical solution of the universal problem. The inset shows the comparable variations for $X\leqslant 0$ . The dotted curve depicts the far-field solution (6.19)–(6.20), transformed using (6.33) for $X=-10$ .

In terms of the $(X,Y)$ coordinates, the local far-field solution (6.19)–(6.20) represents a diffusing Gaussian at large negative $X$ , centred about $Y=X^2/2$ with peak value $(-16\pi X)^{-1/2}$ . This explains the pulse-like behaviour observed in the inset of figure 2, as illustrated there for $X=-10$ .

Both the refined behaviour (6.34) and the numerical results for $\tilde G_{-2}$ at $Y=0$ are required in the asymptotic evaluation of the effective viscosity, which is performed next.

7. Evaluation of $V$

We can now evaluate $V$ at large $ P{\kern-1pt}e$ , going beyond (4.3). Given (5.16), the integral in (3.8) only needs to be taken up to a small angle past the equator $\theta =\pi /2$ . Recalling (see (6.2)) that the angular extent of the equatorial region is $O( P{\kern-1pt}e^{-1/3})$ , we may replace (3.8) by

(7.1) \begin{align} {V}( P{\kern-1pt}e) = \frac {3}{ P{\kern-1pt}e}\int _0^{\frac {\pi }{2}+\varDelta } g(1,\theta ; P{\kern-1pt}e)\sin \theta \cos \theta \, \mathrm {d} \theta , \end{align}

where the choice

(7.2) \begin{align} P{\kern-1pt}e^{-1/3}\ll \varDelta \ll 1 \end{align}

ensures that the entire contribution of the equatorial region is covered.

On the upstream hemisphere $0 \lt \theta \lt \pi /2$ , the boundary-layer solution is given by (cf. (5.34))

(7.3) \begin{align} g(1,\theta ; P{\kern-1pt}e) = P{\kern-1pt}e \frac {\cos \theta }{2} + \left (1 + \frac {1}{2\cos ^2\theta }\right ) + O( P{\kern-1pt}e^{-1}). \end{align}

Using (6.21)–(6.22) and the equatorial-region solution (6.26), we find that near the equator

(7.4) \begin{gather} g(1,\theta ; P{\kern-1pt}e) = P{\kern-1pt}e^{2/3} \tilde G_{-2}\left (X= P{\kern-1pt}e^{1/3}({\pi }/{2}-\theta ),Y=0\right ) + O(1). \end{gather}

Our goal is to evaluate the integral appearing in (7.1) using these two complementary approximations. The $\textrm {ord}( P{\kern-1pt}e)$ contribution to that integral is obtained as $ P{\kern-1pt}e/6$ from integrating the leading term in (7.3) over the upstream hemisphere (5.14). This contribution reproduces the leading-order approximation (4.3).

In going beyond $\textrm {ord}( P{\kern-1pt}e)$ we expect two $\textrm {ord}(1)$ contributions to the integral. The first originates in the $\textrm {ord}(1)$ correction in (7.3), integrated over the upstream hemisphere. The second arises from the leading $\textrm {ord}( P{\kern-1pt}e^{2/3})$ term in the equatorial-region expansion (7.4): indeed, we note that $\cos \theta = \textrm {ord}( P{\kern-1pt}e^{-1/3})$ there, as is the extent of the integration domain.

Due to the diverging term in the $\textrm {ord}(1)$ correction in (7.3), the associated hemispherical integral diverges towards the equator; this divergence is cut off by the equatorial-region solution: see indeed (6.34). We therefore anticipate (Hinch 1991) that the $\textrm {ord}(1)$ contribution incorporates logarithms of $ P{\kern-1pt}e$ .

Employing the cutoff-like parameter (cf. (7.2))

(7.5) \begin{align} P{\kern-1pt}e^{-1/3}\ll \chi \ll 1, \end{align}

we split the integral appearing in (7.1) at $\theta = \pi /2-\chi$ . Thus, we write it as

(7.6) \begin{align} \left (\int _0^{\frac {\pi }{2}-\chi } + \int _{\frac {\pi }{2}-\chi }^{\frac {\pi }{2}+\varDelta }\right )g(1,\theta ; P{\kern-1pt}e)\sin \theta \cos \theta \, \mathrm {d} \theta \stackrel {\text {def}}{=} I_{\textit{hemisphere}} + I_{\textit{equator}} . \end{align}

We now proceed to evaluate the integral to $\textrm {ord}(1)$ in the limit $ P{\kern-1pt}e\gg 1$ . The result must be independent of both $\varDelta$ and $\chi$ .

Using the boundary-layer approximation (7.3) we readily find the ‘hemispherical’ contribution:

(7.7) \begin{align} I_{\textit{hemisphere}}\sim P{\kern-1pt}e\frac {1 - \sin ^3\chi }{6} + \frac {\cos ^2\chi - \ln\sin \chi }{2}, \end{align}

which becomes, upon making use of (7.5),

(7.8) \begin{align} I_{\textit{hemisphere}} \sim P{\kern-1pt}e \frac {1 - \chi ^3}{6} + \frac {1 - \ln \chi }{2}. \end{align}

The asymptotic error in (7.8) is algebraically small, i.e. smaller than some negative power of $ P{\kern-1pt}e$ . We identify the above-mentioned $ P{\kern-1pt}e/6$ term, which reproduces (4.3). It is accompanied by the $\textrm {ord}(1)$ term $1/2$ , as well as intermediate terms (note that $ P{\kern-1pt}e\, \chi ^3\gg 1$ ) that depend upon $\chi$ . These must be cancelled out by the ‘equatorial’ contribution.

Using (6.21), the equatorial contribution is

(7.9) \begin{align} I_{\textit{equator}} \sim \int _{-\textit { Pe}^{1/3}\varDelta }^{\textit { Pe}^{1/3}\chi } \tilde G_{-2}(X,0) X \, \mathrm {d}{X}. \end{align}

Since the attenuation towards the wake is attained exponentially fast (recall (6.19)), we simply replace $ P{\kern-1pt}e^{1/3}\varDelta$ by $\infty$ . The upper integration limit, on the other hand, must be handled with care. Since $\tilde G_{-2}(X,0)$ actually diverges as $X/2$ at large $X$ (see (6.30)), it is clearly forbidden to replace $ P{\kern-1pt}e^{1/3}\chi$ by $\infty$ . We accordingly subtract the asymptotic behaviour $X/2$ . But this is insufficient, since the $\textrm {ord}(X^{-2})$ correction to $X/2$ (see (6.34)) would result in a logarithmically divergent integral. We therefore need to subtract the two-term far-field behaviour (cf. (6.34)):

(7.10) \begin{align} \tilde G_{-2}(X,0) \sim \frac {X}{2} + \frac {1}{2X^2} \quad \text {as} \ X \to \infty . \end{align}

To avoid a spurious singularity at $X=0$ we subtract (7.10) only for $X\gt 1$ . We therefore obtain

(7.11) \begin{align} I_{\textit{equator}} & \sim \int _{-\infty }^{\textit { Pe}^{1/3}\chi } \left [ \tilde G_{-2}(X,0) - \left (\frac {X}{2} + \frac {1}{2X^2}\right )H(X-1)\right ]X \, \mathrm {d}{X} \nonumber \\ & + \int _1^{\textit { Pe}^{1/3} \chi } \left (\frac {X^2}{2} + \frac {1}{2X}\right )\, \mathrm {d}{X} , \end{align}

wherein $H$ is the Heaviside step function.

With the above form, the integrand of the first integral decays faster than $1/X$ as $X\to \infty$ . It is therefore legitimate to replace the endpoint $ P{\kern-1pt}e^{1/3}\chi$ of that integral by $\infty$ . Evaluating the second integral explicitly we obtain

(7.12) \begin{align} I_{\textit{equator}} \sim \mathcal J+ \frac { P{\kern-1pt}e\,\chi ^3 - 1}{6} + \frac {1}{2}\ln \big( P{\kern-1pt}e^{1/3} \chi \big ), \end{align}

wherein

(7.13) \begin{align} \mathcal J= \int _{-\infty }^\infty \left [ \tilde G_{-2}(X,0) - \left (\frac {X}{2} + \frac {1}{2X^2}\right )H(X-1)\right ]X \, \mathrm {d}{X} \end{align}

is a pure number which is defined by the universal problem governing $\tilde G_{-2}$ . Making use of the numerical solution of that problem yields $\mathcal J=-0.1591$ .

Adding up (7.8) and (7.12), we conclude that

(7.14) \begin{align} I_{\textit{hemisphere}} + I_{\textit{equator}} \sim \frac { P{\kern-1pt}e + \ln P{\kern-1pt}e +6\mathcal J+2}{6}. \end{align}

where, as required, the dependence upon $\chi$ has cancelled out. Substitution into (7.1) yields

(7.15) \begin{align} {V}( P{\kern-1pt}e) \sim \frac {1}{2} + \frac {\ln P{\kern-1pt}e + 6\mathcal J+2}{2 P{\kern-1pt}e} \quad \text {as} \ P{\kern-1pt}e\to \infty . \end{align}

8. Transition layer

The jump at $\varpi =1$ between (5.4) and (5.16) necessitates a transition layer which emerges at the equatorial region and propagates downstream to $z\lt 0$ . Since the transition layer is detached from the sphere $r=1$ , it does not directly affect the microrheology – recall (3.8). Nonetheless, its analysis is provided here for completeness.

The transition layer widens by diffusion, so on an $O(1)$ axial length scale we expect the layer to reach an $O( P{\kern-1pt}e^{-1/2})$ thickness. We accordingly employ the stretched coordinate $\rho$ , defined by

(8.1) \begin{align} \varpi = 1 + P{\kern-1pt}e^{-1/2}\rho . \end{align}

For convenience we also employ $\tau =-z$ which is positive downstream of the exclusion sphere. Defining

(8.2) \begin{align} g(\varpi ,z; P{\kern-1pt}e) = h(\rho ,\tau ; P{\kern-1pt}e), \end{align}

equation (3.9) becomes

(8.3) \begin{align} \frac {\partial h}{\partial \tau } =\frac {\partial ^2h}{\partial \rho ^2} + P{\kern-1pt}e^{-1/2} \frac {\partial h}{\partial \rho } + O( P{\kern-1pt}e^{-1}h). \end{align}

In addition to (8.3), $h$ needs to satisfy the matching conditions

(8.4a,b) \begin{align} \lim _{\rho \to \infty }h = 1, \quad \lim _{\rho \to -\infty }h = 0, \end{align}

representing the approach to (5.4) and (5.16), respectively.

Like in the equatorial region, the governing equation (8.3) is parabolic and requires an initial condition. This comes from matching upstream to the large- $t$ behaviour of the equatorial region, but some progress can also be made by instead considering the integral balance:

(8.5) \begin{align} \int _{-\infty }^{\infty } \left [ h-H(\rho )- P{\kern-1pt}e^{-1}\frac {\partial h}{\partial \tau }\right ] \left(1+ P{\kern-1pt}e^{-1/2}\rho \right) \,\mathrm {d}\rho =\frac {1}{2} P{\kern-1pt}e^{1/2} \quad \forall \tau \gt 0, \end{align}

which we derive in Appendix A. It is evident from (8.5) that $h=\textrm {ord}( P{\kern-1pt}e^{1/2})$ . We therefore write

(8.6) \begin{align} h(\rho ,t; P{\kern-1pt}e) \sim P{\kern-1pt}e^{1/2}\,h_{-1} (\rho ,\tau ). \end{align}

At leading order, $\textrm {ord}( P{\kern-1pt}e^{1/2})$ , we find from (8.3)–(8.5) that

(8.7a,b,c) \begin{align} \frac {\partial h_{-1}}{\partial \tau } = \frac {\partial ^2h_{-1}}{\partial \rho ^2}, \quad \lim _{\rho \to \pm \infty }h_{-1} = 0, \quad \int _{-\infty }^\infty h_{-1}(\rho ,\tau ) \,\mathrm {d}{\rho } = \frac {1}{2} , \end{align}

for all $\tau \gt 0$ . Note that the normalisation condition (8.7 c) constitutes a continuation of (6.18).

A solution of (8.7), proportional to the Green’s function of the one-dimensional heat equation, is

(8.8) \begin{align} h_{-1} = \frac {\mathrm {e}^{-(\rho -\varPhi )^2/4\tau }}{2(4\pi \tau )^{1/2}}, \end{align}

where $\varPhi$ is an arbitrary constant. Asymptotic matching of (8.6) and (8.8) with the equatorial-region expansion (6.11) using (6.19)–(6.20) yields $\varPhi =0$ .

While (8.3) and (8.5) appear to suggest that the leading-order correction to (8.6) is $\textrm {ord}(1)$ , matching with the equatorial region actually triggers asymptotically larger homogeneous ‘eigensolutions’. Indeed, a detailed analysis (not shown here) reveals an $\textrm {ord}( P{\kern-1pt}e^{1/3})$ leading-order correction, which can be interpreted as an asymptotically small shift $\varPhi = O( P{\kern-1pt}e^{-1/6})$ of the Green’s function (8.8); this, in turn, originates from an $\textrm {ord}(1)$ shift of $y$ in (6.19) – recall indeed the slight mismatch in the inset of figure 2.

The analysis in this section revolved around the interaction between transport on the $O(1)$ axial length scale – corresponding to the size of that sphere – and on the equatorial length $O( P{\kern-1pt}e^{-1/3})$ . In addition, the solution (8.8) introduces the much longer $O( P{\kern-1pt}e)$ restoration length where the transition layer has grown to $O(1)$ thickness thus filling the wake and reinstating the outer solution (5.4).

9. Comparison with a finite-difference solution

For the purpose of verification, we have constructed a numerical solution of the full equations (3.4)–(3.6) for various values of $ P{\kern-1pt}e$ using a finite-difference scheme. In order to resolve the boundary-layer structures at large $ P{\kern-1pt}e$ , we use a combination of two non-uniform grids as shown in figure 1 (bottom). A spherical grid is used upstream and in a small downstream region, with radial refinement near the surface of the exclusion sphere and azimuthal refinement near the poles and equator. A cylindrical grid is used downstream away from the exclusion sphere, with axial refinement near the equator and radial refinement near $\varpi = 1$ . Interpolation is used across the interface between the two grids. The resulting equations are solved in Matlab using the built-in linear solver.

The numerical scheme is second order in the spatial resolution; with $1\,\%$ relative resolution (i.e. 20 times finer than the sample grid in figure 1) we have verified (by comparison with a finer grid) that the relative error is less than $10^{-4}$ for $ P{\kern-1pt}e$ up to $10^6$ . With that relative error, the $\textrm {ord}(1)$ corrections to the $\textrm {ord}( P{\kern-1pt}e)$ leading-order $g$ values are accurately reproduced for $ P{\kern-1pt}e \lesssim 10^{4}$ . In the subsequent comparison, we did not need to go beyond $ P{\kern-1pt}e=100$ .

Using the numerically generated solution, $V( P{\kern-1pt}e)$ is calculated using (3.8). The results are plotted in figure 3. They are in excellent agreement with the large- $ P{\kern-1pt}e$ asymptotic approximation (7.15).

Figure 3. Rescaled viscosity increment as a function of $ P{\kern-1pt}e$ . Solid: numerical calculation. Dashed: two-term approximation (7.15). Dotted: crude approximation, obtained by overlooking the $\textrm {ord}(1)$ constant in (7.14).

10. Concluding remarks

In classical macrorheology, there is an interest in understanding the nonlinear dependence upon the shear magnitude, and in particular the dependence of the effective viscosity upon the Péclet number under strong shear (Brady & Morris 1997; Morris & Katyal 2002). Our interest is in the analogous limit in the comparable microrheological problem. Using a detailed analysis of the microstructure at large Péclet numbers, we have derived a closed-form approximation for the variation of the viscosity increment with the Péclet number.

In our asymptotic analysis, we have employed throughout the Fraenkel interpretation of matched asymptotic expansions (Fraenkel 1969), where terms that differ by a logarithm of the expansion parameter are considered to be of the same asymptotic order. With that interpretation, our key result (7.15) constitutes a two-term approximation. With the asymptotic error being algebraically small, the excellent agreement between (7.15) and the viscosity increment predicted by the numerical scheme, as depicted in figure 3, is hardly surprising.

Beyond the direct contribution to microrheology, the present analysis may be of fundamental interest in boundary-layer theory – in particular, the structure of boundary layers about spherical surfaces. For ‘conventional’ hydrodynamic advection, associated with an imposed uniform flow, such layers typically break down at the rear stagnation point, morphing into a narrow wake. This is the case for the problem of solute transport at large Péclet numbers (Levich 1962) and vorticity transport about bubbles at large Reynolds numbers (Moore 1963; Harper 1972). The present problem is somewhat non-conventional, as the boundary layer breaks down at the equator, resulting in a broad wake. While there are known examples in the literature of equatorial breakdown (Stewartson 1958; Schnitzer, Frankel & Yariv 2013), they are associated with certain fore–aft symmetries that are absent in uniform-flow advection. In the present problem, the ‘early’ breakdown may be traced back to the neglect of hydrodynamic interactions. With an advection term which is uniform throughout, the equatorial breakdown is a consequence of the blunt violation of the no-flux condition.

There are two natural follow-ups of the present contribution. The first involves the calculation of the velocity fluctuations – a tensor that was calculated by SB05 at leading order. The second is concerned with the calculation of the longitudinal and transverse force-induced diffusivities. In their analysis of that problem, Zia & Brady (2010) employed the physical model of SB05, as used here, where hydrodynamic interactions are neglected. Based upon the leading-order boundary-layer solution for a vectorial ‘displacement’ field, Zia & Brady (2010) obtained the leading-order diffusivities at small and large $ P{\kern-1pt}e$ . In the latter limit, both the longitudinal and transverse diffusivities are proportional to $ P{\kern-1pt}e$ . An analysis in the spirit of the present one could provide the (presumably $\textrm {ord}(1)$ ) corrections in that asymptotic limit.

Other future directions, which have to do with the interactions between the probe and the bath particles, appear more challenging. The first involves extending the present work to incorporate hydrodynamic interactions (Khair & Brady 2006). The second has to do with the validity of the underlying hard-sphere approximation. That approximation implicitly assumes that the range $d$ of colloidal interaction is small compared with the size of the exclusion sphere, $d\ll \ell$ , whereby the interaction only appears as a no-flux boundary condition (preventing colloids from effectively passing through each other). In large- $ P{\kern-1pt}e$ investigations (both in SB05 and the present contribution) the restriction is stricter, as $d$ must be small compared with the boundary-layer thickness $\ell / P{\kern-1pt}e$ . If this is not the case, the boundary-layer analysis must be revised.