Dynamics and instabilities of the free boundary of a two-dimensional dry active nematic aggregate

The dynamics of a two-dimensional aggregate of active rod-shaped particles in the nematic phase with a free boundary is considered theoretically. The aggregate is in contact with a hard boundary at y = 0, a free boundary at y = H(x, t), and in the x-direction the aggregate is of infinite extension. By assuming fast relaxation of the director field, we find instabilities driven by the coupling between the deformation of the free boundary and the active stress in parameter regimes where bulk systems are stable. For a contractile aggregate, when the particles are aligned parallel to the boundaries, the system is unstable in the long wavelengths at any strength of contractility for any H, and the critical wavelength increases as H increases; when the particles are aligned perpendicular to the boundaries, the system acquires a finite-wavelength instability at a critical active stress whose strength decreases as H increases. The behavior for an aggregate with steady-state particle density ρ s , strength of active stress χ, bulk modulus ρ s β, and particles aligned perpendicular to the boundaries can be mapped to one with active stress strength − χ, bulk modulus ρ s (β − χ), and particles aligned parallel to the boundaries. The stability of an extensile aggregate can be obtained from the analysis for contractile aggregates through this mapping as well, even though the corresponding physical mechanisms for the instabilities are different. In the limit H → ∞ , the free boundary is unstable for any contractile or extensile systems in the long-wavelength limit.


