Shaken and stirred: Random organization reduces viscosity and dissipation in granular suspensions

See allHide authors and affiliations

Science Advances  30 Mar 2018:
Vol. 4, no. 3, eaar3296
DOI: 10.1126/sciadv.aar3296


The viscosity of suspensions of large (≥10 μm) particles diverges at high solid fractions due to proliferation of frictional particle contacts. Reducing friction, to allow or improve flowability, is usually achieved by tuning the composition, either by changing particle sizes and shapes or by adding lubricating molecules. We present numerical simulations that demonstrate a complementary approach whereby the viscosity divergence is shifted by driven flow tuning, using superimposed shear oscillations in various configurations to facilitate a primary flow. The oscillations drive the suspension toward an out-of-equilibrium, absorbing state phase transition, where frictional particle contacts that dominate the viscosity are reduced in a self-organizing manner. The method can allow otherwise jammed states to flow; even for unjammed states, it can substantially decrease the energy dissipated per unit strain. This creates a practicable route to flow enhancement across a broad range of suspensions where compositional tuning is undesirable or problematic.


Densely packed suspensions arise widely in industry and manufacturing, where reliable, predictable, and prescribable flow properties are essential (1). A major limiting factor in their processability is the very steep increase of viscosity upon increasing the volume fraction of solid material φ toward the jamming transition (2). This is particularly evident in the non-Brownian regime (particle size, ≥10 μm), where frictional particle contact interactions reduce the jamming density φm (35) and increase dissipation in process flows, resulting in high energy costs.

Empirical strategies that reduce the viscosity and/or net dissipation include tuning the physical properties of the particles—for example, their size, shape, and polydispersity—or modifying their interactions through chemical additives known as plasticizers, emulsifiers, or friction modifiers. These lubricate interparticle contacts and reduce the suspension viscosity by raising φm. Often, though, end-use requirements leave little room for maneuver in the formulation. In calcium phosphate cements for bone injection (6, 7), for example, chemistry and biology both constrain the use of molecular additives. There is, therefore, a practical need for methods of dense suspension flow control that do not require changes to formulation.

Two recent experiments suggest a possible route toward this goal, achieving driven viscosity reduction by superimposing an oscillatory cross shear (OCS) on a primary desired flow. Using OCS, Blanc et al. (8) demonstrated a twofold increase in the sedimentation velocity of an intruder in a granular suspension of rate-independent rheology, whereas Lin et al. (9) measured a two-decade viscosity drop in one of shear-thickening rheology. The latter effect was argued to be a consequence of the fragility of shear-induced particle contacts (9, 10), suggesting that good flowability might be achieved only when the OCS is sufficiently fast to keep the microstructure in a load-incompatible state (10, 11). In this limit, the reduction in primary flow viscosity (unless this is infinite) might easily be outweighed by the high energy cost of implementing fast OCS, particularly for rate-independent suspensions (8) whose primary viscosity drop is much less than in shear-thickening ones (9). More generally, it is not clear how far the benefits of OCS depend on the underlying suspension rheology: The short-range repulsions that prevent frictional particle contacts at low stresses in thickening suspensions (1214) may or may not play a major role during OCS-assisted viscosity reduction.

Here, we present numerical simulations showing that the viscosity drop induced by OCS is generic to suspensions with friction-dominated stress. This includes noninertial flows of most dense suspensions of super-micrometer–sized particles (5). The transverse flow oscillations directly inhibit particle contacts without requiring short-range repulsions, enhancing lubrication and shifting φm to higher values. Consequently, the viscosity reduction increases with increasing φ, so that near jamming the saving in primary flow dissipation outweighs the cost of OCS at any primary flow rate, giving a net reduction in the energy expended per unit strain in the primary direction. We then show that the reduction in particle contacts stems from an OCS-induced “random organization” mechanism (1517); this leads us to an enhanced version of the flow protocol that can reduce the dissipation further. Guided by these results, we argue that driven viscosity control should extend flowability and reduce the associated energy cost across a broad class of materials, including slurries, muds, cement, and other immersed granular systems.


We study a suspension of nearly hard, athermal spheres subject to short-range hydrodynamic and contact interactions with static friction coefficient μs as described in Materials and Methods below. This numerical model [and similar ones (18, 19)] is known to yield accurate predictions for the rheology of non-Brownian hard sphere suspensions. The suspension shows rate-independent rheology, well described under steady simple shear by the viscous number formalism [see the study by Boyer et al. (20) and the Supplementary Materials]. A snapshot of the simulated system is shown in Fig. 1A (i).

