{'fraction_non_alphanumeric': 0.03523267596850885, 'fraction_numerical': 0.009821498168212643, 'mean_word_length': 4.544223566695477, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Active and driven hydrodynamic crystals

Motivated by the experimental ability to produce monodisperse particles in microfluidic devices, we study theoretically the hydrodynamic stability of driven and active crystals. We first recall the theoretical tools allowing to quantify the dynamics of elongated particles in a confined fluid. In this regime hydrodynamic interactions between particles arise from a superposition of potential dipolar singularities. We exploit this feature to derive the equations of motion for the particle positions and orientations. After showing that all five planar Bravais lattices are stationary solutions of the equations of motion, we consider separately the case where the particles are passively driven by an external force, and the situation where they are self-propelling. We first demonstrate that phonon modes propagate in driven crystals, which are always marginally stable. The spatial structure of the eigenmodes depend solely on the symmetries of the lattices, and on the orientation of the driving force. For active crystals, the stability of the particle positions and orientations depends not only on the symmetry of the crystals but also on the perturbation wavelengths and on the crystal density. Unlike unconfined fluids, the stability of active crystals is independent of the nature of the propulsion mechanism at the single particle level. The square and rectangular lattices are found to be linearly unstable at short wavelengths provided the volume fraction of the crystals is high enough. Differently, hexagonal, oblique, and face-centered crystals are always unstable. Our work provides a theoretical basis for future experimental work on flowing microfluidic crystals.


Introduction
The dynamics of passive suspensions is a field with a long history in physical hydrodynamics. Much effort has been devoted to understand, e.g., the origin of fluctuations in the sedimentation of spheres under gravity as well as instabilities in suspensions of elongated fibers (see reviews in [1,2] and references therein). More recently, a significant experimental [3,4,5] and theoretical [6,7,8] research effort has focused on the dynamics of active suspensions where instead of having particles driven by an external field (e.g. gravity), one considers the dynamics and interactions of self-propelled synthetic or biological swimmers. In this case, the interplay of activity and hydrodynamic interactions leads to long-wavelength instabilities [9,10].
Most of the past work on passive (driven) and active suspensions has focused on instabilities and fluctuating behavior in three-dimensional systems. However, over the last ten years microfluidics has offered a number of simple and effective solutions to produce and manipulate large ensemble of highly monodisperse microparticles, prone to form crystal structure in quasi-two dimensional channels [11]. For driven particles, these technological advances have motivated, for example, the study of the nonlinear dynamics of finite flowing crystals [12,13], phonons in one dimensional microfluic-droplet crystals [14] and flow-ing lattices of bubbles [15]. In the case of active particles, these fabrication methods could be extended to selfpropelled catalytic colloids [16,17] or reactive droplets [18].
Motivated by these advances, we take in this paper an approach contrasting with the traditional study of disordered suspensions and consider the dynamics of confined driven and active hydrodynamic crystals. We first develop a formalism to study theoretically position and orientation instabilities for flowing discrete suspensions under confinement. In the case of driven particles, we demonstrate formally that all crystals are marginally stable and study in detail the eigenmodes of deformation for all five two-dimensional Bravais lattices. For active particles, we show that square and rectangular crystals are linearly unstable at short wavelengths provided the volume fraction of the crystals is high enough. Differently, hexagonal and oblique (resp. face-centered) crystals are always unstable for long-(resp. short-) wavelength perturbations. In contrast with past work on three-dimensional swimmers suspensions, the stability of confined active crystals is found to be independent of the pusher vs. puller nature of the actuation of individual active particles [19]. 2 Theoretical setup

Particle crystal in a confined fluid
We start by describing the theoretical framework we use to quantify the large scale dynamics of both active and driven microfluidic crystals. We focus our study on the case of identical particles living in quasi-bidimensional fluids, as sketched in Fig. 1. The fluid is Newtonian and has a homogeneous thickness h in the z-direction, comparable to the size of the particles. Our formalism will be valid both for thin films lying on a solid substrate (with one free surface and one no-slip wall), as well as for microfluidic geometries where the fluid is confined between two parallel plates. The particles can be either axisymmetric or anisotropic and are organized in two-dimensional crystal, see Fig. 1. If a denotes the typical lattice spacing of the crystal and b the typical extent of the particle in the (x, y) plane, we consider in this paper the dynamics in the dilute limit, e.g. a b. In this limit, each particle i is appropriately modeled as a pointwise body characterized by its in-plane position, R i (t) ≡ (x i (t), y i (t)), and its inplane orientation,p i (t), wherep i is a unit vector making an angle θ i (t) with thex-axis. Having microscopic systems in mind, we neglect the particle inertia and work in the limit of zero Reynolds number. In this Hele-Shaw setup, it is a classical result that the fluid flow is potential [20]. The z-averaged fluid velocity, V and the z-averaged pressure, P , are therefore related by where G = αh 2 /η; here η is the fluid viscosity, and α = 1/3 for a thin film, and α = 1/12 for a shallow microchannel. Together with incompressibility, ∇ · V = 0, Eq. 1 determines the fluid flow and stress away from the particles.
We henceforth consider either swimmers moving along their principal axisp i , or passive particles driven by a uniform force field oriented along the x-direction (gravitational, electrostatic, magnetic,...). In all cases, the speed of an isolated particle in a quiescent fluid is constant and denoted U 0 . In addition to their individual dynamics, particles also follow the surrounding flow, and the equation of motion for particle i thus reads whereq =p i for swimmers,q =x for driven particles, and µ is a non-dimensional mobility coefficient [14,21]. Passive tracers have µ = 1. Conversely, for thick particles, the friction against the solid wall(s) can significantly reduce the advection speed, which is smaller that the local fluid velocity, and thus 0 < µ < 1. In principle, µ should be a tensor for anisotropic particles but for simplicity we consider only particles which are weakly anisotropic and thus µ is assumed to remain a scalar 1 . In addition to a change in their velocity, anisotropic particles experience hydrodynamic torques which favor an orientation along the local elongation axis of the flow. This classical hydrodynamic result, which can also be anticipated from symmetry arguments, leads to the so-called Jeffery's orbits [22]. As the flow is irrotational (potential flow), the orientational dynamics reduces to where E is the strain rate tensor, E = 1 2 ∇V + (∇V) T , and γ ≥ 0 is a rotational mobility coefficient which is nonzero for anisotropic particles and zero for axisymmetric bodies.

Long-range hydrodynamic interactions
As a particle located at R(t) moves in the fluid, it induces a far-field velocity, denoted v(r − R), at position r. A given particle i responds to the flow induced by all the other particles in the crystal, and is therefore advected at velocity V(R i ) = µ j =i v(R i − R j ). In this section we provide a quantitative description of the far-field hydrodynamic coupling, v, between the particles; for the sake of clarity, we separate the case of passive and active particles.

Hydrodynamic interactions between driven particles
In the driven case, each particle of the crystal is subject to a constant external force, f = fx, which results in a farfield perturbation which we denote v 1 and is the Green's function of Eq. 1. In three-dimensional flows, the response to a force monopole is known as a Stokeslet, and decays spatially as ∼ 1/r. In our quasi-2D geometries, solid walls act as momentum sinks and screen algebraically the Stokeslet contribution, which then decays as v 1 ∼ 1/r 2 and takes the functional form of a potential source dipole, as shown in [23,24]. In addition, the particles have a finite size and their advection by the surrounding fluid is hindered by the lubrication forces induced by the confining walls (even in the absence of external driving). Due to incompressibility, any relative motion with respect to the fluid results in another algebraic far field contribution, v 2 . As shown, e.g. in [14], v 2 has also the form of a potential dipole with the same spatial decay, v 2 ∼ 1/r 2 . (We note that in unbounded fluids, this potential contribution scales as 1/r 3 and is thus subdominant with respect to the flow induced by a pointwise force, which decays as 1/r). Therefore, in confined flows, the two contributions have the same form [14,23] and the overall flow disturbance, where r = |r| and the dipole strength, σ, is the sum of the two contributions, σ = Ab 2 Gf + Bb 2 U 0 , where A and B are two dimensionless shape factors (I is the identity tensor). The symmetry of the streamlines for this flow field are illustrated in Fig. 2 (left).

Hydrodynamic interactions between active swimmers
By definition swimmers do not require an external force to propel themselves. The stress distribution on the surface of a self-propelled particle has thus, at least, the symmetry of a force dipole [25]. The canonical theoretical setup used to describe (dilute) suspensions of swimmer is to consider an ensemble of such force-dipoles as all other multipolar contributions to the far field are subdominant in an unbounded fluid [9,19,26]. However, as mentioned above, confinement results in an algebraic screening of the hydrodynamic interactions. In the quasi-2D geometry at the center of our paper a force dipole decays spatially as ∼ 1/r 3 , a contribution which is therefore subdominant compared to the ∼ 1/r 2 potential dipole arising from incompressibility (similarly to driven particles) [23]. For active swimmers, the far-field flow disturbance has thus also the symmetry of a potential source dipole, the difference with the passive case being that the dipole direction is now the swimmer orientation (Fig. 2, right). For a swimmer orientated alongp, we obtain a flow given by v s (r,p) = σ 2πr 2 (2rr − I) ·p, with the dipole strength σ = Bb 2 U 0 (B is the same shape factor as in Eq. 4). We therefore see that, that in confined fluids, the usual distinction between pushers and pullers swimmers (contractile and extensile), which is at the heart of qualitatively different behaviors in unconfined fluids [7,8], is irrelevant. The magnitude and sign of the induced dipolar flow is solely set by that of the swimming speed, irrespective of the microscopic swimming mechanism.
In summary, Eqs. 2, 3, and either Eq. 4 (in the driven case) or 5 (active case) fully prescribe the dynamics of the discrete particle positions and orientations. As noted above, the main difference between active and passive particles concerns the orientation of the dipolar flow field: the orientation is slaved to the swimmer direction for active particles whereas it is constant and aligned along the xdirection for driven particles (the difference is further illustrated in Fig. 2). We will show in the next sections that this distinction markedly impacts the large-scale crystal dynamics.

Are hydrodynamic crystals stationary?
When addressing the dynamics of an ordered phase, the first important question is whether this phase does correspond to a stationary state. We focus here on the five planar Bravais crystals, which encompass all possible symmetries for bidimensional mono-atomic crystals (see Fig. 3).
Let us first consider the case of crystals composed of driven axisymmetric particles. The equations of motion reduce to where the R i 's belong to one of the Bravais crystals from Fig. 3. The lattice structure is conserved provided that By definition, all crystals are invariant upon translation along (R i − R j ), which readily implies that the two sums are equal, and therefore that any driven crystal made of axisymmetric particles is stationary structure To extend this result to driven crystals composed of anisotropic particles we first recall that all the Bravais lattices are invariant upon the parity transformation r → −r. Moreover, as v d (r) is invariant upon this transformation whereas the sign of the gradient operator is reversed, we see that the strain rate tensor constructed from a superposition of potential source dipoles, Eq. 4, transforms according to E(−r) = −E(r). This implies that for any Bravais Fig. 3. Geometry of the five planar Bravais lattices. Anisotropic cells are characterized by the ratio, , between the two cell dimensions. The angle β is the tilt angle of the oblique and hexagonal cells. For each particle labeled "0" we also display and number all nearest neighbors.
crystal, E(R i ) is identically zero anywhere on the lattice. It follows from the equation for the orientational dynamics, Eq. 3, that ∂ tpi = 0 for all the particles. In conclusion, in the driven case, both the crystalline structure and the particle orientations remain stationary (in other words, the crystals are fixed points of the dynamical system).
It is straightforward to generalize the above results to swimmer crystals. The equation of motion for the positions, Eq. 6, is given by Obviously, ∂ t (R i − R j ) cannot be zero if thep i are not all identical. Therefore, the crystal structure cannot be conserved if the initial orientation of the particles is not uniform -in such cases the crystal would "melt". For uniform orientations, say alongx, Eqs. 6 and 7 are identical, and so is the equation for the orientational dynamics since v s (r,p i ) = v d (r). We are thus left with the same problem as in the driven case, which implies that the structure of the crystals is conserved as long as the particles all swim along the same direction.

Driven hydrodynamic crystals are marginally stable
We start by investigating in this section the linear stability of the five Bravais crystals with respect to perturbations in both the position and the orientation of the particles, with a special focus on the experimentally-relevant square and hexagonal lattices. Anticipating on our results, we note that the geometrical classification in term of the Bravais lattices might not necessarily be relevant to the dynamics of flowing crystals.
In order to proceed we make use of two additional assumptions. Firstly, we consider the case of particles uniformly aligned along thex-axis prior to the perturbations, as depicted in Fig. 3. Secondly we assume that the driving force is aligned with one of the principal axes of the crystal. Our following study can be easily extended to a more general setup, but this would make the formula and the discussions much more tedious.
We denote δR i and δp i ∼ θ iŷ the infinitesimal perturbations of the particle positions and orientations respectively, so that R i → R i + δR i , andp i →x + θ iŷ . Using the property that E = 0 for dipoles organized into a Bravais lattice (as discussed in the previous section), and after some algebra, the linearization of the equations of motion, Eqs. 3 and 6, yields and, where R ij = R i − R j , and δR ij = δR i − δR j . Eqs. 8 and 9 dictate the dynamics of the elementary excitations in the frame where the unperturbed crystal is stationary. We note that the direction of the crystal translation is, in general, different from the driving direction. We now exploit the symmetries of the dipolar interactions. Inspecting the flow given by Eq. 4, we deduce Using these relations, we look for plane waves solutions, (δX i , δY i , θ i ) ≡ (δX, δY, θ) exp(iωt − iq · R i ). By doing so, we obtain a linear-stability system, which we write in the generic form where the coefficients of the stability matrix M are We readily deduce from the matrix structure that a perturbation in orientations only would not induce any change in the crystal conformation. This is a direct consequence of the dipolar coupling between the particles, v d , which is only a function of the driving force direction and not of the particle orientation. On the contrary, perturbations in the position of a particle modify both position and orientation. In addition, perturbations in orientation only neither relax, grow or propagate. As the third column of the matrix M is always 0, this implies that it will always admit the eigenvalue ω 0 = 0, associated to the pure-orientation eigenmode (0, 0, 1).
The other two eigenvalues of the M matrix are ω ± = ± M 2 1 + M 2 2 . Exploiting again the fact that all Bravais lattices are invariant upon parity transformation to write By definition R ij = −R ji , and due to the dipolar symmetry of the hydrodynamic interaction, we have ∂ . Using these two equalities in Eq. 15, we infer that , and therefore M 1 is a real number. Using the same method, and the identity ∂ , one can show that M 2 is real as well. Therefore, for any symmetry of the crystal, the pulsations of the plane waves, ω ± , are real. In other words, for any Bravais lattice the crystal structure of driven particles is dynamically marginally stable.
Notably, the dipole strength σ can be eliminated from the equations of motion Eqs. 8 and 9 by rescaling the timescale. Therefore, the linear stability of the monocrystals is a purely geometrical problem. The corresponding eigenmodes do not depend on the translational speed U 0 , but only on the orientation, and on the symmetries of the lattice.
Interestingly, we see that phonons propagate with the pulsations ω ± , despite the fact that particles have no inertia and that no potential forces couple the particle displacements. This seemingly counterintuitive result generalizes the experimental observations made by Beatus and coworkers in [14] where they revealed that sound modes propagate along 1D-droplet crystals flowing in quasi-2D microchannels. These results are, importantly, specific to the quasi-2D geometry, which is relevant for numerous microfluidic and thin films applications. In unbounded fluids, the change in the symmetry of the hydrodynamic interactions results in the destabilization of the crystal structure as shown theoretically and experimentally [27].
Below, we derive the dispersion relation for each of the five Bravais crystals, with a special attention given to the case of square and hexagonal lattices

Square lattice
In order to compute the coefficients of the M-matrix analytically we now make a nearest neighbor approximation. In a similar context, this approximation has proven to yield qualitatively correct results for unbounded fluids [27]. We introduce a reference particle labelled as 0. The four nearest neighbors in the square crystal are labeled as Fig. 4. Normalized dispersion relation for the square lattice plotted from Eq. 17 (ω + only), for 2µσ/πa 2 = 1, and a = 1.
1, 2, 3, and 4 ( Fig. 3). In this geometry, we easily compute the coefficient of the M matrix as This dispersion relation is plotted in Fig. 4.
To gain insight into the propagating modes, we focus on the large-scale (long wavelength) response of the crystals. Expanding Eq. 16 at leading order in the wave vector amplitude for q → 0, we find We notice that in this small-q limit, the orientation and the position degrees of freedom are totally decoupled. The three eigenvalues are ω 0 = 0, and ω ± = ±2µσq/(πa 2 ). The two non-trivial modes are non-dispersive and propagate with a constant "sound-velocity" c ± = ±2µσ/(πa 2 ), which increases with the magnitude of the hydrodynamic coupling.
To understand physically how the excitations propagate, we focus on two specific cases. Let us first consider longitudinal perturbations, q = qx. The mode ω − is here associated with the eigenvector (0, 1, 0). It corresponds to shear waves which propagate in the direction opposite to the driving, as illustrated in Fig. 5A. The second sound mode (ω + ) corresponds to compression waves along the x-axis propagating in the driving direction, see Fig. 5B. The corresponding eigenvector is (1, 0, 0). For excitations propagating in the direction transverse to the driving, q = qŷ, the eigenmodes couple the displacements along the two principal axes of the crystal. The mode ω − is associated with the eigenvector (1, 1, 0). It corresponds to the superposition of a compression mode in theŷ direction, in phase with a shear in thex direction. The second mode (ω + ), with eigenvector (−1, 1, 0), is a combination of a dilation in theŷ direction, which propagate in antiphase with a shear wave in thex direction.
To close, we note that the dispersion relation of the phonons remains unchanged if the driving force is not aligned with one of the principal axes of the crystal, although in that case the form of the eigenmodes is more complex.

Hexagonal lattice
We now consider the case of the hexagonal lattice. The main technical difference with the square lattice is that the reference particle 0 has now six nearest neighbors, see Fig. 3. Repeating the same procedure as above, we compute the coefficients of the stability matrix and obtain where ), Notably, the upper-left 2 × 2 sub-bloc of M is diagonal. As a consequence, an excitation of the position along one direction (x or y) induces no net displacement in the transverse direction. The dispersion relations, plotted in Fig. 6, is given by ) . (21) We see from Eq. 21 that there exist two specific orientations of the wavevectors for which no excitation propagates. For perturbations making angles equal to π/6, (resp. π/2) with the x-axis, we obtain M 1 = 0 (resp. M 1 = M 3 = 0); in these cases, the matrix is not diagonalizable and the only solutions of Eq. 19 is ω = 0. Phonons therefore do not propagate in those two directions.
To illustrate the difference in the dynamics between the hexagonal and the square crystals we consider the behavior in the long wavelength limit. Expanding Eq. 19 at leading order in q, we obtain and the eigenvalues takes the form ω ± = ± µσ 4π |3q x q 2 y − q 3 x |. We infer from this formula that hydrodynamic crystals having an hexagonal symmetry are "softer" than square crystals. When q → 0, the speed of sound goes to 0 as q 2 and even the large-wavelength phonons are dispersive. Furthermore, the sound modes couple the displacements and the orientation of the particles.
To convey a more intuitive picture, we again focus on two specific directions of propagation. We first consider longitudinal perturbations along the first principal axis of the crystal, q = qx. As above, we find that one of the eigenvector correspond to a pure compression along the xaxis, (1, 0, 0). The second eigenvector, (0, −iq/9, 1), mixes shear and orientational waves (bending modes) in quadrature, and depends explicitly on q. Such a coupling was not observed for the square lattice. A second simple case concern the excitations propagating along the second principal direction, namely q = cos(π/3)x + sin(π/3)ŷ. Here, Fig. 7. Normalized dispersion relations for the rectangular, oblique and centered rectangular lattices (ω + only), for 2µσ/πa 2 = 1, and a = 1. Rectangular: = 0.8. Oblique: = 0.9 and β = (π/2) − 0.2. Centered rectangular: = 0.9. the two eigenmodes mix the particle displacements (in only one of the two directions, since x and y cannot couple) and their orientation. They are given by (0, −2iq/9, 1) and (−2iq/(9 √ 3), 0, 1) and correspond to ω − and ω + respectively.

