Persistent structures in a three-dimensional dynamical system with flowing and non-flowing regions

Mixing of fluids and mixing of solids are both relatively mature fields. In contrast, mixing in systems where flowing and non-flowing regions coexist remains largely unexplored and little understood. Here we report remarkably persistent mixing and non-mixing regions in a three-dimensional dynamical system where randomness is expected. A spherical shell half-filled with dry non-cohesive particles and periodically rotated about two horizontal axes generates complex structures that vary non-trivially with the rotation angles. They result from the interplay between fluid-like mixing by stretching-and-folding, and solids mixing by cutting-and-shuffling. In the experiments, larger non-mixing regions predicted by a cutting-and-shuffling model alone can persist for a range of protocols despite the presence of stretching-and-folding flows and particle-collision-driven diffusion. By uncovering the synergy of simultaneous fluid and solid mixing, we point the way to a more fundamental understanding of advection driven mixing in materials with coexisting flowing and non-flowing regions.

T he goal of mixing is to rearrange initially segregated matter into states where the constituent elements are homogeneously distributed. In fluids, where the elements are atoms or molecules, mixing at low Reynolds numbers can be achieved by the stretching-and-folding of chaotic flows combined with thermal diffusion which drives mixing at the smallest length scales. In bulk solids composed of macroscopic (athermal) particles, the particles can be deliberately rearranged, as in the cutting-and-shuffling of a deck of cards. Other rheological materials-Bingham fluids for example-fall between these two extremes. While mixing of fluids and mixing of solids have long histories and are relatively mature fields, little is understood about mixing in three-dimensions (3D) when solid and flowing regions coexist. For example, in yield stress materials constituent elements move together as a solid where local stresses are low, but flow in relative motion where the yield stress is exceeded. Common examples of yield stress materials include paint, wet concrete, polymer mixtures,(1) and granular materials, e.g., sand. Understanding 3D mixing in such materials is critical in many domains where variations in local concentration can be disastrous, including the pharmaceutical industry,(2) composite materials, (3) and concrete manufacturing. Also noteworthy is the lottery with the largest payout in the world, the Spanish Christmas Lottery, where 100,000 wooden balls are mixed by rotating a 2 m diameter sphere about a horizontal axis with the expectation of randomness. Our results show that this expectation cannot be taken for granted.
To study the interaction of mixing by stretching-and-folding with mixing by cutting-and-shuffling, we consider a geometrically simple 3D system with localized flow -a spherical tumbler half-filled with a dry granular material and rotated alternately about orthogonal horizontal axes by angles (θz, θx) beyond the repose angle of the granular material, β. In tumblers, particles flow in a relatively thin layer at the free surface, while below the surface, particles move together in solid body rotation about the rotation axis. Using X-ray imaging, we mapped the location of a 4 mm diameter spherical tracer particle in a D(= 2R) = 14 cm diameter spherical tumbler half filled with 2 mm diameter glass spheres after each iteration for numerous protocols [i.e., (θz, θx) pairs]. The larger diameter tracer particle flows on the free surface and is subsequently deposited in the bed near the tumbler wall.
Results and Discussion. Consider first the itinerary of a tracer particle during 500 iterations of the (57 o , 57 o ) protocol, see Fig. 1(a). After each iteration the tracer particle was alternately displaced between three distinct regions as indicated by the gray lines. Figure 1(b) shows a second 500 iteration experiment in which the tracer particle moved periodically between the same three regions (blue points) and between a second set of period-3 regions (red points) as well as aperiodically (black points). Similar period-3 non-mixing regions are observed in the "nearby" protocols (54 o , 54 o ) and (60 o , 60 o ).
To better understand the existence of these non-mixing and mixing regions, we consider the predictions of a standard advection based continuum model (4) for tumbler flow, in which, for rotation about any axis (the z-axis here with x in the streamwise direction and y normal to the free surface), the non-dimensionalized velocity field u = (u, v, w) is piecewise defined such that the flowing layer (0 ≥ y ≥ −δ) velocity is u fl = ((δ + y)/ 2 , xy/δ, 0) and the bulk (y < −δ) solid body rotation velocity is u b = (y, −x, 0). The lenticular flowing layer interface with the bulk is located at ω/γ is the maximal dimensionless flowing layer thickness at the center of the sphere (x = z = 0) for shear rateγ and angular rotation velocity ω. All variables are dimensionless-lengths (x, y, z, r, δ) are normalized by the tumbler radius R and rotation period T is normalized by 1/ω. This flowing layer (FL) model, which includes stretching characteristic of chaotic flows (5), is parametrized by , which is set to 0.15 based on the particles used in the experiments (see Methods section). To characterize mixing in the model, blue passive tracer points are seeded at the intersection of the flowing layer boundary δ with a hemispherical subshell having normalized radius r = 0.9 * before the first rotation, while red points are seeded in the same location after a half iteration (i.e., rotation about the z-axis only). In Fig. 1(c), the Poincaré (stroboscopic) map of tracer points advected by the model displays uniform mixing throughout the domain except for six empty regions (A1-A3 and B1-B3). These non-mixing regions correspond to the regions in experiments where the tracer particle lingers, as indicated by their mapping onto Figs. 1(a,b).
Although the FL model captures the main features of the experiments shown in Figs. 1(a,b), deeper insight into why persistent structures form is gained by focusing exclusively on solids mixing. This is accomplished by taking the δ → 0 limit (an infinitely thin flowing layer) in which tracers instantaneously jump from the point they reach the free surface to a downstream point symmetric about the midpoint of the free surface and, consequently, undergo only solid body rotation. Since there is no flowing layer and, consequently, no shear, the biaxial mixing protocol in the δ → 0 limit corresponds to a radially invariant hemispherical domain with mixing dynamics that are equivalent to slicing the hemisphere into pieces that are rearranged and then reassembled into a hemisphere again. This type of transformation is called a piecewise isometry (PWI). (6) Similar to the FL model, the boundaries formed by the slicing after a rotation about a single axis are used as * Radius r = 0.9 was used to mimic tracer particle position in experiments instead of a larger radius (0.95 ≤ r ≤ 0.97) as elliptic regions at larger r > 0.9 reveal more intricate structures that are not apparent in experiments likely due to particle size effects and collisional diffusion. As the intricate elliptic regions at r > 0.9 occur at the same location on the hemispherical shell as elliptic regions at r = 0.9, we use the elliptic orbits at r = 0.9 for comparison to experiment for visual clarity. initial conditions for the tracer points whose trajectories form a subset of the exceptional set,(7) which is the skeleton for the transport dynamics in PWI systems. Here, blue tracers are seeded where the domain is cut by the action of the θz rotation, while red tracers are seeded where the domain is cut by the θx rotation, similar to the initial conditions used in the FL model. The mixing mechanism for a piecewise isometry is simply cutting-and-shuffling in analogy with mixing a deck of cards. (8) Similar to the FL model, iteration of the PWI model for the (57 o , 57 o ) protocol ( Fig. 1(d)) generates open regions avoided by the tracers, known as cells. These cells vary in size and periodicity across the domain, with certain areas dominated by a particular color of tracer particle. The largest circular cells have the lowest periodicity. The similarity in size and location between the largest cells in the PWI model, the elliptical domains in the FL model ( Fig. 1(c)), and the non-mixing regions in the experiments ( Figs. 1(a,b)) is clear. This agreement suggests that the periodic regions observed in experiment for the (57 o , 57 o ) protocol result from the solids mixing structure generated by cutting-and-shuffling.
The periodicity of non-mixing structures depends on the protocol, for example consider the results for the (90 o , 90 o ) protocol shown in Fig. 2. In experiment ( Fig. 2(a)), a tracer particle initially seeded in the A1 region, cycles between the A1 and A2 regions and between the B1 and B2 regions for a total of 263 and 46 iterations, respectively. Compared to the period-3 non-mixing regions under the (57 o , 57 o ) protocol, the tracer escapes more frequently from the period-2 non-mixing regions. The non-mixing regions A1-A2 and B1-B2 correspond to the white regions avoided by the tracer points in the FL model ( Fig. 2(b)), while in the PWI model ( Fig. 2(c)) non-mixing regions occupy the entire domain since the entire hemisphere returns to its initial condition every two iterations. As for the (57 o , 57 o ) protocol, the finite-depth flowing layer converts a portion of the PWI model's non-mixing regions into mixing regions.
For the half-full spherical tumbler mixed by alternating rotations about orthogonal axis used here, persistent non-mixing period-n regions always appear in pairs, a manifestation of the relationship between the full and half iteration structures. To demonstrate, the sets of non-mixing regions shown in Fig. 2(b) swap positions upon an additional half-cycle of the mixing protocol ( Fig. 2(d)). Thus, the B-set of non-mixing regions is the half-period offset of the A-set of non-mixing regions.
Both PWI and FL models capture the large scale persistent mixing and non-mixing structures observed in experiments for the protocols we have examined so far, but they do not describe the transitions of the tracer particle between islands seen in experiments, e.g., Fig. 1(b) and Fig. 2(a). This is because neither model includes diffusion, which in granular flows is driven by particle-particle and particle-wall collisions. The collisional diffusion coefficient in flowing macroscopic particles scales as d 2γ where d is the diameter of the particles in the tumbler andγ is the shear rate(9); in our tumbler, diffusion displaces a particle on average by about its diameter per flowing layer pass.(10) For the (57 o , 57 o ) protocol which produces ∼ 1 flowing layer pass per iteration, the root-meansquare displacement after 500-iterations (d √ 500) is 4.5 cm (approximately the tumbler radius). Based on the size of the non-mixing regions under the (57 o , 57 o ) protocol, we expect a mean residency time of about 350 iterations in the regions, which is on the order of the observed residence times of at least 500 iterations in Fig. 1(a), 326 iterations for the period-3 A regions in Fig. 1(d), and 101 iterations for the period-3 B regions in in Fig. 1(d). Consequently, particles in experiments "leak" out of non-mixing regions with a residence time that decreases as the square root of the size of the non-mixing regions.
Other protocols can produce completely mixed domains. For example under the rotationally asymmetric (75 o , 60 o ) protocol (Fig. 3), the tracer particle explores most of the domain in experiments ( Fig. 3(a)), while the FL ( Fig. 3(b)) and PWI (Fig. 3(c)) models generate mixing regions that completely and nearly completely, respectively, fill the domain. Mixing is not guaranteed by an asymmetric protocol, i.e., θz = θx (see Extended Data Fig. 1(a)), and, conversely, a symmetric protocol, i.e., θz = θx, does not guarantee the existence of non-mixing regions.
Under the three protocols shown in Figs. 1-3, when nonmixing regions appear in experiments, they are present in the FL model and correspond with the largest cells in the PWI model. In contrast, under the (45 o , 45 o ) protocol (Fig. 4), a mixing barrier emerges that is not captured by the PWI model. A tracer particle in two experiments ( Fig. 4(a-b)) follows two different extended finger-like period-3 structures that are separated by a leaky barrier to particle transport and together cover the entire domain. These interdigitated period-3 structures are readily apparent in the FL model (Fig. 4(c)). Red and blue tracers each dominate half of the domain with an elongated non-mixing region in each "finger," which roughly corresponds to the tracer particle positions from experiment. The region dominated by the blue tracers appears to have only two "fingers", because the third "finger" is mostly contained in the flowing layer, which is not visible in the view shown in Fig. 4(c). The blue edges of this "finger" are evident at the periphery of the domain as well as the periphery in the experiments ( Fig. 4(b)). For the equivalent half iteration structures in the FL model, the red tracer dominated regions alternate with the regions dominated by blue tracers with one of the red "fingers" mapped to the flowing layer. To delineate the two regions, we followed a tracer point in the FL model located between the red and blue tracer dominated regions to produce the grey points in Fig. 4(a-b). The path of this tracer suggests a mixing barrier that wraps around the entire domain. Unlike the PWI model predictions for protocols shown in Figs. 1-3, the predicted structure under the (45 o , 45 o ) protocol ( Fig. 4(d)) is less clearly related to the experiment and the FL model as it lacks the typical large cells that manifest as non-mixing regions in the experiment. Instead the PWI model predicts large arrowhead-like features consisting of multiple small cells with transport barriers between red and blue arrowheads. The colors and positions of the arrowhead features correspond to the "finger" features in the experiment and FL model.
The structures predicted by the FL model are in close agreement with experimental observations regardless of the protocol, while the structures predicted by the PWI model are not always observed. For example, the smaller non-mixing cells in the PWI model for the four protocols we examine above are not seen in experiments. To better understand the relationship between the non-mixing structures predicted by the PWI and FL models, particularly for the (45 o , 45 o ) protocol, we consider the influence of flowing layer thickness. Physically, flowing layer thickness increases with rotation rate ω and decreases with shear rateγ, since = ω/γ = δ(0, 0)/R. (11) Additionally, the flowing layer is typically thicker than ∼ 10 particle diameters in experiments. As their boundaries are set by the boundary of the flowing layer during the interchange of rotation axes, see Fig. 6(a-b). As a result, the periodicity of the non-mixing region is determined by the rotation angle, and periodicity increases as rotation angle decreases because the number of flowing layer passes for a non-mixing region between actions decreases. In contrast, the non-mixing regions that appear only for > 0 periodically land entirely within the flowing layer at the end of a rotation. Their boundaries, unlike persistent non-mixing regions, are not necessarily coincident with the flowing layer boundary at the end of a rotation. In the example shown in Fig. 6(c-d), the non-mixing regions in the bulk are long and thin and their boundaries contact the flowing layer boundary when two of the trio are in the bulk (A2 and A3 at 500 iterations, and B2 and B3 at 500-1/2 iterations) while the other region lands in the flowing layer. Additionally, the boundary of the mixing barrier between red and blue "fingers" maps approximately to the boundary of the flowing layer, but is less clearly defined than the non-mixing region boundaries. Based on this study of granular material in a spherical tumbler under a bi-axial mixing protocol, we expect other 3D dynamical systems with coexisting solid and fluid regions to exhibit similar complex mixing structures that result from the distinct and competing mixing processes associated with fluid stretching-and-folding and solid cutting-and-shuffling. A deeper understanding of these hybrid-mixing systems based on the structures determined by piecewise isometry theory is expected when the spatial extent of flowing regions is limited or when the flows are weakly shearing.

