Spherical particle sedimenting in weakly viscoelastic shear flow

We consider the dynamics of a small spherical particle driven through an unbounded viscoelastic shear flow by an external force. We give analytical solutions to both the mobility problem (velocity of forced particle) and the resistance problem (force on fixed particle), valid to second order in the dimensionless Deborah and Weissenberg numbers, which represent the elastic relaxation time of the fluid relative to the rate of translation and the imposed shear rate. We find a shear-induced lift at $O({\rm Wi})$, a modified drag at $O({\rm De}^2)$ and $O({\rm Wi}^2)$, and a second lift that is orthogonal to the first, at $O({\rm Wi}^2)$. The relative importance of these effects depends strongly on the orientation of the forcing relative to the shear. We discuss how these forces affect the terminal settling velocity in an inclined shear flow. We also describe a new basis set of symmetric Cartesian tensors, and demonstrate how they enable general tensorial perturbation calculations such as the present theory. In particular this scheme allows us to write down a solution to the inhomogenous Stokes equations, required by the perturbation expansion, by a sequence of algebraic manipulations well suited to computer implementation.


I. INTRODUCTION
In this paper we consider the mobility of a small spherical particle driven through an unbounded viscoelastic shear flow by an external force. In a Newtonian fluid the velocity of the particle is determined by the balance between the Stokes drag and the external force, and it is unaffected by the shear flow because of the linearity of Stokes equations. But in a viscoelastic fluid the disturbance flow around the particle interacts non-linearly with the shear flow to induce viscoelastic stresses. As a consequence the mobility depends non-linearly on the forcing and the shear flow.
A viscoelastic shear flow can reduce the terminal velocity of a sphere when the applied shear flow is perpendicular to gravity [1][2][3][4]. This so-called cross-shear flow is a model system for transport of particles in vertical cracks induced by hydraulic fracturing [5]. Experiments by van den Brule and Gheissary [1] first demonstrated that a cross-shear flow strongly reduces the settling velocity, and that fluid elasticity is the dominant mechanism. Recently, numerical simulations by Padhy et al. [3,6] verified an increased drag on a sphere translating through a cross-shear flow, and showed that the experimental observation is explained by a combination of viscoelasticity and the effects of the nearby walls in the experiment. Calculations by Housiadas and Tanner [2,7] demonstrate that the drag is increased also in an unbounded viscoelastic cross-shear flow.
These studies concern settling in a cross-shear flow in which the gravity acts along the vorticity axis. In this case the physical system is invariant under a 180 • rotation around the vorticity axis. This symmetry was exploited in both the analytical and numerical calculations to reduce the number of variables [2,3,7]. In particular the only relevant force is the drag force, and the particle only rotates around the vorticity axis. The effect of the shear on settling is substantial in this symmetrical case, and this fact raises new questions: How are the dynamics affected when gravity acts at an angle to the vorticity? Are there additional forces and torques when the symmetry is broken? How does the drag change as the angle between gravity and flow vorticity changes?
In this paper we calculate the effect of an unbounded shear flow on the terminal particle velocity for any orientation of the external force relative to the shear. We must consequently abandon the simplifications of the symmetrical case. We must allow for lift forces in the flow-shear plane, and for rotation around any axis. Further, since viscoelasticity is a non-linear effect, it is not possible to construct the general result as a linear combination of results for two independent directions. Therefore we must solve the general perturbation problem for the flow velocity u and viscoelastic stress tensor Π around a translating and rotating particle in a shear flow. Our calculation relies on a perturbation theory for weak elasticity, valid to second order in the Deborah and Weissenberg numbers. These dimensionless numbers relate the elastic relaxation time of the fluid relative to the rate of translation and the shear rate.
Brunn [8,9] and Vishnampet and Saintillan [10] considered the first order of this problem, and both found lateral migration, although their detailed results do not agree with each other. The first part of our calculation is an independent check of their results, which we return to in Section IV A. Housiadas and Tanner [11] and D'Avino et al. [12] considered the angular velocity of the sphere in the absence of an external force, and Leslie and Tanner [13] calculated the drag on a sphere in absence of shear flow. These two results coincide with our theory in their respective limits.
We solve the problem in tensorial form. Our solution does not refer to any coordinate representation such as spherical coordinates. Since the governing equations and boundary conditions are tensorial in nature, this substantially simplifies the calculations. All steps of the calculation are algebraic, and therefore well suited to computer implementation. To achieve this we introduce a new basis set of symmetric, rank-n Cartesian tensors. We describe these tensors and how to calculate with them in some detail in Section III, because we expect that they will be useful for treating other problems too.
We present two related calculations. The first is the mobility problem, where we impose an external force F ext on the particle, and compute the resulting particle velocity v(F ext ). This corresponds directly to the experimental protocol of for example van den Brule and Gheissary [1], where they release a sphere in a cylindrical Couette device and measure the steady settling velocity. The other question is the resistance problem, where we prescribe the particle velocity v, and compute the resulting force F (v) exerted by the fluid on the particle. This approach corresponds to the calculations by Housiadas and Tanner [7] and numerical simulations by Padhy et al. [6]. The two are related, because given the solution to the mobility problem v(F ext ), and the solution of the resistance problem F (v), it must hold that F (v) = −F ext . In this paper we solve both the mobility problem and the resistance problem for a freely rotating spherical particle in an unbounded viscoelastic shear flow, with no restriction on the direction of v or F ext relative to the shear.
The rest of this paper is organized as follows. In Section II we describe the problem, give the governing equations, and describe how we apply the Lorentz reciprocal theorem. In Section III we explain our algebraic solution of the inhomogeneous Stokes equation in terms of Cartesian tensors, and summarise their algebraic properties. We summarise our calculation and give the final result in Section IV. We discuss the results and conclude in Section V.