Rectangular,oblique, and face-centered lattices
The lattice geometries corresponding to the rectangular, the oblique and the face-centered lattices are shown in Fig. 3. To derive the eigenmodes we restrict our analysis to calculations with four nearest neighbors, an assumption which restrains the number of crystals for which our calculations are correct (weakly anisotropic and weakly tilted lattices only, as sketched in Fig. 3). It is straightforward to proceed mathematically as in the two previous cases and derive the two sound modes, ω ± , propagating. The results for the dispersion relation are plotted in Fig. 7. In the small q-limit, these phonons always propagate in a dispersive manner. As q goes to zero, the sound velocities reach a constant value which depends on the orientation of the propagation due to the crystal anisotropy.

Response of driven hydrodynamic crystals to finite amplitude perturbations
Before closing this section, we make a final remark regarding the stability of all the five Bravais crystals with respect to finite amplitude perturbations. We start by a simple observation on the relationship between the various lattices. A rectangular lattice corresponds to a square lattice transformed upon a finite homogeneous stretching. An oblique lattice is obtained by stretching and shearing a square lattice. The hexagonal lattice is an oblique lattice with a tilt angle of π/3. Finally, a face-centered lattice is obtained from a rectangular lattice by applying a shear modulated at the highest possible wavelength (q = 2π/a). As all these structures are stationary, and marginally stables at the linear level, we can deduce that any finite amplitude deformation corresponding to an homogeneous shear, or stretch, of the crystal would also be a marginal perturbation: their growth rate would be zero. The same conclusion also holds for rectangular crystals deformed by the specific high-q shear that would transform them into a face-centered lattice.

