Waves, Algebraic Growth and Clumping in Sedimenting Disk Arrays

We study experimentally the Stokesian sedimentation ($Re \sim 10^{-4} $) of a one dimensional lattice of disks in a quasi-two-dimensional geometry with the trajectory of the centres of the disks lying in a plane. We induce initial positional perturbations over a configuration in which the disks are uniformly spaced with their separation vectors and normals aligned, and perpendicular to gravity. Over a range of interparticle separations and density perturbation wavenumbers, we find two regimes of behaviour: (i) a wave-like excitation in orientation and number-density, which grows in amplitude to disrupt the lattice at long times and (ii) a clumping instability resembling that of spheres, but modified by orientations. We construct the equations of motion for displacements and orientations of the lattice using pairwise hydrodynamic interactions and show that the nonnormal nature of the stability matrix is responsible for a new dynamics, in the form of wave-like motion modulated by algebraic growth of perturbations. Linear stability analysis demarcates, in the plane of wavenumber and lattice spacing, the phase boundary between the regimes of nonmodally growing wave-like modes and clumping instability, in quantitative agreement with our experiments.


INTRODUCTION
The collective settling of particles in viscous fluids is a classic and notoriously difficult problem in the physics of strongly interacting driven systems. In the Stokesian limit of Reynolds number Re → 0, inertia is negligible, viscous forces dominate, and a settling particle creates a flow that decays slowly with distance r as 1/r [1][2][3][4]. Furthermore, particles in most natural and industrial settings are not spheres; the hydrodynamics of settling couples their rotational and translational degrees of freedom [5][6][7]. The separation vector of a pair of Stokesian settling spheroids is either time-periodic or asymptotically diverging [8][9][10][11], in sharp contrast to that of two spheres, which is constant [12,13]. Many-body sedimentation [14][15][16][17] thus acquires qualitatively new dynamical features associated with particle geometry; understanding these new features for the many-spheroid problem is the central theme of this work. Concentration fluctuations have been studied in a statistical framework for uniform suspensions of apolar [18][19][20] and polar [21] spheroids. In this article, we study experimentally and theoretically the unexplored case of a sedimenting array of orientable but apolar objects. The array breaks translation invariance and thus offers a reference structure, and consequent dynamics [22][23][24][25][26] , distinct from those of the uniform suspension. In addition, the presence of a lattice clarifies the connection between particle-level interactions and longwavelength collective phenomena such as waves and instabilities.
An array of sedimenting Stokesian spheres with only hydrodynamic interactions is subject to a well-known clumping instability at all wavenumbers [22,23]. This Crowley instability can be understood by the composition of two-body interactions: (i) a locally denser region in the array settles faster due to reduced drag force F D , (ii) the resulting local tilt of the array leads to a lateral drift force F LC , acting along the line joining their centers [see Figure 1  peting mechanism that might stabilize the array, we ask the general question: is a sedimenting lattice of oriented objects stable?
Our experiments find two distinct dynamical regimeswavelike excitations with algebraic growth, with no counterpart in sphere arrays, and linearly unstable modes - separated by a stability boundary in the plane of perturbation wavenumber q and lattice spacing d. Theoretically, we show that the long-wavelength dynamics of an array of uniaxial fore-aft symmetric objects, with long dimension initially aligned vertically, contains terms that compete with Crowley's [22] clumping instability. Explicit construction of the equations of motion for Stokesian sedimenting spheroids, at the level of pair hydrodynamic interactions, determines the values of coefficients in the above coarse-grained theory, and accounts for the experimentally observed behaviour in the q-d plane.
The mode structure of the linearised nearest-neighbour theory, in the limit of thin disks, compares remarkably well with the experimentally measured frequency of the waves. An unusual feature observed in our experiments and in the numerical solution of the far-field equations is FIG. 3. Clumping instability: Overlapped frames (a) t =0, 72, 135 seconds, and (b) t= 0, 120, 255 seconds. In (a) initial horizontal perturbation wavenumber q = 2π/λ = π/2d [same as Fig (2b)] and lattice spacing d = 1.875a, exhibits an unstable case, as line of centers force F LC dominates over orientational drift F θ . In (b) d = 2.5a and qd = π/4, exhibits an unstable case where F LC dominates over F θ , leading to coarsening of the lattice followed by clumping in non-linear regime. The trajectories of disks at nodal points is given by red dashed lines. Clockwise and counterclockwise rotation of disks, depicted as blue and red circular arrows respectively, do not change colour along the trajectory of the disks, in contrast with the wave shown in Fig (2b the nonmodal [27][28][29][30] growth of perturbations in the neutrally stable regime. We term this growth "nonmodal" since it occurs even when all modes of the stability operator are neutral or decaying. The underlying reason is that our dynamical matrix A is nonnormal, i.e., AA † = A † A (where the dagger represents the adjoint). Once the perturbation amplitude due to this nonmodal growth is large enough, nonlinearities can be triggered. We further predict (i) a universal scaled form for the stability boundary, valid across all axisymmetric apolar shapes and (ii) the form of the initial perturbation that leads to maximum transient growth at each point in the neutrally stable regime in the q-d plane.

II. EXPERIMENTS
Our experiments were conducted with disks of radius a=0.4 cm and thickness 1 mm, 3D printed (Form-Labs SLA) with stereolithography using resin of density 1.164 g cm −3 , settling in Silicone oil of density 0.98 cm −3 and kinematic viscosity 5000cSt leading to a typical Reynolds number Re ∼ 10 −4 . The particles are released in a one-dimensional array from the top of a quasitwo-dimensional glass container of width (x-direction) = 225 a , height (z-direction) = 112.5 a , and depth (ydirection) = 12.5 a. As shown in Figure 2 (a), the disks were initialised such that the surface normal K is perpendicular to gravity and parallel to the line joining the disk centres. This was achieved by first placing the disks in slots separated with a centre-to-centre spacing of 0.625a within a frame centred along the depth of the container. The clearance of the disks in the slot sets a precision of 0.0625a in the horizontal position, a deviation of up to 1.8 • in orientation from the vertical, and negligibly small differences in initial vertical positions. Both this frame and the disks are already submerged in the fluid to suppress air bubbles. The disks are then ejected from the slots at the same time with a comb whose teeth fit the slots in the array. The centres of the discs and their surface normals lie in the central (x, z) plane of the experimental geometry for much of their trajectory.
Our reference state is an array of disks with uniform spacing achieved by choosing slots with separation d. On top of this initially uniform lattice, we impose horizontal positional perturbations u x (t = 0) at a wavenumber q by displacing disks appropriately to the right or left slot [see Appendix Figure 6]. The initial perturbations were measured to be sinusoidal with good accuracy, despite the discrete nature of the horizontal displacement. The manual process of poking the disks out of the slots leads to random errors in initial orientations, which contributes to an error of 0.06a in imposed positional perturbation [see Appendix Figure 7].
Images were taken at 1/3 frame per second using a Nikon D750 D-SLR camera. The positions and orientations of the disks were tracked by fitting ellipses to every disk for each image frame. The centroid and angle of the ellipse give the centre positions and orientation of the disks respectively, with precisions of 0.02a and 0.5 deg. The time-dependent amplitude of the positional and orientation perturbation (u x , u z , θ) were measured by fitting a sine wave to the measured particle displacements relative to the reference lattice in the co-moving frame.

III. TWO REGIMES OF DYNAMICS
As we vary lattice spacing, d, and perturbation wavevector q, we experimentally observe two distinct regimes of dynamics, as depicted in Figure 2(b) and 3.
(i) Waves of orientation coupled with number density fluctuations -for the initial condition in Figure 2 (b), we see that the disks at the density nodes hardly rotate, while the orientation and position of disks at the antinodes vary sinusoidally with time [see Appendix figure 7]. Qualitatively, these wave dynamics may be explained by a composition of drag reduction, horizontal glide and mutual rotation as discussed in Figure 1(d). Disks in regions of high number density fall faster than those in less dense regions due to reduced drag. The translational degree of freedom couples with rotations such that the disks in the dense region spread out due to orientational glide, stabilising the lattice. This mechanism is characterized by change in sign of the rotation of disks at the antinodal points, which leads to waves [see Supplementary video 2]. This wave is eventually disrupted [see Supplementary video 3], due to experimental imprecision in the initial orientations [see Appendix Figure 7] that excites nonlinearities via nonmodal growth of perturbations, as shown later in this article.
(ii) Clumping instability decorated with orientations -a different type of dynamics is observed for the initial conditions in Figure 3 where the perturbation quickly sharpens at the displacement nodes, or the high density regions. Just as in the Crowley instability of spheres, the dense regions fall faster due to reduced drag, and the vertical perturbation u z increases. The orientation acts to spread out and rarefy the dense regions, but this effect is suppressed by the line of centers force leading to a clumping instability. The rotation of the antinodal points does not change sign, in contrast with the wave-like regime. Figure 3 (a) depicts a marginally unstable case where the initial horizontal perturbation neither grows nor decays substantially whereas in Figure 3 (b) the horizontal perturbations grows to make dense region more dense [see Supplementary video 4 & 5]. Later in the article, we show experimentally the regime of each of these two types of dynamics by varying initial conditions in the q-d space.
Theoretically, we go beyond our qualitative explanation above at two levels. First, we understand our experimental observations using symmetries of the equations of motion for displacements and orientations, in the continuum limit of our system. Second, to determine the phenomenological coefficients in the symmetry-based equations for Stokesian sedimentation, we construct the dynamical equations of the lattice using pairwise addition of forces and torques resulting from the hydrodynamic interactions. We then show that the linearized dispersion relation of our theory compares quantitatively well with our experiments, while long time nonlinear instabilities can be understood by numerical investigations of the farfield equations of motion.

IV. SEDIMENTING SPHEROID LATTICE: SYMMETRY-BASED CONTINUUM THEORY
We construct the "hydrodynamic" equations of motion of a drifting lattice of oriented objects, in the limit of no inertia, by writing the most general form of the mobility tensor (defined by velocity = mobility × force) allowed by the symmetries of the system, to leading order in a gradient expansion. We find thus that the dynamical response of a lattice of orientable particles when perturbed about a suitable reference state contains terms that can compete with the clumping instability of isotropic particles [22,24]. We discuss the structure of the resulting wavelike modes.
The configurations of a periodic lattice of uniaxial objects are characterized, in a coarse-grained Eulerian description, by the displacement field u of the lattice and the orientation field K defined by the mean local alignment of the particle axes. For our geometry [see Figure 4 (a)] K = (cos θ, 0, sin θ), with θ equivalent to θ + π because the particles are fore-aft symmetric. The equations of motion for u and K, in the presence of a gravitational driving force F, must obey the following symmetries: • Translational invariance • Rotational invariance in the subspace perpendicular to gravity • Symmetry under inversion of orientations, K → −K The mobility cannot depend directly on u due to translational invariance, but dependence on ∇u, K and ∇K is allowed: where M and N are the translational and rotational mobilities respectively, and P ≡ I − KK is the projector transverse to the unit vector K. The other symmetries further constrain the allowed form of translational and rotational mobilities, leading, at lowest order in gradients, for a one-dimensional lattice along x, in a comoving frame, to Here, λ i , α and γ depend on F and the parameters governing the mobilities in (1) and (2). Substituting K = (cos θ, 0, sin θ) in (3)- (5) and linearizing about θ = 0, the state where the particle axes are along x (Fig. 2a) leads, for disturbances with frequency ω and wavenumber q, to the dispersion relations with elasticity contributing to (5) and (6) at order q 2 . For α → 0 the linearized equations for the translational degrees of freedom (u x , u z ) are independent of K and reduce to those of the LR model [24], with wavelike modes or an instability depending on the sign of λ 1 λ 2 [24][25][26]. u affects K through the one-way coupling governed by γ.
For α = 0, translation and rotation are coupled, and the presence of αγ in the dispersion relation opens up the possibility of linearly stable wavelike dynamics even for λ 1 λ 2 < 0. The linearized dynamics about the state where K is vertical corresponds to changing the sign of α in (6).
For a system of sedimenting particles this means that the array is stable either with horizontal orientations or vertical orientations, but not both. Similar considerations arise in principle for the stability and dynamics [26] of driven flux lattices in thin slabs of type-II superconductors if the cross-sections of the flux lines are non-circular.

V. SEDIMENTING SPHEROID LATTICE: PAIR HYDRODYNAMIC INTERACTIONS
We now go beyond symmetry considerations, and explicitly construct the equations of motion for a settling lattice based on single-particle motion and addition of pairwise forces and torques at each particle position. We develop the theory for an array of spheroids, of eccentricity e = 1 − b 2 /a 2 , where a and b are the semi-major and semi-minor axes respectively. In the limit of e → 1, an oblate spheroid approaches a disk shape, as in our experiments. We consider hydrodynamic interactions to leading order in a/r, where r is the separation between two particles. The ingredients of array dynamics are: (i) Lateral drift of a single particle -An isolated settling spheroid cannot rotate, thanks to Stokesian time-reversal symmetry, but drifts horizontally with velocity sin 2θ (7) [2,7] where F is its buoyant weight, µ is the dynamic viscosity of the fluid and the mobility α is a function of the eccentricity. Figure 4(a) shows a schematic of a portion of our array, in which the orientation vector K n of the nth particle is defined as a unit vector along the minor (major) axis for an oblate (prolate) spheroid. The angle θ n is measured from the vertical as shown.
(ii) Mutual drag reduction -Two particles at finite separation fall faster than an isolated one, due to the addition of the flow fields generated by each Stokes monopole [2,22]. In the far-field approximation, the increased ver- x ,u n x ) and orientation perturbation θ n of the n th disk which interacts hydrodynamically with the neighbours n − i, i = 1, 2, .. . For prolate spheroids the orientation vector K n is rotated by π/2 from the one shown. (b) The phase diagram with the stable regime shown in blue and unstable regime in red. The experimental data points (circles) are coloured blue or red by measuring whether the density autocorrelation grows or decays, which is shown in the inset to (b) on the top right, for some representative stable (blue) and unstable (red) points in the q − d plane. The phase boundary predicted by the linear theory with nearest-neighbour interaction is shown for oblate spheroids with eccentricity approximating the experimental thickness, e = 0.9922 (solid line) and and elliposid of zero thickness e = 1 (dashed line). tical velocity to leading order in a/r is (iii) Horizontal drift -The flow generated by the neighbouring particle gives rise to a force along the line joining the centers of the two particles [2,22], which leads to a horizontal component of velocity to leading order in a/r. (iv) Mutual rotational coupling -The presence of a neighbouring particle generates a velocity field of nonzero vorticity, which to leading order in a/r gives a rotationθ We combine ingredients (i) to (iv) to build the dynamics of the array of spheroids. A.
Mode structure for oblate and prolate spheroids We consider an infinite one-dimensional reference lattice along the x-axis of uniformly spaced lattice points with spacing d and falling in the −z direction. As shown in Figure 4(a), we consider identical spheroids with orientation θ n , and centroids at a small displacement (u n x , u n z ) measured from each lattice point, where the superscript n stands for the n th particle. In the mean settling frame, pairwise addition of forces and torques on the n th particle due to hydrodynamic interactions with the (n + l) th and (n − l) th particles, for l = 1, 2, 3...∞, gives the equation of motion of the n th particle as In a quasi two-dimensional geometry the dynamics can be approximated by a nearest-neighbour treatment, where the n th particle interacts hydrodynamically only with the n + 1 th and n − 1 th particle. We nondimensionalise equations (11) -(13) using the lattice separation d and time scale T * = µd 2 /F , and perturb the angle θ = 0 + δθ. Linearising the equations and fourier transforming gives the equationẊ q = A(q)X q , where X q = (u x q , u y q , δθ q ) is the spatial fourier transform of the perturbations with wavenumber q along x, with a nonnormal dynamical matrix: We return to the interesting consequences of the nonnormality of A(q) later. For now, we substitute the translational mobility function [3,7,11] for oblate spheroids: α(e) = 9 − 6e 2 tan −1 e/ √ 1 − e 2 − 9e √ 1 − e 2 /8e 3 , and for prolate spheroids: α(e) = 9 − 3e 2 ln [(e + 1)/(1 − e)] − 18e /16e 3 , which gives the mode structure with two branches around ω = 0 for each: Oblate Spheroids: Prolate Spheroids: In the limit of e → 0, these dispersion relations for both prolate and oblate spheroids correctly reduce to iω ± = ±| sin(q)|/4π, which is just the Crowley instability for spheres. For e = 0, defining the nondimensional quantitỹ d ≡ 2dα(e)/3a, gives a universal condition for stability: so thatd = q cos 2 /2 defines the stability boundary in thed-q plane, separating the regime of kinematic waves (blue) from the clumping instability (red) as shown in the phase diagram of Figure 4(b). This prediction agrees well with our experimental data shown by the red and blue circles, where we have initialised the lattice at those points in thed-q plane. The outcome of any given experiment is identified as being wave-like or clumping by considering the early stages of the time-dependence of ρ(t)δρ(t = 0) , normalized amplitude of the density autocorrelation, which is measured by projecting the particle number density ρ(t) = N m=1 δ(x − x m (t))/N , on the initial density fluctuation δρ(t = 0) of the lattice. We obtained δρ(t = 0) by fitting a sine to the initial horizontal displacement perturbation u x (t = 0), and shifting in phase by π/2. This is shown in the inset to Figure 4(b), where some curves increase in amplitude, and others decay. At later times, even in the wave-like regime, the perturbation becomes very non-sinusoidal, as nonlinear effects become prominent.
More specifically, the limit of disks with zero thickness (e → 1 for oblate spheroids), produces the mode struc-ture shown in Figure 5(a): which gives neutrally stable modes when the lattice spacing d > 8a cos 2 (q/2)/π and clumping instability otherwise. This prediction is compared with experimental data for the frequency in Figure 5(c) for various q and d.
We show solutions corresponding both to zero thickness, as well as for the ellipsoid with 2a and 2b corresponding to the diameter and thickness of our disks. In the long wavelength limit q → 0, (18) reduces to the dispersion relation (6) predicted by that we showed in the previous section based on symmetry arguments. The mobility coefficients for disks are determined to be: λ 1 λ 2 = −1/16π 2 and αγ = d/128aπ.
The limit for needles of zero thickness (prolate spheroids with e → 1) is not well defined, but the dispersion relation for rods of small thickness 2b and length 2a, to leading order in 1 − e is: B.

Waves and nonmodal behaviour of disks
For disks, the eigenfunctions of the dynamical matrix (14) can be used to construct the solution for experimental initial conditions. The eigenvectors corresponding to FIG. 5. Spectrum:(a) The branches of the spectrum. The absolute value of the complex frequency is shown as a function of the nondimensional wavenumber q and spacing d. The blue surface is for purely real frequencies (the neutrally stable modes) and the red surface is for purely imaginary frequencies (clumping instability). The boundary between the stable and unstable regimes is shown by the magenta curve dc. (b) The log of ratio G/G0, of maximum amplitude G of the nonmodal perturbation to the initial amplitude G0 depicted for the neutrally stable regime in the q * d − d/a plane. (c) Measured non-dimensional frequency ωT * (symbols) for various wavenumbers q plotted against the lattice spacing d/a. The curves are predicted by our nearest-neighbour theory for oblate spheroids, solid lines for the experimental eccentricity e = 0.9922 and dashed lines for the limit of zero thickness e = 1. (d) The experimental amplitude for the angle (red patch), horizontal perturbation (blue patch) and vertical perturbation (yellow patch) compared against theory (solid lines) and far-field simulation of oblate spheroids in limit of zero thickness (dashed lines), for the stable case of q = π/2 and d/a = 3.75. The extent of the shaded region shows the corresponding error in measurement of amplitude. (e) The nonmodal growth plotted for the stable case q = π/2 and d/a = 3.75 in experiment (blue), compared with the simulation for the initial condition that gives maximum gain (red) in t = 14T * , where G0 is the initial amplitude.
The time-dependence predicted here is in good agreement with the experimental data shown in Figure 5  The non-modal nature of the dynamical matrix (14) gives rise to non-orthogonal eigenvectors, resulting in non-modal growth of perturbations. So even in the 'stable' regime of the phase diagram, perturbations show transient algebraic growth. If the transient amplitude is large enough, nonlinear growth takes over, as in our experiment. In our far-field solutions on the other hand, we have the facility to reduce the initial amplitude so much that despite transient growth the system remains linear. To quantify the nonmodal growth, we calculate the norm of exp[As] for all times s, and calculate the maximum amplitude G max attained by the perturbation over all s, which is finite in the stable regime and depends on wavenumber q and lattice spacing d [see Fig  5 (b)]. We observe significant nonmodal growth for the neutrally stable mode of experimental perturbations [see Fig 5 (e)]. The singular value decomposition of exp(As) provides the initial condition which gives maximum nonmodal growth, which we compare against the experimental perturbation using far-field simulation of oblate spheroids in the limit of disk e → 1. Further, our numerical study of the far-field equations [see Appendix D] shows that the observed disruption of the lattice in the stable regime results from amplification of the experimental noise in the initial orientations [see Supplementary videos 3 & 7].
We see thus that even in the regime where the orienta-tional degree of freedom defeats the Crowley mechanism, and linear stability predicts waves, transient growth ultimately triumphs. An array of sedimenting spheroids is thus disrupted at all q and d. In our numerical study with periodic boundary conditions we are able to observe the waves and delay the onset of nonlinearity by reducing the amplitude of the initial perturbations, unlike in the experiments where there are inevitable imprecisions in the initial conditions.

VI. CONCLUSIONS
We have shown the existence of two regimes of dynamical behaviour, as a function of lattice spacing and perturbation wavenumber, in the Stokesian sedimentation of a one-dimensional lattice of disks. One of these is an extension of Crowley's clumping instability [22] to nonspherical particles. The second is a hitherto unknown state of orientation and displacement waves, where the drift and mutual interaction of the disks overcome the clumping instability. This competition between orientation and clumping in spheroids is related to an effect predicted for polar particles [21], and suggests a new consideration that must be included in the statistical theory of Koch and Shaqfeh [18].
The wave-like regime is unusual in that we predict, and observe experimentally, large transient growth that ultimately destabilizes the lattice, through nonlinear effects arising from initial experimental error in release. Thus, the lattice is nonlinearly unstable over the entire q − d plane, but due to two very different mechanisms. The fact that our calculation, and the accompanying numerics, allow us to capture both the mode-structure and the growth of perturbation amplitude, reassure us that our numerical model can in future be used to gain a comprehensive understanding of the unstable regime, and of other lattice configurations.
We close with an interesting formal observation: equations (11) and (13) governing the dynamics of the settling array are formally identical to those for the displacement and momentum-density fields respectively of a momentum-conserving lattice of masses and springs. Conservation of the total "momentum" n θ n in (13) and the resulting gapless character (ω → 0 as q → 0) of the modes is a consequence of the orientationindependence of the gravitational energy of the apolar disks. Objects with polar shape would have a preferred orientation in a gravitational field [6,31,32] and hence a damping of the "momentum" in (13). The centre-to-centre distance between adjacent slots is δ and the lattice spacing of reference lattice is d. Initial horizontal positional perturbation ux are shown by fitting a sine function (blue curve) to the measured perturbations (red dots) for the wavelengths λ = 4d, 6d, 8d and 12d.
FIG. 7. Error in initial orientation: (a) The error is quantified by 0.5×(maximum -minimum) of the initial angles in any given run. This measure of the width in release angle is plotted here as a function of lattice spacing for qd = π/2, π/3, π/4 and π/6. (b) a combined distribution of, initial angles of several experimental runs for a fixed qd = π/2 and d = 3.75a.
ure. This perturbation was made as close to sinusoidal as possible within the constrain of discretization imposed by the slots [see Fig (6)]. After arranging the disks across the total length of the release mechanism of 80 cm, the disks were poked out gently while the whole mechanism was submerged roughly 3 cm below the surface of the fluid, to avoid any bubbles. After the disks were released out of the slots, they were measured to have a random orientation error sitting on the spatial perturbation which we imposed [see Fig(7)]. This angular error corresponds to an error in horizontal spatial perturbation in u x of ±0.03 cm. This error in release plays a crucial role in disrupting the lattice at late times in the linearly stable, but transiently growing, regime.

Measuring frequency
The wave nature is evident in the dynamics of the orientations and positional perturbations which exhibit a quarter cycle of the wave with reasonable accuracy, before the non-linear instabilities kick in via an algebraic growth of perturbations, in contrast with the exponential growth of perturbations in the unstable regime. The positions and orientations of the disks are measured at each frame every 3 seconds by fitting an ellipse around each disk. The reference lattice is constructed by measuring the largest node to node distance in the initial condition. It is assumed that this reference lattice settles down vertically with the mean settling speed of the lattice and the orientation and positional perturbations are measured for the particles corresponding to each lattice point. A sine wave is fitted in the perturbation with specified wavenumber at each frame and the amplitude of the fitted wave is measured at each time step [see Fig(8)]. The residual of this fit gives error in frequency measurements as the amplitude for u x and θ is plotted as a function of time. The gradient expansion of the translational mobility M(∇u, K, ∇K) and rotational mobility N(∇u, K, ∇K), to leading orders in gradients gives Here M 0 is the mobility of the undistorted lattice and such a term in not allowed in N due to symmetry under K → −K; and P ≡ I − KK is the projector transverse to the unit vector K. In the first term of (B2), is the Levi-Civita tensor. Retaining only those terms that are allowed by the symmetries, leads to the "hydrodynamic" equations for the displacement field u and orientation field K in one dimension x by dropping z derivatives (3) -(5).

APPENDIX C: WAVE SOLUTIONS FOR SPHEROIDS
The eigenfunctions corresponding to the eigenvalues (λ 1 , λ 2 , λ 3 ) = (0, −iω, iω), where ω ≡ ω + is form (15)- (16), are given by v 1 , v 2 and v 3 respectively giving the solution as a real part here X n = (u n x , u n y , δθ n q ) . For a ≡ (a 1 , a 2 , a 3 ), (C2) becomes X n (t) = B · a, where in the stable regime −ωπ csc 2 q 2 sin(qn − ωt) ωπ csc 2 q 2 sin(qn + ωt) 2dπα(e) 3a csc(q) sin(qn) 1 2 cot q 2 sin(qn − ωt) 1 2 cot q 2 sin(qn + ωt) cos(qn) cos(qn − ωt) cos(qn + ωt) The coefficients a can be determined from the initial condition, a = B −1 · X n | t=0 . Our experimental initial condition is X n (t = 0) = ( sin(qn), 0, 0) , making a = sin 2 (q/2) 2ωπ (0, −1, 1) , which gives the wave solution in stable regime Note that the dependence on eccentricity of the spheroids enters through ω from (15)  To understand the non-linear dynamics of disks in (x, z) plane, we numerically analyse the equations of motion for spheroids in the limiting case of disks e → 1, to leading order in O(a/r), by pairwise addition of hydrodynamic interactions using the method of reflections [11]. We simulate the following equations for positions (x n , z n ) and orientations θ n of the n th spheroid with using fourth order Runge-Kutta method: In the nearest-neighbour approximation, number of interacting neighbours N truncates the spatial summation over m, making implementation of periodic boundaries straightforward. Note that the above equations are non-dimensionalized using length scale d and time scale µd 2 /F . Also, the initial conditions are such that the orientation vector of all the spheroids lie in the (x, z) plane and hence the resulting trajectories are confined to the same plane y = 0. Five spheres of diameter 0.6 cm prepared in an array perturbed around an equally spaced configuration with an interparticle spacing 1.5 ± 0.1 cm. The initial perturbation is like that of Fig1 (b) with amplitude 0.25 ± 0.05 cm. Trajectories of the nodes of this perturbation is shown in red. The three-sphere dynamics at later times is expected to be chaotic [33].

Video 2: Linearly stable wavelike mode
Initial sinusoidal perturbation with d = 3.75a, qd = π/2, and amplitude 0.625a; and with trajectory of nodes shown by the dashed red lines. We zoom in on a region where initial errors in release were small, which shows a half cycle of the wavelike oscillation in orientations and positions. More details are in Fig.2(b).

Video 3: Disruption of waves at late times
This video shows the late-time dynamics for the same q and d as in Video 2. Transient algebraic growth of the perturbations leads to nonlinear effects that disrupt the array.

Video 5: Clumping dynamics at late times.
Late time clumping behaviour of the perturbation with d = 1.875a, qd = π/6 and amplitude 0.625a. Trajectories of nodes are shown in red color. 6. Video 6: Numerical study of wave-like regime Numerical integration of the non-dimensionalised farfield equations [see Appendix D] with initial sinusoidal perturbation of q = π/2. The interaction is cut-off beyond 1.5d, such that only nearest neighbours interact hydrodynamically. The region shown here is the same size as the experimental container, scaled by lattice spacing. The initial conditions is the same as in the experiment of Video 2, albeit with periodic boundary condition and no experimental error in initial condition [see The initial condition is the same as in Video 6, but we add a random error in the initial orientations uniformly randomly distributed between ±8 deg, to reflect the measured experimental initial conditions of Video 3 [see Fig 7].