A. Equation of motion and dimensionless parameters
We consider the steady-state motion of a spherical particle of radius a, suspended in a viscoelastic fluid and subject to an external force F ext . For concreteness we may think of the gravitational force F ext = 4πa 3 (ρ p − ρ f )g/3. The particle moves with center-of-mass velocity v, and rotates with angular velocity ω. Far away from the particle the flow is a simple shear flow where Ω is half the flow vorticity, and the symmetric tensor S is the rate of strain. In a simple shear flow the vorticity and strain are related by SΩ = 0 and 2|Ω| 2 = Tr SS, in contrast to a general linear flow. We work in dimensionless variables. The length scale is given by the particle radius a. The time scale is given by the reciprocal of the imposed shear rate s = √ 2Tr SS, which also determines the scale of u ∞ to sa. The particle and disturbance flow velocities are nondimensionalized by the characteristic flow velocity v c past the particle. In the resistance problem, v c is simply the magnitude |v| of the imposed velocity v. In the mobility problem we estimate the characteristic speed by v c = F ext /aµ, related to the terminal velocity in Stokes flow under an external force of magnitude F ext . Here µ is the total viscosity, defined precisely in conjunction with the constitutive equations below. Stresses are made dimensionless by v c µ/a, and forces by v c µa. In the remainder of this paper all quantities are dimensionless: t = st, r = r/a, S = S/s, Ω = Ω/s, F = F /(v c µa), and so forth. We drop the primes since all quantities are dimensionless.
It follows that there are two dimensionless parameters that govern this problem, corresponding to the translational and rotational motion of the particle compared to the relaxation time λ of the viscoelastic fluid. The Deborah number De = λv c /a is associated with the time scale of convective flow over the particle size. The Weissenberg number Wi = λs is associated with the shear rate. The ratio α = Wi/De measures the relative importance of the imposed shear to the translational motion.
The perturbation theory in this paper is valid in the limit De 1 and Wi 1. This implies that the elastic part of the fluid relaxes quickly relative to the rate at which it is deformed by the moving particle and the shear flow. In the following we expand in De, and treat the ratio α as an O(1)-quantity. But in the end we give the result in terms of De and Wi. We note that the choice between 1/s and a/v c for the characteristic timescale is arbitrary up to factors of α = Wi/De, which we assume to be O(1).
We neglect the effects of fluid inertia. This requires that the viscous relaxation time of the fluid is shorter than the elastic relaxation time λ, so that we can neglect effects of inertio-elastic coupling. More precisely the particle Reynolds number Re p = ρ f v c a/µ De. This condition is equivalent to ρ f a 2 /µ λ. We write down the dimensionless governing equations for with respect to a frame moving with the steady centerof-mass velocity v. In this frame the fluid pressure p and velocity u satisfy The viscoelastic stress tensor Π is modeled by the steady Oldroyd-B constituitive equations [14]. They describe a suspension of elastic dumbbells, which is one of the simplest models of an elastic polymer suspension that exhibits a normal stress difference in a shear flow. The equations are The parameter µ r = µ p /(µ s + µ p ) is the relative contribution to the total viscosity from the elastic polymers, relative to the solvent viscosity µ s . We denote the total viscosity µ = µ s + µ p . The flow problem in Eqns. (2) and (3) is completed by the no-slip boundary condition on the particle surface S p and that it approaches u ∞ as |r| → ∞: The force and torque on the particle are given by We do not consider any external torques in this paper and so T ext = 0. We compute the forces on the particle, or the resulting velocities, with the Lorentz reciprocal theorem [15]. We outline how to apply it to our problem in the following Section.

