1. Introduction
It has been suggested that the diurnal vertical migration of marine organisms (termed the largest migration on Earth; see Martin et al. Reference Martin2020) could mix the oceans and is likely as important as the well-known energy sources due to winds and tides in contributing towards the meridional overturning circulation (Dewar et al. Reference Dewar, Bingham, Iverson, Nowacek, St Laurent and Wiebe2006; Katija & Dabiri Reference Katija and Dabiri2009; Nawroth & Dabiri Reference Nawroth and Dabiri2014; Houghton et al. Reference Houghton, Koseff, Monismith and Dabiri2018). This ‘biogenic mixing’ hypothesis concerns the motion of both passive particles and active swimmers at low (e.g. bacteria) to order-unity (e.g. copepods) Reynolds numbers (
$Re$). Such motion results in intricate flow patterns on scales comparable to and larger than the individual particles/swimmers, with the nature of the patterns being governed by various factors such as their swimming characteristics, fluid physical properties, interactions with other particles/swimmers, etc. (Katija Reference Katija2012; Houghton et al. Reference Houghton, Koseff, Monismith and Dabiri2018). An important element in the aforesaid hypothesis is the amount of fluid displaced by a single entity, either active or passive. This so-called drift volume can, in principle, act as a source of available potential energy which can then couple to motion on scales much larger than the swimmer size. The drift volume may be calculated by considering an initial plane of material points with the particle infinitely far upstream to begin with. The plane gets deformed as the particle translates towards it, and eventually to downstream infinity. For the initial plane being infinite in extent, the volume enclosed between it and the final deformed one is defined as the ‘total drift volume’ (Darwin Reference Darwin1953; Lighthill Reference Lighthill1956). For a sphere in potential flow, with a fore–aft symmetric disturbance velocity field of 
$O(1/r^3)$ (
$r$ being the radial distance from the sphere's centre; see figure 1a), the total drift volume is half the sphere's volume, this being equal to the added mass divided by the fluid density (Lighthill Reference Lighthill1956). In the opposite limit of a viscous fluid ambient (
$Re = 0$), the slower decaying 
${O}(1/r)$ disturbance velocity field (see figure 1b) causes the total drift volume to diverge over any finite time interval. Thus, it becomes necessary to define a ‘partial drift volume’, as the volume enclosed between initial and final material planes of a finite spatial extent, and with the initial plane located at a finite distance away from the sphere (Eames, Belcher & Hunt Reference Eames, Belcher and Hunt1994; Chisholm & Khair Reference Chisholm and Khair2018). Accounting for the effects of weak fluid inertia (see figure 1c) results in a faster 
${O}(1/r^2)$ source-flow-like decay of the velocity field, in almost all directions, at distances larger than the Oseen length (
$l_o = {O}(aRe^{-1})$). However, the original 
${O}(1/r)$ decay persists within a paraboloidal wake region behind the sphere, and as a result, despite there being no finite-time divergence, the drift volume does diverge linearly in the infinite time limit for any non-zero 
$Re$ (Subramanian Reference Subramanian2010; Chisholm & Khair Reference Chisholm and Khair2017).
Figure 1. A schematic (not to scale) of the disturbance flow fields due to a translating sphere (of radius 
$a$) in different physical scenarios. (a) Potential flow, (b) Stokes flow, (c) in the presence of weak fluid inertia and (d) with a weak background stratification. Here, 
$l_o = {O}(a Re^{-1})$, the Oseen length, is the scale at which inertial forces become comparable to viscous ones, and 
$l_c = {O}(a Ri_v^{-1/3})$, the stratification screening length, is the scale at which buoyancy forces become comparable to viscous ones (for large 
$Pe$). The focus of the current work is the flow field beyond 
$l_c$ in (d). The definitions of the non-dimensional numbers Re, Ri and Pe are as given in the fourth paragraph of § 1.
The divergence of the drift volume over a finite time interval, or for infinite time, led to the suggestion, by Katija & Dabiri (Reference Katija and Dabiri2009), of the aforementioned diurnal migration being a biogenic source of oceanic mixing. Although this migration is known to affect the transport of nutrients and oxygen, and the entrapment of CO
$_2$ in the ocean interior (Nowicki, DeVries & Siegel Reference Nowicki, DeVries and Siegel2022), its contribution to oceanic mixing has been proposed to be negligible based on various physical arguments (Visser Reference Visser2007; Subramanian Reference Subramanian2010; Wagner et al. Reference Wagner, Young and Lauga2014; Shaik & Elfring Reference Shaik and Elfring2024), with supporting evidence from numerical simulations (Wang & Ardekani Reference Wang and Ardekani2015). Apart from emphasizing the smallness of the scale at which individual marine organisms (zooplankton) input energy, the arguments rely on two points: (a) consideration of an active instead of a passive particle in the observations of Katija & Dabiri (Reference Katija and Dabiri2009), and (b) neglect of the stable stratification of the oceanic ambient by the said authors. Both of these should lead to a faster decay of the disturbance velocity field. Within a non-interacting framework, where one superposes single-particle contributions, this in turn should lead to a convergent drift volume. In being applicable to both passive and active particles, the ambient stratification is a more important factor than particle activity (Subramanian Reference Subramanian2010).
Mixing efficiencies in stratified fluids have been found to be very low for small Péclet numbers (
$Pe$), corresponding to a rapidly diffusing stratifying agent; in this limit, the energy injected by the particles/swimmers is primarily dissipated, rather than contributing to mixing or an increase in potential energy (Visser Reference Visser2007; Wagner et al. Reference Wagner, Young and Lauga2014). At the other extreme, that is, for large 
$Pe$, quantifying the extent of biogenic mixing requires a study of the drift volume mentioned above. This in turn requires an examination of the disturbance flow field induced by a particle translating vertically in a viscous stably stratified medium. The dominant contributions to the drift volume arise from fluid motion at distances from the particle/swimmer that are large compared with its size (Varanasi & Subramanian Reference Varanasi and Subramanian2022), and thus, the interest is in the far-field characteristics of the disturbance velocity in the parameter regime relevant to the oceanic scenario (see figure 1d).
Particle motion in a stratified fluid ambient is characterized by three non-dimensional numbers: the Reynolds number (
$Re = \rho U a/\mu$), the viscous Richardson number (
$Ri_v = \gamma a^3 g/\mu U$) and the Péclet number (
$Pe = U a/D$). Here, 
$U$ is the particle translational velocity, 
$a$ is the particle size (radius), 
$\rho$ is an appropriate reference fluid density, 
$\mu$ is the fluid viscosity, 
$\gamma$ is the absolute value of the background density gradient, g is acceleration due to gravity and 
$D$ may be taken as the salt diffusivity, keeping in mind the oceanic scenario. The viscous Richardson number is related to the Froude number (
$Fr = U/a\sqrt {g\gamma /\rho }$) as 
$Ri_v = Fr^{-2}Re$. In the vicinity of the particle, viscous forces are dominant, yielding the familiar 
${O}(1/r)$ Stokesian decay of the flow field, as already mentioned in the context of the drift-volume discussion above. A length scale (
$l_c$, measured from the particle), termed the ‘stratification screening length’, can be defined as the distance where the decaying viscous forces become comparable to buoyancy forces (Ardekani & Stocker Reference Ardekani and Stocker2010), leading to a departure from the 
${O}(1/r)$ decay above. The stratification screening lengths (
$l_c$) for small and large 
$Pe$ are 
${O}(a(Ri_vPe)^{-1/4})$ and 
${O}(aRi_v^{-1/3})$, respectively. These are analogous to the well-known Oseen length (
$l_o = {O}(aRe^{-1})$), where the departure from the Stokesian rate of decay occurs due to viscous and inertia forces becoming comparable in magnitude. The existence of both inertia and buoyancy forces, and their competition with viscous forces, leads to two flow regimes: (a) the ‘Stokes-stratification regime’ where 
$l_c \ll l_o$, and buoyancy effects dictate the flow behaviour at distances of order and larger than 
$l_c$, and (b) the ‘inertia-stratification regime’ where 
$l_c \gg l_o$, and inertial effects influence the flow much before buoyancy does.
For a migrating marine bacterium (for instance, Pseudoalteromonas haloplanktis), the aforementioned non-dimensional numbers are 
$Re \sim {O}(10^{-4})$, 
$Ri_v \sim {O}(10^{-12})$ and 
$Pe \sim {O}(10^{-2})$ (Wagner et al. Reference Wagner, Young and Lauga2014), and the bacterium therefore corresponds to the small-
$Pe$ limit of the Stokes-stratification regime. In contrast, typical zooplankton (for instance, copepods), with 
$Re \sim {O}(10^{-1})$, 
$Ri_v \sim {O}(10^{-8})$ and 
$Pe \sim {O}(10^2)$ (Kunze et al. Reference Kunze, Dower, Beveridge, Dewey and Bartlett2006; Varanasi & Subramanian Reference Varanasi and Subramanian2022, henceforth VS22), correspond to the large-
$Pe$ limit of the inertia-stratification regime. Although both scenarios are relevant for oceanic biomass, the latter is more important due to the abundance of zooplankton, and their participation in the diel vertical migration (Martin et al. Reference Martin2020; Wang & Ardekani Reference Wang and Ardekani2015). In figure 2(a), we present typical sizes and speeds of different classes of marine particulate matter and swimming organisms considered in the previous literature, along with estimates of the Reynolds and viscous Richardson numbers shown using blue and red contours. Figure 2(b) displays the same information on a plot whose coordinates are the Stokes-stratification screening length (ordinate) and the Oseen length (abscissa). In addition, in table 1, we have organized the earlier efforts in the literature by identifying the parameter space examined in each of these, along with the particular focus of the work. From this table, it is evident that, while there have been numerous efforts studying particle/drop motion in stratified fluids, the focus, especially for large-
$Pe$, had until recently been on the drag (Zvirin & Chadwick Reference Zvirin and Chadwick1975; Mehaddi et al. Reference Mehaddi, Candelier and Mehlig2018; Lee et al. Reference Lee, Fouxon and Lee2019) and torque corrections (Varanasi et al. Reference Varanasi, Marath and Subramanian2022) arising from the ambient stratification; with fewer works having examined the nature of fluid motion and the associated distortion of the isopycnals. Recent analytical efforts (Shaik & Ardekani Reference Shaik and Ardekani2020b; Varanasi & Subramanian Reference Varanasi and Subramanian2022) have studied the flow field due to a translating particle in the large-
$Pe$ limit, in an attempt to understand the resulting fluid drift. The current study has an analogous motivation, and again concerns the disturbance velocity and density fields induced by a vertically translating passive particle in the large-
$Pe$ limit of the Stokes-stratification regime (see figure 2b).
Figure 2. (a) Typical sizes and speeds of different groups/classes of marine particulate matter: Pseudoalteromonas haloplanktis (Wagner, Young & Lauga Reference Wagner, Young and Lauga2014), Chlamydomonas reinhardtii (More & Ardekani Reference More and Ardekani2020), marine snow (relevant range of values indicated by a purple box, see Turner Reference Turner2002), Volvox carteri (Drescher et al. Reference Drescher, Leptos, Tuval, Ishikawa, Pedley and Goldstein2009), Diaptomus minutus (Malkiel et al. Reference Malkiel, Sheng, Katz and Strickler2003), Daphnia magna (Noss & Lorke Reference Noss and Lorke2014), Artemia salina (Houghton et al. Reference Houghton, Koseff, Monismith and Dabiri2018), Euphasia pacifica (Kunze et al. Reference Kunze, Dower, Beveridge, Dewey and Bartlett2006) and Mastigias sp. (Katija & Dabiri Reference Katija and Dabiri2009). The blue and red lines correspond to constant Reynolds and viscous Richardson numbers, respectively. Here, 
$U \propto a^{-1}$ and 
$U\propto a^3$ along the blue and red lines, respectively. Note that U and a are characteristic speed and size, of the marine particulate matter, respectively. (b) A classification of the marine particulate matter in (a), in dimensionless terms, on a parameter plane comprising the ratio of the stratification screening length to particle size (
$l_c/a$) and the ratio of Oseen length scale to the particle size (
$l_o/l_a$); refer to third paragraph of § 1 for definitions of 
$l_c$ and 
$l_o$.
Table 1. The parameter regimes considered in the previous literature. Here, 
$Re$, 
$Ri$ and 
$Pe$ are the Reynolds number, viscous Richardson number and Péclet number, respectively, and are defined in the third paragraph of § 1. The flow regime considered by each study is given in the fifth column: the Stokes-stratification regime corresponds to 
$l_c \ll l_o$ and the inertia-stratification regime to 
$l_c \gg l_o$. Here, 
$l_c$ is the stratification screening length and 
$l_o$ is the Oseen length (see para 2, § 1). The focus of the different studies can be classified broadly into the flow-field description, drag calculation, drift-volume calculation and mixing efficiency estimation. The nature of the work is also indicated via superscripts 
$A$ (analytical), 
$N$ (numerical) or 
$E$ (experimental).
The disturbance flow field due to a point force in a stratified medium was first studied by List (Reference List1971) and later, for a point force and a point dipole, by Ardekani & Stocker (Reference Ardekani and Stocker2010), with both studies pertaining to the small-
$Pe$ limit of the Stokes-stratification regime. The point force in List (Reference List1971) represented a momentum jet, while it was taken to be a vertically translating particle in the latter study. In making the point-force approximation, the studies focused on the far field and followed the Fourier transform approach pioneered by Childress (Reference Childress1964) and Saffman (Reference Saffman1965). List (Reference List1971) used contour integration to simplify the inverse transform integral, whereas Ardekani & Stocker (Reference Ardekani and Stocker2010) used a fast Fourier transform technique for numerical integration. The resulting flow field was found to be fore–aft symmetric with a series of horizontal recirculating cells induced by buoyancy forces, and that acted to inhibit vertical motion of fluid parcels. In studying the stability of the disturbance flow field of a sphere translating in a viscous stratified ambient to turbulent fluctuations, Fouxon & Leshansky (Reference Fouxon and Leshansky2014) also derived expressions for the far-field velocity disturbance at small 
$Pe$. They too used contour integration to obtain the final integral expression for the streamfunction. Further, by solving the simplified integral along the stagnation streamline, the authors obtained a far-field algebraic decay of 
${O}(1/|z|^9)$.
Zhang et al. (Reference Zhang, Mercier and Magnaudet2019) studied the drag enhancement due to stratification, on a vertically translating sphere, using numerical simulations, and also presented a brief description of the flow field for different 
$Re$, and for small and large 
$Pe$. Their results showed a fore–aft asymmetry of the large-
$Pe$ flow field, caused by the convection of density perturbations. Lee et al. (Reference Lee, Fouxon and Lee2019) considered the flow field generated by both stationary and vertically translating spheres in a stratified ambient. For the latter case, expressions for the disturbance fields, at small 
$Pe$, were derived based on a Lorentz-type reciprocal identity involving an auxiliary problem where the sphere was replaced by singular forcings in the Stokes and convection–diffusion equations. The authors also numerically computed streamline and isopycnal patterns, albeit over domains of a limited spatial extent. While these patterns pertained to a range of 
$Re$, 
$Pe$ and 
$Ri_v$ (see table 2), and did show the emergence of a fore–aft asymmetry with an increase in 
$Pe$, the latter increase was accompanied by an analogous increase in 
$Re$, owing to the Prandtl number being fixed. As a result, one cannot isolate the asymmetry induced due to convection of density perturbations alone (as would be the case in the Stokes-stratification regime). Shaik & Ardekani (Reference Shaik and Ardekani2020b) examined the disturbance flow field and the partial drift volume, induced by a vertically translating sphere and a spherical swimmer at large 
$Pe$, using the fast Fourier transform technique originally employed in Ardekani & Stocker (Reference Ardekani and Stocker2010). The disturbance flow field again exhibited a fore–aft asymmetry, but rather surprisingly, appeared to decay at a slower rate in the upstream direction. The partial drift volume for the passive sphere, although finite, was found to be several times greater than the sphere's volume.
Table 2. The large-
$Pe$ wake region asymptotes; the wake region has a self-similar structure with the similarity variable defined by 
$\tilde {\eta } = \tilde {z}/\tilde {r}_t^{2/5}$.
Recently, VS22 studied features of the disturbance velocity and density fields due to a sphere (modelled as a point force) moving vertically through a stably stratified medium, in the absence of inertia, over a range of 
$Pe$. The authors numerically evaluated the inverse transform integrals to obtain the streamlines and isopycnals as a function of 
$\beta _\infty = Ri_v^{1/3}Pe^{-1}$, a dimensionless parameter that is the ratio of the large-
$Pe$ stratification screening length (
$a Ri_v^{-1/3}$) to the convective screening length (
$aPe^{-1}$). The authors also analytically examined the structure of the disturbance fields at distances much larger than the stratification screening length, identifying both a horizontal wake (as suggested by the recirculating cells in the earlier efforts of Zhang et al. (Reference Zhang, Mercier and Magnaudet2019) and Shaik & Ardekani Reference Shaik and Ardekani2020b) and a vertical jet behind the sphere, directed opposite to the sphere's motion, for 
$Pe = \infty$. While buoyancy forces caused the velocity field to decay rapidly in almost all directions, a much slower 
${O}(1/r)$ decay was observed within the jet region close to the rear stagnation streamline. Identification of the jet led to the emergence of a new screening length for finite 
$Pe$ (termed the ‘secondary screening length’, 
$l_s \sim {O}(aRi_v^{-1/2}Pe^{1/2})$), beyond which density diffusion started to smear out the jet and the ‘reverse-Stokeslet’ decay above gave way to an exponential decay of the disturbance velocity and density fields. It is worth noting that the slower decay of the disturbance velocity along the rear stagnation streamline, until distances of 
$O(l_s)$, is in contrast to Shaik & Ardekani (Reference Shaik and Ardekani2020b) above, who found the decay downstream to be faster. The latter is likely due to the limited spatial extent of the computational domain considered – figure 1(b) in Shaik & Ardekani (Reference Shaik and Ardekani2020b) shows the downstream axial velocity only until the first zero crossing (i.e. until 
$z/a \approx 0.9l_c$); VS22 also provided estimates for the total drift volume. A discussion of the literature considering drops, active/passive particles and the drag induced by their motion has been provided in the recent review by More & Ardekani (Reference More and Ardekani2023).
While VS22 was one of the first efforts to characterize the asymptotic structure of the far-field velocity disturbance at large 
$Pe$, herein, we show that this characterization is nevertheless incomplete, and that the disturbance fields have a richer asymptotic structure. The reverse jet identified in VS22 is not an isolated feature, but merely the central part of a configuration of concentric near-vertical recirculating cells that decrease in thickness, and increase in number, with increasing downstream distance. Identification of this columnar structure leads to the emergence of a tertiary screening length, corresponding to the length of the said cells. In VS22, the velocity along the rear stagnation streamline transitioned from a Stokeslet to a reverse Stokeslet, and then to an exponential decay, for any non-zero 
$\beta _\infty$. Herein, we show that the tertiary screening length leads to a further transition from the exponential to an eventual algebraic (
$1/z^7$) decay; this algebraic decay being the same as that at sufficiently long distances along the front stagnation streamline, albeit with a different numerical prefactor. Finally, although VS22 presented streamline and isopycnal patterns for a range of 
$\beta _\infty$, the spatial extents of these plots were limited. The limitation was quite severe for the smallest 
$\beta _\infty$ values, due to convergence issues in the numerical integration technique. For example, for 
$\beta _\infty = 10^{-5}$, isopycnals were only provided in a domain with dimensions of 16
$l_c$ and 5
$l_c$ in the vertical and horizontal directions, respectively; the flow field for 
$\beta _\infty = 0$ was not provided for the same reason. As already mentioned above, the nature of the 
$\beta _\infty = 0$ flow field at large distances is particularly important for the drift-volume calculation. As part of the current effort, we simplify the inverse transform integrals using contour integration, enabling us to obtain the flow field and isopycnals over much larger distances, and over a wide range of 
$Pe$, including 
$Pe = \infty \,(\beta _\infty = 0)$. The salient features of the flow field based on the analysis of VS22, and those emerging from the analysis in the following sections, are summarized in the pair of schematics in figure 3. It should be noted that, with the far-field restriction, and thence, the treatment of the translating particle as a point force at leading order, the results of both VS22 and the current study are agnostic to particle shape. Thus, the schematics in figure 3 remain valid for anisotropic particles, provided one replaces the velocity vector by the (point) force vector; the two being collinear for spheres. The flow field at distances of order the particle size is sensitive to shape anisotropy. Features of the flow field around an arbitrarily oriented translating spheroidal particle, for small 
$Pe$, have been discussed in Varanasi et al. (Reference Varanasi, Marath and Subramanian2022).
Figure 3. Schematics of the flow field induced by a sphere translating vertically in a viscous density-stratified ambient. Panel (a) is based on VS22, and (b) is based on the present study. In both figures, the primary screening length appears as a dashed circle with radius 
$l_c \sim {O}(aRi_v^{-1/3})$, and diffusion smears out transverse density variations, downstream of the sphere, at distances beyond the secondary screening length, 
$l_s \sim {O}(aRi_v^{-1/2}Pe^{1/2})$. The tertiary screening length, 
$l_t \sim {O}(aRi_v^{-1/2}Pe^{1/2}[\zeta + \frac {13}{4}\ln {\zeta } + ({13^2}/{4^2})({\ln \zeta }/{\zeta })])$ with 
$\zeta = \frac {1}{2}\ln ({\sqrt {{\rm \pi} }Ri_v^{-1}Pe^3}/{2160})$, marks the length of the columnar structure on the right.
In § 2, we first provide the governing equations, and then, the (formal) integral expressions for the velocity and density fields resulting from the linearized Fourier-transformed governing equations. Next, we discuss the reduction of these expressions to double integrals using contour integration, and the numerical method used to evaluate the double integrals; the actual simplified expressions, used for evaluation purposes, are given in Appendix A. In § 3.1, we briefly present the main results of VS22 which provide the context for the investigation here. In § 3.2, using asymptotic evaluations of the integrals, valid at large distances behind the translating sphere, we move beyond the findings of VS22 and characterize the columnar structure downstream of the translating sphere, providing expressions for its radial extent, the number of cells within, and the tertiary screening length (
$l_t$). The columnar structure ends at a distance of 
$O(l_t)$, leading to the disturbance fields transitioning from exponential to algebraically decaying functions of the downstream distance. The significance of the various screening lengths are illustrated via plots of the streamlines and isco-pycnals, in the downstream region, over a range of small 
$\beta _\infty$ values, including 
$\beta _\infty = 0$. In § 3.3, similar to VS22, we present the streamline and isopycnal patterns for a wide range of 
$\beta _\infty$, but over a much larger domain – the domain extends to 25
$l_c$ in the radial direction, and to 
$\max (l_t,100l_c)$ in the vertical (axial) direction. We present a summary of our main findings in § 4, and end with (i) a pair of figures, each with two halves, that help contrast the streamline and isopycnal patterns for large and small 
$\beta _\infty$; (ii) figures that show preliminary results from ongoing drift-volume calculations based on the velocity field obtained here.
2. Governing equations and the disturbance fields
We consider a sphere of radius 
$a$ settling with velocity 
$U$ across a uniformly (stably) stratified fluid, that is, with 
${\rm d}\rho /{\rm d}z = -\gamma$ (
$\gamma$ being a positive constant in the absence of the sphere). The non-dimensional governing equations, with the Boussinesq approximation, are given by (Zvirin & Chadwick Reference Zvirin and Chadwick1975; Ardekani & Stocker Reference Ardekani and Stocker2010; Varanasi & Subramanian Reference Varanasi and Subramanian2022)
\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = 0 \end{gather}
\begin{gather}Re[\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}] ={-}\boldsymbol{\nabla} p + \nabla^2\boldsymbol{u} - Ri_v\rho_f\hat{\boldsymbol{e}}_3, \end{gather}
\begin{gather}1-w+\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\rho_f = \dfrac{1}{Pe}\nabla^2\rho_f, \end{gather}
in a reference frame translating with the sphere. Here, 
$\rho _f$ is the density disturbance field, and 
$\boldsymbol {u}$ is the velocity field, 
$w$ being its vertical component. The length, velocity and density scales used to non-dimensionalize the relevant variables in (2.1)–(2.3) are 
$a$, 
$U$ and 
$\gamma a$, respectively. The no-slip and no-flux boundary conditions are
\begin{equation} \boldsymbol{u} = 0\quad\textrm{and}\quad\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\nabla}\rho_f ={-}\boldsymbol{n}\boldsymbol{\cdot}\hat{\boldsymbol{e}}_3\quad \textrm{at}\ r = 1, \end{equation}and the far-field boundary conditions are
\begin{equation} \boldsymbol{u} \to \hat{\boldsymbol{e}}_3\,(w \to 1), \quad p \to p_\infty, \quad \textrm{and}\quad \rho_f \rightarrow 0\quad \textrm{for}\ r \to \infty, \end{equation}
where 
$r = |\boldsymbol {x}|$, 
$\boldsymbol {x}$ is the position relative to the sphere's centre. The non-dimensional parameters in (2.1)–(2.3) have already been defined in the introduction section. In writing (2.1)–(2.5), we assume a quasi-steady state to have developed in the region of interest surrounding the translating sphere.
Our primary focus is on large 
$Pe$ in the inertialess limit (
$Re = 0$). For 
$Pe\to \infty$, the diffusive term in (2.3) is asymptotically small compared with the convective one. Hence, using 
$\boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {\nabla }\rho _f\sim (1-w)$, along with the Stokes equations, the large-
$Pe$ stratification screening length can be derived as 
$l_c \sim {O}(aRi_v^{-1/3})$. Using this as the characteristic length scale, and further, restricting our attention to the far-field (the outer region) where the sphere appears as a point force, (2.1)–(2.3) can be written in the form
\begin{gather} \tilde{\boldsymbol{\nabla}}\boldsymbol{\cdot}\tilde{\boldsymbol{u}} = 0, \end{gather}
\begin{gather}-\alpha_\infty\dfrac{\partial\tilde{\boldsymbol{u}}}{\partial\tilde{z}} ={-}\tilde{\boldsymbol{\nabla}} \tilde{p} + \tilde{\nabla}^2\tilde{\boldsymbol{u}} - [\tilde{\rho}_f+6{\rm \pi}\delta(\tilde{\boldsymbol{r}})]\hat{\boldsymbol{e}}_3, \end{gather}
\begin{gather}-\hat{\boldsymbol{e}}_3\boldsymbol{\cdot}\tilde{\boldsymbol{u}}+\dfrac{\partial \tilde{\rho}_f}{\partial\tilde{z}} = \beta_\infty\tilde{\nabla}^2\tilde{\rho}_f. \end{gather}
Here, 
$\tilde {\boldsymbol {r}} = Ri_v^{1/3}\boldsymbol {r}$ is the outer-region coordinate, 
$\tilde {\boldsymbol {u}} = Ri_v^{-1/3}(\boldsymbol {u}-\hat {\boldsymbol {e}}_3)$, 
$\tilde {p} = Ri_v^{-2/3}(p-p_\infty )$ and 
$\tilde {\rho }_f = \rho _f$. The delta-function forcing on the right-hand side of (2.7) replaces the no-slip condition above, and the system (2.6)–(2.8) is only required to satisfy far-field decay conditions. There are three screening lengths associated with the dimensionless parameters in (2.1)–(2.3) – the inertial screening length 
$aRe^{-1}$, the stratification screening length 
$aRi_v^{-1/3}$ and the screening length 
$aPe^{-1}$ associated with the convection of density perturbations (for small 
$Pe$). The ratios of these taken pairwise, viz. 
$\alpha _\infty = {Re}/{Ri_v^{1/3}}$ and 
$\beta _\infty = {Ri_v^{1/3}}/{Pe}$, are the parameters in (2.6)–(2.8). We study the nature of the disturbance flow field and isopycnals, as a function of 
$\beta _\infty$, assuming 
$\alpha _\infty = 0$.
Fourier transforming (2.6)–(2.8), one obtains the disturbance fields in terms of the following Fourier integrals:
\begin{gather} \tilde{\boldsymbol{u}}(\tilde{\boldsymbol{r}}) = \dfrac{-3}{4{\rm \pi}^2}\displaystyle\int\dfrac{({\rm i}k_3+\beta_\infty k^2)k^2\left(\hat{\boldsymbol{e}}_3-\dfrac{k_3\boldsymbol{k}}{k^2}\right)}{({\rm i}k_3+\beta_\infty k^2)k^4+k_t^2}{\rm e}^{{\rm i}\boldsymbol{k}\boldsymbol{\cdot}\tilde{\boldsymbol{r}}}\,{\rm d}\boldsymbol{k}, \end{gather}
\begin{gather}\tilde{\rho}_f(\tilde{\boldsymbol{r}}) = \dfrac{-3}{4{\rm \pi}^2}\displaystyle\int\dfrac{k_t^2}{({\rm i}k_3+\beta_\infty k^2)k^4+k_t^2}{\rm e}^{{\rm i}\boldsymbol{k}\boldsymbol{\cdot}\tilde{\boldsymbol{r}}}\,{\rm d}\boldsymbol{k}. \end{gather}
Here, 
$k^2 = k_1^2+k_2^2+k_3^2$ and 
$k_t^2 = k_1^2+k_2^2$. Exploiting axisymmetry with respect to the vertical direction (the direction of particle translation) in physical space and the absence of a swirl, one may define a Stokes streamfunction 
$\tilde {\psi }_s$, with the radial (
$\tilde {u}_r$) and axial (
$\tilde {u}_z$) velocity components related to 
$\tilde {\psi }_s$ as 
$\tilde {u}_r = -({1}/{\tilde {r}_t})({\partial \tilde {\psi }_s}/{\partial \tilde {z}})$ and 
$\tilde {u}_z = ({1}/{\tilde {r}_t})({\partial \tilde {\psi }_s}/{\partial \tilde {r}_t})$ in cylindrical coordinates (
$\tilde {r}_t$,
$\tilde {z}$). The expression for the Stokes streamfunction can then be written as
\begin{equation} \tilde{\psi}_s(\tilde{\boldsymbol{r}}) = \dfrac{3\tilde{r}_t {i}}{4{\rm \pi}^2}\displaystyle\int\dfrac{({\rm i}k_3+\beta_\infty k^2)k_2}{({\rm i}k_3+\beta_\infty k^2)k^4+k_t^2}{\rm e}^{{\rm i}\boldsymbol{k}\boldsymbol{\cdot}\tilde{\boldsymbol{r}}}\,{\rm d}\boldsymbol{k}. \end{equation} The triple integrals in 
$\boldsymbol {k}$ in (2.9)–(2.11) are first reduced to double integrals, which can then be evaluated numerically to obtain the velocity (or Stokes streamfunction) and density fields. This is done using spherical polar coordinates 
$(k,\theta, \phi )$ in Fourier space with the polar axis along 
$\hat {\boldsymbol {e}}_3$. VS22 first evaluated the integral over 
$\phi$, but this led to an integrand with an oscillatory dependence on 
$k$, in turn leading to numerical difficulties. To circumvent this issue, we first evaluate the 
$k$-integral using contour integration. The resulting double integral over 
$\theta$ and 
$\phi$ has finite integration limits, and the associated integrand is easier to evaluate. The detailed expressions for 
$\tilde {\psi }_s$ and density 
$\tilde {\rho }_f$, obtained in this manner, are given in Appendix A. We perform the numerical integrations appearing in these expressions ((A11), (A13), (A21) and (A23)) using the in-built integration routine ‘Clenshaw–Curtis rule’ in Mathematica, and vary the ‘Working Precision’ option to obtain the desired precision of the result. Because of the smallness of the result at large 
$\tilde {z}$, the variable precision arithmetic in Mathematica is well suited to this problem.
3. Towards a complete flow picture
In this section, we first present the important results of VS22 (in § 3.1), which helps set the stage for reporting our findings starting from § 3.2. In § 3.2, we derive expressions for the boundary of, and the number of recirculating cells within, the vertical columnar structure that arises behind the translating sphere in the non-diffusive (
$\beta _\infty = 0$) limit. Next, for non-zero 
$\beta _\infty$, we identify an additional length scale, the tertiary screening length, beyond which the columnar structure ceases to exist. In § 3.3, we put together the complete picture using streamline and isopycnal patterns as a function of 
$\beta _\infty$.
3.1. Varanasi & Subramanian (Reference Varanasi and Subramanian2022)
For 
$Re \ll Ri_v^{1/3}$ and 
$Pe \gg 1$, corresponding to the large-
$Pe$ limit of the Stokes-stratification regime, VS22 evaluated the 
$\phi$-integral in (2.10) and (2.11) to obtain the following double integrals for the Stokes streamfunction and density fields (see VS22-3.4 and VS22-3.2):
\begin{gather} \tilde{\psi}_s = \displaystyle \dfrac{-3\tilde{r}_t^2}{4\sqrt{\tilde{r}_t^2+\tilde{z}^2}} + \dfrac{3\tilde{r}_t}{2{\rm \pi}}\int_0^\infty {\rm d}k \int_0^{\rm \pi} {\rm d}\theta \dfrac{\sin^4\theta {\rm J}_1\left(k\tilde{r}_t\sin\theta\right)\exp({{\rm i}k\tilde{z}\cos\theta})}{k \left({\rm i}k^3\cos\theta + \beta_\infty k^4 +\sin^2\theta\right)}, \end{gather}
\begin{gather}\tilde{\rho}_f = \displaystyle \dfrac{-3}{2{\rm \pi}}\int_0^\infty {\rm d}k \int_0^{\rm \pi} {\rm d}\theta \dfrac{k^2\sin^3\theta {\rm J}_0\left(k\tilde{r}_t\sin\theta\right)\exp({{\rm i}k\tilde{z}\cos\theta})}{\left({\rm i}k^3\cos\theta + \beta_\infty k^4 +\sin^2\theta\right)}. \end{gather}
The first term in (3.1) is the Stokeslet contribution written out separately to facilitate numerical integration of the second term. The integrands in (3.1) (second term) and (3.2) decay as 
$1/k^{11/2}$ and 
$1/k^{5/2}$ for 
$k \to \infty$, for any non-zero 
$\beta _\infty$; the slower decay of the density integrand increases the difficulty of evaluation at larger spatial distances. The presence of 
$\tilde {r}_t$ and 
$\tilde {z}$ in the arguments of the Bessel and complex exponential functions implies an increasingly oscillatory integrand at large spatial distances, again contributing to an increased difficulty of calculation.
The need to retain the full complex exponential in the integrands of (3.1) and (3.2) is consistent with a fore–aft asymmetry of the disturbance fields. This was evident in the different numbers of recirculating cells in front of, and behind the translating sphere, in the streamline and isopycnal patterns obtained from numerically evaluating (3.1) and (3.2) – refer to figure 2 of VS22, and to figures 6 and 8 below. Since the method used made numerical integration difficult at large distances, the authors also conducted an asymptotic analysis of the disturbance fields on scales larger than 
${O}(a Ri_v^{-{1}/{3}})$. In the wake region where the flow is predominantly horizontal, the approximation 
$k_t \ll k_3$ was used (see (2.10) and (2.11)), rendering one of the integrals readily evaluable via contour integration. The remaining one-dimensional (1-D) integral revealed a self-similar wake structure with a width that grows as 
$\tilde {z} \sim \tilde {r}_t^{2/5}$. Evaluating this integral led to a quantitative characterization of the fore–aft asymmetry, and in particular, the different rates of decay of the disturbance fields in the upstream and downstream regions outside the wake; these algebraic asymptotes have been given in table 2.
The authors used the complementary approximation of the flow being predominantly vertical, corresponding to 
$k_t\gg k_3$, to examine what they termed the ‘reverse jet’ behind the sphere. The resulting 1-D integral permitted analytical evaluation only along the rear stagnation streamline 
$(\tilde {r}_t = 0, \tilde {z} > 0)$, giving the axial velocity and density disturbance fields in terms of modified Bessel functions of the second kind, of orders 1 (
$\tilde {u}_z \approx 3\beta _\infty ^{1/2}{\rm K}_1[2\beta _\infty ^{1/2}\tilde {z}]$) and 0 (
$\tilde {\rho }_f \approx -3{\rm K}_0[2\beta _\infty ^{1/2}\tilde {z}]$), respectively. The small-argument limiting form of 
${\rm K}_1$ gave rise to the reverse-Stokeslet decay 
$(\tilde {u}_z \sim 3/2\tilde {z}$), with its large-argument limiting form leading to an exponentially decaying velocity field (
$\tilde {u}_z \sim \exp (- 2\beta _\infty ^{1/2} \tilde {z})/\tilde {z}^{1/2})$. The transition between these decay regimes occurred at 
$\tilde {z} \sim {O}(\beta _\infty ^{-1/2})$ which, in dimensional terms, corresponds to the secondary screening length 
$l_s \sim {O}(aRi_v^{-1/2}Pe^{1/2})$. While the large-argument form of 
${\rm K}_0$ again led to an exponential decay of the density field (
$\tilde {\rho }_f \sim \exp (- 2\beta _\infty ^{1/2} \tilde {z})/\tilde {z}^{1/2}$), interestingly, the small-argument limit revealed a logarithmic behaviour; the density disturbance is logarithmically singular along the rear stagnation streamline for 
$\beta _\infty = 0$. The exponential decay for both 
$\tilde {u}_z$ and 
$\tilde {\rho }_f$ arises from diffusive effects that start to smear out the jet.
3.2. The vertical columnar structure behind the translating sphere
As mentioned in the previous subsection, VS22 obtained analytical forms for the disturbance fields along the rear stagnation streamline, allowing them to infer contrasting behaviour on either side of the secondary screening length. For non-zero 
$\tilde {r}_t$ and 
$\tilde {z} > 0$, the 1-D integrals for the Stokes streamfunction and the density field, with the aforementioned jet approximation (
$k_t \gg k_3$), are given by
\begin{gather} \tilde{\psi}_s= \displaystyle 3\tilde{r}_t\int_0^\infty {\rm d}k \dfrac{{\rm J}_1(k \tilde{r}_t)\exp({-\tilde{z}(\beta_\infty k_t^2 + 1/k_t^2)})}{k^4}, \end{gather}
\begin{gather}\tilde{\rho}_f={-}\displaystyle 3\int_0^\infty {\rm d}k \dfrac{{\rm J}_0(k \tilde{r}_t)\exp({-\tilde{z}(\beta_\infty k_t^2 + 1/k_t^2)})}{k}. \end{gather} While the structure of the flow field away from the rear stagnation streamline can be understood from numerically evaluating (3.3), as was done in VS22 (see figure 7 therein), important features are better understood using an analytical expression. The integral in (3.3) can be evaluated in closed form for 
$\beta _\infty = 0$, and for 
$\tilde {r}_t\tilde {z}^{1/2} \gg 1$, using a steepest-descent method (Bender & Orszag Reference Bender and Orszag2013); details of this calculation are given in Appendix B.1. One obtains the following expression for the Stokes streamfunction:
\begin{align} \tilde{\psi}_s|_{\beta_\infty = 0} &= 3\tilde{r}_t\displaystyle\int_0^\infty {\rm d}k \dfrac{{\rm J}_1(k \tilde{r}_t){\rm e}^{-\tilde{z}/k^2}}{k^4}, \end{align}
\begin{align} &\approx \dfrac{\sqrt{3}}{\sqrt[3]{2}}\left(\tilde{r}_t^{4/3}\tilde{z}^{{-}4/3}\right)\exp\left(-\dfrac{3}{2\sqrt[3]{4}} \left(\tilde{r}_t\tilde{z}^{1/2}\right)^{2/3}\right)\sin\left(\dfrac{3\sqrt{3}}{2\sqrt[3]{4}}\left(\tilde{r}_t\tilde{z}^{1/2}\right)^{2/3}+\dfrac{\rm \pi}{3}\right). \end{align}
The sinusoidal term in (3.6) is indicative of a cellular structure with the flow direction alternating between adjacent cells, and with the exponential pre-factor determining the manner of decay of these oscillations from the central to the peripheral cells. The cell widths may be correlated to the distance between successive zero crossings of the sine, and therefore scale as 
$\tilde {z}^{-1/2}$. The argument, 
$(\tilde {r}_t\tilde {z}^{1/2})^{2/3}$, of the sine implies that, at a fixed 
$\tilde {z}$, the cell widths, when projected onto the horizontal plane, increase with increasing 
$\tilde {r}_t$; on the other hand, a given cell becomes increasingly columnar (that is, approaches a vertical orientation) with increasing 
$\tilde {z}$. Only the central portion of the above columnar structure, corresponding to radial distances less than that of the first zero crossing, was identified in VS22 as the reverse jet. Although formally valid only for 
$\tilde {r}_t \tilde {z}^{1/2} \gg 1$, the location of the first zero crossing obtained from (3.6), given by 
$\tilde {r}_t\tilde {z}^{1/2} = (10{\rm \pi} \sqrt [3]{4}/9\sqrt {3})^{3/2}$, nevertheless gives a good approximation for the reverse-jet boundary.
In figure 4(a), the axial velocity field obtained from the steepest-descent approximation, (3.6), is compared with a numerical evaluation of the original inverse transform integral in (A8) using (A13) for 
$\beta _\infty = 0$. The comparison is as a function of 
$\tilde {z}$ for different 
$\tilde {r}_t$, starting from 
$\tilde {z} = 100$. The latter choice ensures that 
$\tilde {r}_t\tilde {z}^{1/2}$ remains large for the chosen non-zero 
$\tilde {r}_t$ values, and (3.6) therefore remains valid; for 
$\tilde {r}_t$ = 0, (3.6) remains inapplicable regardless of 
$\tilde {z}$, but the decay is known to have a reverse-Stokeslet character in this case. The numerical results, for non-zero 
$\tilde {r}_t$, are seen to agree with the analytical prediction up to a critical 
$\tilde {z}$ (that decreases with increasing 
$\tilde {r}_t$), there being a subsequent departure due to the numerical results transitioning to an algebraic decay for larger 
$\tilde {z}$ (implying the absence of recirculating cells at these distances). This algebraic decay is the same as the far-field asymptote obtained using the wake approximation, and is independent of 
$\tilde {r}_t$ – see table 2. In other words, for a given 
$\tilde {r}_t \neq 0$, the disturbance velocity field shifts from the oscillatory behaviour given by the jet approximation in (3.6), to an algebraic decay given by the far-field asymptote (
$\tilde {u}_z = 3240/\tilde {z}^7$) in table 2. Interestingly, the sum of the aforementioned approximations turns out to match very well with the numerical results over the entire range of 
$\tilde {z}$, and we exploit this fact to find the boundary of the columnar structure comprising the recirculating cells. This boundary is where the flow changes from being predominantly vertical (for smaller 
$\tilde {z}$) to predominantly horizontal (for larger 
$\tilde {z}$), and is identified by equating the amplitude of the oscillations in (3.6) to the algebraic asymptote. That is,
\begin{equation} \dfrac{1620 \tilde{r}_t^2}{\tilde{z}^7} = \dfrac{\sqrt{3}}{\sqrt[3]{2}}\dfrac{(\tilde{r}_t\tilde{z}^{1/2})^{4/3}}{\tilde{z}^2}\exp\left(-\dfrac{3}{2\sqrt[3]{4}} (\tilde{r}_t\tilde{z}^{1/2})^{2/3}\right). \end{equation}
The critical 
$\tilde {r}_t$ corresponding to the boundary can be obtained by solving (3.7) perturbatively, in the limit 
$\tilde {r}_t\tilde {z}^{1/2} \gg 1$, after taking a logarithm on both sides. One obtains (see Appendix B.2 for details)
\begin{equation} \tilde{r}_J = \dfrac{1}{\sqrt{\tilde{z}A^3}}\left[\ln(\tilde{z}^6/B) - \ln\left(\dfrac{1}{A}\left(\ln(\tilde{z}^6/B)\right)\right)\right]^{3/2}, \end{equation}
to second order in the perturbation expansion, with 
$A = 3/2\sqrt [3]{4}$ and 
$B = 1620\sqrt [3]{2}/\sqrt {3}$. The radial distance to the boundary (
$\propto \tilde {z}^{-1/2}(\ln \tilde {z})^{3/2}$ at leading logarithmic order), at a fixed downstream distance, is seen to be logarithmically larger than the width of an individual cell (
$\propto \tilde {z}^{-1/2}$), which points to the number of recirculating cells 
$N(\tilde {z})$ being a logarithmically increasing function of 
$\tilde {z}$. Here, 
$N$ can be obtained by equating the argument of the sinusoidal term in (3.6) to 
$(N+1){\rm \pi}$, recognizing that the number of zero crossings is one more than the cell number. Substituting 
$\tilde {r}_J$ from (3.8) leads to
\begin{equation} N(\tilde{z}) = \dfrac{\sqrt{3}}{{\rm \pi}\sqrt{A}}\left[\ln(\tilde{z}^6/B) - \ln\left(\dfrac{1}{A}\ln(\tilde{z}^6/B) \right)\right]^{3/2}-\dfrac{2}{3}, \end{equation}again to second order in the logarithmic expansion. As expected, (3.9) increases logarithmically with downstream distance.
Figure 4. The downstream axial velocity disturbance (
$\tilde {u}_z$) and streamline pattern for 
$\beta _\infty = 0$. Panel (a) shows 
$|\tilde {u}_z|$ plotted as a function of 
$\tilde {z}\,({>}0)$. The numerical results (bright-coloured solid curves), obtained from evaluating the inverse transform integral in (2.11), are compared withthe asymptotic approximation derived from (3.6) (pale-coloured dashed curves with dots added for better differentiation) for different non-zero 
$\tilde {r}_t$. The reverse-Stokeslet decay for 
$\tilde {r}_t = 0$ is shown as a solid grey line, while the (common) algebraic asymptote, 
$3240/\tilde {z}^7$, appears as a dashed grey line. Panel (b) shows the streamline pattern for 
$\tilde {z} > 1$, depicting the columnar structure behind the translating sphere. The boundary of the columnar structure, which is the red curve connecting the ends of the cell boundaries (black curves) above the wake region (
$\tilde {z} \gtrsim 100$), is given by (3.8).
Expressions (3.8) and (3.9) serve to characterize the columnar structure downstream that, despite narrowing down with increasing 
$\tilde {z}$, nevertheless encloses an increasing number of recirculating cells. This is depicted via the downstream portion of the streamline pattern (for 
$\tilde {z} \geqslant 1$) for 
$\beta _\infty = 0$ in figure 4(b), again obtained from a numerical integration of (2.11) using (A8) and (A13). The red curve in this figure corresponds to (3.8), and evidently demarcates the columnar structure from the predominantly horizontal streamline pattern outside. The black curves in the figure correspond to the zero crossings of the sine function in (3.6), obtained by equating its argument to 
$(m+1){\rm \pi} \,(m=1,2,3,\ldots )$; this leads to 
$\tilde {r}_t\tilde {z}^{1/2} = [(m{\rm \pi} +{2{\rm \pi} }/{3})({2\sqrt [3]{4}}/{3\sqrt {3}})]^{3/2}$.
For 
$\beta _\infty = 0$, the axial velocity along the rear stagnation streamline exhibits the reverse-Stokeslet behaviour for all distances (much) greater than the primary screening length (
$\tilde {z} \gg 1$), and this is indicated by the solid grey line in figure 4(a). For any non-zero 
$\tilde {r}_t$, however, there is an eventual transition to the 
$O(1/\tilde {z}^7)$ wake asymptote, corresponding to the dashed grey line in the said figure. In contrast, for 
$\beta _\infty \neq 0$, the axial velocity along the rear stagnation streamline changes from the reverse-Stokeslet decay to an exponential decay beyond the secondary screening length (
$\tilde {z} \sim \beta _\infty ^{-1/2}$), and then at still larger distances, reverts to the aforementioned wake asymptote. The plot of the axial velocity along the rear stagnation streamline in figure 5(a) shows the aforementioned pair of transitions for different small 
$\beta _\infty$. If we take the second transition to define a tertiary screening length (
$l_t$), then the exponential decay is essentially an intermediate asymptotic regime acting to connect the 
$O(1/\tilde {z})$-asymptote at distances smaller than 
$l_s$, to the 
$O(1/\tilde {z}^7)$-asymptote at distances larger than 
$l_t$. An estimate of 
$l_t$ can be obtained by equating the large-argument form of 
${\rm K}_1(2\beta _\infty ^{1/2}\tilde {z})$, that characterizes the exponential decay phase, to the algebraic wake asymptote. That is, 
$3\beta _\infty ^{1/2}{\rm K}_1(2\beta _\infty ^{1/2}\tilde {z}_t) = {3240}/{\tilde {z}_t^7}$ with 
$\tilde {z}_t \beta _\infty ^{1/2} \gg 1$, 
$\beta _\infty \ll 1$; here, 
$\tilde {z}_t = l_t/l_c$ is the axial location corresponding to the tertiary screening length in units of 
$l_c$. Using 
$\lim _{\tilde {z}_t\beta _\infty ^{1/2} \gg 1} {\rm K}_1(2\beta _\infty ^{1/2}\tilde {z}_t) = \sqrt {{\rm \pi} }\exp (- 2\beta _\infty ^{1/2} \tilde {z})/(4\beta _\infty ^{1/2}\tilde {z})^{1/2}$, one may solve the resulting transcendental equation in a manner similar to 
$\tilde {r}_J$ above, by taking logarithms on both sides. This gives
\begin{equation} l_t = l_c\beta_\infty^{{-}1/2}\left(\zeta + \dfrac{13}{4}\ln\zeta + \dfrac{13^2}{4^2}\dfrac{\ln\zeta}{\zeta}\right), \end{equation}
to third order in the small parameter 
$\zeta ^{-1}$, where 
$\zeta = \frac {1}{2}\ln ({\sqrt {{\rm \pi} }}/{2160\beta _\infty ^3}) \gg 1$; see Appendix C. Using the original dimensionless parameters with 
$l_c = aRi_v^{-1/3}$, one has 
$l_t \sim {O}(aRi_v^{-1/2}Pe^{1/2}\ln [\sqrt {{\rm \pi} }Ri_v^{-1}Pe^3/2160] )$ to leading logarithmic order. We have verified that the three-term expansion in (3.10) is necessary for a precise estimate of the second exponential-to-algebraic transition mentioned above; the first term alone leads to a significant underestimate. In figure 5(a), the tertiary screening lengths, obtained from (3.10), are indicated by (vertical) dashed grey lines, and correlate very well with the transition from an intermediate exponential to the eventual algebraic decay. The secondary screening lengths, given by 
$l_s = 0.5623\beta _\infty ^{-1/2}l_c$ and marking the beginning of the exponential decay phase, are also shown. The numerical prefactor in the said expression is obtained by finding the difference between the axial velocities for non-zero 
$\beta _\infty$ and 
$\beta _\infty = 0$, along the rear-stagnation streamline. This difference increases with 
$\tilde {z}$ to begin with, owing to the onset of diffusion effects, attains a maximum and then decreases. The latter decrease is because the behaviour for large 
$\tilde {z}$ is dictated by the reverse-Stokeslet decay for 
$\beta _\infty = 0$, the contribution due to the non-zero 
$\beta _\infty$ becoming exponentially small. The location of the aforementioned maximum gives the numerical prefactor.
Figure 5. The absolute value of the axial velocity (
$|\tilde {u}_z|$) (a) along the rear stagnation streamline; and (b) along the third (continuous), fifth (dashed) and seventh (dotted) zero crossings (of the Stokes streamfunction), for different 
$\beta _\infty$, including 
$\beta _\infty = 0$. For all non-zero 
$\beta _\infty$, the plots depict the transition from a reverse-Stokeslet decay (
$3/2\tilde {z}$) to the far-field asymptote (
$3240/\tilde {z}^7$) at large 
$\tilde {z}$. The secondary and tertiary screening lengths, marking the beginning and end of the exponential decay phase, are shown for each 
$\beta _\infty$.
Figure 5(b) shows the behaviour of the axial velocity field, as a function of distance, along curves corresponding to the third, fifth and seventh zero crossings of 
$\psi _s$; this is done for different small 
$\beta _\infty$, including 
$\beta _\infty = 0$. In the latter case, the loci of the zero crossings are the black curves shown earlier in figure 4(b). The first zero crossing corresponds to the rear stagnation streamline which serves as the axis of the downstream columnar structure. The flow within the columnar structure reverses direction from one cell to the next, and the zero crossing numbers above correspond therefore to the ‘centrelines’ of every alternate cell, starting from the central reverse jet, with flow directed away from the sphere. Interestingly, for a given 
$\beta _\infty$, the decay has a similar character for the different cells, although its nature is essentially different for zero and non-zero 
$\beta _\infty$. For 
$\beta _\infty = 0$, 
$\tilde {u}_z$ decays as the inverse of the distance along the centreline for all three cells shown, although with smaller numerical pre-factors for the peripheral ones. For 
$\beta _\infty \neq 0$, the inverse-centreline-distance decay only occurs until the secondary screening length, and an exponential decay ensues thereafter. Note that all curves in figure 5(b) start from a finite downstream distance, with this distance being larger for the peripheral cells – this is consistent with the nature of the columnar structure for 
$\beta _\infty = 0$ as evident from the streamline pattern in figure 4(b) above; and for non-zero 
$\beta _\infty$ as shown in § 3.3 (see figures 5 and 7 below). Note also that all the curves for non-zero 
$\beta _\infty$ terminate at the tertiary screening length, as must be the case. Since the latter length scale diverges for 
$\beta _\infty \rightarrow 0$, all three black lines in figure 5(b) continue to 
$\tilde {z} = \infty$.
3.3. The complete streamline and isopycnal patterns
Although different flow features at large distances from the sphere have been predicted using asymptotic methods, as illustrated in § 3.2, the extent to which these approximations are applicable can only be ascertained by a full numerical integration. Towards this end, we now present the complete streamline and isopycnal patterns from a numerical evaluation of the original inverse transform integrals in (2.9) and (2.10) using (A8) and (A16).
Figure 6 shows the streamline patterns around the translating sphere (a point force at the origin) over a range of 
$\beta _\infty$, from the diffusion-dominant limit (
$\beta _\infty = 10$; figure 6a) up until the convection-dominant one (
$\beta _\infty = 10^{-4}$; figure 6f); black arrows attached to select streamlines in each pattern indicate the flow direction. The portions of the patterns in the rather limited central rectangular region, demarcated by black dashed lines in each of figures 6(a)–6(f), were the ones obtained earlier by VS22 (see figures 8 and 9 therein). The increasing importance of convective effects, accompanying a decrease in 
$\beta _\infty$, leads to a progressive departure from the fore–aft symmetry that characterizes the streamline pattern in the limit 
$\beta _\infty \rightarrow \infty$. One manifestation of this asymmetry is a significant reduction in the number of horizontal recirculating cells upstream of the translating sphere (corresponding to negative 
$\tilde {z}$) – from 
$7$ for 
$\beta _\infty = 10$, to 3 for 
$\beta _\infty = 1$, to 
$1$ for 
$\beta _\infty < O(10^{-2})$. The radial extent of all streamline patterns shown is 
$25 aRi_v^{-1/3}$ which is sufficient for 
$\beta _\infty \lesssim 10^{-1}$ in the sense that the number of (horizontal) recirculating cells should remain unchanged for 
$\tilde {r}_t \gtrsim 25$. This constancy of cell number is consistent with the self-similar structure of the large-
$Pe$ horizontal wake at these radial distances; the wake grows as 
$\tilde {z} \propto \tilde {r}_t^{2/5}$ for large 
$Pe$, as shown in VS22, with its structure being a function of the similarity variable 
$\tilde {z}/\tilde {r}_t^{2/5}$. For 
$\beta _\infty \lesssim 10^{-1}$, the number of recirculating cells at larger distances therefore corresponds to the (fixed) number of zero crossings of the similarity solution, when plotted as a function of the similarity variable above. Further, the far-field asymptotes of the large-
$Pe$ similarity solution, given in table 2, give the decay rates outside of the wake recirculating cells, in both the upstream and downstream directions; for the Stokes streamfunction, the upstream and downstream asymptotes are given by 
$1620\tilde {r}_t^2/\tilde {z}^7$ and 
$-1080\tilde {r}_t^2/\tilde {z}^7$, respectively.
Figure 6. Streamline patterns in an axial (half) plane for 
$\beta _\infty = (a)$ 
$10$, (b) 
$1$, (c) 
$10^{-1}$, (d) 
$10^{-2}$, (e) 
$10^{-3}$ and (f) 
$10^{-4}$. The sphere is located at 
$(\tilde {r}_t,\tilde {z}) \equiv (0,0)$, with 
$\tilde {r}_t$ and 
$\tilde {z}$ being scaled by the primary screening length 
$l_c = {O}(aRi_v^{-1/3})$. Here, the dense blue bands mark the zero crossings of the Stokes streamfunction. The streamline patterns within the central regions, enclosed by the black dashed rectangles, are those originally obtained by VS22. The black continuous curves in (a,f) are wake boundaries corresponding to the small- and large
$-Pe$ similarity solutions, respectively.
The wake grows differently for small 
$Pe$, as 
$\tilde {z} \propto \tilde {r}_t^{1/3}$ for distances larger than 
$a(Ri_vPe)^{-1/4}$. An implication of the different screening lengths in the small and large-
$Pe$ limits, for the streamline patterns in figure 6, is that, for a fixed 
$\tilde {r}_t$, the effective size of the domain decreases as 
$\beta _\infty ^{-{1}/{4}}$ for 
$\beta _\infty > 1$ (
$\beta _\infty ^{{1}/{4}}$ being the ratio of the two screening lengths). As a result, for a given choice of the maximum 
$\tilde {r}_t$, the domain will become small enough for sufficiently large 
$\beta _\infty$, in the sense of the numerically generated streamline pattern not connecting to the pattern consistent with the far-field wake-similarity solution. This insufficiency is starting to be evident in figure 6(a), for 
$\beta _\infty = 10$, where the black curves (
$\tilde {z} \propto \tilde {r}_t^{{1}/{3}}$) denoting the self-similar growth of the small-
$Pe$ wake are seen to depart noticeably from the contours that mark the final zero crossings in the upstream and downstream directions. The analogous curves for large 
$Pe$ (
$\tilde {z} \propto \tilde {r}_t^{{2}/{5}}$) exhibit a better match with the zero crossing contours in figure 6(f).
Starting from 
$\beta _\infty = 10^{-2}$, corresponding to figure 6(d), one sees the emergence of an additional set of nearly vertical recirculating cells behind the sphere. Only a limited portion of the central cell in figure 6(f), for 
$\beta _\infty = 10^{-4}$, falls within the domain accessed in VS22 (the dashed black rectangle, as mentioned above), which gave the impression of an isolated rearward jet-like structure. The larger domains in figures 6(d)–6(f) make it evident that this reverse jet is only the central portion of an ensemble of concentrically arranged nearly vertical cells. With a decrease in 
$\beta _\infty$ from 
$10^{-2}$ to 
$10^{-4}$, the number of vertical recirculating cells in the columnar structure increases from 
$4$ to 
$8$ over the axial extent (
$100 l_c$) of the domain examined.
To better understand the structure and spatial extent of the downstream columnar structure in figures 6(d)–6(f), we focus next on the streamline patterns for the smaller 
$\beta _\infty$ values (
${=}10^{-3}$, 
$10^{-4}$, 
$10^{-5}$ and 
$0$), and on a much longer region downstream of the sphere (
$\tilde {z}$ extending up to 
$10^4$). The depiction of the patterns in figures 7(a)–7(d) is now via a semi-log plot, with the 
$\tilde {z}$-axis on a logarithmic scale. For 
$\tilde {z} \lesssim 20$, the streamline patterns in figures 7(a)–7(d) show a largely identical pattern of horizontal recirculating cells, implying that the wake structure has converged to its limiting form for 
$\beta _\infty \rightarrow 0$. In contrast, the emergence of vertical recirculating cells with increasing 
$\tilde {z}$, and the associated pattern of streamlines, continues to be sensitively dependent on 
$\beta _\infty$ even in the said limit. Figure 7(d), for 
$\beta _\infty = 0$, shows a continuous increase in the number of recirculating cells with increasing 
$\tilde {z}$. The radial extents of the individual cells decrease as 
$\tilde {z}^{-1/2}$, with the number of cells increasing in a manner consistent with the prediction given by (3.9) – from 11 cells at 
$\tilde {z} = 100$ to 24 cells at 
$\tilde {z} = 10^4$. For non-zero 
$\beta _\infty$, as the effects of density diffusion start to become important at the secondary screening length, 
$l_s = {O}(l_c\beta _\infty ^{-1/2})$, the recirculating cell boundaries begin to deviate from those for 
$\beta _\infty = 0$, with the cells, especially the ones closer to the rear stagnation streamline, becoming vertical and having a nearly constant width at larger distances. This change due to density diffusion occurs at a smaller downstream distance for the larger 
$\beta _\infty$ values, as a result of which the recirculating cells are wider in these cases. All the recirculating cells for non-zero 
$\beta _\infty$ end at 
$\tilde {z} \sim O(l_t)$ with 
$l_t$ given by (3.10). It is worth reiterating that both the secondary and tertiary screening lengths diverge in the non-diffusive limit (
$\beta _\infty = 0$). Figure 7(e) superposes the columnar structures for different 
$\beta _\infty$, validating the aforementioned significance of the secondary and tertiary screening lengths.
Figure 7. Streamline patterns in an axial (half) plane, and over an extended downstream region (
$1 < \tilde {z} < 10^4$), showing the columnar structure for 
$\beta _\infty = (a)$ 
$10^{-3}$, (b) 
$10^{-4}$, (c) 
$10^{-5}$ and (d) 
$0$. The sphere is at 
$(\tilde {r}_t,\tilde {z}) \equiv (0,0)$, with 
$\tilde {r}_t$ and 
$\tilde {z}$ being scaled using the primary screening length 
$l_c = {O}(aRi_v^{-1/3})$. The dense blue bands mark the zero crossings of the Stokes streamfunction; the region enclosed by black dashed rectangles correspond to the downstream patterns accessed in VS22. The inset in (b) is a zoomed-in version of the red dashed rectangle, showing fixed points at the end of recirculating cells that enable the transition from the cellular to a largely horizontal flow with increasing downstream distance. (e) Features of the columnar structures, for different non-zero 
$\beta _\infty$, compared with the limiting case of 
$\beta _\infty = 0$; horizontal dashed and dot-dashed lines indicate secondary(
$l_s$) and tertiary (
$l_t$) screening lengths, respectively.
Figure 8 shows colour plots (in red) of the perturbation density (
$\tilde {\rho }_f$), on which are superimposed blue isopycnal contours. The isopycnal patterns correspond to the same 
$\beta _\infty$ values as in figure 6 and have the same spatial extents; dashed black rectangles again mark out the domains accessed previously in VS22, with the size of these regions reducing starting from 
$\beta _\infty =10^{-2}$, owing to numerical convergence issues already mentioned in the introduction. Similar to the streamline patterns, there exists a set of horizontal cells comprising the wake region in each of figures 8(a)–8(f), that is largely insensitive to 
$\beta _\infty$ in the limit 
$\beta _\infty \rightarrow 0$ (compare figure 8e–f); and a set of vertical cells comprising the downstream columnar structure that emerge for 
$\beta _\infty \lesssim O(10^{-2})$, and that continue to lengthen and increase in number even as 
$\beta _\infty$ decreases to zero. The decay of the density perturbation outside the aforementioned cellular regions is algebraic, and given in table 2. The colour plots in figure 8 also help identify the buoyant envelope of fluid (bright red) that surrounds the sphere (the point force at the origin). While this envelope remains localized for larger 
$\beta _\infty$ (see figure 8a–c), it starts to get stretched out in the downstream direction with decreasing 
$\beta _\infty$, indicative of the sphere dragging along a column of buoyant fluid behind it.
Figure 8. Isopycnal patterns in an axial (half) plane at 
$\beta _\infty = (a)$ 
$10$, (b) 
$1$, (c) 
$10^{-1}$, (d) 
$10^{-2}$, (e) 
$10^{-2}$ and (f) 
$10^{-3}$. The sphere is located at the origin 
$(\tilde {r}_t,\tilde {z}) = (0,0)$. Here, 
$\tilde {r}_t$ and 
$\tilde {z}$ are non-dimensionalized using the primary screening length 
$l_c = {O}(aRi_v^{-1/3})$, and the dense bands correspond to the zero crossings (i.e. 
$\tilde {\rho }_f = 0$) of the density. The portions of the patterns in the region enclosed by the black dashed rectangles correspond to those obtained by VS22.
The structure of the buoyant tail mentioned above is better seen in figures 9(a)–9(d) where, on account of the domain extending until the tertiary screening length in the downstream direction, one can also see the columnar structure that surrounds the buoyant tail, in its entirety, down until 
$\beta _\infty = 10^{-5}$; features of this structure are analogous to those seen in the streamline patterns in figure 7. The elongated buoyant tail above, with a length of 
$O(l_s)$ is one of the most important manifestations of the fore–aft asymmetry that develops in the isopycnal patterns with decreasing 
$\beta _\infty$. Further, as mentioned earlier, the maximum negative value of density in this region should diverge logarithmically along the entire rear stagnation streamline for 
$\beta _\infty \rightarrow 0$, as may also be inferred from the expression for the density field obtained within the jet approximation (
$\tilde {\rho }_f = -3{\rm K}_0[2\beta _\infty ^{1/2}\tilde {z}]$; see VS22).
Figure 9. Isopycnal patterns in an axial half (plane) (and positive 
$\tilde {z}$) showing the downstream columnar structure at 
$\beta _\infty = (a)$ 
$10^{-3}$, (b) 
$10^{-4}$, (c) 
$10^{-5}$ and (d) 
$0$. The sphere is located at the origin 
$(\tilde {r}_t,\tilde {z}) = (0,0)$. Both 
$\tilde {r}_t$ and 
$\tilde {z}$ are non-dimensionalized using the primary screening length 
$l_c = {O}(aRi_v^{-1/3})$. Here, the dense bands correspond to the zero crossings of perturbation density (i.e. 
$\tilde {\rho }_f = 0$). The isopycnal patterns in the regions enclosed by the black dashed rectangles are those obtained by VS22. (e) The columnar structures at various non-zero 
$\beta _\infty$, compared with the case of 
$\beta _\infty = 0$. The horizontal dashed and dot-dashed lines indicate secondary (
$l_s$) and tertiary (
$l_s$) screening lengths for different 
$\beta _\infty$.
4. Conclusions
In this paper, we have considered the disturbance flow field due to a sphere translating vertically in a viscous stratified fluid ambient, the emphasis being on the convection-dominant limit (
$Pe \gg 1$) within the Stokes-stratification regime; the latter regime is where buoyancy forces are dominant over inertial forces, and therefore first balance viscous forces beyond a (primary) screening length (
$a Ri_v^{-1/3}$). VS22 showed that the flow field in this limit, at large distances from the sphere, consists of a set of horizontal recirculating cells (the wake) surrounding the sphere in the equatorial plane, and a vertical jet in the rear, directed opposite to the sphere's translation. Interestingly, the authors speculated on the existence of additional recirculating cells in the region downstream, and the possibility of the reverse jet being part of a more elaborate structure behind the translating sphere (see bottom of p. 18 in VS22). Our results, both analytical and numerical, confirm this speculation, and help quantitatively characterize this more elaborate columnar structure downstream of the sphere. We show that, although this structure extends to downstream infinity in the absence of diffusion (
$\beta _\infty = 0$), for any non-zero 
$\beta _\infty$, it has a finite length, which we term the tertiary screening length, and that is 
${O}[Ri_v^{-1/2}Pe^{1/2}\ln (Ri_v^{-1}Pe^3)]$ to leading logarithmic order. The analysis in VS22 revealed an exponential decay of the disturbance fields for any non-zero 
$\beta _\infty$ that, at sufficiently large distances, would have led to the downstream influence of the sphere being smaller than its upstream influence; a feature that defies intuition (at least one based on homogeneous fluids at finite Reynolds numbers). However, the aforementioned tertiary screening mechanism intervenes, so to speak, at the right distance, the result being that the aforementioned exponential decay transitions to an algebraic one at larger distances. Further, the rate of algebraic decay is the same as that upstream, albeit with a larger numerical prefactor (
$3240$ as opposed to 
$2160$). Thus, the resulting downstream influence is still larger than that upstream, although the difference between the two is only a constant of order unity (the ratio 
$3240/2160 = 1.5$). A more profound upstream–downstream asymmetry will have to await consideration of inertial forces.
The aforementioned results set up the groundwork for a detailed drift-volume calculation, the significance of which was highlighted in the introduction. The more rapid decay of the velocity field induced by stratification implies that the drift volume is finite. VS22 argued that the drift volume must be 
${O}(a^3Ri_v^{-2/3})$ for 
$Ri_v \rightarrow 0$ in the Stokes-stratification regime, although determining the numerical prefactor requires a detailed calculation. The calculation methodology is depicted via the schematic in figure 10(a) (see also Varanasi Reference Varanasi2022) where, instead of initializing the material plane far upstream of the sphere, we choose it to be the plane 
$z = 0$ containing the sphere's centre. The pathlines of different material points are then obtained by integrating backward (over the interval 
$[0,-T]$) and forward (over the interval 
$[0,T]$) in time, until a time 
$T$ large enough for the displacement to converge to a finite value. For a sufficiently dense ensemble of initial material points, one may use the drift displacements to determine the volume enclosed between the initial and final material planes for the forward (backward) integration, which is termed the ‘downstream (upstream) drift volume’. The total drift volume is the sum of upstream and downstream drift components. Figures 10(c) and 10(d) showcase the contrasting nature of drift displacements, as a function of time, for two points whose initial locations are identified in figure 10(b): (a) point 1 is a material element initially at 
$(\tilde {r}_t,z) \equiv (0.13,0)$, and that traverses the downstream columnar structure on its way past the sphere; (b) point 2 is a material element at 
$(\tilde {r}_t,0) = (10,0)$ which only encounters the wake region. In the backward integration, both material points encounter only the rapidly decaying flow in the upstream wake, leading to a quick saturation of the drift displacements (the blue curves in figure 10c,d). In the forward integration, point 2 passes through a longer sequence of recirculating cells in the downstream wake, leading to an oscillatory response of the drift displacement prior to saturation (the red curve in figure 10d). Point 1, when traversing the columnar structure, experiences the reverse-Stokeslet-like (
$1/\tilde {z}$) decay close to the rear stagnation streamline (see figure 5), which leads to a logarithmic increase in the drift displacement (the red curve in figure 10c). This increase results in an eventual change in sign of the drift displacement (relative to the homogeneous case), leading to a contribution that has the character of a reflux. Our preliminary results for the drift volume show that the upstream and downstream components conform to the VS22 estimate above, but the total drift volume is smaller owing to cancellation induced by the reflux nature of the downstream component. A detailed study of the drift displacement characteristics, and the drift-volume calculation, will be presented in a separate communication.
Figure 10. (a) Schematic of the drift-volume calculation methodology, illustrating the total drift-volume evaluation as the sum of upstream and downstream components. (b) The initial locations of two material points (black stars) superimposed on the streamline pattern for 
$\beta _\infty = 0$; the red lines indicate the disturbance flow field encountered by these points during the sphere's translation. The non-dimensional drift displacements (
$z/a$) of (c) point 1 and (d) point 2, plotted as a function of time (
$tU/a$), for 
$Ri_v = 10^{-4}$. The red curve is the downstream drift displacement resulting from the forward integration, and the blue curve denotes the upstream displacement resulting from the backward integration; the corresponding dashed curves are the drift displacements in the absence of stratification. Black dashed lines mark the time when the points cross the stratification screening length 
$l_c$.
It is worth ending our investigation with figures 11(a) and 11(b) – the first one containing streamline patterns and the second one containing isopycnals – that help contrast the structure of the disturbance fields in the diffusion-dominant and convection-dominant limits. This contrast is achieved by having each of these figures consist of two halves, one corresponding to 
$\beta _\infty \gg 1$, and the other to 
$\beta _\infty \ll 1$. It should be noted that the diffusion-dominant halves in the said figures differ from the versions presented earlier (figures 6a and 8a) because, in the interest of a fair comparison, the length scale used in figures 11(a) and 11(b) is now the small-
$Pe$ screening length of 
${O}[a(Ri_vPe)^{-1/4}]$; the ratio of this length scale to 
$aRi_v^{-1/3}$ is 
$\beta _\infty ^{1/4}$, and the large-
$\beta _\infty$ halves in figures 11(a) and 11(b) are rescaled versions of figures 6(a) and 8(a), obtained by multiplying the axes with 
$\beta _\infty ^{1/4}$.
Figure 11. (a) Streamline and (b) isopycnal patterns in an axial (half) plane for 
$\beta _\infty = 10^{-5}$ (for 
$-14<\tilde {r}_t<0$) and 
$10$ (for 
$0<\bar {r}_t<14$); only the downstream region, corresponding to 
$\tilde {z} > 1$, is shown. Please note that the radial and axial coordinates for 
$\beta _\infty = 10^{-5}$ are scaled by the large-
$Pe$ stratification screening length of 
${O}[a(Ri_v)^{-1/3}]$, whereas those for 
$\beta _\infty = 10$ are scaled by the small-
$Pe$ screening length of 
${O}[a(Ri_vPe)^{-1/4}]$.
Acknowledgements
We acknowledge the use of the computing resources at HPCE, IIT Madras. R.P. acknowledges fruitful discussions with A. Varanasi during the early stages of this work. R.P. and A.R. acknowledge support from IIT Madras for its support of the ‘Geophysical Flows Lab’ research initiative under the Institute of Eminence framework.
Declaration of interests
The authors report no conflicts of interest.
Appendix A. Simplification of the disturbance flow and density fields
A.1. The Stokes streamfunction
The expression for the Stokes streamfunction, obtained via an inverse Fourier transform, is as given in (2.11). For the sake of continuity, it is repeated here
\begin{equation} \tilde{\psi}_s(\tilde{\boldsymbol{r}}) = \dfrac{3\tilde{r}_t i}{4{\rm \pi}^2}\displaystyle\int\dfrac{({\rm i}k_3+\beta_\infty k^2)k_2}{({\rm i}k_3+\beta_\infty k^2)k^4+k_t^2}{\rm e}^{{\rm i}\boldsymbol{k}\boldsymbol{\cdot}\tilde{\boldsymbol{r}}}\,{\rm d}\boldsymbol{k}. \end{equation}For ease of numerical evaluation, the Stokeslet contribution can be separated out first in a manner similar to (3.4) of VS22, in which case (A1) takes the form
\begin{equation} \tilde{\psi}_s(\tilde{\boldsymbol{r}}) = \dfrac{3\tilde{r}_t i}{4{\rm \pi}^2}\displaystyle\int\left(\dfrac{k_2}{k^4} -\dfrac{k_t^2k_2}{({\rm i}k_3+\beta_\infty k^2)k^4+k_t^2} \right) {\rm e}^{{\rm i}\boldsymbol{k}\boldsymbol{\cdot} \tilde{\boldsymbol{r}}}\,{\rm d}\boldsymbol{k}. \end{equation}
Here, the first term within brackets corresponds to the Stokeslet contribution and can be evaluated analytically. The three-dimensional integral in the second term needs to be evaluated numerically; however, it can first be reduced to a two-dimensional integral by a series of simplifications. Using spherical polar coordinates (
$k,\theta,\phi$) in Fourier space for the second term, with the polar axis along 
$\boldsymbol {1}_3$, (A2) can be written as
\begin{equation} \tilde{\psi}_s(\tilde{r}_t,\tilde{z}) = \displaystyle -\dfrac{3\tilde{r}_t^2}{4\sqrt{\tilde{r}_t^2+\tilde{z}^2}} -\dfrac{3\tilde{r}_t i}{4{\rm \pi}^2} \int_0^{\rm \pi}{\rm d}\theta\int_{-{\rm \pi}}^{\rm \pi}{\rm d}\phi \int_0^{\infty}{\rm d}k \dfrac{\sin^4\theta\sin\phi\exp({\rm i}k\bar{\delta}_1)}{k\left({\rm i}k^3 \cos\theta + \beta_{\infty} k^4 + \sin^2\theta\right)}, \end{equation}
where 
$k_t = \sqrt {k_1^2+k_2^2} = k\sin \theta$, 
$k_3 = k\cos \theta$, 
$k_2 = k\sin \theta \sin \phi$ and 
$\bar {\delta }_1 = \tilde {r}_t\sin \theta \sin \phi + \tilde {z}\cos \theta$. Here, the first term is the Stokeslet contribution in physical space. The second term can be simplified further based on angular symmetry, resulting in
\begin{align} \tilde{\psi}_s
&={-}\displaystyle
\dfrac{3\tilde{r}_t^2}{4\sqrt{\tilde{r}_t^2+\tilde{z}^2}} -
\dfrac{3\tilde{r}_t i}{2{\rm \pi}^2} \int_0^{{\rm \pi}/2}{\rm
d}\theta\int_{0}^{{\rm \pi}/2}{\rm d}\phi \nonumber\\ &\quad\times \int_{-\infty}^{\infty}{\rm d}k
\dfrac{\sin^4\theta\sin\phi\left(\exp({\rm
i}k\bar{\delta}_1)-\exp(-{\rm
i}k\delta_1)\right)}{k\left({\rm i}k^3 \cos\theta +
\beta_{\infty} k^4 + \sin^2\theta\right)},
\end{align}
with 
$\delta _1 = \tilde {r}_t\sin \theta \sin \phi - \tilde {z}\cos \theta$.
As mentioned in § 2, we first evaluate the 
$k$-integral in (A4) using contour integration. However, the contour should be carefully chosen considering the signs of 
$\delta _1$ and 
$\bar {\delta }_1$ in the arguments of the exponentials, so as to satisfy Jordan's lemma. Here, we present the steps for positive 
$\tilde {z}$; the expressions for negative 
$\tilde {z}$ can be derived similarly. For positive 
$\tilde {z}$, 
$\bar {\delta }_1$ in (A4) is always positive, whereas 
$\delta _1$ can be either positive for 
$\theta < \theta _c$ or negative for 
$\theta > \theta _c$; 
$\theta _c({=}\arctan (\tilde {z}/\tilde {r}_t\sin \phi ))$ being the critical polar angle at which 
$\delta _1 = 0$. So, for the term involving the first exponential and for the term involving the second exponential when 
$\delta _1<0$, a contour in the upper half of the complex
$-k$ plane, as shown in figure 12(a) (‘case I’), should be chosen. When 
$\delta _1>0$, the contour for the term involving the second exponential should be chosen in the lower half of the complex
$-k$ plane, as shown in figure 12(b) (‘case II’). Therefore, one can write
where the left-hand side denotes the integral corresponding to case I in (A4). The integral along the curve 
$C_{R2}$ has been set to zero due to Jordan's lemma (see Lemma 4.2.2 in Ablowitz & Fokas Reference Ablowitz and Fokas2003). For case II, 
$C_{L1}, C_{L2}, C_{R1}, C_{R2}$ and I in (
) should be replaced by 
$C_{L3}, C_{L4}, C_{R3}, C_{R4}$ and II, respectively. Also, note that the signs of the contour integrals (
$\oint _I$ or 
$\oint _{II}$), and the integrals along the curves 
$C_{R1}$ and 
$C_{R3}$, depend on sense (clockwise/anticlockwise) in which the contours/curves are traversed; see figures 12(a) and 12(b). The contribution from the curves 
$C_{R1}$ and 
$C_{R3}$ can be evaluated (using theorem 4.3.1 in Ablowitz & Fokas Reference Ablowitz and Fokas2003) to be 
$-{\rm i}{\rm \pi} \sin ^2\theta \sin \phi$ and 
${\rm i}{\rm \pi} \sin ^2\theta \sin \phi$, respectively.
Figure 12. Contours of integration used in solving the 
$k$-integral in (A4). Here, (a) shows a contour in the upper half of the complex
$-k$ plane, used when the argument of the complex exponentials in (A4) is positive (case I), whereas (b) shows a contour in the lower half of the complex
$-k$ plane used when the argument of the complex exponentials in (A4) is negative (case II). The parameters 
$C_{L1}-C_{L4}$ are integrals along the real line, 
$R$ is the radius of the curves 
$C_{R2}$ and 
$C_{R4}$ and 
$\epsilon$ is the radius of curves 
$C_{R1}$ and 
$C_{R3}$. Here, 
$\kappa _0$ to 
$\kappa _4$ are zeros of the denominator in the integrand of (A4).
To evaluate the contour integral that remains, the residues at the simple poles of the integrand need to be calculated (see (4.1.10) in Ablowitz & Fokas Reference Ablowitz and Fokas2003). These correspond to the zeros of the denominator of the integrand excluding the origin, and therefore, to the four roots of the quartic polynomial 
$\beta _\infty k^4 + {\rm i}k^3\cos \theta +\sin ^2\theta$. Figure 13 shows the behaviour of the roots in the complex
$-k$ plane, as a function of 
$\theta$ (
$0\leqslant \theta \leqslant {\rm \pi}/2$), for various 
$\beta _\infty$. Two roots exist in the upper half of the complex
$-k$ plane and are located symmetrically on either side of the positive imaginary axis. The remaining two roots exist in the lower half of the complex
$-k$ plane; they can be purely imaginary or can lie symmetrically on either side of the negative imaginary axis, depending on 
$\theta$; note that one of the roots starts at 
$-\beta _\infty ^{-1} i$ at 
$\theta = 0$, as shown in figure 13(b). We number the roots 
$\kappa _1-\kappa _4$ in an anticlockwise sense, so roots 
$\kappa _1$ and 
$\kappa _2$ lie in the upper half, while 
$\kappa _3$ and 
$\kappa _4$ lie in the lower half of the complex
$-k$ plane.
Figure 13. The behaviour of (a) first and second roots, and (b) third and fourth roots of the quartic polynomial (
$\beta _\infty k^4+{\rm i}k^3\cos \theta + \sin ^2\theta$), as a function of 
$\theta$ and for various 
$\beta _\infty$, in the complex
$-k$ plane. Here, 
$\theta$ varies from 
$0$ to 
${\rm \pi} /2$ and increases in the direction shown by the arrows. The continuous and dashed lines correspond to root 1 (3) and 2 (4) in panel (a) (panel b), respectively.
After contour integration, the integral involving the fist exponential, in (A4), can be written as
\begin{equation} \displaystyle \dfrac{3\tilde{r}_t}{2 {\rm \pi}} \int_0^{{\rm \pi}/2}{\rm d}\theta\int_{0}^{{\rm \pi}/2}{\rm d}\phi \left[2 \sum_{m=1}^{2} \dfrac{\sin^4\theta\sin\phi\exp({\rm i}\kappa_m\bar{\delta}_1)}{\kappa_m^3(3{\rm i}\cos\theta+4\beta_\infty \kappa_m)} + \sin^2\theta\sin\phi\right], \end{equation}
where the first term in the brackets is due to the residues associated with 
$\kappa _1$ and 
$\kappa _2$, and the second term is due to the integral along 
$C_{R1}$. The integral in (A4), involving the second exponential, can be written as
\begin{align} &-\dfrac{3\tilde{r}_t}{2{\rm \pi}} \int_{0}^{{\rm \pi}/2}{\rm d}\phi \left(\int_0^{\theta_c}{\rm d}\theta \sum_{m=1}^{2} - \int_{\theta_c}^{{\rm \pi}/2}{\rm d}\theta \sum_{m=3}^{4} \right)\dfrac{2\sin^4\theta\sin\phi\exp(-{\rm i} \kappa_m\delta_1)}{\kappa_m^3(3{\rm i}\cos\theta+4\beta_\infty \kappa_m)} \nonumber\\ &\quad -\dfrac{3\tilde{r}_t}{2{\rm \pi}} \int_{0}^{{\rm \pi}/2}{\rm d}\phi\left(\int_0^{\theta_c}{\rm d}\theta -\int_{\theta_c}^{{\rm \pi}/2}{\rm d}\theta \right) \sin^2\theta\sin\phi, \end{align}
where the first term is due to the roots of the quartic polynomial, and the second term is due to the integral along 
$C_{R1}$ and 
$C_{R3}$. The integrals corresponding to 
$C_{R1}$ and 
$C_{R3}$ have been given above, and taken together, will cancel the Stokeslet contribution in (A4). Finally, one can express the Stokes streamfunction as
\begin{equation} \tilde{\psi}_s = \dfrac{3\tilde{r}_t}{2 {\rm \pi}} (T_{\psi 1}+T_{\psi 2}), \end{equation}
where, for 
$\tilde {z}>0$,
\begin{equation} \left.\begin{gathered}
T_{\psi 1}(\tilde{r}_t, \tilde{z}) = \int_{0}^{{\rm \pi}/2}{\rm
d}\phi \left[\int_0^{{\rm \pi}/2}{\rm d}\theta
\sum_{m=1}^{2}\dfrac{2\sin^4\theta\sin\phi\exp({\rm
i}\kappa_m\bar{\delta}_1)}{\kappa_m^3(3{\rm
i}\cos\theta+4\beta_\infty \kappa_m)}\right]\\ T_{\psi
2}(\tilde{r}_t, \tilde{z}) ={-} \int_{0}^{{\rm \pi}/2}{\rm
d}\phi \left[\left(\int_0^{\theta_c}{\rm d}\theta
\sum_{m=1}^{2} \!-\! \int_{\theta_c}^{{\rm \pi}/2}{\rm d}\theta
\sum_{m=3}^{4}
\right)\dfrac{2\sin^4\theta\sin\phi\exp(-{\rm
i}\kappa_m\delta_1)}{\kappa_m^3(3{\rm
i}\cos\theta+4\beta_\infty \kappa_m)} \right].
\end{gathered}\right\}
\end{equation}
Using similar arguments as above, the expressions for 
$\tilde {z}<0$ can also be written as
\begin{equation} \left.\begin{gathered} T_{\psi 1}(\tilde{r}_t, \tilde{z}) = \int_{0}^{{\rm \pi}/2}{\rm d}\phi \left[ \int_0^{{\rm \pi}/2}{\rm d}\theta \sum_{m=3}^{4}\dfrac{2\sin^4\theta\sin\phi\exp(-{\rm i} \kappa_m\delta_1)}{\kappa_m^3(3{\rm i}\cos\theta+4\beta_\infty \kappa_m)} \right] \\ T_{\psi 2}(\tilde{r}_t, \tilde{z}) ={-}\int_{0}^{{\rm \pi}/2}{\rm d}\phi \left[\left(\int_0^{\theta_c}{\rm d}\theta \sum_{m=3}^{4} - \int_{\theta_c}^{{\rm \pi}/2}{\rm d}\theta \sum_{m=1}^{2} \right)\dfrac{2\sin^4\theta\sin\phi\exp({\rm i}\kappa_m\bar{\delta}_1)}{\kappa_m^3(3{\rm i}\cos\theta+4\beta_\infty \kappa_m)}\right]. \end{gathered}\right\} \end{equation}
The 
$\phi$-integral in the above expressions can be evaluated analytically only for those terms with the 
$\theta$-integral limits independent of 
$\phi$ i.e. for 
$T_{\psi 1}$ in (A9) and in (A10), leading in these cases to a single 
$\theta -$integral in terms of Bessel and Struve functions as given in (A11). Finally, by using the variables 
$(x,y) = (\sin \phi,\sin \theta )$, expressions (A9)–(A10) can be written as
\begin{equation} \left.\begin{gathered}
T_{\psi 1}(\tilde{r}_t, \tilde{z}) =
{\rm \pi}\displaystyle\int_0^1{\rm
d}y(1\!-\!y^2)^{{3}/{2}}\sum_{m=p}^{p+1}\dfrac{\left({\rm i} s
{\rm J}_1(\tilde{r}_t
\kappa_m\sqrt{1\!-\!y^2})+{\rm H}_{{-}1}(\tilde{r}_t
\kappa_m\sqrt{1\!-\!y^2})\right){\rm e}^{{\rm i}
\kappa_m\tilde{z}y}}{\kappa_m^3(3{\rm i}y+4\beta_\infty
\kappa_m)},\\ T_{\psi 2}(\tilde{r}_t, \tilde{z}) =
2\displaystyle\int_0^1{\rm d}\kern0.7pt x\left[\int_0^{y_c}{{\rm d}
y}\sum_{m=q}^{q+1}-\int_{y_c}^1{\rm
d}y\sum_{m=p}^{p+1}\right]\dfrac{(1-y^2)^{{3}/{2}}{\rm
e}^{-{\rm i}\kappa_m\delta}}{\kappa_m^3(3{\rm
i}y+4\beta_\infty \kappa_m)},
\end{gathered}\right\}
\end{equation}
with 
$\delta = s \tilde {r}_t \sqrt {1-x^2} \sqrt {1-y^2} - \tilde {z} y$ and 
$y_c = \left.\tilde {r}_t\sqrt {1-x^2}\right /\sqrt {\tilde {z}^2+\tilde {r}_t^2(1-x^2)}$. Here, 
$\textrm {J}_1$ and 
$H_{-1}$ in (A11) are, respectively, the Bessel and Struve functions of the first kind and of orders 
$1$ and 
$-1$ (Abramowitz & Stegun Reference Abramowitz and Stegun1968). The indices, that appear in the summations in (A11), are defined as 
$(p,q,s) \equiv (1,3,1)$ for 
$\tilde {z} > 0$ and 
$(p,q,s) \equiv (3,1,-1)$ for 
$\tilde {z}<0$; the indices number the roots in terms of their position in the complex
$-k$ plane, taken in an anticlockwise manner starting from the first quadrant.
The limit 
$\beta _\infty \to 0$ is a singular one when the quartic polynomial in 
$k$ degenerates to a cubic one, with one of the aforementioned four roots receding to infinity. As mentioned earlier, root 
$\kappa _4(\theta = 0) = -\beta _\infty ^{-1} i$. Therefore, for 
$\beta _\infty \to 0$, 
$\kappa _4$ approaches 
$-\infty i$, making the evaluation of the residue at 
$\kappa _4$ difficult. So, we set 
$\beta _\infty = 0$ in (A4), reducing the order of the polynomial in the denominator by one. The resulting expression can be simplified as
\begin{align} \tilde{\psi}_s &={-}\dfrac{3\tilde{r}_t^2}{4\sqrt{\tilde{r}_t^2+\tilde{z}^2}} + \dfrac{3\tilde{r}_t}{{\rm \pi}^2} \int_0^{{\rm \pi}/2}{\rm d}\theta\int_{0}^{{\rm \pi}/2}{\rm d}\phi \int_{0}^{\infty}{\rm d}k \nonumber\\ &\quad \times \dfrac{\sin^4\theta\sin\phi\left(\sin^2\theta(\sin k\bar{\delta}_1 + \sin k \delta_1) - k^3\cos\theta(\cos k\bar{\delta}_1 - \cos k\delta_1)\right)}{k(k^6\cos^2\theta + \sin^4\theta)}. \end{align}
In evaluating the 
$k$-integral for the terms including 
$\sin k\bar {\delta }_1$ and 
$\sin k\delta _1$, we use identity (3.738.1) from Gradshteyn & Ryzhik (Reference Gradshteyn and Ryzhik2007). For the terms including 
$\cos k\bar {\delta }_1$ and 
$\cos k\delta _1$, we use identity (3.738.2). Please note that the said identities are applicable only when the arguments of 
$\sin k\bar {\delta }_1$ (or 
$\sin k\delta _1$) and 
$\cos k\bar {\delta }_1$ (or 
$\cos k\delta _1$) are positive. Therefore, when 
$\delta _1$ (
$\bar {\delta }_1$) is negative, for positive (negative) 
$\tilde {z}$, it should be replaced with 
$-|\delta _1|$ (
$-|\bar {\delta }_1|$). So one obtains, for example, 
$\sin k\delta _1 = \sin (-k|\delta _1|) = -\sin k|\delta _1|$ and 
$\cos k\delta _1 = \cos (-k|\delta _1|) = \cos k|\delta _1|$. Subsequent simplifications lead to the following alternate expressions for 
$T_{\psi 1}$, and 
$T_{\psi 2}$ for 
$\beta _\infty = 0$:
\begin{equation} \left.\begin{gathered}
T_{\psi 1}(\tilde{r}_t, \tilde{z}) =
\dfrac{1}{3}\displaystyle \int_0^1{\rm d}\kern0.7pt
x\int_0^1{\rm d}y\sqrt{1-y^2}\left(3-\varPsi_1
(\bar{\delta})\varTheta(\tilde{z})-\varPsi_2(\delta)\varTheta(-\tilde{z})\right),\\
T_{\psi 2}(\tilde{r}_t, \tilde{z}) =
\dfrac{1}{3}\displaystyle \int_0^1{\rm d}\kern0.7pt x
\left[\vphantom{\left.-\displaystyle\int_{y_c}^1{\rm
d}y\sqrt{1-y^2}\left(3-\varPsi_1(\delta)\varTheta(\tilde{z})-\varPsi_2(|\bar{\delta}|)
\varTheta(-\tilde{z})\right)\right]}\int_0^{y_c}{{\rm d}
y}\sqrt{1-y^2}\left(3-\varPsi_1(\bar{\delta})
\varTheta(-\tilde{z})-\varPsi_2(\delta)\varTheta(\tilde{z})\right)\right.\\
\quad\quad \left.-\displaystyle\int_{y_c}^1{\rm
d}y\sqrt{1-y^2}\left(3-\varPsi_1(\delta)\varTheta(\tilde{z})-\varPsi_2(|\bar{\delta}|)
\varTheta(-\tilde{z})\right)\right],
\end{gathered}\right\}
\end{equation}where
\begin{equation} \left.\begin{gathered} \varPsi_1(\varDelta) = 4\exp\left(-\frac{\varDelta}{2}\sqrt[3]{\frac{1-y^2}{y^2}}\right)\cos\left(\frac{\varDelta \sqrt{3}}{2}\sqrt[3]{\frac{1-y^2}{y^2}}\right), \\ \varPsi_2(\varDelta) = 2\exp\left(-\varDelta\sqrt[3]{\frac{1-y^2}{y^2}}\right)\\ \bar{\delta} = s \tilde{r}_t \sqrt{1-x^2} \sqrt{1-y^2} + \tilde{z} y, \quad {\delta} = s \tilde{r}_t \sqrt{1-x^2} \sqrt{1-y^2} - \tilde{z} y, \end{gathered}\right\} \end{equation}
and 
$\varTheta (\tilde {z})$ is the Heaviside function. Other variables, including 
$s$ in particular, have the same definitions as in (A11).
A.2. The density disturbance field
For arbitrary 
$\beta _\infty$, the inverse Fourier transform integral for the density perturbation is given in (2.10)
\begin{equation} \tilde{\rho}_f(\tilde{\boldsymbol{r}}) = \dfrac{-3}{4{\rm \pi}^2}\displaystyle\int\dfrac{k_t^2}{({\rm i}k_3+\beta_\infty k^2)k^4+k_t^2}{\rm e}^{{\rm i}\boldsymbol{k}\boldsymbol{\cdot}\tilde{\boldsymbol{r}}}\,{\rm d}\boldsymbol{k}. \end{equation}Again, use of spherical polar coordinates in Fourier space leads to
\begin{align} \tilde{\rho}_f &={-} \displaystyle \dfrac{3}{2{\rm \pi}^2} \int_0^{{\rm \pi}/2}{\rm d}\theta\int_{0}^{{\rm \pi}/2}{\rm d}\phi \int_{-\infty}^{\infty}{\rm d}k \dfrac{\sin^3\theta k^2\left(\exp({\rm i}k\bar{\delta}_1)+\exp(-{\rm i}k\delta_1)\right)}{\left({\rm i}k^3 \cos\theta + \beta_{\infty} k^4 + \sin^2\theta\right)} \nonumber\\ &={-}\dfrac{3}{2{\rm \pi}}(T_{\rho1} + T_{\rho2}). \end{align}
The 
$k$-integral here can be performed using contour integration similar to the earlier subsection. An important difference is that the contour need not be indented to exclude the origin, i.e. the small semi-circular curve (
$C_{R1}$ and 
$C_{R3}$), as there is no pole at 
$k=0$. Therefore, the integration contour is a semi-circle either in the upper quadrants or the lower quadrants, depending on the sign of 
$\delta _1$ and 
$\bar {\delta }_1$. For 
$\tilde {z}>0$, after the 
$k-$integration, one obtains
\begin{gather} T_{\rho1}(\tilde{r}_t, \tilde{z}) = 2{\rm i} \displaystyle \int_0^{{\rm \pi}/2}{\rm d}\theta\int_{0}^{{\rm \pi}/2}{\rm d}\phi \sum_{m=1}^{2}\dfrac{\sin^3\theta \exp({\rm i}\kappa_m\bar{\delta}_1)}{\left(3{\rm i}\cos\theta+4\beta_\infty \kappa_m\right)}, \end{gather}
\begin{gather}T_{\rho2}(\tilde{r}_t, \tilde{z}) = 2{\rm i} \displaystyle \int_{0}^{{\rm \pi}/2}{\rm d}\phi \left[\int_0^{\theta_c}{\rm d}\theta \sum_{m=1}^{2}-\int_{\theta_c}^{{\rm \pi}/2}{\rm d}\theta \sum_{m=3}^{4}\right]\dfrac{\sin^3\theta \exp(-{\rm i}\kappa_m\delta_1)}{\left(3{\rm i}\cos\theta+4\beta_\infty \kappa_m\right)}. \end{gather}
For 
$\tilde {z}<0$,
\begin{gather} T_{\rho1}(\tilde{r}_t,\tilde{z}) ={-}2{\rm i}\displaystyle \int_0^{{\rm \pi}/2}{\rm d}\theta\int_{0}^{{\rm \pi}/2}{\rm d} \phi \sum_{m=3}^{4}\dfrac{\sin^3\theta \exp(-{\rm i}\kappa_m\delta_1)}{\left(3{\rm i}\cos\theta+4\beta_\infty \kappa_m\right)}, \end{gather}
\begin{gather}T_{\rho2}(\tilde{r}_t,\tilde{z}) ={-}2{\rm i} \displaystyle \int_{0}^{{\rm \pi}/2}{\rm d}\phi \left[\int_0^{\theta_c} {\rm d}\theta \sum_{m=3}^{4}-\int_{\theta_c}^{{\rm \pi}/2}{\rm d}\theta \sum_{m=1}^{2}\right] \dfrac{\sin^3\theta \exp({\rm i}\kappa_m\bar{\delta}_1)}{\left(3{\rm i}\cos\theta+4\beta_\infty, \kappa_m\right)}. \end{gather}
Similar to the case for the Stokes streamfunction, the 
$\phi$-integral in (A17) and (A19) can be evaluated in terms of Bessel and Struve functions. Finally, with a change of variables from 
$(\theta,\phi )$ to 
$(x,y) = (\sin \theta,\sin \phi )$, one obtains
\begin{equation} \left.\begin{gathered}
T_{\rho 1}(\tilde{r}_t, \tilde{z}) = {\rm \pi}{\rm i}
\displaystyle\int_0^1{\rm
d}y(1\!-\!y^2)\sum_{m=p}^{p+1}\left[\dfrac{\left[s {\rm
J}_0(\kappa_m\tilde{r}_t\sqrt{1\!-\!y^2})+ {\rm
i}{\rm H}_0(\kappa_m\tilde{r}_t\sqrt{1\!-\!y^2})\right]{\rm e}^{{\rm
i}\kappa_m\tilde{z}y}}{3{\rm i}y+ 4\beta
\kappa_m}\right],\\ T_{\rho 2}(\tilde{r}_t, \tilde{z}) = 2
{\rm i} \displaystyle\int_0^1\dfrac{{\rm d}\kern0.7pt x}{\sqrt{1-x^2}}\left[-\int_0^{y_c}{{\rm
d}y}\sum_{m=q}^{q+1}+\int_{y_c}^1{\rm
d}y\sum_{m=p}^{p+1}\right]\dfrac{s(1-y^2){\rm e}^{-{\rm
i}\kappa_m\delta}}{3{\rm i}y+ 4\beta \kappa_m}.
\end{gathered}\right\}
\end{equation}
Here, 
$\textrm {J}_0$ and 
${\rm H}_{0}$ are, respectively, the Bessel and Struve functions of the first kind and order 
$0$ (Abramowitz & Stegun Reference Abramowitz and Stegun1968). The definitions for different variables (
$\delta, y_c, \kappa _m$) and indices (
$\kern0.7pt p,q,s$) are as given in (A21). For 
$\beta _\infty = 0$, as mentioned in the previous subsection, we set 
$\beta _\infty$ to zero in (A16) and write
\begin{align} \tilde{\rho}_f &={-} \dfrac{3}{{\rm \pi}^2} \int_0^{{\rm \pi}/2}{\rm d}\theta\int_{0}^{{\rm \pi}/2}{\rm d}\phi \int_{0}^{\infty}{\rm d}k\, k^2\sin^3\theta\nonumber\\ &\quad \times\dfrac{\sin^2\theta(\cos k\bar{\delta}_1 + \cos k\delta_1) + k^3 \cos \theta(\sin k\bar{\delta}_1 - \sin k\delta_1)}{k^6\cos^2\theta + \sin^4\theta}. \end{align}
The 
$k$-integral is then evaluated using identities (3.738.1) and (3.738.2) from Gradshteyn & Ryzhik (Reference Gradshteyn and Ryzhik2007), as done in the previous subsection. Finally, one can write
\begin{equation} \left.\begin{gathered} T_{\rho 1}(\tilde{r}_t, \tilde{z}) = \dfrac{1}{3} \int_0^1{\rm d}\kern0.7pt x\int_0^1{\rm d}y\dfrac{(1-y^2)}{y\sqrt{1-x^2}}\left(\varPsi_1(\bar{\delta})\varTheta(\tilde{z})- \varPsi_2(\delta)\varTheta(-\tilde{z})\right),\\ T_{\rho 2}(\tilde{r}_t, \tilde{z}) = \dfrac{1}{3}\int_0^1{\rm d}\kern0.7pt x \left[\int_0^{y_c}{{\rm d} y}\dfrac{(1-y^2)}{y\sqrt{1-x^2}}\left(\varPsi_1(\bar{\delta})\varTheta(-\tilde{z})-\varPsi_2(\delta) \varTheta(\tilde{z})\right)\right.\\ \qquad\qquad \left.+\int_{y_c}^1{\rm d}y\dfrac{(1-y^2)}{y\sqrt{1-x^2}}\left(\varPsi_1(\delta)\varTheta(\tilde{z})- \varPsi_2(|\bar{\delta}|)\varTheta(-\tilde{z})\right)\right], \end{gathered}\right\} \end{equation}
where the definitions for variables 
$\varPsi _1, \varPsi _2, \varTheta, \delta$ and 
$\bar {\delta }$ are as described in (A13).
Appendix B. Calculating the boundary of the downstream columnar structure
In order to find an expression for the boundary of the columnar structure, we first derive a closed-form expression for the Stokes streamfunction under the jet approximation, that is, the limit 
$k_3 \ll k_t$, corresponding to nearly vertical flow. In this limit, (A1) simplifies to
\begin{equation} \tilde{\psi}_s = \dfrac{3\tilde{r}_t i}{4{\rm \pi}^2}\int \dfrac{({\rm i}k_3+\beta_\infty k_t^2)k_2}{({\rm i}k_3 + \beta_\infty k_t^2)k_t^4 + k_t^2}{\rm e}^{{\rm i}\boldsymbol{k}\boldsymbol{\cdot}\tilde{\boldsymbol{r}}}\,{\rm d}\boldsymbol{k}. \end{equation}
Using cylindrical coordinates in Fourier space – (
$k_t,\phi,k_3$) facilitates the evaluation of the 
$k_3$-integral using contour integration (VS22). Subsequently, performing the 
$\phi$-integral results in the following 1-D integral in terms of 
$k_t$:
\begin{equation} \tilde{\psi}_s = 3\tilde{r}_t\int_0^\infty {\rm d}k_t \dfrac{{\rm J}_1(\tilde{r}_t k_t)\exp({-z(\beta_\infty k_t^2 + 1/k_t^2)})}{k_t^4}. \end{equation}
We now evaluate this integral in closed form in the non-diffusive limit (
$\beta _\infty = 0$).
B.1. Steepest-descent method for the far-field jet approximation integral
For 
$\beta _\infty = 0$, using a change of variables 
$\lambda = \sqrt {z}/k_t$, the integral in (B2) can be written as
\begin{equation} \frac{3\tilde{r}_t}{\tilde{z}^{3/2}}\int_0^\infty {\rm d}\lambda\,\lambda^2 {\rm J}_1\left(\frac{\tilde{r}_t\tilde{z}^{1/2}}{\lambda}\right)\exp(-\lambda^2). \end{equation}
In the limit 
$\tilde {r}_t\tilde {z}^{1/2} \gg 1$, the dominant contribution to the integral comes from 
$\lambda \sim {O}(1)$, owing to the decaying exponential in the integrand. Therefore, 
$(\tilde {r}_t\tilde {z}^{1/2})/\lambda \gg 1$, and one can approximate the Bessel function in terms of its large-argument trigonometric representation 
$\textrm {J}_1(x) = (\sin x-\cos x)/\sqrt {{\rm \pi} x}$, whence (B3) simplifies to
\begin{equation} \frac{3(1+{\rm i})(\tilde{r}_t\tilde{z}^{1/2})^{1/2}}{2\sqrt{\rm \pi}\tilde{z}^{2}}\int_{-\infty}^\infty {\rm d}\lambda\,\lambda^{5/2} \exp\left(-\lambda^2+{\rm i}\dfrac{\tilde{r}_t\tilde{z}^{1/2}}{\lambda}\right). \end{equation}
A further change in variables, 
$\hat {\lambda } = \lambda (\tilde {r}_t\tilde {z}^{1/2})^{-1/3}$, leads to
\begin{equation} \frac{3(1+{\rm i})\tilde{r}_t\tilde{z}^{1/2}}{2\sqrt{\rm \pi}\tilde{z}^{2}}\int_{-\infty}^\infty {\rm d}\hat{\lambda}\,\hat{\lambda}^{5/2} \exp\left((\tilde{r}_t\tilde{z}^{1/2})^{2/3}g(\hat{\lambda})\right),\quad \textrm{where } g(\hat{\lambda}) = \left(-\hat{\lambda}^2+\dfrac{{\rm i}}{\hat{\lambda}}\right). \end{equation}
This integral can be determined using the steepest-descent method, with the large parameter being 
$(\tilde {r}_t\tilde {z}^{1/2})^{2/3}$. Accordingly, the integration path along the real axis is deformed onto the appropriate constant phase contour (
$\textrm {Im}(g(\hat {\lambda }))$ is a constant) in the complex-
$\hat {\lambda }$ plane. Next, the saddle points on this steepest-descent contour are identified. These points correspond to 
$\textrm {Re}(g(\hat {\lambda }))$ attaining a maximum, and for 
$\tilde {r}_t\tilde {z}^{1/2} \gg 1$, the dominant contribution to the integral arises from the neighbourhood of the saddle point (Bender & Orszag Reference Bender and Orszag2013). Substituting 
$\hat {\lambda } = \hat {\lambda }_r + \textrm {i}\hat {\lambda }_i$ in 
$g(\hat {\lambda })$ above and writing the real and imaginary parts separately yields
\begin{equation} g(\hat{\lambda}_r, \hat{\lambda}_i) = \left(-\hat{\lambda}_r^2+\hat{\lambda}_i^2+\dfrac{\hat{\lambda}_i}{\hat{\lambda}_r^2+\hat{\lambda}_i^2}\right)+{i}\left(\dfrac{\hat{\lambda}_r}{\hat{\lambda}_r^2+\hat{\lambda}_i^2} - 2 \hat{\lambda}_r\hat{\lambda}_i\right), \end{equation}where the constant phase curves are given by
\begin{equation} {\rm Im}(g(\hat{\lambda}_r,\hat{\lambda}_i)) = \left(\dfrac{\hat{\lambda}_r}{\hat{\lambda}_r^2+ \hat{\lambda}_i^2} - 2 \hat{\lambda}_r\hat{\lambda}_i\right) = {\rm const}. \end{equation}
The saddle points corresponding to 
$\textrm {d}g(\hat {\lambda })/\textrm {d}\hat {\lambda } = 0$ are given by
\begin{equation} \hat{\lambda}_n = \dfrac{1}{\sqrt[3]{2}}\exp\left(\dfrac{-{\rm i}{\rm \pi}}{3}\left(\dfrac{1}{2}+2n\right)\right), \quad\textrm{where } n \in [0,1,2]. \end{equation}
The constant phase curves that pass through the saddle points can now be drawn by replacing the 
$\textrm {const}.$ in (B5) with 
$\textrm {Im}(g(\hat {\lambda }_n))$. That is,
\begin{equation} \left(\dfrac{\hat{\lambda}_r}{\hat{\lambda}_r^2+\hat{\lambda}_i^2} - 2 \hat{\lambda}_r\hat{\lambda}_i\right) = {\rm Im}(g(\hat{\lambda}_n)), \quad\textrm{where } n \in [0,1,2], \end{equation}
which are depicted in figure 14. We only consider the red and blue curves (corresponding to saddle points 1 and 2) to deform the original integration path (i.e. 
$\hat {\lambda } \in [-\infty, \infty ]$) into; the green curve that passes through saddle point 3 is a steepest-ascent curve, and hence, cannot be used. In other words,
\begin{equation} \int_{-\infty}^\infty {\rm d}\hat{\lambda}\,\hat{\lambda}^{5/2} \exp\left((\tilde{r}_t\tilde{z}^{1/2})^{2/3}g(\hat{\lambda})\right) = \left[\int_R + \int_B \right] {\rm d}\hat{\lambda}\,\hat{\lambda}^{5/2} \exp\left((\tilde{r}_t\tilde{z}^{1/2})^{2/3}g(\hat{\lambda})\right), \end{equation}
where 
$R$ and 
$B$ indicate the red and blue steepest-descent contours, respectively. Next, the integrals are approximated close to saddle points 1 and 2, respectively; as already mentioned, these regions will have the largest contribution. This is done by substituting 
$\hat {\lambda } = \hat {\lambda }_n + \varLambda$ and doing a Taylor series expansion of the integrand, for small 
$\varLambda$, whence one obtains
\begin{align}
&\left(\exp({(\tilde{r}_t\tilde{z}^{1/2})^{2/3}g(\hat{\lambda}_1)})\hat{\lambda}_1^{5/2}+
\exp({(\tilde{r}_t\tilde{z}^{1/2})^{2/3}g(\hat{\lambda}_2)})\hat{\lambda}_2^{5/2}\right)\nonumber\\ &\quad\times \int_{-\infty}
^{\infty} {\rm
d}\varLambda\exp\left({-}3(\tilde{r}_t\tilde{z}^{1/2})^{2/3}\varLambda^2\right),
\end{align}
at the leading order in 
$\varLambda$. Evaluating the Gaussian integral above yields
\begin{equation} \displaystyle \left(\exp({(\tilde{r}_t\tilde{z}^{1/2})^{2/3}g(\hat{\lambda}_1)}) \hat{\lambda}_1^{5/2}+\exp({(\tilde{r}_t\tilde{z}^{1/2})^{2/3}g(\hat{\lambda}_2)}) \hat{\lambda}_2^{5/2}\right)\sqrt{\dfrac{\rm \pi}{3\left(\tilde{r}_t\tilde{z}^{1/2}\right)^{2/3}}}. \end{equation}
Substituting 
$\hat {\lambda }_1$ and 
$\hat {\lambda }_2$ in the expression above, one obtains the Stokes streamfunction in the limit of 
$\tilde {r}_t\tilde {z}^{1/2} \gg 1$, as given in (3.6).
Figure 14. Constant phase curves in the complex-
$\hat {\lambda }$ plane from expression (B7). The blue, red and green curves (both dashed and continuous) correspond to saddle points 1, 2 and 3, respectively; however, we consider only the continuous curves. The red and blue continuous curves are the steepest-descent curves, whereas the green curve is the steepest-ascent curve.
B.2. Boundary of the columnar structure
As described in § 3.2, the boundary separating the columnar structure, and the far-field region described by the algebraic wake asymptote, is given by (see (3.8))
\begin{equation} \dfrac{1620\tilde{r}_t^2}{\tilde{z}^7} = \dfrac{\sqrt{3}}{\sqrt[3]{2}}\dfrac{(\tilde{r}_J\tilde{z}^{1/2})^{4/3}}{\tilde{z}^2}\exp\left(-\dfrac{3}{2\sqrt[3]{4}} (\tilde{r}_J\tilde{z}^{1/2})^{2/3}\right). \end{equation}This can be written as
\begin{equation} \dfrac{B}{\tilde{z}^6} = \dfrac{\exp\left({-}A \xi\right)}{\xi}, \end{equation}
where 
$A =3/2\sqrt [3]{4}$, 
$B = 1620\sqrt [3]{2}/\sqrt {3}$ are constants and 
$\xi = (\tilde {r}_J\tilde {z}^{1/2})^{2/3}$. Taking a logarithm on both sides gives
\begin{equation} \ln\left(\tilde{z}^6/B\right) = A\xi + \ln\xi. \end{equation}
For the columnar structure 
$\tilde {z} \gg 1$, and the left-hand side in (B15) is a logarithmically large quantity. Assuming 
$\xi \gg 1$ in anticipation, the first term on the right-hand side is dominant over the second. Therefore, at leading (logarithmic) order, one obtains
\begin{equation} \xi \sim L_1 \colon = \dfrac{1}{A}\left[\ln\left(\dfrac{\tilde{z}^6}{B}\right)\right]. \end{equation}
In order to obtain a better approximation, one writes 
$\xi = L_1 + L_2$, where 
$L_1$ is given in (B16) and the correction 
$L_2 \ll L_1$. Substituting 
$\xi = L_1 + L_2$ in (B15) gives
\begin{equation} \ln\left(\dfrac{\tilde{z}^6}{B}\right) = A L_1+ A L_2 + \ln(L_1+L_2). \end{equation}
The first term in the right-hand side cancels the left-hand side due to (B16). Next, the third term in (B17) can be written as 
$\ln (L_1+L_2) = \ln L_1 + \ln (1+L_2/L_1) \approx \ln L_1 + L_2/L_1$. Therefore, (B17) simplifies to
\begin{equation} A L_2 ={-} \ln L_1 - \dfrac{L_2}{L_1}. \end{equation}
Since 
$L_2/L_1 \ll L_2$, one obtains
\begin{equation} L_2 ={-}\dfrac{1}{A}\ln\left(L_1\right). \end{equation}
Using 
$L_1$ and 
$L_2$ from (B16) and (B19), one obtains
\begin{equation} \xi = L_1+L_2 = \dfrac{1}{A}\left[\ln\left(\dfrac{\tilde{z}^6}{B}\right) - \ln\left(\dfrac{1}{A}\ln\left(\dfrac{\tilde{z}^6}{B}\right)\right)\right], \end{equation}
to second order. Substituting 
$A, B$ and 
$\xi$ in (B20) results in an explicit expression for the radial extent of the columnar structure boundary (
$\tilde {r}_J$) in terms of 
$\tilde {z}$, given in (3.8). Figure 15(a) shows a comparison of the boundary calculated from 
$L_1$ alone (red curve), and using 
$L_1+L_2$ (black dashed curve), with that obtained from the streamline plot (blue curve). The two-term approximation is better than the single-term one, and shows an exact match with the numerically obtained boundary.
Figure 15. (a) Comparison of the columnar structure boundary obtained in Appendix B.2 with the boundary drawn from ends of different recirculating cells in the streamline plot of 
$\beta _\infty = 0$ (blue curve). The red and black dashed curves are from the single and two term asymptotic solution from (B16) and (B20), respectively. (b) Comparison of the tertiary screening length obtained in Appendix C and that from numerically solving (C1) (blue curve). The red, green and black curves are from the single-, two- and three-termasymptotic solutions given in (C4), (C6) and (C7), respectively.
Appendix C. Tertiary screening length
As described in § 3.2, we equate the large-argument form of the modified Bessel function expression of the axial velocity (along the rear stagnation streamline) with the expected far-field decay asymptote to obtain
\begin{equation} 3\beta_\infty^{1/2}\sqrt{\frac{\rm \pi}{2}}\dfrac{\exp\left({-}2\beta_\infty^{1/2}\tilde{z}_t\right)}{\left({-}2\beta_\infty^{1/2}\tilde{z}_t\right)^{1/2}} = \dfrac{3240}{\tilde{z}_t^7}. \end{equation}
This equation can be simplified by setting 
$\chi = \beta _\infty ^{1/2}\tilde {z}$ and 
$C = 2160\beta _\infty ^{3}/\sqrt {{\rm \pi} }$, so (C1) takes the form
\begin{equation} \exp({-}2 \chi) = \dfrac{C}{\chi^{13/2}}. \end{equation}Taking a logarithm on both sides of (C2) gives
\begin{equation} \ln \left(\dfrac{1}{C}\right) = 2\chi - \dfrac{13}{2}\ln(\chi). \end{equation}
In the convection-dominant limit (
$\beta _\infty \ll 1$), the left-hand side of (C3) is a logarithmically large quantity. Assuming 
$\chi \gg 1$ therefore, one finds that the first term in the right-hand side is dominant over the second; a reasonable assumption the secondary screening length is given by 
$l_s \sim {O}(al_c\beta _\infty ^{-1/2})$, and 
$\chi \gg 1$ implies that the tertiary screening length is asymptotically larger in relation. Therefore, one obtains, to the leading (logarithmic) order,
\begin{equation} \chi \sim M_1 \colon = \dfrac{1}{2}\ln\left(\dfrac{1}{C}\right). \end{equation}
Similar to Appendix B.2, a correction to 
$\chi$ can be found by substituting 
$\chi = M_1 + M_2$, where 
$M_1$ is given by (C4), in (C3). Grouping terms of the same order yields
\begin{equation} \left(\ln\left(\dfrac{1}{C}\right) - 2M_1\right) = \left(2M_2 - \dfrac{13}{2}\ln M_1\right) - \dfrac{13}{2}\dfrac{M_2}{M_1}. \end{equation}
In obtaining the expression above, 
$\ln (M_1+M_2)$ due to the second term in (C3) is written as 
$\ln M_1 + \ln (1+M_2/M_1) \approx \ln M_1 + M_2/M_1$. One then finds
\begin{equation} M_2 = \tfrac{13}{4}\ln M_1. \end{equation}
Figure 15(b) shows a comparison of the tertiary screening length calculated from 
$M_1$ alone (red curve), and using 
$M_1+M_2$ (green curve), with that obtained from numerically solving (C1) (blue curve). In contrast to 
$\tilde {r}_J$ above, there is a small deviation from the numerical results even relative to the two-term approximation. A better approximation is therefore obtained by considering 
$\chi = M_1 + M_2 + M_3$ and substituting in (C2). This results in
\begin{equation} \left(\ln\left(\dfrac{1}{C}\right)-2M_1\right) = \left(2M_2-\dfrac{13}{2}\ln M_1\right) + 2M_3 - \dfrac{13}{2}\ln\left(1+\dfrac{M_2+M_3}{M_1}\right). \end{equation}
Again, the left-hand side and the first two terms in the right-hand side cancel due to (C4) and (C6), respectively. Further, noticing that 
$M_3/M_1 \ll M_3$, one finds
\begin{equation} M_3 = \dfrac{13}{4}\dfrac{M_2}{M_1}. \end{equation}
Finally, substituting 
$M_1$, 
$M_2$ and 
$M_3$, yields
\begin{equation} \chi = M_1 + M_2 + M_3 = \left(\dfrac{1}{2}\ln\left(\dfrac{1}{C}\right)+\dfrac{13}{4}\ln\left(\dfrac{1}{2}\ln\left(\dfrac{1}{C}\right)\right) + \dfrac{169}{16}\dfrac{\ln\left(\dfrac{1}{2}\ln\left(\dfrac{1}{C}\right)\right)}{\dfrac{1}{2}\ln\left(\dfrac{1}{C}\right)} \right). \end{equation}The tertiary screening length calculated from the expression above (black dashed curve) matches well with that from numerically solving (C1) (blue curve), as shown in figure 15(b).
