2.1. Problem formulation

A rigid sphere (radius $a$ ) is held fixed in an unbounded solution (viscosity $\mu$ ) which is subject to an otherwise homogeneous concentration gradient of magnitude $G$ in the direction of gravity. The sphere is impermeable to the solute, resulting in modification to the gradient. We assume a linear dependence of the mass density $\rho '$ upon the solute concentration $c'$ , denoting the proportionality coefficient by $\alpha$ . (Hereafter, the prime indicates dimensional field variables.) Without loss of generality, we write

(2.1) \begin{equation} \rho ' - \rho '_0=\alpha (c'-c'_0), \end{equation}

where $c'_0$ and $\rho '_0$ denote the ambient concentration and density values at the sphere centre. Measuring the position vector $\boldsymbol{x}'$ from that centre, the ambient concentration profile is

(2.2) \begin{equation} c'_\infty (x') = c'_0 - G x', \end{equation}

wherein $x' = \hat {\boldsymbol{{i}}} \boldsymbol{\cdot } \boldsymbol{x}'$ in which $\hat {\boldsymbol{{i}}}$ is a unit vector anti-parallel to gravity. The resulting ambient stratification is obtained from (2.1) as $\rho '_\infty (x')=\rho '_0-\alpha Gx'$ , while the associated hydrostatic profile is

(2.3) \begin{equation} p'_\infty (x') = - \rho '_0 g x'+ \frac {1}{2}\alpha g G {x'}^2, \end{equation}

wherein $g$ is the acceleration due to gravity.

Our interest here is in the steady state. The concentration $c'$ is then governed by: (i) the advection–diffusion equation,

(2.4) \begin{equation} \boldsymbol{u}'\boldsymbol{\cdot }{\nabla }' c' = \kappa {{\nabla }'}^2c', \end{equation}

where $\kappa$ is the molecular diffusivity of the solute and ${\nabla }'$ denotes the gradient operator with respect to $\boldsymbol{x}'$ ; (ii) the no-flux condition at the sphere boundary, representing its impermeability; and (iii) the imposed-gradient condition,