B. Reciprocal theorem
In the context of perturbation theory, the Lorentz reciprocal theorem relates certain integral quantities such as the force or torque to order n + 1, given the detailed flow solution to only order n [15,16]. For non-Newtonian flows in particular, the method has for example been used to calculate lateral migration [17], and orbit drift of non-spherical particles [18]. In this Section we state the theorem as applicable to our problem and notation.
We denote the Newtonian part of the stress by σ N = −pI + (∇u + (∇u) T ), and the "extra stress" σ E = µ r (Π − (∇u + (∇u) T )), so that σ = σ N + σ E . The flow equation of motion (2a) is therefore The Lorentz reciprocal theorem for an arbitrary Stokes flow (ũ,σ) and the flow defined in Eq.
(2) reads [15] Here V is any volume outside the particle, and S denotes the surfaces bounding V . The vector dS = ndS, where n is the surface normal pointing out of V . Using σ = σ N + σ E , it follows from (8) that In the first two surface integrals in (9) we identify the hydrodynamic force and torque on the particle, as given in Eq. (6). We takeũ to be the Stokes flow around a spherical particle translating with velocityṽ and rotating with angular velocityω in an otherwise quiescent fluid. We write this auxiliary flow asũ = M vṽ + M ωω , and we also know thatF = −6πṽ andT = −8πω. Finally, upon taking the size of the volume V to infinity, and inserting the boundary conditions (4) into Eq.
The surface integrals 'at infinity' do not contribute, and the remaining surface integrals are only over the particle surface. This is well known for the disturbance quantities, say σ E (u) − σ E (u ∞ ), because the integrands decay faster than 1/r 2 [16]. The potentially problematic terms are those from σ E (u ∞ ) that are independent of r. But M v is an even function of r, so that surface integral vanishes by symmetry. On the other hand, M ω integrated over the sphere is an antisymmetric tensor that vanishes upon contraction with the symmetric stress tensor. Becauseṽ andω may be chosen arbitrarily, we have two separate theorems for the force and torque: In this paper we do not consider external torques, and therefore T = 0 implies The reciprocal theorem may be used to solve either the resistance problem, or the mobility problem. For the resistance problem we take F ext = 0 and use Eq. (11). For the mobility problem we require the total force F = 0 and find v = 1 6π The integrands in Eqs. (11)(12)(13)(14) are functions of the yet unknown σ E . In Section IV we evaluate these integrals by calculating σ E perturbatively. As a final note, Eq. (11) is equivalent to the integral theorem Ho and Leal [17] used to compute the lateral drift of a spherical particle in wall-bounded flow. Their Eq. (2.22) follows from our Eq. (11) because However, we evaluate the two integral contributions in Eq. (11) separately, because they give the contributions from two different physical mechanisms. The surface integral represents the extra polymer stress acting directly on the particle surface (see Eq. (61) and (67) in Section IV). The volume integral represents the indirect effect that the polymer stress modifies the flow field, which in turn modifies the viscous stress on the particle.