Introduction
Active soft matters are known to behave very differently from their passive counterparts or systems driven by external field [1]. This is because the mechanical energy transformed from chemical energy at the scale of the constituent elements produces large-scale complex moving patterns. For example, active nematics and active polar systems easily become unstable because of the positive feedbacks between the active stress and the deformation of the director field [1,2]. As a result, in many actrive nematic systems topological defects are continuously produced and annihilated [3]. These large-scale movements were predicted theoretically in wet active nematic and polar systems [4][5][6], and experimentally verified in systems ranging from aggregates of biological cells [7] to collective motion in motor-filament systems [8].
The soft boundaries of active soft matter systems also show interesting dynamics that cannot be found in boundaries of passive matter. One early example is that membranes of crawling cells exhibit quite peculiar fluctuations due to the dynamics of the underlying active cytoskeleton [9,10]. Theoretically, systems like cell lamellipodium and thin fluid films containing swimming bacteria were modeled as wet active contractile polar systems bounded by a surface in [11]. Interestingly, theoretical studies on these systems predicted instabilities which come from the interplay between surface deformation, elastic and active stress [11][12][13]. The dynamics of tissue-tissue interface due to the coupling of cell proliferation and the mechanical 2. A two-dimensional dry active nematic aggregate with a free boundary Consider a system composed of active non-polar rod-shaped particles on the xy-plane. The passive interaction between these particles makes them aggregate. The resulting aggregate is in the nematic phase, contains N particles, spanning from a rigid wall located at y = 0 to a free boundary at y = H(x, t). That is, the particle density ρ(x, y, t) is high in the region y < H(x, t), and rapidly decays as y crosses H(x, t), at y > H(x, y) the density is negligible. In this work we are interested in the limit when the width of y ≈ H(x, t) region is small compared to all other important length scales, therefore the free boundary is modeled as a sharp line 4 .The extent of the aggregate in the x-direction, L, is much greater than H, and it is convenient to neglect the effect of the boundaries in the xdirection. The director fieldˆ( ) n x y t , , with |ˆ| = n 1 describes the local particle orientation, and the velocity field of this aggregate is denoted by ( )  v x y t , , . In the absence of non-equilibrium activities, the aggregate is in the nematic phase. Schematic cartoons of the system are illustrated in figure 1.
Since the system is active, the steady state and its stability should be determined from the dynamics. First, the system under consideration is compressible and the density of the particles obeys a continuity equation Due to the drag exerted by the substrate, the momentum of the particles is not conserved. When inertia is negligible, the system follows the force balance condition [30] · ( Here s ij p is the passive part of the stress tensor, it can be constructed from the cost of passive interaction free energy due to a change in the particle density [30] 5 .  x v is the drag between the particles and the substrate and ξ is the drag coefficient. The last term on the right-hand side of the equation comes from the active stress of the particles. From symmetry consideration, the active stress s ij a is related to the tensor º á ñ Q nn ij i j of the system [31] ( ) where n i is the ith component ofn, i = x, y, and á ñ n n i j is the average of n i n j over a small region. When χ < 0 the active stress is extensile; when χ > 0 the active stress is contractile. The viscous stress is neglected because the drag force  x v is the dominant dissipative force in the long-wavelength limit. That is, there is no hydrodynamic interaction between the particles. Finally, the evolution of the broken symmetry variablen is described by [1] Figure 1. Schematic cartoons of an aggregate of active non-polar rod-shaped particles on a flat surface (xy plane) with a hard boundary at y = 0 and a free boundary at y = H(x, t). In the absence of nonequilibrium activities, the system is in the nematic phase and H = H s is flat. The dashed lines represent the director field in the aggregate. In (a) the particles align with the upper and lower boundaries, in (b) the particles are perpendicular to these boundaries. 4 The full effective free energy for the passive interaction can be expressed as F = F ρ + F Q + F c , where F ρ depends on particle density only, to the lowest order in spatial gradients it takes the form . f (ρ) has minima at two different densities, one corresponds to the isotropic gaseous phase and the other one is the condensed nematic phase. For example, ( ) , where a has a dimension of area, f r r = Max and χ Flory ? 1 such that the aggregate is in the two-phase regime, far from the critial point. ζ in the density gradient term of F ρ is related to the width of the interface between the gaseous and condensed phases. F c contains terms that couple ρ with the tensor Q ij ≡ 〈n i n j 〉, F Q are terms that depend only on Q ij . The effect of F c and F Q is to make the high density phase nematic. 5 In analogy to model-H phase ordering kinetics, s ij p is related to the density variation through Since we consider a dry system, the only driving force for director dynamics is the elastic forces due to the deformation of the director field, it is related to the derivative of elastic free energy F d , and Γ is the corresponding kinetic coefficient. Since twist deformation is not allowed for two-dimensional nematics, the elastic free energy of the aggregate with one-constant approximation is simply At y = 0 the boundary is not penetrable, therefore the y-component of velocity The drag in the system mainly comes from the substrate and we assume that the y = 0 boundary does not exert extra drag, this means that there is no boundary condition for v x at y = 0. The director field satisfieŝ |ˆ(ˆ) = = n x y y 0 when the particles are aligned parallel (perpendicular) to the boundaries. At y = H(x, t) we have a free boundary. Letˆ( ) N x t , andˆ( ) T x t , be the unit normal and tangential vector of the free boundary, respectively. The director field satisfiesˆ|ˆ(ˆ) ( , and the only nonvanishing component of the active stress is when the particles are aligned parallel (perpendicular) to the boundaries. For both alignment conditions, the condition of zero shear stress at the free boundary is automatically satisfied, and the condition for the normal stress across the free boundary is [33], where P ext is the external (two-dimensional) pressure applied on the free boundary, ò p = 0(1) when the particles are aligned parallel (perpendicular) to the boundaries, σ is the line tension of the free boundary, and  k is the curvature of the free boundary. In this article we are interested in particles which spontaneously form condensed nematic aggregate even when P ext = 0. Finally, there is a kinematic boundary condition v y | y=H = ∂ t H + v x | y=H ∂ x H [34].
With mass continuity equation (1), momentum balance equation (2), relation between the active stress and the nematic order parameter equation (3), dynamics of the director field equation (4), boundary conditions for the velocity field at y = 0, boundary conditions for the stress tensor at y = H, and the kinematic boundary condition at y = H, the dynamics of the system can be solved.
In the steady state, the free surface is flat. When the particles are aligned parallel to the boundaries, the director field isˆ= n x s , and Q xx = 1, all other Q ij ʼs vanish. On the other hand, when the particles are aligned perpendicular to the boundaries, the director field isˆ= n y s , and Q yy = 1, all other Q ij ʼs vanish. Since the steady state is expected to be translationally invariant in the x-direction, it is easy to see from equations (1)-(2) that in the steady state the velocity field  v s vanishes, and the density ρ s is a constant determined by the continuity of normal stress on the free boundary For an aggretate spanning a length L in the x-direction, H s , the steady state magnitude of H, is determined by H s = N/ρ s L.

Linearized dynamics
Suppose a small perturbation is applied to the free boundary at t < 0, the system deviates from its steady-state, and H(x, t) becomes H s + h(x, t). This induces a velocity field ( )  v x y t , , , and the other variables are also perturbed. That is, , , , , , , ,a n d , , . 8

s s s
When the perturbation is turned off at t = 0, these variables evolve with time according to equations (1)-(2)-(4). The steady state is stable in the absence of nonequilibrium activities [35], and we are interested in how the active processes in the aggregate affect the stability of the steady state. For a perturbation with wave vectorqx, the displacement of the free boundary is ( ) The elastic free energy per unit length associated with the deformation of the free boundary is σq 2 |h q (t)| 2 /2, where σ is the energy per unit length of the free boundary. Because of this deformation, the particles in the system do not align perfectly. For aggregates whose particles are oriented parallel (perpendicular) to the surface at H, let the angle betweenn andx (ŷ ) be θ(x, y, t), the elastic free energy of the director field can be written as ( ) [32]. Throughout this article, we assume that the director field relaxes very fast such that at any instance the spatial arrangement of θ minimizes After dropping terms of higher order and terms of higher spatial derivatives, the passive stress s ij p in the linearized dynamics can be expressed as where ( ) s ij p s is the passive stress in the steady state, β is related to the two-dimensional bulk modulus ρβ of the aggregate at ρ = ρ s . Note that by dropping terms of higher spatial derivatives, we neglect the contribution of passive nematic elasticity to the momentum equation. The only contribution of the director field to the dynamics is through the active contractile/extensile stress described by s ij a .
3.1. Particles aligned parallel to the boundaries When the particles align parallel to the boundaries, | q = ¶ = h y H x s , and θ| y=0 = 0. Solving equation (9) with these boundary conditions, one obtains For small θ, the tensor Q ij becomes The remaining variables v, ρ, and h obey the following boundary conditions. First is the kinematic boundary condition that connects the surface dynamics to the velocity field at the free boundary. Next, the boundary condition for the normal stress at y = H connects the change of normal stress across the free boundary and the line tension σ, Furthermore, the normal component of the velocity field vanishes at the rigid boundary y = 0, (14) and keep the leading terms, the continuity equation, momentum equation, and the boundary condition for normal stress at y = H(x, t) become whereQ ij is the traceless part of Q ij , and repeated indices are summed.

Particles aligned perpendicular to the boundaries
When the particles align perpendicular to the boundaries at y = 0 and y = H, let θ(x, y, t) be the angle betweenn andŷ shown in figure 1(b) such that the boundary conditions for θ are also θ| y=H = ∂ x h and θ| y=0 = 0. As a result, equation (11) also describes the director field of a perturbed aggregate but Q ij takes a different form, Now the active stress contributes to the normal stress on the surface at y = H(x, t), (20) and keep the leading terms, we obtain the momentum equation and the boundary condition for the normal stress at It appears that the dynamics of such an aggregate can be obtained by a simple mapping from the dynamics of an aggregate whose particles are aligned parallel to the boundaries. Let equations (21)-(22) become the same as equations (17) and (18) except c¢ and b¢ replace χ and β. For convenience, we will use a dimensionless parameterc c b º to characterize the strength of active stress, under this mapping˜˜˜( Equations (23)- (24) indicate that a perpendicular-aligned contractile (extensile) aggregate with bulk modulus ρ s β is mapped to a parallel-aligned extensile (contractile) aggregate with bulk modulus ρ s (β − χ). Intuitively, a contractile aggregate actively compresses itself along the directorn, and pushes in the direction perpendicular ton; on the other hand an extensile aggregate does the opposite. This is the physical origin of the mathematical relations equations (23)- (24). These relations simplify the analysis a lot: once parallel-aligned aggregates are studied, we can obtain results for perpendicular-aligned aggregates by a simple mapping.

Characteristic length and characteristic time
From equations (16)-(17)-(18), a characteristic length l 0 ≡ σ/ρ s β can be identified as the crossover length from compressibility-driven relaxation at lengths small compared to l 0 to line-tension-driven relaxation at lengths large compared to l 0 . There is a characteristic time scale t x s º l 0 0 3 associated with the relaxation time of the system at wavelength ∼l 0 in the passive (χ = 0) limit.
It is possible to estimate the magnitudes of ł 0 and τ 0 for dry active nematics made of biological cells or artificial active micron-sized particles. For a monolayer of biological cells, line tension of the aggregate σ ∼ γl, where γ ∼ 1 − 10 mN/m [36] is the typical magnitude of the surface tension of biological cells and~ l (μm) is the thickness of the cell monolayer. The two-dimensional bulk modulus ρ s β ∼ Gl, where G ∼ 10 2 − 10 3 Pa is the bulk modulus of a typical cell [37]. Putting these numbers together one finds tht ł 0 = σ/ρ s β ∼ 1 − 100 μ m. Furthermore, the drag coefficient between a cell and the substrate is ξ ∼ 100 Pa s/μm 2 [38] and this gives the characteristic time t x s =  l s 10 0 0 3 5 , a time scale one expected for non-motile cells. On the other hand, for artificial active micron-sized particles, . This leads to l 0 ∼ 1 − 10 μm. The drag coefficient ξ can be estimated by letting the active force per unit area χρ s /l and drag force per unit area ξv to be the same order of magnitude, and letting χ ∼ 10 k B T, l ∼ 1μm for typical particle size and v 1μm/s for characteristic velocity field in the system. From these numbers one finds that ξ ∼ 0.1 − 1k B T s/μm 4  Note that there is another crossover length h x = l w where η is the viscosity of the active aggregate, on lengths shorter than l w the system should be considered as 'wet' [1]. Whether ql 0 1 behavior of our model can be observed in an experiment depends on whether l w < l 0 or l 0 > l w . For example, a simple estimate from [22] shows that in certain monolayer cell sheets one likely has l w > l 0 . On the other hand, two-dimensional aggregates of human melanocytes can be considered as 'dry' [39].

Dispersion relations and instabilities
In this section, the solutions for the linearized equations of motion with the boundary conditions discussed in the previous section are presented. The situation of finite H s is discussed here, some intermediate steps of the calculation are presented in appendix A and appendix B, and some interesting results in the limit H s /l 0 → ∞ are presented in appendix C.

When particles are aligned parallel to the boundaries
For aggregates with director field parallel to the boundaries, θ is determined by equation (11). The density and momentum equations (16)- (17) can be expressed as First, note that the system collapses whenc -< 1 0 unless higher order terms in the passive stress tensor are included. For simplicity we will only consider the situation whenc < 1. Furthermore, the elastic stress of the director field acts as source terms of equations (27)- (28), and the solution contains a homogeneous part and a particular part. That is, The superscript h denotes the homogeneous part of the solution, we present the calculation for this part in appendix A. The superscript p denotes the particular solution. For a perturbation with fixed wave vectorqx and frequency ω, the particular solution can be obtained by substituting . 29 . . In general γ cannot be solved analytically, but it is possible to find the solution graphically and numerically. As shown in appendix B, γ can only be real or pure imaginary. Since equation (A3) indicates that the decay rate for pure imaginary γ is higher than that for real γ, we focus on the slowest mode here and leave the solutions with pure imaginary γ to appendix B. By solving equation (32) numerically for real γ and substitute the resulting γ into equation (A3), the dispersion relation is obtained. Figure 2(a) shows that the decay rate decreases as the system becomes more contractile, and figure 2(b) shows that the decay rate decreases as H s decreases. These results can be understood qualitatively from the following consideration. When particles are aligned parallel to the boundaries, there is no boundary along the direction of active stress. As the line tension of the boundary at y = H tries to bring the system back to the steady state, contractile active stress in the high-density regions tries to recruit more particles. This mechanism is against the effect of line tension, hence the decay rate decreases as the aggregate becomes more contractile. On the other hand, extensile active stress tries to push particles away from high-density regions. This effect helps to restore the steady state and increases the decay rate.
Since one can also see that the system is unstable at long wavelengths, figure 3 shows how H s andc affect the stability of the steady state for contractile aggregates. Asc increases or H s /l 0 decreases, the aggregate becomes less stable as the range of unstable wavelengths increases. This indicates that active contractility and the boundary condition at y = 0 work together to drive the instability. Note that, as discussed in appendix C, the system is unstable at sufficiently small ql 0 even in the H s → ∞ limit.

When particles are aligned perpendicular to the boundaries
Applying the parallel-perpendicular mapping equations (23)- (24) to the numerical solutions of equations (A3)-(32) gives the dispersion relation for systems with director field perpendicular to the boundaries.
As figure 4(a) shows, when ql 0 is not small the relaxation rate increases as the system becomes more contractile. This is due to the same mechanism which we discussed for particles aligned parallel to the boundaries, here use the mapping provided by equations (23)-(24) we obtain different trends comparing to figure 2.
The inset of figure 4(a) shows that when ql 0 becomes sufficiently large, the relaxation rate becomes independent of the strength of active stress. When ql 0 is large, density relaxation is the main driving force that brings the system back to the steady state. Since δρ is periodic in the x-direction, and the active stress for particles aligned perpendicular to the boundaries is in the y-direction, it has a very small effect on density relaxation, thus in this limit the system behaves as a passive one.
Note that, from figure 4(b), for contractile aggregates there is a finite-wavelength instability for strongly contractile aggregates in the ql 0 < 1 region. Since figure 4 shows the dispersion relation for aggregates with only one H s (H s /l 0 = 30). The stability for contractile aggregates is better seen from the phase diagrams shown in    figure 5(a) shows that for given H s /l 0 , a finite-wavelength instability develops at sufficiently largẽ c. However, asc  1 the system becomes stable again. The critical magnitude ofc and the wavenumber of the first unstable mode both decrease as H s /l 0 increases. In fact, in appendix B we show that when H s → ∞ the instability becomes a long-wavelength one. Next, figure 5(b) shows that for givenc > 0, a finite-wavelength instability can be found when H s becomes sufficiently large.

Discussion: mechanisms for the instabilities
The long-wavelength and finite-wavelength instabilities found in the dispersion relations worth further discussion. Since the finite-wavelength instability shown in figure 5 occurs at large H s /l 0 and small ql 0 , and in this region contractile aggregates aligned parallel to the boundaries also has a long-wavelength instability (figure 3), in this section we provide the physical picture of the instabilities in this region.
In the limit χ → 0, the system is passive, since ql 0 = 1, the relaxation is driven by the line tension of the free boundary. From equations (32)-(A3), the relaxation rate scales as This gives us the characteristic time scales due to passive processes. The contribution of active stress to the dynamics of the free boundary can be obtained from the following consideration. First, the momentum equations (27)-(28) for aggregates with director field aligned parallel to the boundaries tell us that for a contractile aggregate, a splay deformation generates an active stress that pushes the particles. For example, a positive ∂ y θ generates an active stress˜cq cH qh H s s that induces a positive v x . Since in our system ∂ y θ is periodic in the x-direction, figure 6(a) shows that the resulting particle movement enhances the deformation of the free boundary and the splay of the director field. On the contrary, a bend deformation generates an active stress that suppresses the deformation. For example, the upward blue arrows in figure 6(a) shows that a positive ∂ x θ generates an active stress˜c q cq qh 2 that induces a positive v y in regions where h is convex. These two active contributions compete with each other, and affect the stability of the system. When ( ) <  qH 1 s the splay-induced destabilizing effect dominates, it produces a horizontal flux gradient that leads to a contributionc ¶ṽ qh H x x s 2 to the growth of h(x, t). From equation (32),the resulting decay rate of the free boundary becomes The first term comes from the decay due to passive processes at ( ) <  qH 1 s , now modified by active stress. The second term comes from the destabilizing drive from active stress induced by splay deformation. This equation clearly shows that the system is unstable in the long wavelength limit, and the range of unstable wavelengths becomes smaller as H s increases. This explains the features in the phase diagrams plotted in figure 3.
On the other hand, for contractile aggregates with director field perpendicular to the boundaries, splay deformation is associated with ∂ x θ ∼ q 2 h and bend deformation is associated with ∂ y θ ∼ qh/H s . Figure 6(b) shows that again splay deformation tends to enhance fluctuation and bend deformation tends to suppress deformation. As a result, one expects that there is no instability when ( ) <  qH 1 s . Indeed, applying parallel-toperpendicular mapping (note that , the decay rate can be derived from equation (32) with the parallelto-perpendicular mapping. By taking ql 0 = 1, the result is Here the first term is the passive decay in ( ) >  qH 1 s region modified by active stress, and the second term is the contribution from splay-induced destabilizing effect. It is clear that the system becomes unstable when the contractility is sufficiently strong, and the instability is at the small-q part of the ( ) >  qH 1 s region. As a result, contractile aggregates aligned perpendicular to the boundaries have a finite-wavelength instability, and the unstable wavelengths are in the region . This is what the phase diagrams in figure 5 show. Notice that as H s increases the region of ( ) >  qH 1 s extends toward smaller q, as appendix C shows eventually the system has a long-wavelength instability as H s → ∞ .
From the mapping rules equations (23)-(24) between extensile and contractile aggregates, extensile aggregates are also unstable. However, the physical mechanisms for the instabilities in extensile aggregates are different. This is because for extensile aggregates splay deformations are suppressed and bend deformations are enhanced by active stress. The phase diagrams for extensile aggregates are not included in section 4.1 and section 4.2 for simplicity, but they can be obtained by numerically applying the parallel/perpendicular mapping from the dispersion relations for contractile aggregates. The resulting phase diagrams look qualitatively similar to the corresponding phase diagrams for contractile aggregates except the dimensionless activity parameterc is mapped to˜(˜) c c --1 . Finally, it is interesting to take our model to the incompressible limit. By expanding equation (32) around 1/β = 0, keeping the leading terms, and taking parallel-to-perpendicular mapping when needed, one obtains  where −(+) is for aggregates whose director field is aligned parallel (perpendicular) to the boundaries. In this case, when ( )   qH 1 s active stress has a negligible contribution to the dispersion relation, and the system is always stable. When ( )   qH 1 s , for aggregates aligned parallel (perpendicular) to the boundaries, only contractile (extensile) systems become unstable in the long wavelength limit.

Summary
A minimal model is presented to study the dynamics and stability of a two-dimensional dry active nematic aggregate. From this model, we identified the characteristic lengths of the system, the mechanisms that drive the system back to the steady state, and the mechanisms that drive the system to various instabilities. There are three characteristic lengths in the system. The first one, 1/q, is the wavelength of the perturbation on the free boundary. The second one, l 0 = σ/ρ s β, characterizes the crossover from compressibility-dominated relaxation (ql 0 1) to relaxation driven by the line tension of the free boundary (ql 0 1). The third characteristic length H s is the steady-state extension of the aggregation in the y-direction, it affects the relative degree of splay and bend induced by the boundary conditions imposed at both boundaries. From these three characteristic lengths we have two dimensionless parameters ql 0 and H s /l 0 . Together with a dimensionless parameterc c b = for the strength of active stress, they determine the behavior of the aggregate. Our analysis shows that the interplay between the active stress and the spatial distribution of the director field further changes the dynamics of the system. Both contractile and extensile aggregates can become unstable because when a perturbation on the free boundary induces a director field deformation, the resulting active stress can push particles to further enhance the displacement of the free boundary. Such boundary-fluctuationinduced director field deformation depends on the boundary condition for the director field and system size, thus we find rich phase diagrams for the stability of the system. Although the splay-induced instability for contractile systems and bend-induced instability for extensile systems have been discussed for bulk systems before [1,3], in our model we work in the limit of fast relaxation of director field and the instabilities cannot happen without deformation of the free boundary, therefore the instabilities we found are resulting from different physical mechanisms.
The theoretical models for the dynamics of active matter boundaries can be classified by the dimensionality of the free boundary. The free surface of a wet active polar/nematic film on a substrate was first considered by Sankararaman and Ramaswamy [11] and later extended in [12]. The director field was assumed to be parallel to the surface and the focus was the dynamics of the director in the plane parallel to the air/liquid surface. For the dynamics of this system in the direction normal to the free surface, calculations [13] show that when the active stress on the free surface is properly included, only when the viscosity of the passive fluid on top of the active film is sufficiently large, contractile aggregates have finite wavelength instabilities and extensile aggregates are always unstable in the long-wavelength limit. These instabilities are driven by the coupling of splay (for contractile systems) and bend (for extensile systems) to the active stress. Because of hydrodynamic interaction and different geometry, the instabilities found in [11][12][13] are quite different from our two-dimensional dry active nematics.
The contact line dynamics of active thin films spreading on a flat surface is more closely related to the system considered in this article. For example, models of cell crawling often treat the cell or cell fragment as a twodimensional active fluid confined by a closed boundary [19,20]. Although drag between the active fluid and the substrate is present, the relative small size of the system makes it necessary for [19,20] to include viscous dissipation. Furthermore, the shape of the boundary is intimately related to the distribution of the director field and the motion of the active fluid domain. As the onset of motion is associated with a change of the cell shape and the director field pattern in the cell, these models form a special class of dynamic free boundary of active matter because the onset of movement is associated with broken rotational symmetry. It is interesting to point out that in another cell crawling model [21], the thickness of the cell is finite, translational invariance is imposed on the direction normal to the crawling direction. This geometry is somewhat closer to [11][12][13] except here the active fluid has contact lines at the front and rear ends of the cell. The most interesting discovery of [21] is that moving state can be found even for tractionless cells. This moving state is similar to the spontaneous flow transition for active polar/nematic fluid [5,6], but here the active fluid has a free surface that is not parallel to the substrate near cell ends.
Spreading cell monolayers, comparing to the aforementioned systems, are most similar to the systems that we considered in this article. The model considered in [17] has a wet incompressible active matter made of selfpropelling polar cells in two-dimensions. They found that within their model there is no fingering instability, and they speculated that fingering only happens when the system is compressible. This is quite interesting, as we also find no finite-wavelength instability when incompressibility is imposed to our dry system. A recent theoretical model [22] is more complicated than [17]. This model has a two-dimensional compressible wet active fluid on a substrate with cell polarity decays as the distance from the free boundary increases. Fingering instability is found in this model, and the result agrees with experimental observations. This work serves as an example of the fact that once details are included, active matter hydrodynamic theory can describe important features of the collective dynamics of real biological systems.
Comparing with the past theoretical works on the dynamics of the free boundaries of active matter, our model is sufficiently simple such that the stability of the free boundary can be analyzed and solved exactly, and the mechanisms for the instabilities can be clearly identified. Our future work is to add more details to our model, and study how the new details change the physics of the results we presented here. For example, when the system becomes unstable, defects should appear and affect the evolution of the system. To describe defect dynamics, one needs to allow the magnitude of nematic order parameter to vary, and it is necessary to solve the full nonlinear model numerically. It will be interesting to see how the new length scale associated with the magnitude of nematic order parameter affect the evolution of our system.

Acknowledgments
LSL acknowledges the hospitality of Prof. Shigeyuki Komura during his visit to TMU in 2019. LSL and HYC are grateful for the financial support from MOST Taiwan (grant number MOST 109-2112-M-008-019-) and support from NCTS Taiwan.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).

Appendix A. Solutions of the homogeneous equations
Since we consider a perturbation with wave vectorqx, we look for solution of the following form, . A1 Here c. c means complex conjugate. k(q, ω) is a complex function to be solved. Substituting equation (A1) into the homogeneous part of equations (26)-(28) leads to a matrix equation , we find the following equations from the real and imaginary parts of equation (32) By tracing each term during the derivation, one can identify the last terms on the right-hand sides of equations (B1)-(B2) as contribution from the particular solution.
To solve a and b, first consider the solutions a 0 , b 0 for equations (B1)-(B2) without the last terms, . Similar to what we did when solving a 0 , b 0 , we first look for solutions of equation (B1) with a = 0 graphically. Figure 9 shows that there are also infinite number of solutions with pure imaginary γ, as there is a one to one correspondence between these solutions and the solutions for b 0 .
The physical meaning of these modes can be understood by noting that when a = 0, k in equation (A1) is k = ± qb. This tells us that each b corresponds to a bulk mode picked up by the boundary conditions. Furthermore, equation (A3) tells us that, for given q, modes with a = 0 (pure imaginary γ) decay faster than modes with b = 0 (real γ).
Since we are interested in the slow dynamics of the system, modes with pure imaginary γ are ignored. From now on we look for solutions with real γ, and it satisfies equation (32).
Appendix C. Large H s limit C.1. Particles aligned parallel to the boundaries When the particles are allowed to change their orientations, for givenc in the large H s limit, in general γ has to be solved numerically from equation (32). After that − iω is solved by inverting equation (A3). The result is plotted in figure 10. Figure 10(a) shows that, like aggregates with fixed particle orientationˆ= n x, the relaxation rate decreases (increases) as the system becomes more contractile (extensile). Furthermore, figures 10(b)(c) show that when H s → ∞ the system becomes unstable in the long-wavelength limit for any contractile and extensile aggregate. The phase diagram in the large-H s limit is shown in figure 11(a).
The phase boundary of this long-wavelength instability can be found analytically. First, use the fact that from equation (A3), when ω = 0˜( Substituting this result back to equation (C1), one finds two critical magnitudes ofc for given ql 0 . The system becomes unstable whenc exceeds these two critical magnitudes,