Fig. 1 Viscosity and dissipation reduction under superimposed primary and oscillatory flows.

(A) (i) Simulation snapshot showing primary (blue) and cross shear (red) flow directions. (ii and iii) Example flow paths explored for different values of parameters ωpri and ω (values given in insets) with δ = 0. (B) Contour map showing viscosity in primary flow direction as a function of ωpri and ω, for δ = 0 and φ = 0.55. (C) Viscosity as a function of oscillation rate ωγ/γ̇ at various volume fractions φ with amplitude γ = 1% and ωpri = 0, the SO protocol. (D) Viscosity divergence as a function of φ under steady shear (SS) and high-frequency SO with friction coefficient μs = 1. Inset: Difference between SS and SO viscosities. (E) Viscosity divergences for particles with lower friction coefficient μs show diminishing viscosity reduction. (F) Dissipation per unit strain W (rescaled by the φ-dependent steady shear dissipation WSS) as a function of oscillation rate for the same simulations as in (C), for (i) φ = 0.54 and (ii) φ = 0.57, showing contributions in xy and zy. Green areas in (i) and (ii) highlight the region in which both viscosity and dissipation reduction are achieved. (iii) Dissipation in the (ωγ/γ̇, φ) plane, highlighting in white the region for which dissipation may be reduced by at least 5% with SO.

Manipulating suspension viscosity using superimposed oscillations

From the argument that fragility makes contact stresses in suspensions susceptible to driven perturbations (9), it follows that the addition of any arbitrary oscillating flow might lead to viscosity reduction. This hypothesis is in line with experimental (21, 22) and theoretical (23) works that propose applied and endogenous noise, respectively, as sources of opening and closing granular contacts and consequent unjamming. To test this, we first explore a generalization of OCS comprising primary steady shear with rate γ̇ and superimposed oscillatory shears in both the primary and cross shear directions, leading to an overall strain in xy as γpri(t)=γsin(ωprit+δ)+γ̇t and in zy as γOCS(t) = γ sin (ωt). For simplicity, we keep γ = 1% in each case, which Lin et al. (9) found to be an optimal amplitude for viscosity reduction. The remaining dimensionless control parameters are then ωpriγ/γ̇, ωγ/γ̇, and the phase shift δ. This protocol gives strain paths such as those illustrated in Fig. 1A (ii and iii). A characteristic viscosity is computed as ηr=σxy/ηγ̇ averaged over 10/γ̇ time units, where η is the solvent viscosity and σxy is the xy component of the stress. We find that δ has very little effect on the viscosity (see the Supplementary Materials), and present a contour map of ηr in the (ωpriγ/γ̇, ωγ/γ̇) plane at δ = 0 and φ = 0.55 in Fig. 1B. At fixed δ, viscosity minima are obtained as ωpriγ/γ̇0 and ωγ/γ̇6. In this limit, that is, with cross shear oscillations only, we obtain a viscosity drop comparable to that for a thickening suspension (9). This suggests that OCS—by keeping frictional particle contacts open—effectively brings the suspension to a low-friction state. We similarly find a maximal rate of viscosity reduction when ωγ/γ̇ is close to unity. Contrary to the hypothesis made above, however, our results show that the orientation of the oscillatory flow is crucial: At this strain amplitude, any oscillatory component along the primary flow direction makes no useful contribution to improving flowability. If the shear is constrained to a single direction and the amplitude of the oscillations is very small compared to the primary flow, then the net displacements of the particles over large strains are, for rate-independent flow, the same as for steady shear. This is not the case when the oscillations are applied transverse to the primary flow.

Simple OCS: Viscosity reduction using transverse oscillations

In what follows, we therefore revert to the purely transverse case with ωpri = 0, leading to γpri(t)=γ̇t, and γOCS(t) = γ sin (ωt) (see Fig. 2A, inset), hereafter called the “simple OCS” (SO) protocol. (This is to distinguish it from an alternative protocol introduced below.) In Fig. 1C, we report the viscosity ηr under this protocol at γ = 1% as a function of the reduced frequency ωγ/γ̇, and in Fig. 1D we compare, as a function of volume fraction φ, the steady shear viscosity (obtained when ωγ/γ̇=0) to the limiting viscosity under SO (obtained when ωγ/γ̇10). The viscosity drop increases rapidly with φ, reaching a decade at φ = 0.56 and actually diverging between φ = 0. 57 and 0. 58 (Fig. 1D, inset). This reveals that as well as reducing the viscosity, the effect of SO is to slightly delay the jamming transition from φm ≈ 0.58, for steady shear to φm,SO ≈ 0.60, at ωγ/γ̇=10. Although small in absolute terms, shifts of jamming by a couple of percent can have dramatic consequences for formulation and processing (24), as discussed further below.