Hydrodynamic stability of active crystals
We now move on to investigate the linear stability of active swimmer crystals. To do so, we use the same theoretical framework as in the previous section. The swimmers selfpropel along one of the principal axes of the crystals. We also recall that in this active case, the swimming direction is slaved to the particle shape and so is the dipolar flow (Fig. 2). Following the same strategy as in the case of driven particles, we first establish the linearized equations of motion. Combining Eqs. 3, 5 and 7 we obtain and where E s yx is the (y, x) component of the strain rate tensor associated with the dipolar perturbation induced by the swimmer located at j, namely v s (R ij ,p j =x).
Two important remarks can be made at this point. First we see that the stability equations now depend explicitly on the swimming speed of the particles. In addition, contrary to driven lattices, the stability of the swimmer crystals depends on the particle shape through γ. Therefore, we discuss below isotropic and anisotropic swimmers separately.

Isotropic swimmers
Isotropic particles correspond to γ = 0. Their dynamic equation are significantly simplified as Eq. 24 is now trivial and the swimmer orientation remains constant. Since the flow is irrotational, no hydrodynamic torque (from vorticity) is present to modify the orientations of the particles. As in the previous section on driven suspensions, we look for plane waves solutions, from which we infer the form of the stability matrix M. This matrix is analogous to the one defined in Eq. 10 but takes here a slightly different structure where M 1 and M 2 are given by Eq. 11, and Eq. 12 respectively. The two new coefficients are Independently of the crystal symmetry, we see that the eigenvalues of the above matrix are identical to the p Fig. 8. Unstable position/orientation mode for a square lattice of active particles if the volume fraction is high enough (p < 1). The mode is a compression along the ydirection out of phase with a splay perturbation of the particles orientation and leads to the formation of short wavelength bands. one we found for driven crystals, ω 0 = 0, and ω ± = ± M 2 1 + M 2 2 . Therefore, crystals composed of isotropic swimmers are marginally stable and phonons propagate with the same dispersion relations as in driven lattices, albeit with different eigenmodes.