III. METHOD: T-TENSORS
In this Section we introduce a basis set of symmetric rank-n Cartesian tensors T nl i1i2..in . Each of these basis tensors is a linear combination of spherical harmonics Y m l with a particular value of the angular momentum quantum number l, but different modes m. Therefore the Cartesian tensors share many useful properties with the spherical harmonics. For example, surface integrals vanish unless l = 0, tensors with different values of l are orthogonal with respect to integration over the unit sphere, and they have known Fourier transforms.
Because of their direct relation to the spherical harmonics, the T -tensors are an alternative basis for Lamb's general solution for Stokes flow [15]. But as explained in detail below, a rank-n T -tensor is also closely related to the rank-n polyadr i1ri2 ..r in of a unit vectorr. Together with the radial functions 1/r m these polyads are the building blocks of the familiar multipole expansion for Stokes flow, for example the Stokeslet δ ij /r +r irj /r, or the rotlet ε ijkrj /r 2 . Therefore the T -tensors stand as a new alternative between Lamb's general solution in spherical coordinates, and the Cartesian multipole expansion. Although any calculation may in principle be performed in any of these representations, we found that the basis described here is suitable for implementation in computer algebra. In particular it enables us to write down particular solutions to inhomogenous Stokes equations in tensorial form, without any explicit coordinate representation, and without explicitly solving differential equations.
In this Section we use index notation to avoid any ambiguity. When appropriate we use the vector notation T nl , remembering that T nl is rank n and symmetric in all indices.

A. Definition
We consider the rank-n polyadr i1ri2 ..r in of a unit vectorr. Any given polyad is a smooth function defined on the sphere, and may be expanded in the spherical harmonicŝ For our present purposes it is not necessary to calculate the expansion coefficients c nlm i1i2..in , but we may deduce the following important fact. The left hand side is a polynomial of order n, and every term in the sum on the right hand side is a polynomial of order l, so we conclude that c nlm i1i2..in = 0 if l > n. Thereforê We define the tensor T nl i1i2..in as the inner sum in Eq. (17), so that and therefore by constructionr We complete the definition by T 00 = 1.