Fig. 2 Revealing random organization at work during oscillatory shear.

(A) Origin of the viscosity drop at φ = 0.54 with SO. The contact stress contribution is strongly suppressed with increasing oscillation rate. Inset: Schematic of the SO strain profile. (B) The cumulative pair correlation function G(h) under steady shear and ωγ/γ̇=10, demonstrating a room-making process. Inset: Around 20 cycles are needed to minimize the number of particle contacts. (C) Proportion of particles following an irreversible trajectory under successive periods of oscillatory shear (in the absence of primary shear) as a function of the number of cycles, starting from presheared configurations at several volume fractions. The steady decrease of the irreversibility is a signature of random organization.

This shift of jamming under SO naturally raises the question of the sensitivity of the viscosity reduction to particle friction, which may stem from, for example, surface roughness (25). It is expected that in the absence of ordering, which we do not observe in our binary system, the random close packing density φRCP ≈ 0.64 is an upper limit for both φm and φm,SO. Furthermore, it is established that the jamming point φm approaches φRCP as surface friction μs is decreased (26, 27). Because φm < φm,SO, it follows that φm < φm,SO < φRCP and the window between φm and φm,SO consequently vanishes in the limit of low friction (as μs → 0). The performance of SO, thus, diminishes as friction decreases. We demonstrate this in Fig. 1E in the limits of steady shear and SO with ωγ/γ̇=10. As a result, suspensions of rough particles, which are typically the most problematic in terms of processing (18), are best placed to benefit from driven flow control. In the context of friction-driven shear thickening of colloids, this result thus confirms that SO can be successful in reducing the viscosity of a thickened sample, as demonstrated by Lin et al. (9), but that it would fail to reduce the viscosity of a nonthickened sample, that is, one at which the applied stress lies below the onset stress (5, 12, 28, 29).

SO-enabled reductions in energy dissipation

The ability to control suspension viscosity during flow is itself desirable for mitigating instabilities (30) and, for example, when pumps are desired to operate within narrow bounds. Often, though, rheological tuning has a somewhat different objective: to minimize the energy cost of processing. For 0.58 < φ < 0.60, OCS triumphs: It permits flow at finite dissipation rates not otherwise possible. Below jamming (φ < φm) however, its benefits are less obvious. The energy dissipation is given per unit volume and per unit primary strain as W=limT(0Tdtσ:γ̇)/(γpri(T)γpri(0)). Figure 1F (i to iii) shows this quantity (rescaled by the φ-dependent steady shear dissipation WSS) as a function of oscillation rate for our SO protocol, separating out the primary (σxyγ̇) and cross flow (σzyγ.OCS) contributions. The primary dissipation decreases in line with the viscosity, but the direct cost of the cross shear increases as (ωγ/γ̇)2. Summing these, we identify oscillation rates for which W is usefully decreased, highlighted green in Fig. 1F (i and ii) and outlined in white in Fig. 1F (iii). This operating window, although it grows as the density approaches φm, remains narrow at lower densities. We show below that it can be extended significantly by a simple modification to the oscillatory protocol.

Random organization drives the viscosity reduction

