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 Reference Hinch1991) 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 (Reference Khair and Brady2006) 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}

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 Reference Hinch1991) 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).