Anisotropic swimmers
We now explore the richer phenomenology arising from swimmer anisotropy. Generic results cannot be established in a framework as general as in the isotropic case. We proceed with the calculation under the nearest neighbors approximation, and deal with the five Bravais lattices separately.

Square lattice
To establish the linear stability of the square crystal we compute all the coefficients of the M matrix using Eqs. 23 and 24 and obtain with where we introduced the dimensionless number Unlike the driven case, the stability matrix for active particles is not characterized solely by the geometry of the lattice. The parameter p quantifies the relative magnitude of the swimming speed and the dipolar advection velocity induced by a neighboring particle. Recall that σ is itself a function of U 0 , and of the particle shape, and σ ≡ BU 0 b 2 , where B is a shape factor of order 1. Therefore, p scales as p ∼ µ −1 (a/b) 2 . Large values of p correspond to the dilute limit, a b, in which our far field approach is expected to be quantitatively correct. Small values of p correspond to a dense crystal, for which our model should capture the essential physical features. The presence of p in the matrix M means that the crystal stability now strongly depends on the particle volume fraction.
As even in the large-p limit the eigenvalues of M take a quite complex form, we proceed to consider the short and long wavelength excitations separately. In the limit q → 0, the matrix M has again three real eigenvalues, corresponding to three propagating modes with frequencies ω 0 = 2γσq x /(πa 2 ), and ω ± = ±2µσq/(πa 2 ). The mode ω 0 is a combination of phonons and orientation waves, whereas ω ± are the phonon modes we found for driven crystals (Fig. 5).
In the high-q limit (small wavelengths) the phenomenology is markedly different. For waves vectors of the edge of the Brillouin zone, q x = 0 and q y = π/a, we find ω 0 = 0 as well as two non-trivial modes, ω ± = ± 2σ πa 3 6γµ(p − 1). Importantly, ω ± are either real or pure imaginary numbers depending on the magnitude of p. In principle, p > 1 for dilute crystals, and therefore the modes ω ± correspond again to a combination of phonons and orientation waves. However, we can expect our results to hold at a qualitative level for more concentrated systems, for which p < 1. In such a case, the hydrodynamic coupling destroys the square crystal structure. Specifically, the ω − mode is unstable. It correspond to the eigenvector (δX, δY, θ) = (0, −i πa 3 6γσ ω − , 1), which combines a compression along the y-axis and splay distortions of the particle orientation. In this strong hydrodynamic coupling limit, the square crystal evolves to form short wavelength bands aligned with the average swimming direction, as sketched in Fig. 8.
At second order in q → 0 and given that p is small enough, the eigenvalues ω ± have a non-zero imaginary part which scales as q 2 . These eigenvalues correspond to the roots of a 3rd-order polynomial, which has no analytical solution. Therefore, we proceed to a numerical investigation of the short wavelength dynamics of the crystal. We compute numerically the eigenvalues of the matrix M for all q's and 0 < p < 3. We find that the square crystals are indeed always unstable for p < 1. In addition the wavenumbers q x = 0 and q y = π/a correspond to the most unstable mode as shown in Fig. 9 for p = 1/2. Fig. 9. Normalized growth rate (iω − ) of the (q x , q y ) modes of the square lattice of active particles for p = 1/2. The parameters are 2µσ/πa 2 = 1, a = 1 and γ/µ = 1.