The proposed modification exploits mechanistic insights, gleaned from our simulations, into how OCS promotes flowability. To gain these insights, we start by decomposing the viscosity into its hydrodynamic and frictional particle contact contributions, revealing that at φ = 0.54, the stress is dominated by friction for any oscillatory frequency (Fig. 2A). Significantly, the effect of the cross shear oscillations is to decrease this frictional part while leaving the hydrodynamic part unchanged. The loss of friction parallels the shift of jamming to higher φ (Fig. 1D), indicative of a shift from rolling to sliding contacts as ωγ/γ̇ is increased (5, 29). Defining interparticle gaps hij = 2(rijaiaj)/(ai + aj) with rij as the center-to-center distance between particles i and j with radii ai and aj, respectively, we compute G(h), the average number of neighbors around a particle separated at most by h (Fig. 2B). The loss of frictional particle contacts occurs by a “room-making” process, whereby the mean distance between nearest neighbors increases. Consequently, starting with a presheared sample, there is a gradual decrease of particle contacts over O(10) cycles after SO startup (Fig. 2B, inset). Strikingly, room-making does just enough to hinder the stress-generating contacts. This is strongly reminiscent of “random organization,” whereby application of oscillatory shear to a suspension of hard particles drives collective self-organization, leading to configurations that minimize the number of particle contacts generated per cycle (1517). Below a critical volume fraction, the system evolves to an “absorbing state” for which configuration invariance is ensured under further oscillations. Although first elaborated for dilute suspensions, a similar scenario applies at higher density where the absorbing state is defined not by absence of collisions but by absence of plastic rearrangements (31, 32).

Flow-induced random organization offers a natural explanation for the viscosity decrease upon increasing oscillation frequency. At high frequency, the “primary” and “secondary” labels respectively assigned to steady shear and OCS are misnomers. In fact, we have a steady transverse flow that weakly perturbs an oscillatory flow, for which the random organization effect is well established. This makes room around particles, thus decreasing the contact stress, whereas the steady shear slowly consumes this room and simultaneously initiates new particle contacts. At finite primary flow, the absorbing state can never be reached, but its proximity allows particles to avoid frictional particle contacts at densities where these would otherwise cause large viscosities or jamming.

The largest φ at which absorbing states are obtained locates a nonequilibrium phase transition, where the self-organization process is maximized (15, 16, 33). To confirm the role of random organization, we determine the location of this transition in our system. Starting from a presheared configuration, we apply an oscillatory shear at γ = 1% with no primary flow and measure the fraction of particles following irreversible, “active” trajectories, which, for these purposes, we define as those whose net displacement after a cycle stays below a threshold of 10−5ai. For φ ≤ 0.58, this quantity approaches zero, indicating that the system evolves toward absorbing states (Fig. 2C). For an amplitude γ = 1%, the absorbing state transition density is thus estimated as 0.58 < φ < 0.59. Because in practice the transition is cut off by the primary flow, our precise definition of activity is not crucial here, although a more inclusive one (for example, counting all particles that make frictional contact at any point during the cycle) would give a lower estimate for the transition. Nonetheless, our results indicate that the random organization effect is a strong one throughout the density range where our SO protocol is effective.

Alternating OCS: Separated flow phases reduce dissipation further

In the SO protocol above, particle contacts are eliminated by applying OCS concurrently with the desired primary shear. A relatively high energy cost arises from the need to have sufficiently fast oscillations to ensure that random organization can compete with the restoration of frictional particle contacts caused by the primary shear. If this is indeed the mechanism, though, there is no strict requirement that we perform these flows concurrently. Instead, we can use alternating intervals of OCS without primary shear and of primary shear without OCS. The former eliminates particle contacts; the latter restores them, but not before a finite strain has been achieved. The cross shear dissipation can, in principle, be reduced to zero by having long intervals of very slow oscillations, creating an optimization scenario different to that of SO.

We therefore now test a new flow protocol [“alternating OCS” (AO)] that alternates an interval of n periods of oscillation during a time αT with γ̇pri=0 and γOCS(t) = γsin (ωt) for ω = 2πn/(αT), with an interval of primary shear during a time (1 − α)T with γ̇pri=γ̇/(1α) and γ̇OCS=0. Guided by our result in Fig. 1B, the oscillations are applied transverse to the primary flow, where we anticipate that their efficacy will be maximized. Averaged over one cycle T, the shear rate in the primary direction is γ̇. The primary shear strain during each cycle is Γ, that is, T=Γ/γ̇.

In the case of rate-independent dynamics as simulated above, the microstructure depends on the strain path only (sketched in Fig. 3A, inset), not the rate at which it is followed. As a consequence, the viscosity depends on n, γ, and Γ, but not on α (as long as 0 < α < 1). We define the relative viscosity as ηr=σxy/(ηγ̇pri). This viscosity is averaged over the intervals of pure primary shear during the AO protocol, measured over a time period covering 30 strain units in the primary direction. It is reported as a function of n in Fig. 3A for φ = 0.56 and γ = Γ = 1%. The viscosity drops rapidly with n, and, remarkably, even n = 1 is sufficient to achieve a viscosity reduction of ≈ 96% in this unjammed system. This viscosity reduction is already larger than that achieved with SO, as seen by comparing the relative reductions in Fig. 3B (AO) and the inset of Fig. 1D (SO). Finally, we find that the viscosity drop is maximized with AO [as with SO (9)] for γ, Γ ≈ 1 to 5%, whereas for larger γ, interparticle gaps close during cycles, allowing frictional particle contacts and a rapid viscosity increase. A viscosity transient for the AO protocol is given in the Supplementary Materials.