B. Properties of the T -tensors
Symmetry. From the expansion (19) it is clear that any T -tensor is symmetric in all indices. By definition only tensors with l ≤ n are non-zero: Further, the polynomial in the left hand side of Eq. (19) has parity (−1) n under inversion ofr, and every term on the right hand side has parity (−1) l . Therefore T nl is non-zero only if both n and l are even, or if both n and l are odd.
Integrals. Any two T -tensors are orthogonal with respect to integrals over the unit sphere S, because the spherical harmonics enjoy this property: It follows that Cartesian rank. Taking a trace, i.e. contracting any two indices, of a T -tensor lowers its rank n by two: This follows from Eq. (19) and the orthogonality property (21). An important consequence is that T ll i1i2..i l is traceless: Conversely, for l ≤ n − 2, we raise the Cartesian rank n by The parenthesis in Eq. (25) contains n(n − 1)/2 terms, one for each unique pairing of the n indices. We have not proven Eq. (25) for general values of n and l, but it is straightforward to work out the cases l = n − 2, l = n − 4, and so on, by taking a trace of Eq. (25) and using Eqns. (23) and (24) repeatedly. We have checked all values of n and l that are used in our calculations in this paper. Dimensionality. The tensor T ll is a rank-l Cartesian tensor. In general it could have 3 l unique elements (in three spatial dimensions). In contrast, there are only 2l + 1 spherical harmonics Y m l of degree l, corresponding to the values m = −l... l. But we have shown that T ll is symmetric and traceless. A symmetric tensor of rank l has (l + 1)(l + 2)/2 unique elements, and it has l(l − 1)/2 unique traces that we require to vanish. These conditions leave exactly 2l + 1 degrees of freedom for a symmetric and traceless rank-l Cartesian tensor. This is the reason we refer to the T -tensors as a basis set. The coefficients c llm i1i2..i l in Eq. (18) are the elements of the 'rotation matrix' between the two basis sets T ll and Y m l . We claim that this transformation is unitary for a certain choice of normalization of the spherical harmonics. In other words, it is in fact a proper rotation. However, we have not proven this for general values of l, but we confirmed that it is true up to l = 8 by brute force calculation of c llm from the definition Eq. (18), see Appendix A.
Multiplication. The product of two T -tensors follows directly from (19) as a recurrence relation: Here i and j are short-hand for i 1 ..i n and j 1 ..j n . We use Eq. (26) in practical calculations, but there is also a largely unexplored connection to quantum angular momentum algebra and Clebsch-Gordan coefficients. In particular, it can be shown (Appendix A) that for some coupling tensor A independent ofr. Relation to polyads. We convert any polyadic expression into T -tensors by replacingr i → T 11 i , and applying Eq. (26) until no products remain. Conversely any T -tensor is expressed as a polyadic by recursively using and Eq. (25). The first few tensors are Differentiation. In order to calculate using only algebraic manipulations on the T -tensors, we must know how the gradient operator ∇, where ∇ i = ∂/∂r i , acts on them. We will show that the action of ∇ is to 'scatter' a tensor of degree l into a linear combination of tensors with degrees l − 1 and l + 1 (given in Eq. (34)). We will briefly describe how to compute the coefficients of this linear combination, for any l.
Consider the differential operator T nl i (∇), defined by taking the polynomial r l T nl i (r) and replacing the components of r with the partial derivatives ∂/∂r i . Hobson's theorem on differentiation [19,20] explains that any such differential operator built from a harmonic polynomial acts on radial functions in a particularly simple way. In our case we use his result to find Therefore the general J-th order derivative is The product T N J i (∇) T nl i (∇) is given by Eq. (27), and the general formula follows from Eq. (30). In this paper we only consider first order derivatives which correspond to J = 1, because ∂/∂r i = T 11 i (∇). For J = 1 the general formula (32) and Eq. (27) give which becomes, using Eq. (30) and Eq. (31) , To evaluate Eq. (34) in our computer program we compute the product r m−1 T 11 T nl ≡ αr m−1 T l−1,l−1 +βr m−1 T l+1,l+1 using (26), and replace the coefficients by α → (m + l + 1)α, and β → (m − l)β. Fourier transform. We use a symmetric convention for the Fourier transform: The Fourier transform of r m T nl i (r) follows directly from that of the functions r m Y µ l (θ, ϕ) given in Ref. [21]. For all values of m and l that appear in the present calculation When m = l + 2j, r m T nl (r) is a polynomial in the components of r, and its Fourier transform is the Dirac delta function and its derivatives. The case m = −(l + 3) − 2j is more complicated, involving logarithms [21]. Neither of these cases arise in this paper.

C. Particular solution for the inhomogenous Stokes equation
Consider the inhomogenous Stokes problem where we assume that f i is a linear combination of T -tensors. The Fourier transform of Eq. (39) is where k = |k|. This algebraic equation is solved by wherek i ≡ k i /k is a unit vector. In terms of T -tensors, the Fourier space Green's function is The procedure to find the solution u i is therefore In this paper we never need an explicit expression for the pressure p, but if needed it is computed in the analogous way from Eq. (41).

IV. RESULTS
In this Section we give the solutions to both the mobility problem (Section IV A) and the resistance problem (Section IV B) for a freely rotating spherical particle in an unbounded viscoelastic shear flow, with no restriction on the direction of v or F ext relative to the shear.