Hexagonal lattice
When the lattice has hexagonal symmetry, the structure of M is somewhat simplified and we obtain in this case with, Similarly to the square lattice, M does depend on the relative magnitude of the hydrodynamic coupling through p. The general form of the eigenvalues is too complex to yield an intuitive picture. However, the salient features correspond to small wave vectors. In this limit q → 0, the eigenvalues of M are ω 0 = 0, and ω ± = ± σ πa 2 9 2 γµp(q 2 x − q 2 y ). For all values of p, there exists therefore an infinite number of unstable modes growing at a rate |ω ± |. They correspond to perturbations in the position of the particle propagating along a direction making an angle comprised between [±π/4, ±3π/4] with respect to x. This behavior is illustrated in Fig. 10, where we show the variations of the imaginary part of ω ± in the (q x , q y ) plane for p = 5. We note that the most unstable mode is again a combination of compression along the y-axis and splay-like instability of the particle orientation.

Rectangular crystals
The behavior for the rectangular lattices is very similar to what we found for square lattice (within the limits of the Fig. 10. Normalized growth rate (iω − ) of the (q x , q y ) modes of the hexagonal lattice of active particles for p = 5. The parameters are 2µσ/πa 2 = 1, a = 1 and γ/µ = 1.
nearest neighbor approximation). These crystals are all stable for long wavelengths but can display short wavelength instabilities. Denoting the aspect ratio of the lattice cell (Fig. 3) a numerical diagonalization of the stability matrix reveals that again the most unstable mode lies on the edge of the Brillouin zone in the y direction. The associated eigenvalue is Hence, there exists a critical value p c = 1 2 (1 + −2 ) such that the crystal destabilizes for p < p c ; note that p c is a decreasing function of the aspect ratio which plateaus at p = 1/2. Dilute crystals corresponding to high p-values are stable and display phonons modes. Note that, similarly to our observation on the rectangular driven lattices, the latter result implies that dilute swimmer crystals with a rectangular symmetry are marginally stable with respect to finite amplitude stretch deformations.