Fig. 3 Optimizing dissipation reduction using an alternating OCS (AO) protocol.

(A) Primary viscosity as a function of n, the number of cross shear oscillations per phase of the AO protocol, with γ = Γ = 1% and φ = 0.56. Inset: Schematic of the AO strain profile. (B) Relative viscosity drop as a function of volume fraction φ, comparing steady shear and AO. (C) Energy dissipation for the cases n = 1, 2, 3, 5, and 10 as a function of the fraction of time spent oscillating α. (D) Dissipation in the (ωγ/γ̇, φ) plane rescaled by the φ-dependent steady shear dissipation, highlighting the region for which dissipation is reduced by at least 5% with AO. Outlined in black is the analogous region from Fig. 1F (iii) for comparison. (E) Comparison of the minimal dissipation obtained with the AO and SO protocols across a range of volume fractions φ. Inset: The ratio between the minimal dissipation rates achievable for AO and SO, WminAO and WminSO, respectively.

Although α has no role in setting the viscosity, it is a crucial parameter when it comes to the dissipation, which depends on the deformation rate. In Fig. 3C, we show the work per unit primary strain W (rescaled by the φ-dependent steady shear dissipation WSS) as a function of α for several values of n, for φ = 0.56. For n = 1 and for α values between 0.58 and 0.9, the overall dissipation is reduced compared to steady shear, reaching a reduction of around 45% at α ≈ 0.8. Dissipation is minimized when n = 1 for all φ, as the shear rate shoots up quickly with increasing n, swamping any further viscosity reduction achieved for n > 1. To ease the comparison with the SO protocol, we define an oscillatory frequency for AO as ω = 2πnγ/Γ (that is, a given frequency corresponds to the same number of cross shear oscillations per unit strain in the primary direction for SO and AO). In Fig. 3D, we show a map of the relative dissipation W/WSS in the (φ, ω) plane, showing a wider area of reduced dissipation with AO compared to SO [Fig. 1F (iii)]. The AO protocol therefore has clear advantages particularly in avoiding the need to precisely tune the frequency of the driving oscillations. We finally present in Fig. 3E a comparison of the reduction in energy dissipation achieved by the AO and SO protocols as a function of the volume fraction φ. For φ ≲ 0.54, there is no further gain with AO compared to SO, but for larger volume fractions, where the viscosity reduction performance of AO is markedly superior (Fig. 3B), there is an improvement in AO over SO. We quantify this improvement in Fig. 3E (inset), giving the ratio of the minimal dissipation for AO and SO (WminAO/WminSO). This shows that close to the steady shear jamming volume fraction, the improvement of AO over SO can reach almost 40%, suggesting that AO is likely the protocol of choice for dissipation reduction in very concentrated frictional suspensions.

Concluding remarks

Our results show that nonsteady deformation protocols can lead to substantial viscosity and energy dissipation reductions in any friction-dominated suspension flow. The strategy is applicable for most flows involving granular suspensions and related systems in which frictional particle contacts bear most of the stress in steady shear, including Brownian suspensions under very large stresses (34). Because of their simplicity, our protocols, or ones like them, might be readily implemented as precision unblockers and flow controllers in industrial devices, such as extruders or mixers, or as dissipation regulators for active granular damping (35). In particular, an extrusion nozzle might be fitted with an internal coaxial cylindrical actuator that oscillates about its axis with a protocol specified to maximize flowability according to our present results. Moreover, these implementations might be applied not only to minimize viscosities but also to regulate them against a desired set point. We tested such a protocol numerically with good success (see the Supplementary Materials). From a fundamental point of view, the relation to random organization opens new research directions. For example, it suggests that protocols other than oscillatory flow that lead to a similar absorbing phase transition (36) might also be good candidates for driven flow enhancement in complex fluids. It also suggests an unexpected link between rheological properties and hyperuniformity (37, 38).