A. The mobility problem
Here we consider a particle moving under the effect of an external force. The particle velocity v is a function of De and Wi to be determined, and to that end we require that the total force F = 0. We proceed with a regular perturbation expansion in De: At each order Π (k) is given by an algebraic equation, and u (k) by an inhomogenous Stokes equation (except for the lowest order, which is homogenous). To lowest order De 0 we have from Eq. (3) and therefore σ E(0) = 0, and therefore from Eqs. (14) and (13) we have To order De 0 the flow satisfies At order De 1 Eq. (3) gives where u (0) is known, but u (1) is still unknown. Consequently The reciprocal theorem (14) gives At this order the shear flow and the disturbance from the external forcing interact to create a lateral drift perpendicular to both the direction of forcing and the vorticity. The drift arises from both the extra stress on the particle surface, and the viscous stress induced by the viscoelastic medium, in proportion 2 : 3. The lateral drift Eq. (50) was first calculated by Brunn [8]. Our Eq. (50) agrees with his when accounting for the erratum [9] and letting his κ 11 0 = −2κ 22 0 , which corresponds to the second-order fluid limit of the Oldroyd-B model. The mobility derived in Ref. 10 is different from Eq. (50). In particular they report a contribution proportional to SF ext . The O(Wi) contribution to the angular velocity vanishes, in agreement with all the previous results [8][9][10].
With σ E(1) , v (1) , and ω (1) given by Eqs. (49-51) we can write down the inhomogenous Stokes problem for u (1) : subject to First we compute a particular solution u (1)p (r) as explained in Section III C. The flow field u (1)p satisfies the inhomogenous Eq. (52), but not the boundary conditions Eq. (53). We next solve for a Stokes flow u (1)h that satisfies the homogenous equation and the boundary conditions By construction u (1) = u (1)h + u (1)p . At order De 2 we have from Eq. (3) with where u (0) , u (1) , Π (0) and Π (1) are all known. The reciprocal theorem at order De 2 gives We first consider the angular velocity ω (2) . The angular velocity due to the surface integral in Eq. (59) vanishes, so the induced viscous stress alone explains the rotation rate at this order. The full result for the angular velocity, to second order in De and Wi, takes the form The numerically largest contribution is the O(Wi 2 ) slowdown of the rotation around vorticity. This contribution agrees with an earlier analytical result [11] that also explains numerical simulations [12]. The O(De 2 ) contribution shows a coupling between the external force and rotation rate, through the strain. It is numerically small, but may be important because it describes a rotation around another axis than Ω. For the particle velocity v (2) the surface integral of the extra stress on the particle evaluates to We see that this contribution may affect the velocity along F ext , but in particular it gives another lateral drift in a direction perpendicular to the O(Wi) lateral drift calculated above. Next we evaluate the volume integral in (58), and this gives the final result for v to second order in De and Wi: The induced viscous stress, given by the volume integral, contributes to the same terms as the extra stress on the surface shown in Eq. (61). In addition there is yet another velocity proportional to SSF ext , and a component in the direction of F ext . The second order contribution to the velocity directly proportional to F ext consists of one term proportional to De 2 , and one proportional to Wi 2 . The velocity along F ext increases as De increases, but the numerical prefactor is small. The important result is that the velocity decreases with increasing shear rate, as observed in experiment [1]. We discuss this effect for an inclined shear flow in Section V. In the next Section we solve the resistance problem that can be directly compared with earlier calculations for the cross-shear flow.