Oblique and face-centered lattices
The dynamics of active crystals having oblique, or face centered symmetries are much more complex. We here briefly highlight some interesting large scale properties, and comment on the stability of these structures.
The stability matrix of the oblique crystals take a simple for the global modes only, q = 0, yet it reveals an original dynamics. Indeed for q = 0 we get where β is the inclination of the lattice cells and their aspect ratio (Fig. 3). Recall that the nearest neighbor scheme restrains our analysis to weakly tilted and weakly anisotropic lattices. Beyond the ω 0 = 0 mode, the other two eigenvalues are nonzero, and we obtain For weakly tilted lattices β ≈ π/2, so that the frequencies are purely imaginary, and thus q = 0 modes are unstable. Note that this result does not contradict the stationarity of the structure. Indeed, the orientation field is here unstable, thereby inducing a coupled translation of the lattice, as swimmers rotate.
Conversely, in the small-q limit, the face-centered lattices are marginally stable for any amplitude of the hydrodynamic coupling, and phonons and orientation waves propagate in a non dispersive manner. For small wavelengths however, and looking specifically at the combination (q x = 0, q y = π/ a), we see that the eigenmodes corresponding to compression along the y-direction coupled to distortions of the orientation grow exponentially at a rate 2σ π( a) 3 √ 3pγµ. This last results implies that facecentered swimmer lattices are unstable for any amplitude of the hydrodynamic coupling.