We simulated the trajectories of athermal, noninertial particles using a minimal model that comprises short-ranged hydrodynamic lubrication and frictional surface contacts. Our simulations comprise O(103) particles with size ratio 1:1.4 in a periodic box. For a particle pair with positions x1 and x2 and translational and rotational velocities U1, U2 and Ω1, Ω2, respectively, in a background flow described at x1 by U(x1) = Ex1 + Ω × x1, the hydrodynamic forces F1h, F2h and torques Γ1h, Γ2h are given by (3941)[F1hF2hΓ1hΓ2h]=RLub[U(x1)U1U(x2)U2ΩΩ1ΩΩ2EE]+RStokes[U(x1)U1U(x2)U2ΩΩ1ΩΩ2](1)The matrices RLub and RStokes follow our earlier description (28), with the scalar resistances therein comprising only the leading short-ranged diverging contributions, following the study by Ball and Melrose (42). The hydrodynamic stress contribution for particle 1 resulting from its pairwise interaction with particle 2, with force F1h and particle-particle vector r is given by Sh=12(F1hrT+(F1h)Tr).

The leading terms of RLub diverge according to 1/h as particles 1 and 2 approach, with h as the surface-surface distance. Following experimental evidence that lubrication layers break down in suspensions under large stress (43), and, equivalently, for large particles (5), we used a minimum hmin = 0.001a (with a the smaller particle radius), below which hydrodynamic forces are regularized and particles may come into contact. For a particle pair with contact overlap δ and center-center unit vector n, we computed the contact force and torque according to (44)F1c=knδnktu(2a)Γ1c=a1kt(n×u)(2b) where u represents the incremental tangential displacement, reset at the initiation of each contact. kn and kt are stiffnesses and a1 is the radius of particle 1. The tangential force component is restricted by a Coulomb friction coefficient μs, such that |ktu| ≤ μsknδ. For larger values of |ktu|, contacts enter a sliding regime. The contact stress contribution is given by Sc=F1crT for particle-particle vector r and pairwise force F1c. The stress tensor is σ=2ηE+1V(Sh+Sc), where η is the suspending fluid viscosity and the sums are over all relevant pairwise interactions. Throughout the main text, we focused on the shear component in the primary flow direction, σxy.

Trajectories were computed from the above forces using two equivalent schemes. In the first, contact and hydrodynamic forces and torques were summed on each particle [according to the study by Radhakrishnan (45)] and the trajectory was updated according to Newtonian dynamics [using LAMMPS (46)], ensuring the Stokes number (ργ̇a2/η for particle density ρ, suspending fluid viscosity η, and shear rate γ̇) remained ≪ 1 to approximate overdamped conditions. We also set 2γ̇a/kn/(2ρa)<105 to approximate hard spheres. In the second, per-particle forces were explicitly set to zero and the velocities were computed to balance contact and hydrodynamic forces and torques, ensuring strictly inertia-free flow (12, 28). The numerical model generates results consistent with μ(J)-rheology as predicted by the experimental work of Boyer et al. (20) (see the Supplementary Materials).


Supplementary material for this article is available at

section S1. Demonstration that our numerical model predicts μ(J) rheology

section S2. Influence of primary-cross shear phase difference

section S3. Time evolution of the relative viscosity for the SO and AO protocols

section S4. Setting the viscosity using a proportional flow controller

fig. S1. Comparison of the results from the present model to viscous number rheology.

fig. S2. Viscosity maps analogous to Fig. 1B but with varying phase shift δ.

fig. S3. Time evolution of the viscosity for the SO and AO protocols.

fig. S4. Time evolution of the relative viscosity under proportional control in the SO protocol.

This is an open-access article distributed under the terms of the Creative Commons Attribution license, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


Acknowledgments: C.N. thanks R. Radhakrishnan for sharing his derivation of the hydrodynamic forces. Funding: C.N. acknowledges an Engineering and Physical Sciences Research Council Impact Acceleration Account and grant EP/N025318/1 and, subsequently, the Maudslay-Butler Research Fellowship at Pembroke College, Cambridge, for financial support. M.E.C. is funded by the Royal Society. Author contributions: All authors planned the work, interpreted the results, and wrote the manuscript. C.N. and R.M. conducted the simulations. Competing interests: The authors declare that they have no competing interests. Data and materials availability: Codes required to generate the results in this article are given at and All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors

Stay Connected to Science Advances

Navigate This Article