Materials and Methods
Experiments. Experiments were conducted using d = 1.89 ± 0.09mm diameter soda-lime glass beads (SiLiglit Deco Beads, Sigmund Lindner GmbH, Germany) in a half-filled D = 14-cm diameter acrylic spherical tumbler. A density matched (ρ = 2.5 g cm −1 ) dtracer = 4-mm diameter X-ray opaque tracer particle constructed from two 3D-printed plastic hemispherical shells with a Pb-Sn solder sphere in the center was used to visualize the flow under X-ray imaging. The larger diameter of the tracer particle caused the particle to flow at the surface of the flowing layer and near the tumbler wall in the solid bed (bold path in Extended Data Fig. 2). The tumbler was mounted on an apparatus that performs the biaxial protocol while being X-ray imaged. For each action (U or W, in Extended Data Fig. 2), the sphere was rotated about a single axis for protocol angle (θz or θx) by three wheels driven by a motor mounted on the turntable at rotation speed ω = 2.6 rpm. Since particles only flow continuously once the free surface reaches the dynamic angle of repose β with respect to the horizontal, † the tumbler was slowly rotated to the angle of repose prior to each action. After each action, the tumbler was rotated in the reverse direction to ensure the free surface was horizontal before switching the rotation axis. The apparatus reoriented the drive wheels for the next rotation by raising the tumbler off the wheels, rotating the turntable to the next tumbling axis, and then lowering the tumbler back down onto the wheels.
To minimize accumulated error from executing the protocol over large numbers of iterations, a position sensor was used to ensure that the turntable returned to the original axis after each iteration. Additionally, a thin X-ray opaque fiducial marker (lead tape) was mounted on the interior wall of the spherical tumbler to track any tumbling deviations. With these measures in place, the system has an angular displacement error of less than 1 o per protocol iteration. The errors are not systematic and, hence, average to zero.
Components of the spherical tumbler apparatus within the Xray beam path consist of X-ray transparent materials (aluminum and plastic) to ensure a clear image suitable for particle tracking. Lead tape was also affixed to the turntable to reduce variations in X-ray beam attenuation due to variation in path length through the particle bed.
A high speed grayscale camera (Point Grey BlackFly 12A2M) mounted directly beneath the output image port of the X-ray image intensifier (see schematic in Extended Data Fig. 3) provides a subsurface bottom view of the tumbler. Images were acquired prior to each tumbling action: the first frame, when the free surface was horizontal and subsequent frames while the tumbler was slowly rotating to the angle of repose. The direction and magnitude of tracer particle displacement in these frames identified whether the particle was at the free surface or in the bulk near the tumbler wall. The image distortion due to image intensifier and camera optics was corrected using the MATLAB Image Processing Toolbox function imwarp. The geometric transformation required by the function imwarp was a 4th order polynomial model generated from applying the function fitgeotrans to a Cartesian hole pattern (diameter 3.5 mm with 5.1 mm spacing) image. The corrected hole pattern image was used to generate the matrix transformation from the hole pattern's pixel spacing to the physical grid dimensions. To track the tracer particle and the fiducial indicator, each image was divided by a background image generated from the average of all iterations from a particular run to reveal the two features of interest. The tracer particle and the fiducial indicator were distinguishable from each other by their eccentricity and were tracked automatically using 2D feature finding MATLAB algorithms developed by the Kilfoil group ‡ using methods from Crocker et al.(12) † Since the particle bed was stationary at the onset of flow, the free surface was at the static angle of repose βs. There was a small avalanche at the start of flow where the free surface first relaxed to the dynamic angle of repose β < βs ‡ The repository for the code is hosted at http://people.umass.edu/kilfoil/tools.php and they remain there until they reach the flowing layer again. The angle of the flowing free surface is β with respect to gravity, g. Three mean particle paths are shown for a two-dimensional slice for each action, with the bold path (near the free surface and closest to tumbler wall) being a representative path for the large X-ray opaque tracer particle.
Flowing layer (FL) model. The flowing layer (FL) model (4,13) assumes that the flow is primarily two-dimensional in the streamwise direction and confined to a thin lenticular flowing layer with a constant depthwise shear rateγ for each tumbling action U and W (Extended Data Fig. 2). Using a Cartesian coordinate system with its origin at the center of radius R spherical tumbler, rotation is clockwise about the z-axis and x-axis at a rotation speed ω for an angular displacement θz and θx for U and W actions, respectively. For convenience, all variables are dimensionless-length scales (x, y, z, r, and flowing layer thickness δ) are normalized by R and rotation period T is normalized by 1/ω. The interface between the flowing layer and the bulk is given by δ(x, z) = √ 1 − x 2 − z 2 , where = ω/γ is the maximal dimensionless flowing layer thickness (at the center x = z = 0). The velocity u = (u, v, w) is piecewise defined for the flowing layer and the bulk. For the U action, the flowing layer (y ≥ −δ) velocity is u fl = ((δ + y)/ 2 , xy/δ, 0) and the bulk (y < −δ) is in solid body rotation with velocity profile u b = (y, −x, 0). Analogously, the velocity field for the W action is obtained by interchanging x-and z-components. Particle positions are generated by integrating the velocity field using the Runge-Kutta (RK4) method in the flowing layer and the semi-implicit Euler method in the bulk. (4,13) The dimensionless time step is ∆t/T = 5 · 10 −5 to ensure stability. A flowing layer thickness at the midpoint of the flowing layer (δ(0, 0)) of = 0.15 matches the experiments (ω = 2.6 rpm) based on the tracer particle flowing layer passage time. To investigate the impact of the flowing layer thickness, conditions with 0 ≤ ≤ 0.20 were also simulated. For context, previous studies on quasi-2D flows reported flowing layer thickness of ≈ 0.1,(11) depending on the particle to tumbler diameter ratio d/D.
Stroboscopic maps (also known as Poincaré sections or discrete time maps) of 500 iterations were used to investigate mixing and nonmixing behavior and to provide direct comparison to experiments. The stroboscopic maps utilized tracers seeded at the intersection of the flowing layer boundary and the r = 0.9 hemispherical shell at zero iterations (blue points) and at the first half iteration (red points). Regions avoided by the tracer particles correspond to elliptic regions. The outer boundary orbits of these elliptic regions were obtained by finely seeding initial conditions near the boundary in successive stroboscopic maps until the boundary orbits were extracted. The periodicity of these elliptic regions was determined by tracking tracers seeded in these regions. Other larger elliptic orbits that wrap around the bulk and the flowing layer (e.g., grey points in Fig. 4(a-b) and in Extended Data Fig. 1(a)) were extracted by seeding a tracer in the neighbourhood of the orbit.
Piecewise Isometry (PWI) model.. The PWI model is applicable in the infinitely thin flowing layer (ITFL) limit ( = 0) of the FL model. (14)(15)(16) Like the continuum model, stroboscopic maps were used to investigate mixing and non-mixing regions. Because the flow dynamics are radially invariant with an ITFL, tracer positions are mapped by their angular displacements on the hemispherical shell. When tracers reach the free surface, they are instantaneously reflected across the ITFL in the streamwise direction. The initial positions of the passive tracers in the stroboscopic maps were selected to lie along the isometry partitions at zero iterations (blue points) and the first half iteration (red points) and mapped for 10,000 iterations to assure convergence.