Conclusion
In this paper we considered theoretically the dynamics and stability of both driven and active crystals. With a geometry of elongated particles under confinement we derived the dynamical system quantifying the time-evolution of the particle positions and orientations and showed that all five planar Bravais lattices are stationary solutions of the equations of motion. In the case of particles passively driven by an external force we formally demonstrated that all five lattices are always marginally stable. The phonons modes do not depend on the magnitude of the driving force but solely on the orientation and on the symmetries of the lattices. We detailed the spatial structure of the eigenmodes in the square and hexagonal geometry.
In the separate case where the particles are actively self-propelling we showed that the stability of the particle positions and orientations depends not only on the symmetry of the crystals but also on the perturbation wavelengths and the volume fraction of the crystal. We obtained that the square and rectangular lattices are linearly unstable at short wavelengths provided the volume fraction of the crystals is high enough. As a difference, hexagonal, oblique, and face-centered crystals are always unstable.
The results of our work can be compared with past theoretical studies. In the driven case, planar crystalline arrangements, were shown to be hydrodynamically unstable in a three-dimensional fluid at long wavelengths [27]. The results in our paper demonstrate that confinement of the crystals, which algebraically screens hydrodynamic interactions between the particles, leads to a qualitatively different behaviors and all lattices solely support phonon modes.
In the active case, previous work demonstrated the presence of long wavelength instabilities in orientation, density and stress (see [7,8] and references therein). In this past work, aligned suspensions for both pusher and puller swimmers were shown to be unstable in the dilute regime, and so are isotropic suspensions of pushers [10] whereas isotropic puller suspensions, which are linearly stable at zero volume fraction, were shown numerically to be unstable at high volume fraction [28]. In our paper, again because of hydrodynamic screening, the stability characteristics of confined active crystals were found to be independent of the pusher vs. puller nature of the self-propelled particle -the only flow singularity dictating hydrodynamic interactions in this case is the potential flow dipole whose sign is set by the swimming direction only.