2. A sedimenting spheroid in a linearly stratified ambient: the generalized reciprocal theorem formulation

The torque acting on a spheroid, sedimenting in a stably stratified fluid ambient, is derived in the following using the generalized reciprocal theorem (see Kim & Karrila Reference Kim and Karrila1991; Dabade et al. Reference Dabade, Marath and Subramanian2015; Dabade, Marath & Subramanian Reference Dabade, Marath and Subramanian2016). The theorem relates two pairs of stress and velocity fields, and may be stated in the form:

(2.1)\begin{equation} \displaystyle\int_{S_p} \sigma^{(2)}_{ij} u'^{(1)}_i n_j \,\textrm{d}S - \displaystyle\int_{S_p} \sigma^{(1d)}_{ij} u^{(2)}_i n_j \,\textrm{d}S =\displaystyle\int \frac{\partial \sigma_{ij}^{(1d)}}{\partial x_j} u^{(2)}_i \,\textrm{d}V, \end{equation}

where $S_p$ denotes the surface of the spheroid, and with the neglect of the surface integrals at infinity, the volume integral on the right-hand side of (2.1) is over the unbounded fluid domain external to the spheroid. In (2.1), the pair $({\boldsymbol \sigma }^{(1d)},{\boldsymbol u}'^{(1)})$ denotes the dynamic stress and velocity fields associated with the problem of interest, namely a torque-free spheroid sedimenting under gravity in an ambient linearly stratified medium for small Reynolds ($Re$) and viscous Richardson ($Ri_v$) numbers, with ${\boldsymbol u}'^{(1)}$ corresponding to the lab reference frame with a quiescent far-field ambient (the prime indicates that the velocity field has a disturbance-like character in this reference frame, and decays away in the far-field). The Reynolds and viscous Richardson numbers measure the relative importance of inertial and buoyancy forces relative to viscous forces, respectively, and the aforementioned smallness of these parameters corresponds therefore to the case where inertia and stratification act as weak perturbing influences about a leading-order Stokesian approximation; note, however, that the result for the stratification torque obtained in § 4.1 is an exception to this general assumption, in that it continues to be valid for finite values of the Richardson number. The Reynolds and Richardson numbers are defined later in this section when writing down the non-dimensional system of governing equations; the precise definition of the dynamic stress field, ${\boldsymbol \sigma }^{(1d)}$, is also provided at the same place.

The pair $({\boldsymbol \sigma }^{(2)},{\boldsymbol u}^{(2)})$, that defines the test problem in (2.1), corresponds to the stress and velocity fields associated with inertialess ($Re =0$) rotation of the same spheroid, about an axis orthogonal to its axis of symmetry, in a homogeneous and otherwise quiescent ambient with the same (assumed constant) viscosity as the medium in the actual problem. The equations governing the test problem may be written as

(2.2)\begin{gather} \frac{\partial u^{(2)}_i}{\partial x_i} =0, \end{gather}

(2.3)\begin{gather}\mu \frac{\partial^2 u^{(2)}_{i}}{\partial x_j^2} - \frac{\partial p^{(2)}}{\partial x_i} =0, \end{gather}

with the boundary condition ${\boldsymbol u}^{(2)} = {\boldsymbol \varOmega }^{(2)} \wedge {\boldsymbol x}$ on $S_p$, ${\boldsymbol \varOmega }^{(2)}$ being the angular velocity of the spheroid in the test problem, and far-field decay conditions for both ${\boldsymbol u}^{(2)}$ and $p^{(2)}$. Note that the test velocity (${\boldsymbol u}^{(2)}$) and stress (${\boldsymbol \sigma }^{(2)}$) fields decay as $O({1}/{r^2})$ and $O({1}/{r^3})$, respectively, $r$ being the distance away from the spheroid; this decay, together with the decaying (dynamic) stress field in the problem of interest, justifies the neglect of the surface integrals at infinity in (2.1). Use of the aforementioned surface boundary condition in (2.1) leads to

(2.4)\begin{equation} \displaystyle\int_{S_p} \sigma^{(2)}_{ij} u'^{(1)}_i n_j \,\textrm{d}S - \varOmega^{(2)}_j \displaystyle\int_{S_p} \epsilon_{ijk} x_k \sigma^{(1d)}_{il} n_l \,\textrm{d}S =\displaystyle\int \frac{\partial \sigma^{(1d)}_{ij}}{\partial x_j} u^{(2)}_i \,\textrm{d}V, \end{equation}

where the second integral on the left-hand side in (2.4) now denotes the torque due to the dynamic stress field ${\boldsymbol \sigma }^{(1d)}$. We postpone further simplification of (2.4) until after we define the pair $({\boldsymbol \sigma }^{(1d)},{\boldsymbol u}'^{(1)})$.

As mentioned previously, problem 1 corresponds to an arbitrarily oriented spheroid, sedimenting under the action of a gravitational force $F\hat {\boldsymbol g}$, in an ambient fluid that is linearly stratified (along $\hat {\boldsymbol g}$) in the absence of the fluid motion induced by the spheroid. The unit vector $\hat {\boldsymbol g}$ is aligned along gravity, with $g$ denoting the magnitude of the gravitational acceleration, and $F = ({4{\rm \pi} }/{3})Lb^2\Delta \rho g (({4{\rm \pi} }/{3})L^2b \Delta \rho g)$ denoting the buoyant weight for a prolate (oblate) spheroid. Here, $L$ and $b$ are the semi-major and semi-minor axes of the spheroid, with $\kappa = L/b$ and $b/L$ being the aspect ratios of prolate and oblate spheroids, respectively; thus, $\kappa >1$ and $\kappa <1$ for the prolate and oblate cases. The density difference that enters the buoyant weight above is $\Delta \rho = \rho _s -\rho ^{(1)}_\infty ({\boldsymbol x}_c)$, with $\rho _s$ being the density of the spheroid (assumed homogeneous), and $\rho ^{(1)}_\infty ({\boldsymbol x}_c) = \rho _0$ being the ambient fluid density at the centre of the spheroid. The latter simplification arises because of the linear stratification and the fore–aft symmetry of the spheroid, both of which imply that the weight of the equivalent stratified spheroidal fluid blob that gives the buoyant force, within an Archimedean interpretation, is the same as the weight of a homogeneous fluid blob with density equal to the ambient value at the spheroid centre. In a lab-fixed reference frame, the ambient density field in problem 1 may be written in the form:

(2.5)\begin{equation} \rho^{(1)}_\infty({\boldsymbol x}^L) =\rho_0 + \gamma x^L_i \hat{g}_i, \end{equation}

where ${\boldsymbol x}^L$ denotes the position vector in laboratory coordinates with the spheroid centre as the origin, and $\gamma > 0$ is the constant density gradient that characterizes the stable ambient stratification. The calculations for the torque are, however, best done in a reference frame translating with the spheroid where a quasi-steady state is assumed to prevail at leading order. The latter assumption is motivated by the asymptotically weak rotation of the sedimenting spheroid in the limit $Re, Ri_v \ll 1$. The precise condition for the quasi-steady state assumption to hold depends on $Pe$, being more restrictive for large $Pe$, and is stated later alongside the results for the spheroid angular velocity for small and large $Pe$, obtained in the following.

The ambient density in the particle-fixed reference frame takes the form:

(2.6)\begin{equation} \rho^{(1)}_\infty({\boldsymbol x}) =\rho_0 +\gamma ( x_i + U_i t)\hat{g}_i, \end{equation}

${\boldsymbol x}$ being the position vector in the new reference frame. In (2.6), ${\boldsymbol U}$ is the spheroid settling velocity, and related to the force ($F\hat {\boldsymbol g}$) via a mobility tensor that is a known function of the spheroid aspect ratio $\kappa$. In terms of the spheroid orientation vector ${\boldsymbol p}$, one may write ${\boldsymbol U} = ({1}/{\mu L}) [X_A^{-1}{\boldsymbol p}{\boldsymbol p}+ Y_A^{-1}({\boldsymbol I} - {\boldsymbol p}{\boldsymbol p})] \boldsymbol {\cdot } (F\hat {\boldsymbol g})$, $X_A(\kappa )$ and $Y_A(\kappa )$ being the non-dimensional axial and transverse translational resistance functions. The aspect ratio dependence of these functions is well known (see Kim & Karrila Reference Kim and Karrila1991), and is given in Appendix A for convenient reference. Note that the ambient density at the centre of the spheroid (${\boldsymbol x} = 0$) is given by $\rho _0 + \gamma (U_i\hat {g}_i)t$, the time dependence arising from the spheroid translation. The equations of motion for problem 1, within a Boussinesq framework where the fluid density multiplying the inertial terms is taken as a constant $\rho _0$ (say), may be written as

(2.7)\begin{gather} \frac{\partial u^{(1)}_i}{\partial x_i} =0, \end{gather}

(2.8)\begin{gather}\mu \frac{\partial^2 u^{(1)}_{i}}{\partial x_j^2} - \frac{\partial p^{(1)}}{\partial x_i} = \rho_0 u^{(1)}_j\frac{\partial u^{(1)}_i}{\partial x_j} - \rho^{(1)}g_i, \end{gather}

(2.9)\begin{gather}\frac{\partial \rho^{(1)}}{\partial t} + u^{(1)}_j\frac{\partial \rho^{(1)}}{\partial x_j} =D \nabla^2 \rho^{(1)}, \end{gather}

where $D$ is the diffusivity of the stratifying agent (Candelier, Mehaddi & Vauquelin Reference Candelier, Mehaddi and Vauquelin2014; Mehaddi et al. Reference Mehaddi, Candelier and Mehlig2018; Shaik & Ardekani Reference Shaik and Ardekani2020). One now defines the perturbation density ($\rho '^{(1)}$) via $\rho ^{(1)} = \rho _0 + \gamma (x_i + U_i t)\hat {g}_i + \rho '^{(1)}$. Next, using the scales $U=F/(\mu L X_A)$ for the velocity, $L$ for the length, $\mu U/L$ for the pressure and $\gamma L$ for $\rho '^{(1)}$, one obtains the following system of non-dimensional equations:

(2.10)\begin{gather} \frac{\partial u^{(1)}_i}{\partial x_i} =0, \end{gather}

(2.11)\begin{gather}\frac{\partial^2 u^{(1)}_{i}}{\partial x_j^2} - \frac{\partial p^{(1)}}{\partial x_i} + \frac{\rho_0 g L^2}{\mu U}\hat{g}_i + Ri_v (\hat{U}_j t +x_j )\hat{g}_j\hat{g}_i= Re \,u^{(1)}_j\frac{\partial u^{(1)}_i}{\partial x_j} - Ri_v\rho'^{(1)}\hat{g}_i, \end{gather}

(2.12)\begin{gather}u^{(1)}_j\frac{\partial \rho'^{(1)}}{\partial x_j} + (\hat{U}_j + u^{(1)}_j)\hat{g}_j =\frac{1}{Pe} \nabla^2 \rho'^{(1)}, \end{gather}

where $Re = \rho _0UL/\mu$, $Ri_v = \gamma L^3g/(\mu U)$ and $Pe=UL/D$ are the Reynolds, viscous Richardson and Péclet numbers, respectively; note that $Ri_v = Re/Fr^2$, where $Fr$ is the Froude number, and the usual measure of the importance of stratification in the inviscid limit (Turner Reference Turner1973). In (2.10)–(2.12), we continue to use the same notation for the dimensionless fields for simplicity. The velocity fields, in the lab reference frame used in (2.1), and in the particle-fixed reference frame adopted in (2.10)–(2.12) are related as ${\boldsymbol u}'^{(1)} = (\hat {\boldsymbol U} + {\boldsymbol u}^{(1)})$, with $\hat {\boldsymbol U}= [p_i p_j+X_A/Y_A(\delta _{ij}-p_ip_j)]\hat {g}_j$, now being a dimensionless vector along the direction of settling, and $- \hat {\boldsymbol U}$ therefore being the far-field ambient flow in the particle-fixed frame; note that $\hat {\boldsymbol U}$ is not a unit vector for an arbitrarily oriented spheroid, and reduces to one only for a spheroid aligned with gravity. Thus, the combination $(\hat {\boldsymbol U} + {\boldsymbol u}^{(1)}) \boldsymbol {\cdot } \hat {\boldsymbol g}$ in (2.12) denotes the convection of the (constant) base-state density gradient by the component of the disturbance velocity field (${u'}_3^{(1)}$) along gravity. Finally, the time dependence of ${\boldsymbol u}^{(1)}$ in (2.11), and that of $\rho '^{(1)}$ in (2.12) in particular, that arise from the (slow) rotation of the spheroid, have been neglected owing to the quasi-steady state assumption made in (2.10)–(2.12); the time dependence of the density multiplying the inertial terms, on account of spheroid translation, has already been neglected within the Boussinesq approximation.

One now defines a disturbance pressure field ($p'^{(1)}$) via $p^{(1)} = p_0^{(1)} + p'^{(1)}$ with

(2.13)\begin{equation} \frac{\partial p_0^{(1)}}{\partial x_i} =\frac{\rho_0 g L^2}{\mu U}\hat{g}_i + Ri_v (\hat{U}_j t +x_j) \hat{g}_j\hat{g}_i, \end{equation}

so that $p^{(1)}_0$ defines the baseline hydrostatic contribution arising from the ambient linear stratification. Having incorporated the baseline hydrostatic variation in $p^{(1)}_0$, one may write the governing equations above in terms of the disturbance velocity, pressure and density fields as follows:

(2.14)\begin{gather} \frac{\partial u'^{(1)}_i}{\partial x_i} =0, \end{gather}

(2.15)\begin{gather}\frac{\partial \sigma^{(1d)}_{ij}}{\partial x_j} = Re (-\hat{U}_j+ u'^{(1)}_j)\frac{\partial u'^{(1)}_i}{\partial x_j} - Ri_v \rho'^{(1)}\hat{g}_i, \end{gather}

(2.16)\begin{gather}(-\hat{U}_j + u'^{(1)}_j)\frac{\partial \rho'^{(1)}}{\partial x_j} ={-}u'^{(1)}_j\hat{g}_j + \frac{1}{Pe} \nabla^2\rho'^{(1)}. \end{gather}

where the left-hand side of (2.11) has been written in terms of the dynamic stress field ${\boldsymbol \sigma }^{(1d)}$ defined by ${\boldsymbol \sigma }^{(1d)} = -p'^{(1)}{\boldsymbol I} + ( {\boldsymbol \nabla } {\boldsymbol u}'^{(1)} + {\boldsymbol \nabla } {{\boldsymbol u}'^{(1)}}^{\dagger} )$. Thus, one has the relation ${\boldsymbol \sigma }^{(1)} = -p_0^{(1)} {\boldsymbol I} + {\boldsymbol \sigma }^{(1d)}$ between the total (${\boldsymbol \sigma }^{(1)}$) and the dynamic stress fields of problem 1.

Assuming the spheroid in problem 1 to rotate with an angular velocity ${\boldsymbol \varOmega }^{(1)}$, one has the boundary condition ${\boldsymbol u}'^{(1)} = {\boldsymbol \varOmega }^{(1)} \wedge {\boldsymbol x}$ on $S_p$. Using this in the first surface integral in (2.4), and substituting the divergence of the dynamic stress from (2.15) in the volume integral in (2.1), one obtains

(2.17)\begin{align} \varOmega^{(1)}_j {\mathcal{L}}_j^{(2)} -\varOmega^{(2)}_j {\mathcal{L}}_j^{\sigma(1)d} =Re\displaystyle\int u^{(2)}_i \, (-\hat{U}_j+ u'^{(1)}_j)\frac{\partial u'^{(1)}_i}{\partial x_j} \,\textrm{d}V - Ri_v \displaystyle\int \rho'^{(1)} \hat{g}_i\, u_{i}^{(2)}\,\textrm{d}V, \end{align}

where $\boldsymbol {\mathcal {L}}^{\sigma (1)d}$ now denotes the torque contribution due to the dynamic stress ${\boldsymbol \sigma }^{(1d)}$. Now, the particle in problem 1 is torque-free. In light of the above relation between ${\boldsymbol \sigma }^{(1)}$ and ${\boldsymbol \sigma }^{(1d)}$, the total torque ($\boldsymbol {\mathcal {L}}^{(1)}$) may be written as $\boldsymbol {\mathcal {L}}^{(1)} = \boldsymbol {\mathcal {L}}^{\sigma (1)d} + \boldsymbol {\mathcal {L}}^{\sigma (1)s}$, where the dynamic torque component $\boldsymbol {\mathcal {L}}^{\sigma (1)d}$ includes both inertia and stratification-induced contributions, whereas $\boldsymbol {\mathcal {L}}^{\sigma (1)s}$ is the hydrostatic contribution owing to the pressure field $p^{(1)}_0$ associated with the linearly varying density field of the stably stratified ambient, and defined by (2.13). Thus, $\boldsymbol {\mathcal {L}}^{(1)} = 0 \Rightarrow \boldsymbol {\mathcal {L}}^{\sigma (1)d}=-\boldsymbol {\mathcal {L}}^{\sigma (1)s}$, and the relation involving the spheroid angular velocity in problem 1 takes the following form:

(2.18)\begin{align} \varOmega^{(1)}_j {\mathcal{L}}_j^{(2)}= Re\int u_{i}^{(2)} (-\hat{U}_j+ u'^{(1)}_j)\frac{\partial{u'^{(1)}_i}}{\partial x_j} \,\textrm{d}V - \left[\varOmega^{(2)}_j {\mathcal{L}}_j^{\sigma(1)s} + Ri_v \int \rho'^{(1)} \hat{g}_i\, u_{i}^{(2)}\,\textrm{d}V \right], \end{align}

where

(2.19)\begin{equation} {\mathcal{L}}_i^{\sigma(1)s} ={-}\epsilon_{ijk} \displaystyle\int_{S_p} p_0^{(1)} x_j n_k \,\textrm{d}S, \end{equation}

with $p_0^{(1)}$ being defined in (2.13). As the buoyancy force in a homogeneous ambient acts through the centre of the spheroid, only the linearly varying term in (2.13) contributes to the hydrostatic torque, which may therefore be written as

(2.20)\begin{equation} {\mathcal{L}}_i^{\sigma(1)s} ={-}Ri_v\, \epsilon_{ijk} \displaystyle\int_{S_p} \frac{1}{2}(x_l\hat{g}_l)^2x_j n_k \,\textrm{d}S, \end{equation}

the contribution above remaining the same regardless of the choice of reference frame (${\boldsymbol x}$ or ${\boldsymbol x}^L$). On substitution of the above expression for ${\mathcal {L}}_k^{\sigma (1)s}$, and using the relation ${\boldsymbol u}^{(2)}_i = {\boldsymbol U}^{(2)}_{ij}{\boldsymbol \varOmega }^{(2)}_j$ (on account of the linearity of the Stokes equations), the second-order tensor ${\boldsymbol U}^{(2)}_{ij}$ being known in closed form (see Dabade et al. Reference Dabade, Marath and Subramanian2015, and § 3), (2.18) takes the form

(2.21)\begin{gather} \varOmega^{(1)}_j {\mathcal{L}}_j^{(2)}= \varOmega^{(2)}_k\left\{ Re\int U_{jk}^{(2)} \, (-\hat{U}_l+ u'^{(1)}_l)\frac{\partial{u'^{(1)}_j}}{\partial x_l} \,\textrm{d}V \right.\nonumber\\ \left. - Ri_v \left[-\frac{1}{2} \epsilon_{klm} \displaystyle\int_{S_p} (x_j\hat{g}_j)^2 x_l n_m \,\textrm{d}S+ \int \rho'^{(1)} \hat{g}_j\, U_{jk}^{(2)}\,\textrm{d}V \right] \right\}, \end{gather}

Again, on account of linearity, one may write the torque on the rotating spheroid in the test problem, in the form $\boldsymbol {\mathcal {L}}^{(2)} =[X_C {\boldsymbol p}{\boldsymbol p} + Y_C({\boldsymbol I} - {\boldsymbol p}{\boldsymbol p})] \boldsymbol {\cdot } \boldsymbol {\varOmega }^{(2)}$, where $X_C(\kappa )$ and $Y_C(\kappa )$ are the non-dimensional axial and transverse rotational resistance functions, and are known functions of $\kappa$ (Kim & Karrila Reference Kim and Karrila1991) whose expressions are given in Appendix A. By symmetry, the sedimenting spheroid cannot spin about its axis regardless of its orientation, and therefore without loss of generality, the test problem can be taken as that of a transverse rotation in the inertialess limit (${\boldsymbol \varOmega }^{(2)} \boldsymbol {\cdot } {\boldsymbol p} = 0$), in which case the test torque–angular velocity relation takes the simpler form $\boldsymbol {\mathcal {L}}^{(2)} = Y_C \boldsymbol {\varOmega }^{(2)}$. Finally, accounting for the fact that the test angular velocity ${\boldsymbol \varOmega }^{(2)}$ can point in an arbitrary direction in a plane perpendicular to ${\boldsymbol p}$, one obtains the following relation for the spheroid angular velocity in problem 1:

(2.22)\begin{gather} \varOmega^{(1)}_i= \frac{1}{Y_C}\left\{ Re\int U_{ji}^{(2)} \, (-\hat{U}_l+ u'^{(1)}_l)\frac{\partial{u'^{(1)}_j}}{\partial x_l} \,\textrm{d}V + Ri_v \left[ \epsilon_{ilm} \displaystyle\int_{S_p}\frac{1}{2} (x_j\hat{g}_j)^2 x_l n_m \,\textrm{d}S \right. \right. \nonumber\\ -\left. \left. \displaystyle\int \rho'^{(1)} \hat{g}_j\, U_{ji}^{(2)}\,\textrm{d}V \right] \right\}. \end{gather}

As a settling spheroid in a homogeneous ambient must retain its initial orientation in the Stokes limit on account of reversibility, expectedly, the rotation of the spheroid, as given by (2.22), arises due to the combined (weak) effects of fluid inertia and the ambient stratification. The first term within the curly brackets in (2.22) corresponds to the inertial torque, whereas the second and third terms which have been grouped together (within square brackets) correspond to the hydrostatic and hydrodynamic components of the stratification torque, respectively. The hydrostatic torque only involves knowledge of the ambient density field, and is easily evaluated. The inertial torque has a regular character in that the dominant contributions to the $O(Re)$ volume integral in (2.22) arise from a volume of $O(L^3)$ around the sedimenting spheroid, and therefore, the integral may again readily be determined at leading order using Stokesian approximations for the velocity fields involved, as has been done in Dabade et al. (Reference Dabade, Marath and Subramanian2015). The evaluation of these two simpler contributions is detailed in the next section. The nature of the hydrodynamic torque arising from the perturbed stratification depends crucially on $Pe$, and this more complicated calculation is given in §§ 4.1 and 4.2, for small and large $Pe$, respectively.

3. The spheroidal angular velocity due to the inertial and hydrostatic torque contributions

The $O(Re)$ inertial angular velocity in (2.22) has recently been calculated for spheroids, both prolate and oblate, of an arbitrary aspect ratio (see Dabade et al. Reference Dabade, Marath and Subramanian2015). Although the analysis in Dabade et al. (Reference Dabade, Marath and Subramanian2015) pertains to the limit $Re \ll 1$, the results have been shown to remain qualitatively valid even for $Re$ of order unity (see Jiang et al. Reference Jiang, Zhao, Andersson, Gustavsson, Pumir and Mehlig2020). As mentioned previously, this angular velocity has a regular character, as may be seen from the convergence of the inertial volume integral in (2.18) based on the leading-order Stokesian estimate for the integrand. As argued in Dabade et al. (Reference Dabade, Marath and Subramanian2015), the inertial acceleration ${\boldsymbol u}^{(1)} \boldsymbol {\cdot } \boldsymbol {\nabla } {\boldsymbol u'}^{(1)} \sim \hat {\boldsymbol U} \boldsymbol {\cdot } \boldsymbol {\nabla } {\boldsymbol u}'^{(1)} \sim O(1/r^2)$ for distances large compared with $L$, or $r \gg 1$ in dimensionless terms, on using ${\boldsymbol u}'^{(1)} \sim O(1/r)$ for the Stokeslet field due to the translating spheroid. The test velocity field ${\boldsymbol u}^{(2)}$ due to the rotating spheroid has the character of a rotlet-cum-stresslet in the far-field, and is therefore $O(1/r^2)$. This leads to an integrand that decays as $\hat {\boldsymbol U} \boldsymbol {\cdot } \boldsymbol {\nabla } {\boldsymbol u}'^{(1)} \boldsymbol {\cdot } {\boldsymbol u}^{(2)} \sim O(1/r^4)$ for $r \gg 1$, implying a convergent volume integral. This volume integral has been evaluated in closed form using spheroidal coordinates in Dabade et al. (Reference Dabade, Marath and Subramanian2015). For the prolate case, the spheroidal coordinates $(\xi,\eta,\phi )$ are defined by the relations: $x_1 + \mathrm {i} x_2 = \textrm {d}\bar {\xi }\bar {\eta } \exp (\textrm {i} \phi )$, $x_3 = \textrm {d}\xi \eta$, with the three axes of the Cartesian system ($\boldsymbol {1}_3$) aligned with the spheroid axis of symmetry. Here, $1 \leq \xi < \infty$, $|\eta | \leq 1$ and $0 \leq \phi < 2{\rm \pi}$, with $\bar {\xi } = (\xi ^2-1)^{{1}/{2}}$ and $\bar {\eta } = (1-\eta ^2)^{{1}/{2}}$. The constant-$\xi$ surfaces correspond to confocal prolate spheroids and the constant-$\eta$ surfaces to confocal two-sheeted hyperboloids, both with the interfoci distance $2d$, and the constant-$\phi$ surfaces are planes passing through the axis of symmetry. The corresponding expressions for the oblate case may be obtained by the substitutions $d \leftrightarrow -\mathrm {i}d$, $\xi \leftrightarrow \mathrm {i}\bar {\xi }$, the constant-$\xi$ and $\eta$ surfaces now being confocal oblate spheroids and single-sheeted hyperboloids, respectively. In either case, the spheroid is the surface $\xi = \xi _0$, its aspect ratio being given by $\kappa = {\xi _0}/{\overline {\xi _0}}$ and ${\bar {\xi }_0}/{\xi _0}$ for prolate and oblate spheroids; thus, the near-spherical limit ($\kappa \rightarrow 1$) for either prolate or oblate spheroids corresponds to $\xi _0 \rightarrow \infty$, whereas the slender fiber ($\kappa \rightarrow \infty$) and flat disk ($\kappa \rightarrow 0$) limits correspond to $\xi _0 \rightarrow 1$. The fluid domain in the volume integrals in (2.22) corresponds to $\xi \geq \xi _0$.

For a prolate spheroid, the actual velocity field ${\boldsymbol u'}^{(1)}$ and the test velocity field tensor ${\boldsymbol U}^{(2)}$ in (2.22), may be written in the form (see Dabade et al. Reference Dabade, Marath and Subramanian2015, Reference Dabade, Marath and Subramanian2016):

(3.1)\begin{equation} {\boldsymbol u'}^{(1)} = \left(\frac{\hat{\boldsymbol{U}}\boldsymbol{\cdot}\boldsymbol{1}_3}{\xi_0 Q_1^1(\xi_0)+Q_0^0(\xi_0)} \right)\boldsymbol{S}_{1,0}^{(3)} +\left(\frac{\hat{\boldsymbol{U}}\boldsymbol{\cdot}\boldsymbol{1}_1}{3 Q_0^0(\xi_0)- \xi_0 Q_1 ^0(\xi_0)} \right)(\boldsymbol{S}_{1,1}^{(3)}-\boldsymbol{S}_{1,-1}^{(3)}), \end{equation}

(3.2)\begin{equation} {\boldsymbol U}^{(2)}=\boldsymbol{1}_2 \left(\frac{d(2\xi_0^2-1)(\boldsymbol{S}_{1,1}^{(2)}-\boldsymbol{S}_{1,-1}^{(2)})}{\left[2 \xi_0 Q_1^0(\xi_0)-\sqrt{\xi_0^2-1}Q_1^1(\xi_0)\right]}+\frac{d\left[\xi_0 Q_1^1(\xi_0)+2 \sqrt{\xi_0^2-1}Q_1^0(\xi_0)\right](\boldsymbol{S}_{2,1}^{(3)}-\boldsymbol{S}_{2,-1}^{(3)})}{Q_2^1(\xi_0)\left[2 \xi_0 Q_1^0(\xi_0)-\sqrt{\xi_0^2-1}Q_1^1(\xi_0)\right]}\right) . \end{equation}

The $\boldsymbol {S}^{(3)}_{t,s}$ and $\boldsymbol {S}^{(2)}_{t,s}$ in (3.1) and (3.2) denote the decaying (biharmonic and harmonic) vectorial solutions of the Stokes equations in spheroidal coordinates, and are given in Appendix C; the $Q_t^s(\xi )$ denote the associated Legendre functions of the second kind. Using these expressions, the volume integration may then be carried out analytically, and the inertial angular velocity ($\boldsymbol {\varOmega }^{(1)I}$) is given by

(3.3)\begin{equation} \varOmega_i^{(1)I}= Re \left[ \frac{F^p_I(\xi_0)X_A}{Y_C Y_A} ( \epsilon_{ijk}\hat{g}_j p_k \, \hat{g}_lp_l)\right],\end{equation}

for prolate spheroids. The corresponding expression for the oblate case may be obtained by the aforementioned substitutions, namely $d \leftrightarrow -\mathrm {i}d$, $\xi _0 \leftrightarrow \mathrm {i}\bar {\xi }_0$ in the dimensional angular velocity, and is given by

(3.4)\begin{equation} \varOmega_i^{(1)I}= Re\left[ \frac{F^o_I(\xi_0)X_A}{Y_C Y_A}(\epsilon_{ijk}\hat{g}_j p_k \hat{g}_lp_l) \right]. \end{equation}

The expressions for $F^p_I(\xi _0)$ and $F^o_I(\xi _0)$, as functions of the spheroid eccentricity ($e=1/\xi _0$), were first obtained by Dabade et al. (Reference Dabade, Marath and Subramanian2015), and are given in Appendix A. The inertial angular velocity given by (3.3) and (3.4) orients sedimenting spheroids broadside-on regardless of $\kappa$. The combination of the aspect-ratio-dependent functions, $F^{p/o}_I(\xi _0) X_A/(Y_CY_A)$, that multiplies $Re(\hat {\boldsymbol g} \boldsymbol {\cdot }{\boldsymbol p})(\hat {\boldsymbol g} \wedge {\boldsymbol p})$, and that determines the $\kappa$-dependence of the inertial angular velocities above, is plotted as a function of the eccentricity in figure 1, for both the prolate and oblate cases. One obtains the expected $O(1/\xi _0^2)$ scaling in the near-sphere limit ($\xi _0 \rightarrow \infty$); at the other extreme($\xi _0 \rightarrow 1$), the inertial angular velocity approaches zero as $O[\ln (\xi _0-1)]^{-1}$ in the slender fiber limit, consistent with viscous slender body theory (Khayat & Cox Reference Khayat and Cox1989; Subramanian & Koch Reference Subramanian and Koch2005), while remaining finite in the limit of a flat disk.

Figure 1.

The functions $({F^p_I(\xi _0)X_A})/({Y_C Y_A})$ and $({F^o_I(\xi _0)X_A})/({Y_C Y_A})$, that characterize the aspect ratio dependence of the inertial contributions to the angular velocities of prolate and oblate spheroids in (3.3) and (3.4), plotted as a function of the spheroid eccentricity.

The hydrostatic component of the stratification torque is also readily evaluated in spheroidal coordinates. For the prolate case, the dimensionless position vector that appears in (2.20) is given by $\boldsymbol {x} = ({\bar {\xi }_0}/{\xi _0})\bar {\eta }(\cos \phi {\boldsymbol 1}_1 + \sin \phi {\boldsymbol 1}_2) + \eta {\boldsymbol 1}_3$, and the unit normal is

(3.5)\begin{equation} \boldsymbol{n} = {\boldsymbol 1}_\xi = \frac{\xi_0\bar{\eta}}{\sqrt{\xi_0^2 -\eta^2}}(\cos \phi {\boldsymbol 1}_1 + \sin \phi {\boldsymbol 1}_2) + \frac{\overline{\xi_0}\eta}{\sqrt{\xi_0^2 -\eta^2}} {\boldsymbol 1}_3. \end{equation}

Using these expressions, and the areal element $\textrm {d}S = h_\eta h_\phi \,\textrm {d}\eta \,\textrm {d}\phi$, with $h_\eta = {\sqrt {\xi _0^2- \eta ^2}}/{\xi _0\bar {\eta }}$ and $h_\phi = {\bar {\xi }_0\bar {\eta }}/{\xi _0}$, one obtains the angular velocity due to the hydrostatic torque as

(3.6)\begin{equation} \varOmega_i^{(1)s} = Ri_v\frac{4{\rm \pi}}{15Y_C}\frac{1-\xi_0^2 }{\xi_0^4}(\epsilon_{ijk} \hat{g}_j p_k)(\hat{g}_lp_l), \end{equation}

for the prolate case, and using the transformations mentioned previously,

(3.7)\begin{equation} \varOmega_i^{(1)s} =Ri_v\frac{4{\rm \pi}}{15Y_C}\frac{\sqrt{\xi_0^2 -1}}{\xi_0^3}(\epsilon_{ijk} \hat{g}_j p_k)(\hat{g}_lp_l), \end{equation}

for the oblate case. The angular velocities given by (3.6) and (3.7) also orient the spheroid broadside-on like the inertial torque above. The aspect-ratio-dependent functions that multiply $Ri_v(\hat {\boldsymbol g} \boldsymbol {\cdot } {\boldsymbol p})(\hat {\boldsymbol g} \wedge {\boldsymbol p})$ in (3.6) and (3.7) are plotted as functions of $\xi _0$ in figure 2. As the hydrostatic torque is only a function of the particle geometry, these aspect ratio functions are algebraically small in both the near-sphere, and the slender fiber and flat-disk limits.

Figure 2.

The aspect-ratio-dependent functions multiplying $Ri_v(\hat {\boldsymbol g} \boldsymbol {\cdot } {\boldsymbol p})(\hat {\boldsymbol g} \wedge {\boldsymbol p})$ in (3.6) and (3.7), that characterize the hydrostatic contributions for prolate and oblate spheroids, plotted as a function of the spheroid eccentricity.

The inertial and hydrostatic angular velocities above have an identical angular dependence, of the form $(\hat {\boldsymbol g} \boldsymbol {\cdot } {\boldsymbol p})(\hat {\boldsymbol g} \wedge {\boldsymbol p})$, one which is easily inferred based on the requirement that the angular velocity be a pseudovector quadratic in $\hat {\boldsymbol g}$ (Dabade et al. Reference Dabade, Marath and Subramanian2015). The dependence implies that the maximum angular velocity occurs midway between the horizontal ($\hat {\boldsymbol g} \boldsymbol {\cdot } {\boldsymbol p} = 0$) and vertical ($\hat {\boldsymbol g} \boldsymbol {\cdot } {\boldsymbol p} = 1$) orientations. The hydrostatic torque arises because the point of action of the upward buoyant force, the centre of mass of the equivalent stratified fluid blob (in the Archimedean interpretation) lies below the geometric centre through which the weight of the spheroid acts vertically downward. The two forces therefore constitute a couple that turns the spheroid broadside-on. The broadside-on nature of the inertial torque is on account of ‘wake-shielding’: the wake associated with the front portion of the spheroid shields the rear, which catches up with the front as a result. As pointed out in Dabade et al. (Reference Dabade, Marath and Subramanian2015), this is not literally true for small $Re$. A signature of the wake arises only on length scales greater than $O(LRe^{-1})$, the Oseen region, in contrast to the scaling arguments above which show that the $O(Re)$ inertial torque arises from fluid inertial forces in a region of $O(L^3)$ (the inner region) around the sedimenting spheroid. Nevertheless, the velocity field in the inner region reflects the asymmetry of the outer Oseen field, and the sense of rotation remains the same for small $Re$. Importantly, the broadside-on nature of the inertial and hydrostatic torques imply that the transition from broadside-on to edgewise settling, observed in the recent experiments (see Mercier et al. Reference Mercier, Wang, Péméja, Ern and Ardekani2020) discussed in the introduction, must depend entirely on the hydrodynamic component of the stratification torque, that is, the second term within square brackets on the right-hand side in (2.22). While the calculation above shows the hydrostatic component to be $O(Ri_v)$, consistent with the nominal order in (2.22), this is not true of the hydrodynamic component. As shown in § 4, the hydrodynamic component scales as $O(Ri_v)$ only for sufficiently small $Pe (\ll Ri_v^{{3}/{5}}$ for $Ri_v \ll 1$; see § 4.1.1), when the dominant contribution to the associated torque integral comes from length scales of $O(L)$ similar to the inertial torque above. In the opposite limit, and for the so-called Stokes stratification regime corresponding to $Re \ll Ri_v^{{1}/{3}}$ (see Mehaddi et al. Reference Mehaddi, Candelier and Mehlig2018; Varanasi & Subramanian Reference Varanasi and Subramanian2021), the dominant contributions to the torque integral arise from much larger length scales of $O(LRi_v^{-{1}/{3}})$, and the hydrodynamic component scales as $O(Ri_v^{{2}/{3}})$, being much larger than the hydrostatic component above.

Before proceeding with the calculation of the hydrodynamic component of the stratification torque, it is worth remarking on the nature of the coupling between the inertial and stratification torque contributions that is not obvious from the formal result (2.22), where they appear as separate additive contributions. On account of the convergent volume integral, the $O(Re)$ inertial angular velocity, as given by (3.3) and (3.4), only involves the Stokesian fields in a homogeneous ambient, and is evidently independent of the ambient stratification. The correction to this leading-order estimate is dependent on the nature of the ambient stratification, however, even within the Boussinesq framework. To see this, we return to the inertial volume integral, and estimate the next correction. Recall that the angular velocities in (3.3) and (3.4) were based on the approximating the volume integral by Stokesian estimates, and the torque contribution at the next order requires one to examine the next term in the small-$Re$ expansion for the velocity field in problem 1. Writing ${\boldsymbol u}'^{(1)} = {\boldsymbol u}'^{(10)}+ Re\,{\boldsymbol u}'^{(11)}$, ${\boldsymbol u}'^{(10)}$ is the Stokesian approximation given by (3.1) and is $O(1/r)$ for $r \gg 1$, whereas ${\boldsymbol u}'^{(11)}$ remains $O(1)$ in the far-field. The latter, of course, implies that the above regular expansion breaks down at length scales of $O(LRe^{-1})$, a manifestation of the singular nature of inertia in an unbounded domain (the so-called Whitehead's paradox; see Leal Reference Leal1992). Provided one assumes buoyancy forces to dominate the inertial ones on scales much smaller than the inertial screening length (of $O(LRe^{-1}))$, the above far-field estimate of ${\boldsymbol u}'^{(11)}$ may still be used to estimate the correction to the leading $O(Re)$ contribution. The $O(Re^2)$ inertial acceleration is now $\hat {\boldsymbol U} \boldsymbol {\cdot } \boldsymbol {\nabla } {\boldsymbol u}'^{(11)} \sim O(1/r)$, and using ${\boldsymbol u}^{(2)} \sim O(1/r^2)$ , the resulting volume integral, at $O(Re^2)$, is logarithmically divergent. The divergence will be cut off at the stratification screening length that is $O(Ri_vPe)^{-{1}/{4}}$ for $Pe \ll 1$ (List Reference List1971; Ardekani & Stocker Reference Ardekani and Stocker2010), and $O(Ri_v^{-{1}/{3}})$ for $Pe \gg 1$ (Mehaddi et al. Reference Mehaddi, Candelier and Mehlig2018; Varanasi & Subramanian Reference Varanasi and Subramanian2021), implying that the next correction to the inertial angular velocity is $O[Re^2\ln (Ri_v Pe)^{-{1}/{4}}]$ for $Pe \ll 1$ and $O(Re^2\ln Ri_v^{-{1}/{3}})$ for $Pe \gg 1$, and is thereby a function of the ambient stratification. For self consistency, one requires that both of the aforementioned stratification screening lengths be less than $O(Re^{-1})$, which translates to the requirement $Ri_v \gg Pe^{-1} Re^4$ for small $Pe$, and for one to be in the aforementioned Stokes stratification regime ($Ri_v^{{1}/{3}} \gg Re$) for large $Pe$.

4. The spheroidal angular velocity due to the hydrodynamic component of the stratification contribution

Owing to the differing character of the hydrodynamic component in the limits of small and large $Pe$, the calculations in these two asymptotic regimes are carried out in separate subsections here. Keeping in mind that $Pe = RePr$, $Pr$ being the Prandtl number, the small $Pe$ case does not necessarily place a restriction on $Pr$ which may either be small or large (although, large $Pr$ imposes a greater restriction on the smallness of $Re$ because $Re$ must now be smaller than $O(Pr^{-1})$). However, the assumption of small $Re$ implies that the large $Pe$ case necessarily requires a large $Pr$ which may be realized in experiments that use salt as a stratifying agent.

4.1. The hydrodynamic stratification torque in the diffusion-dominant limit ($Pe \ll 1$)

In the limit $Pe \ll 1$, one may neglect the convective terms in the advection diffusion equation (2.16), and the density perturbation $\rho '^{(1)}$ in the stratification torque integral therefore arises as a diffusive response to the no-flux condition that must be satisfied on the spheroid surface. As a result, the spheroid acts as a concentration-dipole singularity in the far-field ($r \gg 1$), implying that $\rho '^{(1)}$ must decay as $O(1/r^2)$. As the test velocity field ${\boldsymbol u}^{(2)}$ corresponding to the rotating spheroid also decays as $O(1/r^2)$, the integrand is $O(1/r^4)$ for $r \gg 1$. This decay is the same as that of the inertial integrand estimated above, and the integral for the hydrodynamic stratification torque is therefore convergent for small $Pe$, based on the leading-order diffusive estimates above. Thus, the effects of stratification arise as a regular perturbation for small $Pe$, or said differently, the dominant contribution to the hydrodynamic component of the stratification torque arises from buoyancy forces in a volume of $O(L^3)$ around the sedimenting spheroid. As a consequence and as is shown in the following, for small $Pe$, the hydrodynamic component is $O(Ri_v)$ similar to the hydrostatic component given in (3.6) and (3.7). It turns out that there is, in fact, a more severe constraint on the Péclet number; as explained in § 4.1.1, the $O(Ri_v)$ scaling for the hydrodynamic component holds only when $Pe \ll Ri_v^{{3}/{5}} \ll 1$.

To determine the detailed dependence of the $O(Ri_v)$ hydrodynamic component on $\kappa$, one needs to solve for the density perturbation $\rho '^{(1)}$ which satisfies

(4.1)\begin{equation} \nabla^2 \rho'^{(1)} = 0 \end{equation}

in the fluid domain $\xi \geq \xi _0$. The no-flux boundary condition on the spheroid surface ($\xi = \xi _0$) may be written as ${\boldsymbol 1}_\xi \boldsymbol {\cdot } {\boldsymbol \nabla } \rho '^{(1)} = -{\boldsymbol 1}_\xi \boldsymbol {\cdot } \hat {\boldsymbol g}$ where we have used that ${\boldsymbol n} = {\boldsymbol 1}_\xi$, and the right-hand side of the boundary condition arises from the gradient of the linearly varying ambient density; there is the additional requirement of far-field decay, namely $\rho '^{(1)} \rightarrow 0$ for $\xi \rightarrow \infty$. Note that the linearity of the governing equation (4.1) and the boundary conditions in $\rho '^{(1)}$, and the linear dependence on $\hat {\boldsymbol g}$ of the surface boundary condition above, imply that $\rho '^{(1)}$ must be linear in $\hat {\boldsymbol g}$ at leading order for small $Pe$. From (2.22), the hydrodynamic component of the stratification torque must therefore be quadratic in $\hat {\boldsymbol g}$. This implies that the hydrodynamic torque must have an angular dependence identical to the inertial and hydrostatic contributions, of the form $(\hat {\boldsymbol g} \wedge {\boldsymbol p})(\hat {\boldsymbol g} \boldsymbol {\cdot } {\boldsymbol p}$), with a multiplicative pre-factor that is a function of $\kappa$. Thus, for small $Pe$, the ratio of the hydrostatic and hydrodynamic components of the stratification torque is independent of the spheroid orientation and $Ri_v$, and only a function of $\kappa$.

To solve (4.1), we write down the explicit form of the no-flux boundary condition, in prolate spheroidal coordinates:

(4.2)\begin{equation} \left.\frac{\partial \rho'^{(1)}}{\partial \xi}\right|_{\xi = \xi_0} =\frac{(P_1^1(\eta) \textrm{e}^{\mathrm{i}\phi} - 2 P_1^{{-}1}(\eta)\textrm{e}^{\mathrm{i}\phi} )}{2\bar{\xi}_0} \sin \psi - \frac{P^0_1(\eta)}{\xi_0} \cos \psi, \end{equation}

where $\psi$ denotes the angle between ${\boldsymbol p}$ and ${\boldsymbol g}$; in (4.2), and $P_t^s(\xi )$ is the associated Legendre function of the first kind with $P_1^{-1}(\eta ) =-P_1^1(\eta ) = -({\bar {\eta }}/{2})$ and $P_1^0({\eta })=\eta$. The $\eta$-dependence of the solution is imposed by that of the boundary condition above. This, and the fact that the Laplacian is separable in spheroidal coordinates with the eigenfunctions in $\xi$ and $\eta$ being the associated Legendre functions of the second and the first kind, respectively, points to the following ansatz for $\rho '^{(1)}$:

(4.3)\begin{equation} \rho'^{(1)} = A_{1,1}Q_1^1(\xi)P_1^1(\eta) \textrm{e}^{\mathrm{i}\phi} + A_{1,-1} Q_1^{{-}1}(\xi)P_1^{{-}1}(\eta)\textrm{e}^{-\mathrm{i}\phi} + A_{1,0}Q_1^0(\xi)P_1^0(\eta), \end{equation}

where $Q_1^{-1}(\xi ) =Q_1^1(\xi )/2$ , the latter being defined in Appendix C (see (C7)). Substitution of (4.3) leads to the following expressions for the $A_{i,j}$ in (4.3):

(4.4)\begin{gather} A_{1,1} = \frac{\bar{\xi}_0^2}{4-2\xi_0^2 + 2\xi_0\bar{\xi}_0^2\coth^{{-}1}\xi_0}\sin \psi, \end{gather}

(4.5)\begin{gather}A_{1,-1} ={-}4A_{1,1}, \end{gather}

(4.6)\begin{gather}A_{1,0} =\frac{\bar{\xi}_0^2}{\xi_0(\xi_0 - \bar{\xi}_0^2\coth^{{-}1}\xi_0)}\cos \psi. \end{gather}

The first two terms in (4.3) correspond to the density perturbation induced by a spheroid oriented transversely to gravity, whereas the third term corresponds to that induced by a spheroid aligned with gravity. The expressions for the $A_{i,j}$ above are therefore consistent with the underlying linearity of the problem, in that the perturbation induced by an arbitrarily oriented spheroid is obtained as a superposition of the transverse and longitudinal problems, the factors involved in this superposition being $\sin \psi$ and $\cos \psi$, respectively.

The density perturbation from (4.3), together with (4.4)–(4.6), and the test velocity field from (3.2) are now substituted in the integral for the hydrodynamic stratification torque in (2.22). The volume integration is carried out in the spheroidal coordinate system introduced in § 3, and leads to the following expression for the angular velocity:

(4.7)\begin{equation} \varOmega_i^{(1)d} =Ri_v \frac{F^{p/o}_s(\xi_0)}{Y_C}\, (\epsilon_{ijk} \hat{g}_j p_k) (\hat{g}_lp_l), \end{equation}

with

(4.8)\begin{align} &F^p_s(\xi_0)=2{\rm \pi} {\overline{\xi_0}}^2\nonumber\\ &{ {\frac{ \left(\left(7 \xi_0^5-7 \xi_0^3-2\xi_0\right)+\overline{\xi_0}^2 C_p \left({-}12 \xi_0^4+6 \xi_0^2+\xi_0 C_p \left(3 \xi_0^4+2 \left(\xi_0^5-\xi_0 \right) C_p -6 \xi_0^2+7\right)-2\right)\right)}{15 \xi_0^5 \left(\xi_0-\overline{\xi_0}^2 C_p\right) \left(-\xi_0^2+\overline{\xi_0}^2 \xi_0 C_p+2\right) \left(\left(\xi_0^2+1\right) C_p-\xi_0\right)}}}, \end{align}

for prolate spheroids, where $C_p=\coth ^{-1}\xi _0$. Using the transformation mentioned in § 3, one obtains

(4.9)\begin{align} &F^o_s(\xi_0)=2{\rm \pi} \nonumber\\ &{\frac{ \left(\left(7 \xi_0^4-7 \xi_0^2-2\right) \bar{\xi}_0+2 \left(\xi_0^4-3 \xi_0^2+2\right) \xi_0^4 C_o^3-2 \left(6 \xi_0^4-9 \xi_0^2+4\right) \xi_0^2 C_o+\left(3 \xi_0^4+4\right) \bar{\xi}_0 \xi_0^2 C_o^2\right)}{15 \xi_0^3 \left(\xi_0^2 C_o)-\bar{\xi}_0\right) \left(\left(\xi_0^2-2\right) \left(C_o\right)-\bar{\xi}_0\right) \left(-\xi_0^2+\bar{\xi}_0 \xi_0^2 C_o-1\right)}}, \end{align}

for the oblate case, where $C_o=\cot ^{-1}\bar {\xi }_0$. Figure 3 shows a plot of the aspect-ratio- dependent functions given by (4.8) and (4.9), divided by $Y_C$. The hydrodynamic component of the stratification torque given by (4.8) always orients a prolate spheroid edgewise for small $Pe$. Interestingly, figure 3 shows that (4.9) changes sign below a critical aspect ratio $\kappa _c \approx 0.41\ (e \approx 0.9)$ and, therefore, the hydrodynamic component acts to orient oblate spheroids, with aspect ratios lower than the aforementioned threshold, broadside-on.

Figure 3.

The aspect-ratio-dependent functions, given by (4.8) and (4.9) divided by $Y_C$, that characterize the hydrodynamic stratification-induced contribution for prolate and oblate spheroids in the small $Pe$ limit, plotted as a function of the spheroid eccentricity.

The hydrodynamic stratification torque arises due to the flow associated with the baroclinic source of vorticity, although the reciprocal theorem formulation used here bypasses the explicit calculation of this flow. The vorticity arises from the tilting of the isopycnals to the vertical (the direction of gravity) due to the requirement that they meet the spheroidal surface in a normal orientation, consistent with a no-flux constraint. A sketch of the deformed isopycnals in the plane $\phi = 0,{\rm \pi}$, and the resulting sign of the baroclinic vorticity field ($\propto \boldsymbol {\nabla } \rho \wedge \hat {\boldsymbol g}$) in different regions of the fluid domain, appears in figure 4, for both prolate and oblate spheroids. The baroclinically induced flow has a dipolar character in the Stokes limit, and the relative sizes of the different flow quadrants are set by the pair of singular isopycnals that meet the spheroid surface in the points $S_1$ and $S_2$. This pair separates the isopycnals that do not meet the spheroid surface from those that do. The baroclinically induced flow on account of diffusive isopycnal tilting has been known for a long time, having originally been proposed in the oceanic context where the induced flow has a boundary layer character on account of the dominance of inertia (Phillips Reference Phillips1970; Wunsch Reference Wunsch1970). Such a flow has also been examined in a Stokesian scenario more recently (Anis Alias & Page Reference Anis Alias and Page2017), although only for the case of a horizontal circular cylinder wherein symmetry precludes a torque contribution. That there must be a torque on an inclined spheroid, due to the aforementioned baroclinic flow, is obvious. Although the sense of the torque (broadside-on versus edgewise) is not readily evident, one may nevertheless rationalize the scalings observed for the extreme aspect ratio cases. Figure 3 shows that the angular velocity remains finite in the limit of a flat disk ($\kappa \rightarrow 0$) which is consistent with the isopycnals being perturbed in a volume of $O(L^3)$ around the spheroid in this limit, with the density perturbation being $O(\gamma L)$, and the test velocity field ${\boldsymbol U}^{(2)}$ being $O(1)$ in this region; the points $S_1$ and $S_2$ remain bounded away from the edges of the flat disk. On the other hand, the angular velocity approaches zero as $O(\xi _0-1)$ in the limit of a slender fiber, with the points $S_1$ and $S_2$ now moving towards the ends of the fiber. The dominant contribution continues to come from a volume of $O(L^3)$. However, although ${\boldsymbol U}^{(2)}$ is $O[\ln (\xi _0-1)]^{-1}$, the density perturbation in this region is algebraically small. The slender fiber only perturbs the isopycnals in a thin $O(d^2L)$ shell around itself, with the density perturbation being $O(\gamma d)$ in this region. Further, each cross-section of the fiber acts as a 2D concentration dipole, implying that the density perturbation decays as $O(1/r)$ for $r$ much greater than $d$, and is therefore $O(\gamma d)(d/L)$ for $r \sim O(L)$. Using these estimates, and dividing by the $O[\ln (\xi _0-1)]^{-1}$ rotational resistance for a slender fiber leads to the aforementioned scaling for the fiber rotation due to the hydrodynamic contribution of the stratification torque.

Figure 4.

Baroclinically driven flow, for small $Pe$, that is responsible for the rotation of an (a) oblate and a (b) prolate spheroid, in a stably stratified ambient. The curved blue arrows denote the sense of the baroclinically induced vorticity in the different quadrants of the fluid domain, with vorticities corresponding to anticlockwise and clockwise senses of rotation being denoted by solid and dashed lines; the blue contours denote the deformed isopycnals around each of the spheroids.

4.1.1. A closer look at the $Pe\ll 1$ analysis

Although the $O(Ri_v)$ angular velocity given by (4.7) was said to be valid for small $Pe$, there are, in fact, multiple contributions to the stratification-induced rotation for $Pe \ll 1$; the more detailed arguments herein, and the analysis in Appendix B, yield a precise estimate of the $Pe$ interval of validity. Although the dominant contribution to (4.7) arose from buoyancy forces on length scales of $O(L)$, as we show in § 4.2, the dominant length scale contributing to the hydrodynamic component of the stratification torque changes to $O(L Ri_v^{-{1}/{3}})$ with increasing $Pe$; one therefore expects the singular effect of convection to already be evident for small but finite $Pe$. To see this, we note that, for small $Pe$, in addition to the density perturbation driven by the no-flux condition on the spheroid surface that has been analysed in section (4.1), an independent contribution arises from perturbation of the ambient stratification on much larger length scales due to weak convection effects. To obtain an estimate for this latter torque contribution, we consider the correction to the leading-order density perturbation, now denoted as $\rho '^{(10)}$, in a manner similar to the velocity field examined at the end of the earlier section; thus, one writes $\rho '^{(1)} = \rho '^{(10)} + Pe\,\rho '^{(11)}$, with $\nabla ^2 \rho '^{(11)} \sim {\boldsymbol u}'^{(1)} \boldsymbol {\cdot } \hat {\boldsymbol g}$. With ${\boldsymbol u}'^{(1)} \sim 1/r$ on length scales smaller than $O(LRe^{-1})$, one obtains $\rho '^{(11)} \sim r$. The stratification torque at the next order is proportional to $Ri_v Pe\int \rho '^{(11)} \hat {\boldsymbol g} {\boldsymbol u^{(2)}}\,\textrm {d}V$, which turns out to diverge as $O(r^2)$, on using the above estimate for $\rho '^{(11)}$. Cutting off the divergence at the small-$Pe$ stratification screening length of $O[L(Ri_vPe)^{-{1}/{4}}]$ (List Reference List1971; Ardekani & Stocker Reference Ardekani and Stocker2010) would seem to lead to an $O(Ri_vPe)^{{1}/{2}}$ torque contribution. However, it is shown in Appendix B that this contribution is identically zero, on account of the fore–aft symmetry of the disturbance density field on scales of $[L(Ri_vPe)^{-{1}/{4}}]$ (Ardekani & Stocker Reference Ardekani and Stocker2010; Varanasi & Subramanian Reference Varanasi and Subramanian2021).

The fore–aft asymmetry of the density perturbation, necessary for a non-trivial torque contribution, requires inclusion of the $O(Pe)$ convective term (${\boldsymbol u} \boldsymbol {\cdot } \boldsymbol {\nabla } \rho '$) in (2.12). It is well known that, for small but finite $Pe$, this convective term becomes comparable to the diffusive term on length scales of $O(LPe^{-1})$, the mass/heat transfer analogue of the inertial screening length (Leal Reference Leal1992), and the angular velocity scaling therefore depends on the relative magnitudes of the convective ($LPe^{-1}$) and stratification ($L(Ri_vPe)^{-{1}/{4}}$) screening lengths, which in turn depends on the relative magnitudes of $Pe$ and $Ri_v^{{1}/{3}}$; this criterion being analogous to the classification into Stokes ($Re \ll Ri_v^{{1}/{3}}$) and inertia ($Re \geq Ri_v^{{1}/{3}}$) stratification regimes based on the structure of the large-$Pe$ disturbance velocity field (Mehaddi et al. Reference Mehaddi, Candelier and Mehlig2018; Varanasi & Subramanian Reference Varanasi and Subramanian2021), except that $Pe$ now replaces $Re$. For $Pe \ll Ri_v^{{1}/{3}}$, fore–aft asymmetric buoyancy forces acting on scales of $O[L(Ri_vPe)^{-{1}/{4}}]$ lead to an $O(Ri_v^{{1}/{4}}Pe^{{5}/{4}})$ hydrodynamic torque contribution, an exact expression for which is obtained in Appendix B. In the opposite limit of $Pe \gg Ri_v^{{1}/{3}}$, the dominant contribution to the torque integral arises on scales of $O(Ri_v^{-{1}/{3}})$, and the resulting torque comes out to be $O(Ri_v^{{2}/{3}})$. Thus, the arguments above, and those in Appendix B, show that, for small $Pe$, in addition to the $O(Ri_v)$ torque contribution given by (4.8) and (4.9), there exists a second independent contribution that increases with $Pe$ as $O(Ri_v^{{1}/{4}}Pe^{{5}/{4}})$ for $Pe \ll Ri_v^{{1}/{3}}$, but is independent of $Pe$, being $O(Ri_v^{{2}/{3}})$ for $Ri_v^{{1}/{3}} \ll Pe \ll 1$. This far-field hydrodynamic contribution, arising from a weak convective distortion of the stratified ambient, can evidently exceed the hydrostatic contribution, possibly leading to an edgewise settling regime even for small $Pe$. In light of this additional contribution, the dominance of the $O(Ri_v)$ hydrostatic contribution and the prevalence of broadside-on settling requires $Ri_v \gg Ri_v^{{1}/{4}}Pe^{{5}/{4}}$, which translates to the stricter criterion $Pe \ll Ri_v^{{3}/{5}} \ll 1$, instead of $Pe \ll 1$, as assumed originally.

It is important to point out that the $O(Ri_v)$ scaling in (3.6), (3.7) and (4.7) implies that the associated (dimensional) angular velocities are independent of $U$. Although this must be the case for the hydrostatic contributions, it turns out to be the case for the hydrodynamic component too, at small $Pe$, because the leading-order density perturbation arises as a diffusive response, and is therefore independent of the ambient flow. Despite this $U$-independence, the torque associated with (4.7) does have a hydrodynamic character, in that it still arises from the flow induced by buoyancy forces. Using the $O(\gamma L)$ density perturbation produced by the diffusive response, one obtains a buoyancy-driven velocity scale of $O(\gamma L^3 g/\mu )$, implying a spheroidal angular velocity of $O(\gamma L^2 g/\mu )$; the latter is $O(Ri_v)$ in units of $U/L$, this being the scale used in the reciprocal theorem formulation in § 2. Importantly, the $U$-independence implies that this torque contribution is not necessarily limited to small $Ri_v$. Instead, it is limited by the assumption of a quasi-steady density perturbation set up on scales of $O(L)$ by diffusion alone, which requires an appropriate Péclet number (based on the aforementioned $O(\gamma L^3 g/\mu )$ velocity scale) to be small; an analogously defined Reynolds number must also be small for the baroclinic flow to be obtained from the Stokes equations. These two requirements translate to $Ri_vPe, Ri_v Re \ll 1$, which also ensures that the spheroid rotation may be neglected when deriving the disturbance fields. Thus, (4.7) remains valid even when $Ri_v \sim O(1)$ provided $Re, Pe \ll 1$. Note, however, that a genuine dependence, of the hydrodynamic component of the angular velocity contribution, on $U$ arises due to contributions from the outer region, and an estimate of this contribution was obtained in the preceding paragraph (also see Appendix B).

To summarize then, for sufficiently small $Pe$, all three contributions in (2.22) have a regular character, and therefore, the same dependence on the spheroid orientation, namely $\sin \psi \cos \psi$ with $\psi$ being the angle between ${\boldsymbol p}$ and ${\boldsymbol g}$ as defined above; the $Pe$ interval of validity depends on $Ri_v$, being $Pe \ll 1$ for $Ri_v \sim O(1)$, and $Pe \ll Ri_v^{{3}/{5}}$ for $Ri_v \ll 1$. The inertial contribution is $O(Re)$, whereas both hydrodynamic and hydrostatic components of the stratification contribution are $O(Ri_v)$, with the hydrodynamic component alone acting to orient the spheroid edgewise for prolate spheroids of arbitrary $\kappa$ and oblate spheroids with $\kappa > 0.41$. Therefore, oblate spheroids with $\kappa < 0.41$ will certainly orient broadside-on for $Pe \ll 1$. Further, it is seen from figures 2 and 3 that the hydrodynamic component always remains smaller in magnitude than the hydrostatic one in the edgewise-rotation regime and, therefore, a sedimenting spheroid, either prolate or oblate, is expected to settle in the broadside-on configuration, regardless of $\kappa$, for sufficiently small $Pe$.

4.2. The angular velocity due to the hydrodynamic torque in the convection-dominant limit ($Pe \gg 1$)

In contrast to the small $Pe$ limit examined in the previous section, for $Pe \gg 1$, the dominant contribution to the integral for the stratification-induced hydrodynamic torque in (2.21) comes from length scales of $O(LRi_v^{-{1}/{3}})$, the stratification screening length in the Stokes stratification regime ($Re \ll Ri_v^{{1}/{3}}\ll 1$). To see this, we note from the right-hand side of the advection diffusion equation (2.16) that, for large $Pe$, the density perturbation is driven by the convection of the base-state stratification of order unity by the vertical component of the Stokeslet field, $u'^{(1)}_3$. As $u'^{(1)}_3 \sim O(1/r)$, one finds $\rho '^{(1)} \sim O(1)$ for $r\gg 1$. This, along with the far-field $O(1/r^2)$ decay of the test velocity field, implies that the integrand in the stratification torque integral decays as $O(1/r^2)$, and that the volume integral therefore diverges as $O(r)$. This divergence is expected to be resolved only when the slow $O(1/r)$ decay of the Stokeslet is accelerated by stratification that, for large $Pe$, occurs on length scales of $O(LRi_v^{-{1}/{3}})$. It has recently been shown in Varanasi & Subramanian (Reference Varanasi and Subramanian2021) that, for a sedimenting sphere at large $Pe$, the density and velocity fields are indeed asymptotically small on length scales larger than $O(Ri_v^{-{1}/{3}})$, except within a horizontal wake whose vertical extent grows as $O(r_t^{{2}/{5}})$, $r_t$ being the distance in the plane transverse to gravity, where the density and axial velocity perturbation decay as $r_t^{-{12}/{5}}$ and $r_t^{-{14}/{5}}$, respectively; and a buoyant jet in the rear where $u'^{(1)}_3$ reverses sign, but continues to exhibit an $O(1/r)$ Stokesian decay. Despite the latter slow decay, the asymptotically narrow character of the buoyant jet implies that the torque integral does converge on length scales of $O(LRi_v^{-{1}/{3}})$, and is $O(Ri_v^{-{1}/{3}})$. The pre-factor of $Ri_v$ in front of the integral in (2.21) implies that the torque and the angular velocity scale as $O(Ri_v^{{2}/{3}})$ for $Pe \rightarrow \infty$.

As the dominant contribution to the torque integral comes from length scales much larger than $O(L)$, of $O(LRi_v^{-{1}/{3}})$, the calculation requires one to rewrite the integral involving $\rho '^{(1)}$, in (2.21), in outer coordinates (defined below), with the sedimenting spheroid in problem 1 now regarded as a point force, and the rotating spheroid in the test problem acting as a combination of rotlet and stresslet singularities (see Marath & Subramanian Reference Marath and Subramanian2017). The details of this calculation are provided below; a similar calculation, but for low $Pe$, with the outer region being characterized by a length scale of $O(Ri_vPe)^{-{1}/{4}}$, has been given in Appendix B and the resulting $O(Ri_v^{{1}/{4}} Pe^{{5}/{4}}$) torque contribution was discussed above in § 4.1.1.

Before delving into the large-$Pe$ analysis, it is worth noting that a reciprocal theorem formulation to determine the stratification-induced correction to the force (that would include both drag and lift components for an arbitrarily oriented spheroid) would involve the test problem of a translating spheroid instead. As the test velocity field now decays as $O(1/r)$ in the far-field, this would lead to a stronger $O(r^2)$ divergence of the force integral, leading to a scaling of $O(Ri_v^{-{2}/{3}})$ on truncation of the divergence and, thence, an $O(Ri_v^{{1}/{3}})$ stratification-induced correction to the Stokes drag. Such a correction was originally calculated for a spherical particle by Zvirin & Chadwick (Reference Zvirin and Chadwick1975). A similar calculation was done by Candelier et al. (Reference Candelier, Mehaddi and Vauquelin2014) for the small-$Pe$ regime, the drag contribution from the outer region now being $O(Ri_v Pe)^{{1}/{4}}$; a later effort has connected this outer-region drag calculation across the different asymptotic regimes (Mehaddi et al. Reference Mehaddi, Candelier and Mehlig2018). In a very recent study, Dandekar et al. (Reference Dandekar, Shaik and Ardekani2020) have examined the force and torque acting on an arbitrarily shaped particle sedimenting in a linearly stratified ambient. For anisotropic particles lacking a handedness (that includes the spheroids examined here), the authors find a correction to the force at $O(Ri_v^{{1}/{3}})$ similar to the case of a spherical particle mentioned above efforts, but end up not finding a torque at this order, a result that is not surprising in light of the above scaling arguments which show the torque to be $O(Ri_v^{{2}/{3}})$. Within the framework of the matched asymptotics expansions approach used by the said authors, the $O(Ri_v^{{1}/{3}})$ correction to the drag appears as a response of the particle to an ‘ambient uniform flow’ that is the limiting form of the outer solution in the matching region ($1 \ll r \ll Ri_v^{-{1}/{3}}$); the uniformity of this flow is consistent with the absence of a torque at this order. Incidentally, the existence of an inertial torque induced by a uniform flow suggests a higher $O(Re Ri_v^{{1}/{3}})$ inner-region contribution to the leading $O(Re)$ inertial angular velocity, in response to the aforementioned $O(Ri_v^{{1}/{3}})$ uniform flow that would again orient the spheroid broadside-on; this is in addition to the outer-region correction obtained in § 3, and is dominant in the Stokes stratification regime.

As mentioned above, the stratification torque integral in (2.22) needs to be evaluated in outer coordinates which are related to the coordinates in the particle-fixed reference frame as $\tilde {\boldsymbol x} = Ri_v^{{1}/{3}}{\boldsymbol x}$, so an $O(1)$ change in $\tilde {\boldsymbol x}$ corresponds to ${\boldsymbol x}$ changing by an amount of order the stratification screening length. However, as originally shown by Childress (Reference Childress1964) and Saffman (Reference Saffman1965), a Fourier space approach turns out to be much more convenient for a calculation involving the outer region, and we therefore consider the Fourier transformed equations of continuity and motion, and the advection diffusion equation for the density field, obtained from (2.14)–(2.16), and given by

(4.10)\begin{gather} k_i \hat{u}'^{(1)}_i = 0, \end{gather}

(4.11)\begin{gather}-4{\rm \pi}^2 k^2\hat{u}'^{(1)}_i - 2{\rm \pi} \mathrm{i}k_i \hat{p}'^{(1)} ={-} Ri_v(\hat{\rho}'^{(1)} \hat{g}_i )+ \tilde{F}\hat{g}_i, \end{gather}

(4.12)\begin{gather}2{\rm \pi} \mathrm{i}k_j \hat{U}_j \hat{\rho}'^{(1)} = \hat{u}'^{(1)}_j\hat{g}_j + \frac{1}{Pe}4{\rm \pi}^2 k^2 \hat{\rho}'^{(1)} + \frac{1}{Pe} \textrm{i} 2 {\rm \pi}k_i D^s_i , \end{gather}

for problem 1. Here, we have used the definition $\hat {f}({\boldsymbol k}) = \int \,\textrm {d}\kern0.06em {\boldsymbol x} \exp ({-2{\rm \pi} \mathrm {i} {\boldsymbol k}\boldsymbol {\cdot }{\boldsymbol x}}) f({\boldsymbol x})$ for the Fourier transform, and the sedimenting spheroid has been replaced by a point force, $\tilde {\boldsymbol F}\delta ({\boldsymbol x})$, on the right-hand side of the physical space equations of motion, namely (2.15), with $\tilde {\boldsymbol F} = \tilde {F}\hat {\boldsymbol g}$, $\tilde {F}$ being the non-dimensional buoyant force (in units of $\mu UL$) exerted by the spheroid; the corresponding dimensional expression has been given in § 2 (given that $U$ is itself defined in terms of $F$, this re-scaling of $F$ merely amounts to multiplication by $X_A$). Note that the inertial term in the original equation of motion, (2.15), has now been omitted in (4.11) because, as argued earlier, the leading $O(Re)$ inertial torque is dominated by the inner region, with the outer region contribution being only $O(Re^2\ln Ri_v^{-{1}/{3}})$ (see the end of § 3), and not considered here. Further, on large length scales relevant to the outer region, the spheroid, on account of the no-flux boundary condition at its surface, appears as a concentration-dipole forcing in the advection–diffusion equation. In physical space, this corresponds to a term of the form $\boldsymbol {D}^s\boldsymbol {\cdot }\boldsymbol {\nabla } \delta ({\boldsymbol x})$, and the Fourier transform of this term appears on the right-hand side of (4.12); $\boldsymbol {D}^{s}$ here is the $Pe$-dependent strength of the dipole forcing; this can be neglected for the same reason as the inertial term above, its contribution being asymptotically small compared with that arising from the distortion of the base-state stratification (the first term on the right-hand side of (4.12)).

Setting $Pe =\infty$ in (4.12), one obtains

(4.13)\begin{equation} \hat{\rho}'^{(1)} = \frac{\hat{u}'^{(1)}_j \hat{g}_j}{2{\rm \pi} \mathrm{i} (k_l \hat{U}_l)}. \end{equation}

Using (4.13) in (4.11), and operating on both sides with $(\delta _{ij} - \hat {k}_i\hat {k}_j)$ to eliminate the pressure field, one obtains

(4.14)\begin{equation} 4{\rm \pi}^2 k^2\hat{u}'^{(1)}_i = Ri_v \hat{g}_j(\delta_{ij} - \hat{k}_i\hat{k}_j)\frac{\hat{u}'^{(1)}_l \hat{g}_l}{2{\rm \pi} \mathrm{i} (k_p \hat{U}_p)} + \tilde{F} \hat{g}_j (\delta_{ij} - \hat{k}_i \hat{k}_j). \end{equation}

Contracting with $\hat {\boldsymbol g}$ gives

(4.15)\begin{equation} \hat{u}'^{(1)}_i\hat{g}_i =\frac{\tilde{F}\hat{g}_i\hat{g}_j(\delta_{ij} - \hat{k}_i\hat{k}_j)}{\left[ 4{\rm \pi}^2 k^2 - Ri_v\dfrac{\hat{g}_m \hat{g}_n(\delta_{mn} - \hat{k}_m\hat{k}_n)}{2{\rm \pi} \mathrm{i} (k_l \hat{U}_l)}\right]}, \end{equation}

which, on using in (4.13), yields the following final expression for $\hat {\rho }'^{(1)}$:

(4.16)\begin{equation} \hat{\rho}'^{(1)}({\boldsymbol k}) = \frac{\tilde{F}[1-(\hat{k}_i \hat{g}_i)^2]}{\{8{\rm \pi}^3\mathrm{i}k^3(\hat{k}_n \hat{U}_n) - Ri_v[1-(\hat{k}_m \hat{g}_m)^2]\}}, \end{equation}

which will be used in the Fourier-space torque integral that is defined below.

The test velocity field $u^{(2)}_i =U^{(2)}_{ij}\varOmega ^{(2)}_j$ satisfies the Stokes equations, with the rotating spheroid, on length scales of $O(LRi_v^{-{1}/{3}})$, acting as a force-dipole singularity that includes both stresslet and rotlet contributions. Thus, the equations of motion may be written in the form (see Marath & Subramanian Reference Marath and Subramanian2017):

(4.17)\begin{equation} \frac{\partial ^2 u^{(2)}_i}{\partial x_j^2} - \frac{\partial p^{(2)}}{\partial x_i} = S^{(2)}_{ij} \frac{\partial }{\partial x_j} [\delta({\boldsymbol x})], \end{equation}

where

(4.18)\begin{equation} S^{(2)}_{ij} = B_1[(\epsilon_{ilm}\varOmega^{(2)}_lp_m )p_j + (\epsilon_{jlm}\varOmega^{(2)}_lp_m )p_i] + B_3 \epsilon_{ijk}\varOmega^{(2)}_k \end{equation}

with

(4.19)\begin{gather} B_1 = \frac{8{\rm \pi}}{\xi_0^3({-}3\xi_0 + 3 \coth^{{-}1}\xi_0(1+\xi_0^2))}, \end{gather}

(4.20)\begin{gather}B_3 =\frac{8{\rm \pi}(1-2\xi_0^2)}{\xi_0^3({-}3\xi_0 + 3 \coth^{{-}1}\xi_0(1+\xi_0^2))}, \end{gather}

for prolate spheroids. There is an additional contribution that is neglected in (4.18), on account of the test spheroid rotating about an axis transverse to ${\boldsymbol p}$, that is, because ${\boldsymbol \varOmega }^{(2)} \boldsymbol {\cdot } {\boldsymbol p}=0$. The term proportional to $B_3$ in (4.18) is the rotlet singularity (due to transverse rotation), whereas that involving $B_1$ is the stresslet singularity. Thus, for $\xi _0 \rightarrow \infty$, $B_3 = -4{\rm \pi}$ and $B_1$ is $O(1/\xi _0^2)$, consistent with a rotating sphere acting as a pure rotlet singularity; note that $B_3 =Y_C/2$, the latter being the resistance function mediating the torque–angular velocity relation for transverse rotation defined earlier.

Fourier transforming (4.17), and contracting with the projection operator $({\boldsymbol I} -\hat {\boldsymbol k}\hat {\boldsymbol k})$, one obtains

(4.21)\begin{align} \hat{u}^{(2)}_i ={-} \frac{\mathrm{i}}{2{\rm \pi} k}\{ B_1[(\epsilon_{mqr}\varOmega^{(2)}_qp_r )p_n + (\epsilon_{nqr}\varOmega^{(2)}_qp_r )p_m] + B_3 \epsilon_{mnq}\varOmega^{(2)}_q\}\hat{k}_n(\delta_{im} - \hat{k}_i\hat{k}_m), \end{align}

so the second-order tensor ${\boldsymbol U}^{(2)}$ is given by

(4.22)\begin{equation} U^{(2)}_{ij}({\boldsymbol k}) ={-}\frac{\mathrm{i}}{2{\rm \pi} k}\{ B_1[(\epsilon_{mjr} p_r )p_n + (\epsilon_{njr}p_r )p_m] + B_3 \epsilon_{mnj}\}\hat{k}_n(\delta_{im} - \hat{k}_i\hat{k}_m). \end{equation}

Now, using the convolution theorem, the integral for the angular velocity contribution due to the hydrodynamic component of the stratification-induced torque in (2.21) may be written as

(4.23)\begin{equation} Ri_v \displaystyle\int \rho'^{(1)} \hat{g}_j U^{(2)}_{ji} \,\textrm{d}V = Ri_v \displaystyle\int \,\textrm{d}{\boldsymbol k} \hat{\rho}'^{(1)}({\boldsymbol k}) \hat{g}_j U^{(2)}_{ji}(-{\boldsymbol k}), \end{equation}

where, in applying the convolution theorem, we have assumed the volume integral on the left-hand side of (4.23) to extend over the entire domain, and thereby, neglected the $O(L^3)$ volume of the spheroid. As the dominant contribution arises from length scales of $O(LRi_v^{-{1}/{3}})$, this neglect only amounts to an error of $O(Ri_v)$ in the torque integral, and $O(Ri_v^2)$ in the resulting angular velocity.

Using (4.16) and (4.22) in the Fourier space torque integral in (4.23), and after some simplification which includes defining a rescaled wave vector $2{\rm \pi} {\boldsymbol k}$, one obtains the angular velocity induced by the hydrodynamic stratification torque as

(4.24)\begin{gather} \varOmega_i^{(1)d} =Ri_v\frac{\mathrm{i}\tilde{F}}{8{\rm \pi}^3Y_C} \displaystyle\int \,\textrm{d}{\boldsymbol k} \frac{[1-(\hat{k}_x\hat{g}_x)^2]}{\{\mathrm{i}k^3(\hat{k}_y \hat{U}_y) - Ri_v[1-(\hat{k}_z \hat{g}_z)^2]\}k}\left[ B_1\left\{ (\epsilon_{irj}p_r\hat{g}_j)(\hat{k}_mp_m) \right. \right. \nonumber\\ \left. \left. + (\epsilon_{irj}p_r\hat{k}_j)(\hat{g}_mp_m)- 2 (\hat{k}_m\hat{g}_m)(\hat{k}_jp_j)(\epsilon_{irl}p_r\hat{k}_l) \right\} + B_3 \epsilon_{ijr}\hat{g}_j \hat{k}_r \right], \end{gather}

the terms proportional to $B_1$ and $B_3$ being the stresslet and rotlet-induced torque contributions, respectively. Redefining the new wave vector to be $Ri_v^{-{1}/{3}}{\boldsymbol k}$, so it remains of order unity on length scales of order the stratification screening length (and, thereby, pertains to the outer region in Fourier space), and considering only the real part of the integral above, one obtains

(4.25) \begin{gather} \varOmega_i^{(1)d} =Ri_v^{{2}/{3}} \frac{\tilde{F}}{8{\rm \pi}^3Y_C} \left[B_1 \displaystyle\int \,\textrm{d}{\boldsymbol k}\frac{[1-(\hat{k}_x\hat{g}_x)^2]k^2(\hat{k}_v \hat{U}_v)}{\{k^6(\hat{k}_y \hat{U}_y)^2 +[1-(\hat{k}_z \hat{g}_z)^2]^2\}}[ \{ (\epsilon_{irj}p_r\hat{g}_j)(\hat{k}_mp_m) \right.\nonumber\\+ (\epsilon_{irj}p_r\hat{k}_j)(\hat{g}_mp_m) - 2 (\hat{k}_m\hat{g}_m)(\hat{k}_jp_j)(\epsilon_{irl}p_r\hat{k}_l)\}]\nonumber\\ \left.+ B_3 \displaystyle\int \,\textrm{d}{\boldsymbol k}\frac{[1-(\hat{k}_x\hat{g}_x)^2]k^2(\hat{k}_v \hat{U}_v)}{\{k^6(\hat{k}_y \hat{U}_y)^2 +[1-(\hat{k}_z \hat{g}_z)^2]^2\}}\epsilon_{ijr}\hat{g}_j \hat{k}_r \right], \end{gather}

where the angular velocity due to the hydrodynamic stratification torque finally comes out to be $O(Ri_v^{{2}/{3}})$, as anticipated by the scaling arguments above, and the convergent Fourier integrals in (4.25) are evaluated below; note that the imaginary part of (4.24), neglected in (4.25), may be shown to equal zero by symmetry. Before evaluating the integrals above using a specific coordinate system, we note that the force–velocity relationship for a sedimenting spheroid, for the scalings used here, is given by

(4.26) \begin{equation} \hat{U}_i = \left[p_i p_j + \frac{X_A}{Y_A}(\delta_{ij} - p_i p_j)\right] \hat{g}_j. \end{equation}

Defining the aspect-ratio-dependent resistance ratio $An(\kappa ) = X_A/Y_A$, equation (4.26) may be written as

(4.27)\begin{equation} \hat{U}_i = [(1-An)p_i p_j + An \delta_{ij}] \hat{g}_j. \end{equation}

where $An$ decreases monotonically from unity for a sphere to a minimum of $1/2$ for an infinitely slender prolate spheroid ($\kappa \rightarrow \infty$); on the oblate side, $An$ increases from unity to a maximum of $3/2$ for a flat disk ($\kappa \rightarrow 0$). Both of these may be readily verified based on the expressions given in Appendix A.

To evaluate the above Fourier integrals, we choose a spherical coordinate system with its polar axis along $\hat {\boldsymbol U}$. Interestingly, a numerical evaluation in the alternate and perhaps more natural choice of a $\hat {\boldsymbol g}$-aligned coordinate system turns out to much more involved, with the individual integrals making up the torque diverging as $Pe^{{1}/{2}}$ in the limit $Pe \rightarrow \infty$, the divergence arising due to the buoyant jet mentioned above, and that corresponds to the region $\hat {\boldsymbol k} \boldsymbol {\cdot } \hat {\boldsymbol U} =0,\ k \ll 1$ in Fourier space (see Varanasi & Subramanian Reference Varanasi and Subramanian2021). We have verified that the divergences of the individual contributions cancel out, and the total torque integral is nevertheless convergent and independent of $Pe$ for $Pe \rightarrow \infty$, matching the result obtained below in the $\hat {\boldsymbol U}$-aligned system. In the latter coordinate system, $\hat {\boldsymbol U} = U_m {\boldsymbol 1}_U$, ${\boldsymbol 1}_U$ being the unit vector along $\hat {\boldsymbol U}$; using (4.27), one finds $U_m = \cos ^2 \psi + An^2 \sin ^2 \psi$, $\psi$ being the angle between $\hat {\boldsymbol g}$ and ${\boldsymbol p}$ defined in earlier sections. The unit wave vector is given by $\hat {\boldsymbol k} = \sin \theta \cos \phi {\boldsymbol 1}_{U^{\perp ^1}} + \sin \theta \sin \phi {\boldsymbol 1}_{U^{\perp ^2}} + \cos \theta {\boldsymbol 1}_U$, with ${\boldsymbol 1}_{U^{\perp ^1}}$ chosen to lie in the plane of sedimentation, that is, the plane containing the vectors ${\boldsymbol p}$, $\hat {\boldsymbol U}$ and $\hat {\boldsymbol g}$, so the torque acting on the spheroid points along $\hat {\boldsymbol {g}}\wedge \hat {\boldsymbol {p}}$ =${\boldsymbol 1}_{U^{\perp ^2}}$. With this choice, $\hat {\boldsymbol k} \boldsymbol {\cdot } \hat {\boldsymbol U} =U_m \cos \theta$. We take ${\boldsymbol p} = \cos \psi _U {\boldsymbol 1}_U + \sin \psi _U {\boldsymbol 1}_{U^{\perp ^1}}$, $\psi _U$ being the angle between ${\boldsymbol p}$ and $\hat {\boldsymbol U}$. Noting from (4.27) that $\hat {\boldsymbol U} \boldsymbol {\cdot } {\boldsymbol p} = \hat {\boldsymbol g} \boldsymbol {\cdot } {\boldsymbol p}$, one has $\cos \psi _U = \cos \psi /U_m$. With $\psi _U$ as defined above, $\hat {\boldsymbol k} \boldsymbol {\cdot } {\boldsymbol p} = \cos \psi _U \cos \theta + \sin \psi _U \sin \theta \cos \phi$, and $\hat {\boldsymbol k} \boldsymbol {\cdot } \hat {\boldsymbol g} = An^{-1}[\hat {\boldsymbol k} \boldsymbol {\cdot } \hat {\boldsymbol U} - (1-An) \cos \psi \hat {\boldsymbol k} \boldsymbol {\cdot } {\boldsymbol p}]$. As mentioned above, the torque only has a component along ${\boldsymbol 1}_{U^{\perp ^2}}$ in the chosen coordinate system. As the vectors involved in the integrals in (4.25) are $\hat {\boldsymbol g} \wedge \hat {\boldsymbol k}$, ${\boldsymbol p} \wedge \hat {\boldsymbol g}$ and ${\boldsymbol p} \wedge \hat {\boldsymbol k}$, we have

(4.28)\begin{gather} \epsilon_{2jr}\hat{g}_j\hat{k}_r =\frac{1}{An}[U_m \sin \theta \cos \phi - (1-An)\cos \psi(\cos \psi_U \sin \theta \cos \phi - \sin \psi_U \cos \theta)], \end{gather}

(4.29)\begin{gather}\epsilon_{2jr}p_j\hat{g}_r ={-}\frac{U_m\sin\psi_U}{An}, \end{gather}

(4.30)\begin{gather}\epsilon_{2jr}p_j\hat{k}_r ={-} \sin \psi_U \cos \theta + \cos \psi_U \sin \theta \cos \phi, \end{gather}

which defines all quantities involved in the calculation. Using these expressions in the integrals in (4.25), the $k$-integration is carried out analytically while the remaining two angular integrals are evaluated numerically using Gaussian quadrature. From (4.28)–(4.30), and other quantities defined in the preceding text, the $\psi$-dependence of the large-$Pe$ hydrodynamic angular velocity is seen to be more complicated than the $\cos \psi \sin \psi$ dependence obtained earlier for the inertial and hydrostatic contributions in § 3, and for the hydrodynamic component, in the limit $Pe \ll Ri_v^{{3}/{5}}$, in § 4.1. This is along expected lines because the large-$Pe$ limit examined in this section is a singular perturbation problem, as evident from the outer region contributing at leading order (this aspect is also seen in the outer-region contribution to the torque at low $Pe$, derived in Appendix B). Another feature of the singular character is that, unlike the earlier torque contributions, the large-$Pe$ hydrodynamic stratification torque is in general a non-separable function of $\psi$ and $\kappa$.

Figures 5(a) and 6(a) show plots of the angular velocity ($\varOmega ^{(1)d}$), due to the hydrodynamic component of the stratification torque, versus $\psi$ for prolate and oblate spheroids, respectively. As evident from figure 5(a), for prolate spheroids, the magnitude of the angular velocity is expectedly small in the near-sphere limit, increasing monotonically with $\kappa$ to a (finite) maximum in the limit of a slender fiber. Figures 5(b) shows plots of the angular velocity scaled with the square of the eccentricity ($\xi _0^{-2}$), so as to obtain a collapse in the near-sphere limit. Note that the finite value of the stratification-induced angular velocity in the limit of a slender fiber ($\xi _0 \rightarrow 1$) is in contrast to the $O[\ln (\xi _0-1)]^{-1}$ scaling exhibited by the inertial angular velocity calculated in § 3 (see also Dabade et al. Reference Dabade, Marath and Subramanian2015), and implies that, for fixed $Re$ and $Ri_v$, the stratification torque invariably becomes dominant for large aspect ratios. Figure 6(b) confirms the squared-eccentricity scaling for oblate spheroids with near-unity aspect ratios; expectedly, the angular velocity approaches a finite value in the limit of a flat disk. Importantly, for both prolate and oblate spheroids, the sign of $\varOmega ^{(1)d}$ is such as to rotate the spheroid onto an edgewise orientation.

Figure 5.

The scaled angular velocity ($Ri_v^{-2/3}\varOmega ^{(1)d}$), due to the hydrodynamic component of the stratification torque, for prolate spheroids of different aspect ratios: (a) $Ri_v^{-2/3}\varOmega ^{(1)d}$, for prolate spheroids, as a function of the spheroid inclination with gravity; (b) $Ri_v^{-2/3}\varOmega ^{(1)d}$, normalized by the near-sphere scaling ($1/\xi _0^2$) for prolate spheroids, plotted as a function of the spheroid inclination with gravity.

Figure 6.

The scaled angular velocity ($Ri_v^{-2/3}\varOmega ^{(1)d}$), due to the hydrodynamic component of the stratification torque, for oblate spheroids of different aspect ratios: (a) $Ri_v^{-2/3}\varOmega ^{(1)d}$, for oblate spheroids, as a function of the spheroid inclination with gravity; (b) $Ri_v^{-2/3}\varOmega ^{(1)d}$, normalized by the near-sphere scaling ($1/\xi _0^2$) for oblate spheroids, plotted as a function of the spheroid inclination with gravity.

The non-trivial orientation dependence of the angular velocity referred to in the previous paragraph is also evident from the plots in figures 5(a) and 6(a). For near-unity aspect ratios, the angular velocity curve is nearly symmetric about $\psi = {{\rm \pi} }/{4}$; that the angular dependence in this limit is indeed of the form $\sin \psi \cos \psi$ may be shown based on the fact that for $An \rightarrow 1$, $\psi _U \approx \psi$. The asymmetry about $\psi = {{\rm \pi} }/{4}$ increases as the aspect ratio departs from unity, with the location of the maximum angular velocity moving to $\psi$ greater than and less than ${{\rm \pi} }/{4}$ for prolate and oblate spheroids, respectively, as shown in figure 7. To see the deviation of the angular dependence from the aforementioned simple form more clearly, in figures 8(a) and 8(b) we plot the angular velocity, scaled by the inertial angular velocity which is proportional to $\sin \psi \cos \psi$, again as a function of $\psi$. For near-unity aspect ratios, one obtains a horizontal line, while for both larger and smaller aspect ratios, this renormalized angular velocity asymptotes from one plateau for $\psi \rightarrow 0$ to a second one for $\psi \rightarrow {{\rm \pi} }/{2}$.

Figure 7.

The angle corresponding to the maximum angular velocity, arising from the hydrodynamic component of the stratification torque, plotted as a function of the spheroid aspect ratio (both prolate and oblate spheroids); the inset shows the variation of this angle on a log–log scale, emphasizing the approach to finite values for extreme aspect ratios ($\kappa = 0$ and $\infty$).

Figure 8.

Ratio of the angular velocities due to the hydrodynamic stratification and inertial torques: (a) prolate spheroid; (b) oblate spheroid. The stratification-induced rotation is higher for the near-vertical and near-horizontal orientations for prolate and oblate spheroids, respectively; note that the angular velocity ratio at $\psi = {\rm \pi}/4$ has been subtracted for convenient depiction.

4.2.1. A closer look at the $Pe\gg 1$ analysis

In contrast to the inertial contribution determined in § 3, which was independent of the ambient stratification at leading order, the stratification-induced torque can, in principle, be coupled to inertial forces even in the limit $Re, Ri_v \ll 1$. For sufficiently small $Pe$ ($Pe \ll Ri_v^{{3}/{5}}$ as argued in § 4.1.1), the density perturbation that determines the stratification torque arises from a diffusive response to the no-flux condition on the surface of the sedimenting spheroid, and is therefore independent of the fluid motion. As a result, a non-trivial coupling between stratification and inertia occurs primarily for large $Pe$. As the dominant length scales contributing to the stratification torque in this limit are much larger than $O(L)$, the magnitude of the density perturbation is controlled by the convection of the ambient stratification by the far-field disturbance fluid motion. The nature of this convection is, therefore, dependent on the form of the disturbance velocity field, and this in turn depends on the relative magnitudes of the inertial ($LRe^{-1}$) and stratification ($LRi_v^{-{1}/{3}}$) screening lengths. The calculation of the angular velocity due to the hydrodynamic stratification torque detailed above pertains to the Stokes stratification regime, with $Re \ll Ri_v^{{1}/{3}}$, where the disturbance velocity field directly transitions from the Stokeslet form to a more rapid decay on length scales of $O(LRi_v^{-{1}/{3}})$ (Varanasi & Subramanian Reference Varanasi and Subramanian2021) and is therefore independent of $Re$. In the Stokes stratification regime, therefore, the stratification-induced rotation, both in the limit of small and large $Pe$, is independent of $Re$, and the inertial and stratification angular velocity contributions are additive. This will no longer true when $Re \geq O(Ri_v^{{1}/{3}})$, corresponding to the so-called inertia-stratification regime, in which case the leading-order stratification-induced rotation for large $Pe$ will be a function of $Re/Ri_v^{{1}/{3}}$. In the limit $Ri_v^{{1}/{3}} \ll Re \ll 1$, opposite to that analysed above, the disturbance velocity transitions from an $O(1/r)$ to an $O(1/r^2)$ decay (outside of a viscous wake) across length scales of order the inertial screening length. This leads to the stratification torque integrand decaying as $O(1/r^3)$ for length scales much larger than $O(LRe^{-1})$, and the torque integral in (2.22) continues to exhibit a logarithmic divergence. This (milder) divergence is only eliminated when buoyancy forces become comparable to inertial forces at a secondary screening length that was estimated in Varanasi & Subramanian (Reference Varanasi and Subramanian2021) to be $O(Re/Ri_v)^{{1}/{2}}$. Accounting for the aforementioned cut-off of the logarithmic divergence, the angular velocity arising from the hydrodynamic stratification torque is expected to have a leading $O[Ri_v Re^{-1}\ln (Re/Ri_v^{{1}/{3}})]$ contribution arising from a region between the primary and secondary screening lengths (that is, due to the logarithmic growth for $Re^{-1} \ll r \ll (Re/Ri_v)^{{1}/{2}}$), with logarithmically smaller $O(Ri_v Re^{-1})$ contributions arising from length scales of order the two screening lengths. Assuming this angular velocity to rotate the spheroid towards an edgewise configuration, and equating it to the $O(Re)$ inertial contribution, one obtains $Re \leq Ri_v^2$ for a transition to an edgewise-settling regime. This, however, contradicts the requirement $Re \gg Ri_v^{{1}/{3}}$ characterizing the inertia-stratification regime, implying that the inertial angular velocity contribution remains dominant in this regime. Thus, one concludes that, in the limit of small $Re$ and $Ri_v$, a broadside-on–edgewise transition is possible only in the Stokes stratification regime.

To end this subsection, we again examine the validity of a quasi-steady state assumed in the analysis above for large $Pe$. As already argued in § 3, momentum diffusion occurs asymptotically fast for small $Re$ and, therefore, the quasi-steady assumption used to evaluate the stratification torque integral relies on the time scale for the density disturbance to approach a steady state being much shorter than that characterizing spheroid rotation. The former time scale may be regarded as that required to convect the density perturbation through the $O(LRi_v^{-{1}/{3}})$ stratification screening length, and is therefore $O(L/URi_v^{-{1}/{3}})$. The time scale of rotation is $O(L/URe^{-1})$ or $O(L/URi_v^{-{2}/{3}})$, depending on which of $Re$ or $Ri_v^{{2}/{3}}$ is greater. In either case the time scale for the development of a steady density perturbation is smaller, provided one remains in the Stokes stratification regime $Re \ll Ri_v^{{1}/{3}} \ll 1$. Note that the perturbation density field, for $Pe = \infty$, is expected to be logarithmically singular along the rear stagnation streamline, as has been shown for the case of a spherical particle (see Varanasi & Subramanian Reference Varanasi and Subramanian2021), with this singularity either being resolved on a larger diffusive time scale (that is, due to $Pe$ being regarded as large but finite), or on account of the unsteadiness arising from the rotation of the settling spheroid. However, the convergence of the Fourier integrals involved in the stratification torque above implies that the contribution of the transiently developing region in the immediate neighbourhood of the rear stagnation streamline is irrelevant as far as the leading-order hydrodynamic stratification torque is concerned, and a quasi-steady analysis of this torque remains valid for $Re \ll Ri_v^{{1}/{3}}$.