(2.5) \begin{equation} c' \sim c'_\infty (x') \quad \text{for} \quad |\boldsymbol{x}'|\gg a. \end{equation}

The fluid velocity $\boldsymbol{u}'$ and pressure $p'$ are governed by: (i) the continuity equation,

(2.6) \begin{equation}{\nabla }' \boldsymbol{\cdot }\boldsymbol{u}' = 0, \end{equation}

which has been tacitly employed in (2.4); (ii) the Stokes equation,

(2.7) \begin{equation} {\nabla }' \! p' = \mu {{\nabla }'}^2\boldsymbol{u}' - \hat {\boldsymbol{{i}}} g[\rho '_0+\alpha (c'-c'_0)], \end{equation}

wherein (2.1) has been employed; (iii) no slip at the sphere boundary; and (iv) velocity decay at large distances. Note the use of the Boussinesq approximation in (2.6)–(2.7).

Towards a dimensionless formulation, we normalise all length variables by $a$ ; thus, the dimensionless position vector and gradient operator are $\boldsymbol{x}=\boldsymbol{x}'/a$ and $\boldsymbol{\nabla } = a{\nabla }'$ . We employ spherical polar coordinates $(r,\theta ,\phi )$ with $r=0$ at the sphere centre and the ray $\theta =0$ pointing in the direction of $\hat {\boldsymbol{{i}}}$ . We use the scale $aG$ for normalising the excess concentration, relative to the ambient profile (2.2); we therefore define

(2.8) \begin{equation} c = \frac {c' -c'_\infty (x')}{aG}. \end{equation}

The body force associated with excess density (relative to the ambient density $\rho '_\infty (x')$ ) provides the scale $\alpha G g a^2$ for the modified pressure (relative to hydrostatic profile (2.3)); we therefore define

(2.9) \begin{equation} p =\frac {p' -p'_\infty (x')}{\alpha G g a^2}. \end{equation}

Balancing (excess) body forces with viscous stresses yields the velocity scale

(2.10) \begin{equation} \mathcal{U}=\frac {\alpha G g a^3}{\mu }. \end{equation}

We accordingly define the dimensionless velocity as $\boldsymbol{u} = \boldsymbol{u}'/\mathcal{U}$ .

The dimensionless problem governing the excess concentration $c$ consists of: (i) the advection–diffusion equation,

(2.11) \begin{equation} {\textit {Ra}}\,\boldsymbol{u}\boldsymbol{\cdot }(\boldsymbol{\nabla } c - \hat {\boldsymbol{{i}}})= {\nabla} ^2 c, \end{equation}

where ${\textit {Ra}}=a\,\mathcal{U}/\kappa$ is the Rayleigh number; (ii) the no-flux condition,

(2.12) \begin{equation} \frac {\partial c}{\partial r}=\cos \theta \quad \text{at} \quad r=1; \end{equation}

and (iii) the decay condition,

(2.13) \begin{equation} \lim _{r\to \infty }c= 0. \end{equation}

The dimensionless flow problem consists of: (i) the continuity equation,

(2.14) \begin{equation} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}=0; \end{equation}

(ii) the Stokes equation, subject to the Boussinesq body force,

(2.15) \begin{equation} \boldsymbol{\nabla }\! p = {\nabla} ^2\boldsymbol{u} - \hat {\boldsymbol{{i}}} c; \end{equation}

(iii) the no-slip condition,

(2.16) \begin{equation} \boldsymbol{u}=\boldsymbol{0} \quad \text{at} \quad r=1; \end{equation}

and (iv) velocity decay,

(2.17) \begin{equation} \lim _{r\to \infty }\boldsymbol{u}=\boldsymbol{0}. \end{equation}

The above system uniquely defines $c$ and $\boldsymbol{u}$ for all $\textit {Ra}$ , the single governing parameter; it defines the pressure up to additive constant, which is physically meaningless. Note that the flow is coupled to the concentration by the body force in (2.15), while the concentration is (nonlinearly) coupled to the flow by the advective term in (2.11).

Despite the nonlinearity, the problem possesses two key symmetries. The first is an obvious axial symmetry, which implies that the velocity possesses the form $\boldsymbol{u}=\hat {\boldsymbol{e}}_r u+\hat {\boldsymbol{e}}_\theta v$ (there is no swirl), where both $u$ and $v$ – as well as $c$ and $p$ – are independent of the azimuthal angle $\phi$ . The second symmetry is about the equatorial plane $x=0$ . It is readily verified that $c$ and $v$ are odd functions of $x$ , while $u$ and $p$ are even functions of $x$ . This symmetry implies that the flow cannot impart any hydrodynamic force on the sphere in the vertical direction, suggesting that the preceding analysis equally holds for a freely suspended particle. This aspect is addressed in § 4. Note that the second symmetry ceases to hold in the presence of inertia.

2.2. Paradox and resolution

Our interest is in the limit where $\textit {Ra}$ is small. Since the nonlinear coupling between $c$ and $\boldsymbol{u}$ is then eliminated in (2.11), this limit may appear to be straightforward. Indeed, setting $\textit {Ra}$ to zero we readily find the solution of (2.11)–(2.13) as

(2.18) \begin{equation} c = -\frac {\cos \theta }{2r^2}, \end{equation}

representing a dipole anti-parallel to $\hat {\boldsymbol{{i}}}$ . The paradox lies in the flow engendered by the body force associated with (2.18). Since that force decays slowly, as $r^{-2}$ , it follows from (2.15) that the velocity field attains an (orientation-dependent) $\textrm {ord}(1)$ velocity field at large $r$ , in contradiction to decay requirement (2.17).

Similarly to classical paradoxes in forced convection (Acrivos & Taylor Reference Acrivos and Taylor1962), the present free convection conundrum may be traced to the breakdown of negligible-advection approximation at large distances. Indeed, with $\boldsymbol{u}$ being $\textrm {ord}(1)$ for large $r$ , it follows that the left-hand side of (2.11) is $\textrm {ord}({\textit {Ra}})$ there; given (2.18), the (presumably dominant) right-hand side of (2.11) behaves as $r^{-4}$ . No matter how small $\textit {Ra}$ is, both sides become comparable at remote distances of order ${\textit {Ra}}^{-1/4}$ . The purely diffusive approximation obtained from the neglect of advection in (2.11) is non-uniformly valid.

We resolve the present impasse using matched asymptotic expansions in the spirit of Proudman & Pearson (Reference Proudman and Pearson1957). Thus, the fluid domain is conceptually decomposed into two regions: an ‘inner’ region on the scale of the particle, and a remote ‘outer’ region at distances of order ${\textit {Ra}}^{-1/4}$ . We separately analyse the coupled flow–transport in each of these regions. For our purposes it suffices to consider only the leading-order problem in each region.

The result (2.18), which happens to satisfy the decay condition (2.13), should be interpreted as a leading-order inner approximation. When solving for the flow problem in the inner region we disregard the decay condition (2.17), which applies to the outer region. We therefore allow the leading-order flow to attain $O(1)$ values as $r\to \infty$ . As will become evident, this will facilitate asymptotic matching with an outer approximation.

2.3. Inner flow problem

Forming the curl of (2.15) yields the equation

(2.19) \begin{equation}{\nabla} ^2 \boldsymbol \omega = \hat {\boldsymbol{{i}}}\times \boldsymbol{\nabla } c, \end{equation}

governing the vorticity $\boldsymbol \omega =\boldsymbol{\nabla }\times \boldsymbol{u}$ . Substituting (2.18) gives, at leading order,

(2.20) \begin{equation} {\nabla} ^2 \boldsymbol \omega = -\frac {3\sin \theta \cos \theta }{4r^3}\hat {\boldsymbol{e}}_\phi . \end{equation}

We introduce the Stokes stream function $\psi$ , defined by $\boldsymbol{u}=\boldsymbol{\nabla }\phi \times \boldsymbol{\nabla }\psi$ (Happel & Brenner Reference Happel and Brenner1965):

(2.21) \begin{equation} u = - \frac {1}{r^2\sin \theta }\frac {\partial \psi }{\partial \theta }, \quad v = \frac {1}{r\sin \theta }\frac {\partial \psi }{\partial r}. \end{equation}

With that definition we have

(2.22) \begin{equation} \boldsymbol \omega = \frac {\hat {\boldsymbol{e}}_\phi }{r\sin \theta }\mathrm{E}^2\psi ,\end{equation}

wherein

(2.23) \begin{equation} \mathrm{E}^2 = \frac {\partial ^2}{\partial r^2} + \frac {\sin \theta }{r^2}\frac {\partial }{\partial \theta }\left (\frac {1}{\sin \theta }\frac {\partial }{\partial \theta }\right )\!. \end{equation}

Substitution of (2.22) into (2.20) yields

(2.24) \begin{equation} \mathrm{E}^4\psi = -\frac {3\zeta (1-\zeta ^2)}{2r^2} ,\end{equation}

where $\zeta =\cos \theta$ . Equation (2.24) constitutes the inhomogeneous counterpart of the familiar equation governing $\psi$ in axisymmetric Stokes flow. It is supplemented by the boundary conditions

(2.25) \begin{equation} \psi =\frac {\partial \psi }{\partial r}=0 \quad \text{at} \quad r=1, \end{equation}

obtained from (2.16).

Figure 1.

Streamlines $\psi =\text{const.}$ of the quadrupolar flow on the particle scale, as obtained from (2.26).

A particular integral of (2.24) is $-r^2\zeta (1-\zeta ^2)/16$ . (Note that it results in a non-decaying velocity, as already anticipated.) The complementary solutions to the homogeneous counterpart of (2.24) are well known (Happel & Brenner Reference Happel and Brenner1965). In seeking solutions proportional to $\zeta (1-\zeta ^2)$ we disregard modes that diverge faster than $r^2$ at large $r$ , as they result in an even more singular behaviour. This leads to $\psi =-\zeta (1-\zeta ^2) (r^2+A+Br^{-2})/16$ . (Conditions (2.25) and the requirement that $\psi =O(r^2)$ at large $r$ exclude the superposition of complementary solutions of different angular dependence.) The constants $A$ and $B$ are determined from (2.25) as $-2$ and $1$ , respectively. We conclude that

(2.26) \begin{equation} \psi = \frac {2-r^2-r^{-2}}{16}\sin ^2\theta \cos \theta , \end{equation}

a solution that has been reported in the supplementary information of Camassa et al. (Reference Camassa, Harris, Hunt, Kilic and McLaughlin2019). The associated streamlines are depicted in figure 1. Near the boundary we find

(2.27) \begin{equation} v \sim - \frac {r-1}{2}\sin \theta \cos \theta , \end{equation}

implying velocity towards the poles ( $\theta =0,\pi$ ). This agrees with the intuitive expectation of upward creep along a sloping bottom wall (Phillips Reference Phillips1970).

At large distance, $\psi$ is dominated by the particular integral, with the associated velocity components satisfying

(2.28) \begin{equation} \lim _{r\to \infty }u = \frac {2-3\sin ^2\theta }{16}, \quad \lim _{r\to \infty }v = -\frac {\sin \theta \cos \theta }{8}. \end{equation}

In particular, the velocity approaches an inward radial velocity at the equatorial plane,

(2.29) \begin{equation} \lim _{r\to \infty } \boldsymbol{u}= -\frac {1}{16}\hat {\boldsymbol{e}}_r \quad \text{at} \quad \theta =\pi /2, \end{equation}

and an outward radial velocity on the symmetry axis,

(2.30) \begin{equation} \lim _{r\to \infty } \boldsymbol{u}=\frac {1}{8}\hat {\boldsymbol{e}}_r \quad \text{at} \quad \theta =0,\pi . \end{equation}

2.4. Outer region

Formally speaking, the inner analysis on the scale of the particle corresponds to the limit process ${\textit {Ra}}\to 0$ with $\boldsymbol{x}$ fixed. We now consider the outer limit process:

(2.31) \begin{equation} {\textit {Ra}}\to 0 \quad \text{with} \quad \boldsymbol{X}={\textit {Ra}}^{1/4}\boldsymbol{x} \quad \text{fixed}. \end{equation}

Given the dipole decay of (2.18) and the $\textrm {ord}(1)$ limit attained by $\boldsymbol{u}$ at large $r$ , we anticipate that the concentration is of order ${\textit {Ra}}^{1/2}$ at the outer region, while the velocity is of order unity there. We therefore employ the rescaled outer variables:

(2.32a,b,c) \begin{equation} C = {\textit {Ra}}^{-1/2} c, \quad \boldsymbol{U}=\boldsymbol{u}, \quad P = {\textit {Ra}}^{-1/4} p, \end{equation}

where the pressure scaling follows from the Stokes equation (2.15). We also employ the stretched radial coordinate and gradient operator:

(2.33a,b) \begin{equation} R = {\textit {Ra}}^{1/4} r, \quad \boldsymbol{\nabla }_{\boldsymbol{X}} = {\textit {Ra}}^{-1/4}\boldsymbol{\nabla }. \end{equation}

(Note that $R=|\boldsymbol{X}|$ .)

We are only concerned with the leading-order outer approximation. The concentration is therefore governed by the simplified advection–diffusion equation (cf. (2.11))

(2.34) \begin{equation} \hat {\boldsymbol{{i}}}\boldsymbol{\cdot }\boldsymbol{U} +{\nabla} _{\boldsymbol{X}}^2C=0, \end{equation}

where the first term represents advection of the ambient concentration profile. It also satisfies the local behaviour,

(2.35) \begin{equation} C\sim \frac {1}{2}\hat {\boldsymbol{{i}}}\boldsymbol{\cdot }\boldsymbol{\nabla }_{\boldsymbol{X}} \frac {1}{R}\quad \text{as} \quad R\to 0, \end{equation}

which follows from asymptotic matching with (2.18). We can combine (2.34) and the singular behaviour (2.35) to form the inhomogeneous differential equation

(2.36) \begin{equation} \hat {\boldsymbol{{i}}}\boldsymbol{\cdot }\boldsymbol{U} + {\nabla} _{\boldsymbol{X}}^2C=-2\pi \hat {\boldsymbol{{i}}}\boldsymbol{\cdot }\boldsymbol{\nabla }_{\boldsymbol{X}}\delta (\boldsymbol{X}), \end{equation}

where $\delta$ denotes the 3-D Dirac delta function. (Note that $\boldsymbol{U}$ is regular at the origin.) Equation (2.36) is supplemented by the decay condition,

(2.37) \begin{equation} \lim _{R\to \infty }C=0, \end{equation}

which follows from (2.13). The flow is governed by the continuity and Stokes equations,

(2.38a,b) \begin{equation} \boldsymbol{\nabla }_{\boldsymbol{X}}\boldsymbol{\cdot }\boldsymbol{U}=0, \quad \boldsymbol{\nabla }_{\boldsymbol{X}} P= {\nabla} _{\boldsymbol{X}}^2\boldsymbol{U} - \hat {\boldsymbol{{i}}} C, \end{equation}

together with the decay condition,

(2.39) \begin{equation} \lim _{R\to \infty }\boldsymbol{U}=\boldsymbol{0}. \end{equation}

We employ the 3-D Fourier transform, defined generically as

(2.40) \begin{equation} \mathcal F_3 \{G\}=\int _{\mathbb{R}^3} G(\boldsymbol{X}) \mathrm{e}^{-{\rm i}\boldsymbol{q}\boldsymbol{\cdot }\boldsymbol{X}}\, {\rm d}^3\boldsymbol{X}, \end{equation}

and written $\tilde {G}(\boldsymbol{q})$ . The inverse transform is

(2.41) \begin{equation} \mathcal F_3^{-1} \{\tilde G\} = \frac {1}{(2\pi )^3} \int _{\mathbb{R}^3} \tilde {G}(\boldsymbol{q}) \mathrm{e}^{{\rm i}\boldsymbol{q}\boldsymbol{\cdot }\boldsymbol{X}}\, {\rm d}^3\boldsymbol{q}. \end{equation}

Forming the transform of (2.36) while making use of the decay condition (2.37) yields

(2.42) \begin{equation} \hat {\boldsymbol{{i}}}\boldsymbol{\cdot }\tilde {\boldsymbol{U}} - q^2\tilde C = -2\pi {\mathrm{i}}\,\hat {\boldsymbol{{i}}}\boldsymbol{\cdot }\boldsymbol{q}. \end{equation}

Forming the transform of (2.38) while making use of the decay condition (2.39) yields

(2.43a,b) \begin{equation} \boldsymbol{q}\boldsymbol{\cdot }\tilde {\boldsymbol{U}}=0, \quad \mathrm{i}\boldsymbol{q}\tilde P= -q^2\tilde {\boldsymbol{U}} - \hat {\boldsymbol{{i}}} \tilde C. \end{equation}

Combining (2.43) to eliminate $\tilde P$ we obtain

(2.44) \begin{equation} \hat {\boldsymbol{{i}}}\boldsymbol{\cdot }(\unicode{x1D644}-\hat {\boldsymbol{q}}\hat {\boldsymbol{q}})\tilde C + q^2 \tilde {\boldsymbol{U}}=\boldsymbol{0}, \end{equation}

where $\hat {\boldsymbol{q}}=\boldsymbol{q}/q$ . Substitution into (2.42) yields

(2.45) \begin{equation} \tilde C = 2\pi {\mathrm{i}} \frac {\hat {\boldsymbol{{i}}}\boldsymbol{\cdot }\hat {\boldsymbol{q}} \,q^3}{1-(\hat {\boldsymbol{{i}}}\boldsymbol{\cdot }\hat {\boldsymbol{q}})^2+q^4}. \end{equation}

Finally, we obtain from (2.44) that

(2.46) \begin{equation} \tilde {\boldsymbol{U}} = -2\pi {\mathrm{i}} \, \frac {\hat {\boldsymbol{{i}}}\boldsymbol{\cdot }(\unicode{x1D644}-\hat {\boldsymbol{q}}\hat {\boldsymbol{q}})\hat {\boldsymbol{{i}}}\boldsymbol{\cdot }\boldsymbol{q}}{1-(\hat {\boldsymbol{{i}}}\boldsymbol{\cdot }\hat {\boldsymbol{q}})^2+q^4}. \end{equation}

Forming the inverse transforms of (2.45) and (2.46) yields, respectively,

(2.47) \begin{gather} C = \frac {2\pi \mathrm{i}}{(2\pi )^3} \int _{\mathbb{R}^3} \frac {\hat {\boldsymbol{{i}}}\boldsymbol{\cdot }\hat {\boldsymbol{q}} \,q^3}{1-(\hat {\boldsymbol{{i}}}\boldsymbol{\cdot }\hat {\boldsymbol{q}})^2+q^4} \mathrm{e}^{{\rm i}\boldsymbol{q}\boldsymbol{\cdot }\boldsymbol{X}}\,{\rm d}^3\boldsymbol{q}, \end{gather}

(2.48) \begin{gather} \boldsymbol{U} = -\frac {2\pi \mathrm{i}}{(2\pi )^3} \int _{\mathbb{R}^3} \frac {\hat {\boldsymbol{{i}}}\boldsymbol{\cdot }(\unicode{x1D644}-\hat {\boldsymbol{q}}\hat {\boldsymbol{q}})\hat {\boldsymbol{{i}}}\boldsymbol{\cdot }\hat {\boldsymbol{q}}\,q}{1-(\hat {\boldsymbol{{i}}}\boldsymbol{\cdot }\hat {\boldsymbol{q}})^2+q^4} \mathrm{e}^{{\rm i}\boldsymbol{q}\boldsymbol{\cdot }\boldsymbol{X}}\,{\rm d}^3\boldsymbol{q}. \end{gather}

2.5. Small- $R$ behaviour of the outer fields

When considering the concentration fields in the inner and outer regions, matching has already been enforced by imposing the local behaviour (2.35). To verify that the outer concentration (2.47) indeed satisfies (2.35) we note that the integral in (2.47) is dominated by the large- $q$ contribution at small $R$ . It follows that

(2.49) \begin{equation} C \sim 2\pi \mathrm{i} \hat {\boldsymbol{{i}}}\boldsymbol{\cdot } \mathcal{F}_3^{-1}\left \{ \frac {\boldsymbol{q}}{q^2} \right \} \quad \text{as}\quad R\to 0. \end{equation}

Recalling the generalised Fourier transform

(2.50) \begin{equation} \mathcal{F}_3\left \{\frac {1}{4\pi R}\right \}= \frac {1}{q^2}, \end{equation}

we find that

(2.51) \begin{equation} \mathcal{F}_3^{-1}\left \{ \frac {\boldsymbol{q}}{q^2} \right \} = \frac {1}{4\pi \mathrm{i}} \boldsymbol{\nabla }_{\boldsymbol{X}}\frac {1}{R}. \end{equation}

We have therefore reproduced (2.35).

Consider now the velocity field. With (2.46) dominated by large $q$ we find that

(2.52) \begin{equation} \boldsymbol{U} \sim -2\pi \mathrm{i} \hat {\boldsymbol{{i}}}\boldsymbol{\cdot } \mathcal{F}_3^{-1}\left \{ \frac {(\unicode{x1D644}-\hat {\boldsymbol{q}}\hat {\boldsymbol{q}})\boldsymbol{q}}{q^4} \right \} \boldsymbol{\cdot } \hat {\boldsymbol{{i}}} \quad \text{as}\quad R\to 0. \end{equation}

Recalling the generalised Fourier transforms

(2.53) \begin{equation} \mathcal{F}_3\left \{\frac {R}{8\pi }\right \}= - \frac {1}{q^4}, \quad \mathcal{F}_3\left \{\frac {R^3}{96\pi }\right \}= \frac {1}{q^6}, \end{equation}

we readily obtain

(2.54) \begin{equation} \mathcal{F}_3^{-1}\left \{ \frac {(\unicode{x1D644}-\hat {\boldsymbol{q}}\hat {\boldsymbol{q}})\boldsymbol{q}}{q^4} \right \} = -\frac {1}{8\pi \mathrm{i}} \unicode{x1D644}\boldsymbol{\nabla }_{\boldsymbol{X}}R +\frac {1}{96\pi \mathrm{i}} \boldsymbol{\nabla }_{\boldsymbol{X}}\boldsymbol{\nabla }_{\boldsymbol{X}}\boldsymbol{\nabla }_{\boldsymbol{X}}R^3. \end{equation}

Making use of the derivatives

(2.55) \begin{equation} \boldsymbol{\nabla }_{\boldsymbol{X}}R = \hat {\boldsymbol{e}}_r, \quad \boldsymbol{\nabla }_{\boldsymbol{X}}\boldsymbol{\nabla }_{\boldsymbol{X}}\boldsymbol{\nabla }_{\boldsymbol{X}}R^3 = 3 \big[ \unicode{x1D644} \hat {\boldsymbol{e}}_r + (\unicode{x1D644}\hat {\boldsymbol{e}}_r)^\dagger + \hat {\boldsymbol{e}}_r\unicode{x1D644} - \hat {\boldsymbol{e}}_r\hat {\boldsymbol{e}}_r\hat {\boldsymbol{e}}_r \big ], \end{equation}

wherein $\dagger$ denotes tensor transposition, we obtain from (2.52) that

(2.56) \begin{equation} \boldsymbol{U} \sim \frac {1}{16}\big[ 2\hat {\boldsymbol{{i}}}\hat {\boldsymbol{{i}}}\boldsymbol{\cdot }\hat {\boldsymbol{e}}_r -\hat {\boldsymbol{e}}_r + (\hat {\boldsymbol{{i}}}\boldsymbol{\cdot }\hat {\boldsymbol{e}}_r)^2\hat {\boldsymbol{e}}_r\big ]\quad \text{as}\quad R\to 0. \end{equation}

Recalling that $\hat {\boldsymbol{{i}}}\boldsymbol{\cdot }\hat {\boldsymbol{e}}_r=\cos \theta$ therefore yields

(2.57) \begin{equation} \boldsymbol{U} \sim \frac {1}{16}\big [ \hat {\boldsymbol{e}}_r(2-3\sin ^2\theta ) - 2\hat {\boldsymbol{e}}_\theta \sin \theta \cos \theta \big ]\quad \text{as}\quad R\to 0, \end{equation}

which trivially matches (2.28).

2.6. Simplified form of the outer fields

With the goal of addressing the large- $R$ limit, we wish to derive a simplified form for the outer fields. To that end, we restrict our attention to the symmetry axis $\theta =0$ , where $\boldsymbol{X}= \hat {\boldsymbol{{i}}} R$ .

We begin with the velocity field $\boldsymbol{U}$ . By axial symmetry, it is purely radial on the symmetry axis, $\boldsymbol{U}=\hat {\boldsymbol{{i}}}U$ . We obtain from (2.48) that

(2.58) \begin{equation} U(\theta =0) = -\frac {\mathrm{i}}{(2\pi )^2} \int _{\mathbb{R}^3} \frac {[1-(\hat {\boldsymbol{{i}}}\boldsymbol{\cdot }\hat {\boldsymbol{q}})^2]\hat {\boldsymbol{{i}}}\boldsymbol{\cdot }\hat {\boldsymbol{q}}\,q}{1-(\hat {\boldsymbol{{i}}}\boldsymbol{\cdot }\hat {\boldsymbol{q}})^2+q^4} \exp \big (\mathrm{i} qR \, \hat {\boldsymbol{{i}}}\boldsymbol{\cdot }\hat {\boldsymbol{q}}\big )\,{\rm d}^3\boldsymbol{q}. \end{equation}

We parametrise $\boldsymbol{q}$ using another set of polar coordinates, say $(q,\vartheta ,\varphi )$ , with the zenith angle $\vartheta$ measured from $\hat {\boldsymbol{{i}}}$ . With that parametrisation, $\hat {\boldsymbol{{i}}}\boldsymbol{\cdot }\hat {\boldsymbol{q}}=\cos \vartheta$ . Substitution into (2.58) and recalling that $U$ is real-valued gives

(2.59) \begin{equation} U(\theta =0) = \frac {1}{2\pi } \iint _{(0,\pi )\times (0,\infty )} \frac {q^3\sin ^3\vartheta \cos \vartheta }{\sin ^2\vartheta +q^4} \sin (q R \cos \vartheta ) \, \mathrm{d}\vartheta \, \mathrm{d} q. \end{equation}

Since the integrand is an even function of $q$ and is invariant under the transformation $\vartheta \mapsto \pi -\vartheta$ we may modify the integration domain to the infinite strip $(0,\pi /2)\times (-\infty ,\infty )$ . Writing the double integral as an iterated integral we therefore obtain

(2.60) \begin{equation} U(\theta =0) = \frac {1}{2\pi } \int _0^{\pi /2} I(R,\vartheta ) \sin ^3\vartheta \cos \vartheta \, \mathrm{d}\vartheta , \end{equation}

wherein

(2.61) \begin{equation} I(R,\vartheta ) = \int _{-\infty }^\infty \frac {q^3 \sin (q R \cos \vartheta ) }{\sin ^2\vartheta +q^4} \, \mathrm{d} q. \end{equation}

Using again the even dependence upon $q$ we replace $\sin (q R \cos \vartheta )$ with $\mathrm{e}^{{\rm i} q R \cos \vartheta }/\mathrm{i}$ , transforming $I$ to a one-dimensional exponential Fourier transform. Making use of the change-of-variables $q=\tau \sin ^{1/2}\vartheta$ therefore gives

(2.62) \begin{equation} I(R,\vartheta ) = \frac {1}{\mathrm{i}}\int _{-\infty }^\infty \frac {\tau ^3 \mathrm{e}^{{\rm i}\tau R \sin ^{1/2}\vartheta \cos \vartheta } }{1+\tau ^4} \, \mathrm{d} \tau . \end{equation}

Noting that $\cos \vartheta$ is non-negative in the integration domain of (2.60), we evaluate the integral using standard residue-theorem techniques to obtain

(2.63) \begin{equation} I(R,\vartheta )=\pi \exp \left (-\frac {R\sin ^{1/2}\vartheta \cos \vartheta }{\sqrt {2}}\right ) \cos \left (\frac {R\sin ^{1/2}\vartheta \cos \vartheta }{\sqrt {2}}\right )\!. \end{equation}

The velocity may therefore be written as

(2.64) \begin{equation} U(\theta =0) = \frac {1}{2} \operatorname {Re}\int _0^{\pi /2} \exp \left \{-\frac {R(1-\mathrm{i})\sin ^{1/2}\vartheta \cos \vartheta }{\sqrt {2}}\right \}\sin ^3\vartheta \cos \vartheta \, \mathrm{d}\vartheta . \end{equation}

Consider now the excess concentration $C$ . We obtain from (2.47) that

(2.65) \begin{equation} C(\theta =0) = \frac {\mathrm{i}}{(2\pi )^2} \int _{\mathbb{R}^3} \frac {\hat {\boldsymbol{{i}}}\boldsymbol{\cdot }\hat {\boldsymbol{q}}\,q^3}{1-(\hat {\boldsymbol{{i}}}\boldsymbol{\cdot }\hat {\boldsymbol{q}})^2+q^4} \exp \big(\mathrm{i} qR \, \hat {\boldsymbol{{i}}}\boldsymbol{\cdot }\hat {\boldsymbol{q}}\big )\,{\rm d}^3\boldsymbol{q}. \end{equation}

Making use of the same parametrisation as above and recalling that $C$ is real-valued gives

(2.66) \begin{equation} C(\theta =0) = -\frac {1}{2\pi } \iint _{(0,\pi )\times (0,\infty )} \frac {q^5\sin \vartheta \cos \vartheta }{\sin ^2\vartheta +q^4} \sin (q R \cos \vartheta ) \, \mathrm{d}\vartheta \, \mathrm{d} q, \end{equation}

or, upon changing the integration domain,

(2.67) \begin{equation} C(\theta =0) = -\frac {1}{2\pi } \iint _{(0,\pi /2)\times (-\infty ,\infty )} \frac {q^5\sin \vartheta \cos \vartheta }{\sin ^2\vartheta +q^4} \sin (q R \cos \vartheta ) \, \mathrm{d}\vartheta \, \mathrm{d} q. \end{equation}

We write the above double integral as an iterated integral:

(2.68) \begin{equation} C(\theta =0) = -\frac {1}{2\pi } \int _0^{\pi /2} J(R,\vartheta ) \sin \vartheta \cos \vartheta \, \mathrm{d}\vartheta , \end{equation}

wherein

(2.69) \begin{equation} J(R,\vartheta ) = \int _{-\infty }^\infty \frac {q^5 \sin (q R \cos \vartheta ) }{\sin ^2\vartheta +q^4} \, \mathrm{d} q. \end{equation}

Since the above integral does not converge in the classical sense, it must be interpreted in the generalised sense. Using again the even dependence upon $q$ we replace $\sin (q R \cos \vartheta )$ with $\mathrm{e}^{{\rm i} q R \cos \vartheta }/\mathrm{i}$ , transforming $J$ to a one-dimensional generalised Fourier transform. Making use of the change-of-variables $q=\tau \sin ^{1/2}\vartheta$ therefore gives

(2.70) \begin{equation} J(R,\vartheta ) = \frac {\sin \vartheta }{\mathrm{i}}\int _{-\infty }^\infty \frac {\tau ^5 \mathrm{e}^{{\rm i}\tau \xi } }{1+\tau ^4} \, \mathrm{d} \tau , \end{equation}

where

(2.71) \begin{equation} \xi =R \sin ^{1/2}\vartheta \cos \vartheta . \end{equation}

Making use of the generalised transform $\int _{-\infty }^\infty \mathrm{e}^{{\rm i}\tau \xi } \, \mathrm{d} \tau =2\pi \delta (\xi )$ (Lighthill Reference Lighthill1958), it follows that

(2.72) \begin{equation} \int _{-\infty }^\infty \tau \mathrm{e}^{{\rm i}\tau \xi } \, \mathrm{d} \tau =-2\pi \mathrm{i} \frac {\mathrm{d} \delta }{\mathrm{d} \xi }. \end{equation}

Employing (2.72) together with standard residue-theorem techniques, the integral appearing in (2.70) is obtained as

(2.73) \begin{equation} -2\pi \mathrm{i} \frac {\mathrm{d} \delta }{\mathrm{d} \xi } - \mathrm{i}\pi \mathrm{e}^{-\xi /\sqrt {2}} \sin (\xi /\sqrt {2}). \end{equation}

We conclude that

(2.74) \begin{align} C(\theta =0) &= \int _0^{\pi /2} \left .\frac {\mathrm{d} \delta }{\mathrm{d} \xi }\right |_{\xi =R \sin ^{1/2}\vartheta \cos \vartheta }\sin ^2\vartheta \cos \vartheta \, \mathrm{d}\vartheta \nonumber \\[3pt] &\quad + \frac {1}{2} \operatorname {Im} \int _0^{\pi /2} \exp \left \{-\frac {R(1-\mathrm{i})\sin ^{1/2}\vartheta \cos \vartheta }{\sqrt {2}}\right \} \sin ^2\vartheta \cos \vartheta \, \mathrm{d}\vartheta , \end{align}

or, upon using the properties of the Dirac delta function (Lighthill Reference Lighthill1958),

(2.75) \begin{equation} C(\theta =0) = -\frac {1}{2R^2} + \frac {1}{2} \operatorname {Im} \int _0^{\pi /2} \exp \left \{-\frac {R(1-\mathrm{i})\sin ^{1/2}\vartheta \cos \vartheta }{\sqrt {2}}\right \} \sin ^2\vartheta \cos \vartheta \, \mathrm{d}\vartheta . \end{equation}

2.7. Far-field decay of the outer fields

For $R\to 0$ we find from (2.64) and (2.75) that

(2.76) \begin{equation} U(\theta =0) \sim \frac {1}{8}, \quad C(\theta =0) \sim -\frac {1}{2R^2}, \end{equation}

in agreement with (2.30) and (2.18).

More importantly, the simplified expressions (2.64) and (2.75) allow us to derive asymptotic expansions for both fields at large $R$ , where the integrals are evaluated using a Laplace-type approach. We first address the velocity (2.64). We note that the integral is dominated by two contributions: one from an $\textrm {ord}(R^{-2})$ neighbourhood of $\vartheta =0$ , and the other from an $\textrm {ord}(R^{-1})$ neighbourhood of $\vartheta =\pi /2$ . To obtain the former we make use of the stretched variable $\eta = R(\vartheta /2)^{1/2}$ and find the following leading contribution to the integral:

(2.77) \begin{equation} \frac {32}{R^8} \int _0^\infty \eta ^7 \mathrm{e}^{-(1-{\rm i})\eta } \, \mathrm{d}\eta = \frac {10{\,}080}{R^8}. \end{equation}

The contribution from the neighbourhood of $\vartheta =\pi /2$ requires more effort. Here, we define the stretched variable $\eta$ as $R(\pi -\vartheta )/2^{1/2}$ . Expanding the integrand of (2.64) for $R\to \infty$ (with $\eta$ fixed) up to order $R^{-8}$ and integrating gives the contribution

(2.78) \begin{equation} \frac {\mathrm{i}}{R^2} + \frac {0}{R^4}- \frac {45\mathrm{i}}{R^6} + \frac {5040}{R^8} + O(R^{-10}), \end{equation}

whose real part is $\textrm {ord}(R^{-8})$ – just as was the case with the neighbourhood of $\vartheta =0$ . Adding the contributions and substituting into (2.64) eventually yields

(2.79) \begin{equation} U(\theta =0) \sim \frac {7560}{R^8} \quad \text{as} \quad R\to \infty . \end{equation}

The concentration field is similarly handled. As with (2.64), the integral appearing in (2.75) is dominated by two contributions: one from an $\textrm {ord}(R^{-2})$ neighbourhood of $\vartheta =0$ , and the other from an $\textrm {ord}(R^{-1})$ neighbourhood of $\vartheta =\pi /2$ . The former is readily found as $-240\mathrm{i}/R^6$ , while the latter is

(2.80) \begin{equation} \frac {\mathrm{i}}{R^2} - \frac {3}{R^4}- \frac {120\mathrm{i}}{R^6} + O(R^{-8}). \end{equation}

Substitution into (2.75) reveals that the $\textrm {ord}(R^{-2})$ terms cancel out to yield the faster decay rate

(2.81) \begin{equation} C(\theta =0) \sim -\frac {180}{R^6} \quad \text{as} \quad R\to \infty . \end{equation}

3.1. Inner flow region

Substituting (3.1) into (2.19), while noting that $\boldsymbol \omega = \omega \hat {\boldsymbol{{k}}}$ in two dimensions, we obtain the leading-order vorticity equation (cf. (2.20)):

(3.2) \begin{equation} {\nabla} ^2 \omega = -\frac {\sin 2\theta }{r^2}. \end{equation}

We here employ the 2-D stream function $\psi$ , defined via $\boldsymbol{u}=\boldsymbol{\nabla }\times (\hat {\boldsymbol{{k}}}\psi )$ :

(3.3) \begin{equation} u = \frac {1}{r}\frac {\partial \psi }{\partial \theta }, \quad v = -\frac {\partial \psi }{\partial r}. \end{equation}

Upon substituting $\omega =-{\nabla} ^2\psi$ into (3.2) we obtain the inhomogeneous biharmonic equation (cf. (2.24)):

(3.4) \begin{equation} {\nabla} ^4\psi = \frac {\sin 2\theta }{r^2}, \end{equation}

which constitutes the present counterpart to the biharmonic equation in 2-D Stokes flows. The no-slip condition is still provided by (2.25).

A particular integral of (3.4) is $-({1}/{16}) r^2\ln r \sin 2\theta$ . It corresponds to a velocity that diverges as $r \ln r$ , slightly worse than originally anticipated. Seeking complementary solutions to the biharmonic equation that do not introduce a worse divergence, we have

(3.5) \begin{equation} \psi = F(r) \sin 2\theta , \end{equation}

where $F(r) = ( {1}/{16})( -r^2\ln r + \textit{Ar}^2+B + Dr^{-2})$ . To determine the integration constants all we have are the two conditions (2.25), which give here $A+B+D=0$ and $2A-2D=1$ . Expressing $B$ and $D$ in terms of $A$ therefore yields

(3.6) \begin{equation} F(r) = \frac {1}{16} \left [-r^2\ln r + \textit{Ar}^2 +1/2-2A + (A-1/2)r^{-2}\right ]. \end{equation}

The inner flow problem is therefore underdetermined. In contrast to the 3-D problem, asymptotic matching is required to uniquely determine the inner solution.

The leading-order inner approximations (3.1) and (3.5)–(3.6) are the same as those obtained by Alias & Page (Reference Alias and Page2017). Near the boundary we find that (cf. (2.27))

(3.7) \begin{equation} v \sim - \frac {4A-3}{8}(r-1) \sin \theta \cos \theta . \end{equation}

For $A\gt 3/4$ , this implies velocity towards the ‘poles’ ( $\theta =0,\pi$ ). Given the intuitive expectation of upward creep along a sloping bottom wall (Phillips Reference Phillips1970), we therefore expect $A\gt 3/4$ . Alias & Page (Reference Alias and Page2017) noted that in this case there exists $r_0\gt 1$ such that $F(r_0)=0$ , implying that $r=r_0$ is a circular streamline where $F$ changes sign, with $F(r)\gt 0$ for $1\lt r\lt r_0$ and $F(r)\lt 0$ for $r\gt r_0$ . (No comparable spherical stream surface exists in the 3-D case; see (2.26).) Moreover, they found from (3.6) that

(3.8) \begin{equation} r_0\sim \mathrm{e}^A \quad \text{as} \quad A\to \infty . \end{equation}

At large distances, $\psi \sim -( {1}/{16}) r^2 (\ln r - A)\sin 2\theta$ . Here, we follow the Fraenkel convention (Fraenkel Reference Fraenkel1969), where we do not separate asymptotic orders by logarithms of the expansion parameter; this allows for the application of Van Dyke’s matching rule (Van Dyke Reference Van Dyke1964). The associated velocity components therefore satisfy

(3.9a,b) \begin{equation} u \sim -\frac {r}{8} (\ln r - A)\cos 2\theta , \quad v \sim \frac {r}{8} (\ln r +1/2 - A)\sin 2\theta . \end{equation}

In particular, we find the radial velocity at the (vertical) $x$ axis ( $\theta =0,\pi$ ):

(3.10) \begin{equation} u(\theta =0) \sim -\frac {r}{8} (\ln r -A)\quad \text{as} \quad r\to \infty , \end{equation}

and the radial velocity at the (horizontal) $y$ axis ( $\theta =\pm \pi /2$ ):

(3.11) \begin{equation} u(\theta =\pi /2) \sim \frac {r}{8} (\ln r -A)\quad \text{as} \quad r\to \infty . \end{equation}

3.2. Outer analysis

With the outer region representing $O({\textit {Ra}}^{-1/4})$ distances, the outer limit is again defined by (2.31). Here, the inner behaviour (3.1) and (3.9) suggest that $c$ is $O({\textit {Ra}}^{1/4})$ and $\boldsymbol{u}=O({\textit {Ra}}^{-1/4})$ in the outer region. We therefore replace (2.32) with

(3.12a,b,c) \begin{equation} C = {\textit {Ra}}^{-1/4} c, \quad \boldsymbol{U}={\textit {Ra}}^{1/4}\boldsymbol{u}, \quad P = p, \end{equation}

and note that (2.33) remains valid. The leading-order solute transport in the outer region is still governed by (2.34), but the matching condition (2.35) is replaced by

(3.13) \begin{equation} C\sim -\hat {\boldsymbol{{i}}}\boldsymbol{\cdot }\boldsymbol{\nabla }_{\boldsymbol{X}} \ln {R}\quad \text{as} \quad R\to 0. \end{equation}

When using the delta-function formulation, we find that the combination of (2.34) and (3.13) results in (2.36), with $\delta$ now interpreted as a 2-D delta function. Equation (2.36) is supplemented by the decay condition (2.37). The outer flow is still governed by (2.38) and (2.39).

We now employ the Fourier transform in two dimensions, given generically as (cf. (2.40))

(3.14) \begin{equation} \mathcal F_2 \{G\}=\int _{\mathbb{R}^2} G(\boldsymbol{X}) \mathrm{e}^{-{\rm i}\boldsymbol{q}\boldsymbol{\cdot }\boldsymbol{X}}\,{\rm d}^2\boldsymbol{X}; \end{equation}

writing it as $\tilde {G}(\boldsymbol{q})$ , the inverse transform is (cf. (2.41))

(3.15) \begin{equation} \mathcal F_2^{-1} \{\tilde G\} = \frac {1}{(2\pi )^2} \int _{\mathbb{R}^2} \tilde {G}(\boldsymbol{q}) \mathrm{e}^{{\rm i}\boldsymbol{q}\boldsymbol{\cdot }\boldsymbol{X}}\, {\rm d}^2\boldsymbol{q}. \end{equation}

Since the mathematical form of the leading-order outer problem is the same in two and three dimensions, we again obtain the closed-form expressions (2.45)–(2.46) for the transforms $\tilde C$ and $\tilde {\boldsymbol{U}}$ . Given the $2\pi$ -factor difference in the form of the inverse transform, (2.47)–(2.48) are replaced by

(3.16) \begin{gather} C = \frac {2\pi \mathrm{i}}{(2\pi )^2} \int _{\mathbb{R}^2} \frac {\hat {\boldsymbol{{i}}}\boldsymbol{\cdot }\hat {\boldsymbol{q}} \,q^3}{1-(\hat {\boldsymbol{{i}}}\boldsymbol{\cdot }\hat {\boldsymbol{q}})^2+q^4} \mathrm{e}^{{\rm i}\boldsymbol{q}\boldsymbol{\cdot }\boldsymbol{X}}\,{\rm d}^2\boldsymbol{q}, \end{gather}

(3.17) \begin{gather} \boldsymbol{U} = -\frac {2\pi \mathrm{i}}{(2\pi )^2} \int _{\mathbb{R}^2} \frac {\hat {\boldsymbol{{i}}}\boldsymbol{\cdot }(\unicode{x1D644}-\hat {\boldsymbol{q}}\hat {\boldsymbol{q}})\hat {\boldsymbol{{i}}}\boldsymbol{\cdot }\hat {\boldsymbol{q}}\,q}{1-(\hat {\boldsymbol{{i}}}\boldsymbol{\cdot }\hat {\boldsymbol{q}})^2+q^4} \mathrm{e}^{{\rm i}\boldsymbol{q}\boldsymbol{\cdot }\boldsymbol{X}}\,{\rm d}^2\boldsymbol{q}. \end{gather}

3.3. Small- $R$ behaviour of the outer fields

Asymptotic matching has already been applied to the concentration field by imposing the local behaviour (3.13). To verify that the outer concentration (3.16) indeed satisfies that behaviour we note that the integral in (3.16) is dominated by large $q$ for $R\ll 1$ . We therefore have

(3.18) \begin{equation} C \sim 2\pi \mathrm{i}\hat {\boldsymbol{{i}}} \boldsymbol{\cdot } \mathcal{F}_2^{-1}\left \{\frac {\boldsymbol{q}}{q^2}\right \} \quad \text{as} \quad R\to 0. \end{equation}

Making use of the generalised transform

(3.19) \begin{equation} \mathcal F_2\left \{\frac {\ln R}{2\pi }\right \}=-\frac {1}{q^{2}}, \end{equation}

we find that $ \mathcal{F}_2^{-1}\{\boldsymbol{q}/q^2\} = -(2\pi \mathrm{i})^{-1} \boldsymbol{\nabla }_{\boldsymbol{X}}\ln R$ . We have therefore reproduced (3.13).

It may naively appear that a similar procedure may be applied to the fluid velocity. Assuming that (2.48) is dominated by large $q$ , that would require the inverse transforms of $\boldsymbol{q}/q^4$ and $\boldsymbol{q}\boldsymbol{q}\boldsymbol{q}/q^6$ . This suggests using the generalised transforms

(3.20) \begin{equation} \mathcal F_2\left \{\frac {R^2 \ln R}{8\pi }\right \}=\frac {1}{q^{4}}, \quad \mathcal F_2\left \{\frac {R^4 \ln R}{128\pi }\right \}=-\frac {1}{q^{6}}. \end{equation}

Here, however, one must recall the non-uniqueness of the generalised transforms. For example, the inverse transform of $q^{-4}$ can be taken as

(3.21) \begin{equation} \frac {R^2 \ln R}{8\pi } + \text{a multiple of }R^2, \end{equation}

since $R^2$ is a solution of the biharmonic equation. Recalling that the functions $R^2$ and $R^2\ln R$ must be considered together when performing Van Dyke matching, the above approach is futile. (It is easily verified that the non-uniqueness has no effect in the preceding 3-D analysis.). Without a small- $R$ approximation for the outer velocity, however, we cannot proceed to asymptotic matching with the inner velocity.

3.4. Simplified expressions for the outer velocity

Since our goal is to determine the unknown $A$ in (3.6), it suffices to consider the outer velocity at certain orientations. In what follows, we derive a reduced expression for that velocity at the vertical ray $\theta =0$ , where $\boldsymbol{X}=\hat {\boldsymbol{{i}}} R$ , and at the horizontal ray $\theta =\pi /2$ , where $\boldsymbol{X}=\hat {\boldsymbol{{j}}} R$ . In both cases it follows from the problem symmetry that the velocity is purely radial, $\boldsymbol{U}=\hat {\boldsymbol{e}}_r U$ . The problem symmetry also implies that $U(\theta =0) = U(\theta =\pi )$ and $U(\theta =\pi /2) = U(\theta =-\pi /2)$ .

At the vertical ray, where $\hat {\boldsymbol{e}}_r=\hat {\boldsymbol{{i}}}$ , we obtain from (3.17) that

(3.22) \begin{equation} U(\theta =0) = -\frac {\mathrm{i}}{2\pi } \int _{\mathbb{R}^2} \frac {[1-(\hat {\boldsymbol{{i}}}\boldsymbol{\cdot }\hat {\boldsymbol{q}})^2]\hat {\boldsymbol{{i}}}\boldsymbol{\cdot }\hat {\boldsymbol{q}}\,q}{1-(\hat {\boldsymbol{{i}}}\boldsymbol{\cdot }\hat {\boldsymbol{q}})^2+q^4} \exp \big(\mathrm{i} qR \, \hat {\boldsymbol{{i}}}\boldsymbol{\cdot }\hat {\boldsymbol{q}}\big )\,{\rm d}^2\boldsymbol{q}. \end{equation}

We parametrise $\boldsymbol{q}$ using the polar coordinates $(q,\vartheta )$ defined similarly to the $(r,\theta )$ coordinates, so that $\hat {\boldsymbol{{i}}}\boldsymbol{\cdot }\hat {\boldsymbol{q}}=\cos \vartheta$ . Substitution into (3.22) and recalling that $U$ is real-valued gives

(3.23) \begin{equation} U(\theta =0) = \frac {1}{2\pi }\iint _{(-\pi ,\pi )\times (0,\infty )} \frac {q^2\sin ^2\vartheta \cos \vartheta }{\sin ^2\vartheta +q^4}\sin (q R \cos \vartheta ) \, \mathrm{d}\vartheta \, \mathrm{d} q. \end{equation}

Noting that the integrand is invariant under the transformations $\vartheta \mapsto - \vartheta$ and $\vartheta \mapsto \pi -\vartheta$ we modify the integration domain:

(3.24) \begin{equation} U(\theta =0) = \frac {2}{\pi }\iint _{(0,\pi /2)\times (0,\infty )} \frac {q^2\sin ^2\vartheta \cos \vartheta }{\sin ^2\vartheta +q^4}\sin (q R \cos \vartheta ) \, \mathrm{d}\vartheta \, \mathrm{d} q. \end{equation}

With $\sin \vartheta$ non-negative we can change to the integration variable $\tau = q / \sin ^{1/2}\vartheta$ :

(3.25) \begin{equation} U(\theta =0) = \frac {2}{\pi }\iint _{(0,\pi /2)\times (0,\infty )} \frac {\tau ^2\sin ^{3/2}\vartheta \cos \vartheta }{1+\tau ^4} \sin \big (\tau R \cos \vartheta \sin ^{1/2}\vartheta \big ) \, \mathrm{d}\vartheta \, \mathrm{d} \tau . \end{equation}

Expressing $U$ as an iterated integral gives

(3.26) \begin{equation} U(\theta =0) = \frac {2}{\pi } \int _0^{\pi /2} K\big (R\cos \vartheta \sin ^{1/2}\vartheta \big )\sin ^{3/2}\vartheta \cos \vartheta \, \mathrm{d}\vartheta , \end{equation}

where

(3.27) \begin{equation} K(k) = \int _0^\infty \frac {\tau ^2\sin (k\tau )}{1+\tau ^4} \, \mathrm{d} \tau . \end{equation}

At the horizontal ray, where $\hat {\boldsymbol{e}}_r=\hat {\boldsymbol{{j}}}$ , we obtain from (3.17) that

(3.28) \begin{equation} U(\theta =\pi /2)= \frac {\mathrm{i}}{2\pi } \int _{\mathbb{R}^2} \frac {(\hat {\boldsymbol{{j}}}\boldsymbol{\cdot }\hat {\boldsymbol{q}})(\hat {\boldsymbol{{i}}}\boldsymbol{\cdot }\hat {\boldsymbol{q}})^2\,q}{1-(\hat {\boldsymbol{{i}}}\boldsymbol{\cdot }\hat {\boldsymbol{q}})^2+q^4} \exp \big (\mathrm{i} qR \, \hat {\boldsymbol{{j}}}\boldsymbol{\cdot }\hat {\boldsymbol{q}}\big )\,{\rm d}^2\boldsymbol{q}. \end{equation}

Employing the same parametrisation as before and noting that $\hat {\boldsymbol{{j}}}\boldsymbol{\cdot }\hat {\boldsymbol{q}}=\sin \vartheta$ we have here

(3.29) \begin{equation} U(\theta =\pi /2) = -\frac {1}{2\pi } \iint _{(-\pi ,\pi )\times (0,\infty )} \frac {q^2 \sin \vartheta \cos ^2\vartheta }{\sin ^2\vartheta +q^4}\sin (q R \sin \vartheta ) \, \mathrm{d}\vartheta \, \mathrm{d} q. \end{equation}

Noting that the integrand is again invariant under the transformations $\vartheta \mapsto - \vartheta$ and $\vartheta \mapsto \pi -\vartheta$ we modify the integration domain to obtain

(3.30) \begin{equation} U(\theta =\pi /2) = -\frac {2}{\pi } \iint _{(0,\pi /2)\times (0,\infty )} \frac {q^2\sin \vartheta \cos ^2\vartheta }{\sin ^2\vartheta +q^4} \sin (q R \sin \vartheta ) \, \mathrm{d}\vartheta \, \mathrm{d} q. \end{equation}

Changing to the integration variable $\tau = q / \sin ^{1/2}\vartheta$ gives (cf. (3.26))

(3.31) \begin{equation} U(\theta =\pi /2) = -\frac {2}{\pi } \int _0^{\pi /2} K\big (R\sin ^{3/2}\vartheta \big )\sin ^{1/2}\vartheta \cos ^2\vartheta \, \mathrm{d}\vartheta . \end{equation}

Both (3.26) and (3.31) are expressed in terms of $K$ , defined by the sine transform (3.27). An alternative representation for $K$ is obtained by employing contour integration in the complex $\tau$ plane of the function $\tau ^2 \mathrm{e}^{{\rm i} k \tau }/(1+\tau ^4)$ , with the closed contour consisting of the positive real axis, a circular arc of large radius in the first quadrant and the positive imaginary axis. Making use of the residue theorem and Jordan’s lemma readily yields,

(3.32) \begin{equation} K(k) = K_1(k) - K_2(k), \end{equation}

wherein

(3.33) \begin{equation} K_1(k) = \frac {\pi }{2\sqrt {2}} \operatorname {Im} \left \{ (1+\mathrm{i}) \mathrm{e}^{{\rm i} k (1+{i})/\sqrt {2}} \right \} \end{equation}

is contributed by the pole at $\tau = \mathrm{e}^{\pi \pi /4}$ and

(3.34) \begin{equation} K_2(k) = \int _0^\infty \frac {\tau ^2 \mathrm{e}^{-k\tau }}{1+\tau ^4}\, \mathrm{d}\tau \end{equation}

is the Laplace transform of $\tau ^2/(1+\tau ^4)$ .

In Appendix A we derive asymptotic properties of $K$ and evaluate certain integrals involving $K$ . These are needed next.

3.5. Asymptotic behaviour of the outer velocity

We can now determine the asymptotic behaviour of the vertical velocity (3.26) and the horizontal velocity (3.31) at small $R$ . To that end we make use of (A6), the approximation of $K$ in the limit $k\to 0$ . Substituting (A6) into (3.26) and performing the integration over $\vartheta$ we then obtain

(3.35) \begin{equation} U(\theta =0) \sim -\frac {R}{8} \left (\ln R + \gamma -\frac {5}{8} - \frac {3}{2}\ln 2\right ) \quad \text{as} \quad R\to 0. \end{equation}

Similarly, substituting (A6) into (3.31) and performing the integration over $\vartheta$ yields

(3.36) \begin{equation} U(\theta =\pi /2) \sim \frac {R}{8} \left (\ln R + \gamma -\frac {5}{8} - \frac {3}{2}\ln 2\right )\quad \text{as} \quad R\to 0. \end{equation}

We can also determine the asymptotic behaviour at large $R$ . Consider first the horizontal velocity (3.31) as $R\to \infty$ . A naive attempt to introduce (A7) and obtain what may appear like an $\textrm {ord}(R^{-3})$ approximation results in a divergent integral due to a strong singularity as $\vartheta \searrow 0$ . It is evident that it is necessary to address the local contribution to the integral from the region $\vartheta =\textrm {ord}(R^{-2/3})$ . In fact, a simple estimate reveals that the local contribution is $\textrm {ord}(R^{-1})$ and therefore dominates the $\textrm {ord}(R^{-3})$ global contribution. Defining the stretched variable $\nu = R\vartheta ^{3/2}$ we readily obtain

(3.37) \begin{equation} U(\theta =\pi /2) \sim -\frac {4}{3\pi R} \int _0^\infty K(\nu ) \, \mathrm{d}\nu \quad \text{as} \quad R\to \infty . \end{equation}

Upon substituting (A9) we obtain

(3.38) \begin{equation} U(\theta =\pi /2) \sim -\frac {1}{3 R} \quad \text{as} \quad R\to \infty . \end{equation}

The calculation of the comparable approximation of the vertical velocity (3.26),

(3.39) \begin{equation} U(\theta =0) \sim -\frac {24}{R^5} \quad \text{as} \quad R\to \infty , \end{equation}

requires more effort. It is derived in Appendix B.

When calculating (3.26) and (3.31) numerically we encounter difficulties when employing expression (3.27) due to the need to deal with an oscillatory integrand over a semi-infinite domain. We therefore employ the alternative form (3.32)–(3.34). In figure 2 we plot the variation with $R$ of the vertical velocity (3.26) and horizontal velocity (3.31) together with the corresponding small- $R$ approximations, (3.35) and (3.36), and large- $R$ approximations, (3.39) and (3.38). While approximation (3.38) becomes effective already for $R\approx 10$ , approximation (3.39) is only relevant once $U(\theta =0)$ becomes extremely small (and negative). This is illustrated in the inset of figure 2.

Figure 2.

Variation with $R$ of the outer velocity. Solid: $U(\theta =0)$ , as provided by (3.26), and $U(\theta =\pi /2)$ , as provided by (3.31). Dashed: corresponding small- $R$ approximations, (3.35) and (3.36). Dotted: large- $R$ approximation (3.38). The inset zooms in $U(\theta =0)$ , with the dotted curve portraying the large- $R$ approximation (3.39).

3.6. Back to the inner region

Using asymptotic matching, we may now evaluate the unknown constant $A$ appearing in the inner flow (3.6). While we only know the outer velocity on the vertical and horizontal axes, that is sufficient for matching purposes.

Consider indeed the vertical velocity. We match (3.35), the inner limit of the outer variable, with (3.10), the outer limit of the inner variable. Using (2.33a ) and (3.12b ) reveals that the $R\ln R$ terms trivially match; then, matching of the $R$ terms yields

(3.40) \begin{equation} A = \frac {1}{4}\ln \frac {1}{{\textit {Ra}}}-\gamma +\frac {5}{8}+\frac {3}{2}\ln 2. \end{equation}

We may similarly match (3.36) with (3.11). This reproduces (3.40).

With $A$ known, we have fully determined the inner flow. We now address the issue of a circular streamline. Substitution of (3.40) into (3.8) yields

(3.41) \begin{equation} r_0 \sim \frac {2^{3/2}\mathrm{e}^{5/8-\gamma }}{{\textit {Ra}}^{1/4}}. \end{equation}

It therefore follows that the streamline is predicted to be positioned at distances where the inner approximation underlying it is inapplicable. The impossibility of such a dividing streamline becomes evident from the need for asymptotic matching. This is best illustrated via use of an intermediate variable, say $\varrho ={\textit {Ra}}^\beta r$ with $0\lt \beta \lt 1/4$ . Considering the limit ${\textit {Ra}}\to 0$ with $\varrho$ fixed we obtain from (3.10) that

(3.42) \begin{equation} u(\theta =0) \sim \frac {\varrho }{8{\textit {Ra}}^\beta } \left [ (1/4-\beta )\ln \frac {1}{{\textit {Ra}}} -\ln \varrho -\gamma +\frac {5}{8}+\frac {3}{2}\ln 2\right ]\!. \end{equation}

With the first term in the brackets dominating in the intermediate limit, $u(\theta =0)$ is positive in the matching region. The prediction of a circular streamline, where the radial velocity changes sign, becomes irrelevant.

To illustrate the flow topology, consider the uniform approximation for the variation with $r$ of the radial velocity (recall (2.33a ) and (3.12b )):

(3.43) \begin{equation} u(r) + {\textit {Ra}}^{-1/4} U({\textit {Ra}}^{1/4}r) - \text{overlap}. \end{equation}

For a given $\textit {Ra}$ , this variation may be evaluated for the vertical velocity, where the overlap expression is provided by (3.10), and the horizontal velocity, where the overlap expression is provided by (3.11). In figure 3 we portray the uniform approximations in both directions for ${\textit {Ra}}=0.01$ . At this value we find from (3.41) that $r_0\approx 9.38$ . It is evident that no velocity reversal takes place in that neighbourhood.

Figure 3.

Uniform approximation at ${\textit {Ra}}=0.01$ to the radial velocity at both $\theta =0$ and $\theta =\pi /2$ .

3.7. Comparison with existing literature

As already mentioned, Alias & Page (Reference Alias and Page2017) analysed the 2-D problem in a different limit process, where the Grashof and Rayleigh numbers are comparably small. The leading-order problems in the inner regions are obviously identical. Moreover, the leading-order outer problems are also identical, since inertia terms are still absent.

The key difference between their work and the present one, associated with the outer analysis, is rather technical: Alias & Page (Reference Alias and Page2017) employ a stream-function formulation in their outer analysis, which is avoided in favour of an invariant notation in the present work. Consequently, Alias & Page (Reference Alias and Page2017) employed a Fourier transform only in the vertical direction, as opposed to the present 2-D transforms (3.14). Their derivation of the inverse transforms therefore becomes rather cumbersome. Comparison of equation (4.29) of Alias & Page (Reference Alias and Page2017) with the present (3.36) or (3.35) reveals that the constant $C_1$ in Alias & Page (Reference Alias and Page2017) is simply given by $\gamma -5/8 - (3/2)\ln 2$ . In their scheme, Alias & Page (Reference Alias and Page2017) chose to calculate it numerically (see their (4.32)).

In other aspects, the present 2-D analysis is complementary to that of Alias & Page (Reference Alias and Page2017). In contrast to the orientation-specific velocity approximations in the present work, Alias & Page (Reference Alias and Page2017) obtained large- $R$ integral expression (4.40) for the stream function, valid for all orientations. While that expression cannot be evaluated in exact form, its form suggests the consideration of a ‘horizontal’ jet-type structure whose scaling (in the present notation) is $X=O(Y^{1/3})$ . In that regime, Alias & Page (Reference Alias and Page2017) obtained a similarity representation for their large- $R$ expression in terms of a special function defined by Moore & Saffman (Reference Moore and Saffman1969). For $X\gg Y^{1/3}$ that representation reduces to a simple closed-form approximation, given by their (4.45a). Since it holds in the vicinity of the $X$ axis, it may be used to provide a large- $R$ velocity approximation on that axis. The resulting approximation is in agreement with the present (3.39). Approximation (4.45a) of Alias & Page (Reference Alias and Page2017) does not apply in the vicinity of the $Y$ axis, where the present analysis predicts (3.38).

The 2-D problem was also analysed by Thomas & Camassa (Reference Thomas and Camassa2023). The authors realised the two-scale structure in the limit ${\textit {Ra}}\to 0$ , but did not solve it using standard matched asymptotic expansions. Rather, they employ an Oseen-type linearisation in the entire fluid domain for small $\textit {Ra}$ , solving the linearised equations using Green functions. Their methodology resembles small-Reynolds-number analyses of the Oseen equation (Goldstein Reference Goldstein1929). It is interesting to note that Thomas & Camassa (Reference Thomas and Camassa2023), who were apparently unaware of Alias & Page (Reference Alias and Page2017), predicted (see their (4.8)) the same $R^{-5}$ decay rate for the radial velocity as that discovered by Alias & Page (Reference Alias and Page2017), and concluded that it breaks down at the horizontal axis. Their analysis there predicted the same $R^{-1}$ decay rate as in the present (3.38).