B. The resistance problem
The calculation for this problem is very similar to that of the mobility problem, so we omit most details. Here we consider a freely rotating sphere moving at velocity v through a shear flow, and calculate the resulting hydrodynamic force F and angular velocity ω. The crucial differences to the mobility problem are that F ext = 0, and v is a constant, independent of De.
The zeroth order problem is the corresponding Stokes problem, which determines u (0) and σ E(1) . The reciprocal theorem (13) gives The first order flow problem has a different boundary condition to that of the mobility problem, because v is independent of De. The first order equations are subject to The reciprocal theorem (11) gives The contributions from the surface and volume integrals are similar to those of the mobility problem. Specifically, After evaluating the volume integrals in Eq. (66), we find for the resistance problem with velocity v As expected, the expression for the resistance force (68) is similar to the expression for the mobility velocity (62) with F ext replaced by −6πv. However, they differ in a term Wi 2 µ 2 r Ω × Ω × v. This difference arises because of the lateral force at O(Wi) for the following reason. In the mobility problem the particle is allowed to relax this lateral hydrodynamic force by drifting sideways. But in the resistance problem we essentially force the fluid, through the boundary conditions, with the lateral force required to keep the particle moving with the prescribed velocity v. This forcing, or lack thereof, at O(Wi) is what gives the differing term at O(Wi 2 ). Upon substitution of the mobility velocity Eq. (62) into the expression for the resistance force Eq. (68) we find F = −F ext to second order in De 2 and Wi 2 , as advertised in the Introduction.
The drag term proportional to De 2 is the drag in absence of shear. This term was first calculated by Leslie and Tanner [13] in the context of sedimentation in a quiescent fluid. Our coefficient matches theirs when µ r = 1, and their = β = 0. For the cross-shear flow, v is parallel to Ω, so that Ω × v = Sv = 0. For this case only the first and third terms on the right-hand side of Eq. (68) remain. This expression agrees with the analytical result of Housiadas and Tanner [7], and therefore with Padhy et al. [3] as shown in their Fig. 7. In this example the external force lies in the plane spanned by the vorticity axis and the flow direction. This situation corresponds to settling between two far-apart shearing walls, parallel to the walls, but where the shearing is at an angle to gravity. We denote by ' the angle between F ext and vorticity.
FIG. 1. The inclined shear flow geometry discussed in Section V. In this example the external force lies in the plane spanned by the vorticity axis and the flow direction. This situation corresponds to settling between two far-apart shearing walls, parallel to the walls, but where the shearing is at an angle to gravity. We denote by ϕ the angle between F ext and vorticity.

