Geometric frustration induces the transition between rotation and counterrotation in swirled granular media
Lisa Lee, John Paul Ryan, Yoav Lahini, Miranda Holmes-Cerfon, Shmuel M. Rubinstein
GGeometric frustration induces the transition between rotation and counterrotation inswirled granular media
Lisa Lee, John Paul Ryan, Yoav Lahini, Miranda Holmes-Cerfon, ∗ and Shmuel M. Rubinstein ∗ John A. Paulson School of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138 Department of Computer Science, Cornell University, Ithaca, NY 14850 Raymond and Beverly Sackler School of Physics and Astronomy, Tel Aviv University Courant Institute of Mathematical Sciences, New York University, New York, NY 10012 (Dated: December 11, 2018)Granular material in a swirled container exhibits a curious transition as the number of particlesis increased: at low densities the particle cluster rotates in the same direction as the swirlingmotion of the container, while at high densities it rotates in the opposite direction. We investigatethis phenomenon experimentally and numerically using a co-rotating reference frame in which thesystem reaches a statistical steady-state. In this steady-state the particles form a cluster whosetranslational degrees of freedom are stationary, while the individual particles constantly circulatearound the cluster’s center of mass, similar to a ball rolling along the wall within a rotating drum.We show that the transition to counterrotation is friction-dependent. At high particle densities,frictional effects result in geometric frustration which prevents particles from cooperatively rollingand spinning. Consequently, the particle cluster rolls like a rigid body with no-slip conditions on thecontainer wall, which necessarily counterrotates around its own axis. Numerical simulations verifythat both wall-disc friction and disc-disc friction are critical for inducing counterrotation.
From hurricanes to bacterial swarms, the emergence ofsystem-scale circulation from local interactions and lo-cal driving is a phenomenon exhibited on many scalesand in many different physical systems. In 2D turbu-lence, vorticity at the small injection scale may cascade tolarger and larger scales, stabilizing into a single system-scale vortex [1, 2]. Analogous behavior is observed inan active fluid of spinners flowing though a lattice ofannular channels and driven by a magnetic field, whichbreaks time reversal symmetry. This active liquid de-velops sound modes that propagate along the boundary,generating global circulation [3]. Finally, dense suspen-sions of self-propelled bacteria in confinement also lead tothe spontaneous formation of stable circulation along thecontainer walls [4, 5]. In all of these systems, rotation isactively or passively injected locally, and interactions be-tween the local units ultimately lead to global circulationof the entire system.There is no need to go to very complex systems toobserve the nontrivial emergence of system-size circula-tion. A handful of marbles swirled in a teacup exhibitssimilar dynamics. When there are only a few marblesin the container, they form a line that rolls along thecontainer wall at the frequency of the circular transla-tions ("snake" mode [6]). When a couple more marblesare added, they form a cluster that sloshes periodicallyaround the container, its individual marbles repeatedlycramming against the wall at the outer edges of the swirlbefore flowing freely through the rest of the container.This sloshing marble ensemble rotates about its own cen-ter of mass in the same direction as the container, muchlike wine swirling in a glass. When even more marbles areadded, the monolayer of marbles will continue to period-ically slosh around the container, but, curiously, above a critical marble density the pack reverses to counterrota-tion, or to rotating in the direction opposite to that ofthe container [7]. This transition to counterrotation isstartling, because the angular momentum of the marblesin their own reference frame changes sign with the addi-tion of a couple marbles, and is ultimately opposite thatof the hand which forces them. The transition is not re-stricted to marbles in a teacup, but can be observed in avariety of shaken containers filled with granular material,from vibration mills for grinding industrial materials [8]to baby rattles.Such a simple and ubiquitous phenomenon demandsan explanation, yet, while it has been studied empirically[6, 7, 9–11]; there is still no clear understanding of theminimal ingredients necessary for it to occur. We proposeto understand this phenomenon via analogies to a singlerigid body and a fluid in a swirled container. A singlerigid body in a swirled container with no-slip boundaryconditions will roll commensurately on the wall (akin toone gear inside another), rotating about its own centerof mass in the opposite direction as the container. Onesees similar behavior for a pancake rolling on the edgeof a swirled frying pan. In contrast, a rigid body withperfect-slip boundary conditions (as if the pancake werein a bath of melted butter) would receive no torque fromthe boundary, and hence experience no overall sense ofrotation. Finally, a fluid will rotate in the same directionas the container is swirled, also to varying degrees basedon its slip conditions with the container wall. Thus therotation behavior of a single body in a swirled containerdepends on its own rigidity as well as its slip interac-tions with the wall. For granular systems composed ofmany interacting bodies, the fluidity and boundary in-teractions are not material constants but are emergent a r X i v : . [ c ond - m a t . s o f t ] D ec dynamic properties [12–16], determined by the frictionbetween the particles as well as the friction with the con-tainer wall.Here we experimentally and numerically investigatethe dynamics of swirling particles with the goal of iden-tifying the minimal physical factors that facilitate thetransition from rotation to counterrotation as the num-ber of particles is increased. We show that the criti-cal control parameter for counterrotation is the particlecluster’s effective slip with the wall, which arises fromfriction. Although the individual particles’ friction coef-ficient with the wall does not depend on the number ofparticles, the effective boundary condition of the particlecollective does and can change from perfect slip to perfectstick as the number of particles is increased. Our exper-imental observations suggest that friction between theparticles prevents densely packed particles from coopera-tively rolling and spinning, causing the particle cluster torigidify and roll on the wall as a whole, ultimately result-ing in counterrotation. Finally, we numerically simulatethe swirling system and find that if either interparticlefriction or particle-wall friction is set to zero, the systemrotates at the same frequency at all densities and nevertransitions to counterrotation, verifying that friction iscritical for counterrotation.Experimentally, the transition from rotation to coun-terrotation of swirling particles is observed upon aRotoMix orbital table performing circular translations,without rotation, of angular velocity ω = 11 . rad/secand amplitude 1.15 cm. The mounted circular container,5.1 cm in diameter, with a vertical edge is partially filledwith a monolayer of N (ranging from 24 to 48) plasticspheres (mass m = . g and diameter 6 mm) and imagedfrom above using a stationary Sony RX100 IV camera at960 fps, as shown for two typical images in Fig 1a. Wecall the collection of particles a “cluster,” and measureits average angular velocity (cid:36) by the average angular ve-locity of the particles about the cluster’s center of mass: (cid:36) = T − N (cid:80) Tt =2 (cid:80) Ni =1 (cid:36) i ( t ) . Here T is the total num-ber of frames and (cid:36) i ( t ) is particle i ’s instantaneous an-gular velocity about the cluster’s center of mass at time t : (cid:36) i ( t ) = t (arctan( y ( t,i ) − y CM x ( t,i ) − x CM ) − arctan( y ( t − ,i ) − y CM x ( t − ,i ) − x CM )) .Here x CM and y CM refer to the x- and y-coordinates ofthe center of mass of all particles.Under these conditions, the average angular velocity (cid:36) crosses zero and transitions between rotation and coun-terrotation at roughly 36 particles, as shown in Fig 1band the SI video. For N > , the cluster counterrotates( (cid:36) < ), and for even higher values of N , the averageangular velocity reaches a minimum and then increasesslightly for the largest values of N .Further analyzing the data in the lab frame is a chal-lenge since the particles experience a rotating time-dependent force from the circular translations of theswirling container; thus the particles slosh around the Figure 1. a: Raw image of particles in a swirling container,and the same container at a later time after it has trans-lated through roughly half its circular trajectory. The whitedotted line represents the outermost points of the containerduring its trajectory. The container does not rotate in the labframe. b: The mean angular velocity ( (cid:36) ) of the particle clus-ter about its own center of mass for different particle counts N . As N increases the cluster transitions from rotation tostalling to counterrotation. Black dashed line indicates theangular velocity of a theoretical single particle or pancake,with perimeter matching that of the cluster, if it were per-fectly rolling along the wall. c: The M-frame rotates withthe container such that the point (left, labeled "M") on thecontainer’s boundary furthest from the center of swirling isalways positioned at the bottom in the M-frame (middle andright). As a result, the walls of the container in the M-frameappear to be rotating clockwise, as depicted by the black ar-rows. Also shown are the point S center of swirling, and thepoint C center of the container. The rightmost M-frame dia-gram depicts the centrifugal force field (F cent , brown) as wellas the Coriolis force (F Cor , red) associated with a particle(green) with indicated velocity ( v , blue). container with no apparent steady-state (SI Video).Therefore, it is illuminating to analyze the data in aframe of reference where the external forces are steady[10]. Consider a frame which rotates at the angular veloc-ity ω of the container, about an axis at the center S of theswirling orbit. In this frame of reference the container’stranslational velocity is zero and it rotates around its owncenter C at a constant angular speed − ω . This frame ofreference is equivalent to centering our camera above thecenter of the container and rotating it such that S is sta-tionary, keeping the boundary point furthest from S atthe bottom; Hereafter we refer to this point as the pointM, as shown in Fig 1C. In this frame of reference there aretwo external forces acting on all particles: the centrifugalforce, which points radially away from S with magnitudeat a point x equal to mω | x − S | , and the Coriolis force,which points perpendicularly to the right of the parti-cle’s velocity vector v , with magnitude m | v | ω . Withinthe container, the centrifugal force pushes the particlesoutward to the boundary, and at most of the boundarythe particles are then pushed towards point M. At theboundary, the particles are also dragged clockwise dueto friction with the moving container wall. We call theoriginal frame of reference the lab frame, and the rotatingframe of reference the M-frame. The latter terminologyis borrowed from Kumar et al [10], which used this frameof reference to study granular particles on the edge of aswirling cylinder.In the M-frame, our system appears to be at steadystate (SI Videos and Fig 2a), characterized by a steadyparticle density pattern and a steady circulation patternwhen averaged over typical fluctuations, as shown in Fig-ures 2b and 2c. Visually, the dynamics resemble particlesin a rotating drum, where the wall rotates at a steady ve-locity and the beads are uniformly pulled downwards. Inthe rotating drum, gravity aggregates the beads at thebottom, where the wall then drags the bead cluster up-wards until it liquefies and the beads avalanche back tothe bottom of the drum, continuing the cycle. In the M-frame, the force of gravity is replaced by the centrifugalforce, both tending to push particles toward point M atthe bottom. Additionally, the M-frame introduces theCoriolis force that the rotating drum lacks, resulting inthe beads tending to the left of the container.In the M-frame, all particles follow clockwise trajec-tories within the container while passing between twodistinct spatial regions: a dense, “pinned” region (solid-like), and a sparse, “loose” region (liquid-like), as shownin Fig 2a. In the pinned region, the particles are packed inlayers against the moving wall and move as a rigid struc-ture, whereas in the loose region the particles performless constrained trajectories as they cross the container,as shown in Fig 2b. When a particle in the loose regiontransitions to the pinned region, it undergoes multiplecollisions before settling into a trajectory parallel to thecontainer’s edge. Increasing N also increases the num-ber of collisions a particle experiences, causing it to enterthe pinned region earlier in its cycle and subsequently in-creasing the size of the pinned region. For all values of N , particles exit the pinned region and detach from therigid structure at the same location. Once loose, a par-ticle experiences centrifugal (0-460 cm/s ) and Coriolis(0-200 cm/s ) forces and accelerates as it traverses theloose region. Multiple collisions at the end of the looseregion pin the particle and this cycle continues.Individual particles circulate between these two re-gions and on average the particle cluster rotates in thesame clockwise direction as the container in the M-frame,shown for two typical examples in Fig 2b and 2c. Impor-tantly, the cluster’s angular velocity in the M-frame, (cid:36) M ,is related to its angular velocity in the lab frame as (cid:36) M = (cid:36) − ω . (1)Therefore, counterrotation in the lab frame ( (cid:36) < ) cor-responds to the particle cluster rotating faster than thedish in the M-frame ( | (cid:36) M | > ω ). Figure 2. a: Sample image in the M-frame, with pinned(red) and loose (blue) region particle centers labeled. A sam-ple angular slice interrogation area is highlighted in green.The angular slice slides around the entire container duringanalysis; the major angular positions θ are denoted on theedge of the container. b: Density histogram of the particlesfor a rotation case at low N = 28 . The white arrows denotethe average local particle velocity deviation from the velocityof the underlying container. c: Same as b, but for counterro-tation at high N = 48 . The coherent pinned region (top lefton the container) grows in size with increasing N , while theloose region shrinks in size. The particles perform clockwisetrajectories around the container while passing between thepinned and loose regions in a cyclic manner. At this point it is insightful to consider a minimalmodel for our swirling system: a single rigid ball of radius R ball swirled in a container of radius R cont (SI Video). Asteady-state for the ball in the M-frame occurs when it isnear the bottom-left of the container (near M), rotatingat constant angular velocity (cid:36) b with the torque alongthe frictional boundary balancing the centrifugal force.If we define a dimensionless parameter γ = | (cid:36) b | ω , thenby (1), γ < corresponds to rotation in the lab frameand γ > corresponds to counterrotation. The value of γ depends on the slip conditions between the ball andthe container wall, which interpolate between two lim-iting cases: perfect no-slip (strong friction) and perfectslip (no friction). When friction is strong, the ball rolls Figure 3. a: Average particle angular velocity (cid:104) (cid:36) Mi ( θ ) (cid:105) i aboutthe cluster’s center of mass, as a function of θ and N . In theloose region, (cid:36) i ( θ ) increases and decreases in a quantitativelysimilar manner for all values of N , with the only differencebeing the location on the dish at which this peak occurs.However, in the pinned region, the average value of | (cid:36) i ( θ ) |consistently increases with increasing N . Plotted in black arethe average locations at which particles enter the pinned andloose regions. Note that the apparent decrease of |(cid:104) (cid:36) Mi ( θ ) (cid:105) i | in the transition regions is an artifact of the ooblong shapeof the particle cluster. The transition regions occur at higherdistances from the center, but the linear velocity of the par-ticles during the transition regions does not change much. b:Average angular velocity of particles circling the dish, normal-ized with respect to the angular velocity of the wall about thecluster’s center of mass. This is done by dividing (cid:104) (cid:36) i ( θ ) (cid:105) i bythe angular velocity of the boundary around the cluster’s cen-ter of mass, a function of θ since the center of mass is not at C .At low N , a normalized velocity <1 in the pinned region in-dicates that those particles are traveling slower than the dishwall. As N is increased, the normalized velocity of particlesin the pinned region approaches 1, meaning the particles aremoving at the velocity of the wall. commensurately on the wall of the container without slip-ping, and γ = γ c = R cont R ball . Since R cont > R ball , the ballthus counterrotates in the lab frame. Conversely, whenthere is no friction, the ball cannot rotate at all due toperfect slip boundary conditions and γ = 0 < , corre-sponding to rotation. In between the two extremes, thedegree of wall-slip can be measured by a slip parameter s = γγ c = (cid:36) M ω R ball R drum , which varies from s = 0 for perfectslip conditions, to s = 1 for perfect no-slip conditions.While the single-particle model captures much of why aswirled particle cluster can either rotate or counterrotate,the details are more subtle since an ensemble of particlesrarely behaves exactly like a single rigid object. Even inthe M-frame, where the dynamics are in steady state, theinstantaneous angular velocity of any single particle (cid:36) i ( t ) depends on where it is in the dish, as shown in Fig 3a. Itis therefore more appropriate to refer to a particle’s in- Figure 4. a: The slip parameter s of the particles in thepinned region approaches 1 as N is increased, correspondingto the pinned region particles moving more coherently withthe moving wall. Additionally, the slip parameter s of rough-ened particles along the wall of the pinned region is closer to1 (non-slip conditions) than for smoother particles, meaningthat roughened particles also move more coherently with themoving wall in the pinned region. b: The local rotation of in-dividual particles along the wall in the pinned area decreaseswith increasing N , a result of increasing frustration betweenthe particles. Inset shows snapshot of a particle with over-laid Particle Image Velocimetry (PIV) vectors. We use PIVon the surface features of the particles to determine generalparticle rotation about any axis (blue) and spinning about anaxis going into the dish through the particle’s center (orange).c: Particles with roughened surfaces, and therefore increasedfriction, transition to counterrotation at a lower N than forsmoother particles. stantaneous velocity as (cid:36) i ( θ ) . The change in (cid:36) i ( θ ) with θ can largely be characterized by the loose and pinnedregions. In the loose region, the value of (cid:36) i ( θ ) increasesand decreases in a quantitatively similar manner for allvalues of N , with the only difference being the locationon the dish at which this peak occurs. However, in thepinned region, the average value of | (cid:36) i ( θ ) | consistentlyincreases with increasing N . Therefore, the pinned re-gion is the major contributor to the different behaviorsseen across varying N ’s. Comparison of (cid:36) i ( θ ) in thepinned region to the velocity of the moving boundary re-veals that the particles in the pinned region lag behindthe moving boundary for low values of N , but approachthe speed of the boundary at high N , as shown in Fig3b. Therefore, the average angular velocity of the clus-ter is dictated by the interaction of the particles in thepinned region with the container’s wall, in analogy to theslip-condition s of the single-particle model.Indeed it is possible to define an analogous slip param-eter for the swirling cluster, s = γγ c = (cid:36) M ω P ball P cont , wherethe radii of the rigid ball and container have been re-placed by their respective perimeters, P ball and P cont , toaccount for the deformability of the cluster. Accordingly,the slip parameter of the pinned particles in our systemincreases from approximately 0.5 at N = 25 to nearly 1at the highest values of N , as shown in Fig 4a. Similarly,our experimentally swirled cluster behaves increasinglysimilarly to a model particle with no-slip boundary con-ditions as N is increased, as shown in Fig 1b, confirmingthat counterrotation is due to increasingly no-slip bound-ary conditions.The effective slip parameter s is set by the dynam-ics of the individual particles along the wall within thepinned region, but what determines the dynamics in thepinned region? At low N , the particles at the wall arefree to spin and roll locally, falling behind the movingwall and effectively reducing the particle cluster’s fric-tion with the container. This effective slip of the clusteron the wall corresponds to s (cid:28) , and the cluster ro-tates in the lab frame. As N increases, the particles inthe pinned region pack together more tightly, and frictionbetween the particles dominates. When interparticle fric-tion is strong between two contacting particles, they mustspin about their own axes in opposite directions, dictat-ing antiferromagnetic-like interactions for spinning. Forsix-fold packing, as we often see in the pinned region,this results in geometric frustration, prohibiting any ofthe particles from rolling or spinning freely at high N ,as shown in Fig 4b. Similarly, particles in strong contactcannot advance as a tight single-file due to geometricfrustration, and can only roll on the ground side-by-side.The inability to individually roll causes the particles toeffectively stick to the container wall, resulting in lesseffective wall-slip ( s ≈ ). Therefore, the increased fric-tional effects at high N cause the dense particle cluster toroll, or rather treadmill, on the container wall, resultingin counterrotation.If friction is indeed what drives counterrotation, onewould expect that increased frictional effects achievedby alternative means would also promote counterrota-tion. Indeed, sandpaper-roughened particles transitionto counterrotation at N = 28 as compared to N = 36 forsmooth particles, as shown in Fig 4c, confirming the im-portance of friction to counterrotation. Furthermore, thecalculated slip parameter s for the roughened particles isconsistently closer to 1 than for the smooth particles, asshown in Fig 4a, affirming that counterrotation is asso-ciated with increasingly no-slip boundary conditions.Further experimentally testing the importance of fric-tion to counterrotation is a challenge since it is difficultto systematically fine tune or completely eliminate thefriction. We therefore turn to numerical simulations,which offer the unique advantage of adjusting physicalconstants that are impossible to change experimentally.We numerically simulate a system where N two-dimensional discs with radius r = 1 are swirled in a cir- cular container with radius R = 8 . , as shown in Fig 5a.The container is translated around a polygonal path with30 sides and amplitude A = 0 . , approximating a circlewhile allowing particle-wall collisions to be solved ana-lytically. The behavior of the particles is simulated usingan event-driven method, with the particles’ linear andangular velocities updated every collision and otherwisedetermined via Newton’s equations. Such a method al-lows us to exactly solve the dynamics up to floating-pointprecision, with the minimal number of parameters andingredients in the model. Collisions are perfectly elasticin the normal direction, and in the tangential directionare subject to frictional impulses derived from Coulomb’slaw, using a coefficient of friction µ d for particle-particlecollisions and µ w for particle-wall collisions [17] (see SIfor Methods).The average angular velocity (cid:36) for the simulationsis qualitatively similar to those measured in the exper-iments, exhibiting a rotation-counterrotation transitionwith increasing N , as shown in Fig 5b. The M-framedensity histograms and relative angular velocities arealso qualitatively similar to those of the experiment, asshown in Fig 5c and Fig 5d. The qualitative similar-ity between the experiments and simulations is striking,given that the simulations are significantly simplified –notably, there is no friction with the substrate, no three-dimensional rolling effects, and no normal damping dur-ing collisions, so the discs never actually stick to eachother. Therefore, these additional factors cannot be crit-ical for observing the transition.If the transition from rotation to counterrotation is in-deed driven by friction both between the discs and withthe container, then decreasing those frictions would re-quire more discs to achieve the same effective fricton,and therefore the counterrotation transition should oc-cur at a higher N . Indeed, decreased disc-disc and disc-wall friction cause the counterrotation transition to oc-cur at higher N , as shown in Fig 5b. Furthermore, com-pletely turning off friction should eliminate the transitionto counterrotation. We test these two cases separately.We first test the importance of boundary friction byeliminating disc-wall friction ( µ w = 0 ), resulting in africtionless boundary that serves only to contain the discsvia hard-core-like elastic interactions. The discs may stillfrustrate each others’ abilities to spin about their ownaxes, but without wall friction the discs are not encour-aged to rotate commensurately with the wall. In the labframe, the cluster of discs behaves as a solid-like unitsloshing around the container without any internal indi-vidual spinning (see SI videos). The angular velocity ofthe cluster therefore remains at that of the container re-gardless of N , as shown in Fig 5b. In the M-frame, thediscs form a rim at the boundary of the container withvery little internal movement (SI video). The steady-state of these particles arises from the centrifugal forcepushing particles to the wall, much like sand settling un- Figure 5. a: Snapshot of the simulated discs. b: Theparticle cluster transitions from rotation to counterrotationas N is increased when friction is present. This rotation-counterrotation transition point occurs at higher N whenfriction is decreased. When either disc-disc or disc-wall fric-tion is completely eliminated, the system never transitionsto counterrotation. Here "High friction" is µ d = µ w = 1 . ,"Medium friction" is µ d = µ w = 0 . , and "Low friction" is µ d = µ w = 0 . . When one friction is completely turned off,the other friction is set to 1.0. c: Density histogram of thediscs for a rotating ( N = 28 ) and counterrotating ( N = 48 )case. d: Quiver plots showing the average local disc veloc-ity deviation from the container for a rotating ( N = 28 ) andcounterrotating ( N = 48 ) case. der gravity. However, without additional forcing from thewall, the angular velocity of the cluster cannot change,so the system never transitions to counterrotation, veri-fying that disc-wall friction is critical to induce counter-rotation.Next, we separately test the importance of disc-discfriction by eliminating it ( µ d = 0 ) and bringing backdisc-wall friction. In this system, the outer discs are ac-celerated by frictional collisions with the wall and spinabout their own axes quickly. However, with no disc-disc friction, none of this spinning is transferred to theinner discs except through particle exchanges, so individ- ual discs spin with no coherence. The cluster of discs isloose and gas-like, unlike the rigid body appearance ofthe cluster when all friction is present (SI video). Withno disc-disc friction the average angular velocity remainsclose to the container velocity, except at very high N where it drops slightly. The system never transitions tocounterrotation, verifying that disc-disc friction is alsocritical to induce counterrotation.In conclusion, we perform experiments and numeri-cal simulations of swirled granular media and identifythe minimal ingredients for the particles and their in-teractions necessary to induce counterrotation. We usea frame of reference in which the dynamics of swirledgranular media are at steady state, the M-frame. In thisframe of reference we introduce a minimal model, wherethe granular cluster is replaced by a single ball within arotating drum. This system can capture the observed dy-namics by only tuning the amount of slip s that the ballexperiences with the wall of the rotating drum. When theball experiences strong slip on the wall ( s (cid:28) ), the sys-tem rotates in the lab frame, while minimal slip ( s ≈ )corresponds to counterrotation in the lab frame. In thegranular system, s (cid:28) is the result of individual par-ticles falling behind the moving wall while rolling freelyalong the wall and the ground, acting as bearings forthe particle cluster against the wall. This is generallytrue at low N when the effective friction is minimal, andcorresponds to rotation in the lab frame. On the otherhand, s ≈ occurs when the individual particles alongthe wall move with the wall due to high packing den-sities and geometric frustration. This is generally truefor densely packed particles at high N . Particle-particlefriction prevents closely-packed particles from spinningand rolling freely, converting the particle ensemble into asolid-like cluster. At the same time, particle-wall frictioncauses that solid-like cluster of particles to stick to thewall, causing the overall cluster to roll commensuratelyalong the wall and counterrotate, much like a pancake ina swirling pan.Our investigation has considered the dynamics of dis-crete particles but it would be interesting to model thesystem with continuum equations, coupling internal spin-ning to a continuum notion of vorticity. Such equationshave modeled related systems [18–20] and could perhapsgive insight into the sensitive interplay between pressure,vorticity, friction, and external forcing that leads to coun-terrotation. While our system is more complicated thanthe others that have been studied as it has no additionalsymmetries that lead to simplified equations, one couldstill study these equations numerically in the M-framewhere the external forces are stationary.In our swirling system, the large scale translation lo-cally drives individual particles to roll on the ground andspin on their neighbors and the boundary. The indi-vidual particles then interact with each other via fric-tional collisions, causing the motion of the particles toeventually coalesce into system-size rotation or counter-rotation. In addition to the swirling container and ro-tating drum systems, the behavior associated with in-teracting, individually-driven particles can give rise tophase changes and bifurcations in other physical systems.Several studies have shown that interactions between ac-tively moving or rotating objects, both self-propelled andsystem-scale driven, can lead to collective angular mo-mentum changes and rich phase behaviors [4, 5, 19–24].Specifically, in our system the interaction is particle fric-tion and geometric frustration. Particle-particle rollingfrustration is likely responsible for the rolling-jammingtransition, resulting in an apparent discontinuous jumpin effective friction when sheared layers of spherical mar-bles exceed a critical thickness [25]. The transition tosolid-like behavior in these systems is driven by the in-ability of contacting particles to co-spin or roll in-line,distinct from solid-like behavior that is achieved via jam-ming [26, 27] or rigidity percolation [28].Finally, our observation that the transition to counter-rotation can be manipulated by roughening the particlesinspires a speculative but interesting analogy with thetransition to turbulence in pipe flow, where a roughnessdependent transition has also been observed [29]. How-ever, unlike pipe turbulence, it is not clear if the tran-sition to counterrotation of swirling particles exhibits atrue critical phenomenon.We are grateful to Tadashi Tokieda for introducingus to this phenomenon. This work was supported bythe NSF (DMR-1420570). M. H.-C. and J. P. R. weresupported by US Department of Energy, Office of Sci-ence, Office of Advanced Scientific Computing Research,Applied Mathematics Program under Award No. DE-SC0012296. M. H.-C and J. P. R. thank Leif Ristroph forprocuring tabletop experimental materials. S. M. R. andM.H.-C acknowledge support from the Alfred P. SloanFoundation. ∗ To whom correspondence should be addressed. E-mail:[email protected], [email protected][1] M. G. Shats, H. Xia, and H. Punzmann, Phys. Rev. E , 046409 (2005).[2] H. Xia, D. Byrne, G. Falkovich, and M. Shats, NaturePhysics , 321 (2011).[3] A. Souslov, B. C. van Zuiden, B. Bartolo, and V. Vitelli,Nature Physics , 1091 (2017).[4] H. Wioland, E. Lushi, and R. E. Goldstein, Proceedingsof the National Academy of Sciences , 9733 (2014).[5] E. Lushi, H. Wioland, and R. E. Goldstein, New Journalof Physics , 9733 (2014).[6] A. Feltrup, K. Huang, C. Kruelle, and I. Rehberg, Eur.Phys. J. Special Topics , 19 (2009).[7] M. A. Scherer, K. Kötter, M. Markus, E. Goles, andI. Rehberg, Phys. Rev. E , 4069 (2000).[8] T. Yokoyama, K. Tamura, H. Usui, and G. Jimbo, in Proceedings of the Eigth European Symposium on Com-minution, Stockholm, 1994 (Elsevier, 1996) pp. 413–424.[9] M. A. Scherer, V. Buchholtz, T. Pöschel, and I. Rehberg,Phys. Rev. E , R4560 (1996).[10] D. Kumar, N. Nitsure, S. Bhattacharya, and S. Ghosh,Proceedings of the National Academy of Sciences ,11443 (2015).[11] M. A. Scherer, T. Mahr, A. Engel, and I. Rehberg, Phys.Rev. E , 6061 (1998).[12] O. Zik and J. Stavans, EPL (Europhysics Letters) ,255 (1991).[13] S. Fauve, S. Douady, C. Laroche, and O. Thual, PhysicaScripta , 250 (1989).[14] H. M. Jaeger, S. R. Nagel, and R. P. Behringer, Rev.Mod. Phys. , 1259 (1996).[15] D. Bi, J. Zhang, B. Chakraborty, and R. P. Behringer,Nature , 355 (2011).[16] T. S. Majmudar, M. Sperl, S. Luding, and R. P.Behringer, Phys. Rev. Lett. , 058001 (2007).[17] Y. Wang and M. T. Mason, J. Appl. Mech. , 635(1992).[18] J. S. Dahler and L. E. Scriven, Nature , 36 (1961).[19] J. C. Tsai, F. Ye, J. Rodriguez, J. P. Gollub, and T. C.Lubensky, Phys. Rev. Lett. , 241 (2005).[20] B. C. van Zuiden, J. Paulose, W. T. M. Irvine, D. Bar-tolo, and V. Vitelli, Proc. Natl. Acad. Sci. , 12919(2016).[21] N. H. P. Nguyen, D. Klotsa, M. Engel, and S. C. Glotzer,Phys. Rev. Lett. , 075701 (2014).[22] K. Yeo, E. Lushi, and P. M. Vlahovska, Phys. Rev. Lett. , 188301 (2015).[23] M. Workamp, G. Ramirez, K. E Daniels, and J. Dijks-man, Soft Matter (2018).[24] F. Moisy, J. Bouvard, and W. Herreman, EPL (Euro-physics Letters) , 34002 (2018).[25] C. Marone, B. M. Carpenter, and P. Schiffer, Phys. Rev.Lett. , 248001 (2008).[26] I. Peters, S. Majumdar, and H. Jaeger, Nature , 214(2016).[27] P. Coussot, N. Roussel, S. Jarny, and H. Chanson,Physics of Fluids , 011704 (2005).[28] S. Alexander, Physics Reports , 65 (1998).[29] N. Goldenfeld, Phys. Rev. Lett. , 044503 (2006).[30] D. Frenkel and B. Smit, Understanding Molecular Sim-ulation , From Algorithms to Applications (AcademicPress, 2001).[31] Y. Wang and M. T. Mason, J. Appl. Mech. , 635(1992).[32] J. B. Keller, J. Appl. Mech. , 1 (1986).[33] L. D. Landau and E. M. Lifshitz, Mechanics (Butterworth-Heinemann, 1976).
Supplementary Information for “Geometric frustration induces the transition betweenrotation and counterrotation in swirled granular media”
VIDEO
Included:Experimental lab frame: rotating, stall, counterrotatingExperimental M-frame: rotating, stall, counterrotatingNumerical lab frame all friction: rotating, counterrotatingNumerical M-frame all friction: rotating, counterrotatingNumerical lab frame no disc friction: rotating, counterrotatingNumerical M-frame no disc friction: rotating, counterrotatingNumerical lab frame no wall friction: rotating, counterrotatingNumerical M-frame no wall friction: rotating, counterrotatingMinimal model M-frame: no slip and slipMinimal model lab frame: no slip and slip
CALCULATION OF LOCATION THAT PARTICLES BECOME PINNED AND FREE (FIG 3)Becoming pinned
Particles experience multiple collisions while transitioning from the loose to the pinned region. The angular locationon the dish at which particles become pinned is therefore determined by the average angle at which these collisionsoccur. A collision event is assigned to a particle p at time t if the particle experiences a sufficiently sharp turn in itstrajectory, calculated using a threshold on the second derivative of its coordinates x ( t, p ) and y ( t, p ) . Becoming free
The angular location on the dish at which particles become free is determined by the angular position at whichthe particle no longer travels parallel to the dish edge. To determine this, we calculate the best fit slope for theparticle’s coordinates between the preceding 7 frames and the following 7 frames. When this slope deviates from theperpendicular to the particle’s radial vector on the dish (or the particle’s azimuthal vector) by more than 10 degrees,we mark this particle as becoming “free.” Mathematically: θ free = allframes-7 (cid:88) t =7 allparticles (cid:88) p =1 θ ( t, p ) ∗ boolean ( | θ ( t, p ) − arctan( − x ( t, p ) − x dish center y ( t, p ) − y dish center ) | > , where θ ( t, p ) = (cid:80) t +7 i = t − ( x ( i, p ) − ¯ X )( y ( i, p ) − ¯ Y ) (cid:80) t +7 i = t − ( x ( i, p ) − ¯ X ) , and ¯ X = (cid:80) t +7 i = t − x ( i, p )15 , ¯ Y = (cid:80) t +7 i = t − y ( i, p )15 . CALCULATION OF PARTICLE FRUSTRATION (FIG 4B)
Particle Image Velocimetry (using MATLAB’s PIVlab) is run on the particles within the pinned region. The PIVinterrogation windows are set such that a total of 16 velocity vectors is calculated within each particle. The velocityvalues ( (cid:126)v
PIV ) are then used to calculated the following:
Average velocity
Average velocity is the average of the velocity vector magnitudes:Average velocity = all vectors (cid:88) | (cid:126)v PIV | . Average particle spin
Average particle spin takes into account only the azimuthal velocity, calculated by averaging the dot product ofeach PIV vector with its azimuthal vector (the azimuthal vector of the coordinates at which the PIV arrows havebeen calculated, taking the middle of the particle as the origin):Average particle spin = all vectors (cid:88) i ( (cid:126)v PIV ,i · azimuthal vector at point i | azimuthal vector at point i | ) . NUMERICAL SIMULATION METHODOverview
Our simulations consider N identical discs of radius R disc in a circular container of radius R container . We treat thecontainer as a disc with a different radius so the number of discs in the simulation is N + 1 . The discs are labelled , , , . . . , N , where 0 indexes the container. The discs have mass m and moment of inertia I . The container isassumed to have infinite mass and moment of inertia. The state of the system at any point in time is given by thepositions of the discs and container { x i } Ni =0 , x i ∈ R , the velocities of the discs and container { v i } Ni =0 , v i ∈ R , andthe angular velocities of the discs { ω i } Ni =1 (the container does not rotate.)The container is “swirled,” translated along a path which is topologically a loop, at frequency ω . We translateat constant speed along a path which is a regular polygon with M sides and amplitude A (distance from center ofpolygon to vertex.) This path approximates a circle but allows us to solve analytically for disc-container collisions.The dynamics of the discs are assumed to satisfy Newton’s equations everywhere except at collisions, so the discstravel with constant linear and angular velocity. Notably, they do not experience friction with the floor. When discscollide with each other or the container, their linear and angular velocities are updated according to a “collision law,”a rule for updating velocities which is consistent with certain physical principles. We choose our disc-disc collisions tobe perfectly elastic in the normal direction, and in the tangential direction apply an impulse consistent with kineticCoulomb friction (we ignore static friction.) Our disc-container collisions are similar but we consider the container tohave infinite mass and moment of inertia, so its velocity doesn’t change after a collision. More details on the collisionlaw are given in section below.The dynamics are implemented using an event-driven method. We consider three kinds of events:1. Two discs collide;2. A disc collides with the container;3. The container reaches a vertex of its polygonal path.The method proceeds by repeating the following steps.(i) Calculate the times of all future events;(ii) Determine the time of the event that happens the soonest;(iii) Determine the positions and velocities of the system at that time and then update them using the rule associatedwith the event.When none of the events occurs we may solve analytically for the positions and velocities of all the discs. When anevent occurs we must change the velocities using the rule associated with that event. Finding the next event
In this section we describe how the times of future events are calculated, considering each type of event in turn.1. We calculate the collision time t ij for each pair of discs ( i, j ) . For some disc i with position x i at the currenttime, the position at any time increment t from the current event (if there have been no collisions) is p i ( t ) = x i + tv i . (S1)A pair of discs ( i, j ) in isolation will collide when the distance between their centers equals the sum of theirradii: | p i ( t ) − p j ( t ) | = (2 R disc ) . (S2)This is a quadratic equation in t with positive leading coefficient, whose roots give the collision times. If bothroots are negative, then the discs will never collide on their present trajectories (though they may have collidedin the past). If the roots have different signs, then the current distance between centers is less than R disc ,and so the discs are overlapping; note that this should never happen if the simulation has been implementedcorrectly. If the roots are both positive, then the lesser of them is the time that the discs will collide. If theroots are not real, then the discs will never collide. Therefore we can set t ij to be the lesser of the roots if it isreal and positive, or t ij = ∞ otherwise.We make one small change to the above to account for numerical errors. It is sometimes the case that two discsthat have just collided are predicted to collide next by the program (i.e., their solutions to (S2) are real andpositive.) In theory, the parabola in t from (S2) should only have one root at the origin, but numerical errorsmay cause the roots to both be real and positive. When this has happened, the predicted next time of collisionhas always been less than − . To prevent these errors from being mistaken for actual collision predictions,we require that a predicted time of collision be greater than − , otherwise we ignore that prediction. Theminimal time between collisions is around − and − with and without friction, respectively.2. For collisions between discs and the container, the collision times are calculated in a similar way as they are forpairs of discs except now we use the quadratic equation | p i ( t ) − p ( t ) | = | R container − R disc | (S3)and choose the greater of the roots. As long as the discs always stay within the container, this quadratic willalways have real roots of different signs (except in the case that the container and the disc have exactly thesame linear velocity.)3. Calculating the time of the next update to the container’s velocity is trivial, since the container reaches a vertexof its polygonal path at fixed time increments.Our current method calculates the collision times of all pairs of discs, which is an O ( n ) calculation. The methodcould be accelerated by using cell lists, making it an O ( n ) calculation [30]. Our method was fast enough for thenumbers of discs under consideration but simulating large numbers of discs would require cell lists. Collision law
We now describe how to update the velocities of the discs after a collision. Our approach for relating tangentialimpulses to normal impulses is based on the calculations in [31], although our equations are simpler since we consideronly discs, and [31] considers arbitrary rigid bodies. These calculations can be justified asymptotically if the timeover which the collision takes place is much shorter than a typical evolution time [32].Suppose discs i and j have just collided so their boundaries are exactly tangent, as in Figure S1. Let ˆ n = ( x j − x i ) / | x j − x i | be the unit vector normal to the point of contact, pointing from x i to x j , and let ˆ r = − ( x j − x i ) ⊥ / | x j − x i | be the unit tangent vector at the contact. Here ( x, y ) ⊥ = ( − y, x ) . For generality we let the discs have masses m i , m j and moments of inertia I i , I j .Let v i , v j , ω i , ω j be the velocities of the discs exactly when they contact but before the collision law has beenimplemented, and let v (cid:48) i , v (cid:48) j , ω (cid:48) i , ω (cid:48) j be the velocities after the collision law. Let v normi ( j ) , v tangi ( j ) be the normal and disc j disc i ω j ω i ˆ n ˆ r Figure S1. Sketch of two discs at the point of contact. Vectors ˆ n , ˆ r are the normal and tangent directions as described in thetext. tangential components of each disc’s initial linear velocity, i.e. v normi = v i · ˆ n , v tangi = v i · ˆ r , and similarly for v j . Thevelocity components are updated by increments ∆ v normi ( j ) , ∆ v tangi ( j ) , ∆ ω i giving updated velocities v (cid:48) i = ( v normi + ∆ v normi )ˆ n + ( v tangi + ∆ v tangi )ˆ r , ω (cid:48) i = ω i + ∆ ω i , and similarly for v (cid:48) j .The change in normal velocities is calculated assuming an elastic collision. Because we consider only discs, thenormal force always points toward the centers of mass and does not mix with the tangential direction, so an elasticcollision is equivalent to conserving energy and momentum conservation in the normal direction. (This is related tothe equivalence of Newton’s law of restitution, and the Poisson hypothesis, for discs, see [31].) The increments are ∆ v normi = 2 m j m i + m j ( v normj − v normi )∆ v normj = 2 m i m i + m j ( v normi − v normj ) = − m i m j ∆ v normi . We wish to derive tangential increments from the Coulomb friction law for kinetic friction. Suppose the collisionis not instantaneous, but rather takes place over a time interval of length τ . Let F n ( t ) , F f ( t ) be the normal andtangential components of the force acting on disc j during the collision; F f ( t ) is the friction force. The forces actingon disc i are equal and opposite, by Newton’s third law. The Coulomb friction law says that [33] F f ( t ) = s µ | F n ( t ) | , (S4)where µ is the coefficient of kinetic friction. The above holds provided the relative velocity at the point of contactis nonzero; after it hits zero F f ( t ) = 0 and the relative velocity is zero for the remainder of the collision. The sign s ∈ {− , } is chosen to ensure that friction reduces magnitude of the relative velocity of the points in contact. It iscomputed from the initial relative velocity v rel (the velocity of the contact point on disc j minus the velocity of thecontact point on disc i ) as s = − sgn( v rel ) , v rel = v tangj − v tangi + R disc ( ω i + ω j ) . (S5)If we knew the normal force F n ( t ) and the length τ of the collision, then we could update the tangential velocitiesusing Newton’s second law by integrating (S4). We don’t know either F n ( t ) or τ (and they probably depend sensitivelyon the details of the collision), but, crucially, we know the total impulse in the normal direction, and this determinesthe total tangential impulse via (S4). Indeed, following [31], we integrate (S4) over the interval [0 , τ ] to obtain J f = s µJ n , (S6)where J f = (cid:82) τ F f ( t ) dt , J n = (cid:82) τ F n ( t ) dt are the impulses in the tangential and normal directions (for disc-disccollisions, by the choice of axes we have F n ( t ) ≥ .) Equation (S6) is valid provided the relative velocities neverchange sign; if they do the tangential impulse must be truncated as described below.The impulse in the normal direction ˆ n for disc j is the change in the normal component of momentum, J n = m j ∆ v normj . (S7)The magnitude of the maximum impulse in the tangential direction is therefore J maxf = µm j ∆ v normj . (S8)The actual tangential impulse is J f = s min( J maxf , J stickf ) , (S9)where J stickf is the magnitude of the impulse required to make the relative tangential velocity between discs exactlyzero (the superscript indicates the discs “stick” to each other in the tangential direction.) A formula for J stickf will bederived in (S11) below.Given J f , the tangential increments are computed from the change in tangential and angular momenta: m j ∆ v tangj = J f , m i ∆ v tangi = − J f , I j ∆ ω j = R disc J f , I i ∆ ω i = R disc J f The tangential and angular velocities are updated as v (cid:48) ,tangi = v tangi − J f m i , ω (cid:48) i = ω i + R disc J f I i v (cid:48) ,tangj = v tangj + J f m j , ω (cid:48) j = ω j + R disc J f I j (S10)The truncated impulse magnitude J stickf is the impulse required to make the relative tangential velocity betweenpoints in contact equal zero. It is found by solving equation v tangj + J stickf m j − (cid:16) v tangi − J stickf m i (cid:17) + R disc (cid:16) ω j + R disc J stickf I j + ω i + R disc J stickf I i (cid:17) = 0 , resulting in J stickf = (cid:12)(cid:12)(cid:12) v tangj − v tangi + R disc ( ω i + ω j ) (cid:12)(cid:12)(cid:12) m i + m j + R disc I i + R disc I j . (S11)Processing collisions with the container is similar, except that we set m = I = ∞ . In this case, the initialincrement is ∆ v normi = 2( v normj − v normi )∆ v normj = 0 Since ω = 0 we do not need to worry about the change in sign in some of the formulas because the collision is interior.We make one final remark, that one can modify this method to allow for non-elastic behavior in the normaldirection. However, it is then possible for discs to stick to each other, after which one has to track the dynamics ofthe cluster. Collision laws between clusters can become arbitrarily complicated as one solves for the propagation ofimpulse between the discs in each cluster. Since the dynamics we seek to understand do not seem to require nonelasticcollisions we did not consider adding this additional complexity. SIMULATION PARAMETERS
The container’s path is a polygon with M = 30 sides, and the period of the swirl is 11.866 units of time. Thecontainer moves at one unit of distance per unit of time, so the radius of the polygon formed by the container’s pathis 1.892. Since the container’s radius is 8.6, this makes the ratio of the radius of the container’s path to the radius ofthe container equal to 0.22, which matches the experimental setup.The following are the default parameters used in the plots and videos.Number of discs N variable; typically 28-48Mass of disc m I R disc R container M ω − Disc/Disc Friction Parameter µ d µ w N between and , ofthe total angular velocity of the discs about their center of mass over iterations.For Figure 5c, the M-frame was divided evenly into a 100x100 grid of cells, and the amount of times that a disc’scenter was found inside a cell was measured for N = 28 and . Cells that are colored red were found to contain disccenters very often.For Figure 5d, the M-frame was divided evenly into a 20x20 grid of cells, and the average linear velocity of discsfound in each of the cells was calculated for N = 28 and48