V. SUMMARY AND DISCUSSION
We have derived analytical results for the linear and angular velocities of a particle driven through a viscoelastic shear flow by an external force, valid to second order in De and Wi, given in Eqs. (62) and (60).
We found three qualitatively different corrections to the predicted velocity in a Newtonian fluid. First, at O(Wi) there is a drift proportional to Ω × F ext , that is, perpendicular to the forcing and vorticity. Second, the resulting velocity along the forcing is modified at O(De 2 ) and O(Wi 2 ). The numerical prefactor of the De 2 -contribution is small, so in practice only the O(Wi 2 ) effect is important. These terms correspond to the effect of the imposed shear flow. Third, at O(Wi 2 ) there is yet another lateral drift, perpendicular to the first one. Even for Wi = 0.5 this second drift may be as strong as the O(Wi)-drift, but it points in another direction. The relative importance of these three effects depends strongly on the direction of external forcing relative to the orientation of the shear flow.
There are two mechanisms that contribute to these corrections. First, the extra stress acts directly on the particle, giving a force and a torque. Secondly, the extra stress acts on the fluid, which modifies the flow and indirectly gives a force and a torque via the viscous and pressure terms. The lateral drift at O(Wi) is a combination of these two mechanisms [Eq. (51)]. The decreased velocity of a sphere sedimenting in a cross-shear flow, however, is due to the indirect increase of viscous stress, and the direct contribution from the extra stress vanishes [Eq. (61)]. This observation is in qualitative agreement with the observations of numerical simulations [3]. When the forcing is at an angle to the vorticity vector, the correction is typically a combination of the two mechanisms.
The angular velocity around the vorticity slows down at O(Wi 2 ), in agreement with earlier results [11,12]. But at O(De 2 )we also find a coupling between the strain and translation that induces a rotation around the axis F ext ×SF ext . The prefactor is quite small, but the effect could be important because it describes rotation around another axis than Ω.
Settling in inclined shear flow. Eq. (62) is valid for any orientation of the forcing relative to the shear flow, described for instance by two angles relative to the vorticity and flow directions. In the remainder of this discussion we focus on the concrete example of a particle settling under gravity, F ext = mg, with the particular set of orientations so that gravity lies in the plane spanned by the vorticity axis and the flow direction, see Fig. 1. This situation corresponds to settling between two far-apart shearing walls, parallel to the walls, but where the shearing is at an angle to gravity. We denote by ϕ the angle between g and the vorticity, see Fig. 1.
We show the resulting settling velocity as a function of Wi in Fig. 2, for ϕ = 0, 45 • , and 90 • . When ϕ = 0 we recover the cross-shear result, and the lateral drift vanishes. As the angle of inclination increases the settling velocity increases, diminishing the shear-induced drag increase described for ϕ = 0 in Refs. [1,3,7]. When gravity acts along the flow direction, ϕ = 90 • , the settling velocity is almost the same as that given by Stokes law, only slightly higher. The direction of the O(Wi) lateral drift is along the −ŷ direction, see Fig. 1. For finite ϕ and small Wi this drift is the dominant feature of the particle velocity. But even for larger Wi the magnitude of this drift is comparable to the reduction in settling velocity when ϕ = 45 • . The additional O(Wi 2 ) drift is in the third independent direction, given by F ext ×ŷ (Fig. 1). For Wi ≈ 0.5 it is comparable in magnitude to both the reduction in settling velocity and the O(Wi) drift.
The direction of the O(Wi) lateral drift can be understood by considering how the elastic fluid is stretched in the vicinity of the sphere. In Fig. 3 we show the trace of the first order elastic stress tensor, Tr Π (1) (r), around the sphere. This trace indicates how strongly the dumbells are stretched by the lowest-order Newtonian flow. The Figure shows a center cross-section of the particle, viewed along the direction of external forcing. In the cross-shear flow, ϕ = 0, the stretching is a complicated function of the spatial variables, but perfectly symmetric around the sphere (top row in Fig. 3). Therefore there is no net force on the particle. But as the external forcing is tilted, the particle is forced to move along the flow direction. Now the particle surface moves opposite to the undisturbed flow on one side, and along the undisturbed flow on the other. This asymmetry results in different stretching of the dumbbells on the two sides, as illustrated for ϕ = 45 • in the second row of Fig. 3. This stress contributes to the particle drift both directly, by forcing the particle surface, and indirectly by forcing the suspending fluid and thereby inducing additional viscous drag.
Method. In this paper we also introduced the tensors T nl i1i2..in (Section III). These tensors are a basis suitable for symbolic calculations of tensorial quantities in spherical geometry. In particular they allow us to write down solutions to inhomogenous Stokes equations in tensorial form, without any explicit coordinate representation, and without explicitly solving differential equations. The calculation in this paper demonstrates the power of our method for treating tensorial equations such as the coupled rank-2 constitutive Eq. (3) and the rank-1 flow Eq. (2). Nevertheless, there are many open questions regarding the algebraic properties of the T -tensors. Most importantly, we have shown that the product T l1l1 T l2l2 is given by a linear combination of T JJ with |l 1 − l 2 | ≤ J ≤ l 1 + l 2 , analogous to the product of two spherical harmonics (see Appendix A). But further work must be done to determine the properties of the coefficients in this linear combination, to determine the general expression for differentiation, and to prove the general case of Eq. (25).
Our tensor formalism can also be extended to other geometries. Nearly spherical geometry can be treated by perturbation theory. Other geometries can be treated by the method of images [22,23]. For example, the flow around a spheroid in unbounded flow is given by a finite distribution of multipoles [23]. However, the radial functions are no longer simply r m , but integrals I n m = c −c ξ n /|r − ξn| m dξ, where n is the direction of the spheroid and c is a shape-dependent constant [23,24]. Their algebraic properties must be worked out in order to use the formulae in Section III e.g. for the Fourier transform. For wall interactions, or other many-center problems, it is possible to derive a translation theorem that expresses |r − r | m T nl (r − r ) as an infinite series of |r| m T nl (r) and |r | m T nl (r ) [20], which restores the linearity of the problem. (A9)